program main
  use sep
  use base_solver_mod
  use wemva_struct_mod
  use wemva_kerpix_incore_mod
  use wemva_offdso_incore_mod
  use wemva_objdso_incore_mod
  use triangle2

  implicit none

  type(img_grid)                              :: igrid
  type(vel_grid)                              :: vgrid
  type(wemva_ftag)                            :: ftag_
  type(axis)                                  :: axk, ayk, azk, aiter
  real            , dimension(:), allocatable :: bslw, dimg, dslw, dspl, temp
  integer                                     :: rect1, rect2, rect3

  real                                        :: bvel_aver
  double precision                            :: ax, bx, cx
  double precision                            :: fa, fb, fc
  double precision                            :: tol, xmin, fmin
  integer                                     :: niter, iter, ITMAX, stat
  character(len=6)                            :: lnsrch

  character(len=256)                          :: mmov_

  mmov_ = "mmov"

  call sep_init(SOURCE)
  call getch_add_string("head=/dev/null")
  call getch_add_string("noheader=y")
  call set_no_putch()
  call sep_begin_prog()
  
  call from_param("lnsrch", lnsrch, "brent")
  call from_param("ITMAX" , ITMAX , 10     )

  call from_param("niter", niter, 5) 
  call from_param("rect1", rect1, 5)
  call from_param("rect2", rect2, 5)
  call from_param("rect3", rect3, 5)

  call from_param("oxk", axk%o, 0.)
  call from_param("dxk", axk%d, 1.)
  call from_param("nxk", axk%n, 1 )

  call from_param("oyk", ayk%o, 0.)
  call from_param("dyk", ayk%d, 1.)
  call from_param("nyk", ayk%n, 1 )

  call from_param("ozk", azk%o, 0.)
  call from_param("dzk", azk%d, 1.)
  call from_param("nzk", azk%n, 1 )


  stat = wemva_kerpix_incore_init(vgrid, igrid, ftag_)
  stat = wemva_offdso_incore_init(igrid%ax%n ,igrid%ay%n ,igrid%az%n ,igrid%ahx%n,igrid%ahy%n,igrid%ahz%n, &
                                  igrid%ahx%o,igrid%ahy%o,igrid%ahz%o,igrid%ahx%d,igrid%ahy%d,igrid%ahz%d)
  stat = base_spln3d_init(vgrid%ax,vgrid%ay,vgrid%az,axk,ayk,azk)
  call triangle2_init(rect1,rect3,vgrid%ax%n,vgrid%az%n)
 
  allocate(dimg(igrid%ax%n*igrid%ay%n*igrid%az%n*igrid%ahx%n*igrid%ahy%n*igrid%ahz%n))
  allocate(temp(igrid%ax%n*igrid%ay%n*igrid%az%n*igrid%ahx%n*igrid%ahy%n*igrid%ahz%n))
  allocate(bslw(vgrid%ax%n*vgrid%ay%n*vgrid%az%n))
  allocate(dslw(vgrid%ax%n*vgrid%ay%n*vgrid%az%n))
  allocate(dspl(axk%n*ayk%n*azk%n))

  dimg = 0.
  temp = 0.
  bslw = 0. 
  dslw = 0.
  dspl = 0.

  call sseek(ftag_%bslw_, 0, 0)
  call sep_read(bslw, ftag_%bslw_)

  call auxinout(ftag_%dimg_)
  call wemva_histo_imag(igrid, ftag_%dimg_)
  call wemva_fio_imag_init(igrid, ftag_%dimg_)

  aiter%o = 1
  aiter%d = 1
  aiter%n = niter

  call auxinout(mmov_)

  call houaxis(mmov_, vgrid%ax, 1)
  call houaxis(mmov_, vgrid%ay, 2)
  call houaxis(mmov_, vgrid%az, 3)
  call houaxis(mmov_, aiter   , 4)

  tol = 1.d-8

  do iter=1, niter
     write(0,*)"***********************"
     write(0,*)"iter = ", iter
     write(0,*)"***********************"
     write(0,*)"dI = D bI"

!    compute the perturbed image
!    ---------------------------
     stat = wemva_objdso_incore_oper(.false.,.false.,bslw,dimg)        !!.. dI = D bI

     write(0,*)"D'D bI"
!    compute D'D bI used for back projection
!    ---------------------------------------
     stat = wemva_offdso_incore_oper(.true. ,.false.,temp,dimg)        !!.. D'D bI
     call sseek(ftag_%dimg_, 0, 0)
     call sep_write(temp, ftag_%dimg_)

     write(0,*)"F'D'D bI"
!    compute the gradient
!    --------------------
     stat = wemva_kerpix_incore_oper(.true. ,.false.,dslw,temp)        !!.. F'D'D bI

!    smooth the gradient with B-spline basis function
!    ------------------------------------------------
     stat = base_spln3d_oper(.true. , .false., dspl, dslw)   
     stat = base_spln3d_oper(.false., .false., dspl, dslw)

!    searching direction is minus the gradient
!    -----------------------------------------
     dslw =  -dslw/maxval(abs(dslw))
     call sseek(ftag_%dslw_, 0, 0)
     call sep_write(dslw, ftag_%dslw_)

     temp = 0.

!    set up the initial brackets
!    ---------------------------
     write(0,*) "line search begins..."
     bvel_aver = sum(1./bslw)/size(bslw)
     !ax = 0.
     !bx = + 100/(bvel_aver+100)/bvel_aver

     ax = 0.
     bx = 0.00002
!    bracketing
!    ----------
     write(0,*)"bracketing..."
     call base_mnbrak_genfnc(base_objfnc, wemva_objdso_incore_oper, bslw, temp, dslw, ax, bx, cx, fa, fb, fc)
     write(0,*)"ax,bx,cx",ax,bx,cx
     write(0,*)"fa,fb,fc",fa,fb,fc

     write(0,*)"Golden/Brent section..."

!    line search with Golden or Brent section method
!    --------------------------------------
     select case(lnsrch)
       case ("golden")
            call base_golden_genfnc(base_objfnc, wemva_objdso_incore_oper, bslw, temp, dslw, ax, bx, cx, tol, fmin, xmin)
       case ("brents")
            call base_brents_genfnc(base_objfnc, wemva_objdso_incore_oper, bslw, temp, dslw, ax, bx, cx, tol, fmin, xmin, ITMAX)
     end select
     write(0,*)"ax,bx,cx",ax,bx,cx
     write(0,*)"xmin = ", xmin, "fmin = ", fmin

!    update the background slowness
!    ------------------------------
     bslw = bslw + xmin*dslw
     call sseek(ftag_%bslw_, 0, 0)
     call sep_write(bslw, ftag_%bslw_)

     call sep_write(bslw, mmov_)

  end do


  
end program
