program ref_2_dip
  use sep3d_struct_mod
  use sep
  type(sep3d) :: slow_s,ref_s
  real, allocatable :: refs(:,:,:),dip(:,:)
  logical, allocatable :: mask(:,:)
  integer :: iref,ipt,i,ibeg,iend,ia1,ia2
  real    :: f,x
  logical :: found

  call sep_init(SOURCE)
  call init_sep3d("in",slow_s,"INPUT")
  call init_sep3d("refs",ref_s,"INPUT")
  allocate(refs(2,ref_s%n(1),ref_s%n(2)))
  allocate(dip(slow_s%n(1),slow_s%n(2)))
  allocate(mask(slow_s%n(1),slow_s%n(2)))
  call sreed("refs",refs,size(refs)*4)
  
  !interpolate reflectors to slowness size
  mask=.false.
  dip=-1.
  do iref=1,ref_s%n(2)
    do ipt=1,slow_s%n(2)
      x=(ipt-1)*slow_s%d(2)+slow_s%o(2)
      if(x <= refs(1,1,iref)) then 
        ibeg=1
      else if(x >= refs(1,ref_s%n(1),iref)) then 
        ibeg=ref_s%n(1)-1
      else
        i=1; found=.false.
        do while(i< ref_s%n(1) .and. .not. found) 
          if(x>=refs(1,i,iref) .and. x < refs(1,i+1,iref)) then
            found=.true.
          else
            i=i+1
          end if
        end do
        if (.not. found) then
          write(0,*) "check",x  ,iref
          call seperr("internal error")
        end if
        ibeg=i
      end if
      i=(refs(2,ibeg,iref)-slow_s%o(1))/slow_s%d(1)+1.5
      mask(i,ipt)=.true.
      dip(i,ipt)=(refs(2,ibeg+1,iref)-refs(2,ibeg,iref))/&
          (refs(1,ibeg+1,iref)-refs(1,ibeg,iref))
    end do
  end do

  !interpolate the depths
  do i2=1,slow_s%n(2)
    ibeg=1
    do while(.not. mask(ibeg,i2)) 
       ibeg=ibeg+1
    end do
    dip(1:ibeg,i2)=dip(ibeg,i2)
    iend=slow_s%n(2)
    iend=51
    do while(.not.mask(iend,i2)) 
      iend=iend-1
    end do
    dip(iend:,i2)=dip(iend,i2)
    ia1=ibeg
    do while(ia1 /= iend) 
      ia2=ia1+1
      do while(.not.mask(ia2,i2).and. ia2 < iend) 
        ia2=ia2+1
      end do
      do i=ia1,ia2
        f=real(i-ia1)/real(ia2-ia1)
        dip(i,i2)=(1.-f)*dip(ia1,i2)+f*dip(ia2,i2)
      end do 
      ia1=ia2
    end do
  end do
  call srite("out",dip,size(dip)*4)
end program
