#source en minuscule
set source $1
set srcMaj $2
set date   $3
set fcycle $4
set lcycle $5
set nframes $6
maxcycle = ${lcycle}+1

set path "/sps/baoradio/AmasNancay/JEC/${srcMaj}"

#first script executed
clearscript image
clearscript doImage

#mergeimg should be executed before timeevol, fwhnInt
clearscript mergeimg
clearscript timeevol
clearscript fwhmInt

#script independant from mergeimg
clearscript timeevol2
#######################
defscript image
#
#compute (ON-OFF)/OFFn with OFF computed from the first 30s of each image
set cycle  $1
openppf ${path}/img_${date}_${source}_cycle${cycle}.ppf 
del ich0
del zCh0
objaoper img${cycle} slicexz 0 ich0
c++exec TMatrix<r_4>zCh0(ich0(Range(0),Range::all()),false); \
for (sa_size_t i=1;i<50;i++){zCh0+=ich0(Range(i),Range::all());} \
zCh0/=50.; KeepObj(zCh0);\
for (sa_size_t i=0;i<ich0.NRows();i++){ich0(Range(i),Range::all())-=zCh0;ich0(Range(i),Range::all()).Div(zCh0);}
#
rename ich0 ich${cycle}0
rename zCh0 tsys${cycle}0

del ich1
del zCh1
objaoper img${cycle} slicexz 1 ich1
c++exec TMatrix<r_4>zCh1(ich1(Range(0),Range::all()),false); \
for (sa_size_t i=1;i<50;i++){zCh1+=ich1(Range(i),Range::all());} \
zCh1/=50.; KeepObj(zCh1);\
for (sa_size_t i=0;i<ich1.NRows();i++){ich1(Range(i),Range::all())-=zCh1;ich1(Range(i),Range::all()).Div(zCh1);}
#
rename ich1 ich${cycle}1
rename zCh1 tsys${cycle}1

endscript
#######################
###################
defscript doImage

for ic ${fcycle}:${maxcycle}
echo "process ${ic}"
image $ic
end
endscript


#######################
defscript mergeimg
#
del ich0Tot
del ich1Tot
del ncycles
c++exec TMatrix<r_4>ich0Tot(ich${fcycle}0(Range(0,${nframes}),Range::all()),false); \
	TMatrix<r_4>ich1Tot(ich${fcycle}1(Range(0,${nframes}),Range::all()),false); \
	TVector<r_4>ncycles(1); ncycles(0)=1; \
	KeepObj(ich0Tot); KeepObj(ich1Tot); KeepObj(ncycles);

fp1cycle = ${fcycle}+1

for ic ${fp1cycle}:${maxcycle}
  c++exec ich0Tot+=ich${ic}0(Range(0,${nframes}),Range::all()); \
  	  ich1Tot+=ich${ic}1(Range(0,${nframes}),Range::all()); \
  	  ncycles(0)+=1;
end

c++exec ich0Tot/=ncycles(0); ich1Tot/=ncycles(0);

setaxesatt 'font=helvetica,bold,14 fixedfontsize'
newwin 
disp ich0Tot 'lut=lin,-0.1,3.5 colbr128 showcmap=top'

newwin 
disp ich1Tot 'lut=lin,-0.1,3.5 colbr128 showcmap=top'
endscript
#######################
defscript timeevol
#
#intensite en ft du temps avec image cumulee
del rg0Tot
del rg1Tot
#1405,1415 MHz => range [5080,5405]
c++exec TMatrix<r_4>tmp0(ich0Tot(Range::all(),Range(5080,5405)),false); \
	TVector<r_4>rg0Tot(ich0Tot.NRows()); \
	for(sa_size_t i=0;i<tmp0.NRows();++i){ \
	  double mean,sigma; \
	  MeanSigma(tmp0(Range(i),Range::all()),mean,sigma); \
	  rg0Tot(i) = mean; \
	} \
	KeepObj(rg0Tot); \
	TMatrix<r_4>tmp1(ich1Tot(Range::all(),Range(5080,5405)),false); \
	TVector<r_4>rg1Tot(ich1Tot.NRows()); \
	for(sa_size_t i=0;i<tmp1.NRows();++i){ \
	  double mean,sigma; \
	  MeanSigma(tmp1(Range(i),Range::all()),mean,sigma); \
	  rg1Tot(i) = mean; \
	} \
	KeepObj(rg1Tot);

newwin
setaxesatt 'grid'
graphicatt '' 
#90sec = milieu du cycle = maximum
#donne 7.56kHz
#273 canaux en 170sec => 8.22kHz
plot2d rg0Tot (n*170)/273 val 1 'blue cpts notit nsta'
plot2d rg1Tot (n*170)/273 val 1 'same red cpts notit nsta'
settitle "${srcMaj} Drift Scan [1405,1415]MHz  all cycles Ch 0 (blue) Ch 1 (red)"
setaxelabels 't(sec)' 'I (a.u)'
endscript

#######################
defscript fwhmInt
#
#spectre de frequence canal en temps [120,146] FWHM.
#
del s0
del s1
c++exec TMatrix<r_4>tmp0(ich0Tot(Range(120,146),Range::all()),false); \
	TVector<r_4>s0(ich0Tot.NCols()); \
	for(sa_size_t i=0;i<tmp0.NCols();++i){ \
	  double mean,sigma; \
	  MeanSigma(tmp0(Range::all(),Range(i)),mean,sigma); \
	  s0(i) = mean; \
	} \
	KeepObj(s0); \
	TMatrix<r_4>tmp1(ich1Tot(Range(120,146),Range::all()),false); \
	TVector<r_4>s1(ich1Tot.NCols()); \
	for(sa_size_t i=0;i<tmp1.NCols();++i){ \
	  double mean,sigma; \
	  MeanSigma(tmp1(Range::all(),Range(i)),mean,sigma); \
	  s1(i) = mean; \
	} \
	KeepObj(s1);
