%; Time-stamp: <2001-10-30 15:54:45 dph> 
%; MIT Directory: ~dph/d1/CXO_Data/AR_Lac/Sl/
%; File: Contin_sim.sl
%; Author: D. Huenemoerder
%; Original version: 2001.10.29
%;====================================================================
% version: 1.0
%
% purpose:  use aped_multi, aped_contin, and fakeit to examine true
%	    continuum vs estimated continuum for multi-Temp model.
%
%

% control variables

%  to make filecopy, set this, then execute the plot blocks w/ a
%  do_filecopy  test, to make gif or ps files
%
if ( not is_defined("do_filecopy") ) { do_filecopy=0 ;};

%
%


%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PART 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SETUP data, responses:
% 
%    follow Contin_tutorial through "plasma(aped);"
%    http://space.mit.edu/CXC/ISIS/contrib/Contin_fit/Contin_tutorial.html
%    http://space.mit.edu/CXC/ISIS/contrib/Contin_fit/Contin_tutorial.sl
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PART 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Establish plasma model, multi-T, based on AR Lac ad hoc model:
% (see http://space.mit.edu/CXC/ISIS/contrib/Aped_fit/Aped_model_demo.html)
%

evalfile("aped_fit_models.sl");

% from def_mdl.sl  --- ad hoc AR Lac DEM:

variable mdl_temps = [
    3.16228e+05,   3.98107e+05,   5.01187e+05,   6.30958e+05,   7.94328e+05,
    1.00000e+06,   1.25893e+06,   1.58489e+06,   1.99526e+06,   2.51189e+06,
    3.16228e+06,   3.98107e+06,   5.01187e+06,   6.30958e+06,   7.94328e+06,
    1.12202e+07,   1.41254e+07,   1.77828e+07,   2.23872e+07,   2.81838e+07,
    3.54814e+07,   4.46684e+07,   5.62341e+07,   7.07946e+07,   8.91251e+07,
    1.12202e+08,   1.41254e+08,   1.77828e+08,   2.23872e+08,   2.81838e+08 
    ];

  variable mdl_wgts = [
    3.257e-07,   5.016e-06,   6.083e-06,   5.357e-06,   4.835e-06,
    5.477e-06,   4.997e-06,   6.666e-06,   1.780e-05,   3.002e-05,
    4.921e-05,   8.800e-05,   1.088e-04,   3.620e-04,   5.124e-04,
    4.618e-04,   3.576e-04,   4.747e-04,   6.833e-04,   7.939e-04,
    5.566e-04,   2.740e-04,   1.361e-04,   4.348e-05,   6.811e-06,
    5.981e-07,   5.107e-08,   5.419e-09,   3.427e-10,   1.232e-11
    ];

% renorm gets wgts of: 
% 5.886142e-06 9.681320e-05 9.030718e-05 5.425298e-04 1.966264e-03
% 8.345779e-03 1.234879e-02 4.951805e-03 1.230903e-04 9.793368e-08

% vs input
% 3.257000e-07 5.357000e-06 4.997000e-06 3.002000e-05 1.088000e-04
% 4.618000e-04 6.833000e-04 2.740000e-04 6.811000e-06 5.419000e-09


variable elem_abund =[1.1, 0.8, 0.8, 0.4, 0.60, 0.80, 1.1, 0.35, 0.5, 0.5];
variable elem =[Al, Ar, Ca, Fe, Mg, N, Ne, O, S, Si];

% truncate to 10 points for faster evaluation:
%
t=mdl_temps[[0:29:3]]; w=mdl_wgts[[0:29:3]];
aped_multi_setup(t,w*3);  % factor 3 for having taken 1/3 points
aped_set_abund(elem, elem_abund);

group_data(3,0);
xnotice(3,1.7,25);
ignore([1,2,4]);
eval_counts;

plot_data_counts(3,2);oplot_model_counts(3,3);
xnotice(3,10,15);
freeze([11:20]);
renorm_counts;

xrange(1.7,25);
xnotice(3,1.7,25); eval_counts;
plot_data_counts(3,2);oplot_model_counts(3,3);


% read actual norms into new w array:
%
wn = @w;    %  copy w array
for (i=1; i<11; i++) wn[i-1] = get_par(i);


% notice all arrays, evaluate model
%
group_data(1,0); group_data(2,0); group_data(3,0); group_data(4,0);
xnotice([3,4], 1.7,25);  xnotice([1,2], 1.7, 17);
eval_counts;

% read the model counts (model folded through response)
%
fhm=get_model_counts(1);  fhp=get_model_counts(2);
fmm=get_model_counts(3);  fmp=get_model_counts(4);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PART 3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Establish plasma continuum model, multi-T, based on above "fit":
%

% define continuum model, evaluate, and retrieve into an array:
%
evalfile("aped_contin.sl");
fun_c = "aped_contin(1)";
for ( i=2; i<=length(wn); i++ )
  fun_c = fun_c + " + aped_contin(" + string(i) + ")";

fit_fun(fun_c);	  % establish the model

% set the params:
%
for (i=0; i<length(wn); i++)
{
  set_par(2*i+1, wn[i], 1, 0, wn[i]*100);
  set_par(2*i+2, t[i],  1, 1.e4, 3.e8);
}

eval_counts;

% read the continuum model counts (model folded through response)
%
chm=get_model_counts(1);  chp=get_model_counts(2);
cmm=get_model_counts(3);  cmp=get_model_counts(4);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PART 4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Compute fake data, based on AR Lac responses and exposure time:
%

% re-establish aped_multi model, but w/ new weights:
%
aped_multi_setup(t,wn);


% set up fake data histograms by assigning arfs and rmfs.
% use histogram indices 5-8 for HEG -1, 1, MEG -1, 1:
% (order of assigns is important!)
%
for (i=0; i < length(fgarfs_heg); i++)
{
  assign_rmf(rmf_heg, i+1 +4);
  assign_arf(i+1,i+1 +4, 1);
}

for (i=0; i < length(fgarfs_meg); i++ )
{
  assign_rmf( rmf_meg, nheg + i+1 +4 );
  assign_arf(nheg + i+1, nheg + i+1 +4, 1 );
}


ignore([1:4]) ;   % don't work w/ real observation.
xnotice([7,8], 1.7,25);  xnotice([5,6], 1.7, 17);
fakeit;


%%%%%%%%%%%%%%%%%%%%%%%%% PART 5 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Now do continuum line-free-region fitting on fake data,
% then compare to model continuum.
%
% --- this follows Contin_tutorial.sl ---
% (http://space.mit.edu/CXC/ISIS/contrib/Contin_fit/Contin_tutorial.html)


fit_fun("aped_contin(1)");               % model is thermal plasma continuum

% set some likely defaults:
                 freeze
%       n  value /thaw  min   max
set_par(1, 0.004, 0,  	0,    0.1);
set_par(2, 2.3e7, 0,  	1.e4, 5.e8);


% Group data coarsely, since this is a continuum fit, and we would
%	like good statistics.

meg_group_factor = 8;   heg_group_factor = 2*meg_group_factor;

% redefine order indices to map to the loaded histograms:

mm = 7; mp = 8;   hm = 5;  hp = 6;

group_data(mm, meg_group_factor); group_data(mp, meg_group_factor); 
group_data(hm, heg_group_factor); group_data(hp, heg_group_factor); 

%%%%%%%%%%%%%%%%%%%%%  
%
% First, use only line-free regions.  Function applies ignore region table.

line_free_regs([hm,hp,mm,mp]);

% For first pass, single high-T component, only notice the
% short-wavelength region, and fit both HEG and MEG:

ignore([hm,hp,mm,mp], 10, 25); 

list_data ;                           % check
yrange(0); xrange(1,8);
renorm_counts;
plot_data_counts(mm,1);  oplot_data_counts(mp,2);
oplot_data_counts(hm,3);  oplot_data_counts(hp,4);


% Do a fit, list the params, compute confidence:
fit_counts; list_par;

oplot_model_counts(mm,5);  oplot_model_counts(mp,6);
oplot_model_counts(hm,7);  oplot_model_counts(hp,8);


% evaluate 90% confidence intervals:
%
(clo_norm, chi_norm) = conf(1);   ()=printf("%% Conf norm = %f , %f\n", clo_norm, chi_norm);
(clo_T, chi_T) = conf(2);         ()=printf("%% Conf T = %e , %e\n", clo_T, chi_T);

% Conf norm = 0.020154 , 0.022441
% Conf T = 2.131334e+07 , 2.460589e+07


%  Filecopy:
%

if (do_filecopy) 
{
  pid = dup_plot("Contin_fit_0.gif/gif",p1);   resize(6*2.54, 0.6);
}  

plot_data_counts(mm,1);  oplot_data_counts(mp,2);
oplot_data_counts(hm,3);  oplot_data_counts(hp,4);
oplot_model_counts(mm,5);  oplot_model_counts(mp,6);
oplot_model_counts(hm,7);  oplot_model_counts(hp,8);

if (do_filecopy) 
{
  close_plot(pid); window(p1);
}  


%%%%%%%%%%%%%%%%%%%%%%%%% PART 5b %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%  If satisfied, freeze this, and add a new component, and
%%%%%%%%%%  fit long-wavelength region where low-T has more affect. 

freeze([1,2]);
fit_fun("aped_contin(1) + aped_contin(2)");   % insert another T component

set_par(3, 0.004, 0,  	0,    0.1);           % set initial params
set_par(4, 2.3e6, 0,  	1.e4, 5.e8);
list_par;

ignore([hm,hp]);                    % no good region for HEG - ignore it.
notice([mm,mp]);		    % notice the MEG, then apply regions
line_free_regs([mm,mp]);	    % line free,
ignore([mm,mp], 1, 8);              % ignore short waves.

yrange(0); xrange(8,25);
plot_data_counts(mm,1);    oplot_data_counts(mp,2);
renorm_counts;

fit_counts; list_par;

oplot_model_counts(mm,5);  oplot_model_counts(mp,6);

(clo_norm2, chi_norm2) = conf(3);   ()=printf("%% Conf norm = %f , %f\n", clo_norm2, chi_norm2);
(clo_T2, chi_T2) = conf(4);         ()=printf("%% Conf T = %e , %e\n", clo_T2, chi_T2);

% Conf norm = 0.002271 , 0.012189
% high T fails ---> don't need high-T continuum component.


%%%%%%%%%%%%%%%%%%%%%  PART 5c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
%  set continuum fit - 1 or 2 components, 
%
% evaluate the model on the full MEG grid and
% save store it for use.

fit_fun("aped_contin(1)");          % reset to component 1.

save_par("Contin_sim_contin.par"); % save par file to disk.

% notice full region

group_data(mm, 0); group_data(mp, 0); 
group_data(hm, 0); group_data(hp, 0); 
xnotice([hm,hp], 1.7, 17);   xnotice([mm,mp], 1.7, 25);

eval_counts;             %   so fold the model:

% look at it:
%
yrange(0);xrange(5,7);plot_data_counts(mm,2);oplot_model_counts(mm,3);


% read the fit continuum model counts (model folded through response)
%
chm_fit=get_model_counts(hm);  chp_fit=get_model_counts(hp);
cmm_fit=get_model_counts(mm);  cmp_fit=get_model_counts(mp);


%%%%%%%%%%%%%%%%%%%%%  PART 6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
%  Examine the model true continuum,  fit continuum, 
%              model counts spectrum, fake counts spectrum
%
%
% variables:  model spectrum:       fhm      fhp      fmm      fmp
%             fake  spectrum: hists 5	     6	      7	       8
%	      model continuum:      chm      chp      cmm      cmp
%	      fit continuum:	    chm_fit  chp_fit  cmm_fit  cmp_fit

% Define bin boundary arrays, for convenience; MEG and HEG
%
variable xm1=fmm.bin_lo, xm2=fmm.bin_hi, xh1=fhm.bin_lo, xh2=fhm.bin_hi;


% Device choices:

froot = "Continuum_simulations";
ps_dev =".ps/cps";
gif_dev = ".gif/gif";         % pick a suffix - edit for successive gifs.

dev = ps_dev;       % pick  a device

if (do_filecopy) 
{
  pid = dup_plot( froot + dev, p1 );   resize(10*2.54, 0.6);
}  ;


%%% Plot model continuum and fit continuum:

multiplot([3,1]);

%% -1 orders:

label("Wavelength [\\A]",      "Counts/bin",
      "Continuum Simulation: Model(red), Fit(green) (MEG,HEG -1)");
yrange(0);xrange(1.7,25);
hplot(xm1,xm2,cmm.value,2);ohplot(xm1,xm2,cmm_fit.value,3);
ohplot(xh1,xh2,chm.value,2);ohplot(xh1,xh2,chm_fit.value,3);

label("Wavelength [\\A]",      "Residuals",
      "Continuum Simulation: Model(red), Fit(green) (MEG,HEG -1)");
yrange();
hplot(xm1,xm2,cmm_fit.value-cmm.value,2);
ohplot(xh1,xh2,chm_fit.value-chm.value,2);


%% +1 orders

label("Wavelength [\\A]",      "Counts/bin",
      "Continuum Simulation: Model(red), Fit(green) (MEG,HEG +1)");
yrange(0);
hplot(xm1,xm2,cmp.value,2);ohplot(xm1,xm2,cmp_fit.value,3);
ohplot(xh1,xh2,chp.value,2);ohplot(xh1,xh2,chp_fit.value,3);
label("Wavelength [\\A]",      "Residuals",
      "Continuum Simulation: Model(red), Fit(green) (MEG,HEG +1)");
yrange();
hplot(xm1,xm2,cmp_fit.value-cmp.value,2);
ohplot(xh1,xh2,chp_fit.value-chp.value,2);


multiplot(1);

%%%%%%%%%% Plot fake data and continuua:

% ad hoc plotting utility function:
define pfunc(xmm, h, s, x1, x2, cf, cm)
{
  xrange(xmm[0], xmm[1]);
  splot_data_counts(h, s, 2);
  ohplot(x1, x2, cf.value, 3);
  ohplot(x1, x2, cm.value, 1);
  plot_group( brightest(20, where( wl( xmm[0], xmm[1] ) ) ) );
}


% MEG -1
%
label("Wavelength [\\A]",      "Counts/bin",
      "Fake Data and Contin: Data(red), CFit(green) CMdl(other) (MEG,-1)");

xmm=[1.7,6]; pfunc(xmm, mm, 0.01, xm1, xm2, cmm_fit, cmm);
xmm=[6, 12]; pfunc(xmm, mm, 0.01, xm1, xm2, cmm_fit, cmm);
xmm=[12,18]; pfunc(xmm, mm, 0.01, xm1, xm2, cmm_fit, cmm);
xmm=[18,25]; pfunc(xmm, mm, 0.01, xm1, xm2, cmm_fit, cmm);
xmm=[20,25]; pfunc(xmm, mm, 0.01, xm1, xm2, cmm_fit, cmm);


% MEG +1
%
label("Wavelength [\\A]",      "Counts/bin",
      "Fake Data and Contin: Data(red), CFit(green) CMdl(other) (MEG,+1)");

xmm=[1.7,6]; pfunc(xmm, mp, 0.01, xm1, xm2, cmp_fit, cmp);
xmm=[6,12]; pfunc(xmm,  mp, 0.01, xm1, xm2, cmp_fit, cmp);
xmm=[12,18]; pfunc(xmm, mp, 0.01, xm1, xm2, cmp_fit, cmp);
xmm=[18,25]; pfunc(xmm, mp, 0.01, xm1, xm2, cmp_fit, cmp);
xmm=[20,25]; pfunc(xmm, mp, 0.01, xm1, xm2, cmp_fit, cmp);


% HEG -1
%
label("Wavelength [\\A]",      "Counts/bin",
      "Fake Data and Contin: Data(red), CFit(green) CMdl(other) (HEG,-1)");

xmm=[1.7,5]; pfunc(xmm, hm, 0.005, xh1, xh2, chm_fit, chm);
xmm=[5,9]; pfunc(xmm,  hm, 0.005, xh1, xh2, chm_fit, chm);
xmm=[9,13]; pfunc(xmm, hm, 0.005, xh1, xh2, chm_fit, chm);
xmm=[13,17]; pfunc(xmm, hm, 0.005, xh1, xh2, chm_fit, chm);


% HEG +1
%
label("Wavelength [\\A]",      "Counts/bin",
      "Fake Data and Contin: Data(red), CFit(green) CMdl(other) (HEG,+1)");

xmm=[1.7,5]; pfunc(xmm, hp, 0.005, xh1, xh2, chp_fit, chp);
xmm=[5,9];  pfunc(xmm, hp, 0.005, xh1, xh2, chp_fit, chp);
xmm=[9,13]; pfunc(xmm, hp, 0.005, xh1, xh2, chp_fit, chp);
xmm=[13,17]; pfunc(xmm, hp, 0.005, xh1, xh2, chp_fit, chp);


%%% Plot Fake data, folded model, folded fit continuum:

% MEG -1
%
label("Wavelength [\\A]",      "Counts/bin",
      "Fake data (red), Folded Model (green) Fit Contin (other) (MEG,-1)");

xmm=[1.7,6]; pfunc(xmm,mm,0.005, xm1, xm2, fmm, cmm_fit);
xmm=[6, 12]; pfunc(xmm,mm,0.005, xm1, xm2, fmm, cmm_fit);
xmm=[9,10]; pfunc(xmm,mm,0.005, xm1, xm2, fmm, cmm_fit);
xmm=[12,18]; pfunc(xmm,mm,0.005, xm1, xm2, fmm, cmm_fit);
xmm=[13,14]; pfunc(xmm,mm,0.005, xm1, xm2, fmm, cmm_fit);
xmm=[18,25]; pfunc(xmm,mm,0.005, xm1, xm2, fmm, cmm_fit);
xmm=[20,25]; pfunc(xmm,mm,0.005, xm1, xm2, fmm, cmm_fit);
xmm=[21.2,22.4]; pfunc(xmm,mm,0.005, xm1, xm2, fmm, cmm_fit);

% MEG +1
%
label("Wavelength [\\A]",      "Counts/bin",
      "Fake data (red), Folded Model (green) Fit Contin (other) (MEG,+1)");

xmm=[1.7,6]; pfunc(xmm,mp,0.005, xm1, xm2, fmp, cmp_fit);
xmm=[6, 12]; pfunc(xmm,mp,0.005, xm1, xm2, fmp, cmp_fit);
xmm=[9,10];  pfunc(xmm,mp,0.005, xm1, xm2, fmp, cmp_fit);
xmm=[12,18]; pfunc(xmm,mp,0.005, xm1, xm2, fmp, cmp_fit);
xmm=[13,14]; pfunc(xmm,mp,0.005, xm1, xm2, fmp, cmp_fit);
xmm=[18,25]; pfunc(xmm,mp,0.005, xm1, xm2, fmp, cmp_fit);
xmm=[20,25]; pfunc(xmm,mp,0.005, xm1, xm2, fmp, cmp_fit);
xmm=[21.2,22.4]; pfunc(xmm,mp,0.005, xm1, xm2, fmp, cmp_fit);


% HEG -1
%
label("Wavelength [\\A]",      "Counts/bin",
      "Fake data (red), Folded Model (green) Fit Contin (other) (HEG,-1)");

xx=[1.7,5]; pfunc(xx,hm,0.0025, xh1, xh2, fhm, chm_fit);
xx=[5, 9]; pfunc(xx,hm,0.0025, xh1, xh2, fhm, chm_fit);
xx=[9, 13]; pfunc(xx,hm,0.0025, xh1, xh2, fhm, chm_fit);
xx=[9, 10]; pfunc(xx,hm,0.0025, xh1, xh2, fhm, chm_fit);
xx=[12, 12.5]; pfunc(xx,hm,0.0025, xh1, xh2, fhm, chm_fit);
xx=[13,17]; pfunc(xx,hm,0.0025, xh1, xh2, fhm, chm_fit);

% HEG +1
%
label("Wavelength [\\A]",      "Counts/bin",
      "Fake data (red), Folded Model (green) Fit Contin (other) (HEG,+1)");

xx=[1.7,5];pfunc(xx,hp,0.005, xh1, xh2, fhp, chp_fit);
xx=[5, 9]; pfunc(xx,hp,0.0025, xh1, xh2, fhp, chp_fit);
xx=[9, 13];pfunc(xx,hp,0.0025, xh1, xh2, fhp, chp_fit);
xx=[9,10]; pfunc(xx,hp,0.0025, xh1, xh2, fhp, chp_fit);
xx=[12, 12.5];pfunc(xx,hp,0.0025, xh1, xh2, fhp, chp_fit);
xx=[13,17];pfunc(xx,hp,0.0025, xh1, xh2, fhp, chp_fit);


if (do_filecopy) 
{
  close_plot(pid); window(p1);
}  ;