PROGRAM planet2 ! d = 2 solar system with major and minor planet DIM x(2),y(2),vx(2),vy(2),ratio(2) LIBRARY "csgraphics" CALL initial(x(),y(),vx(),vy(),t,GM,ratio(),dt,nshow) CALL output(x(),y(),t) LET counter = 0 DO CALL Euler(x(),y(),vx(),vy(),t,GM,ratio(),dt) LET counter = counter + 1 IF mod(counter,nshow) = 0 then CALL output(x(),y(),t) END IF LOOP until key input END SUB initial(x(),y(),vx(),vy(),t,GM,ratio(),dt,nshow) LET t = 0 LET GM = 4.0*pi*pi ! gravitational constant LET dt = 0.001 ! time step (yrs) LET nshow = 20 ! # of iterations between output ! planet one initial conditions LET x(1) = 2.52 LET y(1) = 0 LET vx(1) = 0 LET vy(1) = sqr(GM/x(1)) ! mass of planet 2 is 0.04 mass of sun LET ratio(1) = 0.04*GM ! planet two initial conditions LET x(2) = 5.24 LET y(2) = 0 LET vx(2) = 0 LET vy(2) = sqr(GM/x(2)) ! mass of planet 1 is 0.001 mass of sun LET ratio(2) = -0.001*GM ! set up window CALL compute_aspect_ratio(2*x(2),xwin,ywin) SET WINDOW -xwin,xwin,-ywin,ywin END SUB SUB Euler(x(),y(),vx(),vy(),t,GM,ratio(),dt) ! compute distance between planets 1 and 2 LET dx = x(2) - x(1) LET dy = y(2) - y(1) LET dr2 = dx*dx + dy*dy LET dr = sqr(dr2) LET dr3 = dr2*dr FOR i = 1 to 2 ! sum over planets LET r2 = x(i)*x(i) + y(i)*y(i) LET r = sqr(r2) LET r3 = r2*r LET ax = -GM*x(i)/r3 + ratio(i)*dx/dr3 LET ay = -GM*y(i)/r3 + ratio(i)*dy/dr3 LET vx(i) = vx(i) + ax*dt LET vy(i) = vy(i) + ay*dt LET x(i) = x(i) + vx(i)*dt LET y(i) = y(i) + vy(i)*dt NEXT i LET t = t + dt END SUB SUB output(x(),y(),t) ! plot orbits IF t = 0 then SET COLOR "red" BOX CIRCLE -0.1,0.1,-0.1,0.1 ! sun at origin FLOOD 0,0 END IF SET COLOR "blue" PLOT POINTS: x(1),y(1) ! planet one SET COLOR "green" PLOT POINTS: x(2),y(2) ! planet two END SUB