program Source2P

  !use sep
  use source2p_mod

  implicit none
  integer :: ix,iz,iw,ihy,npx,npy,nw,nhx,nsx,nt,nt_pad,nshots
  integer :: ipx,plane_nsx,nx,nhy,nsy,isy,yint,iyint,isx,ihx
  real :: opx,dpx,opy,dpy,ow,dw,pi,ohx,dhx,osx,dsx
  real :: xc(180),zc(180),dt,to,iw_min,iw_max,ohy,dhy
  real :: px,ox,plane_osx,plane_dsx,ex,dx,esx,dsy,osy
  real :: ehx,minx,maxx

  complex, allocatable:: writeout(:,:,:,:),out(:,:),in(:,:,:,:)

  call sep_init()
  call getch_add_string("head=/dev/null")
  call getch_add_string("noheader=y")
  call set_no_putch()

  !! PW parameters
  call from_param("opx",opx,0.)
  call from_param("dpx",dpx,0.001)
  call from_param("npx",npx,1)
  call from_param("yint",yint,1)
  call from_param("nshots",nshots,100)
  !! INPUT Variables
  call from_aux("Shot","o1",ohx)
  call from_aux("Shot","d1",dhx)
  call from_aux("Shot","n1",nhx)
  call from_aux("Shot","n2",nhy)
  call from_aux("Shot","o2",ohy)
  call from_aux("Shot","d2",dhy)
  call from_aux("Shot","o3", ow)
  call from_aux("Shot","d3", dw)
  call from_aux("Shot","n3", nw)
  call from_aux("Shot","o4",osx)
  call from_aux("Shot","d4",dsx)
  call from_aux("Shot","n4",nsx)
  call from_aux("Shot","n5",nsy)
  call from_aux("Shot","o5",osy)
  call from_aux("Shot","d5",dsy)

  nsx=nshots

  write(0,*) '-------------------------'
  write(0,*) 'WELCOME TO S2P'
  write(0,*) ' '
  write(0,*) 'YOUR INPUT AXES: '
  write(0,*) '     AXIS-1: HX',nhx,dhx,ohx
  write(0,*) '     AXIS-2: HY',nhy,dhy,ohy
  write(0,*) '     AXIS-3: OW',nw,dw,ow
  write(0,*) '     AXIS-4: SX',nsx,dsx,osx
  write(0,*) '     AXIS-5: SY',nsy,dsy,osy
  write(0,*) ' '

  !! FIND PW OFFSET and NUMBERS
  plane_osx = osx;  
  plane_dsx = dsx;  
  plane_nsx = nsx
  ox = osx+ohx
  if (ox .ge. osx ) ox = osx-ohx
  dx = dhx

  !! . . Calculate End of the data set
  esx = osx+(nsx-1.)*dsx
  ehx = ohx+(nhx-1.)*dhx

  !! . . See if want to truncate in x
  call from_param("minx",minx,ox )
  call from_param("maxx",maxx,ehx)

  if (minx .gt. ox  ) ox = minx
  if (maxx .lt. ehx ) ehx= maxx

  !! . . Calculate plane-wave line 
  ex  = esx + ehx
  nx  = (ex-ox)/dx + 1.5

  write(0,*) 'Startline:',ox
  write(0,*) 'Num in line:', nx
  write(0,*) 'Endline: ',ex 

  write(0,*) 'YOUR OUTPUT AXES: '
  write(0,*) '     AXIS-1: NX',nx,dx,ox
  write(0,*) '     AXIS-2: NY',nhy*yint,dhy/dble(yint),ohy
  write(0,*) '     AXIS-3: OW',nw,dw,ow
  write(0,*) '     AXIS-4: PX',npx,dpx,opx
  write(0,*) '     AXIS-5: SY',nsy,dsy,osy

  !! OUTPUT VARIABLES
  call auxputch("n1","i",nx,"Plane")
  call auxputch("o1","r",ox,"Plane")
  call auxputch("d1","r",dx,"Plane")

  call auxputch("n2","i",nhy*yint,"Plane")
  call auxputch("o2","r",ohy,"Plane")
  call auxputch("d2","r",dhy/(1.*yint),"Plane")

  call auxputch("n3","i",nw,"Plane")
  call auxputch("o3","r",ow,"Plane")
  call auxputch("d3","r",dw,"Plane")

  call auxputch("n4","i",npx,"Plane")
  call auxputch("o4","r",opx,"Plane")
  call auxputch("d4","r",dpx,"Plane")

  call auxputch("n5","i",nsy,"Plane")
  call auxputch("o5","r",osy,"Plane")
  call auxputch("d5","r",dsy,"Plane")
  call auxputch("esize","i",8,"Plane")
  allocate( out(nx,nw),in(nhx,nhy,nw,nsx),writeout(nx,yint*nhy,nw,npx) )

  SAILLINES: do isy=1,nsy

     !! Read in shot
     call sreed("Shot",in(:,:,:,1),8*nhx*nhy*nw) 

     !! Spread to all shot locations

     do isx=2,nsx
        !$OMP PARALLEL DO PRIVATE(iw,ihy,ihx) SHARED(isx)
        do iw=1,nw
           do ihy=1,nhy
              do ihx=1,nhx
                 in(ihx,ihy,iw,isx) = in(ihx,ihy,iw,1)
              end do
           end do
        end do
        !$OMP END PARALLEL DO
     end do
   
     !! Zero WF
     !$OMP PARALLEL DO PRIVATE(ix,ihy,iw,ipx) 
     do ipx=1,npx
        do iw=1,nw
           do ihy=1,nhy*yint
              do ix=1,nx
                 writeout(ix,ihy,iw,ipx)=0.
              end do
           end do
        end do
     end do
     !$OMP END PARALLEL DO
     
     !! Phase-encode in px direction for all crossline offsets
     COFFSETS: do ihy = 1,nhy
        
        if (maxval(abs(in(:,ihy,:,:))) .lt. 0.000000001) then
        write(0,*) 'Working on X-line cable: ',ihy,' of ', nhy

        else
        write(0,*) 'Working on X-line cable: ',ihy,' of ', nhy
           do ipx = 1,npx

              out=0. 

              px = (ipx-1.)*dpx + opx
              
              call shottoplane1p(out,in(:,ihy,:,:),ox,dx,nx,px,ow,dw,nw,ohx,dhx,nhx,osx,dsx,nsx)
              
              do iw=1,nw
                 do ix=1,nx
                    writeout(ix,(ihy-1)*yint+1,iw,ipx) = out(ix,iw)
                 end do
              end do
           
 !             write(0,*) ipx,ihy,maxval(abs(out))
           end do
           
        end if
     end do COFFSETS
     
     do ipx=1,npx 
        write(0,*) 'WRITING OUT ipx: ',ipx,' of ',npx
        do iw=1,nw
           call srite("Plane",writeout(:,:,iw,ipx),8*nx*nhy*yint)
!           write(0,*) ipx,iw,maxval(abs(writeout(:,:,iw,ipx)))
        end do
!        do iw=1,nw
!              call sep_write(writeout(:,:,iw,ipx),"Plane",nx*nhy*yint,8)
!               call srite("Plane",writeout(:,ihy,iw,ipx),8*nx)
 !       end do
     end do

  end do SAILLINES

  deallocate( out,in,writeout )
     
end program Source2P
