	PROGRAM SETUP
C
C=========================================================================
C
C	PROGRAM "SETUP" PROMPTS FOR PROCESSING PARAMETERS THAT ARE
C	USED TO SET ARRAY DIMENSIONS IN THE MAIN PROCESSING SUBROUTINE
C	CALLED "POLAR_REF_D".  DIFFERENT MACHINES WILL HAVE DIFFERENT
C	PROCESSING SIZE LIMITS WHICH CAN BE ADJUSTED BY CHANGING THE
C	VALUES IN THE "PARAMETER" STATEMENT.
C
C=========================================================================
C
	COMPLEX aBLOCK(512,512),bBLOCK(512,512)

	COMPLEX  dummy(512)

        real yr1(512), yr2(512)
        real yr3(512), yd(512)

        real total_time
        real npr_time(1), fft_time(1)
        real npr_trans(1), fft_trans(1)
        real total_npr, total_fft
        real read_time, write_time
      
        integer i,k
        integer loop
        real p_time, t_time
C        real timer_stop
C        external timer_stop

C
C       for direct access file
C          
C           Alpha/OSF-1 counts number of words, each complex = 2 words
C           Sun4/SunOS  counts number of bytes, each complex = 8 bytes
C
C

      numrow = 500
      numcol = 397

      print *, 'starting the program'

	OPEN(7, NAME = 'targs_4.dat', STATUS = 'UNKNOWN',
     1	RECL = numcol*2, ACCESS = 'DIRECT',
     2	FORM = 'UNFORMATTED', iostat=ios, err= 110)

      print *, 'input file opened'
      

	OPEN(8, NAME = 'image_4.dat', STATUS = 'UNKNOWN',
     1	RECL = 512*2, ACCESS = 'DIRECT',
     2	FORM = 'UNFORMATTED')

      print *, 'output file opened'

      call timer_start()

      do k=1,512
        dummy(k) = (0.0, 0.0)
      enddo

         do i= 1, numrow
           call DREADR(7,dummy, numcol, i)
           do k=1,numcol
             bblock(k,i) = dummy(k)
           enddo
         enddo

      read_time = timer_stop()

      call read_vec(yr1,yr2,yr3,yd)

      
      print *, 'input read'

      do loop = 1, 1
         
C        bblock = cblock

        print *, '*****   iteration ', loop, '  starting *****'

        print *, '      calling polar reformatter '

        call polar_ref(ablock,bblock,yr1,yr2,yr3,yd,p_time,t_time)

        npr_time(loop) = p_time

        npr_trans(loop) = t_time

        print *, '      calling 2D FFT '

        call fft2d(ablock,bblock,p_time,t_time)

        fft_time(loop) = p_time

        fft_trans(loop) = t_time

      enddo

      print *, 'writing output'

      call timer_start()

      do i= 1, 512
         do k=1, 512
           dummy(k) = bblock(k,i)
         enddo

         call DWRITR(8,dummy, 512, i)
      enddo

      write_time = timer_stop()

      print *, 'read time = ', read_time, ' write time=',write_time
      print *, ' '
      print *, ' iter  npr  npr_trans  fft fft_trans '
      print *, ' '

      total_npr = 0.0
      total_fft = 0.0
      total_time = 0.0

      do i=1,1
        print *, i, npr_time(i), npr_trans(i),fft_time(i),fft_trans(i)
        total_npr  = total_npr  + npr_time(i) + npr_trans(i)
        total_fft  = total_fft  + fft_time(i) + fft_trans(i)
      enddo

        total_time = total_time + total_npr + total_fft

        print *, ' '
        print *, ' total     ', total_npr, total_fft, total_time
        print *, ' '
        print *, ' avg   =   ',  total_time/1
        go to 999

110    print *, 'input open failed ios = ', ios
       go to 999
999    continue
        stop
        end


        subroutine polar_ref(ablock,bblock,yr1,yr2,yr3,yd,pt,tt)

	COMPLEX aBLOCK(512,512),bBLOCK(512,512)
        real yr1(512), yr2(512)
        real yr3(512), yd(512)
        real pt,tt

	COMPLEX  POLAR(512)
	COMPLEX  TEMP(512)
	COMPLEX  dummy(512)
	REAL     FLT_R(512*2),DELTAR(512)
	REAL     PHI(512),CR(512),PROSCAL(512)
	INTEGER  ISOR(512),IHOR(512),RECT_R,RECT_C,POL_R
	INTEGER  ISOC(512)

        real RX,RY,RZ
        integer nopol, TBLSIZ
        integer L,nq,ixmaxr,ixmaxc,ifx,ify
        real    q,RC,R0,brf_r,brf_c
        integer i,j,row
        real THETA,FS,F0,FDOT
        integer ISTEER

      pt = 0.0
      tt = 0.0

      call timer_start()

      call init_part1(RX,RY,RZ,pol_r,nopol,rect_r,rect_c)

      do i=1,512
         dummy (i) = (0.0,0.0)
      enddo

      do i=1,512
      do j=1,512
         ablock (i,j) = (0.0,0.0)
      enddo
      enddo



C	      PROJECT -   SETS PARAMETERS FOR PROJECTION INTO DESIRED
C	                  SLANT PLANE.  ALSO, IF BEAM STEERING IS
C	                  CHOSEN, FOLDS THIS INFORMATION INTO THE PROJECTION.


	CALL PROJECT(PHI,THETA,PROSCAL,DELTAR,nopol,VPMAG,ISTEER,
     *  FS,F0,FDOT,RX,RY,RZ,yr1,yr2,yr3,yd)

        call init_part2(L,q,RC,R0,nq,ixmaxr,ixmaxc,brf_r,brf_c,
     *         ifx,ify,pol_r)


	DELTAY=REAL(POL_R-1)/REAL(RECT_R-1)*VPMAG !SPACING 1 Y-DIM. OF RECT. GRID 
						!1 FRACTIONS OF SAMPLES.
C
C
C	SET AN UPPER LIMIT ON THE SIZE OF THE RANGE INTERPOLATION FILTER.
C

	PMAX=0.
	DO  I=1,500	!GET THE MAXIMUM PROSCAL(*) TO USE 1 LEN_R CALC.
	   PMAX=AMAX1(PMAX,PROSCAL(I))
	   IF(PMAX.EQ.PROSCAL(I))  IMAX=I
	END DO

	LEN_R=2*IXMAXR*DELTAY*BRF_R*PMAX+1 !GOOD FOR 2*IXMAXR*DELTAY*BRF<=100
				    !& THETA_MAX<=8 DEG

	PI=4.*ATAN(1.0)
	R1=R0


	DO  NUMBLK=1,500 ! 

C        pdo (numblk=1:500)
C        pinout bblock(:,numblk)

           I=NUMBLK	!KEEP TRACK OF POLAR REC. NO.           

               do ix1 = 1, pol_r
                  polar(ix1) = bblock(ix1,numblk)
               enddo

C            polar = bblock(:,numblk)


C	      APPLY A BEAM STEERING CORRECTION, IF REQUESTED.

           isteer = 0
	      IF(ISTEER.EQ.1) CALL HEIFER(POLAR,POL_R,RX,RY,RZ,FS,F0,FDOT)
C
	      PH=PHI(I)*PI/180.
C
C	      ADJUST THE RANGE SAMPLING DISTANCE TO REFLECT THE ANGLE
C	      OF THE POLAR RECORD.
C
	      DELTA=DELTAY/COS(PH)
      
	      XMAX=IXMAXR*DELTA*PROSCAL(I)*BRF_R
              

C        FLT_R(NQ,LEN_R) got set in MAKFILT()

	      CALL MAKFILT(FLT_R,LEN_R,NQ,BRF_R,DELTA,XMAX,TBLSIZ)

C
C	      COMPUTE THE QUANTIZED LOCATIONS FOR ONE ENTIRE RANGE
C	      RECORD, THE FILTER POINTERS, AND THE INPUT DATA POINTERS.
C

C       ISOR and IHOR got set in LOCATE_R()

	      CALL LOCATE_R(PHI(I),DELTAY,RECT_R,R1,DELTAR(I),ISOR,POL_R,
     *        IHOR,NQ,VPMAG,PROSCAL(I))


C        print *, 'step 3: numblk =  ', numblk
C
C	      PERFORM THE RANGE INTERPOLATION.
C

	      CALL INTERP_R(POLAR,POL_R,dummy,nopol,RECT_R,
     *        TBLSIZ,FLT_R,NQ,ISOR,IHOR,PHI(I),DELTAY,IXMAXR,BRF_R,
     *        PROSCAL(I))

              do ix1 = 1, 500
                bblock(ix1,numblk) = dummy(ix1)
              enddo

C               bblock(:,numblk) = dummy
	END DO

C        end pdo

        pt = timer_stop()

