C********************************************************************** SUBROUTINE CALIB (CTS, COUNTS, NP, L0, TEXP, W0, DW, NW, # WVDISP, PXDISP, NDISP, # WVPHOT, FCPHOT, NPHOT, # Y2DISP, LAMBDA, FLUX, STATUS) C********************************************************************** *23456789 123456789 123456789 123456789 123456789 123456789 123456789 12 C C Version 2.0, 11 June 1992 - modified by W. Hack C - pass as variables the dispersion and phot. arrays C - photometrically correct the pixels THEN re-bins them C - dispersion curve limited to 75 elements C - photometric calibration curve limited to 2750 elements C C Version 1.0, October 9, 1989 C C Photometrically corrects raw 1D FOC objective prism spectra C from counts per pixel row to flux (ergs/s/cm2/Ang) per (equidistant) C wavelength bin. C C Input: C C REAL CTS(I) = real array containing summed counts in pixel row I=1,NP C REAL COUNTS(I) = CTS, but photometrically corrected; passed to rebin C C INTEGER NP = length of raw spectrum (normally NP=512) C [Note: raw spectrum is assumed to have opposite sense of L coordinate, C i.e., wavelength decreasing with increasing I. C C REAL L0 = L coordinate (in pixels) of undispersed image, measured C from L coordinate of start first pixel of raw spectrum in CTS C C [i.e., if start L=0, f/96 FUV Prism: L0= -89.5 for POS TARG 0,-7.6 C f/96 NUV Prism: L0=+156.0 for POS TARG 0,-2.2 C f/48 FUV Prism: L0= +83.3 for POS TARG 0,-7.6 C f/48 NUV Prism: L0=+226.5 for POS TARG 0,-1.3] C C REAL TEXP = exposure time of image in seconds C C REAL W0 = wavelength (in Ang) of first pixel in re-sampled spectrum C C REAL DW = binwidth (in Ang) of resampled spectrum C C INTEGER NW = length of re-sampled spectrum C C Output: C C REAL LAMBDA(NW) = array of central wavelengths of pixels of re-sampled C spectrum C REAL FLUX(NW) = array of fluxes (ergs/s/cm2/Ang) of re-sampled spectrum C C INTEGER STATUS = zero if OK C C Calls: C C SPLINE, LAMPIX, PHOTFC, REBIN, LLINT C C Calibration files: C C F96_FUV_DISP.DAT, F96_FUV_PHOT.DAT C or F96_NUV_DISP.DAT, F96_NUV_PHOT.DAT C or F48_FUV_DISP.DAT, F48_FUV_PHOT.DAT C or F48_NUV_DISP.DAT, F48_NUV_PHOT.DAT C C********************************************************************** INTEGER NP, NW, NDISP, NPHOT REAL CTS(NP), COUNTS(NP), LAMBDA(NW), FLUX(NW) REAL L0, TEXP, W0, DW REAL WVDISP(NDISP), PXDISP(NDISP), WVPHOT(NPHOT), FCPHOT(NPHOT) REAL Y2DISP(NDISP) INTEGER STATUS INTEGER I, J, NSAMPL REAL LAMPIX, PHOTFC REAL WAVE, WAVELO, WAVEHI, FAC STATUS = 0 C photometrically calibrate raw spectrum pixel by pixel DO 20 I = 1, NP WAVEHI = LAMPIX (L0 - FLOAT(I-1), WVDISP, PXDISP, Y2DISP, # NDISP, STATUS) IF (STATUS .NE. 0) RETURN WAVELO = LAMPIX (L0 - FLOAT(I), WVDISP, PXDISP, Y2DISP, # NDISP, STATUS) IF (STATUS .NE. 0) RETURN NSAMPL = INT((WAVEHI - WAVELO) / 5.) + 3 FAC = 0. DO 10 J = 1, NSAMPL WAVE = WAVELO + (WAVEHI - WAVELO) / 1 FLOAT(NSAMPL-1) * FLOAT(J-1) FAC = FAC + PHOTFC (WAVE, WVPHOT, FCPHOT, NPHOT) 10 CONTINUE FAC = FAC / FLOAT(NSAMPL) COUNTS(I) = CTS(I) * FAC / TEXP / DW 20 CONTINUE C rebin image CALL REBIN (COUNTS, NP, L0, W0, DW, NW, # WVDISP, PXDISP, NDISP, # WVPHOT, FCPHOT, NPHOT, # Y2DISP, LAMBDA, FLUX, STATUS) RETURN END C********************************************************************** SUBROUTINE REBIN (CTS, NP, L0, W0, DW, NW, # WVDISP, PXDISP, NDISP, # WVPHOT, FCPHOT, NPHOT, # Y2DISP, LAMBDA, FLUX, STATUS) C********************************************************************** C C Version 1.0, October 9 1989 C C Rebins raw 1D FOC objective prism spectra from counts per pixel row C to counts per (equidistant) wavelength bin. C The resampling is photon-conserving. C C Input: C C REAL CTS(I) = real array containing summed counts in pixel row I=1,NP C C INTEGER NP = length of raw spectrum (normally NP=512) C [Note: raw spectrum is assumed to have opposite sense of L coordinate, C i.e., wavelength decreasing with increasing I. C C REAL L0 = L coordinate (in pixels) of undispersed image, measured C from L coordinate of start first pixel of raw spectrum in CTS C C [i.e., if start L=0, f/96 FUV Prism: L0= -89.5 for POS TARG 0,-7.6 C f/96 NUV Prism: L0=+156.0 for POS TARG 0,-2.2 C f/48 FUV Prism: L0= +83.3 for POS TARG 0,-7.6 C f/48 NUV Prism: L0=+226.5 for POS TARG 0,-1.3] C C REAL W0 = wavelength (in Ang) of first pixel in re-sampled spectrum C C REAL DW = binwidth (in Ang) of resampled spectrum C C INTEGER NW = length of re-sampled spectrum C C Output: C C REAL LAMBDA(NW) = array of central wavelengths of pixels of re-sampled C spectrum C REAL FLUX(NW) = array of counts in re-sampled spectrum C C INTEGER STATUS = zero if OK C C Calls: C C SPLINE, SPLINT, PIXLAM C C Calibration files: C C F96_FUV_DISP.DAT C or F96_NUV_DISP.DAT C or F48_FUV_DISP.DAT C or F48_NUV_DISP.DAT C C********************************************************************** INTEGER NP, NW, NDISP, NPHOT REAL CTS(NP), LAMBDA(NW), FLUX(NW) REAL L0, W0, DW REAL WVDISP(NDISP), PXDISP(NDISP), WVPHOT(NPHOT), FCPHOT(NPHOT) REAL Y2DISP(NDISP) INTEGER STATUS INTEGER ILEFT, IRIGHT, I, J, K, ILAST REAL PLEFT, PRIGHT, PLAST, FLEFT, FRIGHT, PIXLAM STATUS = 0 C check requested wavelength limits * PIX0 = -W0 * W0 = LAMPIX (PIX0 + L0, STATUS) * IF (STATUS .NE. 0) RETURN * LMIN = LAMPIX (-512. + L0, STATUS) * IF (STATUS .NE. 0) RETURN * IF (W0.LT.LMIN) W0 = LMIN * LMAX = LAMPIX (L0, STATUS) * IF (STATUS .NE. 0) RETURN * IF (W0 + NW * DW.GT.LMAX) NW = INT((LMAX-W0) / DW) C create re-binned spectrum I = 0 PRIGHT = L0 - PIXLAM (W0-DW/2., # WVDISP, PXDISP, Y2DISP, NDISP, STATUS) IF (STATUS .NE. 0) RETURN IRIGHT = INT(PRIGHT) + 1 100 I = I + 1 IF (I.GT.NW) GOTO 200 LAMBDA(I) = W0 + DW * FLOAT(I-1) PLEFT = L0 - PIXLAM (LAMBDA(I) + DW/2., # WVDISP, PXDISP, Y2DISP, NDISP, STATUS) ILEFT = INT(PLEFT) + 1 IF (STATUS .NE. 0) RETURN IF (ILEFT .LT. IRIGHT) THEN C wavelength pixel spans two or more FOC pixels IF (ILEFT .LT. 1 .OR. ILEFT .GT. NP) THEN FLEFT = 0. ELSE FLEFT = (FLOAT(ILEFT) - PLEFT) * CTS(ILEFT) END IF IF (IRIGHT .LT. 1 .OR. IRIGHT .GT. NP) THEN FRIGHT = 0. ELSE FRIGHT = (PRIGHT - FLOAT(IRIGHT-1)) * CTS(IRIGHT) END IF FLUX(I) = FLEFT + FRIGHT IF (IRIGHT-ILEFT .GT. 1) THEN DO 110 J = MAX (ILEFT+1, 1), MIN (IRIGHT-1, NP) FLUX(I) = FLUX(I) + CTS(J) 110 CONTINUE ENDIF ELSE C wavelength bin must lie within one FOC pixel - find edge J = 0 120 PLAST = PLEFT ILAST = ILEFT J = J + 1 IF ((I+J).GT.NW) GOTO 150 LAMBDA(I+J) = W0 + DW * FLOAT(I+J-1) PLEFT = L0 - PIXLAM (LAMBDA(I+J) + DW/2., # WVDISP, PXDISP, Y2DISP, NDISP, STATUS) IF (STATUS .NE. 0) RETURN ILEFT = INT(PLEFT) + 1 IF (ILEFT.EQ.IRIGHT) GOTO 120 150 CONTINUE IF (IRIGHT .GE. 1 .AND. IRIGHT .LE. NP) THEN DO 170 K = 0, J-1 FLUX(I+K) = CTS(IRIGHT) * (PRIGHT - PLAST) / FLOAT(J) 170 CONTINUE ENDIF IF ((I+J).GT.NW) GOTO 200 IF (ILEFT .LT. 1 .OR. ILEFT .GT. NP) THEN FLEFT = 0. ELSE FLEFT = (FLOAT(ILEFT) - PLEFT) * CTS(ILEFT) END IF IF (ILAST .LT. 1 .OR. ILAST .GT. NP) THEN FRIGHT = 0. ELSE FRIGHT = (PLAST - FLOAT(ILAST-1)) * CTS(ILAST) END IF FLUX(I+J) = FLEFT + FRIGHT IF (ILAST-ILEFT.GT.1) THEN DO 180 K = MAX (ILEFT+1, 1), MIN (ILAST-1, NP) FLUX(I+J) = FLUX(I+J) + CTS(K) 180 CONTINUE ENDIF I = I + J ENDIF PRIGHT = PLEFT IRIGHT = ILEFT GOTO 100 200 CONTINUE RETURN END C********************************************************************** REAL FUNCTION PIXLAM (LAMBDA, WVDISP, PXDISP, # Y2DISP, NDISP, STATUS) C********************************************************************** C C Version 1.0, October 9 1989 C C Input: REAL LAMBDA = wavelength in Ang C C Output: REAL PIXLAM = offset from un-dispersed image in pixels C [Note: offset measured NEGATIVE in +L direction] C C INTEGER STATUS = zero if OK C C Calls: SPLINT (SPLINE assumed called from calling program) C C Calibration files: (assumed loaded from calling program) C C F96_FUV_DISP.DAT C or F96_NUV_DISP.DAT C or F48_FUV_DISP.DAT C or F48_NUV_DISP.DAT C C********************************************************************** REAL LAMBDA INTEGER NDISP REAL WVDISP(NDISP), PXDISP(NDISP) REAL Y2DISP(NDISP) INTEGER STATUS REAL OFF CALL SPLINT (WVDISP, PXDISP, Y2DISP, NDISP, LAMBDA, OFF, STATUS) PIXLAM = OFF RETURN END C********************************************************************** REAL FUNCTION LAMPIX (OFFSET, WVDISP, PXDISP, Y2DISP, # NDISP, STATUS) C********************************************************************** C C Version 1.0, October 9 1989 C C Input: REAL OFFSET = offset from un-dispersed image in pixels C [Note: offset measured NEGATIVE in +L direction] C C Output: REAL LAMPIX = wavelength in Ang C [Note: LAMPIX returned as first or last entry in C DISP.DAT table if out of range] C C INTEGER STATUS = zero if OK C C Calls: PIXLAM (assumed initialized from calling program) C C Calibration files: (assumed loaded from calling program) C C F96_FUV_DISP.DAT C or F96_NUV_DISP.DAT C or F48_FUV_DISP.DAT C or F48_NUV_DISP.DAT C C********************************************************************** REAL OFFSET INTEGER NDISP REAL WVDISP(*), PXDISP(*), Y2DISP(*) INTEGER STATUS REAL WAVE1, WAVE2, WAVMID, OFFTST, PIXLAM INTEGER J INTEGER K C Initial values. STATUS = 0 LAMPIX = 0. C find bracketing values in PXDISP array IF (OFFSET.LE.PXDISP(1))THEN LAMPIX = WVDISP(1) RETURN ELSE IF (OFFSET.GT.PXDISP(NDISP))THEN LAMPIX = WVDISP(NDISP) RETURN ELSE C J = 2 C DO WHILE ((OFFSET.GE.PXDISP(J)).AND.(J.LT.NDISP)) C J = J + 1 C END DO J = 2 DO 10 K = 2, NDISP-1 IF (OFFSET .GE. PXDISP(J)) THEN J = J + 1 ELSE GO TO 20 ENDIF 10 CONTINUE 20 CONTINUE WAVE1 = WVDISP(J-1) WAVE2 = WVDISP(J) ENDIF C check for good offset value OFFTST = PIXLAM (WAVE1, WVDISP, PXDISP, Y2DISP, NDISP, # STATUS) - OFFSET IF (STATUS .NE. 0) RETURN WAVMID = PIXLAM (WAVE2, WVDISP, PXDISP, Y2DISP, NDISP, # STATUS) - OFFSET IF (STATUS .NE. 0) RETURN IF (OFFTST * WAVMID.GT.0) THEN C 'Bad offset value in LAMPIX' STATUS = 1 RETURN ENDIF C check sense of slope IF (OFFTST.GT.0) THEN WAVMID = WAVE2 WAVE2 = WAVE1 WAVE1 = WAVMID ENDIF C find LAMPIX through bisection DO 30 J = 1, 250 WAVMID = (WAVE1 + WAVE2) / 2. OFFTST = PIXLAM (WAVMID, WVDISP, PXDISP, Y2DISP, NDISP, # STATUS) - OFFSET IF (STATUS .NE. 0) RETURN IF (OFFTST.LT.0.) THEN WAVE1 = WAVMID ELSE WAVE2 = WAVMID ENDIF IF ((ABS(OFFTST).LT.0.005) # .OR. (ABS(WAVE1 - WAVE2).LT.0.01)) THEN LAMPIX = WAVMID RETURN ENDIF 30 CONTINUE C 'No convergence in LAMPIX' STATUS = 2 LAMPIX = 0. RETURN END C********************************************************************** REAL FUNCTION PHOTFC (LAMBDA, WVPHOT, FCPHOT, NPHOT) C********************************************************************** C C Version 1.0, October 9 1989 C modified July 1, 1994 by PEH: C Extrapolate using a constant end value. Previously llint was C linearly extrapolating, which resulted in huge values. C C Interpolate within photometric conversion curve FCPHOT. WVPHOT is the C independent variable array, and FCPHOT is the dependent variable array. C C Input: REAL LAMBDA = wavelength in Ang C REAL WVPHOT = array of wavelengths C REAL FCPHOT = photometric conversion factor array C INTEGER NPHOT = size of arrays WVPHOT and FCPHOT. C C Output: REAL PHOTFC = photometric conversion factor C [ergs/cm2 per count] C C Calls: LLINT C C Calibration files: (assumed loaded from calling program) C C F96_FUV_PHOT.DAT C or F96_NUV_PHOT.DAT C or F48_FUV_PHOT.DAT C or F48_NUV_PHOT.DAT C C********************************************************************** REAL LAMBDA, WVPHOT(*), FCPHOT(*) INTEGER NPHOT REAL FAC IF (LAMBDA .LE. WVPHOT(1)) THEN FAC = FCPHOT(1) ELSE IF (LAMBDA .GE. WVPHOT(NPHOT)) THEN FAC = FCPHOT(NPHOT) ELSE CALL LLINT (WVPHOT, FCPHOT, NPHOT, LAMBDA, FAC) END IF PHOTFC = FAC RETURN END C********************************************************************** SUBROUTINE LLINT (X, Y, NPTS, XP, YP) C********************************************************************** C C Log-linear interpolation routine C REAL X, Y, XP, YP, X1, X2, Y1, Y2, A0, A1 INTEGER NPTS, J, J1, J2 INTEGER K DIMENSION X(NPTS), Y(NPTS) C Interpolate. IF (XP.LE.X(2))THEN J1 = 1 J2 = 2 ELSE IF (XP.GE.X(NPTS-1))THEN J1 = NPTS - 1 J2 = NPTS ELSE C J = 1 C DO WHILE ((XP.GE.X(J)).AND.(J.LT.NPTS)) C J = J + 1 C END DO J = 1 DO 10 K = 1, NPTS-1 IF (XP .GE. X(J)) THEN J = J + 1 ELSE GO TO 20 ENDIF 10 CONTINUE 20 CONTINUE IF (J.EQ.1)THEN J1 = J J2 = J + 1 ELSE IF (J.EQ.NPTS)THEN J1 = J - 1 J2 = J ELSE J1 = J - 1 J2 = J ENDIF ENDIF X1 = X(J1) Y1 = LOG(Y(J1)) X2 = X(J2) Y2 = LOG(Y(J2)) A1 = (Y1 - Y2) / (X1 - X2) A0 = Y1 - X1 * A1 YP = EXP(A0 + A1 * XP) RETURN END C********************************************************************** SUBROUTINE SPLINT (XA, YA, Y2A, N, X, Y, STATUS) C********************************************************************** C C Cubic spline interpolation routine C straight from Press et al. "Numerical Recipies" C except for extra argument STATUS C [Note: SPLINE must be called first and 2nd derivative array Y2=Y2A C passed to SPLINT] C Output: C INTEGER STATUS = zero of OK C DIMENSION XA(N), YA(N), Y2A(N) INTEGER STATUS KLO = 1 KHI = N 1 IF (KHI-KLO.GT.1) THEN K = (KHI + KLO) / 2 IF (XA(K).GT.X)THEN KHI = K ELSE KLO = K ENDIF GOTO 1 ENDIF H = XA(KHI) - XA(KLO) IF (H.EQ.0.) THEN C 'Bad XA input.' STATUS = 3 RETURN ENDIF A = (XA(KHI) - X) / H B = (X - XA(KLO)) / H Y = A * YA(KLO) + B * YA(KHI) + * ((A**3 - A) * Y2A(KLO) + (B**3 - B) * Y2A(KHI)) * (H**2) / 6. STATUS = 0 RETURN END C********************************************************************** * Note that U has been added to the calling sequence. This is scratch space. SUBROUTINE SPLINE (X, Y, N, YP1, YPN, U, Y2) C********************************************************************** C C Cubic spline 2nd derivative calculator C straight from Press et al. "Numerical Recipies" C C [Note: SPLINE must be called first and 2nd derivative array Y2=Y2A C passed to SPLINT] C DIMENSION X(N), Y(N), Y2(N), U(N) IF (YP1 .GT. 0.99E30) THEN Y2(1) = 0. U(1) = 0. ELSE Y2(1) = -0.5 U(1) = (3. / (X(2) - X(1))) * ((Y(2) - Y(1)) / 1 (X(2) - X(1)) - YP1) ENDIF DO 11 I = 2, N-1 SIG = (X(I) - X(I-1)) / (X(I+1) - X(I-1)) P = SIG * Y2(I-1) + 2. Y2(I) = (SIG-1.) / P U(I) = (6. * 1 ((Y(I+1) - Y(I)) / (X(I+1) - X(I)) - 2 (Y(I) - Y(I-1)) / (X(I) - X(I-1))) / 3 (X(I+1) - X(I-1)) - SIG * U(I-1)) / P 11 CONTINUE IF (YPN .GT. 0.99E30) THEN QN = 0. UN = 0. ELSE QN = 0.5 UN = (3. / (X(N) - X(N-1))) * (YPN - (Y(N) - Y(N-1)) / 1 (X(N) - X(N-1))) ENDIF Y2(N) = (UN - QN * U(N-1)) / (QN * Y2(N-1) + 1.) DO 12 K = N-1, 1, -1 Y2(K) = Y2(K) * Y2(K+1) + U(K) 12 CONTINUE RETURN END