program DIM2_
  !! Analytical image perturbation
  !! Paul Sava (paul@sep.stanford.edu)
  !! 6/20/01

  !! complex input/output
  !! input  is  z- hx- hy- mx- my ( image)
  !! output is  z- hx- hy- mx- my (dimage)

  use sep
  use utilities
  use dimfft_
  use dimdrv2_
  use dimwgh2_

  implicit none

  integer :: stat
  logical :: verb

  type(paxis) :: az,ahx,ahy,amx,amy
  integer     :: iz,ihx,ihy,imx,imy

  !! z-hx-hy-mx-my
  complex, dimension(:,:,:,:,:), pointer :: cdata

  call sep_init()
  call from_param("verb",verb,.false.)

  call hinaxis(az ,1); az%l ="z";  call houaxis(az ,1)
  call hinaxis(ahx,2); ahx%l="hx"; call houaxis(ahx,2)
  call hinaxis(ahy,3); ahy%l="hy"; call houaxis(ahy,3)
  call hinaxis(amx,4); amx%l="mx"; call houaxis(amx,4)
  call hinaxis(amy,5); amy%l="my"; call houaxis(amy,5)

  call sep_close()

  if(verb) write(0,*) "fft init"
  call dimfft_init(az,ahx,ahy,amx,amy)
  if(verb) write(0,*) "drv init"  
  call dimdrv2_init(az,ahx,ahy,amx,amy)
  if(verb) write(0,*) "wgh init"
  call dimwgh2_init(az,ahx,ahy,amx,amy)
  !----------------------------------------------------------------
  if(verb) write(0,*) "allocate"
  allocate(cdata(az%n,ahx%n,ahy%n,amx%n,amy%n), stat=stat)
  if(stat /=0) call seperr("cant allocate cdata")
  cdata = 0.
  !----------------------------------------------------------------
  if(verb) write(0,*) "read"
  call sep_read(cdata)

  !! conjugate for compatibility with WEI
  cdata = conjg(cdata)
  write(0,*) minval(cabs(cdata)),maxval(cabs(cdata))

  if(verb) write(0,*) "drv"
  call dimdrv2(cdata) !! dR/dkz
  write(0,*) minval(cabs(cdata)),maxval(cabs(cdata))

  if(verb) write(0,*) "fft"
  call dimfft(cdata) !! FFT
  write(0,*) minval(cabs(cdata)),maxval(cabs(cdata))

  if(verb) write(0,*) "wgh"
  call dimwgh2(cdata) !! weight
  write(0,*) minval(cabs(cdata)),maxval(cabs(cdata))

  if(verb) write(0,*) "ift"
  call dimift(cdata) !! IFT
  write(0,*) minval(cabs(cdata)),maxval(cabs(cdata))

  !! conjugate for compatibility with WEI
  cdata = conjg(cdata)

  cdata(   1,:,:,:,:) = cdata(   1+1,:,:,:,:)
  cdata(az%n,:,:,:,:) = cdata(az%n-1,:,:,:,:)

  if(verb) write(0,*) "write"
  call sep_write(cdata)
  !----------------------------------------------------------------
  deallocate(cdata, stat=stat)
  if(stat /=0) call seperr("cant deallocate cdata")
  !----------------------------------------------------------------
  call exit(0)
end program DIM2_










