fortran code do not cycle

Question

i try to solve this simple three body problem with the following code:

           Program Main
           Implicit real*8 (A-H,O-Z)
           real*8 ome,mu, rho, R
           duepi=8*datan(1.d0)
           ome=1
           mu=0.001


           T_per=duepi/ome


           rho=0.1
           R=1.0
           N_step=100

c Open the file              
             OPEN(unit=11, file="prova1.txt")
c Nested do loops
             do iy0=1,100
               do iP0=1,100
c Calc value for 0
                 y0 = real(iy0)/100.
                 x0 = 0
c Calc value for py0
                 py0 = real(iP0)/100.
                 px0 = 0

           x=x0
           y=y0
           px=px0
           py=py0


           dt=T_per/N_step

           E0=H(x,y,px,py)

           k_max=100*N_step

           k=0
           t=0
           errh=0

c---------
c  start integration loop
c--------
           do k=1,k_max
           call sym4(x,y,px,py,dt)
           E= H(x,y,px,py)
           errh=abs(E-E0)
           t=k*dt

           enddo



           do k=1,k_max
           call sym4(x,y,px,py,-dt)
           E= H(x,y,px,py)
           errh=abs(E-E0)
           t=t-dt

           enddo

        write(11,*) y0, py0, errh 
               enddo ! iP0
             enddo ! iy0
        close(11)

           end



           subroutine sym1(x,y,px,py,dt)
           Implicit real*8 (A-H,O-Z)
c
           call f(x,y,fx,fy)

           pxnew=px+dt*fx
           pynew=py+dt*fy


           xnew=x+dt*pxnew
           ynew=y+dt
c
           x=xnew
           y=ynew
           px=pxnew
           py=pynew

           end

           subroutine sym1_B(x,y,px,py,dt)
           Implicit real*8 (A-H,O-Z)
c

           xnew=x+dt*px
           ynew=y+dt

           call f(xnew,ynew,fxnew,fynew)
           pxnew=px+dt*fxnew
           pynew=py+dt*fynew 
c
           x=xnew
           y=ynew
           px=pxnew
           py=pynew

           end



           subroutine  f(x,y,fx,fy)
           Implicit real*8 (A-H,O-Z)
           real*8 ome,mu,rho,R


            fx = ((1-mu)*(rho+x))/((rho*rho+2*rho*x+y*y)**(1.5)) -
     &           (mu*(R+x))/((R**2-2*R*x+x*x+y*y)**(1.5))
            fy = ((1-mu)*(rho+x))/((rho**2+x**2+2*rho*x+y**2)**(1.5))/
     &       + (mu*y)/((R**2-2*R*x+x**2+y**2)**(1.5))   

           return
           end

           real*8 function H(x,y,px,py)
           Implicit real*8 (A-H,O-Z)
           real*8 ome,mu,rho,R

c          h=px*px/2.d0+ py +(1+eps*cos(ome*y))*x*x/2

c           h=px*px/2.d0+ py -(1+eps*cos(ome*y))*cos(x) 

    r12 = sqrt( ( (x*cos(ome*y)-y*sin(ome*y))+rho*cos(ome*y) )**2
     &        + ( x*sin(ome*y)+y*cos(ome*y) + rho*sin(ome*y) )**2 )
        r13 = sqrt( ( (x*cos(ome*y)-y*sin(ome*y))-R*cos(ome*y) )**2
     &        + ( x*sin(ome*y)+y*cos(ome*y) - R*sin(ome*y) )**2 )       

        h=(px**2+py**2)/2.d0 - (1-mu)/r12 - mu/r13 + py

           return
           end


           subroutine sym2(x,y,px,py,dt)
           Implicit real*8 (A-H,O-Z)

           call f(x,y,fx,fy)


           xnew= x+ px*dt +    fx*dt**2/2.d0   
           ynew= y+ dt                         ! così   è giusto

           call f(xnew,ynew,fxnew,fynew) 
           pxnew= px+ dt*(fx+fxnew )/2.d0
           pynew= py+ dt*(fy+fynew )/2.d0

           x=xnew
           y=ynew
           px=pxnew
           py=pynew

           end


           subroutine sym4(x,y,px,py,dt)
           Implicit real*8 (A-H,O-Z)
           sq2=2**(1.d0/3.d0)
           alpha= 1.d0/(2-sq2)
           beta= sq2/(2-sq2)
           dt1= dt*alpha
           dt2=-dt*beta
           call sym2(x,y,px,py,dt1)
           call sym2(x,y,px,py,dt2)
           call sym2(x,y,px,py,dt1)
           return
       end

the code calls sympletic integrators and solve the 3body problem. But when i try to run it, there aren't compiling errors, the output.txt file show only the initial grid and not the errh this column give me only NaN, can someone help me?

Is maybe an initial condition errors (strange initial conditions for velocities or positions, omega is wrong ...)?

Answer 1

As stated already in the comments, your program has many code style problems (wrong intentation, ...) as well as not using best Fortran practises (e.g. implicit none). Your problem can be solved by trivial usage of a debugger.

An obvious problem is that you are using uninitialized variables in the functions:

1. subroutine f(x,y,fx,fy) : rho, mu, R in calculation of fx and fy, which produces NaN
2. function H() : ome, rho, mu, ...

and similarly on other places