%
% file: contam_twolevel.i
%

% Create the "two-level" contamintation model in ISIS...
%
% Thanks to Dave H for acis_contam.sl, etc.
%
% 10/10/03, 12/2/03 - dd
%

% - - - - -
% Get the aciscontam model...
()=evalfile("/nfs/cxc/h1/dph/libisis/acis_contam.sl");

% - - - - -
variable ltime, llines, hc;

hc = 12.3985;

% Set it for the nominal contamination at ~ 12/20/2002 (E0102 II, obsid 3828)
ltime = 2002.97;


fit_fun("aciscontam(1)");

% Time since 1998.0
set_par(1, (ltime-1998.0)*365.25*24.0*3600.0);

% Relative amounts of the absorbing elements:
set_par(2, 1.0);
set_par(3, 1.0);
set_par(4, 1.0);

list_par;

% Some variable needed...
% (declare them here so this file and be .load'ed or evalfile'd.)
variable w1, w2, y, ytl, vps;

% Get the function values into an array
(w1,w2)=linear_grid(1.0, 50.0, 2048);
y=eval_fun(w1,w2);

% - - - Modified version of aciscontam - - - -
%
% Now for the modified version - the sum of two thicknesses
% Use the nominal curve above as the basis for this...
%   and 
% Modify slightly the O and F amounts to give the same edge jumps:
set_par("aciscontam(1).OK", 0.92);
set_par("aciscontam(1).FK", 0.90);
%
fit_fun(" (0.8*(aciscontam(1))^0.94 + 0.2*(aciscontam(1))^2.80) ");
%
% - - - - - - - - - - - - - - - - - - - - - - -

list_par;

% and extract this to another array
ytl=eval_fun(w1,w2);


% - - - - - Make a Plot - - - -
% Open a postscript output
vps = open_plot("twolevel_isis.ps/vcps", 1, 2);
window(vps);

label("Wavelength (A)","Transmission, Ratio",
  "Ratio(green) of Contamarf(black) and Two-level version(red)");
xlin;
ylin;
xrange(0.,50.);
yrange(0.,1.1);

% plot the original contamination transmission:
hplot(w1,w2,y);

% and the two-level version:
ohplot(w1,w2,ytl);

% and plot their ratio...
ohplot(w1,w2,ytl/y);


% and show where the ECS ~700 eV is:
oplot(hc/0.7*[1.,1.],[0.0,1.1],blue);
% and 670 eV and 730:
oplot(hc/0.670*[1.,1.],[0.0,1.1],blue);
oplot(hc/0.730*[1.,1.],[0.0,1.1],blue);

% Close the postscript plot
close_plot(vps);



% Print out the ratio at some values:
vmessage("");
vmessage(" Fluffium transmission at t = "+string(ltime));
vmessage("");
variable ii;
ii=0; while ( ii < length(y) ) {
    vmessage("%6.2f   %6.4f  ",
	0.5*(w1[ii]+w2[ii]),
	ytl[ii]/y[ii]  );
  ii=ii+50;
  }
% last point too
ii = length(w1)-1;
vmessage("%6.2f   %6.4f  ",
	0.5*(w1[ii]+w2[ii]),
	ytl[ii]/y[ii]  );
vmessage("");