-1
(평면) 3 체 문제를 해결하기 위해 RK-4 차 순차 방법을 사용하려고합니다. 나는 4 * N 1 차 ODE의 집합을 가지고 있는데, 여기서 N은 객체의 수이다. 4는 각 객체에 대해 {x, y}와 같은 위치 벡터와 {vx, vy}와 같은 속도 벡터가 있기 때문입니다.FORTRAN : dummy 인수
A dummy argument with an explicit INTENT(OUT) declaration is not given an explicit value
내가하는 서브 루틴의 모든 의도 선언을 제거하려고 :
불행하게도이 프로그램은 나에게 내가 이해할 수없는 파생 서브 루틴에 경고를주고, 제대로 작동하지 않습니다 코드에서 문제를 해결하려고 노력하십시오. 그 후 경고는 나타나지 않지만 여전히 프로그램은 계속 잘못 작동합니다.MODULE constant
!----------------------------
! Initial Condition
!----------------------------
REAL*8, PARAMETER :: m0 = 1.9891e30 ! mass sun (kg)
REAL*8, PARAMETER :: m1 = 5.9736e24 ! mass earth (kg)
REAL*8, PARAMETER :: m2 = 1e20 ! mass of test particle (kg)
REAL*8, PARAMETER :: G = 6.6738e-11 ! G (Nm**2/kg**2)
END MODULE constant
PROGRAM bodies_rk
use constant
implicit none
integer :: n_step,n_periods,k_max,k
real*8 :: ome,dt,duepi,m_per,t_per
real*8 :: x0,y0,z0,vx0,vy0,vz0
real*8 :: x1,y1,z1,vx1,vy1,vz1
real*8 :: x2,y2,z2,vx2,vy2,vz2
real*8 :: a_step,h,t,phi,covx,covy,covz,comx,comy,comz
real*8 :: r0,r1,r2
real*8, DIMENSION(12):: x
!--------------------------------
! Test run for Runge-Kutta
!--------------------------------
duepi=8.d0*DATAN(1.d0)
ome=1.d0
t_per=duepi/ome
n_step=1000 ! step per period
n_periods=5000 ! periods
!! position of the Sun !! (in m)
x0 = 0.d0 ; y0 = 0.d0 ; z0 = 0.d0
!! position of the Earth !!
x1 = 1.496e11 ; y1 = 0.d0 ; z1 = 0.d0 !
!! position of the test-particle !!
x2 = 3e12 ; y2 = 0.d0 ; z2 = 0.d0
r0 = sqrt(x0**2.d0+y0**2.d0+z0**2.d0)
r1 = sqrt(x1**2.d0+y1**2.d0+z1**2.d0)
r2 = sqrt(x2**2.d0+y2**2.d0+z2**2.d0)
!! velocity of the Sun !! (in m/s)
vx0 = 0.d0 ; vy0 = 0.d0; vz0 = 0.d0
!! velocity of the Earth !!
vx1 = 0.d0; vy1 = 29.8e3 ; vz1 = 0.d0 !
!! velocity of the test-particle !!
vx2 = 0.d0; vy2 = sqrt((G*(m0+m1+m2))/r2) ; vz2 = 0.d0 !
!! time step size and number of time step !!
a_step=dfloat(n_step)
dt=t_per/a_step
k_max=n_periods*n_step
!! move to center of Mass and velocity !!
comx = (m0*x0+m1*x1+m2*x2)/(m0+m1+m2)
comy = (m0*y0+m1*y1+m2*y2)/(m0+m1+m2)
comz = (m0*z0+m1*z1+m2*z2)/(m0+m1+m2)
covx = (m0*vx0+m1*vx1+m2*vx2)/(m0+m1+m2)
covy = (m0*vy0+m1*vy1+m2*vy2)/(m0+m1+m2)
covz = (m0*vz0+m1*vz1+m2*vz2)/(m0+m1+m2)
x0 = x0 - comx ; x1 = x1 - comx ; x2 = x2 - comx
y0 = y0 - comy ; y1 = y1 - comy ; y2 = y2 - comy
z0 = z0 - comz ; z1 = z1 - comz ; y2 = y2 - comy
vx0 = vx0 - covx ; vx1 = vx1 - covx ; vx2 = vx2 - covx
vy0 = vy0 - covy ; vy1 = vy1 - covy ; vx2 = vx2 - covx
vz0 = vz0 - covz ; vz1 = vz1 - covz ; vx2 = vx2 - covx
!! declare the vector of position for the objects !!
x(1)=x0
x(2)=vx0
x(3)=y0
x(4)=vy0
x(5)=x1
x(6)=vx1
x(7)=y1
x(8)=vy1
x(9)=x2
x(10)=vx2
x(11)=y2
x(12)=vy2
t=0.d0
DO k=1,k_max ! -------------------------->
CALL kutta(x,t,dt)
IF(MOD(k,n_step) == 0) THEN
WRITE(23,*) x(5),x(7) ! test print: orbit of the Earth
END IF
END DO ! <---------------------
END PROGRAM bodies_rk
SUBROUTINE kutta(x,t,dt)
implicit none
REAL*8, INTENT(OUT) :: x(12)
REAL*8, INTENT(IN OUT) :: t
!----------------------------------------------------------
! Subroutine Runge-Kutta 4th Explicit
!----------------------------------------------------------
REAL*8 :: k1(12), k2(12), k3(12), k4(12)
real*8 :: told,dt,h
integer :: j
REAL*8, DIMENSION(12) ::f
told=t
call phi(t,x(1),x(2),x(3),x(4), &
& x(5),x(6),x(7),x(8), &
& x(9),x(10),x(11),x(12), &
& f(1),f(2),f(3),f(4), &
& f(5),f(6),f(7),f(8), &
& f(9),f(10),f(11),f(12))
DO j=1,12
k1(j)=f(j)
END DO
h=0.5d0*dt
call phi(t+h,x(1)+h*k1(1),x(2)+h*k1(2),x(3)+h*k1(3),x(4)+h*k1(4), &
& x(5)+h*k1(5),x(6)+h*k1(6),x(7)+h*k1(7),x(8)+h*k1(8), &
& x(9)+h*k1(9),x(10)+h*k1(10),x(11)+h*k1(11),x(12)+h*k1(12), &
& f(1),f(2),f(3),f(4),f(5),f(6),f(7),f(8),f(9),f(10),f(11),f(12))
DO j=1,12
k2(j)=f(j)
END DO
call phi(t+h,x(1)+h*k2(1),x(2)+h*k2(2),x(3)+h*k2(3),x(4)+h*k2(4), &
& x(5)+h*k2(5),x(6)+h*k2(6),x(7)+h*k2(7),x(8)+h*k2(8), &
& x(9)+h*k2(9),x(10)+h*k2(10),x(11)+h*k2(11),x(12)+h*k2(12), &
& f(1),f(2),f(3),f(4),f(5),f(6),f(7),f(8),f(9),f(10),f(11),f(12))
DO j=1,12
k3(j)=f(j)
END DO
call phi(t+dt,x(1)+dt*k3(1),x(2)+dt*k3(2),x(3)+dt*k3(3),x(4)+dt*k3(4), &
& x(5)+dt*k3(5),x(6)+dt*k3(6),x(7)+dt*k3(7),x(8)+dt*k3(8), &
& x(9)+dt*k3(9),x(10)+dt*k3(10),x(11)+dt*k3(11),x(12)+dt*k3(12), &
& f(1),f(2),f(3),f(4),f(5),f(6),f(7),f(8),f(9),f(10),f(11),f(12))
DO j=1,12
k4(j)=f(j)
END DO
DO j=1,12
x(j)=x(j)+(k1(j)+2.d0*k2(j)+2.d0*k3(j)+k4(j))*dt/6.d0
END DO
t=t+dt
RETURN
END SUBROUTINE kutta
subroutine phi(t,x0,vx0,y0,vy0, &
& fx0,fvx0,fy0,fvy0, &
& x1,vx1,y1,vy1, &
& fx1,fvx1,fy1,fvy1, &
& x2,vx2,y2,vy2, &
& fx2,fvx2,fy2,fvy2)
use constant
implicit none
REAL*8, INTENT(IN OUT) :: t
REAL*8, INTENT(IN) :: x1,x0,x2
REAL*8, INTENT(IN) :: vx1,vx0,vx2
REAL*8, INTENT(IN) :: y1,y0,y2
REAL*8, INTENT(IN) :: vy1,vy0,vy2
REAL*8, INTENT(OUT) :: fx1,fx0,fx2
REAL*8, INTENT(OUT) :: fvx1,fvx0,fvx2
REAL*8, INTENT(OUT) :: fy1,fy0,fy2
REAL*8, INTENT(OUT) :: fvy1,fvy0,fvy2
real*8 :: r0,r1,r2
!----------
! Define the components of the derivative of the vectorial field to solve:
! dx/dt=fx(t,x,vx,y,vy) dvx/dt=fvx(t,x,vx,y,vy)
! dy/dt=fy(t,x,vx,y,vy) dvy/dt=fvy(t,x,vx,y,vy)
!----------
r0=SQRT(((x2-x1)**2.d0+(y2-y1)**2.d0)**3.d0)
r1=SQRT(((x0-x2)**2.d0+(y0-y2)**2.d0)**3.d0)
r2=SQRT(((x1-x0)**2.d0+(y1-y0)**2.d0)**3.d0)
fx0=vx0
fy0=vy0
fvx0 = -G*((m1*(x0-x1))/r2**3.d0)-G*((m2*(x0-x2))/r1)
fvy0 = -G*((m1*(y0-y1))/r2**3.d0)-G*((m2*(y0-y2))/r1)
fx1=vx1
fy1=vy1
fvx0 = -G*((m2*(x1-x0))/r0**3.d0)-G*((m0*(x1-x2))/r2)
fvy0 = -G*((m2*(y1-y0))/r0**3.d0)-G*((m0*(y1-y2))/r2)
fx2=vx2
fy2=vy2
fvx2 = -G*((m0*(x2-x0))/r1**3.d0)-G*((m1*(x2-x1))/r0)
fvy2 = -G*((m0*(y2-y0))/r1**3.d0)-G*((m1*(y2-y1))/r0)
RETURN
END subroutine phi
코드의 버전 난 그냥 서브 루틴 내부의 INTENT을 삭제하지 마십시오에서 : 다음은 코드입니다. 아마도 INTENT 경고는 코드의 잘못된 출력과 관련이 없지만 적어도 문제 중 하나를 해결할 수 있습니다.
고마워요.
당신은 내 대답을 반영하기 위해 질문을 편집함으로써 방금 경고의 원인을 제거했다는 사실을 알고 있습니까? 현재 버전에서는 경고가 더 이상 나타나지 않습니다! –
실제로 그것을 고치기 때문에 정정하십시오. : D 그래서이 게시물을 닫아야합니다. 원하는 경우 모두가 프로그램을 사용할 수 있기 때문에이 문제를 해결합니다. –
글쎄, 그게 어떻게 작동하지 않습니다 ... 정답이 주어지면 게시물은 닫히지 않습니다! 그들이 문을 닫았다면 다른 사람들이 당신의 질문으로부터 어떻게 혜택을받을 수 있습니까? –