% -*- mode: SLang; mode: fold -*-
%; Time-stamp: <2008-08-30 15:06:49 dph>
%; Directory:  ~dph/h3/Analysis/He_like_emis/
%; File:       t_he_modifier.sl
%; Author:     David P. Huenemoerder <dph@space.mit.edu>
%; Orig. version: 2008.08.19
%;========================================
% version:  1.0.0
%          
% purpose: test he_modifier.sl
%          

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
#iffalse

Perform some tests using the He-triplet line emissivity modifier
function provided by he_modifier.sl and it's required table,
He_like_xyz_norm_coef.fits.  These provide density-dependent
emissivity modifiers for use with the APED (1.3.1).

There are three parts to these tests:

1 - Evaluation of the model spectra vs density

2 - Evaluation of line fluxes and ratios vs density

3 - Fitting an observed spectral region with a density dependent model.


#endif
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Set up:

_auto_declare = 1 ; 

plasma(aped);
require( "he_modifier.sl" ) ; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
#iffalse

1 - evaluate spectra for the He-like triplets of C, N, O, Ne, Mg, Si,
    and S for a grid of densities.  We will use isothermal plasma
    models with a temperature near the peak emissivity of each ion.

#endif
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


use_thermal_profile; % so we see can some line width in models


% Define an APED single temperature plasma model called "Aped_1", and
% use it as a fit function, using the emissivity modifier (see
% create_aped_line_modifier and create_aped_fun for details)

create_aped_fun( "Aped_1", default_plasma_state ) ; 
fit_fun( "Aped_1( 1, He_triplets(1) )" );


% Define a function to help automate evaluation of the model in
% spectral regions containing the various triplets:

define he_flux( wmin, wmax, t, n )
{
    % wmin, wmax: min and max wavelengths for the triplet region.
    % t: temperature
    % n: grid of densities

    variable w1, w2, y, i ;

    ( w1, w2 ) = linear_grid( wmin, wmax, 512 ) ;

    variable n_dens = length( n ) ; 
    y = Struct_Type[ n_dens ] ;

    variable z = struct{ bin_lo, bin_hi, value } ; 

    set_par( "Aped_1(1).temperature", t, 1, 0, 0 );

    for (i=0; i<n_dens; i++)
    {
	set_par( "He_triplets(1).density", n[i], 1, 0, 0 );
	y[ i ] = @z ;
	y[ i ].value = eval_fun( w1, w2 ) ;
	y[ i ].bin_lo = @w1 ;
	y[ i ].bin_hi = @w2 ; 
    }
    return ( y ) ; 
}

% Define a density grid over the range of interest

ngrid = 10^[8:16] ; 

pid = open_plot("p_he_triplets_flux.ps/vcps",1,4); resize(16,1.4);

vp = get_outer_viewport;
vp.xmin=0.12; vp.xmax=0.98; vp.ymin=0.12; vp.ymax=0.88;
set_outer_viewport(vp);

xlabel("Wavelength");
ylabel(  latex2pg( "Flux [photons/cm^2/s/A]" ) );
xlin; xrange;
%ylin;  yrange(0);
ylog;  yrange;
charsize( 1.4 ) ; 
plot_bin_density;
color(1);
set_line_width( 5 ) ; 

title("C V" ) ;  
f = he_flux( 40.7, 41.5, 10^6.0, ngrid ) ; 
hplot( f[0] ) ;
array_map( Void_Type, &ohplot, f[[1:]] );

title( "N VI" );
f = he_flux( 29.05, 29.58, 10^6.2, ngrid ) ; 
hplot( f[0] ) ;
array_map( Void_Type, &ohplot, f[[1:]] );

title(  "O VII" ) ;
f = he_flux( 21.76, 22.15, 10^6.3, ngrid ) ; 
hplot( f[0] ) ;
array_map( Void_Type, &ohplot, f[[1:]] );

title ( "Ne IX" ) ; 
f  = he_flux( 13.54, 13.72, 10^6.6, ngrid ) ; 
hplot( f[0] ) ;
array_map( Void_Type, &ohplot, f[[1:]] );

title(  "Mg XI" ) ; 
f = he_flux( 9.22, 9.33, 10^6.8, ngrid ) ; 
hplot( f[0] ) ;
array_map( Void_Type, &ohplot, f[[1:]] );

