program pw_source
 use sep
 use igrad1

 implicit none

 type axis
    integer 					:: n
    real    					:: o, d
 end type axis

 real, allocatable				:: vel(:,:),velw(:,:),source(:),source2d(:,:)
 real, allocatable				:: dsource(:)
 real						:: xmin,xmax,dx,theta
 integer					:: nx,nwav2d,i,j,stat,k
 integer					:: nwav

 type(axis)					:: wav

 call sep_init()
 call sep_begin_prog()

 call from_aux("source","n1",wav%n); call from_aux("source","d1",wav%d); call from_aux("source","o1",wav%o)

 call from_param("xmin", xmin)
 call from_param("xmax", xmax)
 call from_param("dx", dx)
 call from_param("theta",theta)

 nx = nint((xmax-xmin)/dx)
 theta=30

!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  allocate(dsource(wav%n),source(wav%n))

  call sreed("source", source, 4*wav%n)
  dsource=0
  stat = igrad1_lop (.true.,.false.,dsource,source)
  deallocate(source)
  

  write(0,*) maxval(dsource)

  nwav2d=nint(nx*dx/5.5*sind(theta))
  write(0,*) nwav2d

  call to_history("n1",nwav2d); call to_history("d1",wav%d); call to_history("o1",wav%o) 
  call to_history("n2",nx); call to_history("d2",dx); call to_history("o2",xmin)  
  call to_history("n3",1); call to_history("d3",1); call to_history("o3",1)


  allocate(source2d(nwav2d,nx))

  source2d=0

  do i=1,nx
     j=int(((i-1)*dx/5500*sind(theta))/wav%d)+1
!     write(0,*) j
     do k=1,nwav
        source2d(j+k-1,i)=dsource(k)
     end do
  end do

  call srite("out",source2d,4*nx*nwav2d)

  deallocate(dsource,source2d)

 call sep_end_prog()
end program
