% ** THIS SCRIPT ASSUMES YOU HAVE PROPERLY LOADED MY .isisrc FILES **

% IMPORTANT NOTE:  In a *script*, all variables must first be declared
% (unless one first sets: _auto_declare =1;), and commands are ended 
% with a semi-colon.  ISIS run *interactively* does *not* need the 
% 'variable' declaration, and can be set so that it automatically 
% appends the semi-colon for you (Isis_Append_Semicolon=1;). I default 
% my .isisrc to Isis_Append_Semicolon=0;  I.e., I require the semi-colon 
% to be used (easier for cutting and pasting from scripts, I think).

% Load radio data

variable radio_data = read_radio("radio.data",50);   % My Function
variable radio_id = load_radio(radio_data);          % My Function

% Load pca data, add 0.5% systematics across the board, group to a S/N
% of 5.5 (~30 counts/bin), and notice data in the 3-22 keV range

variable pca_id = load_data("pca.pha");           % ISIS intrinsic
set_sys_err_frac(pca_id,Double_Type[129]+0.005);  % ISIS intrinsic
group(pca_id;min_sn=5.5,bounds=3.,unit="kev");    % ISIS intrinsic
notice_values(pca_id,3.,22.;unit="kev");          % ISIS intrinsic

% Load the HEXTE data, group it to a minimum signal-to-noise of 20
% per bin, and notice data in the 18-200 keV range

variable hexte_id = load_data("hxt.pha");
group(hexte_id;min_sn=20.,bounds=18.,unit="kev");             
notice_values(hexte_id,18.,200.;unit="kev"); 

% Create an extended user grid, so we can accurately evaluate 
% the reflection convolution model.

define grid_hook(id,s)                               
{
    switch(id)
    {
       case(1):
       s.bin_lo = _A(10^[-9.:-7.:0.01]);
       s.bin_hi = make_hi_grid(s.bin_lo);
    }
    {
       s.bin_lo = _A(10^[-0.5:3:0.002]);
       s.bin_hi = make_hi_grid(s.bin_lo);
    }
    return s;
}

set_eval_grid_method (USER_GRID, [1:3], &grid_hook); % ISIS intrinsic

% For those datasets where reflection matters, make sure to evaluate
% the grid at all energies.  This and the above is a generalized
% version of the XSPEC command 'extend'.  This setting of the kernels
% is, of course, *not* required if you are *not* using a convolution
% model that requires model information beyond the boundaries of the
% nominal grid

set_kernel (2, "std;eval=all");                      % ISIS intrinsic
set_kernel (3, "std;eval=all");                      % ISIS intrinsic

% Define the model to be a broken power law (to get the radio), plus
% diskline, plus reflection. model() is my alias of the ISIS intrinsic
% function fit_fun().  Note: Almost all ISIS functions can be broken
% up onto multiple lines without worry.  Declaring a string variable,
% as is being done in the model call, is the rare exception. (Although
% one can write the string in multiple substrings that are added in
% the funciton call.)

model("reflect(1,constant(Isis_Active_Dataset)*phabs(1)*highecut(1)*(bknpower(1)+diskline(1)))");     

% Set the model parameters.  You can do this either with edit_par, or
% with set_par, which I've put a wrapper around into newpar() (again,
% trying to recover some of the familiar XSPEC syntax).
% set_par/newpar will take either parameter indices *or* strings for
% the parameter names

newpar(2,0.6,-1);
newpar(3,50,0,10,200);
newpar(4,150,0,1,1000);
newpar(5,230,0,0,1.e4);
newpar(6,0.829,-1,0.7,1.0);  % To save computational time, freeze radio slope
newpar(7,5.e-4,0,1.e-6,10);
newpar(8,1.7,0,1.4,2);
newpar(9,2.e-3,0,0,0.1);
newpar(10,6.4,-1);
newpar(11,-3,0,-5,-1);
newpar(12,6,-1);
newpar(13,100,0);
newpar(14,45,-1);
newpar(15,0.4,0,0,2);
newpar(19,0.707,-1);
newpar(21,0.98,0,0.5,1.5);

% Freeze the model constants

freeze([1,20]);           % The opposite is thaw! freeze() & thaw()
			  % will also take parameter strings instead or indices

% The "() = " below are to capture return flags.  Useful in scripts for checking
% success/failure, but they can easily be skipped on the command line.

() = eval_counts;         % ISIS intrinsic; evaluate the model to see how close we are
() = renorm_counts;       % ISIS intrinsic; Just fit normalizations at first
() = fit;                 % My own alias for fit_counts();
 
% Set up plot defaults.

% Use thick lines and larger characters and points, since they look 
% better on screen:

nice_width;                            % My own plotting function to set widths
charsize(1.05); point_size(1.5);       % ISIS intrinsics

% Using my plot functions with just data indices will get you sensible
% defaults, including X & Y ranges set by all of the data sets.
% Otherwise, you can choose plot options by qualifiers, or via
% structure variables.  The latter might be better if you are doing
% multiple plots with just a few minor changes each time.  Here I
% demonstrate both.

% Plot the PCA data. 

xlog; ylog;                  % ISIS intrinsics

plot_counts(2;xrange={3.,22},yrange={1.01e4,4.3e5,-3,8},res=2); % My function

% Make some slightly nicer pgplot colors (a window has to be active to
% run this, so when printing, this has to be run each time a ps plot
% is opened).

pg_color;                    % My function

% Use the global variable popt, set-up in my .isisrc, to set plot
% options

fancy_plot_unit("kev");
popt.xrange={2.e-9,180}; popt.yrange={6.e-3,20,-4,3.99};
popt.dsym = {18,4,6}; popt.dcol = {8,17,4};
popt.rsym = {18,4,6}; popt.rcol = {8,17,4};
popt.res=1; popt.power=2;

% In ISIS, there is no difference between "plotting" and "interactive
% plotting".  Just use the same commands, but write to a postscript
% file instead (which, by the way, you won't have to rotate by 270
% degrees, if you use the vcps device)

variable id = open_print("refmod_kev.ps/vcps");  % My own variation of
						 % ISIS intrinsic
						 % open_plot();
keynote_size; nice_width;   % My own wrappers for plot width & size
plot_unfold({1,2,3},popt);  % My function

close_print(id,"gv");       % My own variation of ISIS intrinsic
			    % close_plot();

% Let's plot the same thing with Hz instead

fancy_plot_unit("hz");      % My variation on ISIS plot_unit()

% Just change X & Y ranges to accomodate

popt.xrange={7.e8,4.e19}; popt.yrange={6.e-3,20,-4,3.99};

id = open_print("refmod_hz.ps/vcps");
keynote_size; nice_width;   % My own wrappers for plot width & size
plot_unfold({1,2,3},popt);  % My function

close_print(id,"gv");

% Let's just zoom in on the RXTE data

fancy_plot_unit("kev");

popt.xrange={3.,180}; popt.yrange={6.e-3,0.13,-4,3.99};
popt.dsym = {4,6}; popt.dcol = {17,4};
popt.rsym = {4,6}; popt.rcol = {17,4};

id = open_print("refmod_pca_kev.ps/vcps");
keynote_size; nice_width;   % My own wrappers for plot width & size
plot_unfold({2,3},popt);    % My function

close_print(id,"gv");

% Let's be wacky and plot the RXTE data in Angstroms, 
% using Counts/Angstrom instead, and do ratio residuals.

fancy_plot_unit("a");
popt.xrange = {_A(180),_A(3.)};    % _A: ISIS intrinsic to convert A <-> keV
popt.yrange = {NULL,NULL,0.9,1.19};
popt.res=3;

id = open_print("refmod_cnts_pca_a.ps/vcps");
keynote_size; nice_width;
plot_data({2,3},popt);    % My function

close_print(id,"gv");

% Reflection is pretty slow, and we're not getting super wonderful
% fits.  Switch to a doubly broken power law, which is faster and
% works better.

model("constant(Isis_Active_Dataset)*phabs(1)*highecut(1)*(bkn2pow(1)+diskline(1))");     

% ISIS remembers all the best fit values from the previous run, so only set
% the new ones.  But, given that we've dumped reflection, tweak the highecut model, too

newpar(3,40);
newpar(4,150);
newpar(5,230);
newpar(6,0.829,-1);             % Start with radio slope frozen
newpar(7,4.e-4,0,1.e-6,1.e-2);
newpar(8,1.68);
newpar(9,11.,0,9.,12.);
newpar(10,1.55);

() = renorm_counts;
() = fit;

thaw(6);                         % Let the radio slope go after initial fit
() = fit;

% Let's just zoom in on the RXTE data

fancy_plot_unit("kev");

popt.xrange={3.,180}; popt.yrange={6.e-3,0.13,0.9,1.19};

id = open_print("bkn2pow_pca_kev.ps/vcps");
keynote_size; nice_width;

plot_unfold({2,3},popt);

close_print(id,"gv");

% Now let's have some fun with functions of parameters.  In the above
% fits, sometimes the diskline model gets unhappy, as the outer radius
% wanders to lower values than the inner radius.  Watch how we avoid
% that.

% The constant(10) model doesn't change the model values, since it's
% multiplied by 0.

