%; Time-stamp: <2007-05-22 16:08:00 dph>
%; Directory:  ~dph/h3/CXC/ISIS/tgplot_demo/
%; File:       isis_demo_line_ratios.sl -> demo_02.sl
%; Author:     David P. Huenemoerder <dph@space.mit.edu>
%; Orig. version: 2007.05.11
%;========================================
% version:  1
%          
% purpose: isis demo
%          
%          
%
% Lying Line Ratios
%
% 
% 1 - Define a multi-thermal  plasma model 
% 2 - Load HETGS ARFs and RMFs for first orders,  set exposure, assign responses.
% 3 - Make fake data
% 4 - Define fitting function for O7 He-L\beta and O7 H-L\alpha
%     Fit each line, compute 1\sigma confidence
% 5 - Evaluate theoretical isothermal ratio of O7 HeLb/HLa
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  shorten some function names:
%
alias( "plot_data_counts",   "pdc") ;
alias( "oplot_data_counts",  "opdc");
alias( "plot_model_counts",  "pmc") ;
alias( "oplot_model_counts", "opmc");

% for more pleasing display:
%
putenv("PGPLOT_BACKGROUND=white");
putenv("PGPLOT_FOREGROUND=black");



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 0 - create some windows for plots.
%
variable p1, p2, p3 ; 

p3 = open_plot( "/xwin" ) ;        resize( 18,1.4) ; erase ; 
p2 = open_plot( "/xwin", 1, 2 ) ;  resize( 17,1.4) ; erase ;
p1 = open_plot( "/xwin", 1, 3 ) ;  resize( 16,1.4) ; erase ;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 1 - Define a plasma model

message( "\nDefine a multi-thermal APED plasma model............\n" ) ; 

plasma( aped ) ;   % load the database

variable logtgrid = [ 5.8 : 8.3 : 0.1 ] ; 
variable tgrid = 10^logtgrid ; 
variable norms ; 

% make a 2-peak emd, peaks at 6.5, 7.2, and a hot tail,,,
%
norms = eval_fun2( "gauss", logtgrid, logtgrid+0.1, [0.6, 6.3, 0.15 ] )
      + eval_fun2( "gauss", logtgrid, logtgrid+0.1, [1.3, 7.0, 0.3 ] )
      + eval_fun2( "gauss", logtgrid, logtgrid+0.1, [0.1, 7.7, 0.1 ] );

norms /= sum( norms ) ; 

variable VEM = 1.e53 ; 
variable dpc = 30.0 ; 

norms *=  1.e-14 * VEM / ( dpc^2 * 1.1965e38  ) ; 

variable Plas = default_plasma_state ; 

Plas.norm = norms;
Plas.temperature = tgrid ; 

% alter the Ne and Fe abundances:

Plas.elem = [ Ne, Fe ] ;  
Plas.elem_abund = [ 1.5, 0.5 ] ; 

create_aped_fun( "plas", Plas ) ; 

fit_fun("plas(1)");

window(p1);

define p1a()
{
    charsize(1.33);
    xrange;
    yrange(0);
    xlabel( "Log T" );
    ylabel( "Plasma Norm" );
    set_line_width( 5 ) ;
    title("Emission measure normalizations for a plasma model");
    plot( logtgrid, norms ) ; 
}

p1a;
plot_pause;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 2 - Define and load ARFs, RMFs

message( "\nLoad responses and assign to fake data...............\n") ; 

variable farf = "data_01/evt2"
              + ["HEG_-1", "HEG_1", "MEG_-1", "MEG_1"]
              + "_garf.fits" ;

variable frmf = "data_01/HA_"
              + ["m1_h", "p1_h", "m1_m", "p1_m"]
              + ".rmf" ;

variable a = array_map( Integer_Type, &load_arf, farf ) ; 
variable r = array_map( Integer_Type, &load_rmf, frmf ) ; 

