% zo_1d_analysis.i 
%
% Read in the isis-format zo 1D histogram files
% from Capella, Mkn 421, and MARX capella simulation
% and scale and plot the distributions.
% 1/5/00 dd
%
% Generalize to study focus/LRF vs time...
% 1/12/00 dd

delete_data([1:6]);

% get_data variables:
variable dds_obs_out_dir, data_set_name;

% - - - - -
% more variables ...

% Simulation variables:
variable simm, simp, csimm, csimp;

% Variables to scale the observations, etc.
variable cap_fact, sim_fact, streak_fact;

% variables to keep average background count in
variable background;

% variables for the result of get_data
variable gdlo, gdhi, gdval, gderr;

% and pixel scales
variable meg_angpix, heg_angpix, use_angpix;
meg_angpix = 0.011124491;
heg_angpix = 0.0055625373;

% - - - - -

% Default selection of counts or flux here...
use_counts;


% Load some zero-order profiles:
%
% For comparison:
%  Capella 1318  zero-order
%  Capella 1318  streak
()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
	"obsid_1318_Sky_zo_1d_isis.dat");
()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
	"obsid_1318_Sky_strk_1d_isis.dat");

%******************************************
% Observation under measurement:
%  TBD test obs  zero-order
%  TBD test obs  streak
% Set cap_fact so that Observation and Capella
% zeroth-order wings agree
%******************************************
% Data sets:
% - - -
%%data_set_name = "Capella (obsid 1235)";
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_1235_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_1235_Sky_strk_1d_isis.dat");
%%cap_fact = 0.5;
% - - -
%%data_set_name = "HR1099 (obsid 62538)";
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_62538_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_62538_Sky_strk_1d_isis.dat");
%%cap_fact = 6.0;
% - - -
%%data_set_name = "HR1099 (obsid 1252)";
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_1252_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_1252_Sky_strk_1d_isis.dat");
%%cap_fact = 0.9;
% - - -
%%data_set_name = "SS 433 (obsid 106)";
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_106_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_106_Sky_strk_1d_isis.dat");
%%cap_fact = 0.8;
% - - -
%%data_set_name = "Capella (obsid 1103)";
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_1103_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_1103_Sky_strk_1d_isis.dat");
%%cap_fact = 1.4;
% - - -
%%data_set_name = "Capella (obsid 1318)";
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_1318_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_1318_Sky_strk_1d_isis.dat");
%%cap_fact = 1.01;
% - - -
%%data_set_name = "AB Dor (obsid 16)";
%%()=load_data("/spectra/d2/Fmeo/hak_16/Obs_out/" + 
%%	"obsid_16_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/Fmeo/hak_16/Obs_out/" + 
%%	"obsid_16_Sky_strk_1d_isis.dat");
%%cap_fact = 0.8;
% - - -
%%data_set_name = "Q 0836 (obsid 1450)";
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_1450_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_1450_Sky_strk_1d_isis.dat");
%%cap_fact = 0.5;
% - - -
%%data_set_name = "II Peg (obsid 1451)";
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_1451_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_1451_Sky_strk_1d_isis.dat");
%%cap_fact = 1.3;
% - - -
%%data_set_name = "4U 1636 (obsid 105)";
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_105_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_105_Sky_strk_1d_isis.dat");
%%cap_fact =3.0;
% - - -
%%data_set_name = "Cyg X-3 (obsid 1456)";
%%()=load_data("/spectra/d2/Fmeo/hak_1456/Obs_out/" + 
%%	"obsid_1456_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/Fmeo/hak_1456/Obs_out/" + 
%%	"obsid_1456_Sky_strk_1d_isis.dat");
%%cap_fact = 4.0;
% - - -
%%data_set_name = "Mkn 421 (obsid 457)";
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_457_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_457_Sky_strk_1d_isis.dat");
%%cap_fact = 1.0/3.8;
% - - -
%%data_set_name = "PKS 2149-306 (obsid 1481)";
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_1481_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_1481_Sky_strk_1d_isis.dat");
%%cap_fact = 0.2;
% - - -
%%data_set_name = "Lambda And (obsid 609)";
%%()=load_data("/spectra/d2/Fmeo/hak_609/Obs_out/" + 
%%	"obsid_609_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/Fmeo/hak_609/Obs_out/" + 
%%	"obsid_609_Sky_strk_1d_isis.dat");
%%cap_fact = 1.0;
% - - -
%%data_set_name = "Cyg X-3,II (obsid 1456b)";
%%()=load_data("/spectra/d2/Fmeo/hak_1456p2/Obs_out/" + 
%%	"obsid_1456_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/Fmeo/hak_1456p2/Obs_out/" + 
%%	"obsid_1456_Sky_strk_1d_isis.dat");
%%cap_fact =2.0;
% - - -
%%data_set_name = "3C 273 (obsid 459)";
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_459_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_459_Sky_strk_1d_isis.dat");
%%cap_fact = 2.0;
% - - -
%%data_set_name = "Cyg X-1 (obsid 1511)";
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_1511_Sky_zo_1d_isis.dat");
% no streak (all streak?!) because CC mode
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_1511_Sky_zo_1d_isis.dat");
%%cap_fact = 5.0;
% - - -
%%data_set_name = "UX Ari (obsid 605)";
%%()=load_data("/spectra/d2/Fmeo/hak_605/Obs_out/" + 
%%	"obsid_605_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/Fmeo/hak_605/Obs_out/" + 
%%	"obsid_605_Sky_strk_1d_isis.dat");
%%cap_fact = 0.8;
% - - -
%%data_set_name = "delta Ori (obsid 639)";
%%()=load_data("/spectra/d2/Fmeo/hak_639/Obs_out/" + 
%%	"obsid_639_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/Fmeo/hak_639/Obs_out/" + 
%%	"obsid_639_Sky_strk_1d_isis.dat");
%%cap_fact = 0.05;
% - - -
%%data_set_name = "NGC 3783 (obsid 373)";
%%()=load_data("/spectra/d2/Fmeo/hak_373/Obs_out/" + 
%%	"obsid_373_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/Fmeo/hak_373/Obs_out/" + 
%%	"obsid_373_Sky_strk_1d_isis.dat");
%%cap_fact = 1.0;
% - - -
%%data_set_name = "NGC 5548 (obsid 837)";
%%()=load_data("/spectra/d2/Fmeo/hak_837/Obs_out/" + 
%%	"obsid_837_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/Fmeo/hak_837/Obs_out/" + 
%%	"obsid_837_Sky_strk_1d_isis.dat");
%%cap_fact = 0.8;
% - - -
%%data_set_name = "LMC X-3 (obsid 92,N002)";
%%()=load_data("/spectra/d2/Fmeo/hak_92/Obs_out/" + 
%%	"obsid_92N002_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/Fmeo/hak_92/Obs_out/" + 
%%	"obsid_92N002_Sky_strk_1d_isis.dat");
%%cap_fact = 1.0;
% - - -
%%data_set_name = "Cir X-1 (obsid 706)";
% NOTE: exposures are 1.7s; streak is from 165 rows
% put the following into code below:
% For obsid 706 frame=1.7s; streak rows = 165:
%%streak_fact = 1700.0 / (41.0 * ((165.0)/1022.0)) ;
%
%%()=load_data("/spectra/d2/Fmeo/hak_706/Obs_out/" + 
%%	"obsid_706_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/Fmeo/hak_706/Obs_out/" + 
%%	"obsid_706_Sky_strk_1d_isis.dat");
%%cap_fact = 3.5;
% - - -
%%data_set_name = "Capella (obsid 57)";
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_57_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/GtoHAK/Obs_out/" + 
%%	"obsid_57_Sky_strk_1d_isis.dat");
%%cap_fact = 1.0;
% - - -
%%data_set_name = "Cen X-3 (obsid 705)";
%%()=load_data("/spectra/d2/Fmeo/hak_705/Obs_out/" + 
%%	"obsid_705_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/Fmeo/hak_705/Obs_out/" + 
%%	"obsid_705_Sky_strk_1d_isis.dat");
%%cap_fact = 1.0;
% - - -
%%data_set_name = "HD 68273 (obsid 629)";
%%()=load_data("/spectra/d2/Fmeo/hak_629/Obs_out/" + 
%%	"obsid_629_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/Fmeo/hak_629/Obs_out/" + 
%%	"obsid_629_Sky_strk_1d_isis.dat");
%%cap_fact = 0.3;
% - - -
data_set_name = "Zeta Pup (obsid 640)";
()=load_data("/spectra/d2/Fmeo/hak_640/Obs_out/" + 
	"obsid_640_Sky_zo_1d_isis.dat");
