%% Written (badly!) by C.Cabanac
%% v1.0 2th July 2010
() = evalfile("isis_fancy_plots.sl");
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

public define msg(str_array)
{
   () = printf("\n");
   foreach(str_array)
   {
      variable str = ();
      () = printf(" %s\n", str);
   }
   () = printf("\n");
   return;
}

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



public define plot_counts_sub_comp( )
{
 if( _NARGS !=1 and _NARGS !=2 and  _NARGS !=3 and  _NARGS !=4 )
{
variable str;
 str=[" plot_counts_sub_comp(id,close,resid,type)  ;","",
"								",
" This routine plot all the sub component of a model (in green),",
" together with the whole background spectrum (in red) using the",
" ISIS routine 'plot_counts'.",
" Inputs:","",
" * id			: data index ",
" * close		: specify if the current plot window has to",
" be closed (CLOSE=1) or not (CLOSE!=1) before plotting. It is  ",
" advised to set it to 1 if a plot has already been generated in",
" the current active window.					",
" * resid		: specify the type of residuals to plot   ",
" (value=0, 1, 2). See plot_counts for further explanations.	",
" * type		: specify the type of plot (value=0, 1, 2).",
" It is the same as the keyword POWER int the plot_counts 	",
" routine.							"];
msg(str);
return;
 }
variable id,close,resid,type;
if (_NARGS==3) 
{
(id,close,resid)=();
type = 1;
} else if (_NARGS==2)
{
(id,close)=();
type = 1;
resid=0;
} else if (_NARGS==1)
{
(id)=();
type = 1;
resid=0;
close=1;
} else {
(id,close,resid,type)=();
}
if (close==1) close_plot;
% Save parameters as they initially are
save_par("/tmp/initial.par");
%Note the index of the parameters where the norm of each component is.
variable npar = get_num_pars ();
variable nmod = length(get_fun_components());
variable parmod_norm_name = String_Type [npar];
variable parmod_name = String_Type [npar];

variable i=0;
variable v, w, w_tmp, rnk;

i=0;
while (i<=npar-1){
v=get_params[i];
parmod_name[[i]]=v.name;
parmod_norm_name[[i]]=strchop(parmod_name[i],'.',0)[1];
i++;
}

w=where(parmod_norm_name=="norm");

% Let's make a 3 panel plot, where the top is counts, the middle is
% "unfolded" spectra, and the bottom is residuals with and without
% the lines.
% 3 panel plot, with vertical ratios of 5:5:2
if (resid !=0) multiplot([5,1]);
% pgplot is much happier if you just work on 1 pane at a time.
% Start with the top most pane.
if (resid !=0) mpane(1);
% For single pane plots, the oplt=1 plot option works fine. For
% multi-panel plots, we have to tell pgplot to *not* automatically
% jump to the next pane after doing a plot. This global variable
% does that for my routines.
no_reset=1;
% Uncomment this line for specifying the Units you want to appear on the graph.
% If it has already been defined, do not uncomment it.
%xlog; ylog; Plot_Unit("psd_rms");

% The following lines of code set each but one components normalizations to 0, and plot it in green.

i=0;
while (i<=nmod-1) 
{

if (i==0) 
{ 
w_tmp=w[[1:nmod-1]]+1;
} else if (i==(nmod-1)) {
w_tmp=w[[0:nmod-2]]+1;
} else {
w_tmp=[w[[0:i-1]], w[[i+1:nmod-1]]]+1;
eval_counts;
}


set_par(w_tmp,0); () = eval_counts;

if (i==0)
{
 plot_counts(id;dsym=-4,dcol=1,decol=1,mcol=3,oplt=0,power=type);
} else {
plot_counts(id;dsym=-4,dcol=1,decol=1,mcol=3,oplt=1,power=type);
}
load_par("/tmp/initial.par"); () = eval_counts;

i++;
}

%%% Plot the whole model normally in red
load_par("/tmp/initial.par"); () = eval_counts;
popt.mcol=2;
plot_counts(id;dsym=-4,dcol=1,decol=1,mcol=2,oplt=1,power=type);
% Now do the same thing with the residuals in the lower
% panel;
if (resid !=0) mpane(2);
% Make sure we use the same x-range as above (y-range is different, of course...)


% The next plots after this one we want to go back to "normal" behavior.
no_reset=0;
load_par("/tmp/initial.par"); () = eval_counts;
if (resid !=0) plot_residuals(id;res=resid,rsym=0,rcol=1,recol=1);


}