%
% pick an exposure time and set the ARF exposure (the one that matters):
%
variable t_exp = 100.e3 ; 
array_map( Void_Type, &set_arf_exposure, a, t_exp );


define p1b()
{
    xrange(1.5, 25 ) ; 
    yrange(0);
    xlabel("Wavelength");
    ylabel(  latex2pg( "Effective Area [cm^2]" ));
    title( "HETG Arfs (MEG: +1 black, -1 red;  HEG: -1 green, +1 blue)" ) ; 

    hplot(get_arf(a[3]));
    ohplot(get_arf(a[2]));
    ohplot(get_arf(a[0]));
    ohplot(get_arf(a[1]));
}

p1b;
plot_pause;


%
% Define data histograms for fake data (preserve any existing indices):
%
variable h = [1:4] ; 
if ( all_data != NULL )
{
    h = h + max( all_data ) ; 
}

%
% assign responses to empty histograms to use for fake data:
%
array_map( Void_Type, &assign_rsp, a, r, h );

%
% set the data exposure 
%
array_map( Void_Type, &set_data_exposure, a, t_exp );

%
% list everything:
%
 message( "\n List of data, arfs, rmfs................\n" ) ; 
%
list_data;
list_arf;
list_rmf;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 3 - make fake data:
%
% Notice all the data bins, 
%  set grouping to 1 (moot if 1st time),
%  and generate simulated counts:
%

message("\n Generate simulated counts.......................\n");

group_data( h, 1 ) ; 
xnotice( h ) ;
fakeit;

variable slbl = line_label_default_style ; 


define p1c()
{
    set_line_width( 5 ) ; 
    title("Simulated thermal spectrum, HEG, MEG 1st orders") ; 

    pdc(  h[2] ) ; 
    opdc( h[3] ) ; 
    opdc( h[0] ) ; 
    opdc( h[1] ) ; 

    % label a few lines:
    %
    slbl.char_height = 1.5 ; 
    plot_group( brightest( 10, where( wl( 1.5, 10 ) ) ), 1, slbl ) ; 
    plot_group( brightest( 10, where( wl( 12, 25 ) ) ), 1, slbl ) ; 

    message("\nThe bright lines labled are... \n");
    page_group( brightest( 10, where( wl( 1.5, 10 ) ) ) ) ; 
    page_group( brightest( 10, where( wl( 12, 25 ) ) ) ) ; 
}

p1c;
plot_pause ; 

% replot to a file:
%
variable pid = open_plot("demo_02_f1.ps/vcps", 1, 3);
resize( 6*2.54, 1.4 );
p1a; p1b; p1c;
close_plot(pid);
%

% save some line indices for labeling later:
%
variable lox = where( trans(O,7,13,1) or trans( O, 8, 4, 1 ) ) ; 
slbl.label_type  = 1 ; 
slbl.angle = 45 ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 4 - Define fitting function for O7 and O8 lines:
%     gaussian + poly (folded through response)
%

message("\n Fit the O VII 18.6A and O VIII 19A lines.");
message("Compare their ratio to the theoretical one.");
message("Evaluate theoretical temperature dependence of the ratio.");

%
%  plot the region to be fit, and the plasma model
%
window( p2 ) ; 

define p2a()
{
    set_line_width( 5 ) ; 
    title("Plasma Model: fake data, model, residuals");
    yrange(0); 
    xrange( 18.55, 19.05 );
    rplot_counts( h[2] ) ; 

    plot_group( lox, 14, slbl ) ; 
}

p2a;
plot_pause; 

pid = open_plot("demo_02_f2.ps/vcps", 1, 2);
resize( 6*2.54, 1.4 );
p2a;

window(p2);

message("\nFit the oxygen lines................................... \n");

% Fit the O7HeLb 18.6 line:
%

set_fit_method( "subplex" ) ; 

fit_fun("gauss(1) + poly(1)");

