C===========================================================================
C Copyright (C) 1995-2005 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 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence 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
      PROGRAM FOURI
C
C++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  program FOURI		version 1.00	881028
C  K. Banse                  	ESO - Garching
C
C.KEYWORDS
C  Fourier transformation
C
C.PURPOSE
C  calculate the Fourier transformation of given real + imaginary frame
C
C.ALGORITHM
C  use the FFT algorithm as in E.O. Brigham: "The Fast Fourier Transform",
C  Prentice-Hall, Inc., 1974, page 164
C
C.INPUT/OUTPUT
C  the following keywords are used:
C
C  P1/C/1/120			name of real data frame
C  P2/C/1/120			name of imaginary data frame
C                               if IN_B is "+", only real input data
C  P3/C/1/120			name of frame holding real part of FT
C  P4/C/1/120			name of frame holding imaginary part of FT
C
C  CACTIO/C/1/3			'???' = normal way
C				'INV' = inverse FFT
C				'POW' = return also power spectrum
C				'FRE' = use frequencies 
C				        (center shifted to the middle)
C				'FPO' = FRE + POW combined
C				'FIN' = FRE + INV combined
C  TMPBUF/C/1/3			(1:1) = spare
C				(2:2) = I for inverse, P for power stuff
C				(3:3) will be set to Y, if expansion
C				      to a power2 frame took place or not
C	
C  P5/C/1/120			if action = POWER or FPOWER,
C				name of power spectrum
C
C.VERSIONS
C 
C 050629	last modif
C	
C--------------------------------------------------
C
      IMPLICIT NONE
C
      INTEGER    IAV,M,N,NAXISA,NAXISB,NN,NOPIX
      INTEGER    IMNOA,IMNOB,IMNOC,IMNOD,IMNOP
      INTEGER*8  PNTRA,PNTRB,PNTRC,PNTRD,PNTRP
      INTEGER*8  DI,DR,CT,DT,TS,TC,SWAP1,SWAP2
      INTEGER    IMNODI,IMNODR,IMNOTC,IMNOTS
      INTEGER    SIZE,DSTAT,STAT
      INTEGER    NLOG(3),NPIXA(3),NPIXB(3),NPIXF(3),EXTFLG,POW2(22)
      INTEGER    LOW(3)
      INTEGER    UNI(1),NULO,MADRID(1)
      INTEGER    EC,EL,ED
C
      CHARACTER*120 WFRAME,FRAMEA,FRAMEB,FRAMEC,FRAMED,FRAMEP
      CHARACTER     IDENT*72,CUNIT*48,ACTION*3,TMP*3
C
      DOUBLE PRECISION STARTA(3),STARTB(3),STEPA(3),STEPB(3)
      DOUBLE PRECISION DBUF(6)
C
      REAL      CUTS(4),RC
C
      COMMON  /VMR/ MADRID
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
C
      DATA      FRAMEA /' '/,    FRAMEB /' '/
      DATA      FRAMEC /' '/,    FRAMED /' '/
      DATA      IDENT  /' '/,    CUNIT  /' '/
      DATA      ACTION /' '/
      DATA      NPIXA  /3*1/,    NPIXB  /3*1/,    NPIXF  /3*1/
      DATA      STARTA /3*0.D0/, STARTB /3*0.D0/
      DATA      STEPA  /3*1.D0/, STEPB  /3*1.D0/
      DATA      EXTFLG /0/
      DATA      POW2   /1,2,4,8,16,32,64,128,256,512,1024,2048,4096,
     +                  8192,16384,32768,65536,131072,262144,524288,
     +                  1048576,2097152/
      DATA      RC     /0.0/
      DATA      LOW    /3*1/     NLOG   /3*0/
      DATA      CT     /1/,      DT     /1/
C
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
C  get into MIDAS environment
      CALL STSPRO('FOURI')
C
C  get action
      CALL STKRDC('CACTIO',1,1,3,IAV,ACTION,UNI,NULO,STAT)
      CALL UPCAS(ACTION,ACTION)
      CALL STKRDC('TMPBUF',1,1,3,IAV,TMP,UNI,NULO,STAT)
C
C  get data frames involved
      CALL STKRDC('P1',1,1,120,IAV,FRAMEA,UNI,NULO,STAT)
      CALL STKRDC('P2',1,1,120,IAV,FRAMEB,UNI,NULO,STAT)
      CALL STKRDC('P3',1,1,120,IAV,FRAMEC,UNI,NULO,STAT)
      CALL STKRDC('P4',1,1,120,IAV,FRAMED,UNI,NULO,STAT)
C
C  map real input frame
      CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +            3,NAXISA,NPIXA,STARTA,STEPA,
     +            IDENT,CUNIT,PNTRA,IMNOA,STAT)
C
      IF (NAXISA.GT.2) 
     +   CALL STETER(1,'only 2-dim frames supported currently...')
C
C  get imaginary part
      IF (FRAMEB(1:1).NE.'+') THEN
         CALL STIGET(FRAMEB,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +               3,NAXISB,NPIXB,STARTB,STEPB,
     +               IDENT,CUNIT,PNTRB,IMNOB,STAT)
         IF (NAXISB.NE.NAXISA) GOTO 9000
C
         DO 200, N=1,NAXISA                      !frames have to be identical
            IF (NPIXA(N).NE.NPIXB(N)) GOTO 9000
            IF (ACTION(1:1).EQ.'F') THEN
               IF (ABS(STEPA(N)-STEPB(N)).GT.10.D-15) GOTO 9000
               IF (ABS(STARTA(N)-STARTB(N)).GT.10.D-15) GOTO 9000
            ENDIF
200      CONTINUE
      ENDIF
C
C  for inverse FFT check descriptor FFT_SPECIAL with original START, STEP
C  and eventual NPIX before expansion
C
      DBUF(1) = STARTA(1)
      DBUF(2) = STARTA(2)
      DBUF(3) = STEPA(1)
      DBUF(4) = STEPA(2)
      DBUF(5) = DBLE(NPIXA(1))
      DBUF(6) = DBLE(NPIXA(2))