title ( "Si XIII" );
f = he_flux( 6.68, 6.75, 10^6.7, ngrid ) ; 
hplot( f[0] ) ;
array_map( Void_Type, &ohplot, f[[1:]] );

title ( "S XV" ) ;  
f = he_flux( 5.055, 5.11, 10^7.2, ngrid ) ; 
hplot( f[0] ) ;
array_map( Void_Type, &ohplot, f[[1:]] );

charsize( 1.0); set_line_width(1);
color(1);
xylabel( 5.06, 0.2, "Note: DB limitation: little or no change in input data for i-line");

close_plot(pid); window(1);

ylin; xlin; xrange; yrange;
plot_bin_integral;

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
#iffalse

2 - Evaluate the line fluxes and ratios vs density.  Plot the f/i
    ratio, which is the primarily density-sensitive ratio.
    We will again use a temperature near the peak emissivity for each
    ion. 

#endif
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% Define some functions to facilitate repetive operations:
%
define rif_flux1( id, t, n )
{
    % id = APED line id (scalar)
    % t = temperature (scalar)
    % n = density array

    variable n_dens = length( n ) ; 
    variable f = Double_Type[ n_dens ] ; 

    set_par( "Aped_1(1).temperature", t, 1, 0, 0 );

    variable line_inf = line_info( id ) ;
    variable w = line_inf.lambda ; 

    variable i ; 
    for (i=0; i<n_dens; i++)
    {
	set_par( "He_triplets(1).density", n[ i ], 1, 0, 0 );
	() = eval_fun( w-0.05, w+0.05 );
	f[ i ] = line_info( id ).flux ; 
    }
    return ( f ) ; 
}


define rif_flux( atom, t, n )
{
    % atom = element (the atomic number, or a string element name)
    % t    = temperature
    % n    = grid of densities

    if (typeof( atom ) == String_Type ) atom = eval( atom ) ; 
    variable id_r, id_i, id_f ;  % APED line ids.

    id_r = where( trans( atom, atom-1,     7, 1 ) ) ;  % w:   resonance line
    id_i = where( trans( atom, atom-1, [6,5], 1 ) ) ;  % x,y: intercombination
    id_f = where( trans( atom, atom-1,     2, 1 ) ) ;  % z:   forbidden

    variable fr, fi, ff ;

    use_delta_profile;  % for faster evaluation

    fr = rif_flux1( id_r[0], t, n ) ;
    ff = rif_flux1( id_f[0], t, n ) ;
    fi = rif_flux1( id_i[0], t, n ) + rif_flux1( id_i[1], t, n ) ; 

    variable g, r, l ;

    g = (ff + fi) / fr ;  %  "G" ratio, = (f+i)/r
    r = ff / fi ;         %  "R" ratio, = f/i

    g[where(g != g)] = 0.0 ;      % remove /0 - bad for pgplot
    r[where(fi != fi)] = 0.0 ;

    % package into struct and return:

    variable z = struct{ a, t, n, fr, ff, fi, g, r } ;
    z.a = atom ;
    z.t = t ;
    z.n = @n ;
    z.fr = @fr ;
    z.fi = @fi ;
    z.ff = @ff ;
    z.g = @g ;
    z.r = @r ;

    return( z ) ; 
}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Define a density grid and evaluate the line fluxes and ratios vs
% density for each ion.  Overplot all the f/i ratios vs density.

nn = 10^[ 7.0 : 15.0 : 0.2] ;

fc  = rif_flux( C,  10^6.0, nn);
fn  = rif_flux( N,  10^6.2, nn);
fo  = rif_flux( O,  10^6.3, nn);
fne = rif_flux( Ne, 10^6.6, nn);
fmg = rif_flux( Mg, 10^6.8, nn);
fsi = rif_flux( Si, 10^7.0, nn);
fs  = rif_flux( S,  10^7.2, nn);

pid = open_plot("p_he_triplets_fi.ps/vcps"); resize(16,0.6);

vp = get_outer_viewport;
vp.xmin=0.12; vp.xmax=0.98; vp.ymin=0.12; vp.ymax=0.88;
set_outer_viewport(vp);


xrange( 1.e7, 1.e15);  xlog;
yrange( 0.8, 13 ) ;    ylog;

xlabel("Density"); 
ylabel("R = f/i"); 
title("");

charsize(1.4);
set_line_width(7);

