c     roi -- patch mode interferometric processor with prati filter 
c     
      program roi

      integer*4    nnn, mmm, overlap, ndiv, nbytesmax, ranfftiq, ranfftov
      integer*4    ranfft
      parameter    (nnn = 4096,mmm = 5120, overlap = 0)
      parameter    (nbytesmax = 16384, ndiv = 4)
      parameter    (ranfftov = 16384, ranfftiq = 8192)

      complex*8    trans1(nnn,mmm),trans2(nnn,mmm) 
      complex*8    data1(nnn),data2(nnn)
      complex*8    outarray(nnn/4,mmm),outline(mmm)
      complex*8    amps(nnn/4,mmm),outlinep(mmm)
      complex*8    ref1(ranfftov+1),ref2(ranfftov+1)
      complex*8    ref(ranfftov+1), refn2(ranfftov+1)
      real*4       unpacki1(0:255),unpackq1(0:255)
      real*4       unpacki2(0:255),unpackq2(0:255)
      real*4	   azres
      real*4       spec(256,3),kx01, kx02
      real*8       fd1, fdd1, fddd1, fd2, fdd2, fddd2
      real*8       r001, r002, r01, r02, r, delr, ht1, ht2, wavl, re, gm
      real*8       slpr, intr, slpa, inta, coefa(2)
      real*8       sloper, interr, slopea, intera
      real*8       dsloper, dinterr, dslopea, dintera
      character*64 name_in_1,name_in_2,name_out,name_outp, name_8lk 
      character*1  deskew, iqflip, srm
      integer*4    ideskew, iflip
      integer*4    FDSC1, FDSC2, FDSC3, FDSC4, FDSC5, FDSC6
      data         gm/0.398600485042958863d15/

      save         trans1,trans2,outarray,amps

c     **********  input parameters  **********
c     
c     enter input parameters
c     
      print *
     $     ,'ROI - Repeat Orbit Interferometric processor'
      print *
c     
      print *,'First input data file ? '
      read(*,'(a)')name_in_1
      print *,'Second input data file ? '
      read(*,'(a)')name_in_2
      print *,'Output data file ? '
      read(*,'(a)')name_out
      print *,'Output amplitudes file ? '
      read(*,'(a)')name_outp
      print *,'8lk output file ? '
      read(*,'(a)')name_8lk
      print *,'debug flag? (0) '
      read(*,*)iflag
      print *,'How many input bytes per line ?  '
      read(*,*)nbytes
      print *,'How many good bytes per line, including header ? '
      read(*,*)ngood
      print *,'First line to read ? (start at 0) '
      read(*,*)ifirstline
      print '(a,$)',' enter # of range input patches: '
      read(*,*)npatches
      print '(a,$)',' First sample pair to use (start at zero; 206 ERS) ? ' 
      read(*,*)ifirst
      print '(a,$)',' Number of valid points in azimuth ? '
      read(*,*)na_valid
      print *,'Deskew the image ? '
      read(*,'(a1)')deskew
      print *,'Caltone % of sample rate) ? '
      read(*,*) caltone1, caltone2
      print *,'Start range bin, number of range bins to process ? '
      read(*,*) isave, nlinesaz
      print '(a,$)',' Delta azimuth, range pixels for second file: '
      read(*,*)iazdelta,iradelta
      print *, 'Image 1 Doppler centroid quad coefs (Hz/prf)? '
      read(*,*) fd1, fdd1, fddd1
      print *, 'Image 2 Doppler centroid quad coefs (Hz/prf)? '
      read(*,*) fd2, fdd2, fddd2
      print *, '1 = use file 1 doppler, 2 = file 2, 3 = avg '
      read(*,*) iusedopp
      print *, 'Earth Radius (m)? '
      read(*,*) re
      print *, 'Body Fixed S/C velocities 1,2 (m/s)? '
      read(*,*) vel1, vel2
      print *, 'Spacecraft height 1,2 (m)? '
      read(*,*) ht1, ht2
      print *, 'Range of first sample in raw data file 1,2 (m)? '
      read(*,*) r001, r002
      print *, 'PRF 1,2 (pps)? '
      read(*,*) prf1, prf2
      print *, 'i/q means, i1,q1, i2,q2? '
      read(*,*) xmi1, xmq1, xmi2, xmq2
      print *, 'Flip i/q or offset video (y/n flip, o offset video)? '
      read(*,'(a)')iqflip
      print *, 'Desired azimuth resolution (m)? '
      read(*,*) azres
      print *, 'Number of azimuth looks? '
      read(*,*) nlooks
      print *, 'Range sampling rate (Hz)? '
      read(*,*) fs
      print *, 'Chirp Slope (Hz/s)? '
      read(*,*) slope
      print *, 'Pulse Duration (s)? '
      read(*,*) pulsedur
      print *, 'Chirp extension points? '
      read(*,*) nextend
      print *, 'Secondary range migration correction (y/n)? '
      read(*,'(a)') srm
      print *, 'Radar Wavelength (m)? '
      read(*,*) wavl
      print *, 'Range Spectral Weighting (1.=none, 0.54=Hamming) ? '
      read(*,*) rhww
      print *,'For the following, positive saves first file low freqs' 
      print '(a,$)',' Fraction of range bandwidth to remove: '
      read(*,*)pctbw, pctbwaz
      print *, 'linear resampling coefs:  sloper, intr, slopea, inta'
      read(*,*)sloper,interr,slopea, intera
      print *, 'linear resampling deltas: dsloper, dintr, dslopea, dinta'
      read(*,*)dsloper,dinterr,dslopea, dintera