C
C	   TRANSPOSED
C


c     transpose

      call timer_start()

      call idxperm(ablock,bblock)

      tt = timer_stop()

      call timer_start()
C
C	SET UP FOR CROSS-RANGE INTERPOLATION.  DETERMINE THE DISTANCE FROM
C	AN X-AXIS ORIGIN ON THE RECTANGULAR GRID TO RANGE-CORRECTED SAMPLES
C	ON EACH HORIZONTAL LINE OF THE FULL GRID.
C
	PI=4.*ATAN(1.)

	TH=PI*(THETA/2.)/180.		!MAKE POS. THIS TIME SO XORIG IS POS.
	XORIG=(R1+POL_R-1)*VPMAG*TAN(TH)	!X-AXIS ORIGIN OF FULL RECT. GRID
					!THIS IS THE DISTANCE LEFT OF THE
					!CENTER LINE, 1 SAMPLES, FOR THE
					!DATA FARTHEST 1 RANGE.
	DELTAX=(2.*XORIG)/REAL(RECT_C-1) !SPACING 1 X-DIMENSION OF RECT. GRID 
					!1 FRACTIONS OF SAMPLES.
					!NOTE THAT 2*XORIG IS THE X-DIM LENGTH
					!OF THE FULL GRID. 
C
C	PROCESS ONE HORIZONTAL LINE OF THE FULL GRID AT A TIME.  REMEMBER,
C	WE HAVE JUST DONE A TRANSPOSE.  THERE ARE NOW "RECT_R" LINES AND
C	"NOPOL" COLUMNS.  THE RESULT WILL BE PUT INTO BACK INTO THE MEMORY
C	BLOCK UNTIL THE BLOCK IS FULL, THEN WRITTEN TO A DISK FILE.
C


	XMAX=IXMAXC*DELTAX*BRF_C

C
C	INITIALIZE TEMP FOR USE 1 ZERO PADDING.
C
	DO  K=1,IFX
	   TEMP(K)=(0.,0.)
	END DO
C

C
C	DECIDE WHERE TO ZERO-PAD SUBPATCH 1 X-DIMENSION (CRS-RNG)
C
C
C	ZERO-PAD SOME INITIAL RECORDS OF THE OUTPUT FILE.
C

C
C	CROSS-RANGE CORRECT THE SAR DATA WITHIN RANGE RECORDS BOUNDED BY
C	THE SUB-PATCH OF INTEREST.  FILL THE MAIN MEMORY BLOCK WITH
C	TRANSPOSED DATA, INTERPOLATE ALL CROSS-RANGE RECORDS, THEN WRITE
C	2 THE ENTIRE BLOCK.
C
           nstr = rect_c
           klen = rect_r
           l = rect_c

C
C	   PROCESS RECORDS 1 THE BLOCK.
C
C** CONCURRENT COMPILER DIRECTIVE - BEGINNING
           ibeg = 1
           iend = 397
C           nstr = 502
         
C	   pdo (row=ibeg:iend)
C           pinout ablock(:, row)
          
	   DO  II=IBEG,IEND
         
              row = II
              do ix1 = 1, 512
                dummy(ix1) = ablock(ix1,row)
              enddo

C                 dummy = ablock(:,row)

C              I = II
	      I = row
	      RNG=R1*VPMAG+(I-1)*DELTAY	!RANGE TO HORIZONTAL LINE OF INTEREST

	      CALL LOCATE_CR(CR,nopol,PHI,THETA,DELTAX,DELTAY,RNG,
     *        XORIG,ISOC,RECT_C)

C              if (i.eq. 1) then
C                print *, 'inside pdo'
C                do ix1 = 1, 512
C                  print *, 'dummy(',ix1,') =', dummy(ix1)
C                enddo
C              endif

	      CALL INTERP_CR(CR,nopol,TEMP,RECT_C,ISOC,
     *        DELTAX,dummy,BRF_C,XMAX)

C
C	      COPY RESULT BACK INTO BLOCK FOR LATER OUTPUT.
C

	      DO  K=1,RECT_C
	         aBLOCK(K,row)=TEMP(K)
	      END DO

C              ablock(:,row) = TEMP

	   END DO

C           end pdo

           pt = pt + timer_stop()
C
C	   WRITE 2 THE BLOCK TO FINAL FILE.
C

C

        return

	END
C
C


        SUBROUTINE init_part1(RX,RY,RZ,pol_r,nopol,rect_r,rect_c)
        real RX,RY,RZ
        integer pol_r,nopol,rect_r,rect_c

        RX=0.
	RY=0.
	RZ=0.

        pol_r = 397

        nopol = 500

        rect_r = 397
        rect_c = 500
      
        return
        end

        SUBROUTINE init_part2(L,q,RC,R0,nq,ixmaxr,ixmaxc,brf_r,brf_c,
     *         ifx,ify,pol_r)
        integer L,nq,ixmaxr,ixmaxc,ifx,ify,pol_r
        real    q,RC,R0,brf_r,brf_c
       
	L = 500

        q = 20.1342

	RC=Q*POL_R	!CONVERT TO PATCH CENTER DIST., SAMPLES
	R0=RC-REAL(POL_R)/2.

        nq = 8

        ixmaxr = 8  
        ixmaxc = 8  

        brf_r = 1.04
        brf_c = 1.04
  
        ifx = 500
        ify = 397

        return
        end


        subroutine read_vec(yr1,yr2,yr3,yd)
        real yr1(512),yr2(512),yr3(512),yd(512)
        integer J

        open(9,file='targs_4_vec.dat',iostat=ios, err = 997,
     *            status = 'unknown')

	DO I=1,500

          READ(9,*,END=999)J,yr1(i),yr2(i),yr3(i),yd(i)

C         if (i .le. 10)  then
C           print *, J,yr1(i),yr2(i),yr3(i),yd(i)
C         endif
       enddo

       go to 1000
997    continue
       print *, 'open vec file failed'
       go to 1000
999    continue
       print *, 'unexpected end of file'
1000   continue
       close(9)
       return
       end


	SUBROUTINE PROJECT(PHI,ANGLE,PROSCAL,DELTAR,nopol,VPMAG,ISTEER,
     *  FS,F0,FDOT,RX,RY,RZ,yr1,yr2,yr3,yd)

        integer nopol,ISTEER
        real ANGLE,FS,F0,FDOT
	REAL  R(3),PROSCAL(500),PHI(500)
	REAL  IPRIME(3),JPRIME(3),KPRIME(3),DELTAR(500)
	REAL  IDP(3),JDP(3),KDP(3),REF_F(3),VF(3),VP(3)
	REAL  R20(3),R80(3)
        REAL  RX,RY,RZ
        real yr1(512), yr2(512)
        real yr3(512), yd(512)

C
C**- PROJECT.F  PROJECT THE SAR RANGE RECORDS INTO DESIRED PROC. PLANE.
C
C	CONCEPT TEST PROGRAM.  IMPLEMENTS DENNIS'S FORMULATION FOR THE
C	COORDINATE TRANSFORMATION.
C
C	PHI(I)		ANGLE THE POLAR RECORD MAKES WITH THE PATCH
C			CENTERLINE, 1 DEGREES.
C	ANGLE		ANGULAR EXTENT OF THE PATCH, 1 DEGREES.
C	PROSCAL(I)	MAGNITUDE OF THE PROJECTION SCALE FACTOR FOR EACH
C			RANGE RECORD.  IF SLANT PLANE OF COLLECTED DATA
C			= PROJECTED PLANE, THEN PROSCAL(I)=1.
C	500		NO. POLAR RANGE RECORDS TO PROCESS.
C	VPMAG		PROJECTION OF CENTER RANGE RECORD INTO PROC. PLANE.
C	ISTEER		BEAM STEERING PARAMETER.  "0" MEANS NO STEERING.
C			"1" MEANS STEERING.
C	FS		RADAR A/D SAMPLE RATE (HZ) ALONG THE RANGE SWATH.
C	F0		RADAR CENTER FREQUENCY (HZ).
C	FDOT		RADAR FM RATE (HZ/SEC).
C	(RX,RY,RZ)	VECTOR COMPONENTS POINTING FROM THE CURRENT
C			PATCH CENTER TO THE NEW PATCH CENTER (1 METERS).
C
C
C        fileptr = yuppy_open("targs_4_vec.dat", "r+")
        
C        open(9,file='targs_4_vec.dat',iostat=ios, err = 997,
C     *            status = 'unknown')

        kdp(1) = 0
        kdp(2) = 0
        kdp(3) = 1

	CALL VNORM(KDP)
	ICOUNT=0
	MID=500/2	!DEFINE THE "MIDDLE" RANGE REC. AS CENTERLINE
	M20=NINT(.2*500)
	M80=NINT(.8*500)

	DO I=1,500
