% file: mxs_abund_pkg.sl
% 
% Stand-alone routines for determining plasma properties from a given
% abundance set.   Reference abundances are either "angr" or "wilm".
% 
% mxs = MIT X-ray Spectroscopy
% 
% --- Summary
%
% Public variables (values are set by the routines)
%
%     Used_numfrac[Z]  - to-be-used array of number fractions of ions Z
%             AofZ[Z]  - array of mass of ion Z in amu.
%     AG89_numfrac[Z]  - "angr" array of number fractions of ions Z
%     Wilm_numfrac[Z]  - "wilm" array of number fractions of ions Z
%           
% Public routines
%
%        mxs_abund_set(abund_str);   Select "angr" or "wilm" abund.s to
%                                     be used (i.e., put into Used_numfrac[Z].)
%                                 
%     (...) = get_par_xabunds(model_str);
%     [...] = calc_QofZ (xzs, TcieK);
%       val = calc_amu_per_ion(xzs, xabunds);
%       val = calc_es_per_ion(xzs, xabunds, TcieK);
%       val = calc_amu_per_particle(xzs, xabunds, TcieK);
%       val = calc_es_per_H(xzs, xabunds, TcieK);
%     [...] = calc_QofZ (xzs, TcieK);
%
%    
% --- Useage and Examples
%
% Load the routines, etc (Assuming this file is on the isis load path):
%     .load mxs_abund_calcs
%     ;
%           or
%     evalfile("mxs_abund_calcs");
%
% Set the desired abundance reference set to use:
%     mxs_abund_set("wilm");
%   
% Get the elements and (relative) abundances from an XSPEC model:
%     variable xzs, xabunds;
%     fit_fun("vgnei(1)");       
%     (xzs, xabunds) = get_par_xabunds ("vgnei(1)");
%
% Public routines to calculate various useful scalar values:
%  
%     val = calc_amu_per_ion(xzs,xabunds);
%     val = calc_es_per_ion(xzs,xabunds,1.e7);
%     val = calc_amu_per_particle(xzs,xabunds,1.e7);
%     val = calc_es_per_H(xzs,xabunds,1.e7);
%
% Get the array of electrons-per-ion for each element given
% the plasma (equivalent CIE) Temperature in K.  This routine
% gives the user access to these values which are calculated and
% used internally for several of the scalar functions above.)
% 
%     TcieK = 1.e7;
%     es_per_ion_array = calc_QofZ (xzs, TcieK);
%
%
% - - -
%
% 29-Jan-2010 version from $SCI/CasA/isis_100120/mxs_abund_calcs.sl
%
% 21-Nov-2012 Updated to include "wilm" abundances and ability to
%             choose between "angr" and "wilm"s, mxs_abund_set().
%
% - - -


% - - -
%
% The Atomic mass of the elements/ions
%
public variable AofZ = Float_Type[93];
% and fill them...
AofZ[1]  =  1.0;  % H
AofZ[2]  =  4.0;  % He
AofZ[3]  =  6.0;  % Li
AofZ[4]  =  8.0;  % Be
AofZ[5]  = 10.0;  % B
AofZ[6]  = 12.0;  % C
AofZ[7]  = 14.0;  % N
AofZ[8]  = 16.0;  % O
AofZ[9]  = 19.0;  % F
AofZ[10] = 20.0;  % Ne
AofZ[11] = 23.0;  % Na
AofZ[12] = 24.0;  % Mg
AofZ[13] = 27.0;  % Al
AofZ[14] = 28.0;  % Si
AofZ[15] = 31.0;  % P
AofZ[16] = 32.0;  % S
AofZ[17] = 35.0;  % Cl
AofZ[18] = 40.0;  % Ar
AofZ[19] = 39.0;  % K
AofZ[20] = 40.0;  % Ca
AofZ[21] = 45.0;  % Sc
AofZ[22] = 48.0;  % Ti
AofZ[23] = 51.0;  % V
AofZ[24] = 52.0;  % Cr
AofZ[25] = 55.0;  % Mn
AofZ[26] = 56.0;  % Fe
AofZ[27] = 59.0;  % Co
AofZ[28] = 59.0;  % Ni
% - - -


% - - -
%
% The "angr" abundances, AG89_numfrac(z) :
%
% Define the AG89 abundances by number up to z=92 :
public variable AG89_numfrac = Float_Type[93];
% and fill it:
AG89_numfrac[1] = 1.0;
AG89_numfrac[2] = 0.097723722;
AG89_numfrac[6] = 0.00036307805;
AG89_numfrac[7] = 0.00011220185;
AG89_numfrac[8] = 0.00085113804;
AG89_numfrac[10] = 0.00012302688;
AG89_numfrac[11] = 2.1379621e-06;
AG89_numfrac[12] = 3.8018940e-05;
AG89_numfrac[13] = 2.9512092e-06;
AG89_numfrac[14] = 3.5481339e-05;
AG89_numfrac[15] = 2.8183829e-07;
AG89_numfrac[16] = 1.6218101e-05;
AG89_numfrac[17] = 3.1622777e-07;
AG89_numfrac[18] = 3.6307805e-06;
AG89_numfrac[19] = 1.3182567e-07;
AG89_numfrac[20] = 2.2908677e-06;
AG89_numfrac[21] = 1.2589254e-09;
AG89_numfrac[22] = 9.7723722e-08;
AG89_numfrac[23] = 1.0000000e-08;
AG89_numfrac[24] = 4.6773514e-07;
AG89_numfrac[25] = 2.4547089e-07;
AG89_numfrac[26] = 4.6773514e-05;
AG89_numfrac[27] = 8.3176377e-08;
AG89_numfrac[28] = 1.7782794e-06;
% NORMALIZE it to 1.0 ion:
AG89_numfrac = AG89_numfrac/sum(AG89_numfrac);
% - - -


