C ***************************** EWCL.F ******************************** C A VERSION OF EWIMH.FOR, CONVERTED TO RUN AS A CL-SCRIPT. GRAPHICS C ARE NON-INTERACTIVE 'IRAF/GRAPH'. BASED ON EWIMH2.FOR (DEC90). C ACL JUL94 C THIS VERSION EMPLOYS SIGMA-CLIPPING ON THE CONTINUUM BANDS. C AUTOMATICALLY REMOVES EFFECTS OF COSMIC RAYS, BAD COLUMNS, C AND MAYBE EVEN THE OCCASIONAL DEEP ABSORBTION LINE. ACL MAR96 C C Semi-Automatic Equivalent Width Computation Program with Plotting C Written by Taft Armandroff, December 1984 C VAX Version July 1987 by Gary DaCosta C Modified for use in IRAF environment by Andy Layden, November 1990 C & for variably-placed, self-selected continuum windows, AL Dec 1990 C CHARACTER FILE*15,LABEL*20,FEAT*6,CHA*1 LOGICAL*1 BEGIN DIMENSION FLX(4096),WAV(4096) INTEGER IAXLEN(7) COMMON/CTFXWV/NX,WPC,FLX,WAV DATA C/2.998E+5/ C Open file containing continuum windows, "ewimh.win" WRITE(6,8) 8 FORMAT('EWIMH: opening feature window file `ewimh.win`') OPEN(9,FILE='ewimh.win',STATUS='OLD') C First number in file gives number of windows to measure READ(9,*) NWINDOWS C Open output log: just append onto file "ewimh.log" WRITE(6,9) 9 FORMAT('EWIMH: appending to output logfile `ewimh.log`') OPEN(8,FILE='ewimh.log',ACCESS='APPEND') C Format is: C n feature line[A] rv[km/s] snr ew[A] {wav(AVC1) wav(AVC2)} C Chose an image and prompt for its radial velocity WRITE(6,30) 30 FORMAT(/,' Enter image name (eg obj102f): ',$) READ(5,100)FILE 100 FORMAT(A) WRITE(6,157) 157 FORMAT(' Enter the star`s radial velocity [km/s]: ',$) READ(5,*) RV C Get image information from header and read in spectrum using IMFORT CALL IMOPEN(FILE,1,IM,IER) IF(IER.NE.0) GOTO 1010 CALL IMGSIZ(IM,IAXLEN,NAXIS,IDTYP,IER) IF(IER.NE.0) GOTO 1010 IF(IDTYP.EQ.3) FEAT='short' IF(IDTYP.EQ.6) FEAT='real' CALL IMGKWC(IM,'title',LABEL,IER) IF(IER.NE.0) GOTO 1010 CALL IMGKWR(IM,'crval1',W0,IER) IF(IER.NE.0) GOTO 1010 CALL IMGKWR(IM,'cdelt1',WPC,IER) IF(IER.NE.0) GOTO 1010 CALL IMGKWI(IM,'np2',NX,IER) IF(IER.NE.0) GOTO 1010 CALL IMGL1R(IM,FLX,IER) IF(IER.NE.0) GOTO 1010 CALL IMCLOS(IM,IER) IF(IER.NE.0) GOTO 1010 WRITE(6,41) LABEL,FEAT,W0,WPC,NX 41 FORMAT(/,' TITLE: ',a20,' DATA TYPE: ',a5,' W0=',f8.2, & ' WPC=',f6.4,' Npts=',i5) C Compute wavelengths at each pixel number WAV(I) from W0 and WPC. DO 51 I=1,NX WAV(I) = W0+WPC*FLOAT(I-1) 51 CONTINUE C Write header line to EWIMH.LOG WRITE(8,42) FILE,LABEL 42 FORMAT(/,'FILE: ',A15,5X,'TITLE: ',A20) C Write header for screen output. WRITE(6,44) 44 FORMAT(/,' n feature line[A] rv[km/s] snr ', & 'ew[A] {Nsig wav(AVC1) wav(AVC2)}') C MAIN LOOP starts here: loop through program NWINDOWS times, C once for each window to be measured. NOTE: NWINDOWS should C match number of plotting iterations in EWCL.CL script. DO 123 IWIND=1,NWINDOWS READ(9,121) I,XLA,XLB,C1A,C1B,C2A,C2B,WF,FEAT C WRITE(6,121) I,XLA,XLB,C1A,C1B,C2A,C2B,WF,FEAT 121 FORMAT(I2,7F9.3,1X,A6) C 2 modes are available for continuum placement. Mode 0: use specified C values of C1A,C1B,C2A,C2B for continuum windows. Mode 1: C1A and C2A C give the starting wavelengths from which begin local searches for the C maximum average flux within a continuous interval of C1B,C2B Angstroms. C These average fluxes and thier average wavelengths are used to define C the continuum. Mode 0 is the default, Mode 1 is triggered if CB < CA. C C Pick mode based on continuum window values in EWIMH.WIN IF(C1B.LT.C1A) THEN MODE = 1 ELSE MODE = 0 ENDIF C Correct bandpasses for star's radial velocity RV C1A=C1A*(1.+RV/C) IF(MODE.EQ.0) C1B=C1B*(1.+RV/C) C2A=C2A*(1.+RV/C) IF(MODE.EQ.0) C2B=C2B*(1.+RV/C) XLA=XLA*(1.+RV/C) XLB=XLB*(1.+RV/C) C Compute max and min X values for plot window PLMIN=FLOAT(NINT(C1A)-50) PLMAX=FLOAT(NINT(C2B)+50) IF(MODE.EQ.1) PLMAX=FLOAT(NINT(C2A)+70) YMAX=0. C If MODE = 0, compute average fluxes in continuum bands. IF(MODE.EQ.0) THEN CALL M0CONT(C1A,C1B,AVC1) CALL M0CONT(C2A,C2B,AVC2) ELSE C If MODE = 1: Find continuum bands and compute average fluxes and wavlens. C Also, return final left and right continuum boundaries as C_A and C_B. CALL M1CONT(C1A,C1B,dsig1,AVC1) CALL M1CONT(C2A,C2B,dsig2,AVC2) ENDIF C Compute the SLOPE and y-intercept CONST of the pseudo-continuum line C between the two continuum bands. AVW1=(C1A+C1B)/2. AVW2=(C2A+C2B)/2. SLOPE=(AVC2-AVC1)/(AVW2-AVW1) CONST=AVC1-SLOPE*AVW1 C Compute signal to noise ratio in each cont bandpass and their weighted mean CALL COMPSNR(CONST,SLOPE,C1A,C1B,SNR1,NSNR1) CALL COMPSNR(CONST,SLOPE,C2A,C2B,SNR2,NSNR2) SNR = (SNR1*FLOAT(NSNR1)+SNR2*FLOAT(NSNR2))/FLOAT(NSNR1+NSNR2) C Open the GRAPH-ready file, using my dorkey but effective programming! IF(IWIND.EQ.1) CHA = '1' IF(IWIND.EQ.2) CHA = '2' IF(IWIND.EQ.3) CHA = '3' IF(IWIND.EQ.4) CHA = '4' IF(IWIND.EQ.5) CHA = '5' IF(IWIND.EQ.6) CHA = '6' IF(IWIND.EQ.7) CHA = '7' IF(IWIND.EQ.8) CHA = '8' IF(IWIND.EQ.9) CHA = '9' OPEN(10,FILE='tmp.'//CHA,STATUS='NEW') C This loop computes the number of pixels IPL and the maximum y-value YMAX C in the plot window, and dumps the spectrum into the GRAPH file. IPL=0 YMAX=0. DO 200 I=1,NX IF(WAV(I).LT.PLMIN.OR.WAV(I).GT.PLMAX)GOTO 200 IPL=IPL+1 WRITE(10,*) WAV(I),FLX(I) IF(IPL.EQ.1)J=I IF(FLX(I).GT.YMAX) YMAX=FLX(I) IF(WAV(I).GT.C2B.AND.WAV(I).GT.PLMAX)GOTO 210 200 CONTINUE 210 YMAX=1.1*YMAX C Draw continuum and feature band limits on the plot. C In GRAPH file, move (outside of frame) from end of spectrum to below C C2B, and draw C2B. WRITE(10,*) PLMAX,FLX(I) WRITE(10,*) PLMAX,-1. WRITE(10,*) C2B,-1. WRITE(10,*) C2B,2. C Move (outside of frame) from C2B to C2A, and draw C2A. WRITE(10,*) C2A,2. WRITE(10,*) C2A,-1. C Move (outside of frame) from C2A to XLB, and draw XLB. Make a C double line between XLA and XLB, to visually define feature window. WRITE(10,*) XLB,-1. WRITE(10,*) XLB,0.2 WRITE(10,*) XLA,0.2 WRITE(10,*) XLA,0.21 WRITE(10,*) XLB,0.21 WRITE(10,*) XLB,2. C Move (outside of frame) from XLB to XLA, and draw XLA. WRITE(10,*) XLA,2. WRITE(10,*) XLA,-1. C Move (outside of frame) from XLA to C1B, and draw C1B. WRITE(10,*) C1B,-1. WRITE(10,*) C1B,2. C Move (outside of frame) from C1B to C1A, and draw C1A. WRITE(10,*) C1A,2. WRITE(10,*) C1A,-1. C Now plot the continuum line (first move from CA1 to PLMIN). 587 WRITE(10,*) PLMIN,-1. WRITE(10,*) PLMIN,PLMIN*SLOPE+CONST WRITE(10,*) PLMAX,PLMAX*SLOPE+CONST C And close up the plot file. CLOSE(10) C Do some number counting... remnant of Taft's version. BEGIN=.FALSE. DO 280 K=J+1,J+IPL-1 IF(WAV(K).LT.XLA.OR.WAV(K).GT.XLB)GOTO 280 IF(.NOT.BEGIN) ISTART=K BEGIN=.TRUE. IEND=K 280 CONTINUE C Following Code is Taft's original numerical integration scheme. XLSUM = 0. 561 DO 290 K=ISTART,IEND+1 A=1.-(FLX(K)/(WAV(K)*SLOPE+CONST)) IF(I.EQ.ISTART) A=A*(WAV(K+1)-XLA)/WPC IF(I.EQ.IEND+1) A=A*(XLB-WAV(K-1))/WPC XLSUM=XLSUM+A 290 CONTINUE IMID=(ISTART+IEND)/2 WIDTH=XLSUM*(WAV(IMID)-WAV(IMID-1)) C 2002june19: a new bit so user can grade quality of ew measurement C WRITE(6,468) C 468 FORMAT(/,'Grade the quality of that EW (0=bad 4=best): ',$) C READ(5,*) IGRADE C And write everything to screen and into EWIMH.LOG IF(MODE.EQ.0) THEN WRITE(6,310) IWIND,FEAT,WF,RV,SNR,WIDTH WRITE(8,310) IWIND,FEAT,WF,RV,SNR,WIDTH ELSE WRITE(6,311) IWIND,FEAT,WF,RV,SNR,WIDTH, & dsig1,AVW1,dsig2,AVW2 WRITE(8,311) IWIND,FEAT,WF,RV,SNR,WIDTH, & dsig1,AVW1,dsig2,AVW2 ENDIF 310 FORMAT(1X,I2,3X,A6,3X,F9.3,3X,F5.0,3X,F6.2,3X,F5.2) 311 FORMAT(1X,I2,3X,A6,3X,F9.3,3X,F5.0,3X,F6.2,3X,F5.2,3X, & f4.1,2x,F9.3,3X,f4.1,2x,F9.3) C Go to next window 123 CONTINUE GOTO 1001 C Error exit from IMFORT subroutines. 1010 CALL IMEMSG(IER,LINE) WRITE(6,1005) LINE 1005 FORMAT(' IMFORT ERROR: ',A65) C Close it up and take it home 1001 CLOSE(9) CLOSE(8) C2000 CALL FINITT(0,700) STOP END C------------------------------------------------------------------------- C Subroutine to compute average flux within a continuum window via MODE 0. C CA and CB are the blue and red edges of the window in wavelength units. C C NOTE: THIS VERSION USES SIGMA-CLIPPING TO REJECT BAD PIXELS. C SUBROUTINE M0CONT(CA,CB,AVC) DIMENSION FLX(4096),WAV(4096) integer ibad(4096) LOGICAL*1 BEGIN COMMON/CTFXWV/NX,WPC,FLX,WAV SUM=0. POINTS=0. nitter = 10 dsig = 3. BEGIN=.FALSE. DO 190 I=1,NX IF(WAV(I).LT.CA.OR.WAV(I).GT.CB)GOTO 190 IF(.NOT.BEGIN) ISTART=I BEGIN=.TRUE. IEND=I 190 CONTINUE c set values of ibad (bad pix rejection flag) for each pixel to zero do 191 i=1,nx ibad(i) = 0 191 continue avc = 1. sigma = 9999. irej = 0 do 193 j=1,nitter points = 0. sum = 0. ssq = 0. DO 192 I=ISTART,IEND+1 FRACT=1. c if the pixel is rejected, skip it if(ibad(i).eq.1) goto 192 IF(I.EQ.ISTART) FRACT=(WAV(I)-CA)/WPC IF(I.EQ.IEND+1) FRACT=(CB-WAV(I-1))/WPC POINTS=POINTS+FRACT SUM=SUM+FRACT*FLX(I) ssq = ssq+fract*flx(i)*flx(i) c check to see if the pixel should be rejected in the next itter if(abs(flx(i)-avc).gt.dsig*sigma) then ibad(i) = 1 irej = irej+1 endif C write(6,*) i,wav(i),flx(i),ibad(i),sum,ssq 192 CONTINUE AVC=SUM/POINTS sigma = sqrt((ssq-sum*sum/points)/points) C write(6,*) j,avc,sigma,irej 193 continue 195 RETURN END C ------------------------------------------------------------------------ C Subroutine to compute average flux within a continuum window via MODE 1. C Starting at the first pixel (ISTART) with a wavelength greater than CA, C compute the average flux (AVFLX) between WAV(ISTART) and WAV(ISTART)+CB. C Step along 1 pixel and compute the average flux between WAV(ISTART+1) C and WAV(ISTART+1)+CB. Continue until WAV(I) > WAV(ISTART)+CMAX. Then C return the maximum AVFLX as AVC. Also return final starting and ending C window wavelengths as CA and CB (ie overwrite original values!). C SUBROUTINE M1CONT(CA,CB,dsig,AVC) DIMENSION FLX(4096),WAV(4096) integer ibad(4096) LOGICAL*1 BEGIN COMMON/CTFXWV/NX,WPC,FLX,WAV nitter = 5 c use a 2-sigma clip for CaK, since have few pixels (1-2 bad pix c are a much larger fraction of the Ntot of pix than for Hlines) dsig = 2. C Find the pixel numbers at CA (ISTART) and CA+CMAX (IEND). C NOTE: change CMAX to suit region of spectrum being searched. CMAX=25. BEGIN=.FALSE. DO 10 I=1,NX IF(WAV(I).LT.CA.OR.WAV(I).GT.(CA+CMAX)) GOTO 10 IF(.NOT.BEGIN) ISTART=I BEGIN=.TRUE. IEND=I 10 CONTINUE C Start the search. Outer loop advances the window 1 pixel at a time. C Inner loops compute average flux within the window using sigma-clip. AVC=0. ISTOP=IEND+INT(CB/WPC)+1 DO 30 ISTEP=ISTART,IEND POINTS=0. SUM=0. WSTART=WAV(ISTEP-1) WEND=WSTART+CB c set values of ibad (bad pix rejection flag) for each pixel to zero do 191 i=1,nx ibad(i) = 0 191 continue c loop to iterate through the sigma rejection... do it dumb... just c do nitter iterations. avflx = 1. sigma = 9999. irej = 0 do 193 j=1,nitter points = 0. sum = 0. ssq = 0. c loop to sum up good pixels within this itter DO 20 i=ISTEP,ISTOP FRACT=1. c if the pixel is rejected, skip it if(ibad(i).eq.1) goto 20 IF(WAV(i).GT.WEND) THEN FRACT=(WEND-WAV(i-1))/WPC IF(FRACT.LT.0.) GOTO 21 ENDIF POINTS=POINTS+FRACT SUM=SUM+FRACT*FLX(i) ssq = ssq+fract*flx(i)*flx(i) c check to see if the pixel should be rejected in the next itter if(abs(flx(i)-avflx).gt.dsig*sigma) then ibad(i) = 1 irej = irej+1 endif C write(6,*) i,wav(i),flx(i),ibad(i),sum,ssq 20 CONTINUE 21 AVFLX=SUM/POINTS sigma = sqrt((ssq-sum*sum/points)/points) C write(6,*) i,avflx,sigma,irej 193 continue C write(6,*) istep,istop,avflx,sigma,irej c if this value of avflx is higher than the previous ones, c make this the new best flux point (higher still implies c better, since itter SHOULD have removed all CR effects). IF(AVFLX.GT.AVC) THEN AVC=AVFLX CATEMP=WSTART CBTEMP=WEND ENDIF 30 CONTINUE CA=CATEMP CB=CBTEMP RETURN END C ----------------------------------------------------------------------- C Subroutine to compute the signal-to-noise ratio (and number of pixels C involved) in a continuum band, given upper and lower bounds, and slope C and intercept of pseudo-continuum line. Uses only whole pixels. C SUBROUTINE COMPSNR(CONST,SLOPE,CA,CB,SNR,NSNR) DIMENSION FLX(4096),WAV(4096) integer ibad(4096) COMMON/CTFXWV/NX,WPC,FLX,WAV c sigma-clipped version ... use dsig intermed between H and Ca defaults nitter = 5 dsig = 2.5 c set values of ibad (bad pix rejection flag) for each pixel to zero do 191 i=1,nx ibad(i) = 0 191 continue ave = 1. RMS = 9999. irej = 0 do 193 j=1,nitter SUM = 0. SUMSQ = 0. TOTFLX = 0. NSNR = 0 DO 10 I=1,NX IF(WAV(I).GT.CB) GOTO 10 if(ibad(i).eq.1) goto 10 IF(WAV(I).GT.CA) THEN DIFF = WAV(I)*SLOPE+CONST-FLX(I) SUM = SUM+DIFF SUMSQ = SUMSQ+DIFF*DIFF TOTFLX = TOTFLX+FLX(I) NSNR = NSNR+1 c check to see if the pixel should be rejected in the next itter if(abs(diff).gt.dsig*RMS) then ibad(i) = 1 irej = irej+1 endif C write(6,*) i,wav(i),flx(i),ibad(i),sum,ssq ENDIF 10 CONTINUE AVE = TOTFLX/FLOAT(NSNR) RMS = SQRT((SUMSQ-SUM*SUM/FLOAT(NSNR))/FLOAT(NSNR)) SNR = AVE/RMS C write(6,*) j,ave,rms,irej 193 continue RETURN END