source: BAORadio/AmasNancay/trunk/imgdrift.pic @ 638

#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}"

clearscript image
#######################
defscript image
#
set cycle  $1
openppf ${path}/img_${date}_${source}_cycle${cycle}.ppf 
del ich0
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.; \
for (sa_size_t i=0;i<ich0.NRows();i++){ich0(Range(i),Range::all())-=zCh0;}
#
rename ich0 ich${cycle}0

del ich1
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.; \
for (sa_size_t i=0;i<ich1.NRows();i++){ich1(Range(i),Range::all())-=zCh1;}
#
rename ich1 ich${cycle}1

endscript
#######################

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

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'

#intensite en ft du temps
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 2cycles Ch 0 (blue) Ch 1 (red)"
setaxelabels 't(sec)' 'I (a.u)'

#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 2cycles Ch 0 (blue) Ch 1 (red)"
setaxelabels 'Freq. (MHz)' 'I (a.u)'