C        COPYRIGHT 1985 UNIVERSITY OF GLASGOW                           0001.000
C        FAST FOURIER TRANSFORM PROGRAM FFT                             0002.000
C        GX PACKAGE  APRIL 1985   K.W. MUIR                             0003.000
C                                                                       0004.000
      PROGRAM NFFTKWM                                                   0005.000
C     MODIFICATION OF PROGRAM BY L. F. TEN EYCK                         0006.000
C     REVISED BY KWM SEP 1981. EXPAND HAS BEEN REWRITTEN TO PROVIDE     0007.000
C     A WIDER RANGE OF MAP TYPES AND TO USE THE MODEL FILE.             0008.000
C     THE MAP IS NOW WRITTEN IN SECTIONS.                               0009.000
C     PART 1. THIS SECTION IS MAINLY CONCERNED WITH INPUT AND OUTPUT.   0010.000
C     ------------------------------------------------------------------0011.000
C     READ OUTPUT TAPE FROM MIRAGE &  EXPAND DATA TO SPACE GROUP P1     0012.000
C     BEFORE CALLING FOURIER TRANSFORM ROUTINES                         0013.000
      CHARACTER COMPID(72)                                              0014.000
      CHARACTER*8 LABEL(200),FTYPE(13)                                  0015.000
      INTEGER CN(200),CONTNT(13)                                        0016.000
      COMMON X(30000)                                                   0017.000
      COMMON /MODELC/COMPID,LABEL,FTYPE                                 0018.000
      COMMON /MODEL/ NATOM,NR,NT,NTYPE,NFTYPE(200),P(3,200),            0019.000
     1CN,SOF(200),UIJ(6,200),EP(3,200),ESOF(200),EUIJ(6,200),           0020.000
     2WAVEL,CELL(6),ECELL(6),R(24,3,4),T(3,4),ICENT,                    0021.000
     3SFAC(14,13),CONTNT,SCALE(4),ETA,NQ,EXTNCT                         0022.000
      COMMON /MAP/IMAP(8000)                                            0023.000
      COMMON /REFL/ IH(100),G(2,100)                                    0024.000
      COMMON /PERIF/ NMAP, LMAP,LT,LO,LSCRA,TWOPI,LM,LPP                0025.000
      COMMON/CONST/NEX,KUSER1,DTOR,LC,LIN,NP(3)                         0026.000
      COMMON /TRAN/ RMAX,RMIN,RSCALE,FZERO   ,MH(3),LEVEL,IT4,XLIM(3)   0027.000
      COMMON/TAPE/NOSET,PARA(6),NSYM,NCENT,LATT,NPIC,TS(3,24),IS(2,3,24)0028.000
     1,DUMMY(49)                                                        0029.000
C                                                                       0030.000
C     PRELIMINARY SECTION        DEFINE ALL NECESSARY PROGRAM PARAMETERS0031.000
C                                                                       0032.000
C     SET LOGICAL NUMBERS FOR PERIPHERALS                               0033.000
C     LC = CARD READER                                                  0034.000
      LC = 5                                                            0035.000
C     LSCRA = FFT SCRATCH TAPE                                          0036.000
      LSCRA = 8                                                         0037.000
C     LT = LINE PRINTER                                                 0038.000
      LT = 7                                                            0039.000
C     LIN = INPUT FILE FROM MIRAGE                                      0040.000
      LIN = 11                                                          0041.000
C     LO = FILE OF EXPANDED DATA ( MAY BE THE SAME AS LMAP )            0042.000
      LO = 13                                                           0043.000
C     LMAP = OUTPUT FILE FOR INTERPRETATION                             0044.000
      LMAP = 13                                                         0045.000
C     PROMPT DEVICE                                                     0046.000
      LPP=6                                                             0047.000
C     MODEL FILE                                                        0048.000
      LM=9                                                              0049.000
C     DEFINE CONSTANTS                                                  0050.000
C     KUSER1 = DIMENSION OF ARRAYS IHKL,E,IPHAZ,EP = MAXIMUM NUMBER     0051.000
C     OF INDEPENDENT REFLECTIONS WHICH CAN BE INPUT  (KUSER1 .GE. NUMB) 0052.000
C     DIMENSION OF X = 3 * KUSER1                                       0053.000
      KUSER1 = 10000                                                    0054.000
      NSIZE = 3 * KUSER1 - 1                                            0055.000
      DTOR = ATAN(1.0) / 45.0                                           0056.000
      TWOPI = 360.0 * DTOR                                              0057.000
      CALL EXPAND                                                       0058.000
      CALL PP1(X,NSIZE,NP(1),NP(2),NP(3))                               0059.000
      STOP                                                              0060.000
      END                                                               0061.000
C     ----------------------------------------------------------------  0062.000
C     MAKE NULL TAPE FOR SEARCH                                         0063.000
      SUBROUTINE ERROR                                                  0064.000
      COMMON /PERIF/ NMAP, LMAP,LT,LO,LSCRA,TWOPI,LM,LPP                0065.000
      I = 0                                                             0066.000
      NDIFF = -1                                                        0067.000
      WRITE (LMAP) (I,J=1,20),NDIFF,(I,J=1,226)                         0068.000
      IF(I)1,2,1                                                        0069.000
    2 STOP                                                              0070.000
    1 RETURN                                                            0071.000
      END                                                               0072.000