c     make nlinesaz at most mmm


      nlinesaz = min (nlinesaz,mmm)

c     determine whether iq or offset video
      if(iqflip.eq.'O')iqflip='o'
      if(iqflip.eq.'o')then
         print *,'Offset video format assumed.'
         ranfft=ranfftov
      else
         print *,'I/Q sampling format assumed.'
         ranfft=ranfftiq
      end if

c     compute the azimuth parameters

      dxsamp1 = vel1/prf1
      dxsamp2 = vel2/prf2
      print *, 'actual orbital pulse spacing ', dxsamp1, dxsamp2

c     need to compute effective velocity for chirp rate due to orbit
c     curvature.

      vel1 = vel1 * sqrt(re/(re+ht1))
      vel2 = vel2 * sqrt(re/(re+ht2))

      dxsamp1 = vel1/prf1
      dxsamp2 = vel2/prf2
      delr=299792456./fs/2.

      print *
      print *,'First  input data file ',name_in_1
      print *,'Second input data file ',name_in_2
      print *,'Output interferogram   ',name_out
      print *,'Output amplitudes      ',name_outp
      print *,'First line to read in file 1 (start at 0) ', ifirstline 
      print *,'# of range input patches ', npatches
      print *,'First sample pair to use ',ifirst
      print *,'Number of valid points in azimuth ',na_valid
      print *,'Deskew the image? ',deskew
      print *,'First range bin to save in file 1 ', isave
      print *,'Number of range bins to process ', nlinesaz
      print *,'Delta azimuth for second file ',iazdelta
      print *,'Delta range pixels for second file: ',iradelta
      print *, 'Caltone location 1,2 ',caltone1, caltone2
      print *, 'Doppler centroid quad coef 1 (Hz/prf) ', fd1, fdd1, fddd1
      print *, 'Doppler centroid quad coef 2 (Hz/prf) ', fd2, fdd2, fddd2
      print *, 'Using Doppler flag ', iusedopp
      print *, 'Effective S/C Body fixed velocity 1 (m/s) ', vel1
      print *, 'Effective S/C Body fixed velocity 2 (m/s) ', vel2
      print *, 'Effective Azimuth sample spacing 1 (m) ', dxsamp1
      print *, 'Effective Azimuth sample spacing 2 (m) ', dxsamp2
      print *, 'Earth Radius (m) ', re
      print *, 'Spacecraft height 1 (m) ', ht1
      print *, 'Spacecraft height 2 (m) ', ht2
      print *, 'Range of first pixel in range compressed file 1 (m) ',
     $     r001 
      print *, 'Range of first pixel in range compressed file 2 (m) ',
     $     r002 
      print *, 'PRF 1 (pps) ', prf1
      print *, 'PRF 2 (pps) ', prf2
      print *, 'i/q means, i1,q1, i2,q2 ', xmi1, xmq1, xmi2, xmq2
      print *, 'Flip i/q ', iqflip
      print *, 'Desired azimuth resolution (m) ', azres
      print *, 'Number azimuth looks (m) ', nlooks
      print *, 'Range sampling rate (Hz) ', fs
      print *, 'Chirp Slope (Hz/s) ', slope
      print *, 'Pulse Duration (s) ', pulsedur
      print *, 'Chirp extension ', nextend
      print *, 'Secondary Range Correction ', srm
      print *, 'Radar Wavelength (m) ', wavl
      print *, 'Range Spectral Weighting ', rhww
      print *, 'Fractional bandwidth to remove ', pctbw, pctbwaz
      print *, 'linear resampling coefs:  sloper, intr, slopea, inta',
     .     sloper,interr,slopea, intera
      print *, 'linear resampling deltas: dsloper, dintr, dslopea, dinta'
     $     ,dsloper,dinterr,dslopea, dintera
      print *

      t0 = secnds(0.0)          ! start timer

      pi = 4.*atan(1.)
      pi2 = 2. * pi
      theta1=0.
      theta2=0.