()=load_data("/spectra/d2/Fmeo/hak_640/Obs_out/" + 
	"obsid_640_Sky_strk_1d_isis.dat");
cap_fact = 0.3;
% - - -
%%data_set_name = "Algol (obsid 604)";
%%()=load_data("/spectra/d2/Fmeo/hak_604/Obs_out/" + 
%%	"obsid_604_Sky_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/Fmeo/hak_604/Obs_out/" + 
%%	"obsid_604_Sky_strk_1d_isis.dat");
%%cap_fact = 2.0;

%******************************************

% default is for cannonical set of four HETG spectra:
megm=1; megp=2; hegm=3; hegp=4;
% modify colors
cmegm=red; cmegp=orange; chegm=yellow; chegp=green;

% Also load simulation

% Simulation: aTDET w/background
%%()=load_data("/spectra/d2/GtoHAK/MARX_out/" + 
%%	"sim_capella_aTDET_zo_1d_isis.dat");
%%()=load_data("/spectra/d2/GtoHAK/MARX_out/" + 
%%	"sim_capella_aTDET_strk_1d_isis.dat");
% Set sim_fact so that sim_capella and Capella STREAKs agree
% 1/0.66 factor from time apodization in sim aTDET analysis
%%sim_fact = 1.0/0.66;