C 
      IF (TMP(2:2).EQ.'I') THEN
         CALL STECNT('GET',EC,EL,ED)
         CALL STECNT('PUT',1,0,0)             !test, if FFT_SPECIAL is there
         CALL STDRDD(IMNOA,'FFT_SPECIAL',1,6,IAV,DBUF,UNI,NULO,DSTAT)
         CALL STECNT('PUT',EC,EL,ED)
         IF (DSTAT.EQ.0) THEN
            STARTA(1) = DBUF(1)
            STARTA(2) = DBUF(2)
            STEPA(1) = DBUF(3)
            STEPA(2) = DBUF(4)
         ENDIF
      ELSE
         DSTAT = 1                            !indicate error...
      ENDIF
C
C  calculate log of dimension(s)
      DO 400, M=1,NAXISA
         NN = NPIXA(M)
         NPIXF(M) = NN
         DO 300, N=1,22
            IF (NN.LE.POW2(N)) THEN
               IF (NN.LT.POW2(N)) THEN
                  NPIXF(M) = POW2(N)
                  EXTFLG = 1
               ENDIF
               NLOG(M) = N - 1
               GOTO 400                           !jump out of loop
            ENDIF
300      CONTINUE
         CALL STETER(3,'x- or y-dimension larger than 2097152 ...')
400   CONTINUE
C
C  calculate size
      IF (NAXISA.EQ.1) THEN                     !max. 2-dim frame...
         SIZE = NPIXF(1)
      ELSE
         SIZE = NPIXF(1) * NPIXF(2)
         CALL STFXMP(2*SIZE,D_R4_FORMAT,CT,STAT)
         DT = CT + SIZE
      ENDIF
C
C  see, if we have to expand real input frame
      IF (EXTFLG.EQ.1) THEN
         CALL STFCRE('wworrkfr',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +               SIZE,IMNOC,STAT)
         CALL STFMAP(IMNOC,F_X_MODE,1,SIZE,IAV,PNTRC,STAT)
         CALL OPFC('FC=',0.0,RC,MADRID(PNTRC),SIZE,RC,STAT)    !fill with zero
         CALL COPWND(MADRID(PNTRA),NPIXA,MADRID(PNTRC),NPIXF,
     +               LOW,LOW,NPIXA)                   !and insert input frame
         CALL STFCLO(IMNOA,STAT)
         IMNOA = IMNOC
         PNTRA = PNTRC
C
         IF (FRAMEB(1:1).NE.'+') THEN              !also expand imaginary file
            CALL STFCRE('wworrkfi',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,
     +                  SIZE,IMNOC,STAT)
            CALL STFMAP(IMNOC,F_X_MODE,1,SIZE,IAV,PNTRC,STAT)
            CALL OPFC('FC=',0.0,RC,MADRID(PNTRC),SIZE,RC,STAT)   !fill with zero
            CALL COPWND(MADRID(PNTRB),NPIXA,MADRID(PNTRC),NPIXF,
     +                  LOW,LOW,NPIXA)                   !and insert 
            CALL STFCLO(IMNOB,STAT)
            IMNOB = IMNOC
            PNTRB = PNTRC
         ENDIF
      ENDIF
C
C  update start, step for forward FFT
      IF (TMP(2:2).NE.'I') THEN
         DO 2400, N=1,NAXISA
            STEPA(N) = ABS( 1.D0 / (DBLE(NPIXF(N)) * STEPA(N)) )
            IF (ACTION(1:1).EQ.'F') THEN
               STARTA(N) = - ( (NPIXF(N)/2) * STEPA(N) )
            ELSE
               STARTA(N) = 0.0D0
            ENDIF
2400     CONTINUE
      ENDIF
C
C  create + map result frames
      IF (EXTFLG.EQ.1) THEN
         WFRAME(1:) = 'exp_fftr '
      ELSE
         WFRAME(1:) = FRAMEC(1:)
      ENDIF
      CALL STIPUT(WFRAME,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +            NAXISA,NPIXF,STARTA,STEPA,
     +            IDENT,CUNIT,PNTRC,IMNOC,STAT)
      IF (EXTFLG.EQ.1) THEN
         WFRAME(1:) = 'exp_ffti '
      ELSE
         WFRAME(1:) = FRAMED(1:)
      ENDIF
      CALL STIPUT(WFRAME,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +            NAXISA,NPIXF,STARTA,STEPA,
     +            IDENT,CUNIT,PNTRD,IMNOD,STAT)
C
C  allocate memory for working buffer DREAL + DIMAG 
C  as well as for sine + cosine look-up tables (all double prec. data)
C
      NOPIX = MAX(NPIXF(1),NPIXF(2))
      CALL STFCRE('wworkkdr',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
     +            NOPIX,IMNODR,STAT)
      CALL STFMAP(IMNODR,F_X_MODE,1,NOPIX,IAV,DR,STAT)
      CALL STFCRE('wworkkdi',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
     +            NOPIX,IMNODI,STAT)
      CALL STFMAP(IMNODI,F_X_MODE,1,NOPIX,IAV,DI,STAT)
C 
      CALL STFCRE('wworkkts',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
     +            NOPIX,IMNOTS,STAT)
      CALL STFMAP(IMNOTS,F_X_MODE,1,NOPIX,IAV,TS,STAT)
      CALL STFCRE('wworkktc',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
     +            NOPIX,IMNOTC,STAT)
      CALL STFMAP(IMNOTC,F_X_MODE,1,NOPIX,IAV,TC,STAT)
C 
C  and allocate space for swap tables
      CALL STFXMP(NOPIX,D_I4_FORMAT,SWAP1,STAT)
      SWAP2 = SWAP1 + NOPIX/2