c     compute range reference function

      fwidth=slope*pulsedur
      if(iqflip.eq.'o')then
         npts=2.*fs*pulsedur
         ts=1./(2.*fs)
      else
         npts=fs*pulsedur
         ts=1./fs
      end if
      if(mod(npts,2) .eq. 0) npts=npts+1
      print *,'Pulse length in points, time/sample in us: ',npts,ts*1.e6
      fd1 = fd1 * prf1
      fdd1 = fdd1 * prf1
      fddd1 = fddd1 * prf1
      fd2 = fd2 * prf2
      fdd2 = fdd2 * prf2
      fddd2 = fddd2 * prf2

      print *, 'Doppler coefficients referenced to range bin index: '
      print *, 'fd (Hz)              ', fd1, fd2
      print *, 'd fd/dr (Hz/s)       ', fdd1, fdd2
      print *, 'd^2 fd/dr^2 (Hz/s/s) ', fddd1,fddd2
      
c     modify dopplers as desired in iusedopp

      call moddopp(fd1,fdd1,fddd1,fd2,fdd2,fddd2,iusedopp,iradelta)

c     reference doppler coefficients to range

      call radopp(fd1,fdd1,fddd1,r001,delr)
      call radopp(fd2,fdd2,fddd2,r002,delr)

      print *, 'Doppler coefficients referenced to absolute range: '
      print *, 'fd (Hz)              ', fd1, fd2
      print *, 'd fd/dr (Hz/s)       ', fdd1, fdd2
      print *, 'd^2 fd/dr^2 (Hz/s/s) ', fddd1,fddd2
      
      if(iusedopp .eq. 0) then
         print *, 'Doppler check: ',fd1+fdd1*r001+fddd1*r001*r001, 
     .        fd2+fdd2*r002+fddd2*r002*r002 
      else
         r = r001 + float(iradelta)*delr
         print *, 'Doppler check: ',fd1+fdd1*r001+fddd1*r001*r001, 
     .        fd2+fdd2*r+fddd2*r*r 
      endif

