;+
; NAME: EPSILON
;	
; PURPOSE:
;	To calculate the epsilon coefficient for a filter acting 
;       on a modelled spectrum.
; CATEGORY:
;	
; CALLING SEQUENCE:
;       EPSILON,SFILENAME,NCW,TAU,EPS,SPLOT=splot,LIMITS=limits
;
; INPUTS:
;       SFILENAME: A string variable giving the full path and
;                  filename of a synthetic spectrum.
;
;       NCW: A 3 digit integer, giving the center wavlenght of 
;            the filter in angstroms.
;
; KEYWORD PARAMETERS:
;       SPLOT:  If set, create a PostScript plot file showing the
;               results.
;
;       LIMITS: Optionally provides the wavelength limits (nm) over 
;               which the integrations are to be performed. If this
;               keyword is specified, no plotting is done either to the
;               screen or to a PostScript file.
;
; OUTPUTS:
;       EPS: The epsilon coefficient.
;
;       TAU: The peak transmission coefficient.
;
; COMMON BLOCKS:
;       None.
;
; SIDE EFFECTS:
;       None.
;
; RESTRICTIONS:
;       None known.
;
; PROCEDURE:
;       First, the transmission coefficient tau, is extracted from
;       the file ~steele/PAC/filter/filt.trans. Next the modelled spectrum
;       is recovered from the file specified by SFILENAME. A relative 
;       transmission curve corresponding to the center wavelength, NWC, 
;       is then built. The wavelength intervals over which the filter
;       and the spectral curves are defined, are expanded so that they
;       are equal and span a wavelength interval that contains them both.
;       The curves are set to zero where they were not previously defined
;       to accomodate the possible change in the interval over which they
;       are now defined. The mouse is then used to defined the limits of
;       integration that will be used to calculate epsilon.
;
; MODIFICATION HISTORY:
;       Written June 1994, by T A Oliynyk.	
;
;       June 23, 1994: The wavelengths from the modeled spectrum
;                      can now have intervals different from 0.02 
;                      nanometers.
;                      The limits of integration can now be specified
;                      by using the mouse.
;       June 24, 1994: If the wavelength intervals of the filter and
;                      spectral curve are not compatable, then the 
;                      user will be asked to redefine a compatable
;                      interval.
;-

PRO epsilon,sfilename,ncw,tau,eps,splot=splot,limits=limits

IF N_Params() NE 4 THEN BEGIN
    Doc_Library,'epsilon'
    Return
    ENDIF

IF N_Elements(limits) GT 0 THEN BEGIN
    noplot=1
    limits=limits(Sort(limits))
    ENDIF ELSE noplot=0
    
@isitdos

; Get all peak transmission coefficients.
rdcols,filtrans,n,ncws,taus

; Get the modelled spectrum.
check=RFindFile(sfilename,Count=nch)
IF nch NE 1 THEN BEGIN
    Message,'Modelled spectrum not found!',/Informational
    Return
    ENDIF
rdcols,sfilename,n,wl,inty

; select proper peak transmission coefficient.
l=WHERE(ncws EQ ncw)
tau=taus(l(0))
wl=wl/10  ; change units to nanometers.
int=Median(wl-Shift(wl,1))
int=Round(int*100.)/100.

; Get the absolute filter transmission curve.
filtrans,ncw,int,cw,rpb,hw,tw,w1,rm1

tol=int/10.
; Find smallest and largest wavelength of ...
minw1=MIN(w1,MAX=maxw1) ; filter curve, and ...
minwl=MIN(wl,MAX=maxwl) ; spectral curve.

; Ensure the synthetic spectrum and filter transmission are on
; compatible wavelength scales.
i=0
j=0
interval=String(int,Format='(f5.3)')
; Print,Abs(minw1-minwl)/int,Long(Abs(minw1-minwl)/int)
; Stop
WHILE (Abs(minw1-minwl)/int-Long(Abs(minw1-minwl)/int)) GT 0.0001 DO BEGIN
    j=1
    IF i GT 0 THEN Print, 'Not a valid interval choose again!'
    IF i EQ 0 THEN BEGIN 
        Print,'The wavelength intervals for the filter and the spectral'
        PRINT,'curves do not match. Choose a new interval, smaller than'
        PRINT,interval+', such that it will divide ' $
             +String(ABS(minw1-minwl),Format='(f8.3)')+' evenly.'
        ENDIF 
    PRINT,FORMAT='($,"Enter new interval")'
    READ,int
    i=1
    ENDWHILE

; Interpolate the synthetic spectrum and filter transmission onto a
; common wavelength scale, using splines.
IF j EQ 1 THEN BEGIN
    w1n=FINDGEN((1/int)*(maxw1-minw1)+1)*int+minw1
    wln=FINDGEN((1/int)*(maxwl-minwl)+1)*int+minwl
    rm1=Spl_Interp(w1,rm1,Spl_Init(w1,rm1),w1n)
    inty=Spl_Interp(wl,inty,Spl_Init(wl,inty),wln)>0.0
    wl=wln
    w1=w1n
    ; Find smallest and largest wavelengths of ...
    minw1=MIN(w1,MAX=maxw1)     ; the filter transmission, and ...
    minwl=MIN(wl,MAX=maxwl)     ; the synthetic spectrum.
    tol=int/10.
    ENDIF
    
