Fortran Initial condition distribution

Question

i have a program and i want to implement this with a distribution of initial condition (x0=0.1,0.2,...,1 ; and the same for p0). So i want a 10x10 (or nxn) grid of points. I suppose that i must use a double DO cycle:

do x0=0.1,1.0,0.1
 do p0=0.1,1.0,0.1 
 ...
 enddo
enddo

Also i want to read all the different output in a single file (so i have an array of n*3 where n is the number of points).

           Program Main
           Implicit real*8 (A-H,O-Z)
           common ome,eps
           duepi=8*datan(1.d0)
           ome=sqrt(2.d0)

           T_per=duepi/ome


           eps=0.81
           N_step=100

           x0=0.4   
           y0=0
           px0=1.0  
           py0=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 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
           if(mod(k,N_step).eq.0) then
           xex=x0*cos(t)+px0*sin(t)
           pxex=-x0*sin(t)+px0*cos(t)
           yex=y0*cos(t)+py0*sin(t)
           pyex=-y0*sin(t)+py0*cos(t)
           err=sqrt((x-xex)**2+(p-pex)**2+(y-yex)**2+(py-pyex)**2)
           endif 
           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
           if(mod(k,N_step).eq.0) then
           xex=cos(t)+sin(t)
           pex=-sin(t)+cos(t)
           err=sqrt( (x-xex)**2+(p-pex)**2) 
           endif 
           enddo

c--------
c if i put the OPEN-WRITE command here i have the last value of the two do cycle but i
c whant to write EVERY last values for every different initial conditions WITHOUT
c  OVER-WRITE the file.   
c--------

           OPEN(unit=11, file="prova0.txt")
           write(11,'(I3,8g15.6)') k, x0, px0, errh 

           end


           subroutine  f(x,y,fx,fy)
           Implicit real*8 (A-H,O-Z)
           common ome,eps

           fx=-x*(1+eps*cos(ome*y))             
           fy= eps* ome*sin(ome*y)*x*x/2.d0    


           return
           end

           real*8 function H(x,y,px,py)
           Implicit real*8 (A-H,O-Z)
           common ome,eps

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

           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                        

           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

Thanks a lot for your help !

Answer 1

Maybe like this:

1. Open the file before the loops
2. Do the nested loops and write at the end of the loop
3. Close the file

           Program Main
             Implicit real*8 (A-H,O-Z)
             common ome,eps
             duepi=8*datan(1.d0)
             ome=sqrt(2.d0)

             T_per=duepi/ome


             eps=0.81
             N_step=100

c Open the file              
             OPEN(unit=11, file="prova0.txt")
c Nested do loops
             do iX0=1,10
               do iP0=1,10
c                 x0=0.4   
c Calc value for x0
                 x0 = real(iX0)/10.
                 y0=0
c                 px0=1.0  
c Calc value for px0
                 px0 = real(iP0)/10.
                 py0=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 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
                 if(mod(k,N_step).eq.0) then
                 xex=x0*cos(t)+px0*sin(t)
                 pxex=-x0*sin(t)+px0*cos(t)
                 yex=y0*cos(t)+py0*sin(t)
                 pyex=-y0*sin(t)+py0*cos(t)
                 err=sqrt((x-xex)**2+(p-pex)**2+(y-yex)**2+(py-pyex)**2)
                 endif 
                 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
                 if(mod(k,N_step).eq.0) then
                 xex=cos(t)+sin(t)
                 pex=-sin(t)+cos(t)
                 err=sqrt( (x-xex)**2+(p-pex)**2) 
                 endif 
                 enddo

c--------
c Write results
                 write(11,'(I3,8g15.6)') k, x0, px0, errh 
               enddo ! iP0
             enddo ! iX0
c Close file
             close(11)
           end
c ... omitted remaining program 

---

Please note that there are additional issues with your code:

• p is never initialized in

err=sqrt((x-xex)**2+(p-pex)**2+(y-yex)**2+(py-pyex)**2)

and

err=sqrt( (x-xex)**2+(p-pex)**2) 

• The format of your write statement is insufficient: k might be as high as 10000 but you try to represent it with three digits! Try using

write(11,*) k, x0, px0, errh 

• There seem to be quite a number of similar variables like pxex <-> pex. In your question you ask for p0 which actually does not appear in your code. Maybe it would be a good idea to start by using implicit none and manually specifying your variables... BTW: indentation does help ;-)