program homework4_2

implicit none
real :: v(2,50001),r(2,50001),rad(50000),tstep,sunmass,theta(50000),t(50000)
integer i
call Orbit(v,r,tstep,sunmass,theta,t,i) !should contain theta and t in order to return the value
call Store(theta,t,i)
end



!Orbit subroutine 
subroutine Orbit(v,r,tstep,sunmass,theta,t,i)
implicit none
real :: v(2,50001),r(2,50001),tstep,sunmass,pi,rad(50000),theta(50000),t(50000),a,minr
integer i,m

m=1
a=0.002
tstep=0.0001
sunmass=1
pi=acos(-1.0)


r(1,1)=0.30966
r(2,1)=0
v(1,1)=0
v(2,1)=12.4577
minr=0.30966
theta=0
t=0

!the iterativeprocess
do i=1,50000
rad(i)=sqrt(r(1,i)**2+r(2,i)**2)
if(rad(i)<minr) then
minr=rad(i)
endif
v(1,i+1)=v(1,i)-4*pi**2*sunmass*r(1,i)*(1+a/rad(i)**2)/sqrt(r(1,i)**2+r(2,i)**2)**3*tstep
v(2,i+1)=v(2,i)-4*pi**2*sunmass*r(2,i)*(1+a/rad(i)**2)/sqrt(r(1,i)**2+r(2,i)**2)**3*tstep
r(1,i+1)=r(1,i)+v(1,i+1)*tstep
r(2,i+1)=r(2,i)+v(2,i+1)*tstep
enddo



!to determine the aphelion and perihelion
do i=1,50000-2,1
if(rad(i+1)>rad(i)) then 
if(rad(i+1)>rad(i+2)) then
if(rad(i+1)>minr) then  !it is necessary to add this condition rad(i+1)>minr
t(m)=i*tstep
theta(m)=180/pi*atan(r(2,i+1)/r(1,i+1))
m=m+1
endif
endif
endif
enddo
return
end


!Store subroutine
subroutine Store(theta,t,i)
implicit none
real :: theta(50000),t(50000)
integer i
open(1,file='precession.txt')
do i=1,20
write(1,20) t(i),theta(i)
20 format(1x,1p,2(e12.5,2x))
enddo
close(1)
return 
end