C      icnts = yuppy_fscanf(fileptr,"%d %f %f %f %f",
C     $     J,R(1),R(2),R(3),DELTAR(I))

C      READ(9,*,END=999)J,R(1),R(2),R(3),DELTAR(I)
       J=I
       R(1) = yr1(i)
       R(2) = yr2(i)
       R(3) = yr3(i)
       DELTAR(i) = yd(i)

C      if (i .le. 10)  then
C        print *, J,R(1),R(2),R(3),DELTAR(I)
C      endif
C           R(1) = R1(i)
C           R(2) = 25980.8
C           R(3) = 15000.0
C           DELTAR(i) = 0.0

C	   READ(16,*,END=999)J,R(1),R(2),R(3),DELTAR(I)
C
C	   SUBTRACT THE BEAM STEERING COMPONENT.
C
	   R(1)=R(1)-RX
	   R(2)=R(2)-RY
	   R(3)=R(3)-RZ
	   ICOUNT=ICOUNT+1
	   IF(I.EQ.MID)  THEN
	      R0X=R(1)
	      R0Y=R(2)
	      R0Z=R(3)
	   ELSEIF(I.EQ.M20) THEN
	      R20(1)=R(1)
	      R20(2)=R(2)
	      R20(3)=R(3)
	   ELSEIF(I.EQ.M80) THEN
	      R80(1)=R(1)
	      R80(2)=R(2)
	      R80(3)=R(3)
	   END IF
	END DO

	NMAX=ICOUNT

        kprime(1) = 0
        kprime(2) = 0
        kprime(3) = 0


	IF(KPRIME(1).EQ.0..AND.KPRIME(2).EQ.0..AND.KPRIME(3).EQ.0.)THEN
C**	COMPUTE NORMAL TO 20% AND 80% POINTING VECTORS
C**	THIS NORMAL WILL DEFINE THE PROCESSING PLANE (DEFAULT CONDITION).
	CALL CROSS(R20,R80,KPRIME)
	CALL VNORM(KPRIME)
C**	MAKE SURE KPRIME(3)>0.
		IF(KPRIME(3).LT.0.)THEN
		   KPRIME(1)=-KPRIME(1)
		   KPRIME(2)=-KPRIME(2)
		   KPRIME(3)=-KPRIME(3)
		END IF
	ELSE
		CALL VNORM(KPRIME)
	END IF

	CALL DOTPROD(KPRIME,KDP,B)

C	DEPRESS=ACOSD(B)
C       ACOSD is not available in FX

        DEPRESS=30.0

	S=SQRT(R0X**2+R0Y**2+R0Z**2)
	R0X=R0X/S
	R0Y=R0Y/S
	R0Z=R0Z/S
C
C	COMPUTE REFERENCE AXES 1 FOCUS PLANE.
C
	REF_F(1) = R0X
	REF_F(2) = R0Y
	REF_F(3) = R0Z
	CALL PROJ_F(REF_F,KDP,JDP)
	CALL VNORM(JDP)
	CALL CROSS(JDP,KDP,IDP)
C
C	COMPUTE REFERENCE AXES 1 PROCESSING PLANE.
C
	CALL INTERSECT(JDP,KDP,KPRIME,JPRIME)
	CALL VNORM(JPRIME)
	CALL CROSS(JPRIME,KPRIME,IPRIME)
C
C	WRITE 2 ANSWERS.
C
C	WRITE(6,*)
C	WRITE(6,*)'PROCESSING PLANE AXES:  '
C	WRITE(6,*)'i'' = ',IPRIME(1),'i + ',IPRIME(2),'j + ',IPRIME(3),'k'
C	WRITE(6,*)'j'' = ',JPRIME(1),'i + ',JPRIME(2),'j + ',JPRIME(3),'k'
C	WRITE(6,*)'k'' = ',KPRIME(1),'i + ',KPRIME(2),'j + ',KPRIME(3),'k'
C
C	WRITE(6,*)
C	WRITE(6,*)'FOCUS PLANE AXES:  '
C	WRITE(6,*)'i = ',IDP(1),'i + ',IDP(2),'j + ',IDP(3),'k'
C	WRITE(6,*)'j = ',JDP(1),'i + ',JDP(2),'j + ',JDP(3),'k'
C	WRITE(6,*)'k = ',KDP(1),'i + ',KDP(2),'j + ',KDP(3),'k'
C
C	NOW READ DATA FILE CONTAINING THE RANGE SWATH POINTING VECTORS
C
C	WRITE(6,*)
C	WRITE(6,*)'PROJECTING VECTORS INTO PROCESSING PLANE AND COMPUTING'
C	WRITE(6,*)'PHI AND |V|.'
	PI=4.*ATAN(1.)
	RD=180./PI
	MID=NMAX/2		!DEFINE MIDDLE RECORD
C	WRITE(6,*)'PROJECTED RESULTS:'
	INTERVAL=NMAX/40
	DO I=1,NMAX
C	   READ(16,*)J,R(1),R(2),R(3),DUMMY
C	   icnt = yuppy_fscanf("%d %f %f %f %f",J,R(1),R(2),R(3),DUMMY)
           J = I
           R(1) = yr1(i)
           R(2) = yr2(i)
           R(3) = yr3(i)
           DUMMY = 0.0

	   R(1)=R(1)-RX
	   R(2)=R(2)-RY
	   R(3)=R(3)-RZ
C	   UNIT NORMAL POINTING VECTOR FOR PULSE I.
	   CALL VNORM(R)
C	   PROJECT FROM ACQUISITION PLANE TO FOCUS PLANE.
	   CALL PROJ_F(R,KDP,VF)
C	   SEE WHERE PROJECTION INTERSECTS PROCESSING PLANE.
	   CALL INTERSECT(VF,KDP,KPRIME,VP)
C	   GET THE MAGNITUDE OF THIS VECTOR.
	   CALL VECMAG(VP,S)
	   PROSCAL(I) = 1 / S
C	   DOT VECTOR INTO I AND J AXES, THEN COMPUTE ANGLE FROM ATAN(*,*).
	   CALL DOTPROD(VP,JPRIME,VPDJP)
	   CALL DOTPROD(VP,IPRIME,VPDIP)
	   PHI(I) = ATAN2(VPDIP,VPDJP) * RD
	   IF(I.EQ.MID)  THEN
	      VPMAG=S			!MAG. OF PROJ. CENTER RECORD
	      VPTHETA=PHI(I)		!ANG.OF PROJECTED PATCH CENTERLINE
	      VPINDX=I			!RECORD NO. OF PROJECTED CENTERLINE
	   END IF
CCCC	   IF(MOD(I,INTERVAL).EQ.1) WRITE(6,*)I,PROSCAL(I),PHI(I)
	END DO
C
C	WE WON'T CLOSE THE "POINTING VECTOR" FILE YET.  IF BEAM STEERING IS 
C	REQUIRED, THE FILE WILL HAVE TO BE READ AGAIN JUST BEFORE
C	RANGE INTERPOLATION.  SIMPLY REWIND THE FILE.
C
C	REWIND(UNIT=16)
C
	ANGLE=ABS(PHI(1)-PHI(NMAX))	!ANGULAR EXTENT OF PATCH
C	WRITE(6,*)
C	WRITE(6,*)'ANGULAR EXTENT (IN DEG.) OF POLAR DATA = ',ANGLE
CC	WRITE(7,*)'ANGULAR EXTENT (IN DEG.) OF POLAR DATA = ',ANGLE
C	WRITE(6,*)'CENTER LINE (IN DEG.) = ',VPTHETA
C	WRITE(6,*)'RELATIVE START OF RECT. GRID = ',VPMAG
C	WRITE(7,*)'CENTER LINE (IN DEG.) = ',VPTHETA
C

       go to 1000
C997    continue
C       print *,'Open file failed, ios = ', ios
C       go to 1000
C999    continue
C       print *, 'end of file : ios = ', ios
1000   continue
C       close(9)
       RETURN

C	STOP
C
	END

C
C
	SUBROUTINE PROJ_F(VEC_IN,UNORM_F,VEC_OUT)
	REAL VEC_IN(3),UNORM_F(3),VEC_OUT(3),B(3)
C
C	GIVEN THE INPUT VECTOR VEC_IN, VEC_OUT IS THE VECTOR PROJECTED
C	INTO THE PLANE DEFINED BY THE UNIT NORMAL, UNORM_F.
C
C	THIS SUBROUTINE IS USED TO PROJECT A VECTOR INTO THE FOCUS PLANE.
C
	CALL DOTPROD(VEC_IN,UNORM_F,A)
	CALL VSCMULT(UNORM_F,A,B)
	CALL VSUB(VEC_IN,B,VEC_OUT)
	RETURN
	END
