program source2d
 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
 integer					:: nhpad,nto,nwav2

 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)

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

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

  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)

  nwav2d=nint(nx*dx/5.5)

  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)

  nx=501;theta=30

  allocate(source2d(nwav2d,nx))

  do i=1,nx
     j=nint(((i-1)*dx/5500*sind(theta))/wav%d)
     source2d(j:j+nwav-1,i)=dsource(:)
  end do

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

  deallocate(velw,dsource,source2d)

 call sep_end_prog()
end program
