%; Time-stamp: <2002-03-05 16:57:48 dph> 
%; File: Fold_fit.sl
%; Author: D. Huenemoerder
%; Original version: 2002.03.04 
%;====================================================================
% version: 1
%
% purpose:
%   example: fake a spectrum; fit a line group via rmf*arf folding.


variable wmin=6.5, wmax=6.9;

% For demonstration, only do small range.  If you want to fit other
% lines, open up the range here so the fake data is for the whole
% observable spectrum.

xnotice([1,2], wmin, wmax);  
xnotice([3,4], wmin, wmax);

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

list_data;
list_arf;
list_rmf;

% limit the plot range to the fit region
yrange(0);
xrange(wmin,wmax);

% fit HEG and MEG simultaneously, over the fit region
ignore(all_data);
xnotice([1:4], wmin, wmax);

% look at the fake data and the folded model:
%
plot_data_counts( 3, red );  oplot_model_counts( 3, green );


% define the fit function as a constant + 3 gaussians.
% In order to fit unresolved lines, the gaussian sigmas will be frozen
% very small. 
%
fit_fun( "poly(1) + gauss(1) + gauss(2) + gauss(3)");

% polynomial: make a constant term only:
%
set_par(1, 0.001, 0, 0, 0.01 );
set_par(2,0,1);
set_par(3,0,1);

% Guess some gaussian amplitudes:
%
set_par(4, 0.0001,  0, 0, 0.001);
set_par(7, 0.00002, 0, 0, 0.001);
set_par(10, 0.0001, 0, 0, 0.001);


% use the trans_list assoc array to set line wavelengths:

variable w1= line_info( (where(eval(trans_list["Si13HeLar"])))[0]).lambda,
         w2= line_info( (where(eval(trans_list["Si13HeLai"])))[0]).lambda,
	 w3= line_info( (where(eval(trans_list["Si13HeLaf"])))[0]).lambda;


variable dw = 0.03;  % a range for line fit limits.

% set the initial positions:
%
set_par(5, w1,  0, w1-dw, w1+dw);
set_par(8, w2,  0, w2-dw, w2+dw);
set_par(11, w3, 0, w3-dw, w3+dw);

% set the gaussian sigmas small, and freeze them:

set_par(6, 0.001, 1);
set_par(9, 0.001, 1);
set_par(12, 0.001, 1);

% renormalize (depends on params being designated as "norms" --- done
% in the fit function).

renorm_counts;



fit_counts;  % At last - do the fit.

list_par;    % list result to screen


% evaluate and format confidence limits:

()=printf( "#\n");
()=printf( "# Conf(n,%d) [0=>68%%, 1=>90%%, 2=>99%%]\n",  1);
()=printf( "#\n");
()=printf( "# Conf par %d: %12.4e %12.4e\n", 1, conf(1));
()=printf( "# Conf par %d: %12.4e %12.4e\n", 4, conf(1));
()=printf( "# Conf par %d: %12.4e %12.4e\n", 7, conf(1));
()=printf( "# Conf par %d: %12.4e %12.4e\n", 10, conf(1));
()=printf( "# Conf par %d: %12.4e %12.4e\n", 5, conf(1));
()=printf( "# Conf par %d: %12.4e %12.4e\n", 8, conf(1));
()=printf( "# Conf par %d: %12.4e %12.4e\n", 11, conf(1));


% make some plots; one to the screen; that and more to a file:

window(p1);

label("Wavelength [\\A]",
      "Flux [couts/bin]",
      "Fake DEM model, Si XIII He-like triplet");

plot_data_counts(3,2);
oplot_model_counts(3,3);
oplot_data_counts(1,5);
oplot_model_counts(1,6);


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

plot_data_counts(3,2);
oplot_model_counts(3,3);
oplot_data_counts(1,5);
oplot_model_counts(1,6);

plot_data_counts(4,2);
oplot_model_counts(4,3);
oplot_data_counts(2,5);
oplot_model_counts(2,6);

color(8);
rplot_counts(4);

close_plot(pid); window(p1);


% write params to file; then open file and append conf limits:
%
variable f_par = "Example_Si13fir.par";
save_par( f_par );

variable fp = fopen( f_par, "a" );

()=fprintf(fp, "#\n");
()=fprintf(fp, "# Conf(n,%d) [0=>68%%, 1=>90%%, 2=>99%%]\n",  1);
()=fprintf(fp, "#\n");

()=fprintf(fp, "# Conf par %d: %12.4e %12.4e\n", 1, conf(1));
()=fprintf(fp, "# Conf par %d: %12.4e %12.4e\n", 4, conf(1));
()=fprintf(fp, "# Conf par %d: %12.4e %12.4e\n", 7, conf(1));
()=fprintf(fp, "# Conf par %d: %12.4e %12.4e\n", 10, conf(1));
()=fprintf(fp, "# Conf par %d: %12.4e %12.4e\n", 5, conf(1));
()=fprintf(fp, "# Conf par %d: %12.4e %12.4e\n", 8, conf(1));
()=fprintf(fp, "# Conf par %d: %12.4e %12.4e\n", 11, conf(1));

fclose(fp);

%%% all done.