C
C
	SUBROUTINE INTERSECT(VEC_IN,UNORM_IN,UNORM_OUT,VEC_OUT)
	REAL VEC_IN(3),UNORM_IN(3),UNORM_OUT(3),VEC_OUT(3),TEMP(3)
C
C	THIS ROUTINE TAKES THE VECTOR VEC_IN THAT LIES 1 THE PLANE
C	DEFINED BY UNORM_IN AND PROJECTS IT 1 THE DIRECTION OF UNORM_IN
C	UNTIL IT INTERSECTS THE PLANE DEFINED BY UNORM_OUT.  VEC_OUT
C	IS THE VECTOR THAT NOW LIES 1 THE PLANE DEFINED BY UNORM_OUT.
C
	CALL DOTPROD(VEC_IN,UNORM_OUT,A)
	CALL DOTPROD(UNORM_IN,UNORM_OUT,B)
	IF(B.EQ.0.)THEN
C		WRITE(6,*)'*** CHOSEN VECTORS ARE ORTHOGONAL ***'
C		WRITE(6,*)'INTERSECTION CALCULATIONS IMPOSSIBLE...'
		STOP
	END IF
	X=A/B
	CALL VSCMULT(UNORM_IN,X,TEMP)
	CALL VSUB(VEC_IN,TEMP,VEC_OUT)
	RETURN
	END
C
C
C****************************************************************************
C	VECTOR OPERATIONS SUBROUTINES , I.E. DOT PRODUCT, CROSS PRODUCT, ETC.
C****************************************************************************
C
C
	SUBROUTINE DOTPROD(A,B,C)
C	COMPUTES C=A.B
	DIMENSION A(*),B(*)
	C=0.
	DO I=1,3
	C=C+A(I)*B(I)
	END DO
	RETURN
	END
C
	SUBROUTINE VSCMULT(A,S,C)
C	COMPUTES C=S*A
	DIMENSION A(*),C(*)
	DO I=1,3
	C(I)=S*A(I)
	END DO
	RETURN
	END
C
	SUBROUTINE VADD(A,B,C)
C	COMPUTES C=A+B
	DIMENSION A(*),B(*),C(*)
	DO I=1,3
	C(I)=A(I)+B(I)
	END DO
	RETURN
	END
C
	SUBROUTINE VSUB(A,B,C)
C	COMPUTES C=A-B
	DIMENSION A(*),B(*),C(*)
	DO I=1,3
	C(I)=A(I)-B(I)
	END DO
	RETURN
	END
C
	SUBROUTINE VSCADD(A,S,B)
C	COMPUTES B=A+S
	DIMENSION A(*),B(*)
	DO I=1,3
	B(I)=A(I)+S
	END DO
	RETURN
	END
C
	SUBROUTINE CROSS(A,B,C)
C	COMPUTES C=A X B
	DIMENSION A(*),B(*),C(*)
	C(1)=A(2)*B(3)-A(3)*B(2)
	C(2)=A(3)*B(1)-A(1)*B(3)
	C(3)=A(1)*B(2)-A(2)*B(1)
	RETURN
	END
C
	SUBROUTINE VECMAG(A,AMAG)
C	COMPUTES SQRT(A.A)
	DIMENSION A(*)
	CALL DOTPROD(A,A,C)
	AMAG=SQRT(C)
	RETURN
	END
C
	SUBROUTINE VNORM(A)
	DIMENSION A(*)
	CALL VECMAG(A,S)
	DO I=1,3
	A(I)=A(I)/S
	END DO
	RETURN
	END
C***********************************************************
C
	SUBROUTINE HEIFER(A,N,RX,RY,RZ,FS,F0,FDOT)
	REAL*8 PI,TP
	COMPLEX A(*)
C
C	PROGRAMMER: DENNIS C. GHIGLIA
C	DATE:       NOV. 2,1990
C	MODIFIED:   JUN. 1991 TO FIX SAMPLE TIMING ALONG THE PULSE TO
C		    BE SYMMETRIC ABOUT 2R0/C.  SEE D.C.G. NOTES FOR
C		    EXPLANATION.
C
C	**** NOTE:  THE UNITS USED 1 THIS PROGRAM MUST BE CONSISTENT
C		    WITH THE UNITS ASSOCIATED WITH THE PARAMETERS 1
C		    THE ARGUMENT LIST.  FOR EXAMPLE, IF THE UNITS
C		    1 THE ARGUMENT LIST ARE METERS, HZ, ETC. THEN
C		    THE SPEED OF LIGHT USED 1 THIS SUBROUTINE MUST BE
C		    1 METERS PER SECOND.  ****
C
C	THIS SUBROUTINE PERFORMS A COMPLEX MULTIPLICATION OF
C	THE COMPLEX VECTOR "A" BY EXP(-j*phase(n)) TO INVOKE STEERING
C	OF THE SAR PATCH CENTER TO ANOTHER POINT DEFINED BY THE
C	VECTOR COMPONENTS (RX,RY,RZ) POINTING FROM THE CURRENT
C	PATCH CENTER TO THE NEW DESIRED PATCH CENTER.
C
C	"A" IS THE COMPLEX 1D ARRAY CONTAINING A RANGE SWATH OF
C		PHASE HISTORY DATA ALREADY PHASE COMPENSATED TO PATCH
C		CENTER.
C
C	"R" IS A 3-COMPONENT ARRAY CONTAINING THE POINTING VECTOR
C		COMPONENTS FROM THE CURRENT PATCH CENTER TO THE RADAR
C		FOR THE RANGE SWATH OF INTEREST (METERS).
C
C	(RX,RY,RZ) ARE THE VECTOR COMPONENTS POINTING FROM THE CURRENT
C		PATCH CENTER TO THE NEW PATCH CENTER (METERS).
C
C	"N" IS THE NUMBER OF SAMPLES PER RANGE SWATH TO BE PROCESSED.
C
C	"FS" IS THE A/D SAMPLE RATE (HZ) ALONG THE RANGE SWATH.
C
C	"F0" IS THE RADAR CENTER FREQUENCY (HZ).
C
C	"FDOT" IS THE RADAR FM RATE (HZ/SEC).
C
C
C
	C=3.0E8	!SPEED OF LIGHT 1 M/SEC

C	PI=4.*DATAN(1.D0)
C       *** we seem to have problem with datan(1.)

        pi = 22.0/7.0
        
	TP=2.*PI
C
C	READ NEXT POINTING VECTOR FROM "POINTING VECTOR" FILE.
C
C
	RETURN
C
C	HANDLE FILE READING ERRORS.
C
C
	END
C
C***********************************************************************
C
C
	SUBROUTINE MAKFILT(FLT_R,LEN_R,NQ,BRF,DELTA,XMAX,TBLSIZ)
	REAL  FLT_R(NQ,LEN_R)
	INTEGER  TBLSIZ
C
C	FILENAME:    MAKFILT.FOR
C	LANGUAGE:    FORTRAN 77
C	PROGRAMMER:  DENNIS C. GHIGLIA AND GARY A. MASTIN
C	COMPANY:     SANDIA NATIONAL LABORATORIES
C	             ALBUQUERQUE, NM
C	DATE:        FEBRUARY 1990
C
C	"FLT_R" WILL CONTAIN A FINITE NUMBER OF RANGE INTERPOLATION
C	FILTERS.  THIS FINITE NUMBER IS THE NUMBER OF QUANTIZATION
C	LOCATIONS BETWEEN ANY TWO INPUT SAMPLES WHERE AN INTERPOLATED
C	OUTPUT SAMPLE MAY LIE.
C
C	
C	DETERMINE THE BANDWIDTH REDUCTION FACTOR FOR A FILTER OF 
C	A GIVEN SAMPLE LENGTH AND A SPECIFIED SAMPLE REDUCTION FACTOR.
C
C
	TBLSIZ=LEN_R*NQ
C
	PI=4.*ATAN(1.)
C
C	COMPUTE THE FILTER.
C
C	**** SINC INTERPOLATOR ****
C
	ARG1=PI/BRF
	ARG4=PI/XMAX
	DO  N=1,NQ
C
	DQ=REAL(N-1)/NQ
	NS=INT(XMAX-DQ)
	NF=INT(XMAX+DQ)
C*********************************
C	ADDED ERROR HANDLER ON 2/11/91.
C
	IF(NS+NF+1.GT.LEN_R)  THEN
	   WRITE(6,*)'FILTER LENGTH EXCEEDS LEN_R = ',LEN_R,' IN MAKFILT.'
	   WRITE(6,*)'*** THIS IS A FATAL ERROR. ***'
	   STOP
	END IF
C********************************
C
	   DO  M=1,NS+NF+1