; Expand the wavelength interval over which the filter and the spectral
; curves are defined, to an interval which contains both the original
; intervals. In addition, set the filter and the spectral curves to
; zero for any new wavelengths added by the new interval.
IF ABS(minw1-minwl) GT tol THEN BEGIN
 IF minw1 LT minwl THEN BEGIN
  sind=WHERE(ABS(w1-minwl) LT tol )
  filli=FLTARR(sind(0))
  fillw=w1(0:sind(0)-1) 
  wl=[fillw,wl]
  inty=[filli,inty]
 ENDIF ELSE BEGIN
  sind=WHERE(ABS(wl-minw1) LT tol)
  filli=FLTARR(sind(0))
  fillw=wl(0:sind(0)-1)
  w1=[fillw,w1]
  rm1=[filli,rm1] 
 ENDELSE
ENDIF
IF ABS(maxw1-maxwl) GT tol THEN BEGIN
 IF maxw1 GT maxwl THEN BEGIN
  eind=WHERE(ABS(w1-maxwl) LT tol )
  eind=eind+1
  n1=N_ELEMENTS(w1)
  n1=n1-eind(0)
  fillw=w1(eind(0):*) 
  filli=FLTARR(n1)
  wl=[wl,fillw]
  inty=[inty,filli]
 ENDIF ELSE BEGIN
  eind=WHERE(ABS(wl-maxw1) LT tol )
  eind=eind+1
  nl=N_ELEMENTS(wl)
  nl=nl-eind(0)
  fillw=wl(eind(0):*) 
  filli=FLTARR(nl)
  w1=[w1,fillw]
  rm1=[rm1,filli]
 ENDELSE
ENDIF

maxi=MAX(inty)

t=STRING(tau,FORMAT='(f6.3)')
fn=STRTRIM(ncw,2)
up=MAX(w1,MIN=dn)
rng=up-dn

IF NOT noplot THEN BEGIN
    ; plot the filter and spectral curve over the expanded interval.
    plot,w1,inty/maxi $
        ,CHARSIZE=1.5,XTITLE='Wavelength (nm)' $
        ,YTITLE='Transm. / Rel. Brightness',YRANGE=[0,1.1] $
        ,TITLE='Filter Transmission'
    oplot,w1,rm1,THICK=2.0
    ENDIF

; set up integration interval.
dw=Median(w1-SHIFT(w1,1))
IF N_Elements(limits) GT 0 THEN BEGIN
    r=Where((limits(0) LE w1) AND (w1 LE limits(1)))
;    Print,limits,w1(r(0)),w1(r(N_Elements(r)-1))
    ENDIF ELSE BEGIN
    WShow,0,1
    r=get_roi(w1,rm1,/nohighlight)
    ENDELSE
    
; perform required integrations
; integraten=TOTAL(inty(r)*rm1(r))*dw
integraten=Int_Tabulated(w1(r),inty(r)*rm1(r),/Double)
; integrated=TOTAL(inty(r))*dw
integrated=Int_Tabulated(w1(r),inty(r),/Double)

IF NOT noplot THEN BEGIN
    ; set up vertical bars on plot to show the limits of integration.
    x1=MIN(w1(r),MAX=x2)
    yb=MAX(rm1([r(0),r(N_ELEMENTS(r)-1)]))+0.02
    yt=0.09+yb
    plots,[x1,x1],[yb,yt],THICK=4.0
    plots,[x2,x2],[yb,yt],THICK=4.0
    ENDIF

; calculate epsilon.
eps=integraten/integrated
e=STRING(eps,FORMAT='(f6.4)')

IF (KEYWORD_SET(splot) AND (NOT noplot)) THEN BEGIN
    lw=MIN(w1(r),MAX=uw)
    lw=STRING(lw,FORMAT='(f6.1)')
    uw=STRING(uw,FORMAT='(f6.1)')
    PRINT,'Enter the spectrum type and characteristic temperature.'
    spec='' & temp='' & trial=''
    PRINT, FORMAT='($,"Spectrum: ")'
    READ,spec
    PRINT, FORMAT='($,"Temperature: ")'
    READ,temp
    PRINT, FORMAT='($,"Trial #: ")'
    READ,trial
    psname=tproot+dd+String(ncw,Format='(I3,"c")')+dd $
          +spec+temp+trial+'.ps'
    psopen,psname
    PLOT,w1,inty/maxi,CHARSIZE=1.5,XTITLE='Wavelength (nm)' $
        ,YTITLE='Transm. / Rel. Brightness' $
        ,TITLE='Filter Transmission ' $
        ,YRANGE=[0,1.1]
    OPLOT,w1,rm1,THICK=3.0
    XYOUTS,up-0.05*rng,0.98,spec+' '+temp+' K',ALIGNMENT=1.0
    XYOUTS,up-0.05*rng,0.93,'Peak filter trans.: '+t,ALIGNMENT=1.0
    XYOUTS,up-0.05*rng,0.87,'Epsilon value: '+e,ALIGNMENT=1.0
    XYOUTS,up-0.05*rng,0.81,'Intg. limits: ['+lw+','+uw+']',ALIGNMENT=1.0

    plots,[x1,x1],[yb,yt],THICK=4.0
    plots,[x2,x2],[yb,yt],THICK=4.0

    psclose
    ENDIF

Return
END