C
C  now do it in a subroutine
      IF (FRAMEB(1:1).NE.'+') THEN                  !imaginary part exists...
         IF (TMP(2:2).EQ.'I') THEN
            CALL DFFTB(MADRID(PNTRA),MADRID(PNTRB),
     +                 MADRID(PNTRC),MADRID(PNTRD),
     +                 NPIXF,ACTION,MADRID(DR),MADRID(DI),NLOG,POW2,
     +                 MADRID(TS),MADRID(TC),
     +                 MADRID(CT),MADRID(DT),
     +                 MADRID(SWAP1),MADRID(SWAP2))
         ELSE
            CALL DFFTF(MADRID(PNTRA),MADRID(PNTRB),
     +                 MADRID(PNTRC),MADRID(PNTRD),
     +                 NPIXF,ACTION,MADRID(DR),MADRID(DI),NLOG,POW2,
     +                 MADRID(TS),MADRID(TC),
     +                 MADRID(CT),MADRID(DT),
     +                 MADRID(SWAP1),MADRID(SWAP2))
         ENDIF
         CALL STFCLO(IMNOB,STAT)
      ELSE                                          !just real data
         IF (TMP(2:2).EQ.'I') THEN
            CALL DRFFTB(MADRID(PNTRA),MADRID(PNTRC),
     +                  MADRID(PNTRD),NPIXF,ACTION,
     +                  MADRID(DR),MADRID(DI),NLOG,POW2,
     +                  MADRID(TS),MADRID(TC),
     +                  MADRID(CT),MADRID(DT),
     +                  MADRID(SWAP1),MADRID(SWAP2))
         ELSE
            CALL DRFFTF(MADRID(PNTRA),MADRID(PNTRC),
     +                  MADRID(PNTRD),NPIXF,ACTION,
     +                  MADRID(DR),MADRID(DI),NLOG,POW2,
     +                  MADRID(TS),MADRID(TC),
     +                  MADRID(CT),MADRID(DT),
     +                  MADRID(SWAP1),MADRID(SWAP2))
         ENDIF
      ENDIF
C 
C  get rid of all the allocated temporary space
      CALL STFCLO(IMNODR,STAT)
      CALL STFCLO(IMNODI,STAT)
      CALL STFCLO(IMNOTS,STAT)
      CALL STFCLO(IMNOTC,STAT)
C 
      CALL STFCLO(IMNOA,STAT)
      CALL STFOPN(FRAMEA,D_OLD_FORMAT,0,F_IMA_TYPE,IMNOA,STAT)
      CALL DSCUPT(IMNOA,IMNOC,' ',STAT)
      CALL STDWRD(IMNOC,'FFT_SPECIAL',DBUF,1,6,UNI,STAT)
C
C  if we also want power spectrum (absolute value), do it now...
      IF (TMP(2:2).EQ.'P') THEN
         CALL STKRDC('P5',1,1,120,IAV,FRAMEP,UNI,NULO,STAT)
         IF (EXTFLG.EQ.1) THEN
            WFRAME(1:) = 'exp_pow '
         ELSE
            WFRAME(1:) = FRAMEP(1:)
         ENDIF
         CALL STIPUT(WFRAME,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +               NAXISA,NPIXF,STARTA,STEPA,
     +               IDENT,CUNIT,DR,IMNODR,STAT)
         CALL POWER(MADRID(PNTRC),MADRID(PNTRD),MADRID(DR),SIZE)
         CALL DSCUPT(IMNOA,IMNODR,' ',STAT)
         CALL STDWRD(IMNODR,'FFT_SPECIAL',DBUF,1,6,UNI,STAT)
C 
         IF (EXTFLG.EQ.1) THEN
            CALL STIPUT(FRAMEP,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +                  NAXISA,NPIXA,STARTA,STEPA,
     +                  IDENT,CUNIT,PNTRP,IMNOP,STAT)
            CALL COPWND(MADRID(DR),NPIXF,MADRID(PNTRP),NPIXA,
     +                  LOW,LOW,NPIXA) 
            SIZE = NPIXA(1)*NPIXA(2)
            CALL MYMX(MADRID(PNTRP),SIZE,CUTS(1))
         ELSE
            CALL MYMX(MADRID(DR),SIZE,CUTS(1))
         ENDIF
         CUTS(3) = CUTS(1)
         CUTS(4) = CUTS(2)
         CALL STDWRR(IMNODR,'LHCUTS',CUTS,1,4,UNI,STAT)
         CALL STFCLO(IMNODR,STAT)
         IF (EXTFLG.EQ.1) THEN                     !also update extracted frame
            CALL DSCUPT(IMNOA,IMNOP,' ',STAT)
            CALL STDWRR(IMNOP,'LHCUTS',CUTS,1,4,UNI,STAT)
            CALL STDWRD(IMNOP,'FFT_SPECIAL',DBUF,1,6,UNI,STAT)
            CALL STFCLO(IMNOP,STAT)
         ENDIF
      ENDIF
C 
C  for inverse FFT, check DBUF(5,6) to see if we have to extract frames
      IF ((TMP(2:2).EQ.'I').AND.(DSTAT.EQ.0)) THEN
         NLOG(1) = NINT(DBUF(5))
         NLOG(2) = NINT(DBUF(6))
         IF ((NLOG(1).NE.NPIXA(1)) .OR. (NLOG(2).NE.NPIXA(2))) THEN
            FRAMEC(1:) = 'orig_r '
            FRAMED(1:) = 'orig_i '
            NPIXA(1) = NLOG(1)
            NPIXA(2) = NLOG(2)
            EXTFLG = 2
         ENDIf
      ENDIF
C 
C  if necessary, extract result frames
      IF (EXTFLG.GE.1) THEN
         CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +               NAXISA,NPIXA,STARTA,STEPA,
     +               IDENT,CUNIT,DR,IMNODR,STAT)
         CALL COPWND(MADRID(PNTRC),NPIXF,MADRID(DR),NPIXA,
     +               LOW,LOW,NPIXA) 
         CALL STFCLO(IMNOC,STAT)
         IMNOC = IMNODR
         PNTRC = DR
         CALL DSCUPT(IMNOA,IMNOC,' ',STAT)
         CALL STDWRD(IMNOC,'FFT_SPECIAL',DBUF,1,6,UNI,STAT)
C 
         CALL STIPUT(FRAMED,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +               NAXISA,NPIXA,STARTA,STEPA,
     +               IDENT,CUNIT,PNTRB,IMNOB,STAT)
         CALL COPWND(MADRID(PNTRD),NPIXF,MADRID(PNTRB),NPIXA,
     +               LOW,LOW,NPIXA) 
         CALL STFCLO(IMNOD,STAT)
         IMNOD = IMNOB
         PNTRD = PNTRB
         SIZE = NPIXA(1)*NPIXA(2)
         IF (EXTFLG.EQ.1) THEN
            TMP(3:3) = 'Y'                   !tell fft.prg about the extension
         ELSE
            TMP(3:3) = 'Z'                   !tell fft.prg about the extraction
         ENDIF
         CALL STKWRC('TMPBUF',1,TMP,1,3,UNI,STAT)
      ENDIF
