program Hessang2
  use lint2
  use bin2
  use igrad2
  use laplac2
  use solver_smp_mod
  use solver_prc_mod
  use solver_reg_mod
  use dottest
  use icai1
  use polydiv
  use helicon
  use cgstep_mod
  use createhelixmod
  use lapfac
  use mad_routines
  use sep


  implicit none
  integer                            :: n1,n2,m1,m2,na,stat1,stat2,method
  integer                            :: niter,id
  integer                            :: pef1,pef2   
  integer                            ::center1,center2,gap1,gap2,nfilt1,nfilt2
  integer                            :: center(2), gap(2), pef_axis(2), data_axis(2)
  real                               :: o1,o2,d1,d2,xmin,xmax,ymin,ymax
  real                               :: eps,eps_fac
  real,   dimension(:),pointer       :: data_1D,data_1D_grad,data_1D_ones,model_1D,fold_1D,model0,initsol
  real,   dimension(:),pointer       :: wt_in,wt_conv,wt,t_d,pef_in,res,res1,res1_grad,res2,res3,data_2
  real,   dimension(:),pointer       ::trend_1D
  real,   dimension(:),allocatable   ::wt_conv2,data_detrend_1D,mask_out
  real,    dimension( :,:),pointer   :: data ,xy,data_detrend
  logical                            :: norm,reg,bin
  type(filter)                       :: aa,peff

  call initpar()

  call from_history("n1",n1)
  call from_history("n2",n2)
  call from_param("m1",m1,200)
  call from_param("m2",m2,200)
  call from_param("method",method,1)
  call from_param("niter",niter,20)
  call from_param("eps",eps,.01)
  call from_param("norm",norm,.false.)
  call from_param("bin",bin,.false.)
  allocate( data(n1,n2),data_1D(n1),data_1D_ones(n1),model_1D(m1*m2),fold_1D(m1*m2),xy(n1,2))
  allocate(wt_in(n1),wt_conv(n1),wt(n1),wt_conv2(n1),model0(m1*m2),initsol(m1*m2))
  allocate(res(n1),res1(n1),res2(n1),res3(n1),res1_grad(n1),data_2(n1))
  allocate(data_detrend(n1,n2),data_detrend_1D(n1),trend_1D(m1*m2),mask_out(m1*m2))
  call sreed("in", data, 4*n1*n2) 

  allocate(t_d(2)) ! Find wt by convolving the trackends and bad data with 1,1
  t_d(1)=1.
  t_d(2)=1.
  call icai1_init(t_d,1)
  stat2 =  icai1_lop(.false.,.false.,wt_in,wt_conv) !lhcalc.  d*=Grad d;
  wt_conv2=wt_conv  
  wt=0.
  where (wt_conv2.eq.0.0) wt=1.


  t_d(1)=1
  t_d(2)=-1
  call icai1_init(t_d,1) !reinitialize the convolution with 1,-1

  call madchain_init(wt, n1)
  call lint2_init(m1,m2,o1,d1,o2,d2,xy )
  call bin2_init(m1,m2,o1,d1,o2,d2,xy )


  select case(method)
  case (1)     !  L*d , bin the data
     if (bin) then  
        stat1=bin2_lop(.true.,.false.,model_1D,data_1D)
     else
        stat1=lint2_lop(.true.,.false.,model_1D,data_1D)
     end if
     if (norm) then !normalize by fold 
        data_1D_ones=1.
        if (bin) then
           stat1=bin2_lop(.true.,.false.,fold_1D,data_1D_ones)
        else
           stat1=lint2_lop(.true.,.false.,fold_1D,data_1D_ones)
        end if
        where (fold_1D.ne.0.) model_1D=model_1D/abs(fold_1D) 
        end if
        if  (exist_file("data_detrend")) then
           call detrend(model_1D,m1,m2,trend_1D) !detrend and remove average
           if (bin) then  
              stat1=bin2_lop(.false.,.false.,trend_1D,data_detrend_1D)
           else
              stat1=lint2_lop(.false.,.false.,trend_1D,data_detrend_1D)
           end if

           data_detrend=data
           data_detrend(:,3)=data(:,3)-data_detrend_1D  !remove trend from data
           call to_history("n1",n1 ,"data_detrend")
           call to_history("n2",n2 ,"data_detrend")
           call sep_write(data_detrend,"data_detrend")
        end if
        if (exist_file("maskdat")) then
           mask_out=0.
           where (fold_1D.ne.0.) mask_out=1.
           call to_history ("n1",  m1,"maskdat")
           call to_history ("n2",  m2,"maskdat")
           call sep_write (mask_out,"maskdat")
           call to_history ("esize",4,"maskdat")
        end if
     case(5)    ! W d/dt [L (1/H) p - d] 
        eps_fac=.0001
        na=15
        aa=lapfac2(eps_fac,m1,50)
        call polydiv_init(m1*m2,aa)
        if (exist_file("m0")) then    ! if starting guess is provided
           call sreed("m0",model0,4*m1*m2)
           call helicon_init(aa)
           stat1=helicon_lop(.false.,.false.,model0,initsol)  ! Hh=p
        else
           initsol=0.0
        end if
        if (bin) then  
           call solver_prc(model_1D, data_1D,bin2_lop,polydiv_lop,cgstep,m1*m2&
           &,niter,eps,verb=.true.,Wop=madchain_exec,p0=initsol,resd=res)
        else
           call solver_prc(model_1D, data_1D,lint2_lop,polydiv_lop,cgstep,m1*m2&
           &,niter,eps,verb=.true.,Wop=madchain_exec,p0=initsol,resd=res)
        end if
        call cgstep_close()
     case(6)    ! W d/dt [L (1/A) p - d] where A is a 2D pef estimated on the model
        call from_aux("pef","n1",pef1)
        call from_aux("pef","n2",pef2)
        call from_param("center1",center1)
        call from_param("center2",center2)
        call from_param("a1",nfilt1)
        call from_param("a2",nfilt2)
        call from_param("gap1",gap1,0)
        call from_param("gap2",gap2,0)
        call from_param("reg",reg,.FALSE.)
        allocate(pef_in(pef1*pef2))
        call sreed("pef",pef_in,4*pef1*pef2)
        center(1)=center1
        center(2)=center2
        gap(1)=gap1
        gap(2)=gap2
        data_axis(1)=m1
        data_axis(2)=m2
        pef_axis(1)=nfilt1
        pef_axis(2)=nfilt2
        peff =  createhelix(data_axis,center, gap, pef_axis)
        peff%flt =pef_in
        !     allocate(peff%mis(n1))
        call polydiv_init(m1*m2,peff)
        if (exist_file("m0")) then    ! if starting guess is provided
           call sreed("m0",model0,4*m1*m2)
           call helicon_init(aa)
           stat1=helicon_lop(.false.,.false.,model0,initsol)  ! Hh=p
        else
           initsol=0.0
        end if
        if (bin) then
           call solver_prc(model_1D, data_1D,bin2_lop,polydiv_lop,cgstep,m1*m2&
           &,niter,eps,verb=.true.,Wop=madchain_exec,p0=initsol,resd=res)
        else
           call solver_prc(model_1D, data_1D,lint2_lop,polydiv_lop,cgstep,m1*m2&
           &,niter,eps,verb=.true.,Wop=madchain_exec,p0=initsol,resd=res)
        end if

        call cgstep_close()
        end if
     end select


     call sep_write (model_1D)
     call to_history ("n1",  m1)
     call to_history ("n2",  m2)
     call to_history ("o1",  o1)
     call to_history ("o2",  o2)
     call to_history ("d1",  d1)
     call to_history ("d2",  d2)
     call to_history ("esize",4)

     call madchain_close()
     deallocate( data,data_1D,model_1D,xy,wt_in,wt_conv,wt,wt_conv2,model0,initsol,res1,res2)
     deallocate(res,res3,res1_grad,data_2,mask_out,trend_1D,data_detrend,data_detrend_1D)
     deallocate(data_1D_ones,fold_1D)
   end program Madmain
