set user     "AST"
set source   "Abell1205"
set srclower "abell1205"
set date     "20110415"
set cycle    "8"
set mode     "Off"
set toppath  "/sps/baoradio/AmasNancay/${user}"
set path     "${toppath}/${source}/dab"
set ppf      "0"

# Read DAB calibration in Jansky (Voie E / PKS1127-14)
#-----------------------------------------------------
n/read nt_dabjy ${toppath}/DAB_Jy.dat " " dabfreq,0 dabjy,1
newwin 1 1
n/pl nt_dabjy.dabjy%dabfreq 1 ! "marker=star,4 magenta notit"
settitle "DAB Voie E / PKS1127-14. Juin 2011"


# Read DAB calibration in A.U. (from BAORadio data)
#--------------------------------------------------
openppf ${path}/dabchrom_${date}_${srclower}_${mode}.ppf
ls

## Baseline
del b0 b1 b0f b1f
objaoper /home/base row 0 b0
objaoper /home/base row 1 b1
objaoper /home/basefilt row 0 b0f
objaoper /home/basefilt row 1 b1f

## First DAB
del dab01 dab11 dab01f dab11f dab01n dab11n
objaoper /home/firstdab row 0 dab01
objaoper /home/firstdab row 1 dab11
objaoper /home/firstdabfilt row 0 dab01f
objaoper /home/firstdabfilt row 1 dab11f
objaoper /home/firstdelta row 0 dab01d
objaoper /home/firstdelta row 1 dab11d

## Second DAB
del dab02 dab12 dab02f dab12f dab02n dab12n
objaoper /home/seconddab row 0 dab02
objaoper /home/seconddab row 1 dab12
objaoper /home/seconddabfilt row 0 dab02f
objaoper /home/seconddabfilt row 1 dab12f
objaoper /home/seconddelta row 0 dab02d
objaoper /home/seconddelta row 1 dab12d

## Mean(delta(first,second DAB)) per channel
c++exec TVector<r_4> meandab0d(dab01d,false); KeepObj(meandab0d); \
        TVector<r_4> meandab1d(dab11d,false); KeepObj(meandab1d); \
        meandab0d+=dab02d; meandab0d/=2.; \
        meandab1d+=dab12d; meandab1d/=2.;

## Rebin BAO and calibrated data in 0.5Hz/bin to perform division
exptovec vf nt_dabjy dabfreq
exptovec vjy nt_dabjy dabjy
c++exec TVector<r_4> dumm0(vf), dumm1(vjy);\
        TVector<r_4> vcalfreq=dumm0.SubVector(Range(320,854)).TransposeSelf();\
        KeepObj(vcalfreq);\
        TVector<r_4> vcaljy=dumm1.SubVector(Range(320,854)).TransposeSelf();\
        KeepObj(vcaljy);

del invec0 invec1 invec2 invec3 val0 val1
cp meandab0d invec0
cp meandab1d invec1
cp vcaljy invec2
cp vcalfreq invec3
ls

c++compile interpol
c++link interpol.so dointerpol
call dointerpol

c++exec TVector<r_8> calcoeff0f1400a1420=calcoeff0.SubVector(Range(300,340)); \
	KeepObj(calcoeff0f1400a1420);\
        TVector<r_8> calcoeff1f1400a1420=calcoeff1.SubVector(Range(300,340)); \
	KeepObj(calcoeff1f1400a1420);
projh1d h0 calcoeff0f1400a1420 val
projh1d h1 calcoeff1f1400a1420 val

## Save Mean(delta(first,second DAB)) per channel into a ppf file
#if ( ppf == "1" ) then
#  del daboff
#  c++exec TMatrix<r_4> daboff(2,dab01.NElts()); KeepObj(daboff); \
#          daboff.Row(0)=meandab0d; daboff.Row(1)=meandab1d;
#  saveppf daboff "${toppath}/${source}/daboff_${date}_${srclower}_calibcycle8.fits"
#endif 

# Plot results
#--------------
graphicatt "xylimits=1250,1500,0,3"

## Baseline
newwin 1 1
plot2d b0 (n/8192)*250+1250 val n>0 "blue cpts nstat notit"
plot2d b1 (n/8192)*250+1250 val n>0 "same red cpts nstat notit"
settitle "Baseline raw $source $date $mode Ch0 blue Ch1 red"

newwin 1 1
plot2d b0f (n/8192)*250+1250 val n>0 "blue cpts nstat notit"
plot2d b1f (n/8192)*250+1250 val n>0 "same red cpts nstat notit"
settitle "Baseline filt $source $date $mode Ch0 blue Ch1 red"

## First DAB
#newwin 1 1
#plot2d dab01 (n/8192)*250+1250 val n>0 "blue cpts nstat notit"
#plot2d dab11 (n/8192)*250+1250 val n>0 "same red cpts nstat notit"
#settitle "First DAB raw $source $date $mode Ch0 blue Ch1 red"

newwin 1 1
plot2d dab01f (n/8192)*250+1250 val n>0 "blue cpts nstat notit"
plot2d dab11f (n/8192)*250+1250 val n>0 "same red cpts nstat notit"
settitle "First DAB filt $source $date $mode Ch0 blue Ch1 red"