C
C  get dynamic range
      CALL MYMX(MADRID(PNTRC),SIZE,CUTS(1))
      CUTS(3) = CUTS(1)
      CUTS(4) = CUTS(2)
      CALL STDWRR(IMNOC,'LHCUTS',CUTS,1,4,UNI,STAT)
      CALL MYMX(MADRID(PNTRD),SIZE,CUTS(1))
      CUTS(3) = CUTS(1)
      CUTS(4) = CUTS(2)
      CALL STDWRR(IMNOD,'LHCUTS',CUTS,1,4,UNI,STAT)
C
C  That's it folks...
      CALL STSEPI
C 
C  here if real + imaginary input frames don't match
9000  CALL STETER(2,'input frames do not have same dimensions...')
      END

      SUBROUTINE DFFTF(A,B,C,D,NPIX,ACT,DREAL,DIMAG,NLOG,POW2,
     +                 TSIN,TCOS,CT,DT,SWAPI,SWAPO)
C
      IMPLICIT NONE
C
      INTEGER    INDX,N,NDIM,NX,NY,MM
      INTEGER    XFROM,NSIZE(2),NWORK
      INTEGER    KMM,MPOW(21)
      INTEGER    NPIX(2),NLOG(*),POW2(*)
      INTEGER    SWAPI(*),SWAPO(*),SWALIM
C
      REAL      A(*),B(*),C(*),D(*),CT(*),DT(*)
C
      DOUBLE PRECISION   DREAL(*),DIMAG(*),TSIN(*),TCOS(*)
      DOUBLE PRECISION   DFACT
      DOUBLE PRECISION   DX,DY
C
      CHARACTER*1  ACT
C
C  init
      NDIM = NPIX(1)
      DFACT = 1.D0 / DBLE(FLOAT(NDIM))
C
C  compute array MPOW:            2**(NLOG(1)-1), ..., 2**1, 2**0
      MM = NLOG(1) 
      DO 200, N=1,NLOG(1)
         MPOW(N) = POW2(MM)
         MM = MM - 1
200   CONTINUE
C
C  generate sine + cosine look-up table in x-dimension
      CALL DSCTB(TSIN,TCOS,NDIM,NLOG(1),SWAPI,SWAPO,SWALIM)
C
C  here we handle 1-dim FFT
C
      IF (NPIX(2).EQ.1) THEN
         IF (ACT.NE.'F') THEN
            DO 400, NX=1,NDIM
               DREAL(NX) = DBLE(A(NX))                  !copy A to DREAL
               DIMAG(NX) = DBLE(B(NX))                  !copy B to DIMAG
400         CONTINUE
C
            CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(1),MPOW,
     +               TSIN,TCOS,DFACT,C,D,SWAPI,SWAPO,SWALIM)        !do the FFT
         ELSE
C
C  here for frequency business (1-dim case)
C
            MM = 0
            DO 1000, NX=1,NDIM
               DX = DBLE(A(NX))
               DY = DBLE(B(NX))
               IF (MM.EQ.0) THEN
                  DREAL(NX) = DX                 !copy A*(-1)**(NX+NY) to DREAL
                  DIMAG(NX) = DY                 !copy B*(-1)**(NX+NY) to DIMAG
                  MM = 1
               ELSE
                  DREAL(NX) = -DX                !copy A*(-1)**(NX+NY) to DREAL
                  DIMAG(NX) = -DY                !copy B*(-1)**(NX+NY) to DIMAG
                  MM = 0
               ENDIF
1000        CONTINUE
C
            CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(1),MPOW,
     +               TSIN,TCOS,DFACT,C,D,SWAPI,SWAPO,SWALIM)         !do the FFT
         ENDIF
         RETURN
      ENDIF
C
C
C  -------------- here we handle 2-dim FFT -------------------
C
      NSIZE(1) = 128
      NSIZE(2) = 256
C
C  branch, if frequency scaling is wanted
      XFROM = 1
      IF (ACT.EQ.'F') GOTO 6000
      NWORK = NDIM - 1
C
C  do FFT on each line
      DO 5800, NY=1,NPIX(2)
         N = 1
         DO 5400, NX = XFROM,XFROM+NWORK
            DREAL(N) = DBLE(A(NX))                  !copy A to DREAL
            DIMAG(N) = DBLE(B(NX))                  !copy B to DIMAG
            N = N + 1
5400     CONTINUE
C
         CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW,
     +            TSIN,TCOS,DFACT,C(XFROM),D(XFROM),
     +            SWAPI,SWAPO,SWALIM)                    !do the FFT
C
         XFROM = XFROM + NDIM
5800  CONTINUE
C
      GOTO 7000
C
C  here for frequency business (2-dim case)
6000  KMM = 0
      DO 6800, NY=1,NPIX(2)
         N = 1
         MM = KMM
         NWORK = XFROM - 1
C
         DO 6400, NX = 1,NDIM
            INDX = NWORK + NX
            DX = DBLE(A(INDX))
            DY = DBLE(B(INDX))
            IF (MM.EQ.0) THEN
               DREAL(N) = DX                  !copy A*(-1)**(NX+NY) to DREAL
               DIMAG(N) = DY                  !copy B*(-1)**(NX+NY) to DIMAG
               MM = 1
            ELSE
               DREAL(N) = -DX                 !copy A*(-1)**(NX+NY) to DREAL
               DIMAG(N) = -DY                 !copy B*(-1)**(NX+NY) to DIMAG
               MM = 0
            ENDIF
            N = N + 1
6400     CONTINUE
C
         CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW,
     +            TSIN,TCOS,DFACT,C(XFROM),D(XFROM),
     +            SWAPI,SWAPO,SWALIM)                    !do the FFT

C
         XFROM = XFROM + NDIM
         KMM = 1 - KMM                        !alternate between 0 and 1
6800  CONTINUE
C
C  we now do another FFT along the columns
7000  CALL LINCOL(C,NPIX,NSIZE,CT)            !change lines -> columns
      CALL LINCOL(D,NPIX,NSIZE,DT)
      NDIM = NPIX(2)
      DFACT = DFACT / DBLE(FLOAT(NDIM))
C
C  compute array MPOW in y-dim:            2**(NLOG(2)-1), ..., 2**1, 2**0
      IF (NLOG(2).NE.NLOG(1)) THEN
         MM = NLOG(2)
         DO 7400, N=1,NLOG(2)
            MPOW(N) = POW2(MM)
            MM = MM - 1