newwin
setaxesatt 'grid'
graphicatt 'xylimits=1250,1500,0,2' 
plot2d s0 (n/8192)*250+1250 val n>0 'blue cpts notit nsta'
plot2d s1 (n/8192)*250+1250 val 1 'same red cpts notit nsta'
settitle "${srcMaj} Drift Scan all cycles Ch 0 (blue) Ch 1 (red)"
setaxelabels 'Freq. (MHz)' 'I (a.u)'
endscript
#######################
#######################
defscript timeevol2

#append all images 
del ch0evol
del ch1evol
c++exec sa_size_t nrows=(${nframes}+1)*(${lcycle}-${fcycle}+1); \
	sa_size_t ncols=ich${fcycle}0.NCols(); \
	TMatrix<r_4>ch0evol(nrows,ncols); TMatrix<r_4>ch1evol(nrows,ncols); \
	KeepObj(ch0evol); KeepObj(ch1evol);

for ic ${fcycle}:${maxcycle}
  c++exec sa_size_t curRow=(${ic}-${fcycle})*(${nframes}+1); sa_size_t nextm1CurRow=curRow+(${nframes}+1)-1; \
	  cout << "ic, curRow, next : "<<${ic}<<" "<<curRow<<" "<<nextm1CurRow<<endl; \
	  ch0evol(Range(curRow,nextm1CurRow),Range::all())=ich${ic}0(Range(0,${nframes}),Range::all()); \
	  ch1evol(Range(curRow,nextm1CurRow),Range::all())=ich${ic}1(Range(0,${nframes}),Range::all());
end

#intensite en ft du temps avec image cumulee
del i0evol
del i1evol
#1405,1415 MHz => range [5080,5405]
c++exec TMatrix<r_4>tmp0(ch0evol(Range::all(),Range(5080,5405)),false); \
	TVector<r_4>i0evol(ch0evol.NRows()); \
	for(sa_size_t i=0;i<tmp0.NRows();++i){ \
	  double mean,sigma; \
	  MeanSigma(tmp0(Range(i),Range::all()),mean,sigma); \
	  i0evol(i) = mean; \
	} \
	KeepObj(i0evol); \
	TMatrix<r_4>tmp1(ch1evol(Range::all(),Range(5080,5405)),false); \
	TVector<r_4>i1evol(ch1evol.NRows()); \
	for(sa_size_t i=0;i<tmp1.NRows();++i){ \
	  double mean,sigma; \
	  MeanSigma(tmp1(Range(i),Range::all()),mean,sigma); \
	  i1evol(i) = mean; \
	} \
	KeepObj(i1evol);



newwin
setaxesatt 'grid'
graphicatt '' 
plot2d i0evol n  val 1 'blue cpts notit nsta'
plot2d i1evol n  val 1 'same red cpts notit nsta'
settitle "${srcMaj} Drift Scan [1405,1415]MHz all cycles Ch 0 (blue) Ch 1 (red)"
setaxelabels 't(a.u)' '(I-OFF)/OFF (a.u)'


#Tsys en ft du temps
c++exec sa_size_t nrows=(${lcycle}-${fcycle}+1); \
	sa_size_t ncols=ich${fcycle}0.NCols(); \
	TMatrix<r_4>Tsys0evol(nrows,ncols); TMatrix<r_4>Tsys1evol(nrows,ncols); \
	KeepObj(Tsys0evol); KeepObj(Tsys1evol);

for ic ${fcycle}:${maxcycle}
  c++exec sa_size_t curRow=(${ic}-${fcycle}); \
	  cout << "ic, curRow, next : "<<${ic}<<" "<<curRow<<endl; \
	  Tsys0evol(Range(curRow),Range::all())=tsys${ic}0(Range(0),Range::all()); \
	  Tsys1evol(Range(curRow),Range::all())=tsys${ic}1(Range(0),Range::all());
end

del T0evol
del T1evol
#1405,1415 MHz => range [5080,5405]
c++exec TMatrix<r_4>tmp0(Tsys0evol(Range::all(),Range(5080,5405)),false); \
	TVector<r_4>T0evol(Tsys0evol.NRows()); \
	for(sa_size_t i=0;i<tmp0.NRows();++i){ \
	  double mean,sigma; \
	  MeanSigma(tmp0(Range(i),Range::all()),mean,sigma); \
	  T0evol(i) = mean; \
	} \
	KeepObj(T0evol); \
	TMatrix<r_4>tmp1(Tsys1evol(Range::all(),Range(5080,5405)),false); \
	TVector<r_4>T1evol(Tsys1evol.NRows()); \
	for(sa_size_t i=0;i<tmp1.NRows();++i){ \
	  double mean,sigma; \
	  MeanSigma(tmp1(Range(i),Range::all()),mean,sigma); \
	  T1evol(i) = mean; \
	} \
	KeepObj(T1evol);

newwin
setaxesatt 'grid'
graphicatt '' 
plot2d T0evol n+${fcycle}  val 1 'blue cpts notit nsta'
plot2d T1evol n+${fcycle}  val 1 'same red cpts notit nsta'
settitle "${srcMaj} Drift Scan [1405,1415]MHz all cycles Ch 0 (blue) Ch 1 (red)"
setaxelabels 'cycle num.' 'OFF/Gain (a.u)'

	
endscript