 program rtm_SO_pw

 use sep
 use ext_SO_pw_mod
! use fff_mod
 use sep3d_struct_mod
 use igrad1

 implicit none

 real, allocatable				:: pwsource(:,:)
 real, allocatable				:: source(:),vel(:,:),velw(:,:),us(:,:,:)
 real, allocatable				:: dsource(:),dsource2d(:,:),x(:)
 real						:: xmax,zmax,oto,theta
 integer					:: is,it,ix,iz,ixw,iwav,stat,ih,ixm,i,idelt
 integer					:: nto,nwav2

 type(axis)					:: xm,xwrk,z,t,wav,vx,vz,ang
 type(sep3d)					:: axgr
 

 call sep_init()
 call sep_begin_prog()

 call from_history("n1",t%n); call from_history("d1",t%d); call from_history("o1",t%o)

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

 call from_aux("vel","n1",vz%n); call from_aux("vel","d1",vz%d); call from_aux("vel","o1",vz%o)
 call from_aux("vel","n2",vx%n); call from_aux("vel","d2",vx%d); call from_aux("vel","o2",vx%o)

 call from_param("xmin", xm%o, vx%o)
 call from_param("xmax", xmax, vx%o+(vx%n-1)*vx%d)
 call from_param("zmax", zmax)

 call from_param("angmin", ang%o)
 call from_param("d_ang", ang%d)

 ang%n=int(2*abs(ang%o)/ang%d)+1 
 z%n=int((zmax-vz%o)/vz%d)+1
 z%o=vz%o
 z%d=vz%d
 xm%n=int((xmax-xm%o)/vx%d)+1
 xm%d=vx%d
 xwrk=xm

 call to_history("n1",z%n); call to_history("d1",z%d); call to_history("o1",z%o)
 call to_history("n2",xm%n); call to_history("d2",xm%d); call to_history("o2",xm%o)  
 call to_history("n3",xso%n); call to_history("d3",xso%d); call to_history("o3",xso%o)  
 call to_history("n4",1); call to_history("d4",1); call to_history("o4",1)

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

  if (xm%n>vx%n) call seperr("Image domain is bigger than the velocity model in x - pad velocity (use extent = 1)")

  allocate(source(wav%n),vel(vz%n,vx%n),velw(z%n,xm%n))
  allocate(dsource(wav%n),us(z%n,xwrk%n,3))

  nto=t%n;oto=t%o
  nwav2=int(abs(wav%o)/wav%d)+1
  t%n=t%n+nwav2
  t%o=wav%o

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

  w = shift_init(x, z, t)

  allocate(dsource2d(wav%n,xm%n))
  do i=1,xm%n
     dsource2d(:,i)=dsource
  end do
  deallocate(dsource)

!  allocate (C(w%n,xm%n)); C=0! nt/2 frequencies, no zero freq
!  C = rft_fwd_axis1(dsource2)

  call sreed("vel", vel, 4*vz%n*vx%n)
  ixi=int(abs(vx%o-xm%o)/vx%d)+1
  ixf=ixi+xm%n-1
  velw=vel(1:z%n,ixi:ixf)
  deallocate(vel)

  allocate(x(xm%n)) 
  do i=1,xm%n;x(i)=(i-1)*xm%d+xm%o;end do
  idelt1=maxval(int(abs(sind(ang%o))*x/velw(1,:)/t%d))+1
  deallocate(x)
  allocate(pwsource(idelt1,xm%n))

  call extr_init(xm,z,t,wav,velw)
!
! Dip loop - Propagates plane source wavefield
!
  do ip=1,ang%n

     pwsource=0; theta=(ip-1)*ang%d+ang%o
     do i=1,xm%n
        idelt=int(sind(theta)*(i-1)*xm%d/velw(1,i)/t%d)+1
        pwsource(idelt1+idelt+1:idelt1+idelt+nwav,i)=dsource
     end do

     us=0.;

     call init_sep3d(axgr, "OUTPUT","FLOAT", "REGULAR", ctag="green")
     call axis_allocate(axgr,3)
     axgr%o(1)=z%o; axgr%n(1)=z%n; axgr%d(1)=z%d
     axgr%o(2)=xwrk%o; axgr%n(2)=xwrk%n; axgr%d(2)=xwrk%d
     axgr%o(3)=oto; axgr%n(3)=nto; axgr%d(3)=t%d
     call sep3d_set_sep3d(axgr)
     call sep3d_write_description("green",axgr)

     do it=2,t%n-1
	 if (mod(it,250).eq.0) write(0,*) "p = ", (ip-1)*p%d+p%o,"(",ip,")"," from #",p%n," time = ",it," from #",t%n
	 if(it <= wav%n) us(1,:,1:2)=us(1,:,1:2)+dsource2d(:,it-1:it)
	 us(:,:,3)=0
	 call sht_ext(us)
	 if (it>=nwav2) call srite("green",us(:,:,3),4*z%n*xwrk%n)
	 us(:,:,1:2)=us(:,:,2:3);
     end do

     call sep3df_delete(axgr,"green")
     call init_sep3d("green",axgr,"SCRATCH")

     deallocate(pwsource)

  end do

  deallocate(velw,us)
  deallocate(dsource2d)

 call sep_end_prog()
end program