7400     CONTINUE
C
C  and generate sine + cosine look-up table in y-dimension
         CALL DSCTB(TSIN,TCOS,NDIM,NLOG(2),SWAPI,SWAPO,SWALIM) 
      ENDIF
C
      XFROM = 1
      NWORK = NDIM - 1
C
      DO 8000, NY=1,NPIX(1)
         N = 1
         DO 7700, NX=XFROM,XFROM+NWORK
            DREAL(N) = DBLE(CT(NX))
            DIMAG(N) = DBLE(DT(NX))
            N = N + 1
7700     CONTINUE
C
         CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(2),MPOW,
     +            TSIN,TCOS,DFACT,CT(XFROM),DT(XFROM),
     +            SWAPI,SWAPO,SWALIM)                    !do the FFT
C
         XFROM = XFROM + NDIM
8000  CONTINUE
C
C  finally transpose back again
      CALL TRAPO(CT,DT,NPIX,NSIZE,C,D)
C
      RETURN
      END

      SUBROUTINE DFFTB(A,B,C,D,NPIX,ACT,DREAL,DIMAG,NLOG,POW2,
     +                 TSIN,TCOS,CT,DT,SWAPI,SWAPO)
C
      IMPLICIT NONE
C
      INTEGER    N,NDIM,NX,NY,MM
      INTEGER    XFROM,NSIZE(2),NWORK
      INTEGER    NPIX(2),NLOG(*),POW2(*)
      INTEGER    KMM,MPOW(21)
      INTEGER    SWAPI(*),SWAPO(*),SWALIM
C
      REAL      A(*),B(*),C(*),D(*),CT(*),DT(*)
C
      DOUBLE PRECISION      DREAL(*),DIMAG(*),TSIN(*),TCOS(*)
      DOUBLE PRECISION      DFACT
C
      CHARACTER*1      ACT
C
C  init
      NDIM = NPIX(1)
      DFACT = 1.D0 / DBLE(FLOAT(NDIM))
C
C  compute array MPOW:            2**(NLOG(1)-1), ..., 2**1, 2**0
      MM = NLOG(1)
      DO 200, N=1,NLOG(1)
         MPOW(N) = POW2(MM)
         MM = MM - 1
200   CONTINUE
C
C  generate sine + cosine look-up table in x-dimension
      CALL DSCTB(TSIN,TCOS,NDIM,NLOG(1),SWAPI,SWAPO,SWALIM)
C
C  do 1-dim inverse FFT
      IF (NPIX(2).EQ.1) THEN
         IF (ACT.NE.'F') THEN
            DO 400, NX=1,NDIM
               DREAL(NX) = DBLE(A(NX))            !copy A to DREAL
               DIMAG(NX) = DBLE(-B(NX))           !copy -B to DIMAG
400         CONTINUE
C
            CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW,
     +               TSIN,TCOS,DFACT,C,D,SWAPI,SWAPO,SWALIM)        !do the FFT
         ELSE
C
C  frequency business...
            DO 600, NX=1,NDIM
               DREAL(NX) = DBLE(A(NX))            !copy A to DREAL
               DIMAG(NX) = DBLE(-B(NX))           !copy -B to DIMAG
600         CONTINUE
C
            CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW,
     +               TSIN,TCOS,DFACT,C,D,SWAPI,SWAPO,SWALIM)        !do the FFT
C 
            MM = 0
            DO 1000, NX=1,NDIM
               IF (MM.EQ.0) THEN
                  MM = 1
               ELSE
                  C(NX) = - C(NX)
                  D(NX) = - D(NX)
                  MM = 0
               ENDIF
1000        CONTINUE
         ENDIF
C
         RETURN
      ENDIF
C
C
C  ------------------------- do 2-dim inverse FFT -------------------------
C
      NSIZE(1) = 128
      NSIZE(2) = 256
C
C  handle each line
      XFROM = 1
      NWORK = NDIM - 1
C
      DO 5800, NY=1,NPIX(2)
         N = 1
         DO 5400, NX=XFROM,XFROM+NWORK
            DREAL(N) = DBLE(A(NX))                  !copy A to DREAL
            DIMAG(N) = DBLE(-B(NX))                 !copy -B to DIMAG
            N = N + 1
5400     CONTINUE
C
         CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW,
     +            TSIN,TCOS,DFACT,C(XFROM),D(XFROM),
     +            SWAPI,SWAPO,SWALIM)                    !do the FFT
         XFROM = XFROM + NDIM
5800  CONTINUE
C
C  handle each column
      CALL LINCOL(C,NPIX,NSIZE,CT)
      CALL LINCOL(D,NPIX,NSIZE,DT)
      NDIM = NPIX(2)
CC      DFACT = DFACT * DBLE(FLOAT(NDIM))
      DFACT = 1.D0
C
C  compute array MPOW in y-dim:            2**(NLOG(2)-1), ..., 2**1, 2**0
      IF (NLOG(2).NE.NLOG(1)) THEN
         MM = NLOG(2)
         DO 6000, N=1,NLOG(2)
            MPOW(N) = POW2(MM)
            MM = MM - 1
6000     CONTINUE
C
C  and generate sine + cosine look-up table in y-dimension
         CALL DSCTB(TSIN,TCOS,NDIM,NLOG(2),SWAPI,SWAPO,SWALIM) 
      ENDIF
C
      XFROM = 1
      NWORK = NDIM - 1
C
      IF (ACT.NE.'F') THEN
         DO 6800, NY=1,NPIX(1)
            N = 1
            DO 6400, NX=XFROM,XFROM+NWORK
               DREAL(N) = DBLE(CT(NX))
               DIMAG(N) = DBLE(DT(NX))
               N = N + 1
6400        CONTINUE
C
            CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(2),MPOW,
     +               TSIN,TCOS,DFACT,CT(XFROM),DT(XFROM),
     +               SWAPI,SWAPO,SWALIM)                    !do the FFT
            XFROM = XFROM + NDIM
6800     CONTINUE
C
      ELSE
C
C  frequency business...
         KMM = 0
         DO 7800, NY=1,NPIX(1)
            N = 1
            DO 7400, NX=XFROM,XFROM+NWORK
               DREAL(N) = DBLE(CT(NX))
               DIMAG(N) = DBLE(DT(NX))
               N = N + 1
