# Script pour lire un fichier FITS avec les donnees FITS FFT_SHRINK et faire des plots de 
# controle
#  Appel : ckshrink.pic Directory NumeroFITS 
delobjs *
set DIR $1
set OSIG signal$2
set file $DIR/$OSIG.fits
echo '----> ckshrink: Opening file ' $file
readfits $file

listobjs 
newh2d C1ri -128.5 128.5 257 -128.5 128.5 257
newh2d C2ri -128.5 128.5 257 -128.5 128.5 257
newh1d C1exp -128.5 128.5 257
newh1d C2exp -128.5 128.5 257
newh1d C1val -128.5 128.5 257
newh1d C2val -128.5 128.5 257

newvec specC1 1024 
newvec specC2 1024 

rename $OSIG sigdum

c++exec \
  int C1OFFSET=24;  \
  int C2OFFSET=2072; \
  cout << " ===== Signal Analysis ===== START " << C1OFFSET << " , " << C2OFFSET << endl; \
  for(int j=0; j<sigdum.NRows(); j++)  \
    for(int i=0; i<sigdum.NCols(); i++)  \
      sigdum(j,i) = (sigdum(j,i)<127.5)?sigdum(j,i):sigdum(j,i)-256.; \
  for(int j=0; j<sigdum.NRows(); j++) { \
    C1exp.Add(sigdum(j,C1OFFSET+2046)); \
    C2exp.Add(sigdum(j,C2OFFSET+2046)); \
    for(int i=C1OFFSET; i<C1OFFSET+2042; i+=2) { \
        specC1((i-C1OFFSET)/2) += sqrt(sigdum(j,i)*sigdum(j,i)+sigdum(j,i+1)*sigdum(j,i+1));  \
        C1ri.Add(sigdum(j,i),sigdum(j,i+1));  \
        C1val.Add(sigdum(j,i));  C1val.Add(sigdum(j,i+1));  \ 
    }  \
    for(int i=C2OFFSET; i<C2OFFSET+2042; i+=2) { \
        specC2((i-C2OFFSET)/2) += sqrt(sigdum(j,i)*sigdum(j,i)+sigdum(j,i+1)*sigdum(j,i+1));  \
        C2ri.Add(sigdum(j,i),sigdum(j,i+1));  \
        C2val.Add(sigdum(j,i));  C2val.Add(sigdum(j,i+1));  \ 
    }  \
  }  \
  specC1 /= (double)sigdum.NRows();  \
  specC2 /= (double)sigdum.NRows();  
#
# 	
echo ' --- processing done, start plotting ... '

mv sigdum $OSIG

setaxesatt 'font=helvetica,bold,16 fixedfontsize minorticks'
zone 3 2 
disp C1exp 'font=helvetica,roman,14 blue' 
disp C1val 'font=helvetica,roman,14 blue' 
#disp C1ri 'font=helvetica,roman,14 blue'  
plot2d C1ri x y val>0 "notit nsta blue"
disp C2exp 'font=helvetica,roman,14 red' 
disp C2val 'font=helvetica,roman,14 red' 
#disp C2ri 'font=helvetica,roman,14 red' 
plot2d C2ri x y val>0 "notit nsta red"

newwin 
zone 1 2 
plot2d specC1 n*250/4096+125 val 1 'blue cpts font=helvetica,roman,14 nsta notit'
settitle 'Spec Voie 1 ' ' ' ' blue font=helvetica,bold,16 '
plot2d specC2 n*250/4096+125 val 1 'red cpts font=helvetica,roman,14 nsta notit'
settitle 'Spec Voie 2 ' ' ' ' red font=helvetica,bold,16 '