C @(#)old_fftetc.for	19.1 (ES0-DMD) 02/25/03 13:25:06
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
      SUBROUTINE DFILL(DATA,NX,NY,OUT)
C
C Fill up an array in the way required by the Numerical Recipes
C FFT routines
C
C Note this routine has a 'fudge factor' of 1 pixel to align coordinate
C systems correctly.
C
      IMPLICIT NONE
                   
      INTEGER NX,NY,IP,I,J
      DOUBLE PRECISION DATA(NX,NY)
      DOUBLE PRECISION OUT(NX*NY*2)
                                   
      IP=1
          
      IF(NY.GT.1) THEN
         DO I=1,NX
            OUT(IP)=0.0D0
            IP=IP+2
         ENDDO
              
         DO J=1,NY-1
            OUT(IP)=0.0D0
            IP=IP+2
            DO I=1,NX-1
               OUT(IP-1)=0.0D0
               OUT(IP)=DATA(I,J)
               IP=IP+2
            ENDDO
          ENDDO
       ELSE
         DO I=1,NX
            OUT(IP)=DATA(I,1)
            OUT(IP+1)=0.0D0
            IP=IP+2
         ENDDO
      ENDIF
           
      RETURN
      END
         
      LOGICAL FUNCTION POWER2(IV)
                                 
      IMPLICIT NONE
                   
      INTEGER IV
      REAL V
            
      V=FLOAT(IV)
      V=ALOG10(V)/ALOG10(2.0)
                             
      IF(ABS(V-FLOAT(NINT(V))).LT.1.0E-5) THEN
         POWER2=.TRUE.
      ELSE
         POWER2=.FALSE.
      ENDIF
           
      RETURN
      END
         
      SUBROUTINE FILL(DATA,NX,NY,OUT)
C
C Fill up an array in the way required by the Numerical Recipes
C FFT routines
C
C Single precision version
C
C Note this routine has a 'fudge factor' of 1 pixel to align coordinate
C systems correctly.
C
      IMPLICIT NONE
                   
      INTEGER NX,NY,IP,I,J
      REAL DATA(NX,NY)
      REAL OUT(NX*NY*2)
                       
      IP=1
          
      IF(NY.GT.1) THEN
         DO I=1,NX
            OUT(IP)=0.0
            IP=IP+2
         ENDDO
              
         DO J=1,NY-1
            OUT(IP)=0.00
            IP=IP+2
            DO I=1,NX-1
               OUT(IP-1)=0.0
               OUT(IP)=DATA(I,J)
               IP=IP+2
            ENDDO
          ENDDO
       ELSE
         DO I=1,NX
            OUT(IP)=DATA(I,1)
            OUT(IP+1)=0.0
            IP=IP+2
         ENDDO
      ENDIF
           
      RETURN
      END
         
      SUBROUTINE DCONV(DATA,NX,NY,WORK,PSFFFT,OUT,FLAG)
C
C Convolve a one or two dimensional REAL array with a PSF which is supplied as
C a double precision FFT. Uses Numerical Recipes
C
C Richard Hook, ST-ECF, October 1991
C Double precision version, March 1992
C 1D version, August 1992
C
      IMPLICIT NONE
                   
      INTEGER I,J,N,IP,IQ
      DOUBLE PRECISION A,B,C,D,FAC
      INTEGER NX,NY                ! ARRAY SIZE
      INTEGER NDIMS,DIMS(2),NX2,NY2
      INTEGER FLAG                 ! ROTATED OR NOT?
      REAL*8 DATA(NX,NY)             ! INPUT DATA ARRAY
      REAL*8 OUT(NX,NY)
      REAL*8 PSFFFT(NX*NY*2)         ! FFT OF PSF
      REAL*8 WORK(NX*NY*2)           ! WORKSPACE FOR FFT OF DATA
      REAL*8 DFLAG
                  
      DFLAG=DBLE(FLAG)
      IF(NY.EQ.1) THEN
         NDIMS=1
      ELSE
         NDIMS=2
      ENDIF
           
C Copy the data into the real part of the working array
      IP=1
      DO J=1,NY
         DO I=1,NX
            WORK(IP)=DATA(I,J)
            WORK(IP+1)=0.0D0
            IP=IP+2
         ENDDO
      ENDDO
           
C Take the forward FFT of the data array
      DIMS(1)=NX
      DIMS(2)=NY
                
      CALL DFOURN(WORK,DIMS,NDIMS,+1)
                                     
C Do a complex multiplication of the FFT of the data and of the PSF
C
C If the rotation flag is set (to -1) we multiply the imaginary part
C of the FFT by -1
      IP=1
      IQ=2
      DO N=1,NX*NY
         A=WORK(IP)
         B=WORK(IQ)
         C=PSFFFT(IP)
         D=DFLAG*PSFFFT(IQ)
         WORK(IP)=A*C-B*D
         WORK(IQ)=A*D+B*C
         IP=IP+2
         IQ=IQ+2
      ENDDO
           
C Take the inverse FFT
      CALL DFOURN(WORK,DIMS,NDIMS,-1)
                                     
C Copy the result into the output array
C and 'quadrant swap':
      FAC=DFLOAT(NX*NY)
      IP=1
          
      NX2=MAX(1,NX/2)
      NY2=MAX(1,NY/2)
                     
      IF(NY.GT.1) THEN
       DO J=1,NY2
          DO I=1,NX2
             OUT(I+NX/2,J+NY/2)=WORK(IP)/FAC
             IP=IP+2
          ENDDO
          IP=IP+NX
       ENDDO
            
       IP=NX+1
       DO J=1,NY2
          DO I=1,NX2
             OUT(I,J+NY/2)=WORK(IP)/FAC
             IP=IP+2
          ENDDO
          IP=IP+NX
       ENDDO
            
       IP=NX*NY+1
       DO J=1,NY2
          DO I=1,NX2
             OUT(I+NX/2,J)=WORK(IP)/FAC
             IP=IP+2
          ENDDO
          IP=IP+NX
       ENDDO
            
       IP=NX*NY+NX+1
       DO J=1,NY2
          DO I=1,NX2
             OUT(I,J)=WORK(IP)/FAC
             IP=IP+2
          ENDDO
          IP=IP+NX
       ENDDO
      ELSE
          
C 1D case
       DO I=1,NX2
         OUT(I+NX2,1)=WORK(IP)/FAC
         IP=IP+2
       ENDDO
            
       DO I=1,NX2
         OUT(I,1)=WORK(IP)/FAC
         IP=IP+2
       ENDDO
      ENDIF
           
      RETURN
      END
      SUBROUTINE DFOURN(DATA,NN,NDIM,ISIGN)
C
C Numerical Recipes FFT routine - modified to use entirely
C DOUBLE PRECISION throughout.
C
      IMPLICIT NONE
C
      REAL*8 WR,WI,WPR,WPI,WTEMP,THETA,TEMPR,TEMPI
      INTEGER NN(NDIM),NTOT,IDIM,NDIM,NPREV,NREM
      INTEGER IP1,IP2,IP3,I2REV,I3REV,N
      INTEGER I1,I2,I3,K1,K2,ISIGN,IBIT,IFP1,IFP2
      DOUBLE PRECISION DATA(*)
      NTOT=1
      DO 11 IDIM=1,NDIM
        NTOT=NTOT*NN(IDIM)
11    CONTINUE
      NPREV=1
      DO 18 IDIM=1,NDIM
        N=NN(IDIM)
        NREM=NTOT/(N*NPREV)
        IP1=2*NPREV
        IP2=IP1*N
        IP3=IP2*NREM
        I2REV=1
        DO 14 I2=1,IP2,IP1
          IF(I2.LT.I2REV)THEN
            DO 13 I1=I2,I2+IP1-2,2
              DO 12 I3=I1,IP3,IP2
                I3REV=I2REV+I3-I2
                TEMPR=DATA(I3)
                TEMPI=DATA(I3+1)
                DATA(I3)=DATA(I3REV)
                DATA(I3+1)=DATA(I3REV+1)
                DATA(I3REV)=TEMPR
                DATA(I3REV+1)=TEMPI
12            CONTINUE
13          CONTINUE
          ENDIF
          IBIT=IP2/2
1         IF ((IBIT.GE.IP1).AND.(I2REV.GT.IBIT)) THEN
            I2REV=I2REV-IBIT
            IBIT=IBIT/2
          GO TO 1
          ENDIF
          I2REV=I2REV+IBIT
14      CONTINUE
        IFP1=IP1
2       IF(IFP1.LT.IP2)THEN
          IFP2=2*IFP1
          THETA=ISIGN*6.28318530717959D0/(IFP2/IP1)
          WPR=-2.D0*DSIN(0.5D0*THETA)**2
          WPI=DSIN(THETA)
          WR=1.D0
          WI=0.D0
          DO 17 I3=1,IFP1,IP1
            DO 16 I1=I3,I3+IP1-2,2
              DO 15 I2=I1,IP3,IFP2
                K1=I2
                K2=K1+IFP1
                TEMPR=WR*DATA(K2)-WI*DATA(K2+1)
                TEMPI=WR*DATA(K2+1)+WI*DATA(K2)
                DATA(K2)=DATA(K1)-TEMPR
                DATA(K2+1)=DATA(K1+1)-TEMPI
                DATA(K1)=DATA(K1)+TEMPR
                DATA(K1+1)=DATA(K1+1)+TEMPI
15            CONTINUE
16          CONTINUE
            WTEMP=WR
            WR=WR*WPR-WI*WPI+WR
            WI=WI*WPR+WTEMP*WPI+WI
17        CONTINUE
          IFP1=IFP2
        GO TO 2
        ENDIF
        NPREV=N*NPREV
18    CONTINUE
      RETURN
      END
