%
% file: contam_llines.i
%

% Plot the ~700 eV L-line Tau vs time based on
% the "two-level" model...
%
% (adapted from contam_twolevel.i)
%
% 10/15/03, 11/14/03 - dd
%

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

% - - - - -
variable ltimes, ltaus, cktaus;

variable ltime, lline, ckhi, cklo, hc;

hc = 12.3985;

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

% loop over the time
ltime = 2000.5;

ltimes=[ltime];
ltaus=[0.0];
cktaus=[0.0];

do
  {
%  -----

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);
%
% and redefine to use the two-level model:
fit_fun(" (0.80*(aciscontam(1))^0.94 + 0.20*(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): " +
	" Time = "+string(ltime));
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 695 eV and 705:
oplot(hc/0.695*[1.,1.],[0.0,1.1],blue);
oplot(hc/0.705*[1.,1.],[0.0,1.1],blue);

% Close the postscript plot
%%close_plot(vps);


% - - - - - Evaluate the transmission at 700 eV vs Time - - - - -

(w1,w2)=linear_grid(hc/0.705, hc/0.695, 15);
lline = eval_fun(w1,w2);

% average it
lline = sum(lline)/length(lline);

oplot([hc/0.705,hc/0.695], lline*[1.,1.],orange);

% Average Tau at L lines:
print(-1.0*log( lline ));


% - - - Evaluate the C-K optical depth - - -

(w1,w2)=linear_grid(42.0, 44.0, 100);
ckhi = eval_fun(w1,w2);
%%cklo = min(ckhi);
%%ckhi = max(ckhi);
%%ckhi/cklo;

cklo = ckhi[0];
ckhi = ckhi[length(ckhi)-1];

% - - - - Save these values

ltaus = [ltaus, -1.0*log(lline)];
ltimes = [ltimes, ltime];
cktaus = [cktaus, log(ckhi/cklo)];

ltime = ltime + 0.2;


%  -----
  }

  while ( ltime < 2004.0 );


% - - - End of looping - - -

% Plot the predicted C-K depth and the ~700 eV depth

% Open a postscript output
%%%vps = open_plot("tau700_contamarfECS.ps/vcps", 1, 2);
vps = open_plot("tau700_twolevelECS.ps/vcps", 1, 2);
window(vps);

label("Time (years)","Tau",
	"Tau at 700 eV: the Two-level Model (green) and the Fit to ECS data (blue)");
xrange; yrange;      
xrange(1999.0, 2004.0);
yrange(0.0,0.7);

plot(ltimes[[1:]], ltaus[[1:]], green);

% scale the C-K values
%%oplot(ltimes, cktaus/10.0, red);


% And overplot the approx fit to ECS data:
variable diyr, tauinf, t0, tefold, t, model_tau;
diyr = 365.25;
tauinf = 0.589;
t0=1999.0+239.0/diyr;
tefold=761./diyr;

model_tau = tauinf * ( 1.0 - exp(-1.0*(ltimes-t0)/tefold) );

oplot(ltimes, model_tau, blue);

% And overplot Herman's C-K depth
%%oplot([1999.7, 2003.5],[0.62,2.33]/10.0, orange);

close_plot(vps);