! Programa para determinar orbitas en base a tres posiciones
! utilizando el metodo de Gauss - Encke - Merton.
!
! El programa asume que las cordenadas heliocentricas de la
! Tierra son eclipticas, pero permite corregir por el baricentro
! Tierra+Luna y por topocentro. 
! Las tres posiciones del objeto pueden ser tanto en ascension 
! recta y declinacion como en longitud  y latitud ecliptica.
! Se asume J(2000.0) y la fecha debe ser expresada en UT.
!
! El archivo de entrada GEM.IN contiene:
!
! CT-UT(seg),flag baricentro(1 0), flag ecuatoriales(1 0)
! fecha(yyyy mm dd.d), coord#1, coord#2, coor. eclip. helio. Tierra, codigo
! ...
!
! donde el flag baricentro indica si se debe corregir del baricentro
! Tierra-Luna al centro de la Tierra, y el flag ecuatoriales indica si 
! coord#1 y coord#2  son ecuatoriales o no.
! Coord#1 y coord#2 pueden ser asc. recta (HH MM SS.S) y declinacion
! (sGG MM SS.S), o long. y lat. eclipticas (GG.GGGGG). El codigo se refiere
! al codigo del lugar de observacion.
!
! El programa guarda la siguiente informacion en las siguientes
! matrices y vectores:
!
! pos(i,3)  = para cada posicion "i", guarda DJ, ar y dec (o long y lat)
! hti(i,3)  = para cada "i", coor. eclipticas heliocentricas de la Tierra 
! rh(i,3)   = para cada "i", componentes del versor de posicion
! tra(3,3)  = en cada columna, los versores xi, eta, y zeta de la transf.
! rho(3,3)  = rh transformada
! rad(i,3)  = para cada "i", componentes vector heliocentrico
! cti(3,3)  = hti transformada
! 
! El programa trabaja en DOBLE PRECISION. Se debe compilar con 
! las subrutinas en CURSO-MEC-1.F90 o posterior
!
! rgh - Julio 2008
! 
program gem
implicit none
  integer, parameter ::d=selected_real_kind(15,200)
  integer ::ibar,iecu,i,j,icont,ifg
  integer ::yy,mm,arh,arm,deg,dem
  real(d), parameter ::gk=0.01720209895_d
  real(d) ::ctut,pi,obl
  real(d) ::dd,ars,des,w1,w2
  real(d) ::t1,t3,c1,c3,ff,fg
  real(d), dimension(3) ::cbar,stacor,top,ybar
  real(d), dimension(3) ::gd,ogd,hd,tt
  real(d),dimension(3,3) ::pos,hti,rh,tra,rho,rad,cti
  real(d),dimension(7) ::coo
  real(d),dimension(8) ::ele
  character(40) ::namobs
  character(3) ::cod,icod
  character ::sn

  pi=4._d*atan(1._d)
  obl=23.43929111_d/180._d*pi

! lee el archivo de entrada
!
  open(7,file='gem.in',status='old')
  read(7,*)ctut,ibar,iecu

! lee observaciones para las tres fechas
!
  do i=1,3
    select case(iecu)

! las posiciones son ecuatoriales
!
      case(1)
       read(7,*)yy,mm,dd,arh,arm,ars,deg,dem,des,(hti(i,j),j=1,3),cod
       pos(i,2)=((ars/60._d+dble(arm))/60._d+dble(arh))/12._d*pi
       pos(i,3)=((des/60._d+dble(dem))/60._d+dble(abs(deg)))* &
         sign(1._d,dble(deg))/180._d*pi

! las posiciones son eclipticas
!
      case(0)
       read(7,*)yy,mm,dd,(pos(i,j),j=2,3),(hti(i,j),j=1,3),cod
       do j=2,3
         pos(i,j)=pos(i,j)/180._d*pi
       enddo
    end select

! calcula dia juliano y corrige por diferencia CT-UT
! (guardo DJ para subrutina topo)
!
    call jdate(yy,mm,dd,pos(i,1))
    dd=pos(i,1)
    pos(i,1)=pos(i,1)+ctut/3600._d/24._d

! encuentra las correcciones a los vectores debido
! al baricentro
!
    cbar=0._d
    if(ibar == 1)call baric(cbar,pos(i,1))

