program UnLMO
  use SEP
  implicit none

  integer                              :: n1,n2,n3,m1,m2,pad2,tpfact,xpfact,ix,iw,jj
  real                                 :: o1,o2,p1,d1,d2,xloc,angle,w0,w,dw
  real                                 :: avel,pi,sig1,dx
  real, dimension(:,:,:), allocatable  :: data,modl
  real, dimension(:,:)  , allocatable  :: tmp
  complex, dimension(:,:), allocatable :: cdata,cmodl

  call sep_init()
  call from_history("n1",n1               ) !
  call from_history("n2",n2               ) !
  call from_history("n3",n3               )
  call from_history("o1",o1               ) !
  call from_history("o2",o2               ) !
  call from_history("p1",p1               ) !
  call from_history("d1",d1               )
  call from_history("d2",d2               )
  call from_param  ("avel",avel,6.2       )          
  tpfact=2  !! Padding factor to each side
  xpfact=3
  m1=pad2(n1)
  m2=pad2(n2)
  allocate( data(n1,n2,n3),  modl(m1,n2,n3) )
  allocate( cdata(m1,n2),cmodl(m1,n2), tmp(m1,n2) )
  data=0.;      modl=0.;
  cdata=(0.,0.);cmodl=(0.,0.);

  !! . . Definitions
  pi = 3.14159265; sig1 = +1.;
  w0=-pi/d1;  dw=2.*pi/(m1*d1);

  !! . . Read in File
  call sep_read( data )
  write(0,*) 'Using p1 slowness of: ',p1
  do jj=1,n3
     cdata(1:n1,:)=cmplx(data(1:n1,:,jj),0.)
     call ft1axis(0,sig1,m1,n2,cdata)
     do ix=1,n2
        xloc=(ix-1)*d2+o2
!        write(0,*) p1*(n2-ix)*d2
        do iw=1,m1
           w=w0+(iw-1)*dw 
           if (p1 .lt. 0 ) then 
              cmodl(iw,ix)=cdata(iw,ix)*exp(cmplx(0.,-w*p1*( (n2-ix)*d2 ) ) )
           else
              cmodl(iw,ix)=cdata(iw,ix)*exp(cmplx(0.,w*p1*xloc))
           end if
        end do
     end do
     call ft1axis(1,sig1,m1,n2,cmodl)
     tmp=real(cmodl)
     modl(:,:,jj)=tmp
  end do

  call to_history("n1",m1     )
  call to_history("n2",n2     )
  call to_history("o1",o1     )
  call to_history("o2",o2     )
  call hclose()
  call sep_write(modl)
  deallocate(data,cdata,modl,cmodl)
end program UnLMO
