C @(#)sfft.for	19.1 (ES0-DMD) 02/25/03 13:30:48
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  SUBROUTINE FFT           VERSION 1.0               JAN. 30, 1983
C  IEEE TRANS AU 17
C 
C FAST FOURIER ANALYSIS SUBROUTINE (FFA)
C THIS PROGRAM IS TAKEN FROM IEEE TRANS. AU 17 (1969) 138 
C THIS SUBROUTINE REPLACES THE REAL VECTOR B(K), (K=1,2,...,N)
C WITH ITS FINITE DISCRETE FOURIER TRANSFORM X(J) WHICH IS RETURNED 
C SUCH THAT X(0) = B(1), X(N/2) = B(2), AND X(J) = B(2*J+1)+I*B(2*J+2), 
C (I.E. X(1) = B(3) + I B(4), X(2) = B(5) + I B(6), ETC.) THE 
C SUBROUTINE IS CALLED AS FFA(M,B) WHERE N = 2**M AND B IS AN 
C N TERM REAL ARRAY. A REAL-VALUED, RADIX 8 ALGORITHM IS USED 
C WITH IN-PLACE REORDERING AND THE TRIG FUNCTIONS ARE COMPUTED AS 
C NEEDED. 
C
C     SINGLE PRECISION VERSION
C 
C---------------------------------------------------------------------
      SUBROUTINE FFA(M,B)
      IMPLICIT NONE
      INTEGER  M
      REAL     B(1)
C
      INTEGER  N, NN
      INTEGER  N8POW
      INTEGER  INT, IT
C
      N = 2**M
      N8POW = M/3 
C 
C DO A RADIX 2 OR RADIX 4 ITERATION FIRST IF ONE IS REQUIRED
C 
      IF (M-N8POW*3-1) 30,20,10
   10 NN = 4
      INT = N/NN
      CALL R4TR(INT,B(1),B(INT+1),B(2*INT+1),B(3*INT+1))
      GO TO 40
   20 NN = 2
      INT = N/NN
      CALL R2TR(INT,B(1),B(INT+1))
      GO TO 40
   30 NN = 1
C 
C PERFORM RADIX 8 ITERATIONS
C 
   40 IF (N8POW) 70,70,45
   45 DO 50 IT = 1,N8POW
      NN = NN*8 
      INT = N/NN
      CALL R8TR(INT,NN,B(1),B(INT+1),B(INT*2+1),B(3*INT+1), 
     1 B(4*INT+1),B(5*INT+1),B(6*INT+1),B(7*INT+1),B(1),B(INT+1), 
     2 B(2*INT+1),B(3*INT+1),B(4*INT+1),B(5*INT+1),B(6*INT+1),
     3 B(7*INT+1))
   50 CONTINUE
C 
C PERFORM IN-PLACE REORDERING 
C 
      CALL ORD1(M,B)
      CALL ORD2(M,B)
   70 RETURN
      END
C--------------------------------------------------------------------
      SUBROUTINE FFS(M,B) 
C 
C FAST FOURIER SYNTHESIS SUBROUTINE (FFS) 
C THIS SUBROUTINE SYNTHESIZES THE REAL VECTOR B(K), WHERE 
C K = 1,2,...N. THE INITIAL FOURIER COEFFICIENTS ARE PLACED IN
C THE B ARRAY WITH THE DC AND HIGHEST HARMONIC COEFFICIENTS 
C STORED IN B(1) AND B(2) RESPECTIVELY AND ALL OTHERS LOCATED 
C SUCH THAT THE JTH HARMONIC IS STORED AS B(2*J+1) + I*B(2*J+2).
C B IS THE N TERM REAL ARRAY DISCUSSED ABOVE. 
C------------------------------------------------------------------- 
      IMPLICIT  NONE
      INTEGER   M
      REAL      B(1)
C
      INTEGER   N, NN
      INTEGER   N8POW
      INTEGER   INT, IT
C
      N = 2**M
      N8POW = M/3 
      IF (N8POW) 15,15,5 
C 
C REORDER THE INPUT FOURIER COEFFICIENTS
C 
    5 CALL ORD2(M,B)
      CALL ORD1(M,B)
C 
C PERFORM THE RADIX 8 ITERATIONS
C 
      NN = 8*N
      DO 10 IT=1,N8POW
      NN = NN/8 
      INT = N/NN
      CALL R8SYN(INT,NN,B(1),B(INT+1),B(2*INT+1),B(3*INT+1),
     1 B(4*INT+1),B(5*INT+1),B(6*INT+1),B(7*INT+1),B(1),
     2 B(INT+1),B(2*INT+1),B(3*INT+1),B(4*INT+1),B(5*INT+1),
     3 B(6*INT+1),B(7*INT+1)) 
   10 CONTINUE
C 
C DO A RADIX 2 OR RADIX 4 ITERATION IF ONE IS REQUIRED
C 
   15 IF (M-N8POW*3-1) 40,30,20
   20 INT = N/4 
      CALL R4SYN(INT,B(1),B(INT+1),B(2*INT+1),B(INT*3+1)) 
      GO TO 40
   30 INT = N/2 
      CALL R2TR(INT,B(1),B(INT+1))
   40 RETURN
      END 
C---------------------------------------------------------------------
      SUBROUTINE R2TR(INT,B0,B1)
C 
C RADIX 2 ITERATION SUBROUTINE
C 
      IMPLICIT   NONE	
      INTEGER   INT
      REAL      B0(1),B1(1) 
C
      INTEGER   K
      REAL      T
C
      DO 100 K=1,INT
      T = B0(K)+B1(K) 
      B1(K) = B0(K)-B1(K) 
      B0(K) = T 
  100 CONTINUE
      RETURN
      END 
C---------------------------------------------------------------------
      SUBROUTINE R4TR(INT,B0,B1,B2,B3)
C 
C RADIX 4 ITERATION SUBROUTINE
C 
      IMPLICIT  NONE
      INTEGER   INT
      REAL      B0(1),B1(1),B2(1),B3(1) 
C
      INTEGER   K
      REAL      R0, R1
C
      DO 200 K=1,INT
      R0 = B0(K)+B2(K)
      R1 = B1(K)+B3(K)
      B2(K) = B0(K)-B2(K) 
      B3(K) = B1(K)-B3(K) 
      B0(K) = R0+R1 
      B1(K) = R0-R1 
  200 CONTINUE
      RETURN
      END 
C---------------------------------------------------------------------
      SUBROUTINE R4SYN(INT,B0,B1,B2,B3) 
C 
C RADIX 4 SYNTHESIS SUBROUTINE
C 
      IMPLICIT  NONE
      INTEGER   INT
      REAL      B0(1),B1(1),B2(1),B3(1) 
C
      INTEGER   K
      REAL      T0, T1, T2, T3
C
      DO 200 K=1,INT
      T0 = B0(K)+B1(K)
      T1 = B0(K)-B1(K)
      T2 = B2(K)+B2(K)
      T3 = B3(K)+B3(K)
      B0(K) = T0+T2 
      B2(K) = T0-T2 
      B1(K) = T1+T3 
      B3(K) = T1-T3 
  200 CONTINUE
      RETURN
      END 
C-------------------------------------------------------------------
      SUBROUTINE R8TR(INT,NN,BR0,BR1,BR2,BR3,BR4,BR5,BR6,BR7
     1 ,BI0,BI1,BI2,BI3,BI4,BI5,BI6,BI7)
C 
C RADIX 8 ITERATION SUBROUTINE
C 
      IMPLICIT  NONE
      INTEGER   INT
      INTEGER   NN
      REAL      BR0(1),BR1(1),BR2(1),BR3(1),BR4(1), 
     1          BR5(1),BR6(1),BR7(1),BI0(1),BI1(1),BI2(1),BI3(1),
     2          BI4(1),BI5(1),BI6(1),BI7(1)
      INTEGER   L(15)
C
      INTEGER   INT8
      INTEGER   J, J0, JLAST
      INTEGER   JI, JL, JR
      INTEGER   J1, J2, J3, J4, J5, J6, J7, J8, J9, J10
      INTEGER   J11, J12, J13, J14
      INTEGER   JTHET 
      INTEGER   K0, KL
      INTEGER   K, KS
      INTEGER   L1, L2, L3, L4, L5, L6, L7, L8, L9, L10
      INTEGER   L11, L12, L13, L14, L15
      REAL      ARG
      REAL      C1, C2, C3, C4, C5, C6, C7
      REAL      C22
      REAL      PR, PI
      REAL      PIOVN, P7
      REAL      S1, S2, S3, S4, S5, S6, S7
      REAL      S22
      REAL      T0, T1, T2, T3, T4, T5, T6, T7
      REAL      TH2
      REAL      TR0, TR1, TR2, TR3, TR4, TR5, TR6, TR7
      REAL      TI0, TI1, TI2, TI3, TI4, TI5, TI6, TI7
C
      EQUIVALENCE (L15,L(1)),(L14,L(2)),(L13,L(3)),(L12,L(4)),
     1 (L11,L(5)),(L10,L(6)),(L9,L(7)),(L8,L(8)),(L7,L(9)), 
     2 (L6,L(10)),(L5,L(11)),(L4,L(12)),(L3,L(13)),(L2,L(14)),
     3 (L1,L(15)) 
C
      L(1) = NN/8 
      DO 10 K=2,15
      IF(L(K-1)-2) 11,12,13 
   11 L(K-1) = 2
   12 L(K) = 2
      GO TO 10
   13 L(K) = L(K-1)/2 
   10 CONTINUE
      PIOVN = 3.14159265358979D0/DBLE(NN) 
      P7 = 0.707106781186548D0
      C22 = 0.923879532511287D0
      S22 = 0.382683432365090D0
      JI = 3
      JL = 2
      JR = 2
      DO 70 J1=2,L1,2 
      DO 70 J2=J1,L2,L1 
      DO 70 J3=J2,L3,L2 
      DO 70 J4=J3,L4,L3 
      DO 70 J5=J4,L5,L4 
      DO 70 J6=J5,L6,L5 
      DO 70 J7=J6,L7,L6 
      DO 70 J8=J7,L8,L7 
      DO 70 J9=J8,L9,L8 
      DO 70 J10=J9,L10,L9 
      DO 70 J11=J10,L11,L10 
      DO 70 J12=J11,L12,L11 
      DO 70 J13=J12,L13,L12 
      DO 70 J14=J13,L14,L13 
      DO 70 JTHET=J14,L15,L14 
      TH2 = JTHET-2 
      IF(TH2) 71,71,76
   71 DO 72 K=1,INT 
      T0 = BR0(K)+BR4(K)
      T1 = BR1(K)+BR5(K)
      T2 = BR2(K)+BR6(K)
      T3 = BR3(K)+BR7(K)
      T4 = BR0(K)-BR4(K)
      T5 = BR1(K)-BR5(K)
      T6 = BR2(K)-BR6(K)
      T7 = BR3(K)-BR7(K)
      BR2(K) = T0-T2
      BR3(K) = T1-T3
      T0 = T0+T2
      T1 = T1+T3
      BR0(K) = T0+T1
      BR1(K) = T0-T1
      PR = P7*(T5-T7) 
      PI = P7*(T5+T7) 
      BR4(K) = T4+PR
      BR7(K) = T6+PI
      BR6(K) = T4-PR
      BR5(K) = PI-T6
   72 CONTINUE
      IF(NN-8) 70,70,73 
   73 K0 = INT*8+1
      KL = K0+INT-1 
      DO 75 K=K0,KL 
      PR = P7*(BI2(K)-BI6(K)) 
      PI = P7*(BI2(K)+BI6(K)) 
      TR0 = BI0(K)+PR 
      TI0 = BI4(K)+PI 
      TR2 = BI0(K)-PR 
      TI2 = BI4(K)-PI 
      PR = P7*(BI3(K)-BI7(K)) 
      PI = P7*(BI3(K)+BI7(K)) 
      TR1 = BI1(K)+PR 
      TI1 = BI5(K)+PI 
      TR3 = BI1(K)-PR 
      TI3 = BI5(K)-PI 
      PR = TR1*C22-TI1*S22
      PI = TI1*C22+TR1*S22
      BI0(K) = TR0+PR 
      BI6(K) = TR0-PR 
      BI7(K) = TI0+PI 
      BI1(K) = PI-TI0 
      PR = -TR3*S22-TI3*C22 
      PI =  TR3*C22-TI3*S22 
      BI2(K) = TR2+PR 
      BI4(K) = TR2-PR 
      BI5(K) = TI2+PI 
      BI3(K) = PI-TI2 
   75 CONTINUE
      GO TO 70
   76 ARG = TH2*PIOVN 
      C1 = COS(ARG) 
      S1 = SIN(ARG) 
      C2 = C1*C1-S1*S1
      S2 = C1*S1+C1*S1
      C3 = C1*C2-S1*S2
      S3 = C2*S1+S2*C1
      C4 = C2*C2-S2*S2
      S4 = C2*S2+C2*S2
      C5 = C2*C3-S2*S3
      S5 = C3*S2+S3*C2
      C6 = C3*C3-S3*S3
      S6 = C3*S3+C3*S3
      C7 = C3*C4-S3*S4
      S7 = C4*S3+S4*C3
      INT8 = INT*8
      J0 = JR*INT8+1
      K0 = JI*INT8+1
      JLAST = J0+INT-1
      DO 77 J=J0,JLAST
      KS = K0+J-J0
      DO 77 K=KS,KS 
      TR1 = BR1(J)*C1-BI1(K)*S1 
      TI1 = BR1(J)*S1+BI1(K)*C1 
      TR2 = BR2(J)*C2-BI2(K)*S2 
      TI2 = BR2(J)*S2+BI2(K)*C2 
      TR3 = BR3(J)*C3-BI3(K)*S3 
      TI3 = BR3(J)*S3+BI3(K)*C3 
      TR4 = BR4(J)*C4-BI4(K)*S4 
      TI4 = BR4(J)*S4+BI4(K)*C4 
      TR5 = BR5(J)*C5-BI5(K)*S5 
      TI5 = BR5(J)*S5+BI5(K)*C5 
      TR6 = BR6(J)*C6-BI6(K)*S6 
      TI6 = BR6(J)*S6+BI6(K)*C6 
      TR7 = BR7(J)*C7-BI7(K)*S7 
      TI7 = BR7(J)*S7+BI7(K)*C7 
      T0 = BR0(J)+TR4 
      T1 = BI0(K)+TI4 
      TR4 = BR0(J)-TR4
      TI4 = BI0(K)-TI4
      T2 = TR1+TR5
      T3 = TI1+TI5
      TR5 = TR1-TR5 
      TI5 = TI1-TI5 
      T4 = TR2+TR6
      T5 = TI2+TI6
      TR6 = TR2-TR6 
      TI6 = TI2-TI6 
      T6 = TR3+TR7
      T7 = TI3+TI7
      TR7 = TR3-TR7 
      TI7 = TI3-TI7 
      TR0 = T0+T4 
      TI0 = T1+T5 
      TR2 = T0-T4 
      TI2 = T1-T5 
      TR1 = T2+T6 
      TI1 = T3+T7 
      TR3 = T2-T6 
      TI3 = T3-T7 
      T0 = TR4-TI6
      T1 = TI4+TR6
      T4 = TR4+TI6
      T5 = TI4-TR6
      T2 = TR5-TI7
      T3 = TI5+TR7
      T6 = TR5+TI7
      T7 = TI5-TR7
      BR0(J) = TR0+TR1
      BI7(K) = TI0+TI1
      BI6(K) = TR0-TR1
      BR1(J) = TI1-TI0
      BR2(J) = TR2-TI3
      BI5(K) = TI2+TR3
      BI4(K) = TR2+TI3
      BR3(J) = TR3-TI2
      PR = P7*(T2-T3) 
      PI = P7*(T2+T3) 
      BR4(J) = T0+PR
      BI3(K) = T1+PI
      BI2(K) = T0-PR
      BR5(J) = PI-T1
      PR = -P7*(T6+T7)
      PI = P7*(T6-T7) 
      BR6(J) = T4+PR
      BI1(K) = T5+PI
      BI0(K) = T4-PR
      BR7(J) = PI-T5
   77 CONTINUE
      JR = JR+2 
      JI = JI-2 
      IF (JI-JL) 78,78,70
   78 JI = 2*JR-1 
      JL = JR 
   70 CONTINUE
      RETURN
      END 
C-----------------------------------------------------------------
      SUBROUTINE R8SYN(INT,NN,BR0,BR1,BR2,BR3,BR4,BR5,BR6,BR7 
     1 ,BI0,BI1,BI2,BI3,BI4,BI5,BI6,BI7)
C 
C RADIX 8 SYNTHESIS SUBROUTINE
C

      IMPLICIT  NONE
      INTEGER   INT
      INTEGER   NN
      REAL      BR0(1),BR1(1),BR2(1),BR3(1),BR4(1),
     1          BR5(1),BR6(1),BR7(1),BI0(1),BI1(1),BI2(1),BI3(1),
     2          BI4(1),BI5(1),BI6(1),BI7(1)
      INTEGER   L(15)
C
      INTEGER   INT8
      INTEGER   J, J0, JLAST
      INTEGER   JI, JL, JR
      INTEGER   J1, J2, J3, J4, J5, J6, J7, J8, J9, J10
      INTEGER   J11, J12, J13, J14
      INTEGER   JTHET
      INTEGER   K0, KL
      INTEGER   K, KS
      INTEGER   L1, L2, L3, L4, L5, L6, L7, L8, L9, L10
      INTEGER   L11, L12, L13, L14, L15
      REAL      ARG
      REAL      C1, C2, C3, C4, C5, C6, C7
      REAL      C22
      REAL      PR, PI
      REAL      PIOVN, P7
      REAL      P7TWO
      REAL      RI, RR
      REAL      S1, S2, S3, S4, S5, S6, S7
      REAL      S22
      REAL      T0, T1, T2, T3, T4, T5, T6, T7
      REAL      TH2
      REAL      TR0, TR1, TR2, TR3, TR4, TR5, TR6, TR7
      REAL      TI0, TI1, TI2, TI3, TI4, TI5, TI6, TI7
      REAL      T8, TT0, TT1
      REAL      TTR6, TTR7
C 
      EQUIVALENCE (L15,L(1)),(L14,L(2)),(L13,L(3)),(L12,L(4)),
     1 (L11,L(5)),(L10,L(6)),(L9,L(7)),(L8,L(8)),(L7,L(9)), 
     2 (L6,L(10)),(L5,L(11)),(L4,L(12)),(L3,L(13)),(L2,L(14)),
     3 (L1,L(15)) 
C
      L(1) = NN/8 
      DO 10 K=2,15
      IF (L(K-1)-2) 11,12,13 
   11 L(K-1) = 2
   12 L(K) = 2
      GO TO 10
   13 L(K) = L(K-1)/2 
   10 CONTINUE
      PIOVN = 3.14159265358979D0/DBLE(NN) 
      P7 = 0.707106781186548D0
      P7TWO = 1.414213562373095D0
      C22 = 0.923879532511287D0
      S22 = 0.382683432365090D0
      JI = 3
      JL = 2
      JR = 2
      DO 70 J1=2,L1,2 
      DO 70 J2=J1,L2,L1 
      DO 70 J3=J2,L3,L2 
      DO 70 J4=J3,L4,L3 
      DO 70 J5=J4,L5,L4 
      DO 70 J6=J5,L6,L5 
      DO 70 J7=J6,L7,L6 
      DO 70 J8=J7,L8,L7 
      DO 70 J9=J8,L9,L8 
      DO 70 J10=J9,L10,L9 
      DO 70 J11=J10,L11,L10 
      DO 70 J12=J11,L12,L11 
      DO 70 J13=J12,L13,L12 
      DO 70 J14=J13,L14,L13 
      DO 70 JTHET=J14,L15,L14 
      TH2 = JTHET-2 
      IF (TH2) 71,71,76
   71 DO 72 K=1,INT 
      T0 = BR0(K)+BR1(K)
      T1 = BR0(K)-BR1(K)
      T2 = BR2(K)+BR2(K)
      T3 = BR3(K)+BR3(K)
      T4 = BR4(K)+BR6(K)
      T6 = BR7(K)-BR5(K)
      T5 = BR4(K)-BR6(K)
      T7 = BR7(K)+BR5(K)
      PR = P7*(T7+T5) 
      PI = P7*(T7-T5) 
      TT0 = T0+T2 
      TT1 = T1+T3 
      T2 = T0-T2
      T3 = T1-T3
      T4 = T4+T4
      T5 = PR+PR
      T6 = T6+T6
      T7 = PI+PI
      BR0(K) = TT0+T4 
      BR1(K) = TT1+T5 
      BR2(K) = T2+T6
      BR3(K) = T3+T7
      BR4(K) = TT0-T4 
      BR5(K) = TT1-T5 
      BR6(K) = T2-T6
      BR7(K) = T3-T7
   72 CONTINUE
      IF (NN-8) 70,70,73 
   73 K0 = INT*8+1
      KL = K0+INT-1 
      DO 75 K=K0,KL 
      T1 = BI0(K)+BI6(K)
      T2 = BI7(K)-BI1(K)
      T3 = BI0(K)-BI6(K)
      T4 = BI7(K)+BI1(K)
      PR = T3*C22+T4*S22
      PI = T4*C22-T3*S22
      T5 = BI2(K)+BI4(K)
      T6 = BI5(K)-BI3(K)
      T7 = BI2(K)-BI4(K)
      T8 = BI5(K)+BI3(K)
      RR = T8*C22-T7*S22
      RI = -T8*S22-T7*C22 
      BI0(K) = (T1+T5)+(T1+T5)
      BI4(K) = (T2+T6)+(T2+T6)
      BI1(K) = (PR+RR)+(PR+RR)
      BI5(K) = (PI+RI)+(PI+RI)
      T5 = T1-T5
      T6 = T2-T6
      BI2(K) = P7TWO*(T6+T5)
      BI6(K) = P7TWO*(T6-T5)
      RR = PR-RR
      RI = PI-RI
      BI3(K) = P7TWO*(RI+RR)
      BI7(K) = P7TWO*(RI-RR)
   75 CONTINUE
      GO TO 70
   76 ARG = TH2*PIOVN 
      C1 = COS(ARG) 
      S1 = -SIN(ARG)
      C2 = C1*C1-S1*S1
      S2 = C1*S1+C1*S1
      C3 = C1*C2-S1*S2
      S3 = C2*S1+S2*C1
      C4 = C2*C2-S2*S2
      S4 = C2*S2+C2*S2
      C5 = C2*C3-S2*S3
      S5 = C3*S2+S3*C2
      C6 = C3*C3-S3*S3
      S6 = C3*S3+C3*S3
      C7 = C3*C4-S3*S4
      S7 = C4*S3+S4*C3
      INT8 = INT*8
      J0 = JR*INT8+1
      K0 = JI*INT8+1
      JLAST = J0+INT-1
      DO 77 J=J0,JLAST
      KS = K0+J-J0
      DO 77 K=KS,KS 
      TR0 = BR0(J)+BI6(K) 
      TI0 = BI7(K)-BR1(J) 
      TR1 = BR0(J)-BI6(K) 
      TI1 = BI7(K)+BR1(J) 
      TR2 = BR2(J)+BI4(K) 
      TI2 = BI5(K)-BR3(J) 
      TR3 = BI5(K)+BR3(J) 
      TI3 = BI4(K)-BR2(J) 
      TR4 = BR4(J)+BI2(K) 
      TI4 = BI3(K)-BR5(J) 
      T0 = BR4(J)-BI2(K)
      T1 = BI3(K)+BR5(J)
      TR5 = P7*(T0+T1)
      TI5 = P7*(T1-T0)
      TR6 = BR6(J)+BI0(K) 
      TI6 = BI1(K)-BR7(J) 
      T0 = BR6(J)-BI0(K)
      T1 = BI1(K)+BR7(J)
      TR7 = -P7*(T0-T1) 
      TI7 = -P7*(T1+T0) 
      T0 = TR0+TR2
      T1 = TI0+TI2
      T2 = TR1+TR3
      T3 = TI1+TI3
      TR2 = TR0-TR2 
      TI2 = TI0-TI2 
      TR3 = TR1-TR3 
      TI3 = TI1-TI3 
      T4 = TR4+TR6 
      T5 = TI4+TI6 
      T6 = TR5+TR7 
      T7 = TI5+TI7 
      TTR6 = TI4-TI6
      TI6 = TR6-TR4 
      TTR7 = TI5-TI7
      TI7 = TR7-TR5 
      BR0(J) = T0+T4
      BI0(K) = T1+T5
      BR1(J) = C1*(T2+T6)-S1*(T3+T7)
      BI1(K) = C1*(T3+T7)+S1*(T2+T6)
      BR2(J) = C2*(TR2+TTR6)-S2*(TI2+TI6) 
      BI2(K) = C2*(TI2+TI6)+S2*(TR2+TTR6) 
      BR3(J) = C3*(TR3+TTR7)-S3*(TI3+TI7) 
      BI3(K) = C3*(TI3+TI7)+S3*(TR3+TTR7) 
      BR4(J) = C4*(T0-T4)-S4*(T1-T5)
      BI4(K) = C4*(T1-T5)+S4*(T0-T4)
      BR5(J) = C5*(T2-T6)-S5*(T3-T7)
      BI5(K) = C5*(T3-T7)+S5*(T2-T6)
      BR6(J) = C6*(TR2-TTR6)-S6*(TI2-TI6) 
      BI6(K) = C6*(TI2-TI6)+S6*(TR2-TTR6) 
      BR7(J) = C7*(TR3-TTR7)-S7*(TI3-TI7) 
      BI7(K) = C7*(TI3-TI7)+S7*(TR3-TTR7) 
   77 CONTINUE
      JR = JR+2 
      JI = JI-2 
      IF (JI-JL) 78,78,70
   78 JI = 2*JR-1 
      JL = JR 
   70 CONTINUE
      RETURN
      END 
C---------------------------------------------------------------------
      SUBROUTINE ORD1(M,B)
C 
C IN-PLACE REORDERING SUBROUTINES 
C 
      IMPLICIT  NONE
      INTEGER   M
      REAL      B(1)
C
      INTEGER   K, KL
      INTEGER   N
      INTEGER   J
      REAL      T
C
      K = 4 
      KL = 2
      N = 2**M
      DO 94 J=4,N,2 
      IF (K-J) 92,92,91 
   91 T = B(J)
      B(J) = B(K) 
      B(K) = T
   92 K = K-2 
      IF (K-KL) 93,93,94
   93 K = 2*J 
      KL = J
   94 CONTINUE
      RETURN
      END 
C---------------------------------------------------------------------
      SUBROUTINE ORD2(M,B)
C
      IMPLICIT NONE
      INTEGER  M
      REAL     B(1)
C
      INTEGER  IJ
      INTEGER  J1, J2, J3, J4, J5, J6, J7, J8, J9, J10
      INTEGER  J11, J12, J13, J14, J15, J16, J17, J18
      INTEGER  L1, L2, L3, L4, L5, L6, L7, L8, L9, L10
      INTEGER  L11, L12, L13, L14, L15, L16, L17, L18, L19
      INTEGER  JI
      INTEGER  L(20)
      INTEGER  N
      INTEGER  K
      REAL     T
C
      EQUIVALENCE (L19,L(1)),(L18,L(2)),(L17,L(3)),(L16,L(4)),
     1 (L15,L(5)),(L14,L(6)),(L13,L(7)),(L12,L(8)),(L11,L(9)),
     2 (L10,L(10)),(L9,L(11)),(L8,L(12)),(L7,L(13)),(L6,L(14)), 
     3 (L5,L(15)),(L4,L(16)),(L3,L(17)),(L2,L(18)),(L1,L(19)) 
C
      N = 2**M
      L(1) = N
      DO 101 K=2,M
  101 L(K) = L(K-1)/2 
      DO 102 K=M,19 
  102 L(K+1) = 2
      IJ = 2
      DO 103 J1=2,L1,2
      DO 103 J2=J1,L2,L1
      DO 103 J3=J2,L3,L2
      DO 103 J4=J3,L4,L3
      DO 103 J5=J4,L5,L4
      DO 103 J6=J5,L6,L5
      DO 103 J7=J6,L7,L6
      DO 103 J8=J7,L8,L7
      DO 103 J9=J8,L9,L8
      DO 103 J10=J9,L10,L9
      DO 103 J11=J10,L11,L10
      DO 103 J12=J11,L12,L11
      DO 103 J13=J12,L13,L12
      DO 103 J14=J13,L14,L13
      DO 103 J15=J14,L15,L14
      DO 103 J16=J15,L16,L15
      DO 103 J17=J16,L17,L16
      DO 103 J18=J17,L18,L17
      DO 103 JI =J18,L19,L18
      IF (IJ-JI) 108,103,103
 108  T = B(IJ-1) 
      B(IJ-1) = B(JI-1) 
      B(JI-1) = T 
      T = B(IJ) 
      B(IJ) = B(JI) 
      B(JI) = T 
  103 IJ = IJ+2 
      RETURN
      END 