! busca informacion sobre el lugar de observacion
!
    top=0._d
    open(8,file='observ',status='old')
    icod='   '
    do
      read(8,'(a3,1x,f8.4,1x,f7.5,1x,f8.5,1x,a40)',iostat=ifg)icod, &
       stacor,namobs
      if((cod == icod).or.(ifg /= 0))exit
    enddo

    close(8)

! detiene el programa si no hay datos para ese 
! lugar de observacion
!
    if(cod /= icod)then
      write(*,'(/a)')'No hay datos para el observatorio pedido.'
      stop
    else  
      write(*,'(a,i1,a,a3,1x,a40)')'Obs.#',i,' realizada en :',icod,namobs
      stacor(1)=stacor(1)/180._d*pi

! encuentra las correcciones para llevar los vectores
! al topocentro
!
      call topo(top,stacor,dd)
    endif

! suma las correcciones del baricentro y del topocentro a
! las eclipticas heliocentricas de la Tierra, y rota los
! vectores si se ingresaron coordenadas ecuatoriales para 
! las observaciones
!
    hti(i,1:3)=hti(i,1:3)+cbar+top

    if(iecu == 1)call rotvect(hti(i,1:3),-obl,1)    
  enddo
  close(7)

! comienza el calculo de matrices para realizar
! la transformacion geometrica de Cunningham
!

! obtiene las componentes de los versores para las
! tres posiciones
!
  do i=1,3
    rh(i,1)=cos(pos(i,3))*cos(pos(i,2))
    rh(i,2)=cos(pos(i,3))*sin(pos(i,2))
    rh(i,3)=sin(pos(i,3))
  enddo


! calcula los nuevos ejes para la transformacion
!
  tra=0._d

! versor xi en columna 1
!
  tra(1:3,1)=rh(1,1:3)

! versor eta en columna 2
!
  tra(1:3,2)=(rh(3,1:3)-rh(1,1:3)*dot_product(rh(3,1:3),rh(1,1:3))) &
    /sqrt(1._d-(dot_product(rh(3,1:3),rh(1,1:3)))**2)

! versor zeta en columna 3
!
  tra(1,3)=tra(2,1)*tra(3,2)-tra(3,1)*tra(2,2)
  tra(2,3)=tra(3,1)*tra(1,2)-tra(1,1)*tra(3,2)
  tra(3,3)=tra(1,1)*tra(2,2)-tra(2,1)*tra(1,2)

! ahora transforma los versores de las observaciones
!
  rho=matmul(rh,tra)

! control del valor de nu_2
!
  write(*,'(a,g11.4,a)',advance='no')'Valor de nu_2: ',rho(2,3), &
   ' Continua(S/N)?'
  read(*,*)sn
  if((sn /= 'S').and.(sn /='s'))stop

  write(*,'(/a,a)')' #      Fecha 1         Fecha 2         ', &
    'Fecha 3         Ajuste'

! transforma las coordenadas heliocentricas 
! de la Tierra
!
  cti=matmul(hti,tra)

! calcula los intervalos de tiempo y los valores
! iniciales de c1 y c3
!
  do i=1,3
    tt(i)=pos(i,1)
  enddo
  
  c1=(tt(3)-tt(2))/(tt(3)-tt(1))      
  c3=-(tt(1)-tt(2))/(tt(3)-tt(1))

! aqui se inicia la iteracion
!
  gd=0._d
  ogd=0._d
  icont=0

  bucle: do
    t1=gk*(tt(1)-tt(2))
    t3=gk*(tt(3)-tt(2))

!calcula distancias geocentricas
!
    w1=-(cti(2,3)-c1*cti(1,3)-c3*cti(3,3))/rho(2,3)
    w2=(w1-gd(2))**2
    gd(2)=w1
    w1=(gd(2)*rho(2,2)-c1*cti(1,2)+cti(2,2)-c3*cti(3,2))/(c3*rho(3,2))
    w2=w2+(w1-gd(3))**2
    gd(3)=w1
    w1=(gd(2)*rho(2,1)-c3*gd(3)*rho(3,1)-c1*cti(1,1)+cti(2,1)-c3*cti(3,1))/c1
    w2=w2+(w1-gd(1))**2
    gd(1)=w1