7400        CONTINUE
C
            CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(2),MPOW,
     +               TSIN,TCOS,DFACT,CT(XFROM),DT(XFROM),
     +               SWAPI,SWAPO,SWALIM)                    !do the FFT
C
            MM = KMM
            DO 7600, NX=XFROM,XFROM+NWORK
               IF (MM.EQ.0) THEN
                  MM = 1
               ELSE
                  CT(NX) = - CT(NX)
                  DT(NX) = - DT(NX)
                  MM = 0
               ENDIF
7600        CONTINUE
C
            XFROM = XFROM + NDIM
            KMM = 1 - KMM                     !alternate between 0 and 1
7800     CONTINUE
C
      ENDIF
C
C  finally transpose back again
      CALL TRAPO(CT,DT,NPIX,NSIZE,C,D)
C
      RETURN
      END

      SUBROUTINE DRFFTF(A,C,D,NPIX,ACT,DREAL,DIMAG,NLOG,POW2,
     +                  TSIN,TCOS,CT,DT,SWAPI,SWAPO)
C
      IMPLICIT NONE
C
      INTEGER    INDX,N,NDIM,NX,NY,MM
      INTEGER    XFROM,NSIZE(2),NWORK
      INTEGER    NPIX(2),NLOG(*),POW2(*)
      INTEGER    KMM,MPOW(21)
      INTEGER    SWAPI(*),SWAPO(*),SWALIM
C
      REAL      A(*),C(*),D(*),CT(*),DT(*)
C
      DOUBLE PRECISION      DREAL(*),DIMAG(*),TSIN(*),TCOS(*)
      DOUBLE PRECISION      DFACT
      DOUBLE PRECISION      DX
C
      CHARACTER*1  ACT
C
C  init
      NDIM = NPIX(1)
      DFACT = 1.D0 / DBLE(FLOAT(NDIM))
C
C  compute array MPOW:            2**(NLOG(1)-1), ..., 2**1, 2**0
      MM = NLOG(1)
      DO 200, N=1,NLOG(1)
         MPOW(N) = POW2(MM)
         MM = MM - 1
200   CONTINUE
C
C  generate sine + cosine look-up table in x-dimension
      CALL DSCTB(TSIN,TCOS,NDIM,NLOG(1),SWAPI,SWAPO,SWALIM)
C
C  here we handle 1-dim FFT
C
      IF (NPIX(2).EQ.1) THEN
         IF (ACT.NE.'F') THEN
C
            DO 400, NX=1,NDIM
               DREAL(NX) = DBLE(A(NX))                  !copy A to DREAL
               DIMAG(NX) = 0.D0
400         CONTINUE
C
            CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(1),MPOW,
     +               TSIN,TCOS,DFACT,C,D,SWAPI,SWAPO,SWALIM)        !do the FFT
         ELSE
C
C  here for frequency business
            MM = 0
            DO 1000, NX=1,NDIM
               DX = DBLE(A(NX))
               IF (MM.EQ.0) THEN
                  DREAL(NX) = DX                 !copy A*(-1)**(NX+NY) to DREAL
                  MM = 1
               ELSE
                  DREAL(NX) = -DX                !copy A*(-1)**(NX+NY) to DREAL
                  MM = 0
               ENDIF
               DIMAG(NX) = 0.D0
1000        CONTINUE
C
            CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(1),MPOW,
     +               TSIN,TCOS,DFACT,C,D,SWAPI,SWAPO,SWALIM)        !do the FFT
         ENDIF
C
         RETURN
      ENDIF
C
C
C  ---------------------- here we handle 2-dim FFT ---------------------
C
      NSIZE(1) = 128
      NSIZE(2) = 256
C
C  branch, if frequency scaling is wanted
      XFROM = 1
      IF (ACT.EQ.'F') GOTO 6000
      NWORK = NDIM - 1
C
C  do FFT on each line
      DO 5800, NY=1,NPIX(2)
         N = 1
         DO 5400, NX = XFROM,XFROM+NWORK
            DREAL(N) = DBLE(A(NX))                  !copy A to DREAL
            DIMAG(N) = 0.D0
            N = N + 1
5400     CONTINUE
C
         CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW,
     +            TSIN,TCOS,DFACT,C(XFROM),D(XFROM),
     +            SWAPI,SWAPO,SWALIM)                    !do the FFT
         XFROM = XFROM + NDIM
5800  CONTINUE
C
      GOTO 7000
C
C  here for frequency business
6000  KMM = 0
      DO 6800, NY=1,NPIX(2)
         N = 1
         MM = KMM
         NWORK = XFROM - 1
C
         DO 6400, NX = 1,NDIM
            INDX = NWORK + NX
            DX = DBLE(A(INDX))
            IF (MM.EQ.0) THEN
               DREAL(N) = DX                  !copy A*(-1)**(NX+NY) to DREAL
               MM = 1
            ELSE
               DREAL(N) = -DX                 !copy A*(-1)**(NX+NY) to DREAL
               MM = 0
            ENDIF
            DIMAG(N) = 0.D0
            N = N + 1
6400     CONTINUE
C
         CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW,
     +            TSIN,TCOS,DFACT,C(XFROM),D(XFROM),
     +            SWAPI,SWAPO,SWALIM)                    !do the FFT
C
         KMM = 1 - KMM
         XFROM = XFROM + NDIM
6800  CONTINUE
C
C  we now do another FFT along the columns
7000  CALL LINCOL(C,NPIX,NSIZE,CT)            !change lines -> columns
      CALL LINCOL(D,NPIX,NSIZE,DT)
      NDIM = NPIX(2)
      DFACT = DFACT / DBLE(FLOAT(NDIM))
C
C  compute array MPOW in y-dim:            2**(NLOG(2)-1), ..., 2**1, 2**0
      IF (NLOG(2).NE.NLOG(1)) THEN
         MM = NLOG(2)
         DO 7400, N=1,NLOG(2)
            MPOW(N) = POW2(MM)
            MM = MM - 1
7400     CONTINUE
C
C  and generate sine + cosine look-up table in y-dimension
         CALL DSCTB(TSIN,TCOS,NDIM,NLOG(2),SWAPI,SWAPO,SWALIM) 
      ENDIF
C
      XFROM = 1
      NWORK = NDIM - 1