c     update r001 and r002 to reflect the offsets given in isave,
c     iradelta, and nextend

      
      r01 = r001 + (isave-nextend-1) * delr
      r02 = r002 + (isave-nextend-1+iradelta) * delr

      print *, 'updated start ranges (m) ',r01, r02

      r1  = r01 + delr*(nlinesaz-1)
      r2  = r02 + delr*(nlinesaz-1)
      print *,'far range 1,2 (m) ', r1,r2
      npfin1 = r1*wavl/(2.0*azres*dxsamp1)+2
      npfin2 = r2*wavl/(2.0*azres*dxsamp2)+2
      print *,'filter points in far range ', npfin1,npfin2

      a21 = -2.0*pi/(dxsamp1*float(nnn)) ! for deskew
      a22 = -2.0*pi/(dxsamp2*float(nnn)) ! for deskew
      a41 = wavl/(2.0*azres*dxsamp1) ! for np
      a42 = wavl/(2.0*azres*dxsamp2) ! for np
      print *,'coefficient for filter points a4 ',a41, a42
      print *, 'near range chirp rate ', -2. * vel1**2/(wavl*r01), -2.
     $     * vel2**2/(wavl*r02) 
      print *, 'far  range chirp rate ', -2. * vel1**2/(wavl*r1), -2.
     $     * vel2**2/(wavl*r2) 
      
      theta1 = theta1*pi/180.0  !use only if deskew desired
      theta2 = theta2*pi/180.0  !use only if deskew desired

      rhww1 = 1.0-rhww

      k=0
      if(pctbw.ge.0.0)then
         k1start=abs(pctbw)*npts
         k1end=npts
         k2start=0
         k2end=npts-abs(pctbw)*npts
      else
         k1start=0
         k1end=npts-abs(pctbw)*npts
         k2start=abs(pctbw)*npts
         k2end=npts
      end if
      
      print *,'npts, k1 start, end, k2 start, end'
      print *,npts,k1start,k1end,k2start,k2end
      
      do i=-npts/2,npts/2
         k=k+1
         t=i*ts
         if(iqflip.eq.'o')then
            phase = pi*slope*t*t+pi*fs*t
            if(k.ge.k1start.and.k.le.k1end)
     +           ref1(i+npts/2+1)=cmplx(cos(phase),0.)
            if(k.ge.k2start.and.k.le.k2end)
     +           ref2(i+npts/2+1)=cmplx(cos(phase),0.)
            ref(i+npts/2+1)=cmplx(cos(phase),0.)
         else
            phase = pi*slope*t*t
            if(k.ge.k1start.and.k.le.k1end)
     +           ref1(i+npts/2+1)=cmplx(cos(phase),sin(phase))
            if(k.ge.k2start.and.k.le.k2end)
     +           ref2(i+npts/2+1)=cmplx(cos(phase),sin(phase))
            ref(i+npts/2+1)=cmplx(cos(phase),sin(phase))
         end if
      end do
      
      do i=1,npts
         win1 = (rhww-rhww1*cos((2.*pi*float(i-1))/float(npts-1)))
         ref(i)=ref(i)*win1
         ref1(i)=ref1(i)*win1
         ref2(i)=ref2(i)*win1
      end do

      
      if(nextend .gt. 0) then
         do i = 1 , npts
            k = i - nextend
            if(k .le. 0) k = k + ranfft
            ref(k) = ref(i)
            ref1(k) = ref1(i)
            ref2(k) = ref2(i)
         end do
         do i = 1 , nextend	
       	    k = npts - nextend+i
            if(k .le. 0) k = k + ranfft
            ref(k) = cmplx(0.,0.)
            ref1(k) = cmplx(0.,0.)
            ref2(k) = cmplx(0.,0.)
         end do
      end if
      if(nextend .lt. 0) then
         do i = npts, 1, -1
            k = i - nextend
            if(k .gt. ranfft) k = k - ranfft
            ref(k) = ref(i)
            ref1(k) = ref1(i)
            ref2(k) = ref2(i)
         end do
         do i = 1 , -nextend
            ref(i) = cmplx(0.,0.)
            ref1(i) = cmplx(0.,0.)
            ref2(i) = cmplx(0.,0.)
         end do
      end if

      print *,'reference calculated in time domain.'
      
c     init transform lengths

      call cfft1d(nnn,ref,0)
      call cfft1d(ranfft/2,ref,0)
      call cfft1d(ranfft,ref,0)
      do k=5,16
         len=2**k
         call cfft1d(len,ref,0)
      end do

c     calculate fft of range reference function

      call cfft1d(ranfft,ref,-1)
      call cfft1d(ranfft,ref1,-1)
      call cfft1d(ranfft,ref2,-1)
c     zero out dc and caltone location
      icaltone1=nint(caltone1*ranfft)
      icaltone2=nint(caltone2*ranfft)
      print *,'Caltone bins: ',icaltone1, icaltone2

      do i = 1, 6
         wgt = 0.5 - 0.5 * cos((i-1)/5.*pi)
         ref(i) = ref(i) * wgt
         ref1(i) = ref1(i) * wgt
         ref2(i) = ref2(i) * wgt
         if(iqflip.eq.'o')then
            ref(i+ranfft/2) = ref(i+ranfft/2) * wgt
            ref1(i+ranfft/2) = ref1(i+ranfft/2) * wgt
            ref2(i+ranfft/2) = ref2(i+ranfft/2) * wgt
            ref(ranfft/2+1-i) = ref(ranfft/2+1-i) * wgt
            ref1(ranfft/2+1-i) = ref1(ranfft/2+1-i) * wgt
            ref2(ranfft/2+1-i) = ref2(ranfft/2+1-i) * wgt
         end if
         ref(i+icaltone1) = ref(i+icaltone1) * wgt
         ref1(i+icaltone1) = ref1(i+icaltone1) * wgt
         ref2(i+icaltone2) = ref2(i+icaltone2) * wgt
         k=icaltone1+1-i
         k2=icaltone2+1-i
         if(k.le.0)k=k+ranfft
         if(k2.le.0)k2=k2+ranfft
         ref(k) = ref(k) * wgt
         ref1(k) = ref1(k) * wgt
         ref2(k2) = ref2(k2) * wgt
         if(iqflip.eq.'o')then
            k=i+ranfft-icaltone1
            k2=i+ranfft-icaltone2
            if(k.gt.ranfft)k=k-ranfft
            if(k2.gt.ranfft)k2=k2-ranfft
            ref(k) = ref(k) * wgt
            ref1(k) = ref1(k) * wgt
            ref2(k2) = ref2(k2) * wgt
            ref(ranfft-icaltone1+1-i) = ref(ranfft-icaltone1+1-i) * wgt
            ref1(ranfft-icaltone1+1-i) = ref1(ranfft-icaltone1+1-i) * wgt
            ref2(ranfft-icaltone2+1-i) = ref2(ranfft-icaltone2+1-i) * wgt
         end if
         ref(ranfft+1-i) = ref(ranfft+1-i) * wgt
         ref1(ranfft+1-i) = ref1(ranfft+1-i) * wgt
         ref2(ranfft+1-i) = ref2(ranfft+1-i) * wgt
      end do