C
	      ARG2=(M-1-NS-DQ)/DELTA
	      ARG3=ARG1*ARG2
	      IF(ARG3.EQ.0.)  THEN
	         FLT_R(N,M)=1.0
	      ELSE
	         FLT_R(N,M)=(SIN(ARG3)/(ARG3))*(.5+.5*COS(ARG4*ARG2))
	      END IF
	   END DO
	END DO

C
C	**** TWO POINT LINEAR INTERPOLATOR
C
C	DO  N=1,NQ

C	DQ=REAL(N-1)/NQ
C	NS=INT(1.-DQ)	!SETS XMAX=1. BY DEFAULT
C	NF=INT(1.+DQ)	!SETS XMAX=1. BY DEFAULT

C	   DO  M=1,NS+NF+1
C	      ARG=ABS(M-1-NS-DQ)
C	      FLT_R(N,M)=1.-ARG
C	   END DO
C	END DO
C
C	**** UNIFORM INTERPOLATOR ****
C
C	DO  N=1,NQ
C
C	DQ=REAL(N-1)/NQ
C	NS=INT(XMAX-DQ)
C	NF=INT(XMAX+DQ)
C
C	   DO  M=1,NS+NF+1
C	         FLT_R(N,M)=1.0
C	   END DO
C	END DO

	RETURN
	END
C
C
C
	SUBROUTINE LOCATE_R(ANG,DELTA,LEN,R1,DELTAR,ISO,LENSIG,IHO,
     *  NQ,VC,PROSCAL)
C	INTEGER ISO(LENSIG),IHO(LEN)
	INTEGER ISO(LEN),IHO(LEN)	!JAN. 1991
C
C**- LOCATE_R.FOR  LOCATE DATA POINTS FOR EACH RANGE INTERPOLATION FILTERING.
C
C	FILENAME:    LOCATE_R.FOR
C	LANGUAGE:    FORTRAN 77
C	PROGRAMMER:  DAVID K. JOHNSON & GARY A. MASTIN
C	COMPANY:     SANDIA NATIONAL LABORAOTORIES
C	             ALBUQUERQUE, NM
C	DATE:        JANUARY 1990
C	             MODIFIED CALLING ARGUMENTS, 2-OF-PLANE PROJECTION.
C	             DOUBLE PRECISION "RNG" ON FEB. 21, 1991
C
C
C	1) LOCATE THE EXACT POINT BETWEEN ANY TWO INPUT POINTS WHERE
C	THE DESIRED INTERPOLATION POINT SHOULD BE, 2) QUANTIZE THAT
C	LOCATION TO THE NEAREST OF "NQ" QUANTIZATION LOCATIONS,
C	3) FIND THE CENTRAL INPUT VALUE OF THE FINITE-POINT INTERPOLATION
C	THAT IS TO BE DONE, AND 4) SELECT THE FILTER NUMBER FOR THAT
C	INTERPOLATION.
C
C	CONVERT "PHI" 1 DEGREES TO RADIANS FOR CALCULATING THE COSING
C	OF "PHI" TO DETERMINE THE RECTANGULAR POSITION.
C
	PI=4.*ATAN(1.0)			!Machine precision of PI
	PHI=ANG*PI/180.			!Convert PHI in degrees to radians
C
C	BEGIN LOOP TO DETERMINE THE QUANTIZATION POINT (ISO) AND THE
C	FILTER TO BE USED (IHO).
C
C	"DELTA" IS THE SPACING 1 RANGE OF THE RECTANGULAR GRID.
C

	DO I=1,LEN
C
C	   INCREMENT 1 RANGE ON RECTANGULAR GRID.
C
C	   RNG=R1+(I-1)*DELTA
C	   ******* CHANGED 2/21/91 *****
	   RNG=DBLE(R1)*DBLE(VC)+DBLE(I-1)*DBLE(DELTA)
C	   *****************************
C
C	   COMPUTE WHERE HORIZONTAL LINE OF GRID CROSSES POLAR RECORD.
C
C	   R=ABS(RNG/COS(PHI))
C
C	   SUBTRACT THE RANGE TO THE NEAR-EDGE OF THE PATCH.
C
C	   R=R*VC/VMAG-R0
C	   ******* CHANGED 2/21/91 *******
	   R=DABS(RNG/DCOS(DBLE(PHI)))*PROSCAL-DBLE(R1)-DBLE(DELTAR)
C	   *******************************
C
C
C	   QUANTIZE LOCATION?
C
	   REM=AMOD(R,1.0)
C
C	   NINT(REM*NQ) IS AN INTEGER EQUAL TO THE NEAREST QUANTIZATION PT.
C	   IHO(I) POINTS TO THE FILTER TO BE USED, HAVING VALUES BETWEEN
C	   0 AND NQ-1.
C
	   IVAL=NINT(REM*NQ)
C
C	   ISO(I) IS THE INTEGER VALUE OF "R" AFTER QUANTIZATION.
C
	   ISO(I)=INT(R)+IVAL/NQ
	   IHO(I)=MOD(IVAL,NQ)
	END DO
C
	RETURN
	END
C
C
	SUBROUTINE LOCATE_CR(CR,NUMSMP,PHI,THETA,DELTAX,DELTAY,RNG,
     *  XORIG,ISO,RECT_C)
C	PARAMETER  (IOFF=64)
	INTEGER  RECT_C
	INTEGER  ISO(RECT_C)
	REAL  CR(NUMSMP),PHI(NUMSMP)

C
C	FILENAME:    LOCATE_CR.FOR
C	LANGUAGE:    FORTRAN 77
C	PROGRAMMER:  GARY A. MASTIN
C	COMPANY:     SANDIA NATIONAL LABORATORIES
C	             ALBUQUERQUE, NM
C	DATE:        FEB. 1990
C
C	LOCATE THE CROSS-RANGE SAMPLE POINTS FOR ONE HORIZONTAL LINE
C	OF THE OUTPUT GRID.  DETERMINE THE FIRST SAMPLE POINT TO USE
C	1 INTERPOLATION FOR EACH OUTPUT POINT.
C
C	NUMSMP		NUMBER OF CROSS-RANGE SAMPLES AFTER RANGE CORRECTION.
C	CR(*)		ARRAY CONTAINING LOCATION 1 CROSS-RANGE OF SAMPLES
C			AFTER RANGE CORRECTION.
C	IOFF		WHEN SEARCHING FOR THE ELEMENT TO PUT "CR(*)"
C			WHICH IS CLOSEST TO THE DESIRED OUTPUT POINT, WE
C			WISH TO LIMIT THE SEARCH AREA 1 ORDER TO SAVE TIME
C			THIS IS DONE BY LIMITING THE SEARCH FROM THE POINT
C			CLOSEST TO THE LAST DESIRED PT. TO A POINT "IOFF"
C			SAMPLES FROM THE LAST CLOSEST POINT.
C			NOTE:  IF "THETA" IS MORE THAN A FEW DEGREES, 
C			"IOFF" MAY NEED TO BE INCREASED.
C	PHI(K)		ANGLE (DEG.) BETWEEN Kth POLAR RECORD AND
C			THE CENTER LINE.  THIS IMPLIES THAT THE RECORDS
C			NEED NOT BE EQUALLY SEPARATED 1 ANGLE. 
C	THETA		ANGULAR EXTENT (DEG.) OF POLAR DATA.
C	DELTAX		X-DIMENSION (CROSS-RANGE) SAMPLE SPACING.
C	DELTAY		Y-DIMENSION (RANGE) SAMPLE SPACING.
C	RNG		NUMBER OF SAMPLES TO THE FULL GRID HORIZONTAL LINE
C			OF INTEREST.
C	XORIG		X-AXIS (CROSS-RANGE) ORIGIN OF THE FULL RECT. GRID.
C	ISO(K)		THE Kth ELEMENT OF ISO POINTS TO THE ELEMENT NUMBER
C			OF "CR(*)" WHICH IS THE CENTRAL SAMPLE PT.
C			TO USE 1 INTERPOLATION TO GET THE Kth OUTPUT PT.
C	RECT_C		NO. OF CROSS-RANGE ELEMENTS 1 FULL RECT. GRID.
C
C	INITIALIZATION.
C
	PI=4.*ATAN(1.)
C	IOFF=16*NINT(DELTAX)	!OLD STATEMENT
	IOFF=NINT(16*DELTAX)	!NEW STATEMENT 9-4-91
	DO  J=1,NUMSMP
	   ANG=PHI(J)*PI/180.
C
C	   HORIZ. DISTANCE FROM THE CENTERLINE TO THE DESIRED RANGE-CORRECTED
C	   POINT ON A SPECIFIED CROSS-RANGE RECORD.
C
	   DST=RNG*TAN(ANG)
C
C	   CROSS-RANGE POSITION OF SAMPLE WITH RESPECT TO "XORIG".
C
	   CR(J)=XORIG+DST