% Simulation: Sky, no background
()=load_data("/spectra/d2/GtoHAK/MARX_out/" + 
	"Capella_Sky0.35_Sky_zo_1d_isis.dat");
()=load_data("/spectra/d2/GtoHAK/MARX_out/" + 
	"Capella_Sky0.35_Sky_strk_1d_isis.dat");
% Set sim_fact so that sim_capella and Capella STREAKs agree
% Factor of 1.0 should be close ...
sim_fact = 1.0;

% define simulation indices and colors
simm=5; simp=6;
csimm=purple; csimp=mediblue;

%--------------- remove the average background level --------------------

% Loop over data sets!
message("");
fdata=1;
while ( fdata < 7 ) {
  rsout = region_sum(fdata, 1.0,50.5);
  background = rsout.sum;
  rsout = region_sum(fdata, 149.5,199.0);
  background = rsout.sum + background;
  % get the data set
  (gdlo, gdhi, gdval, gderr) = get_data(fdata);
  % the sum was over 1/2 of the bins
  background = background/(0.5*length(gdval));
  message("hist index "+string(fdata)+" background cnts/bin = "+
	string(background));
  gdval = gdval - background;
  put_data(fdata, gdlo, gdhi, gdval, gderr);

  fdata = fdata + 1;
  }
message("");


%--------------- set the data colors --------------------
set_data_color(megm,cmegm);
set_data_color(megp,cmegp);
set_data_color(hegm,chegm);
set_data_color(hegp,chegp);
set_data_color(simm,csimm);
set_data_color(simp,csimp);

%
% Process the 1D profiles (in pixels) to be offset to
% realistic Angstroms...
%

use_angpix = meg_angpix;


% Normalize the streak to the same scale as the zeroth-order
% Ratio of zeroth-order exposure time to streak exposure time,
% include fact that the HAK streak excludes ~240 rows:
streak_fact = 3200.0 / (41.0 * ((1022.-240.0)/1022.0)) ;
% For obsid 706 frame=1.7s; streak rows = 165:
%%streak_fact = 1700.0 / (41.0 * ((165.0)/1022.0)) ;

% Capella
% Zero-order
(gdlo, gdhi, gdval, gderr) = get_data(megm);
gdlo = 15.0+use_angpix*(gdlo-100.0);
gdhi = 15.0+use_angpix*(gdhi-100.0);
gdval = cap_fact * gdval;
gderr = cap_fact * gderr;
put_data(megm, gdlo, gdhi, gdval, gderr);
% Streak, scale the counts
(gdlo, gdhi, gdval, gderr) = get_data(megp);
gdlo = 15.0+use_angpix*(gdlo-100.0);
gdhi = 15.0+use_angpix*(gdhi-100.0);
gdval = cap_fact * gdval * streak_fact;
gderr = cap_fact * gderr * streak_fact;
put_data(megp, gdlo, gdhi, gdval, gderr);

