PROGRAM three_body_problem_in_rotating_frame_by_4_rungekutta !--変数の宣言をする--------------------------------------------------- implicit none integer :: p,q,i,j real(8) :: x,y,vx,vy,fx,fy real(8) :: xh,yh,vxh,vyh,fxh,fyh real(8) :: xn,yn,vxn,vyn,fxn,fyn real(8) :: x1,y1,x2,y2,r1,r2,xx,yy real(8) :: m1,m2,mu,mu1,mu2 real(8) :: w,t,dt,maxrevo,Cj,Cj0,U,forbid,PI real(8) :: kx(4),ky(4),lx(4),ly(4) real(8) :: w1,w2,Tk1,Tk2 !--書き込むfileの準備------------------------------------------------- open(unit=100, file="4orbit.dat", status="unknown") open(unit=200, file="locate.dat", status="unknown") open(unit=300, file="4jacobi.dat", status="unknown") open(unit=400, file="forbidden.dat", status="unknown") write(400,*) " " !-forbidden.dat の中身のクリア !--定数の値を決める--------------------------------------------------- w = 1.0d0 PI = 3.1415926535898d0 !--中心星と惑星の質量比を求める----------------------------------- m1 = 1.9891d33 m2 = 1.8986d30 mu = m2/(m1+m2) mu1 = 1.0d0-mu mu2 = mu !--中心星と惑星の位置を決める------------------------------------- x1 = -mu2 y1 = 0.0d0 x2 = mu1 y2 = 0.0d0 write(*,*) "factor" write(*,*) "mu1=",mu1,"mu2=",mu2 write(*,*) " x1=",x1," y1=",y1 write(*,*) " x2=",x2," y2=",y2 write(*,*) write(200,*) x1,y1 write(200,*) x2,y2 !--初期条件----------------------------------------------------------- !--粒子の初期位置および初速度を与える----------------------------- x = 1.4960d0 / 5.2026d0 y = 0.0d0 vx = 0.0d0 vy = 1.57d0 write(*,*) "initial condition" write(*,*) " x =",x," y =",y write(*,*) " vx=",vx," vy=",vy write(*,*) !--jacobi定数から定まる粒子の運動の禁止領域を求める------------------- !--jacobi定数の計算----------------------------------------------- r1 = sqrt( (x-x1)**2.0d0 + y**2.0d0 ) r2 = sqrt( (x-x2)**2.0d0 + y**2.0d0 ) Cj = w**2.0d0/2.0d0*(x**2.0d0+y**2.0d0) + mu1/r1 + mu2/r2 - & (vx**2.0d0+vy**2.0d0)/2.0d0 Cj0 = Cj write(*,*) "Cj=",Cj !--xとyを-2から2までloopさせる------------------------------------ do i = 0,200 do j = 0,200 xx = -2.0d0+dble(i)*2.0d-2 yy = -2.0d0+dble(j)*2.0d-2 r1 = sqrt( (xx-x1)**2.0d0 + yy**2.0d0 ) r2 = sqrt( (xx-x2)**2.0d0 + yy**2.0d0 ) !--ある点(x,y)でのポテンシャル U とjacobi定数の差をとる---- U = w**2.0d0/2.0d0*(xx**2.0d0+yy**2.0d0) + mu1/r1 + mu2/r2 forbid = U - Cj !--差が0以下である(粒子の速度が0となる)場合、その点を書き出す- if (forbid .le. 0.0d0) then write(400,*) xx,yy,forbid end if end do !-j end do !-i !--軌道計算----------------------------------------------------------- !--時間が maxrevo となるまで計算する------------------------------ t = 0.0d0 maxrevo = 2.0d0*PI * 1.0d0 do while (t .lt. maxrevo) !--進める時間 dt を決める-------------------------------------- !--現在の位置でそれぞれの天体の回りを円運動するときの周期の-- !--1/1000と1.0d-3を比較して、最も小さいものを dt とする------ dt = 1.0d-3 r1 = sqrt( (x-x1)**2.0d0 + y**2.0d0 ) r2 = sqrt( (x-x2)**2.0d0 + y**2.0d0 ) w1 = sqrt( mu1 / r1**3.0d0 ) Tk1 = 2.0d0*PI/w1 w2 = sqrt( mu2 / r2**3.0d0 ) Tk2 = 2.0d0*PI/w2 dt = min(dt,Tk1/1.0d3,Tk2/1.0d3) !--4 次のrunge-kutta法------------------------------------------ !--現在の位置(x,y,t)での傾き(速度) kx(1),ky(1) を求める------- !--同様に(x,y,vx,vy,t)での速度の傾き(力) lx(1),ly(1) を求める- !--力 fx,fy は最後の Subroutine で計算している---------------- call f(x,y,vx,vy,fx,fy) kx(1) = vx ky(1) = vy lx(1) = fx ly(1) = fy !--(x,y,t)から傾き(1)で時間を dt/2.0 だけ進める--------------- !--進んだ先(xh, yh, t+dt/2.0)での傾き(2)を求める-------------- !--力についても同様に求める----------------------------------- xh = ???????????????? yh = ???????????????? vxh = ???????????????? vyh = ???????????????? call f(xh,yh,vxh,vyh,fxh,fyh) kx(2) = vxh ky(2) = vyh lx(2) = fxh ly(2) = fyh !--(x,y,t)から傾き(2)で時間を dt/2.0 だけ進める--------------- !--進んだ先(xh, yh, t+dt/2.0)での傾き(3)を求める-------------- !--力についても同様に求める----------------------------------- xh = ???????????????? yh = ???????????????? vxh = ???????????????? vyh = ???????????????? call f(xh,yh,vxh,vyh,fxh,fyh) kx(3) = vxh ky(3) = vyh lx(3) = fxh ly(3) = fyh !--(x,y,t)から傾き(3)で時間を dt だけ進める------------------- !--進んだ先(xn, yn, t+dt)での傾き(4)を求める------------------ !--力についても同様に求める----------------------------------- xn = ?????????????? yn = ?????????????? vxn = ?????????????? vyn = ?????????????? call f(xn,yn,vxn,vyn,fxn,fyn) kx(4) = vxn ky(4) = vyn lx(4) = fxn ly(4) = fyn !--ここまでで求めた4つの傾きの平均で時間 dt 進める------------ x = ??????????????????? y = ??????????????????? vx = ??????????????????? vy = ??????????????????? write(100,*) x,y !--jacobi定数の計算--------------------------------------------- r1 = sqrt( (x-x1)**2.0d0 + y**2.0d0 ) r2 = sqrt( (x-x2)**2.0d0 + y**2.0d0 ) Cj = (w**2.0d0)/2.0d0*(x**2.0d0+y**2.0d0) + mu1/r1 + mu2/r2 - & (vx**2.0d0+vy**2.0d0)/2.0d0 write(300,*) t,Cj,abs(Cj-Cj0) t = t + dt end do !-while !--------------------------------------------------------------------- close(unit=100,status="keep") close(unit=200,status="keep") close(unit=300,status="keep") close(unit=400,status="keep") CONTAINS !-力を求めるSubroutine------------------------------------------------ !--本文中で call f(1,2,3,4,answer1,answer2) と書かれたとき、------- !--1,2,3,4 を X,Y,Vx,Vy へ代入して Fx,Fy を計算し、それらを-------- !--answer1,answer2 として返す-------------------------------------- SUBROUTINE f(X,Y,Vx,Vy,Fx,Fy) implicit none real(8) :: X,Y,Vx,Vy,Fx,Fy r1 = sqrt( (X-x1)**2.0d0 + Y**2.0d0 ) r2 = sqrt( (X-x2)**2.0d0 + Y**2.0d0 ) Fx = -mu1*(X+mu2)/(r1**3.0d0) - mu2*(X-mu1)/(r2**3.0d0) + & 2.0d0*w*Vy + w**2.0d0*X Fy = -mu1*Y/(r1**3.0d0) - mu2*Y/(r2**3.0d0) - & 2.0d0*w*Vx + w**2.0d0*Y END SUBROUTINE f !--------------------------------------------------------------------- END PROGRAM