module steepest_mod

use sep
use OW_grabvel_areal
use OW_exec_time
use OW_parms_areal
use OW_taper_areal
use OW_wei_areal
use lloyd_vel
use base_spline_mod
use triangle2_mod
use cgstep_mod0
use solver_tiny_mod
implicit none
real, private, allocatable             :: vel0(:),taper(:),pick1(:)
integer, private                       :: iter, rect1, rect2,rect3,realz

contains 

!---------------------------------------------------------------------
subroutine fgeval_init(vel0_in,taper_in,iter_in,pick1_in)
 real                        :: vel0_in(:),taper_in(:)
 real,optional               :: pick1_in(:)
 integer                     :: iter_in
 iter=iter_in
 allocate(vel0(img%xm%n*nz),taper(kxm%n))
 vel0=vel0_in
 taper=taper_in
 if (present(pick1_in)) then
    allocate(pick1(img%xm%n))
    pick1=pick1_in
 end if
 if (nreal.ne.1) then
    realz=mod(iter,nreal)
 else
    realz=1
 end if

end subroutine
!---------------------------------------------------------------------
subroutine fgeval_close()
 if (allocated(vel0))   deallocate(vel0)
 if (allocated(taper)) deallocate(taper)
 if (allocated(pick1)) deallocate(pick1)
end subroutine
!---------------------------------------------------------------------
subroutine feval(fun,R_image)

real, optional                                :: R_image(:,:,:)
real                                          :: start, finish, start1, finish1
integer                                       :: ni,i,ishot,j,k
real                                          :: deltaI(kxm%n,nsecdim,img%z%n),illum(kxm%n,img%z%n),vel(img%z%n*kxm%n)
real                                          :: B_image(img%xm%n,nsecdim,img%z%n)
real, allocatable                             :: vref(:,:,:)
complex                                       :: recv(rec%h%n,rec%w%n),source(rec%h%n,sou%w%n),source1(sou%w%n)
double precision                              :: fun

deltaI=0.
B_image=0.
vel=0
illum=0
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Computing bkg.image
!
       if (extrap(1:4).eq."pspi") then
          if(lloyd) then
            allocate(vref(nref,2,img%z%n))
            call lloyd_vel1(vel0,vref)
          end if
       end if
       start1=0.;finish1=0.;
       B_image=0.

       call sseek("rec",0,0)
       call sseek("sou",0,0)
       call owexec_time(start1)
       if (not(areal)) then
          source1=0;
          source=0;
          call sreed("sou",source1,8*sou%w%n)
          source(abs(floor(rec%h%o/rec%h%d)) + 1,:)=source1
       end if
       do ishot=1,rec%xd%n
          start=0.;finish=0. 
          call owexec_time(start)
          if (areal) then
             source=0;
             call sreed_window("sou",5,(/sou%h%n,1,sou%w%n,sou%xd%n,  nreal/),&
                                       (/sou%h%n,1,sou%w%n,       1,      1/),&
                                       (/      0,0,      0, ishot-1,realz-1/),&
                                       (/      1,1,      1,      1,       1/),&
                                       8,source) 
!             call sreed("sou",source,8*sou%h%n*sou%w%n)
             vel=vel0
             ni=1
          else
             call owgrab_vel(vel0,vel,.false.,ni,ishot)
          end if
          recv=0;deltaI=0;illum=0.
          if (areal) then
             call sreed_window("rec",5,(/sou%h%n,1,sou%w%n,sou%xd%n,  nreal/),&
                                       (/sou%h%n,1,sou%w%n,       1,      1/),&
                                       (/      0,0,      0, ishot-1,realz-1/),&
                                       (/      1,1,      1,      1,       1/),&
                                       8,recv) 
          else
             call sreed("rec",recv,8*rec%h%n*rec%w%n)
          end if
          if (lloyd) then
             call owwei_init()
             call owwei(taper,recv,source,vel,deltaI,illum,vref)
          else
             call owwei_init()
             call owwei(taper,recv,source,vel,deltaI,illum)
          end if
          B_image(ni:ni+kxm%n-1,:,:)=B_image(ni:ni+kxm%n-1,:,:)+deltaI

          call owexec_time(finish)
          write(3,111) "EXTRAPOLATED WAVEFIELDS = ",ishot," OUT OF ",rec%xd%n," IN ",(finish-start)/60," min"
       end do

       call owexec_time(finish1)
       if(present(R_image)) then
          call to_history("n4",iter,"bkgimg")
          call srite("bkgimg",B_image,4*img%z%n*nsecdim*img%xm%n)
       end if
       write(0,112) "COMPUTED BACKGROUND IMAGE IN ",(finish1-start1)/60," min OF ITERATION # ", iter
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Computing DSO
!
      call owdso(.false. , B_image)
      call owdso(.true.  , B_image)

      if (allocated(pick1)) then
         do i=1,img%z%n
         do j=1,img%sh%n
         do k=1,img%xm%n
            if ((i-1)*img%z%d+zmin < int(pick1(k))) B_image(k,j,i)=0
         end do
         end do
         end do
      end if
      if (present(R_image)) then
         R_image=B_image
         call to_history("n4",iter,"Dimage")
         call srite("Dimage",R_image,4*img%z%n*nsecdim*img%xm%n)
      end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!  Computing the objective function
