您所在位置：主题：大虾 求调试 求纠正
神奇的UNO
[专家分：0] 发布于 2012-04-21 21:27:00
下面是一段模拟天体运动的代码 （最多4个天体 数据写在一个叫input.be 里面） 为什么总是出现段错误。。。
input.be 里面的数据是：
0.00000000E+01
0.00000000E+01
0.00000000E+01
0.00000000E+01
1.00000000E+01
0.00000000E+01
1.00000000E+01
-1.00000000E+01
0.00000000E+01
3.00000000E-06
0.00000000E+01
5.36000000E+01
-0.42500000E+01
0.00000000E+01
9.55000000E-04
3.47500000E+02
0.00000000E+01
0.00000000E+01
0.02960000E+01
1.00000000E-14
program Bewegung
implicit none
!-----------------------------------------------------
integer :: n_max
integer :: n
integer :: i,j,k
real :: t0,te,delt
REAL, DIMENSION(:,:) , allocatable :: r
REAL, DIMENSION(:,:) , allocatable :: x
REAL, DIMENSION(:,:) , allocatable :: y
REAL(kind = 8), DIMENSION(:,:) , allocatable :: v
REAL(kind = 8), DIMENSION(:,:) , allocatable :: F
REAL(kind = 8), DIMENSION(:,:) , allocatable :: Fx
REAL(kind = 8), DIMENSION(:,:) , allocatable :: Fy
REAL(kind = 8), DIMENSION(:) , allocatable :: m
REAL(kind = 8), DIMENSION(:) , allocatable :: r1
REAL(kind = 8), DIMENSION(:) , allocatable :: r2
REAL(kind = 8) :: t
REAL(kind = 8) :: F1
REAL(kind = 8) :: F2
!-----------------------------------------------------
print*,'Geben Sie bitte N'
read(5,'(I10)')n_max
if ( n_max .le. 1) then
print*,'Mindest 2 Koeper ==> exit'
stop
end if
print*,'Geben Sie bitte Anfang von Zeit t0'
read(5,'(E16.8)') t0
if ( t0 .le. -1d-06) then
print*,'Zeit soll groesser als 0 ==> exit'
stop
end if
print*,'Geben Sie bitte Ende von Zeit te'
read(5,'(E16.8)') te
if ( te .le. t0) then
print*,'Zeit soll groesser als 0 ==> exit'
stop
end if
print*,'Geben Sie bitte Aenderung von Zeit te'
read(5,'(E16.8)') delt
if ( delt .le. 1.0d-06) then
print*,'Zeit soll groesser als 0 ==> exit'
stop
end if
print*,' Delta t = ',delt
!------------------------------------------------------
allocate(r(n_max,2),v(n_max,2), Fx(n_max,n_max), Fy(n_max,n_max), m(n_max))
allocate(x(n_max,n_max),y(n_max, n_max),F(n_max,2),r1(n_max),r2(n_max))
do n=1, n_max
r1(n)=0.0
r2(n)=0.0
end do
open(unit=10, FILE='input.be', STATUS='UNKNOWN')
do n=1, n_max
print*,'Zeit'
do i=1, 2
read(10,'(E16.8)') r(n,i)
end do
do i=1, 2
read(10,'(E16.8)') v(n,i)
END do
read(10,'(E16.8)') m(n)
END do
close(unit=10,status='KEEP')
!-----------------------------------------------------
!CALL initF()
do i=1, n_max-1
print*,'Zeit1'
Fx(i,i)=0.0d+00
Fy(i,i)=0.0d+00
x(i,i)=0.0d+00
y(i,i)=0.0d+00
F(i,1)=0.0d+00
F(i,2)=0.0d+00
print*,'Zeit1'
do j=i+1, n_max
x(i,j)=r(i,1)-r(j,1)
y(i,j)=r(i,2)-r(j,2)
Fx(i,j)=(m(i)*m(j)/((sqrt((x(i,j)**2 + y(i,j)**2)))**3))*x(i,j)
Fy(i,j)=(m(i)*m(j)/((sqrt((x(i,j)**2 + y(i,j)**2)))**3))*y(i,j)
Fx(j,i)=-Fx(i,j)
Fy(j,i)=-Fx(i,j)
end do
print*,'Zeit2'
end do
!----------------------------------------------------
do i=1, n_max
do j=1,n_max
F(i,1)=Fx(i,j)+F(i,1)
F(i,2)=Fy(i,j)+F(i,2)
end do
end do
!-----------------------------------------------------
t=t0
do while ( t .le. te )
do i=1,n_max
!CALL update_r(r,v,F,delt,m,i,n_max)
r(i,1)=r(i,1)+delt*v(i,1)+F(i,1)*delt**2/(2*m(i))
r(i,2)=r(i,2)+delt*v(i,2)+F(i,2)*delt**2/(2*m(i))
!-----------------------------------------------------
F1=F(i,1)
F2=F(i,2)
!CALL update_F(m,Fx,Fy,r,i,n_max)
do j=1, (n_max-1)
print*,"t= ",t
Fx(i,i)=0.0d+00
Fy(i,i)=0.0d+00
x(i,i)=0.0d+00
y(i,i)=0.0d+00
do k=j+1, n_max
x(i,j)=r(i,1)-r(j,1)
y(i,j)=r(i,2)-r(j,2)
Fx(i,j)=(m(i)*m(j)/((sqrt((x(i,j)**2 + y(i,j)**2)))**3))*x(i,j)
Fy(i,j)=(m(i)*m(j)/((sqrt((x(i,j)**2 + y(i,j)**2)))**3))*y(i,j)
Fx(j,i)=-Fx(i,j)
Fy(j,i)=-Fx(i,j)
end do
end do
do j=1, n_max
do k=1,n_max
F(i,1)=Fx(i,j)+F(i,1)
F(i,2)=Fy(i,j)+F(i,2)
end do
end do
!-----------------------------------------------------
!CALL update_v(F,F1,F2,m,v)
v(i,1)=v(i,1) + (F(i,1)+F1)*delt/(2.0d+00*m(i))
v(i,2)=v(i,2) + (F(i,2)+F2)*delt/(2.0d+00*m(i))
print*,' ', t, r(i,1), r(i,2), F(i,1), F(i,2), v(i,1), v(i,2)
print*,'Zeit1'
end do
t=t+delt
end do
1000 format(' ',e14.8,' ',e14.8,' ',e14.8,' ',e14.8,' 'e14.8,' ',e14.8,' ',e14.8)
!-----------------------------------------------------
stop
END PROGRAM Bewegung
input.be 里面的数据是：
0.00000000E+01
0.00000000E+01
0.00000000E+01
0.00000000E+01
1.00000000E+01
0.00000000E+01
1.00000000E+01
-1.00000000E+01
0.00000000E+01
3.00000000E-06
0.00000000E+01
5.36000000E+01
-0.42500000E+01
0.00000000E+01
9.55000000E-04
3.47500000E+02
0.00000000E+01
0.00000000E+01
0.02960000E+01
1.00000000E-14
program Bewegung
implicit none
!-----------------------------------------------------
integer :: n_max
integer :: n
integer :: i,j,k
real :: t0,te,delt
REAL, DIMENSION(:,:) , allocatable :: r
REAL, DIMENSION(:,:) , allocatable :: x
REAL, DIMENSION(:,:) , allocatable :: y
REAL(kind = 8), DIMENSION(:,:) , allocatable :: v
REAL(kind = 8), DIMENSION(:,:) , allocatable :: F
REAL(kind = 8), DIMENSION(:,:) , allocatable :: Fx
REAL(kind = 8), DIMENSION(:,:) , allocatable :: Fy
REAL(kind = 8), DIMENSION(:) , allocatable :: m
REAL(kind = 8), DIMENSION(:) , allocatable :: r1
REAL(kind = 8), DIMENSION(:) , allocatable :: r2
REAL(kind = 8) :: t
REAL(kind = 8) :: F1
REAL(kind = 8) :: F2
!-----------------------------------------------------
print*,'Geben Sie bitte N'
read(5,'(I10)')n_max
if ( n_max .le. 1) then
print*,'Mindest 2 Koeper ==> exit'
stop
end if
print*,'Geben Sie bitte Anfang von Zeit t0'
read(5,'(E16.8)') t0
if ( t0 .le. -1d-06) then
print*,'Zeit soll groesser als 0 ==> exit'
stop
end if
print*,'Geben Sie bitte Ende von Zeit te'
read(5,'(E16.8)') te
if ( te .le. t0) then
print*,'Zeit soll groesser als 0 ==> exit'
stop
end if
print*,'Geben Sie bitte Aenderung von Zeit te'
read(5,'(E16.8)') delt
if ( delt .le. 1.0d-06) then
print*,'Zeit soll groesser als 0 ==> exit'
stop
end if
print*,' Delta t = ',delt
!------------------------------------------------------
allocate(r(n_max,2),v(n_max,2), Fx(n_max,n_max), Fy(n_max,n_max), m(n_max))
allocate(x(n_max,n_max),y(n_max, n_max),F(n_max,2),r1(n_max),r2(n_max))
do n=1, n_max
r1(n)=0.0
r2(n)=0.0
end do
open(unit=10, FILE='input.be', STATUS='UNKNOWN')
do n=1, n_max
print*,'Zeit'
do i=1, 2
read(10,'(E16.8)') r(n,i)
end do
do i=1, 2
read(10,'(E16.8)') v(n,i)
END do
read(10,'(E16.8)') m(n)
END do
close(unit=10,status='KEEP')
!-----------------------------------------------------
!CALL initF()
do i=1, n_max-1
print*,'Zeit1'
Fx(i,i)=0.0d+00
Fy(i,i)=0.0d+00
x(i,i)=0.0d+00
y(i,i)=0.0d+00
F(i,1)=0.0d+00
F(i,2)=0.0d+00
print*,'Zeit1'
do j=i+1, n_max
x(i,j)=r(i,1)-r(j,1)
y(i,j)=r(i,2)-r(j,2)
Fx(i,j)=(m(i)*m(j)/((sqrt((x(i,j)**2 + y(i,j)**2)))**3))*x(i,j)
Fy(i,j)=(m(i)*m(j)/((sqrt((x(i,j)**2 + y(i,j)**2)))**3))*y(i,j)
Fx(j,i)=-Fx(i,j)
Fy(j,i)=-Fx(i,j)
end do
print*,'Zeit2'
end do
!----------------------------------------------------
do i=1, n_max
do j=1,n_max
F(i,1)=Fx(i,j)+F(i,1)
F(i,2)=Fy(i,j)+F(i,2)
end do
end do
!-----------------------------------------------------
t=t0
do while ( t .le. te )
do i=1,n_max
!CALL update_r(r,v,F,delt,m,i,n_max)
r(i,1)=r(i,1)+delt*v(i,1)+F(i,1)*delt**2/(2*m(i))
r(i,2)=r(i,2)+delt*v(i,2)+F(i,2)*delt**2/(2*m(i))
!-----------------------------------------------------
F1=F(i,1)
F2=F(i,2)
!CALL update_F(m,Fx,Fy,r,i,n_max)
do j=1, (n_max-1)
print*,"t= ",t
Fx(i,i)=0.0d+00
Fy(i,i)=0.0d+00
x(i,i)=0.0d+00
y(i,i)=0.0d+00
do k=j+1, n_max
x(i,j)=r(i,1)-r(j,1)
y(i,j)=r(i,2)-r(j,2)
Fx(i,j)=(m(i)*m(j)/((sqrt((x(i,j)**2 + y(i,j)**2)))**3))*x(i,j)
Fy(i,j)=(m(i)*m(j)/((sqrt((x(i,j)**2 + y(i,j)**2)))**3))*y(i,j)
Fx(j,i)=-Fx(i,j)
Fy(j,i)=-Fx(i,j)
end do
end do
do j=1, n_max
do k=1,n_max
F(i,1)=Fx(i,j)+F(i,1)
F(i,2)=Fy(i,j)+F(i,2)
end do
end do
!-----------------------------------------------------
!CALL update_v(F,F1,F2,m,v)
v(i,1)=v(i,1) + (F(i,1)+F1)*delt/(2.0d+00*m(i))
v(i,2)=v(i,2) + (F(i,2)+F2)*delt/(2.0d+00*m(i))
print*,' ', t, r(i,1), r(i,2), F(i,1), F(i,2), v(i,1), v(i,2)
print*,'Zeit1'
end do
t=t+delt
end do
1000 format(' ',e14.8,' ',e14.8,' ',e14.8,' ',e14.8,' 'e14.8,' ',e14.8,' ',e14.8)
!-----------------------------------------------------
stop
END PROGRAM Bewegung