module OW_wemva_feval

use sep 
use OW_type
use OW_grabvel_areal
use OW_exec_time
use OW_parms_areal
use OW_taper_areal
use OW_wei_areal
!use OW_imageIO
use OW_tim2ang
use solver_lbfgsb_mod
use base_spline_mod
use fx2d
use lloyd_vel

implicit none

real, private, allocatable                            :: pick1(:)

contains
subroutine feval_init(pick1_in)
  real            :: pick1_in(:)
  allocate(pick1(img%xm%n))
  pick1=pick1_in
end subroutine
!================================================
subroutine feval_close()
  deallocate(pick1)
end subroutine
!================================================
subroutine feval(vel0,fun,g,numg)

real                                          :: start, finish, start1, finish1, start2, finish2
real                                          :: a1,b1,vmin,vmax,stat1,v0,r,p_up,p_lo
integer                                       :: ixs,ino,ni,i,st,ishot,t3,iz,j,k,iter,ix,ish,titer,i1,i2,i3,minfun(1)
real, allocatable                             :: deltaI(:,:,:),scatS(:),deltaS(:)
real, allocatable                             :: R_image(:,:,:),B_image(:,:,:)
real, allocatable                             :: vref(:,:,:)
complex, allocatable                          :: recv(:,:),source(:,:),source1(:)
logical                                       :: cond,bounds,cbounds

 real, allocatable                            :: vel(:),taper(:)   !, l(:), u(:)
 double precision, optional                   :: g(:)
 integer, optional                            :: numg
 double precision                             :: fun,vel0(:)
 double precision, allocatable                :: gspl(:)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 start=0.;finish=0.
 call owexec_time(start)
 if (axkn.ne.99999 .or. azkn.ne.99999) then
    stat1 = base_spln2d_init(img%xm,img%z,axk,azk)
    allocate(gspl(axk%n*azk%n));gspl=dble(0.)
!    stat1 = base_spln3d_init(img%xm,img%ym,img%z,axk,ayk,azk)
 end if

 if (extrap(1:4).eq."pspi") then
    if(lloyd) then
      allocate(vref(nref,2,img%z%n))
      call lloyd_vel1(sngl(vel0),vref)
    end if
 end if

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Memory allocation
!
!allocate(l(nmax));l=dble(0.)
!allocate(u(nmax));u=dble(0.)
allocate(deltaS(img%z%n*kxm%n));deltaS=0.
allocate(scatS(img%z%n*img%xm%n));scatS=0.
allocate(deltaI(kxm%n,nsecdim,img%z%n));deltaI=0.
allocate(R_image(img%xm%n,nsecdim,img%z%n));R_image=0.
allocate(B_image(img%xm%n,nsecdim,img%z%n));B_image=0.
allocate(vel(img%z%n*kxm%n),taper(kxm%n))
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Set up the taper
!
 call owtaper(tapertr,taper)

! call from_param("bounds",bounds,.false.)
! if (bounds) then
!    call from_param("cbounds",cbounds,.false.)
!    if (cbounds) then
!       call from_param("vmin",vmin)
!       call from_param("vmax",vmax)
!       l=dble(1/vmax)
!       u=dble(1/vmin)
!    else
!       call from_param("p_up",p_up)
!       call from_param("p_lo",p_lo)
!       if (p_up >= 1.) p_up=1+p_up/100
!       if (p_lo >= 1.) p_lo=1-p_lo/100
!       l=p_lo*x
!       u=p_up*x
!    end if
!    nbd=2
! else
!    nbd=0
! end if 

! if (allocated(pick1)) then
!    call oweditvel(pick1,x,vel)
! end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Computing bkg.image
!
       start1=0.;finish1=0.;
       if (allocated(R_image)) deallocate(R_image) ;allocate(R_image(img%xm%n,nsecdim,img%z%n));R_image=0.
       if (allocated(B_image)) deallocate(B_image) ;allocate(B_image(img%xm%n,nsecdim,img%z%n));B_image=0.
       B_image=0.;R_image=0.;

       call sseek("rec",0,0)
       call sseek("sou",0,0)
       call owexec_time(start1)

       if (not(areal)) then
          if (allocated(source1)) deallocate(source1); allocate(source1(sou%w%n))
          source1=0;
          if (allocated(source)) deallocate(source); allocate(source(rec%h%n,sou%w%n))
          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
             if (allocated(source)) deallocate(source); allocate(source(sou%h%n,sou%w%n))
             source=0;
             call sreed("sou",source,8*sou%h%n*sou%w%n)
             vel=sngl(vel0)
             ni=1
          else
             call owgrab_vel(sngl(vel0),vel,.false.,ni,ishot)
          end if
          if (allocated(recv))   deallocate(recv);   allocate(recv(rec%h%n,rec%w%n))
          recv=0;deltaI=0
          call sreed("rec",recv,8*rec%h%n*rec%w%n)
          if (lloyd) then
             call owwei_init()
             call owwei(taper=taper,recv=recv,source=source,vel=vel,image=deltaI,sref=vref)
          else
             call owwei_init()
             call owwei(taper=taper,recv=recv,source=source,vel=vel,image=deltaI)
          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)
       write(3,112) "COMPUTED BACKGROUND IMAGE IN ",(finish1-start1)/60," min OF ITERATION # ", iter
       if (present(g)) then
          call to_history("n4",numg,"bkgimg")
          call srite("bkgimg",B_image,4*img%z%n*nsecdim*img%xm%n)
       end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Computing DSO
