      PROGRAM CYLPRBC
C
C   The rationalized mks system of units is used throughout.
C
      PARAMETER (NPH = 360, IC2T = 1024)
C
C-  Definitions of symbolic names in parameter statement
C-         IC2T = Column length of zero padded array (Z-direction)
C-         NPH = Number of points in PHI direction of processed data
C
      COMPLEX DATA
C
      CHARACTER IDI*40, IDW*70, IDAYHR*10, filein*10
      CHARACTER*100 idp1p, idp1z, idp2p, idp2z
      CHARACTER*80 HEAD, IDF(2), ID1P, ID1Z, ID2P, ID2Z, IDC
      CHARACTER*1 ISCAN, ICOMP
C
      DIMENSION NN(2), NM(1), WORK(nph*32)
      DIMENSION DATA(NPH*IC2T), TNG(nph*32), TNGR(nph*32)
      DIMENSION RPHI1(nph*2), RPHI2(nph*2), RZ1(nph*2), RZ2(nph*2)
      DIMENSION DAT1(nph*2), DAT2(nph*2), TNG1(nph*2), TNG2(nph*2)

      COMMON/FARFLD/ FREQ, ELMAX, ELMIN, DELZ, DELP, FPHI, FPHIMIN,
     1        PHSHFT, EMAX, NMAX, IPHIM, JMAX, JMIN, IZMAX, IPHMAX,
     1        GAINP, PHIPMAX, THETPMAX, GAINZ, PHIZMAX, THETZMAX, STR
      COMMON /DDID /IDI
      COMMON/FRAME/ NOFRAME
      COMMON /PLOT/ IPLT
C
      EQUIVALENCE (TNG(1), TNGR(1)),
     +            (DATA(1), RPHI1(1)), (DATA(nph*2+1), RPHI2(1)),
     +            (DATA(nph*4+1), RZ1(1)), (DATA(nph*6+1), RZ2(1)),
     +            (DATA(nph*8+1), DAT1(1)), (DATA(nph*10+1), DAT2(1)),
     +            (DATA(nph*12+1), TNG1(1)), (DATA(nph*14+1), TNG2(1)),
     +            (DATA(nph*16+1), TNG(1))

      PI=4.*ATAN(1.0)
      DTR = PI/180.
      RTD = 180./PI
      CEE = .2997925
      ISF = 0
      LUC = 0
      PIX2 = 2.*PI
      NOFRAME = 0
C
C-READ AND PRINT ---
C-           NRDATA = Number of rows (points per scan) of data
C-           NCDATA = Number of columns (scans) of data
C-             NPHI = Number of phi points to complete a 360 degree scan
C-           IC2TON = Size of zero filled array for the z-direction
C-             DELZ = Near field measurment spacing, Z direction (cm)
C-             FREQ = Near field measurment frequency in GHz
C-             DIAM = Diameter of test antenna
C-             DIST = Antenna - horn separation (cm)
C-              STR = Reflection coefficient of test antenna  
C-          AOBLOG1 = Ratio of input to output, PROBE 1 (nominal Z)
C-          AOBLOG2 = Ratio of input to output, PROBE 1 (nominal PHI)
C-         
C
      PRINT *, ' Input file? '
      READ(5, '(a)') filein
      iplt = 'Y'
      WRITE (6, '(a)') ' Do you want screen plots? (Y,n) '
      READ(5, '(a)') iplt
      OPEN(1,FORM = 'FORMATTED',FILE = filein)
      OPEN(6,FORM = 'FORMATTED',FILE = 'CYLOUT')
      READ  (1, 1680) IDI
      PRINT 1670, IDI
C
 101  READ (1, *, END = 540) NRDATA, NCDATA, MIDROW, MIDCOL, NRMOD,
     +                     NCMOD, NPHI, IC2TON
      PRINT 1500, NRDATA, NCDATA, MIDROW, MIDCOL, NRMOD, NCMOD,
     +                     NPHI, IC2TON 
C
      READ (1, *) DELZ, FREQ, DIAM, DIST, STR, AOBLOG1, AOBLOG2
      PRINT 1510, DELZ, FREQ, DIAM, DIST, STR, AOBLOG1, AOBLOG2
C
      READ (1, *) AMIN, AMAX, ELMIN, ELMAX
      PRINT 1510, AMIN, AMAX, ELMIN, ELMAX
C
      READ (1, *) ISCAN, ICOMP, ZINIT, PHINIT, RELPHS, DLPFF
      PRINT 1710, ISCAN, ICOMP, ZINIT, PHINIT, RELPHS, DLPFF
C
      READ (1, 1680) IDF(1)
      PRINT 1670, IDF(1)
C
      READ (1, 1680) IDF(2)
      PRINT 1670, IDF(2)
C
      READ (1, 1680) ID1P
      PRINT 1670, ID1P
C
      READ (1, 1680) ID1Z
      PRINT 1670, ID1Z
C
      READ (1, 1680) ID2P
      PRINT 1670, ID2P
C
      READ (1, 1680) ID2Z
      PRINT 1670, ID2Z
