%; Time-stamp: <2002-03-04 17:03:06 dph> 
%; MIT Directory: ~dph/d1/CXO_Data/AR_Lac/Sl
%; File: get_line_fluxes.sl
%; Author: D. Huenemoerder
%; Original version: 2001.01.05
%;====================================================================
% version: 1.1    2002.03.04 add line names argument to get_line_fluxes
% version: 1
%
% purpose: utilities to compute line fluxes, given a plasma model.
%          Line info is given by
%              Trans_List_define.sl, which defines an assoc array, used via:
%
% () = evalfile("Sl/Trans_List_define.sl");
%
% examples:
%
% trans_list["Fe25HeLa"];
%
% page_group(where(eval(trans_list["Fe25HeLa"])));
%
% g=(where(eval(trans_list["Fe25HeLa"])));
%
% w = (line_info(g[0])).lambda;
%

% we will assume that a model has been defined, so that
% eval_fun() returns the model flux array.
% To get flux in al line, we will :

%  f = sum( eval_fun( xlo, xhi ) );
%  
% which assumes that aped_fit_models.sl has been eval'd, and the model
% is aped_multi or aped_multi_rv.  These models have a flag for
% evaluating line-only models (in isis 0.9.70 or greater).
%

%+                  pre-define names; function defs to follow.
  define wave_mean();
  define find_logtmax();
%-

%+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  define get_line_fluxes( )
  {

 % trans_list is an associative array which maps line mnemonics, such
 % as "O8HLa" or "Fe17w17.01" to a function string which evaluates to
 % a list of aped line indices.
 % One definition of a trans_list is in Trans_List_define.sl

    variable keys, trans_list;

    if ( _NARGS == 2 )
    {
       keys = ();
       trans_list = ();
    }
    else 
    {
      trans_list = ();
      keys = assoc_get_keys( trans_list );
    }
    
    variable nlines = length( keys );

    variable f, w, ii ;

    variable result = struct { flux, wave, name, LogTmax };
    result.flux = Float_Type[nlines];
    result.wave = Float_Type[nlines];
    result.name = String_Type[nlines];
    result.LogTmax = Float_Type[nlines];
  
%  specify a temperature grid for finding LogT_max
%
    variable dw = 0.05, 
             wstep = 0.005,
	     t=10.^[5.5:8.5:0.1];

%
% Set a flag so that binning the plasma model spectrum w/
% model_spectrum() evaluates only the line-component.

    Aped_Filter.type = MODEL_LINES;  

% loop over all features in the trans_list[]
%
    for (ii=0; ii<nlines; ii++)
    {
      variable g = where( eval( trans_list[ keys[ii] ] ) );

      w = wave_mean( g );  % mean wavelength of the bump; if the bump
                           % were a complex blend and broad, this
                           % might be the wrong thing to use for
                           % defining the center of the region bounded
                           % by +-dw.

% LogT_max is tabulated in the aped files, but ISIS <=0.9.71 does not
% retrieve this.  So do it empirically.  This is also safer if the
% bump is actually a blend (though the max may not be meaningful if
% logT_max of components makes the emissivity bi-modal)

      result.LogTmax[ii] = find_logtmax(t,g);   

% define the region of interest.  This is in model-space, so +-0.05A
% is probably plenty.

      variable x = [w-dw : w+dw : wstep];


% set the model_spectrum() filter to the lines of interest:

      Aped_Filter.lines = g;

% Evaluate the model flux of the lines.  This assumes that the
% fit_fun() has been defined.

      f = sum( eval_fun( x, x+wstep ) );

      result.flux[ii] = f;
      result.wave[ii] = w;
      result.name[ii] = keys[ii];

% DEBUG print each line as it is processed...

      ()=printf("%8.3f %6.2f %10.3e %20s\n",
	 result.wave[ii],
	 result.LogTmax[ii],
	 result.flux[ii],
	 result.name[ii]);
    } 


    Aped_Filter.type = NULL;   % reset, so model_spectrum() returns line+continuum.

    return result;

  }
%- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    

% UTILITY functions

%+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  define wave_mean( g )
  {
    variable ng = length(g),
             ii,
	     wmean = 0.0;

    for (ii=0; ii<ng; ii++)
    {
      wmean += ( line_info( g[ii] ) ).lambda;
    }

    wmean /= ng;

    return wmean;

  }
%- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  define find_logtmax( t, g )
  {
    variable eps=Float_Type[length(t)],
             LogTmax;

    eps = line_em( g, t );
    LogTmax =  log10 ( t[ where( eps == max(eps) ) ] );

    if (length(LogTmax) > 1)      LogTmax = LogTmax[0];

    return LogTmax[0];

  }
%- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%