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
C     This version uses different HEADER KEYWORDS than the old 1994 version.
C     Apparently the old version worked with 1D spectra.  Newer spectra
C     extracted with APALL have 3 axes with multiple spectra (the one you)
C     want should be in position [1]), and so some keywords changed.  You 
C     may have to change them again for the keywords in your spectra: see 
C     the comments with "***" in the "Get image information" section, and 
C     recompile (in IRAF, fc).  This version does not need to be "pre-
C     treated" by inserting the keyword "np2 = nnnn" into the header.
C               ACL Oct2014

       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
C   Gets image size, number of axes (3 for apall), and data type ***
       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'
C   Gets image title and writes in into chararcter string LABEL ***
       CALL IMGKWC(IM,'title',LABEL,IER)
       IF(IER.NE.0) GOTO 1010
C   Gets wavelength of first pixel and writes it into real# W0 ***
       CALL IMGKWR(IM,'crval1',W0,IER)
       IF(IER.NE.0) GOTO 1010
C   Gets the wavelength dispersion and writes it into real* WPC ***
C       CALL IMGKWR(IM,'cdelt1',WPC,IER)    # Old pre-apall?
       CALL IMGKWR(IM,'cd1_1',WPC,IER)
       IF(IER.NE.0) GOTO 1010
C   Gets the number of pixels in the spectrum, writes to integer NX ***
C       CALL IMGKWI(IM,'np2',NX,IER)        # Old pre-apall?
       CALL IMGKWI(IM,'naxis1',NX,IER)
       IF(IER.NE.0) GOTO 1010
C   Gets the flux as a function of pixel, writes to real# FLX
       CALL IMGL1R(IM,FLX,IER)
       IF(IER.NE.0) GOTO 1010
C   Closes the image
       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
