%; Time-stamp: <2002-03-04 17:07:50 dph> 
%; File: Run_model.sl
%; Author: D. Huenemoerder
%; Original version: 2002.03.04
%;====================================================================
% version: 1
%
% purpose:
%   example simulations, via fakeit
%   evaluates multi-thermal aped plasma on heg and meg grids.
%   does for unabsorbed and absorbed.


ignore(all_data);
xnotice([1,2],1.7,17);
xnotice([3,4],1.7,25);

Define_model;          % defines unabsorbed APED multi-thermal plasma model.
                       % "list_par;"  will show the model.

vmessage("...............Creating fake data, unabsorbed case...................\n");
fakeit;                % evaluate over noticed bins, and add poisson noise

plot_data_counts( 3, red ); oplot_model_counts( 3, green );

% get the simulated spectrum and the model into variables, to save
% them.
%
cm_a = get_data_counts(3).value + get_data_counts(4).value;
ch_a = get_data_counts(1).value + get_data_counts(2).value;
mm_a = get_model_counts(3).value + get_model_counts(4).value;
mh_a = get_model_counts(1).value + get_model_counts(2).value;
fm_a = get_model_flux(3).value;


% add an absorption component from the xspec models:
% "aped_multi" is the name of the plasma model, and can be seen via
% "list_par;"
%
import("xspec");
fit_fun( "wabs(1) * aped_multi(1)" );
set_par(1,1.0);    % set n_H to 1.0e+22

vmessage("...............Creating fake data, absorbed case...................\n");
fakeit;            % generate model and simulated spectrum.

list_data;
list_arf;
list_rmf;

% get the data into arrays.
%
cm_b = get_data_counts(3).value + get_data_counts(4).value;
ch_b = get_data_counts(1).value + get_data_counts(2).value;
mm_b = get_model_counts(3).value + get_model_counts(4).value;
mh_b = get_model_counts(1).value + get_model_counts(2).value;
fm_b = get_model_flux(3).value;

%%%%%%%%%%%%%%%%%%%%%%%%%% PLOTS

%   next line isn't necessary, since it is default.  But is example of
%   how to specify (could use kev, for exampe).
%
window(p1);

label("Wavelength [\\A]",
  "Flux [counts/0.005A]",
  "                                 Example 100ks fakeit's: MEG+-1, DEM");
plot_unit("angstrom");
custom_plot( 1.5, 25.5, xm1, xm2, cm_a, cm_b, fm_a, fm_b, mm_a, mm_b, T_exp);

pid = dup_plot( "Example_spec_100ks_1.ps/cps",  p1);
resize(10*2.54,0.6);


custom_plot( 1.5, 25.5, xm1, xm2, cm_a, cm_b, fm_a, fm_b, mm_a, mm_b, T_exp);
custom_plot( 1.5, 6, xm1, xm2, cm_a, cm_b, fm_a, fm_b, mm_a, mm_b, T_exp);
custom_plot( 6, 7,      xm1, xm2, cm_a, cm_b, fm_a, fm_b, mm_a, mm_b, T_exp);
custom_plot( 7, 9,      xm1, xm2, cm_a, cm_b, fm_a, fm_b, mm_a, mm_b, T_exp);
custom_plot( 9, 11,     xm1, xm2, cm_a, cm_b, fm_a, fm_b, mm_a, mm_b, T_exp);
custom_plot( 11, 13,    xm1, xm2, cm_a, cm_b, fm_a, fm_b, mm_a, mm_b, T_exp);
custom_plot( 13, 15,    xm1, xm2, cm_a, cm_b, fm_a, fm_b, mm_a, mm_b, T_exp);
custom_plot( 15, 18.8,  xm1, xm2, cm_a, cm_b, fm_a, fm_b, mm_a, mm_b, T_exp);
custom_plot( 18.8, 19.3,xm1, xm2, cm_a, cm_b, fm_a, fm_b, mm_a, mm_b, T_exp);
custom_plot( 21,25,     xm1, xm2, cm_a, cm_b, fm_a, fm_b, mm_a, mm_b, T_exp);

close_plot(pid);  window(p1);


% evaluate some line fluxes; definitions are from Trans_List_define.sl

ff = get_line_fluxes( trans_list, ["Fe25HeLa", "Ne10HLa", "O8HLa"] );






%%%%%%%%%%%%%%
%%%%%%%%%%%%%%% IMAGING MODE !!!

variable HC = 12.39854;

load_arf("Ard/I3_aim_arf.fits");
load_rmf("Ard/I3_pi_rmf.fits");

set_arf_exposure(5, T_exp);                    
assign_arf(5,5);
assign_rmf(3,5);
ignore([1:4]);

Define_model;
fakeit;

d1 = get_data_counts(5);
m1 = get_model_flux(5);

fit_fun( "wabs(1) * aped_multi(1)" );
set_par(1,1.0);
fakeit;

d2 = get_data_counts(5);
m2 = get_model_flux(5);

%%%%

e1 = HC / d1.bin_hi;
e2 = HC / d1.bin_lo;


%%%%%%%%%%%%%%%%%%%%%%%%%% sim, abs or not, acis-i space.

window(p2); 

label("Energy [keV]",
  "Flux [counts/bin]",
  "Example 100ks fakeits: ACIS-I; DEM");

pid = dup_plot( "Example_spec_100ks_2.ps/cps",  p2);
resize(10*2.54,0.6);

xlog; xrange(0.5,10);
ylog; yrange(1, 10000);

hplot(e1, e2, d1.value, 2);        % unabsorbed
color(2);xylabel(log10(0.52), log10(6000), "n\\dH\\u = 0.0");

ohplot(e1, e2, d2.value, 3);       % absorbed
color(3); xylabel(log10(0.52), log10(100), "n\\dH\\u = 1.0e+22");

close_plot(pid);  window(p2);

% draw to screen...

xlog; xrange(0.5,10);
ylog; yrange(1, 10000);
hplot(e1, e2, d1.value, 2);        % hot, unabsorbed
color(2);xylabel(log10(0.52), log10(6000), "n\\dH\\u = 0.0");

ohplot(e1, e2, d2.value, 3);       % hot, absorbed
color(3); xylabel(log10(0.52), log10(100), "n\\dH\\u = 1.0e+22");


%%%%%%%%%%%%%%%%%%%%%%% MODEL plot in acis-i space...

window(p3); 

label("Energy [keV]",
  "Flux [counts/bin]",
  "Example 100ks fakeits: ACIS-I; DEM");

pid = dup_plot( "Example_spec_100ks_3.ps/cps",  p3);
resize(10*2.54,0.6);

xlog; xrange(0.5,10);
ylog;yrange(1.0e-8, 10);
plot_unit("kev");
plot_model_flux(5,2);

p_groups(HC/10,HC/5);
p_groups(HC/5,HC/2);
p_groups(HC/2,HC/1);
p_groups(HC/1,HC/0.5);

close_plot(pid);  window(p3);


% draw to screen...

window(p3);  

xlog; xrange(0.5,10);
ylog;yrange(1.0e-8, 10);
plot_unit("kev");
plot_model_flux(5,2);

plot_model_flux(5,2);
p_groups(HC/10,HC/5);
p_groups(HC/5,HC/2);
p_groups(HC/2,HC/1);
p_groups(HC/1,HC/0.5);

window(p1);