c     refn2 is ref without the secondary range correction phase
c     here for debugging

      do i = 1 , ranfft
         refn2(i) = ref(i)
      end do

      if(srm .eq. 'y' .or. srm .eq. 'Y') then

c     include secondary range migration correction to the
c     chirp transform

         delw = pi2 * fs / ranfft

         rc1 = r01 + delr*nlinesaz/2
         rc2 = r02 + delr*nlinesaz/2
         kx01 = pi2*(fd1 + fdd1*rc1 + fddd1*rc1**2)/vel1
         kx02 = pi2*(fd2 + fdd2*rc2 + fddd2*rc2**2)/vel2
         omega0 = pi2 * 299792456. / wavl
         print *, 'delw =    ',delw
         print *, 'rc =      ',rc1,rc2
         print *, 'kx0 =     ',kx01,kx02
         print *, 'omega0 =  ',omega0
         phase1 = - 0.25 * rc1 * (wavl / pi2) * 
     .        kx01**2 * ((omega0+ranfft/2.*delw)/omega0)**2 
         phase2 = - 0.25 * rc2 * (wavl / pi2) * 
     .        kx02**2 * ((omega0+ranfft/2.*delw)/omega0)**2 
         print *, 'phaser =   ',phase1, phase2
         phase1 = - 0.25 * rc1 * (wavl / pi2) * kx01**2
         phase2 = - 0.25 * rc2 * (wavl / pi2) * kx02**2
         print *, 'phase0 =   ',phase1,phase2
         phase1 = - 0.25 * rc1 * (wavl / pi2) * 
     .        kx01**2 * ((omega0-ranfft/2.*delw)/omega0)**2 
         phase2 = - 0.25 * rc2 * (wavl / pi2) * 
     .        kx02**2 * ((omega0-ranfft/2.*delw)/omega0)**2 
         print *, 'phasel =   ',phase1,phase2
         
         do i = 0 , ranfft/2
            omega = omega0 + (i-4096)* delw
            phase1 = - 0.25 * rc1 * (wavl / pi2) * 
     .           kx01**2 * (omega/omega0)**2 
            phase2 = - 0.25 * rc2 * (wavl / pi2) * 
     .           kx02**2 * (omega/omega0)**2 
            ref (i+1) = ref (i+1) * cmplx(cos(phase1),sin(phase1))
            ref1(i+1) = ref1(i+1) * cmplx(cos(phase1),sin(phase1))
            ref2(i+1) = ref2(i+1) * cmplx(cos(phase2),sin(phase2))
         end do
         do i = 0 , ranfft/2
            omega = omega0 - i* delw
            phase1 = - 0.25 * rc1 * (wavl / pi2) * 
     .           kx01**2 * (omega/omega0)**2 
            phase2 = - 0.25 * rc2 * (wavl / pi2) * 
     .           kx02**2 * (omega/omega0)**2 
            ref (ranfft+1-i) =ref (ranfft+1-i) * cmplx(cos(phase1),sin(phase1))
            ref1(ranfft+1-i) =ref1(ranfft+1-i) * cmplx(cos(phase1),sin(phase1))
            ref2(ranfft+1-i) =ref2(ranfft+1-i) * cmplx(cos(phase2),sin(phase2))
         end do

      end if

c     scale reference for channel gain, conjugate

      gcal=1./ranfft
      do i=1,ranfft
         refn2(i)=conjg(refn2(i))*gcal
         ref(i)=conjg(ref(i))*gcal
         ref1(i)=conjg(ref1(i))*gcal
         ref2(i)=conjg(ref2(i))*gcal
      end do

