Fortran beginner - выходной файл показывает NaN

Question

Я написал простую программу Fortran для трех задач тела, используя алгоритм Эйлера-Ричардсона. По какой-то причине выходные файлы не дают ничего, кроме NaN. Будет ли проблема решена с помощью подпрограмм или функций?

PROGRAM Threebody 
IMPLICIT none 



!************** Variable declarations *************** 

!real(8), parameter :: G = 6.6738D-11 

INTEGER :: IOSTAT, dt, t_step, io_error, i 
DOUBLE PRECISION:: st, g, m0, m1, m2, force0x, force0y, force1x, force1y, force2x, force2y, n0, n1, n2, x0,v0,x1,v1,x2,v2,y0,w0,y1 
DOUBLE PRECISION:: w1,y2,w2 

m0 = 1.d0        ! the masses of the three bodies 
m1 = 1.d0 
m2 = 1.d0 

g = 9.80d0        ! m/s^2 s due to gravity 

t_step = 1 

y0 = 0.d0       
x0 = 0.d0 
v0 = 0.d0        
w0 = 0.d0 

y1 = 0.d0       
x1 = 1.d0 
v1 = 2.d0        
w1 = 1.118d0 

y2 = 0.d0       
x2 = -1.d0 
v2 = -1.118d0        
w2 = 0.d0 


dt = 1 ! time step 

OPEN(unit=5, file='out.txt', status='replace',action='write', IOSTAT=io_error) 
OPEN(unit=6, file='out1.txt', status='replace',action='write', IOSTAT=io_error) 
OPEN(unit=7, file='out2.txt', status='replace',action='write', IOSTAT=io_error) 





! particle no.1 
DO i=0,1000,t_step       ! 1000 is the total time 

st=i!/10.d0 

n0 = sqrt(((x2 - x1)**2 + (y2 - y1)**2)**3) 

force0x = (- m1*((x0 - x1)/n2)) - (m2*((x0 - x2)/n1)) 

x0 = x0 + (w0 + 0.5*force0x*st)*st 
w0 = w0 + force0x*st 

force0y = (- m1 * ((y0 - y1)/n2)) - (m2 * ((y0 - y2)/n1)) 

y0 = y0 + (v0 + 0.5*force0y*st)*st 
v0 = v0 + force0y*st 


WRITE(5,*) i*t_step, x0, y0, n0 

END DO 




! particle no.2 
DO i=0,1000,dt    

st=i 

n1 = sqrt(((x0 - x2)**2 + (y0 - y2)**2)**3) 

force1x = (- m2 * ((x1 - x0)/(n0))) - m0 * ((x1 - x2)/n2) 

x1 = x1 + (w1 + 0.5*force1x*st)*st 
w1 = w1 + force1x*st 

force1y = (- m2 * (y1 - y0)/n0) - (m0 * ((y1 - y2)/n2)) 

y1 = y1 + (v1 + 0.5*force1y*st)*st 
v1 = v1 + force1y*st 


WRITE(6,*) i*t_step, x1, y1 
END DO 




! particle no.3 
DO i=0,1000,dt 

st=i    

n2 = sqrt (((x1 - x0)**2 + (y1 - y0)**2)**3) 

force2x = (- m0 * ((x2 - x0)/n1)) - (m1 * ((x2 - x1)/n0)) 

x2 = x2 + (w2 + 0.5*force2x*st)*st 
w2 = w2 + force2x*st 

force2y = (- m0 * ((y2 - y0)/n1)) - (m1 * ((y2 - y1)/n0)) 

y2 = y2 + (v2 + 0.5*force2y*st)*st 
v2 = v2 + force2y*st 


WRITE(7,*) i*t_step, x2, y2, n2 
END DO 


CLOSE(unit=5) 
CLOSE(unit=6) 
CLOSE(unit=7) 

END PROGRAM Threebody 

Answer 1

Компиляция с -Wall в gfortran дает подсказку:

f.f90:56:0: warning: ‘n1’ may be used uninitialized in this function [-Wmaybe-uninitialized] 
     force0x = (- m1*((x0 - x1)/n2)) - (m2*((x0 - x2)/n1)) 

n1 и n2 не инициализируются, таким образом они могут содержать что-либо. В моем случае 6,95e-310, вызывая force0x, дают NaN.

Чтобы отладить это, вы можете распечатать переменные на экране, чтобы проверить, когда они становятся NaN. Для этого не используйте единицы ниже 10, потому что они могут быть зарезервированы. Я думаю, что Unit 5 является стандартным выходом, поэтому, повторно используя его, вы ничего не можете распечатать на экране.