 program rtm

!(xm,xs,xh,z,t,v,wav,it0,npd)

! function img = rtm(xm,xs,xh,z,t,v,shot,wave)
!       xm --> x output domain
!       xs --> x shot positions
!       xh --> offset vector
!       v  --> must have (z,xm) dimension
!       ds and dxm must be equal
!

 use sep
 use sht_ext_mod 
! use imaging_mod
 use igrad1

 implicit none

 real, allocatable				:: source(:),data(:,:),vel(:,:),us(:,:,:)
 real, allocatable				:: dsource(:),ddata(:,:)   !,:),ur(:,:,:),img(:,:)
 real						:: xmax
 integer					:: is,it,ix,iz,ixw,iwav,stat,ih
 integer					:: nhpad,nto

 type(axis)					:: xs,xh,xm,xwrk,z,t,wav,vx
 

 call sep_init()

 call from_history("n1",t%n); call from_history("d1",t%d); call from_history("o1",t%o)
 call from_history("n2",xh%n); call from_history("d2",xh%d); call from_history("o2",xh%o)  
 call from_history("n3",xs%n); call from_history("d3",xs%d); call from_history("o3",xs%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",z%n); call from_aux("vel","d1",z%d); call from_aux("vel","o1",z%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("nhpad", nhpad, 100)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! To do: 01 - Different cmp, shot, receiver spacing
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 xm%n=int((xmax-xm%o)/vx%d)+1
 xwrk%n=2*nhpad+xh%n
 xwrk%d=vx%d
 xm%d=vx%d
 ixw=nhpad-int(xh%o/xh%d)+1

 call to_history("n3",t%n); call to_history("d3",t%d); call to_history("o3",t%o)
 call to_history("n2",xm%n); call to_history("d2",xm%d); call to_history("o2",xm%o)  
 call to_history("n1",z%n); call to_history("d1",z%d); call to_history("o1",z%o)


 write(0,*) "**************************************************"
 write(0,*) "***************** Parameters *********************"
 write(0,*) " "
 write(0,*) " ns = ",xs%n, " ds = ",xs%d," os = ",xs%o
 write(0,*) " nt = ",t%n, " dt = ",t%d," nh = ",xh%n," dh = ",xh%d," oh = ",xh%o
 write(0,*) " nwav = ",wav%n, " dwav = ",wav%d," owav = ",wav%o
 write(0,*) " nvx = ",vx%n," dvx = ",vx%d," ovx = ",vx%o," nvz = ",z%n," dvz = ",z%d," ovz = ",z%o
 write(0,*) " xmin = ",xm%o," xmax = ",xmax," nhpad = ", nhpad
 write(0,*) " nxm = ",xm%n," nwrk = ",xwrk%n," dwrk = ",xwrk%d," ixw = ",ixw
 write(0,*) "**************************************************"
 write(0,*) "**************************************************"

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

  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),data(t%n,xh%n),vel(z%n,vx%n))
  allocate(dsource(wav%n),ddata(t%n,xh%n))

  nto=t%n
  t%n=t%n+int(abs(wav%o)/wav%d)
  t%o=wav%o

  allocate(us(z%n,xwrk%n,3))!,ur(z%n,xwrk%n,3),img(z%n,xm%n))
!  allocate(us(z%n,xwrk%n,t%n))

open(unit=2,file='lixo')
  call sreed("source", source, 4*wav%n)
  dsource=0
  stat = igrad1_lop (.true.,.false.,dsource,source)
  call sreed("vel", vel, 4*z%n*vx%n)

  call extr_init(ixw,xs,xh,xwrk,xm,z,t,wav,nhpad)

  do is=1,xs%n
!      img=0;
      us=0;data=0;ddata=0
!      call sreed("in", data, 4*nt*nh)
!      do ih=1,xh%n
!         stat = igrad1_lop (.true.,.false.,ddata(:,ih),data(:,ih))
!      end do
!      us(1,ixw,1:wav%n)=dsource
!      ur(1,:,)=ddata()
      do it=2,t%n-1

         if (it <= wav%n)	us(1,ixw,1:2) = us(1,ixw,1:2) + dsource(it-1:it)
!write(0,*) "passei01"
!	 if(it <= wav%n) us(1,ixw,1:wav%n)=us(1,ixw,1:wav%n)+dsource
	 call sht_ext(is,it,vel,us)
write(0,*) "passei02"
!         call rec_ext(is,it,vel,ddata(t%n+1-it,:),ur)
!         do ix=1:nh+2*npd
!             do iz=1:nz
!                 img(iz,ix)=img(iz,ix)+us(iz,ix,nt-it+1)*ur(iz,ix,nt-it+1)
!             end do
!         end do
!         if (it>=int(abs(wav%o)/wav%d)) call srite("out",us(:,:,it),4*z%n*xwrk%n)
         if (it>=int(abs(wav%o)/wav%d)) call srite("out",us(:,:,3),4*z%n*xwrk%n)
!write(0,*) "passei03"
      end do
!     ixm=floor(((is-1)*ds+os+oh-npd*dh-xmin)/dxm)+1
!     image(:,ixm:ixm+nh-1+2*npd,is)=image(:,ixm-1:ixm+nh-2+2*npd,is)+img(:,:)

  end do
close (2)
end program
