program RCabcd2_
  use sep
  use util_

  implicit none

  complex, pointer :: rays(:,:)
  real,    pointer :: x(:,:),z(:,:)
  real,    pointer :: h1(:,:),h2(:,:)
  real,    pointer :: m1(:,:),m2(:,:)
  real,    pointer :: a(:,:),b(:,:),c(:,:),d(:,:)

  type(myaxis) ::at,ag
  integer      ::it,ig,j
  real         ::gx,gz

  real         :: eps,peps

  call sep_init()

  call hinaxis(at,1)
  call hinaxis(ag,2)

  call to_history("esize",4)
  call to_history("n3",4)

  call from_param("peps",peps,0.1)
  call sep_close()

  allocate(rays(at%n,ag%n))
  allocate(  h1 (at%n,ag%n),h2 (at%n,ag%n))
  allocate(  m1 (at%n,ag%n),m2 (at%n,ag%n))
  allocate(   x(at%n,ag%n))
  allocate(   z(at%n,ag%n))
  allocate(   a(at%n,ag%n),b(at%n,ag%n),c(at%n,ag%n),d(at%n,ag%n))

  call sep_read(rays)
  z=aimag(rays)
  x= real(rays)

  !! h1
  do it=1,at%n-1
     do ig=1,ag%n
        gx = (x(it+1,ig)-x(it,ig))/at%d
        gz = (z(it+1,ig)-z(it,ig))/at%d
        h1(it,ig) = sqrt(gx*gx+gz*gz)
     end do
  end do
  h1(at%n,:) = h1(at%n-1,:)

  !! h2
  do it=1,at%n
     do ig=1+1,ag%n-1
        gx = (x(it,ig+1)-x(it,ig-1))/2/ag%d
        gz = (z(it,ig+1)-z(it,ig-1))/2/ag%d
        h2(it,ig) = sqrt(gx*gx+gz*gz)
     end do
  end do
  h2(:,1)    = h2(:,2)
  h2(:,ag%n) = h2(:,ag%n-1)
  h2(1,:)    = h2(2,:)

  !! m1 = d(h2/h1)/dt = ( d(h2)/dt * h1 - h2 * d(h1)/dt )/ h1^2 
  do it=1,at%n-1
     do ig=1,ag%n
        m1(it,ig) &
        = h1(it,ig)* (h2(it+1,ig)-h2(it,ig))/at%d &
        - h2(it,ig)* (h1(it+1,ig)-h1(it,ig))/at%d
     end do
  end do
  m1(at%n,:) = m1(at%n-1,:)

  !! m2 = d(h1/h2)/dg = ( d(h1)/dt * h2 - h1 * d(h2)/dt )/ h2^2 
  do it=1,at%n
     do ig=1+1,ag%n-1
        m2(it,ig) &
        = h2(it,ig)* (h1(it,ig+1)-h1(it,ig-1))/2/ag%d &
        - h1(it,ig)* (h2(it,ig+1)-h2(it,ig-1))/2/ag%d
     end do
  end do
  m2(:,1)    = m2(:,2)
  m2(:,ag%n) = m2(:,ag%n-1)

!!$  eps=peps*(maxval(abs(h1))+minval(abs(h1)))/2.
!!$  where(abs(h1)<eps)
!!$     h1=eps
!!$!     h1=h1+eps-abs(h1)
!!$  end where
  
  eps=peps*(maxval(abs(h2))+minval(abs(h2)))/2.
  where(abs(h2)<eps)
     h2=eps
!     h2=h2+eps-abs(h2)
  end where

  a= 1/h1**2
  d= 1/h2**2

!  where(h1/=0.)
!     a= 1/h1**2
!  end where
!  where(h2/=0.)
!     d= 1/h2**2
!  end where


!  where(h1/=0. .and. h2/=0.)
     b=m1*a/h1/h2
!  end where
!  where(h1/=0. .and. h2/=0.)
     c=m2*d/h1/h2
!  end where

  call sep_write(a)
  call sep_write(b)
  call sep_write(c)
  call sep_write(d)

  call exit(0)
end program RCabcd2_