c     save spectra of reference functions for checking

      if(iand(iflag,16) .eq. 16) then
         do i=1,256
            sum1=0.
            sum2=0.
            sum3=0.
            do k=0,31
               sum1=sum1+cabs(ref(i*32-31+k))**2
               sum2=sum2+cabs(ref1(i*32-31+k))**2
               sum3=sum3+cabs(ref2(i*32-31+k))**2
            end do
            spec(i,1)=sum1
            spec(i,2)=sum2
            spec(i,3)=sum3
         end do
         open(31,file='rangespecs.out')
         do i=1,256
            write(31,*) i,(spec(i,k),k=1,3)
         end do
         close(31)
         
         open(31,file='fullrangespecs.out')
         do i=1,ranfft/2
            write(31,*) i,cabs(ref(i)),atan2(aimag(ref(i)),real(ref(i)))
     $           ,cabs(refn2(i)),atan2(aimag(refn2(i)),real(refn2(i))) 
         end do
         close(31)
      end if
c     
c     open input files
c     
      FDSC1 = initdk(lun,name_in_1)
      FDSC2 = initdk(lun,name_in_2)
c     
c     open output files
c     
      fdsc4 = initdk(lun,name_out)
      fdsc5 = initdk(lun,name_outp)
      if(iand(iflag,16) .eq. 16) fdsc6 = initdk(lun,name_8lk)

      print *,'Files opened.'

c     offset into valid data in a patch

      ifs = (nnn-na_valid)/2

      iflip = 0
      if(iqflip .eq. 'Y' .or. iqflip .eq. 'y') iflip = 1

c     load the unpacking array

      if(iflip .eq. 0) then
         do i=0,255
            unpacki1(i)=float(i)-xmi1
            unpackq1(i)=float(i)-xmq1
            unpacki2(i)=float(i)-xmi2
            unpackq2(i)=float(i)-xmq2
         end do
      else
         do i=0,255
            unpacki1(i)=float(i)-xmq1
            unpackq1(i)=float(i)-xmi1
            unpacki2(i)=float(i)-xmq2
            unpackq2(i)=float(i)-xmi2
         end do
      end if

      ideskew = 0
      if(deskew .eq. 'y' .or. deskew .eq. 'Y') ideskew = 1

c     
c     begin loop to range process data
c     
      
      do ipatch=1,npatches
         
c     get offsets for the given patch
         
         slpr = sloper+(ipatch-1)*dsloper
         intr = interr+(ipatch-1)*dinterr
         slpa = slopea+(ipatch-1)*dslopea
         inta = intera+(ipatch-1)*dintera
         print *,'Patch ',ipatch,' resampling inputs: ',slpr,intr,slpa
     $        ,inta

         coefa(1) = slpa
         coefa(2) = inta

c     convert to function of range instead of pixel

         intr = intr - slpr*r001/delr
         slpr = slpr/delr

         inta = inta - slpa*r001/delr
         slpa = slpa/delr

         print *, 'Converted Resampling inputs ', slpr,intr,slpa,inta
c     
c     compress channels
         
         print *
         print *,' range compressing channels.....'
         
         print *
 100     format(' processing line: ',i5,f10.3,' seconds.')

c     do file 1

c     irec=loop+ifirstline+(ipatch-1)*na_valid
         irec=ifirstline+(ipatch-1)*na_valid
         ifrst= ifirst+isave-1

         if(iqflip.eq.'o')then
            call rcov(FDSC1,nnn,nlinesaz,trans1,unpacki1,unpackq1,
     $           ref1,irec,ifrst,nbytes,ngood)
         else
            call rciq(FDSC1,nnn,nlinesaz,trans1,unpacki1,unpackq1,
     $           ref1,irec,ifrst,nbytes,ngood,iflip)
         end if


c  write out rc data if iflag contains 32
         if(ipatch.eq.1 .and. (iand(iflag,32) .eq. 32))then
            print *, 'writing first channel range compressed data '
            FDSC3 =initdk(lun,'im.rc1 ')
            nwr = iowrit(FDSC3,trans1,8*nnn*nlinesaz)
            call closedk(lun,fdsc3)
         end if

         t1=secnds(t0)
         print *,' first channel range processing elapsed time: ',t1,'sec'

c     do the same for channel 2

         irec=ifirstline+iazdelta+(ipatch-1)*na_valid
         ifrst=ifirst+isave-1+iradelta

         if(iqflip.eq.'o')then
            call rcov(FDSC2,nnn,nlinesaz,trans2,unpacki2,unpackq2,
     $           ref2,irec,ifrst,nbytes,ngood)
         else
            call rciq(FDSC2,nnn,nlinesaz,trans2,unpacki2,unpackq2,
     $           ref2,irec,ifrst,nbytes,ngood,iflip)
         end if

