program main
  use sep
  use off2ang_fd_mod
  implicit none


  real, dimension(:), pointer :: modl(:), data(:)
  integer                     :: nz, nhx, nax, stat
  real                        :: oz, ohx  
  real                        :: dz, dhx, dax
  real                        :: axmin, axmax 
  real, parameter             :: PI = 3.14159265359
  logical, parameter          :: T = .true. , F = .false.
  call initpar()

  call from_history("n1", nz ); call from_history("d1", dz ); call from_history("o1", oz )
  call from_history("n2", nhx); call from_history("d2", dhx); call from_history("o2", ohx)  

  call from_param("amin", axmin, -60.)
  call from_param("amax", axmax, 60.)

  nax = nhx
  if (nax == 1) then
      dax = 0
  else
      dax = (axmax - axmin) / (nax - 1)
  end if

  allocate(modl(nz*nhx), data(nz*nhx))

  call sreed("in", data, 4*nz*nhx)
  
  call off2ang_fd_init(dz, dhx, nz, nhx, axmin, axmax)

  stat = off2ang_fd_oper(T, F, modl, data)
  !call solver_tiny(m=modl, d=data, Fop=off2ang_sd_oper, stepper=cgstep, niter=50)


  call to_history("n1", nz); call to_history("d1", dz); call to_history("o1", oz)
  call to_history("n2", nax); call to_history("d2", dax); call to_history("o2", axmin) 

  call srite("out", modl, 4*nz*nhx)


end program