program homework4_1

implicit none
real :: v(2,10001),r(2,10001),kenetic(10001),potential(10001),energy(10001),tstep,sunmass
integer i
call Orbit(v,r,kenetic,potential,energy,tstep,sunmass,i)
call Store(v,r,kenetic,potential,energy,i)
end



!Orbit subroutine 
subroutine Orbit(v,r,kenetic,potential,energy,tstep,sunmass,i)
implicit none
real :: v(2,10001),r(2,10001),kenetic(10001),potential(10001),energy(10001),tstep,sunmass,pi,eccentricity,period,maxr,minr,rad,vel
integer i

tstep=0.01
sunmass=1
pi=acos(-1.0)
maxr=1
minr=1

r(1,1)=1
r(2,1)=0
v(1,1)=0
v(2,1)=4.5
kenetic(1)=0.5*(v(1,1)**2+v(2,1)**2)
potential(1)=-4*pi**2
energy(1)=kenetic(1)+potential(1)



!the iterativeprocess
do i=1,10000
v(1,i+1)=v(1,i)-4*pi**2*sunmass*r(1,i)/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)/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
vel=sqrt(v(1,i+1)**2+v(2,i+1)**2)
rad=sqrt(r(1,i+1)**2+r(2,i+1)**2)
kenetic(i+1)=0.5*vel**2
potential(i+1)=-4*pi**2/rad
energy(i+1)=kenetic(i+1)+potential(i+1)

!to determine the aphelion and perihelion
if(rad>maxr) then
maxr=rad
endif
if(rad<minr) then 
minr=rad
endif
enddo

!calculate the eccentricity and the period
eccentricity=1-2*minr/(maxr+minr)
print *,' eccentricity is ',eccentricity
period=sqrt((maxr/2+minr/2)**3/sunmass)
print *, 'period is ', period 
print *, 'maxr is ', maxr

end


!Store subroutine
subroutine Store(v,r,kenetic,potential,energy,i)
implicit none
real :: v(2,10001),r(2,10001),kenetic(10001),potential(10001),energy(10001)
integer i
open(1,file='orbit.txt')
do i=1,10000
write(1,20) r(1,i),r(2,i),v(1,i),v(2,i),kenetic(i),potential(i),energy(i)
20 format(1x,1p,7(e12.5,2x))
enddo
close(1)
return 
end