program SLS_
  !! Paul Sava (paul@sep.stanford.edu)
  !! 1/09/03

  !! real input/output
  !! input  is  z- mx       - other
  !! output is  z- mx - dip - other 

  use sep
  use util_
  use slo_
  use dottest_
  use solver_tiny_
  use cgstep_

  implicit none

  integer :: stat
  logical :: verb
  logical :: adj,add,dot,inv
  integer :: niter

  type(myaxis) :: a1,a2,a3,aa,al
  integer      :: i1,i2,i3,ia
  real, pointer :: image(:), slant(:)

  call sep_init()
  call from_param("verb",verb,.false.)
  call from_param("adj",adj,.false.)
  call from_param("add",add,.false.)
  call from_param("dot",dot,.false.)
  call from_param("inv",inv,.false.)

  if(dot) adj=.true.
  if(inv) adj=.true.

  call hinaxis(a1,1); a1%l="z"; call houaxis(a1,1)
  call hinaxis(a2,2); a2%l="x"; call houaxis(a2,2)

  if(adj) then
     call hinaxis(a3,3); a3%l="a3"; call houaxis(a3,4)
     call from_param("na",aa%n)
     call from_param("da",aa%d)
     call from_param("oa",aa%o)
     aa%l="a"
     call houaxis(aa,3)
  else
     call hinaxis(a3,4); a3%l="a3"; call houaxis(a3,3)
     aa%l="a"
     call hinaxis(aa,3)
  end if

  call from_param("nl",al%n,1)
  call from_param("dl",al%d)
  call from_param("ol",al%o)
  aa%l="l"

  if(inv) then
     call from_param("niter",niter,10)
  end if

  call sep_close()
  call slo_init(a1,a2,aa,al,verb)

  if(verb) then
     call printaxis(a1)
     call printaxis(a2)
     call printaxis(a3)
     call printaxis(al)
     call printaxis(aa)
  end if

  !----------------------------------------------------------------
  if(verb) write(0,*) "allocate"
  allocate(image(a1%n*a2%n     ),stat=stat)
  if(stat /=0) call seperr("cant allocate image"); image=0.
  allocate(slant(a1%n*a2%n*aa%n),stat=stat)
  if(stat /=0) call seperr("cant allocate slant"); slant=0.
  !----------------------------------------------------------------
  if(dot) then
     if(verb) write(0,*) "dottest"
     call dottest(slo_lop,slant,image)
     if(verb) write(0,*) "dottest OK"
  else
     if(inv) then
        if(verb) write(0,*) "inversion"
        call sep_read(image)
        call solver_tiny(slant,image,slo_lop,cgstep,niter)
        call sep_write(slant)
        if(verb) write(0,*) "inversion OK"
     else
        if(adj) then
           if(verb) write(0,*) "adjoint"
           call sep_read(image) 
           stat=slo_lop(adj,add,slant,image)
           call sep_write(slant)
           if(verb) write(0,*) "adjoint OK"
        else
           if(verb) write(0,*) "forward"
           call sep_read(slant)
           stat=slo_lop(adj,add,slant,image)
           call sep_write(image)
           if(verb) write(0,*) "forward OK"
        end if
     end if
  end if
  !----------------------------------------------------------------
  call slo_close()
  if(verb) write(0,*) "deallocate"
  deallocate(image, stat=stat)
  if(stat /=0) call seperr("cant deallocate image")
  deallocate(slant, stat=stat)
  if(stat /=0) call seperr("cant deallocate slant")
  !----------------------------------------------------------------
  call exit(0)
end program SLS_










