-
program kepler integer i ,j ,k !========반복문 변수 double precision c,d,e,f double precision c1,d1,e1,f1 !===== 배열 보조 변수 double precision :: l1 , l2 double precision cc,dd,ee,ff ,cc1,dd1,ee1,ff1 double precision al,m1,m2,pn,np,ts !==========================입력 변수를 넣는 변수 double precision rl !=========================================== 두 입자사이의 거리 double precision dt,a(10) !================================time scale !========================================================== double precision cx1(10000) double precision cy1(10000) double precision cx2(10000) double precision cy2(10000) !================ coordinate of x,y double precision vx1(10000) double precision vy1(10000) !============================= coordinate of x,y !!================================ veolocity double precision vy2(10000) double precision vx2(10000) !================================= double precision ax2(10000) double precision ax1(10000) double precision ay1(10000) double precision ay2(10000) double precision a2x1(100000) double precision a2y1(100000) double precision a2x2(100000) double precision a2y2(100000) double precision a3x1(10000) double precision a3y1(10000) double precision a3x2(10000) double precision a3y2(10000) double precision l !=========================== cutting time double precision a4x1(10000) double precision a4x2(10000) double precision a4y1(100000) double precision a4y2(100000) !============================ open(unit=3,file='intit.txt',status='old') read(3,*) a(1:7) m1=a(1) m2=a(2) al=a(3) e=a(4) np=a(5) ts=a(6) pn=a(7) close(3) !cx1(1)=-m2/m1*(al*e+al) cx1(1)=1 cx2(1)=-1 l1=al+al*e !cx2(1)=(al+al*e) cy1(1)=0 cy2(1)=0 vx1(1)=0 vx2(1)=0 vy1(1)=-0.555 vy2(1)=0.555 !vy1(1)=m1/m2*((m1+m2)*(2/((1+e)*al)-1/al))**(0.5)/(1+m2/m1) !vy2(1)= -((m1+m2)*(2/((1+e)*al)-1/al))**(0.5)/(1+m2/m1) !===================================이 위만 다시 실행해볼것! !vy1(1)=-(((1-e**2)*(m1*m2)*(m1*m2))*al)**(0.5)/(l1*m1) !vy2(1)=(((1-e**2)*(m1*m2)*(m1*m2))*al)**(0.5)/(l2*m2) dt=0.005 do i=1,9000,1 cc=vx1(i) dd=vx2(i) ee=vy1(i) ff=vy2(i) !======================= c=cx1(i) d=cy1(i) e=cx2(i) f=cy2(i) !============================== rl=abs(((c-e)*(c-e)+(d-f)*(d-f))**(0.5)) ax1(i)=-m2*(c-e)/(rl**3) ax2(i)=-m1*(e-c)/(rl**3) ay1(i)=-m2*(d-f)/(rl**3) ay2(i)=-m1*(f-d)/(rl**3) a2x1(i)=(-m2*((cc-dd)/rl**3)+m2*3*((cc-dd)*(c-e)+(ee-ff)*(d-f))*(c-e)/(rl**5)) a2x2(i)=(-m1*((dd-cc)/rl**3)+m1*3*((dd-cc)*(e-c)+(ee-ff)*(d-f))*(e-c)/(rl**5)) a2y1(i)=(-m2*((ee-ff)/rl**3)+m2*3*((ee-ff)*(d-f)+(cc-dd)*(c-e))*(d-f)/(rl**5)) a2y2(i)=(-m1*((ff-ee)/rl**3)+m1*3*((ff-ee)*(d-f)*(f-d)+(dd-cc)*(c-e)*(f-d))/(rl**5)) !============고침 vx1(i+1)=vx1(i)+ax1(i)*dt+a2x1(i)*(dt**2)/2 vx2(i+1)=vx2(i)+ax2(i)*dt+a2x2(i)*(dt**2)/2 vy1(i+1)=vy1(i)+ay1(i)*dt+a2y1(i)*(dt**2)/2 vy2(i+1)=vy2(i)+ay2(i)*dt+a2y2(i)*(dt**2)/2 cx1(i+1)=cx1(i)+vx1(i)*dt+ax1(i)*(dt**2)/2+a2x1(i)*(dt**3)/6 cx2(i+1)=cx2(i)+vx2(i)*dt+ax2(i)*(dt**2)/2+a2x2(i)*(dt**3)/6 cy1(i+1)=cy1(i)+vy1(i)*dt+ay1(i)*(dt**2)/2+a2y1(i)*(dt**3)/6 cy2(i+1)=cy2(i)+vy2(i)*dt+ay2(i)*(dt**2)/2+a2y2(i)*(dt**3)/6 !===================================================================== cc1=vx1(i+1) dd1=vx2(i+1) ee1=vy1(i+1) ff1=vy2(i+1) !======================= c1=cx1(i+1) d1=cy1(i+1) e1=cx2(i+1) f1=cy2(i+1) !============================== rl=((c1-e1)**2+(d1-f1)**2)**0.5 ax1(1+i)=-m2*(c1-e1)/(rl**3) ax2(i+1)=-m1*(e1-c1)/(rl**3) ay1(i+1)=-m2*(d1-f1)/(rl**3) ay2(i+1)=-m1*(f1-d1)/(rl**3) a2x1(i+1)=(-m2*((cc1-dd1)/rl**3)+m2*3*((cc1-dd1)*(c1-e1)+(ee1-ff1)*(d1-f1))*(c1-e1)/(rl**5)) a2x2(i+1)=(-m1*((dd1-cc1)/rl**3)+m1*3*((dd1-cc1)*(e1-c1)+(ee1-ff1)*(d1-f1))*(e1-c1)/(rl**5)) a2y1(i+1)=(-m2*((ee1-ff1)/rl**3)+m2*3*((ee1-ff1)*(d1-f1)+(cc1-dd1)*(c1-e1))*(d1-f1)/(rl**5)) a2y2(i+1)=(-m1*((ff1-ee1)/rl**3)+m1*3*((ff1-ee1)*(d1-f1)*(f1-d1)+(dd1-cc1)*(c1-e1)*(f1-d1))/(rl**5)) a3x1(i+1)=12*(ax1(i)-ax1(i+1))/dt**3+6*(a2x1(i)+a2x1(i+1))/dt**2 a3x2(i+1)=12*(ax2(i)-ax2(i+1))/dt**3+6*(a2x2(i)+a2x2(i+1))/dt**2 a3y1(i+1)=12*(ay1(i)-ay1(i+1))/dt**3+6*(a2y1(i)+a2y1(i+1))/dt**2 a3y2(i+1)=12*(ay2(i)-ay2(i+1))/dt**3+6*(a2y2(i)+a2y2(i+1))/dt**2 !========================a3 3계 미분!!!!! a4x1(i+1)=-6*(ax1(i)-ax1(i+1))/dt**2-2*(2*a2x1(i)+a2x1(i+1))/dt a4x2(i+1)=-6*(ax2(i)-ax2(i+1))/dt**2-2*(2*a2x2(i)+a2x2(i+1))/dt a4y1(i+1)=-6*(ay1(i)-ay1(i+1))/dt**2-2*(2*a2y1(i)+a2y1(i+1))/dt a4y2(i+1)=-6*(ay2(i)-ay2(i+1))/dt**2-2*(2*a2y2(i)+a2y2(i+1))/dt !==================================================================== a4y a4x 2계 미분 vx1(i+1)=vx1(i+1)+a4x1(i+1)*dt**3/6+a3x1(i+1)*dt**4/24 vx2(i+1)=vx2(i+1)+a4x2(i+1)*dt**3/6+a3x2(i+1)*dt**4/24 vy1(i+1)=vy1(i+1)+a4y1(i+1)*dt**3/6+a3y1(i+1)*dt**4/24 vy2(i+1)=vy2(i+1)+a4y2(i+1)*dt**3/6+a3y2(i+1)*dt**4/24 cx1(i+1)=cx1(i+1)+a4x1(i+1)*dt**4/24+a3x1(i+1)*dt**5/120 cx2(i+1)=cx2(i+1)+a4x2(i+1)*dt**4/24+a3x2(i+1)*dt**5/120 cy1(i+1)=cy1(i+1)+a4y1(i+1)*dt**4/24+a3y1(i+1)*dt**5/120 cy2(i+1)=cy2(i+1)+a4y2(i+1)*dt**4/24+a3y2(i+1)*dt**5/120 !cx1(i+1)=cx1(i)+(c-e)+(cc-dd)*dt+(c-e)/(rl**3)*(dt**2)+(-(cc-dd)/rl**3+3*(cx1(i)-cx2(i))*(c-e)*(c-e)/(rl**5))*(dt**3)/6 !cy1(i+1)=cy1(i)+(d-f)+(ee-ff)*dt+(d-f)/(rl**3)*(dt**2)+(-(cc-dd)/rl**3+3*(cy1(i)-cy2(i))*(d-f)*(d-f)/(rl**5))*(dt**3)/6 !cx2(i+1)=cx2(i)+(c-e)+(dd-cc)*dt+(e-c)/(rl**3)*(dt**2)+3*(-cx1(i)+cx2(i))*(c-e)*(c-e)/(rl**5)*(dt**3)/6 !cy2(i+1)=cy2(i)+(d-f)+(ff-ee)*dt+(f-d)/(rl**3)*(dt**2)+3*(-cy1(i)+cy2(i))*(d-f)*(d-f)/(rl**5)*(dt**3)/6 !vx1(i+1)=vx1(i)+(cc-dd)+(c-e)/(rl**3)*(dt)+(cx1(i)-cx2(i))*(c-e)*(c-e)/(rl**5)*(dt**2)/2 !vx2(i+1)=vx2(i)+(cc-dd)+(d-f)/(rl**3)*(dt)+(-cx1(i)+cx2(i))*(c-e)*(c-e)/(rl**5)*(dt**2)/2 !vy1(i+1)=vy1(i)+(ee-ff)+(e-c)/(rl**3)*(dt)+(cy1(i)-cy2(i))*(d-f)*(d-f)/(rl**5)*(dt**2)/2 !vy2(i+1)=vy2(i)+(ee-ff)+(f-d)/(rl**3)*(dt)-(cy1(i)-cy2(i))*(d-f)*(d-f)/(rl**5)*(dt**2)/2 end do open(unit=3,file='ksp.txt',status='replace') do i=1,8000 write(3,*) cx1(i), cy1(i) end do close(3) end program kepler