!
      call owdso(.false. , B_image)
      if (nxfilt.ne.99999) then
         stat1 = fx2d_lop (img%z%n,nsecdim,img%xm%n,15,150,5,int(1.5*img%z%n/(5+1)),&
                 int(1.5*img%xm%n/(50-5+1)),20.,B_image,R_image)
         B_image=R_image
      end if
      call owdso(.true.  , B_image)
      if (allocated(R_image)) deallocate(R_image)
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      !  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
                if (exist_file("horv0")) then
                   if ((i-1)*img%z%d+zmin < int(pick1(k))) B_image(k,j,i)=0
                end if
                fun = fun + B_image(k,j,i)*B_image(k,j,i)
             end do
          end do
       end do
       if (present(g)) then
          call to_history("n4",numg,"Dimage")
          call srite("Dimage",B_image,4*img%z%n*nsecdim*img%xm%n)
       end if

    if (present(g)) then
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       !
       ! Computing the gradient
       !

       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
          if (allocated(source1)) deallocate(source1); allocate(source1(sou%w%n))
          source1=0;
          if (allocated(source)) deallocate(source); allocate(source(rec%h%n,sou%w%n))
          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
             if (allocated(source)) deallocate(source); allocate(source(sou%h%n,sou%w%n))
             source=0;
             call sreed("sou",source,8*sou%h%n*sou%w%n)
             vel=sngl(vel0)
             ni=1
          else
             call owgrab_vel(sngl(vel0),vel,.false.,ni,ishot)
          end if
             
          if (allocated(recv))   deallocate(recv);   allocate(recv(rec%h%n,rec%w%n))
          recv=0;deltaS=0
          call sreed("rec",recv,8*rec%h%n*rec%w%n)

          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
          call owexec_time(start2)
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          ! B-spline representation of grad
          gspl=dble(0.)
          stat1 = base_spln2d_oper(.true. , .false., gspl, g)
          stat1 = base_spln2d_oper(.false. ,.false., gspl, g)
          call owexec_time(finish2)
          !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
       end if

       write(0,*) "main => ",maxval(g),minval(g)
       write(3,112) "SMOOTHED GRADIENT IN ",(finish2-start2)/60

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! importante verificar
!
       if (allocated(pick1)) then
          call oweditvel(pick1,g)
       end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       call owexec_time(finish1)
       write(3,112) "COMPUTED THE GRADIENT IN ",(finish1-start1)/60," min OF ITERATION # ", iter
       call to_history("n4",numg,"grad")
       call srite("grad",sngl(g),4*img%z%n*img%xm%n)
    end if

 111 format(a26,i4,a8,i4,a4,f6.2,a4)
 112 format(a29,f10.2,a20,i4)
 113 format(2(a8,i4),a8,a5,a8,e15.6,2(a8,I5))
! if (allocated(l)) deallocate(l)
! if (allocated(u)) deallocate(u)
! if (allocated(x)) deallocate(x)
! if (allocated(nbd)) deallocate(nbd)
 if (allocated(source)) deallocate(source)
 if (allocated(recv)) deallocate(recv)
 if (allocated(deltaS)) deallocate(deltaS)
 if (allocated(scatS)) deallocate(scatS)
 if (allocated(deltaI)) deallocate(deltaI)
 if (allocated(R_image)) deallocate(R_image)
 if (allocated(gspl)) deallocate(gspl)
 if (allocated(B_image)) deallocate(B_image)
 if (allocated(vel)) deallocate(vel)
 if (allocated(vref)) deallocate(vref)
! if (allocated(taper)) deallocate(taper)
! if (allocated(vel0)) deallocate(vel0)

 close(2)
 close(3)
11 format(a8,i4,a11,i5,a8,i4,2(a16,f20.10))

end subroutine
end module