!
       fun=dble(0.)
       do i=1,img%z%n
          do j=1,img%sh%n
             do k=1,img%xm%n
                fun = fun + B_image(k,j,i)*B_image(k,j,i)
             end do
          end do
       end do
       if (allocated(vref)) deallocate(vref)       

 111 format(a26,i4,a8,i4,a4,f6.2,a4)
 112 format(a29,f10.2,a20,i4)
end subroutine
!---------------------------------------------------------------------
 subroutine geval(g,B_image)

 real                                          :: B_image(:,:,:)
 real                                          :: start, finish, start1, finish1, start2, finish2
 integer                                       :: ni,i,ishot,j,k,stat1
 real                                          :: deltaS(img%z%n*kxm%n),scatS(img%z%n*img%xm%n),vel(img%z%n*kxm%n)
 real, allocatable                             :: vref(:,:,:)
 complex                                       :: recv(rec%h%n,rec%w%n),source(rec%h%n,sou%w%n),source1(sou%w%n)
 double precision                              :: g(:),gspl(img%z%n*img%xm%n)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Computing the gradient
!
       if (extrap(1:4).eq."pspi") then
          if(lloyd) then
            allocate(vref(nref,2,img%z%n))
            call lloyd_vel1(vel0,vref)
          end if
       end if
       start1=0.;finish1=0.;g=dble(0.)

       call sseek("rec",0,0)
       call sseek("sou",0,0)
       call owexec_time(start1)
       if (not(areal)) then
          source1=0;
          source=0;
          call sreed("sou",source1,8*sou%w%n)
          source(abs(floor(rec%h%o/rec%h%d)) + 1,:)=source1
       end if

       start1=0.;finish1=0.;scatS=0.
       call owexec_time(start1)
       do ishot=1,rec%xd%n
          start=0.;finish=0. 
          call owexec_time(start)
          if (areal) then
             source=0;
             call sreed_window("sou",5,(/sou%h%n,1,sou%w%n,sou%xd%n,  nreal/),&
                                       (/sou%h%n,1,sou%w%n,       1,      1/),&
                                       (/      0,0,      0, ishot-1,realz-1/),&
                                       (/      1,1,      1,      1,       1/),&
                                       8,source) 
!             call sreed("sou",source,8*sou%h%n*sou%w%n)
             vel=vel0
             ni=1
          else
             call owgrab_vel(vel0,vel,.false.,ni,ishot)
          end if
          recv=0;deltaS=0

          if (areal) then
             call sreed_window("rec",5,(/sou%h%n,1,sou%w%n,sou%xd%n,  nreal/),&
                                       (/sou%h%n,1,sou%w%n,       1,      1/),&
                                       (/      0,0,      0, ishot-1,realz-1/),&
                                       (/      1,1,      1,      1,       1/),&
                                       8,recv) 
          else
             call sreed("rec",recv,8*rec%h%n*rec%w%n)
          end if
          call owwemva_init(ishot)
          if (lloyd) then
             call owwemva(.true.,taper,recv,source,vel,deltaS,B_image(ni:ni+kxm%n-1,:,:),vref)
          else
             call owwemva(.true.,taper,recv,source,vel,deltaS,B_image(ni:ni+kxm%n-1,:,:))
          end if
          if (areal) then
             scatS=scatS+deltaS
          else
             do j=ni,ni+kxm%n-1;do k=1,img%z%n
                scatS((k-1)*img%xm%n+j)=scatS((k-1)*img%xm%n+j)+deltaS((k-1)*kxm%n+j-ni+1)
             end do;end do
          end if

          call owexec_time(finish)
          write(3,111) "EXTRAPOLATED WAVEFIELDS = ",ishot," OUT OF ",rec%xd%n," IN ",(finish-start)/60," min"
       end do
       g=dble(-scatS)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Smoothing the gradient
       write(0,*) "main => ",maxval(g),minval(g)
       if (axkn.ne.99999 .or. azkn.ne.99999) then
          axk%n=axk%n*(2**(iter-1))
          axk%d=axk%d/(2**(iter-1))
          if (axk%d < dxkmin) then
             axk%d = dxkmin
             axk%n = nxkmin
          end if
          stat1 = base_spln2d_init(img%xm,img%z,axk,azk)
!          stat1 = base_spln3d_init(img%xm,img%ym,img%z,axk,ayk,azk)
       end if

!       if (exist_file("horv0")) call oweditvel(pick1,g)

       if (axkn.ne.99999 .or. azkn.ne.99999) then
          call owexec_time(start2)
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          ! B-spline representation of grad
          gspl=dble(0.)
          stat1 = base_spln2d_oper(.true. , .false., gspl, g)
          g=dble(0.)
          stat1 = base_spln2d_oper(.false. ,.false., gspl, g)
          call owexec_time(finish2)
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
       end if