color(1);
plot(  nn, fc.r ) ; 
oplot( nn, fn.r ) ; 
oplot( nn, fo.r ) ; 
oplot( nn, fne.r ) ; 
oplot( nn, fmg.r ) ; 
oplot( nn, fsi.r ) ; 
oplot( nn, fs.r,8 ) ; 

set_line_width( 5 ) ; 
color(1); xylabel( 1.5e7,  10, "C");
color(2); xylabel( 2.e7,  5.5, "N");
color(3); xylabel( 3.e7,  4.2, "O");
color(4); xylabel( 5.e7,  3.3, "Ne");
color(5); xylabel( 5.e11, 3.1, "Mg");
color(6); xylabel( 1.e13, 2.0, "Si");
color(8); xylabel( 1.e14, 1.5, "S");

close_plot(pid); window(1);
xlin; ylin;

%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
#iffalse

3 - Apply the model to an observed HETGS spectrum.

Test the density-dependent He-triplet models on the observed spectrum
of TW Hya.  We will use OBSID 7435, a 150 ks HETGS observation, which
is long enough for a detection of the O VII forbidden line.  For
simplicity, we will only use the MEG -1 order from a single
observation.  The O VII region is very clean and we will not have to
consider blends from other ions, such as in the Ne IX region.

We will fit the density-dependent plasma model to the triplet region
and evaluate the confidence limits for the density.  We will then
evaluate and save the models for the best fit and the 90% confidence
limits and then overplot them.

#endif


% point to the data; (default to ./Data)

pdir = "./Data" ; 
% pdir = getenv("HOME")+"/d1/CXO_Data/TGCat_Test_Data/obs_7435/P_01"$ ; 


% Load the data; Load the responses and assign them; list the loaded set:

load_data( "$pdir/pha2"$, 9 ) ;      % MEG -1 counts spectrum
load_arf(  "$pdir/meg_-1.arf"$ ) ;   % MEG -1 ARF (effective area)
load_rmf(  "$pdir/meg_-1.rmf"$ ) ;   % MEG -1 RMF (line profile and area factor)

assign_rsp( 1, 1, 1 ) ; 

list_data;
list_arf;
list_rmf;

% Define an isothermal plasma model with Solar photosperic abundances;
% specify it as the fit function and set the parameters to some
% initial guesses.  (Note that the global density parameter of the
% plasma state has no effect for APED 1.3.1.)

p             = default_plasma_state;
p.norm        = 0.001 ; 
p.temperature = 10^6.6;

create_aped_fun( "Aped_1", p ) ; 

fit_fun( "Aped_1( 1, He_triplets(1) )" );

set_par( "He_triplets(1).density", 1.e10,  0, 1.e8,   1.e14 );
set_par( "Aped_1(1).norm",         0.001,  0, 1.e-6,  0.1 );
set_par( "Aped_1(1).temperature",  10^6.6, 1, 10^6.2, 10^6.9 ); %frozen
set_par( "Aped_1(1).redshift",     -2.e-4, 0, -0.001, 0.001 );


% To speed model evaluation, make sure we are not computing the
% thermal plus turbulent profile:

use_delta_profile;


% Bin the data a bit, and set the notice ranges to cover the O VII
% triplet and a little continuum around each line:

group_data( 1, 2 );
xnotice( 1, 21.5, 22.2 );
ignore( 1, 21.68, 21.75 ) ; 
ignore( 1, 21.88, 22.04 ) ; 


% Adjust the normalization term, only and inspect:

renorm_counts;

color(1);
ylin; xlin;
xrange( 21.5, 22.2 ) ;
yrange(0);
set_line_width(3);
charsize( 1.3 ) ; 

rplot_counts(1);


% Now do a fit with norm, redshift, and density free.  Use the subplex
% method for better results.  Since there are <10 counts/bin in the i
% and f lines, use the Cash statistic.  We'll leave temperature frozen
% for the moment - we know roughly what it should be.

set_fit_method("subplex");
set_fit_statistic("cash");

fit_counts;
rplot_counts(1);

% The redshift is not an interesting parameter for these purposes, so
% we can freeze it now for faster minimization:

freeze( "Aped_1(1).redshift" ) ; 


% Evaluate the 90% confidence limits on the density and then use the
% limits to evaluate the model counts for the high and low limits.
% Also evaluate the model counts for the best fit and for a density
% well below the critical density for O VII, to see the low-density
% limit on the lines.