C
      DO 8000, NY=1,NPIX(1)
         N = 1
         DO 7700, NX=XFROM,XFROM+NWORK
            DREAL(N) = DBLE(CT(NX))
            DIMAG(N) = DBLE(DT(NX))
            N = N + 1
7700     CONTINUE
C
         CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(2),MPOW,
     +            TSIN,TCOS,DFACT,CT(XFROM),DT(XFROM),
     +            SWAPI,SWAPO,SWALIM)                    !do the FFT
         XFROM = XFROM + NDIM
8000  CONTINUE
C
C  finally transpose back again
      CALL TRAPO(CT,DT,NPIX,NSIZE,C,D)
C
      RETURN
      END

      SUBROUTINE DRFFTB(A,C,D,NPIX,ACT,DREAL,DIMAG,NLOG,POW2,
     +                  TSIN,TCOS,CT,DT,SWAPI,SWAPO)
C
      IMPLICIT NONE
C
      INTEGER    N,NDIM,NX,NY,MM
      INTEGER    XFROM,NSIZE(2),NWORK
      INTEGER    NPIX(2),NLOG(*),POW2(*)
      INTEGER    KMM,MPOW(21)
      INTEGER    SWAPI(*),SWAPO(*),SWALIM
C
      REAL      A(*),C(*),D(*),CT(*),DT(*)
C
      DOUBLE PRECISION      DREAL(*),DIMAG(*),TSIN(*),TCOS(*)
      DOUBLE PRECISION      DFACT
C
      CHARACTER*1      ACT
C
C  init
      NDIM = NPIX(1)
      DFACT = 1.D0 / DBLE(FLOAT(NDIM))
C
C  compute array MPOW:            2**(NLOG(1)-1), ..., 2**1, 2**0
      MM = NLOG(1)
      DO 200, N=1,NLOG(1)
         MPOW(N) = POW2(MM)
         MM = MM - 1
200   CONTINUE
C
C  generate sine + cosine look-up table in x-dimension
      CALL DSCTB(TSIN,TCOS,NDIM,NLOG(1),SWAPI,SWAPO,SWALIM)
C
C  do 1-dim inverse FFT
      IF (NPIX(2).EQ.1) THEN
         IF (ACT.NE.'F') THEN
            DO 400, NX=1,NDIM
               DREAL(NX) = DBLE(A(NX))            !copy A to DREAL
               DIMAG(NX) = 0.D0
400         CONTINUE
C
            CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW,
     +               TSIN,TCOS,DFACT,C,D,SWAPI,SWAPO,SWALIM)        !do the FFT
         ELSE
C
C  frequency business...
            DO 600, NX=1,NDIM
               DREAL(NX) = DBLE(A(NX))            !copy A to DREAL
               DIMAG(NX) = 0.D0
600         CONTINUE
C
            CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW,
     +               TSIN,TCOS,DFACT,C,D,SWAPI,SWAPO,SWALIM)        !do the FFT
C
            MM = 0
            DO 1000, NX=1,NDIM
               IF (MM.EQ.0) THEN
                  MM = 1
               ELSE
                  C(NX) = - C(NX)
                  D(NX) = - D(NX)
                  MM = 0
               ENDIF
1000        CONTINUE
C
         ENDIF
         RETURN
      ENDIF
C
C
C  ------------------------- do 2-dim inverse FFT -------------------------
C
      NSIZE(1) = 128
      NSIZE(2) = 256
C
C  handle each line
      XFROM = 1
      NWORK = NDIM - 1
C
      DO 5800, NY=1,NPIX(2)
         N = 1
         DO 5400, NX=XFROM,XFROM+NWORK
            DREAL(N) = DBLE(A(NX))                  !copy A to DREAL
            DIMAG(N) = 0.D0
            N = N + 1
5400     CONTINUE
C
         CALL DFFT(0,DREAL,DIMAG,NDIM,NLOG(1),MPOW,
     +            TSIN,TCOS,DFACT,C(XFROM),D(XFROM),
     +            SWAPI,SWAPO,SWALIM)                    !do the FFT
C
         XFROM = XFROM + NDIM
5800  CONTINUE
C
C  handle each column
      CALL LINCOL(C,NPIX,NSIZE,CT)
      CALL LINCOL(D,NPIX,NSIZE,DT)
      NDIM = NPIX(2)
      DFACT = DFACT * DBLE(FLOAT(NDIM))
C
C  compute array MPOW in y-dim:            2**(NLOG(2)-1), ..., 2**1, 2**0
      IF (NLOG(2).NE.NLOG(1)) THEN
         MM = NLOG(2)
         DO 6000, N=1,NLOG(2)
            MPOW(N) = POW2(MM)
            MM = MM - 1
6000     CONTINUE
C
C  and generate sine + cosine look-up table in y-dimension
         CALL DSCTB(TSIN,TCOS,NDIM,NLOG(2),SWAPI,SWAPO,SWALIM) 
      ENDIF
C
      XFROM = 1
      NWORK = NDIM - 1
C
      IF (ACT.NE.'F') THEN
         DO 6800, NY=1,NPIX(1)
            N = 1
            DO 6400, NX=XFROM,XFROM+NWORK
               DREAL(N) = DBLE(CT(NX))
               DIMAG(N) = DBLE(DT(NX))
               N = N + 1
6400        CONTINUE
C
            CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(2),MPOW,
     +               TSIN,TCOS,DFACT,CT(XFROM),DT(XFROM),
     +               SWAPI,SWAPO,SWALIM)                    !do the FFT
C
            XFROM = XFROM + NDIM
6800     CONTINUE
C
      ELSE
C
C  frequency business...
         KMM = 0
         DO 7800, NY=1,NPIX(1)
            N = 1
            DO 7400, NX=XFROM,XFROM+NWORK
               DREAL(N) = DBLE(CT(NX))
               DIMAG(N) = DBLE(DT(NX))
               N = N + 1
7400        CONTINUE
C
            CALL DFFT(1,DREAL,DIMAG,NDIM,NLOG(2),MPOW,
     +               TSIN,TCOS,DFACT,CT(XFROM),DT(XFROM),
     +               SWAPI,SWAPO,SWALIM)                    !do the FFT
C
            MM = KMM
            DO 7600, NX=XFROM,XFROM+NWORK
               IF (MM.EQ.0) THEN
                  MM = 1
               ELSE
                  CT(NX) = - CT(NX)
                  DT(NX) = - DT(NX)
                  MM = 0
               ENDIF