C
C
C
C-The next line read is a control.
C
C-  ICM  = 1  If converting inches to centimeters.
C-  IAMP = 1  If printing amplitude. (not currently implemented)
C-  IPHS = 1  If printing phase. (not currently implemented)
C- MXAMP = 1  If printing and plotting far-field cuts through main
C-            beam maximum (otherwise cuts are through El=Az=0).
C-  INRM = 1  If data normalized before printing.
C-  IBNG = 1  If printing Bng's (not currently implemented)

      READ (1, *) ICM, IAMP, IPHS, MXAMP, INRM, IBNG
      PRINT 1750, ICM, IAMP, IPHS, MXAMP, INRM, IBNG
C
      IF(NPH*IC2T .LT. NPHI*IC2TON)  THEN
        STOP ' Array DATA is not dimensioned large enough, change NPH an
     +d or IC2T and recompile'
      ENDIF

      DELP = 360./NPHI
      DELZ = DELZ*.01
      DIST = DIST*.01
      DIAM = DIAM*.01
C
      IF(IPLT .EQ. 1)   LUC = 1
      IF (ICM .EQ. 1)  THEN
        DELZ = DELZ*2.54
        DIST = DIST*2.54
        DIAM = DIAM*2.54
      ENDIF
C
      AOB1 = 10. ** (AOBLOG1/20.)
      AOB2 = 10. ** (AOBLOG2/20.)
C
      HEAD = IDF(1)

      IF(NRMOD .LT. 1)  NRMOD = 1
      IF(NCMOD .LT. 1)  NCMOD = 1

C
C     PARAMETERS DESCRIBING THE TEST ANTENNA
C
C     A = minimum radius of cylinder circumscribing the test antenna.
C
C     NMAX is the highest order of Hankel function used in the cylindrical
C     expansion.
C
      IR2TON = 360./DELP + .001
      IR2X2 = 2*IR2TON
      IF (IC2TON .EQ. 0)  IC2TON = NCDATA
C
      EZ0=1.0
      A = DIAM/2.
      WAVLGTH = CEE/FREQ
      XK=2.*PI/WAVLGTH
      XKA=XK*A
      NMX=XKA + 10.
C
C
C
C     PARAMETERS DESCRIBING PROBE
C
C     S01 is the ratio of the amplitude of the probe output B0 to the
C     amplitude of the incident transverse electric field (polarization
C     matched) at the mouth of the ideal dipole probe.   That is,
C     S01 = ABS(B0/ET).  S01 is the receiving element S01(K = 0) of the
C     receiving matrix for the probe defined by Kerns.  Here S01 is
C     chosen as 1.0.
C
C     RHO0 is the distance from the center of the scan cylinder to the
C     mouth of the probe.
C
      S01=1.0
      RHO0 = DIST
C
C
C
C     PARAMETERS DESCRIBING SCAN
C
C     DELZ = distance between scan points along the axis of the cylinder.
C     DELP = angle in degrees between scan points around the cylinder.
C     ZINIT, ZFINAL = Z-coordinates respectively of the initial and final
C           Z points with nonzero amplitude.
C     PHINIT, PHIFIN = PHI-coordinates respectively of the initial and final
C           PHI points (in degrees) with nonzero amplitude.
C
      IDDPHI = DELP + .001
      DPHI = DELP*DTR
      IPHINIT = PHINIT/DELP + .001

      NN(1) = IR2TON
      NN(2) = IC2TON
      IN=NN(1)

      IF(ISCAN .EQ. 'Z')  THEN
        PHICTR = PHINIT + (MIDCOL - 1)*DELP/FLOAT(NCMOD)
        ZCTR = ZINIT + (MIDROW - 1)*DELZ/FLOAT(NRMOD)
        NZPTS = NRDATA
        NPHIPTS = NCDATA
        ZFINAL = ZINIT + (NZPTS - 1)*DELZ/FLOAT(NRMOD)
        PHIFIN = PHINIT + (NPHIPTS - 1)*DELP/FLOAT(NCMOD)
      ELSE
        PHICTR = PHINIT + (MIDROW - 1)*DELP/FLOAT(NRMOD)
        ZCTR = ZINIT + (MIDCOL - 1)*DELZ/FLOAT(NCMOD)
        NZPTS = NCDATA
        NPHIPTS = NRDATA
        ZFINAL = ZINIT + (NZPTS - 1)*DELZ/FLOAT(NCMOD)
        PHIFIN = PHINIT + (NPHIPTS - 1)*DELP/FLOAT(NRMOD)
      ENDIF



C
C     Parameters describing far-field EZ
C
C     The far field EZ is printed between 180+FPHIMIN and 180+FPHIMAX
C     degrees for IPHIM points, therefore separated by increments FPHI.
C     MM is the total number of computed far-field points from 0-360 degrees.
      FPHIMIN = AMIN
      FPHIMAX = AMAX
      IPHIM = (AMAX - AMIN)/DLPFF + 1.001
      FPHI = DLPFF
      MM=360./FPHI+.001
      NM(1) = MM
      MMM1=MM-1
      PHSHFT = AMOD(PHINIT,DELP)
C
      PRINT 5
      PRINT 125,(180.+FPHIMIN),(180.+FPHIMAX)
      PRINT 615, FPHI,MM
      NROWX2 = NRDATA * 2
      LU = 15
      LUOUT = 20
