%; Time-stamp: <2002-05-08 11:39:10 dph> 
%; MIT Directory: ~dph/libisis
%; File: aped_fit_models.sl
%; Author: D. Huenemoerder
%; Original version: 2001.10.23
%;====================================================================
%
% purpose:  Define isis user models for aped thermal plasmas.
%
% version 3.1 2002.04.03  added aped_fit - basic model, one component.
%
% version: 3.0 2002.01.06 add global variable and switch to support 
%                         qualifier of model_spectrum, and use for
%                         evaluating a model for specific lines or
%                         contin. 
%
% version: 2.0 2001.11.29 include rv version of model
%                fixed bug in aped_multi_fit, w/ array counting/indexing.
% version: 1.1 2001.10.29 added is_a_norm flags to
%			  add_slang_function() call
% version: 1.0
%
%
   public variable Aped_Filter = struct {type, lines};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
provide("apedx_fit");
provide("aped_multi_fit");
provide("aped_multi_init");
provide("aped_multi_setup");
provide("aped_set_abund");

provide("aped_rv_fit");
provide("aped_multi_rv_fit");
provide("aped_multi_rv_init");
provide("aped_multi_rv_setup");
provide("aped_multi_2rv_setup");

%
%
define aped_multi_fit();       % multple T isis model
define aped_multi_init();      % defines the model for given # components
define aped_multi_setup();     % establishes fit function, sets the
			       %  default model parameters
define aped_set_abund();       % sets/resets local static abundance arrays.
%
% ditto, but w/ radial velocity as a parameter (single for all components.) 
%
define aped_multi_rv_fit();       % multple T isis model
define aped_multi_rv_init();      % defines the model for given # components
define aped_multi_rv_setup();     % establishes fit function, sets the
			       %  default model parameters


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

%
% define local variables for use in storing alternative abundance
% mixes for the plasma models.  These are not used as parameters of a
% fit, but for ad hoc tweaks of abundances. See aped_set_abund();
%
%
static variable Aped_Elem_List,  % an array of Z values; e.g [Fe, Ne] or [26, 10]
		Aped_Abund_List; % an array of factores relative to
				 % cosmic, e.g., [0.1, 2.0];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
define apedx_fit ( xlo, xhi, par )
{
%
% PURPOSE:  isis user model, APED thermal plasma; simple interface to
% the plasma state structure.
%	    
% xlo, xhi = grid of wavelength bins;
% 
% params are 
%
%    norm = 1
%    temperature = 1e+07
%    vturb = 0
%    redshift = 0
%

  variable ii,                        	     % loop counter
	   y;                         	     % return value

  variable mdl = default_plasma_state();  % plasma model definition

% NOTE: density, absorption not currently supported....

    mdl.norm = par[0];     
    mdl.temperature = par[1];
    mdl.vturb = par[2];              % given in km/s?
    mdl.redshift = par[3]/2.99792458e5;   % param is in velocity, km/s.

    if ( __is_initialized( &Aped_Elem_List ) )   % apply different abundances, if set
    {
        mdl.elem = Aped_Elem_List;
        mdl.elem_abund = Aped_Abund_List;
    }

  define_model(mdl);                      % define the model (use list_model; to see).

   switch (Aped_Filter.type)
   {
       case MODEL_LINES:
       if (Aped_Filter.lines == NULL)
         y = model_spectrum (xlo, xhi, MODEL_LINES);
       else
         y = model_spectrum (xlo, xhi, MODEL_LINES, Aped_Filter.lines);
   }
   {
       case MODEL_CONTIN:
       y = model_spectrum (xlo, xhi, MODEL_CONTIN);
   }
   {
       case NULL:
       y = model_spectrum (xlo, xhi);
   }
   
  return y;

}  
%-

variable params = ["norm", "temp", "vturb", "vrad"];
add_slang_function( "apedx", params );




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

% xlo, xhi = grid of wavelength bins;
% 

  variable ii,                        	     % loop counter
	   y,                         	     % return value
           ncomps = length(par)/2,    	     % # model components
	   pstate = default_plasma_state();  % plasma model definition

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

  for (ii = 0; ii < ncomps; ii++ )
  {
    mdl[ii] =             @pstate;
    mdl[ii].temperature = temps[ii];
    mdl[ii].norm =        norms[ii];

    if ( __is_initialized( &Aped_Elem_List ) )   % apply different abundances, if set
    {
        mdl[ii].elem = Aped_Elem_List;
        mdl[ii].elem_abund = Aped_Abund_List;
    }

  }

  define_model(mdl);                      % define the model (use list_model; to see).

   switch (Aped_Filter.type)
   {
       case MODEL_LINES:
       if (Aped_Filter.lines == NULL)
         y = model_spectrum (xlo, xhi, MODEL_LINES);
       else
         y = model_spectrum (xlo, xhi, MODEL_LINES, Aped_Filter.lines);
   }
   {
       case MODEL_CONTIN:
       y = model_spectrum (xlo, xhi, MODEL_CONTIN);
   }
   {
       case NULL:
       y = model_spectrum (xlo, xhi);
   }
   
  return y;

}  
%-

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
define aped_multi_init ( ncomps )
{
% Dynamically define a 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, "norm_"+string(i)];

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

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

}
%-



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