c  write out rc data if iflag contains 32
         if(ipatch.eq.1 .and. (iand(iflag,32) .eq. 32))then
            print *, 'writing second channel range compressed data '
            FDSC3 =initdk(lun,'im.rc2 ')
            nwr = iowrit(FDSC3,trans2,8*nnn*nlinesaz)
            call closedk(lun,fdsc3)
         end if

         t1=secnds(t0)
         print *,' second channel range processing elapsed time: ',t1,'sec'


         t1=secnds(t0)
         print *,' range processing elapsed time: ',t1,'sec'

c     transform lines
         
         print *, 'transforming lines ',t1, ' sec'

         do i=1,nlinesaz
            call cfft1d(nnn,trans1(1,i),-1)
         end do

         do i=1,nlinesaz
            call cfft1d(nnn,trans2(1,i),-1)
         end do

         t1=secnds(t0)
         print *, 'transformed lines ',t1, ' sec'

         
         if(ipatch.eq.1 .and. (iand(iflag,8) .eq. 8))then
            print *, 'writing azimuth spectrum '
            FDSC3 =initdk(lun,'im.rcas1 ')
            nwr = iowrit(FDSC3,trans1,8*nnn*nlinesaz)
            call closedk(lun,fdsc3)
         end if

c     start the range migration correction
         
         print *,'start range migration correction.'
         
         nl=nlinesaz
         r=r01
         call RMpatch(trans1,0.d0,0.d0,
     .        nnn,nl,nls,r,delr,wavl,vel1,fd1,fdd1,fddd1,prf1,ideskew)



         nl = nlinesaz
         r = r02
         call RMpatch(trans2,
     .        slpr,intr,
     .        nnn,nl,nls,r,delr,wavl,vel2,fd2,fdd2,fddd2,prf2,ideskew)


         
         if(ipatch.eq.1 .and. (iand(iflag,4) .eq. 4))then
            print *, 'writing range migrated data '
            FDSC3 =initdk(lun,'im.rcrm1 ')
            nwr = iowrit(FDSC3,trans1,8*nnn*nlinesaz)
            call closedk(lun,fdsc3)
         end if
         if(ipatch.eq.2 .and. (iand(iflag,4) .eq. 4))then
            print *, 'writing range migrated data '
            FDSC3 =initdk(lun,'im.rcrm2 ')
            nwr = iowrit(FDSC3,trans1,8*nnn*nlinesaz)
            call closedk(lun,fdsc3)
         end if
         
         
c     multiply by reference and inverse transform lines
         
         t1=secnds(t0)
         print *, 'inverse transforming lines ',t1, ' sec'

         nl=nlinesaz
         r = r01
         call ACpatch(trans1,nnn,nl,r,delr,
     .        wavl,vel1,fd1,fdd1,fddd1,prf1,ht1,theta1,npfin1,
     .        a21,a41,0.d0,0.d0,ideskew,na_valid)
         
c     r = r02 + (line-1)*delr
         r = r02
         call ACpatch(trans2,nnn,nl,r,delr,
     .        wavl,vel2,fd2,fdd2,fddd2,prf2,ht2,theta2,npfin2,
     .        a22,a42,slpa,inta,ideskew,na_valid)
         t1=secnds(t0)
c     print *, line,nl, t1, 'sec '
c     end do

         print *,'Range-Doppler done',t1,' sec'
         
         if(ipatch.eq.1 .and. (iand(iflag,16) .eq. 16))then
            
c     take looks

            do line = 1 , nlinesaz, 2
               do i = 1 , na_valid/8
                  data1(i) = cmplx(0.,0.)
                  data2(i) = cmplx(0.,0.)
               end do
               do j = 1 , 2
                  do i=1,na_valid/8
                     power=0.
                     do k=0,7
                        power = power+
     .                       cabs(trans1(ifs+i*8-7+k,line+j-1))**2
                     end do
                     data1(i)=data1(i) + power
                     power=0
                     do k=0,7
                        power = power+
     .                       cabs(trans2(ifs+i*8-7+k,line+j-1))**2
                     end do
                     data2(i)=data2(i)+power
                  end do
               end do
               do i = 1 , na_valid/8
                  data1(i) = cmplx(sqrt(real(data1(i))),0.)
                  data2(i) = cmplx(sqrt(real(data2(i))),0.)
               end do
               nwr = iowrit(fdsc6,data1,na_valid)
               nwr = iowrit(fdsc6,data2,na_valid)
            end do
         end if