C
      OPEN (11, FORM = 'UNformatted')
      OPEN (12, FORM = 'UNformatted')
      OPEN (13, FORM = 'UNformatted')
      OPEN (14, FORM = 'UNformatted')
      READ (11) IDP1P
      READ (IDP1P, '(A70, 3I10)') IDW, IC2TON, NTot, NGAM

      ntx2 = ntot*2

      PRINT 1535, IDW, IC2TON, NTX2, NGAM
C
      NTOT = NTX2/2
      IGMAX = NGAM/2
      MMLEAST = NTOT
      NMAX = NTOT/2
      NMAXP1=NMAX+1
      IF(NMX .NE. NMAX)  THEN
        print *, 'nmx = ', nmx, '  nmax = ', nmax
        STOP ' Antenna parameters specified for this run do not agree wi
     +th those used in calculating probe coefficients.'
      ENDIF

      IF(NPHI .LT. 2*NMAX+1)  THEN
        STOP ' Not enough data points in phi to calculate the required n
     +umber of modes.'
      ENDIF

      PRINT 1800, IC2TON, NTOT, NGAM, IGMAX, MMLEAST, NMAX
 1800 FORMAT (1X, 8I10)
C
      READ (12) IDP1Z
      READ (13) IDP2P
      READ (14) IDP2Z
C
      CALL TIMER
C-
C- READ IN AND TRANSFORM PHI AND Z-COMPONENT OF INPUT DATA
C-
      DO 100 ICMP = 1, 2
        LUOUT = 21
        IF(ICMP .EQ. 2)  ICOMP = 'X'
        IF (ICOMP .EQ. 'P' .OR. ICOMP .EQ. 'X')  LUOUT = 22
C-
C- FIND PROPER NEAR-FIELD DATA FILE
C-
        IDAYHR = IDF(ICMP)(71:80)
        print *, idayhr
        OPEN(LU, FILE = IDAYHR, FORM = 'UNformatted')
        CALL FINDID(IDC, IDAYHR, LU, 5)
C-
C- Read input data, transpose if necessary and add phase shift to
C- eliminate necessity for rearranging.
C-
        CALL GTCLARY(LU, 0, DATA, NRDATA, NCDATA, NRMOD, NCMOD, IR2TON,
     1   IC2TON, NMAX, IPHINIT, WORK, ISCAN)

        CLOSE(LU)
C-
C- Perform two-dimensional FFT.
C-
        CALL FOURT(DATA, NN, 2, -1, +1, WORK)
C-
C- Write transformed output to disk file.
C-
        CALL WRTTNG(data, nphi, ic2ton, igmax, mmleast, luout)

  100   CONTINUE
      CLOSE(15)
C-
C- Calculate cylindrical mode coefficients for test antenna.
C-
      
      CALL DCUPL(RPHI1, RPHI2, RZ1, RZ2, DAT1, DAT2, TNG1, TNG2,
     1           AOB1, AOB2, MMLEAST, NGAM)
C-
C- Calculate far field from cylindrical coefficients.
C-
      close(11)
      close(12)
      close(13)
      close(14)
      close(21)
      close(22)

      CALL FAREVAL(TNG, TNGR, WORK, NM, NGAM, IC2TON)

      CALL timer
  
