%; 2001.06.30 dph
%; plot letgs geom test results.
%;
%;          This is a "cut-and-paste" file, not an executable script!
%;
%; file "-1" is old geom: telD1999-07-23geomN0003.fits
%; file "-2" is new geom: telD1999-07-23geomN0004B.fits

sldir=getenv("HOME")+ "/libisis/";
()=evalfile(sldir+"gsmooth.sl");
()=evalfile(sldir+"presid.sl");

fn1 = "./Pha/hrcf01248-1_pha2.fits.gz";
fn2 = "./Pha/hrcf01248-2_pha2.fits.gz";

load_data(fn1);
load_data(fn2);

p1=open_plot("/xwin");resize(18);
label("Wavelength [A]", "Flux [counts/bin]", "Capella (ObsID 1248):  Test geom file");
plot_bin_integral;

cm1=get_data_counts(1);
cp1=get_data_counts(2);

ym1=cm1.value;
yp1=cp1.value;
x1=cm1.bin_lo;
x2=cm1.bin_hi;
(sx1,sx2,sym1)=gsmooth(x1,x2,ym1,0.0125);
(sx1,sx2,syp1)=gsmooth(x1,x2,yp1,0.0125);

cm2=get_data_counts(3);
cp2=get_data_counts(4);
ym2=cm2.value;
yp2=cp2.value;
(sx1,sx2,sym2)=gsmooth(x1,x2,ym2,0.0125);
(sx1,sx2,syp2)=gsmooth(x1,x2,yp2,0.0125);

plasma(aped);       % load for line-id's.
load_model("Tbl/Capella_1T.tbl"); % LogT=6.8 - capella max.
f=model_spectrum(x1,x2);        % evaluate a model so "brightest" works.

%%%%%%%%%%%%%%%% SETUP is done.  Make plots...................

% for output to a file, select one of these:
%+
psid=dup_plot("Cmp_geom_letgs.ps/cps",p1);resize(18,0.6);

% make one as gif for html doc:
gifid=dup_plot("Cmp_geom_letgs.gif/gif",p1);resize(18,0.6);
%-
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

xmin=12; xmax=20; ymin=-50; ymax=1400; nlines=20;
xrange(xmin,xmax); yrange(ymin,ymax);
hplot(x1,x2,ym1,2);
ohplot(x1,x2,ym2,3);
ohplot(x1,x2,ym2-ym1,4);
plot_group(brightest(nlines,where(wl(xmin,xmax))),1);
color(1);xylabel(xmin+1,ymax-(ymax-ymin)/10.,"Minus orders");

close_plot(gifid);window(p1);              %. optional


hplot(x1,x2,yp1,2);
ohplot(x1,x2,yp2,3);
ohplot(x1,x2,yp2-yp1,4);
plot_group(brightest(nlines,where(wl(xmin,xmax))),1);
color(1);xylabel(xmin+1,ymax-(ymax-ymin)/10.,"Minus orders");


xmin=103;xmax=110;ymin=-10;ymax=250;nlines=2;
nxrange(xmin,xmax);yrange(ymin,ymax);
hplot(sx1,sx2,sym1,2);
ohplot(sx1,sx2,sym2,3);
ohplot(sx1,sx2,sym2-sym1,4);
plot_group(brightest(2,where(wl(xmin,xmax))),1);
color(1);xylabel(xmin+1,ymax-(ymax-ymin)/10.,"Minus orders");

hplot(sx1,sx2,syp1,2);
ohplot(sx1,sx2,syp2,3);
ohplot(sx1,sx2,syp2-syp1,4);
plot_group(brightest(2,where(wl(xmin,xmax))),1);
color(1);xylabel(xmin+1,ymax-(ymax-ymin)/10.,"Plus orders");

