%; Time-stamp: <2002-03-12 14:41:51 dph> 
%; MIT Directory: ~dph/libisis
%; File: apec_contin_fit.sl
%; Author: D. Huenemoerder
%; Original version: 2001.10.19
%;====================================================================
%; version 2.0  2002.03.12; fix regrid -> rebin
% version: 1.0
%%%%%%%%%%%%%%%%%%%%%%%%
%
% Name: aped_contin_fit( xlo, xhi, par )
%
% Purpose: define an apec model which computes the apec continuum
%	   spectrum for 1 temperature component.  Does not
%	   support continuua by ion or for variable abundances
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Name: fast_aped_contin_fit( xlo, xhi, par )
%
% Purpose: Use a global continuum model variable for
%	   a re-normalizable continuum shape model.  If the norm is
%	   frozen, then it is an absolute continuum model.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Name: fast_aped_contin_init(nhist)
%
% Purpose: Read a continuum model (from a fit) and store into a global
%	   variable for use by fast_aped_contin_fit().
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
static variable Contin_X_lo, Contin_X_hi, Contin_Value;
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
define fast_aped_contin_init (nhist)
{  %+ define
   % ASSUMES that a continuum model has been evaluated for histogram
   % nhist.  This should be called with an MEG histogram, in order to
   % provide enough bandpass for both HEG and MEG.
   % Resolution is not very important, since it is a continuum.

% Setup is not done here.  First, the equivalent of the following
% should be done:
%%   fit_fun("aped_contin(1) + aped_contin(2) + ... + aped_contin(n)");
%%   set_par(...);
%%   fit_counts;
%%   % etc.
%% Prepare to grab the final continuum array:
%%   ignore([lo:hi]);
%%   notice(hist);
%%   group_data(nhist,0);     % undo any grouping or binning
%%   rebin_data(nhist,0);
%%   eval_flux;		    % re-evaluate flux, unbinned.

% Here we  get the continuum array and set global variables:
   variable y;
   y=get_model_flux(nhist);
   Contin_Value = y.value;
   Contin_X_lo = y.bin_lo;
   Contin_X_hi = y.bin_hi;
%
}   %- define

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Now we simply regrid the Contin_Value to the specified grid,
%% allowing a re-normalization:
%
define fast_aped_contin_fit (xlo, xhi, par)
{  %+ define
   variable y =
      rebin (xlo, xhi, Contin_X_lo, Contin_X_hi, Contin_Value) * par[0];
   return y;
}   %- define


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define an aped continuum function.
% Proper use of this requires ignoring line features.
%
define aped_contin_fit ( xlo, xhi, par )
{ %+ define
  % xlo, xhi = grid of wavelength bins;
  % par[0] = norm;              10^{-14} * VEM/ (4piD^2)
  %				10^{-14}~ int{dV n_e n_h} / (4piD^2)
  % par[1] = electron temperature [K];

  variable y;
  variable KNORM=1.e14;

  y = (get_contin(xlo, xhi, par[1])).pseudo;
  y = (y + (get_contin(xlo, xhi, par[1])).true) * par[0]*KNORM;
    
  return y;

} %- define

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define an aped multi-T continuum function.
% Proper use of this requires ignoring line features.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
define aped_multi_contin_init ( ncomps )
{
% Dynamically define a continuum model w/ ncomps temperature compents.

  variable i, func_name, T_params, Norm_params, params;

  params = "norm_1";

  for (i=2; i <= ncomps; i++)
    params = [params, sprintf("norm_%d",i)];

  for (i=1; i <= ncomps; i++)
    params = [params, sprintf("Temperature_%d", i)];

  add_slang_function( "aped_multi_contin", params, [1:ncomps] );

}
%-

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+

define aped_multi_contin_setup(Temps, Norms)
{
  variable i;

  aped_multi_contin_init( length(Temps) );

  fit_fun("aped_multi_contin(1)");

  for ( i = 0; i < length(Temps); i++)
  {
    set_par( i+1, Norms[i], 0, 0, 100*Norms[i] );
    set_par( i+length(Temps)+1, Temps[i], 0, 1.e4, 7.95e8 );
  }

  list_par;
}  

%-

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
define aped_multi_contin_fit ( xlo, xhi, par )
{
%
% PURPOSE:  isis user model, APED thermal plasma continuum, multi-temperature,
%	    cosmic abundances, low density limit.

% xlo, xhi = grid of wavelength bins;
% 
  variable ii,                        	     % loop counter
	   y, y_tmp,                         % return value
           ncomps = length(par)/2;    	     % # model components
  variable KNORM=1.e14;

  variable norms=par[ [0:ncomps-1] ];
  variable temps=par[ [ncomps:2*ncomps-1] ];	   	   

  y = Float_Type[ length(xlo) ];
  for (ii = 0; ii < ncomps; ii++ )
  {
    y_tmp = get_contin( xlo, xhi, temps[ii] );
    y = y + (y_tmp.pseudo + y_tmp.true) * norms[ii];
  }
  return y*KNORM;
}  
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



add_slang_function("aped_contin", ["norm","T"]);
add_slang_function ("fast_aped_contin", ["norm"]);