## Second DAB
#newwin 1 1
#plot2d dab02 (n/8192)*250+1250 val n>0 "blue cpts nstat notit"
#plot2d dab12 (n/8192)*250+1250 val n>0 "same red cpts nstat notit"
#settitle "Second DAB raw $source $date $mode Ch0 blue Ch1 red"

newwin 1 1
plot2d dab02f (n/8192)*250+1250 val n>0 "blue cpts nstat notit"
plot2d dab12f (n/8192)*250+1250 val n>0 "same red cpts nstat notit"
settitle "Second DAB filt $source $date $mode Ch0 blue Ch1 red"

## Delta=DAB-baseline + Calibrated DAB in Jansky
graphicatt "xylimits=1100,1600,0,30"
newwin 1 1
n/pl nt_dabjy.dabjy%dabfreq 1 ! "marker=star,4 turquoise notit nstat"
plot2d dab01d (n/8192)*250+1250 val*18 n>0 "same blue cpts nstat notit"
plot2d dab11d (n/8192)*250+1250 val*22 n>0 "same red cpts nstat notit"
settitle "First DAB delta $source $date $mode"
setaxelabels "Freq (MHz)" "Intensity" "bigboldfont"
addline 1420 4 1430 4 "turquoise line=solid,3"
addline 1420 3 1430 3 "blue line=solid,3"
addline 1420 2 1430 2 "red line=solid,3"
addtext 1440 4 "DAB Voie E (Jy)" "bigfont"
addtext 1440 3 "Ch0*18 (A.U.)" "bigfont"
addtext 1440 2 "Ch1*22 (A.U.)" "bigfont"

newwin 1 1
n/pl nt_dabjy.dabjy%dabfreq 1 ! "marker=star,4 turquoise notit nstat"
plot2d dab02d (n/8192)*250+1250 val*18 n>0 "same blue cpts nstat notit"
plot2d dab12d (n/8192)*250+1250 val*22 n>0 "same red cpts nstat notit"
settitle "Second DAB delta $source $date $mode"
setaxelabels "Freq (MHz)" "Intensity" "bigboldfont"
addline 1420 4 1430 4 "turquoise line=solid,3"
addline 1420 3 1430 3 "blue line=solid,3"
addline 1420 2 1430 2 "red line=solid,3"
addtext 1440 4 "DAB Voie E (Jy)" "bigfont"
addtext 1440 3 "Ch0*18 (A.U.)" "bigfont"
addtext 1440 2 "Ch1*22 (A.U.)" "bigfont"

## Delta=Mean(DAB-baseline) + Calibrated DAB in Jansky
graphicatt "xylimits=1100,1600,0,30"
newwin 1 1
n/pl nt_dabjy.dabjy%dabfreq 1 ! "marker=star,4 turquoise notit nstat"
plot2d meandab0d (n/8192)*250+1250 val*18 n>0 "same blue cpts nstat notit"
plot2d meandab1d (n/8192)*250+1250 val*23 n>0 "same red cpts nstat notit"
settitle "Mean DAB delta $source $date $mode"
setaxelabels "Freq (MHz)" "Intensity" "bigboldfont"
addline 1420 4 1430 4 "turquoise line=solid,3"
addline 1420 3 1430 3 "blue line=solid,3"
addline 1420 2 1430 2 "red line=solid,3"
addtext 1440 4 "DAB Voie E (Jy)" "bigfont"
addtext 1440 3 "Ch0*18 (A.U.)" "bigfont"
addtext 1440 2 "Ch1*23 (A.U.)" "bigfont"

graphicatt ""
newwin 1 1
plot2d meandab0d (n/8192)*250+1250 val n>0 "blue cpts nstat notit"
plot2d meandab1d (n/8192)*250+1250 val n>0 "same red cpts nstat notit"
settitle "Raw MeanDAB Ch0(blue) Ch1(red)"

## Rebinned DAB spectra
graphicatt ""
newwin 1 1
plot2d val0 (n/500)*250+1250 val 1 "cpts blue notit nstat"
plot2d val1 (n/500)*250+1250 val 1 "cpts same red notit nstat"
settitle "BAO DAB spectra Ch0 (blue) Ch1(red)"
setaxelabels "Freq. (MHz)" "I (a.u.)" "bigboldfont"

newwin 1 1
plot2d val2 (n/500)*250+1250 val 1 "cpts blue notit nstat"
settitle "Calibrated DAB spectra"
setaxelabels "Freq. (MHz)" "I (Jy)" "bigboldfont"

## Calibration coefficients Jy/a.u. in the whole bandwidth
newwin 1 1
plot2d calcoeff1 (n/500)*250+1250 val 1 "cpts red notit nstat"
plot2d calcoeff0 (n/500)*250+1250 val 1 "cpts same blue notit nstat"
settitle "Calibration coefficients spectrum Ch0(blue) Ch1(red)"
setaxelabels "Freq. (MHz)" "Coeff Jy/a.u."  "bigboldfont"

## Calibration coefficients Jy/a.u. in [1400,1420]MHz
newwin 1 1
disp h0 "notit"
settitle "CalCoeff (Jy/a.u.) Ch0"
newwin 1 1
disp h1 "notit"
settitle "CalCoeff (Jy/a.u.) Ch1"