C.IDENTIFICATION C C Program IUECONCAT version 1.1 C A. Fuente INSA/VILSPA May 21, 97. C C.KEYWORDS C C TABLE, ORDERS, IUE HIGH DISPERSION, C C.PURPOSE C C Rebin IUE High Dispersion Extracted Spectra & connect overlapping regions C C.ALGORITHM C C C.INPUT C C Keys: IN_A/C/1/60 input table C OUT_A/C/1/60 output table C COLW/C/1/60 wavelength col C COLO/C/1/60 order col C COLF/C/1/60 flux col C COLE/C/1/60 epsilon/quality col C CAME/I/1/1 camera number C FLAG/C/1/1 selection method flag C INPUTR/R/1/3 start,end,stepsize of the output table C C.OUTPUT C C C----------------------------------------------------------- C C PROGRAM IUECONCAT IMPLICIT NONE INTEGER ICOL,MAX,MAXORD PARAMETER (ICOL = 4) PARAMETER (MAX = 70000) PARAMETER (MAXORD = 125) C INTEGER COL(ICOL) INTEGER ADDR(ICOL) INTEGER INDEX(MAX) CHARACTER*60 TABIN CHARACTER*60 TABOUT CHARACTER*60 COLS(ICOL) INTEGER NCAM CHARACTER*1 FLAG,CAME REAL LMIN(MAXORD),LMAX(MAXORD) C CHARACTER*80 MSG INTEGER N REAL DIMS(3),START,STEP REAL RX(MAX) REAL RY(MAX) INTEGER REPS(MAX) INTEGER MASK INTEGER I,IACT,KUN,KNUL,ISTAT,Z INTEGER NCOL,NRO,NSC,NAC,NAR,NP INTEGER MADRID INTEGER TID C COMMON /VMR/MADRID(1) C DATA MSG /'*** ERROR (IUECONCAT):'/ C C ... get into MIDAS C CALL stspro('IUECONCAT') C C ... get input parameters C CALL STKRDC('IN_A',1,1,60,IACT,TABIN,KUN,KNUL,ISTAT) CALL STKRDC('OUT_A',1,1,60,IACT,TABOUT,KUN,KNUL,ISTAT) CALL STKRDC('COLW',1,1,60,IACT,COLS(1),KUN,KNUL,ISTAT) CALL STKRDC('COLO',1,1,60,IACT,COLS(2),KUN,KNUL,ISTAT) CALL STKRDC('COLF',1,1,60,IACT,COLS(3),KUN,KNUL,ISTAT) CALL STKRDC('COLE',1,1,60,IACT,COLS(4),KUN,KNUL,ISTAT) CALL STKRDC('FLAG',1,1,1,IACT,FLAG,KUN,KNUL,ISTAT) CALL STKRDC('CAME',1,1,1,IACT,CAME,KUN,KNUL,ISTAT) READ (CAME,FMT='(I1)') NCAM CALL STKRDR('INPUTR',1,3,IACT,DIMS,KUN,KNUL,ISTAT) C C ... init input table C CALL TBTOPN(TABIN,2,TID,ISTAT) IF (ISTAT.NE.0) GO TO 1000 CALL TBIGET(TID,NCOL,NRO,NSC,NAC,NAR,ISTAT) IF (ISTAT.NE.0) GO TO 1000 CALL TBCMAP(TID,0,MASK,ISTAT) IF (ISTAT.NE.0) GO TO 1000 C C ... find wavelength column i=1 C ... find order column i=2 C ... find flux column i=3 C ... find epsilon column i=4 C DO I=1,ICOL CALL TBCSER(TID,COLS(I),COL(I),ISTAT) IF (ISTAT.NE.0) GO TO 1000 IF (COL(I).EQ.-1) GO TO 1000 CALL TBCMAP(TID,COL(I),ADDR(I),ISTAT) IF (ISTAT.NE.0) GO TO 1000 ENDDO C C ... compute valid wavelength range for each order C C CALL ORDRANGE(NCAM,NRO,MADRID(ADDR(1)),MADRID(ADDR(2)), & MADRID(ADDR(3)),FLAG,LMIN,LMAX,Z) IF (Z.EQ.0) THEN WRITE(MSG,FMT='(A)')'*** ERROR: flux column not calibrated' CALL STTPUT(MSG,ISTAT) GOTO 1000 ENDIF C C ... select and sort the points C CALL ORDSEL(NRO,MADRID(ADDR(1)),MADRID(ADDR(2)), & LMIN,LMAX,INDEX,NP) C C ... rebin the table C START = DIMS(1) STEP = DIMS(3) N = (DIMS(2) - DIMS(1)) / STEP C CALL REBIN(START,STEP,N,MADRID(ADDR(1)),MADRID(ADDR(3)), & MADRID(ADDR(4)),INDEX,NP,RX,RY,REPS) C C ... write output table C CALL WRITEOUT(TABOUT,N,RX,RY,REPS,COLS(4)) C C ... end C CALL TBTCLO(TID, ISTAT) C C ... end C 1000 IF (ISTAT.NE.0) THEN CALL STETER(ISTAT,MSG) ENDIF 1002 CALL STSEPI ('IUECONCAT') END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE ORDSEL (NROW,WAVE,ORDER,LMIN,LMAX,INDEX,NP) C C------------------------------------------------------------------------------ C IMPLICIT NONE INTEGER NROW REAL WAVE(NROW) INTEGER ORDER(NROW) REAL LMIN(*),LMAX(*) INTEGER INDEX(*) INTEGER NP INTEGER I C C NP = 0 DO I = 1, NROW IF (WAVE(I).GE.LMIN(ORDER(I)).AND. & WAVE(I).LT.LMAX(ORDER(I))) THEN NP = NP + 1 INDEX(NP) = I ENDIF ENDDO C C ... sort the index in wavelength ascending order C CALL INDEXX(NP,WAVE,INDEX) C RETURN C END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE ORDRANGE(CAM,NROW,WAVE,ORDER,FLUX,FLAG,LMIN,LMAX,Z) C C--------------------------------------------------------------------------- C a. fuente Jul. 17, 1991. C define wavelength regions for overlapping orders: C note that the order number m, decreases with wavelength C C order m min(m) ---------------------max(m) C order m-1 min(m-1)----------------------max(m-1) C C ^ ^ ^ C L C R C C flag = L eft cut wavelength = min(m-1) C = R ight = max(m) C = C enter = (max(m) + min(m-1))/2. C = O ptimal (default) C SWP: the order ends have worse quality C LW : the order beginings have worse quality C = E chelle define wavelength regions for order m as C [lc(m+1/2),lc(m-1/2)] C coef. from NASA IUE Newsletters, July 82, n. 19, pag. 37-47 C Ake (1982) C IMPLICIT NONE INTEGER CAM,NROW REAL WAVE(*) INTEGER ORDER(*) REAL FLUX(*) CHARACTER*1 FLAG REAL LMIN(*),LMAX(*) INTEGER Z C REAL ZERO PARAMETER (ZERO = 0.0) C REAL MIN(125),MAX(125),DELTA INTEGER LAST_ORD(4),FIRST_ORD INTEGER M,MMIN,MMAX,STAT,I,NEXT REAL K1(4),K2(4),K3(4),CO,CN CHARACTER*80 RECORD,HD,HD2 C DATA FIRST_ORD /125/ DATA LAST_ORD /72,72,66,66/ DATA K1 /230036.0 ,230036.0 ,138827.0,0/ DATA K2 / 15.3456 ,15.3456 ,-27.426 ,0/ DATA K3 /-0.050638,-0.050638,0.165883,0 / IF (FLAG.EQ.'R'.OR.FLAG.EQ.'r') THEN CO = 1.00 CN = 1.00 ELSEIF (FLAG.EQ.'L'.OR.FLAG.EQ.'l') THEN CO = 0.00 CN = 0.00 ELSEIF (FLAG.EQ.'C'.OR.FLAG.EQ.'c') THEN CO = 0.50 CN = 0.50 ELSE C ! default optimal cuts IF (CAM.EQ.3) THEN C ! SWP CO = 0.30 CN = 0.23 ELSEIF (CAM.EQ.1.OR.CAM.EQ.2) THEN C !LWP & LWR CO = 0.75 CN = 0.65 ENDIF ENDIF DO M=FIRST_ORD, LAST_ORD(CAM), -1 MIN(M) = 5000.0 MAX(M) = 0.0 ENDDO C Z = 0 DO I=1,NROW IF (FLUX(I).NE.ZERO) THEN Z = Z + 1 IF (WAVE(I).GE.MAX(ORDER(I))) MAX(ORDER(I)) = WAVE(I) IF (WAVE(I).LE.MIN(ORDER(I))) MIN(ORDER(I)) = WAVE(I) ENDIF ENDDO C LMIN(FIRST_ORD) = MIN(FIRST_ORD) LMAX(LAST_ORD(CAM)) = MAX(LAST_ORD(CAM)) C IF (FLAG.EQ.'E'.OR.FLAG.EQ.'e') THEN DO M = FIRST_ORD, LAST_ORD(CAM), -1 MMIN = FLOAT(M) + 0.5 MMAX = FLOAT(M) - 0.5 LMIN(M) = (K1(CAM) / MMIN) + K2(CAM) + (K3(CAM) * MMIN) LMAX(M) = (K1(CAM) / MMAX) + K2(CAM) + (K3(CAM) * MMAX) ENDDO ELSE DO M=FIRST_ORD, LAST_ORD(CAM)+1, -1 NEXT = M-1 IF (MAX(M).GT.MIN(NEXT)) THEN C ! Overlapping DELTA = MAX(M) - MIN(NEXT) LMAX(M) = MIN(NEXT) + CO*DELTA LMIN(NEXT) = MIN(NEXT) + CN*DELTA ELSE LMAX(M) = MAX(M) LMIN(NEXT) = MIN(NEXT) ENDIF ENDDO ENDIF C HD =' Order Lambda_min Lambda_max Cut_min Cut_max ' HD2=' ----- ---------- ---------- ---------- ----------' C C ! Write in the MIDAS log too CALL STTPUT(HD,STAT) CALL STTPUT(HD2,STAT) DO M=FIRST_ORD,LAST_ORD(CAM), -1 IF (MIN(M).LT.5000) THEN WRITE(RECORD,1001) M,MIN(M),MAX(M),LMIN(M),LMAX(M) 1001 FORMAT(1X,I5,2X,F10.3,2X,F10.3,2X,F10.3,2X,F10.3) CALL STTPUT(RECORD,STAT) ENDIF ENDDO C RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE WRITEOUT(TABOUT,NP,WAVE,FLUX,EPSI,EPSLAB) C C----------------------------------------------------------------------------- C IMPLICIT NONE CHARACTER*60 TABOUT INTEGER NP REAL WAVE(*) REAL FLUX(*) INTEGER EPSI(*) C CHARACTER*(*) EPSLAB CHARACTER*16 LABEL(3) CHARACTER*16 UNIT(3) CHARACTER*8 FMT(3) INTEGER DATYPE(3) INTEGER IROW,C,TID,STAT,I C INCLUDE 'MID_INCLUDE:ST_DEF.INC' INCLUDE 'MID_INCLUDE:ST_DAT.INC' C DATA LABEL/'WAVELENGTH', 'FLUX', 'EPSILON'/ DATA UNIT /'ANGSTROM', 'ERG/CM2/A', 'UNITLESS'/ DATA FMT /'F8.3', 'E12.4', 'I5'/ DATYPE(1) = D_R4_FORMAT DATYPE(2) = D_R4_FORMAT DATYPE(3) = D_I4_FORMAT C C ... create the output file and columns C LABEL(3) = EPSLAB(1:16) CALL TBTINI(TABOUT,F_TRANS,F_IO_MODE,5,NP,TID,STAT) IF (STAT.NE.0) RETURN C DO I=1,3 CALL TBCINI(TID,DATYPE(I),1,FMT(I),UNIT(I),LABEL(I),C,STAT) ENDDO C C ... write columns C DO IROW=1,NP CALL TBEWRR(TID,IROW,1,WAVE(IROW),STAT) ENDDO C DO IROW=1,NP CALL TBEWRR(TID,IROW,2,FLUX(IROW),STAT) ENDDO C DO IROW=1,NP CALL TBEWRI(TID,IROW,3,EPSI(IROW),STAT) ENDDO C CALL TBTCLO(TID,STAT) C RETURN C END C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE REBIN(START,STEP,N,X,Y,EPS,INDEX,NP,RX,RY,REPS) C C---------------------------------------------------------------------- C IMPLICIT NONE REAL START REAL STEP INTEGER N REAL X(*) REAL Y(*) INTEGER EPS(*) INTEGER INDEX(*) INTEGER NP REAL RX(*) REAL RY(*) INTEGER REPS(*) C INTEGER GOODEPS PARAMETER (GOODEPS = 0) REAL TRESHOLD PARAMETER (TRESHOLD = 1/10.) C ! fraction of bad pixles into a bin C ! to be marked as bad bin C INTEGER I,K,PIX,BADPIX,TINULL,MINEPS,MAXEPS REAL TRNULL,STEP2 DOUBLE PRECISION SY,RXMIN,RXMAX,TDNULL C CALL TBMNUL(TINULL,TRNULL,TDNULL) C C ... C STEP2 = STEP / 2. C DO I=1,N RX(I) = START + STEP * (I-1) RY(I) = TRNULL REPS(I) = TINULL ENDDO C I = 1 RXMIN = START - STEP2 1000 CONTINUE IF (X(INDEX(I)).LE.RXMIN) THEN I = I + 1 IF (I.GT.NP) RETURN GOTO 1000 ENDIF C DO K = 1, N C! RXMIN = RX(K) - STEP2 RXMAX = RX(K) + STEP2 PIX = 0 BADPIX = 0 SY = 0 MINEPS = GOODEPS MAXEPS = GOODEPS 2000 CONTINUE IF (X(INDEX(I)).LE.RXMAX) THEN PIX = PIX + 1 SY = SY + Y(INDEX(I)) IF (EPS(INDEX(I)).LT.GOODEPS) BADPIX = BADPIX + 1 IF (MINEPS.GT.EPS(INDEX(I))) MINEPS = EPS(INDEX(I)) IF (MAXEPS.LT.EPS(INDEX(I))) MAXEPS = EPS(INDEX(I)) I = I + 1 IF (I.GT.NP) RETURN GOTO 2000 ENDIF IF (PIX.GT.0) THEN RY(K) = SY / PIX IF (BADPIX.GT.(PIX*TRESHOLD)) THEN REPS(K) = MINEPS ELSE REPS(K) = MAXEPS ENDIF ENDIF ENDDO C RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C SUBROUTINE INDEXX(N,ARRIN,INDX) C C--------------------------------------------------------------------------- C Sorts the index to an array ARRIN in ascending order C ARRIN values are not changed. C C From Numerical Recipes pag. 233 C IMPLICIT NONE INTEGER N REAL ARRIN(*) INTEGER INDX(N) C INTEGER L,IR,INDXT,I,J REAL Q L=N/2+1 IR=N 10 CONTINUE IF (L.GT.1) THEN L=L-1 INDXT=INDX(L) Q=ARRIN(INDXT) ELSE INDXT=INDX(IR) Q=ARRIN(INDXT) INDX(IR)=INDX(1) IR=IR-1 IF (IR.EQ.1) THEN INDX(1)=INDXT RETURN ENDIF ENDIF I=L J=L+L 20 IF (J.LE.IR) THEN IF (J.LT.IR) THEN IF (ARRIN(INDX(J)).LT.ARRIN(INDX(J+1))) J=J+1 ENDIF IF (Q.LT.ARRIN(INDX(J))) THEN INDX(I)=INDX(J) I=J J=J+J ELSE J=IR+1 ENDIF GO TO 20 ENDIF INDX(I)=INDXT GO TO 10 END