%
% 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, 1/17/04 - dd
%

% - - - - -
% Get the aciscontam model...
()=evalfile("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 = 2004.4;  %tau 700 eV = 0.60  for twolevel model
ltime = 2002.97; %           = 0.48
%ltime = 2000.6;  %           = 0.24
%ltime = 1999.8;  %           = 0.12


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:
% 
% 20% thick 0.92 version:  (*** used to create N9998 ...)
%%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) ");
%
% 50-50 0.885 version:
%%set_par("aciscontam(1).OK", 0.81);
%%set_par("aciscontam(1).FK", 0.81);
%%fit_fun(" (0.5*(aciscontam(1))^0.805 + 0.5*(aciscontam(1))^2.0) ");
%
% 20-80 0.885 version:
%%set_par("aciscontam(1).OK", 0.95);
%%set_par("aciscontam(1).FK", 0.95);
%%fit_fun(" (0.8*(aciscontam(1))^0.95 + 0.2*(aciscontam(1))^3.9) ");
%
% 15-85 0.885 version:   (*** used for N9997 ...)
%   the IDL contam_twolevel.pro parameters were:
% thkfill = 0.15
% thkfac = 5.7
% thnfac = 0.97
%
set_par("aciscontam(1).OK", 0.97);
set_par("aciscontam(1).FK", 0.90);
fit_fun(" (0.85*(aciscontam(1))^0.97 + 0.15*(aciscontam(1))^5.7) ");
%

% - - - - - - - - - - - - - - - - - - - - - - -

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.8,1.05);

% 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);

% Use Alexey's ECS line energies:
%
% weights for transmission for these are: 0.05, 0.47, 0.48
oplot(hc/0.618*[1.,1.],[0.0,1.1],blue);
oplot(hc/0.682*[1.,1.],[0.0,1.1],blue);
oplot(hc/0.7056*[1.,1.],[0.0,1.1],blue);
%

% 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 (only) 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+200;
  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("");


% Use weighted value at the three energies:
variable lline17, lline18, lline20, lline;
(w1,w2)=linear_grid(hc/0.7069, hc/0.7047, 5);
lline = eval_fun(w1,w2);
lline17 = sum(lline)/length(lline);

(w1,w2)=linear_grid(hc/0.6834, hc/0.6813, 5);
lline = eval_fun(w1,w2);
lline18 = sum(lline)/length(lline);

(w1,w2)=linear_grid(hc/0.6317, hc/0.6194, 5);
lline = eval_fun(w1,w2);
lline20 = sum(lline)/length(lline);


% weighted average of the three values
lline = 0.48 * lline17 + 0.47 * lline18 + 0.05 * lline20;

% Average Tau at L lines:
print("---------------------------");
print(ltime);
print("total tau(weighted) = "+string(-1.0*log( lline )));