c     if first time through, save complex images
         
         if(ipatch.eq.1 .and. (iand(iflag,1) .eq. 1))then
            print *,'writing channel 1 range doppler image'
            fdsc3 = initdk(lun,'trans1.dat ')
            nwr = iowrit(fdsc3,trans1,nlinesaz*8*nnn)
            call closedk(lun,fdsc3)
         end if
         if(ipatch.eq.1 .and. (iand(iflag,2) .eq. 2))then
            print *,'writing channel 2 range doppler image'
            fdsc3 = initdk(lun,'trans2.dat ')
            nwr = iowrit(fdsc3,trans2,nlinesaz*8*nnn)
            call closedk(lun,fdsc3)
            print *, ' written ',nwr
         end if
         
c     interfere  the patches


c  if nlooks greater than 1, take looks

         if(nlooks.gt.1)then

         call multilook(trans1,trans2,nlinesaz,nnn,nlooks,
     .        outarray,amps)
         print *,' writing patch out...'
c     write in range line format
         do line=ifs/nlooks+1,(ifs+na_valid)/nlooks
            do k=1,nlinesaz
               outline(k)=outarray(line,k)
               outlinep(k)=amps(line,k)
            end do
            nwr = iowrit(fdsc4,outline,nlinesaz*8)
            nwr = iowrit(fdsc5,outlinep,nlinesaz*8)
         end do
         
         else

         print *,' writing single look data out...'

c  if bit = 64 set in flag, write out individual 1-look complex files
            if(iand(iflag,64).eq.0)then

            do line=ifs+1,(ifs+na_valid)
               do k=1,nlinesaz
                  outline(k)=trans1(line,k)*conjg(trans2(line,k))
                  outlinep(k)=cmplx(cabs(trans1(line,k)),cabs(trans2(line,k)))
               end do
               nwr = iowrit(fdsc4,outline,nlinesaz*8)
               nwr = iowrit(fdsc5,outlinep,nlinesaz*8)
            end do

            else

            do line=ifs+1,(ifs+na_valid)
               do k=1,nlinesaz
                  outline(k)=trans1(line,k)
                  outlinep(k)=trans2(line,k)
               end do
               nwr = iowrit(fdsc4,outline,nlinesaz*8)
               nwr = iowrit(fdsc5,outlinep,nlinesaz*8)
            end do

            end if

         end if

         






         
         
         t1=secnds(t0)
         print *,' patch finished, time =', t1, ' seconds.'
      end do                    !end patch loop
      
      t1=secnds(t0)
      print *,' elapsed time =', t1, ' seconds.'
      end
        
      subroutine moddopp(fd1, fdd1, fddd1, fd2, fdd2, fddd2,
     .      iusedopp, iradelta)

      integer*4 iusedopp, iradelta
      real*8 fd1, fdd1, fddd1, fd2, fdd2, fddd2
      real*8 temp1, temp2, temp3

      if(iusedopp .eq. 0) then
         print *, 'using doppler as is '
      elseif(iusedopp .eq. 1) then
         print *, 'referencing ch 2 dopp to ch 1'
         fd2 = fd1 - float(iradelta) * fdd1 + float(iradelta)**2
     .        * fddd1
         fdd2 = fdd1 - 2.d0 * float(iradelta) * fddd1
         fddd2 = fddd1
      elseif(iusedopp .eq. 2) then
         print *, 'referencing ch 1 dopp to ch 2'
         fd1 = fd2 + float(iradelta) * fdd2 + float(iradelta)**2
     .        * fddd2
         fdd1 = fdd2 + 2.d0 * float(iradelta) * fddd2
         fddd1 = fddd2
      elseif(iusedopp .eq. 3) then
         print *, 'averaging dopplers'
         temp1 = fd2 + float(iradelta) * fdd2 + float(iradelta)**2
     .        * fddd2
         temp2 = fdd2 + 2.d0 * float(iradelta) * fddd2
         temp3 = fddd2
         fd1   = (fd1 + temp1)/2.d0
         fdd1  = (fdd1 + temp2)/2.d0
         fddd1 = (fddd1 + temp3)/2.d0
         fd2 = fd1 - float(iradelta) * fdd1 + float(iradelta)**2
     .        * fddd1
         fdd2 =  fdd1 - 2.d0 * float(iradelta) * fddd1
         fddd2 = fddd1
      end if
      return
      end

      subroutine radopp(fd, fdd, fddd, r, del)

      real*8 	fd, fdd, fddd, r, del, temp1, temp2, temp3

      temp1 = fd - fdd * (r/del) + fddd * (r/del)**2
      temp2 = fdd/del - 2.d0 * fddd*(r/del)/del
      temp3 = fddd / del**2

      fd   = temp1
      fdd  = temp2
      fddd = temp3

      return
      end

      include 'multilook.f'
c      include 'resample.f'







