 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 ext_mod 
 use sep3d_struct_mod
 use igrad1

 implicit none

 logical					:: forward
 real, allocatable				:: source(:),data(:,:),vel(:,:),us(:,:,:),ws(:,:,:)
 real, allocatable				:: dsource(:),ddata(:,:),ur(:,:,:),img(:,:),image(:,:,:)
 real						:: xmax,zmax
 real						:: oto,dto
 integer					:: is,it,ix,iz,ixw,iwav,stat,ih,ixm
 integer					:: nhpad,nto,nwav2

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

 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",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("nhpad", nhpad, 100)
 call from_param("zmax", zmax)
 call from_param("forward",forward,.false.)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! To do: 01 - Different cmp, shot, receiver spacing
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 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
 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("n4",xs%n,"green1"); call to_history("d4",xs%d,"green1"); call to_history("o4",xs%o,"green1")
! call to_history("n3",t%n,"green1"); call to_history("d3",t%d,"green1"); call to_history("o3",t%o,"green1")
! call to_history("n2",xwrk%n,"green1"); call to_history("d2",xwrk%d,"green1"); 
! call to_history("n1",z%n,"green1"); call to_history("d1",z%d,"green1"); call to_history("o1",z%o,"green1")
 
 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",xs%n); call to_history("d3",xs%d); call to_history("o3",xs%o)  

 if(exist_file("wrec")) then
    call to_history("n4",xs%n,"wrec"); call to_history("d4",xs%d,"wrec"); call to_history("o4",xs%o,"wrec")
    call to_history("n3",t%n,"wrec"); call to_history("d3",t%d,"wrec"); call to_history("o3",-t%o,"wrec")
    call to_history("n2",xwrk%n,"wrec"); call to_history("d2",xwrk%d,"wrec"); call to_history("o2",xwrk%o,"wrec")  
    call to_history("n1",z%n,"wrec"); call to_history("d1",z%d,"wrec"); call to_history("o1",z%o,"wrec")
 end if

 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),vel(z%n,vx%n))
  allocate(dsource(wav%n),ddata(t%n,xh%n))
  allocate(image(z%n,xm%n,xs%n))
  allocate(data(t%n,xh%n))
  allocate(ws(z%n,xwrk%n,1))

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

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

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

  call sreed("vel", vel, 4*z%n*vx%n)

  call extr_init(ixw,xs,xh,xwrk,xm,z,t,wav,nhpad)
!
! Shot loop
!
  do is=1,xs%n
      us=0;ddata=0;
      data=0
      call sreed("in", data, 4*t%n*xh%n)
!!
!! 1st derivative of the shot - accounts for phase
!!
      do ih=1,xh%n
         stat = igrad1_lop (.true.,.false.,ddata(:,ih),data(:,ih))
      end do
      xwrk%o=(is-1)*xs%d+xs%o+xh%o-nhpad*xwrk%d
      ixm=floor((xwrk%o-xm%o)/xm%d)+1

      write(0,*) ixw,xwrk%o,xm%o,ixm

!      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)=dto
!!      axgr%o(4)=xs%o; axgr%n(4)=xs%n; axgr%d(4)=xs%d
!      call sep3d_set_sep3d(axgr)
!      call sep3d_write_description("green",axgr)

!!
!! Propagates source
!!

!     do it=2,t%n-1
!	 if(it <= wav%n) us(1,ixw,1:2)=us(1,ixw,1:2)+dsource(it-1:it)
!	 call sht_ext(is,vel,us)
!	 if (it>=nwav2) call srite("green",us(:,:,3),4*z%n*xwrk%n)
!	 us(:,:,1:2)=us(:,:,2:3);us(:,:,3)=0
!	 if (mod(it,250).eq.0) write(0,*) "shot = ", is," from ",xs%n," time = ",it," from ",t%n
!     end do

!     if(forward) stop

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

!
! Back-propagates receivers
!
!     call to_history("o2",xwrk%o,"green1")  
     ur=0;img=0
     do it=t%n-nwav2,t%n-nto,-1
	if (mod(it,250).eq.0) write(0,*) "receivers from shot = ", is," from ",xs%n," time = ",it," from ",t%n
	ur(1,nhpad+1:nhpad+xh%n,2:3) = &
		ur(1,nhpad+1:nhpad+xh%n,2:3) + transpose(ddata(it-1:it,:))
        call rec_ext(is,it,vel,ur)
        if(exist_file("wrec")) then
		if (it>=nwav2) then
		   call srite("wrec",ur(:,:,1),4*z%n*xwrk%n)
		end if
	end if

!
! Applies the imaging condition
!
!        ws=0
!	if (mod())
!        call sreed_window("green",3,(/z%n,xwrk%n,t%n/),(/z%n,xwrk%n,1/),(/0,0,it-1/),(/1,1,1/),4,ws)
!        call imaging(ws(:,:,1),ur(:,:,1),img)
	ur(:,:,2:3)=ur(:,:,1:2);ur(:,:,1)=0

!        call srite("green1",ws(:,:,1),4*z%n*xwrk%n)
!
! Sum migrated images
!
!!!!!!!! Verificar os indices na imagem

!        image(:,ixm:ixm+xwrk%n-1,is)=image(:,ixm:ixm+xwrk%n-1,is)+img

     end do

!     call sep3df_delete(axgr,"green")

  end do

!  call srite("out",image,4*z%n*xm%n*xs%n)

!  deallocate(ddata,img,vel,dsource,ur,ws)
!  deallocate(data,image,us)

end program