7600        CONTINUE
C
            XFROM = XFROM + NDIM
            KMM = 1 - KMM
7800     CONTINUE
C
      ENDIF
C
C  finally transpose back again
      CALL TRAPO(CT,DT,NPIX,NSIZE,C,D)
C
      RETURN
      END

      SUBROUTINE DFFT(IFACT,XREAL,XIMAG,N,NU,MPOW,TSIN,TCOS,DFACT,
     +                XR,XI,SWAPI,SWAPO,SWLIM)
C
C  FAST FOURIER TRANSFORM ALGORITHM
C  K. Banse        version 2.10      850625
C                          2.20      860217
C                          2.40      870109
C                          2.50      921022
C
C  IFACT - factor flag:  1 - yes, multiply by DFACT
C                     0 - no, do not multiply
C  XREAL - double precision real part of input array
C  XIMAG - double precision imaginary part of input array
C  N     - dimension of input array
C  NU    - LOG(N) i.e.: N  =  2**NU
C  MPOW  - array with powers of 2 from 1 - N
C  TSIN  - array precomputed sine values
C  TCOS  - array precomputed cosine values
C  DFACT - double precision constant for multiplication (IFACT = 1)
C  XR    - real part of XREAL
C  XI    - imag part of XIMAG
C
      IMPLICIT NONE
C
      INTEGER    IFACT,N,NU
      INTEGER    N1,MPL
      INTEGER    I,K,K1,K1N2,L
      INTEGER    N2,LP
      INTEGER    MPOW(*)
      INTEGER    SWAPI(*),SWAPO(*),SWLIM
C
      DOUBLE PRECISION      XREAL(*),XIMAG(*),TSIN(*),TCOS(*),DFACT
      DOUBLE PRECISION      TREAL,TIMAG,TCC,TSS
C
      REAL      XR(*),XI(*)
      REAL      TR
C
C  initialize
      N2 = MPOW(1)                        ! the same as: N2 = N/2
      N1 = N - 1
C
C  handle vertical arrays
      DO 500, L=1,NU
         K = 0
         MPL = MPOW(L)
C
50       DO 200, I=1,N2
            LP = K/MPL + 1
            TCC = TCOS(LP)
            TSS = TSIN(LP)
            K = K + 1
            K1N2 = K + N2
            TREAL = XREAL(K1N2)*TCC + XIMAG(K1N2)*TSS
            TIMAG = XIMAG(K1N2)*TCC - XREAL(K1N2)*TSS
            XREAL(K1N2) = XREAL(K) - TREAL
            XIMAG(K1N2) = XIMAG(K) - TIMAG
            XREAL(K) = XREAL(K) + TREAL
            XIMAG(K) = XIMAG(K) + TIMAG
200      CONTINUE
C
         K = K + N2
         IF (K.LT.N1) GOTO 50
C
         N2 = N2/2                  !set up variables for next array
500   CONTINUE
C
C  convert from double precision to single + multiply (if IFACT = 1)
      IF (IFACT.EQ.1) THEN
         DO 1000, K=1,N
            XR(K) = SNGL(XREAL(K) * DFACT)
            XI(K) = SNGL(XIMAG(K) * DFACT)
1000     CONTINUE
      ELSE
         DO 1200, K=1,N
            XR(K) = SNGL(XREAL(K))
            XI(K) = SNGL(XIMAG(K))
1200     CONTINUE
      ENDIF
C
C  unscramble final vector
      DO 2000, K=1,SWLIM 
         K1 = SWAPI(K)
         I = SWAPO(K)
         TR = XR(K1)
         XR(K1) = XR(I)
         XR(I) = TR
         TR = XI(K1)
         XI(K1) = XI(I)
         XI(I) = TR
2000  CONTINUE
C
C  that's it folks...
      RETURN
      END

      SUBROUTINE DSCTB(TSIN,TCOS,NDIM,NU,TIN,TOUT,SWLIM)
C
C subroutine to generate sine and cosine look-up tables
C and also to generate swap tables
C
      IMPLICIT NONE
C
      INTEGER    IBITR,NDIM,NU,K,M,IBR
      INTEGER    TIN(NDIM),TOUT(NDIM),SWLIM
C
      DOUBLE PRECISION      TSIN(*),TCOS(*),ARG,ARG1
C
      ARG1 = 6.2831852D0 / DBLE(FLOAT(NDIM))
      SWLIM = 0
C
      DO 1000, K=1,NDIM
         M = K - 1
         IBR = IBITR(M,NU)
         ARG = ARG1 * IBR
         TSIN(K) = SIN(ARG)
         TCOS(K) = COS(ARG)
C
         IF (IBR.GT.M) THEN                !fill swap tables already
            SWLIM = SWLIM + 1
            TIN(SWLIM) = K 
            TOUT(SWLIM) = IBR + 1
         ENDIF
1000  CONTINUE
C
      RETURN
      END

      SUBROUTINE TRAPO(A,B,NPIX,NSIZE,C,D)
C      
      IMPLICIT NONE
C      
      INTEGER    NPIX(2),NSIZE(2)
      INTEGER    NPAX(2)
C      
      REAL       A(*),B(*),C(*),D(*)
C      
      NPAX(1) = NPIX(2)
      NPAX(2) = NPIX(1)
      CALL LINCOL(A,NPAX,NSIZE,C)
      CALL LINCOL(B,NPAX,NSIZE,D)
C      
      RETURN
      END

      FUNCTION IBITR(J,NU)
C      
      IMPLICIT NONE
C      
      INTEGER    IBITR,J,NU
      INTEGER    I,J1,J2
C      
      J1 = J
      IBITR = 0
      DO 1000 I=1,NU
         J2 = J1/2
         IBITR = IBITR*2 + (J1-2*J2)
         J1 = J2
1000  CONTINUE
C      
      RETURN
      END

      SUBROUTINE POWER(A,B,C,NDIM)
C      
      IMPLICIT NONE
C      
      INTEGER    N,NDIM
C      
      REAL       A(NDIM),B(NDIM),C(NDIM)
C      
      DO 1000 N=1,NDIM
         C(N) = SQRT(  ( A(N)*A(N) )  + ( B(N)*B(N) )  )
1000  CONTINUE
C      
      RETURN
      END