C***********************************************************************
C
C	   DEBUG.  UN-COMMENT TO WRITE AN OUTPUT FILE FOR PLOTTING THE
C	   RANGE-CORRECTED DATA ON A RECT. GRID.
C
C	   WRITE(10,*)CR(J),RNG	!EXISTING X, EXISTING Y
C
C***********************************************************************
	END DO
C
C	IF PHI(1) IS LESS THAN PHI(NUMSMP), THEN WE WILL SELECT
C	ISO VALUES FROM LEFT TO RIGHT.  HOWEVER, IF PHI(1)
C	IS GREATER THAN PHI(NUMSMP),  WE WILL SELECT ISO VALUES
C	FROM RIGHT TO LEFT.  THIS DISTINCTION IS VERY IMPORTANT
C	WHEN WE TRY TO LIMIT THE SEARCH AREA OF CR(J) TO SELECT
C	AN ISO VALUE.
C
	IF(PHI(1).LT.PHI(NUMSMP))  THEN
	   IDIR=1	!LEFT TO RIGHT
	ELSE
	   IDIR=0	!RIGHT TO LEFT
	END IF
C
C	NOW SELECT IOS(*).
C	IT POINTS TO THE CENTRAL INPUT SAMPLE TO USE 1 FILTERING.
C
	ISTRT=1		!ELEMENT NO. OF CR(*) WHICH IS THE START OF
			!OUR SEARCH FOR THE ELEMENT OF CR(*) CLOSEST
			!TO THE DESIRED OUTPUT POINT OF THE GRID.
	IFIN=NUMSMP	!WORST CASE FINAL ELEMENT OF CR(*) TO SEARCH.
	DO  I=1,RECT_C
	   XD=(I-1)*DELTAX	!DESIRED OUTPUT GRID LOCATION.
C***********************************************************************
C
C	   DEBUG.  UN-COMMENT TO WRITE AN OUTPUT FILE FOR PLOTTING THE
C	   FULLY CORRECTED DATA ON A RECTANGULAR GRID.
C
C	   WRITE(11,*)XD,RNG
C
C***********************************************************************
	   RMIN=NUMSMP*DELTAX	!AN ADEQUATELY LARGE VALUE.
	   DO  J=ISTRT,IFIN
	      DIST=ABS(XD-CR(J))
	      IF(DIST.LE.RMIN)  THEN
	         RMIN=DIST
	         JSAV=J
	      END IF
	   END DO
C
	   ISO(I)=JSAV		!CENTRAL SAMPLE POINT TO USE.
	   IF(IDIR.EQ.1)  THEN	!LEFT TO RIGHT
	      ISTRT=JSAV		!LIMIT THE SEARCH BOUNDS.
	   ELSE
	      ISTRT=JSAV-IOFF	      
	      IF(ISTRT.LT.1)  ISTRT=1
	   END IF
	   IFIN=ISTRT+IOFF	!SET LAST SEARCH PT. REASONABLY LARGE.
	   IF(IFIN.GT.NUMSMP)  IFIN=NUMSMP
	END DO
C
	RETURN
	END
C
C
	SUBROUTINE INTERP_R(INPUT,LENIN,OUTPUT,ISIZE,LENOUT,iTBLSIZ,
     *  FILTER,NQ,ISO,IHO,ANG,DELTA,IXMAX,BRF,PROSCAL)
	INTEGER ISO(LENIN),IHO(LENOUT),iTBLSIZ
C	ISIZE IS PROBABLY MUCH LARGER THAN LENOUT
	REAL    FILTER(iTBLSIZ)
	COMPLEX INPUT(LENIN),OUTPUT(ISIZE)
C
C
C	FILENAME:    INTERP_R.FOR
C	LANGUAGE:    FORTRAN 77
C	PROGRAMMER:  GARY A. MASTIN
C	COMPANY:     SANDIA NATIONAL LABORATORIES
C	             ALBUQUERQUE, NM
C	DATE:        JANUARY 1990
C	             MODIFIED CALLING ARGUMENTS 2/11/91.
C
C
C	RQ		RANGE (MINUS RANGE TO NEAR-EDGE) 1 SAMPLES TO THE
C			QUANTIZATION POINT APPROX. THE EXACT LOCATION OF THE
C			DESIRED INTERPOLATION POINT.
C	XMAX		ONE-SIDED LENGTH OF THE FILTER 1 TERMS OF INPUT
C			SAMPLES.
C	XL		LEFT-HAND LIMIT OF THE FILTER SUPPORT 1 TERMS OF
C			INPUT SAMPLES.
C	XR		RIGHT-HAND LIMIT OF THE FILTER SUPPORT 1 TERMS OF
C			INPUT SAMPLES.
C	ISTRT		FIRST INPUT SAMPLE TO BE USED 1 INTERPOLATION.
C	IFIN		LAST INPUT SAMPLE TO BE USED 1 INTERPOLATION.
C			NOTE:  THE FULL LENGTH OF THE INTERPOLATION FILTER
C			VARIES WITH "PHI", BUT IS ALWAYS DEFINED BY
C			ISTRT AND IFIN.
C	BRF		FILTER BANDWIDTH REDUCTION FACTOR ABOVE THE NYQUIST
C			LIMIT FOR THE GIVEN OUTPUT SAMPLE SPACING.  FOR
C			EXAMPLE, BRF=1.04 IMPLIES THAT THE NYQUIST BANDWIDTH
C			WILL BE REDUCED BY 1/1.04.
C	PROSCAL		PROJECTION SCALING FACTOR FOR 2-OF-PLANE 
C			CORRECTION
C
C	INITIALIZE THE OUTPUT RECORD.
C

	DO  I=1,LENOUT
	   OUTPUT(I)=(0.,0.)
	END DO
C
C	DO THE INTERPOLATION.
C
	PI=4.*ATAN(1.)
	PHI=ANG*PI/180.
        XMAX=IXMAX*ABS((DELTA*PROSCAL)/COS(PHI))*BRF
	DO  K=1,LENOUT	!THERE WILL BE "LENOUT" SAMPLE POINTS OUTPUT
        NS=INT(XMAX-REAL(IHO(K))/NQ)
        NF=INT(XMAX+REAL(IHO(K))/NQ)
        LENFLT=NS+NF+1

C        if ((fxcellid .eq. 0) .and. (k.eq. 1)) then
C          print *, 'ns = ', ns, ' nf= ', nf, ' lenflt = ', lenflt
C          print *, 'iso( ', k, ')= ', iso(k), ' lenin = ', lenin
C        endif

	   IF(ISO(K).GE.1.AND.ISO(K).LE.LENIN)  THEN
	      ISTRT=ISO(K)-NS
	      IFIN=ISO(K)+NF
C	      WRITE(7,*)
C	      WRITE(7,*)'OUTPUT SAMPLE POINT IS ',K
C	      WRITE(7,*)'ISO(',K,') = ',ISO(K)
C	      WRITE(7,*)'IHO(',K,') = ',IHO(K)
C	      WRITE(7,*)'XMAX = ',XMAX
C	      WRITE(7,*)'ISTRT = ',ISTRT,' IFIN = ',IFIN
C
C	      DON'T ALLOW USE OF POINTS OUTSIDE THE SAMPLED DOMAIN FOR INTERP.
C

	      IF(ISTRT.GE.1.AND.IFIN.LE.LENIN)  THEN
C
C	         THIS IS THE WEIGHTING (INTERPOLATION) FUNCTION FOR THE Kth PT.
C   
	         RN=0.
	         DO  J=1,LENFLT	!THERE ARE "LENFLT" WEIGHTINGS PER OUTPUT PT.
	            INPT=ISTRT+J-1
	            IFLT=IHO(K)+NQ*(J-1)+1		!FILTER SAMPLE NUMBER
C	            WRITE(7,*)'OUTPUT(',K,')=OUTPUT(',K,')+INPUT(',INPT,
C     *              ')*FILTER(',IFLT,')'
C	            WRITE(7,*)'OUTPUT(',K,')=',OUTPUT(K),'+',INPUT(INPT),
C     *              '*',FILTER(IFLT)
	            RN=RN+FILTER(IFLT)
	            OUTPUT(K)=OUTPUT(K)+INPUT(INPT)*FILTER(IFLT)
C	            WRITE(7,*)'OUTPUT(',K,') = ',OUTPUT(K)
	         END DO
	         OUTPUT(K)=OUTPUT(K)/RN
	      END IF
	   END IF
	END DO
C
	RETURN
	END
C
C
	SUBROUTINE INTERP_CR(CR,nopol,OUTPUT,iRECT_C,ISTRT,
     *  DELTAX,INPUT,BRF,XMAX)
	INTEGER  iRECT_C,ISTRT(iRECT_C)
	COMPLEX  INPUT(500),OUTPUT(512)
	REAL     CR(500)
        integer rect_c