set_par("gauss(1).area",    1.e-4, 0, 0.0,    0.01  ) ;
set_par("gauss(1).center", 18.627, 0, 18.61,  18.64 ) ; 
set_par("gauss(1).sigma",   0.001, 1, 0.0001, 0.02  ) ; 

set_par("poly(1).a0", 3.e-4, 0, 0.0, 1.e-2);
set_par("poly(1).a1", 0, 1 ) ; 
set_par("poly(1).a2", 0, 1 ) ; 

ignore(all_data); 
xnotice( h[2], 18.55, 18.9 ) ; 
() = eval_counts;
() = fit_counts( ) ; 

%
message("compute confidence in flux........");
%
variable c0, c1, c2 ; 
c0  = get_par( "gauss(1).area" ) ; 
(,) = conf( "gauss(1).area", 1 ) ; 

%
% re-fit, since conf() often finds a better fit:
%
() = fit_counts( ) ; 
c0  = get_par( "gauss(1).area" ) ; 
(c1, c2) = conf( "gauss(1).area", 0 ) ; 

message("\n O VII 18.6 fit.................\n");
list_par;


% 
message("\n Fit the O8HLa 19A line; since it is a doublet, set sigma free\n");
%
fit_fun("gauss(2)+poly(1)");

set_par("gauss(2).area", 1.e-3, 0, 0.0, 0.01);
set_par("gauss(2).center", 18.967, 0, 18.94, 18.99 ) ; 
set_par("gauss(2).sigma", 0.004, 0, 0.0001, 0.02 ) ;  % doublet

freeze( "poly(1).a0" ) ;   % use same as prev

xnotice( h[2], 18.85, 19.05 ) ; 

() = fit_counts;

(,) = conf( "gauss(2).sigma"); 

% and re-fit after first conf():
%
() = fit_counts;

variable c0a, c1a, c2a ; 
c0a = get_par( "gauss(2).area" ) ; 
(c1a, c2a) = conf( "gauss(2).area", 1 ) ; 

message("\nO VIII 19A fit................... \n");

list_par;


%
% Evaluate full model over region:
%

fit_fun( "gauss(1) + gauss(2) + poly(1)" ) ; 

xnotice( h[2], 18.55, 19.05);
() = eval_counts;

%
% plot gaussian fit, model, and residuals:
%
define p2b()
{
    title("Parametric Model: data, model, residuals");
    rplot_counts( h[2] ) ; 
}

p2b;
plot_pause ; 

% replot to file:
window(pid);
p2b; 
close_plot(pid);
%





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 5- evaluate line ratios; observed, and theoretical 
%    (plasma model and isothermal) ; 

message("\nEvaluate ratios, observed and theoretical.\n");
message("\nCompute theoretical isothermal temperature from the ratio.......... \n");

%
% "observed" ratio, from the fit:
%
variable r = c0 / c0a ; 
variable dr = r * sqrt ( ((c2-c1)/2./c0)^2 + ( (c2a-c1a)/c0a)^2 ) ; 

% evaluate the theoretical isothermal ratio:
%
variable l1, l2, rth ; 

l1 = where( trans( O, 7, 13, 1 ) ) ;      %  13->1 is HeLb in APED notation.
l2 = where( trans( O, 8, [4,3], 1 ) ) ;   %  [4,3]->1 is HLa doublet APED notation.
rth = ratio_em( l1, l2, tgrid ) ; 

% hack 1st two bad points...
rth[0] = rth[2];
rth[1] = rth[2];

% for better interpolation of temperature, make finer grid:
variable xx, yy ; 
xx = [min(logtgrid) : max(logtgrid) : 0.01 ] ; 
yy = interpol( xx, logtgrid, rth ) ; 

%
% also evaluate model feature fluxes:
%
variable fHeLb, fHLa; 

% need to re-establish the plasma model and evaluate it:
% fit_fun("plas(1)"); () = eval_counts;

