! Metodos Numericos y Simulaciones en Astrofisica
! FCEFN - UNSJ, 2010
!
! ejemplo 5d
!
! Metodo de Euler para problemas de valor inicial
!
! problema de los dos cuerpos usando diferencias centrales
!
!
program eje5d
implicit none
  integer, parameter ::p=selected_real_kind(16,200)
  integer :: nro
  real(p) :: pi,hh,time,ee,rad,vv,ene
  real(p),dimension(3) :: vel0,vec0
  real(p),dimension(3) :: vel1,dumm
  real(p),dimension(3) :: vel,vec,dvdt
  
  pi=4._p*atan(1._p)
  
! datos iniciales: objeto en perihelio y a=1
!
  write(*,*)'Ingrese excentricidad: '
  read(*,*)ee
  
  vec=0._p
  vec(1)=1._p*(1._p-ee)
  vel=0._p
  vel(2)=sqrt(2._p/vec(1)-1._p)

! calcula el paso
!
  write(*,*)'Ingrese nro. puntos: '
  read(*,*)nro
  hh=2._p*pi/float(nro)
  
! valor inicial para t
!
  time=0._p
  open(7,file='eje5d.out')
  write(*,*)time
  write(7,*)time,vec,vel
  
! valor inicial para la posicion
!
  rad=sqrt(vec(1)**2+vec(2)**2+vec(3)**2)
  dvdt=-vec/rad**3
  vec0=vec-hh*vel
  vel0=vel-hh*dvdt
  
! inicia bucle de calculo
!
  do while(time < 20._p*pi)
    time=time+hh
	
! calcula valor de posicion
!
    dumm=vec
    vec=2._p*vec-vec0+hh**2*dvdt
	vec0=dumm
	
! calcula valor de la velocidad
!
    dumm=vel
	vel=vel0+2._p*hh*dvdt
    vel0=dumm

    rad=sqrt(vec(1)**2+vec(2)**2+vec(3)**2)
    dvdt=-vec/rad**3	
	
! calcula energia
!
    vv=sqrt(vel(1)**2+vel(2)**2+vel(3)**2)
	ene=0.5_p*vv**2-1._p/rad
    write(*,*)time
    write(7,*)time,vec,vel,ene
  enddo
  close(7)
  
end program eje5d
