C @(#)old_fftetc.for 17.1.1.1 (ES0-DMD) 01/25/02 17:15:07 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