!!
!! This program is to compute weights based on the diagonal of the Hessian
!! computed in ODCIG and mapped into angle by multiparametric OFF2ANG.
!!
!! The input datasets 
!! - dataset1 = ADCIG resulting from the transformation of an ODCIG previously 
!!              multiplied by diag(Hessian).
!! - dataset2 = ADCIG obtained with no previous multiplication 
!! 
!! They must be sorted according to:
!!
!!          z  -  cmp  -  angle
!!        dim1   dim2     dim3
!!
!!          
!!
!!

  program Hessang

  use sep
  use cgstep_mod
  use adj_mod
  use solver_prc_mod
  use hesang
  use laplac2 

!  use cartesian
!  use helix
!  use nhelix
!  use nhelicon
!  use nchelicon
!  use nchconest
!  use lapfac

  implicit none
  
  character*20					:: date,tempo
  integer					:: niter
  integer                         		:: i1,i2,i3
  integer					:: n1,n2,n3
  real                            		:: o1,o2,o3,d1,d2,d3,eps
  real,dimension(:,:),allocatable 		:: dado		!,ang1, ang2
  real,dimension(:),allocatable   		:: diagon, ang1, ang2, res
!  type (filter)            			:: aa

  call date_and_time(date,tempo)
  write(0,*) date,'  ',tempo

  call sep_init()
  call sep_begin_prog()
  !! get info from history file
  call from_history("n1",n1)
  call from_history("n2",n2)
  call from_history("n3",n3)
  call from_history("d1",d1)
  call from_history("d2",d2)
  call from_history("d3",d3)
  call from_history("o1",o1)
  call from_history("o2",o2)
  call from_history("o3",o3)

  !! get info from parameter
  call from_param("niter",niter,200)
  call from_param("eps",eps,1000.)

  !! write info out to history file
  call to_history("n1",n1)
  call to_history("n2",n2)
  call to_history("n3",n3)
  call to_history("d1",d1)
  call to_history("d2",d2)
  call to_history("d3",d3)
  call to_history("o1",o1)
  call to_history("o2",o2)
  call to_history("o3",o3)

!! allocate working space
  allocate(ang1(n1*n2),ang2(n1*n2),diagon(n1*n2),res(n1*n2),dado(n1,n2))
!!
write(0,*) "***Allocate***" 
!! compute the helix version of the laplacian regularizer
!eps_fac=.01
!aa=lapfac2(eps_fac,n1,15)
!

call laplac2_init(n1,n2)

!call sreed_window("in",3,(/n1,n2,n3/),(/n1,n2,1/),(/0,0,i3-1/),(/1,1,1/),4,ang1)
!call sreed_window("adcig2",3,(/n1,n2,n3/),(/n1,n2,1/),(/0,0,i3-1/),(/1,1,1/),4,ang2)
write(0,*) "***Read***"

do i3=1,n3
  ang1=0.;ang2=0.;diagon=0.
  call sep_read(ang1,"in")
  call sep_read(ang2,"adcig2")
  call random_number(res)
  res = res*0.01
  if (mod(i3,1).eq.0) write(0,*) "angle  ", i3, " of ", n3

  call hesang_init(n1,n2,ang1)
write(0,*) "***Init***"

  call solver_prc(m=diagon, d=ang2, Fop=hesang_lop,Sop=laplac2_lop,stepper=cgstep, &
	nSop=n1*n2,niter=niter,eps=eps,verb=.false.,resd=res)

write(0,*) "***Solver***"

  call hesang_close()

!  dado = reshape (diagon, (/n1,n2/))

!  call srite_window("out",3,(/n1,n2,n3/),(/n1,n2,1/),(/0,0,i3-1/),(/1,1,1/),4,dado)
  call sep_write(diagon,"out")

end do

call cgstep_close()

call date_and_time(date,tempo)
write(0,*) date,'  ',tempo

!! clean up
deallocate(dado,ang1,ang2,diagon)

call sep_end_prog()

end program Hessang