define aped_multi_setup(Temps, Norms)
{
  variable i;

  aped_multi_init( length(Temps) );

  fit_fun("aped_multi(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_set_abund()       % sets local static abundance arrays.
{
  variable  elems;

% example usage:
% aped_set_abund([Fe, Ne], [0.1, 2.0]);   % change abunds of iron and  neon
% aped_set_abund([Fe, Ne]);   % reset abunds of iron and  neon to  cosmic 
% aped_set_abund();           % reset current list to cosmic.

  switch (_NARGS)
  {
    case 0:           % reset list to cosmic (1.0)
    if (__is_initialized (&Aped_Abund_List) )
      Aped_Abund_List = Aped_Abund_List*0 + 1;     % set to cosmic;
  }
  {
    case 1:   % have element list, so set them to cosmic.
    elems=();
    variable i;
    for ( i = 0;  i < length(elems); i++)
      Aped_Abund_List[where(Aped_Elem_List == elems[i])] = 1.;    % reset list to 1
  }
  {
    case 2:
    Aped_Abund_List = ();
    Aped_Elem_List = ();
  }
  {
    vmessage("Bad arguments.  USAGE: aped_set_abund([elems [,abunds]]);");
  }
}
%-



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+  Version w/ radial velocity delta.
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
define aped_multi_rv_fit ( xlo, xhi, par )
{
%
% PURPOSE:  isis user model, APED thermal plasma, multi-temperature,
%	    cosmic abundances, low density limit, w/ velocity shift 

% xlo, xhi = grid of wavelength bins;
% 

  variable ii,                        	     % loop counter
	   y,                         	     % return value
           ncomps = (length(par)-1) / 2,     % # model temperature components
	   pstate = default_plasma_state();  % plasma model definition

  variable rv = par[0];               % radial velocity, km/s
  variable mdl = Struct_Type[ncomps];
  variable norms = par[ [1 : ncomps] ];  % normalizations, 1.e-14*EM/4piD^2
  variable temps=par[ [1+ncomps : 2*ncomps]];  % temperatures, K

  for (ii = 0; ii < ncomps; ii++ )
  {
    mdl[ii] =             @pstate;
    mdl[ii].temperature = temps[ii];
    mdl[ii].norm =        norms[ii];
    mdl[ii].redshift =    rv / 2.99792458e5;

    if ( __is_initialized( &Aped_Elem_List ) )   % apply different abundances, if set
    {
        mdl[ii].elem = Aped_Elem_List;
        mdl[ii].elem_abund = Aped_Abund_List;
    }

  }

  define_model(mdl);                      % define the model (use list_model; to see).

   switch (Aped_Filter.type)
   {
       case MODEL_LINES:
       if (Aped_Filter.lines == NULL)
         y = model_spectrum (xlo, xhi, MODEL_LINES);
       else
         y = model_spectrum (xlo, xhi, MODEL_LINES, Aped_Filter.lines);
   }
   {
       case MODEL_CONTIN:
       y = model_spectrum (xlo, xhi, MODEL_CONTIN);
   }
   {
       case NULL:
       y = model_spectrum (xlo, xhi);
   }

  return y;

}  
%-

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%+
define aped_multi_rv_init ( ncomps )
{
% Dynamically define a model w/ ncomps temperature components.

  variable i, func_name, T_params, Norm_params, params;

  params = "vrad";

  for (i=1; i <= ncomps; i++)
    params = [params, "norm_"+string(i)];

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

  add_slang_function( "aped_multi_rv", params, [2:ncomps+1] );

}
%-



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

define aped_multi_rv_setup(Vrad, Temps, Norms)
{
  variable i;

  aped_multi_rv_init( length(Temps) );

  fit_fun("aped_multi_rv(1)");

  set_par( 1, Vrad, 0, Vrad-600, Vrad+600);

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

%  list_par;
}  

%-


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

define aped_multi_2rv_setup(Vrads, Temps, Norms1, Norms2)
{
  variable i, ncomps;

  ncomps = length(Temps);

  aped_multi_rv_init( ncomps );

  fit_fun("aped_multi_rv(1)+aped_multi_rv(2)");

  set_par( 1, Vrads[0], 0, Vrads[0]-600, Vrads[0]+600);

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


  set_par( ncomps*2+2, Vrads[1], 0, Vrads[1]-600, Vrads[1]+600);

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

%  list_par;
}  

%-