% - - -
%
% The "wilm" abundances, Wilm_numfrac(z).
% The initial values below are from column "ISM^d" of Table 2 of
% Wilms, Allen & McCray (2000) 0.0 for untabulated values.
% 
% Define the Wilm abundances by number up to z=92 :
public variable Wilm_numfrac = -10.0+Float_Type[93];
% and fill it:
Wilm_numfrac[1] = 12.0;
Wilm_numfrac[2] = 10.99;
Wilm_numfrac[6] = 8.38;
Wilm_numfrac[7] = 7.88;
Wilm_numfrac[8] = 8.69;
Wilm_numfrac[10] = 7.94;
Wilm_numfrac[11] = 6.16;
Wilm_numfrac[12] = 7.40;
Wilm_numfrac[13] = 6.33;
Wilm_numfrac[14] = 7.27;
Wilm_numfrac[15] = 5.42;
Wilm_numfrac[16] = 7.09;
Wilm_numfrac[17] = 5.12;
Wilm_numfrac[18] = 6.41;
%%Wilm_numfrac[19] = 0.0;
Wilm_numfrac[20] = 6.20;
%%Wilm_numfrac[21] = 0.0;
Wilm_numfrac[22] = 4.81;
%%Wilm_numfrac[23] = 0.0;
Wilm_numfrac[24] = 5.51;
Wilm_numfrac[25] = 5.34;
Wilm_numfrac[26] = 7.43;
Wilm_numfrac[27] = 4.92;
Wilm_numfrac[28] = 6.05;
% Convert to non-log values:
Wilm_numfrac = 10.0^(Wilm_numfrac-12.0);
% zero the unspecified ones:
Wilm_numfrac[where(Wilm_numfrac < 1.e-20)]=0.0;
% NORMALIZE it to 1.0 ion:
Wilm_numfrac = Wilm_numfrac/sum(Wilm_numfrac);
% - - -


% - - -
% Define the 
%     Used_numfrac[Z]  - to-be-used array of number fractions of ions Z
% Default is "angr":
% 
public variable Used_numfrac = Float_Type[93];
Used_numfrac = 1.0*AG89_numfrac;
public variable Used_abund_str="angr";

% Routine to set the desired abundance set
public define mxs_abund_set ( abund_str )
{
   Used_numfrac[0]=-1.0;
   if ((abund_str == "angr") or (abund_str == "AG89"))
     {
	Used_abund_str = "angr";
	Used_numfrac = 1.0*AG89_numfrac;
     }
   if ((abund_str == "wilm") or (abund_str == "Wilm"))
     {
	Used_abund_str = "wilm";
	Used_numfrac = 1.0*Wilm_numfrac;
     }
   % catch invalid inputs:
   if (Used_numfrac[0] < 0.0)
     {
	message(" *** mxs_abund_set: Unrecognized abund string: "+abund_str);
	Used_numfrac[0]=0.0;
     }
}
% - - -


% - - -
% Simple command-line help:
public define mxs_abund_help ()
{
   message("");
   message("    mxs_abund_pkg      version: 26 Nov 2012");
   message("");
   message(" Public arrays:");
   message("   Used_numfrac[Z]  <-- filled with \""+Used_abund_str+"\" values");
   message("   AofZ[Z], AG89_numfrac[Z], Wilm_numfrac[Z]");
   message(" Public routines:");
   message("   mxs_abund_set(\"angr\"); % Set Used_numfrac to angr or wilm ");   
   message("   (xzs,xabunds) = get_par_xabunds(\"vpshock(1)\"); ");
   message("   val = calc_amu_per_ion(xzs, xabunds); ");
   message("   val = calc_es_per_ion(xzs, xabunds, TcieK); ");
   message("   val = calc_amu_per_particle(xzs, xabunds, TcieK); ");
   message("   val = calc_es_per_H(xzs, xabunds, TcieK); ");
   message("   array = calc_QofZ(xzs, TcieK); ");
   message("");
}
% - - -

