SUBROUTINE SPECPL C C This FORTRAN program extracts any number of 1-D spectra of C point sources from a 2-D image and its statistical error C image using the Lucy-Hook two channel restoration method. C The PSF spectrum image is required with the same input C dimensions as the data. C C The 1-D spectrum and statistical error + data quality are C output as a table. Each separate source has its own table C specified by a root name with an extension _n C The output background spectrum, either with the point C source spectra removed or included, together with error C and data quality image, is also output. C The errors on the extracted spectrum and fitted C background can be formed by Monte Carlo trials. C Data quality is handled through input and output for all C files except the PSF image. C C Written by J. R. Walsh, ST-ECF, ESO (jwalsh@eso.org) September 1998 - C May 2000 C IMPLICIT NONE INTEGER I,J INTEGER STAT INTEGER STAT1 INTEGER STAT2 INTEGER STAT3 INTEGER STAT4 INTEGER DTYPE INTEGER NAXIS1 INTEGER NAXIS2 INTEGER NAXIS3 INTEGER NAXIS4 INTEGER N1 INTEGER N2 INTEGER E1 INTEGER E2 INTEGER D1 INTEGER D2 INTEGER P1 INTEGER P2 INTEGER NELEM INTEGER ILEN INTEGER MEMI(1) INTEGER IMDSCR1 INTEGER IMDSCR2 INTEGER IMDSCR3 INTEGER IMDSCR4 INTEGER IMDSCR5 INTEGER IMDSCR6 INTEGER IMDSCR7 INTEGER TDSCR INTEGER PDSCR INTEGER NST INTEGER NST1 INTEGER BX INTEGER COLID1 INTEGER COLID2 INTEGER COLID3 INTEGER NROWS INTEGER NCOLS INTEGER DIMEN(7) INTEGER DATYP(7) INTEGER CDENT(7) INTEGER ODIMEN(7) INTEGER NITER INTEGER ICSUM INTEGER ICSTEP INTEGER INTYPE INTEGER NMONTE INTEGER ISEED INTEGER BX INTEGER DQLIM INTEGER DQINP INTEGER IMINP INTEGER ERRINP INTEGER DQINP INTEGER PSFINP INTEGER CENCX INTEGER CENCY INTEGER CENM INTEGER CENT1 INTEGER CENT INTEGER NST1 REAL MEMR(1) REAL PIST REAL LAMST REAL DELLAM REAL RCENT REAL OFSTOP REAL RVAR REAL GW REAL AG INTEGER WI1 INTEGER WR1 INTEGER WR2 INTEGER WD1 INTEGER WD2 INTEGER WD3 INTEGER WD4 INTEGER WD5 INTEGER WD6 INTEGER WD7A INTEGER WD7B INTEGER WD7 INTEGER WD8 INTEGER WD9 INTEGER WD10 INTEGER WD11 INTEGER WD12 INTEGER WD13 INTEGER WD14 INTEGER WD15 INTEGER WD16 INTEGER WD17 INTEGER WD18 INTEGER WD19 INTEGER WD20 INTEGER WD21 INTEGER WD22 INTEGER WD23 INTEGER WD24 INTEGER WD25 INTEGER WD26 INTEGER WD27 INTEGER WD28 INTEGER WD29 INTEGER WD30 INTEGER WD31 INTEGER WD32 INTEGER WD33 INTEGER POSPEC INTEGER WAV INTEGER POERR INTEGER PODQ INTEGER POISPEC INTEGER POIERR INTEGER POIDQ INTEGER IMBACK INTEGER IMBERR INTEGER IMBDQ DOUBLE PRECISION MEMD(1) CHARACTER*1 PSMETH CHARACTER*1 POMETH CHARACTER*24 IMANAM CHARACTER*24 ERRNAM CHARACTER*24 DQNAM CHARACTER*24 BAKIMA CHARACTER*24 BAKERR CHARACTER*24 BAKDQ CHARACTER*24 PSFNAM CHARACTER*24 INTERP CHARACTER*24 POITAB CHARACTER*24 ROONAM CHARACTER*32 TABNAM CHARACTER*20 COLNAM(4) CHARACTER*20 COLUNIT(4) CHARACTER*20 COLFMT(4) CHARACTER*132 OUTEXT LOGICAL MEMB(1) LOGICAL LEXIST LOGICAL NFLAG LOGICAL LNOERR LOGICAL LNODQ LOGICAL LMETH LOGICAL LPOBA LOGICAL LVERB LOGICAL TEXIST COMMON/MEM/MEMD EQUIVALENCE (MEMI,MEMR,MEMD,MEMB) C C Get the name of the input 2-D spectral image to be restored C CALL UCLGST('inpima',IMANAM,STAT) C C Open the file and get the dimensions C 10 CALL UIMOPN(IMANAM,1,IMDSCR1,STAT1) IF (STAT1.EQ.0) THEN CALL UIMGID(IMDSCR1,DTYPE,NAXIS1,DIMEN,STAT2) IF (NAXIS1.NE.2) THEN WRITE(OUTEXT,11) 11 FORMAT(' Input image must be 2-dimensional') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF N1=DIMEN(1) N2=DIMEN(2) C C Check that the dimensions are EVEN, if not abort C IF (MOD(N2,2).NE.0) THEN CALL UMSPUT('! Image Y dimension must be even', : 1,0,STAT) GO TO 990 ENDIF C C Allocate dynamic memory for the 2-D input image C of dimensions N1 by N2 (real number) C NELEM=N1*N2 CALL UDMGET(NELEM,6,IMINP,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,12) 12 FORMAT(' Unable to assign memory for internal 2-D array') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF IF (STAT2.EQ.0) THEN CALL UIGS2R(IMDSCR1,1,N1,1,N2,MEMR(IMINP),STAT3) IF (STAT3.NE.0) THEN WRITE(OUTEXT,13) IMANAM 13 FORMAT(' Error reading data from file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF ELSE WRITE(OUTEXT,14) IMANAM 14 FORMAT(' Error reading file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF ELSE WRITE(OUTEXT,15) IMANAM 15 FORMAT(' Failed to open file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Attempt to read the header parameters for the reference pixel, C wavelength start and increment. If succesful create a wavelength C array. If no wavelength array created, use pixels numbers. C CALL UHDGSR(IMDSCR1,'CRPIX1',PIST,STAT2) CALL UHDGSR(IMDSCR1,'CRVAL1',LAMST,STAT2) CALL UHDGSR(IMDSCR1,'CD1_1',DELLAM,STAT3) C C Allocate dynamic memory for the 1-D wavelength array C of length N1 (double precision) C CALL UDMGET(N1,7,WAV,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,16) 16 FORMAT(' Unable to assign memory for internal 1-D array') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL WAVCRE(N1,PIST,LAMST,DELLAM,MEMD(WAV)) C C Get the name of the input 2-D error image C 20 CALL UCLGST('errima',ERRNAM,STAT) IF (ERRNAM(:1).EQ.' ') THEN LNOERR=.TRUE. E1=1 E2=1 GO TO 30 ENDIF LNOERR=.FALSE. C C Open the file and get the dimensions C CALL UIMOPN(ERRNAM,1,IMDSCR2,STAT1) IF (STAT1.EQ.0) THEN CALL UIMGID(IMDSCR2,DTYPE,NAXIS2,DIMEN,STAT2) IF (NAXIS2.NE.2) THEN WRITE(OUTEXT,11) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF E1=DIMEN(1) E2=DIMEN(2) IF (E1.NE.N1.OR.E2.NE.N2) THEN WRITE(OUTEXT,21) ERRNAM 21 FORMAT(' Error image not same size as data image ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Allocate dynamic memory for the 2-D input error image C of dimensions E1 by E2 (real number) C NELEM=E1*E2 CALL UDMGET(NELEM,6,ERRINP,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,12) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF IF (STAT2.EQ.0) THEN CALL UIGS2R(IMDSCR2,1,E1,1,E2,MEMR(ERRINP),STAT3) IF (STAT3.NE.0) THEN WRITE(OUTEXT,13) ERRNAM CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF ELSE WRITE(OUTEXT,14) ERRNAM CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF ELSE WRITE(OUTEXT,15) ERRNAM CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Get the name of the input 2-D data quality image C 30 CALL UCLGST('dqima',DQNAM,STAT) IF (DQNAM(:1).EQ.' ') THEN LNODQ=.TRUE. D1=1 D2=1 DQLIM=2 GO TO 50 ENDIF LNODQ=.FALSE. C C Open the file and get the dimensions C CALL UIMOPN(DQNAM,1,IMDSCR3,STAT1) IF (STAT1.EQ.0) THEN CALL UIMGID(IMDSCR3,DTYPE,NAXIS3,DIMEN,STAT2) IF (NAXIS3.NE.2) THEN WRITE(OUTEXT,11) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF D1=DIMEN(1) D2=DIMEN(2) IF (D1.NE.N1.OR.D2.NE.N2) THEN WRITE(OUTEXT,11) DQNAM 31 FORMAT(' Data quality image not same size as data image ', :A24) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Allocate dynamic memory for the 2-D input data quality image C of dimensions D1 by D2 (integer) C NELEM=D1*D2 CALL UDMGET(NELEM,4,DQINP,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,12) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF IF (STAT2.EQ.0) THEN CALL UIGS2I(IMDSCR3,1,D1,1,D2,MEMI(DQINP),STAT3) IF (STAT3.NE.0) THEN WRITE(OUTEXT,13) DQNAM CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF ELSE WRITE(OUTEXT,14) DQNAM CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF ELSE WRITE(OUTEXT,15) DQNAM CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Get the data quality limit such that values with data quality less C than this value will be considered GOOD. This value is also the C fill value for the point sources table C 50 CALL UCLGSI('dqlim',DQLIM,STAT) C C Get the name of the input 2-D PSF image to match the C spectrum image C CALL UCLGST('psfima',PSFNAM,STAT) C C Open the file and get the dimensions C CALL UIMOPN(PSFNAM,1,IMDSCR4,STAT1) IF (STAT1.EQ.0) THEN CALL UIMGID(IMDSCR4,DTYPE,NAXIS4,DIMEN,STAT2) IF (NAXIS4.NE.2) THEN WRITE(OUTEXT,11) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF P1=DIMEN(1) P2=DIMEN(2) C C Dimensions of PSF image must be identical to spectrum C image C IF (P1.NE.N1.OR.P2.NE.N2) THEN WRITE(OUTEXT,52) 52 FORMAT(' PSF and data images must be same size') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Allocate dynamic memory for the 2-D input error image C of dimensions E1 by E2 (real number) C NELEM=P1*P2 CALL UDMGET(NELEM,6,PSFINP,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,12) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF IF (STAT2.EQ.0) THEN CALL UIGS2R(IMDSCR4,1,P1,1,P2,MEMR(PSFINP),STAT3) IF (STAT3.NE.0) THEN WRITE(OUTEXT,13) PSFNAM CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF ELSE WRITE(OUTEXT,14) PSFNAM CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF ELSE WRITE(OUTEXT,15) PSFNAM CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UIMCLO(IMDSCR4,STAT) C C Get the name of the table file lisiting the reference X C and Y positions and slope wrt X to compute the Y positions C as a function of the X position of the point sources C whose spectra are required C CALL UCLGST('poitab',POITAB,STAT) CALL UTTACC(POITAB,LEXIST,STAT) IF (.NOT.LEXIST) THEN WRITE(OUTEXT,212) POITAB 212 FORMAT(' Table file ',A24,' does not exist') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Open the table file, get the number of rows and C columns C CALL UTTOPN(POITAB,1,PDSCR,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,211) POITAB 211 FORMAT(' Failed to open table file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UTPGTI(PDSCR,22,NCOLS,STAT1) ! No. of columns CALL UTPGTI(PDSCR,21,NROWS,STAT1) ! No. of rows C C Check that number of columns is at least three C (reference X, Y and slope) else exit C IF (NCOLS.LT.3) THEN WRITE(OUTEXT,213) POITAB 213 FORMAT(' Less than 3 columns in table file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Allocate dynamic memory for the 1-D arrays C of point source reference X, reference Y, slope C with respect to X and calculated Y centre for C all X channels (all real) C CALL UDMGET(NROWS,6,CENCX,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,16) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(NROWS,6,CENCY,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,16) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(NROWS,6,CENM,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,16) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(NROWS,6,CENT1,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,16) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UTCNUM(PDSCR,1,COLID1,STAT1) ! Column ID for column 1 CALL UTCNUM(PDSCR,2,COLID2,STAT1) ! Column ID for column 2 CALL UTCNUM(PDSCR,3,COLID3,STAT1) ! Column ID for column 3 C C Read each row value of first column into array CENCX C J=1 DO I=1,NROWS,1 CALL UTRGTR(PDSCR,COLID1,1,I,RCENT,NFLAG,STAT1) IF (STAT1.EQ.0) THEN IF (.NOT.NFLAG) THEN CALL TCOPUT(RCENT,NROWS,J,MEMR(CENCX)) J=J+1 ENDIF ENDIF ENDDO C C Read each row value of first column into array CENCY C J=1 DO I=1,NROWS,1 CALL UTRGTR(PDSCR,COLID2,1,I,RCENT,NFLAG,STAT1) IF (STAT1.EQ.0) THEN IF (.NOT.NFLAG) THEN CALL TCOPUT(RCENT,NROWS,J,MEMR(CENCY)) J=J+1 ENDIF ENDIF ENDDO C C Read each row value of first column into array CENM C J=1 DO I=1,NROWS,1 CALL UTRGTR(PDSCR,COLID3,1,I,RCENT,NFLAG,STAT1) IF (STAT1.EQ.0) THEN IF (.NOT.NFLAG) THEN CALL TCOPUT(RCENT,NROWS,J,MEMR(CENM)) J=J+1 ENDIF ENDIF ENDDO CALL UTTCLO(PDSCR,STAT) C C Determine the number of point source spectra in the Y C range 2 - N2-2 of the input spectral image, NST, and the C first source fully above the 2nd row, NST1 C CALL POILIM(NROWS,MEMR(CENCX),MEMR(CENCY),MEMR(CENM), : MEMR(CENT1),N1,N2,NST,NST1) IF (NST.EQ.0) THEN WRITE(OUTEXT,217) 217 FORMAT(' No point source spectra in range of input image') CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMFRE(CENT1,6,STAT) C C Allocate space for the array of the centres of the point C sources actually in range of the image for all X channels CALL UDMGET(N1*NST,6,CENT,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,16) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Get the method for fitting the PSF to determine its C centre position C C for centroid; G for Gaussian C CALL UCLGST('psfmeth',PSMETH,STAT) C C Check that value of PSFMETH is C or G C LMETH=.FALSE. IF (PSMETH.EQ.'C'.OR.PSMETH.EQ.'c') THEN LMETH=.TRUE. ENDIF IF (PSMETH.EQ.'G'.OR.PSMETH.EQ.'g') THEN LMETH=.TRUE. ENDIF IF (.NOT.LMETH) THEN CALL UMSPUT('! Invalid value for psfmeth ',1,0,STAT) GO TO 990 ENDIF C C Determine whether the centre positions of the point sources C are to be determined by cross correlation or given exactly by C the values in the table POIFIT (C or E) C CALL UCLGST('posmeth',POMETH,STAT) C C Check that value of POMETH is C or E C LMETH=.FALSE. IF (POMETH.EQ.'C'.OR.POMETH.EQ.'c') THEN LMETH=.TRUE. ENDIF IF (POMETH.EQ.'E'.OR.POMETH.EQ.'e') THEN LMETH=.TRUE. ENDIF IF (.NOT.LMETH) THEN CALL UMSPUT('! Invalid value for posmeth ',1,0,STAT) GO TO 990 ENDIF C C Enter the number of columns to sum for the data array for C computing the cross correlation with the PSF C CALL UCLGSI('icsum',ICSUM,STAT) ICSTEP=ICSUM C C Enter the number of iterations for the Lucy restoration C CALL UCLGSI('niter',NITER,STAT) C C Enter the stopping value for the objective function C CALL UCLGSR('objfun',OFSTOP,STAT) C C Get the interpolation type for the shifted PSF's C CALL UCLGST('interpol',INTERP,STAT) IF (INTERP.EQ.'nearest') THEN INTYPE=1 ELSE IF (INTERP.EQ.'linear') THEN INTYPE=2 ELSE IF (INTERP.EQ.'poly3') THEN INTYPE=3 ELSE IF (INTERP.EQ.'poly5') THEN INTYPE=4 ELSE IF (INTERP.EQ.'spline3') THEN INTYPE=5 ELSE CALL UMSPUT('! Invalid interpolant specified', : 1,0,STAT) GO TO 990 ENDIF C C Get the fractional contribution for the entropy term C CALL UCLGSR('fracent',RVAR,STAT) AG=DBLE(RVAR) C C Get the sub-sampling factor for the PSF image C CALL UCLGSI('psfb',BX,STAT) C C Verify PSF subsampling factor is sensible C IF (BX.LT.1) THEN CALL UMSPUT('! Invalid PSF blocking factors', : 1,0,STAT) GO TO 990 ENDIF C C Get the smoothing kernel sigma for the floating prior C CALL UCLGSR('skernel',RVAR,STAT) GW=DBLE(RVAR) C C Get the number of Monte Carlo trials to perform for C error estimation C CALL UCLGSI('ntrial',NMONTE,STAT) IF (NMONTE.LE.1) THEN NMONTE=1 LNOERR=.TRUE. ! Even if errors available ENDIF IF (LNOERR.AND.(NMONTE.GT.1)) THEN WRITE(OUTEXT,221) 221 FORMAT('! Multiple trials not possible without error array') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Get the seed for the random number generator C CALL UCLGSI('seed',ISEED,STAT) C C Get the switch to determine if the background is output C or the background + point source (TRUE=point sources + C background; FALSE=background only) C CALL UCLGSB('poi_back',LPOBA,STAT) C C Get the reporting level for the running of the C Lucy-Hook two channel restoration C CALL UCLGSB('verbose',LVERB,STAT) C C Allocate dynamic memory for the 1-D arrays required for C storage and workspace for the two channel Lucy Hook restoration C First integer arrays C CALL UDMGET(N2,4,WI1,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,227) 227 FORMAT(' Unable to assign memory for internal 1-D ', : 'integer array') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Now real arrays C CALL UDMGET(N2,6,WR1,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,231) 231 FORMAT(' Unable to assign memory for internal 1-D ', : 'real array') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(P2,6,WR2,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,231) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Now double precision arrays C CALL UDMGET(N2,7,WD1,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) 233 FORMAT(' Unable to assign memory for internal 1-D ', : 'double precision array') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2*2,7,WD2,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2*2,7,WD3,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(NST,7,WD4,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD5,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD6,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2*2,7,WD7A,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2*2,7,WD7B,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(NST,7,WD7,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD8,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD9,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD10,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD11,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(P2,7,WD12,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(NST,7,WD13,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(NST,7,WD14,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(P2,7,WD15,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(P2*2,7,WD16,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD17,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD18,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD19,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD20,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD21,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(P2,7,WD22,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD23,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD24,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(P2*NST,7,WD25,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,235) 235 FORMAT(' Unable to assign memory for internal 2-D ', : 'double precision array') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD26,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2*2,7,WD27,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,235) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD28,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD29,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(NST,7,WD30,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD31,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(NST,7,WD32,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N2,7,WD33,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,233) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Allocate dynamic memory for the output 2-D arrays C of point source flux, error and data quality C 250 NELEM=N1*NST CALL UDMGET(NELEM,6,POSPEC,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,251) 251 FORMAT(' Unable to assign memory for output 2-D ', : 'real array') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(NELEM,6,POERR,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,251) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(NELEM,4,PODQ,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,252) 252 FORMAT(' Unable to assign memory for internal 2-D ', : 'integer array') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Allocate dynamic memory for the output 2-D images C of background flux, error and data quality C 260 NELEM=N1*N2 CALL UDMGET(NELEM,6,IMBACK,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,261) 261 FORMAT(' Unable to assign memory for output 2-D ', : 'real array') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF NELEM=E1*E2 CALL UDMGET(NELEM,6,IMBERR,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,261) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF NELEM=D1*D2 CALL UDMGET(NELEM,4,IMBDQ,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,262) 262 FORMAT(' Unable to assign memory for output 2-D ', : 'integer array') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Call the subroutine to perform the CPLUCY restoration C of the input spectrum y-section by section. At each C set the point source flux into the output spectrum C vector with the errors and data quality C c call umsput(' into popcpl ',1,0,stat) CALL POPCPL(N1,N2,MEMD(WAV),MEMR(IMINP),E1,E2,MEMR(ERRINP), : D1,D2,MEMI(DQINP),LNOERR,LNODQ,MEMR(PSFINP),NROWS, : MEMR(CENCX),MEMR(CENCY),MEMR(CENM),MEMR(CENT),NST,NST1, : PSMETH,POMETH,ICSTEP,NITER,OFSTOP,INTYPE,AG,BX,GW,NMONTE, : ISEED,DQLIM,LVERB,MEMI(WI1),MEMR(WR1),MEMD(WD1),MEMD(WD2), : MEMD(WD3),MEMD(WD4),MEMD(WD5),MEMD(WD6),MEMD(WD7A), : MEMD(WD7B),MEMD(WD7),MEMD(WD8),MEMD(WD9),MEMD(WD10), : MEMD(WD11),MEMD(WD12),MEMD(WD13),MEMD(WD14),MEMD(WD15), : MEMD(WD16),MEMD(WD17),MEMD(WD18),MEMD(WD19),MEMD(WD20), : MEMD(WD21),MEMD(WD22),MEMD(WD23),MEMD(WD24),MEMD(WD25), : MEMD(WD26),MEMR(WR2),MEMD(WD27),MEMD(WD28),MEMD(WD29), : MEMD(WD30),MEMD(WD31),MEMD(WD32),MEMD(WD33), : MEMR(POSPEC),MEMR(POERR),MEMI(PODQ),LPOBA,MEMR(IMBACK), : MEMR(IMBERR),MEMI(IMBDQ)) c call umsput(' Really outo popcpl ',1,0,stat) C C Free up some of the allocated dynamic memory C CALL UDMFRE(IMINP,6,STAT) CALL UDMFRE(ERRINP,6,STAT) CALL UDMFRE(DQINP,4,STAT) CALL UDMFRE(PSFINP,6,STAT) CALL UDMFRE(CENCX,6,STAT) CALL UDMFRE(CENCY,6,STAT) CALL UDMFRE(CENM,6,STAT) CALL UDMFRE(CENT,6,STAT) CALL UDMFRE(WI1,4,STAT) CALL UDMFRE(WR1,6,STAT) CALL UDMFRE(WR2,6,STAT) CALL UDMFRE(WD1,7,STAT) CALL UDMFRE(WD2,7,STAT) CALL UDMFRE(WD3,7,STAT) CALL UDMFRE(WD4,7,STAT) CALL UDMFRE(WD5,7,STAT) CALL UDMFRE(WD6,7,STAT) CALL UDMFRE(WD7,7,STAT) CALL UDMFRE(WD8,7,STAT) CALL UDMFRE(WD9,7,STAT) CALL UDMFRE(WD10,7,STAT) CALL UDMFRE(WD11,7,STAT) CALL UDMFRE(WD12,7,STAT) CALL UDMFRE(WD13,7,STAT) CALL UDMFRE(WD14,7,STAT) CALL UDMFRE(WD15,7,STAT) CALL UDMFRE(WD16,7,STAT) CALL UDMFRE(WD17,7,STAT) CALL UDMFRE(WD18,7,STAT) CALL UDMFRE(WD19,7,STAT) CALL UDMFRE(WD20,7,STAT) CALL UDMFRE(WD21,7,STAT) CALL UDMFRE(WD22,7,STAT) CALL UDMFRE(WD23,7,STAT) CALL UDMFRE(WD24,7,STAT) CALL UDMFRE(WD25,7,STAT) CALL UDMFRE(WD26,7,STAT) CALL UDMFRE(WD27,7,STAT) CALL UDMFRE(WD28,7,STAT) CALL UDMFRE(WD29,7,STAT) CALL UDMFRE(WD30,7,STAT) CALL UDMFRE(WD31,7,STAT) CALL UDMFRE(WD32,7,STAT) CALL UDMFRE(WD33,7,STAT) C C Set the arrays for the column names C COLNAM(1)='WAVELENGTH' COLNAM(2)='FLUX' COLNAM(3)='STATISTICAL_ERROR' COLNAM(4)='DATA_QUALITY' COLUNIT(1)='Angstroms' COLUNIT(2)='Erg/s/cm**2/Angstrom' COLUNIT(3)='Erg/s/cm**2/Angstrom' COLUNIT(4)=' ' COLFMT(1)='25.16g' COLFMT(2)='15.7g' COLFMT(3)='15.7g' COLFMT(4)='I4' DATYP(1)=7 ! Double Precision Wavelength DATYP(2)=6 ! Real flux DATYP(3)=6 ! Real flux error DATYP(4)=4 ! Integer data quality C C Allocate dynamic memory for the three 1-D arrays C for the point source flux, error and data quality C for writing out to the output table columns C CALL UDMGET(N1,6,POISPEC,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,303) 303 FORMAT(' Unable to assign memory for internal 1-D ', : 'real array') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N1,6,POIERR,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,303) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMGET(N1,4,POIDQ,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,305) 305 FORMAT(' Unable to assign memory for internal 1-D ', : 'integer array') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C For each extracted point source spectrum within the range of C the input image write out a table file of the wavelengths, C the total flux in the point source, the error and the data C quality. Keep the same index numbers for the extracted C spectra as in the input table C CALL UCLGST('Rootab',ROONAM,STAT) DO I=1,NST,1 CALL LISNAM(ROONAM,I+NST1-1,TABNAM) C C Open a table file to list the wavelengths, the total flux C in the point source, the error and the data quality C ILEN=INDEX(TABNAM,' ') -1 CALL UTTACC(TABNAM,TEXIST,STAT) IF (TEXIST) THEN WRITE(OUTEXT,307) TABNAM(:ILEN) 307 FORMAT(' Error writing table. ',A,' already exists') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 330 ENDIF CALL UTTINN(TABNAM,TDSCR,STAT1) CALL UTPPTI(TDSCR,3,N1,STAT2) CALL UTCDEF(TDSCR,COLNAM,COLUNIT,COLFMT,DATYP,4,CDENT,STAT3) IF (STAT1.EQ.0.AND.STAT2.EQ.0.AND.STAT3.EQ.0) THEN CALL UTTCRE(TDSCR,STAT) ELSE WRITE(OUTEXT,311) TABNAM 311 FORMAT(' Failed to open table file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 330 ENDIF C C Copy the relevant section of the 2-D arrays of point C source spectrum, error and data quality to 1-D holding C arrays for writing out C CALL OUT2T1(N1,NST,I,MEMR(POSPEC),MEMR(POISPEC), : MEMR(POERR),MEMR(POIERR),MEMI(PODQ), : MEMI(POIDQ)) C C Put the wavelength, flux, statistical error and data quality C values for the extracted point source into the columns of the C table C CALL UTCPTD(TDSCR,CDENT(1),1,N1,MEMD(WAV),STAT1) IF (STAT1.NE.0) THEN WRITE(OUTEXT,313) TABNAM 313 FORMAT(' Failed to copy wavelengths to table file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) ENDIF CALL UTCPTR(TDSCR,CDENT(2),1,N1,MEMR(POISPEC),STAT1) IF (STAT1.NE.0) THEN WRITE(OUTEXT,315) TABNAM 315 FORMAT(' Failed to copy fluxes to table file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) ENDIF CALL UTCPTR(TDSCR,CDENT(3),1,N1,MEMR(POIERR),STAT1) IF (STAT1.NE.0) THEN WRITE(OUTEXT,317) TABNAM 317 FORMAT(' Failed to copy errors to table file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) ENDIF CALL UTCPTI(TDSCR,CDENT(4),1,N1,MEMI(POIDQ),STAT1) IF (STAT1.NE.0) THEN WRITE(OUTEXT,319) TABNAM 319 FORMAT(' Failed to copy data quality to table file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) ENDIF CALL UTTCLO(TDSCR,STAT) 330 CONTINUE ENDDO C C Open a new file for the fitted background image and write IMBACK C into it C 400 CALL UCLGST('bakima',BAKIMA,STAT) ODIMEN(1)=N1 ODIMEN(2)=N2 CALL UIMCRE(BAKIMA,6,2,ODIMEN,IMDSCR5,STAT1) IF (STAT1.EQ.0) THEN CALL UHDCPY(IMDSCR1,IMDSCR5,STAT2) CALL UIPS2R(IMDSCR5,1,N1,1,N2,MEMR(IMBACK),STAT3) IF (STAT3.NE.0) THEN WRITE(OUTEXT,411) BAKIMA 411 FORMAT(' Error writing data to file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UIMCLO(IMDSCR5,STAT4) ELSE WRITE(OUTEXT,412) BAKIMA 412 FORMAT(' Failed to open data file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Open a new file for the fitted background image error map C and write the array IMBERR into it (if input errors) C 430 IF (LNOERR) THEN GO TO 460 ENDIF CALL UCLGST('bakerr',BAKERR,STAT) CALL UIMCRE(BAKERR,6,2,ODIMEN,IMDSCR6,STAT1) IF (STAT1.EQ.0) THEN CALL UHDCPY(IMDSCR2,IMDSCR6,STAT2) CALL UIPS2R(IMDSCR6,1,N1,1,N2,MEMR(IMBERR),STAT3) IF (STAT3.NE.0) THEN WRITE(OUTEXT,431) BAKERR 431 FORMAT(' Error writing error data to file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UIMCLO(IMDSCR6,STAT4) ELSE WRITE(OUTEXT,432) BAKERR 432 FORMAT(' Failed to open error file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Open a new file for the fitted background image data quality C map and write the array IMBDQ into it (if input data quality) C 460 IF (LNODQ) THEN GO TO 900 ENDIF CALL UCLGST('bakdq',BAKDQ,STAT) CALL UIMCRE(BAKDQ,6,2,ODIMEN,IMDSCR7,STAT1) IF (STAT1.EQ.0) THEN CALL UHDCPY(IMDSCR3,IMDSCR7,STAT2) CALL UIPS2I(IMDSCR7,1,N1,1,N2,MEMI(IMBDQ),STAT3) IF (STAT3.NE.0) THEN WRITE(OUTEXT,461) BAKDQ 461 FORMAT(' Error writing data quality to file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UIMCLO(IMDSCR7,STAT4) ELSE WRITE(OUTEXT,462) BAKDQ 462 FORMAT(' Failed to open data quality file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Close the files left open for copying of headers C 900 CALL UIMCLO(IMDSCR1,STAT) IF (.NOT.LNOERR) THEN CALL UIMCLO(IMDSCR2,STAT) ENDIF IF (.NOT.LNODQ) THEN CALL UIMCLO(IMDSCR3,STAT) ENDIF C C Free up the allocated dynamic memory C CALL UDMFRE(POSPEC,6,STAT) CALL UDMFRE(POERR,6,STAT) CALL UDMFRE(PODQ,4,STAT) CALL UDMFRE(POISPEC,6,STAT) CALL UDMFRE(POIERR,6,STAT) CALL UDMFRE(POIDQ,4,STAT) CALL UDMFRE(IMBACK,6,STAT) CALL UDMFRE(IMBERR,6,STAT) CALL UDMFRE(IMBDQ,4,STAT) GO TO 999 990 OUTEXT=' Program failed. No output written' CALL UMSPUT(OUTEXT,1,0,STAT) 999 END SUBROUTINE POPCPL(N1,N2,XIN,IMAG,E1,E2,ERRIM,D1,D2,DQIM, : LNOERR,LNODQ,PSFIM,NS,CENCX,CENCY,CENM, : CENT,NST,NST1,PSFMETH,POSMETH,ICSTEP, : NITER,OFSTOP,INTERP,RAG,BX,RGW,NTRIAL, : IDUM,DQLIM,VERBOSE,YDQ,WORKR1,WORKD1, : WORKD2,WORKD3,WORKD4,WORKD5,WORKD6, : WORKD7,WORKD8,FLUXVAR,PSBVAR,X,Y,YNY, : PSF,POS,FLUX,PSFP,FPFFT,RP,CH,CF,CFP, : SD,PHB,DPSB,RU,PSFS,PSFPNT,PSFR,PSFFFT, : PHS,PHSM,FLUXMT,PSBMT,FLUXO,PSB,SPEC, : SPECERR,SPECDQ,LPOBA,BACK,BACKER,BACKDQ) C C This subroutine copies each y-section from the input image C IMAG to a 1-d array and performs a two channel restoration C to the resulting 1-d image with the background spatial C variation controlled by a resolution kernel. The same C y-section of the PSF image is also copied to a 1-D array. C The position of the centre of the PSF is determined either C by finding the simple centroid or fitting a Gaussian. The C positions of the spectra are then determined by cross- C correlation with the PSF or as given by the value in CENT. C For the cross correlation, if the peak is < 5 sigma, then C 10 y-channels are summed and the cross-correlation repeated; C if no peak still found the 100 y-channels are summed. The C point sources provide one channel for the restoration. The C background over the whole image provides the other channel. C The read-out noise is RON and the ADU/electron conversion C factor is CADU. C The point source flux is returned as a vector SPEC C together with the error, which may be calculated from C multiple trials, and the data quality. C IMPLICIT NONE INTEGER N1 ! 1st dimension of input array IMAG INTEGER N2 ! 2nd dimension of input array IMAG INTEGER E1 ! 1st dimension of input array ERRIM INTEGER E2 ! 2nd dimension of input array ERRIM INTEGER D1 ! 1st dimension of input array DQIM INTEGER D2 ! 2nd dimension of input array DQIM INTEGER NITER ! Number of iterations of the Lucy-Hook restoration INTEGER ICSTEP ! Number of X-sections to add in forming the profile for cross-correlation INTEGER INTERP ! Interpolation code: 1=Nearest; 2=Linear; 3=Poly3; 4=Poly5; 5=Spline3 INTEGER NS ! Number of point source spectra in input arrays INTEGER NST ! Number of point sources to extract in range of IMAG INTEGER NST1 ! First point source in input arrays in range of IMAG INTEGER BX ! Sum-sampling factor (Y direction) of the PSF with respect to the data INTEGER NTRIAL ! Number of Monte Carlo trials for error estimation of restoration INTEGER IDUM ! Value of seed for random number generator INTEGER DQLIM ! Limiting data quality; value < has GOOD quality INTEGER DQIM(D1,D2) ! Input array of data quality for IMAG INTEGER YDQ(D2) ! Holding array for data quality for each X-section INTEGER SPECDQ(N1,NS) ! Output spectra of the NS restored point sources INTEGER BACKDQ(D1,D2) ! Output 2-D spectrum of restored background REAL IMAG(N1,N2) ! Input array of 2D spectral image to be restored REAL ERRIM(E1,E2) ! Input statistical error image corresponding to IMAG REAL PSFIM(N1,N2) ! Input image of Slit Spread Function for IMAG REAL CENCX(NS) ! Array of fiducial X position for calculating point souce positions REAL CENCY(NS) ! Array of Y position for calculating point souce positions REAL CENM(NS) ! Array of slope of Y positions of point sources as function of X REAL CENT(N1,NST) ! Array of derived Y positions of point sources in range of image REAL OFSTOP ! Value of objective function for convergence criterion REAL RAG ! Value of multiplier for background regularization entropic term REAL RGW ! Sigma of gaussian for background smoothing kernel REAL WORKR1(N2) ! Work array REAL PSFR(N2) ! Holding array for PSF for shifting REAL SPEC(N1,NST) ! Output array of the NS extracted spectra REAL SPECERR(N1,NST) ! Output array of the NS extracted spectra statistical error REAL BACK(N1,N2) ! Output restored background image REAL BACKER(E1,E2) ! Output restored background statistical error DOUBLE PRECISION XIN(N1) ! Input array of centres of X channels DOUBLE PRECISION WORKD1(N2) ! Work array for DOUBLE PRECISION WORKD2(N2*2) ! Work array for DOUBLE PRECISION WORKD3(N2*2) ! Work array for DOUBLE PRECISION WORKD4(NST) ! Work array for DOUBLE PRECISION WORKD5(N2) ! Work array for DOUBLE PRECISION WORKD6(N2) ! Work array for DOUBLE PRECISION WORKD7(N2*2) ! Work array for DOUBLE PRECISION WORKD8(N2*2) ! Work array for DOUBLE PRECISION FLUXVAR(NST) DOUBLE PRECISION PSBVAR(N2) DOUBLE PRECISION X(N2) DOUBLE PRECISION Y(N2) DOUBLE PRECISION YNY(N2) DOUBLE PRECISION PSF(N2) DOUBLE PRECISION POS(NS) DOUBLE PRECISION FLUX(NST) DOUBLE PRECISION PSFP(N2) DOUBLE PRECISION FPFFT(N2*2) DOUBLE PRECISION RP(N2) DOUBLE PRECISION CH(N2) DOUBLE PRECISION CF(N2) DOUBLE PRECISION CFP(N2) DOUBLE PRECISION SD(N2) DOUBLE PRECISION PHB(N2) DOUBLE PRECISION DPSB(N2) DOUBLE PRECISION RU(N2) DOUBLE PRECISION PSFS(N2,NST) DOUBLE PRECISION PSFPNT(N2) DOUBLE PRECISION PSFFFT(N2*2) DOUBLE PRECISION PHS(N2) DOUBLE PRECISION PHSM(N2) DOUBLE PRECISION FLUXMT(NST) DOUBLE PRECISION PSBMT(N2) DOUBLE PRECISION FLUXO(NST) DOUBLE PRECISION PSB(N2) CHARACTER*1 PSFMETH ! Key to method to use for PSF finding CHARACTER*1 POSMETH ! Key to method for applying star positions LOGICAL LNOERR ! Flag to indicate errors available (no errors set TRUE) LOGICAL LNODQ ! Flag to indicate data quality available (no DQ set TRUE) LOGICAL VERBOSE ! Flag to indicate whether to report iteration progress LOGICAL LPOBA C C Local variables C INTEGER I,J,II,K,L INTEGER MUNIT INTEGER LMI INTEGER N11 INTEGER IDUM INTEGER NB INTEGER D21 INTEGER DIMS(2) INTEGER NNEG INTEGER BY INTEGER NTR INTEGER IL INTEGER IR INTEGER ISUM INTEGER ICENPSF INTEGER IBGW INTEGER NITPC INTEGER ISTAT REAL VAL REAL RAN1 REAL NORDEV DOUBLE PRECISION AG DOUBLE PRECISION GW DOUBLE PRECISION DOFLIM DOUBLE PRECISION OVERFL DOUBLE PRECISION VMIN DOUBLE PRECISION PSFPOS DOUBLE PRECISION PSFPK DOUBLE PRECISION PSFWID DOUBLE PRECISION PSFCEN DOUBLE PRECISION PSFWHM DOUBLE PRECISION OBACK DOUBLE PRECISION YSUM DOUBLE PRECISION CCSHIFT DOUBLE PRECISION WID DOUBLE PRECISION BDFAC DOUBLE PRECISION OFLU CHARACTER*96 TEXT LOGICAL LTOT LOGICAL LOVER LOGICAL LMONIT LOGICAL LFAIL LOGICAL LWAR AG=DBLE(RAG) GW=DBLE(RGW) DOFLIM=DBLE(OFSTOP) OVERFL=1.0D30 ! If FLUXMT larger than this then multiple trial results not summed c************* cC cC Set up unit number for monitoring iteration progress, the cC X channel to monitor and the number of the star whose flux to be cC monitored cC c MUNIT=16 c OPEN(UNIT=MUNIT,NAME='monit1.dat',TYPE='NEW',FORM='FORMATTED') c WRITE(MUNIT,11) c11 FORMAT(' NIT , Log.Lik , Entrop , Obj Fun , Mean BG , Flux ,', c :' Total RMS') c WRITE(MUNIT,12) RAG,RGW c12 FORMAT(' Fracent ',F7.3,' Skernel ',F7.3) c LMI=128 c**************** N11=1 LMI=0 C C If performing multiple trials initialize random C number generator C IF (NTRIAL.GT.1) THEN VAL=RAN1(-IDUM) ENDIF C C Create an X array with real column numbers (using the IRAF C convention for pixel centres) C DO I=1,N2,1 X(I)=DBLE(I) ENDDO C C For the floating prior get the smoothing kernel sigma C and create an image for it. Then take the FFT C NB=1 D21=1 DIMS(1)=N2 DIMS(2)=NB CALL FGAUSS(PSFP,N2,NB,GW) CALL DFILL(PSFP,N2,NB,FPFFT) CALL DFOURT(FPFFT,DIMS,1,-1,0,WORKD2) C C Call POICENN to set the Y centres of the NST profiles C for all X channels I over the Y channel range 2 to N2-2. C CALL POICENX(NS,CENCX,CENCY,CENM,N1,NST,NST1,CENT) C C For each x-section of the data and PSF 2-D arrays, C copy the values to 1-D arrays. Loop over the number C of Monte Carlo trials. Determine the channel C corresponding to the maximum of the PSF C OBACK=1.0D0 ! Safe value for mean background estimate BY=1 DO I=1,N1,1 C C Initialize the output variance arrays C IF (.NOT.LNOERR) THEN DO K=1,NST,1 FLUXVAR(K)=0.0D0 ENDDO DO K=1,N2,1 PSBVAR(K)=0.0D0 ENDDO ENDIF C C Copy I'th X-section of the PSF array to a 1-D array C DO J=1,N2,1 PSF(J)=DBLE(PSFIM(I,J)) ENDDO C C Check for negative values, if there are any set C them to zero. Check that the total count in the C PSF is normalised to 1.0. If not normalise. If C total count =0 then skip over this X-section C CALL ZAPTOT1(N2,PSF,I,VERBOSE,LTOT) IF (.NOT.LTOT) THEN WRITE(TEXT,'(''! Warning, X-sect '',I6,''; sum of PSF'', : '' zero. No restoration'')') I CALL UMSPUT(TEXT,1,0,ISTAT) DO K=1,NST,1 FLUXO(K)=0.0D0 ENDDO DO J=1,N2,1 PSB(J)=0.0D0 ENDDO GO TO 200 ENDIF C C Determine rough position of peak, peak value and C approximate FWHM of the PSF C CALL PROFPR(N2,X,PSF,PSFPOS,PSFPK,PSFWID) C C Determine exact position of PSF centre using C either simple centroid (PSFMETH=C) or Gaussian fit C (PSFMETH=G) C CALL PSFCNT(N2,X,PSF,PSFPOS,PSFPK,PSFWID,PSFMETH, : WORKD1,WORKD5,PSFCEN,PSFWHM) ICENPSF=NINT(PSFCEN) C C Fill the PSF internal version (for the background C convolutions) with a block-summed version of the input one C CALL PSFBLOK(N2,NB,PSF,BX,BY,PSFCEN,ICENPSF,INTERP, : WORKR1,WORKD1,N2,NB,PSFPNT) c do j=1,n2,1 c write(text,305) I,J,PSFPNT(J) c305 format(' I:J,PSFPNT(J)',2(I5),1x,E12.5) c call umsput(text,1,0,istat) c enddo C C Check for negative values, if there are any set C them to zero. Check that the total count in the C PSF is normalised to 1.0. If not normalise C CALL ZAPTOT1(N2,PSFPNT,I,.FALSE.,LTOT) C C Form the FFT of the block-summed PSF C CALL DFILL(PSFPNT,N2,NB,PSFFFT) CALL DFOURT(PSFFFT,DIMS,1,-1,0,WORKD2) C C Form a real version of the PSF image to be used for C shifting the individual PSF's C CALL COREAL(PSF,PSFR,N2,NB) 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 IF (POSMETH.EQ.'C'.OR.POSMETH.EQ.'c') THEN 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 ENDIF C C If the position of the point sources are given C exactly (POSMETH=e), then set the array of point source C positions POS C IF (POSMETH.EQ.'E'.OR.POSMETH.EQ.'e') THEN DO K=1,NST,1 POS(K)=DBLE(CENT(I,K)) ENDDO ENDIF C C Perform restoration over number of trials C NTR=1 DO L=1,NTRIAL,1 C C Copy I'th X-section of the data array to a 1-D array. C Copy the range of X-sections from I-ICSTEP/2 to I+ICSTEP/2 C for cross correlation. C DO J=1,N2,1 IF (LNOERR) THEN Y(J)=DBLE(IMAG(I,J)) ELSE C C Errors available so multiple trials possible. Set Y value C to data value if first trial, otherwise take a Gauusian C error distributed value. C IF (L.EQ.1) THEN Y(J)=DBLE(IMAG(I,J)) ELSE Y(J)=DBLE(NORDEV(IDUM,IMAG(I,J),ERRIM(I,J))) ENDIF ENDIF C C Set the data qaulity for this X-section, if available C IF (.NOT.LNODQ) THEN YDQ(J)=DQIM(I,J) ENDIF IF (POSMETH.EQ.'C'.OR.POSMETH.EQ.'c') THEN YSUM=0.0D0 ISUM=0 DO II=IL,IR,1 IF (LNODQ) THEN IF (LNOERR) THEN YSUM=YSUM+DBLE(IMAG(II,J)) ISUM=ISUM+1 ELSE IF (L.EQ.1) THEN YSUM=YSUM+DBLE(IMAG(II,J)) ELSE YSUM=YSUM+ : DBLE(NORDEV(IDUM,IMAG(II,J),ERRIM(II,J))) ENDIF ISUM=ISUM+1 ENDIF ELSE c write(text,664) ii,j,dqim(II,J) c664 format('ii,j,dqim(II,J)',3(1x,I5)) c call umsput(text,1,0,istat) IF (DQIM(II,J).LT.DQLIM) THEN IF (LNOERR) THEN YSUM=YSUM+DBLE(IMAG(II,J)) ISUM=ISUM+1 ELSE IF (L.EQ.1) THEN YSUM=YSUM+DBLE(IMAG(II,J)) ELSE YSUM=YSUM+ : DBLE(NORDEV(IDUM,IMAG(II,J),ERRIM(II,J))) ENDIF ISUM=ISUM+1 ENDIF ENDIF ENDIF ENDDO IF (ISUM.GT.0) THEN YNY(J)=YSUM/DBLE(ISUM) ELSE YNY(J)=0.0D0 ENDIF ENDIF ENDDO C C Check that the data array Y is non negative. If non-negative C values found set to zero C CALL ZAPNEG(Y,N2,NB,NNEG,VMIN,0.0D0) IF (NNEG.GT.0) THEN IF (VERBOSE) THEN WRITE(TEXT,'(''! Warning, X-sect '',I6,''; data'', : '' image contains'',I7,'' negative values,'', : '' set to zero.'')') I,NNEG CALL UMSPUT(TEXT,1,0,ISTAT) C WRITE(TEXT,'(''--Minimum value: '',D20.10)') VMIN C CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF ENDIF C C Find the total in the data array Y, taking account C of the data quality if applicable, and check greater C than zero. If not skip over this X-section C CALL TOTALDQ(N2,NB,Y,D2,D21,YDQ,LNODQ,DQLIM,YSUM) IF (YSUM.LE.0.0D0) THEN WRITE(TEXT,'(''! Warning, X-sect '',I6,''; sum of data'', : '' zero. No restoration'')') I CALL UMSPUT(TEXT,1,0,ISTAT) DO K=1,NST,1 FLUXO(K)=0.0D0 ENDDO DO J=1,N2,1 PSB(J)=0.0D0 ENDDO GO TO 200 ENDIF C C If POSMETH=C then determine the position of the point C sources by cross correlation with the PSF C IF (POSMETH.EQ.'C'.OR.POSMETH.EQ.'c') THEN NITPC=20 ! Number of iterations of LS Gaussian fitting of Cross-Corr peak DO K=1,NST,1 CALL POSCORR(N2,X,YNY,PSF,DBLE(CENT(I,K)),PSFCEN,PSFWHM, : NITPC,WORKD1,WORKD5,WORKD6,WORKD7,WORKD8, : WORKD2,2*N2,WORKD3,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) POS(K)=DBLE(CENT(I,K)) ELSE POS(K)=CCSHIFT + PSFCEN ENDIF ENDDO ENDIF C C Estimate the fluxes in the lines above the local continuum C WID=PSFWHM/DBLE(BX) ! Estimate of width of profile IBGW=INT(WID) ! No. of channels to sum in forming left and right mean continuum BDFAC=1.5D0 ! Factor * FWHM for ofset position of background from peak DO K=1,NST,1 CALL FLUEST(N2,Y,NST,POS,K,WID,BDFAC,IBGW,OFLU) FLUX(K)=OFLU ENDDO IF (VERBOSE) THEN WRITE(TEXT,'('' '')') CALL UMSPUT(TEXT,1,0,ISTAT) WRITE(TEXT,'(''--Working on X-channel: '',I5)') I CALL UMSPUT(TEXT,1,0,ISTAT) c IF (L.EQ.1) THEN c WRITE(TEXT,'(''--PSF: Cen,FWHM '',2(1x,F9.4))') c : PSFCEN,PSFWHM c CALL UMSPUT(TEXT,1,0,ISTAT) c ENDIF DO K=1,NST,1 WRITE(TEXT,'(''--Spectrum No. '',I4,'' at '', : F9.4,'' flux '',E11.4E2)') K+NST1-1,POS(K),FLUX(K) CALL UMSPUT(TEXT,1,0,ISTAT) ENDDO ENDIF IF (I.EQ.LMI) THEN LMONIT=.TRUE. ELSE LMONIT=.FALSE. ENDIF 100 CALL CPLUCY(N2,Y,D2,D21,YDQ,LNODQ,DQLIM,N2,NITER, : DOFLIM,BX,INTERP,AG,POS,FLUX,YSUM,OBACK, : PSFCEN,ICENPSF,FPFFT,RP,CH,CF,CFP,SD, : DPSB,RU,WORKD1,WORKD2,WORKD4,NST,N11, : N2,PSFS,PSFPNT,PSFR,PSFFFT,PHB,PHS, : PHSM,FLUXMT,PSBMT,VERBOSE,MUNIT, : LMONIT) C C If not performing multiple trials or first of multiple C trials copy point source flux and background holding C arrays to output arrays. Else form variance on point C source flux and background C IF (L.EQ.1) THEN DO K=1,NST,1 FLUXO(K)=FLUXMT(K) ENDDO IF (LPOBA) THEN DO K=1,N2,1 PSB(K)=PHSM(K) ENDDO ELSE DO K=1,N2,1 PSB(K)=PSBMT(K) ENDDO ENDIF NTR=NTR+1 ELSE C C Check all values of FLUXMT aren't overflowed. If all C values are below OVERFL then sum the variance on the C point source fluxes and the background array C CALL FMTCHK(NST,FLUXMT,OVERFL,LOVER) IF (LOVER) THEN DO K=1,NST,1 FLUXVAR(K)=FLUXVAR(K) + (FLUXMT(K)-FLUXO(K))**2. ENDDO DO K=1,N2,1 IF (LPOBA) THEN PSBVAR(K)=PSBVAR(K) + (PSBMT(K)-PHSM(K))**2. ELSE PSBVAR(K)=PSBVAR(K) + (PSBMT(K)-PSB(K))**2. ENDIF ENDDO NTR=NTR+1 ENDIF ENDIF ENDDO C C Copy the flux of the restored point source to the output C spectrum at channel I and the restored background to the C output background array. Copy the input errors and data C quality to the output background arrays if no mutliple C trials, else copy the rms on the mean of the multiple trials C to the errors arrays of the point source flux and the C background C 200 DO K=1,NST,1 SPEC(I,K)=REAL(FLUXO(K)) IF (FLUXO(K).LE.0.0D0) THEN SPECERR(I,K)=0.0 SPECDQ(I,K)=DQLIM ELSE IF (NTR.EQ.2) THEN SPECERR(I,K)=0.0 SPECDQ(I,K)=0 ELSE IF (FLUXVAR(K).LE.0.0D0) THEN SPECERR(I,K)=0.0 SPECDQ(I,K)=DQLIM ELSE SPECERR(I,K)= : REAL(SQRT(FLUXVAR(K)/DBLE(NTR-1))) SPECDQ(I,K)=0 ENDIF ENDIF ENDIF ENDDO DO J=1,N2,1 BACK(I,J)=REAL(PSB(J)) IF (.NOT.LNOERR) THEN BACKER(I,J)=REAL(SQRT(PSBVAR(J)/DBLE(NTR-1))) ENDIF IF (.NOT.LNODQ) THEN BACKDQ(I,J)=DQIM(I,J) ENDIF ENDDO ENDDO CLOSE(UNIT=16) 999 END SUBROUTINE FGAUSS(IMAGE,NX,NY,GW) C C Fill the centre of an empty image with a 2d circular C Gaussian of sigma=GW. We go out as far as 5 times sigma. C IMPLICIT NONE INTEGER NX INTEGER NY DOUBLE PRECISION IMAGE(NX,NY) DOUBLE PRECISION GW C C Local Variables C INTEGER I,J DOUBLE PRECISION R2 DOUBLE PRECISION S5S DOUBLE PRECISION T DOUBLE PRECISION DX S5S=25.0D0*GW*GW T=0.0D0 DO J=1,NY,1 DO I=1,NX,1 DX=DBLE(I-NX/2) R2=DX*DX IF (R2.LT.S5S) THEN IMAGE(I,J)=EXP(-R2/(2.0D0*GW*GW)) T=T+IMAGE(I,J) ELSE IMAGE(I,J)=0.0D0 ENDIF ENDDO ENDDO IF (T.NE.0.0D0) THEN DO J=1,NY,1 DO I=1,NX,1 IMAGE(I,J)=IMAGE(I,J)/T ENDDO ENDDO ENDIF 99 END SUBROUTINE ZAPTOT1(N1,ARRY,IS,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 If LREP is set true then the number of negative values C and totals are reported. C The total is then determined. If the total is not C normalised to 1.0, then the values are renormalised C to 1.0. If the total is zero, then LZERO is set false. C IMPLICIT NONE INTEGER N1 ! 1st dimension of input array ARRY INTEGER IS ! X-section from which ARRY derives, being checked DOUBLE PRECISION ARRY(N1) ! Input array to be zapped and renormalised LOGICAL LREP ! If TRUE then report no. negative values and totals per X-section LOGICAL LZERO ! If sum of ARRY is <=zero, then set FALSE C C Local variables C 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 N11=1 CALL ZAPNEG(ARRY,N1,N11,NNEG,VMIN,0.0D0) IF (NNEG.GT.0) THEN IF (LREP) THEN WRITE(TEXT, : '(''! Warning, X-sect '',I6,''; PSF '', : '' had '',I7,'' negative values, set'', : '' to zero.'')') IS,NNEG CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF ENDIF C C Get the total count in array ARRY. If zero then set LZERO to C FALSE and exit. Check that the total count of ARRY is C normalised to 1.0. If not normalise. C CALL TOTAL(ARRY,N1,N11,PT) IF (PT.LE.0.0D0) THEN IF (LREP) THEN WRITE(TEXT, : '(''! Warning, X-sect '',I6,''; PSF zero '')') IS CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF LZERO=.FALSE. GO TO 99 ENDIF LZERO=.TRUE. IF (DABS(PT-1.0D0).GT.1.0D-7) THEN CALL MULC(ARRY,N1,N11,1.0D0/PT,ARRY) IF (LREP) THEN WRITE(TEXT, : '(''! Warning, X-sect '',I6,''; PSF totals '', : F15.7,'' - renormalising.'')') IS,PT CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF ENDIF 99 END SUBROUTINE TOTALDQ(N1,N2,DATA,D1,D2,DQ,LNODQ,DQLIM,TOT) C C This subroutine adds up all the elements in an image C with data qaulity DQ< DQLIM C IMPLICIT NONE INTEGER N1 ! 1st dimension of input array DATA INTEGER N2 ! 2nd dimension of input array DATA INTEGER D1 ! 1st dimension of input array DQ INTEGER D2 ! 2nd dimension of input array DQ INTEGER DQ(D1,D2) ! Data quality array. A value < DQLIM signifies GOOD data INTEGER DQLIM ! Limiting value for GOOD data DOUBLE PRECISION DATA(N1,N2) DOUBLE PRECISION TOT LOGICAL LNODQ ! Flag to signify if data quality used. Set TRUE if NOT used C C Local variables C INTEGER I,J TOT=0.0D0 DO J=1,N2 DO I=1,N1 IF (.NOT.LNODQ) THEN IF (DQ(I,J).LT.DQLIM) THEN TOT=TOT+DATA(I,J) ENDIF ELSE TOT=TOT+DATA(I,J) ENDIF ENDDO ENDDO 99 END SUBROUTINE PROFPR(N,X,Y,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 DOUBLE PRECISION X(N) DOUBLE PRECISION Y(N) DOUBLE PRECISION POS DOUBLE PRECISION PKHT DOUBLE PRECISION WID C C Local variables C INTEGER I INTEGER IPOS INTEGER LEFT INTEGER RIGHT C C Find the peak value and its position in the array Y C PKHT=-1.0D30 IPOS=N/2 DO I=1,N,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.N) THEN DO I=IPOS+1,N,1 IF (DABS(Y(I)).LT.0.50D0*DABS(PKHT)) THEN RIGHT=I GO TO 50 ENDIF ENDDO ELSE RIGHT=N ENDIF 50 WID=X(RIGHT)-X(LEFT)-1.0D0 IF (WID.LT.1.0D0) THEN WID=1.0D0 ENDIF 99 END SUBROUTINE PSFCNT(N,WAV,ARRY,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 DOUBLE PRECISION WAV(N) DOUBLE PRECISION ARRY(N) DOUBLE PRECISION POS DOUBLE PRECISION PK DOUBLE PRECISION WID DOUBLE PRECISION WHT(N) DOUBLE PRECISION FIT(N) DOUBLE PRECISION CEN DOUBLE PRECISION FWHM 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 PSFBLOK(PX,PY,PSF,BX,BY,CENI,ICENI,ITYPE,WORKR,WORKD, : NX,NY,PSFC) C C This subroutine shifts the input (subsampled by BX) PSF C array PSF to channel PX/2 and blocked-sums the subsampled PSF C to the output array PSFC C IMPLICIT NONE INTEGER PX ! 1st dimension of the input array PIN INTEGER PY ! 2nd dimension of the input array PIN INTEGER BX ! Subsampling factor in 1st dimension of PIN INTEGER BY ! Subsampling factor in 2nd dimension of PIN INTEGER ICENI ! Pixel number of peak of PSF profile INTEGER ITYPE ! Interpolation code: 1=Nearest; 2=Linear; 3=Poly3; 4=Poly5; 5=Spline3 INTEGER NX ! 1st dimension of the output array PSF INTEGER NY ! 2nd dimension of the output array PSF REAL WORKR(PX,PY) ! ! Work array for storing the shifted PSF DOUBLE PRECISION PSF(PX,PY) ! Input array of subsampled PSF DOUBLE PRECISION CENI ! Position of centre of PSF profile DOUBLE PRECISION WORKD(PX,PY) ! Work array for storing the shifted PSF DOUBLE PRECISION PSFC(NX,NY) ! Output array of block-summed PSF C C Local variables C INTEGER I,J,K INTEGER II DOUBLE PRECISION CENO c integer stat c character*128 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 to be centred at A position C so that the block summed PSF is at channel PX/2 C CENO=DBLE(PX/2) CALL SHIFYT(PX,PY,PSF,CENI,CENO,ITYPE,WORKR,WORKD) C C Block sum the data into the larger pixels C DO J=1,PY,1 II=NX/2 - PX/BX/2 + 1 DO I=1,PX/BX,1 DO K=1,BX,1 IF (II.GT.1.AND.II.LT.NX) THEN PSFC(II,J)=PSFC(II,J) + WORKD((I-1)*BX+K,J) ENDIF ENDDO II=II+1 ENDDO ENDDO 99 END SUBROUTINE SHIFYT(P1,P2,PSF,PCEN,CENR,ITYPE,YSR,PSFSH) C C This subroutine shifts each X-section of the array C PSF in the Y direction so that the peak is at position C CEN. 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 REAL YSR(P1) ! Work array required for interpolation function DOUBLE PRECISION PSF(P1,P2) ! Input PSF image DOUBLE PRECISION PCEN ! Position of peak of profile in PSF DOUBLE PRECISION CENR ! Desired centre of shifted PSF array DOUBLE PRECISION PSFSH(P1,P2) ! Output shifted PSF array centred at CEN C C Local variables C INTEGER I,J INTEGER NY REAL X1 REAL Y1 REAL MRIEVL c integer stat c character*80 text NY=1 Y1=1.0 C C Initialize output array C DO J=1,P2,1 DO I=1,P1,1 PSFSH(I,J)=0.0D0 ENDDO ENDDO DO J=1,P2,1 C C Convert the input YSF to real in preparation for shifting C DO I=1,P1,1 YSR(I)=REAL(PSF(I,J)) ENDDO C C Shift the PSF by CENR-CENP C DO I=1,P1,1 X1=REAL(DBLE(I)-(CENR-PCEN)) IF (X1.GE.1.0.AND.X1.LE.REAL(P1)) THEN PSFSH(I,J)=DBLE(MRIEVL(X1,Y1,YSR,P1,NY,P1,ITYPE)) ENDIF ENDDO ENDDO 99 END SUBROUTINE COREAL(IN,OUT,NX,NY) C C Convert an image from double to real C IMPLICIT NONE INTEGER NX INTEGER NY REAL OUT(NX,NY) DOUBLE PRECISION IN(NX,NY) C INTEGER I,J DO J=1,NY,1 DO I=1,NX,1 OUT(I,J)=REAL(IN(I,J)) ENDDO ENDDO END SUBROUTINE POICENX(NW,CENCX,CENCY,CENM,N1,NST,NST1,CEN) C C This subroutine sets the real centres of the NST profiles C for all X channels which lie in the range of the data C defined as the values in the range NST1 to NST1+NST-1 of C the NW input values C IMPLICIT NONE INTEGER NW ! Number of elements in input arrays INTEGER N1 ! Number of X values for centres of profiles INTEGER NST ! Number of profiles within Y channel range 2 to N2-2 INTEGER NST1 ! First valid point source at Y channel >=2 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(N1,NST) ! Output centres of profiles for channel I C C Local variables C INTEGER I,J,JJ C C Compute the Y centre of the profiles for channel I. C Only include values in the range 2 - N2-2 in the array C CEN. Flag the first valid point source above Y channel 2 C DO I=1,N1,1 DO J=1,NST,1 JJ=J+NST1-1 CEN(I,J)=(REAL(I)-CENCX(JJ))*CENM(JJ) + CENCY(JJ) ENDDO ENDDO 99 END SUBROUTINE CPLUCY(NX,PHT,D1,D2,DQ,LNODQ,DQLIM,PX,NITER, : OFLIM,BX,INTERP,AG,POS,FLUX,TSUM,OLBAC, : CENPSF,ICPSF,FPFFT,RP,CH,CF,CFP,SD,DPSB, : RU,WORKD1,WORKD2,CS,NS,NSU,NPX,PSFS, : PSFPNT,PSFR,PSFFFT,PHB,PHS,PHSM,FLUXO,PSB, : VERBOSE,MUNIT,LMONIT) C C Lucy point-sources+smooth extended objects restoration. C Original written by R. Hook, based on ideas and code by C Leon Lucy C C Copied from cplucy.f Version 0.2 and turned into a subroutine C for 1-D data and PSF's by J. R. Walsh, December 1997 C IMPLICIT NONE INTEGER NX ! Dimension of input signal array PHT INTEGER D1 ! 1st dimension of input array DQ INTEGER D2 ! 2nd dimension of input array DQ INTEGER DQ(D1,D2) ! Input data qaulity array: value < DQLIM signifies GOOD INTEGER DQLIM ! Limiting value for GOOD data quality INTEGER PX ! 1st dimension of PSF array INTEGER NITER ! No. of iterations for restoration cycles INTEGER BX ! X subsampling factor for PSF array INTEGER INTERP ! Interpolation code: 1=Nearest; 2=Linear; 3=Poly3; 4=Poly5; 5=Spline3 INTEGER ICPSF ! Pixel number of centre of PSF profile INTEGER NS ! Number of 'stars' (spectra) to fit INTEGER NSU ! Index number of star to write to monitoring output INTEGER NPX ! X-dimension of shifted PSF array INTEGER MUNIT ! Unit number of output file for reporting monitoring of iterations REAL PSFR(PX) ! Real array of PSF image PINPNT DOUBLE PRECISION PHT(NX) ! Input data array to be restored DOUBLE PRECISION OFLIM ! Stopping value for objective function being maximized during iterations DOUBLE PRECISION AG ! Fractional contribution of the entropy term DOUBLE PRECISION POS(NS) ! Input array of positions of point source spectra DOUBLE PRECISION FLUX(NS) ! Input flux of NS point source spectra DOUBLE PRECISION TSUM ! Total flux in input data array PHT DOUBLE PRECISION OLBAC ! Current guesimate of mean background per pixel DOUBLE PRECISION CENPSF ! Position of centre of PSF profile DOUBLE PRECISION FPFFT(NX*2) ! FFT of floating prior array DOUBLE PRECISION RP(NX) ! Scaling array DOUBLE PRECISION CH(NX) ! Array of background convolved with smoothing kernel DOUBLE PRECISION CF(NX) ! Array of convolved correction factor to background DOUBLE PRECISION CFP(NX) ! Increment array DOUBLE PRECISION SD(NX) ! Array of the partial derivative of the correction to the background DOUBLE PRECISION PHB(NX) ! Convolved background array DOUBLE PRECISION DPSB(NX) ! Background correction array DOUBLE PRECISION RU(NX) ! Ratio array DOUBLE PRECISION PSFS(NPX,NS) ! 2-d array of shifted PSF's DOUBLE PRECISION PSFPNT(NX) ! PSF array - block summed to same pixel size as data & centred at NX/2 DOUBLE PRECISION PSFFFT(NX*2) ! FFT of block summed and NX/2 centred PSF (PSFPNT) DOUBLE PRECISION WORKD1(NX) ! Work space array for subroutine PSFBLOK DOUBLE PRECISION WORKD2(NX*2) ! Work space array for subroutine DCONV DOUBLE PRECISION CS(NS) ! Correction factor array for point sources DOUBLE PRECISION PHS(NX) ! Sources array DOUBLE PRECISION PHSM(NX) ! Sources + background array DOUBLE PRECISION FLUXO(NS) ! Output flux of point source spectra DOUBLE PRECISION PSB(NX) ! Output background array LOGICAL LNODQ ! Flag to signify if data quality used. Set TRUE if NOT used LOGICAL VERBOSE ! If true then report all verbose statements LOGICAL LMONIT ! If true then write monitoring quantities to output file C C Local variables C INTEGER I INTEGER ISTAT INTEGER NY INTEGER PY INTEGER N INTEGER NPY INTEGER BY INTEGER NNEG INTEGER IRMAX INTEGER JRMAX INTEGER FLAG DOUBLE PRECISION EBACK DOUBLE PRECISION SMB DOUBLE PRECISION RMSRES DOUBLE PRECISION RESMAX DOUBLE PRECISION SMS DOUBLE PRECISION TT DOUBLE PRECISION XMU DOUBLE PRECISION XMU1 DOUBLE PRECISION XMU2 DOUBLE PRECISION HH DOUBLE PRECISION SS DOUBLE PRECISION QQ DOUBLE PRECISION OFLAST DOUBLE PRECISION VMIN DOUBLE PRECISION LV DOUBLE PRECISION XL CHARACTER*128 TEXT LOGICAL LOFLIM C C Data and PSF only 1-d in X C NY=1 PY=1 NPY=1 BY=1 C C Initialise the maximised quantity (objective function) to C something v. small C LV=-1.0D30 C C Initialize the array of the ratio of the data to convolution C CALL FILCON(RU,NX,NY,0.0D0) C C Generate the shifted and rebinned PSF's for the NS point C sources from the subsampled (real) PSF C CALL EXPSFS1(NX,NY,PSFR,CENPSF,ICPSF,BX,BY,NS,POS,INTERP, : NPX,NPY,NS,PSFS) c do i=1,nx,1 c write(text,300) i,psfs(i,1),psfs(i,2) c300 format(' I,PSFS(I)s',I5,1x,2(1x,E12.5)) c call umsput(text,1,0,istat) c enddo C C Prepare the first estimate - just a flat image with C the flux not in the point soruces C CALL FILFLAT(NX,NY,TSUM,NS,FLUX,OLBAC,PSB,EBACK) IF (VERBOSE) THEN WRITE(TEXT,'( : ''--Setting background in first estimate to: '', : F12.3)') EBACK CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF C C Put the guesses for the star fluxes in the fit array C DO I=1,NS,1 FLUXO(I)=FLUX(I) ENDDO C C Now start the iterative procedure C OFLAST=1.0D20 ! Huge value for beginning increase in objective function LOFLIM=.FALSE. DO N=1,NITER,1 IF(VERBOSE) THEN CALL UMSPUT(' ',1,0,ISTAT) WRITE(TEXT,'(''# Starting iteration '',I5, : ''/'',I5,''-----------'')') N,NITER CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF C C Check the background for non-negativity C CALL ZAPNEG(PSB,NX,NY,NNEG,VMIN,1.0D0) IF (NNEG.GT.0) THEN IF (VERBOSE) THEN WRITE(TEXT,'(''! Background image '', : ''contains '',I7,'' negative values, set to unity'')') : NNEG CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF ENDIF C C Get the total in the background, excluding points with C bad data quality if applicable, and check that background C and star flux correctly normalised C CALL TOTALDQ(NX,NY,PSB,D1,D2,DQ,LNODQ,DQLIM,TT) c write(text,115) n,tt,tsum c115 format(' iter, TT,TSUM',I5,2(1x,F12.5)) c call umsput(text,1,0,istat) IF (VERBOSE) THEN WRITE(TEXT,'(''--Flux distribution, Stars: '', : F6.2,''% Background: '',F6.2,''%.'')') : 100.0D0*(1.0D0-(TT/TSUM)),100.0D0*(TT/TSUM) CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF IF (DABS((TT/TSUM)-1.0D0).GT.1.0D-10) THEN IF (VERBOSE) THEN WRITE(TEXT,'( : ''--Sum of background image is '',D20.10)') TT CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF ENDIF C C Convolve the background image using FFTs C CALL DCONV(PSB,NX,NY,WORKD2,PSFFFT,PHB,1) C C If there is a stars list, create the convolutions directly C IF (NS.GT.0) THEN FLAG=+1 CALL CONVLSDQ(NPX,NS,PSFS,FLUXO,D1,D2,DQ,LNODQ, : DQLIM,FLAG,PHS,CS) ENDIF C C Add the background and stars components C CALL ADDIM(PHB,PHS,NX,NY,PHSM) C C Calculate and report some information on the residuals C between the data and the restoration to assist C closeness of fit assessment, Points with BAD C data quality are not considered, if applicable C CALL RESINF(NX,NY,PHSM,PHT,D1,D2,DQ,LNODQ,DQLIM, : RMSRES,RESMAX,IRMAX,JRMAX) IF (VERBOSE) THEN WRITE(TEXT,'(''--RMS residual: '',D10.4, : '' Max residual of '',D10.4, : '' at ('',I4,'')'')') : RMSRES,RESMAX,IRMAX CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF C C Now convolve the background with the smoothing C kernel to get the default solution C CALL DCONV(PSB,NX,NY,WORKD2,FPFFT,CH,1) C C Calculate the sum of the continuum component, taking C account of the data quality if applicable C CALL TOTALDQ(NX,NY,PSB,D1,D2,DQ,LNODQ,DQLIM,SMB) C C And the stars C SMS=0.0D0 DO I=1,NS,1 SMS=SMS+FLUXO(I) ENDDO c IF (VERBOSE) THEN c WRITE(TEXT,'('' Background sum: '', c : G13.5,'' Stars sum: '',G13.5,'' Total: '',G13.5)') c : SMB,SMS,SMB+SMS c CALL UMSPUT(TEXT,1,0,ISTAT) c ENDIF C C Divide the data by the convolution taking account of data C quality if applicable C CALL DIVIDQ(NX,NY,PHT,PHSM,D1,D2,DQ,LNODQ,DQLIM, : VERBOSE,RU) C C For the floating prior, calculate the RP array C by division of PSB by CH C CALL DIVIDQ(NX,NY,PSB,CH,D1,D2,DQ,LNODQ,DQLIM, : VERBOSE,RP) C C Now do the correlations to get the correction factors C First the FFT one C CALL DCONV(RU,NX,NY,WORKD2,PSFFFT,CF,-1) C C Then the default background image C CALL DCONV(RP,NX,NY,WORKD2,FPFFT,CFP,-1) C C and now the direct one using the small PSFs C FLAG=-1 CALL CONVLSDQ(NPX,NS,PSFS,FLUXO,D1,D2,DQ,LNODQ, : DQLIM,FLAG,RU,CS) C C Calculate the entropy and likelihood and fill the C derivatives array taking account of the data quality C if applicable C CALL GETOF(NX,NY,PSB,PHSM,PHT,CH,D1,D2,DQ,LNODQ, : DQLIM,VERBOSE,SD,CFP,SS,HH) C C Calculate the objective function C QQ=HH+AG*SS C C Report some diagnostics on the Log Likelihood, the Entropy C and the value of the Objective Function C IF (VERBOSE) THEN WRITE(TEXT,'(''--Log.Lik: '',E12.5,'' Ent: '',E12.5, : '' Obj.Func: '',E14.6)') : HH,SS,QQ CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF C C Write monitoring quantities to output file C IF (LMONIT) THEN CALL MONIT(MUNIT,N,HH,SS,QQ,SMB,NX,NS,FLUXO,NSU,RMSRES) ENDIF IF (QQ.LT.LV) THEN IF (VERBOSE) THEN WRITE(TEXT, : '(''! Warning, objective function fell by '',D15.7)') : LV-QQ CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF ELSE IF (VERBOSE.AND.N.GT.1) THEN WRITE(TEXT, : '(''--Objective function increased by '',D15.7)') QQ-LV CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF C C Stop iterations if value of QQ-LV for this and previous C iteration both below OFLIM C IF ((QQ-LV).LT.OFLIM.AND.OFLAST.LT.OFLIM) THEN LOFLIM=.TRUE. ELSE OFLAST=QQ-LV ENDIF ENDIF LV=QQ C C Calculate Lagrange multiplier C XMU1=-TSUM/(SMB+SMS) XMU2=-AG*SS/(SMB+SMS) XMU=XMU1+XMU2 C C Calculate the correction for the background C CALL GTDPSB(NX,NY,PSB,CF,SD,DPSB,AG,XMU) XL=1.0D0 C C Now apply the correction to the background image C CALL BCORRDQ(NX,NY,PSB,DPSB,XL,D1,D2,DQ,LNODQ,DQLIM) C C and finally to the stars C DO I=1,NS,1 FLUXO(I)=FLUXO(I)*(1.0D0 + (XL*(CS(I)+XMU))) IF (VERBOSE) THEN WRITE(TEXT, : '('' Star # '',I4,'' Posn: '',F8.2,'' Flux: '', : G12.4,'' Orig: '',G12.4)') I,POS(I),FLUXO(I),FLUX(I) CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF ENDDO C C Stop iterations if reached limit on value of objective C function C IF (LOFLIM) THEN GO TO 90 ENDIF ENDDO 90 CONTINUE 99 END SUBROUTINE FMTCHK(N,FLUX,OVER,LOVER) C C This subroutine checks all the N 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 DOUBLE PRECISION FLUX(N) DOUBLE PRECISION OVER LOGICAL LOVER C INTEGER I LOVER=.TRUE. DO I=1,N,1 IF (FLUX(I).GT.OVER) THEN LOVER=.FALSE. ENDIF GO TO 99 ENDDO 99 END SUBROUTINE FILFLAT(N1,N2,TOT,NP,FLUX,OBACK,IMOUT,BACK) C C This subroutine fills the array IMOUT with the C estimated background value BACK determined as the C difference between the total count and the sum in C the NS point source fluxes FLUX. If the background C estimate is <= zero then it is replaced by OBACK, C else if this is zero it is set to 1.0 C IMPLICIT NONE INTEGER N1 ! 1st dimension of output array INTEGER N2 ! 2nd dimension of output array INTEGER NP ! Number of points sources (dimension of array FLUX DOUBLE PRECISION TOT ! Total flux (stars+background) in input array FLUX DOUBLE PRECISION FLUX(NP) ! Array of input flux values of stars DOUBLE PRECISION OBACK ! Input/output current background estimate DOUBLE PRECISION IMOUT(N1,N2) ! Output array of set background DOUBLE PRECISION BACK ! Output mean background C C Local variables C INTEGER I,J DOUBLE PRECISION SUMS DOUBLE PRECISION DIFF C C Get total in stars C SUMS=0.0D0 DO I=1,NP,1 SUMS=SUMS+FLUX(I) ENDDO C C Find difference between total and stars and fill C IMOUT with this value. Don't allow background C estimate to be less than 1.0D-30. If background valid C then update OBACK for next X channel C DIFF=TOT-SUMS BACK=DIFF/DBLE(N1*N2) ! Mean background IF (BACK.LT.1.0D0) THEN BACK=OBACK IF (BACK.LT.1.0D-20) THEN BACK=1.0D0 ENDIF ELSE OBACK=BACK ENDIF DO J=1,N2,1 DO I=1,N1,1 IMOUT(I,J)=BACK ENDDO ENDDO 99 END SUBROUTINE EXPSFS1(NX,NY,PSF,CEN,ICEN,BX,BY,NS,POS, : INTYPE,PX,PY,NP,PSFS) C C This subroutine extracts, block averages by factor BX C and shifts the subsampled PSF array PSF, centred at CEN, C for the NS point sources to the positions given by CEN. C 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 array PSFS consists of NS shifted 1-D PSF C array sections C IMPLICIT NONE INTEGER NX ! 1st dimension of input array PSF INTEGER NY ! 2nd dimension of input array PSF INTEGER ICEN ! Integer value of centre of PSF profile INTEGER BX ! PSF subsampling factor in 1st dimension INTEGER BY ! PSF subsampling factor in 2nd dimension INTEGER PX ! 1st dimension of output array PSFS INTEGER PY ! 2nd dimension of output array PSFS INTEGER NP ! 3rd dimension of output array PSFS INTEGER NS ! 1st dimension of input array POS INTEGER INTYPE ! Key to type of interpolation for shifting PSF REAL PSF(NX,NY) ! Input array of the subsampled PSF DOUBLE PRECISION CEN ! Position of the centre of the subsampled PSF DOUBLE PRECISION PSFS(PX,NP) ! Output array of the shifted and binned PSF's DOUBLE PRECISION POS(NS) ! Input array of the centres of the point sources C C Local variables C INTEGER I INTEGER II INTEGER N REAL OFFI REAL AI REAL X REAL Y REAL MRIEVL DOUBLE PRECISION T c integer istat c character*128 text C C Bin the sub-sampled PSF by a factor BX and shift to the C positions of the NS point sources C OFFI=(REAL(BX)-1.0)/2.0/REAL(BX) ! Offset so that the big pixels are sampled uniformely and not shifted DO N=1,NS,1 C C Initialize the output section for the shifted PSF C DO I=1,PX,1 PSFS(I,N)=0.0D0 ENDDO Y=1.0 DO AI=1.0-OFFI,REAL(PX)+OFFI,1.0/REAL(BX) C C Calculate the shifted X value for the interpolation C c X=REAL(BX)*(AI-REAL(ICEN+1)) + REAL(ICEN) - c : REAL(BX)*(REAL(POS(N)) - REAL(NINT(POS(N)))) ! Original X=REAL(BX)*(AI-REAL(ICEN+1)) + REAL(CEN) - ! REAL(ICEN) (was 128.0) : REAL(BX)*(REAL(POS(N)) - REAL(NINT(POS(N)))) ! Original IF (X.GE.1.0.AND.X.LE.REAL(NX)) THEN C C Set the PSF down at the correct channel position for this C point source C II=NINT(AI) + NINT(POS(N))-ICEN-1 IF (II.GE.1.AND.II.LE.PX) THEN PSFS(II,N)=PSFS(II,N) + : DBLE(MRIEVL(X,Y,PSF,NX,NY,NX,INTYPE)) ENDIF ENDIF ENDDO C C Zap any negative values in the N'th section of the C output array C DO I=1,PX,1 IF (PSFS(I,N).LT.0.0D0) THEN PSFS(I,N)=0.0D0 ENDIF ENDDO C C Normalize the the N'th section of the shifted PSF C array C T=0.0D0 DO I=1,PX,1 T=T+PSFS(I,N) ENDDO DO I=1,PX,1 PSFS(I,N)=PSFS(I,N)/T ENDDO ENDDO 99 END SUBROUTINE CONVLS(NX,NP,PSFS,FLUX,FLAG,IMOUT,CS) C C This subroutine convolves (FLAG=+1) or correlates (FLAG=-1) C a set of point source PSF's (PSFS) by an array of their C flux (FLUX). The output array is either the array of the C sum of the flux for the convolved PSFs (IMOUT) or the C correlated flux of each point source (CS) C IMPLICIT NONE INTEGER NX ! 1st dimension of input array PSFS INTEGER NP ! 2nd dimension of input array PSFS (=1st dimension of array FLUX) INTEGER FLAG ! FLAG to indicate action: +1 = convolve; -1 = convolve DOUBLE PRECISION PSFS(NX,NP) ! Input array of the NP shifted 1-D PSFs DOUBLE PRECISION FLUX(NP) ! Input flux of the NP point sources DOUBLE PRECISION IMOUT(NX) ! Output array of the NP flux-convolved PSF's DOUBLE PRECISION CS(NP) ! Output flux of the NP correlated point sources C C Local variables C INTEGER I INTEGER N C C FLAG=1, then convolve the PSF for each point source C with tits flux C IF (FLAG.EQ.1) THEN C C Initialise the output array C DO I=1,NX,1 IMOUT(I)=0.0D0 ENDDO C C Add in the scaled stars C DO N=1,NP,1 DO I=1,NX,1 IMOUT(I)=IMOUT(I) + PSFS(I,N)*FLUX(N) ENDDO ENDDO ELSE C C FLAG=-1, then correlate the PSF of each point source C with the array IMOUT C DO N=1,NP,1 C C Initialise the output array C CS(N)=0.0D0 C C Sum the product of the stars and the image supplied C DO I=1,NX CS(N)=CS(N)+IMOUT(I)*PSFS(I,N) ENDDO ENDDO ENDIF 99 END SUBROUTINE CONVLSDQ(NX,NP,PSFS,FLUX,D1,D2,DQ,LNODQ,DQLIM, : FLAG,IMOUT,CS) C C This subroutine convolves (FLAG=+1) or correlates (FLAG=-1) C a set of point source PSF's (PSFS) by an array of their C flux (FLUX). The output array is either the array of the C sum of the flux for the convolved PSFs (IMOUT) or the C correlated flux of each point source (CS) C IMPLICIT NONE INTEGER NX ! 1st dimension of input array PSFS INTEGER NP ! 2nd dimension of input array PSFS (=1st dimension of array FLUX) INTEGER D1 ! 1st dimension of input array DQ INTEGER D2 ! 2nd dimension of input array DQ INTEGER DQ(D1,D2) ! input data qaulity array: value < DQLIM signifies GOOD INTEGER DQLIM ! Limiting value for GOOD data quality INTEGER FLAG ! FLAG to indicate action: +1 = convolve; -1 = convolve DOUBLE PRECISION PSFS(NX,NP) ! Input array of the NP shifted 1-D PSFs DOUBLE PRECISION FLUX(NP) ! Input flux of the NP point sources DOUBLE PRECISION IMOUT(NX) ! Output array of the NP flux-convolved PSF's DOUBLE PRECISION CS(NP) ! Output flux of the NP correlated point sources LOGICAL LNODQ ! Flag to signify if data quality used. Set TRUE if NOT used C C Local variables C INTEGER I INTEGER NY INTEGER N NY=1 C C FLAG=1, then convolve the PSF for each point source C with its flux C IF (FLAG.EQ.1) THEN C C Initialise the output array C DO I=1,NX,1 IMOUT(I)=0.0D0 ENDDO C C Add in the scaled stars C DO N=1,NP,1 DO I=1,NX,1 IF (.NOT.LNODQ) THEN IF (DQ(I,NY).LT.DQLIM) THEN IMOUT(I)=IMOUT(I) + PSFS(I,N)*FLUX(N) ENDIF ELSE IMOUT(I)=IMOUT(I) + PSFS(I,N)*FLUX(N) ENDIF ENDDO ENDDO ELSE C C FLAG=-1, then correlate the PSF of each point source C with the array IMOUT C DO N=1,NP,1 C C Initialise the output array C CS(N)=0.0D0 C C Sum the product of the stars and the image supplied C DO I=1,NX,1 IF (.NOT.LNODQ) THEN IF (DQ(I,NY).LT.DQLIM) THEN CS(N)=CS(N)+IMOUT(I)*PSFS(I,N) ENDIF ELSE CS(N)=CS(N)+IMOUT(I)*PSFS(I,N) ENDIF ENDDO ENDDO ENDIF 99 END SUBROUTINE ADDIM(A,B,NX,NY,OUT) C C Add an image to another C IMPLICIT NONE INTEGER NX INTEGER NY DOUBLE PRECISION A(NX,NY) DOUBLE PRECISION B(NX,NY) DOUBLE PRECISION OUT(NX,NY) INTEGER I,J DO J=1,NY DO I=1,NX OUT(I,J)=A(I,J)+B(I,J) ENDDO ENDDO END SUBROUTINE RESINF(N1,N2,IM1,IM2,D1,D2,DQ,LNODQ,DQLIM, : RMSRES,RESMAX,IRMAX,JRMAX) C C This subroutine calculates the RMS residual between C the two images IM1 and IM2, taking account of the data C quality, if applicable, and the also reports the C largest residual and its pixel position C IMPLICIT NONE INTEGER N1 ! 1st dimension of input array IM1 INTEGER N2 ! 2nd dimension of input array IM1 INTEGER D1 ! 1st dimension of input array DQ INTEGER D2 ! 2nd dimension of input array DQ INTEGER DQ(D1,D2) ! Data quality array. A value < DQLIM signifies GOOD data INTEGER DQLIM ! Limiting value for GOOD data INTEGER IRMAX ! X pixel position of maximum residual in array IM1 INTEGER JRMAX ! Y pixel position of maximum residual in array IM1 DOUBLE PRECISION IM1(N1,N2) ! First input array to be compared DOUBLE PRECISION IM2(N1,N2) ! Second input array to be compared DOUBLE PRECISION RMSRES ! Output root mean square between IM1 and IM2 DOUBLE PRECISION RESMAX ! Maximum residual between IM1 and IM2 LOGICAL LNODQ ! Flag to signify if data quality used. Set TRUE if NOT used C C Local variables C INTEGER I,J INTEGER ISUM DOUBLE PRECISION RES DOUBLE PRECISION T T=0.0D0 RESMAX=0.0D0 ISUM=0 DO J=1,N2 DO I=1,N1 IF (.NOT.LNODQ) THEN IF (DQ(I,J).LT.DQLIM) THEN RES=IM1(I,J)-IM2(I,J) T=T+RES*RES ISUM=ISUM+1 IF (DABS(RES).GT.DABS(RESMAX)) THEN RESMAX=RES IRMAX=I JRMAX=J ENDIF ENDIF ELSE RES=IM1(I,J)-IM2(I,J) T=T+RES*RES ISUM=ISUM+1 IF (DABS(RES).GT.DABS(RESMAX)) THEN RESMAX=RES IRMAX=I JRMAX=J ENDIF ENDIF ENDDO ENDDO IF (ISUM.GT.1) THEN RMSRES=DSQRT(T/DBLE(ISUM)) ENDIF 99 END SUBROUTINE DIVIDQ(N1,N2,IN1,IN2,D1,D2,DQ,LNODQ,DQLIM, : VERBOSE,OUT) C C This subroutine divides array IN1 by IN2 taking account of C the data quality, if applicable. C If the denominator is zero the output is set to zero. C The number of attempts to divide by zero is reported C IMPLICIT NONE INTEGER N1 ! 1st dimension of input array IM1 INTEGER N2 ! 2nd dimension of input array IM1 INTEGER D1 ! 1st dimension of input array DQ INTEGER D2 ! 2nd dimension of input array DQ INTEGER DQ(D1,D2) ! Data quality array. A value < DQLIM signifies GOOD data INTEGER DQLIM ! Limiting value for GOOD data DOUBLE PRECISION IN1(N1,N2) ! First input array to be divided (numerator) DOUBLE PRECISION IN2(N1,N2) ! Second input array to be divided (denominator) DOUBLE PRECISION OUT(N1,N2) ! Output array of ratio IN1/IN2 LOGICAL LNODQ ! Flag to signify if data quality used. Set TRUE if NOT used LOGICAL VERBOSE ! If true then report all verbose statements C C Local variables C INTEGER I,J INTEGER NZERO INTEGER ISTAT DOUBLE PRECISION LAST CHARACTER*80 TEXT NZERO=0 C C Divide array IN1 by IN2 taking account of the data quality C if appropriate. If data quality set then set bad pixel C values of the output array to the last good value rather C than zero. Monitor number of attempts to devide by zero. C LAST=1.0D0 DO J=1,N2 DO I=1,N1 IF (.NOT.LNODQ) THEN IF (DQ(I,J).LT.DQLIM) THEN IF (IN2(I,J).EQ.0.0D0) THEN OUT(I,J)=0.0D0 NZERO=NZERO+1 ELSE OUT(I,J)=IN1(I,J)/IN2(I,J) LAST=OUT(I,J) ENDIF ELSE OUT(I,J)=LAST ENDIF ELSE IF (IN2(I,J).EQ.0.0D0) THEN OUT(I,J)=0.0D0 NZERO=NZERO+1 ELSE OUT(I,J)=IN1(I,J)/IN2(I,J) ENDIF ENDIF ENDDO ENDDO C C Report number of attempts to divide by zero C IF (NZERO.GT.0) THEN IF (VERBOSE) THEN WRITE(TEXT, : '(''! '',I6,'' attempts to divide by zero.'')') NZERO CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF ENDIF 99 END SUBROUTINE DIVIDE(N1,N2,IN1,IN2,VERBOSE,OUT) C C This subroutine divides array IN1 by IN2 taking account of C the data quality, if applicable. C If the denominator is zero the output is set to zero. C The number of attempts to divide by zero is reported C IMPLICIT NONE INTEGER N1 ! 1st dimension of input array IM1 INTEGER N2 ! 2nd dimension of input array IM1 DOUBLE PRECISION IN1(N1,N2) ! First input array to be divided (numerator) DOUBLE PRECISION IN2(N1,N2) ! Second input array to be divided (denominator) DOUBLE PRECISION OUT(N1,N2) ! Output array of ratio IN1/IN2 LOGICAL VERBOSE ! If true then report all verbose statements C C Local variables C INTEGER I,J INTEGER NZERO INTEGER ISTAT CHARACTER*80 TEXT NZERO=0 DO J=1,N2 DO I=1,N1 IF (IN2(I,J).EQ.0.0D0) THEN OUT(I,J)=0.0D0 NZERO=NZERO+1 ELSE OUT(I,J)=IN1(I,J)/IN2(I,J) ENDIF ENDDO ENDDO IF (NZERO.GT.0) THEN IF (VERBOSE) THEN WRITE(TEXT, : '(''! '',I6,'' attempts to divide by zero.'')') NZERO CALL UMSPUT(TEXT,1,0,ISTAT) ENDIF ENDIF 99 END SUBROUTINE GETOF(NX,NY,PSB,PHSM,PHT,CH,D1,D2,DQ, : LNODQ,DQLIM,VERBOSE,SD,CFP,SS,HH) C C This subroutine computes the entropy (SS) and C Log Likelihood terms (HH) and fills the the partial C derivatives array SD. Data qaulity is taken into account C in computing the entropy and Log Likelihood C IMPLICIT NONE INTEGER NX ! First dimension of input arrays INTEGER NY ! Second dimension of input arrays INTEGER D1 ! 1st dimension of input array DQ INTEGER D2 ! 2nd dimension of input array DQ INTEGER DQ(D1,D2) ! Data quality array. A value < DQLIM signifies GOOD data INTEGER DQLIM ! Limiting value for GOOD data DOUBLE PRECISION PSB(NX,NY) ! Array of restored background DOUBLE PRECISION PHSM(NX,NY) ! Array of restored point sources and background DOUBLE PRECISION PHT(NX,NY) ! Data array to be restored DOUBLE PRECISION CH(NX,NY) ! Array of background convolved with smoothing kernel DOUBLE PRECISION SD(NX,NY) ! Array of the partial derivative of the correction to the background DOUBLE PRECISION CFP(NX,NY) ! Array of the convolved correction to the background DOUBLE PRECISION SS ! Value of entropy DOUBLE PRECISION HH ! Value of Log Likelihood LOGICAL LNODQ ! Flag to signify if data quality used. Set TRUE if NOT used LOGICAL VERBOSE ! If true then report all verbose statements C C Local variables C INTEGER I,J INTEGER ISTAT DOUBLE PRECISION ARG C C Initialize entropy and Log likelihood C SS=0.0D0 HH=0.0D0 C C Compute the entropy and Lok Likelihood and set the entropy C decrement array SD, taking account of the data quality if C appropriate C DO J=1,NY,1 DO I=1,NX,1 SD(I,J)=-1.0D0+CFP(I,J) IF (.NOT.LNODQ) THEN IF (DQ(I,J).LT.DQLIM) THEN IF (PSB(I,J).LE.0.0D0.OR.CH(I,J).LE.0.0D0) THEN GO TO 11 ENDIF ARG=PSB(I,J)/CH(I,J) IF (ARG.LE.0.0D0) THEN IF (VERBOSE) THEN CALL UMSPUT('! Warning: ARG negative',1,0,ISTAT) ENDIF ELSE SS=SS-PSB(I,J)*DLOG(ARG) SD(I,J)=SD(I,J)-DLOG(ARG) ENDIF ENDIF ELSE IF (PSB(I,J).LE.0.0D0.OR.CH(I,J).LE.0.0D0) THEN GO TO 11 ENDIF ARG=PSB(I,J)/CH(I,J) IF (ARG.LE.0.0D0) THEN IF (VERBOSE) THEN CALL UMSPUT('! Warning: ARG negative',1,0,ISTAT) ENDIF ELSE SS=SS - PSB(I,J)*DLOG(ARG) SD(I,J)=SD(I,J) - DLOG(ARG) ENDIF ENDIF C C Now calculate Log Likelihood C 11 IF (.NOT.LNODQ) THEN IF (DQ(I,J).LT.DQLIM) THEN IF (PHT(I,J).GT.0.0D0) THEN IF (PHSM(I,J).LE.0.0D0) THEN IF (VERBOSE) THEN CALL UMSPUT('! Warning: PHSM negative',1,0,ISTAT) ENDIF ELSE HH=HH + PHT(I,J)*DLOG(PHSM(I,J)) ENDIF ENDIF ENDIF ELSE IF (PHT(I,J).GT.0.0D0) THEN IF (PHSM(I,J).LE.0.0D0) THEN IF (VERBOSE) THEN CALL UMSPUT('! Warning: PHSM negative',1,0,ISTAT) ENDIF ELSE HH=HH + PHT(I,J)*DLOG(PHSM(I,J)) ENDIF ENDIF ENDIF ENDDO ENDDO 99 END SUBROUTINE GTDPSB(N1,N2,PSB,CF,SD,DPSB,AG,XMU) C C This subroutine calculates the unaccelerated increment C DPSB in the background image PSB C IMPLICIT NONE INTEGER N1 ! First dimension of input arrays INTEGER N2 ! Second dimension of input arrays DOUBLE PRECISION PSB(N1,N2) ! Array of restored background DOUBLE PRECISION CF(N1,N2) ! Array of convolved correction factor to background DOUBLE PRECISION SD(N1,N2) ! Array of the partial derivative of the correction to the background DOUBLE PRECISION DPSB(N1,N2) ! Array of unaccelerated increment to background array DOUBLE PRECISION AG ! Fraction contribution of the entropy term DOUBLE PRECISION XMU ! Lagrange multiplier for increment to background C C Local variables C INTEGER I,J DO J=1,N2,1 DO I=1,N1,1 DPSB(I,J)=PSB(I,J)*(CF(I,J)+AG*SD(I,J)+XMU) ENDDO ENDDO 99 END SUBROUTINE BCORRDQ(NX,NY,PSB,DPSB,XL,D1,D2,DQ,LNODQ,DQLIM) C C This subroutine applies the additive correction XL*DSPB C to the restored background PSB. If applicable data quality C is taken into account. C XL is the acceleration factor. C IMPLICIT NONE INTEGER NX ! First dimension of the input arrays INTEGER NY ! Second dimension of the input arrays INTEGER D1 ! 1st dimension of input array DQ INTEGER D2 ! 2nd dimension of input array DQ INTEGER DQ(D1,D2) ! Data quality array. A value < DQLIM signifies GOOD data INTEGER DQLIM ! Limiting value for GOOD data DOUBLE PRECISION PSB(NX,NY) ! Input/output array of the restored background estimate DOUBLE PRECISION DPSB(NX,NY) ! Input array of the additive correction to the background DOUBLE PRECISION XL ! Acceleration factor LOGICAL LNODQ ! Flag to signify if data quality used. Set TRUE if NOT used C C Local variables C INTEGER I,J DOUBLE PRECISION LAST C C Amend array PSB by the value of DPSB*XL, taking account C of the data quality if appropriate. If data quality set C then set bad pixel values of the output array to the last C good value rather than zero. C DO J=1,NY,1 DO I=1,NX,1 IF (.NOT.LNODQ) THEN IF (DQ(I,J).LT.DQLIM) THEN PSB(I,J)=PSB(I,J)+XL*DPSB(I,J) LAST=PSB(I,J) ELSE PSB(I,J)=LAST ENDIF ELSE PSB(I,J)=PSB(I,J)+XL*DPSB(I,J) ENDIF ENDDO ENDDO 99 END SUBROUTINE MONIT(MUNIT,NIT,HH,SS,QQ,SUMB,NX,NS,FLUX,NF,RMS) C C This subroutine writes out the following quantities per column C to the file open on unit UNIT. C Iteration number NIT C Value of Log likelihood HH C Value of entropy SS C Value of objective function QQ C Mean value of fit to background (SUMB/NX) C Point source flux for profile NF FLUX(NF) C Rms on fit to background and point sources RMS C IMPLICIT NONE INTEGER MUNIT INTEGER NIT INTEGER NX INTEGER NS INTEGER NF DOUBLE PRECISION RMS DOUBLE PRECISION HH DOUBLE PRECISION SS DOUBLE PRECISION QQ DOUBLE PRECISION SUMB DOUBLE PRECISION FLUX(NS) WRITE(MUNIT,20) NIT,HH,SS,QQ,SUMB/DBLE(NX),FLUX(NF),RMS 20 FORMAT(1X,I5,1X,6(1X,E14.7)) 99 END