! Metodos Numericos y Simulaciones en Astrofisica
! FCEFN - UNSJ, 2013
!
! calculo de N-cuerpos usando el corrector de hermite 
! avanza las particulas en bloques temporales jerarquicos
!
! el archivo de incio (HERMITE.IN) debe contener:
!
! - tiempo inicial y final (yr)
! - paso temporal (yr)
! - numero de pasos para guardar informacion (int)
! - cte. de amortiguamiento (UA)
! - factor temporal eta 
! - numero de objetos (int)
! - codigo (int) y masa (masa solar) para obj 1
! - coord. X, Y, Z (UA) y vel. VX, VY,VZ (UA/yr) para obj 1
! - idem para obj 2
!   etc.
! ...
! ...
!
! La unidad de distancia es la unidad astronomica
! La unidad de masa es la masa solar
! La unidad de tiempo es el a#o sidereo
!
! El numero maximo de objetos es de 100
!
! se utiliza un COMMON para acelerar la transferencia
! de variables a las subrutinas
!
!
program hermite
implicit none
  integer, parameter ::p=selected_real_kind(16,200)

  integer :: ii,jj,kk,ll
  integer :: ncuer,pasos
  real(p) :: tt0,ttf,tnext,dt,tmin,vk,time
  real(p) :: gg,eps,eta,smax,ene,ss
  integer, dimension(100) :: cod,prox
  real(p), dimension(3) :: ff,ff1,ff2,ff3,off,off1
  real(p), dimension(100) :: masa,dtpar,t0
  real(p),dimension(100,3) :: pos,vel,pos0,vel0,ace,ace1
  real(p), dimension(100,3) :: posx,velx
  
  common pos,vel,masa,ncuer,gg,eps

! parametros iniciales
!
! si quiere cambiar las unidades debe modificar 
! la cte. de gravitacion en la siguiente linea
!
  gg=(0.01720209895_p*365.25_p)**2
  
! lee archivo de entrada
!
  open(7,file='hermite.in')
  read(7,*)tt0,ttf
  read(7,*)dt
  read(7,*)pasos
  read(7,*)eps
  read(7,*)eta
  read(7,*)ncuer

! lee codigos, masas, posiciones y velocidad
! para los NCUER objetos
!
  do ii=1,ncuer
    read(7,*)cod(ii),masa(ii)
    read(7,*)(pos(ii,jj),jj=1,3),(vel(ii,jj),jj=1,3)
  enddo
  close(7)

! calcula los pasos temporales iniciales para cada objeto,
! inicializa el t0, encuentra el paso minimo y calcula
! la fuerza inicial y su derivada
!
  tmin=1.e10_p
  smax=0._p
  do ii=1,ncuer
    call fuerzas1(ii,ff,ff1)
    if(sqrt(ff(1)**2+ff(2)**2+ff(3)**2) /= 0._p)then
      dtpar(ii)=sqrt(eta/sqrt(ff(1)**2+ff(2)**2+ff(3)**2))

! se asegura que los pasos individuales sean menores o iguales
! al ingresado en el archivo de inicio
!
      if(dtpar(ii) > dt) dtpar(ii)=dt
      smax=max(dtpar(ii),smax)
    endif

! guarda valor inicial de la aceleracion y su primer derivada
!
    ace(ii,1:3)=ff
    ace1(ii,1:3)=ff1
  enddo
  write(*,*)'smax: ', smax,' dt: ',dt
  pause

! encuantra para cada objeto el paso correcto en factores de
! 1/2 para facilitar la integracion
!
  do ii=1,ncuer
    vk=2._p+int(log(smax/dtpar(ii))/log(2._p))
    dtpar(ii)=smax/2._p**(vk-1._p)
    tmin=min(tt0+dtpar(ii),tmin)

! en t0 van los tiempos en que fue avanzada la particula
! por ultima vez
!
    t0(ii)=tt0
  enddo

! abre archivo de salida
!
  open(7,file='hermite.out')

  kk=pasos
  pos0=pos
  vel0=vel
  tnext=tt0

! bucle de calculo hasta alcanzar el tiempo 
! de finalizacion
!
  do while(tnext <= ttf)

  write(*,*)tnext,ttf,tmin
  
! decide si escribe en el archivo de salida 
!
    if(kk == pasos)then
      kk=1
      call energia(ene)
      write(7,*)tnext,ncuer,ene
      do ii=1,ncuer
        write(7,*)cod(ii),(pos0(ii,jj),jj=1,3),(vel0(ii,jj),jj=1,3)
      enddo
    else
      kk=kk+1
    endif

    time=0._p
    tnext=tnext+dt
    
! avanza el bloque de particulas que tienen el menor paso
! 
    do while(time < tnext)

! bloque de particulas a avanzar hasta tmin y guarda su indice
! en el array PROX
!
      ll=0
      do ii=1,ncuer
        if(t0(ii)+dtpar(ii) == tmin)then
          ll=ll+1
          prox(ll)=ii
        endif
      enddo
      time=tmin

! predictor de posiciones y velocidades a bajo orden
!
! se asegura que los arrays POS y VEL contengan los
! valores correctos para calcular las nuevas posiciones y
! velocidades
!
      pos=pos0
      vel=vel0
      do ii=1,ncuer