% Capella SIMULATION
% Zero-order
(gdlo, gdhi, gdval, gderr) = get_data(simm);
gdlo = 15.0+use_angpix*(gdlo-100.0);
gdhi = 15.0+use_angpix*(gdhi-100.0);
gdval = cap_fact * sim_fact * gdval;
gderr = cap_fact * sim_fact * gderr;
put_data(simm, gdlo, gdhi, gdval, gderr);
% Streak, scale the counts
(gdlo, gdhi, gdval, gderr) = get_data(simp);
gdlo = 15.0+use_angpix*(gdlo-100.0);
gdhi = 15.0+use_angpix*(gdhi-100.0);
gdval = cap_fact * sim_fact * gdval * streak_fact;
gderr = cap_fact * sim_fact * gderr * streak_fact;
put_data(simp, gdlo, gdhi, gdval, gderr);

% Obs being measured
% Zero-order
(gdlo, gdhi, gdval, gderr) = get_data(hegm);
gdlo = 15.0+use_angpix*(gdlo-100.0);
gdhi = 15.0+use_angpix*(gdhi-100.0);
gdval = gdval;
gderr = gderr;
put_data(hegm, gdlo, gdhi, gdval, gderr);
% Streak, scale the counts
(gdlo, gdhi, gdval, gderr) = get_data(hegp);
gdlo = 15.0+use_angpix*(gdlo-100.0);
gdhi = 15.0+use_angpix*(gdhi-100.0);
gdval = gdval * streak_fact;
gderr = gderr * streak_fact;
put_data(hegp, gdlo, gdhi, gdval, gderr);


% Notice only the core of the streak histograms
fdata=2;
while ( fdata < 7 ) {
  xnotice(fdata, 15.0-2.0*use_angpix, 15.0+2.0*use_angpix);
  fdata = fdata + 2;
  }


% Plot the data...
LABEL_BY_DEFAULT=0;
window(p1);
plot_bin_density;
ylog;
yrange(2.0, 7.5);
xrange(15.0-10.0*use_angpix, 15.0+10.0*use_angpix);
color(white);

% Capella
plot_data(megp); 
oplot_data(megm);
color(cmegp);
xylabel(14.92, 6.4, "Capella Streak");
color(cmegm);
xylabel(14.92, 6.2, "Capella Zeroth-order");
color(white);
xylabel(15.04, 6.4, "streak factor = "+string(streak_fact) );

% Observation
oplot_data(hegm);
color(chegm);
xylabel(14.92, 5.8, data_set_name+" Zeroth-order");
% Observation streak
oplot_data(hegp);
color(chegp);
xylabel(14.92, 6.0, data_set_name+" Streak");
%
color(white);
xylabel(15.04, 6.2, "Cap factor = "+string(cap_fact) );


data_xtitle="Pseudo Wavelength (" +
	string(use_angpix) + " A per pixel)";

label(data_xtitle,"log Counts/Pseudo-A",
	data_set_name+": Zeroth Order 1D Projections");

% simulation
while (0 > 1) {
plot_data(simm);
oplot_data(simp);
color(white);
xylabel(15.04, 6.0, "simulation factor = "+string(sim_fact) );
color(csimp);
xylabel(14.985, 3.6, "Sim_1318 Streak");
color(csimm);
xylabel(14.985, 3.4, "Sim_1318 Zeroth-order");
color(white);
xylabel(15.04, 6.4, "streak factor = "+string(streak_fact) );
}

%
% Fit the observation's zo and streak
%
% Do it by manually pasting these lines...
%

while (0 > 1) {

% Zeroth order
fdata = hegm;
cfdata = chegm;
()=evalfile("fit_line.i");

% Streak
fdata = hegp;
cfdata = chegp;
()=evalfile("fit_line.i");


}

LABEL_BY_DEFAULT=1;