
      
      pi=3.141592
      id1=99
      
      ix1=21
      iy1=11
      iz1=101
      dx=ix1-1
      dy=iy1-1
      dz=iz1-1
      flx=2.0
      fly=5.0
      rr=2.5
      r0=20.0
      open(id1,file="M3DLat_samp04.dat",status="new")
      write(id1,*) 1
      write(id1,*) ix1,iy1,iz1
      do 10 iz=1,iz1
         tz=2.0*pi/(iz1-1)*(iz-1)
      do 10 ix=1,ix1
         x=flx/(ix1-1)*(ix-1)-flx*0.5
      do 10 iy=1,iy1
         y=fly/(iy1-1)*(iy-1)-fly*0.5
         r=sqrt(x*x+y*y)
         if (r.eq.0.0) then
            thi=0.0
         else
            thi=acos(x/r)
            if (y.ne.0.0) thi=thi*abs(y)/y
         end if
         
         x1=r*cos(thi+tz*rr)-r0
         y1=r*sin(thi+tz*rr)
         z1=x1*sin(tz)
         x1=x1*cos(tz)
         p1=z1
         write(id1,*) x1,y1,z1,p1
   10 continue
      
      close(id1)
      stop
      end
      
      