!!
!       call triangle2_init(10*rect1,10*rect3,img%xm%n,img%z%n)
!       if (exist_file("horv0"))  call oweditvel(pick1,g)
!       call owexec_time(start2);scatS=0
!       stat1 = triangle2_lop(.true.,.false.,scatS,sngl(g))
!       g=dble(scatS)
!       call owexec_time(finish2)
!       if (exist_file("horv0"))  call oweditvel(pick1,g)
!       call triangle2_close()
       write(0,*) "main => ",maxval(g),minval(g)

       write(0,112) "SMOOTHED GRADIENT IN ",(finish2-start2)/60
       if (exist_file("horv0")) call oweditvel(pick1,g)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       call owexec_time(finish1)
       write(0,112) "COMPUTED THE GRADIENT IN ",(finish1-start1)/60," min OF ITERATION # ", iter
       call to_history("n4",iter,"grad")
       call srite("grad",sngl(g),4*img%z%n*img%xm%n)

       if (allocated(vref)) deallocate(vref)
 111 format(a26,i4,a8,i4,a4,f6.2,a4)
 112 format(a29,f10.2,a20,i4)
 end subroutine
!---------------------------------------------------------------------
 subroutine A2g(deltaS,R_image)

 real                                          :: R_image(img%xm%n,nsecdim,img%z%n)
 real                                          :: start, finish, start1, finish1, start2, finish2
 integer                                       :: ni,i,ishot,j,k,stat1
 real                                          :: deltaS(:),deltaI(kxm%n,nsecdim,img%z%n),vel(img%z%n*kxm%n)
 real, allocatable                             :: vref(:,:,:)
 complex                                       :: recv(rec%h%n,rec%w%n),source(rec%h%n,sou%w%n),source1(sou%w%n)

       R_image=0
       if (extrap(1:4).eq."pspi") then
          if(lloyd) then
            allocate(vref(nref,2,img%z%n))
            call lloyd_vel1(vel0,vref)
          end if
       end if
       start1=0.;finish1=0.

       call sseek("rec",0,0)
       call sseek("sou",0,0)
       call owexec_time(start1)
       if (not(areal)) then
          source1=0;
          source=0;
          call sreed("sou",source1,8*sou%w%n)
          source(abs(floor(rec%h%o/rec%h%d)) + 1,:)=source1
       end if

       start1=0.;finish1=0.;deltaI=0.
       call owexec_time(start1)
       do ishot=1,rec%xd%n
          start=0.;finish=0. 
          call owexec_time(start)
          if (areal) then
             source=0;
             call sreed("sou",source,8*sou%h%n*sou%w%n)
             vel=vel0
             ni=1
          else
             call owgrab_vel(vel0,vel,.false.,ni,ishot)
          end if
             
          recv=0;deltaI=0
          call sreed("rec",recv,8*rec%h%n*rec%w%n)
          call owwemva_init(ishot)
          if (lloyd) then
             call owwemva(.false.,taper,recv,source,vel,deltaS,deltaI,vref)
          else
             call owwemva(.false.,taper,recv,source,vel,deltaS,deltaI)
          end if
          R_image(ni:ni+kxm%n-1,:,:)=R_image(ni:ni+kxm%n-1,:,:)+deltaI
          call owexec_time(finish)
       end do
       if (allocated(vref)) deallocate(vref)

end subroutine
!---------------------------------------------------------------------
 subroutine ow_smoovel(x,vel0,pick1)

 real, optional             :: pick1(:),vel0(:)
 double precision           :: x(:),ex,exs
! double precision           :: gspl(axk%n*azk%n)
 real                       :: start2,finish2,slow(img%z%n*img%xm%n)
 real, allocatable          :: gspl(:)
 integer                    :: stat1,i,j,ix,iz
 
       call triangle2_init(rect1,rect3,img%xm%n,img%z%n)
       if (present(pick1)) call oweditvel(pick1,x,vel0)
       call owexec_time(start2)
       stat1 = triangle2_lop(.true.,.false.,slow,sngl(x))
       x=dble(slow)
       call owexec_time(finish2)
       if (present(pick1)) call oweditvel(pick1,x,vel0)
       call triangle2_close()

end subroutine
!---------------------------------------------------------------------
 subroutine owtri_init(rect1_in,rect2_in,rect3_in)
 integer           :: rect1_in,rect2_in,rect3_in
   rect1=rect1_in
   rect2=rect2_in
   rect3=rect3_in
 end subroutine
!---------------------------------------------------------------------
real function norm(x,pick1)
  double precision                  :: x(:)
  real,optional                     :: pick1(:)
  integer                           :: i,k
  norm=0.
  if (present(pick1)) then
     do i=1,img%xm%n
        do k=1,img%z%n
           if ((k-1)*img%z%d+img%z%o > pick1(i)) then
              norm=norm+x((k-1)*img%xm%n+i)**2
           end if
        end do
     end do
     norm=sqrt(norm)
  else
     norm = sqrt(sum(x*x))
  end if
end function
!---------------------------------------------------------------------
end module
