!!
!! 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 irls_mod
!  use solver_irls_mod
  use adj_mod
  use identity_mod
  use laplac2 
  use igrad1
  use cgstep_mod
  use solver_wprc_mod

  implicit none
  
  character*20					:: date,tempo
  integer					:: niter
  integer                         		:: i1,i2,i3,i,j
  integer					:: n1,n2,n3
  real                            		:: o1,o2,o3,d1,d2,d3,eps
!  real,dimension(:,:),allocatable 		:: ang_1, ang_2
  real,dimension(:),allocatable   		:: diagon, ang1, ang2,data

  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("n1",n1)
  call from_param("n2",n2)
  call from_param("n3",n3)
  call from_param("d1",d1)
  call from_param("d2",d2)
  call from_param("d3",d3)
  call from_param("o1",o1)
  call from_param("o2",o2)
  call from_param("o3",o3)

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

  !! 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(ang_1(n1,n2),ang_2(n1,n2),ang1(n1*n2),ang2(n1*n2),diagon(n1*n2),data(n1*n2))
  allocate(ang1(n1*n2),ang2(n1*n2),diagon(n1*n2),data(n1*n2))
!!

call laplac2_init(n1,n2)
call identity_init(n1,n2)

do i3=1,n3
  ang1=0.;ang2=0.;diagon=0.
  call sreed("adcig1",ang1,n1*n2*4)
  call sreed("adcig2",ang2,n1*n2*4)
!  call sep_read(ang_1,"adcig1")
!  call sep_read(ang_2,"adcig2")

  if (mod(i3,1).eq.0) write(0,*) "angle  ", i3, " of ", n3

  data=ang2*ang1/(ang1*ang1+abs(sum(ang1)/size(ang1)));
 
!  call solver_wprc(m=diagon, d=data, Fop=identity_lop,Sop=identity_lop,stepper=cgstep, &
!	nSop=n1*n2,Wdiag=l1,niter=niter,nfreq=2,eps=eps,verb=.true.)
  call solver_wprc(m=diagon, d=data, Fop=identity_lop,Sop=igrad1_lop,stepper=cgstep, &
	nSop=n1*n2,Wdiag=cauchy,niter=niter,nfreq=2,eps=eps,verb=.true.)
!  call solver_wprc(m=diagon, d=data, Fop=identity_lop,Sop=laplac2_lop,stepper=cgstep, &
!	nSop=n1*n2,Wdiag=l1,niter=niter,nfreq=2,eps=eps,verb=.true.)

diagon=data-diagon
  call sep_write(diagon,"out")
!  call sep_write(data,"out")

end do

call cgstep_close()

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

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

call sep_end_prog()

end program Hessang