(nlo, nhi) = conf( "He_triplets(1).density"); 

message("%% nlo=$nlo  nhi=$nhi"$);
%% nlo=1.67967e+11  nhi=1.58274e+12;   Log:    11.2, 12.2

nfit = get_par( "He_triplets(1).density" );  % get a copy of the best fit param


xnotice(1, 21.5, 22.2 ) ; % notice all bins in region for nicer plots

% Get 90% high-limit spectrum:

set_par( "He_triplets(1).density", nhi ); 
eval_counts;
yhi = get_model_counts( 1 ) ; 


% Get 90% low-limit spectrum:

set_par( "He_triplets(1).density", nlo ); 
eval_counts;
ylo = get_model_counts(1);


% Get low-density limit spectrum:

set_par( "He_triplets(1).density", 1.e8 ) ; 
eval_counts;
ymin = get_model_counts(1);


% Restore best fit and plot:

set_par( "He_triplets(1).density", nfit ) ; 
eval_counts;
rplot_counts(1);


% Plot the best fit model and overplot the additional models:

pid = open_plot( "o7_models.ps/vcps" ); resize( 16, 0.6);

xlabel("Wavelength");
ylabel("Model counts/bin");
title("O VII model comparisons");
set_line_width(7);
plot_model_counts( 1, 1 ) ;  % black
set_line_width(4);
ohplot( ylo, 2 ) ;           % blue
ohplot( yhi, 2 ) ;           % green
set_line_width(1);
ohplot( ymin, 15 ) ;         % gray
oplot_model_counts( 1, 1 ) ;  % black

close_plot( pid ) ; window( 1 ) ; 


% Also try fitting the temperature
% ( the G-ratio is temperature sensitive: G = (f+i)/r )
% and compute the confidence limits:

thaw( 3 ) ;

xnotice( 1, 21.5, 22.2 );
ignore( 1, 21.68, 21.75 ) ; 
ignore( 1, 21.88, 22.04 ) ; 
fit_counts;
rplot_counts(1);

log10(get_par(3));  % 6.65  vs our ad hoc model, 6.6
log10(get_par(1));  % 11.8  vs prev fit of 11.7

conf(3);   %  3.2e6, 5.3e6        6.5, 6.7
conf(1);   %  1.8e11, 2.6e12      11.2, 12.4

save_par( "o7.par" ) ; 

pid = open_plot( "o7_fit.ps/vcps"); resize(16,0.6);

xnotice(1, 21.5, 22.2 ) ; 
xrange( 21.5, 22.2 ) ; 
eval_counts;
rplot_counts(1);

close_plot(pid);  window(1);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
#iffalse

For a more detailed look at the T - n_e solution, we can compute
confidence maps.  We will define linear grids in temperature and
density and compute the confidence contours of the statistic.  (If you
want logarithmic grids, see the help for conf_map_counts).  This can
be quite time consuming, so it may take some experimentation with
coarse grids to determine grid ranges before adopting the final grid.

#endif
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Reset the noticed bins here:

exclude( all_data ) ; 
include( 1 ) ; 
group_data( 1, 4 );
xnotice( 1, 21.50, 22.2 );
ignore(  1, 21.68, 21.75 ) ; 
ignore(  1, 21.88, 22.04 ) ; 

load_par( "o7.par" ) ;   % reload the best fit model


% Define the grids in the temperature and density parameters:
% For subplex with 40x40 the time is about 1200 s on a 2.6 GHz AMD
% Athlon.   Since it takes a while to compute, we will save the
% confidence map to a FITS file.

gT = conf_grid( "Aped_1(1).temperature",  2.0e6,  6.0e6, 32 );
gN = conf_grid( "He_triplets(1).density", 1.0e8,  4.e12, 32 ) ; 

set_fit_method( "subplex" ) ;
set_fit_statistic( "cash" ) ; 

tic; cTN = conf_map_counts( gT, gN ) ;  toc;  

save_conf( cTN, "confmap_o7.fits");


pid = open_plot("confmap_o7.ps/vcps");  resize(16, 1.0);

xrange;
yrange(0);
charsize( 1.1 ) ; 
xlabel( "T" ) ; 
ylabel( "n_e" ) ; 
title ( "O VII" ) ; 
set_line_width( 7 );
plot_conf( cTN ) ; 

close_plot(pid);