C
C**- INTERP_CR.FOR  INTERPOLATION OF UNEQUALLY SPACED CROSS-RANGE SAMPLES 
C
C	FILENAME:    INTERP_CR
C	LANGUAGE:    FORTRAN 77
C	PROGRAMMER:  GARY A. MASTIN
C	COMPANY:     SANDIA NATIONAL LABORATORIES
C	             ALBUQUERQUE, NM
C	DATE:        FEB. 1990
C
C	IMPORTANT ARRAYS, VARIABLES, ETC.:
C	CR(*)		ARRAY CONTAINING LOCATIONS 1 CROSS-RANGE OF SAMPLES
C			AFTER RANGE CORRECTION.
C	500		NUMBER OF CROSS-RANGE SAMPLES AFTER RANGE CORRECTION.
C			(I.E. BEFORE CROSS-RANGE INTERPOLATION)
C	RECT_C		NUMBER OF CROSS-RANGE SAMPLES AFTER CROSS-RANGE INTERP.
C			IS SAME AS NUMBER OF VERTICAL LINES 1 FULL RECT. GRID.
C	DELTAX		X-DIMENSION (CROSS-RANGE) SAMPLE SPACING.
C	ISTRT(K)	THE Kth ELEMENT OF ISTRT POINTS TO THE ELEMENT NUMBER
C			OF "CR(*)" WHICH IS THE CENTRAL POINT
C			TO USE 1 INTERPOLATION TO GET THE Kth OUTPUT PT.
C	RECT_C		NO. OF CROSS-RANGE ELEMENTS 1 FULL RECT. GRID.
C	OUTPUT(*)	CONTAINS THE SAR SAMPLES EQUALLY SPACED ON A LINE
C			OF THE FULL RECT. GRID.
C	INPUT(*)	CONTAINS THE UNEQUALLY SPACED COMPLEX SAR AMPLITUDES
C			ON A LINE OF THE FULL RECT. GRID.
C	BRF		BANDWIDTH REDUCTION FACTOR.
C
C	INITIALIZE.
C
        rect_c = irect_c

	DO  I=1,RECT_C
	   OUTPUT(I)=(0.,0.)
	END DO
C
C	INTERPOLATE.
C
	PI=4.*ATAN(1.)
	ARG1=PI/BRF
	ARG4=PI/XMAX

C        if (fxcellid .eq. 0) then
C          do i=121,130
C            print *, 'istrt(',i,')= ', istrt(i)
C          enddo
C        endif

	DO  I=1,RECT_C
	   XD=(I-1)*DELTAX
	   RNORM=0.

	   IF((XD+XMAX.GT.AMAX1(CR(1),CR(500))).OR.(XD-XMAX.LT.
     *     AMIN1(CR(1),CR(500))))  THEN
	      OUTPUT(I)=(0.,0.)
	   ELSE
C
C	      MOVE CENTER TO RIGHT.
C

	      DO  J=ISTRT(I),500
	          IF(ABS(CR(J)-XD).LE.XMAX)  THEN
	            ARG2=(CR(J)-XD)/DELTAX
	            ARG3=ARG1*ARG2
	            IF(ARG3.EQ.0.0)  THEN
	               FILTER=1.
	            ELSE
	               FILTER=SIN(ARG3)/(ARG3)*(.5+.5*COS(ARG4*ARG2))
	            END IF
	            RNORM=RNORM+FILTER

C                  print *,"i=",i," j=",j,"  rnorm=",rnorm
C                  print *,"output(i)=",output(i)," input(j)=",input(j)

	            OUTPUT(I)=OUTPUT(I)+FILTER*INPUT(J)
	         ELSE
C	            PREVENT DEAD CYCLING.
	            GO TO 200
	         END IF
	      END DO
200	      CONTINUE

C
C	      MOVE CENTER TO LEFT.
C
	      DO  J=ISTRT(I)-1,1,-1
	         IF(ABS(CR(J)-XD).LE.XMAX)  THEN
	            ARG2=(CR(J)-XD)/DELTAX
	            ARG3=ARG1*ARG2
	            IF(ARG3.EQ.0.0)  THEN
	               FILTER=1.
	            ELSE
	               FILTER=SIN(ARG3)/(ARG3)*(.5+.5*COS(ARG4*ARG2))
	            END IF
	            RNORM=RNORM+FILTER
	            OUTPUT(I)=OUTPUT(I)+FILTER*INPUT(J)
	         ELSE
	            GO TO 400
	         END IF
	      END DO
400	      CONTINUE

	      IF(RNORM.EQ.0.0)  THEN
	         OUTPUT(I)=(0.,0.)
	      ELSE
	         OUTPUT(I)=OUTPUT(I)/RNORM
	      END IF

	   END IF
	END DO
C
	RETURN
	END
C

        subroutine fft2d(block,tblock,pt,tt)

	COMPLEX  BLOCK(512,512),WINDOW(512),TBLOCK(512,512)

	COMPLEX  dummy(512)
        real pt,tt

        integer i
        

      pt = 0.0
      tt = 0.0

      call timer_start()

      call init_all(numrow,numcol,MINP,M,NINP,N,
     *                      IDIR,MWSIZE,NWSIZE,IOPT,S1,NBAR)


C
C
C*********************************************************************
C
C	PERFORM 1-D TRANSFORM, WRITE RESULTS 2 TO BLOCK TRANSPOSE
C	TEMPORARY FILE.
C
C*********************************************************************
C
C	CALL CFFT1D(M,BLOCK,0)	!INITIALIZE FOR FFT.

           do i=1,512
             window(i) = (0.0,0.0)
           enddo
	   CALL GEN_TAYLOR(S1,NBAR,MWSIZE,WINDOW)

C** BEGIN CONCURRENCY LOOP

	   DO  II=1,NINP
C           pdo ii=1,512
C           pin block(:,II)

              I = II
              do ix = 1, 512
                dummy(ix) = BLOCK(ix,II)
              enddo


C	      TAYLOR WEIGHT.
C
C	      DO  KK=1,MINP
C	         BLOCK(kk,II)=BLOCK(KK,II)*WINDOW(KK)
C	      END DO
C
	      DO  KK=1,MINP
	         dummy(kk)=dummy(kk)*WINDOW(KK)
	      END DO


C	      PERFORM THE 1-D FAST FOURIER TRANSFORM.
C

C	      IF(IDIR.EQ.-1)  CALL MUXIT(BLOCK(1,II),M,N,I)	!FORWARD
C	      CALL CFFT1D(M,BLOCK(1,II),IDIR)

	      IF(IDIR.EQ.-1)  CALL MUXIT(dummy,M,N,I)	!FORWARD
	      CALL CFFT1D(M,dummy,IDIR)

              do ix = 1, 512
                 BLOCK(ix,II) = dummy(ix)
              enddo
             
	   END DO
C           endpdo

        pt = timer_stop()

C** END CONCURRENCY LOOP

      call timer_start()

      call idxperm(tblock,block)

      tt = timer_stop()

 
C	   CALL TRSP(BLOCK,TBLOCK)

C*************************************************************************
C
C	READ DATA BACK INTO BLOCK 1 TRANSPOSED FORMAT.  PERFORM
C	1-D FFT. WRITE 2 FINAL DATA.
C
C*************************************************************************
C
C	GET READY TO LOAD A MAIN MEMORY BLOCK.  SET "BLKELM" SO THE
C	DATA ARE READY.
C
C	CALL CFFT1D(N,tBLOCK,0)	!INITIALIZE FOR FFT.

           call timer_start()

           do i=1,512
             window(i) = (0.0,0.0)
           enddo

	   CALL GEN_TAYLOR(S1,NBAR,NWSIZE,WINDOW)

	   DO  II=1,512

C           pdo II=1,512
C           pin tBLOCK(:,II)
     
              I = II

              do ix =1, 512
                dummy(ix) = tBLOCK(ix,II)
              enddo

C	      DO  KK=1,NINP
C	         tBLOCK(KK,II)=tBLOCK(KK,II)*WINDOW(KK)
C	      END DO

	      DO  KK=1,NINP
	         dummy(KK)=dummy(KK)*WINDOW(KK)
	      END DO

C	      CALL CFFT1D(N,tBLOCK(1,II),IDIR)
C	      IF(IDIR.EQ.1)  CALL MUXIT(tBLOCK(1,II),N,M,I)	!INVERSE


	      CALL CFFT1D(N,dummy,IDIR)
	      IF(IDIR.EQ.1)  CALL MUXIT(dummy,N,M,I)	!INVERSE

              do ix =1, 512
                tBLOCK(ix,II) = dummy(ix)
              enddo
	   END DO
C	   endpdo

           pt = pt + timer_stop()

C** END CONCURRENCY LOOP

        return
	END