C
C     Printout of the main parameters
C
      PRINT 20, Freq, Diam, STR
      PRINT 12
      PRINT 15, Nzpts, Nphipts, Delz, Delp, Zinit, Phinit, 
     +          Zfinal, Phifin, Rho0, Nmax, Wavlgth, WAVLGTH/2., 
     +          Wavlgth*90./(A*Pi)
      PRINT 16, idf(1)(71:80), idf(2)(71:80), idp1p, idp1z, idp2p, idp2z
      IF (NMAXP1 .GT. NN(1))  THEN
        PRINT 85
        PRINT 86
        PRINT 87
      ENDIF
      PRINT 122
      PRINT 123, 'Az component gain = ', GAINP, PHIPMAX, THETPMAX
      PRINT 123, 'El component gain = ', GAINZ, PHIZMAX, THETZMAX
      IF(mxamp .EQ. 1)  THEN
        print 630, 'Far field pattern cuts are through main beam peak.'
      ELSE 
        print 630, 'Far field pattern cuts are through El=0 and Az=0.'
        iphmax = 0
      ENDIF
      NPHIM = IPHIM
      IF(inrm .NE. 1) emax = 0.
      CALL FFPRNT(DATA, IC2TON, HEAD, amax, amin, ngam)

      GO TO 101
  540 CONTINUE
      endfile 6
      STOP ' End of CYLPRBC'

   5  FORMAT(' THETA = 90 DEGREES'//)
    6 FORMAT(/////)
   12 FORMAT (///, T28, 'NEAR-FIELD SCAN PARAMETERS'//)
   15 FORMAT (T12, 'Nz = ', I4, T52, 'Nphi = ', I4,/,
     +      T8, 'DELTAz = ', F8.3, ' meters',
     +      T48, 'DELTAphi = ', F7.2, ' degrees',/,
     +      T8, 'Zfirst = ', F8.3, ' meters',
     +      T48, 'PHIfirst = ', F7.2,' degrees'/,
     +      T9, 'Zlast = ', F8.3, ' meters',
     +      T49, 'PHIlast = ', F7.2, ' degrees',/,
     +      T26, 'Scan Radius = ', F8.4, ' meters',//,
     +      T30, '(Derived Parameters)'//,
     +      T33, 'Nmax = ', I4,/
     +      T27, 'Wavelength = ', F9.4, ' meters',/
     +      T27, 'DELTAz max = ', F9.4, ' meters',/
     +      T25, 'DELTAphi max = ', F7.2, ' degrees')
   16 FORMAT(//T30, 'Near-field Data Files',//T25, 'Probe 1 - ', a10,/
     +       T25, 'Probe 2 - ', a10,//, T28,
     +       'Probe Coefficient File IDs',
     +       //, 1x, a, /,1x,  a, /, 1x, a, /, 1x, a, /)
   20 FORMAT ('1', ///, T32, 'ANTENNA PARAMETERS',//,
     +      T29, 'Frequency = ', F7.3, ' GHz',/,
     +      T30, 'Diameter = ', F7.3, ' meters',/,
     +      T16, 'Reflection Coefficient = ', F7.3)
   40 FORMAT (1X ,2E15.5,5X,2E15.5,I12)
   50 FORMAT (' TIME IN SECONDS IS',F12.5)
   65 FORMAT ('1',13X,'BNG2P',30X,'BNG2M',18X,'ABS VALUE OF ORDER' )
   70 FORMAT (8X,'AMPL',7 X,'PHASE(RADIANS)',10X,'AMPL',7X,'PHASE(RADIAN
     1S)'/)
   85 FORMAT(60X,'NOTE THAT NMAXP1 IS LARGER THAN NN(1).  THUS THE ')
   86 FORMAT(60X,'PROGRAM WILL PROBABLY RUN BUT THE BNG2P(N) AND')
   87 FORMAT(60X,'BNG2M(N) WILL BE ERRONEOUS FROM N=NN(1)+1 TO NMAXP1.')
  122 FORMAT(///, T33, 'FAR-FIELD PARAMETERS'///)
  123 FORMAT(T5, A25, F8.3, ' dB at Phi = ', f9.3, ' and Theta = ',
     1      F9.3, ' Degrees')
  125 FORMAT (' THE FAR-FIELD IS PRINTED FROM',F12.4,'  TO',F12.4,'  DEG
     1REES'/)
  160 FORMAT (///' EZMAX =', F12.3,'  DB   AT PHI =', F12.3,'  DEGREES')
  170 FORMAT (1H1)
  451 FORMAT('   IN .LT. 2*NMAX + 1')
  510 FORMAT('1 PRINTOUT OF REAL AND IMAGINARY PART OF DATA ARRAY'//)
  511 FORMAT('1 PRINTOUT OF REAL AND IMAGINARY PART OF DATA1 ARRAY'//)
  515 FORMAT (1X,3(E15.5,E15.5,I8))
  516 FORMAT(1X,8E14.4)
  615 FORMAT( ' FAR-FIELD INCREMENT IN PHI =',F12.5,' DEGREES, SO NUMBER
     1 OF FAR-FIELD POINTS FROM 0-360 DEGREES =',I8//)
  625 FORMAT('1 EXECUTION STOPPED BECAUSE MM LESS THAN (2 NMAX+1).')
  630 FORMAT (/,t20, a)
 1500 FORMAT(BZ,8(6X,I4))
 1510 FORMAT(BZ,8F10.5)
 1520 FORMAT(5F10.4, 10X, 4I5)
 1530 FORMAT(1X, 8A10, 2I10/)
 1535 FORMAT(1X, A70, 3I10/)
 1540 FORMAT(1X,'SKIPPED EOF',I5)
 1670 FORMAT(1X, A)
 1680 FORMAT(A)
 1710 FORMAT(BZ,2(9X,A1),6F10.4)
 1720 FORMAT(1X,2(9X,A1),6F10.4)
 1740 FORMAT(BZ,80I1)
 1750 FORMAT(1X, 80I1)
 2000 FORMAT(' ID(5) .NE. IDAYHR')
      END


      SUBROUTINE GTCLARY(LU, IEOF, DATA, NROW, NCOL, NRMOD, NCMOD, MD,
     1 ND, NMAX, IPHI0, STOR, ISCAN)
C-
C- THIS ROUTINE READS AN ENTIRE ARRAY OF CYLINCRICAL NEAR-FIELD DATA
C- AND APPLIES THE APPROPRIATE PHASE CORRECTION SO THAT REARRANGEMENT
C- OF THE TRANSFORMED DATA IS UNNECESSARY.  DATA ARE TRANSPOSED IF
C- THEY WERE TAKEN WITH SCANS IN THE Z-DIRECTION.
C-
      DIMENSION STOR(*)

      COMPLEX DATA(MD, ND)
C
      CHARACTER*1 ISCAN
C
      PI = 4.*ATAN(1.)
      PIX2 = PI*2.
      DTR = PI/180.
      C1 = PIX2*NMAX/MD
      C2 = PI*(ND - 2)/ND
C
C
      IF (LU .EQ. 0) RETURN
C
C-BUFFER IN NCOL COLUMNS OF DATA ON LU, SEPARATE AMPLITUDE AND PHASE,
C-CONVERT FROM POLAR TO RECTANGULAR COORDINATES AND STORE ADJACENTLY
C-(AS COMPLEX) IN ARRAY DATA              (1 DEGREE=.01745329252 RADIAN)
C
C- IF DATA WAS TAKEN IN THE Z-DIRECTION, THEN DATA ARRAY IS FILLED ROW
C- BY ROW

C- IF THE DATA WAS TAKEN IN THE PHI-DIRECTION, THEN THE DATA ARRAY IS
C- FILLED COLUMN BY COLUMN
C
C- ZERO OUT ARRAY DATA.
C
      DO 50 IR = 1, MD
      DO 50 IC = 1, ND
   50 DATA(IR, IC) = (0., 0.)
C
      IF (ISCAN .NE. 'Z')  THEN
C
        PRINT 1530
C
        NR = NROW/NRMOD
        IF(MOD(NROW, NRMOD) .NE. 0)  NR = NR + 1
        SUMAMP = 0.0
        IRFST = IPHI0
        ICFST = (ND - NCOL/NCMOD)/2
        IF(MOD(NCOL, NCMOD) .NE. 0)  ICFST = ICFST - 1
        IF(ICFST .LT. 1)  ICFST = 1
        ICOL = ICFST
C
        IF (NCOL .LT. 1) GO TO 150
        DO 140 IC = 1,NCOL
          ITST = MOD(IC + 1, NCMOD)
C
          CALL RDCOLU(STOR, NROW * 2, LU)
          IF(ITST .EQ. 0)  THEN
          ICOL = ICOL + 1
C  
          IF (NROW .LT. 1) GO TO 130
          DO 110 IR = 1, NR
 
            IDTA = 1 + (IR - 1)*NRMOD
            IDPNR = IDTA + NROW
            IROW = IR + IRFST
            IRPNR = IR + NROW
            TAMP = STOR(IDTA)
            ANGLE = STOR(IDPNR)*DTR + C1*(IROW - 1) + C2*(ICOL - 1)
            SUMAMP = SUMAMP + TAMP * * 2
            IRPNR1 = IRPNR - 1
C      IF(STOR(IRPNR).EQ.STOR(IRPNR1))PRINT1560,STOR(IRPNR),IROW,ICOL
C     IF(IROW .EQ. 1) PRINT 1520, TAMP,ANGLE
C     IF(IROW .EQ. NROW) PRINT 1520,TAMP,ANGLE
C
            IRX2 = IROW*2
            IRX2M1 = IRX2 - 1
            DATA(IROW, ICOL) = TAMP*CMPLX(COS(ANGLE), SIN(ANGLE))
C
 110      CONTINUE
 130      CONTINUE
        ENDIF
 140    CONTINUE
        PRINT 1580, NCOL, NROW*2, LU
 150    CONTINUE
C
        PRINT 1570, SUMAMP
C
        RETURN
C
      ENDIF
C
      PRINT 1540
C
      SUMAMP = 0.0
      NR = NROW/NRMOD
      IRFST = IPHI0
      IROW = IRFST
      ICFST = (ND - NROW/NRMOD)/2
      IF (MOD(NROW, NRMOD) .NE. 0)  THEN
        NR = NR + 1
        ICFST = ICFST - 1
      ENDIF
      IF (ICFST .LT. 1)   ICFST = 1

      IF (NCOL .LT. 1)  GO TO 220
      DO 210 IR = 1, NCOL

        ITST = MOD(IR + 1, NCMOD)

        CALL RDCOLU(STOR, NROW*2, LU)
C
        IF (ITST .EQ. 0)  THEN
        IROW = IROW + 1
        IRX2 = IROW*2
        IRX2M1 = IRX2 - 1

        IF (NR .LT. 1)   GO TO 200
        DO 190 IC = 1, NR

          IDTA = 1 + (IC - 1)*NRMOD
          IDPNR = IDTA + NROW
          ICOL = IC + ICFST

          TAMP = STOR(IDTA)

          SUMAMP = SUMAMP + TAMP * * 2
          ANGLE = STOR(IDPNR)*DTR + C1*(IROW - 1) + C2*(ICOL - 1)
C
          DATA(IROW, ICOL) = TAMP*CMPLX(COS(ANGLE), SIN(ANGLE))
C
 190    CONTINUE
 200  CONTINUE
      ENDIF
 210  CONTINUE
      PRINT 1580, NCOL, NROW*2, LU
 220  CONTINUE
      PRINT 1570, SUMAMP
C
      RETURN
1500  FORMAT  (I10)
1510  FORMAT  ('0' ,2F12.5 )
 1530 FORMAT(20X, 'ANTENNA SCANNED IN THE PHI-DIRECTION, SO DATA NOT
     1TRANSPOSED IN GTCLARY')
 1540 FORMAT(20X, 'ANTENNA SCANNED IN THE Z-DIRECTION, SO DATA WAS
     1TRANSPOSED IN GTCLARY')
1550  FORMAT  (5X, I5, 3F12.5)
1560  FORMAT  (10X,F10.3,'AT ROW NO.',I5,'AT COL NO.',I5)
1570  FORMAT (' SUM OF SQUARES OF INPUT AMPLITUDES = ', F20.4,//)
1580  FORMAT (/// 2X, I5, 5X, I5, ' WORD RECORDS IN FROM UNIT ', I5///)
1590  FORMAT(2X, I5, 5X, I5, ' WORD RECORD IN FROM LU ',I4)
C
      END


      SUBROUTINE DCUPL(RPHI1, RPHI2, RZ1, RZ2, DAT1, DAT2, TNG1, TNG2,
     1 AOB1, AOB2, NROW, NCOL)
C-
C- THIS SUBROUTINE READS PROBE RECEIVING COEFFICIENTS FROM UNITS 11,
C- 12, 13 AND 14 AND THE TRANSFORMED INPUT DATA FROM UNITS 21 AND 22.
C- IT CALCULATES THE TRANSMITTING COEFFICIENTS AND WRITES THEM ON
C- UNITS 23 AND 24.
C-
      COMPLEX RPHI1, RZ1, RPHI2, RZ2, DAT1, DAT2, TNG1, TNG2, detr

      DIMENSION RPHI1(*), RZ1(*), RPHI2(*), RZ2(*),
     +       DAT1(*), DAT2(*), TNG1(*), TNG2(*)
C-
      NRX2 = NROW*2
      DO 10 IGAM = 1, NCOL
        CALL RDCOLU (RPHI1, NRX2, 11)
        CALL RDCOLU (RZ1, NRX2, 12)
        CALL RDCOLU (RPHI2, NRX2, 13)
        CALL RDCOLU (RZ2, NRX2, 14)
        CALL RDCOLU (DAT1, NRX2, 21)
        CALL RDCOLU (DAT2, NRX2, 22)
        DO 20 IPHI = 1, NROW
          DETR = RZ2(IPHI)*RPHI1(IPHI) - RZ1(IPHI)*RPHI2(IPHI)
          IF(CABS(DETR) .NE. 0.)  THEN
            TNG1(IPHI) = (RZ2(IPHI)*DAT1(IPHI)/AOB1 -
     +                    RZ1(IPHI)*DAT2(IPHI)/AOB2)/DETR
            TNG2(IPHI) = (RPHI1(IPHI)*DAT2(IPHI)/AOB2 - 
     +                    RPHI2(IPHI)*DAT1(IPHI)/AOB1)/DETR
          ELSE
            TNG1(IPHI) = (0., 0.)
            TNG2(IPHI) = (0., 0.)
          ENDIF
   20     CONTINUE
        CALL WRTONLU (TNG1, NRX2, 23)
        CALL WRTONLU (TNG2, NRX2, 24)
   10   CONTINUE
      REWIND 23
      REWIND 24
 3000 FORMAT(i10, 2e20.4)
      RETURN
      END


      SUBROUTINE FAREVAL(TNG, TNGR, WORK, NM, NGAM, IC2TON)
C-
C- THIS SUBROUTINE SUMS THE CYLINDRICAL MODES FOR THE THETA AND PHI
C- COMPONENTS TO OBTAIN FAR-FIELD PATTERNS.
C-
      COMMON/FARFLD/ FREQ, ELMAX, ELMIN, DELZ, DELP, FPHI, FPHIMIN,
     1       PHSHFT, EMAX, NMAX, IPHIM, JMAX, JMIN, IZMAX, IPHMAX, 
     1       GAINP, PHIPMAX, THETPMAX, GAINZ, PHIZMAX, THETZMAX, STR
C
      DIMENSION TNG(*), TNGR(*), NM(*), WORK(*)
C
      COMPLEX TNG, EZ, EMIP, COEF
C
      CHARACTER*10 IDENT, head*100
C-
C- MISCELLANEOUS
C-
      PI = 4.*ATAN(1.)
      PIX2 = 2.*PI
      DTR = PI/180.
      RTD = 180./PI
      CEE = .2997925
      Z0 = 376.73
      WAVLGTH = CEE/FREQ
      XK = PIX2/WAVLGTH
      RFACT = 1. - STR*STR
      MM = NM(1)
      NGMID = (NGAM + 1)/2
      IF (NGAM .EQ. IC2TON)   NGMID = IC2TON/2
      NTOT = 2*NMAX + 1
      THMIN = DTR*(90. - ELMAX)
      THMAX = DTR*(90. - ELMIN)
      JMIN = INT(IC2TON*DELZ*COS(THMax)/WAVLGTH + NGMID + .0001) + 1
      JMAX = INT(IC2TON*DELZ*COS(THMin)/WAVLGTH + NGMID + .0001)
      IF (JMIN .LT. 1)  JMIN = 1
      IF (JMAX .GT. NGAM)  JMAX = NGAM
      NSKIP = JMIN - 1
      IDENT = ' EPHIMAX ='
C-
      DO 10 LUIN = 23, 24
      EMAX = 0.
      head(1:7) = 'CYLPRBC'
      if (luin .eq. 23) then
         head(10:20) = 'Az comp'
      else
         head(10:20) = 'El comp'
      endif

      luout = luin + 2
      IF (NSKIP .NE. 0)  THEN
        DO 15 I = 1, NSKIP
          CALL RDCOLU (TNG, 0, LUIN)
   15   CONTINUE
      ENDIF
C-
      ncout = jmax - jmin + 1
      write(head(81:90), '(i10)')  iphim
      write(head(91:100), '(i10)')  ncout
      write(luout)  head(1:100) 

        DO 100 IC = JMIN, JMAX
          CALL RDCOLU (TNG, NTOT*2, LUIN)
          GAMMA = (IC - NGMID)*PIX2/IC2TON/DELZ
          CAPKSQ = XK*XK - GAMMA*GAMMA
          CAPK = SQRT(CAPKSQ)
          EZ = 1.
          EMIP = CMPLX(0., -1.)
          IF (LUIN .EQ. 24)   EZ = EMIP
          COEF2 = DELZ*DELP*DTR/PIX2/PIX2
          COEF = -2.*CAPK*COEF2*EMIP**(NMAX + 1)
C-
          DO 200 L = 1, NTOT
            EZ = EZ*EMIP
            TNG(L) = EZ*COEF*TNG(L)
  200       CONTINUE
          DO 250 L = NTOT-1, MM
  250       TNG(L) = (0., 0.)
C-
          IM = (1. + FPHIMIN/180.)*MM/2. + .0001
          NPHSHFT = IFIX(PHSHFT/FPHI + .0001)
          CALL FOURT(TNG, NM, 1, 1, 1, WORK)
C-
          DO 300 I = 1, IPHIM
  300       TNG(I) = TNG(IM + I - NPHSHFT)
 
          CALL ARAYRTP(TNGR, MM*2, 1)
C-
          DO 350 IPHI = 1, IPHIM
            IPHIX2 = 2*IPHI
            IPHX2M1 = IPHI*2 - 1
            PHI = PI*(1. + FPHIMIN/180.) + FPHI*(IPHI - 1)*DTR
            PHID = PHI*RTD
            COEFF = -NMAX*(PHID - 90.) + (IC - NGMID)*180.*(IC2TON - 2)
     1               /IC2TON
            PTEMP = TNGR(IPHIX2) + COEFF
            PTEMP = AMOD(PTEMP, 360.)
            IF (PTEMP .LT. 0.)   PTEMP = PTEMP + 360.
            TNGR(IPHIX2) = PTEMP
C
C-   FIND MAXIMUM VALUE OF FAR-FIELD
C
            IF (TNGR(IPHX2M1) .GT. EMAX)   THEN
              EMAX = TNGR(IPHX2M1)
              IZMAX = IC
              IPHMAX = IPHI
            ENDIF
  350       CONTINUE
            CALL WRTONLU (TNGR, IPHIM*2, LUOUT)
  100     CONTINUE
          IF (LUIN .EQ. 24)   IDENT = '  EZMAX = '
          GAIN = EMAX*EMAX*2.*PIX2/Z0/RFACT
          GAIN = 10.*ALOG10(GAIN)
          EMAX = 20.*ALOG10(EMAX)
          IF(luin .EQ. 23)  THEN
            emax1 = emax
            izmax1 = izmax
            iphmax1 = iphmax
          ELSE
            emax2 = emax
            izmax2 = izmax
            iphmax2 = iphmax
          ENDIF
          PHIDMAX = 180. + FPHIMIN + FPHI*(IPHMAX - 1)
          THETMAX = WAVLGTH*(IZMAX - NGMID)/DELZ/IC2TON
          THETMAX = ACOS(THETMAX)*RTD
          IF (LUIN .eq. 23) THEN
            GAINP = GAIN
            PHIPMAX = PHIDMAX
            THETPMAX = THETMAX
          ELSE
            GAINZ = GAIN
            PHIZMAX = PHIDMAX
            THETZMAX = THETMAX
          ENDIF
          PRINT 1010, IDENT, GAIN, PHIDMAX, THETMAX
          rewind luout
   10   CONTINUE
      If (emax1 .GT. emax2)  THEN
        emax = emax1
        izmax = izmax1
        iphmax = iphmax1
      ELSE
        emax = emax2
        izmax = izmax2
        iphmax = iphmax2
      ENDIF

      RETURN
 1010 FORMAT(//// A10,     F12.3, ' DB  AT PHI =', F12.3, ' AND THETA ='
     1 , F12.3, '  DEGREES')
      END


      SUBROUTINE FFPRNT(DATAF, IC2TON, HEAD, amax, amin, ngam)

      DIMENSION DATAF(*), WA(1440), WE(1024)
      DIMENSION FFA(1440, 2), FPA(1440, 2), FFE(1024, 2), FPE(1024, 2)
      CHARACTER HEAD*80, xlabel*30, id*100

      COMMON/FARFLD/ FREQ, ELMAX, ELMIN, DELZ, DELP, FPHI, FPHIMIN,
     +         PHSHFT, EMAX, NMAX, IPHIM, JMAX, JMIN, IZMAX, IPHMAX,
     +         GAINP, PHIPMAX, THETPMAX, GAINZ, PHIZMAX, THETZMAX, STR

C-   Plot azimuth patterns

      LUOUTR = 24
      JTOT = JMAX - JMIN + 1
      NGMID = IZMAX - JMIN + 1
      PI = 4.*atan(1.)
      RTD = 180./PI
      CEE = .2997925
      WAVLGTH = CEE/FREQ
      IF(iphmax .EQ. 0)  THEN
        iphmax = (iphim + 1)/2
        ngmid = (jtot + 1)/2
      ENDIF
      IR = 2*IPHMAX - 1
      IRP1 = IR + 1
C
      HEAD(41:50) = '  A-COMP  '
C
      DO 100 I = 1, 2
C
        luout = luoutr + i
        rewind luout
        read (luout) id
        IF (I .EQ. 2)   HEAD(41:50) = '  E-COMP  '
        DO 110 IC = 1, JTOT
          CALL RDCOLU(DATAF, IPHIM*2, LUOUT)
          DLP = (AMAX - AMIN)/FLOAT(IPHIM - 1)
          IF(IC .EQ. NGMID)  THEN
          DO 600 IPHI = 1, IPHIM
            FFA(IPHI, I) = DATAF(2*IPHI - 1)
            IF(FFA(IPHI, I) .NE. 0.)  THEN
              FFA(IPHI, I) = 20.*ALOG10(FFA(IPHI, I)) - EMAX
            ELSE
              FFA(IPHI, I) = -100.
            ENDIF
            FPA(IPHI, I) = DATAF(2*IPHI)
            WA(IPHI) = AMIN + DLP*FLOAT(IPHI - 1)
  600     CONTINUE
        ENDIF

C-   Plot elevation patterns

          JZ = IC
          IF (DATAF(IR) .NE. 0.)  THEN
            FFE(JZ, I) = 20.*ALOG10(DATAF(IR)) - EMAX
          ELSE
            FFE(JZ, I) = -100.
          ENDIF
          FPE(JZ, I) = DATAF(IRP1)
          csth = (jz-1+jmin-(ngam+1)/2)*wavlgth/delz/ic2ton
          we(jz) = 90. - acos(csth)*rtd
  110   CONTINUE
C
          HEAD(51:60) = '  AZ-cut  '
          HEAD(61:70) = '          '
          xlabel = 'Azimuth - degrees'
          CALL CRTPLT3(WA, FFA(1, I), AMAX, AMIN,  0.  ,  -50.   ,
     +               IPHIM, HEAD, 1, 1, 1, xlabel)
          CALL CRTPLT3(WA, FPA(1, I), AMAX, AMIN,  360.  ,    0. ,
     +               IPHIM, HEAD, 1, 1, 1, xlabel)
 
          HEAD(51:60) = '  EL-cut  '
          NTHPTS = JMAX - JMIN + 1

          xlabel = 'Elevation - degrees'
          CALL CRTPLT3(WE, FFE(1, I), ELMAX, ELMIN,     0.,    -100.,
     +                 NTHPTS, HEAD, 1, 1, 1, xlabel)
          CALL CRTPLT3(WE, FPE(1, I), ELMAX, ELMIN,   360.,      0.,
     +                  NTHPTS, HEAD, 1, 1, 1, xlabel)
  100 CONTINUE
      PRINT 1001
      PRINT 615
      PRINT 620, (WE(JZ), FFE(JZ, 1), FPE(JZ, 1),
     +            FFE(JZ, 2), FPE(JZ, 2), JZ = 1, JTOT)
      PRINT 1000
      PRINT 610 
      PRINT 620, (WA(IPHI), FFA(IPHI, 1), FPA(IPHI, 1),
     +            FFA(IPHI, 2), FPA(IPHI, 2), IPHI = 1, IPHIM)

      RETURN
    6 FORMAT(/////)
  610 FORMAT(T28, 'Azimuth Component', T57, 'Elevation Component',//
     +       T9,'Azimuth', T25, 'Amplitude ',
     +       T42, 'Phase', T55, 'Amplitude', T72, 'Phase',/, 
     +       T9, 'Degrees', T29, 'dB', T41, 'Degrees', 
     +       T59, 'dB', T71, 'Degrees',//)
  615 FORMAT(T28, 'Azimuth Component', T57, 'Elevation Component',//
     +       T8,'Elevation', T25, 'Amplitude ',
     +       T42, 'Phase', T55, 'Amplitude', T72, 'Phase',/, 
     +       T9, 'Degrees', T29, 'dB', T41, 'Degrees', 
     +       T59, 'dB', T71, 'Degrees',//)
  620 FORMAT(1X, F15.3, F15.3, F15.2, F15.3, F15.2)
 1000 FORMAT(1H1, T25, ' FAR-FIELD PATTERN - Azimuth Cut'//)
 1001 FORMAT(1H1, T25, ' FAR-FIELD PATTERN - Elevation Cut'//)
      END


        SUBROUTINE WRTTNG(data, nphi, ic2ton, igmax, mmleast, luout)

        DIMENSION DATA(2*NPHI, IC2TON)

        I1 = IC2TON/2 - IGMAX
        I2 = IC2TON/2 + IGMAX
        DO 350 IC = I1, I2
          CALL WRTONLU(DATA(1, IC), 2*MMLEAST, LUOUT)
  350     CONTINUE
        REWIND LUOUT
        RETURN
        END