% - - -
% 
% Routine to get the XSPEC zs and abundances from an isis model component:
%
% Only compile if isis is available:
#ifexists get_par
public define get_par_xabunds ( model_component_string )
{
   % set the nominal XSPEC Zs and their strings
   variable xabund_zs, xabund_strs, xabunds, iz;
   xabund_zs = [1,2,6,7,8,10,12,14,16,18,20,26,28];   
   xabund_strs = ["H","He","C","N","O","Ne","Mg","Si","S","Ar","Ca","Fe","Ni"];
   xabunds = Float_Type [length(xabund_zs)];

   iz=0;
   _for iz (0, length(xabund_zs)-1, 1)
     {
    xabunds[iz] = get_par(model_component_string+"."+xabund_strs[iz]);
     }
      
   return xabund_zs, xabunds;
}
#endif
% - - -


% - - -
%
% The number of electrons per ion, for CIE at temp T in K
%
public define calc_QofZ (zs, TcieK)
{
   variable sel, qfrac = Float_Type[length(zs)];
   % Ideally this would be looked up from tables...
   % Can get a reasonable approximation from a simple
   % z^2-scaled equation:
   %    q(z,T) = z*0.412206*(log10(t*(28.0/z)^2) - 5.17609) ;
   % This equation was determined using the QofZ_exact_CIE.sl
   % script which is included for reference after this definition.
   % Just fill these:
   qfrac = 0.412206*(log10(TcieK*(28.0/zs)^2) - 5.17609) ;
   % clip to 0 - 1 range:
   sel = where(qfrac < 0.0);
   if (length(sel) > 0) qfrac[sel]=0.0;
   sel = where(qfrac > 1.0);
   if (length(sel) > 0) qfrac[sel]=1.0;
   % scale by z and return:
   return zs*qfrac;
}

#iffalse
% The following script -- QofZ_exact_CIE.sl --
% is ignored but left here for reference
% as the source of the approximation in calc_QofZ.
%
% Use apec and calculate the number of free electrons
% vs temp for the elements.

% Initialize the data base:
plasma(aped);

% temp grid:
t = 10.0^[4.0:8.0:0.05];
z = 1;

% setup the plot evaluating the Z=1 curve
qplus1 = 1;
frac = ion_frac(z, qplus1, t);
% put it in a .ps file:
variable PSwin=open_plot("qofz_approx_plot.ps/CPS");
xrange(1.e4,1.e9); xlog; yrange(-0.02,1.02);
label("T * (28/Z)^2  [K]","Q_T(Z) / Z",
      "Average ion charge-state vs Scaled Temperature, Z=1-28");
set_line_width(3);
plot(t,frac,0);  % invisible curve

% Will save the Si curve in this array
si_qoft = 0.0*t;

% loop over zs
iclr=1;
_for z (1,28,1)
{
   
   % reset this
   qoft = 0.0*t;
   
   _for qplus1 (1, z+1, 1)
     {
	frac = ion_frac(z, qplus1, t);
	qoft = qoft + frac*(qplus1 - 1.0);
     }
   
   % show the scaled charge vs scaled temp:
   oplot(t*(28.0/z)^2, qoft/z, iclr);
   color(iclr);
   xylabel(5.e8, z/29.0, string(z));
   iclr = iclr+1; if (iclr == 9) iclr = 1;
   % save the Si:
   if (z == 14) si_qoft = qoft/z;
}

set_line_width(7);
% Show the Si curve on top (use the Si color)
oplot(t*(28.0/14)^2, si_qoft,6);
% Show the simple approximation - a manually adjusted straight line:
oplot([1.5e5,4.e7],[0.0,1.0],1);
% label the Si and approximation
color(1);
xylabel(1.e5,0.9,"Simple approximation (black line)");
% all done with the plot
close_plot(PSwin);
% - - -
% So use this approximation (clipped at 0 and 1 ) :
%
%   Q_T(Z) / Z = 0.412206*(log10(T*(28.0/Z)^2) - 5.17609) ;
%              = 0.412206 * log10( 100 * (T/3.75e6)*(14/Z)^2 )  ;
%   
#endif
% - - -


% - - - 
% Calculating useful plasma scalars:  mu_A, q_A, mu and (n_e/n_H)

% mu_A = average ion mass:
public define calc_amu_per_ion (xzs, xabunds)
{
   return sum(AofZ[xzs]*xabunds*Used_numfrac[xzs])/sum(xabunds*Used_numfrac[xzs]);
}

% q_A = average number of electrons per ion:
public define calc_es_per_ion (xzs, xabunds, TcieK)
{
   return sum(calc_QofZ(xzs,TcieK)*xabunds*Used_numfrac[xzs])/sum(xabunds*Used_numfrac[xzs]);
}

% "mu" = average mass per particle (ions+e's):
public define calc_amu_per_particle(xzs, xabunds, TcieK)
{
   return calc_amu_per_ion(xzs,xabunds)/(1.0+calc_es_per_ion(xzs,xabunds,TcieK));
}

% (n_e/n_H) = average number of electrons per H:
public define calc_es_per_H (xzs, xabunds, TcieK)
{
   variable iH = (where(xzs == 1))[0];
   variable H_per_ion = xabunds[iH]*Used_numfrac[1]/sum(xabunds*Used_numfrac[xzs]);
   return calc_es_per_ion(xzs, xabunds, TcieK)/H_per_ion;
}
% - - -