C
C
        subroutine init_all(numrow,numcol,MINP,M,NINP,N,
     *                      IDIR,MWSIZE,NWSIZE,IOPT,S1,NBAR)
        integer numrow,numcol,MINP,M,NINP,N
        integer IDIR,MWSIZE,NWSIZE,IOPT,NBAR
        real S1

        numrow = 397
        numcol = 500
        MINP = 500
        M=512
        NINP = 397
        N=512

        IDIR=-1
        MWSIZE = 500

        NWSIZE = 397
        IOPT = 4
        S1= -35.0
        NBAR= 5

        return
        end
C
C
C
	SUBROUTINE CFFT1D(N,X,IDIR)
	COMPLEX  X(*)
C
C	MAKE 1-D FFT CALL LOOK LIKE OPTIMIZED ONE ON ALLIANT.
C
	IF(IDIR.EQ.0)  THEN
C	   INITIALIZE.
	   RETURN
	ELSE IF(IDIR.EQ.-1)  THEN
C	   FORWARD FFT.
	   INV=0
	ELSE IF(IDIR.EQ.1)  THEN
C	   INVERSE
	   INV=1
	END IF
	CALL FFT(X,N,INV)
	RETURN
	END
C
C-----------------------------------------------------------------------
C SUBROUTINE:  FFT
C JIM COOLEY'S SIMPLE FFT PROGRAM--USES DECIMATION IN TIME ALGORITHM
C X IS AN N=2**M POINT COMPLEX ARRAY THAT INITIALLY CONTAINS THE INPUT
C AND ON OUTPUT CONTAINS THE TRANSFORM
C THE PARAMETER INV SPECIFIED DIRECT TRANSFORM IF 0 AND INVERSE IF 1
C-----------------------------------------------------------------------
C
      SUBROUTINE FFT(X, N, INV)
      COMPLEX X(1), U, W, T, CMPLX
C
C   X = COMPLEX ARRAY OF SIZE N--ON INPUT X CONTAINS
C       THE SEQUENCE TO BE TRANSFORMED
C       ON OUTPUT X CONTAINS THE DFT OF THE INPUT
C   N = SIZE OF FFT TO BE COMPUTED--N=2**M FOR 1.LE.M.LE.15
C INV = PARAMETER TO DETERMINE WHETHER TO DO A DIRECT TRANSFORM (INV=0)
C       OR AN INVERSE TRANSFORM (INV=1)
C
      M = ALOG(FLOAT(N))/ALOG(2.) + .1
      NV2 = N/2
      NM1 = N - 1
      J = 1
      DO 40 I=1,NM1
        IF (I.GE.J) GO TO 10
        T = X(J)
        X(J) = X(I)
        X(I) = T
  10    K = NV2
  20    IF (K.GE.J) GO TO 30
        J = J - K
        K = K/2
        GO TO 20
  30    J = J + K
  40  CONTINUE
      PI = 4.*ATAN(1.0)
      DO 70 L=1,M
        LE = 2**L
        LE1 = LE/2
        U = (1.0,0.0)
        W = CMPLX(COS(PI/FLOAT(LE1)),-SIN(PI/FLOAT(LE1)))
        IF (INV.NE.0) W = CONJG(W)
        DO 60 J=1,LE1
          DO 50 I=J,N,LE
            IP = I + LE1
            T = X(IP)*U
            X(IP) = X(I) - T
            X(I) = X(I) + T
  50      CONTINUE
          U = U*W
  60    CONTINUE
  70  CONTINUE
      IF (INV.EQ.0) RETURN
      DO 80 I=1,N
C        X(I) = X(I)/CMPLX(FLOAT(N),0.)
        X(I) = X(I) * (1.0 / FLOAT(N))
  80  CONTINUE
      RETURN
      END
C
	SUBROUTINE MUXIT(A,M,N,LINNO)
	COMPLEX A(M)
C
C	SUBROUTINE TO MULTIPLEX THE COMPLEX ARRAY A(M,N)
C	WHERE M AND N ARE POWERS OF TWO.  THIS ROUTINE
C	IS USED TO SHIFT THE DC TERM TO SAMPLE M/2 + 1, N/2 + 1,
C	PRIOR TO A 2D FFT CALL.
C
	JSIGN=1-(MOD(LINNO,2)*2)	!ROW SIGN
	ISIGN=-1			!COLUMN SIGN
	DO I=1,M
	   ISIGN=-ISIGN
	   A(I)=ISIGN*JSIGN*A(I)
	END DO
	RETURN
	END
C
C
C
C
	SUBROUTINE GEN_TAYLOR(S1,NBAR,WSIZE,WINDOW)
	INTEGER  NBAR,WSIZE,P
	REAL  S1,WINDOW(*),FM(10),NUM(10),DEN(10)
	REAL  DEN_PROD,NUM_PROD,PROD1,PROD2
C+
C**- GEN_TAYLOR.F  GENERATE A TAYLOR WINDOW FOR SAR PHASE HISTORY.
C
C	FILENAME:    GEN_TAYLOR.F
C	LANGUAGE:    FORTRAN 77
C	PROGRAMMER:  DAN WAHL, ORG. 9117
C	             IMPLEMENTED HERE BY GARY MASTIN
C	COMPANY:     SANDIA NATIONAL LABORATORIES
C	             ALBUQUERQUE, NM
C	DATE:        APRIL 1991
C
C	IMPORTANT ARGUMENTS, VARIABLES, AND ARRAYS:
C	S1		NUMBER OF db OF SIDE LOBE SUPPRESSION (I.E. 35 db)
C	NBAR		SPECIFIED ALONG WITH SIDE LOBE LEVEL:
C	 		20 db	NBAR = 2
C			25 db 	NBAR = 3
C			30 db	NBAR = 4
C			35 db	NBAR = 5
C			40 db	NBAR = 6
C			45 db	NBAR = 8
C	NBAR	
C-
	PI=4.*ATAN(1.)
C
	B=10.0**(-S1/20.0)
	A=LOG(B + SQRT(B*B-1.0))/PI
C
	SIGP = FLOAT(NBAR)/SQRT(A*A + (FLOAT(NBAR) - 0.5)**2)
C
	DO 10  M=1,NBAR-1
	   DEN_PROD=1.
	   NUM_PROD=1.0
	   DO  20  N=1,NBAR-1
	      NUM(N)=(1.0-((FLOAT(M)**2)/((SIGP*SIGP)*
     *	             (A*A + (FLOAT(N) - 0.5)**2))))
	      NUM_PROD=NUM_PROD*NUM(N)
20	   CONTINUE
	   DO  30  P=1,NBAR-1
	      DEN(P)=(1.0 -(FLOAT(M)/FLOAT(P))**2)
	      IF(DEN(P).EQ.0.0)  DEN(P)=1.0
	      DEN_PROD=DEN_PROD*DEN(P)
30	   CONTINUE
	   FM(M)=0.5*((-1.0)**(M+1))*NUM_PROD/DEN_PROD
10	CONTINUE
C
C	NOW CALCULATE THE WINDOW.
C
	DO 100  I=1,WSIZE
	   T=(I-1.0-FLOAT(WSIZE)/2.0)/FLOAT(WSIZE)
c	   T=(FLOAT(I-1-FLOAT(WSIZE)/2.0))/FLOAT(WSIZE)
	   PROD1=0.0
	   PROD2=0.0
	   DO 110  J=1,NBAR-1
	      PROD1=PROD1+FM(J)*COS(2.0*PI*T*FLOAT(J))
	      PROD2=PROD2+FM(J)
110	   CONTINUE
	   WINDOW(I)=(1.0 + 2.0*PROD1)/(1.0 + 2.0*PROD2)
100	CONTINUE
C
	RETURN
	END



	SUBROUTINE DREADR( CHANNEL, BUFFER, RECORDSIZE, RECORDNUMBER )
	INTEGER CHANNEL, RECORDSIZE, RECORDNUMBER
	complex BUFFER(recordsize)
	READ(CHANNEL, rec=RECORDNUMBER) BUFFER
	RETURN
	END

	SUBROUTINE DWRITR( CHANNEL, BUFFER, RECORDSIZE, RECORDNUMBER )
	INTEGER CHANNEL, RECORDSIZE, RECORDNUMBER
	complex buffer(recordsize)
        WRITE(CHANNEL,rec=RECORDNUMBER) BUFFER
	RETURN
	END

        subroutine idxperm(out,in)
        complex out(512,512), in(512,512)
        integer i,j

        do i = 1,512
           do j=1, 512
             out(i,j) = in(j,i)
           enddo
        enddo

        return
        end

        subroutine timer_start()
        real tarray(2)

        call dtime(tarray)
        return
        end

        function timer_stop()
        real timer_stop
        real tarray(2)

        timer_stop = dtime(tarray)
        return
        end