model("constant(Isis_Active_Dataset)*phabs(1)*highecut(1)*(bkn2pow(1)+diskline(1))+0.*constant(10)"); 

% But it does give us a parameter to tie to other ones.  We use it to get a 
% Delta_R, for difference between inner and outer radius.  
% for a few more examples.

newpar("constant(10).factor",get_par(15)-get_par(14),0,0.1,994);
set_par_fun(15,"_par(14)+_par(17)");

() = renorm_counts;
() = fit;

% We could use the ISIS conf() or vconf() functions to do error bar
% searches.  But we might find lower chi^2 and have to refit.  To that
% end, I have written "error_loop" to do the error bars, and refit if
% necessary, and write out intermediate steps along the way.

% Remind ourselves of the free parameters

list_free;                                         % ISIS intrinsic

(,) = conf_loop(,1,1.e-3;save,prefix="bkn2pow_diskline.");   % ISIS intrinsic

% You know, the diskline is a slow model, and not really doing us a
% lot of good, so let's swap it out for a gaussian, so we can do a fast
% demonstration of error contours.

% NOTE:  ISIS has a gauss model appropriate for fitting lines to gratings 
% data in wavelength space.  So, here's the one instance of where the ISIS
% model name is different than that from XSPEC.

model("constant(Isis_Active_Dataset)*phabs(1)*highecut(1)*(bkn2pow(1)+gaussian(1))");

% Again, ISIS keeps old parameter values, so only define what's new

newpar("gaussian(1).norm",1.e-3);
newpar("gaussian(1).LineE",6.4,-1);
newpar("gaussian(1).Sigma",0.3,0,0,1.5);

() = renorm_counts;
() = fit;

% As before, do conf_loop.  The error_save script uses save_par to save best fit values

(,) = conf_loop(,1,1.e-3;save,prefix="bkn2pow_gaussian.");   

% Let's make a slightly wacky plot: BTU/Acre/s vs. Slug*Sverdrup*Hz/TerraSmoot

add_plot_unit("sshps","btupas";xscale=5.353038e10,yscale=3274.3436,is_energy=1);
fancy_plot_unit("sshps","btupas");

xrange; % Reset xrange to autoscale. (Previously *set* xranges are preserved - I'm
        % not sure if that's a wise idea or not.)

id= open_print("BTU_vs_Stuff.ps/vcps");
keynote_size; nice_width;
plot_unfold({2,3};power=3,dsym={4,6},dcol={17,4},
                  xlabel="\\fs Energy (Slug\\.Sv\\.Hz/TSmoot)",
                  ylabel="\\fs Flux (BTU/Acre/s)");
close_print(id,"gv");

% Let's make a somewhat more practical plot, and do velocity/redshift axes.

fancy_plot_unit("kev","photons");

id= open_print("data_vs_velocity.ps/vcps");
keynote_size; nice_width;

% Reference zero velocity to 6.4 keV
plot_data({2,3};power=3,dsym={4,6},dcol={17,4},vzero=6.4); 

close_print(id,"gv");

id= open_print("data_vs_redshift.ps/vcps");
keynote_size; nice_width;

% Reference zero velocity to 6.4 keV, but use a redshift axis instead
plot_data({2,3};power=3,dsym={4,6},dcol={17,4},vzero=6.4, zaxis=1); 

close_print(id,"gv");

% Now let's do error contours of Gamma_2 vs. Gamma_3.  First, set up a
% 31 X 31 confidence contour grid.

variable px = conf_grid(8,1.67,1.70,31);         % Gamma_2
variable py = conf_grid(10,1.43,1.52,31);        % Gamma_3

% Now do the confidence contours.  NOTE: I have placed a simple
% wrapper around the ISIS intrinsic conf_map_counts() to use the built
% in ISIS functionality of a failure hook.  If a fit fails along the
% way, ISIS temporarily switches from the lmdif fit method to the subplex
% fit method (or visa versa, if I'm using subplex to begin with)

variable cntr = conf_map(px,py);

% Save the results into a fits file.

() = save_conf(cntr,"bkn2pow_gauss_g2_vs_g3.contour");  % ISIS intrinsic; see also load_conf()

% Plot the results

xlin; ylin; xrange; yrange;

id = open_print("bkn2pow_gauss_g2_vs_g3_contour.ps/vcps");
keynote_size; nice_width;

xlabel("\\fr Photon Index, \\gG\\d2");    % ISIS intrinsic
ylabel("\\fr Photon Index, \\gG\\d3");    % ISIS intrinsic
plot_conf(cntr);                          % ISIS intrinsic; see also oplot_conf


close_print(id);