fHeLb = line_info( l1[0] ).flux ; 
fHLa  = line_info( l2[0] ).flux + line_info( l2[1] ).flux ; 

% evaluate line emissivities:
%
variable eHeLb, eHLa ; 

eHeLb = line_em( l1, tgrid ) ; 
eHLa =  line_em( l2, tgrid ) ; 

% locate the temperature limits implied by the line ratio
% confidence interval:
%
variable ltmin, ltmax ; 
ltmin = logtgrid[ where( rth <= ( r+dr )  ) ] [ 0 ] ; 
ltmax = logtgrid[ where( rth <  ( r-dr )  ) ] [ 0 ] ; 

ltmin = xx[ where( yy <= ( r+dr )  ) ] [ 0 ] ; 
ltmax = xx[ where( yy <  ( r-dr )  ) ] [ 0 ] ; 


window( p3 );

define p3a()
{
    multiplot([1,1,1,1] ) ;
    resize(7*2.54, 1.4);
    set_line_width(5);
    charsize(1.0);

    title( sprintf("f(O7)/f(O8) = %5.3f (%5.3f)", r, dr )) ; 

    yrange( 1.e-3, 1 ) ;  ylog ; 
    xrange( 6.1, 7.8 ) ; xlin;
    xlabel( "Log T" ) ; 
    ylabel( "f(O7)/f(O8) [theory]");
    %plot( logtgrid, rth ) ;
    plot( xx, yy ) ;
    color(1);
    xylabel( 6.8, 0.3, "T-dependent line ratio");

    plot_pause ; 


    oplot( [1,1]*ltmin, [1.e-3, 1], 15 ) ; 
    oplot( [1,1]*ltmax, [1.e-3, 1], 15 ) ; 

    color(14);
    xylabel( 6.45, 1.5, "T(r)~6.48", 0, 0.5 ) ; 

    plot_pause ; 


    ylabel("Relative Emissivity");
    ylin; yrange(0, 1.02);
    plot( logtgrid, eHLa/max(eHLa), 1 ) ; 
    oplot( logtgrid, eHeLb/max(eHeLb), 2 ) ; 
    color(1);
    xylabel( 6.7, 0.9, latex2pg( "O VIII HL\\alpha" ) );
    color(2);
    xylabel( 6.65, 0.2, latex2pg( "O VII HeL\\beta" ) );

    oplot( [1,1]*ltmin, [0,1],  15 ) ; 
    oplot( [1,1]*ltmax, [0, 1], 15 ) ; 

    plot_pause ; 


    ylabel("Rel. DEM");
    ylin; yrange(0);
    plot( logtgrid, norms, 1);
    color( 14 ) ; xylabel( 7.2, 6.e-4, "Emission Measure Dist.");
    oplot( [1,1]*ltmin, [0,1]*8.e-4,  15 ) ; 
    oplot( [1,1]*ltmax, [0, 1]*8.e-4, 15 ) ; 


    plot_pause ; 

    ylabel("EM Weighted Emissivity");
    ylin; yrange(0);
    variable k1;
    k1 = norms*eHLa ;
    plot( logtgrid, k1/max(k1), 1 ) ; 
    k1 = norms*eHeLb ;
    oplot( logtgrid, k1/max(k1), 2 ) ; 
    color(1); xylabel( 6.6, 0.9, latex2pg("O VII HL\\alpha"));
    color(2); xylabel( 6.6, 0.3, latex2pg("O VII HeL\\beta"));
    color(14); xylabel( 7, 0.6, "Flux = \int(curve)" ) ; 

    oplot( [1,1]*ltmin, [0,1],  15 ) ; 
    oplot( [1,1]*ltmax, [0, 1], 15 ) ; 

    plot_pause ; 
}

p3a;

% replot to file:
pid = open_plot("demo_02_f3.ps/vcps");
resize( 6*2.54, 1.4 );
p3a;
close_plot(pid);
%