SUBROUTINE SPEPSF C C This FORTRAN program produces a 2-D longslit spectrum of a C Point Spread Function from a number of PSF's at given C wavelengths. C A reference file is read to determine the wavelength range C of the PSF file. A list of wavelengths (must lie within the C the range of wavelengths of the reference file) is read and C then a list of files at the wavelengths. The PSF images are C assembled into a cube. The PSF is integrated over the slit C width (number of pixels in slit determined by slit width C and detector pixel size) and may be displaced from the peak C of each PSF. The set of 1-d PSF's are then assembled into an C array. The centre of the PSF can be specified to be centred C at a selected position in the output image. The PSF is C interpolated in the wavelength direction with a cubic C spline (1-d in the wavelength direction) to produce an image C with wavelength increasing along the X axis and the Y axis C corresponding to the slit length. C C Written by J. R. Walsh, ST-ECF, ESO (jwalsh@eso.org) May 2000 C IMPLICIT NONE INTEGER J INTEGER K INTEGER N1 INTEGER N2 INTEGER STAT INTEGER STAT1 INTEGER STAT2 INTEGER STAT3 INTEGER STAT4 INTEGER STATT INTEGER DIMEN(7) INTEGER DTYPE INTEGER NAXIS1 INTEGER IMDSCR1 INTEGER NAXISI INTEGER NOPSF INTEGER TBSCR INTEGER COLID INTEGER IMDSCR2 INTEGER IMDSCR3 INTEGER IMDSCR4 INTEGER ODIMEN(7) INTEGER P1 INTEGER P2 INTEGER P3 INTEGER P4 INTEGER NL INTEGER NELEM INTEGER WAV INTEGER PWAV INTEGER PWLOG INTEGER HOLD INTEGER PSFCUB INTEGER WIND INTEGER RWKPSF INTEGER RWK1 INTEGER RWK2 INTEGER RWK3 INTEGER DWKPSF INTEGER DWKFIT INTEGER DWK1 INTEGER DWK2 INTEGER DWK3 INTEGER DWK4 INTEGER DWK5 INTEGER DWK6 INTEGER DWK7 INTEGER DWK8 INTEGER DWK9 INTEGER DWK10 INTEGER DWK11 INTEGER DWK12 INTEGER DWK13 INTEGER DWK14 INTEGER DWK15 INTEGER DWK16 INTEGER DWK17 INTEGER IMSPSF INTEGER MEMI(1) REAL PIST REAL LAMST REAL DELLAM REAL SLWID REAL PXWID REAL SLIPIX REAL SUBSAM REAL CENOFF REAL PSFCEN REAL COFF REAL MEMR(1) DOUBLE PRECISION MEMD(1) CHARACTER*24 IMANAM CHARACTER*24 TABNAM CHARACTER*24 GENAME CHARACTER*24 NAME CHARACTER*24 OUTPSF CHARACTER*132 OUTEXT LOGICAL WINC LOGICAL TEXIST LOGICAL LOFLIM LOGICAL MEMB(1) COMMON/MEM/MEMD EQUIVALENCE (MEMI,MEMR,MEMD,MEMB) C C Get the reference data set name and determine the C wavelength extent C CALL UCLGST('refima',IMANAM,STAT) C C Attempt to open the data file 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(' Error. Reference image must be 2-dimensional') CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF N1=DIMEN(1) N2=DIMEN(2) ELSE WRITE(OUTEXT,13) IMANAM 13 FORMAT(' Error. Failed to open file ',A24) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF C C Attempt to read the header parameters for the wavelength C start and increment. If successfull create a wavelength array. C If no wavelengths available, use pixels numbers C 20 CALL UHDGSR(IMDSCR1,'CRPIX1',PIST,STAT2) IF (STAT2.EQ.40) THEN ! Header CRPIX1 not pesent PIST=1.0 ENDIF CALL UHDGSR(IMDSCR1,'CRVAL1',LAMST,STAT2) IF (STAT2.EQ.40) THEN ! Header CRVAL1 not pesent LAMST=1.0 ENDIF CALL UHDGSR(IMDSCR1,'CD1_1',DELLAM,STAT3) IF (STAT2.EQ.40) THEN ! Header CD1_1 not pesent DELLAM=1.0 ENDIF C C Allocate dynamic memory for the 1-D wavelength array C CALL UDMGET(N1,7,WAV,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,21) 21 FORMAT(' Error. Unable to assign memory for internal 1-D array') CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF C C Create a double precision array of the wavelengths C of the X pixels C CALL WAVCRE(N1,PIST,LAMST,DELLAM,MEMD(WAV)) C C Open the table for the list of wavelengths of the PSF's C to be input C 30 CALL UCLGST('psfwavs',TABNAM,STAT) CALL UTTACC(TABNAM,TEXIST,STAT) IF (.NOT.TEXIST) THEN WRITE(OUTEXT,31) TABNAM 31 FORMAT(' Error. Table file does not exist ',A24) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UTTOPN(TABNAM,1,TBSCR,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,32) TABNAM 32 FORMAT(' Error. Table could not be opened ',A24) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF C C Get the number of rows in the table C CALL UTPGTI(TBSCR,21,NOPSF,STAT1) IF (STAT1.NE.0) THEN WRITE(OUTEXT,33) TABNAM 33 FORMAT(' Error. Cannot get no. of rows in table ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) ENDIF C C Allocate dynamic memory for the 1-D arrays of the C wavelength of the PSF's, the array of flags for the C values and the integer channel numbers corresponding C to these wavelengths C CALL UDMGET(NOPSF,6,PWAV,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,21) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(NOPSF,1,PWLOG,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,34) 34 FORMAT(' Error. Unable to assign memory for internal 1-D', : ' Boolean array') CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(NOPSF,6,WIND,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,35) 35 FORMAT(' Error. Unable to assign memory for internal 1-D', : ' real array') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Get the column ID's and read the wavelengths to the C array PWAV and close the table file C CALL UTCNUM(TBSCR,1,COLID,STAT) CALL UTCGTR(TBSCR,COLID,1,NOPSF,MEMR(PWAV),MEMB(PWLOG),STAT) CALL UTTCLO(TBSCR,STAT) C C Get the number of non-zero values in the array PWAV C CALL PWAVCHK(NOPSF,MEMR(PWAV),N1,MEMD(WAV),WINC) IF (.NOT.WINC) THEN WRITE(OUTEXT,37) 37 FORMAT(' Error. All wavelengths not in range of reference', : ' image') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Determine real channel numbers corresponding to wavelengths of C PSF's C CALL WAVIND(N1,MEMD(WAV),NOPSF,MEMR(PWAV),MEMR(WIND),LOFLIM) IF (.NOT.LOFLIM) THEN WRITE(OUTEXT,38) 38 FORMAT(' Error. Wavelengths cannot be matched to pixel', : ' numbers') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF C C Initialize the filename template processing for the C list of PSF images C 100 CALL UCLGST('psfdata',GENAME,STAT) CALL TIMOTP(GENAME,IMDSCR2,STAT) C C Get first data set name and read the values into the C first plane of the cube PSFHOLD C 110 CALL TIMXTP(IMDSCR2,NAME,STATT) IF (STATT.NE.0) THEN WRITE(OUTEXT,111) NAME 111 FORMAT(' Error opening data file ',A24) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF IF (STATT.EQ.0) THEN CALL UIMOPN(NAME,1,IMDSCR3,STAT1) CALL UIMGID(IMDSCR3,DTYPE,NAXISI,DIMEN,STAT2) IF (NAXIS1.NE.2) THEN WRITE(OUTEXT,11) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF P1=DIMEN(1) P2=DIMEN(2) C C Allocate dynamic memory for the 2-D input images C of dimensions P1 by P2 (real number) C NELEM=P1*P2 CALL UDMGET(NELEM,6,HOLD,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,113) 113 FORMAT(' Error. Unable to assign memory for internal', : ' 2-D real array') CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF C C Allocate dynamic memory for the cube of NOPSF input C images of dimensions P1 by P2 (real numbers) C NELEM=NOPSF*P1*P2 CALL UDMGET(NELEM,6,PSFCUB,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,115) 115 FORMAT(' Error. Unable to assign memory for internal', : ' 3-D real array') CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF IF (STAT2.EQ.0) THEN CALL UIGS2R(IMDSCR3,1,P1,1,P2,MEMR(HOLD),STAT3) IF (STAT3.NE.0) THEN WRITE(OUTEXT,117) NAME 117 FORMAT(' Error reading data from file ',A24) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF ELSE WRITE(OUTEXT,119) NAME 119 FORMAT(' Failed to read ',A24) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF ELSE C C Input list empty C WRITE(OUTEXT,121) 121 FORMAT(' Error. No data files read') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UIMCLO(IMDSCR3,STAT) C C Copy the array into the first plane of cube C CALL INTOCUB(P1,P2,MEMR(HOLD),P1,P2,NOPSF,MEMR(PSFCUB),1) C C Loop through the list of data sets reading the values for C each into array HOLD and copying each into a plane C of the 3-D array PSFCUB C 200 J=2 DO K=2,NOPSF,1 CALL TIMXTP(IMDSCR2,NAME,STATT) IF (STATT.EQ.10) THEN WRITE(OUTEXT,211) NAME 211 FORMAT(' Error opening data file ',A24) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF IF (STATT.EQ.-2) THEN GO TO 290 ENDIF IF (STATT.EQ.0) THEN CALL UIMOPN(NAME,1,IMDSCR3,STAT1) CALL UIMGID(IMDSCR3,DTYPE,NAXISI,DIMEN,STAT2) IF (NAXIS1.NE.2) THEN WRITE(OUTEXT,11) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF P3=DIMEN(1) P4=DIMEN(2) IF (P3.NE.P1.OR.P4.NE.P2) THEN WRITE(OUTEXT,212) NAME 212 FORMAT(' File ',A24,' not same size') CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF IF (STAT2.EQ.0) THEN CALL UIGS2R(IMDSCR3,1,P3,1,P4,MEMR(HOLD),STAT3) IF (STAT3.NE.0) THEN WRITE(OUTEXT,117) NAME CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UIMCLO(IMDSCR3,STAT4) CALL INTOCUB(P3,P4,MEMR(HOLD),P1,P2,NOPSF,MEMR(PSFCUB),J) J=J+1 ELSE WRITE(OUTEXT,213) NAME 213 FORMAT(' Failed to read ',A24) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF ENDIF ENDDO 290 IF (J-1.NE.NOPSF) THEN WRITE(OUTEXT,291) 291 FORMAT(' Incompatible number of wavelengths and input PSFs') CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF C C Close down the filename template C 295 CALL TIMCTP(IMDSCR2,STAT) C C Get the width of the slit in arcsec and the pixel C width in arcsec C CALL UCLGSR('sliwid',SLWID,STAT) CALL UCLGSR('pixwid',PXWID,STAT) C C Get the ratio of the reference frame pixels to the C PSF pixels C CALL UCLGSR('subsam',SUBSAM,STAT) IF (PXWID.NE.0.0) THEN SLIPIX=SLWID/PXWID ELSE SLIPIX=SLWID ENDIF C C Get the position of the slit center relative to the C peak of the PSF (offset in the dispersion direction) C in PSF pixels C CALL UCLGSR('cenoff',CENOFF,STAT) IF (PXWID.NE.0.0) THEN COFF=CENOFF/PXWID ELSE COFF=CENOFF ENDIF C C Get the pixel value at which to centre the PSF C (perpendicular to the wavelength direction) C CALL UCLGSR('psfcen',PSFCEN,STAT) C C Allocate dynamic memory for the 2-D and 1-D real and C double precision internal work arrays required for C the interpolation and fitting C NL=P2*NINT(SUBSAM) NELEM=P1*P2 CALL UDMGET(NELEM,6,RWKPSF,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,311) 311 FORMAT(' Unable to assign memory for internal 2-D'// : ' real array') CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF C C Now double precision arrays C CALL UDMGET(NOPSF,7,RWK1,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) 315 FORMAT(' Unable to assign memory for internal 1-D'// : ' double precision array') CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(P1,7,RWK2,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(P2,7,RWK3,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(NL,7,DWK1,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(NL,7,DWK2,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(NL,7,DWK3,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(P1,7,DWK4,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF NELEM=NOPSF*NL CALL UDMGET(NELEM,7,DWKPSF,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,317) 317 FORMAT(' Unable to assign memory for internal 2-D'// : ' double precision array') CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF NELEM=NOPSF*NL CALL UDMGET(NELEM,7,DWKFIT,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,317) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF NELEM=P1*P2 CALL UDMGET(NELEM,7,DWK5,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,317) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(NELEM,7,DWK6,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,317) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(P2,7,DWK7,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(P2,7,DWK8,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(P2,7,DWK9,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(P1,7,DWK10,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(NOPSF,7,DWK11,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(NOPSF,7,DWK12,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(NOPSF,7,DWK13,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(NL,7,DWK14,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(NL,7,DWK15,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(NOPSF,7,DWK16,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF CALL UDMGET(NL,7,DWK17,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,315) CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF C C Allocate dynamic memory for the 2-D output PSF C spectrum image of dimensions N1 by N2 (real number) C NELEM=N1*N2 CALL UDMGET(NELEM,6,IMSPSF,STAT) IF (STAT.NE.0) THEN WRITE(OUTEXT,321) 321 FORMAT(' Unable to assign memory for output 2-D array') CALL UMSPUT(OUTEXT,3,0,STAT) GO TO 990 ENDIF C C Call the subroutine to produce each PSF integrated over the C slit width and then interpolate the PSF's onto the C wavelength grid of the reference image C CALL PSFPRO(P1,P2,NOPSF,MEMR(PSFCUB),MEMR(WIND),SLIPIX,SUBSAM, : COFF,PSFCEN,MEMR(RWKPSF),MEMD(DWKPSF), : MEMD(DWKFIT),MEMD(RWK1),MEMD(RWK2),MEMD(RWK3),NL, : MEMD(DWK1),MEMD(DWK2),MEMD(DWK3),MEMD(DWK4), : MEMD(DWK5),MEMD(DWK6),MEMD(DWK7), : MEMD(DWK8),MEMD(DWK9),MEMD(DWK10),MEMD(DWK11), : MEMD(DWK12),MEMD(DWK13),MEMD(DWK14),MEMD(DWK15), : MEMD(DWK16),MEMD(DWK17),N1,N2,MEMR(IMSPSF)) C C C Free the dynamic memory allocated which is no longer C required C CALL UDMFRE(WIND,6,STAT) CALL UDMFRE(PSFCUB,6,STAT) CALL UDMFRE(RWKPSF,6,STAT) CALL UDMFRE(RWK1,7,STAT) CALL UDMFRE(RWK2,7,STAT) CALL UDMFRE(RWK3,7,STAT) CALL UDMFRE(DWKPSF,7,STAT) CALL UDMFRE(DWKFIT,7,STAT) CALL UDMFRE(DWK1,7,STAT) CALL UDMFRE(DWK2,7,STAT) CALL UDMFRE(DWK3,7,STAT) CALL UDMFRE(DWK4,7,STAT) CALL UDMFRE(DWK5,7,STAT) CALL UDMFRE(DWK6,7,STAT) CALL UDMFRE(DWK7,7,STAT) CALL UDMFRE(DWK8,7,STAT) CALL UDMFRE(DWK9,7,STAT) CALL UDMFRE(DWK10,7,STAT) CALL UDMFRE(DWK11,7,STAT) CALL UDMFRE(DWK12,7,STAT) CALL UDMFRE(DWK13,7,STAT) CALL UDMFRE(DWK14,7,STAT) CALL UDMFRE(DWK15,7,STAT) CALL UDMFRE(DWK16,7,STAT) CALL UDMFRE(DWK17,7,STAT) C C Open an output file for the PSF spectrum image C 400 CALL UCLGST('outpsf',OUTPSF,STAT) ODIMEN(1)=N1 ODIMEN(2)=N2 CALL UIMCRE(OUTPSF,6,2,ODIMEN,IMDSCR4,STAT1) IF (STAT1.EQ.0) THEN CALL UHDCPY(IMDSCR1,IMDSCR4,STAT2) CALL UIPS2R(IMDSCR4,1,N1,1,N2,MEMR(IMSPSF),STAT3) IF (STAT3.NE.0) THEN WRITE(OUTEXT,411) OUTPSF 411 FORMAT(' Error writing data to file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UIMCLO(IMDSCR4,STAT4) ELSE WRITE(OUTEXT,412) OUTPSF 412 FORMAT(' Failed to open data file ',A24) CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 990 ENDIF CALL UDMFRE(IMSPSF,6,STAT) CALL UIMCLO(IMDSCR1,STAT1) GO TO 999 990 CALL UMSPUT(' Error occurred. No output',1,0,STAT) 999 END SUBROUTINE PWAVCHK(NPSF,PWAV,N,WAV,WIN) C C This subroutine does some simple checking of the C correspondence of the wavelengths of the PSF (PWAV) C against the total extent of the reference spectrum C wavelengths (WAV). If WIN is returned false there C is a mismatch C IMPLICIT NONE INTEGER NPSF ! Dimension of input array PWAV INTEGER N ! Dimension of input array WAV REAL PWAV(NPSF) ! Array of wavelengths of the PSFs DOUBLE PRECISION WAV(N) ! Array of the wavelengths LOGICAL WIN ! Flagged false if PWAV does not fall within WAV C C Local variables C INTEGER I,J J=1 DO I=1,NPSF,1 IF (PWAV(I).NE.0.0) THEN PWAV(J)=PWAV(I) J=J+1 ENDIF ENDDO NPSF=J-1 C C Check that the wavelength values in PWAV are in the range C of wavelengths of WAV C DO I=1,NPSF,1 IF (PWAV(I).GT.REAL(WAV(1)).AND.PWAV(I).LT.REAL(WAV(N))) THEN WIN=.TRUE. ELSE WIN=.FALSE. ENDIF ENDDO 99 END SUBROUTINE WAVIND(N,WAV,NW,PWAV,WIND,LOFLIM) C C This subroutine returns a real array, WIND, of C the channel numbers in the array WAV of the C wavelengths in the array PWAV. If for some C reason the pixel numbers cannot be found over C the wavelength range then LOFLIM is set FALSE C IMPLICIT NONE INTEGER N ! Dimension of input array WAV INTEGER NW ! Dimension of input array PWAV REAL PWAV(NW) ! Array of wavelengths of PSFs DOUBLE PRECISION WAV(N) ! Array of wavelengths REAL WIND(NW) ! Output array of positions of PWAV values in WAV LOGICAL LOFLIM ! Logical flag to indicate failure to set pixel numbers C C Local variables C INTEGER I,J REAL LAM LOGICAL LOFF C C For all the wavelengths of the array PWAV, determine the C corresponding pixel number in the array WAV C LOFLIM=.TRUE. DO I=1,NW,1 LAM=PWAV(I) LOFF=.FALSE. DO J=1,N-1,1 IF (LAM.GE.REAL(WAV(J)).AND.LAM.LT.REAL(WAV(J+1))) THEN WIND(I)=REAL(J) + ( LAM-REAL(WAV(J)) )/ :( REAL(WAV(J+1)) - REAL(WAV(J)) ) LOFF=.TRUE. GO TO 20 ENDIF ENDDO 20 CONTINUE ENDDO IF (.NOT.LOFF) THEN LOFLIM=.FALSE. ENDIF 99 END SUBROUTINE INTOCUB(N1,N2,HOLD,N3,N4,NP,CUB,NS) C C This subroutine copies the 2-D array HOLD into C plane NS of the 3-d array CUB C IMPLICIT NONE INTEGER N1 ! 1st dimension of input array HOLD INTEGER N2 ! 2nd dimension of input array HOLD INTEGER N3 ! 1st dimension of updated array CUB INTEGER N4 ! 2nd dimension of updated array CUB INTEGER NP ! 3rd dimension of updated array CUB INTEGER NS ! Plane of CUB in which to copy array HOLD REAL HOLD(N1,N2) ! Input array to copy to cube REAL CUB(N3,N4,NP) ! Cube to update with array HOLD in plane NS C C Local variables C INTEGER I,J DO I=1,N1,1 DO J=1,N2,1 CUB(I,J,NS)=HOLD(I,J) ENDDO ENDDO 99 END SUBROUTINE PSFPRO(N1,N2,NW,CUB,RWIN,WID,SUBSAM,XOFF,CEN, : PSFA,PSFA1D,PSFIT,DWK1,X,Y,NL,YL, : YU,YC,DWK2,DWK3,DWK4,DWK5,DWK6, : DWK7,DWK8,DWK9,DWK10,DWK11,DWK12, : DWK13,DWK14,DWK15,N3,N4,OUTIMA) C C This subroutine takes NW PSF's from the planes of the 3-D C array CUB, integrates the signal in the X-direction across C the slit width and integrates them in the Y direction onto C the pixels, which may be subsampled a factor SUBSAM. The C centre of the region of integration is offset from the C peak of the PSF by XOFF (pixels) and the Y centre of the C PSF is shifted to CEN. The 1-d PSF's at their respective C wavelengths are then used to produce a PSF spectrum by C interpolating in the X-direction over the extent of the C output image. C IMPLICIT NONE INTEGER N1 ! 1st dimension of input array CUB INTEGER N2 ! 2nd dimension of input array CUB INTEGER NW ! 3rd dimension of input array CUB INTEGER N3 ! 1st dimension of output array OUTIMA INTEGER N4 ! 2nd dimension of output array OUTIMA INTEGER NL ! Dimension of input arrays YL and YU REAL CUB(N1,N2,NW) ! Input cube of NW PSF arrays REAL RWIN(NW) ! Array of the X pixel (wavelength) positions of the PSF's REAL WID ! Pixel width in X direction REAL SUBSAM ! Subsampling factor data pixels/PSF pixels REAL XOFF ! X offset of PSF centre from slit centre REAL CEN ! Required centre in Y direction for PSF spectrum REAL PSFA(N1,N2) ! Holding array for 2-d PSF array REAL OUTIMA(N3,N4) ! Output array of PSF spectrum DOUBLE PRECISION X(N1) ! Array for real X pixel values DOUBLE PRECISION Y(N2) ! Array for real Y pixel values DOUBLE PRECISION PSFA1D(NW,NL) ! Holding array for NW 1-d fitted PSF spectra DOUBLE PRECISION PSFIT(NW,NL) ! Holding array for NW 1-d PSF spectra DOUBLE PRECISION YL(NL) ! Array of lower pixel limits of PSF (Y direction) DOUBLE PRECISION YU(NL) ! Array of upper pixel limits of PSF (Y direction) DOUBLE PRECISION YC(NL) ! Array of centres of pixels of PSF (Y direction) DOUBLE PRECISION DWK1(NW) ! Holding array for centres of PSFs DOUBLE PRECISION DWK2(N1) ! Work array DOUBLE PRECISION DWK3(N1,N2) ! Work array DOUBLE PRECISION DWK4(N1,N2) ! Work array DOUBLE PRECISION DWK5(N2) ! Work array DOUBLE PRECISION DWK6(N2) ! Work array DOUBLE PRECISION DWK7(N2) ! Work array DOUBLE PRECISION DWK8(N1) ! Work array DOUBLE PRECISION DWK9(NW) ! Work array DOUBLE PRECISION DWK10(NW) ! Work array DOUBLE PRECISION DWK11(NW) ! Work array DOUBLE PRECISION DWK12(NL) ! Work array DOUBLE PRECISION DWK13(NL) ! Work array DOUBLE PRECISION DWK14(NW) ! Work array DOUBLE PRECISION DWK15(NL) ! Work array C C Local variables C INTEGER I INTEGER NL INTEGER IFAIL INTEGER STAT DOUBLE PRECISION X1 DOUBLE PRECISION X2 CHARACTER*72 OUTEXT C C Form real arrays of the X and Y array indexes C Use the IRAF convention: pixel 1 extends from X=0.5 to 1.5 C DO I=1,N1,1 X(I)=DBLE(I) ENDDO DO I=1,N2,1 Y(I)=DBLE(I) ENDDO C C For each PSF in each plane of CUB, find the peak and extract C the 1-d PSF integrated over the slit of width WID centred at C XOFF relative to PSF peak C DO I=1,NW,1 C C Copy the Ith plane of CUB to a 2-d array C WRITE(OUTEXT,380) I 380 FORMAT(' Working on PSF #',I4) CALL UMSPUT(OUTEXT,3,0,STAT) CALL OUTOCUB(N1,N2,NW,CUB,I,PSFA) C C Find the peak and the pixel range of the slit width and C the arrays of values of the y positions of the reference C pixels (lower edge, centre and upper edge) C CALL SLIEXT(N1,N2,PSFA,X,Y,WID,XOFF,SUBSAM,DWK2,X1,X2,NL, : YL,YU,YC,IFAIL) IF (IFAIL.NE.0) THEN WRITE(OUTEXT,19) I 19 FORMAT(1X,I4,'th PSF image empty. Output empty.') CALL UMSPUT(OUTEXT,1,0,STAT) GO TO 99 ENDIF C C Integrate the area of the PSF over the slit and the C reference pixels and place the 1-d array in the Ith C section of the 2-d array PSFA1D C CALL PSF1D(N1,N2,X,Y,PSFA,X1,X2,NL,YL,YU, : DWK3,DWK4,DWK5,DWK6,DWK7,DWK8,NW,I,PSFA1D) ENDDO C C Interpolate the NW 1-d PSF's in the wavelength direction C onto the 2-D output array C CALL PSF1D2(NW,NL,N1,N2,X,Y,PSFA1D,PSFIT,RWIN,CEN, : SUBSAM,YL,YU,YC,DWK9,DWK1,DWK10,DWK11,DWK12, : DWK13,DWK14,DWK15,N3,N4,OUTIMA) 99 END SUBROUTINE OUTOCUB(N1,N2,NP,CUB,NS,HOLD) C C This subroutine copies the NS'th plane of C the cube CUB to the 2-D array HOLD C IMPLICIT NONE INTEGER N1 ! 1st dimension of input array CUB INTEGER N2 ! 2nd dimension of input array CUB INTEGER NP ! 3rd dimension of input array CUB INTEGER NS ! Plane of CUB to copy to 2-d array REAL CUB(N1,N2,NP) ! Input 3-D array REAL HOLD(N1,N2) ! Output 2-d array sliced from CUB C C Local variables C INTEGER I,J DO I=1,N1,1 DO J=1,N2,1 HOLD(I,J)=CUB(I,J,NS) ENDDO ENDDO 99 END SUBROUTINE SLIEXT(N1,N2,INAR,X,Y,WID,XOFF,SUBSAM,XARR,X1,X2, : NL,YL,YU,YC,IFAIL) C C This subroutine determines the X-limits of the slit C width in pixels and the Y-limits for the reference pixels C in terms of the original pixels. If no peak pixel C is found in array INAR, the subroutine returns with C IFAIL set to 1 C IMPLICIT NONE INTEGER N1 ! 1st dimension of input array INAR INTEGER N2 ! 2nd dimension of input array INAR INTEGER NL ! Input dimensions of array YL and YU INTEGER IFAIL ! Failure indicator if peak not found in PSF REAL INAR(N1,N2) ! Input 2d array REAL WID ! Width of pixel in wavelength REAL XOFF ! Offset of PSF centre from slit centre REAL SUBSAM ! Subsampling factor data/PSF DOUBLE PRECISION X(N1) ! Array of X pixel values for INAR DOUBLE PRECISION Y(N2) ! Array of Y pixel values for INAR DOUBLE PRECISION XARR(N1) ! Holding array for Y direction collapsed image DOUBLE PRECISION X1 ! X position of centroid of INAR DOUBLE PRECISION X2 ! Y position of centroid of INAR DOUBLE PRECISION YL(NL) ! Output array of lower Y limits of reference pixels in the Y in terms of the original pixels DOUBLE PRECISION YU(NL) ! Output array of upper Y limits of reference pixels in the Y in terms of the original pixels DOUBLE PRECISION YC(NL) ! Output array of centres of Y sumsampled reference pixels in t C C Local variables C INTEGER I,J,JJ INTEGER IX INTEGER ISUBSAM c INTEGER STAT DOUBLE PRECISION SUM DOUBLE PRECISION MAXI DOUBLE PRECISION XCEN DOUBLE PRECISION XFWHM DOUBLE PRECISION YY DOUBLE PRECISION DYSAM LOGICAL LFIND c CHARACTER*132 OUTEXT C C Collapse the 2-D array in the Y direction C DO I=1,N1,1 SUM=0.0D0 DO J=1,N2,1 SUM=SUM+DBLE(INAR(I,J)) ENDDO XARR(I)=SUM ENDDO C C Find the peak pixel in the array XARR C LFIND=.FALSE. MAXI=-1.0D30 DO I=1,N1,1 IF (XARR(I).GT.MAXI) THEN MAXI=INAR(I,J) IX=I LFIND=.TRUE. ENDIF ENDDO IF (.NOT.LFIND) THEN IFAIL=1 GO TO 99 ENDIF C C Centroid the peak in the array XARR (i.e. X direction) C CALL CENTCEN(N1,X,XARR,DBLE(IX),XCEN,XFWHM) C C Determine X limits of slit allowing for the offset of the slit C from the centre of the source C X1=XCEN + DBLE(XOFF - 0.50*WID) X2=XCEN + DBLE(XOFF + 0.50*WID) c WRITE(OUTEXT,411) X1,XCEN,X2 c411 FORMAT(' SLIT EXTENT: LEFT, CENTRE, RIGHT ',3(F12.4)) c CALL UMSPUT(OUTEXT,3,0,STAT) C C Form the arrays of the lower and upper limits of the reference C pixels in the Y direction in terms of the original pixels C JJ=1 ISUBSAM=INT(SUBSAM) DYSAM=DBLE(SUBSAM) DO I=1,N2,1 YY=DBLE(I)-0.50D0 DO J=1,ISUBSAM,1 YL(JJ)=YY+(DBLE(J-1)/DYSAM) YU(JJ)=YY+(DBLE(J)/DYSAM) YC(JJ)=YY+(DBLE(J-1)/DYSAM) + (0.50D0/DYSAM) JJ=JJ+1 ENDDO ENDDO 99 END SUBROUTINE PSF1D(N1,N2,X,Y,INAR,X1,X2,NL,YL,YU, : DARR,DARRD,WK1,WK2,WK3,WK4,NW,NS,OUT1D) C C This subroutine forms the 1-d integrated profile through the C 2-d array INAR. The integration is perfomed from X1 to X2 C in the X direction and over the array of limits specified C by the arrays YL and YU. The resulting 1-D array is copied C to the NS'th row of the array OUT1D C IMPLICIT NONE INTEGER N1 ! 1st dimension of array INAR INTEGER N2 ! 2nd dimension of array INAR INTEGER NL ! Number of elements in YL and YU INTEGER NW ! No of PSF's to integrate (1st dimension of ouput array) INTEGER NS ! Indicator for which section of OUT1D to copy to REAL INAR(N1,N2) ! Input 2D array of PSF DOUBLE PRECISION X(N1) ! Array of real X pixel values DOUBLE PRECISION Y(N2) ! Array of real Y pixel values DOUBLE PRECISION X1 ! Lower X limit for integration DOUBLE PRECISION X2 ! Upper X limit for integration DOUBLE PRECISION YL(NL) ! Output array of lower Y limits of reference pixels in the Y in terms of the original pixels DOUBLE PRECISION YU(NL) ! Output array of upper Y limits of reference pixels in the Y in terms of the original pixels DOUBLE PRECISION DARR(N1,N2) ! Holding array for double precision copy of INAR DOUBLE PRECISION DARRD(N1,N2) ! Holding array for array of second derivatives of the interpolating function DOUBLE PRECISION WK1(N2) ! Work array for interpolation DOUBLE PRECISION WK2(N2) ! Work array for interpolation DOUBLE PRECISION WK3(N2) ! Work array for interpolation DOUBLE PRECISION WK4(N1) ! Work array for interpolation DOUBLE PRECISION OUT1D(NW,NL) ! Output array of Ns'th 1-d PSF C C Local variables C INTEGER I,J c INTEGER STAT DOUBLE PRECISION DY1 DOUBLE PRECISION DY2 DOUBLE PRECISION Z c CHARACTER*72 OUTEXT LOGICAL LZERO COMMON/LOW/DY1 COMMON/UPP/DY2 C C Copy input values to double precision C DO I=1,N1,1 DO J=1,N2,1 DARR(I,J)=DBLE(INAR(I,J)) ENDDO ENDDO C C Initialize output array C DO I=1,NL,1 OUT1D(NS,I)=0.0D0 ENDDO C C Form the array of second derivatives of the interpolating C function at the tabulated points for the bicubic spline fit C CALL SPLIE2(N1,N2,Y,DARR,WK1,WK2,WK3,DARRD) C C Integrate the flux across the slit (range X1 to X2) and DY1 to DY2 DO I=1,NL,1 DY1=YL(I) DY2=YU(I) C C Check if not all zeros over region of integration C CALL CHECK0(N1,N2,X1,X2,DY1,DY2,DARR,LZERO) IF (LZERO) THEN CALL QUAD2D(X1,X2,Z,N1,N2,X,Y,DARR,DARRD, : WK1,WK2,WK3,WK4) IF (Z.GT.1.0D-30) THEN OUT1D(NS,I)=Z c WRITE(OUTEXT,135) I,Z c135 FORMAT(' Flux integrated over slit: Ypix, Val ', c :I5,1X,E12.5) c CALL UMSPUT(OUTEXT,3,0,STAT) ENDIF ENDIF ENDDO 99 END SUBROUTINE CHECK0(N1,N2,X1,X2,Y1,Y2,INAR,LZERO) C C This subroutine checks that there are non-zero values C of the array INAR over the x range X1 to X2 and the C y range Y1 to Y2. LZERO is returned true if non-zero C values are found C IMPLICIT NONE INTEGER N1 ! 1st dimension of input array INAR INTEGER N2 ! 2nd dimension of input array INAR DOUBLE PRECISION X1 ! Lower X limit for integration DOUBLE PRECISION X2 ! Upper X limit for integration DOUBLE PRECISION Y1 ! Lower Y limit for integration DOUBLE PRECISION Y2 ! Upper Y limit for integration DOUBLE PRECISION INAR(N1,N2) ! Input 2-d array to be checked LOGICAL LZERO ! Output flag - false if ony zeros found in rgion of INAR C C Local variables C INTEGER I,J INTEGER IX1 INTEGER IX2 INTEGER IY1 INTEGER IY2 IX1=INT(X1) IF (IX1.LT.1) THEN IX1=1 ENDIF IF (IX1.GT.N1) THEN IX1=N1 ENDIF IX2=INT(X2) + 1 IF (IX2.LT.1) THEN IX2=1 ENDIF IF (IX2.GT.N1) THEN IX2=N1 ENDIF IY1=INT(Y1) IF (IY1.LT.1) THEN IY1=1 ENDIF IF (IY1.GT.N2) THEN IY1=N2 ENDIF IY2=NINT(Y2) + 1 IF (IY2.LT.1) THEN IY2=1 ENDIF IF (IY2.GT.N2) THEN IY2=N2 ENDIF LZERO=.FALSE. DO J=IY1,IY2,1 DO I=IX1,IX2,1 IF (INAR(I,J).GT.1.0D-30) THEN LZERO=.TRUE. GO TO 99 ENDIF ENDDO ENDDO 99 END SUBROUTINE PSF1D2(NW,NL,N1,N2,X,Y,PSFA1D,SPLFIT,WIND,CEN, : SUBSAM,YL,YU,YC,XIN,PSMAX,DXARR,DXARRD,DYARR, : DYARRD,DUX,DUY,N3,N4,OUTARR) C C This subroutine takes the NW 1-D profiles of the C slit integrated PSF from the array PSFA1D, interpolates C each one using cubic spline interpolation onto the C subsampled grid point and shifts the profiles so that on C the output grid the peak is at position CEN. C Each 1-D shifted and subsampled array at its X C position in the array WIND is then interpolated and C extrapolated, again using cubic splines, in the C X-direction of the output array OUTARR. C IMPLICIT NONE INTEGER NW ! 1st dimension of input array PSFA1D INTEGER NL ! 2nd dimension of input array PSFA1D INTEGER N1 ! 1st dimension of input array X INTEGER N2 ! 1st dimension of input array Y INTEGER N3 ! 1st dimension of output array OUTARR INTEGER N4 ! 2nd dimension of output array OUTARR REAL WIND(NW) ! Real X pixel values of wavelengths of PSF's REAL CEN ! Required centre of slit in X direction REAL SUBSAM ! Pixel subsampling factor REAL OUTARR(N3,N4) ! Output 2-D spectrum of point source DOUBLE PRECISION X(N1) ! Array of real X pixel values of PSF images DOUBLE PRECISION Y(N2) ! Array of real Y pixel values of PSF images DOUBLE PRECISION PSFA1D(NW,NL) ! Input array of 1-d PSF's for interpolation in wavelength DOUBLE PRECISION SPLFIT(NW,NL) ! Array of 1-d interpolated PSFs DOUBLE PRECISION YL(NL) ! Array of lower limits of Y pixels DOUBLE PRECISION YU(NL) ! Array of upper limits of Y pixels DOUBLE PRECISION YC(NL) ! Array of centres of Y pixels DOUBLE PRECISION XIN(NW) ! Array of X pixel positions of PSF's DOUBLE PRECISION PSMAX(NW) ! Array of positions of centroid of PSF DOUBLE PRECISION DXARR(NW) ! Array of values to interpolate in X direction DOUBLE PRECISION DXARRD(NW) ! Array of second derivatives of the interpolating function DOUBLE PRECISION DYARR(NL) ! Array of values of 1-D PSF for interpolation DOUBLE PRECISION DYARRD(NL) ! Array of second derivatives of the interpolating function DOUBLE PRECISION DUX(NW) ! Holding array for SPLINE (X direction interpolation) DOUBLE PRECISION DUY(NL) ! Hilding array for SPLINE (Y direction interpolation) C C Local variables C INTEGER I,J,K INTEGER IPM INTEGER ISHIFT REAL ZSUM DOUBLE PRECISION PMAX DOUBLE PRECISION PCEN DOUBLE PRECISION PFWHM DOUBLE PRECISION Z DOUBLE PRECISION YY DOUBLE PRECISION XX LOGICAL LVAL c integer stat c character*80 outext C C Initialize the output array C DO I=1,N3,1 DO J=1,N4,1 OUTARR(I,J)=0.0 ENDDO ENDDO C C Set the pixel centres of the subsampled pixels to be C a running number C DO I=1,NL,1 YC(I)=DBLE(I) ENDDO C C Find the peaks of the PSF's in the input array, C by centroiding the non-zero values C DO K=1,NW,1 IPM=NL/2 ! Safe default value PMAX=-1.0D30 DO J=1,NL,1 IF (PSFA1D(K,J).GT.1.0D-30.AND.PSFA1D(K,J).LT.1.0D30) THEN DYARR(J)=PSFA1D(K,J) IF (DYARR(J).GT.PMAX) THEN PMAX=DYARR(J) IPM=J ENDIF ELSE DYARR(J)=0.0D0 ENDIF ENDDO CALL CENTCEN(NL,YC,DYARR,DBLE(IPM),PCEN,PFWHM) PSMAX(K)=PCEN c write(outext,301) K,psmax(K) c301 format(' PSF1D2: K,PSFMAX ',I5,1x,F12.5) c call umsput(outext,1,0,stat) C C For each input PSF fit with a cubic spline, shifting the C Y values so that the peak is at channel CEN. Interpolate C the values to those of the original spectrum image (N3xN4) C C Form the array of second derivatives of the interpolating C function at the tabulated points for the bicubic spline fit C CALL SPLINE(NL,YC,DYARR,2.D30,2.D30,DUY,DYARRD) C C Evalaute the cubic spline on the output Y values C ISHIFT=INT(DBLE(CEN) - PSMAX(1)) c write(outext,302) ISHIFT c302 format(' PSF1D2: ISHIFT ',I5) c call umsput(outext,1,0,stat) DO J=1,NL,1 c YY=YC(J) - (DBLE(CEN) - PSMAX(K)*DBLE(SUBSAM) - c : DBLE(ISHIFT)) YY=YC(J) - (DBLE(CEN) - PSMAX(K) - DBLE(ISHIFT)) IF (YY.GE.YC(1).AND.YY.LE.YC(NL)) THEN CALL SPLINT(NL,YC,DYARR,DYARRD,YY,Z) SPLFIT(K,J)=Z ELSE SPLFIT(K,J)=0.0D0 ENDIF ENDDO ENDDO C C Interpolate in the X-direction across the NW PSF's onto C the X-values of the output array C DO K=1,NW,1 XIN(K)=DBLE(WIND(K)) ENDDO DO J=1,NL,1 C C Check if all values for interpolating array are greater C than some low value (1E-12 of the peak). If not then C don't attempt to interpolate C LVAL=.FALSE. DO K=1,NW,1 IF (DABS(SPLFIT(K,J)).GT.1.0D-12) THEN LVAL=.TRUE. ENDIF ENDDO IF (LVAL) THEN C C Copy the Ith row of SPLFIT to a 1-D double C precision array for interpolation C DO K=1,NW,1 DXARR(K)=SPLFIT(K,J) ENDDO C C Form the array of second derivatives of the C interpolating function at the tabulated points C for the bicubic spline fit C CALL SPLINE(NW,XIN,DXARR,2.D30,2.D30,DUX,DXARRD) C C Evaluate the cubic spline on the original X values C DO I=1,N3,1 XX=DBLE(I) CALL SPLINT(NW,XIN,DXARR,DXARRD,XX,Z) IF ((J+ISHIFT).GE.1.AND.(J+ISHIFT).LE.N4) THEN IF (Z.GT.0.0D0) THEN OUTARR(I,J+ISHIFT)=REAL(Z) ENDIF ENDIF ENDDO ENDIF ENDDO C C Normalise all the X-sections of the output spectrum C to SUM of 1.0 C DO I=1,N3,1 ZSUM=0.0 DO J=1,N4,1 ZSUM=ZSUM+OUTARR(I,J) ENDDO DO J=1,N4,1 OUTARR(I,J)=OUTARR(I,J)/ZSUM ENDDO ENDDO 99 END SUBROUTINE QUAD2D(X1,X2,SS,M,N,XA,YA,XYA,XY2D, : YTMP,UTMP,Y2TMP,YYTMP) C C This is the Numerical Recipes 3-D integration routine C (using Gaussian quadrature) adapted for 2-D. C G is a function of X and Y C IMPLICIT NONE INTEGER M ! 1st dimension input array of XYA INTEGER N ! 2nd dimension input array of XYA DOUBLE PRECISION X1 ! Lower X limit for integration DOUBLE PRECISION X2 ! Upper X limit for integration DOUBLE PRECISION SS ! Value of integral DOUBLE PRECISION XA(M) ! Holding array for X values for XYA DOUBLE PRECISION YA(N) ! Holding array for Y values for XYA DOUBLE PRECISION XYA(M,N) ! Holding array for tabulated function to be interpolated DOUBLE PRECISION XY2D(M,N) ! Holding array for 2nd derivatives of interpolant DOUBLE PRECISION YTMP(N) ! Holding array for 1-d array to be interpolated DOUBLE PRECISION UTMP(N) ! Holding array for interpolating function DOUBLE PRECISION Y2TMP(N) ! Holding array for 1-d array for 2nd derivatives of interpolant DOUBLE PRECISION YYTMP(M) ! Holding array for 1-d array of interpolated values C DOUBLE PRECISION G EXTERNAL G CALL QGAUSX(G,X1,X2,SS,M,N,XA,YA,XYA,XY2D,YTMP, : UTMP,Y2TMP,YYTMP) 99 END SUBROUTINE QGAUSX(FUNT,A,B,SS,M,N,XA,YA,XYA,XY2D, : YTMP,UTMP,Y2TMP,YYTMP) C C This is the Numerical Recipes Gaussian Quadrature C routine C IMPLICIT NONE INTEGER M INTEGER N DOUBLE PRECISION FUNT ! Function name to pass DOUBLE PRECISION A ! All as in QUAD2D DOUBLE PRECISION B DOUBLE PRECISION SS DOUBLE PRECISION XA(M) DOUBLE PRECISION YA(N) DOUBLE PRECISION XYA(M,N) DOUBLE PRECISION XY2D(M,N) DOUBLE PRECISION YTMP(N) DOUBLE PRECISION UTMP(N) DOUBLE PRECISION Y2TMP(N) DOUBLE PRECISION YYTMP(M) EXTERNAL FUNT C INTEGER J DOUBLE PRECISION DX DOUBLE PRECISION XM DOUBLE PRECISION XR DOUBLE PRECISION W(5) DOUBLE PRECISION V(5) SAVE W,V DATA W/0.2955242247D0,0.2692667193D0,0.2190863625D0, :0.1494513491D0,0.0666713443D0/ DATA V/0.1488743389D0,0.4333953941D0,0.6794095682D0, :0.8650633666D0,0.9739065285D0/ XM=0.50D0*(B+A) XR=0.50D0*(B-A) SS=0 DO J=1,5,1 DX=XR*V(J) SS=SS + W(J)*( : FUNT(XM+DX,M,N,XA,YA,XYA,XY2D,YTMP,UTMP,Y2TMP,YYTMP) + : FUNT(XM-DX,M,N,XA,YA,XYA,XY2D,YTMP,UTMP,Y2TMP,YYTMP) ) ENDDO SS=XR*SS 99 END DOUBLE PRECISION FUNCTION G(XX,M,N,XA,YA,XYA,XY2D, : YTMP,UTMP,Y2TMP,YYTMP) C C This is the required function subprogram of the C function to be integrated with respect to X C IMPLICIT NONE INTEGER M INTEGER N DOUBLE PRECISION XX ! Independent variable at which to evaluate G DOUBLE PRECISION XA(M) ! All as in QUAD2D DOUBLE PRECISION YA(N) DOUBLE PRECISION XYA(M,N) DOUBLE PRECISION XY2D(M,N) DOUBLE PRECISION YTMP(N) DOUBLE PRECISION UTMP(N) DOUBLE PRECISION Y2TMP(N) DOUBLE PRECISION YYTMP(M) C DOUBLE PRECISION F DOUBLE PRECISION Y1 DOUBLE PRECISION Y2 DOUBLE PRECISION X DOUBLE PRECISION Y DOUBLE PRECISION SS EXTERNAL F COMMON /XY/ X,Y X=XX CALL QGAUSY(F,Y1(X),Y2(X),SS,M,N,XA,YA,XYA,XY2D, : YTMP,UTMP,Y2TMP,YYTMP) G=SS 99 END SUBROUTINE QGAUSY(FUNT,A,B,SS,M,N,XA,YA,XYA,XY2D, : YTMP,UTMP,Y2TMP,YYTMP) C C This is the Numerical Recipes Gaussian Quadrature C routine C IMPLICIT NONE INTEGER M INTEGER N DOUBLE PRECISION FUNT ! Function name to pass DOUBLE PRECISION A ! All as QUAD2D DOUBLE PRECISION B DOUBLE PRECISION SS DOUBLE PRECISION XA(M) DOUBLE PRECISION YA(N) DOUBLE PRECISION XYA(M,N) DOUBLE PRECISION XY2D(M,N) DOUBLE PRECISION YTMP(N) DOUBLE PRECISION UTMP(N) DOUBLE PRECISION Y2TMP(N) DOUBLE PRECISION YYTMP(M) EXTERNAL FUNT C INTEGER J DOUBLE PRECISION DY DOUBLE PRECISION YM DOUBLE PRECISION YR DOUBLE PRECISION W(5) DOUBLE PRECISION V(5) SAVE W,V DATA W/0.2955242247D0,0.2692667193D0,0.2190863625D0, :0.1494513491D0,0.0666713443D0/ DATA V/0.1488743389D0,0.4333953941D0,0.6794095682D0, :0.8650633666D0,0.9739065285D0/ YM=0.5*(B+A) YR=0.5*(B-A) SS=0 DO J=1,5,1 DY=YR*V(J) SS=SS + W(J)*( : FUNT(YM+DY,M,N,XA,YA,XYA,XY2D,YTMP,UTMP,Y2TMP,YYTMP) + : FUNT(YM-DY,M,N,XA,YA,XYA,XY2D,YTMP,UTMP,Y2TMP,YYTMP) ) ENDDO SS=YR*SS 99 END DOUBLE PRECISION FUNCTION F(YY,M,N,XA,YA,XYA,XY2D, : YTMP,UTMP,Y2TMP,YYTMP) C C This is the required function subprogram of the C function to be integrated with respect to Y C IMPLICIT NONE INTEGER M INTEGER N DOUBLE PRECISION YY ! Independent variable at which to evaluate F DOUBLE PRECISION XA(M) ! All as QUAD2D DOUBLE PRECISION YA(N) DOUBLE PRECISION XYA(M,N) DOUBLE PRECISION XY2D(M,N) DOUBLE PRECISION YTMP(N) DOUBLE PRECISION UTMP(N) DOUBLE PRECISION Y2TMP(N) DOUBLE PRECISION YYTMP(M) C DOUBLE PRECISION FUNC DOUBLE PRECISION X DOUBLE PRECISION Y COMMON /XY/ X,Y Y=YY F=FUNC(X,Y,M,N,XA,YA,XYA,XY2D,YTMP,UTMP,Y2TMP,YYTMP) 99 END DOUBLE PRECISION FUNCTION FUNC(X,Y,M,N,XA,YA,XYA,XY2D, : YTMP,UTMP,Y2TMP,YYTMP) C C This is the required function subprogram which C supplies the `Z' values for given X,Y C IMPLICIT NONE INTEGER M INTEGER N DOUBLE PRECISION X ! 1st independent variable at which to evaluate FUNC DOUBLE PRECISION Y ! 2nd independent variable at which to evaluate FUNC DOUBLE PRECISION XA(M) ! All as QUAD2D DOUBLE PRECISION YA(N) DOUBLE PRECISION XYA(M,N) DOUBLE PRECISION XY2D(M,N) DOUBLE PRECISION YTMP(N) DOUBLE PRECISION UTMP(N) DOUBLE PRECISION Y2TMP(N) DOUBLE PRECISION YYTMP(M) C DOUBLE PRECISION Z CALL SPLIN2(M,N,XA,YA,XYA,XY2D,YTMP,UTMP,Y2TMP,YYTMP,X,Y,Z) FUNC=Z 99 END DOUBLE PRECISION FUNCTION Y1(X) C C This is the required function subprogram which C supplies the lower limit on Y as a function of X C IMPLICIT NONE DOUBLE PRECISION X ! Dummy independent variable C DOUBLE PRECISION LOWER COMMON/LOW/LOWER Y1=X*0.0D0 + LOWER 99 END DOUBLE PRECISION FUNCTION Y2(X) C C This is the required function subprogram which C supplies the upper limit on Y as a function of X C IMPLICIT NONE DOUBLE PRECISION X ! Dummy independent variable C DOUBLE PRECISION UPPER COMMON/UPP/UPPER Y2=X*0.0D0 + UPPER 99 END SUBROUTINE SPLIE2(M,N,X2A,YA,YTMP,Y2TMP,UTMP,Y2A) C C This is the Numerical recipes bicubic spline interpolation routine C C Given an M by N tabulated function YA, and tabulated C independent variable X2A (N values), this routine C constructs one-dimensional natural cubic splines C of the rows of YA and returns the second derivatives in the C array Y2A. C IMPLICIT NONE INTEGER M ! 1st dimension of function to be interpolated INTEGER N ! 2nd dimension of function to be interpolated DOUBLE PRECISION X2A(N) ! Array of Y values of YA DOUBLE PRECISION YA(M,N) ! Tabulated function to be interpolated DOUBLE PRECISION YTMP(N) ! Holding array for 1-d array to be interpolated DOUBLE PRECISION Y2TMP(N) ! Holding array for 1-d array of derivatives of interpolan DOUBLE PRECISION UTMP(N) ! Holding array for interpolating function DOUBLE PRECISION Y2A(M,N) ! Output array of 2nd derivatives of interpolant C INTEGER J INTEGER K DO J=1,M,1 DO K=1,N,1 YTMP(K)=YA(J,K) ENDDO CALL SPLINE(N,X2A,YTMP,2.D30,2.D30,UTMP,Y2TMP) DO K=1,N,1 IF (ABS(Y2TMP(K)).GT.1.0D-30) THEN Y2A(J,K)=Y2TMP(K) ELSE Y2A(J,K)=0.0D0 ENDIF ENDDO ENDDO 99 END SUBROUTINE SPLINE(N,X,Y,YP1,YPN,U,Y2) C C This is the Numerical Recipes Cubic Spline routine to C compute the array of second derivatives of the interpolating C function at the tabulated points. C C Given arrays X and Y of length N containing a tabulated function C i.e. y=f(x), with X1