SUBROUTINE WAVCRE(N,PST,LAMST,DELAM,WAV) C C This subroutine creates the wavelength array C WAV(N) from the reference pixel PST and the C start wavelength LAMST and the wavelength step C DELAM. If LAMST is not set (no wavelength C information) then the wavelength is set to C (real) channel numbers C IMPLICIT NONE INTEGER N REAL PST REAL LAMST REAL DELAM DOUBLE PRECISION WAV(N) C C Local variables C INTEGER I IF (LAMST.LE.1.0) THEN DO I=1,N,1 WAV(I)=DBLE(I) ! IRAF convention ENDDO ELSE DO I=1,N,1 WAV(I)=DBLE(LAMST + (REAL(I)-PST)*DELAM) ENDDO ENDIF 99 END SUBROUTINE POILIM(NW,CENCX,CENCY,CENM,CENT,N1,N2,NL,NS1) C C This subroutine returns the number of point sources NL C within the Y pixel range 2 - N2-1 for the NW profiles for C all X channels 1 - N1. The position of each point source is C calculated from the reference X,Y centre (CENCX,CENCY) and C the slope with respect to the X axis CENM. The first point C source above Y pixel 2 is NS1 C IMPLICIT NONE INTEGER NW ! Number of elements in input arrays INTEGER N1 ! Number of X channels for profile centre calculation INTEGER N2 ! Number of Y channels to contain point sources INTEGER NL ! Number of point soures in Y channel range 2 - N2-1 INTEGER NS1 ! Index number of first point source in input arrays REAL CENCX(NW) ! Input array of reference X value of profiles REAL CENCY(NW) ! Input array of reference Y value of profiles REAL CENM(NW) ! Input array of Y v. X slope of profile centres REAL CENT(NW) ! Centres of profiles for channel I C C Local variables C INTEGER I,J INTEGER ILIM INTEGER FLIM LOGICAL LSET LOGICAL LSET1 C C Compute the Y centre of the profiles for channel I. C Find the first index number of the first point source above C Y=2 for all X channels. Find the index number of the point C source at the highest Y channel just contained within the C range 2 - N2-1 for all X channels C LSET=.FALSE. LSET1=.FALSE. FLIM=0 DO J=1,NW,1 DO I=1,N1,1 CALL POICENI(NW,CENCX,CENCY,CENM,I,CENT) IF (INT(CENT(J)).GT.2.AND.INT(CENT(J)).LT.N2-1) THEN ILIM=J LSET=.TRUE. ELSE LSET=.FALSE. GO TO 30 ENDIF ENDDO 30 IF (LSET) THEN FLIM=ILIM IF (.NOT.LSET1) THEN NS1=ILIM LSET1=.TRUE. ENDIF ENDIF ENDDO NL=FLIM - NS1 + 1 99 END SUBROUTINE TCOPUT(VAL,N,IC,CEN) C C This subroutine places the value VAL into the IC'th element C of the array CEN C IMPLICIT NONE INTEGER N INTEGER IC REAL VAL REAL CEN(N) CEN(IC)=VAL 99 END SUBROUTINE OUT2T1(N1,N2,IS,SPEC12,SPEC11,SPEC22,SPEC21, : ISPEC2,ISPEC1) C C This subroutine copies the IS'th row of the 2-D arrays C SPEC12, SPEC22, ISPEC2 to the 1-D arrays SPEC11, SPEC21 and C ISPEC1 respectively C IMPLICIT NONE INTEGER N1 INTEGER N2 INTEGER IS INTEGER ISPEC2(N1,N2) INTEGER ISPEC1(N1) REAL SPEC12(N1,N2) REAL SPEC22(N1,N2) REAL SPEC11(N1) REAL SPEC21(N1) C INTEGER J DO J=1,N1,1 SPEC11(J)=SPEC12(J,IS) SPEC21(J)=SPEC22(J,IS) ISPEC1(J)=ISPEC2(J,IS) ENDDO 99 END SUBROUTINE LISNAM(ROOT,N,ONAME) C C This subroutine appends _N.tab to the name ROOT and C returns it as ONAME C IMPLICIT NONE INTEGER N CHARACTER*24 ROOT CHARACTER*32 ONAME C C Local variables C INTEGER ILENR INTEGER ILENN CHARACTER*24 CNUM ILENR=INDEX(ROOT,' ') -1 IF (N.LT.10) THEN ILENN=1 CALL NAMNUM(N,ILENN,CNUM) ENDIF IF (N.GE.10.AND.N.LT.100) THEN ILENN=2 CALL NAMNUM(N,ILENN,CNUM) ENDIF IF (N.GE.100.AND.N.LT.1000) THEN ILENN=3 CALL NAMNUM(N,ILENN,CNUM) ENDIF IF (N.GE.1000.AND.N.LT.10000) THEN ILENN=4 CALL NAMNUM(N,ILENN,CNUM) ENDIF ONAME=ROOT(:ILENR)//'_'//CNUM(1:ILENN)//'.tab' 99 END SUBROUTINE NAMNUM(N,ILEN,CNUM) C C This subroutine returns the character value of C N, which has ILEN digits, as CNUM C INTEGER N INTEGER ILEN CHARACTER*24 CNUM C C Local variables C INTEGER N1 INTEGER N2 INTEGER N3 INTEGER N4 CHARACTER*1 CHANUM IF (ILEN.EQ.1) THEN N1=N CNUM(1:2)=CHANUM(N1) ENDIF IF (ILEN.EQ.2) THEN N1=INT(N/10) N2=N - N1*10 CNUM(1:2)=CHANUM(N1) CNUM(2:3)=CHANUM(N2) ENDIF IF (ILEN.EQ.3) THEN N1=INT(N/100) N2=INT((N - N1*100)/10) N3=N - N1*100 - N2*10 CNUM(1:2)=CHANUM(N1) CNUM(2:3)=CHANUM(N2) CNUM(3:4)=CHANUM(N3) ENDIF IF (ILEN.EQ.4) THEN N1=INT(N/1000) N2=INT((N - N1*1000)/100) N3=INT((N - N1*1000 - N2*100)/10) N4=N - N1*1000 - N2*100 - N3*10 CNUM(1:2)=CHANUM(N1) CNUM(2:3)=CHANUM(N2) CNUM(3:4)=CHANUM(N3) CNUM(4:5)=CHANUM(N4) ENDIF 99 END CHARACTER FUNCTION CHANUM(N) C C This function subprogram returns the character string C of the integer number N C INTEGER N IF (N.EQ.0) THEN CHANUM='0' ENDIF IF (N.EQ.1) THEN CHANUM='1' ENDIF IF (N.EQ.2) THEN CHANUM='2' ENDIF IF (N.EQ.3) THEN CHANUM='3' ENDIF IF (N.EQ.4) THEN CHANUM='4' ENDIF IF (N.EQ.5) THEN CHANUM='5' ENDIF IF (N.EQ.6) THEN CHANUM='6' ENDIF IF (N.EQ.7) THEN CHANUM='7' ENDIF IF (N.EQ.8) THEN CHANUM='8' ENDIF IF (N.EQ.9) THEN CHANUM='9' ENDIF 99 END SUBROUTINE MKXYCEN(N1,N2,X,Y,XFIT,YFIT,XYCEN) C C This subroutine forms the array XYCEN which C specifies the real Y values of the lower edges of C the pixels at their XIN midpoints. C C XYCEN gives the mid point positions of lower edges C of pixels on rectified grid - i.e. positions C of asterisks in diagram. C This array defines distortion of grid for code. C C +--- ---+--- ---+--- ---+ C | | | | C pixel |i-1,j | i,j |i+1,j | Y C | | | | C +---*---+---*---+---*---+ C X C The IRAF convention is followed - viz. pixel 1 extends C from 0.5 to 1.5; pixel 2 frm 1.5 to 2.5 ... C For the last Y pixel the upper upper is at N2+0.5 C IMPLICIT NONE INTEGER N1 ! X dimension of output image XYCEN INTEGER N2 ! Y dimension of output image XYCEN DOUBLE PRECISION X(N1) ! Array of the X (wavelength) values of the pixels DOUBLE PRECISION Y(N2) ! Array of float values of centres of Y pixels DOUBLE PRECISION XFIT(N1,N2) ! Array of polynomial coefficients to generate the real X values of the X midpoints DOUBLE PRECISION YFIT(N1,N2) ! Array of polynomial coefficients to generate the real Y values of the X midpoints DOUBLE PRECISION XYCEN(N1,N2) ! Array of Y values of X mid-points of pixels C C Local variables C INTEGER I,J DOUBLE PRECISION CONST DOUBLE PRECISION DIC CONST=1.0D0 DIC=0.50D0 DO J=1,N2,1 DO I=1,N1,1 C XYCEN(I,J)=XFIT(I,J)*X(I) + YFIT(I,J)*Y(I) + DBLE(J-1)*CONST C : + DIC ! IRAF convention XYCEN(I,J)=DIC + DBLE(J-1)*CONST ENDDO ENDDO 99 END SUBROUTINE ZAPTOT(N1,N2,ARRY,DWK,LREP,LZERO) C C This subroutine checks for any negative values in the C array ARRY; if any are found they are set to zero. C It then finds the total in each of the N1 X-sections; C if the total in any X-section is not normalised to 1.0, C then the values in that X-section are renormalised C to 1.0. If any X-section has zero count then the flag C LZERO is set FALSE C IMPLICIT NONE INTEGER N1 ! 1st dimension of input array ARRY INTEGER N2 ! 2nd dimension of input array ARRY DOUBLE PRECISION ARRY(N1,N2) ! Input array to be zapped and renormalised DOUBLE PRECISION DWK(N2) ! Work array for holding X-sections of array ARRY LOGICAL LREP ! If TRUE then report no. negative values and totals per X-section LOGICAL LZERO ! If any X-sections of ARRY are <=0, then set FALSE C C Local variables C INTEGER I,J INTEGER NNEG INTEGER ISTAT INTEGER N11 DOUBLE PRECISION VMIN DOUBLE PRECISION PT CHARACTER*128 TEXT C C First check for negative values, if there are any set C them to zero C CALL ZAPNEG(ARRY,N1,N2,NNEG,VMIN,0.0D0) IF (NNEG.GT.0) THEN IF (LREP) THEN WRITE(TEXT, : '(''! Warning, PSF spectrum had '',I7, : '' negative values, set to zero.'')') NNEG CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF ENDIF C C Check that for each X-section of ARRY, the total C count is normalised to 1.0. If not normalise. C LZERO=.TRUE. N11=1 DO I=1,N1,1 DO J=1,N2,1 DWK(J)=ARRY(I,J) ENDDO CALL TOTAL(DWK,N2,N11,PT) IF (PT.LE.0.0D0) THEN IF (LREP) THEN WRITE(TEXT, : '(''! Warning, X-sect '',I6,''; PSF zero'')') I CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF LZERO=.FALSE. GO TO 50 ENDIF IF (DABS(PT-1.0D0).GT.1.0D-7) THEN CALL MULC(DWK,N2,N11,1.0D0/PT,DWK) DO J=1,N2,1 ARRY(I,J)=DWK(J) ENDDO IF (LREP) THEN WRITE(TEXT, : '(''! Warning, X-sect '',I6,''; PSF totals '', : F15.7,'' - renormalising.'')') I,PT CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF ENDIF 50 CONTINUE ENDDO 99 END SUBROUTINE ZAPNEG(DATA,NX,NY,NNEG,VMIN,VZAP) C C Find negative points in an array and set them to VZAP C Also return the number of negative points and the smallest C value. C IMPLICIT NONE INTEGER NX INTEGER NY INTEGER NNEG DOUBLE PRECISION DATA(NX,NY) DOUBLE PRECISION VMIN DOUBLE PRECISION VZAP C INTEGER I,J NNEG=0 VMIN=1.0D20 DO J=1,NY DO I=1,NX IF (DATA(I,J).LT.0.0D0) THEN NNEG=NNEG+1 IF (DATA(I,J).LT.VMIN) THEN VMIN=DATA(I,J) ENDIF DATA(I,J)=VZAP ENDIF ENDDO ENDDO 99 END SUBROUTINE TOTAL(DATA,NX,NY,TOT) C C Add up all elements in an image C IMPLICIT NONE INTEGER NX INTEGER NY DOUBLE PRECISION DATA(NX,NY) DOUBLE PRECISION TOT INTEGER I,J TOT=0.0D0 DO J=1,NY DO I=1,NX TOT=TOT+DATA(I,J) ENDDO ENDDO 99 END SUBROUTINE MULC(IN1,NX,NY,F,OUT) C C Just multiply one array by a constant C IMPLICIT NONE INTEGER NX INTEGER NY DOUBLE PRECISION IN1(NX,NY) DOUBLE PRECISION OUT(NX,NY) DOUBLE PRECISION F C C Local variables C INTEGER I,J DO J=1,NY DO I=1,NX OUT(I,J)=IN1(I,J)*F ENDDO ENDDO 99 END SUBROUTINE PSFPARN(P1,P2,PSF,N2,Y,WK,POS,PKT,WID) C C This subroutine returns rough parameters of the single C profile in each X-section of the array PSF. The position C of the peak (POS), the peak height (PKT) and the FWHM C (WID) are returned C IMPLICIT NONE INTEGER P1 ! 1st dimension of input array PSF INTEGER P2 ! 2nd dimension of input array PSF INTEGER N2 ! Dimension of input arrays Y and WK DOUBLE PRECISION PSF(P1,P2) ! Input array of N1 profiles to be fitted DOUBLE PRECISION Y(N2) ! Input array of real values of Y channel numbers for array PSF DOUBLE PRECISION WK(N2) ! Work array for 1-D sections of PSF DOUBLE PRECISION POS(P1) ! Output array of rough position of profile peaks DOUBLE PRECISION PKT(P1) ! Output array of rough peak height of profiles DOUBLE PRECISION WID(P1) ! Output array of rough FWHM of profiles C C Local variables C INTEGER I,J DO I=1,P1,1 DO J=1,P2,1 WK(J)=PSF(I,J) ENDDO CALL PROFPAR(N2,Y,WK,P2,POS(I),PKT(I),WID(I)) ENDDO 99 END SUBROUTINE PROFPAR(N,X,Y,M,POS,PKHT,WID) C C This subroutine returns a rough position of the peak C of the single profile in the array Y, of peak height C PKHT and FWHM WID C IMPLICIT NONE INTEGER N ! Dimension of iput arrays INTEGER M ! Extent of data in Y input array DOUBLE PRECISION X(N) ! Array of X pixel position DOUBLE PRECISION Y(N) ! Array of signal values (1-M significant) DOUBLE PRECISION POS ! Output position of found peak DOUBLE PRECISION PKHT ! Output height of found peak DOUBLE PRECISION WID ! Output width of found peak C C Local variables C INTEGER I INTEGER IPOS INTEGER LEFT INTEGER RIGHT PKHT=-1.0D30 IPOS=M/2 ! Guess on position of peak DO I=1,M,1 IF (Y(I).GT.PKHT) THEN PKHT=Y(I) IPOS=I ENDIF ENDDO POS=X(IPOS) C C Find pixels to left and right of peak position where peak C height drops to < 0.5 of peak C 30 IF (IPOS.GT.1) THEN DO I=IPOS-1,1,-1 IF (DABS(Y(I)).LT.0.50D0*DABS(PKHT)) THEN LEFT=I GO TO 40 ENDIF ENDDO ELSE LEFT=1 ENDIF 40 IF (IPOS.LT.M) THEN DO I=IPOS+1,M,1 IF (DABS(Y(I)).LT.0.50D0*DABS(PKHT)) THEN RIGHT=I GO TO 50 ENDIF ENDDO ELSE RIGHT=M ENDIF 50 WID=X(RIGHT)-X(LEFT)-1.0D0 IF (WID.LT.1.0D0) THEN WID=1.0D0 ENDIF 99 END SUBROUTINE PSFCENTN(P1,P2,PSF,N2,Y,POS,PKT,WID,METH,YVAL, : WHT,FIT,CEN,FWHM) C C This subroutine determines the position of the centre of C the profile and its FWHM for the P1 X-sections of input C array PSF at approximate position POS, either by simple C centroid of the peak (METH=C) or by fitting a Gaussian C (METH=G) C IMPLICIT NONE INTEGER P1 ! 1st dimension of input array PSF INTEGER P2 ! 2nd dimension of input array PSF INTEGER N2 ! Dimension of input array Y DOUBLE PRECISION PSF(P1,P2) ! Input array of P1 profiles to be fitted DOUBLE PRECISION Y(N2) ! Input array of real values of Y channel number for array PSF DOUBLE PRECISION POS(P1) ! Input array of rough position of profile peaks DOUBLE PRECISION PKT(P1) ! Input array of rough peak height of profiles DOUBLE PRECISION WID(P1) ! Input array of rough FWHM of profiles DOUBLE PRECISION YVAL(N2) ! Work array for 1-D sections of PSF DOUBLE PRECISION WHT(N2) ! Work array for weights of 1-D section of PSF DOUBLE PRECISION FIT(N2) ! Work array for fit of 1-D section of PSF DOUBLE PRECISION CEN(P1) ! Output array of fitted position of (P1) profiles DOUBLE PRECISION FWHM(P1) ! Output array of fitted FWHM of (P1) profiles CHARACTER*1 METH C C Local variables C INTEGER I,J DO I=1,P1,1 DO J=1,P2,1 YVAL(J)=PSF(I,J) ENDDO C C Take care to zeroise rest of YVAL if N2>P2 C IF (P2.LT.N2) THEN DO J=P2+1,N2,1 YVAL(J)=0.0D0 ENDDO ENDIF CALL PSFCENT(N2,Y,YVAL,P2,POS(I),PKT(I),WID(I),METH, : WHT,FIT,CEN(I),FWHM(I)) ENDDO 99 END SUBROUTINE PSFCENT(N,WAV,ARRY,M,POS,PK,WID,METH,WHT,FIT,CEN,FWHM) C C This subroutine determines the position of the centre of C the profile and its FWHM in the array ARRY at approximate C position POS either by simple centroid of the peak (METH=C) C or by fitting a Gaussian (METH=G) C IMPLICIT NONE INTEGER N ! Dimension of input arrays INTEGER M ! Extent of data in input array ARRY DOUBLE PRECISION WAV(N) ! Input array of X pixel position DOUBLE PRECISION ARRY(N) ! Input array of signal value DOUBLE PRECISION POS ! Input approximate position of profile centre DOUBLE PRECISION PK ! Input approximate height of profile peak DOUBLE PRECISION WID ! Input approximate width of profile DOUBLE PRECISION WHT(N) ! Holding array for weights in fitting DOUBLE PRECISION FIT(N) ! Holding array for fitted profile data DOUBLE PRECISION CEN ! Output improved position of profile centre DOUBLE PRECISION FWHM ! Output improved width of profile CHARACTER*1 METH IF (METH.EQ.'C'.OR.METH.EQ.'c') THEN CALL CENTCEN(N,WAV,ARRY,POS,CEN,FWHM) ENDIF IF (METH.EQ.'G'.OR.METH.EQ.'g') THEN CALL GAUSPOS(N,WAV,ARRY,WHT,FIT,POS,PK,WID,CEN,FWHM) ENDIF 99 END SUBROUTINE CENTCEN(N,X,Y,POS,CENT,FWHM) C C This subroutine determines the simple centroid of the C profile in the array Y at position POS as CENT C and its FWHM. The centroiding is performed over the C range above which the signal is non-zero, else the C whole array Y C IMPLICIT NONE INTEGER N DOUBLE PRECISION X(N) DOUBLE PRECISION Y(N) DOUBLE PRECISION POS DOUBLE PRECISION CENT DOUBLE PRECISION FWHM C C Local variables C INTEGER I INTEGER ILEFT INTEGER IRIGHT INTEGER IH DOUBLE PRECISION BACK DOUBLE PRECISION PK DOUBLE PRECISION SUMX DOUBLE PRECISION SUMY C C Determine the left and right limits about the peak where C the Y values go to zero. Set background as mean of values C at these two limits C ILEFT=1 IRIGHT=N DO I=NINT(POS)-1,1,-1 IF (Y(I).EQ.0.0D0) THEN ILEFT=I GO TO 20 ENDIF ENDDO 20 DO I=NINT(POS)+1,N,1 IF (Y(I).EQ.0.0D0) THEN IRIGHT=I GO TO 30 ENDIF ENDDO 30 BACK=0.50D0*(Y(ILEFT)+Y(IRIGHT)) C C Determine simple centroid C SUMX=0.0D0 SUMY=0.0D0 DO I=ILEFT,IRIGHT,1 SUMX=SUMX+X(I)*(Y(I)-BACK) SUMY=SUMY+(Y(I)-BACK) ENDDO IF (SUMY.NE.0.0D0) THEN CENT=SUMX/SUMY ELSE CENT=POS ENDIF C C Check that value of CENT is within valid limits C IF (CENT.LT.X(1).OR.CENT.GT.X(N)) THEN CENT=0.0D0 FWHM=0.0D0 GO TO 99 ENDIF IH=NINT(CENT) C C Determine an estimate of the FWHM from the half C height positions C PK=Y(IH)-BACK ILEFT=IH-1 IF (ILEFT.LT.1) THEN ILEFT=1 ENDIF IRIGHT=IH+1 IF (IRIGHT.GT.N) THEN IRIGHT=N ENDIF 50 DO I=IH,1,-2 IF ((Y(I)-BACK).LT.(0.5D0*PK).AND. : (Y(I-1)-BACK).LT.(0.5D0*PK)) THEN ILEFT=I GO TO 60 ENDIF ENDDO 60 DO I=IH,N,2 IF ((Y(I)-BACK).LT.(0.5D0*PK).AND. : (Y(I+1)-BACK).LT.(0.5D0*PK)) THEN IRIGHT=I GO TO 70 ENDIF ENDDO 70 FWHM=X(IRIGHT)-X(ILEFT) - 1.0D0 99 END SUBROUTINE GAUSPOS(N,X,Y,WHT,YFIT,POS,PK,WID,CENT,FWHM) C C This subroutine determines the position of the peak C of the profile Y as CENT by fitting a Gaussian to the C line. The Gaussian FWHM is also returned. The C continuum is assumed to be zero. C IMPLICIT NONE INTEGER N ! X dimension of input arrays DOUBLE PRECISION X(N) ! Input array of X pixel values DOUBLE PRECISION Y(N) ! Input array of signal values DOUBLE PRECISION WHT(N) ! Input array of Y array weights DOUBLE PRECISION YFIT(N) ! Holding array for fit to Y DOUBLE PRECISION POS ! Input approximate X position of peak of Y array DOUBLE PRECISION PK ! Input approximate height of peak in Y array DOUBLE PRECISION WID ! Input approximate FWHM of peak in Y array DOUBLE PRECISION CENT ! Output fitted position of X position of peak in Y array DOUBLE PRECISION FWHM ! Output fitted FWHM of peak in Y array C C Local variables C INTEGER I,K INTEGER NTERMS INTEGER NITER INTEGER IFAIL DOUBLE PRECISION FLAMDA DOUBLE PRECISION A(3) DOUBLE PRECISION CHISQ DOUBLE PRECISION OCHISQ DOUBLE PRECISION FRACHI C C Set the weights to 1.0 for unit weighting C DO I=1,N,1 WHT(I)=1.0D0 ENDDO C C Set up parameters for CURFITS to fit the line by a C Gaussian and the continuum by a straight line C using least squares with linearization of the fitting C function C NTERMS=3 FLAMDA=0.001D0 A(1)=PK A(2)=POS A(3)=WID/2.354820D0 C C Call the least squares fitting routine until the fractional C change in Chi-squared is less than some value (set to 1E-6) C NITER=10 ! Plenty of iterations of LS Gaussian fitting FRACHI=1.0D-6 OCHISQ=1.0D-7 ! Small starting value DO K=1,NITER,1 CALL CURFITG(N,X,Y,WHT,NTERMS,A,FLAMDA,YFIT,CHISQ,IFAIL) IF (DABS((CHISQ/OCHISQ)-1.0D0).LT.FRACHI) THEN GO TO 80 ELSE OCHISQ=CHISQ ENDIF ENDDO 80 CENT=A(2) FWHM=A(3)*2.354820D0 99 END SUBROUTINE CURFITG(NPTS,X,Y,WEIGHT,NTERMS,A, : FLAMDA,YFIT,CHISQR,IFAIL) C C Performs a least squares fit to a non-linear function with C a linearization of the fitting function. C The function is a single Gaussian on a zero continuum. C The program is adapted from Bevington, 1968, p.237-9. C C Subprograms called: C FUNCTG(X,N,I,A,NTERMS) - evaluates the fitting function C for the Ith term C GCHISQ(NPTS,Y,NFREE,YFIT) - evaluates C reduced chi-squared for fit C GDERIV(N,X,A,NTERMS,I,DERIV) - evaluates the derivatives C of the fitting function for the Ith term with C respect to each parameter C MATINV(ARRAY,NTERMS,DET) - inverts a symmetric two dimensional C matrix of degree NTERMS and calculates its C determinant C IMPLICIT NONE INTEGER NPTS INTEGER NTERMS INTEGER IFAIL DOUBLE PRECISION X(NPTS) DOUBLE PRECISION Y(NPTS) DOUBLE PRECISION WEIGHT(NPTS) DOUBLE PRECISION A(NTERMS) DOUBLE PRECISION FLAMDA DOUBLE PRECISION YFIT(NPTS) DOUBLE PRECISION CHISQR C C Local variables C INTEGER I,J,K INTEGER NFREE INTEGER STAT DOUBLE PRECISION ALPHA(3,3) DOUBLE PRECISION BETA(3) DOUBLE PRECISION DERIV(3) DOUBLE PRECISION ARRAY(3,3) DOUBLE PRECISION B(3) DOUBLE PRECISION DET DOUBLE PRECISION CHISQ1 DOUBLE PRECISION FUNCTG DOUBLE PRECISION GCHISQ EXTERNAL FUNCTG,GCHISQ C C Get number of degrees and test C NFREE=NPTS-NTERMS IF (NFREE.LE.1) THEN CALL UMSPUT(' Too few points for fit',1,0,STAT) CHISQR=0.0D0 IFAIL=1 GO TO 999 ENDIF IFAIL=1 C C Evaluate ALPHA and BETA matrices C DO J=1,NTERMS,1 BETA(J)=0.0D0 DO K=1,J,1 ALPHA(J,K)=0.0D0 ENDDO ENDDO DO I=1,NPTS,1 CALL GDERIV(NPTS,X,A,NTERMS,I,DERIV) DO J=1,NTERMS,1 BETA(J)=BETA(J) + : (WEIGHT(I)*( Y(I)-FUNCTG(NPTS,X,I,A,NTERMS) )* : DERIV(J)) DO K=1,J,1 ALPHA(J,K)=ALPHA(J,K) + (WEIGHT(I)*DERIV(J)*DERIV(K)) ENDDO ENDDO ENDDO DO J=1,NTERMS,1 DO K=1,J,1 ALPHA(K,J)=ALPHA(J,K) ENDDO ENDDO C C Evaluate chi-squared at starting point C DO I=1,NPTS,1 YFIT(I)=FUNCTG(NPTS,X,I,A,NTERMS) ENDDO CHISQ1=GCHISQ(NPTS,Y,WEIGHT,NFREE,YFIT) C C Invert modified curvature matrix to find new parameters C 71 DO J=1,NTERMS,1 DO K=1,NTERMS,1 IF ((ALPHA(J,J)*ALPHA(K,K)).GT.0.0D0) THEN ARRAY(J,K)=ALPHA(J,K)/DSQRT(ALPHA(J,J)*ALPHA(K,K)) ELSE ARRAY(J,K)=0.0D0 ENDIF ENDDO ARRAY(J,J)=1.0D0 + FLAMDA ENDDO CALL MATINV(ARRAY,NTERMS,DET) DO J=1,NTERMS,1 B(J)=A(J) DO K=1,NTERMS,1 IF ((ALPHA(J,J)*ALPHA(K,K)).GT.0.0D0) THEN B(J)=B(J) + BETA(K)*ARRAY(J,K)/DSQRT(ALPHA(J,J)*ALPHA(K,K)) ELSE B(J)=B(J) ENDIF ENDDO ENDDO C C If chi-squared increased, increase Flamda and retry C DO I=1,NPTS,1 YFIT(I)=FUNCTG(NPTS,X,I,A,NTERMS) ENDDO CHISQR=GCHISQ(NPTS,Y,WEIGHT,NFREE,YFIT) IF (FLAMDA.EQ.0.0D0) THEN GO TO 101 ENDIF IF ((CHISQ1-CHISQR).LT.0.0D0) THEN FLAMDA=10.0D0*FLAMDA GO TO 71 ENDIF C C Evaluate parameters C 101 DO J=1,NTERMS,1 A(J)=B(J) ENDDO FLAMDA=FLAMDA/10.0D0 999 END DOUBLE PRECISION FUNCTION FUNCTG(N,X,IW,A,NTERMS) C C Function subprogram to evaluate a Gaussian of C height A(1), peak position A(2) and FWHM A(3). C Taken from Bevington, 1969, p.214 C IMPLICIT NONE INTEGER N INTEGER IW INTEGER NTERMS DOUBLE PRECISION X(N) DOUBLE PRECISION A(NTERMS) C C Local variables C DOUBLE PRECISION XX DOUBLE PRECISION Z,Z2 XX=X(IW) Z=(XX-A(2))/A(3) Z2=Z*Z IF (Z2.LT.200.D0) THEN FUNCTG=A(1)*DEXP(-Z2/2.0D0) ENDIF 99 END DOUBLE PRECISION FUNCTION GCHISQ(NPTS,Y,WEIGHT,NFREE,YFIT) C C Function subprogram to evaluate reduced chi-squared for fit to C data. C Taken from Bevington, 1969, p. 194. C IMPLICIT NONE INTEGER NPTS INTEGER NFREE DOUBLE PRECISION Y(NPTS) DOUBLE PRECISION WEIGHT(NPTS) DOUBLE PRECISION YFIT(NPTS) C INTEGER I DOUBLE PRECISION CHISQ CHISQ=0.0D0 IF (NFREE.LE.1) THEN GCHISQ=0.0D0 GO TO 99 ENDIF C C Accumulate chi-squared C DO I=1,NPTS,1 CHISQ=CHISQ + ( WEIGHT(I)*((Y(I)-YFIT(I))**2) ) ENDDO C C Divide by number of degrees of freedom C IF (NFREE.GT.1) THEN GCHISQ=CHISQ/DBLE(NFREE) ELSE GCHISQ=0.0D0 ENDIF 99 END SUBROUTINE GDERIV(N,X,A,NTERMS,IW,DERIV) C C Subroutine subprogram to evaluate the derivatives of a C Gaussian of height A(1), peak position A(2) and FWHM C A(3). C Taken from Bevington, 1969, p.241 C IMPLICIT NONE INTEGER N INTEGER IW INTEGER NTERMS DOUBLE PRECISION X(N) DOUBLE PRECISION A(NTERMS) DOUBLE PRECISION DERIV(NTERMS) C C Local variables C DOUBLE PRECISION XX DOUBLE PRECISION Z DOUBLE PRECISION Z2 XX=X(IW) Z=(XX-A(2))/A(3) Z2=Z*Z IF (Z2.GT.200.0D0) THEN DERIV(1)=0.0D0 DERIV(2)=0.0D0 DERIV(3)=0.0D0 ELSE DERIV(1)=DEXP(-Z2/2.0D0) DERIV(2)=A(1)*DERIV(1)*Z/A(3) DERIV(3)=DERIV(2)*Z ENDIF 99 END SUBROUTINE PSFBLKY(P1,P2,PSF,BX,BY,CENI,ITYPE,WORKR,N2,WORKD, : PX,PY,PSFC) C C This subroutine shifts the input (subsampled by BX) PSF C array PSF to channel P2/2 and blocked-sums the subsampled PSF C to the output array PSFC for all X-channels of the input C array C IMPLICIT NONE INTEGER P1 ! 1st dimension of the input array PIN INTEGER P2 ! 2nd dimension of the input array PIN INTEGER BX ! Subsampling factor in 1st dimension of PSF INTEGER BY ! Subsampling factor in 2nd dimension of PSF INTEGER ITYPE ! Interpolation code: 1=Nearest; 2=Linear; 3=Poly3; 4=Poly5; 5=Spline3 INTEGER N2 ! Dimension of work array WORKD INTEGER PX ! 1st dimension of the output array PSF INTEGER PY ! 2nd dimension of the output array PSF REAL WORKR(P2) ! Work array for storing the shifted PSF DOUBLE PRECISION PSF(P1,P2) ! Input array of subsampled PSF DOUBLE PRECISION CENI(P1) ! Array of position of centre of PSF profile for all X-sections DOUBLE PRECISION WORKD(N2) ! Work array for storing the shifted PSF DOUBLE PRECISION PSFC(PX,PY) ! Output array of block-summed PSF C C Local variables C INTEGER I,J,K INTEGER JJ DOUBLE PRECISION CENO c integer istat c character*80 text C C Initialize the output array C DO J=1,PY,1 DO I=1,PX,1 PSFC(I,J)=0.0D0 ENDDO ENDDO C C Shift the array PSF in the Y direction to be centred at C position P1/2 for all X-sections C DO I=1,P1,1 C C Shift the I'th X-section of the array PSF in the Y C direction so that the block summed PSF is at channel C P2/2 C CENO=DBLE(P2/2) CALL SHIFTY(P1,P2,PSF,CENI,CENO,ITYPE,I,WORKR,N2,WORKD) c if (i.eq.128) then c do j=1,n2,1 c if (workd(j).gt.0.0d0) then c write(text,33) j,workd(J) c33 format(' J,WORKD ',I5,1x,F12.5) c call umsput(text,1,0,istat) c endif c enddo c endif C C Block sum the data into the larger pixels C JJ=PY/2 - (P2/2)/BY + 1 DO J=1,PY/BY,1 DO K=1,BY,1 IF (JJ.GT.1.AND.JJ.LT.PY) THEN PSFC(I,JJ)=PSFC(I,JJ) + WORKD((J-1)*BY+K) ENDIF ENDDO c write(text,43) jj,j,psfc(i,jj) c43 format(' JJ,J,PSFC ',2(1x,I5),1x,F12.5) c call umsput(text,1,0,istat) JJ=JJ+1 ENDDO ENDDO 99 END SUBROUTINE SHIFTY(P1,P2,PSF,PCEN,CENR,ITYPE,IS,YSR,N2,PSFSH) C C This subroutine shifts the IS'th X-section of the array C PSF in the Y direction so that the peak is at position C CENR. The interpolation method is specified by ITYPE. The C output shifted PSF is PSFSH C IMPLICIT NONE INTEGER P1 ! 1st dimension of input image PSF INTEGER P2 ! 2nd dimension of input image PSF INTEGER ITYPE ! Interpolation type for rebinning INTEGER IS ! X-section to be shifted INTEGER N2 ! Dimension of output array PSFSH REAL YSR(P1) ! Work array required for interpolation function DOUBLE PRECISION PSF(P1,P2) ! Input PSF image DOUBLE PRECISION PCEN(P1) ! Array of positions of peak of profile in PSF for all X-sections DOUBLE PRECISION CENR ! Desired centre of shifted PSF array DOUBLE PRECISION PSFSH(N2) ! Output shifted PSF array centred at CEN C C Local variables C INTEGER J,JJ INTEGER NY REAL X1 REAL Y1 REAL MRIEVL NY=1 Y1=1.0 C C Initialize output array C DO J=1,N2,1 PSFSH(J)=0.0D0 ENDDO C C Convert the input YSF to real in preparation for shifting C DO J=1,P2,1 YSR(J)=REAL(PSF(IS,J)) ENDDO C C Shift the PSF by CENR-CENP in the output array C DO J=1,P2,1 X1=REAL(DBLE(J)-(CENR-PCEN(IS))+INT(CENR-PCEN(IS))) JJ=J+INT(CENR-PCEN(IS)) IF (X1.GE.1.0.AND.X1.LE.REAL(P2). : AND.JJ.GE.1.AND.JJ.LE.N2) THEN PSFSH(JJ)=DBLE(MRIEVL(X1,Y1,YSR,P2,NY,P2,ITYPE)) ENDIF ENDDO 99 END SUBROUTINE POICENN(NW,CENCX,CENCY,CENM,CENT,N1,NS,NS1,CEN) C C This subroutine returns the real centres of the NS profiles C within the actual range of the spectrum image for all channels C N1, calculated from the reference X,Y centre (CENCX,CENCY) C and the slope with respect to the X axis CENM C IMPLICIT NONE INTEGER NW ! Number of elements in input arrays INTEGER N1 ! Number of X channels for profile centre calculation INTEGER NS ! Number of profiles in Y channel range 2 to N2-2 INTEGER NS1 ! Index number of first profile in input arrays REAL CENCX(NW) ! Input array of reference X value of profiles REAL CENCY(NW) ! Input array of reference Y value of profiles REAL CENM(NW) ! Input array of Y v. X slope of profile centres REAL CENT(NW) ! Centres of profiles for channel I DOUBLE PRECISION CEN(N1,NS) ! Output centres of profiles for all N1 channels C C Local variables C INTEGER I,J C C Compute the Y centre of the profiles for channel I. C DO I=1,N1,1 CALL POICENI(NW,CENCX,CENCY,CENM,I,CENT) DO J=1,NS,1 CEN(I,J)=DBLE(CENT(NS1+J-1)) ENDDO ENDDO 99 END SUBROUTINE POICENI(NW,CENCX,CENCY,CENM,IS,CEN) C C This subroutine returns the real centres of the NW profiles C appropriate for channel IS, calculated from the reference C X,Y centre (CENCX,CENCY) and the slope with respect to the C X axis CENM C IMPLICIT NONE INTEGER NW ! Number of elements in input arrays INTEGER IS ! I channel number to calculate centres REAL CENCX(NW) ! Input array of reference X value of profiles REAL CENCY(NW) ! Input array of reference Y value of profiles REAL CENM(NW) ! Input array of Y v. X slope of profile centres REAL CEN(NW) ! Output centres of profiles for channel I C C Local variables C INTEGER I C C Compute the Y centre of the profiles for channel I C DO I=1,NW,1 CEN(I)=(REAL(IS)-CENCX(I))*CENM(I) + CENCY(I) ENDDO 99 END SUBROUTINE PSFKER(N,Y,CEN,SIG,PSFR,IRL1,IRL2) C C This subroutine forms the 1-D PSF of the resolution kernel C PSFR of sigma SIG centred at CEN in array Y. C The channel limits corresponding to +/- 10 sigma are also C returned as IRLI1,IRLI2 C IMPLICIT NONE INTEGER N ! Dimension of output array PSFR INTEGER IRL1 ! Lower limit of line centre - 10*SIG INTEGER IRL2 ! Upper limit of line centre + 10*SIG DOUBLE PRECISION Y(N) ! Y values of the PSFR pixels DOUBLE PRECISION SIG ! Sigma of Gaussian to produce DOUBLE PRECISION CEN ! Centre of Gaussian peak DOUBLE PRECISION PSFR(N) ! Output 1-D Gaussian PSF C C Local variables C INTEGER I DOUBLE PRECISION XD DOUBLE PRECISION GAUSSF DO I=1,N,1 XD=Y(I)-CEN PSFR(I)=GAUSSF(XD,SIG) ENDDO C C Set limits for extent of PSFR to + & - 10 sigma about peak C IRL1=INT(CEN - 10.0D0*DABS(SIG)) IRL2=NINT(CEN + 10.0D0*DABS(SIG)) IF (IRL1.LT.1) THEN IRL1=1 ENDIF IF (IRL2.GT.N) THEN IRL2=N ENDIF 99 END SUBROUTINE POSCORN(N1,N2,Z,P1,P2,PSF,Y,NW,CENT,ICSTEP,PCEN,PFWHM, : ZS,PSFS,WT,YCOR,YFIT,N22,YC,PSFC,CCROS,TEMP) C C This subroutine cross correlates Z against PSF X-section C by section, fits, in the Y direction, the peak of the cross C correlation nearest CENT by a Gaussian. C If the peak of Z is at a higher pixel value than C for PSF, then the shift is POSITIVE C The accurate value for the position of the profile for C each X-section of Z is returned using the shift of the peak C of the profile from the PSF array. C If the FWHM of the cross correlation was estimated at zero C or the peak height of the cross correlation is less than or C equal to zero, then LFAIL is set true and the input value C of the position is used. C If the cross correlation peak was outside the range C PKPOS-FWHM to PKPOS+FWHM (where FWHM is the estimated C full width at half max) then LWAR is set true and C the input value of the position is used C IMPLICIT NONE INTEGER N1 ! 1st dimension of input array Z INTEGER N2 ! 2nd dimension of input array Z INTEGER P1 ! 1st dimension of input array PSF INTEGER P2 ! 2nd dimension of input array PSF INTEGER NW ! 2nd dimension of input array CENT INTEGER ICSTEP ! No. of X-sections of array Z to sum for cross correlation INTEGER N22 ! Dimension of arrays for FFT DOUBLE PRECISION Z(N1,N2) ! Input array whose profile positions to be determined DOUBLE PRECISION PSF(P1,P2) ! Input PSF array for cross-corr against DOUBLE PRECISION Y(N2) ! Input array of values of Y pixel centres DOUBLE PRECISION CENT(N1,NW) ! Input array of rough positions of (N1) profiles = Output array of exact profile positions DOUBLE PRECISION PCEN(N1) ! Input array of (N1) positions of PSF DOUBLE PRECISION PFWHM(N1) ! Input array of (N1) FWHM's of PSF DOUBLE PRECISION ZS(N2) ! Work array for 1-D X-sections of array Z DOUBLE PRECISION PSFS(N2) ! Work array for 1-D X-sections of array PSF DOUBLE PRECISION WT(N2) ! Work array for 1-D weight for Gaussian fitting DOUBLE PRECISION YCOR(N2) ! Work array for 1-D signal for Gaussian fitting DOUBLE PRECISION YFIT(N2) ! Work array for 1-D fit for Gaussian fitting DOUBLE PRECISION TEMP(N22) ! Work array for storage in FFT DOUBLE COMPLEX YC(N2) ! Work array for 1-D complex copy of Z DOUBLE COMPLEX PSFC(N2) ! Work array for 1-D complex copy of PSF DOUBLE COMPLEX CCROS(N2) ! Work array for 1-D complex cross-correlation C C Local variables C INTEGER I,J,K INTEGER II INTEGER NITPC INTEGER IL INTEGER IR INTEGER ISTAT DOUBLE PRECISION YSUM DOUBLE PRECISION CCSHIFT CHARACTER*80 TEXT LOGICAL LFAIL LOGICAL LWAR LFAIL=.FALSE. LWAR=.FALSE. NITPC=20 ! Number of iterations of LS Gaussian fitting of Cross-Corr peak DO I=1,N1,1 C C Copy the I'th X-section of the PSF to a 1-D array. Take account C that the PSF array Y dimension P2 may not be equal to N2 C DO J=1,N2,1 IF (J.GE.1.AND.J.LE.P2) THEN PSFS(J)=PSF(I,J) ELSE PSFS(J)=0.0D0 ENDIF c if (psfs(J).gt.0.0d0) then c write(text,33) j,psfs(J) c33 format(' J, PSFS(J) ',I5,1x,F12.5) c call umsput(text,1,0,istat) c endif ENDDO C C Set up the lower (IL) and upper (IR) limits of the C X region to sum for determining the position of the C point source(s) by cross correlation C IL=I - ICSTEP/2 IF (IL.LT.1) THEN IL=1 ENDIF IR=IL+ICSTEP-1 IF (IR.GT.N1) THEN IR=N1 ENDIF IF ((IR-IL).LT.ICSTEP-1) THEN IL=N1-(ICSTEP-1) IF (IL.LT.1) THEN IL=1 ENDIF ENDIF C C Copy the range of X-sections from IL to IR of the array Z C to a 1-D array C DO J=1,N2,1 YSUM=0.0D0 DO II=IL,IR,1 YSUM=YSUM+Z(II,J) ENDDO ZS(J)=YSUM c if (zs(J).gt.0.0d0) then c write(text,43) j,zs(J) c43 format(' J, ZS(J) ',I5,1x,F12.5) c call umsput(text,1,0,istat) c endif ENDDO C C Determine the position of the point sources by cross correlation C with the PSF C DO K=1,NW,1 CALL POSCORR(N2,Y,ZS,PSFS,CENT(I,K),PCEN(I),PFWHM(I), : NITPC,WT,YCOR,YFIT,YC,PSFC,CCROS, : 2*N2,TEMP,CCSHIFT,LFAIL,LWAR) IF (LFAIL.OR.LWAR) THEN WRITE(TEXT,41) I 41 FORMAT(' ** Warning. Cross-correlation failed at ', : 'X-section ',I5) CALL UMSPUT(TEXT,1,0,ISTAT) ELSE CENT(I,K)=CCSHIFT + PCEN(I) ENDIF ENDDO ENDDO 99 END SUBROUTINE POSCORR(N2,X,Y,PSF,PKPOS,PCENT,PFWHM,NITER, : WT,YCOR,YFIT,YC,PSFC,CCROS,N22,TEMP, : SHIFT,LFAIL,LWAR) C C This subroutine cross correlates spectrum Y against PSF, C fits, in the X direction, the peak of the cross C correlation nearest PKPOS by a Gaussian. C If the peak of Y is at a higher pixel value than C for PSF, then the shift is POSITIVE C The shift SHIFT of the peak of the profile is returned. C If the FWHM was estimated at zero or the peak of the cross C correlation is less than or equal to zero, then LFAIL is C set true. C If the cross correlation peak was outside the range C PKPOS-FWHM to PKPOS+FWHM (where FWHM is the estimated C full width at half max) then LWAR is set true C IMPLICIT NONE INTEGER N2 ! X dimension of the input arrays INTEGER N22 ! X input of the storage arrays for FFT INTEGER NITER ! Number of iterations for fitting of cross-corr peak by a Gaussian DOUBLE PRECISION PKPOS ! Approximate X position of peak of Y array DOUBLE PRECISION PCENT ! Exact X position of peak of PSF array DOUBLE PRECISION PFWHM ! FWHM of peak of Y array DOUBLE PRECISION X(N2) ! Input array of X pixel values DOUBLE PRECISION Y(N2) ! Input array of spectrum signal values containing peak DOUBLE PRECISION PSF(N2) ! Input array of PSF to be cross-correlated with Y DOUBLE PRECISION WT(N2) ! Holding array for weights for Gaussian fitting DOUBLE PRECISION YCOR(N2) ! Holding array for unwrapped cross-correlation DOUBLE PRECISION YFIT(N2) ! Holding array for Gaussian fit to cross-correl peak DOUBLE PRECISION TEMP(N22) ! Holding array for FFT of spectrum Y DOUBLE PRECISION SHIFT ! Output shift of Y with respect to PSF DOUBLE COMPLEX YC(N2) ! Holding array for Y values for FFT DOUBLE COMPLEX PSFC(N2) ! Holding array for PSF values for FFT DOUBLE COMPLEX CCROS(N2) ! Holding array for cross correlation of Y against PSF LOGICAL LFAIL ! Output failure flag - set TRUE to indicate failure in fitting cross-correl peak LOGICAL LWAR ! Output warning flag - set true to indicate a problem in fitting cross-correl peak C C Local variables C INTEGER J,K c INTEGER STAT DOUBLE PRECISION DPKPOS DOUBLE PRECISION CENT DOUBLE PRECISION RCENT c CHARACTER*72 OUTEXT LOGICAL LWRAP1 LOGICAL LWRAP2 LOGICAL LFAILS LOGICAL LWARS LFAIL=.FALSE. LWAR=.FALSE. c write(outext,1291) PCENT,PFWHM c1291 format(' PCENT,PFWHM',2(F12.5)) c call umsput(outext,1,0,stat) C C Initialize the weight array to unity and zero the cross C correlation array C DO J=1,N2,1 WT(J)=1.0D0 YCOR(J)=0.0D0 ENDDO C C Prepare the data and PSF arrays for Fourrier C transform C DO J=1,N2,1 YC(J)=DCMPLX(Y(J)) PSFC(J)=DCMPLX(PSF(J)) ENDDO C C Cross correlate YC v. PSFC C CALL COREL(N2,YC,PSFC,N22,TEMP,CCROS) C C Unwrap the cross correlation and set centre at C N2/2 +1 C LWRAP1=.FALSE. LWRAP2=.FALSE. DO J=1,N2,1 K=J+N2/2 IF (K.GT.N2) THEN K=K-N2 ENDIF YCOR(K)=DBLE(CCROS(J)) ENDDO RCENT=DBLE(N2/2+1) ! Zero lag at N/2+1 c do j=1,n2,1 c if (ycor(j).gt.0.0d0) then c write(outext,119) J,YCOR(J) c119 format(' J,YCOR ',i5,1x,1(F12.5)) c call umsput(outext,1,0,stat) c endif c enddo C C Fit the line and background of both PSF and data. Take C care if the expected position of the cross-corr peak C is <1, caused by wrapping of CCROS C DPKPOS=PKPOS - PCENT + RCENT ! Expected posn of cross-corr peak IF (DPKPOS.LT.1.0D0) THEN LWRAP1=.TRUE. DPKPOS=DPKPOS+N2 ENDIF IF (DPKPOS.GT.DBLE(N2)) THEN LWRAP2=.TRUE. DPKPOS=DPKPOS-DBLE(N2+1) ENDIF c write(outext,129) PCENT,RCENT,DPKPOS,PFWHM c129 format(' PCENT,RCENT,DPKPOS',4(F12.5)) c call umsput(outext,1,0,stat) C C Fit peak of cross-correlation by a Gaussian using Levenberg Marquardt C CALL GAUSFIN(N2,X,YCOR,DPKPOS,PFWHM,NITER,WT,YFIT,CENT, : LFAILS,LWARS) IF (LFAILS) THEN LFAIL=.TRUE. ENDIF IF (LWARS) THEN LWAR=.TRUE. ENDIF C C Convert the fitted cross correlation peak to a shift C taking account of any wrap of the cross-correlation C about position of zero lag (LWRAP1 set TRUE) or C shifts larger than N2 (LWRAP2) C SHIFT=(CENT-RCENT) IF (LWRAP1) THEN SHIFT=(CENT-RCENT)-DBLE(N2) ENDIF IF (LWRAP2) THEN SHIFT=(CENT+RCENT) - 2.0D0 ENDIF c write(outext,149) SHIFT c149 format(' SHIFT ',1(F12.5)) c call umsput(outext,1,0,stat) c IF (LWARS) THEN c WRITE(OUTEXT,298) I,CENT c298 FORMAT(' X-chan. ',I5, c : '. WARNING - Problem fitting cross-corr. for peak',F8.3) c CALL UMSPUT(OUTEXT,1,0,STAT) c ENDIF 999 END SUBROUTINE GAUSFIN(N,X,Y,PKPOS,PWID,NITER,WHT,YFIT,CENT, : LFAIL,LWAR) C C This subroutine determines the position of the peak C of the profile Y as CENT by fitting a Gaussian to the C line and a straight line to the continuum by non-linear C least squares. C If the FWHM was estimated at zero or the peak of the cross C correlation is less than or equal to zero, then LFAIL is C set true. C If an the least squares fitted cross-corr peak is outside C the range CENT-FWHM to CENT+FWHM then LWAR is set true C IMPLICIT NONE INTEGER N ! X dimension of input arrays INTEGER NITER ! No. of iterations for Levenberg-Marquardt fitting DOUBLE PRECISION X(N) ! Input array of X pixel positions DOUBLE PRECISION Y(N) ! Input array of signal values to fit DOUBLE PRECISION PKPOS ! Estimated pixel position of centre of cross-correl. DOUBLE PRECISION PWID ! Estmated FWHM of peak of cross-correl. DOUBLE PRECISION WHT(N) ! Holding array for weights for fitting DOUBLE PRECISION YFIT(N) ! Holding array for fit to Y array DOUBLE PRECISION CENT ! Output Gaussian fitted X centre of cross-correl. peak LOGICAL LFAIL ! Output flag - set TRUE if approximate FWHM of peak is zero or peak height < 0 LOGICAL LWAR ! Output flag - set TRUE if an error occured in the fitting routine CURFITS C C Local variables C INTEGER K INTEGER IH INTEGER IFWHM INTEGER ILOW INTEGER IHIG INTEGER NTERMS INTEGER IFAIL DOUBLE PRECISION BG DOUBLE PRECISION WFAC DOUBLE PRECISION FWHM DOUBLE PRECISION M DOUBLE PRECISION C DOUBLE PRECISION FLAMDA DOUBLE PRECISION A(5) DOUBLE PRECISION CHISQ DOUBLE PRECISION OCHISQ DOUBLE PRECISION FRACHI c integer stat c character*128 text LFAIL=.FALSE. LWAR=.FALSE. IH=INT(PKPOS) C C Determine a suitable mean background value as the mean C value of the median in the two regions from the peak to C the ends. Calculate the slope and intercept of a straight C line fit to these two median values C CALL BACKCC(N,X,Y,IH,YFIT,BG,M,C) c write(text,121) bg c121 format(' Mean median BG',E12.5) c call umsput(text,1,0,stat) C C Determine centroid of cross-corr profile and FWHM C WFAC=2.0D0 ! Centroid over PKPOS - WFAC*PWID to PKPOS+WFAC*PWID CALL CENTPOS(N,X,Y,PKPOS,PWID,WFAC,BG,CENT,FWHM) c write(text,123) CENT,FWHM c123 format(' CENTPOS: CEN,FWHM ',2(1x,F12.5)) c call umsput(text,1,0,stat) IH=INT(CENT) C C Error condition if FWHM < zero C IF (FWHM.LE.0.0D0) THEN LFAIL=.TRUE. GO TO 99 ENDIF C C Set lower and upper limits to region to consider for C goodness of fit in the linear least sqaures fitting C Use region peak-3*fwhm to peak+3*fwhm C IFWHM=INT(FWHM) ILOW=IH-3*IFWHM IF (ILOW.LT.1) THEN ILOW=1 ENDIF IHIG=IH+3*IFWHM IF (IHIG.GT.N) THEN IHIG=N ENDIF C C Set up parameters for CURFITS to fit the line by a C Gaussian and the continuum by a straight line C using least squares with linearization of the fitting C function C c write(text,201) CENT,(Y(IH)-(M*X(IH)+C)),FWHM,M,C c201 format(' Xpk,Ypk,FWHM,M,C',5(1x,F11.5)) c call umsput(text,1,0,stat) NTERMS=5 FLAMDA=0.001D0 A(1)=Y(IH)-BG A(2)=CENT A(3)=FWHM/2.354820D0 A(4)=M A(5)=C C C Error condition if cross-corr height above continuum < zero C IF (A(1).LE.0.0D0) THEN LFAIL=.TRUE. GO TO 99 ENDIF C C Call the least squares fitting routine until the fractional C change in Chi-squared is less than some value (set to 1.0E-6) C FRACHI=1.0D-6 OCHISQ=1.0D-6 ! Small starting value DO K=1,NITER,1 CALL CURFITS(N,X,Y,WHT,NTERMS,ILOW,IHIG,A,FLAMDA,YFIT, : CHISQ,IFAIL) IF (IFAIL.EQ.0) THEN IF (A(2).LT.(PKPOS-FWHM).OR.A(2).GT.(PKPOS+FWHM)) THEN CENT=PKPOS LWAR=.TRUE. GO TO 99 ENDIF IF (CHISQ.LT.1.0D-6) THEN GO TO 80 ENDIF IF (DABS((CHISQ/OCHISQ)-1.0D0).LT.FRACHI) THEN GO TO 80 ELSE OCHISQ=CHISQ ENDIF ELSE LFAIL=.TRUE. GO TO 99 ENDIF ENDDO 80 CENT=A(2) 99 END SUBROUTINE BACKCC(N,X,Y,IPK,YY,BG,M,C) C C This subroutine returns the mean value of the C background BG from the medians in the two sections C from the peak to the ends. The slope M and intercept C C of the straight line defined by the two median values C are also returned C IMPLICIT NONE INTEGER N INTEGER IPK DOUBLE PRECISION X(N) DOUBLE PRECISION Y(N) DOUBLE PRECISION YY(N) DOUBLE PRECISION BG DOUBLE PRECISION M DOUBLE PRECISION C C C Local variables C INTEGER I,J INTEGER IP INTEGER ILEF INTEGER IRIG DOUBLE PRECISION LMED DOUBLE PRECISION UMED LOGICAL LLEF LOGICAL LRIG LLEF=.TRUE. LRIG=.TRUE. C C Copy the values of Y from 1 to peak channel-2 C to a holding array. Take care of peak near to end C of array C IP=IPK-2 ! Start to left of peak IF (IP.GT.1) THEN J=1 DO I=1,IP,1 YY(J)=Y(I) J=J+1 ENDDO C C Compute median C CALL MDIAN1(IP,YY,LMED) ILEF=(IP-1)/2 IF (ILEF.LT.1) THEN ILEF=1 ENDIF ELSE LLEF=.FALSE. ILEF=1 ENDIF C C Copy the values of Y from peak channel+2 to N C to a holding array. Take care of peak near to end C of array C IP=IPK+2 ! Start to right of peak IF (IP.LT.N) THEN J=1 DO I=IP,N,1 YY(J)=Y(I) J=J+1 ENDDO C C Computer median C CALL MDIAN1(N-IP+1,YY,UMED) IRIG=(N-IP)/2 + IP IF (IRIG.GT.N) THEN IRIG=N ENDIF ELSE IRIG=N LRIG=.FALSE. ENDIF C C Compute the mean of the two medians C IF (.NOT.LLEF) THEN BG=UMED ENDIF IF (.NOT.LRIG) THEN BG=LMED ENDIF IF (LLEF.AND.LRIG) THEN BG=0.50D0*(LMED+UMED) ENDIF C C Compute the slope and intercept of the two background C estimates from a straight line fit C M=(UMED-LMED)/(X(IRIG)-X(ILEF)) C=UMED-M*X(IRIG) 99 END SUBROUTINE MDIAN1(N,X,XMED) C C This is the Numerical Recipes 1st edition subroutine C to compute the rough median of an array X of N numbers. C The array is modified and returned sorted into ascending C order C IMPLICIT NONE INTEGER N DOUBLE PRECISION X(N) DOUBLE PRECISION XMED C C Local variables C INTEGER N2 CALL SORT(N,X) N2=N/2 IF (2*N2.EQ.N) THEN XMED=0.50D0*(X(N2)+X(N2+1)) ELSE XMED=X(N2+1) ENDIF 99 END SUBROUTINE SORT(N,RA) C C This is the Numerical Recipes subroutine SORT to perform the C HEAPSORT algorithm on array RA. C It sorts an array RA of length N into ascending numerical order. C RA is replaced on output by its sorted rearrangement. C IMPLICIT NONE INTEGER N DOUBLE PRECISION RA(N) C C Local variables C INTEGER I,J,L,IR DOUBLE PRECISION RRA L=N/2+1 IR=N 10 CONTINUE IF (L.GT.1) THEN L=L-1 RRA=RA(L) ELSE RRA=RA(IR) RA(IR)=RA(1) IR=IR-1 IF (IR.EQ.1) THEN RA(1)=RRA GO TO 99 ENDIF ENDIF I=L J=L+L 20 IF (J.LE.IR) THEN IF (J.LT.IR) THEN IF (RA(J).LT.RA(J+1)) THEN J=J+1 ENDIF ENDIF IF (RRA.LT.RA(J)) THEN RA(I)=RA(J) I=J J=J+J ELSE J=IR+1 ENDIF GO TO 20 ENDIF RA(I)=RRA GO TO 10 99 END SUBROUTINE CENTPOS(N,X,Y,POS,WID,WFAC,BACK,CENT,FWHM) C C This subroutine determines the simple centroid of the C profile in the array Y at position POS as CENT C and its FWHM. The centroiding is performed over the C range +/- WFAC*WID about the peak position above the C mean background BACK C IMPLICIT NONE INTEGER N ! X dimension of input arrays DOUBLE PRECISION X(N) ! Input array of X pixel values DOUBLE PRECISION Y(N) ! Input array of signal values DOUBLE PRECISION POS ! Input expected position of peak of Y array DOUBLE PRECISION WID ! Input expected FWHM of peak of Y array DOUBLE PRECISION WFAC ! Factor * FWHM over which to search for peak DOUBLE PRECISION BACK ! Median estimate of flat background to Y array DOUBLE PRECISION CENT ! Output centroid position of peak of Y array DOUBLE PRECISION FWHM ! Output FWHM of peak of Y array C C Local variables C INTEGER I INTEGER ILEFT INTEGER IRIGHT INTEGER IH DOUBLE PRECISION LEFT DOUBLE PRECISION RIGHT DOUBLE PRECISION PK DOUBLE PRECISION SUMX DOUBLE PRECISION SUMY c integer stat c character*128 outext C C Determine the range over which to compute the centroid. C Take the limits as +/- WFAC*FWHM C LEFT=POS-WFAC*WID RIGHT=POS+WFAC*WID ILEFT=1 IRIGHT=N DO I=2,N-1,1 IF (LEFT.GT.X(I-1).AND.LEFT.LE.X(I)) THEN ILEFT=I ENDIF IF (RIGHT.GT.X(I).AND.RIGHT.LE.X(I+1)) THEN IRIGHT=I ENDIF ENDDO C C Check WFAC*FWHM limits within array bounds C IF (ILEFT.LT.1) THEN ILEFT=1 ENDIF IF (IRIGHT.GT.N) THEN IRIGHT=N ENDIF C C Determine simple centroid C SUMX=0.0D0 SUMY=0.0D0 DO I=ILEFT,IRIGHT,1 SUMX=SUMX+X(I)*(Y(I)-BACK) SUMY=SUMY+(Y(I)-BACK) ENDDO IF (SUMY.NE.0.0D0) THEN CENT=SUMX/SUMY ELSE CENT=POS ENDIF c write(outext,180) CENT c180 format(' CENT ',F12.5) c call umsput(outext,1,0,stat) C C Check that value of CENT is within valid limits C IF (CENT.LT.X(1).OR.CENT.GT.X(N)) THEN CENT=0.0D0 FWHM=0.0D0 GO TO 99 ENDIF IH=NINT(CENT) C C Determine an estimate of the FWHM from the half C height positions C PK=Y(IH)-BACK ILEFT=IH-1 IF (ILEFT.LT.1) THEN ILEFT=1 ENDIF IRIGHT=IH+1 IF (IRIGHT.GT.N) THEN IRIGHT=N ENDIF 50 DO I=IH,1,-2 IF ((Y(I)-BACK).LT.(0.5D0*PK).AND. : (Y(I-1)-BACK).LT.(0.5D0*PK)) THEN ILEFT=I GO TO 60 ENDIF ENDDO 60 DO I=IH,N,2 IF ((Y(I)-BACK).LT.(0.5D0*PK).AND. : (Y(I+1)-BACK).LT.(0.5D0*PK)) THEN IRIGHT=I GO TO 70 ENDIF ENDDO 70 FWHM=X(IRIGHT)-X(ILEFT) - 1.0D0 99 END SUBROUTINE CURFITS(NPTS,X,Y,WEIGHT,NTERMS,ILEF,IRIG,A, : FLAMDA,YFIT,CHISQR,IFAIL) C C Performs a least squares fit to a non-linear function with C a linearization of the fitting function. C The function is a single Gaussian on a stright line. C The goodness of fit is only considered over the channel C range ILEF to IRIG C The program is adapted from Bevington, 1968, p.237-9. C C Subprograms called: C FUNCTN(X,N,I,A,NTERMS) - evaluates the fitting function C for the Ith term C FCHISQ(NPTS,Y,NFREE,ILEF,IRIG,YFIT) - evaluates C reduced chi-squared for fit over range ILEF to IRIG C FDERIV(N,X,A,NTERMS,I,DERIV) - evaluates the derivatives C of the fitting function for the Ith term with C respect to each parameter C MATINV(ARRAY,NTERMS,DET) - inverts a symmetric two dimensional C matrix of degree NTERMS and calculates its C determinant C IMPLICIT NONE INTEGER NPTS INTEGER NTERMS INTEGER ILEF INTEGER IRIG INTEGER IFAIL DOUBLE PRECISION X(NPTS) DOUBLE PRECISION Y(NPTS) DOUBLE PRECISION WEIGHT(NPTS) DOUBLE PRECISION A(NTERMS) DOUBLE PRECISION FLAMDA DOUBLE PRECISION YFIT(NPTS) DOUBLE PRECISION CHISQR C C Local variables C INTEGER I,J,K INTEGER NFREE INTEGER STAT DOUBLE PRECISION ALPHA(5,5) DOUBLE PRECISION BETA(5) DOUBLE PRECISION DERIV(5) DOUBLE PRECISION ARRAY(5,5) DOUBLE PRECISION B(5) DOUBLE PRECISION DET DOUBLE PRECISION CHISQ1 DOUBLE PRECISION FUNCTN DOUBLE PRECISION FCHISQ EXTERNAL FUNCTN,FCHISQ c character*128 outext C C Get number of degrees of freedom for number of points C in range ILEF to IRIG and test C NFREE=(IRIG-ILEF+1)-NTERMS IF (NFREE.LE.1) THEN CALL UMSPUT(' Too few points for fit',1,0,STAT) CHISQR=0.0D0 IFAIL=1 GO TO 999 ENDIF IFAIL=0 C C Zero ALPHA and BETA matrices C DO J=1,NTERMS,1 BETA(J)=0.0D0 DO K=1,J,1 ALPHA(J,K)=0.0D0 ENDDO ENDDO C C Evaluate ALPHA and BETA matrices C DO I=ILEF,IRIG,1 CALL FDERIV(NPTS,X,A,NTERMS,I,DERIV) DO J=1,NTERMS,1 BETA(J)=BETA(J) + : (WEIGHT(I)*( Y(I)-FUNCTN(NPTS,X,I,A,NTERMS) )* : DERIV(J)) DO K=1,J,1 ALPHA(J,K)=ALPHA(J,K) + (WEIGHT(I)*DERIV(J)*DERIV(K)) ENDDO ENDDO ENDDO DO J=1,NTERMS,1 DO K=1,J,1 ALPHA(K,J)=ALPHA(J,K) ENDDO ENDDO C C Evaluate chi-squared at starting point C DO I=ILEF,IRIG,1 YFIT(I)=FUNCTN(NPTS,X,I,A,NTERMS) ENDDO CHISQ1=FCHISQ(NPTS,Y,WEIGHT,NFREE,ILEF,IRIG,YFIT) C C Invert modified curvature matrix to find new parameters C 71 DO J=1,NTERMS,1 DO K=1,NTERMS,1 IF ((ALPHA(J,J)*ALPHA(K,K)).GT.0.0D0) THEN c write(outext,111) j,k,alpha(j,j),alpha(k,k) c111 format('j,k,alpha(j,j),alpha(k,k)',i4,1x,i4,2(1x,E13.5)) c call umsput(outext,1,0,stat) ARRAY(J,K)=ALPHA(J,K)/DSQRT(ALPHA(J,J)*ALPHA(K,K)) ELSE ARRAY(J,K)=0.0D0 ENDIF ENDDO ARRAY(J,J)=1.0D0 + FLAMDA ENDDO CALL MATINV(ARRAY,NTERMS,DET) DO J=1,NTERMS,1 B(J)=A(J) DO K=1,NTERMS,1 IF ((ALPHA(J,J)*ALPHA(K,K)).GT.0.0D0) THEN B(J)=B(J) + BETA(K)*ARRAY(J,K)/DSQRT(ALPHA(J,J)*ALPHA(K,K)) ELSE B(J)=B(J) ENDIF ENDDO ENDDO C C If chi-squared increased, increase Flamda and retry C DO I=ILEF,IRIG,1 YFIT(I)=FUNCTN(NPTS,X,I,A,NTERMS) ENDDO CHISQR=FCHISQ(NPTS,Y,WEIGHT,NFREE,ILEF,IRIG,YFIT) IF (FLAMDA.EQ.0.0D0) THEN GO TO 101 ENDIF IF ((CHISQ1-CHISQR).LT.0.0D0) THEN FLAMDA=10.0D0*FLAMDA GO TO 71 ENDIF C C Evaluate parameters and errors C 101 DO J=1,NTERMS,1 A(J)=B(J) ENDDO FLAMDA=FLAMDA/10.0D0 999 END DOUBLE PRECISION FUNCTION FUNCTN(N,X,IW,A,NTERMS) C C Function subprogram to evaluate a Gaussian of C height A(1), peak position A(2) and FWHM A(3) on C a linear background Y=A(4)*X + A(5) C Taken from Bevington, 1969, p.214 C IMPLICIT NONE INTEGER N INTEGER IW INTEGER NTERMS DOUBLE PRECISION X(N) DOUBLE PRECISION A(NTERMS) C C Local variables C DOUBLE PRECISION XX DOUBLE PRECISION Z,Z2 XX=X(IW) FUNCTN=A(5) + A(4)*XX Z=(XX-A(2))/A(3) Z2=Z*Z IF (Z2.LT.200.D0) THEN FUNCTN=FUNCTN + A(1)*DEXP(-Z2/2.0D0) ENDIF 99 END DOUBLE PRECISION FUNCTION FCHISQ(NPTS,Y,WEIGHT,NFREE,ILEF, : IRIG,YFIT) C C Function subprogram to evaluate reduced chi-squared for fit to C data over the range ILEF to IRIG only. C Taken from Bevington, 1969, p. 194 C IMPLICIT NONE INTEGER NPTS INTEGER NFREE INTEGER ILEF INTEGER IRIG DOUBLE PRECISION Y(NPTS) DOUBLE PRECISION WEIGHT(NPTS) DOUBLE PRECISION YFIT(NPTS) C INTEGER I DOUBLE PRECISION CHISQ CHISQ=0.0D0 IF (NFREE.LE.1) THEN FCHISQ=0.0D0 GO TO 99 ENDIF C C Accumulate chi-squared C DO I=ILEF,IRIG,1 CHISQ=CHISQ + ( WEIGHT(I)*((Y(I)-YFIT(I))**2) ) ENDDO C C Divide by number of degrees of freedom C FCHISQ=CHISQ/DBLE(NFREE) 99 END SUBROUTINE FDERIV(N,X,A,NTERMS,IW,DERIV) C C Subroutine subprogram to evaluate the derivatives of a C Gaussian of height A(1), peak position A(2) and FWHM C A(3) on a linear background Y=A(4)*X + A(5) C Taken from Bevington, 1969, p.241 C IMPLICIT NONE INTEGER N INTEGER IW INTEGER NTERMS DOUBLE PRECISION X(N) DOUBLE PRECISION A(NTERMS) DOUBLE PRECISION DERIV(NTERMS) C C Local variables C DOUBLE PRECISION XX DOUBLE PRECISION Z DOUBLE PRECISION Z2 XX=X(IW) Z=(XX-A(2))/A(3) Z2=Z*Z IF (Z2.GT.200.0D0) THEN DERIV(1)=0.0D0 DERIV(2)=0.0D0 DERIV(3)=0.0D0 ELSE DERIV(1)=DEXP(-Z2/2.0D0) DERIV(2)=A(1)*DERIV(1)*Z/A(3) DERIV(3)=DERIV(2)*Z ENDIF DERIV(4)=XX DERIV(5)=1.0D0 99 END SUBROUTINE MATINV(ARRAY,NORDER,DET) C C Subroutine subprogram to invert a symmetric matrix and calculate C its determinant C The input matrix ARRAY is replaced by its inverse C Taken from Bevington, 1969, p.302-3 C IMPLICIT NONE INTEGER NORDER DOUBLE PRECISION ARRAY(NORDER,NORDER),DET INTEGER I,J,K,L,JK(6),IK(6) DOUBLE PRECISION AMAX,SAVE DET=1.0D0 DO K=1,NORDER,1 C C Find largest element in rest of matrix C AMAX=0.0D0 21 DO I=K,NORDER,1 DO J=K,NORDER,1 IF (DABS(AMAX)-DABS(ARRAY(I,J)).LE.0.0D0) THEN AMAX=ARRAY(I,J) IK(K)=I JK(K)=J ENDIF ENDDO ENDDO C C Interchange rows and columns to put AMAX in ARRAY(K,K) C IF (AMAX.EQ.0.0D0) THEN DET=0.0D0 GO TO 999 ENDIF I=IK(K) IF ((I-K).LT.0) THEN GO TO 21 ENDIF IF ((I-K).EQ.0) THEN GO TO 51 ENDIF IF ((I-K).GT.0) THEN GO TO 43 ENDIF 43 DO J=1,NORDER,1 SAVE=ARRAY(K,J) ARRAY(K,J)=ARRAY(I,J) ARRAY(I,J)=-SAVE ENDDO 51 J=JK(K) IF ((J-K).LT.0) THEN GO TO 21 ENDIF IF ((J-K).EQ.0) THEN GO TO 61 ENDIF IF ((J-K).GT.0) THEN GO TO 53 ENDIF 53 DO I=1,NORDER,1 SAVE=ARRAY(I,K) ARRAY(I,K)=ARRAY(I,J) ARRAY(I,J)=-SAVE ENDDO C C Accumulate elements of inverse matrix C 61 DO I=1,NORDER,1 IF ((I-K).NE.0) THEN ARRAY(I,K)=-ARRAY(I,K)/AMAX ENDIF ENDDO 71 DO I=1,NORDER,1 DO J=1,NORDER,1 IF ((I-K).EQ.0) THEN GO TO 80 ENDIF IF ((J-K).EQ.0) THEN GO TO 80 ENDIF ARRAY(I,J)=ARRAY(I,J) + ( ARRAY(I,K)*ARRAY(K,J) ) 80 CONTINUE ENDDO ENDDO 81 DO J=1,NORDER,1 IF ((J-K).NE.0) THEN ARRAY(K,J)=ARRAY(K,J)/AMAX ENDIF 90 CONTINUE ENDDO ARRAY(K,K)=1.0D0/AMAX DET=DET*AMAX ENDDO C C Restore ordering of matrix C DO L=1,NORDER,1 K=NORDER - L + 1 J=IK(K) IF ((J-K).LE.0) THEN GO TO 111 ENDIF DO I=1,NORDER,1 SAVE=ARRAY(I,K) ARRAY(I,K)=-ARRAY(I,J) ARRAY(I,J)=SAVE ENDDO 111 I=JK(K) IF ((I-K).LE.0) THEN GO TO 130 ENDIF DO J=1,NORDER,1 SAVE=ARRAY(K,J) ARRAY(K,J)=-ARRAY(I,J) ARRAY(I,J)=SAVE ENDDO 130 CONTINUE ENDDO 999 END SUBROUTINE COREL(N,Y1,Y2,N2,WORK,YC) C C This subroutine correlates the array Y1 and Y2, C writing the copmplex output as YC C IMPLICIT NONE INTEGER N INTEGER N2 DOUBLE PRECISION WORK(N2) DOUBLE COMPLEX Y1(N) DOUBLE COMPLEX Y2(N) DOUBLE COMPLEX YC(N) C C Local variables C INTEGER I INTEGER NN(1) DOUBLE PRECISION NO2 NN(1)=N ! Length of dimensions of arrays Y1 and Y2 C C Form the Fourrier transforms of the two C input arrays C CALL DFOURT(Y1,NN,1,-1,+1,WORK) CALL DFOURT(Y2,NN,1,-1,+1,WORK) NO2=DBLE(N)/2.0D0 ! Normalization factor for FFT C C Correlate the arrays Y1 and Y2; multiply the FT of C Y1 by the complex conjugate of Y2 C DO I=1,N,1 YC(I)=Y1(I)*DCONJG(Y2(I))/NO2 ENDDO C C Inverse transform to get the cross correlation C CALL DFOURT(YC,NN,1,+1,+1,WORK) 99 END SUBROUTINE FLUESTN(N1,N2,Z,NW,CENT,PFWHM,BY,ZS,POS,FLUX) C C This subroutine returns the estimated flux in the C N1 1-D profiles based on an estimate of the mean C background at distances from the peak position. C The profile flux is summed over an extent POS-WID-1 to C POS+WID+1 where WID is set by the FWHM of the PSF, C scaled by the rebinning factor C IMPLICIT NONE INTEGER N1 ! 1st dimension of the input array Z INTEGER N2 ! 2nd dimension of the input array Z INTEGER NW ! 2nd dimension of the input array CENT INTEGER BY ! Subsampling factor for the PSF relative to the data DOUBLE PRECISION Z(N1,N2) ! Input array of the N1 1-D profiles DOUBLE PRECISION CENT(N1,NW) ! Input array of the accurate positions of the (N1) profiles DOUBLE PRECISION PFWHM(N1) ! Input array of the FWHM of the N1 PSF profiles DOUBLE PRECISION ZS(N2) ! Work array for the 1-D X-sections of Z DOUBLE PRECISION POS(NW) ! Work array for the positions of the profiles DOUBLE PRECISION FLUX(N1,NW) ! Output array of the estimated line fluxes for the NW profiles at each N1 positions C C Local variables C INTEGER I,J,K,L INTEGER IBGW DOUBLE PRECISION BDFAC DOUBLE PRECISION WID DOUBLE PRECISION OFLU BDFAC=1.50D0 ! Factor* FWHM for offset position from peak C C Copy each X section to a 1D array. For each point source C estimate the flux in the source above the local background C DO I=1,N1,1 DO J=1,N2,1 ZS(J)=Z(I,J) ENDDO WID=PFWHM(I)/DBLE(BY) IBGW=INT(WID) DO K=1,NW,1 POS(K)=CENT(I,K) ENDDO DO L=1,NW,1 CALL FLUEST(N2,ZS,NW,POS,L,WID,BDFAC,IBGW,OFLU) FLUX(I,L)=OFLU ENDDO ENDDO 99 END SUBROUTINE FLUEST(N,Y,NS,POS,IS,WID,FAC,IBGW,FLUX) C C This subroutine returns the estimated flux in the C profile in the array Y centred at POS(IS) and of FWHM C WID based on an estimate of the mean background C at a distance of FAC*WID of extent IBGW from both C sides of the peak. The flux is summed over an extent C POS-WID-1 to POS+WID+1 channels C IMPLICIT NONE INTEGER N INTEGER NS INTEGER IS INTEGER IBGW DOUBLE PRECISION Y(N) DOUBLE PRECISION POS(NS) DOUBLE PRECISION WID DOUBLE PRECISION FAC DOUBLE PRECISION FLUX C C Local variables C INTEGER I INTEGER ISUM INTEGER LEFT INTEGER RIGHT DOUBLE PRECISION SUM DOUBLE PRECISION YCON DOUBLE PRECISION SMAX c integer istat c character*128 text C C Estimate mean continuum on both sides of line at positions C FAC*FWHM from peak C SUM=0.0D0 ISUM=0 LEFT=NINT(POS(IS) - (FAC*WID)) IF (LEFT.LT.1) THEN LEFT=2 ENDIF IF (LEFT.GT.N) THEN LEFT=N-1 ENDIF RIGHT=NINT(POS(IS) + (FAC*WID)) IF (RIGHT.LT.1) THEN RIGHT=2 ENDIF IF (RIGHT.GT.N) THEN RIGHT=N-1 ENDIF C C Check that LEFT and RIGHT do not lie within WID of another point C source. If so place background region at FAC*WID from obstracting C point source position C DO I=1,NS,1 IF (I.NE.IS) THEN IF (LEFT.GE.INT(POS(I)-WID).AND.LEFT.LE.NINT(POS(I)+WID)) : THEN LEFT=INT(POS(I) - (FAC*WID)) IF (LEFT.LE.1) THEN LEFT=2 ENDIF ENDIF IF (RIGHT.GE.INT(POS(I)-WID).AND.RIGHT.LE.NINT(POS(I)+WID)) : THEN RIGHT=NINT(POS(I) + (FAC*WID)) IF (RIGHT.GE.N) THEN RIGHT=N-1 ENDIF ENDIF ENDIF ENDDO c WRITE(TEXT,'('' Rev. Left,Right '',2(I5))') LEFT,RIGHT c CALL UMSPUT(TEXT,1,0,ISTAT) C C Form the summed flux over the two offset background regions C and calculate the mean C SUM=0.0D0 ISUM=0 DO I=LEFT,LEFT-IBGW,-1 IF (I.GE.1.AND.I.LE.N) THEN SUM=SUM + Y(I) ISUM=ISUM+1 ENDIF ENDDO DO I=RIGHT,RIGHT+IBGW,1 IF (I.GE.1.AND.I.LE.N) THEN SUM=SUM + Y(I) ISUM=ISUM+1 ENDIF ENDDO IF (ISUM.GT.0) THEN YCON=SUM/DBLE(ISUM) ELSE YCON=0.0D0 ENDIF C C Subtract estimated background from signal and sum flux C over extent NINT(POS - WID)-1 to NINT(POS + WID)+1 C SUM=0.0D0 SMAX=0.0D0 DO I=NINT(POS(IS)-WID)-1,NINT(POS(IS)+WID)+1,1 IF (I.GE.1.AND.I.LE.N) THEN SUM=SUM + (Y(I)-YCON) IF ((Y(I)-YCON).GT.SMAX) THEN SMAX=Y(I)-YCON ENDIF ENDIF ENDDO c WRITE(TEXT,'('' ISUM,SMAX,WID,YCON '',I5,3(1x,F11.4))') ISUM, c :SMAX,WID,YCON c CALL UMSPUT(TEXT,1,0,ISTAT) C C Set extimated flux. If SUM < 0 then set flux to C absolute value of WID*peak value or 1.0, since point C source fluxes must be positive C IF (SUM.GT.0.0D0) THEN FLUX=SUM ELSE FLUX=DABS(SMAX)*WID c WRITE(TEXT,'('' ISUM,SMAX,WID,YCON '',I5,3(1x,F11.4))') ISUM, c :SMAX,WID,YCON c CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF IF (FLUX.LT.1.0D0) THEN FLUX=1.0D0 ENDIF c WRITE(TEXT,'('' IS,POS,LEFT,RIGHT,SMAX,FLUX '',I5,1X,F10.4, c :2(1X,I5),2(1x,F10.2))') IS,POS(IS),LEFT,RIGHT,SMAX,FLUX c CALL UMSPUT(TEXT,1,0,ISTAT) 99 END SUBROUTINE PSFMANY(P1,P2,PSF,PCEN,BY,NS,POS,INTYPE, : RWK,N2,DWK,PSFS,ISHPSF) C C This subroutine extracts, bins and shifts the PSFs for C the NS individual point sources from the subsampled PSF C image. The interpolation is done using the standard IRAF C routine MRIEVL which can handle several interpolation C schemes, defined by the value of the INTYPE parameter. C The output is a 3-D array with each plane of 1st by 2nd C dimension for the PSF spectrum of each point source C shifted by the fractional pixel from P2/2. The integer C pixel shift is stored in ISHPSF C IMPLICIT NONE INTEGER P1 ! 1st dimension of input array PSF INTEGER P2 ! 2nd dimension of input array PSF INTEGER BY ! PSF subsampling factor in Y direction INTEGER NS ! 1st dimension of input array POS INTEGER INTYPE ! Key to type of interpolation for shifting PSF INTEGER N2 ! Dimension of input array DWK INTEGER ISHPSF(NS) ! Array of integer pixel shifts of NS PSF's REAL RWK(P2) ! 1-D work array for storing PSF to be shifted REAL PSFS(P1,P2,NS) ! Output array of the NP shifted and binned 2-D PSF's DOUBLE PRECISION PSF(P1,P2) ! Input array of the subsampled PSF DOUBLE PRECISION PCEN(P1) ! Position of the centre of the subsampled PSF for each X-section DOUBLE PRECISION POS(P1,NS) ! Input array of the centres of the point sources for all X-sections DOUBLE PRECISION DWK(P2) ! Holding array for output planes of PSFS C C Local variables C INTEGER I,J INTEGER JJ INTEGER NY INTEGER ICEN INTEGER N REAL OFFI REAL CEN REAL AI REAL X REAL Y REAL MRIEVL DOUBLE PRECISION T c integer istat c character*128 text C C Calculate the offsets so that the big pixels are C sampled uniformely and not shifted C OFFI=(REAL(BY)-1.0)/2.0/REAL(BY) C C Compute the array of the integer pixel centres of the C shifted PSF's for all the point sources C DO I=1,NS,1 ISHPSF(I)=NINT((POS(P1/2,I)-PCEN(P1/2)) ) ! /REAL(BY)) c write(text,23) I,POS(P1/2,I),PCEN(P1/2),ISHPSF(I) c23 format(' I,PCEN(P1/2),ISHPSF(I) ',I5,2(1x,F10.4),1x,I5) c call umsput(text,1,0,istat) ENDDO C C Work through the number of PSF's and the range of X-sections C (1st dimension of PSF image), computing the shifted and C rebinned 1-D PSF for each separate PSF position and writing C each to the 3rd dimension of PSFS C Y=1.00 NY=1 DO N=1,NS,1 DO I=1,P1,1 CEN=REAL(PCEN(I)) ICEN=INT(PCEN(I)) C C Convert the Ith section of the PSF to real in preparation C for shifting. C Initialize the Ith section of the output array C DO J=1,P2,1 RWK(J)=REAL(PSF(I,J)) DWK(J)=0.0D0 ENDDO DO AI=1.0-OFFI,REAL(P2)+OFFI,1.0/REAL(BY) C C Calculate the shifted X value for the interpolation C to bring it to position ISHPSF(N) C CBest X=REAL(BY)*(AI - REAL(POS(I,N)-PCEN(I))) X=REAL(BY)*(AI-REAL(ICEN+1)) + REAL(CEN) - : REAL(BY)*(REAL(POS(I,N)) - REAL(NINT(POS(I,N)))) IF (X.GE.1.0.AND.X.LE.REAL(P2)) THEN JJ=NINT(AI) + NINT(POS(I,N)) - ICEN-1 - ISHPSF(N) c write(text,45) JJ,X c45 format(' JJ,X ',I5,1x,F10.4) c call umsput(text,1,0,istat) IF (JJ.GE.1.AND.JJ.LE.P2) THEN DWK(JJ)=DWK(JJ) + : DBLE(MRIEVL(X,Y,RWK,P2,NY,P2,INTYPE)) c write(text,90) I,AI,X,JJ,DWK(JJ) c90 format(' I,AI,X,JJ,DWK ',I5,2(1x,F11.5),1x,I5,1x,F11.5) c call umsput(text,1,0,istat) ENDIF ENDIF ENDDO C C Zap any negative values in the I'th X-section of the C N'th plane of the output cube C DO J=1,P2,1 IF (DWK(J).LT.0.0D0) THEN DWK(J)=0.0D0 ENDIF ENDDO C C Normalize the PSF to 1.0 for this X-section C T=0.0D0 DO J=1,P2,1 T=T+DWK(J) ENDDO IF (T.GT.0.0D0) THEN DO J=1,P2,1 PSFS(I,J,N)=REAL(DWK(J)/T) ENDDO ENDIF ENDDO ENDDO 99 END SUBROUTINE EFFPSF(N1,N2,PSF,FB,PSFR,IRL1,IRL2,WT,PSFM, : PSFFFT,WK1,PSFE) C C This subroutine calculates the effective PSF PSFE to C be applied to the restored background. The wavelength C dependent PSF is weighted by the current estimate of the C position independent spectrum of the background FB and C convolved with the Gaussian resolution kernel PSFR C C Lucy routine: EFFPSF C IMPLICIT NONE INTEGER N1 ! X dimension of input PSF image INTEGER N2 ! Y dimension of input PSF image INTEGER IRL1 ! Lower limit for 10 sigma extent of PSFR INTEGER IRL2 ! Upper limit for 10 sigma extent of PSFR DOUBLE PRECISION PSF(N1,N2) ! Input PSF image with PSF centred at channel N2/2 DOUBLE PRECISION FB(N1) ! Array of spectrum of background after integrating over the spatial profile DOUBLE PRECISION PSFR(N2) ! Array of the background resolution Gaussian kernel DOUBLE PRECISION WT ! Sum over iterations of the convergence of restored point source fluxes DOUBLE PRECISION PSFM(N2) ! Array for the 1-D PSF DOUBLE PRECISION PSFFFT(N2*2) ! Array of Fourrier transform of PSFR DOUBLE PRECISION WK1(N2*2) ! Work array for DCONV DOUBLE PRECISION PSFE(N2) ! Array of the effective PSF C C Local variables C INTEGER I,J INTEGER NY INTEGER DIMS(2) INTEGER STAT CHARACTER*80 TEXT NY=1 DIMS(1)=N2 DIMS(2)=NY C C Initialize PSFM C DO I=1,N2,1 PSFM(I)=0.0D0 ENDDO C C Compute the PSF as a weighted mean with wavelength C IF (WT.EQ.0.0D0) THEN WRITE(TEXT,29) 29 FORMAT('! Warning, PSF weighting factor iz zero. : Resetting to 1.0') CALL UMSPUT(TEXT,1,0,STAT) WT=1.0D0 ENDIF DO I=1,N1,1 DO J=IRL1,IRL2,1 PSFM(J)=PSFM(J) + FB(I)*PSF(I,J)/WT ENDDO ENDDO C C Prepare the array PSFR for Fourrier transform and perform FFT C CALL DFILL(PSFR,N2,NY,PSFFFT) CALL DFOURT(PSFFFT,DIMS,1,-1,0,WK1) C C Convolve the weighted PSF with the resolution kernel C to obtain the effective PSF C CALL DCONV(PSFM,N2,NY,WK1,PSFFFT,PSFE,1) 99 END SUBROUTINE TRANS(N1,N2,XYCEN,PHIN,LNODQ,D1,D2,DQ,DQLIM, : IN,SUBPX,PHOUT,DQO) C C This subroutine outputs the shifted version PHOUT C of the 1-D input array PHIN from the observed grid C to the rectified grid. The Y values of the lower edges C of the X mid points of the pixels are given by XYCEN. C The computation is done by subdividing the pixels SUBPX C times. C It is assumed that distortion is small - i.e. the C Jacobian is 1.0. If the input point has bad data quality C the output data quality is set to DQLIM C Lucy routine: TRANS + INDJ C IMPLICIT NONE INTEGER N1 ! X dimension of input image INTEGER N2 ! Y dimension of input image INTEGER D1 ! X dimension of input data quality image INTEGER D2 ! Y dimension of input data quality image INTEGER IN ! The X index value at which PHIN is computed INTEGER DQ(N1,N2) ! Data quality array of input image IMA INTEGER DQLIM ! Limiting data quality value for `good' data INTEGER SUBPX ! No. of subpixels per pixels for interpolation INTEGER DQO(N2) ! Data quality array of output image IMA DOUBLE PRECISION XYCEN(N1,N2) ! Array of Y values of X mid-points of pixels DOUBLE PRECISION PHIN(N2) ! Array of counts per pixel in rectified grid DOUBLE PRECISION PHOUT(N2) ! Array of counts per pixel in distorted grid LOGICAL LNODQ ! TRUE if data quality DQ not set C C Local variables C INTEGER J,K INTEGER JJ DOUBLE PRECISION CF DOUBLE PRECISION DYR DOUBLE PRECISION DYK DOUBLE PRECISION YM c integer stat c character*80 text C C Initialize output arrays C DO J=1,N2,1 PHOUT(J)=0.0D0 DQO(J)=0 ENDDO C C Loop through each rectified pixel calculating the SUBPX C contributions on the distorted grid C DO J=1,N2,1 CF=PHIN(J)/DBLE(SUBPX) IF (J.EQ.N2) THEN DYR=DABS(XYCEN(IN,J)-XYCEN(IN,J-1)) ELSE DYR=DABS(XYCEN(IN,J+1)-XYCEN(IN,J)) ENDIF DYK=DYR/DBLE(SUBPX) DO K=1,SUBPX,1 YM=XYCEN(IN,J) + DBLE(K-1)*DYK + 0.50D0*DYK ! IRAF convention C C Find the index JJ of the pixel in XYCEN that contains C Y coordinate YM C DO JJ=1,N2,1 IF (JJ.NE.N2) THEN IF (YM.GE.XYCEN(IN,JJ).AND.YM.LT.XYCEN(IN,JJ+1)) THEN PHOUT(JJ)=PHOUT(JJ) + CF c write(text,61) IN,J,JJ,YM,xycen(IN,JJ) c61 format(' IN,J,JJ,YM,xycen(IN,JJ) ',3(I5),2(F11.4)) c call umsput(text,1,0,stat) GO TO 50 ENDIF ELSE IF (YM.GE.XYCEN(IN,JJ).AND.YM.LT.XYCEN(IN,JJ)+DYR) THEN PHOUT(JJ)=PHOUT(JJ) + CF c write(text,61) IN,J,JJ,YM,xycen(IN,JJ) c call umsput(text,1,0,stat) GO TO 50 ENDIF ENDIF ENDDO 50 IF (.NOT.LNODQ) THEN IF (DQ(IN,J).GE.DQLIM) THEN DQO(JJ)=2*DQLIM ENDIF ENDIF ENDDO ENDDO 99 END SUBROUTINE PRFLB(P1,P2,PSF,N1,N2,ARRI,IN,XYCEN,LNODQ, : D1,D2,DQ,DQLIM,SUBPX,IWK1,WK1,WK2, : WK3,WK4,ARRO) C C This subroutine convolves the IN'th section of the C array PSF with the 1-d array ARRI. Data quality is C taken into account if applicable. The output array is C ARRO C IMPLICIT NONE INTEGER P1 ! X dimension of input PSF image INTEGER P2 ! Y dimension of input PSF image INTEGER N1 ! X dimension of input array ARRI INTEGER N2 ! Y dimension of input array ARRI INTEGER D1 ! X dimension of input data quality array DQ INTEGER D2 ! Y dimension of input data quality array DQ INTEGER IN ! Value of X section of input PSF to convolve INTEGER DQ(D1,D2) ! Data quality array for input image ARRI INTEGER DQLIM ! Limiting data quality value for `good' data INTEGER SUBPX ! No. of subpixels per pixels for interpolation in TRANS INTEGER IWK1(N2) ! Holding array for output data quality from TRANS DOUBLE PRECISION PSF(P1,P2) ! Input PSF image DOUBLE PRECISION ARRI(N2) ! Input image to be convolved with PSF DOUBLE PRECISION XYCEN(N1,N2) ! Array of Y values of X mid-points of pixels DOUBLE PRECISION WK1(N2) ! Holding array for 1-D section of PSF DOUBLE PRECISION WK2(N2) ! Holding array for transformed 1-D section of PSF DOUBLE PRECISION WK3(N2*2) ! Array of Fourrier transform of WK1 DOUBLE PRECISION WK4(N2*2) ! Work array required by DCONV DOUBLE PRECISION ARRO(N2) ! Output convolved array LOGICAL LNODQ ! TRUE if the data quality DQ not set C C Local variables C INTEGER I INTEGER NY INTEGER DIMS(2) NY=1 DIMS(1)=N2 DIMS(2)=NY C C Copy the NI'th section of PSF to array WK1 C DO I=1,P2,1 IF (I.LE.N2) THEN WK1(I)=PSF(IN,I) ENDIF ENDDO CALL TRANS(N1,N2,XYCEN,WK1,LNODQ,D1,D2,DQ,DQLIM,IN,SUBPX, : WK2,IWK1) C C Prepare WK2 for Fourier transform and perform FFT C CALL DFILL(WK2,N2,NY,WK3) CALL DFOURT(WK3,DIMS,1,-1,0,WK4) C C Convolve ARRI with WK2 C CALL DCONV(ARRI,N2,NY,WK4,WK3,ARRO,1) 99 END DOUBLE PRECISION FUNCTION GAUSSF(X,SIG) C C This function returns the value for a Gaussian C at the value X and with sigma SIG of unit area C IMPLICIT NONE DOUBLE PRECISION X DOUBLE PRECISION SIG C C Local variables C DOUBLE PRECISION PIE DOUBLE PRECISION FAC DOUBLE PRECISION NORM PARAMETER (PIE=3.14159265358979324D0) FAC=(X*X)/(SIG*SIG) IF (FAC.LT.138.0D0) THEN ! For non-zero limit on e**-(FAC) of 1.0E-30 NORM=DSQRT(2.0D0*PIE)*SIG GAUSSF=DEXP(-0.50D0*FAC)/NORM ELSE GAUSSF=0.0D0 ENDIF 99 END REAL FUNCTION NORDEV(IDUM,MEAN,SIGMA) C C This is the Numerical Recipes function to return a normally C distributed deviate with mean value MEAN and sigma SIGMA, C using RAN1(IDUM) as a source of uniform deviates C IMPLICIT NONE INTEGER IDUM REAL MEAN,SIGMA INTEGER ISET REAL FAC,GSET,RSQ,V1,V2,RAN1,GASDEV SAVE ISET,GSET DATA ISET/0/ IF (ISET.EQ.0) THEN 10 V1=2.*RAN1(IDUM) - 1.0 V2=2.*RAN1(IDUM) - 1.0 RSQ=V1**2 + V2**2 IF (RSQ.GE.1.0.OR.RSQ.EQ.0.0) THEN GO TO 10 ENDIF FAC=SQRT(-2.*LOG(RSQ)/RSQ) GSET=V1*FAC GASDEV=V2*FAC ISET=1 ELSE GASDEV=GSET ISET=0 ENDIF C C GASDEV has zero mean and unit variance - convert C to a value with MEAN and SIGMA C NORDEV=MEAN + SIGMA*GASDEV 99 END REAL FUNCTION RAN1(IDUM) C C This is the Numerical Recipes random number generator C of Park amd Miller with Bays-Durham shuffle safeguards. C IMPLICIT NONE INTEGER IDUM INTEGER IA,IM,IQ,IR,NTAB,NDIV REAL AM,EPS,RNMX PARAMETER (IA=16807,IM=2147483647,AM=1./IM, :IQ=127773,IR=2836,NTAB=32,NDIV=1+(IM-1)/NTAB, :EPS=1.2E-7,RNMX=1.-EPS) INTEGER J,K,IV(NTAB),IY SAVE IV,IY DATA IV /NTAB*0/, IY /0/ IF (IDUM.LE.0.OR.IY.EQ.0) then IDUM=MAX(-IDUM,1) DO 11 J=NTAB+8,1,-1 K=IDUM/IQ IDUM=IA*(IDUM-K*IQ)-IR*K IF (IDUM.LT.0) THEN IDUM=IDUM+IM ENDIF IF (J.LE.NTAB) THEN IV(J)=IDUM ENDIF 11 CONTINUE IY=IV(1) ENDIF K=IDUM/IQ IDUM=IA*(IDUM-K*IQ)-IR*K IF (IDUM.LT.0) THEN IDUM=IDUM+IM ENDIF J=1+IY/NDIV IY=IV(J) IV(J)=IDUM RAN1=MIN(AM*IY,RNMX) 99 END SUBROUTINE FMTCHKA(N,NS,FLUX,OVER,LOVER) C C This subroutine checks all the N*NS values in the array C FLUX. If any element is found with value greater than C OVER, then LOVER is set FALSE C IMPLICIT NONE INTEGER N ! 1st dimension of array FLUX INTEGER NS ! 2nd dimension of array FLUX DOUBLE PRECISION FLUX(N,NS) ! 2-D array of values to be checked DOUBLE PRECISION OVER ! Bound value to check for in FLUX LOGICAL LOVER ! Logical flag: set FALSE if any element of FLUX > OVER C C Local variables C INTEGER I,J LOVER=.TRUE. DO J=1,NS,1 DO I=1,N,1 IF (FLUX(I,J).GT.OVER) THEN LOVER=.FALSE. ENDIF GO TO 99 ENDDO ENDDO 99 END SUBROUTINE DIVID1(ARR1,ARR2,N1,N2,NA,NB,RAT,FLAG) C C This subroutine divides ARR1(NA,NB) by ARR2(NA,NB) C returning the result as RAT. If a divide by zero C is attempted FLAG is returned FALSE C IMPLICIT NONE INTEGER N1 ! X dimension of input images INTEGER N2 ! Y dimension of input images INTEGER NA ! Index X value of element to operate on INTEGER NB ! Index Y value of element to operate on DOUBLE PRECISION ARR1(N1,N2) ! Numerator input array DOUBLE PRECISION ARR2(N1,N2) ! Denominator input array DOUBLE PRECISION RAT ! Output value of ratio LOGICAL FLAG ! Set to FALSE if denominator value is zero FLAG=.TRUE. IF (ARR2(NA,NB).EQ.0.0D0) THEN FLAG=.FALSE. GO TO 99 ELSE RAT=ARR1(NA,NB)/ARR2(NA,NB) ENDIF 99 END SUBROUTINE INTGPTS(P1,P2,NP,PSFFI,ISHPSF,LNODQ,D1,D2,DQ, : DQLIM,N1,N2,XYCEN,SUBPX,FP,PHYT,PHYX, : PHY,DQO) C C This subroutine computes PHY, the sum over N1 wavelengths C of the contribution of all points sources as a function C of the spatial dimension weighted by the amalgmated PSF's C of the NP point sources. The input point source fluxes are C held as FP(N1,NP) and the PSF of the NP'th point source C as PSFFI(P1,2,NP) C C Lucy routine: INTGPTS C IMPLICIT NONE INTEGER P1 ! X dimension of input array PSFFI INTEGER P2 ! Y dimension of input array PSFFI INTEGER NP ! Z dimension of input array PSFFI INTEGER ISHPSF(NP) ! Array of the fiducial pixel for each NP PSF in PSFFI INTEGER D1 ! X dimension of input data quality image DQ INTEGER D2 ! Y dimension of input data quality image DQ INTEGER DQ(D1,D2) ! Data quality array of input image PHT INTEGER DQLIM ! Limiting data quality value for `good' data INTEGER N1 ! X dimension of output array FP INTEGER N2 ! Y dimension of output array FP INTEGER SUBPX ! No. of subpixels per pixels for interpolation in TRANS INTEGER DQO(N2) ! Data quality array of output image PHY REAL PSFFI(P1,P2,NP) ! Cube of NP PSF's shifted to centre of point source DOUBLE PRECISION FP(N1,NP) ! Array of integrated spectra of NP point sources DOUBLE PRECISION XYCEN(N1,N2) ! Array of Y values of X mid-points of pixels DOUBLE PRECISION PHYT(N2) ! Holding array for PHY computation DOUBLE PRECISION PHYX(N2) ! Holding array for PHY computation in TRANS DOUBLE PRECISION PHY(N2) ! Array of wavelength summed contributions of all point sources LOGICAL LNODQ ! TRUE if data quality DQ not set C C Local variables C INTEGER I,J,K INTEGER JJ C C Initialize output array C DO J=1,N2,1 PHY(J)=0.0D0 ENDDO C C Loop through each X section, summing predicted contribution C of point sources to image on rectified grid. Shift stored C PSF to correct Y position of profile C DO I=1,N1,1 DO J=1,N2,1 PHYT(J)=0.0D0 DO K=1,NP,1 JJ=J - ISHPSF(K) ! - 1 ! - (N2/2) + (P2/2) IF (.NOT.LNODQ) THEN IF (DQ(I,J).LT.DQLIM) THEN IF (JJ.GE.1.AND.JJ.LE.P2) THEN PHYT(J)=PHYT(J) + FP(I,K)*DBLE(PSFFI(I,JJ,K)) ENDIF ENDIF ELSE IF (JJ.GE.1.AND.JJ.LE.P2) THEN PHYT(J)=PHYT(J) + FP(I,K)*DBLE(PSFFI(I,JJ,K)) ENDIF ENDIF ENDDO ENDDO C C Transform the array PHYT to PHY on the observed grid C CALL TRANS(N1,N2,XYCEN,PHYT,LNODQ,D1,D2,DQ,DQLIM,I, : SUBPX,PHYX,DQO) C C Copy sum of all point source PSF contributions to output C array DO J=1,N2,1 PHY(J)=PHY(J)+PHYX(J) ENDDO ENDDO 99 END C FFTETC.F - some useful FFT and convolution routines. C C This file contains three subroutines: C C i) An FFT implementation which is fast and flexible C (even though it is 25 years old). It is called DFOURT. C C ii) A simple piece of code which puts an image into the C correct form to be swallowed by DFOURT. It is DFILL. C C iii) A convolution code which uses the former to convolve a C supplied image with the FFT of a supplied PSF. It is C DCONV. C C The authors of the first of these are given in the comment C block at the start of the code. The other two smaller routines C were written by Richard Hook, Autumn 1992. C C Note that all these routines work fully in DOUBLE PRECISION. C C Note - this is a stripped down version, the full one C with different kinds of FFT and single precisio C versions is in 'bigfftetc.f' C CCC Note-------- C Note - as of October 1993 a bug was fixed in the DFILL C routine. This affects the first row and column and normally has C no effect. "Normally" means cases where the PSF is zero at the C edge of the frame. However if you have been using an earlier version C please change to this one. CCC Note-------- C SUBROUTINE DFILL(DATA,NX,NY,OUT) C C Fill up an array in the way required by the DFOURT C FFT routines C C Supplied: C C DATA - a double precision 2d array. C C NX,NY - integers, the array dimensions. C C Returned: C C OUT - a double precision 1d array in the correct form C to be used by the DFOURT FFT routines. The data values C in DATA are simply put into alternating elements and C the rest set to zero. 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 OUT(IP)=DATA(NX,NY) OUT(IP+1)=0.0D0 IP=IP+2 DO I=1,NX-1 OUT(IP)=DATA(I,NY) OUT(IP+1)=0.0D0 IP=IP+2 ENDDO DO J=1,NY-1 OUT(IP)=DATA(NX,J) OUT(IP+1)=0.0D0 IP=IP+2 DO I=1,NX-1 OUT(IP)=DATA(I,J) OUT(IP+1)=0.0D0 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 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 the DFOURT FFT code. C C Supplied: C C DATA - 2d double precision image array of dimensions (NX,NY) C C NX,NY - integers, the dimensions of this array C C WORK - double precision 1d array of dimension (NX*NY*2) C This is used as workspace for the FFTs. C C PSFFFT - double precision 1d array of dimension (NX*NY*2) C This is the FFT of the PSF as created by the routine C DFOURT. C C FLAG - integer (1 or -1). If this is set to 1 a simple C convolution operation will be performed. If -1 C a correlation will be done instead. This is equivalent C to using a PSF rotated by 180 degrees or multiplying the C imaginary part of the PSF FFT by -1. C C Returned: C C OUT - 2d double precision image array of dimensions (NX,NY) C The result of the convolution. C C Richard Hook, ST-ECF, October 1991 C Double precision version, March 1992 C 1D version, August 1992 C Modified to use DFOURT and improved comments, October 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? DOUBLE PRECISION DATA(NX,NY) ! INPUT DATA ARRAY DOUBLE PRECISION OUT(NX,NY) DOUBLE PRECISION PSFFFT(NX*NY*2) ! FFT OF PSF DOUBLE PRECISION WORK(NX*NY*2) ! WORKSPACE FOR FFT OF DATA DOUBLE PRECISION SCRATCH(40960) DOUBLE PRECISION 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 DFOURT(WORK,DIMS,NDIMS,-1,0,SCRATCH) 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 DFOURT(WORK,DIMS,NDIMS,1,1,SCRATCH) 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 DFOURT(DATA,NN,NDIM,ISIGN,IFORM,WORK) C C Double Precision version of FOURT.F C C Richard Hook, ST-ECF, October 1992 C C-------------------Start of original code (modifed to be more readable) C C June 30, 1983 C C Watch out for INTEGER*4 on VAX!! NN, ISIGN and IFORM must all C be dimensioned INTEGER*4 C C The Fast Fourier Transform in USASI basic FORTRAN C ------------------------------------------------- C C TRANSFORM(J1,J2,...) = SUM(DATA(I1,I2,...)*W1**((I1-1)*(J1-1)) C *W2**((I2-1)*(J2-1))*...), C C where I1 and J1 run from 1 to NN(1) and W1=EXP(ISIGN*2*PI* C SQRT(-1)/NN(1)), etc. there is no limit on the dimensionality C (number of subscripts) of the data array. If an inverse C transform (ISIGN=+1) is performed upon an array of transformed C (ISIGN=-1) data, the original data will reappear, C multiplied by NN(1)*NN(2)*... The array of input data may be C real or complex, at the programmers option, with a saving of C up to forty per cent in running time for real over complex. C (for fastest transform of real data, NN(1) should be even.) C the transform values are always complex, and are returned in the C original array of data, replacing the input data. the length C of each dimension of the data array may be any integer. the C program runs faster on composite integers than on primes, and is C particularly fast on numbers rich in factors of two. C C Timing is in fact given by the following formula. Let NTOT be the C total number of points (real or complex) in the data array, that C is, NTOT=NN(1)*NN(2)*... decompose NTOT into its prime factors, C such as 2**K2 * 3**K3 * 5**K5 * ... Let SUM2 be the sum of all C the factors of two in NTOT, that is, SUM2 = 2*K2. Let SUMF be C the sum of all other factors of NTOT, that is, SUMF = 3*K3+5*K5+.. C The time taken by a multidimensional transform on these NTOT data C point add time = six microseconds), T = 3000 + NTOT*(600+40*SUM2+ C is T = T0 + NTOT*(T1+T2*SUM2+T3*SUMF). On the CDC 3300 (floating C 175*SUMF) microseconds on complex data. C C Implementation of the definition by summation will run in a time C proportional to NTOT*(NN(1)+NN(2)+...). For highly composite NTOT C the savings offered by this program can be dramatic. A one-dimen- C sional array 4000 in length will be transformed in 4000*(600+ C 40*(2+2+2+2+2)+175*(5+5+5)) = 14.5 seconds versus about 4000* C 4000*175 = 2800 seconds for the straightforward technique. C C [Note added by Richard Hook, October 1992: as a rough guide ] C [the time taken for a real and imaginary FFT of a 512*512 ] C [array is about 8s on a SPARCStation 2 ] C C The Fast Fourier algorithm places two restrictions upon the C nature of the data beyond the usual restriction that C the data form one cycle of a periodic function. They are-- C 1. The number of input data and the number of transform values C must be the same. C 2. Considering the data to be in the time domain, C they must be equi-spaced at intervals of DT. Further, the trans- C form values, considered to be in frequency space, will be equi- C spaced from 0 to 2*PI*(NN(I)-1)/(NN(I)*DT) at intervals of C 2*PI/(NN(I)*DT) for each dimension of length NN(I). Of course, C dt need not be the same for every dimension. C C The calling sequence is-- C C CALL DFOURT(DATA,NN,NDIM,ISIGN,IFORM,WORK) C C DATA is the array used to hold the real and imaginary parts C of the data on input and the transform values on output. It C is a multidimensional floating point array, with the real and C imaginary parts of a datum stored immediately adjacent in storage C (such as FORTRAN IV places them). The extent of each dimension C is given in the integer array NN, of length NDIM. ISIGN is -1 C to indicate a forward transform (exponential sign is -) and +1 C for an inverse transform (sign is +). IFORM is +1 if the data and C the transform values are complex. It is 0 if the data are real C but the transform values are complex. If it is 0, the imaginary C parts of the data should be set to zero. As explained above, the C transform values are always complex and are stored in array data. C work is an array used for working storage. It is not necessary C if all the dimensions of the data are powers of two. In this case C it may be replaced by 0 in the calling sequence. Thus, use of C powers of two can free a good deal of storage. If any dimension C is not a power of two, this array must be supplied. It is C floating point, one dimensional of length equal to twice the C largest array dimension (i.e., NN(I) ) that is not a power of C two. Therefore, in one dimension for a non power of two, C work occupies as many storage locations as data. If supplied, C work must not be the same array as data. All subscripts of all C arrays begin at one. C C Example 1. Three-dimensional forward Fourier Transform of a C complex array dimensioned 32 by 25 by 13 in FORTRAN IV. C C DOUBLE PRECISION DATA(32,25,13),WORK(50),NN(3) C COMPLEX DATA C DATA NN/32,25,13/ C DO 1 I=1,32 C DO 1 J=1,25 C DO 1 K=1,13 C 1 DATA(I,J,K)=COMPLEX VALUE C CALL DFOURT(DATA,NN,3,-1,1,WORK) C C Example 2. One-dimensional forward transform of a real array of C length 64 in FORTRAN II. C C DOUBLE PRECISION DATA(2,64) C DO 2 I=1,64 C DATA(1,I)=REAL PART C 2 DATA(2,I)=0. C CALL DFOURT(DATA,64,1,-1,0,0) C C There are no error messages or error halts in this program. The C program returns immediately if ndim or any nn(i) is less than one. C C Program by Norman Brenner from the BASIC program by Charles C Rader (both of MIT Lincoln Laboratory). May 1967. The idea C for the digit reversal was suggested by Ralph Alter (also MIT LL). C This is the fastest and most versatile version of the FFT known C to the author. A program called four2 is available that also C performs the fast fourier transform and is written in usasi basic C fortran. It is about one third as long and restricts the C dimensions of the input array (which must be complex) to be powers C of two. Another program, called four1, is one tenth as long and C runs two thirds as fast on a one-dimensional complex array whose C length is a power of two. C C Reference-- C Fast Fourier Transforms for Fun and Profit, W. Gentleman and C G. Sande, 1966 Fall Joint Computer Conference. C C The work reported in this document was performed at Lincoln Lab- C oratory, a center for research operated by Massachusetts Institute C of Technology, with the support of the U.S. Air Force under C contract af 19(628)-5167. C C-------------End of original comment block. C IMPLICIT NONE C DOUBLE PRECISION DATA(1),WORK(1) DOUBLE PRECISION TWOPI,RTHLF,THETA,WR,WI DOUBLE PRECISION W2R,W2I,W3R,W3I,SUMR,SUMI,DIFR,DIFI DOUBLE PRECISION TEMPR,TEMPI,T2R,T2I,T3R,T3I,WTEMP DOUBLE PRECISION T4R,T4I,WSTPR,WSTPI,THETM,WMINR,WMINI DOUBLE PRECISION TWOWR,SR,SI,OLDSR,OLDSI,STMPR,STMPI DOUBLE PRECISION U1R,U1I,U2R,U2I,U3R,U3I,U4I,U4R INTEGER IFACT(32),NN(1) INTEGER ISIGN,IFORM INTEGER NDIM,NTOT,NHALF,IMAX,IMIN,JMAX INTEGER NP1,NP2,I,J,N,M,NTWO,IF,IDIV,IQUOT,IREM,INON2,NON2P INTEGER ICASE,IFMIN,I1RNG,IDIM,NP0,NPREV INTEGER I1,I2,I3,J1,J2,J3,K1,K2,K3,IPAR,KDIF,NP1TW,MMAX INTEGER NP2HF,I1MAX,NWORK,IFP1,IFP2,I2MAX,LMAX,L INTEGER KMIN,KSTEP,K4,J2MAX,JMIN,J3MAX C Initialise constants TWOPI=6.28318530717959D0 RTHLF=DSQRT(2.0D0)/2.0D0 C Start of executable code IF(NDIM-1)920,1,1 1 NTOT=2 DO 2 IDIM=1,NDIM IF(NN(IDIM))920,920,2 2 NTOT=NTOT*NN(IDIM) C C main loop for each dimension C NP1=2 DO 910 IDIM=1,NDIM N=NN(IDIM) NP2=NP1*N IF(N-1)920,900,5 C C is n a power of two and if not, what are its factors C 5 M=N NTWO=NP1 IF=1 IDIV=2 10 IQUOT=M/IDIV IREM=M-IDIV*IQUOT IF(IQUOT-IDIV)50,11,11 11 IF(IREM)20,12,20 12 NTWO=NTWO+NTWO IFACT(IF)=IDIV IF=IF+1 M=IQUOT GO TO 10 20 IDIV=3 INON2=IF 30 IQUOT=M/IDIV IREM=M-IDIV*IQUOT IF(IQUOT-IDIV)60,31,31 31 IF(IREM)40,32,40 32 IFACT(IF)=IDIV IF=IF+1 M=IQUOT GO TO 30 40 IDIV=IDIV+2 GO TO 30 50 INON2=IF IF(IREM)60,51,60 51 NTWO=NTWO+NTWO GO TO 70 60 IFACT(IF)=M 70 NON2P=NP2/NTWO C C Separate four cases-- C 1. Complex transform C 2. Real transform for the 2nd, 3rd, etc. dimension. method-- C transform half the data, supplying the other half by con- C jugate symmetry. C 3. Real transform for the 1st dimension, n odd. method-- C set the imaginary parts to zero. C 4. Real transform for the 1st dimension, n even. method-- C transform a complex array of length n/2 whose real parts C are the even numbered real values and whose imaginary parts C are the odd numbered real values. separate and supply C the second half by conjugate symmetry. C ICASE=1 IFMIN=1 I1RNG=NP1 IF(IDIM-4)74,100,100 74 IF(IFORM)71,71,100 71 ICASE=2 I1RNG=NP0*(1+NPREV/2) IF(IDIM-1)72,72,100 72 ICASE=3 I1RNG=NP1 IF(NTWO-NP1)100,100,73 73 ICASE=4 IFMIN=2 NTWO=NTWO/2 N=N/2 NP2=NP2/2 NTOT=NTOT/2 I=1 DO 80 J=1,NTOT DATA(J)=DATA(I) 80 I=I+2 C C Shuffle data by bit reversal, since n=2**k. as the shuffling C can be done by simple interchange, no working array is needed C 100 IF(NON2P-1)101,101,200 101 NP2HF=NP2/2 J=1 DO 150 I2=1,NP2,NP1 IF(J-I2)121,130,130 121 I1MAX=I2+NP1-2 DO 125 I1=I2,I1MAX,2 DO 125 I3=I1,NTOT,NP2 J3=J+I3-I2 TEMPR=DATA(I3) TEMPI=DATA(I3+1) DATA(I3)=DATA(J3) DATA(I3+1)=DATA(J3+1) DATA(J3)=TEMPR 125 DATA(J3+1)=TEMPI 130 M=NP2HF 140 IF(J-M)150,150,141 141 J=J-M M=M/2 IF(M-NP1)150,140,140 150 J=J+M GO TO 300 C C Shuffle data by digit reversal for general n C 200 NWORK=2*N DO 270 I1=1,NP1,2 DO 270 I3=I1,NTOT,NP2 J=I3 DO 260 I=1,NWORK,2 IF(ICASE-3)210,220,210 210 WORK(I)=DATA(J) WORK(I+1)=DATA(J+1) GO TO 240 220 WORK(I)=DATA(J) WORK(I+1)=0.0D0 240 IFP2=NP2 IF=IFMIN 250 IFP1=IFP2/IFACT(IF) J=J+IFP1 IF(J-I3-IFP2)260,255,255 255 J=J-IFP2 IFP2=IFP1 IF=IF+1 IF(IFP2-NP1)260,260,250 260 CONTINUE I2MAX=I3+NP2-NP1 I=1 DO 270 I2=I3,I2MAX,NP1 DATA(I2)=WORK(I) DATA(I2+1)=WORK(I+1) 270 I=I+2 C C Main loop for factors of two. C W=EXP(ISIGN*2*PI*SQRT(-1)*M/(4*MMAX)). Check for W=ISIGN*SQRT(-1) C and repeat for W=W*(1+ISIGN*SQRT(-1))/SQRT(2). C 300 IF(NTWO-NP1)600,600,305 305 NP1TW=NP1+NP1 IPAR=NTWO/NP1 310 IF(IPAR-2)350,330,320 320 IPAR=IPAR/4 GO TO 310 330 DO 340 I1=1,I1RNG,2 DO 340 K1=I1,NTOT,NP1TW K2=K1+NP1 TEMPR=DATA(K2) TEMPI=DATA(K2+1) DATA(K2)=DATA(K1)-TEMPR DATA(K2+1)=DATA(K1+1)-TEMPI DATA(K1)=DATA(K1)+TEMPR 340 DATA(K1+1)=DATA(K1+1)+TEMPI 350 MMAX=NP1 360 IF(MMAX-NTWO/2)370,600,600 370 LMAX=MAX0(NP1TW,MMAX/2) DO 570 L=NP1,LMAX,NP1TW M=L IF(MMAX-NP1)420,420,380 380 THETA=-TWOPI*DFLOAT(L)/DFLOAT(4*MMAX) IF(ISIGN)400,390,390 390 THETA=-THETA 400 WR=DCOS(THETA) WI=DSIN(THETA) 410 W2R=WR*WR-WI*WI W2I=2.0D0*WR*WI W3R=W2R*WR-W2I*WI W3I=W2R*WI+W2I*WR 420 DO 530 I1=1,I1RNG,2 KMIN=I1+IPAR*M IF(MMAX-NP1)430,430,440 430 KMIN=I1 440 KDIF=IPAR*MMAX 450 KSTEP=4*KDIF IF(KSTEP-NTWO)460,460,530 460 DO 520 K1=KMIN,NTOT,KSTEP K2=K1+KDIF K3=K2+KDIF K4=K3+KDIF IF(MMAX-NP1)470,470,480 470 U1R=DATA(K1)+DATA(K2) U1I=DATA(K1+1)+DATA(K2+1) U2R=DATA(K3)+DATA(K4) U2I=DATA(K3+1)+DATA(K4+1) U3R=DATA(K1)-DATA(K2) U3I=DATA(K1+1)-DATA(K2+1) IF(ISIGN)471,472,472 471 U4R=DATA(K3+1)-DATA(K4+1) U4I=DATA(K4)-DATA(K3) GO TO 510 472 U4R=DATA(K4+1)-DATA(K3+1) U4I=DATA(K3)-DATA(K4) GO TO 510 480 T2R=W2R*DATA(K2)-W2I*DATA(K2+1) T2I=W2R*DATA(K2+1)+W2I*DATA(K2) T3R=WR*DATA(K3)-WI*DATA(K3+1) T3I=WR*DATA(K3+1)+WI*DATA(K3) T4R=W3R*DATA(K4)-W3I*DATA(K4+1) T4I=W3R*DATA(K4+1)+W3I*DATA(K4) U1R=DATA(K1)+T2R U1I=DATA(K1+1)+T2I U2R=T3R+T4R U2I=T3I+T4I U3R=DATA(K1)-T2R U3I=DATA(K1+1)-T2I IF(ISIGN)490,500,500 490 U4R=T3I-T4I U4I=T4R-T3R GO TO 510 500 U4R=T4I-T3I U4I=T3R-T4R 510 DATA(K1)=U1R+U2R DATA(K1+1)=U1I+U2I DATA(K2)=U3R+U4R DATA(K2+1)=U3I+U4I DATA(K3)=U1R-U2R DATA(K3+1)=U1I-U2I DATA(K4)=U3R-U4R 520 DATA(K4+1)=U3I-U4I KDIF=KSTEP KMIN=4*(KMIN-I1)+I1 GO TO 450 530 CONTINUE M=M+LMAX IF(M-MMAX)540,540,570 540 IF(ISIGN)550,560,560 550 TEMPR=WR WR=(WR+WI)*RTHLF WI=(WI-TEMPR)*RTHLF GO TO 410 560 TEMPR=WR WR=(WR-WI)*RTHLF WI=(TEMPR+WI)*RTHLF GO TO 410 570 CONTINUE IPAR=3-IPAR MMAX=MMAX+MMAX GO TO 360 C C Main loop for factors not equal to two. C W=EXP(ISIGN*2*PI*SQRT(-1)*(J1+J2-I3-1)/IFP2) C 600 IF(NON2P-1)700,700,601 601 IFP1=NTWO IF=INON2 610 IFP2=IFACT(IF)*IFP1 THETA=-TWOPI/DFLOAT(IFACT(IF)) IF(ISIGN)612,611,611 611 THETA=-THETA 612 WSTPR=DCOS(THETA) WSTPI=DSIN(THETA) DO 650 J1=1,IFP1,NP1 THETM=-TWOPI*DFLOAT(J1-1)/DFLOAT(IFP2) IF(ISIGN)614,613,613 613 THETM=-THETM 614 WMINR=DCOS(THETM) WMINI=DSIN(THETM) I1MAX=J1+I1RNG-2 DO 650 I1=J1,I1MAX,2 DO 650 I3=I1,NTOT,NP2 I=1 WR=WMINR WI=WMINI J2MAX=I3+IFP2-IFP1 DO 640 J2=I3,J2MAX,IFP1 TWOWR=WR+WR J3MAX=J2+NP2-IFP2 DO 630 J3=J2,J3MAX,IFP2 JMIN=J3-J2+I3 J=JMIN+IFP2-IFP1 SR=DATA(J) SI=DATA(J+1) OLDSR=0.0D0 OLDSI=0.0D0 J=J-IFP1 620 STMPR=SR STMPI=SI SR=TWOWR*SR-OLDSR+DATA(J) SI=TWOWR*SI-OLDSI+DATA(J+1) OLDSR=STMPR OLDSI=STMPI J=J-IFP1 IF(J-JMIN)621,621,620 621 WORK(I)=WR*SR-WI*SI-OLDSR+DATA(J) WORK(I+1)=WI*SR+WR*SI-OLDSI+DATA(J+1) 630 I=I+2 WTEMP=WR*WSTPI WR=WR*WSTPR-WI*WSTPI 640 WI=WI*WSTPR+WTEMP I=1 DO 650 J2=I3,J2MAX,IFP1 J3MAX=J2+NP2-IFP2 DO 650 J3=J2,J3MAX,IFP2 DATA(J3)=WORK(I) DATA(J3+1)=WORK(I+1) 650 I=I+2 IF=IF+1 IFP1=IFP2 IF(IFP1-NP2)610,700,700 C C Complete a real transform in the 1st dimension, N even, by con- C jugate symmetries. C 700 GO TO (900,800,900,701),ICASE 701 NHALF=N N=N+N THETA=-TWOPI/DFLOAT(N) IF(ISIGN)703,702,702 702 THETA=-THETA 703 WSTPR=DCOS(THETA) WSTPI=DSIN(THETA) WR=WSTPR WI=WSTPI IMIN=3 JMIN=2*NHALF-1 GO TO 725 710 J=JMIN DO 720 I=IMIN,NTOT,NP2 SUMR=(DATA(I)+DATA(J))/2.0D0 SUMI=(DATA(I+1)+DATA(J+1))/2.0D0 DIFR=(DATA(I)-DATA(J))/2.0D0 DIFI=(DATA(I+1)-DATA(J+1))/2.0D0 TEMPR=WR*SUMI+WI*DIFR TEMPI=WI*SUMI-WR*DIFR DATA(I)=SUMR+TEMPR DATA(I+1)=DIFI+TEMPI DATA(J)=SUMR-TEMPR DATA(J+1)=-DIFI+TEMPI 720 J=J+NP2 IMIN=IMIN+2 JMIN=JMIN-2 WTEMP=WR*WSTPI WR=WR*WSTPR-WI*WSTPI WI=WI*WSTPR+WTEMP 725 IF(IMIN-JMIN)710,730,740 730 IF(ISIGN)731,740,740 731 DO 735 I=IMIN,NTOT,NP2 735 DATA(I+1)=-DATA(I+1) 740 NP2=NP2+NP2 NTOT=NTOT+NTOT J=NTOT+1 IMAX=NTOT/2+1 745 IMIN=IMAX-2*NHALF I=IMIN GO TO 755 750 DATA(J)=DATA(I) DATA(J+1)=-DATA(I+1) 755 I=I+2 J=J-2 IF(I-IMAX)750,760,760 760 DATA(J)=DATA(IMIN)-DATA(IMIN+1) DATA(J+1)=0. IF(I-J)770,780,780 765 DATA(J)=DATA(I) DATA(J+1)=DATA(I+1) 770 I=I-2 J=J-2 IF(I-IMIN)775,775,765 775 DATA(J)=DATA(IMIN)+DATA(IMIN+1) DATA(J+1)=0. IMAX=IMIN GO TO 745 780 DATA(1)=DATA(1)+DATA(2) DATA(2)=0. GO TO 900 C C Complete a real transform for the 2nd, 3rd, etc. dimension by C conjugate symmetries. C 800 IF(I1RNG-NP1)805,900,900 805 DO 860 I3=1,NTOT,NP2 I2MAX=I3+NP2-NP1 DO 860 I2=I3,I2MAX,NP1 IMAX=I2+NP1-2 IMIN=I2+I1RNG JMAX=2*I3+NP1-IMIN IF(I2-I3)820,820,810 810 JMAX=JMAX+NP2 820 IF(IDIM-2)850,850,830 830 J=JMAX+NP0 DO 840 I=IMIN,IMAX,2 DATA(I)=DATA(J) DATA(I+1)=-DATA(J+1) 840 J=J-2 850 J=JMAX DO 860 I=IMIN,IMAX,NP0 DATA(I)=DATA(J) DATA(I+1)=-DATA(J+1) 860 J=J-NP0 C C End of loop on each dimension C 900 NP0=NP1 NP1=NP2 910 NPREV=N 920 RETURN END SUBROUTINE DIRCON(DATA,NX,NY,PSF,PNX,PNY,OUT,FLAG) C C DIRCON - direct convolution, not using FFTs. C C First simple version, Richard Hook, June 1993 C IMPLICIT NONE INTEGER NX,NY,PNX,PNY,FLAG DOUBLE PRECISION DATA(NX,NY),PSF(PNX,PNY),OUT(NX,NY) C Local variables INTEGER I,J,N,M,IS,JS DOUBLE PRECISION T DO J=PNY/2,NY-PNY/2 JS=J-PNY/2 DO I=PNX/2,NX-PNX/2 IS=I-PNX/2 T=0.0D0 DO N=1,PNY DO M=1,PNX T=T+DATA(IS+M,JS+N)*PSF(M,N) ENDDO ENDDO OUT(I,J)=T ENDDO ENDDO END SUBROUTINE SIRCON(DATA,NX,NY,PSF,PNX,PNY,OUT,FLAG) C C SIRCON - direct convolution, not using FFTs. C Single precision version. C C First simple version, Richard Hook, June 1993 C IMPLICIT NONE INTEGER NX,NY,PNX,PNY,FLAG REAL DATA(NX,NY),PSF(PNX,PNY),OUT(NX,NY) C Local variables INTEGER I,J,N,M,IS,JS REAL T DO J=PNY/2,NY-PNY/2 JS=J-PNY/2 DO I=PNX/2,NY-PNX/2 IS=I-PNX/2 T=0.0 DO N=1,PNY DO M=1,PNX T=T+DATA(IS+M,JS+N)*PSF(M,N) ENDDO ENDDO OUT(I,J)=T ENDDO ENDDO END SUBROUTINE FILCON(DATA,NX,NY,VAL) C C Fill up a double precision array with a constant value C IMPLICIT NONE INTEGER NX,NY,I,J DOUBLE PRECISION DATA(NX,NY),VAL DO J=1,NY DO I=1,NX DATA(I,J)=VAL ENDDO ENDDO RETURN END SUBROUTINE SILCON(DATA,NX,NY,VAL) C C Fill up a single precision array with a constant value C IMPLICIT NONE INTEGER NX,NY,I,J REAL DATA(NX,NY),VAL DO J=1,NY DO I=1,NX DATA(I,J)=VAL ENDDO ENDDO RETURN END