C @(#)necopt.for 17.1.1.1 (ESO-DMD) 01/25/02 17:51:30 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, C MA 02139, USA. C C Correspondence concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C C+ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 21:13 - 3 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR : M Peron, S. Wolf C C.IDENTIFICATION C program ECHOPT version 1.0 910120 C C.KEYWORDS C C ECHELLE, CASPEC, DATA EXTRACTION C C.PURPOSE C C Extract orders from an echelle spectrum, the order position is defined by C the coefficients of a regression. The slit width and the orders to be C extracted are defined by the keyword INPUTI C The flux is calculated by assigning to each pixel a weighting factor C proportional to the amount of signal in it. C C.ALGORITHM C The algorithm is based on a paper of Koji Mukai ( 1990, Optimal Extraction of C Cross-Dispersed Spectra, Pub. of Astr.Soc. of the Pacific,102:183-189) C AND a paper of Keith Horne: C (1986, Pub. of Astr.Soc. of the Pacific, 98:609-617) C C.SYNTAX C EXTRACT/ORDER input output width,ord1,ord2 ron,g,sigma table coeffs C C the connection file contains C IN_A input FRM (has to be de-biased BEFORE) C OUTPUTC variance image of input data (same dimension as input frame) C OUT_A output FRM C OUT_B weight map FRM C IN_B table with positions of the orders C INPUTC coeffs C INPUTI slit,ord1,ord2, weight_determination_flag C INPUTR ron,g,sigma (ron: readout-noise, g: inverse gain, C sigma: threshold for cosmic ray rejection) C C.HISTORY: C 1999.04.08-SW save weight map, make use of variances, C iterative cosmic ray rejection (negative sigma switches off) C calculates the variance of the extracted spectrum C 1999.09.08-SW the variance of the extracted spectrum is calculated in C any case of INPUTI(4). C bug fix: the variance has to be devided by the squared C number of extracted lines (represented usually by the C keyword SLIT) in order to be consistent with the C extracted orders which are averaged data! The number of C lines is specified by INPUTI(1). C 1999.10.28-SW Do a gaussian fit to the profils (CALCPROFI(), C G_PROF()). C Save the cosmic ray detections in a table if EXTSIGMA>0 C 1999.11.09-SW Prevent division by zero in CALCULx() and CALC_x(). C 1999.11.11-SW Profile normalization removed in VSAMPLE() so that the C variances may be easily derived for the cosmic ray C rejection. C C C 010625 last modif C C*************************************************************************** C PROGRAM ECHOPT C IMPLICIT NONE C INTEGER NDIM PARAMETER (NDIM=2) C INTEGER MADRID, NAC, NAR, IMNO,SCAN(2),BINSIZE INTEGER I,IAV,INDX,IORD,IORD1,IORD2 INTEGER ICOL1,ICOL2 INTEGER ISTAT,NA,NAXISA,NAXISB,NAXISC INTEGER NCOL,NPAR,NROW,NSORT,NCR INTEGER STATUS,TID,MTID,INO,INVO INTEGER NPIXA(3),NPIXB(3),NPIXC(3) INTEGER IPAR(10),KUN,KNUL INTEGER IKIND,IWIDTH,DIMP INTEGER NWY,NWXY,WID,VID,DETWGT C INTEGER*8 PNTRA,PNTRB,PNTRC,PNTRV,PNTVO !pointers ... C REAL RPAR(10),WIDTH REAL PARAMS(12),OFFSET C DOUBLE PRECISION STEPA(3),STEPB(3),STARTA(3),STARTB(3) DOUBLE PRECISION DPAR(40) C CHARACTER*72 IDENA,IDENB,LINE CHARACTER*64 CUNITA,CUNITB,DESNAM CHARACTER*60 FRMIN,FRMOUT,FRMVOUT,FRMWGT,FRMVAR CHARACTER*60 INA1,INA2,INA3,INA4,TABLE CHARACTER*8 NAME,TYPE C COMMON /CRMASK/MTID,NCR INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC' C DATA CUNITB/ + 'FLUX PIXEL ORDER'/ C C initialize system C CALL STSPRO('ECHOPT') CALL STKRDC('IN_A',1,1,60,I,INA1,KUN,KNUL,STATUS) CALL STKRDC('IN_B',1,1,60,I,TABLE,KUN,KNUL,STATUS) CALL STKRDC('OUT_A',1,1,60,I,INA2,KUN,KNUL,STATUS) CALL STKRDC('OUT_B',1,1,60,I,INA3,KUN,KNUL,STATUS) CALL STKRDC('OUTPUTC',1,1,60,I,INA4,KUN,KNUL,STATUS) CALL STKRDC('INPUTC',1,1,8,I,NAME,KUN,KNUL,STATUS) CALL STKRDR('INPUTR',1,3,I,RPAR,KUN,KNUL,STATUS) CALL STKRDI('INPUTI',1,4,I,IPAR,KUN,KNUL,STATUS) CALL STKRDI('SCAN',1,2,I,SCAN,KUN,KNUL,STATUS) CALL STKRDR('OFFSET',1,1,I,OFFSET,KUN,KNUL,STATUS) c CALL STKRDI('BINSIZE',1,1,I,BINSIZE,N1,N1,N1) BINSIZE = 10 !binsize for virtual resampling VSAMPLE(). IORD1 = IPAR(2) IORD2 = IPAR(3) DETWGT = IPAR(4) !1: determine 0: use weight map CALL CLNFRA(INA1,FRMIN,0) !input frame (pixel-pixel-space) CALL CLNFRA(INA2,FRMOUT,0) !output frame (extracted orders) CALL CLNFRA(INA3,FRMWGT,0) !weight map (pixel-full_order-space) CALL CLNFRA(INA4,FRMVAR,0) !variance frame (pixel-pixel-space) CALL CLNTAB(TABLE,TABLE,0) !order table (order positions) C C input parameters in table C CALL TBTOPN(TABLE,F_U_MODE,TID,STATUS) CALL TBIGET(TID,NCOL,NROW,NSORT,NAC,NAR,STATUS) C C C map variance frame C CALL STIGET(FRMVAR,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE, . NDIM,NAXISA,NPIXA,STARTA,STEPA,IDENA, + CUNITA,PNTRV,VID,STATUS) C C map input FRM C CALL STIGET(FRMIN,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, . NDIM,NAXISA,NPIXA,STARTA,STEPA,IDENA, + CUNITA,PNTRA,IMNO,STATUS) PARAMS(1) = IPAR(1) !slitwidth (in spatial direction) PARAMS(2) = RPAR(3) !threshold factor for CR-rejection PARAMS(3) = RPAR(2) !invers gain factor [e-/ADU] PARAMS(4) = RPAR(1) !readout-noise [e-] C C read coeffs C INDX = INDEX(NAME//' ',' ') - 1 DESNAM = NAME(1:INDX)//'C' CALL STDRDC(TID,DESNAM,1,17,4,IAV,TYPE,KUN,KNUL,ISTAT) TYPE = 'MULT' IF (TYPE.EQ.'MULT') THEN DESNAM = NAME(1:INDX)//'I' CALL STDRDI(TID,DESNAM,4,4,IAV,IPAR,KUN,KNUL,ISTAT) NA = (IPAR(3)+1)* (IPAR(4)+1) DESNAM = NAME(1:INDX)//'R' CALL STDRDR(TID,DESNAM,1,4,IAV,PARAMS(5),KUN,KNUL,ISTAT) DESNAM = NAME(1:INDX)//'D' CALL STDRDD(TID,DESNAM,1,NA,IAV,DPAR,KUN,KNUL,ISTAT) C C ... map output FRM C IF (IORD1.EQ.IORD2 .AND. IORD1.EQ.0) THEN IORD1 = PARAMS(7) IORD2 = PARAMS(8) END IF STARTB(1) = STARTA(1) STARTB(2) = IORD1 STEPB(1) = STEPA(1) STEPB(2) = 1. IDENB = IDENA NAXISB = 2 NPIXB(1) = NPIXA(1) NPIXB(2) = IORD2 - IORD1 + 1 FRMVOUT = 'var_'//FRMOUT I = MAX(INDEX(FRMVOUT,' ')-1,1) CALL STIPUT(FRMOUT,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXISB,NPIXB,STARTB,STEPB,IDENB,CUNITB,PNTRB,INO,STATUS) WRITE(LINE,*) 'OUTPUT:',FRMOUT(1:I),' extracted spectrum' CALL STTPUT(LINE,STATUS) c c ... variance of extracted orders c IDENB = 'variance' CALL STIPUT(FRMVOUT,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXISB,NPIXB,STARTB,STEPB,IDENB,CUNITB,PNTVO,INVO,STATUS) WRITE(LINE,*) ' ',FRMVOUT(1:I) + ,' variance of extr. spectrum' CALL STTPUT(LINE,STATUS) C C ... map weight FRM C WIDTH = PARAMS(1) IWIDTH = WIDTH IF (MOD(IWIDTH,2).EQ.0)THEN IKIND=1 !even slit width DIMP = IWIDTH+1 !dimension of the weight profil ELSE IKIND=0 !odd slit width DIMP = IWIDTH+2 !dimension of the weight profil ENDIF NWY = DIMP*(IORD2 - IORD1 + 1) NWXY = NPIXA(1)*NWY STARTB(1) = STARTA(1) STARTB(2) = 1. STEPB(1) = STEPA(1) STEPB(2) = 1. IDENB = 'weight map' NAXISC = 2 NPIXC(1) = NPIXA(1) NPIXC(2) = NWY IF (DETWGT .EQ. 1) THEN CALL STIPUT(FRMWGT,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXISC,NPIXC,STARTB,STEPB, + IDENB,CUNITB,PNTRC,WID,STATUS) C C cosmic ray rejection ONLY if weights should be determined C IF (PARAMS(2).LT.0) THEN WRITE(LINE,*) ' NO cosmic ray rejection' CALL STTPUT(LINE,STATUS) MTID = -1 ELSE C create CR-mask table FRMOUT = 'crm_'//FRMIN I = MAX(INDEX(FRMOUT,' ')-4,1) FRMOUT = FRMOUT(1:I)//'tbl' I = I + 3 CALL TBTINI(FRMOUT,F_TRANS,F_O_MODE,2,100,MTID $ ,ISTAT) CALL TBCINI(MTID,D_I4_FORMAT,1,'I','PIXEL',':X',ICOL1 $ ,ISTAT) CALL TBCINI(MTID,D_I4_FORMAT,1,'I','PIXEL',':Y',ICOL2 $ ,ISTAT) NCR = 0 WRITE(LINE,*) ' Cosmic ray rejection (>',PARAMS(2) . ,'*sigma)' CALL STTPUT(LINE,STATUS) WRITE(LINE,*) ' save detections on: ' $ ,FRMOUT(1:I) CALL STTPUT(LINE,STATUS) ENDIF ELSE CALL STIGET(FRMWGT,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, . NDIM,NAXISC,NPIXC,STARTB,STEPB,IDENB, + CUNITB,PNTRC,WID,STATUS) ENDIF I = MAX(1,INDEX(FRMWGT,' ')-1) WRITE(LINE,*) ' ',FRMWGT(1:I),' weight map' CALL STTPUT(LINE,STATUS) I = MAX(1,INDEX(TABLE,' ')-1) WRITE(LINE,*) ' ',TABLE(1:I),' order table' CALL STTPUT(LINE,STATUS) WRITE(LINE,*) ' OFFSET to order positions:',OFFSET CALL STTPUT(LINE,STATUS) WRITE(LINE,*) ' BINSIZE of virtually resampled data:' $ ,BINSIZE CALL STTPUT(LINE,STATUS) WRITE(LINE,*) ' RON=',PARAMS(4),' [e], GAIN=', $ PARAMS(3),' [e/ADU]' CALL STTPUT(LINE,STATUS) C C extract order C NPAR = 11 PARAMS(9) = IPAR(3) PARAMS(10) = IPAR(4) IF (SCAN(1).EQ.0) THEN SCAN(1) = 1 SCAN(2) = NPIXA(2) ENDIF DO 10 I = IORD1,IORD2 WRITE(LINE,*) 'Order number : ',I CALL STTPUT(LINE,STATUS) PARAMS(11) = I IORD = I - IORD1 + 1 IF (DETWGT .EQ. 1) THEN c do an optimal extraction, save weights and variances on disk. CALL EXTRACT(MADRID(PNTRA),MADRID(PNTRV),NPIXA(1) $ ,NPIXA(2),MADRID(PNTRB),MADRID(PNTVO),NPIXB(1) $ ,NPIXB(2),MADRID(PNTRC),NPIXC(1),NPIXC(2),IORD $ ,NPAR,PARAMS,NA,DPAR,SCAN,OFFSET,BINSIZE,STARTA $ ,STEPA,IKIND,IWIDTH,DIMP) ELSE c do an immediate extraction using the weights of a previous optimal extract. CALL EXTRACT_I(MADRID(PNTRA),MADRID(PNTRV),NPIXA(1) $ ,NPIXA(2),MADRID(PNTRB),MADRID(PNTVO),NPIXB(1) $ ,NPIXB(2),MADRID(PNTRC),NPIXC(1),NPIXC(2),IORD $ ,NPAR,PARAMS,NA,DPAR,SCAN,OFFSET,STARTA,STEPA $ ,IKIND,IWIDTH,DIMP) ENDIF C 10 CONTINUE ELSE CALL STTPUT('Error in Dispersion coefficients',ISTAT) GOTO 20 ENDIF IF (MTID.GT.0) CALL TBTCLO(MTID,ISTAT) !close CR-mask table C CALL DSCUPT(IMNO,INO,' ',STATUS) C C free data C 20 CALL STSEPI END SUBROUTINE EXTRACT(IN,VARI,NPIXA1,NPIXA2,OUT,VOUT,NB1,NB2,WGT $ ,NC1,NC2,INDX,NPAR,PAR,ND,DPAR,SCAN,OFFSET,BINSIZE,STARTA $ ,STEPA,IKIND,IWIDTH,DIMP) C C C ARGUMENTS C IN REAL*4 INPUT IMAGE C VARI REAL*4 VARIANCE IMAGE C NPIXA1 INTG*4 NO. OF SAMPLES C NPIXA2 INTG*4 NO. OF LINES C OUT REAL*4 EXTRACTED ORDERS C VOUT REAL*4 VARIANCE OF EXTRACTED ORDERS C NB1 INTG*4 NO. OF SAMPLES C NB2 INTG*4 NO. OF ORDERS C WGT REAL*4 WEIGHT MAP C NC1 INTG*4 NO. OF SAMPLES C NC2 INTG*4 NO. OF ORDERS * THICKNESS OF ORDERS (profil dimension) C INDX INTG*4 SEQUENTIAL NO. OF THE ORDER C NPAR INTG*4 NO. OF EXTRACTION PARAMETERS C PAR REAL*4 EXTRACTION PARAMETERS C ND INTG*4 NO. OF COEFFS. C DPAR REAL*8 COEFFS. TO DEFINE ORDER POSITION C XMIN REAL*4 MINIMUM VALUE C XMAX REAL*4 MAXIMUMVALUE C IKIND INTG*4 EVEN/ODD FLAG: 1: EVEN, 0: ODD SLIT WIDTH SIZE C IWIDTH INTG*4 SLIT WIDTH (size perpendicular to dispersion in pixel) C DIMP INTG*4 DIMENSION OF THE PROFILE OF THE ORDERS C IMPLICIT NONE C INTEGER NPIXA1, NPIXA2, NB1, NB2, NC1, NC2, NPAR, ND INTEGER IKIND, IWIDTH, DIMP INTEGER NWW, IYLIM, NX, N1, N2, K, L, IX, INDX INTEGER NBY1,ISTAT,MADRID,DIM INTEGER NBY2,BINSIZE,IAC,NBY4 INTEGER FID1,FID2,FID3,FID4,FID5,FID6 INTEGER IR(2),SCAN(2),IXMIN,IXMAX INTEGER IXWIND,IOFF C INTEGER*8 PT1,PT2,PT3,PT4,PT5,PT6 !pointers C REAL WGT(NC1,NC2) REAL OUT(NB1,NB2),VOUT(NB1,NB2) REAL IN(NPIXA1,NPIXA2), VARI(NPIXA1,NPIXA2) REAL PAR(NPAR) REAL RON, CONAD, THRESH REAL XMIN,XMAX,OFFSET C DOUBLE PRECISION DPAR(ND),START(2),STEP(2) DOUBLE PRECISION WIDTH,HWIDTH DOUBLE PRECISION YV DOUBLE PRECISION SKWDX,SKWDY,SKWOF DOUBLE PRECISION DEGREE, STARTA(2), STEPA(2) CHARACTER*64 CUNIT INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC' DATA DEGREE/0.017453292519943295769237D+0/ DATA CUNIT/ + 'INTENSITY PIXEL ORDER'/ C NWW = 50 C C ... define extraction constants C WIDTH = IWIDTH HWIDTH = 0.5*WIDTH SKWDX = 0.0 SKWDY = 1.0 SKWOF = 0.0 XMIN = 0.0 XMAX = XMIN IYLIM = NPIXA2 - IWIDTH NX = SKWOF N1 = NX N2 = NPIXA1 - NX C C ... get RON and CONAD C THRESH = PAR(2) !cosmic ray rejection threshold CONAD = PAR(3) !inverse GAIN: CONAD [e-/ADU] RON = PAR(4) !readout-noise [e-] RON = (RON / CONAD)**2 !-> converted to ADU^2 (variance!) C C ... get ordinate values C K = PAR(9) !Dimension of the coefficient L = PAR(10) !matrix to get order pos. YV = PAR(11) !Current order C C ... get offset within weight buffer for corresponding order C IOFF = (INDX - 1) * DIMP + 1 C C ... loop over the samples C NBY1 = NPIXA1 CALL STFCRE('vdummx',D_R8_FORMAT,F_X_MODE,1,NBY1,FID1,ISTAT) CALL STFMAP(FID1,F_X_MODE,1,NBY1,IAC,PT1,ISTAT) CALL CENT(MADRID(PT1),NPIXA1,DPAR,ND,K,L,YV,IWIDTH,SCAN,OFFSET, 1 IXMIN,IXMAX,STARTA,STEPA) DIM = ((IWIDTH/2)*2+2)*BINSIZE+1 START(1) = 0 START(2) = 0 STEP(1) = 1. STEP(2) = 1./BINSIZE IR(1) = NPIXA1 IR(2) = DIM NBY2 = NPIXA1 * DIM CALL STFCRE('vdummy',D_R4_FORMAT,F_X_MODE,1,NBY2,FID2,ISTAT) CALL STFMAP(FID2,F_X_MODE,1,NBY2,IAC,PT2,ISTAT) CALL STFCRE('mask',D_I2_FORMAT,F_X_MODE,1,NBY2,FID3,ISTAT) CALL STFMAP(FID3,F_X_MODE,1,NBY2,IAC,PT3,ISTAT) IR(2) = 4 !parameter buffer NBY4 = IR(1) * IR(2) !size of the buffer CALL STFCRE('pdummy',D_R4_FORMAT,F_X_MODE,1,NBY4,FID4,ISTAT) CALL STFMAP(FID4,F_X_MODE,1,NBY4,IAC,PT4,ISTAT) IR(1) = DIM IR(2) = 1 !parameter buffer NBY4 = IR(1) * IR(2) !size of the buffer CALL STFCRE('colum',D_R4_FORMAT,F_X_MODE,1,NBY4,FID5,ISTAT) CALL STFMAP(FID5,F_X_MODE,1,NBY4,IAC,PT5,ISTAT) CALL STFCRE('sigma',D_R4_FORMAT,F_X_MODE,1,NBY4,FID6,ISTAT) CALL STFMAP(FID6,F_X_MODE,1,NBY4,IAC,PT6,ISTAT) DO IX=1,NPIXA1 OUT(IX,INDX) = 0 ENDDO CALL VSAMPLE(IN,WGT(1,IOFF),MADRID(PT1),MADRID(PT2), + MADRID(PT3),DIMP,IWIDTH,BINSIZE,NPIXA1,IXMIN, + IXMAX,NPIXA2,DIM,IKIND) IXWIND = 15 !window size of the order chunks CALL CALCPROFI(MADRID(PT2),MADRID(PT3),MADRID(PT4),MADRID(PT5) + ,MADRID(PT6),NPIXA1,IXWIND,IXMIN,IXMAX,DIM,4,RON,CONAD,THRESH + ,BINSIZE,WIDTH) IF (IKIND.EQ.0) THEN CALL CALCULI(IN,VARI,WGT(1,IOFF),MADRID(PT1),MADRID(PT2) + ,MADRID(PT3),OUT,VOUT,DIMP,IWIDTH,BINSIZE,NPIXA1,IXMIN + ,IXMAX,NPIXA2,DIM,NB1,NB2,INDX,PAR,NPAR) ELSE CALL CALCULP(IN,VARI,WGT(1,IOFF),MADRID(PT1),MADRID(PT2), + MADRID(PT3),OUT, + VOUT,DIMP,IWIDTH,BINSIZE,NPIXA1,IXMIN,IXMAX,NPIXA2,DIM, + NB1,NB2,INDX,PAR,NPAR) C ENDIF 99 CALL STFCLO(FID1,ISTAT) CALL STFCLO(FID2,ISTAT) CALL STFCLO(FID3,ISTAT) CALL STFCLO(FID4,ISTAT) CALL STFCLO(FID5,ISTAT) CALL STFCLO(FID6,ISTAT) RETURN C END SUBROUTINE EXTRACT_I(IN,VARI,NPIXA1,NPIXA2 + ,OUT,VOUT,NB1,NB2,WGT,NC1,NC2,INDX,NPAR,PAR + ,ND,DPAR,SCAN,OFFSET,STARTA,STEPA,IKIND,IWIDTH,DIMP) C C Do an immediate extract using the weights of a previous optimal extraction C C ARGUMENTS C IN REAL*4 INPUT IMAGE C VARI REAL*4 VARIANCE IMAGE C NPIXA1 INTG*4 NO. OF SAMPLES C NPIXA2 INTG*4 NO. OF LINES C OUT REAL*4 EXTRACTED ORDERS C VOUT REAL*4 VARIANCE OF EXTRACTED ORDERS C NB1 INTG*4 NO. OF SAMPLES C NB2 INTG*4 NO. OF ORDERS C WGT REAL*4 WEIGHT MAP C NC1 INTG*4 NO. OF SAMPLES C NC2 INTG*4 NO. OF ORDERS * THICKNESS OF ORDERS (profil dimension) C INDX INTG*4 SEQUENTIAL NO. OF THE ORDER C NPAR INTG*4 NO. OF EXTRACTION PARAMETERS C PAR REAL*4 EXTRACTION PARAMETERS C ND INTG*4 NO. OF COEFFS. C DPAR REAL*8 COEFFS. TO DEFINE ORDER POSITION C IKIND INTG*4 EVEN/ODD FLAG: 1: EVEN, 0: ODD SLIT WIDTH SIZE C IWIDTH INTG*4 SLIT WIDTH (size perpendicular to dispersion in pixel) C DIMP INTG*4 DIMENSION OF THE PROFILE OF THE ORDERS C IMPLICIT NONE C INTEGER NPIXA1, NPIXA2, NB1, NB2, NC1, NC2, NPAR, ND INTEGER IKIND, IWIDTH, DIMP INTEGER NWW, K, L, IX, INDX INTEGER NBY1,ISTAT,MADRID INTEGER IAC,FID1 INTEGER SCAN(2),IXMIN,IXMAX,IOFF C INTEGER*8 PT1 C REAL WGT(NC1,NC2) REAL OUT(NB1,NB2),VOUT(NB1,NB2) REAL IN(NPIXA1,NPIXA2), VARI(NPIXA1,NPIXA2) REAL PAR(NPAR) REAL OFFSET C DOUBLE PRECISION DPAR(ND) DOUBLE PRECISION YV DOUBLE PRECISION DEGREE, STARTA(2), STEPA(2) C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC' DATA DEGREE/0.017453292519943295769237D+0/ C NWW = 50 C C ... get ordinate values C K = PAR(9) !Dimension of the coefficient L = PAR(10) !matrix to get order pos. YV = PAR(11) !Current order C C ... get offset within weight buffer for corresponding order C IOFF = (INDX - 1) * DIMP + 1 C C ... loop over the samples C NBY1 = NPIXA1 CALL STFCRE('vdummx',D_R8_FORMAT,F_X_MODE,1,NBY1,FID1,ISTAT) CALL STFMAP(FID1,F_X_MODE,1,NBY1,IAC,PT1,ISTAT) c c get the centers of the orders: stored in yval[] == madrid(pt1) c CALL CENT(MADRID(PT1),NPIXA1,DPAR,ND,K,L,YV,IWIDTH,SCAN,OFFSET, 1 IXMIN,IXMAX,STARTA,STEPA) DO IX=1,NPIXA1 OUT(IX,INDX) = 0 ENDDO IF (IKIND.EQ.0) THEN CALL CALC_I(IN,VARI,WGT(1,IOFF),MADRID(PT1),OUT,VOUT $ ,NPIXA1,NPIXA2,NB1,NB2,DIMP,IXMIN,IXMAX,IWIDTH,INDX) ELSE CALL CALC_P(IN,VARI,WGT(1,IOFF),MADRID(PT1),OUT,VOUT $ ,NPIXA1,NPIXA2,NB1,NB2,DIMP,IXMIN,IXMAX,IWIDTH,INDX) C ENDIF C CALL STFCLO(FID1,ISTAT) C END SUBROUTINE CENT(YVAL,NPIXA1,DPAR,ND,K,L,YV,IWIDTH, + SCAN,OFFSET,IXMIN,IXMAX,STARTA,STEPA) C C Get the centres of the orders at each wavelength C IMPLICIT NONE C INTEGER K,L,IX1,L1,K1,ND,NPIXA1,IWIDTH,IP,IXMIN,IXMAX INTEGER SCAN(2),ITEST C DOUBLE PRECISION YVAL(NPIXA1),DPAR(ND),AX,XVAL,DC,YV DOUBLE PRECISION WW(50),STARTA(2), STEPA(2) DOUBLE PRECISION XI_R8R8, IX_R8R8 C REAL OFFSET C AX = 1 IXMIN = 1 IXMAX = NPIXA1 ITEST = 0 DO 50, IX1 = 1,NPIXA1 IP = 0 DC = 1.D0 YVAL(IX1) = 0.D0 XVAL = IX_R8R8(AX, STARTA(1), STEPA(1)) DO 30, L1 = 0,L IP = IP + 1 WW(IP) = DC YVAL(IX1) = YVAL(IX1)+ WW(IP)*DPAR(IP) DO 20, K1 = 1,K IP = IP + 1 WW(IP) = WW(IP-1)*XVAL YVAL(IX1) = YVAL(IX1) + WW(IP)*DPAR(IP) 20 CONTINUE DC = DC*YV 30 CONTINUE YVAL(IX1) = XI_R8R8(YVAL(IX1),STARTA(2),STEPA(2)) + OFFSET AX = AX+1 IF ( INT(YVAL(IX1)-IWIDTH/2-1).LE.SCAN(1)) THEN IXMIN = IX1 ELSE IF ( INT(YVAL(IX1)+IWIDTH/2+2).GE.SCAN(2).AND. + ITEST.EQ.0) THEN IXMAX = IX1 ITEST = 1 ENDIF 50 CONTINUE RETURN END SUBROUTINE CALCPROFI(FBVAL,MASK,PBUF,COLUM,SIGMA,NPIXA1,IXWIND + ,IXMIN,IXMAX,DIM,PDIM,RON,CONAD,THRESH,BINSIZE,WIDTH) C IMPLICIT NONE C INTEGER NUM,I1,I2,K,JMAX INTEGER DIM,PDIM,I,J,NPIXA1 INTEGER IXMIN,IXMAX,IXWIND,IXWMAX INTEGER BINSIZE,ISTAT,JWIN INTEGER LP(4), MFIT C INTEGER*2 MASK(NPIXA1,DIM),IMASK C REAL YMIN,YMAX REAL FBVAL(NPIXA1,DIM),PBUF(PDIM,NPIXA1) REAL RON,CONAD,THRESH REAL COLUM(DIM),SIGMA(DIM),PARAM(4),CHI2 REAL AMIN1,LOW,HGH C DOUBLE PRECISION WIDTH C DATA LP/1,2,3,4/ DATA MFIT/4/ c c Chop the order (FBVAL) into equaly spaced windows (IXWIND) in c wavelength space. c Do an averaging in the chunks to get averaged order profiles c Fit a gaussian to each of these averaged profiles. c Save the parameters of the gaussian in the first four rows c of FBVAL. c These parameters will finally be fitted in wavelength space, c this will be done by G_PROF(). c IMASK = 1 !mask: 1=valid pixel, 0:invalid LOW = BINSIZE*WIDTH/10.0 !these are empircal data !!! HGH = BINSIZE*WIDTH/3.0 !.. and should be replaced by measurements NUM = (IXMAX - IXMIN) / IXWIND I1 = IXMIN I2 = I1 + IXWIND - 1 DO K=1,NUM YMIN = 1E33 YMAX = -YMIN DO J = 1,DIM !Do the averaging within one window COLUM(J) = 0.0 DO I=I1,I2 COLUM(J) = COLUM(J) + FBVAL(I,J)*MASK(I,J) ENDDO COLUM(J) = COLUM(J) / IXWIND SIGMA(J) = (RON + ABS(COLUM(J))/CONAD) / IXWIND / IXWIND YMIN = AMIN1(COLUM(J),YMIN) !get the minimum IF (COLUM(J) .GT. YMAX) THEN !get maximum JMAX = J YMAX = COLUM(JMAX) ENDIF ENDDO C C Determine the initial parameters for gaussfit JWIN = 2 * BINSIZE PARAM(4) = (YMAX - YMIN) * 0.01 PARAM(1) = YMAX - PARAM(4) PARAM(2) = JMAX IF (JMAX.LT.JWIN .OR. DIM-JMAX.LT.JWIN) THEN IMASK = 0 GOTO 90 ENDIF PARAM(3) = DIM / 8.0 I = 4 CALL GAUSS_FIT(COLUM,SIGMA,DIM,PARAM,I,LP,MFIT,CHI2,ISTAT) IF (ISTAT.LT.0 .OR. + PARAM(3).LT.LOW .OR. PARAM(3).GT.HGH) IMASK = 0 90 I = IXMIN+K-1 PBUF(1,I) = 1.0 DO J=2,4 PBUF(J,I) = PARAM(J) MASK(I,J) = IMASK !mask out if parameters are BAD ENDDO 98 IMASK = 1 I1 = I2 + 1 I2 = I1 + IXWIND - 1 ENDDO c c Determine the gaussian profils at each wavelength. At the beginning c the gaussians have a zero background and an amplitude of 1.0. c IXWMAX = IXMIN + (IXMAX - IXMIN) / IXWIND - 1 CALL G_PROF(FBVAL,MASK,PBUF,COLUM,SIGMA,NPIXA1,IXWMAX,IXWIND,IXMIN + ,IXMAX,DIM,PDIM,RON,CONAD,THRESH,BINSIZE) c c Now determine background and amplitude for fixed c sigma and position of maximum c RETURN END SUBROUTINE G_PROF(FBVAL,MASK,PBUF,COLUM,SIGMA,NPIXA1,IXWMAX,IXWIND + ,IXMIN,IXMAX,DIM,PDIM,RON,CONAD,THRESH,BINSIZE) C IMPLICIT NONE C INTEGER DIM,PDIM,I,J,POLORD,NPIXA1 INTEGER IXWMAX,IXWIND,IXMIN,IXMAX,BINSIZE INTEGER I1,NOUT,ITER,MXITER C INTEGER*2 MASK(NPIXA1,DIM) C REAL FBVAL(NPIXA1,DIM),PBUF(PDIM,NPIXA1) REAL COLUM(DIM),SIGMA(DIM) REAL A(10),B(10),S(10),W(10),CHISQ,X,DY(4) REAL RON,CONAD,THRESH,THRESH2 C CC EXTERNAL FGAUSS_C DATA MXITER/5/ c c Fit the parameters of all the gaussians fitted to the averaged c order profils. As a result you get the parameters of a gaussian c at each wavelength, stored in PBUF[]. The gaussians itself are c saved in FBVAL[]. c c At first determin sigma and position of maximum of each gaussian C THRESH2 = THRESH**2 C I1 = IXMIN + 15 DO J=2,PDIM-1 POLORD = 2 CALL LSORTHO (PBUF,MASK,A,B,S,W,NPIXA1,I1,IXWMAX, 1 PDIM,J,CHISQ,POLORD,IXWIND) DO I = IXMIN,IXMAX CALL POLY(FLOAT(I),PBUF(J,I),A,B,S,W,POLORD) ENDDO ENDDO c The parameters of the gaussians at each wavelength are determined c Now calculate the gaussian profils and save them in FBVAL[]. DO I = IXMIN,IXMAX DO J = 1,DIM MASK(I,J) = 1 !initialize mask ENDDO ITER = 1 101 PBUF(1,I) = 1.0 !amplitude of gaussian PBUF(4,I) = 0.0 !background (additional constant) c The amplitude and the background of the gaussian c may be determined analytically DO J = 1,5 A(J) = 0.0 ENDDO DO J = 1,DIM IF (MASK(I,J).GT.0) THEN X = J CALL XGAUSS_C(X,PBUF(1,I),COLUM(J),DY) SIGMA(J) = RON + ABS(FBVAL(I,J))/CONAD B(1) = COLUM(J)/SIGMA(J) B(2) = 1.0/SIGMA(J) A(1) = A(1) + B(1) * COLUM(J) A(2) = A(2) + B(1) A(3) = A(3) + B(1) * FBVAL(I,J) A(4) = A(4) + B(2) A(5) = A(5) + B(2) * FBVAL(I,J) ENDIF ENDDO PBUF(4,I) = A(1)*A(4) - A(2)*A(2) IF (PBUF(4,I) .NE. 0) THEN PBUF(4,I) = (A(1)*A(5) - A(2)*A(3)) / PBUF(4,I) ELSE PBUF(4,I) = 1E33 WRITE(*,12345) I,A(1),A(2),A(3),A(4),A(5) ENDIF IF (A(1) .NE. 0) THEN PBUF(1,I) = (A(3) - A(2)*PBUF(4,I)) / A(1) ELSE PBUF(1,I) = 1E33 WRITE(*,12346) I,A(1),A(2),A(3),PBUF(4,I) ENDIF c Now the amplitude and the background are known c Determine the final gaussian profil c Do the kappa sigma clipping if enabled (THRESH>0) IF (THRESH.GT.0) THEN NOUT = 0 DO J = 1,DIM IF (MASK(I,J).GT.0) THEN X = J CALL XGAUSS_C(X,PBUF(1,I),COLUM(J),DY) IF (FBVAL(I,J).GT.0) THEN !only interested in CR !! DY(1) = MAX(0.0, FBVAL(I,J) - COLUM(J)) DY(2) = DY(1)*DY(1)/SIGMA(J) IF (DY(2).GT.THRESH2) THEN MASK(I,J) = 0 NOUT = NOUT+1 ENDIF ENDIF ENDIF ENDDO ITER = ITER + 1 IF (NOUT.GT.0.AND.ITER.LT.MXITER) GOTO 101 DO J = 1,DIM X = J CALL XGAUSS_C(X,PBUF(1,I),FBVAL(I,J),DY) ENDDO ELSE DO J = 1,DIM X = J CALL XGAUSS_C(X,PBUF(1,I),FBVAL(I,J),DY) ENDDO ENDIF ENDDO C 12345 FORMAT('(4) BAD profile at X = ',I5,5G13.7) 12346 FORMAT('(1) BAD profile at X = ',I5,4G13.7) C RETURN END SUBROUTINE PROFI(FBVAL,MASK,A,B,S,W,NPIXA1, + IXMIN,IXMAX,DIM,POLORD) C INTEGER DIM,I,J,POLORD,NPIXA1,IXMIN,IXMAX C INTEGER*2 MASK(NPIXA1,DIM) C REAL FBVAL(NPIXA1,DIM) REAL A(10),B(10),S(10),W(10),CHISQ c c old version to determine the order profiles c using a polynomial least squares fit of the resampled c intensity profils along the wavelength in each row of the c order. c DO J=1,DIM POLORD = 2 CALL LSORTH (FBVAL,MASK,A,B,S,W,NPIXA1,IXMIN,IXMAX, 1 DIM,J,CHISQ,POLORD,1) DO I = IXMIN,IXMAX CALL POLY(FLOAT(I),FBVAL(I,J),A,B,S,W,POLORD) ENDDO ENDDO RETURN END SUBROUTINE CALCULI(IN,VARI,PROF,YVAL,FBVAL,MASK,OUT,VOUT,DIMP, 1 IWIDTH,BINSIZE,NPIXA1,IXMIN,IXMAX, 2 NPIXA2,DIM,NB1,NB2,INDX,PAR,NPAR) C IMPLICIT NONE C INTEGER MTID,ISTAT INTEGER NPIXA1,NPIXA2,NB1,NB2,NPAR,I1 INTEGER I,J,K,KMAX,L,IWIDTH,DIM,BINSIZE,INDX,DIMP INTEGER INDY,FIND,IXMIN,IXMAX INTEGER NUM,NBAD,NCR C INTEGER*2 MASK(NPIXA1,DIM) C DOUBLE PRECISION YVAL(NPIXA1),X1,X2,X C REAL YPOS, YYY REAL FBVAL(NPIXA1,DIM),PROF(NPIXA1,DIMP) REAL SUM,SUM1,SUM2 REAL IN(NPIXA1,NPIXA2),VARI(NPIXA1,NPIXA2) REAL OUT(NB1,NB2),VOUT(NB1,NB2) REAL PAR(NPAR),RON,CONAD C CHARACTER*80 LINE C COMMON /CRMASK/MTID,NCR CONAD = PAR(3) !inverse GAIN: CONAD [e-/ADU] RON = PAR(4) !readout-noise [e-] RON = (RON / CONAD)**2 !-> converted to ADU^2 (used for variance!) NBAD = 0 C C Loop through the 'wavelength'-range [IXMIN,IXMAX] C DO I = IXMIN,IXMAX SUM = 0 I1 = INT(YVAL(I)-IWIDTH/2-1) c K = INT(YVAL(I)-IWIDTH/2) K = INT(YVAL(I)-IWIDTH/2.) KMAX = K+DIMP-1 YPOS = YVAL(I) C C Do the linear interpolation of the profile in the 'virtual space' C to get the profile in 'real space'. C At first find the nearest neighbours of X: (X1,X2) C X -> (I,J-K+1) is the profile position of interest. C DO J = K,KMAX-1 X = J + 0.5 c X = INT(YVAL(I)-IWIDTH/2)+J-K +0.5 c X = INT(YVAL(I)-IWIDTH/2.)+J-K +0.5 X1 = YVAL(I)-IWIDTH/2-1 + J-K X2 = X1 + 1./BINSIZE FIND = 0 10 IF (FIND.EQ.0) THEN IF ((X1.LE.X).AND.(X.LE.X2)) THEN FIND = 1 ELSE X1 = X2 X2 = X2 +1./BINSIZE ENDIF GOTO 10 ENDIF YYY = (X1-YVAL(I)+IWIDTH/2+1)*BINSIZE+1 INDY = NINT( YYY ) PROF(I,J-K+1) = FBVAL(I,INDY)*(X2-X)*BINSIZE+ 1 FBVAL(I,INDY+1)*(X-X1)*BINSIZE IF (MASK(I,INDY).EQ.0) THEN NCR = NCR + 1 !count total number of cosmics CALL TBEWRI(MTID,NCR,1,I,ISTAT) CALL TBEWRI(MTID,NCR,2,J,ISTAT) IN(I,J) = PROF(I,J-K+1) ENDIF ENDDO C C Do the extraction C Start with the normalization of the profile C SUM = 0 IF (NINT(YPOS).EQ.INT(YPOS)) THEN DO L = 1,DIMP-1 SUM = SUM + ABS(PROF(I,L)) ENDDO PROF(I,DIMP) = 0 !set to 0 as not used ELSE PROF(I,1) = 0 !set to 0 as not used DO L = 2,DIMP SUM = SUM + ABS(PROF(I,L)) ENDDO ENDIF IF (SUM.NE.0) THEN DO L = 1,DIMP !normalize profile PROF(I,L) = PROF(I,L)/SUM ENDDO ELSE !All PROF() are 0.0 SUM2 = 0 !--> skip OUT(I,INDX) = 0.0 VOUT(I,INDX) = 0.0 NBAD = NBAD+1 GOTO 998 ENDIF C C Recalculate the variance C DO J=K,KMAX L = J-K+1 IF (PROF(I,L).NE.0) 1 VARI(I,J)=RON+ABS(IN(I,K))/CONAD ENDDO C C Now the profile is determined. Extract the spectrum: C SUM2 = 0 !do the extraction ... SUM1 = 0 !... sum over spatial coordinate IF (NINT(YPOS).EQ.INT(YPOS)) THEN SUM1 = PROF(I,1)*IN(I,K)/VARI(I,K)/(DIMP-2)* 1 (NINT(YPOS-0.5)-YPOS+0.5) SUM2 = PROF(I,1)*PROF(I,1)/VARI(I,K) DO J=K+1,KMAX-2 L = J-K+1 SUM1 = SUM1 + PROF(I,L)*IN(I,J)/VARI(I,J)/(DIMP-2) SUM2 = SUM2 + PROF(I,L)*PROF(I,L)/VARI(I,J) ENDDO J = KMAX-1 SUM1 = SUM1+PROF(I,DIMP-1)*IN(I,J)/VARI(I,J)/(DIMP-2)* 1 (0.5+YPOS-NINT(YPOS-0.5)) SUM2 = SUM2+PROF(I,DIMP-1)*PROF(I,DIMP-1)/VARI(I,J) ELSE SUM1 = PROF(I,2)*IN(I,K+1)/VARI(I,K+1)/(DIMP-2)* 1 (1-(YPOS-0.5-NINT(YPOS-0.5))) SUM2 = PROF(I,2)*PROF(I,2)/VARI(I,K+1) DO J= K+2,KMAX-1 L = J-K+1 SUM1 = SUM1 + PROF(I,L)*IN(I,J)/VARI(I,J)/(DIMP-2) SUM2 = SUM2 + PROF(I,L)*PROF(I,L)/VARI(I,J) ENDDO J = KMAX SUM1 = SUM1 +PROF(I,DIMP)*IN(I,J)/VARI(I,J)/(DIMP-2)* 1 (YPOS-0.5-NINT(YPOS-0.5)) SUM2 = SUM2+PROF(I,DIMP)*PROF(I,DIMP)/VARI(I,J) ENDIF NUM = DIMP-2 OUT(I,INDX) = SUM1/SUM2 !extracted spectrum VOUT(I,INDX) = 1.0/SUM2/NUM/NUM !and its variance c To be consistent the variance has to be devided by the square of the c number of extracted lines as the extracted data are averged data! C C Now store the weights in PROF() and register all bad pixels in the table C 998 IF (NINT(YPOS).EQ.INT(YPOS)) THEN c IF (NINT(YVAL(I)).EQ.INT(YVAL(I))) THEN PROF(I,1) = PROF(I,1)/VARI(I,K)/SUM2* 1 (NINT(YPOS-0.5)-YPOS+0.5) DO J=K+1,KMAX-2 L = J-K+1 PROF(I,L) = PROF(I,L)/VARI(I,J)/SUM2 ENDDO J = KMAX-1 PROF(I,DIMP-1) = PROF(I,DIMP-1)/VARI(I,J)/SUM2* 1 (0.5+YPOS-NINT(YPOS-0.5)) ELSE PROF(I,2) = PROF(I,2)/VARI(I,K+1)/SUM2* 1 (1-(YPOS-0.5-NINT(YPOS-0.5))) DO J=K+2,KMAX-1 L = J-K+1 PROF(I,L) = PROF(I,L)/VARI(I,J)/SUM2 ENDDO J = KMAX PROF(I,DIMP) = PROF(I,DIMP)/VARI(I,J)/SUM2* 1 (YPOS-0.5-NINT(YPOS-0.5)) ENDIF ENDDO IF (NBAD.GT.0) THEN WRITE(LINE,*) NBAD,' bad profils found (extracted to 0.0).' CALL STTPUT(LINE,NBAD) ENDIF RETURN END SUBROUTINE CALCULP(IN,VARI,PROF,YVAL,FBVAL,MASK,OUT,VOUT,DIMP, + IWIDTH,BINSIZE,NPIXA1,IXMIN,IXMAX, + NPIXA2,DIM,NB1,NB2,INDX,PAR,NPAR) C INTEGER NPIXA1,NPIXA2,NB1,NB2,NPAR,I1 INTEGER I,J,K,KMAX,L,IWIDTH,DIM,BINSIZE,INDX,DIMP INTEGER INDY,FIND,IXMIN,IXMAX INTEGER NUM,NBAD,MTID,NCR C INTEGER*2 MASK(NPIXA1,DIM) C DOUBLE PRECISION YVAL(NPIXA1),X1,X2,X C REAL FBVAL(NPIXA1,DIM),PROF(NPIXA1,DIMP) REAL SUM,SUM1,SUM2 REAL IN(NPIXA1,NPIXA2),VARI(NPIXA1,NPIXA2) REAL OUT(NB1,NB2),VOUT(NB1,NB2) REAL PAR(NPAR),RON,CONAD C CHARACTER*80 LINE C COMMON /CRMASK/MTID,NCR C CONAD = PAR(3) !inverse GAIN: CONAD [e-/ADU] RON = PAR(4) !readout-noise [e-] RON = (RON / CONAD)**2 !-> convert to [ADU^2] DO I = IXMIN,IXMAX SUM = 0 I1 = INT(YVAL(I)-IWIDTH/2) K = I1 KMAX = K+DIMP-1 C C Do the linear interpolation of the profile in the 'virtual space' C to get the profile in 'real space'. C At first find the nearest neighbours of X: (X1,X2) C X -> (I,J-K+1) is the profile position of interest. C DO J = K,KMAX X = INT(YVAL(I)-IWIDTH/2)+J-K +0.5 X1 = YVAL(I)-IWIDTH/2 +J-K-1 X2 = YVAL(I)-IWIDTH/2 + J-K-1 +1./BINSIZE FIND = 0 10 IF (FIND.EQ.0) THEN IF ((X1.LE.X).AND.(X.LE.X2)) THEN FIND = 1 ELSE X1 = X2 X2 = X2 +1./BINSIZE ENDIF GOTO 10 ENDIF INDY = NINT((X1-YVAL(I)+IWIDTH/2+1)*BINSIZE+1) PROF(I,J-K+1)= FBVAL(I,INDY)*(X2-X)*BINSIZE+ 1 FBVAL(I,INDY+1)*(X-X1)*BINSIZE IF (MASK(I,INDY).EQ.0) THEN NCR = NCR + 1 !count total number of cosmics CALL TBEWRI(MTID,NCR,1,I,ISTAT) CALL TBEWRI(MTID,NCR,2,J,ISTAT) IN(I,J) = PROF(I,J-K+1) ENDIF ENDDO C C Do the extraction C Start with the normalization of the profile C SUM = 0 DO L = 1,DIMP SUM = SUM + PROF(I,L) ENDDO IF (SUM.NE.0) THEN DO L = 1,DIMP !normalize profile PROF(I,L) = PROF(I,L)/SUM ENDDO ELSE !All PROF() are 0.0 SUM2 = 0 !--> skip OUT(I,INDX) = 0.0 VOUT(I,INDX) = 0.0 NBAD = NBAD+1 GOTO 998 ENDIF C C Recalculate the variance C DO J=K,KMAX L = J-K+1 IF (PROF(I,L).NE.0) 1 VARI(I,J)=RON+ABS(IN(I,K))/CONAD c 1 VARI(I,J)=RON+ABS(OUT(I,INDX)*PROF(I,L))/CONAD ENDDO C C Now the profile is determined. Extract the spectrum: C SUM1 = PROF(I,1)*IN(I,K)/VARI(I,K)/(DIMP-1)* 1 (INT(YVAL(I))+1-YVAL(I)) SUM2 = PROF(I,1)*PROF(I,1)/VARI(I,K) DO J=K+1,KMAX-1 L = J-K+1 SUM1 = SUM1 + PROF(I,L)*IN(I,J)/VARI(I,J)/(DIMP-1) SUM2 = SUM2 + PROF(I,L)*PROF(I,L)/VARI(I,J) ENDDO J = K+DIMP-1 SUM1 = SUM1+PROF(I,DIMP)*IN(I,J)/VARI(I,J)/(DIMP-1)* 1 (YVAL(I)-INT(YVAL(I))) SUM2 = SUM2+PROF(I,DIMP)*PROF(I,DIMP)/VARI(I,J) NUM = DIMP-1 OUT(I,INDX) = SUM1/SUM2 !extracted spectrum VOUT(I,INDX) = 1.0/SUM2/NUM/NUM !and its variance c To be consistent the variance has to be devided by the square of the c number of extracted lines as the extracted data are averged data! C C Now store the weights in PROF() and register all bad pixels in the table C 998 PROF(I,1) = PROF(I,1)/VARI(I,K)/SUM2* 1 (INT(YVAL(I))+1-YVAL(I)) DO J=K+1,KMAX-1 L = J-K+1 PROF(I,L) = PROF(I,L)/VARI(I,J)/SUM2 ENDDO J = K+DIMP-1 PROF(I,DIMP) = PROF(I,DIMP)/VARI(I,J)/SUM2* 1 (YVAL(I)-INT(YVAL(I))) ENDDO IF (NBAD.GT.0) THEN WRITE(LINE,*) NBAD,' bad profils found (extracted to 0.0).' CALL STTPUT(LINE,NBAD) ENDIF RETURN END SUBROUTINE CALC_I(IN,VARI,PROF,YVAL,OUT,VOUT 1 ,NPIXA1,NPIXA2,NB1,NB2,DIMP,IXMIN,IXMAX,IWIDTH,INDX) INTEGER I,J,K,KMAX,L INTEGER IXMIN,IXMAX,IWIDTH,INDX,DIMP INTEGER NPIXA1,NPIXA2,NB1,NB2 REAL IN(NPIXA1,NPIXA2),VARI(NPIXA1,NPIXA2) REAL OUT(NB1,NB2),VOUT(NB1,NB2) REAL PROF(NPIXA1,DIMP),PART DOUBLE PRECISION SUM1,SUM2,SUM3 DOUBLE PRECISION YVAL(NPIXA1) DO I = IXMIN,IXMAX c K = INT(YVAL(I)-IWIDTH/2) K = INT(YVAL(I)-IWIDTH/2.0) KMAX = K+DIMP-1 C C Do the extraction C SUM3 = 0 SUM2 = 0 !do the extraction ... SUM1 = 0 !... sum over spatial coordinate IF (NINT(YVAL(I)).EQ.INT(YVAL(I))) THEN PART = NINT(YVAL(I)-0.5)-YVAL(I)+0.5 SUM1 = PROF(I,1)*IN(I,K)*PART SUM2 = PROF(I,1) SUM3 = PROF(I,1)*PROF(I,1)*VARI(I,K)*PART DO J=K+1,KMAX-2 L = J-K+1 SUM1 = SUM1 + PROF(I,L)*IN(I,J) SUM2 = SUM2 + PROF(I,L) SUM3 = SUM3 + PROF(I,L)*PROF(I,L)*VARI(I,J) ENDDO J = KMAX-1 PART = 0.5+YVAL(I)-NINT(YVAL(I)-0.5) SUM1 = SUM1 + PROF(I,DIMP-1)*IN(I,J)*PART SUM2 = SUM2 + PROF(I,DIMP-1) SUM3 = SUM3 + PROF(I,DIMP-1)*PROF(I,DIMP-1)*VARI(I,J)*PART ELSE PART = 1-(YVAL(I)-0.5-NINT(YVAL(I)-0.5)) SUM1 = PROF(I,2)*IN(I,K+1)*PART SUM2 = PROF(I,2) SUM3 = PROF(I,2)*PROF(I,2)*VARI(I,K+1)*PART DO J= K+2,KMAX-1 L = J-K+1 SUM1 = SUM1 + PROF(I,L)*IN(I,J) SUM2 = SUM2 + PROF(I,L) SUM3 = SUM3 + PROF(I,L)*PROF(I,L)*VARI(I,J) ENDDO J = KMAX PART = YVAL(I)-0.5-NINT(YVAL(I)-0.5) SUM1 = SUM1 + PROF(I,DIMP)*IN(I,J)*PART SUM2 = SUM2 + PROF(I,DIMP) SUM3 = SUM3 + PROF(I,DIMP)*PROF(I,DIMP)*VARI(I,J)*PART ENDIF IF (SUM2.NE.0) THEN OUT(I,INDX) = SUM1/SUM2/(DIMP-2) !extracted spect. VOUT(I,INDX) = SUM3/SUM2/SUM2/(DIMP-2)/(DIMP-2) !and its variance ELSE OUT(I,INDX) = 0.0 VOUT(I,INDX) = 0.0 ENDIF c To be consistent the variance has to be devided by the square of the c number of extracted lines as the extracted data are averged data! ENDDO RETURN END SUBROUTINE CALC_P(IN,VARI,PROF,YVAL,OUT,VOUT $ ,NPIXA1,NPIXA2,NB1,NB2,DIMP,IXMIN,IXMAX,IWIDTH,INDX) INTEGER I,J,K,KMAX,L INTEGER IXMIN,IXMAX,IWIDTH,IWIDTH2,INDX,DIMP INTEGER NPIXA1,NPIXA2,NB1,NB2 REAL IN(NPIXA1,NPIXA2),VARI(NPIXA1,NPIXA2) REAL OUT(NB1,NB2),VOUT(NB1,NB2) REAL PROF(NPIXA1,DIMP),PART DOUBLE PRECISION SUM1,SUM2,SUM3 DOUBLE PRECISION YVAL(NPIXA1) IWIDTH2 = IWIDTH/2 DO I = IXMIN,IXMAX SUM = 0 K = INT(YVAL(I)-IWIDTH2) KMAX = K+DIMP-1 C C Do the extraction C PART = INT(YVAL(I))+1-YVAL(I) SUM1 = PROF(I,1)*IN(I,K)*PART SUM2 = PROF(I,1) SUM3 = PROF(I,1)*PROF(I,1)*VARI(I,K)*PART DO J=K+1,KMAX-1 L = J-K+1 SUM1 = SUM1 + PROF(I,L)*IN(I,J) SUM2 = SUM2 + PROF(I,L) SUM3 = SUM3 + PROF(I,L)*PROF(I,L)*VARI(I,J) ENDDO J = K+DIMP-1 PART = YVAL(I)-INT(YVAL(I)) SUM1 = SUM1 + PROF(I,DIMP)*IN(I,J)*PART SUM2 = SUM2 + PROF(I,DIMP) SUM3 = SUM3 + PROF(I,DIMP)*PROF(I,DIMP)*VARI(I,DIMP)*PART IF (SUM2.NE.0) THEN OUT(I,INDX) = SUM1/SUM2/(DIMP-1) !extracted spect. VOUT(I,INDX) = SUM3/SUM2/SUM2/(DIMP-1)/(DIMP-1) !and its variance ELSE OUT(I,INDX) = 0.0 VOUT(I,INDX) = 0.0 ENDIF c To be consistent the variance has to be devided by the square of the c number of extracted lines as the extracted data are averged data! ENDDO RETURN END SUBROUTINE VSAMPLE(IN,PROF,YVAL,FBVAL,MASK,DIMP, 1 IWIDTH,BINSIZE,NPIXA1,IXMIN,IXMAX,NPIXA2,DIM,IKIND) C C Do the sampling in the virtual space. C INTEGER I,J,IWIDTH,DIM,IX,BINSIZE,IND,DIMP,IY INTEGER NPOINT,NPIXA1,NPIXA2,K,IKIND,IXMIN,IXMAX C INTEGER*2 MASK(NPIXA1,DIM) C DOUBLE PRECISION YVAL(NPIXA1),YDVAL,XDVAL C REAL FBVAL(NPIXA1,DIM) REAL PROF(NPIXA1,DIMP),DX,DY REAL IN(NPIXA1,NPIXA2),YDEC,SUM C DO J=1,DIM DO I=IXMIN,IXMAX FBVAL(I,J) = 0 MASK(I,J) = 0 ENDDO ENDDO DO J=1,DIMP DO I=IXMIN,IXMAX PROF(I,J) = 0 ENDDO ENDDO IF (IKIND.EQ.0) THEN DO I = IXMIN,IXMAX SUM = 0 K = INT(YVAL(I)-IWIDTH/2-1) PROF(I,1) = IN(I,K) DO J = K+1,K+DIMP-2 PROF(I,J-K+1) = IN(I,J) SUM = SUM+IN(I,J) ENDDO IF ( NINT(YVAL(I)).EQ.INT(YVAL(I))) THEN SUM = SUM + IN(I,K) ELSE SUM = SUM + IN(I,K+DIMP-1) ENDIF PROF(I,DIMP) = IN(I,K+DIMP-1) c IF (SUM.EQ.0) SUM=1 c DO J = K,K+DIMP-1 c PROF(I,J-K+1) = PROF(I,J-K+1)/SUM c ENDDO ENDDO ELSE DO I=IXMIN,IXMAX SUM = 0 K = INT(YVAL(I)-IWIDTH/2) DO J = K,K+DIMP-1 PROF(I,J-K+1) = IN(I,J) SUM = SUM+IN(I,J) ENDDO c IF (SUM.EQ.0) SUM = 1 c DO J = K,K+DIMP-1 c PROF(I,J-K+1) = PROF(I,J-K+1)/SUM c ENDDO ENDDO ENDIF IF (IKIND.EQ.0) THEN DO J=1,DIM NPOINT =0 DO I= IXMIN,IXMAX XDVAL = I IND = I K = INT(YVAL(I)-IWIDTH/2-1) YDVAL = YVAL(IND)-IWIDTH/2-1 + (J-1.)/BINSIZE IX = XDVAL YDEC = INT(YDVAL)+0.5 DX = XDVAL-IX DY = YDVAL-YDEC IY = INT(YDVAL) IF (DY.GT.-0.5.AND.DY.LE.0) THEN IF ((IY-K).GT.0) THEN FBVAL(I,J) = (1+DY)*PROF(IX,IY-K+1) - 1 DY*PROF(IX,IY-K) MASK(I,J) = 1 NPOINT = NPOINT+1 ENDIF ELSE IF (DY.LT.0.5.AND.DY.GE.0) THEN IF ((IY-K+1).LT.DIMP) THEN FBVAL(I,J) = (1-DY)*PROF(IX,IY-K+1) + 1 DY*PROF(IX,IY-K+2) MASK(I,J) = 1 NPOINT = NPOINT+1 ENDIF ENDIF ENDDO ENDDO ELSE DO J=1,DIM NPOINT =0 DO I= IXMIN,IXMAX XDVAL = I IND = I K = INT(YVAL(I)-IWIDTH/2) YDVAL = YVAL(IND)-IWIDTH/2-1 + (J-1.)/BINSIZE IX = XDVAL YDEC = INT(YDVAL)+0.5 DX = XDVAL-IX DY = YDVAL-YDEC IY = INT(YDVAL) IF (DY.GT.-0.5.AND.DY.LE.0.AND. 1 (IY-K+1).LE.DIMP.AND.(IY-K).GE.0) THEN IF ((IY-K).NE.0.AND.(IY-K).NE.DIMP) THEN FBVAL(I,J) = (1+DY)*PROF(IX,IY-K+1) - 1 DY*PROF(IX,IY-K) ELSE IF((IY-K).EQ.0) THEN FBVAL(I,J) = PROF(IX,IY-K+1) ELSE FBVAL(I,J) = PROF(IX,IY-K) ENDIF MASK(I,J) = 1 NPOINT = NPOINT+1 ELSE IF (DY.LT.0.5.AND.DY.GE.0.AND. 1 (IY-K+1).LE.DIMP.AND.(IY-K).GE.0) THEN IF ((IY-K+1).NE.DIMP.AND.(IY-K+1).NE.0) THEN FBVAL(I,J) = (1-DY)*PROF(IX,IY-K+1) + 1 DY*PROF(IX,IY-K+2) ELSE IF ((IY-K+1).EQ.DIMP) THEN FBVAL(I,J) = PROF(IX,IY-K+1) ELSE FBVAL(I,J) = PROF(IX,IY-K+2) ENDIF MASK(I,J) = 1 NPOINT = NPOINT+1 ENDIF ENDDO ENDDO ENDIF RETURN END SUBROUTINE LSORTHO (DAT,MASK,A,B,S,W,NPIXA1, 1 IXMIN,IXMAX,DIM,LIN,CHISQ,IORD,IXSTP) C IMPLICIT NONE C INTEGER IMAX,IMAXST,IFAUTO,IORD,I,N,I1,IFF,K,K1 INTEGER IFTEST,DIM,LIN,NPIXA1,IXMIN,IXMAX,IXSTP C INTEGER*2 MASK(NPIXA1,DIM) C REAL DAT(DIM,NPIXA1),A(10),B(10),S(10),W(10),P(10) REAL CHISQ,DEL,XNU,F,F95,DCHI,P2 REAL X C DATA IMAXST/10/,IFTEST/2/ C C IMAX=IMAXST IFAUTO=0 IF(IORD .LT. 10) GOTO 110 IORD=9 IFAUTO=1 110 CONTINUE IF(IORD .NE. 0) IMAX=MAX0(IABS(IORD)+1,2) DO 10 I=1,IMAXST W(I)=0. S(I)=0. A(I)=0. B(I)=0. 10 CONTINUE X = IXMIN - IXSTP P(1)=1. DO 20 N=IXMIN,IXMAX IF (MASK(N,LIN).EQ.1) THEN X = X + IXSTP W(1)=W(1)+1. S(1)=S(1)+DAT(LIN,N) A(1)=A(1)+X c A(1)=A(1)+N ENDIF 20 CONTINUE IF (W(1) .EQ. 0) RETURN S(1)=S(1)/W(1) A(1)=A(1)/W(1) I=1 XNU=INT(W(1))-1 30 IFF=1 40 I1=I IF(IMAX .GT. I) I = I+1 X = IXMIN - IXSTP CHISQ=0. DO 70 N=IXMIN,IXMAX X = X + IXSTP IF (MASK(N,LIN).EQ.1) THEN c P(2)=N-A(1) P(2)=X-A(1) IF(I .EQ. 2) GOTO 60 DO 50 K=3,I K1=K-1 P(K)=(X-A(K1))*P(K1)-B(K1)*P(K-2) c P(K)=(N-A(K1))*P(K1)-B(K1)*P(K-2) 50 CONTINUE 60 DEL=DAT(LIN,N)-S(I1)*P(I1) DAT(LIN,N)=DEL CHISQ=CHISQ+DEL*DEL IF(IMAX .LE. I1) GOTO 70 P2=P(I)**2 S(I)=S(I)+DEL*P(I) A(I)=A(I)+X*P2 c A(I)=A(I)+N*P2 W(I)=W(I)+P2 ENDIF 70 CONTINUE IF(IMAX .LE. I1) GOTO 80 A(I)=A(I)/W(I) B(I)=W(I)/W(I1) S(I)=S(I)/W(I) XNU=XNU-1. DCHI=S(I)**2*W(I) IF(DCHI .GE. CHISQ) GOTO 30 F=XNU*DCHI/(CHISQ-DCHI) F95=3.84+(10.+(12.+(30.+105./XNU/XNU)/XNU)/XNU)/XNU IF(F .GT. F95) GOTO 30 IF(IFAUTO.EQ.0) GOTO 30 ! F-TEST OFF XNU=XNU+1. IFF=IFF+1 S(I)=0. IF(IFF .LE. IFTEST) GOTO 40 80 IORD=MIN0(IMAX-1,I1)-IFF+1 RETURN END SUBROUTINE POLY(XFIT,YFIT,A,B,S,W,IORD) C IMPLICIT NONE INTEGER K,IORD REAL A(10),B(10),S(10),W(10),P(10) REAL XFIT,YFIT C P(1)=1. YFIT = S(1) C P(2)=XFIT-A(1) YFIT = YFIT+S(2)*P(2) C DO 10 K=3,IORD+1 P(K)=(XFIT-A(K-1))*(P(K-1))-B(K-1)*P(K-2) YFIT =YFIT+S(K)*P(K) 10 CONTINUE RETURN END SUBROUTINE XGAUSS_C(X,A,Y,DYDA) REAL X,Y,A(4),DYDA(4) DOUBLE PRECISION ARG,EX,FAC Y=0. ARG=(X-A(2))/A(3) EX=EXP(-ARG**2) FAC=A(1)*EX*2.*ARG Y=A(1)*EX + A(4) DYDA(1)=EX DYDA(2)=FAC/A(3) DYDA(3)=FAC*ARG/A(3) DYDA(4)=1.0 RETURN END c c OBSOLETE SUBROUTINES: c c SUBROUTINE LSORTH (FBVAL,MASK,A,B,S,W,NPIXA1, c 1 IXMIN,IXMAX,DIM,LIN,CHISQ,IORD) SUBROUTINE LSORTH (FBVAL,MASK,A,B,S,W,NPIXA1, 1 IXMIN,IXMAX,DIM,LIN,CHISQ,IORD,IXSTP) C IMPLICIT NONE INTEGER IMAX,IMAXST,IFAUTO,IORD,I,N,I1,IFF,K,K1 INTEGER IFTEST,DIM,LIN,NPIXA1,IXMIN,IXMAX,IXSTP INTEGER*2 MASK(NPIXA1,DIM) REAL FBVAL(NPIXA1,DIM),A(10),B(10),S(10),W(10),P(10) REAL CHISQ,DEL,XNU,F,F95,DCHI,P2 REAL X C DATA IMAXST/10/,IFTEST/2/ IMAX=IMAXST IFAUTO=0 IF(IORD .LT. 10) GOTO 110 IORD=9 IFAUTO=1 110 CONTINUE IF(IORD .NE. 0) IMAX=MAX0(IABS(IORD)+1,2) DO 10 I=1,IMAXST W(I)=0. S(I)=0. A(I)=0. B(I)=0. 10 CONTINUE X = IXMIN - IXSTP P(1)=1. DO 20 N=IXMIN,IXMAX IF (MASK(N,LIN).EQ.1) THEN X = X + IXSTP W(1)=W(1)+1. S(1)=S(1)+FBVAL(N,LIN) A(1)=A(1)+X c A(1)=A(1)+N ENDIF 20 CONTINUE IF (W(1) .EQ. 0) RETURN S(1)=S(1)/W(1) A(1)=A(1)/W(1) I=1 XNU=INT(W(1))-1 30 IFF=1 40 I1=I IF(IMAX .GT. I) I = I+1 X = IXMIN - IXSTP CHISQ=0. DO 70 N=IXMIN,IXMAX X = X + IXSTP IF (MASK(N,LIN).EQ.1) THEN c P(2)=N-A(1) P(2)=X-A(1) IF(I .EQ. 2) GOTO 60 DO 50 K=3,I K1=K-1 P(K)=(X-A(K1))*P(K1)-B(K1)*P(K-2) c P(K)=(N-A(K1))*P(K1)-B(K1)*P(K-2) 50 CONTINUE 60 DEL=FBVAL(N,LIN)-S(I1)*P(I1) FBVAL(N,LIN)=DEL CHISQ=CHISQ+DEL*DEL IF(IMAX .LE. I1) GOTO 70 P2=P(I)**2 S(I)=S(I)+DEL*P(I) A(I)=A(I)+X*P2 c A(I)=A(I)+N*P2 W(I)=W(I)+P2 ENDIF 70 CONTINUE IF(IMAX .LE. I1) GOTO 80 A(I)=A(I)/W(I) B(I)=W(I)/W(I1) S(I)=S(I)/W(I) XNU=XNU-1. DCHI=S(I)**2*W(I) IF(DCHI .GE. CHISQ) GOTO 30 F=XNU*DCHI/(CHISQ-DCHI) F95=3.84+(10.+(12.+(30.+105./XNU/XNU)/XNU)/XNU)/XNU IF(F .GT. F95) GOTO 30 IF(IFAUTO.EQ.0) GOTO 30 ! F-TEST OFF XNU=XNU+1. IFF=IFF+1 S(I)=0. IF(IFF .LE. IFTEST) GOTO 40 80 IORD=MIN0(IMAX-1,I1)-IFF+1 RETURN END SUBROUTINE DUMP(X,N1,N2) C IMPLICIT NONE INTEGER N1, N2 REAL X(N1, N2) INTEGER I1, I2, MIN C INTEGER I, J I1 = MIN(N1,5) I2 = MIN(N2,5) C C DO 20 I = 1, I2 C WRITE(*,1000) (X(J,I),J=1,I1) C 20 CONTINUE C RETURN 1000 FORMAT(1X,5E13.4) END