C     ------------------------------------------------------------------0073.000
C     INTERACTIVE DATA INPUT. EXPANSION OF DATA TO FULL HEMISPHERE      0074.000
      SUBROUTINE EXPAND                                                 0075.000
      INTEGER TITLE(68)                                                 0076.000
      CHARACTER COMPID(72),LINE(148)                                    0077.000
      CHARACTER*8 LABEL(200),FTYPE(13),MENU(7),TAG,FMAP(7),ITAPE(3)     0078.000
      INTEGER CN(200),CONTNT(13),FCH(20),FWD(20)                        0079.000
      DIMENSION PC(3),CFOM(2250),EQ(10000),LH(60),LK(60),LL(60),        0080.000
     1EW(60),ED(60),IHKL(10000),E(10000),IPHAZ(10000)                   0081.000
      COMMON IHKL,E,IPHAZ                                               0082.000
      COMMON/MODELC/COMPID,LABEL,FTYPE                                  0083.000
      COMMON/MODEL/NATOM,NR,NT,NTYPE,NFTYPE(200),P(3,200),              0084.000
     1CN,SOF(200),UIJ(6,200),EP(3,200),ESOF(200),EUIJ(6,200),           0085.000
     2WAVEL,CELL(6),ECELL(6),R(24,3,4),T(3,4),ICENT,SFAC(14,13),        0086.000
     3CONTNT,SCALE(4),ETA,NQ,EXTNCT                                     0087.000
      COMMON/REFL/IH(100),G(2,100)                                      0088.000
      COMMON/PERIF/ NMAP,LMAP,LT,LO,LSCRA,TWOPI,LM,LPP                  0089.000
      COMMON/CONST/NEX,KUSER1,DTOR,LC,LIN,NP(3)                         0090.000
      COMMON/TRAN/ RMAX,RMIN,RSCALE,FZERO,MH(3),LEVEL,IT4,XLIM(3)       0091.000
      COMMON/TAPE/NOSET,PARA(6),NSYM,NCENT,LATT,NPIC,TS(3,24),          0092.000
     1IS(2,3,24),TAB(15),KIND(24),I1(3),I2(3),IP5,FSUM,NUMB,NBATCH      0093.000
      EQUIVALENCE  (IPHAZ(1),EQ(1),CFOM(1))                             0094.000
      EQUIVALENCE (IH(1),LH(1)),(IH(61),LK(1)),(G(1,11),LL(1)),         0095.000
     1  (G(1,41),EW(1)),(G(1,71),ED(1))                                 0096.000
      DATA LINE/148*' '/,MENU/'FMAP','FGRID','FSET','FPRINT',           0097.000
     1'FSCALE','FZERO','FREJ'/,FMAP/'DF','WDF','FOBS','FCALC','FSQ',    0098.000
     2'ESQ','EMAP '/,NMEN/7/,ITAPE/'L.S.','CAD4','EMAP '/               0099.000
      LK10=2**10                                                        0100.000
      LK16=2**16                                                        0101.000
      LK7=2**7                                                          0102.000
      LK8=2**8                                                          0103.000
 1000 FORMAT('0CHOOOSE INPUT OPTIONS FROM THE MENU',//                  0104.000
     1' FMAP     DF/FOBS/FCALC/FSQ ',/                                  0105.000
     2' FGRID   <GRID IN ANGSTROMS>',/                                  0106.000
     4' FPRINT  <XMAX> <YMAX> <ZMAX> <<LEVEL>>',/                       0107.000
     5' FSCALE  <DENSITY SCALE>',/                                      0108.000
     6' FZERO   <F(OOO)>',/                                             0109.000
     7' FREJ    < SIG  <A>> < MIN  <B>> < MAX  <C>> < PHASE  <D>>'/     0110.000
     8' MENU    TO OBTAIN A FRESH COPY OF THE MENU'/                    0111.000
     9' <RETURN> STARTS THE CALCULATION'/)                              0112.000
 1001 FORMAT(//' TYPE AN INSTRUCTION')                                  0113.000
 1002 FORMAT(72A1)                                                      0114.000
 1003 FORMAT(/ 1H  ,14X,A5,' MAP TO BE CALCULATED')                     0115.000
 1004 FORMAT(/ 6X,' THE GRID SPACING IS APPROXIMATELY',F6.3,' ANGSTROMS'0116.000
     1)                                                                 0117.000
 1005 FORMAT(/  9X,' PHASE SET',I3,' TO BE USED FOR EMAPS')             0118.000
 1006 FORMAT(/  9X,' THE MAP WILL NOT BE PRINTED')                      0119.000
 1007 FORMAT(/ 9X,' THE MAP WILL BE PRINTED FROM X=0.0 TO X=',F5.3,/,38X0120.000
     1,' Y=0.0 TO Y=',F5.3,/38X,' Z=0.0 TO Z=',F5.3/,                   0121.000
     29X,' POINTS GREATER THAN ',I4,' WILL BE UNDERLINED'/)             0122.000
 1008 FORMAT(/  9X,' THE SCALE FACTOR WILL BE CHOSEN BY THE PROGRAM')   0123.000
 1009 FORMAT(/  12X,' THE SCALE FACTOR IS ',F15.5)                      0124.000
 1010 FORMAT(/  9X,' F(000) WILL BE CHOSEN BY THE PROGRAM')             0125.000
 1011 FORMAT(/ 19X,' F(000)=',F8.2)                                     0126.000
 1012 FORMAT(/ 11X,' A REFLECTION WILL BE REJECTED IF:',/               0127.000
     110X,' (1) ITS INTENSITY IS LESS THAN ',F5.2,' TIMES ITS E.S.D.',/ 0128.000
     210X,' (2) SINTHETA/LAMBDA IS LESS THAN ',F5.3,                    0129.000
     3' OR GREATER THAN ',F5.3,/                                        0130.000
     410X,' (3) ABS(FCALC/FOBS) IS LESS THAN ',F5.3)                    0131.000
 1013 FORMAT(' FFT PROGRAM. VERSION OF JAN 1984',/1H0,72A1,             0132.000
     1//8X,' REFLECTION DATA WILL BE READ FROM A ',A6,' FILE')          0133.000
C                                                                       0134.000
C                                                                       0135.000
    1 NMAP=1                                                            0136.000
      FACTOR=1.0/3.0                                                    0137.000
      NOSET=0                                                           0138.000
      LP=LPP                                                            0139.000
      IT4=0                                                             0140.000
      LEVEL=100                                                         0141.000
      DO 2 I=1,3                                                        0142.000
      MH(I)=0                                                           0143.000
    2 XLIM(I)=0.                                                        0144.000
      RSCALE=-1.0                                                       0145.000
      FZERO=-1.0                                                        0146.000
      FT1=3.0                                                           0147.000
      FT2=0.0                                                           0148.000
      FT3=1.5                                                           0149.000
      FT4=0.0                                                           0150.000
      NPIC=0                                                            0151.000
      ASSIGN 4 TO NJUMP                                                 0152.000
C                                                                       0153.000
C                                                                       0154.000
C                                                                       0155.000
    3 WRITE(LP,1000)                                                    0156.000
    4 WRITE(LP,1001)                                                    0157.000
    5 READ(LC,1002,END=25)(LINE(I),I=2,73)                              0158.000
      CALL PARSE(LINE,NFS,FCH,FWD)                                      0159.000
      IF(NFS.EQ.0)GOTO25                                                0160.000
      CALL AFORMT(LINE,FCH(1),FWD(1),TAG)                               0161.000
      IF(TAG.EQ.'MENU')GOTO3                                            0162.000
      NTAG=NMEN+1                                                       0163.000
      DO 6 I=1,NMEN                                                     0164.000
    6 IF(TAG.EQ.MENU(I))NTAG=I                                          0165.000
      IF(NFS.EQ.1)GOTO7                                                 0166.000
      GOTO(8,11,13,15,18,20,22,4),NTAG                                  0167.000
    7 GOTO(10,12,14,17,19,21,24,4),NTAG                                 0168.000
C                                                                       0169.000
C     FMAP    *******************************************************   0170.000
    8 CALL AFORMT(LINE,FCH(2),FWD(2),TAG)                               0171.000
      DO 9 I=1,7                                                        0172.000
    9 IF(TAG.EQ.FMAP(I))NMAP=I                                          0173.000
   10 WRITE(LP,1003)FMAP(NMAP)                                          0174.000
      GOTO NJUMP                                                        0175.000
C                                                                       0176.000
C     FGRID   ********************************************************  0177.000
   11 CALL FFORMT(LINE,FCH(2),FWD(2),FACTOR)                            0178.000
   12 WRITE(LP,1004)FACTOR                                              0179.000
      GOTO NJUMP                                                        0180.000
C                                                                       0181.000
C     FSET    **********************************************************0182.000
   13 CALL IFORMT(LINE,FCH(2),FWD(2),NOSET)                             0183.000
   14 WRITE(LP,1005)NOSET                                               0184.000
      GOTO NJUMP                                                        0185.000
C                                                                       0186.000
C     FPRINT                                                            0187.000
   15 IF(NFS.LT.4)GOTO17                                                0188.000
      DO 16 I=1,3                                                       0189.000
   16 CALL FFORMT(LINE,FCH(I+1),FWD(I+1),XLIM(I))                       0190.000
      IF(NFS.GE.5)CALL IFORMT(LINE,FCH(5),FWD(5),LEVEL)                 0191.000
      IT4=LT                                                            0192.000
      IF(XLIM(1).LT..0001.AND.XLIM(2).LT..0001.AND.XLIM(3).LT..0001)IT4=0193.000
     10                                                                 0194.000
   17 IF(IT4.EQ.0)WRITE(LP,1006)                                        0195.000
      IF(IT4.NE.0)WRITE(LP,1007)XLIM,LEVEL                              0196.000
      GOTO NJUMP                                                        0197.000
C                                                                       0198.000
C     FSCALE   *********************************************************0199.000
   18 CALL FFORMT(LINE,FCH(2),FWD(2),RSCALE)                            0200.000
   19 IF(RSCALE.LT.0.)WRITE(LP,1008)                                    0201.000
      IF(RSCALE.GE.0.)WRITE(LP,1009)RSCALE                              0202.000
      GOTO NJUMP                                                        0203.000
C                                                                       0204.000
C     FZERO   ********************************************************* 0205.000
   20 CALL FFORMT(LINE,FCH(2),FWD(2),FZERO)                             0206.000
   21 IF(FZERO.LT.0.)WRITE(LP,1010)                                     0207.000
      IF(FZERO.GE.0.)WRITE(LP,1011)FZERO                                0208.000
      GOTO NJUMP                                                        0209.000
C                                                                       0210.000
C     FREJ    **********************************************************0211.000
   22 IF(MOD(NFS,2).NE.1)GOTO24                                         0212.000
      I=2                                                               0213.000
   23 CALL AFORMT(LINE,FCH(I),FWD(I),TAG)                               0214.000
      CALL FFORMT(LINE,FCH(I+1),FWD(I+1),FF)                            0215.000
      IF(TAG.EQ.'SIG')FT1=FF                                            0216.000
      IF(TAG.EQ.'MIN')FT2=FF                                            0217.000
      IF(TAG.EQ.'MAX')FT3=FF                                            0218.000
      IF(TAG.EQ.'PHASE')FT4=FF                                          0219.000
      I=I+2                                                             0220.000
      IF(I.LT.NFS)GOTO 23                                               0221.000
   24 WRITE(LP,1012)FT1,FT2,FT3,FT4                                     0222.000
      GOTO NJUMP                                                        0223.000
C                                                                       0224.000
C     READ THE MODEL FILE IF NOT A EMAP  RUN                            0225.000
   25 IF(NMAP.GT.5)GOTO 33                                              0226.000
      WRITE(LP,2500)                                                    0227.000
 2500 FORMAT(/10X,' PLEASE WAIT. I AM READING THE MODEL FILE'/)         0228.000
      CALL MDLIN(LM)                                                    0229.000
      NCENT=1                                                           0230.000
      IF(ICENT.EQ.1)NCENT=0                                             0231.000
      LATT=1                                                            0232.000
      IF(NT.EQ.2)GOTO 26                                                0233.000
      IF(NT.EQ.3)LATT=7                                                 0234.000
      IF(NT.EQ.4)LATT=6                                                 0235.000
      GOTO 28                                                           0236.000
   26 LATT=5                                                            0237.000
      DO 27 I=1,3                                                       0238.000
   27 IF(ABS(T(I,2)).LT..1)LATT=I+1                                     0239.000
   28 NSYM=NR                                                           0240.000
      DO 30 J=1,NSYM                                                    0241.000
      DO 29 I=1,3                                                       0242.000
      TS(I,J)=R(J,I,4)                                                  0243.000
      DO 29 K=1,2                                                       0244.000
   29 IS(K,I,J)=0                                                       0245.000
      DO 30 I=1,3                                                       0246.000
      K=1                                                               0247.000
      DO 30 L=1,3                                                       0248.000
      IF(ABS(R(J,I,L)).LT.0.1)GOTO 30                                   0249.000
      IS(K,I,J)=ISIGN(L,NINT(R(J,I,L)))                                 0250.000
      K=2                                                               0251.000
   30 CONTINUE                                                          0252.000
      DO 31 I=1,3                                                       0253.000
      PARA(I)=CELL(I)                                                   0254.000
   31 PARA(I+3)=ACOS(CELL(I+3))/DTOR                                    0255.000
      IF(FZERO.GE.0.0)GOTO 200                                          0256.000
      DO 32 I=1,NTYPE                                                   0257.000
      DO 32 J=1,9,2                                                     0258.000
   32 FZERO=FZERO-CONTNT(I)*SFAC(J,I)                                   0259.000
      GOTO 200                                                          0260.000
C                                                                       0261.000
C     READ A EMAP  FILE                                                 0262.000
C                                                                       0263.000
   33 REWIND LIN                                                        0264.000
      WRITE(LP,3333)                                                    0265.000
 3333 FORMAT(/10X,' PLEASE WAIT. I AM READING THE EMAP  FILE')          0266.000
         READ (LIN) NDIFF,TITLE          , NCENT, LATT, NSYM,           0267.000
     1  ((TS(I,J),(IS(K,I,J),K=1,2),I=1,3),J=1,NSYM),                   0268.000
     2 MH, PARA, NPIC, NUMB, (IHKL(I),E(I), I=1,NUMB)                   0269.000
C     CHECK FOR NULL TAPE                                               0270.000
      IF (NDIFF.EQ.(-1)) CALL ERROR                                     0271.000
      IF(NDIFF.EQ.0) GO TO 145                                          0272.000
C     READ FILE PREPARED BY NORMAL FOR WEIGHTED OR DIFFERENCE FOURIER   0273.000
      NUMB=0                                                            0274.000
  115 READ(LIN) LH,LK,LL,EW,ED                                          0275.000
      DO 120 I=1,60                                                     0276.000
      IF(LH(I).LT.(-999)) GO TO 200                                     0277.000
      NUMB=NUMB+1                                                       0278.000
      ISC=1                                                             0279.000
      IF(LH(I).LT.0) ISC=-1                                             0280.000
      IHKL(NUMB)=ISC*(LH(I)*262144+LK(I)*512+LL(I))+256*(512+1)         0281.000
      MH(1)=MAX0(MH(1),IABS(LH(I)))                                     0282.000
      MH(2)=MAX0(MH(2),IABS(LK(I)))                                     0283.000
      MH(3)=MAX0(MH(3),IABS(LL(I)))                                     0284.000
      E(NUMB)=EW(I)                                                     0285.000
      EQ(NUMB)=FLOAT(ISC)*ED(I)                                         0286.000
      IF(NUMB.EQ.KUSER1) GO TO 130                                      0287.000
  120 CONTINUE                                                          0288.000
      GO TO 115                                                         0289.000
  130 WRITE(LT,140) KUSER1                                              0290.000
  140 FORMAT(1H0,25X,5HFIRST,I6,19HREFLEXIONS ACCEPTED)                 0291.000
      GO TO 200                                                         0292.000
C                                                                       0293.000
C     EXAMINE SETS OF PHASES IN ORDER OF COMBINED FIGURE OF MERIT UNLESS0294.000
C     SET IS SPECIFIED BY THE USER                                      0295.000
C                                                                       0296.000
  145 IF(NOSET.GT.0) GOTO 170                                           0297.000
  150 READ (LIN) NIX                                                    0298.000
      IF (NIX .GE. 0) GO TO 150                                         0299.000
      READ (LIN) NSET, (CFOM(J), J=1,NSET)                              0300.000
      NPOS = 1                                                          0301.000
      DO 160 J=2,NSET                                                   0302.000
      IF (CFOM(J) .GT. CFOM(NPOS)) NPOS = J                             0303.000
  160 CONTINUE                                                          0304.000
C     MARK SELECTED SET                                                 0305.000
      CFOM(NPOS) = CFOM(NPOS)-100.0                                     0306.000
      BACKSPACE LIN                                                     0307.000
      WRITE (LIN) NSET, (CFOM(J), J=1,NSET)                             0308.000
C     READ TO FIND SELECTED SET                                         0309.000
  170 REWIND LIN                                                        0310.000
      READ (LIN)                                                        0311.000
      NCS = 0                                                           0312.000
  180 READ (LIN) NO, ABSFOM, PSIZRO, RESID,RN4,RN5,(IPHAZ(I),I=1,NUMB)  0313.000
      NCS = NCS+1                                                       0314.000
      IF(NOSET.EQ.0.AND.NCS.NE.NPOS) GOTO 180                           0315.000
      IF(NOSET.NE.0.AND.NOSET.NE.NO) GOTO 180                           0316.000
      NOSET = NO                                                        0317.000
      DO 185 I=1,NUMB                                                   0318.000
      WT=FLOAT(MOD(IPHAZ(I),LK10))/100.0                                0319.000
      E(I)=E(I)*WT                                                      0320.000
  185 CONTINUE                                                          0321.000
  190 FORMAT (13X,6HNUMSET,7X,7HABS FOM,11X,8HPSI ZERO,11X,5HRESID,     0322.000
     1 11X,'NQEST',11X,'NQINT'/                                         0323.000
     1   9X,I8,F16.4,F18.2,F17.2,2F16.2)                                0324.000
      DO 192 I=1,NUMB                                                   0325.000
      PHASE=DTOR*FLOAT(IPHAZ(I)/LK10)                                   0326.000
      EQ(I)=E(I)*SIN(PHASE)                                             0327.000
      E(I)=E(I)*COS(PHASE)                                              0328.000
  192 CONTINUE                                                          0329.000
C                                                                       0330.000
C     SET DEFAULT VALUES OF PROGRAM PARAMETERS UNLESS SPECIFIED BY USER 0331.000
C                                                                       0332.000
C                                                                       0333.000
C     CALCULATE NUMBER OF GRID POINTS FOR FOURIER MAP                   0334.000
C                                                                       0335.000
  200 DO 280 I=1,3                                                      0336.000
      NP(I)=PARA(I)/FACTOR+0.5                                          0337.000
      IF(I.EQ.3) NP(I)=NP(I)+MOD(NP(I),2)                               0338.000
  250 NTEST=NP(I)                                                       0339.000
      DO 270 J=2,5                                                      0340.000
  260 IF(NTEST.NE.(NTEST/J)*J) GO TO 270                                0341.000
      NTEST=NTEST/J                                                     0342.000
      IF(NTEST.EQ.1) GO TO 280                                          0343.000
      GO TO 260                                                         0344.000
  270 CONTINUE                                                          0345.000
      NP(I)=NP(I)+1                                                     0346.000
      IF(I.EQ.3) NP(I)=NP(I)+1                                          0347.000
      GO TO 250                                                         0348.000
  280 CONTINUE                                                          0349.000
      NTAPIN=1                                                          0350.000
      GOTO(281,281,282,282,283,284,286),NMAP                            0351.000
  281 IF(RSCALE.LT.0.)RSCALE=100.                                       0352.000
      IF(FZERO.LT.0.)FZERO=0.                                           0353.000
      GOTO 287                                                          0354.000
  282 IF(RSCALE.LT.0.)RSCALE=10.                                        0355.000
      IF(FZERO.LT.0.)FZERO=-1.0-FZERO                                   0356.000
      GOTO 287                                                          0357.000
  283 NTAPIN=2                                                          0358.000
  284 NCENT=1                                                           0359.000
      DO 285 I=1,NSYM                                                   0360.000
      DO 285 J=1,3                                                      0361.000
  285 TS(J,I)=0.0                                                       0362.000
      FZERO=0.                                                          0363.000
      IF(NMAP.EQ.5)GOTO 287                                             0364.000
  286 NTAPIN=3                                                          0365.000
      FZERO=0.                                                          0366.000
  287 WRITE(LT,1013)COMPID,ITAPE(NTAPIN)                                0367.000
      I=0                                                               0368.000
      ASSIGN 288 TO NJUMP                                               0369.000
      LP=LT                                                             0370.000
  288 I=I+1                                                             0371.000
      GOTO(10,12,289,17,19,21,290,291),I                                0372.000
  289 IF(NMAP.NE.7)GOTO NJUMP                                           0373.000
      WRITE(LT,190)NO,ABSFOM,PSIZERO,RESID,RN4,RN5                      0374.000
      GOTO NJUMP                                                        0375.000
  290 IF(NTAPIN.NE.3)GOTO 24                                            0376.000
      GOTO NJUMP                                                        0377.000
  291 LP=LPP                                                            0378.000
      DO 292 I=1,3                                                      0379.000
  292 PC(I)=COS(PARA(I+3)*DTOR)                                         0380.000
      CVOL=PARA(1)*PARA(2)*PARA(3)*SQRT(1.0+2.0*PC(1)*PC(2)*PC(3)       0381.000
     1-PC(1)*PC(1)-PC(2)*PC(2)-PC(3)*PC(3))                             0382.000
      RSCALE=RSCALE/CVOL                                                0383.000
      FZERO=FZERO*RSCALE                                                0384.000
      IF(NPIC.EQ.0)NPIC=CVOL/(10.0*NSYM*NT*(NCENT+1))                   0385.000
      WRITE (LT,300) PARA,NP                                            0386.000
  300 FORMAT (// 30X,1HA,9X,1HB,9X,1HC,6X,5HALPHA,5X,4HBETA,4X,5HGAMMA/ 0387.000
     1   3X,20HUNIT CELL PARAMETERS,3F10.3,3F9.1//                      0388.000
     2                                     2X,21HNUMBER OF GRID POINTS, 0389.000
     3  6X,2HNX,8X,2HNY,8X,2HNZ/ 6X,15HIN FOURIER MAP ,3I10)            0390.000
      IF (NP(1) .LE. 200) GO TO 398                                     0391.000
      WRITE (LT,320)                                                    0392.000
  320 FORMAT (1H0,22X,76HNX GREATER THAN 200. INCREASE ARRAY INUT IN SUB0393.000
     1ROUTINE OUTPUT OF FFT PROGRAM)                                    0394.000
      CALL ERROR                                                        0395.000
C                                                                       0396.000
C     GENERATE ONE COMPLETE HEMISPHERE OF REFLEXIONS                    0397.000
C                                                                       0398.000
  398 DO 399 I=1,15                                                     0399.000
      TAB(I)=SIN(FLOAT(30*I)*DTOR)                                      0400.000
  399 CONTINUE                                                          0401.000
      WRITE(LP,397)ITAPE(NTAPIN)                                        0402.000
  397 FORMAT(/10X,' I AM NOW EXPANDING DATA READ FROM THE '             0403.000
     1,A5,' FILE ')                                                     0404.000
      NEX = 0                                                           0405.000
      IUN = 0                                                           0406.000
      FSUM=0.0                                                          0407.000
      FT1=2.0*FT1                                                       0408.000
      IP5=256*(262144+512+1)                                            0409.000
      NBATCH=0                                                          0410.000
      JN=0                                                              0411.000
  400 GOTO(402,401,410),NTAPIN                                          0412.000
  401 READ(LIN,END=801)JAA,JBB,FF,SIGI                                  0413.000
      FOBS=SQRT(FF)                                                     0414.000
      SIG=0.5*SIGI/FOBS                                                 0415.000
      GOTO 403                                                          0416.000
  402 READ(LIN,END=801)JAA,JBB,FF,SIGI,AC,BC,SS                         0417.000
      FOBS=SQRT(FF)                                                     0418.000
      SIG=0.5*SIGI/FOBS                                                 0419.000
      FC=SQRT(AC*AC+BC*BC)                                              0420.000
      DC=FOBS-FC                                                        0421.000
      SS=SQRT(SS)                                                       0422.000
      IF(SS.LT.FT2.OR.SS.GT.FT3)GOTO 402                                0423.000
      IF(ABS(FC/FOBS).LT.FT4)GOTO 402                                   0424.000
      A=AC/FC                                                           0425.000
      B=BC/FC                                                           0426.000
  403 IF(FOBS.LT.FT1*SIG)GOTO 400                                       0427.000
      NSET=JAA/16777216                                                 0428.000
      JAA=JAA-16777216*NSET                                             0429.000
      I1(1)=JAA/65536                                                   0430.000
      JAA=JAA-65536*I1(1)                                               0431.000
      I1(1)=I1(1)-128                                                   0432.000
      I1(2)=JAA/256                                                     0433.000
      I1(3)=JAA-256*I1(2)-128                                           0434.000
      I1(2)=I1(2)-128                                                   0435.000
      JN=JN+1                                                           0436.000
      GOTO(404,405,406,407,408),NMAP                                    0437.000
  404 FF=DC                                                             0438.000
      GOTO 409                                                          0439.000
  405 W=1.0                                                             0440.000
      FF=W*DF                                                           0441.000
      GOTO 409                                                          0442.000
  406 FF=FOBS                                                           0443.000
      GOTO 409                                                          0444.000
  407 FF=FC                                                             0445.000
      GOTO 409                                                          0446.000
  408 A=FOBS*FOBS                                                       0447.000
      B=0.0                                                             0448.000
      GOTO 450                                                          0449.000
  409 A=FF*A                                                            0450.000
      B=FF*B                                                            0451.000
      GOTO 450                                                          0452.000
  410 JN=JN+1                                                           0453.000
      IF(JN.LT.NUMB) GOTO 411                                           0454.000
      JN=NUMB                                                           0455.000
      GOTO 801                                                          0456.000
  411 A=E(JN)                                                           0457.000
      B=EQ(JN)                                                          0458.000
      I1(1)=IHKL(JN)/LK16                                               0459.000
      I1(2)=MOD(IHKL(JN),LK16)/LK8-LK7                                  0460.000
      I1(3)=MOD(IHKL(JN),LK8)-LK7                                       0461.000
      IF(NMAP.NE.6)GOTO 450                                             0462.000
      A=A*A-1.0                                                         0463.000
      B=0.0                                                             0464.000
  450 DO 700 J=1,NSYM                                                   0465.000
      DO 540 I=1,3                                                      0466.000
  540 I2(I)=0                                                           0467.000
      DO 560 I=1,3                                                      0468.000
      DO 550 K=1,2                                                      0469.000
      M=IABS(IS(K,I,J))                                                 0470.000
      IF(M.NE.0) I2(M)=I2(M)+I1(I)*ISIGN(1,IS(K,I,J))                   0471.000
  550 CONTINUE                                                          0472.000
  560 CONTINUE                                                          0473.000
      B1=1.0                                                            0474.000
      JS=1                                                              0475.000
      IF(I2(3)) 590,570,600                                             0476.000
  570 IF(I2(1)) 590,580,600                                             0477.000
  580 IF(I2(2)) 590,600,600                                             0478.000
  590 B1=-1.0                                                           0479.000
      JS=-1                                                             0480.000
  600 KIND(J)=IP5+JS*(262144*I2(1)+512*I2(2)+I2(3))                     0481.000
      DO 610 I=1,3                                                      0482.000
      IF(IABS(I2(I)).GE.NP(I)/2) GO TO 700                              0483.000
      IF(I2(I).GT.MH(I)) MH(I)=I2(I)                                    0484.000
  610 CONTINUE                                                          0485.000
      IF(J.EQ.1) GO TO 630                                              0486.000
      ND=J-1                                                            0487.000
      DO 620 I=1,ND                                                     0488.000
      IF(KIND(J).EQ.KIND(I)) GO TO 700                                  0489.000
  620 CONTINUE                                                          0490.000
  630 NBATCH=NBATCH+1                                                   0491.000
      NEX=NEX+1                                                         0492.000
      NU=0                                                              0493.000
      DO 640 I=1,3                                                      0494.000
  640 NU=NU-I1(I)*INT(TS(I,J)*12.0+0.1)                                 0495.000
      NU=MOD(NU,12)                                                     0496.000
      IF(NU.LE.0) NU=NU+12                                              0497.000
      XC=TAB(NU+3)                                                      0498.000
      XS=TAB(NU)                                                        0499.000
      GR=XC*A-XS*B                                                      0500.000
      GI=(XS*A+XC*B)*B1                                                 0501.000
      FSUM=FSUM+SQRT(GR*GR+GI*GI)                                       0502.000
      G(1,NBATCH)=GR                                                    0503.000
      G(2,NBATCH)=GI                                                    0504.000
      IH(NBATCH)=KIND(J)                                                0505.000
      IF(NBATCH.LT.100) GO TO 700                                       0506.000
      WRITE(LO) IH,G                                                    0507.000
      NBATCH=0                                                          0508.000
  700 CONTINUE                                                          0509.000
  800 GOTO 400                                                          0510.000
  801 IF(NBATCH.NE.0) WRITE(LO) IH,G                                    0511.000
      GOTO(900,900,900,900,802,802,804),NMAP                            0512.000
  802 RSCALE=500.0/FSUM                                                 0513.000
      GOTO 900                                                          0514.000
  804 RSCALE=3000.0/FSUM                                                0515.000
  900 WRITE(LT,901)MH,JN,NEX                                            0516.000
  901 FORMAT(//30X,1HH,9X,1HK,9X,1HL/ 5X,16HMAXIMUM INDICES ,3I10//     0517.000
     1I6,' REFLECTIONS HAVE BEEN EXPANDED TO',I6,' REFLECTIONS'//)      0518.000
      WRITE(LP,8888)                                                    0519.000
 8888 FORMAT(/10X,' I AM NOW CALCULATING THE MAP')                      0520.000
      RETURN                                                            0521.000
      END                                                               0522.000
C     ------------------------------------------------------------------0523.000
      SUBROUTINE PP1 (X, NSIZE, NX, NY, NZ)                             0524.000
      REAL X(NSIZE)                                                     0525.000
      INTEGER  P1, P2, R, SKIP, RECS,  D(5)                             0526.000
      CHARACTER COMPID(72)                                              0527.000
      CHARACTER*8 LABEL(200),FTYPE(13)                                  0528.000
      CHARACTER*4 ITLE(18),JUNK                                         0529.000
      COMMON/MODELC/COMPID,LABEL,FTYPE                                  0530.000
      COMMON /MAP/IMAP(8000)                                            0531.000
      EQUIVALENCE(COMPID(1),ITLE(1))                                    0532.000
      COMMON /PERIF/ NMAP, LMAP,LT,LO,LSCRA,TWOPI,LM,LPP                0533.000
      COMMON /TRAN/ RMAX,RMIN, RSCALE,FZERO,MH(3),LEVEL,IT4,XLIM(3)     0534.000
      COMMON/TAPE/NOSET,PARA(6),NSYM,NCENT,LATT,NPIC,TS(3,24),IS(2,3,24)0535.000
C                                                                       0536.000
      RMIN=0.                                                           0537.000
      RMAX=0.                                                           0538.000
      NSEC2=(NX+2)*NZ                                                   0539.000
      P1 = NSIZE/(2*NY*(MH(3) + 1))                                     0540.000
      P2 = NSIZE/(NX*NZ)                                                0541.000
      IF (P1 .LE. 0 .OR. P2 .LE. 0) GO TO 300                           0542.000
C                                                                       0543.000
C     PASS ONE. TRANSFORM ON K AND WRITE INTERMEDIATE RESULTS ON LSCRA  0544.000
C                                                                       0545.000
      REWIND LSCRA                                                      0546.000
      R = -MH(1)                                                        0547.000
  100 IF (R + P1 .GT. MH(1)) P1 = MH(1) + 1 - R                         0548.000
      CALL READ KL(X, NY, MH(3)+1, P1, R)                               0549.000
      D(1) = 2*(NY*P1*(MH(3)+1))                                        0550.000
      D(2) = 2                                                          0551.000
      D(3) = D(1)                                                       0552.000
      D(4) = D(1)                                                       0553.000
      D(5) = 2*NY                                                       0554.000
      CALL CMPLFT (X(1), X(2), NSIZE, NY, D)                            0555.000
      CALL WRITE Y (X, NY, MH(3)+1, P1, P2)                             0556.000
      R = R + P1                                                        0557.000
      IF (R .LE. MH(1)) GO TO 100                                       0558.000
C                                                                       0559.000
C     END OF PASS 1.  PASS 2 READS BACK THE INTERMEDIATE RESULTS,       0560.000
C     CALCULATES THE TRANSFORMS ON H AND L, AND WRITES OUT THE FINAL    0561.000
C     Y SECTIONS.                                                       0562.000
C                                                                       0563.000
      REWIND LMAP                                                       0564.000
      WRITE(LMAP)ITLE,NMAP,RSCALE,NOSET,PARA,NSYM,NCENT,LATT,NPIC,TS,IS 0565.000
      WRITE(LMAP) NX,NZ,NY                                              0566.000
C                                                                       0567.000
      REWIND LSCRA                                                      0568.000
      SKIP = 0                                                          0569.000
      R = 0                                                             0570.000
      P1 = NSIZE/(2*NY*(MH(3) + 1))                                     0571.000
      RECS = (NY - 1)/P2                                                0572.000
  200 IF (R + P2 .GT. NY) P2 = NY - R                                   0573.000
      CALL READ HL (X, NX, NZ/2, P2, MH(1), MH(3), P1, SKIP, RECS)      0574.000
      IF (R + P2 .LT. NY) REWIND LSCRA                                  0575.000
      SKIP = SKIP + 1                                                   0576.000
C                                                                       0577.000
C     X NOW CONTAINS INTERMEDIATE RESULTS STORED WITH L DOWN THE        0578.000
C     COLUMNS AND H ACROSS THE ROWS.                                    0579.000
C                                                                       0580.000
      D(1) = NX*NZ*P2                                                   0581.000
      D(2) = NZ                                                         0582.000
      D(3) = NZ*NX                                                      0583.000
      D(4) = 2*(MH(3) + 1)                                              0584.000
      D(5) = 2                                                          0585.000
      CALL CMPLFT (X(1), X(2), NSIZE, NX, D)                            0586.000
      D(2) = 2                                                          0587.000
      D(3) = D(1)                                                       0588.000
      D(4) = D(1)                                                       0589.000
      D(5) = NZ                                                         0590.000
      CALL HERMFT (X(1), X(2), NSIZE, NZ/2, D)                          0591.000
      CALL OUTPUT (X, NZ, NX, P2, R, NY,IMAP, NSEC2)                    0592.000
      R = R + P2                                                        0593.000
      IF (R .LT. NY) GO TO 200                                          0594.000
      WRITE(LMAP)RMIN,RMAX                                              0595.000
      WRITE(LT,201)RMIN,RMAX                                            0596.000
  201 FORMAT(/6X,' FUNCTION VALUES IN THE MAP RANGE FROM',F10.2,        0597.000
     1' TO ',F10.2)                                                     0598.000
C                                                                       0599.000
C     END OF JOB                                                        0600.000
C                                                                       0601.000
      RETURN                                                            0602.000
C                                                                       0603.000
  300 P1 = 2*NY*(MH(3)+1)                                               0604.000
      P2 = NX * NZ                                                      0605.000
      WRITE (LT,350) NSIZE, P1, P2                                      0606.000
  350 FORMAT(1H0,8X,37HREQUEST TOO LARGE. SPACE AVAILABLE IS,I10,       0607.000
     1 33H AND SPACE REQUESTED IS LARGER OF,I10,4H AND,I10)             0608.000
      CALL ERROR                                                        0609.000
C     THIS STATEMENT IS INACCESSIBLE AND IS FOR COMPILATION ONLY        0610.000
      RETURN                                                            0611.000
      END                                                               0612.000
C     ------------------------------------------------------------------0613.000
      SUBROUTINE READ KL(X, NY, NZ, NX, HS)                             0614.000
      INTEGER HS, H, HM, HL                                             0615.000
      DIMENSION X(2,NY,NZ,NX)                                           0616.000
      COMMON /REFL/ IH(100), G(2,100)                                   0617.000
      COMMON /PERIF/ NMAP, LMAP,LT,LO,LSCRA,TWOPI,LM,LPP                0618.000
      COMMON /CONST/ NEX                                                0619.000
C                                                                       0620.000
      REWIND LO                                                         0621.000
      NBLOCK = NEX/100+1                                                0622.000
      HM=262144*(HS+NX-1+256)+(512+1)*(2*256-1)                         0623.000
      HL=262144*(HS+256)                                                0624.000
      DO 100 H = 1, NX                                                  0625.000
      DO 100 L = 1, NZ                                                  0626.000
      DO 100 K = 1, NY                                                  0627.000
      X(1,K,L,H)=0.0                                                    0628.000
      X(2,K,L,H)=0.0                                                    0629.000
  100 CONTINUE                                                          0630.000
C                                                                       0631.000
      M=100                                                             0632.000
      DO 500 IB=1,NBLOCK                                                0633.000
      IF(IB.EQ.NBLOCK) M=MOD(NEX,100)                                   0634.000
      IF(M.EQ.0) GO TO 500                                              0635.000
      READ(LO) IH,G                                                     0636.000
      DO 300 J = 1, M                                                   0637.000
      IHJ=IH(J)                                                         0638.000
      IF(IHJ.GT.HM.OR.IHJ.LT.HL) GO TO 300                              0639.000
      H=IHJ/262144-256                                                  0640.000
      L=IHJ-262144*(H+256)                                              0641.000
      K=L/512-256                                                       0642.000
      L=L-512*(K+256)-256                                               0643.000
      NOKO = 0                                                          0644.000
      IF (H.EQ.0.AND.L.EQ.0.AND.K.NE.0) NOKO = NY-K+1                   0645.000
      H = H - HS + 1                                                    0646.000
      IF(K.LT.0)K=NY+K                                                  0647.000
      K = K + 1                                                         0648.000
      L = L + 1                                                         0649.000
      X(1,K,L,H)=G(1,J)                                                 0650.000
      X(2,K,L,H)=G(2,J)                                                 0651.000
      IF (NOKO.EQ.0) GOTO 600                                           0652.000
      X(1,NOKO,L,H)=X(1,K,L,H)                                          0653.000
      X(2,NOKO,L,H)=-X(2,K,L,H)                                         0654.000
  600 CONTINUE                                                          0655.000
  300 CONTINUE                                                          0656.000
  500 CONTINUE                                                          0657.000
      RETURN                                                            0658.000
      END                                                               0659.000
C     ------------------------------------------------------------------0660.000
      SUBROUTINE WRITE Y (X, NY, NZ, NX, SIZE)                          0661.000
      INTEGER SIZE, H, P, Q, R                                          0662.000
      DIMENSION X(2,NY,NZ,NX)                                           0663.000
      COMMON /PERIF/ NMAP, LMAP,LT,LO,LSCRA,TWOPI,LM,LPP                0664.000
C                                                                       0665.000
      P = SIZE                                                          0666.000
      Q = 0                                                             0667.000
  100 R = Q + 1                                                         0668.000
      IF (Q + P .GT. NY) P = NY - Q                                     0669.000
      Q = Q + P                                                         0670.000
      WRITE (LSCRA) (((X(1,K,L,H),X(2,K,L,H), H = 1, NX), K = R, Q), L =0671.000
     11,NZ)                                                             0672.000
      IF (Q .LT. NY) GO TO 100                                          0673.000
      RETURN                                                            0674.000
      END                                                               0675.000
C     ------------------------------------------------------------------0676.000
      SUBROUTINE READ HL (X, NX, NZ, NY, HMAX, LMAX, SIZE, SKIP, RECS)  0677.000
      INTEGER HMAX, SIZE, SKIP ,RECS, H, HL, HU, P, Q                   0678.000
      DIMENSION X(2,NZ,NX,NY)                                           0679.000
C                                                                       0680.000
C     READS IN INTERMEDIATE RESULTS FOR ALL H AND L FOR A BLOCK OF K.   0681.000
C     THE INPUT IS STORED WITH L DOWN THE COLUMNS AND H ON THE ROWS.    0682.000
      COMMON /PERIF/ NMAP, LMAP,LT,LO,LSCRA,TWOPI,LM,LPP                0683.000
C                                                                       0684.000
      LM = LMAX + 1                                                     0685.000
      P = SIZE                                                          0686.000
      HU = NX - HMAX                                                    0687.000
      IF (SKIP .LE. 0) GO TO 200                                        0688.000
      DO 100 Q = 1, SKIP                                                0689.000
      READ (LSCRA)                                                      0690.000
  100 CONTINUE                                                          0691.000
C                                                                       0692.000
C     READ DATA FOR NEGATIVE H                                          0693.000
C                                                                       0694.000
  200 IF (HU + P .GT. NX) GO TO 400                                     0695.000
      HL = HU + 1                                                       0696.000
      HU = HU + P                                                       0697.000
      READ (LSCRA) (((X(1,L,H,K),X(2,L,H,K), H = HL, HU), K = 1, NY), L 0698.000
     1=1,LM)                                                            0699.000
      IF (RECS .LE. 0) GO TO 200                                        0700.000
      DO 300 Q = 1, RECS                                                0701.000
      READ (LSCRA)                                                      0702.000
  300 CONTINUE                                                          0703.000
      GO TO 200                                                         0704.000
C                                                                       0705.000
C     GO PICK UP RECORD WHICH SPANS H = 0 IF NECESSARY                  0706.000
C                                                                       0707.000
  400 IF (HU .NE. NX) GO TO 700                                         0708.000
      HU = 0                                                            0709.000
C                                                                       0710.000
C     READ RECORDS FOR POSITIVE H                                       0711.000
C                                                                       0712.000
  500 IF (HU + P .GT. HMAX + 1) P = HMAX + 1 - HU                       0713.000
      HL = HU + 1                                                       0714.000
      HU = HU + P                                                       0715.000
      READ (LSCRA) (((X(1,L,H,K),X(2,L,H,K), H = HL, HU), K = 1, NY), L 0716.000
     1=1,LM)                                                            0717.000
  550 IF (HU .EQ. HMAX + 1) GO TO 800                                   0718.000
      IF (RECS .LE. 0) GO TO 500                                        0719.000
      DO 600 Q = 1, RECS                                                0720.000
      READ (LSCRA)                                                      0721.000
  600 CONTINUE                                                          0722.000
      GO TO 500                                                         0723.000
C                                                                       0724.000
C     SECTION TO READ A RECORD SPANNING H = 0                           0725.000
C                                                                       0726.000
  700 HL = HU + 1                                                       0727.000
      HU = HU + P - NX                                                  0728.000
      IF (HU .GT. HMAX + 1) HU = HMAX + 1                               0729.000
      READ (LSCRA) (((X(1,L,H,K),X(2,L,H,K), H = HL, NX), (X(1,L,H,K),X(0730.000
     12,L,H,K),H=1,HU),K=1,NY),L=1,LM)                                  0731.000
      GO TO 550                                                         0732.000
C                                                                       0733.000
C     PAD THE ARRAY WITH ZERO WHERE THERE ARE NO DATA, AND FILL IN THE  0734.000
C     -H K 0 RESULTS BY SYMMETRY.                                       0735.000
C                                                                       0736.000
  800 DO 900 H = 2, HU                                                  0737.000
      HL = NX + 2 - H                                                   0738.000
      DO 900 K = 1, NY                                                  0739.000
      X(1,1,HL,K)=X(1,1,H,K)                                            0740.000
      X(2,1,HL,K)=-X(2,1,H,K)                                           0741.000
  900 CONTINUE                                                          0742.000
      HL = HMAX + 2                                                     0743.000
      HU = NX - HMAX                                                    0744.000
      IF (HU .LT. HL) GO TO 920                                         0745.000
      DO 910 L = 1, LM                                                  0746.000
      DO 910 K = 1, NY                                                  0747.000
      DO 910 H = HL, HU                                                 0748.000
      X(1,L,H,K)=0.0                                                    0749.000
      X(2,L,H,K)=0.0                                                    0750.000
  910 CONTINUE                                                          0751.000
C                                                                       0752.000
  920 IF (LM .GE. NZ) GO TO 940                                         0753.000
      P = LM + 1                                                        0754.000
      DO 930 K = 1, NY                                                  0755.000
      DO 930 L = P, NZ                                                  0756.000
      DO 930 H = 1, NX                                                  0757.000
      X(1,L,H,K)=0.0                                                    0758.000
      X(2,L,H,K)=0.0                                                    0759.000
  930 CONTINUE                                                          0760.000
C                                                                       0761.000
C     END OF INPUT                                                      0762.000
C                                                                       0763.000
  940 RETURN                                                            0764.000
      END                                                               0765.000
C     ------------------------------------------------------------------0766.000
      SUBROUTINE OUTPUT (X, NZ, NX, NY, Y, NYT, IMAP,NSEC2)             0767.000
      CHARACTER COMPID(72)                                              0768.000
      CHARACTER*8 LABEL(200),FTYPE(13)                                  0769.000
      INTEGER IMAP(NSEC2)                                               0770.000
      REAL X(NZ,NX,NY)                                                  0771.000
      INTEGER STAR,BLANK,Y,SEC,XL,XU                                    0772.000
      COMMON /MODELC/COMPID,LABEL,FTYPE                                 0773.000
      COMMON /PERIF/ NMAP, LMAP,LT,LO,LSCRA,TWOPI,LM,LPP                0774.000
      COMMON /TRAN/ RMAX,RMIN, RSCALE,FZERO,MH(3),LEVEL,IT4,XLIM(3)     0775.000
      COMMON /TAPE/ INUT(200),LINE(50)                                  0776.000
      DATA STAR/2H**/,BLANK/4H    /,NCOL/36/                            0777.000
C                                                                       0778.000
C     PUT OUT SCALED Y-SECTIONS                                         0779.000
C                                                                       0780.000
      IF(IT4.EQ.0) GO TO 250                                            0781.000
      NXM=MIN0(NX,INT(XLIM(1)*FLOAT(NX))+1)                             0782.000
      NYM=MIN0(NYT,INT(XLIM(2)*FLOAT(NYT))+1)                           0783.000
      NZM=MIN0(NZ,INT(XLIM(3)*FLOAT(NZ))+1)                             0784.000
      ASC=1000.0/FLOAT(NX)                                              0785.000
C                                                                       0786.000
C                                                                       0787.000
  250 DO 900 K = 1, NY                                                  0788.000
      SEC = Y + K - 1                                                   0789.000
      DO 350 J=1,NZ                                                     0790.000
      KK=(J-1)*(NX+2)                                                   0791.000
      DO 330 I=1,NX                                                     0792.000
      XX=X(J,I,K)*RSCALE + FZERO                                        0793.000
      X(J,I,K)=XX                                                       0794.000
      RMIN=AMIN1(RMIN,XX)                                               0795.000
      RMAX=AMAX1(RMAX,XX)                                               0796.000
      KK=KK+1                                                           0797.000
      IMAP(KK)=0                                                        0798.000
      IF(XX.NE.0.0) IMAP(KK)=XX+SIGN(0.5,XX)                            0799.000
  330 CONTINUE                                                          0800.000
      IMAP(KK+1)=IMAP(KK+1-NX)                                          0801.000
      IMAP(KK+2)=IMAP(KK+2-NX)                                          0802.000
  350 CONTINUE                                                          0803.000
      KK=KK+2                                                           0804.000
      WRITE(LMAP)(IMAP(I),I=1,KK)                                       0805.000
      IF (IT4 .EQ. 0) GO TO 900                                         0806.000
      XU = 0                                                            0807.000
  450 XL = XU + 1                                                       0808.000
      XU = XU + NCOL                                                    0809.000
      IF(XU.GT.NXM)XU=NXM                                               0810.000
      DO 460 L=XL,XU                                                    0811.000
      INUT(L)=FLOAT((L-1))*ASC+0.5                                      0812.000
  460 CONTINUE                                                          0813.000
      NEC=1000.*(FLOAT(SEC))/(FLOAT(NYT))+0.5                           0814.000
      WRITE (IT4,480) COMPID,NEC                                        0815.000
  480 FORMAT(1H1,7X,72A1,8X,12HSECTION  Y =,I4/)                        0816.000
      WRITE (IT4,500) (INUT(NN),NN=XL,XU,2)                             0817.000
  500 FORMAT(8X,3HX =,I4,17I6)                                          0818.000
      ILIM = XL + 1                                                     0819.000
      WRITE (IT4,520) (INUT(NN),NN=ILIM,XU,2)                           0820.000
  520 FORMAT(12X,18I6)                                                  0821.000
      WRITE(IT4,530)                                                    0822.000
  530 FORMAT(1H0)                                                       0823.000
      DO 700 I = 1,NZM                                                  0824.000
      LINE(1)=FLOAT((I-1))*1000./FLOAT(NZ)+0.5                          0825.000
      L = 1                                                             0826.000
      DO 600 J = XL, XU                                                 0827.000
      L = L + 1                                                         0828.000
      INUT(L)=BLANK                                                     0829.000
      XX=X(I,J,K)                                                       0830.000
      LINE(L)=0                                                         0831.000
      IF(XX.NE.0.0) LINE(L)=INT(0.1*XX+SIGN(0.5,XX))                    0832.000
      IF (IABS(LINE(L)) .GT. 99) LINE(L) = ISIGN(99,LINE(L))            0833.000
      IF (10 * LINE(L) .GT. LEVEL) INUT(L) = STAR                       0834.000
  600 CONTINUE                                                          0835.000
      WRITE (IT4,630) (LINE(J), J = 1, L)                               0836.000
  630 FORMAT(4H Z =,I4,4H *  ,36I3)                                     0837.000
      WRITE(IT4,660) (INUT(J),J=2,L)                                    0838.000
  660 FORMAT(12X,36(1X,A2))                                             0839.000
  700 CONTINUE                                                          0840.000
      IF(XU.LT.NXM)GO TO 450                                            0841.000
      IF(SEC.GE.NYM) IT4=0                                              0842.000
  900 CONTINUE                                                          0843.000
      RETURN                                                            0844.000
      END                                                               0845.000
C     END OF PART 1 OF EXFFT                                            0846.000
C     ------------------------------------------------------------------0847.000
C     PART 2 OF EXFFT. THE TEN EYCK FFT ROUTINES.                       0848.000
C     INTERACTION WITH PART 1 IS VIA (1) SUBROUTINE CALLS, OR           0849.000
C     (2) COMMON /PERIF/.                                               0850.000
      SUBROUTINE CMPLFT (X, Y, NSIZE, N, D)                             0851.000
      REAL X(NSIZE), Y(NSIZE)                                           0852.000
      INTEGER D(5),PMAX,PSYM,TWOGRP,FACTOR(15),SYM(15),UNSYM(15)        0853.000
C                                                                       0854.000
C     COMPLEX FINITE DISCRETE FOURIER TRANSFORM                         0855.000
C     TRANSFORMS ONE DIMENSION OF MULTI-DIMENSIONAL DATA                0856.000
C     MODIFIED BY L. F. TEN EYCK FROM A ONE-DIMENSIONAL VERSION WRITTEN 0857.000
C     BY G. T. SANDE, 1969.                                             0858.000
C                                                                       0859.000
C     THIS PROGRAM CALCULATES THE TRANSFORM                             0860.000
C               (X(T) + I*Y(T))*(COS(2*PI*T/N) - I*SIN(2*PI*T/N))       0861.000
C                                                                       0862.000
      COMMON /PERIF/ NMAP, LMAP,LT,LO,LSCRA,TWOPI,LM,LPP                0863.000
      P MAX = 5                                                         0864.000
      TWO GRP = 4                                                       0865.000
C                                                                       0866.000
      CALL S R FP (N, P MAX, TWO GRP, FACTOR, SYM, P SYM, UN SYM)       0867.000
      CALL MDFTKD (N, FACTOR, D, X, Y, NSIZE)                           0868.000
      CALL DIPRP (N, SYM, P SYM, UN SYM, D, X, Y, NSIZE)                0869.000
      RETURN                                                            0870.000
      END                                                               0871.000
C     ------------------------------------------------------------------0872.000
      SUBROUTINE S R FP (PTS,PMAX,TWO GRP,FACTOR,SYM,P SYM,UN SYM)      0873.000
C     SYMMETRIZED REORDERING FACTORING PROGRAM                          0874.000
C                                                                       0875.000
      INTEGER PTS,PMAX,TWO GRP,P SYM, FACTOR (10), SYM (10), UN SYM (10)0876.000
      INTEGER PP(14), QQ (7), F,P,P TWO,Q,R                             0877.000
      COMMON /PERIF/ NMAP, LMAP,LT,LO,LSCRA,TWOPI,LM,LPP                0878.000
C                                                                       0879.000
      NEST=14                                                           0880.000
C                                                                       0881.000
      N=PTS                                                             0882.000
      P SYM=1                                                           0883.000
      F=2                                                               0884.000
      P=0                                                               0885.000
      Q=0                                                               0886.000
  100 IF (N.LE.1) GO TO 500                                             0887.000
      DO 200 J=F,PMAX                                                   0888.000
      IF (N.EQ.(N/J)*J) GO TO 300                                       0889.000
  200 CONTINUE                                                          0890.000
      CALL ERROR                                                        0891.000
  300 F=J                                                               0892.000
      N=N/F                                                             0893.000
      IF (N.EQ.(N/F)*F) GO TO 400                                       0894.000
      Q=Q+1                                                             0895.000
      QQ(Q)=F                                                           0896.000
      GO TO 100                                                         0897.000
  400 N=N/F                                                             0898.000
      P=P+1                                                             0899.000
      PP(P)=F                                                           0900.000
      P SYM=P SYM*F                                                     0901.000
      GO TO 100                                                         0902.000
C                                                                       0903.000
  500 R=1                                                               0904.000
      IF (Q.EQ.0) R=0                                                   0905.000
      IF (P.LT.1) GO TO 700                                             0906.000
      DO 600 J=1,P                                                      0907.000
      JJ=P+1-J                                                          0908.000
      SYM(J)=PP(JJ)                                                     0909.000
      FACTOR(J)=PP(JJ)                                                  0910.000
      JJ=P+Q+J                                                          0911.000
      FACTOR(JJ)=PP(J)                                                  0912.000
      JJ=P+R+J                                                          0913.000
      SYM(JJ)=PP(J)                                                     0914.000
  600 CONTINUE                                                          0915.000
  700 IF (Q.LT.1) GO TO 900                                             0916.000
      DO 800 J=1,Q                                                      0917.000
      JJ=P+J                                                            0918.000
      UN SYM(J)=QQ(J)                                                   0919.000
      FACTOR(JJ)=QQ(J)                                                  0920.000
  800 CONTINUE                                                          0921.000
      SYM(P+1)=PTS/P SYM**2                                             0922.000
  900 JJ=2*P+Q                                                          0923.000
      FACTOR(JJ+1)=0                                                    0924.000
      P TWO=1                                                           0925.000
      J=0                                                               0926.000
 1000 J=J+1                                                             0927.000
      IF (FACTOR(J).EQ.0) GO TO 1200                                    0928.000
      IF (FACTOR(J).NE.2) GO TO 1000                                    0929.000
      P TWO=P TWO*2                                                     0930.000
      FACTOR(J)=1                                                       0931.000
      IF (P TWO.GE.TWO GRP) GO TO 1100                                  0932.000
      IF (FACTOR(J+1).EQ.2) GO TO 1000                                  0933.000
 1100 FACTOR(J)=P TWO                                                   0934.000
      P TWO=1                                                           0935.000
      GO TO 1000                                                        0936.000
 1200 IF (P.EQ.0) R=0                                                   0937.000
      JJ=2*P+R                                                          0938.000
      SYM(JJ+1)=0                                                       0939.000
      IF (Q.LE.1) Q=0                                                   0940.000
      UN SYM(Q+1)=0                                                     0941.000
      RETURN                                                            0942.000
      END                                                               0943.000
C     ------------------------------------------------------------------0944.000
      SUBROUTINE DIPRP (PTS, SYM, P SYM, UN SYM, DIM, X, Y, NSIZE)      0945.000
C     DOUBLE IN PLACE REORDERING PROGRAM                                0946.000
C                                                                       0947.000
      REAL X(NSIZE), Y(NSIZE)                                           0948.000
      INTEGER SYM(10), UN SYM(10), DIM(5), PTS, PSYM, DK, P UN SYM, TEST0949.000
      LOGICAL ONE MOD                                                   0950.000
      INTEGER SEP, DELTA, P, P0, P1, P2, P3, P4, P5, SIZE               0951.000
      INTEGER V(14), MODULO(14), S(14), U(14)                           0952.000
      COMMON /PERIF/ NMAP, LMAP,LT,LO,LSCRA,TWOPI,LM,LPP                0953.000
C                                                                       0954.000
      NEST=14                                                           0955.000
C                                                                       0956.000
      NT = DIM(1)                                                       0957.000
      SEP = DIM(2)                                                      0958.000
      P2 = DIM(3)                                                       0959.000
      SIZE = DIM(4) - 1                                                 0960.000
      P4 = DIM(5)                                                       0961.000
      IF (SYM(1).EQ.0) GO TO 500                                        0962.000
      DO 100 J=1,NEST                                                   0963.000
      U(J)=1                                                            0964.000
      S(J)=1                                                            0965.000
  100 CONTINUE                                                          0966.000
      N=PTS                                                             0967.000
      DO 200 J=1,NEST                                                   0968.000
      IF (SYM(J).EQ.0) GO TO 300                                        0969.000
      JJ=NEST+1-J                                                       0970.000
      U(JJ)=N                                                           0971.000
      N=N/SYM(J)                                                        0972.000
      S(JJ)=N                                                           0973.000
  200 CONTINUE                                                          0974.000
  300 JJ=0                                                              0975.000
      L=1                                                               0976.000
      V(1)=1                                                            0977.000
  310 L=L+1                                                             0978.000
      V(L)=V(L-1)                                                       0979.000
  320 IF(L.LT.NEST) GO TO 310                                           0980.000
      N=V(NEST)                                                         0981.000
      JJ=JJ+1                                                           0982.000
      IF (JJ.GE.N) GO TO 400                                            0983.000
      DELTA = (N-JJ)*SEP                                                0984.000
      P1 = (JJ-1)*SEP + 1                                               0985.000
      DO 350 P0 = P1, NT, P2                                            0986.000
      P3 = P0 + SIZE                                                    0987.000
      DO 350 P = P0, P3, P4                                             0988.000
      P5 = P + DELTA                                                    0989.000
      T = X(P)                                                          0990.000
      X(P) = X(P5)                                                      0991.000
      X(P5) = T                                                         0992.000
      T = Y(P)                                                          0993.000
      Y(P) = Y(P5)                                                      0994.000
      Y(P5) = T                                                         0995.000
  350 CONTINUE                                                          0996.000
  400 V(L)=V(L)+S(L)                                                    0997.000
      IF(V(L).LE.U(L)) GO TO 320                                        0998.000
      L=L-1                                                             0999.000
      IF(L.NE.0) GO TO 400                                              1000.000
C                                                                       1001.000
  500 IF (UN SYM(1).EQ.0) GO TO 1900                                    1002.000
       P UN SYM=PTS/P SYM**2                                            1003.000
       MULT=P UN SYM/UN SYM(1)                                          1004.000
       TEST=(UN SYM(1)*UN SYM(2)-1)*MULT*P SYM                          1005.000
       LK=MULT                                                          1006.000
       DK=MULT                                                          1007.000
       DO 600 K=2,NEST                                                  1008.000
       IF (UN SYM(K).EQ.0) GO TO 700                                    1009.000
       LK=LK*UN SYM(K-1)                                                1010.000
       III=UN SYM(K)                                                    1011.000
       JJJ=DK                                                           1012.000
       IK=JJJ/III                                                       1013.000
       U(K)=(LK-IK)*P SYM                                               1014.000
       DK=IK                                                            1015.000
       MODS=K                                                           1016.000
  600 CONTINUE                                                          1017.000
  700 ONE MOD=MODS.LT.3                                                 1018.000
      IF (ONE MOD) GO TO 900                                            1019.000
      DO 800 J=3,MODS                                                   1020.000
      JJ=MODS+3-J                                                       1021.000
      MODULO(JJ)=U(J)                                                   1022.000
  800 CONTINUE                                                          1023.000
  900 MODULO(2)=U(2)                                                    1024.000
      JL=(P UN SYM-3)*P SYM                                             1025.000
      MS=P UN SYM*P SYM                                                 1026.000
C                                                                       1027.000
      DO 1800 J=P SYM,JL,P SYM                                          1028.000
      K=J                                                               1029.000
C                                                                       1030.000
 1000 K=K*MULT                                                          1031.000
      IF (ONE MOD) GO TO 1200                                           1032.000
      DO 1100 I=3,MODS                                                  1033.000
      K=K-(K/MODULO(I))*MODULO(I)                                       1034.000
 1100 CONTINUE                                                          1035.000
 1200 IF (K.GE.TEST) GO TO 1300                                         1036.000
      K=K-(K/MODULO(2))*MODULO(2)                                       1037.000
      GO TO 1400                                                        1038.000
 1300 K=K-(K/MODULO(2))*MODULO(2)+MODULO(2)                             1039.000
 1400 IF (K.LT.J) GO TO 1000                                            1040.000
C                                                                       1041.000
      IF (K.EQ.J) GO TO 1800                                            1042.000
      DELTA = (K-J)*SEP                                                 1043.000
      DO 1600 L=1,P SYM                                                 1044.000
      DO 1500 M=L,PTS,MS                                                1045.000
      P1 = (M+J-1)*SEP + 1                                              1046.000
      DO 1500 P0 = P1, NT, P2                                           1047.000
      P3 = P0 + SIZE                                                    1048.000
      DO 1500 JJ = P0, P3, P4                                           1049.000
      KK = JJ + DELTA                                                   1050.000
      T=X(JJ)                                                           1051.000
      X(JJ)=X(KK)                                                       1052.000
      X(KK)=T                                                           1053.000
      T=Y(JJ)                                                           1054.000
      Y(JJ)=Y(KK)                                                       1055.000
      Y(KK)=T                                                           1056.000
 1500 CONTINUE                                                          1057.000
 1600 CONTINUE                                                          1058.000
 1800 CONTINUE                                                          1059.000
 1900 RETURN                                                            1060.000
      END                                                               1061.000
C     ------------------------------------------------------------------1062.000
      SUBROUTINE MDFTKD (N, FACTOR, DIM, X, Y, NSIZE)                   1063.000
C     MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL DRIVER         1064.000
C                                                                       1065.000
      INTEGER FACTOR(10), DIM(5), F, P, R, S                            1066.000
      REAL X(NSIZE), Y(NSIZE)                                           1067.000
      COMMON /PERIF/ NMAP, LMAP,LT,LO,LSCRA,TWOPI,LM,LPP                1068.000
C                                                                       1069.000
      S = DIM(2)                                                        1070.000
      F = 0                                                             1071.000
      M = N                                                             1072.000
  100 F = F + 1                                                         1073.000
      P = FACTOR(F)                                                     1074.000
      IF (P.EQ.0) RETURN                                                1075.000
      M = M/P                                                           1076.000
      R = M*S                                                           1077.000
      GO TO (100, 200, 300, 400, 500), P                                1078.000
C                                                                       1079.000
  200 CALL R2 CFTK (N, M, X(1), Y(1), X(R+1), Y(R+1), DIM)              1080.000
      GO TO 100                                                         1081.000
C                                                                       1082.000
  300 CONTINUE                                                          1083.000
      CALL R3 CFTK (N, M, X(1), Y(1), X(R+1), Y(R+1), X(2*R+1), Y(2*R+1)1084.000
     ., DIM)                                                            1085.000
      GO TO 100                                                         1086.000
C                                                                       1087.000
  400 CALL R4 CFTK (N, M, X(1), Y(1), X(R+1), Y(R+1), X(2*R+1), Y(2*R+1)1088.000
     ., X(3*R+1), Y(3*R+1), DIM)                                        1089.000
      GO TO 100                                                         1090.000
C                                                                       1091.000
  500 CALL R5 CFTK (N, M, X(1), Y(1), X(R+1), Y(R+1), X(2*R+1), Y(2*R+1)1092.000
     ., X(3*R+1), Y(3*R+1), X(4*R+1), Y(4*R+1), DIM)                    1093.000
      GO TO 100                                                         1094.000
C                                                                       1095.000
      END                                                               1096.000
C     ------------------------------------------------------------------1097.000
      SUBROUTINE R2 CFTK (N, M, X0, Y0, X1, Y1, DIM)                    1098.000
C     RADIX 2 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL        1099.000
C                                                                       1100.000
      INTEGER DIM(5), SIZE, SEP                                         1101.000
      REAL X0(10), Y0(10), X1(10), Y1(10), IS, IU                       1102.000
      LOGICAL FOLD,ZERO                                                 1103.000
      COMMON /PERIF/ NMAP, LMAP,LT,LO,LSCRA,TWOPI,LM,LPP                1104.000
C                                                                       1105.000
      NT = DIM(1)                                                       1106.000
      SEP = DIM(2)                                                      1107.000
      L1 = DIM(3)                                                       1108.000
      SIZE = DIM(4) - 1                                                 1109.000
      K2 = DIM(5)                                                       1110.000
      NS = N*SEP                                                        1111.000
      M2=M*2                                                            1112.000
      FM2 = FLOAT(M2)                                                   1113.000
      M OVER 2=M/2+1                                                    1114.000
      MM2 = SEP*M2                                                      1115.000
C                                                                       1116.000
      FJM1 = -1.0                                                       1117.000
      DO 600 J=1,M OVER 2                                               1118.000
      FOLD=J.GT.1 .AND. 2*J.LT.M+2                                      1119.000
      K0 = (J-1)*SEP + 1                                                1120.000
      FJM1 = FJM1 + 1.0                                                 1121.000
      ANGLE = TWO PI*FJM1/FM2                                           1122.000
      ZERO=ANGLE.EQ.0.0                                                 1123.000
      IF (ZERO) GO TO 200                                               1124.000
      C=COS(ANGLE)                                                      1125.000
      S=SIN(ANGLE)                                                      1126.000
      GO TO 200                                                         1127.000
  100 FOLD=.FALSE.                                                      1128.000
      K0 = (M+1-J)*SEP + 1                                              1129.000
      C=-C                                                              1130.000
  200 DO 500 KK = K0, NS, MM2                                           1131.000
      DO 440 L = KK, NT, L1                                             1132.000
      K1 = L + SIZE                                                     1133.000
      DO 420 K = L, K1, K2                                              1134.000
      RS=X0(K)+X1(K)                                                    1135.000
      IS=Y0(K)+Y1(K)                                                    1136.000
      RU=X0(K)-X1(K)                                                    1137.000
      IU=Y0(K)-Y1(K)                                                    1138.000
      X0(K)=RS                                                          1139.000
      Y0(K)=IS                                                          1140.000
      IF (ZERO) GO TO 300                                               1141.000
      X1(K)=RU*C+IU*S                                                   1142.000
      Y1(K)=IU*C-RU*S                                                   1143.000
      GO TO 420                                                         1144.000
  300 X1(K)=RU                                                          1145.000
      Y1(K)=IU                                                          1146.000
  420 CONTINUE                                                          1147.000
  440 CONTINUE                                                          1148.000
  500 CONTINUE                                                          1149.000
      IF (FOLD) GO TO 100                                               1150.000
  600 CONTINUE                                                          1151.000
      RETURN                                                            1152.000
      END                                                               1153.000
C     ------------------------------------------------------------------1154.000
      SUBROUTINE R3 CFTK (N, M, X0, Y0, X1, Y1, X2, Y2, DIM)            1155.000
C     RADIX 3 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL        1156.000
C                                                                       1157.000
      INTEGER DIM(5), SIZE, SEP                                         1158.000
      REAL X0(10),Y0(10),X1(10),Y1(10),X2(10),Y2(10), I0,I1,I2,IA,IB,IS 1159.000
      LOGICAL FOLD,ZERO                                                 1160.000
      COMMON /PERIF/ NMAP, LMAP,LT,LO,LSCRA,TWOPI,LM,LPP                1161.000
      DATA A/-0.5/, B/0.86602540/                                       1162.000
C                                                                       1163.000
      NT = DIM(1)                                                       1164.000
      SEP = DIM(2)                                                      1165.000
      L1 = DIM(3)                                                       1166.000
      SIZE = DIM(4) - 1                                                 1167.000
      K2 = DIM(5)                                                       1168.000
      NS = N*SEP                                                        1169.000
      M3=M*3                                                            1170.000
      FM3 = FLOAT(M3)                                                   1171.000
      MM3 = SEP*M3                                                      1172.000
      M OVER 2=M/2+1                                                    1173.000
C                                                                       1174.000
      FJM1 = -1.0                                                       1175.000
      DO 600 J=1,M OVER 2                                               1176.000
      FOLD=J.GT.1 .AND. 2*J.LT.M+2                                      1177.000
      K0 = (J-1)*SEP + 1                                                1178.000
      FJM1 = FJM1 + 1.0                                                 1179.000
      ANGLE = TWO PI*FJM1/FM3                                           1180.000
      ZERO=ANGLE.EQ.0.0                                                 1181.000
      IF (ZERO) GO TO 200                                               1182.000
      C1=COS(ANGLE)                                                     1183.000
      S1=SIN(ANGLE)                                                     1184.000
      C2=C1*C1-S1*S1                                                    1185.000
      S2=S1*C1+C1*S1                                                    1186.000
      GO TO 200                                                         1187.000
  100 FOLD=.FALSE.                                                      1188.000
      K0 = (M+1-J)*SEP + 1                                              1189.000
      T=C1*A+S1*B                                                       1190.000
      S1=C1*B-S1*A                                                      1191.000
      C1=T                                                              1192.000
      T=C2*A-S2*B                                                       1193.000
      S2=-C2*B-S2*A                                                     1194.000
      C2=T                                                              1195.000
  200 DO 500 KK = K0, NS, MM3                                           1196.000
      DO 440 L = KK, NT, L1                                             1197.000
      K1 = L + SIZE                                                     1198.000
      DO 420 K = L, K1, K2                                              1199.000
      R0=X0(K)                                                          1200.000
      I0=Y0(K)                                                          1201.000
      RS=X1(K)+X2(K)                                                    1202.000
      IS=Y1(K)+Y2(K)                                                    1203.000
      X0(K)=R0+RS                                                       1204.000
      Y0(K)=I0+IS                                                       1205.000
      RA=R0+RS*A                                                        1206.000
      IA=I0+IS*A                                                        1207.000
      RB=(X1(K)-X2(K))*B                                                1208.000
      IB=(Y1(K)-Y2(K))*B                                                1209.000
      IF (ZERO) GO TO 300                                               1210.000
      R1=RA+IB                                                          1211.000
      I1=IA-RB                                                          1212.000
      R2=RA-IB                                                          1213.000
      I2=IA+RB                                                          1214.000
      X1(K)=R1*C1+I1*S1                                                 1215.000
      Y1(K)=I1*C1-R1*S1                                                 1216.000
      X2(K)=R2*C2+I2*S2                                                 1217.000
      Y2(K)=I2*C2-R2*S2                                                 1218.000
      GO TO 420                                                         1219.000
  300 X1(K)=RA+IB                                                       1220.000
      Y1(K)=IA-RB                                                       1221.000
      X2(K)=RA-IB                                                       1222.000
      Y2(K)=IA+RB                                                       1223.000
  420 CONTINUE                                                          1224.000
  440 CONTINUE                                                          1225.000
  500 CONTINUE                                                          1226.000
      IF (FOLD) GO TO 100                                               1227.000
  600 CONTINUE                                                          1228.000
      RETURN                                                            1229.000
      END                                                               1230.000
C     ------------------------------------------------------------------1231.000
      SUBROUTINE R4 CFTK (N, M, X0, Y0, X1, Y1, X2, Y2, X3, Y3, DIM)    1232.000
C     RADIX 4 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL        1233.000
C                                                                       1234.000
      INTEGER DIM(5), SIZE, SEP                                         1235.000
      REAL X0(10),Y0(10),X1(10),Y1(10),X2(10),Y2(10),X3(10),Y3(10)      1236.000
      LOGICAL FOLD,ZERO                                                 1237.000
      REAL I1,I2,I3,IS0,IS1,IU0,IU1                                     1238.000
      COMMON /PERIF/ NMAP, LMAP,LT,LO,LSCRA,TWOPI,LM,LPP                1239.000
C                                                                       1240.000
      NT = DIM(1)                                                       1241.000
      SEP = DIM(2)                                                      1242.000
      L1 = DIM(3)                                                       1243.000
      SIZE = DIM(4) - 1                                                 1244.000
      K2 = DIM(5)                                                       1245.000
      NS = N*SEP                                                        1246.000
      M4=M*4                                                            1247.000
      FM4 = FLOAT(M4)                                                   1248.000
      MM4 = SEP*M4                                                      1249.000
      M OVER 2=M/2+1                                                    1250.000
C                                                                       1251.000
      FJM1 = -1.0                                                       1252.000
      DO 600 J=1,M OVER 2                                               1253.000
      FOLD=J.GT.1 .AND. 2*J.LT.M+2                                      1254.000
      K0 = (J-1)*SEP + 1                                                1255.000
      FJM1 = FJM1 + 1.0                                                 1256.000
      ANGLE = TWO PI*FJM1/FM4                                           1257.000
      ZERO=ANGLE.EQ.0.0                                                 1258.000
      IF (ZERO) GO TO 200                                               1259.000
      C1=COS(ANGLE)                                                     1260.000
      S1=SIN(ANGLE)                                                     1261.000
      C2=C1*C1-S1*S1                                                    1262.000
      S2=S1*C1+C1*S1                                                    1263.000
      C3=C2*C1-S2*S1                                                    1264.000
      S3=S2*C1+C2*S1                                                    1265.000
      GO TO 200                                                         1266.000
  100 FOLD=.FALSE.                                                      1267.000
      K0 = (M+1-J)*SEP + 1                                              1268.000
      T=C1                                                              1269.000
      C1=S1                                                             1270.000
      S1=T                                                              1271.000
      C2=-C2                                                            1272.000
      T=C3                                                              1273.000
      C3=-S3                                                            1274.000
      S3=-T                                                             1275.000
  200 DO 500 KK = K0, NS, MM4                                           1276.000
      DO 440 L = KK, NT, L1                                             1277.000
      K1 = L + SIZE                                                     1278.000
      DO 420 K = L, K1, K2                                              1279.000
      RS0=X0(K)+X2(K)                                                   1280.000
      IS0=Y0(K)+Y2(K)                                                   1281.000
      RU0=X0(K)-X2(K)                                                   1282.000
      IU0=Y0(K)-Y2(K)                                                   1283.000
      RS1=X1(K)+X3(K)                                                   1284.000
      IS1=Y1(K)+Y3(K)                                                   1285.000
      RU1=X1(K)-X3(K)                                                   1286.000
      IU1=Y1(K)-Y3(K)                                                   1287.000
      X0(K)=RS0+RS1                                                     1288.000
      Y0(K)=IS0+IS1                                                     1289.000
      IF (ZERO) GO TO 300                                               1290.000
      R1=RU0+IU1                                                        1291.000
      I1=IU0-RU1                                                        1292.000
      R2=RS0-RS1                                                        1293.000
      I2=IS0-IS1                                                        1294.000
      R3=RU0-IU1                                                        1295.000
      I3=IU0+RU1                                                        1296.000
      X2(K)=R1*C1+I1*S1                                                 1297.000
      Y2(K)=I1*C1-R1*S1                                                 1298.000
      X1(K)=R2*C2+I2*S2                                                 1299.000
      Y1(K)=I2*C2-R2*S2                                                 1300.000
      X3(K)=R3*C3+I3*S3                                                 1301.000
      Y3(K)=I3*C3-R3*S3                                                 1302.000
      GO TO 420                                                         1303.000
  300 X2(K)=RU0+IU1                                                     1304.000
      Y2(K)=IU0-RU1                                                     1305.000
      X1(K)=RS0-RS1                                                     1306.000
      Y1(K)=IS0-IS1                                                     1307.000
      X3(K)=RU0-IU1                                                     1308.000
      Y3(K)=IU0+RU1                                                     1309.000
  420 CONTINUE                                                          1310.000
  440 CONTINUE                                                          1311.000
  500 CONTINUE                                                          1312.000
      IF (FOLD) GO TO 100                                               1313.000
  600 CONTINUE                                                          1314.000
      RETURN                                                            1315.000
      END                                                               1316.000
C     ------------------------------------------------------------------1317.000
      SUBROUTINE R5 CFTK(N,M,X0, Y0, X1, Y1, X2, Y2, X3, Y3, X4, Y4,DIM)1318.000
C     RADIX 5 MULTI-DIMENSIONAL COMPLEX FOURIER TRANSFORM KERNEL        1319.000
C                                                                       1320.000
      INTEGER DIM(5), SIZE, SEP                                         1321.000
      REAL X0(10),Y0(10),X1(10),Y1(10),X2(10),Y2(10),X3(10),Y3(10),     1322.000
     1 X4(10),Y4(10),I0,I1,I2,I3,I4,IA1,IA2,IB1,IB2,IS1,IS2,IU1,IU2     1323.000
      LOGICAL FOLD,ZERO                                                 1324.000
      COMMON /PERIF/ NMAP, LMAP,LT,LO,LSCRA,TWOPI,LM,LPP                1325.000
      DATA A1/0.30901699/,B1/0.95105652/,A2/-0.80901699/,B2/0.58778525/ 1326.000
C                                                                       1327.000
      NT = DIM(1)                                                       1328.000
      SEP = DIM(2)                                                      1329.000
      L1 = DIM(3)                                                       1330.000
      SIZE = DIM(4) - 1                                                 1331.000
      K2 = DIM(5)                                                       1332.000
      NS = N*SEP                                                        1333.000
      M5=M*5                                                            1334.000
      FM5 = FLOAT(M5)                                                   1335.000
      MM5 = SEP*M5                                                      1336.000
      M OVER 2=M/2+1                                                    1337.000
C                                                                       1338.000
      FJM1 = -1.0                                                       1339.000
      DO 600 J=1,M OVER 2                                               1340.000
      FOLD=J.GT.1 .AND. 2*J.LT.M+2                                      1341.000
      K0 = (J-1)*SEP + 1                                                1342.000
      FJM1 = FJM1 + 1.0                                                 1343.000
      ANGLE = TWO PI*FJM1/FM5                                           1344.000
      ZERO=ANGLE.EQ.0.0                                                 1345.000
      IF (ZERO) GO TO 200                                               1346.000
      C1=COS(ANGLE)                                                     1347.000
      S1=SIN(ANGLE)                                                     1348.000
      C2=C1*C1-S1*S1                                                    1349.000
      S2=S1*C1+C1*S1                                                    1350.000
      C3=C2*C1-S2*S1                                                    1351.000
      S3=S2*C1+C2*S1                                                    1352.000
      C4=C2*C2-S2*S2                                                    1353.000
      S4=S2*C2+C2*S2                                                    1354.000
      GO TO 200                                                         1355.000
  100 FOLD=.FALSE.                                                      1356.000
      K0 = (M+1-J)*SEP + 1                                              1357.000
      T=C1*A1+S1*B1                                                     1358.000
      S1=C1*B1-S1*A1                                                    1359.000
      C1=T                                                              1360.000
      T=C2*A2+S2*B2                                                     1361.000
      S2=C2*B2-S2*A2                                                    1362.000
      C2=T                                                              1363.000
      T=C3*A2-S3*B2                                                     1364.000
      S3=-C3*B2-S3*A2                                                   1365.000
      C3=T                                                              1366.000
      T=C4*A1-S4*B1                                                     1367.000
      S4=-C4*B1-S4*A1                                                   1368.000
      C4=T                                                              1369.000
  200 DO 500 KK = K0, NS, MM5                                           1370.000
      DO 440 L = KK, NT, L1                                             1371.000
      K1 = L + SIZE                                                     1372.000
      DO 420 K = L, K1, K2                                              1373.000
      R0=X0(K)                                                          1374.000
      I0=Y0(K)                                                          1375.000
      RS1=X1(K)+X4(K)                                                   1376.000
      IS1=Y1(K)+Y4(K)                                                   1377.000
      RU1=X1(K)-X4(K)                                                   1378.000
      IU1=Y1(K)-Y4(K)                                                   1379.000
      RS2=X2(K)+X3(K)                                                   1380.000
      IS2=Y2(K)+Y3(K)                                                   1381.000
      RU2=X2(K)-X3(K)                                                   1382.000
      IU2=Y2(K)-Y3(K)                                                   1383.000
      X0(K)=R0+RS1+RS2                                                  1384.000
      Y0(K)=I0+IS1+IS2                                                  1385.000
      RA1=R0+RS1*A1+RS2*A2                                              1386.000
      IA1=I0+IS1*A1+IS2*A2                                              1387.000
      RA2=R0+RS1*A2+RS2*A1                                              1388.000
      IA2=I0+IS1*A2+IS2*A1                                              1389.000
      RB1=RU1*B1+RU2*B2                                                 1390.000
      IB1=IU1*B1+IU2*B2                                                 1391.000
      RB2=RU1*B2-RU2*B1                                                 1392.000
      IB2=IU1*B2-IU2*B1                                                 1393.000
      IF (ZERO) GO TO 300                                               1394.000
      R1=RA1+IB1                                                        1395.000
      I1=IA1-RB1                                                        1396.000
      R2=RA2+IB2                                                        1397.000
      I2=IA2-RB2                                                        1398.000
      R3=RA2-IB2                                                        1399.000
      I3=IA2+RB2                                                        1400.000
      R4=RA1-IB1                                                        1401.000
      I4=IA1+RB1                                                        1402.000
      X1(K)=R1*C1+I1*S1                                                 1403.000
      Y1(K)=I1*C1-R1*S1                                                 1404.000
      X2(K)=R2*C2+I2*S2                                                 1405.000
      Y2(K)=I2*C2-R2*S2                                                 1406.000
      X3(K)=R3*C3+I3*S3                                                 1407.000
      Y3(K)=I3*C3-R3*S3                                                 1408.000
      X4(K)=R4*C4+I4*S4                                                 1409.000
      Y4(K)=I4*C4-R4*S4                                                 1410.000
      GO TO 420                                                         1411.000
  300 X1(K)=RA1+IB1                                                     1412.000
      Y1(K)=IA1-RB1                                                     1413.000
      X2(K)=RA2+IB2                                                     1414.000
      Y2(K)=IA2-RB2                                                     1415.000
      X3(K)=RA2-IB2                                                     1416.000
      Y3(K)=IA2+RB2                                                     1417.000
      X4(K)=RA1-IB1                                                     1418.000
      Y4(K)=IA1+RB1                                                     1419.000
  420 CONTINUE                                                          1420.000
  440 CONTINUE                                                          1421.000
  500 CONTINUE                                                          1422.000
      IF (FOLD) GO TO 100                                               1423.000
  600 CONTINUE                                                          1424.000
      RETURN                                                            1425.000
      END                                                               1426.000
C     ------------------------------------------------------------------1427.000
      SUBROUTINE HERMFT (X, Y, NSIZE, N, DIM)                           1428.000
      REAL X(NSIZE), Y(NSIZE)                                           1429.000
      INTEGER DIM(5), D2, D3, D4, D5                                    1430.000
C                                                                       1431.000
C     HERMITIAN SYMMETRIC FOURIER TRANSFORM                             1432.000
C                                                                       1433.000
      COMMON /PERIF/ NMAP, LMAP,LT,LO,LSCRA,TWOPI,LM,LPP                1434.000
C                                                                       1435.000
      TWO N = FLOAT(2*N)                                                1436.000
      NT = DIM(1)                                                       1437.000
      D2 = DIM(2)                                                       1438.000
      D3 = DIM(3)                                                       1439.000
      D4 = DIM(4) - 1                                                   1440.000
      D5 = DIM(5)                                                       1441.000
C                                                                       1442.000
      DO 100 I0 = 1, NT, D3                                             1443.000
      I1 = I0 + D4                                                      1444.000
      DO 100 I = I0, I1, D5                                             1445.000
      A = X(I)                                                          1446.000
      B = Y(I)                                                          1447.000
      X(I) = A + B                                                      1448.000
      Y(I) = A - B                                                      1449.000
  100 CONTINUE                                                          1450.000
C                                                                       1451.000
      N OVER 2 = N/2 + 1                                                1452.000
      IF (N OVER 2 .LT. 2) GO TO 500                                    1453.000
      DO 400 I0 = 2, N OVER 2                                           1454.000
      ANGLE = TWO PI*FLOAT(I0-1)/TWO N                                  1455.000
      CO = COS(ANGLE)                                                   1456.000
      SI = SIN(ANGLE)                                                   1457.000
      K = (N + 2 - 2*I0)*D2                                             1458.000
      K1 = (I0 - 1)*D2 + 1                                              1459.000
      DO 300 I1 = K1, NT, D3                                            1460.000
      I2 = I1 + D4                                                      1461.000
      DO 200 I = I1, I2, D5                                             1462.000
      J = I + K                                                         1463.000
      A = X(I) + X(J)                                                   1464.000
      B = X(I) - X(J)                                                   1465.000
      C = Y(I) + Y(J)                                                   1466.000
      D = Y(I) - Y(J)                                                   1467.000
      E = B*CO + C*SI                                                   1468.000
      F = B*SI - C*CO                                                   1469.000
      X(I) = A + F                                                      1470.000
      X(J) = A - F                                                      1471.000
      Y(I) = E + D                                                      1472.000
      Y(J) = E - D                                                      1473.000
  200 CONTINUE                                                          1474.000
  300 CONTINUE                                                          1475.000
  400 CONTINUE                                                          1476.000
C                                                                       1477.000
      CALL CMPLFT (X, Y, NSIZE, N, DIM)                                 1478.000
C                                                                       1479.000
  500 RETURN                                                            1480.000
      END                                                               1481.000
C     FINISH OF PART 2 OF EXFFT                                         1482.000
