Skip to content

GitLab

Commit c7924977 authored by Reza ANSARI's avatar Reza ANSARI
Browse files

ajout programme visi2dtacx.cc (remplissage datatable auto-corr + cross-corr... 

ajout programme visi2dtacx.cc (remplissage datatable auto-corr + cross-corr 1H-2H 1H-3H ... 3H-4H) et script de trace placx.pic , et MAJ script shell doproc, Reza 19/05/2015
parent 2f505fb7
Hide whitespace changes
Inline Side-by-side
Showing with 369 additions and 0 deletions
placx.pic 0 → 100644
View file @ c7924977
################## Script placx.pic #################
## Script pour tracer les 8 auto-correlations et 6 cross-correlation 1H-2H ... 3H-4H de PAON4
if ( $# < 1 ) then
echo 'placx.pic/missing argument '
echo 'Usage: placx input_base_name (without .ppf)'
return
endif
set flnm $1
set titre 'PAON4 - transit'
# Pour les axes
setaxesatt 'font=helvetica,bold,16 fixedfontsize minorticks'
# Ouvrir le fichier
openppf $flnm.ppf
set DT $flnm
print $DT
## Facteur de normalisation entre les voies x 1.e-4
## 1 1/3.48 1.06 0.92 1 0.98 0.377 1.32
set FNORM ( 1.e-4 0.287e-4 1.06e-4 0.92e-4 1.e-4 0.96e-4 0.377e-4 1.32e-4 )
## couleurs pour les traces
set COLS ( black red orange brown blue green skylbue forestgreen )
###### Trace de la puissance en bande 1
## On trace p1 en fct du temps
newwin 1 1 800 600
plot2d $DT time/3600 p1*$FNORM[0] 1 " nsta cpts notit $COLS[0] "
### trace de p2 ... p8
for i 1:8
j = i+1
plot2d $NT time/3600 p$j*$FNORM[i] 1 "nsta cpts same notit $COLS[i] "
end
### Axes et titres
setaxelabels 'Time (hours)' 'P1...P8 (arb.units)' 'font=helvetica,bolditalic,16 black'
settitle "$titre Bande 1" ' ' 'font=helvetica,bold,18 black'
###### Trace de la puissance en bande HI 1420+-1MHz
## On trace p1 en fct du temps
newwin 1 1 800 600
plot2d $DT time/3600 p1HI*$FNORM[0] 1 " nsta cpts notit $COLS[0] "
### trace de p2 ... p8
for i 1:8
j = i+1
plot2d $NT time/3600 p$jHI*$FNORM[i] 1 "nsta cpts same notit $COLS[i] "
end
### Axes et titres
setaxelabels 'Time (hours)' 'P1...P8 (arb.units)' 'font=helvetica,bolditalic,16 black'
settitle "$titre @1420 MHz" ' ' 'font=helvetica,bold,18 black'
###### Trace des 6 cross-correlations
set XCOR ( xh1h2 xh1h3 xh1h4 xh2h3 xh2h4 xh3h4 )
set XCORTIT ( 1H-2H 1H-3H 1H-4H 2H-3H 2H-4H 3H-4H )
## On trace p1 en fct du temps
newwin 3 2 900 600
for i 0:6
plot2d $DT time/3600 sqrt(pow($XCOR[i].real,2)+pow($XCOR[i].imag,2)) 1 " nsta cpts notit black "
plot2d $DT time/3600 $XCOR[i].real 1 " nsta cpts notit same blue "
plot2d $DT time/3600 $XCOR[i].imag 1 " nsta cpts notit same red "
setaxelabels 'Time (hours)' 'Cross-Corr (arb.units)' 'font=helvetica,bolditalic,16 black'
settitle "$titre $XCORTIT[i]" ' ' 'font=helvetica,bold,18 black'
end
###########################################################
visi2dtacx.cc 0 → 100644
View file @ c7924977
// Utilisation de SOPHYA pour faciliter les tests ...
#include "sopnamsp.h"
#include "machdefs.h"
/* ----------------------------------------------------------
Projet BAORadio/PAON4 - (C) LAL/IRFU 2008-2015
visi2dtacx: programme de lecture des fichiers matrices de
visibilites de PAON4, remplissage et ecriture d'un
DataTable avec les cross-correlations en fonction du temps
R. Ansari, J.E. Campagne, C. Magneville - LAL/Irfu
---------------------------------------------------------- */
// include standard c/c++
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <iostream>
#include <string>
#include "pexceptions.h"
#include "tvector.h"
#include "fioarr.h"
// #include "tarrinit.h"
#include "ntuple.h"
#include "datatable.h"
#include "histinit.h"
#include "matharr.h"
#include "timestamp.h"
#include <utilarr.h>
// include sophya mesure ressource CPU/memoire ...
#include "resusage.h"
#include "ctimer.h"
#include "timing.h"
// include lecteur de fichiers visibilites
#include "visip4reader.h"
//--------------------------- Fonctions de ce fichier -------------------
int Usage(bool fgshort=true);
// int DecodeArgs(int narg, char* arg[]);
/* --Fonction-- */
int Usage(bool fgea)
{
cout << " --- visi2dtacx.cc : Read PPF files produced by mfacq time-frequency\n" << endl;
if (fgea) cout << " visi2dtacx: Missing or wrong argument ! " << endl;
cout << " Usage: visi2dtacx InPathBAO5 InPathBAO6 Imin,Imax GainPPFFile OutPPFFile Freq1 Freq2 DeltaFreq DeltaIAvg [PrtLev=0] \n"
<< " InPathBAO5/6: path for input vismtx_J_iii.ppf on bao5/6 \n"
<< " Imin,Imax : read files for Imin<=iii<=Imax iii+=step \n"
<< " GainPPFFile: Input gains PPF file name \n"
<< " OutPPFFile: Output PPF file name \n"
<< " Freq1 Freq2: the two frequencies for which cross-correlations are computed [=1375 1410 MHz] \n"
<< " DeltaFreq: frequency band (in MHz) [=1 MHz] \n"
<< " DeltaIAvg: compute average power every DeltaI input vismtx files, def=10 \n"
<< " Note: bande_freq=[freq1,2-deltafreq , freq1,2+deltafreq] \n"<< endl;
/*
if (fgshort) {
cout << " rdvisip4 -h for detailed instructions" << endl;
return 1;
}
*/
return 1;
}
//----------------------------------------------------
//----------------------------------------------------
int main(int narg, char* arg[])
{
if ((narg>1)&&(strcmp(arg[1],"-h")==0)) return Usage(false);
if (narg<10) return Usage(true);
string path5 = arg[1];
string path6 = arg[2];
int Imin,Imax,Istep=1;
sscanf(arg[3],"%d,%d",&Imin,&Imax);
Istep=1;
string gainfile=arg[4];
string outfile=arg[5];
// frequency range definition
double deltafreq=1.;
double freq1=1375.;
double freq2=1410.;
if ((narg>7)&&(strcmp(arg[6],"-")!=0)) freq1=atof(arg[6]);
if ((narg>8)&&(strcmp(arg[7],"-")!=0)) freq2=atof(arg[7]);
if ((narg>6)&&(strcmp(arg[8],"-")!=0)) deltafreq=atof(arg[8]);
// time averaging interval definition, in term of visibility files
int deltaIavg=10;
if ((narg>9)&&(strcmp(arg[9],"-")!=0)) deltaIavg=atoi(arg[9]);
if (deltaIavg<1) deltaIavg=1;
int prtlev=1;
if (narg>10) prtlev=atoi(arg[10]);
cout << " visi2dtacx/Info: Path BAO5:"<<path5<<" BAO6:"<<path6<<"\n"
<< " Imin,max,step="<<Imin<<","<<Imax<<","<<Istep<<" OutFile:"<<outfile<<"\n"
<< " freq1="<<freq1<<" freq2="<<freq2<<" deltafreq="<<deltafreq<<" MHz\n"
// << JFmin<<" <=NumFreq<= "<<JFmax<<" "<<JFminB<<" <=NumFreqB<= "<<JFmaxB
<<" DeltaIAvg="<<deltaIavg<<" PrtLev="<<prtlev<<endl;
//---- calcul des index de bandes de frequences
// pas en frequence du firmware FFT
double deltanufft=250./4096; // 250 MHz en 4096 frequences
sa_size_t DJF=deltafreq/deltanufft;
double freqstart=1250.; // Bande de 1250-1500 MHz
sa_size_t JFmin=(freq1-1250.-deltafreq)/deltanufft;
sa_size_t JFmax=(freq1-1250.+deltafreq)/deltanufft;
sa_size_t JFminB=(freq2-1250.-deltafreq)/deltanufft;
sa_size_t JFmaxB=(freq2-1250.+deltafreq)/deltanufft;
cout << " visi2dtacx/Info: frequency bands ... " <<endl;
cout << " Band 1: "<<freq1<<" +/-"<<deltafreq<<" MHz -> "<<JFmin<<" <=NumFreq<= "<<JFmax<<endl;
cout << " Band 2: "<<freq2<<" +/-"<<deltafreq<<" MHz -> "<<JFminB<<" <=NumFreq<= "<<JFmaxB<<endl;
// Bande de frequence de 2 MHz autour du 21 cm (1420 MHz -> JF=2792)
sa_size_t JFmin21=2776,JFmax21=2808; // +/- 1 MHz autour de 1420 MHz
sa_size_t JFmin21_5=2710,JFmax21_5=2874; // +/- 5 MHz autour de 1420 MHz
cout << " pjHI +/-1 MHz @1420 MHz " << JFmin21<<" <=NumFreq<= "<<JFmax21
<< " pjHI +/-5 MHz @1420 MHz " << JFmin21_5<<" <=NumFreq<= "<<JFmax21_5 << endl;
HiStatsInitiator _inia;
// TArrayInitiator _inia;
int rc = 0;
try {
ResourceUsage resu;
// Numero des lignes des auto-correlations dans la matrice des visibilites
sa_size_t KVAC[8]={0,8,15,21,26,30,33,35};
// Numero des lignes des cross-correlations 1H-2H 1H-3H 1H-4H 2H-3H 2H-4H 3H-4H dans la matrice des visibilites
sa_size_t KVCXHH[6]={1,2,3,9,10,16};
cout << " visi2dtacx/Info: reading input gain file"<<gainfile<<" ...";
Matrix gains;
Vector gn;
{
PInPersist pin(gainfile);
pin>>PPFNameTag("gains")>>gains;
pin>>PPFNameTag("gn")>>gn;
for (int p=0; p<8; p++) cout<<"gn["<<p<<"]="<<gn(p)<<" ";
// on transforme gains en facteur multiplicatif
for(sa_size_t r=0;r<gains.NRows();r++) {
for(sa_size_t c=0; c<gains.NCols(); c++) gains(r,c)=1./gains(r,c);
}
cout<<endl;
}
// Gains pour les cross-correlations g_ij(nu) = sqrt(g_ii(nu) 8 g_jj(nu) )
Matrix gaincx(6, gains.NCols());
for(sa_size_t c=0; c<gains.NCols(); c++) {
gaincx(0,c)=sqrt(gains(0,c)*gains(1,c)); // 1H-2H
gaincx(1,c)=sqrt(gains(0,c)*gains(2,c)); // 1H-3H
gaincx(2,c)=sqrt(gains(0,c)*gains(3,c)); // 1H-4H
gaincx(3,c)=sqrt(gains(1,c)*gains(2,c)); // 2H-3H
gaincx(4,c)=sqrt(gains(1,c)*gains(3,c)); // 2H-4H
gaincx(5,c)=sqrt(gains(2,c)*gains(3,c)); // 3H-4H
}
cout << " DONE"<<endl;
const char* dtnames[53]={"k","day","time",
"p1","p2","p3","p4","p5","p6","p7","p8",
"p1B","p2B","p3B","p4B","p5B","p6B","p7B","p8B",
"p1HI","p2HI","p3HI","p4HI","p5HI","p6HI","p7HI","p8HI",
"p1HI5","p2HI5","p3HI5","p4HI5","p5HI5","p6HI5","p7HI5","p8HI5",
"xh1h2","xh1h3","xh1h4","xh2h3","xh2h4","xh3h4",
"xh1h2B","xh1h3B","xh1h4B","xh2h3B","xh2h4B","xh3h4B",
"xh1h2HI","xh1h3HI","xh1h4HI","xh2h3HI","xh2h4HI","xh3h4HI"};
DataTable dt(384);
dt.AddIntegerColumn(dtnames[0]);
dt.AddIntegerColumn(dtnames[1]);
dt.AddFloatColumn(dtnames[2]);
for(int kk=0; kk<32; kk++)
dt.AddFloatColumn(dtnames[3+kk]);
for(int kk=0; kk<18; kk++)
dt.AddComplexColumn(dtnames[3+32+kk]);
Range freqs(JFmin,JFmax);
Range freqs21(JFmin21,JFmax21);
Range freqs21_5(JFmin21_5,JFmax21_5);
Range freqsB(JFminB,JFmaxB);
vector<string> paths;
paths.push_back(path5);
paths.push_back(path6);
VisiP4Reader vreader(paths, Imin,Imax,Istep,true);
vreader.setPrintLevel(prtlev);
bool fgok=true;
TMatrix< complex<r_4> > vismtx;
TMatrix< complex<r_4> > acsum;
TMatrix< complex<r_4> > cxsum;
double acdt[32]; // les 4*8=32 valeurs d'autocorrelation pour remplissage dans la table
complex<double> cxdt[18]; // les 3*6=18 valeurs de cross-correlations pour remplissage dans la table
TimeStamp dateobs;
TimeStamp dateorg(2015,1,1,12,0,0.); // Date origine 1 jan 2015
double mttag;
int cnt=0, cntnt=0, pcntnt=0;
int I=0;
// Getting en empty row_ptr
DataTableRowPtr rowp = dt.EmptyRowPtr();
//
while (fgok) {
fgok=vreader.ReadNext(vismtx, dateobs, mttag);
if (fgok) {
if (cnt==0) { //resizing matrices for sum of auto-correlations and sum of 6 cross-correlations
acsum.SetSize(8, vismtx.NCols());
cxsum.SetSize(6, vismtx.NCols());
}
if (I==0) { // start filling a new DataTable row
dt.NextRowPtr(rowp);
rowp(0)=(int_4)cnt;
rowp(1)=(int_4)(dateobs.DaysPart()-dateorg.DaysPart());
rowp(2)=(r_4)dateobs.SecondsPart();
if (prtlev>0)
cout << "visi2dtacx/Info: dateobs="<<dateobs<<" SecondsPart()="<<dateobs.SecondsPart()<<endl;
for(int k=0; k<32; k++) acdt[k]=0.;
for(int k=0; k<18; k++) cxdt[k]=complex<r_8>(0.,0.);
acsum = complex<r_4>(0.,0.);
cxsum = complex<r_4>(0.,0.);
}
// sum (integration) along the time axis
for(int k=0; k<8; k++) acsum.Row(k) += vismtx.Row(KVAC[k]); // Les auto-correlations
for(int k=0; k<6; k++) cxsum.Row(k) += vismtx.Row(KVCXHH[k]); // les cross-correlations
I++; // incrementing DeltaTime counter
if (I==deltaIavg) {
//---- On s'occupe d'abord des autocorrelations P1 ... P8
Vector vac(acsum.NCols());
for(int k=0; k<8; k++) {
// correct for gain(nu)
for(sa_size_t c=0; c<vac.Size(); c++) vac(c)=acsum(k,c).real()*gains(k,c);
acdt[k] += vac(freqs).Sum(); // integration en bande 1
acdt[k+8] += vac(freqsB).Sum(); // integration en bande large B
acdt[k+16] += vac(freqs21).Sum(); // integration a +/- 1 MHz @ 1420 MHz
acdt[k+24] += vac(freqs21_5).Sum(); // integration a +/- 5 MHz @ 1420 MHz
}
//---- On s'occupe des cross-correlations 1H-2H ... 3H-4H
TVector< complex<r_8> > vcx(acsum.NCols());
for(int k=0; k<6; k++) {
for(sa_size_t c=0; c<vcx.Size(); c++)
vcx(c)=complex<r_8>((r_8)cxsum(k,c).real(), (r_8)cxsum(k,c).imag())*complex<r_8>(gaincx(k,c), 0.);
cxdt[k] += vcx(freqs).Sum(); // integration en bande 1
cxdt[k+8] += vcx(freqsB).Sum(); // integration en bande large B
cxdt[k+16] += vcx(freqs21).Sum(); // integration a +/- 1 MHz @ 1420 MHz
}
// Moyennes dans la bande 1
double fnorm=(1./(double)deltaIavg)/(double)(JFmax-JFmin+1);
for(int k=0; k<8; k++) acdt[k]*=fnorm;
for(int k=0; k<6; k++) cxdt[k]*=complex<r_8>(fnorm,0.);
// Moyennes dans la bande 2
fnorm=(1./(double)deltaIavg)/(double)(JFmaxB-JFminB+1);
for(int k=0; k<8; k++) acdt[k+8]*=fnorm;
for(int k=0; k<6; k++) cxdt[k+6]*=complex<r_8>(fnorm,0.);
// Moyennes ds +/5 MHz @1420 , exclu +/- 1 MHz centre a 1420
fnorm=(1./(double)deltaIavg)/(double)((JFmax21_5-JFmin21_5)-(JFmax21-JFmin21));
for(int k=0; k<8; k++) acdt[k+24] = (acdt[k+24]-acdt[k+16])*fnorm;
// Moyennes ds +/1 MHz @1420
fnorm=(1./(double)deltaIavg)/(double)(JFmax21-JFmin21+1);
for(int k=0; k<8; k++) acdt[k+16]*=fnorm;
for(int k=0; k<6; k++) cxdt[k+12]*=complex<r_8>(fnorm,0.);
// remplissage datatable
for(int k=0; k<32; k++) rowp(3+k) = (float)acdt[k];
for(int k=0; k<18; k++) rowp(35+k) = complex<float>((float)cxdt[k].real(), (float)cxdt[k].imag());
// ... done
I=0; cntnt++;
}
cnt++;
if ((cntnt>0)&&(cntnt%10==0)&&(cntnt>pcntnt)) {
cout << " visi2dtacx/Info: fill-count="<<cntnt<<" DateObs="<<dateobs<<endl;
pcntnt=cntnt;
}
}
}
cout << " visi2dtacx/Info: count="<<cnt<<" visimtx read "<<endl;
dt.SetShowMinMaxFlag(true);
cout << dt;
cout<<" visi2dtacx/Info: Saving auto-corr=f(time) ntuple to PPF file "<<outfile<<endl;
POutPersist po(outfile);
po<<dt;
// resu.Update();
cout << resu; // Update est fait lors du print
}
catch (PException& exc) {
cerr << " visi2dtacx.cc catched PException " << exc.Msg() << endl;
rc = 77;
}
catch (std::exception& sex) {
cerr << "\n visi2dtacx.cc std::exception :"
<< (string)typeid(sex).name() << "\n msg= "
<< sex.what() << endl;
rc = 78;
}
catch (...) {
cerr << " visi2dtacx.cc catched unknown (...) exception " << endl;
rc = 79;
}
cout << ">>>> visi2dtacx.cc ------- END ----------- RC=" << rc << endl;
return rc;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment