! Metodos Numericos y Simulaciones en Astrofisica
! FCEFN - UNSJ, 2013
!
! calculo de N-cuerpos por leapfrog
! fuerzas por suma directa y paso fijo 
!
! el archivo de incio (NCUERPOS0.IN) debe contener:
!
! - tiempo inicial y final (yr)
! - paso temporal (yr)
! - numero de pasos para guardar informacion (int)
! - 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 viene dado por 
! el parametro CMAX.
!
program ncuerpos0
implicit none
  integer, parameter ::p=selected_real_kind(16,200)
  integer, parameter :: cmax=50

  integer :: ii,jj,kk
  integer :: ncuer,pasos
  real(p) :: tt0,ttf,dt,gg
  integer, dimension(cmax) :: cod
  real(p), dimension(cmax) :: masa
  real(p),dimension(cmax,3) :: pos,vel,vel0,ace
  
! parametros iniciales
!
  gg=(0.01720209895_p*365.25_p)**2

! lee archivo de entrada
!
  open(7,file='ncuerpos0.in')
  read(7,*)tt0,ttf
  read(7,*)dt
  read(7,*)pasos
  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)

! abre archivo de salida y
! guarda posicion y velocidad inicial
!
  open(7,file='ncuerpos0.out')
  write(7,*)tt0,ncuer
  do ii=1,ncuer
    write(7,*)cod(ii),(pos(ii,jj),jj=1,3),(vel(ii,jj),jj=1,3)
  enddo

! avanza medio paso temporal en velocidad para
! iniciar el proceso de calculo por leapfrog
! conservando precision de segundo orden
!
  call fuerzas(pos,masa,ncuer,cmax,gg,ace)
  vel=vel+ace*dt/2._p

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

! calcula posicion en t+dt
!
    pos=pos+dt*vel

! calcula aceleraciones usando las nuevas
! posiciones en t+dt
!
    call fuerzas(pos,masa,ncuer,cmax,gg,ace)

! calcula velocidades en t+3*dt/2
!
    vel0=vel
    vel=vel+ace*dt

! incrementa el tiempo
!
    tt0=tt0+dt

! tiene que escribir en el archivo de salida?
!
    if(kk == pasos)then
      kk=1

! calcula velocidad centrada en t
!
      vel0=(vel+vel0)/2._p

! escribe en el archivo
!
      write(7,*)tt0,ncuer
      write(*,*)tt0,ncuer
      do ii=1,ncuer
        write(7,*)cod(ii),(pos(ii,jj),jj=1,3),(vel0(ii,jj),jj=1,3)
      enddo
    else
      kk=kk+1
    endif

  enddo
  close(7)

end program ncuerpos0

! subrutina para calcular la aceleracion ejercida sobre 
! una particula por las restantes (N-1) particulas del sistema
! calculada por suma directa
!
!
subroutine fuerzas(pos,masa,ncuer,cmax,gg,ace)
implicit none
  integer, parameter ::p=selected_real_kind(16,200)
  integer :: ii,jj
  integer :: cmax,ncuer
  real(p) :: gg
  real(p), dimension(3) :: cooi,cooj,fue
  real(p), dimension(cmax) :: masa
  real(p), dimension(cmax,3) :: pos,ace

! pone los acumuladores a cero
!
  ace=0._p

! doble bucle de calculo
!
  do ii=1,ncuer
    cooi=pos(ii,1:3)
    do jj=ii+1,ncuer
      cooj=pos(jj,1:3)

! vector de posicion relativa de ii respecto de jj
!
      fue=cooi-cooj

! termino de la sumatoria de fuerzas
!
      fue=fue/(fue(1)**2+fue(2)**2+fue(3)**2)**(3._p/2._p)

! fuerza de jj sobre ii
!
      ace(ii,1:3)=ace(ii,1:3)-gg*masa(jj)*fue

! fuerza de ii sobre jj
!
      ace(jj,1:3)=ace(jj,1:3)+gg*masa(ii)*fue
    enddo
  enddo

end subroutine