! recupera las aceleraciones y sus primeras derivadas
!
        ff=ace(ii,1:3)
        ff1=ace1(ii,1:3)
        ss=time-t0(ii)

! posiciones y velocidades mediante un desarrollo de Taylor
! alrededor de t0
!
        pos(ii,1:3)=((ff1/6._p*ss+ff/2._p)*ss+vel0(ii,1:3)) &
          *ss+pos0(ii,1:3)
        vel(ii,1:3)=(ff1/2._p*ss+ff)*ss+vel0(ii,1:3)

      enddo

! ahora calcula las aceleraciones y sus primeras derivadas
! utilizando las nuevas posiciones y velocidades obtenidas
! del predictor
!
      tmin=1.e10_p
      posx=pos
      velx=vel
      do jj=1,ll
        ii=prox(jj)
        call fuerzas1(ii,ff,ff1)
        ss=time-t0(ii)
        t0(ii)=time

! calcula segunda y tercer derivada 
!
        off=ace(ii,1:3)
        off1=ace1(ii,1:3)
        ff2=2._p*(-3._p*(off-ff)-(2._p*off1+ff1)*ss)/ss**2
        ff3=6._p*(2._p*(off-ff)+(off1+ff1)*ss)/ss**3

! aplica el corrector de hermite
!
        pos0(ii,1:3)=posx(ii,1:3)+ff2*ss**4/24._p+ff3*ss**5/120._p
        vel0(ii,1:3)=velx(ii,1:3)+ff2*ss**3/6._p+ff3*ss**4/24._p

! guarda la fuerza y su primera derivada
!
        ace(ii,1:3)=ff
        ace1(ii,1:3)=ff1
   
! calcula el nuevo paso temporal para la particula y el nuevo
! paso minimo
!
        if(sqrt(ff(1)**2+ff(2)**2+ff(3)**2) /= 0._p)then
          dtpar(ii)=sqrt(eta/sqrt(ff(1)**2+ff(2)**2+ff(3)**2))
          if(dtpar(ii) > dt) dtpar(ii)=dt
          vk=2._p+int(log(smax/dtpar(ii))/log(2._p))
          dtpar(ii)=smax/2._p**(vk-1._p)
          tmin=min(t0(ii)+dtpar(ii),tmin)
        endif
      enddo
    enddo
  enddo
  close(7)

end program hermite

! subrutina para evaluar la energia total del sistema.
! los valores de entrada ingresan por el COMMON y la
! salida es en la variable ENE.
!
! 
subroutine energia(ene)
implicit none
  integer, parameter ::p=selected_real_kind(16,200)
  integer :: ii,jj
  integer :: ncuer
  real(p) :: gg,rij2,ene,eps
  real(p), dimension(3) :: fue
  real(p), dimension(100) :: masa
  real(p), dimension(100,3) :: pos,vel

  common pos,vel,masa,ncuer,gg,eps

! pone acumulador a cero
!
  ene=0._p

! doble bucle de calculo
!
  do ii=1,ncuer
    do jj=ii+1,ncuer

! vector de posicion relativa de la particula ii respecto jj
!
      fue=pos(ii,1:3)-pos(jj,1:3)
      ene=ene-gg*masa(ii)*masa(jj)/sqrt(fue(1)**2+fue(2)**2+fue(3)**2 &
         +eps**2)
    enddo
    ene=ene+0.5_p*masa(ii)*(vel(ii,1)**2+vel(ii,2)**2+vel(ii,3)**2)
  enddo
end subroutine 

! subrutina para calcular la aceleracion y su primera
! derivada ejercida sobre una particula por las restantes 
! (N-1) particulas del sistema calculada por suma directa
!
! se ingresa con la particula a calcular y la salida se hace
! en los vectores FF y FF1. Las restantes cantidades ingresan
! por el COMMON (posiciones, velocidades, masas, el numero de 
! particulas, el maximo numero de cuerpos, la cte. de gravitacion, 
! y el parametro de amortiguamiento).
!
!
subroutine fuerzas1(ii,ff,ff1)
implicit none
  integer, parameter ::p=selected_real_kind(16,200)
  integer :: ii,jj
  integer :: ncuer
  real(p) :: gg,rij2,rv,eps
  real(p), dimension(3) :: fue,dfue,ff,ff1
  real(p), dimension(100) :: masa
  real(p), dimension(100,3) :: pos,vel

  common pos,vel,masa,ncuer,gg,eps

! pone los acumuladores a cero
!
  ff=0._p
  ff1=0._p

! bucle de calculo
!
  do jj=1,ncuer

    if(ii /= jj)then

! vectores de posicion y velocidad relativa de la
! particula ii respecto jj
!
      fue=pos(ii,1:3)-pos(jj,1:3)
      dfue=vel(ii,1:3)-vel(jj,1:3)

! cuadrado de la distancia entre particulas y producto r*v
!
      rij2=fue(1)**2+fue(2)**2+fue(3)**2+eps**2
      rv=fue(1)*dfue(1)+fue(2)*dfue(2)+fue(3)*dfue(3)

! fuerza y su primera derivada de jj sobre ii
!
      ff=ff-gg*masa(jj)*fue/rij2**(3._p/2._p)
      ff1=ff1-gg*masa(jj)/rij2**(3._p/2._p)*(dfue- &
          3._p*rv*fue/rij2)
    endif
  enddo
end subroutine
