program Rat
  use sep

  implicit none
  integer                                    :: n1, n2, w1, w2, i1, i2
  integer                                    :: x1, x2, y1, y2
  real, dimension (:,:), allocatable, target :: denom, numer
  real, dimension (:,:), allocatable         :: ratio
  real, dimension (:,:), pointer             :: dwind, nwind
  real                                       :: over, under

  call sep_init ()
  call from_history (n1, n2)
  call from_par ("w1", w1)
  call from_par ("w2", w2)
  call sep_close ()

  allocate (denom (n1, n2), numer (n1, n2), ratio (n1, n2))
  
  call sep_read (numer)
  call sep_read (denom, file="denom")

  ratio = 0.
  do i2 = 1, n2
     if (i2 <= w2/2) then
        y1 = 1
        y2 = 1 + w2
     else if (i2 >= n2 - w2/2) then
        y1 = n2 - w2
        y2 = n2
     else
        y1 = i2 - w2/2
        y2 = i2 + w2/2
     end if
     do i1 = 1, n1        
        if (i1 <= w1/2) then
           x1 = 1
           x2 = 1 + w1
        else if (i1 >= n1 - w1/2) then
           x1 = n1 - w1
           x2 = n1
        else
           x1 = i1 - w1/2
           x2 = i1 + w1/2
        end if
        dwind => denom (x1:x2,y1:y2)
        nwind => numer (x1:x2,y1:y2)
        
        over  = sum (dwind*nwind)
        under = sum (dwind*dwind)
        if (abs(over)  > epsilon (over)   &
        .and.    under > epsilon (under)) &
        ratio (i1,i2) = -180/3.1415926536*atan(over / under)
     end do
  end do

  call sep_write (ratio)

  deallocate (denom, numer, ratio)
end program Rat