close_plot(psid);window(p1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%
% while we're at it, compare +- in file -2;


p2=open_plot("/xwin",1,2,[3,1]);resize(25,0.4);
label("Wavelength [A]", "Flux [counts/bin]", "Capella (ObsID 1248):  Compare +/-");

psid=dup_plot("Cmp_pm_letgs.ps/cps",p2);resize(25,0.4);    %. optional

_pgpage();
xmin=6.5; xmax=8.5; ymin=-0; ymax=300; nlines=10;
presid(p2,[xmin,xmax,ymin,ymax],x1,x2,sym2,syp2);
pane(1);plot_group(brightest(nlines,where(wl(xmin,xmax))),1);


_pgpage();
xmin=12; xmax=15.5; ymin=0; ymax=800; nlines=10;
presid(p2,[xmin,xmax,ymin,ymax],x1,x2,sym2,syp2);
pane(1);plot_group(brightest(nlines,where(wl(xmin,xmax))),1);


_pgpage();
xmin=15.5; xmax=20; ymin=0; ymax=1200; nlines=10;
presid(p2,[xmin,xmax,ymin,ymax],x1,x2,sym2,syp2);
pane(1);
plot_group(brightest(nlines,where(wl(xmin,xmax))),1);

_pgpage();
xmin=21.4; xmax=25.4; ymin=0; ymax=250; nlines=6;
presid(p2,[xmin,xmax,ymin,ymax],x1,x2,sym2,syp2);
pane(1);plot_group(brightest(nlines,where(wl(xmin,xmax))),1);

_pgpage();
xmin=25; xmax=30; ymin=0; ymax=60; nlines=6;
presid(p2,[xmin,xmax,ymin,ymax],x1,x2,sym2,syp2);
pane(1);plot_group(brightest(nlines,where(wl(xmin,xmax))),1);
	
_pgpage();
xmin=40; xmax=45; ymin=0; ymax=120; nlines=12;
presid(p2,[xmin,xmax,ymin,ymax],x1,x2,sym2,syp2);
pane(1);plot_group(brightest(nlines,where(wl(xmin,xmax))),1);

_pgpage();
xmin=45; xmax=55; ymin=0; ymax=120; nlines=12;
presid(p2,[xmin,xmax,ymin,ymax],x1,x2,sym2,syp2);
pane(1);plot_group(brightest(nlines,where(wl(xmin,xmax))),1);

_pgpage();
xmin=55; xmax=75; ymin=0; ymax=120; nlines=12;
presid(p2,[xmin,xmax,ymin,ymax],x1,x2,sym2,syp2);
pane(1);plot_group(brightest(nlines,where(wl(xmin,xmax))),1);

_pgpage();
xmin=93; xmax=95; ymin=0; ymax=320; nlines=4;
presid(p2,[xmin,xmax,ymin,ymax],x1,x2,sym2,syp2);
pane(1);plot_group(brightest(nlines,where(wl(xmin,xmax))),1);

_pgpage();
xmin=103.5; xmax=109; ymin=0; ymax=150; nlines=4;
presid(p2,[xmin,xmax,ymin,ymax],x1,x2,sym2,syp2);
pane(1);plot_group(brightest(nlines,where(wl(xmin,xmax))),1);

_pgpage();
xmin=116; xmax=123; ymin=0; ymax=100; nlines=8;
presid(p2,[xmin,xmax,ymin,ymax],x1,x2,sym2,syp2);
pane(1);plot_group(brightest(nlines,where(wl(xmin,xmax))),1);

_pgpage();
xmin=128; xmax=134; ymin=0; ymax=100; nlines=8;
presid(p2,[xmin,xmax,ymin,ymax],x1,x2,sym2,syp2);
pane(1);plot_group(brightest(nlines,where(wl(xmin,xmax))),1);



gifid=dup_plot("Cmp_pm_letgs.gif/gif",p2);resize(25,0.4);  %. optional

_pgpage();
xmin=5; xmax=140; ymin=0; ymax=1000; nlines=52;
presid(p2,[xmin,xmax,ymin,ymax],x1,x2,sym2,syp2);
pane(1);
plot_group(brightest(nlines,where(wl(xmin,xmax) and el_ion([C,N,O,Mg,Si,S,Ni,Ne]))),1);
s=line_label_default_style();
s.top_frac=0.4;s.bottom_frac=0.3;
plot_group(brightest(nlines,where(wl(xmin,xmax) and el_ion(Fe))),4,s);

close_plot(gifid);window(p2);    %. optional

close_plot(psid);window(p2);    %. optional