! muestra distancias geocentricas
!
    write(*,'(i3,3f16.12,g16.9)')icont,gd,w2

! converge o no?
!
    if((w2 < 1.d-12))then
      write(*,'(/a)',advance='no')'Acepta dist. geocentricas (S/N)?'
      read(*,*)sn
      if((sn == 'S').or.(sn =='s'))exit bucle
      write(*,'(a)',advance='no')'Ingrese dist. geocentricas: '
      read(*,*)gd
      icont=0
      write(*,'(/a,a)')' #      Fecha 1         Fecha 2         ', &
        'Fecha 3         Ajuste'
    else if(icont > 100)then
      write(*,'(/a)')'El calculo no converge.'
      write(*,'(a)',advance='no')'Ingrese dist. geocentricas: '
      read(*,*)gd
      icont=0
      write(*,'(/a,a)')' #      Fecha 1         Fecha 2         ', &
        'Fecha 3         Ajuste'
    endif

! corrige tiempos por aberracion planetaria
!
    do i=1,3
      tt(i)=pos(i,1)-0.005768_d*gd(i)
    enddo

! calcula distancias heliocentricas y componentes
!
    hd=0._d
    do i=1,3
      rad(i,1:3)=gd(i)*rho(i,1:3)+cti(i,1:3)
      do j=1,3
        hd(i)=hd(i)+rad(i,j)**2
      enddo
      hd(i)=sqrt(hd(i))
    enddo

! calcula las relaciones sector/triangulo
!
    w1=sqrt(2._d*(hd(2)*hd(3)+rad(2,1)*rad(3,1)+rad(2,2)*rad(3,2)+ &
      rad(2,3)*rad(3,3)))
    call sectri(hd(2),hd(3),t3,w1,ybar(1))
    w1=sqrt(2._d*(hd(1)*hd(3)+rad(1,1)*rad(3,1)+rad(1,2)*rad(3,2)+ &
      rad(1,3)*rad(3,3)))
    call sectri(hd(1),hd(3),t3-t1,w1,ybar(2))
    w1=sqrt(2._d*(hd(1)*hd(2)+rad(2,1)*rad(1,1)+rad(2,2)*rad(1,2)+ &
      rad(2,3)*rad(1,3)))
    call sectri(hd(1),hd(2),-t1,w1,ybar(3))

! recalcula c1 y c3
!
    c1=ybar(2)/ybar(1)*(tt(3)-tt(2))/(tt(3)-tt(1))
    c3=ybar(2)/ybar(3)*(tt(2)-tt(1))/(tt(3)-tt(1))

    icont=icont+1

  enddo bucle

! encuentra las coordenadas eclipticas geocentricas
! en el sistema no transformado
!
  do i=1,3
   rad(i,1:3)=hti(i,1:3)+gd(i)*rh(i,1:3)
   if(iecu == 1)call rotvect(rad(i,1:3),obl,1)   
  enddo 

! obtiene las funciones F y G
!
  ff=1._d-2._d*(tt(2)-tt(1))**2*gk**2/(w1**2*ybar(3)**2*hd(1))
  fg=(tt(2)-tt(1))/ybar(3)

! calcula las componentes del vector velocidad
!
  do i=5,7
    coo(i)=(rad(2,i-4)-ff*rad(1,i-4))/fg
  enddo

  coo(2:4)=rad(1,1:3)
  coo(1)=tt(1)

! obtiene elementos osculadores
!
  call vectoe(coo,ele,gk**2)

  do i=3,7
    ele(i)=ele(i)*180._d/pi
  enddo  
  
! salida de resultados
!
  write(*,'(/a)')'Elem. Osculadores J(2000.0):'
  write(*,'(a,f9.5)')'a= ',ele(1)
  write(*,'(a,f9.6)')'e= ',ele(2)
  write(*,'(a,f9.5)')'i= ',ele(3)
  write(*,'(a,f9.5)')'W= ',ele(4)
  write(*,'(a,f9.5)')'w= ',ele(5)
  write(*,'(a,f9.5)')'M= ',ele(6)
  write(*,'(a,f9.7)')'n= ',ele(7)
  write(*,'(a,f15.7)')'T= ',ele(8)

end program gem