Program Planet01 c----------------------------------------------------------------- c Motion of a planet around a star c method: Solving a second order ODE with initial conditions c written by: Alex Godunov c last revision: November 2006 c----------------------------------------------------------------- c input from a file (self explanatory) c see file planet1.dat c output ... c to a file named by a user c----------------------------------------------------------------- implicit none real*8 x0, y0, vx0, vy0, dt real*8 t0, tf, r0, v0, E0, Lz0, E, Lz, Torb0 real*8 GM, pi, diff, r, rmax real*8 y(4), relerr, abserr, work(27) real*8 xmin, xmax, ymin, ymax, a, b, ecc integer iflag, iwork(5), ne integer*4 j, jmax character output*12 parameter (pi=3.1415926) parameter (relerr=1.0e-9, abserr=0.0) parameter (jmax=10000) common/const/ GM external planet1 c*** read initial data from a file read 201, output open (unit=7,file=output) read 203, x0 read 203, y0 read 203, vx0 read 203, vy0 read 203, dt read 203, rmax c*** end reading and set initial time to 0.0 t0 = 0.0 c*** calculations in AU units (GM=4*pi**2) GM = 4.0*pi**2 r0 = sqrt(x0**2 + y0**2) v0 = sqrt(GM/r0) Torb0 = 2.0*pi*r0/v0 vx0= vx0*v0 vy0= vy0*v0 c initial energy and angular momentum in E/m units c they must be conserved E0 = 0.5*(vx0**2+vy0**2) - GM/r0 Lz0= x0*vy0 - y0*vx0 write(7,210) write(7,211) t0, x0, y0, vx0, vy0, E0, Lz0 c initial steps to calculate eccentricity c assuptions for calculations are based on initial conditions c initial position (x0,0 ) c initial velocity (0 ,v0) c sun is at (0 ,0 ) xmin = x0 ymax = y0 c rkf45 initial data and conditions for rkf45 and first call c it is very important to call rkf45 for the first time with c iflag = 1 (otherwise the code does not run) ne = 4 iflag = 1 y(1) = x0 y(2) = y0 y(3) = vx0 y(4) = vy0 c*** loop over time till the planet is back in initial point c*** or time is more than 10*T0 (T0 is a period for a circular orbit) j=0 do j=1,jmax tf = t0 + dt call rkf45(planet1,ne,y,t0,tf,relerr,abserr,iflag,work,iwork) if(iflag.eq.7) iflag = 2 E = 0.5*(y(3)**2+y(4)**2) - GM/sqrt(y(1)**2+y(2)**2) Lz= y(1)*y(4) - y(2)*y(3) write(7,211) tf, y(1),y(2),y(3),y(4), E, Lz t0 = tf xmin = min(xmin,y(1)) ymax = max(ymax,y(2)) r = sqrt(y(1)**2+y(2)**2) diff = abs(x0-y(1)) if(diff.le.0.001.or.r.ge.rmax) exit end do a = (x0-xmin)/2.0 b = ymax - y0 ecc = 0.0 if(b.lt.a) ecc = sqrt(1.0-b**2/a**2) write(7,212) ecc 201 format (a12) 202 format (i5) 203 format (f10.4) 210 format(4x,'time',6x,'x',10x,'y',10x,'Vx',9x,'Vy',7x,'E',11x,'Lz') 211 format (f8.2,4f11.3,1pe12.2,1e12.2) 212 format (/,' eccentricity = ',f10.4) pause end subroutine planet1(t, y, yp) c--------------------------------------------- c first and second derivatives for rkf45 c definition of the differential equations c y(1) = x yp(1)=vx=y(3) c y(2) = y yp(2)=vy=y(4) c y(3) = vx yp(3)=d2x/dt2 = - GM*x/r**3 c y(4) = vy yp(4)=d2y/dt2 = - GM*y/r**3 c--------------------------------------------- implicit none Real*8 t, y(4), yp(4), GM, r common/const/ GM yp(1) = y(3) yp(2) = y(4) c equations of motion r = sqrt(y(1)**2+y(2)**2) yp(3) = (-1.0)*GM*y(1)/r**3 yp(4) = (-1.0)*GM*y(2)/r**3 return end