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'