Commit 10062417 authored by OP's avatar OP
Browse files

Merge branch 'master' of gitlab.in2p3.fr:baoradio/AnaPAON4

parents 422ddf17 ca975190
Some GuideLines to tune postFringe
---------------------------------
JEC: 7 Juin 18
o Pour postFringe v3
Running fitphi.pic selon
(cd CasA; spiapp -term -exec fitphi.pic cor_CasA30mai18b ; stty sane)
Le plot des residus des fits permets d'investiguer des problemes
dans la plus part des cas venant de la procedure qui met bout a bout
les segments de phase defilant entre -Pi et Pi.
La routine incriminee est TransformFitData qui drecoit en entree
une collection (GeneralFitData) de points (freq_i, val_i) qui sont
des medianes sur un certain nombres de points d'origine
(voir estimePhase les vecteurs x et y issus de medianReduct)
Donc dans TransformFitData:
1) on cherche les points de transition -Pi->Pi (cas actuel)
ou +Pi->-Pi par la methode des differences adjacentes, puis en
la seuillant par une valeur Min et Max (TUNING POSSIBLE) pour
reconnaitre des points pathologiques que l'on suprime
2) avec le nouveau set de points cleaned on reprocede de nouveau
a des differences adjacentes. Un nouveau seuillage (TUNING POSSIBLE)
permet de determiner les points de transitions (vecteur idy)
3) on controle si il y a des pts de transition aberrants par exemple
en controlant que la phase avant et apres ce point doit faire un saut
suffisant (TUNING POSSIBLE)
4) une fois les points de transition valides on procede
progressivement au calcul des nouvelles phases en "additionnant" -2pi
selon le signe de la pente de variation.
Une fois les nouveaux points (freq_i, phase_i) on procede au fit
polynomial de la relation phase(freq). Cependant, on procede
progressivement en prenant un polynome de degre {1,...., polDegreMax}
avec polDegreMax controle en ligne de commande par (-degre <d>).
Le controle se fait en calculant l'inter-quartile normalise des
valeurs absolues des residus. Le seuil qui accepte un polynome
est iqrn_thr (TUNING POSSIBLE).
Si aucun polynome est satisfaisant alors l'objet PolyPhi<phase> n'est
cree dans le ficher de controle
fOutName = outputDir + "/" + "PhaseParam-" + source + "-" + Num2String(FreqMin) + "-" + Num2String(FreqMax) + "-" + Type +".txt";
#! /bin/bash
tag="$1"
find *${tag}*.eps | awk 'BEGIN { FS="." } {print $1 }' | xargs -I {} convert -background white -flatten {}.eps {}.png
set rootdir ./CasA/
set source $1
setaxesatt "font=times,bold,20 fixedfontsize"
graphicatt "font=times,bold,20 fixedfontsize grid"
set vecAUTO ( 1H 2H 3H 4H )
newwin 2 2
newwin 2 2 900 900
foreach x ( $vecAUTO )
openppf fitAllAuto-${source}-${x}.ppf
openppf ${rootdir}fitAllAuto-${source}-${x}.ppf
mv dataFit dataFit${x}
mv dataFit_func dataFit_func${x}
nt2d dataFit${x} x0 y ex0 ey 1 ! "marker=fcircle,5 nsta ntit"
nt2d dataFit_func${x} x0 y ex0 ey 1 ! "same cpts red nsta ntit"
nt2d dataFit_func${x} x0 y ex0 ey 1 ! "same cpts red thickdashedline nsta ntit"
settitle "${source} : ${x}"
setaxelabels "h.a (rad)" "Auto (a.u)" "font=times,bold,20 fixedfontsize"
end
w2eps ${source}-autofit.eps
set vecCROSS ( 1H2H 1H3H 1H4H 2H3H 2H4H 3H4H )
newwin 3 2
newwin 3 2 1200 900
foreach x ( $vecCROSS )
openppf fitAllCross-${source}-${x}-RE.ppf
openppf ${rootdir}fitAllCross-${source}-${x}-RE.ppf
mv dataFit dataFit${x}
mv dataFit_func dataFit_func${x}
nt2d dataFit${x} x0 y ex0 ey 1 ! "marker=fcircle,5 nsta ntit"
nt2d dataFit_func${x} x0 y ex0 ey 1 ! "same cpts red nsta ntit"
nt2d dataFit_func${x} x0 y ex0 ey 1 ! "same cpts red thickdashedline nsta ntit"
settitle "${source} : ${x}"
setaxelabels "h.a (rad)" "Cross (a.u)" "font=times,bold,20 fixedfontsize"
end
w2eps ${source}-crossfit.eps
......@@ -12,6 +12,7 @@
#include "tarrinit.h"
#include "ppersist.h"
#include "matharr.h"
#include "datime.h"
//astrop
......@@ -48,7 +49,7 @@ typedef ParFitCont_t::const_iterator ParFit_CstIter;
//201: pour extraire le fit a une frequence, 301: pour en plus debbuger
#define DEBUG_THRESH 200
#define DEBUG_THRESH 0
#if DEBUG_THRESH > 200
bool firstDBG = true;
int nloop = 0;
......@@ -237,7 +238,8 @@ Beam::Beam(r_8 waveL, r_8 raSrc, r_8 decSrc, const Now& np, int nPar):
np_(np) {
//position nominale des parabolles ideale
r_8 l0 = -sin(decSrc - np.n_lat);
//JEC 13/6/18 BUG r_8 l0 = -sin(decSrc - np.n_lat);
r_8 l0 = sin(decSrc - np.n_lat);
r_8 m0 = 0.;
r_8 n0 = cos(decSrc - np.n_lat);
......@@ -260,7 +262,8 @@ r_8 Beam::Value(r_8 const xp[], r_8 const* Par) {
r_8 ha = xp[0]; // l'angle horaire de la source (LST - RA_source)
//cos directeur de la direction de la source dans le repere horizontal
r_8 ls = cos(decSrc_)*cos(ha)*sin(np_.n_lat)-sin(decSrc_)*cos(np_.n_lat);
//JEC 13/6/18 BUG r_8 ls = cos(decSrc_)*cos(ha)*sin(np_.n_lat)-sin(decSrc_)*cos(np_.n_lat);
r_8 ls = -cos(decSrc_)*cos(ha)*sin(np_.n_lat)+sin(decSrc_)*cos(np_.n_lat);
r_8 ms = cos(decSrc_)*sin(ha);
r_8 ns = sin(np_.n_lat)*sin(decSrc_)+ cos(np_.n_lat)*cos(decSrc_)*cos(ha);
Vector3d nStar(ls,ms,ns);
......@@ -477,7 +480,8 @@ r_8 Cross::Value(r_8 const xp[], r_8 const* Par) {
r_8 ha = xp[0]; // l'angle horaire de la source (LST - RA_source)
//cos directeur de la direction de la source dans le repere horizontal
r_8 ls = cos(decSrc_)*cos(ha)*sin(np_.n_lat)-sin(decSrc_)*cos(np_.n_lat);
//JEC 13/6/18 BUG r_8 ls = cos(decSrc_)*cos(ha)*sin(np_.n_lat)-sin(decSrc_)*cos(np_.n_lat);
r_8 ls = -cos(decSrc_)*cos(ha)*sin(np_.n_lat)+sin(decSrc_)*cos(np_.n_lat);
r_8 ms = cos(decSrc_)*sin(ha);
r_8 ns = sin(np_.n_lat)*sin(decSrc_)+ cos(np_.n_lat)*cos(decSrc_)*cos(ha);
Vector3d nStar(ls,ms,ns);
......@@ -923,6 +927,7 @@ void fitAutoCrossCorr(){
r_8 decSrc = degrad(decSrcDeg+decSrcMin/60.+decSrcSec/3600.); //J2000
cout << "raSrc (hr): " << raSrc << ", decSrc (rad): " << decSrc << "\n";
......@@ -936,6 +941,9 @@ void fitAutoCrossCorr(){
if(fabs(hourAngleMin)<HATHRESHOLD)
throw ParmError("fitAllCross |hourAngleMin|<HATHRESHOLD");
Now LocalGeoTime = param.LocalGeoTime;
//--------------------
// Prepare numbering & naming stuff
//--------------------
......@@ -965,10 +973,21 @@ void fitAutoCrossCorr(){
vector<r_8> posXCh(Ndish), posYCh(Ndish), posZCh(Ndish);
//Geometrical initial values
posYCh[0] = 0.; posXCh[0] = 0.; posZCh[0] = 0.;
posYCh[1] = 5.993; posXCh[1] = 0.; posZCh[1] = 0.;
posYCh[2] = -4.380; posXCh[2] = -5.996; posZCh[2] = 0.;
posYCh[3] = -4.380; posXCh[3] = 5.995; posZCh[3] = 0.;
posXCh[0] = 0.;
posYCh[0] = 0.;
posZCh[0] = 0.;
posXCh[1] = 0.;
posYCh[1] = 5.993;
posZCh[1] = 0.;
posXCh[2] = -5.996;
posYCh[2] = -4.380;
posZCh[2] = 0.;
posXCh[3] = 5.995;
posYCh[3] = -4.380;
posZCh[3] = 0.;
vector< pair<int, int> > crossPair; //(i,j) numbering starts at 0
......@@ -1105,31 +1124,41 @@ void fitAutoCrossCorr(){
PInPersist fIn(fName);
//Time vector transformation in Hour Angle (MAY BE ALREADY DONE IN THE INPUT FILE)
TArray<r_8> TimeVec;
fIn >> PPFNameTag("TimeVec") >> TimeVec;
TimeVec.Show(cout);
TVector<r_4> HAVec(TimeVec.Size());
//JEC 10/6/18 No more usefull, use RAVec directly
// //Time vector transformation in Hour Angle (MAY BE ALREADY DONE IN THE INPUT FILE)
// TArray<r_8> TimeVec;
// fIn >> PPFNameTag("TimeVec") >> TimeVec;
// TimeVec.Show(cout);
int year0 = param.year0;
int month0 = param.month0;
int iday0 = param.iday0;
Now LocalGeoTime = param.LocalGeoTime;
double lst;
double ha;
for(sa_size_t i=0; i<TimeVec.Size();i++){
r_8 day0 = iday0 +HdecfrHMS(0,0,0)/24.; //noon
day0 += TimeVec(i)/(24.*3600.);
LocalGeoTime.n_mjd = MJDfrDate(day0,month0,year0);
now_lst(&LocalGeoTime,&lst);
ha = hrrad(lst-raSrc);
if(ha<-M_PI) // to use [-PI; PI] convention
ha += 2*M_PI;
HAVec(i) = ha;
// cout << "t , Ra(rad) " << TimeVec(i) << ", " << HAVec(i) << "\n";
}
//RAVec
TArray<r_8> RAVec;
fIn >> PPFNameTag("RAVec") >> RAVec;
TVector<r_4> HAVec(RAVec.Size());
//JEC 9/6/18 seems that xephem routine has bugs
for(sa_size_t i=0; i<RAVec.Size();i++){
HAVec(i) = hrrad(RAVec(i)-raSrc); //hour angle = RA-raSrc hour -> radian
}
// int year0 = param.year0;
// int month0 = param.month0;
// int iday0 = param.iday0;
// double lst;
// double ha;
// for(sa_size_t i=0; i<TimeVec.Size();i++){
// r_8 day0 = iday0 +HdecfrHMS(0,0,0)/24.; //noon
// day0 += TimeVec(i)/(24.*3600.);
// LocalGeoTime.n_mjd = MJDfrDate(day0,month0,year0);
// now_lst(&LocalGeoTime,&lst);
// ha = hrrad(lst-raSrc);
// if(ha<-M_PI) // to use [-PI; PI] convention
// ha += 2*M_PI;
// HAVec(i) = ha;
// // cout << "t , Ra(rad) " << TimeVec(i) << ", " << HAVec(i) << "\n";
// }
//Selection des Angles Horaires
......@@ -2035,32 +2064,36 @@ void rotateVisi() {
TArray<r_8> TimeVec;
fInPhFreeVisi >> PPFNameTag("TimeVec") >> TimeVec;
TVector<r_4> HAVec(TimeVec.Size());
//save Freq in output
fOut << PPFNameTag("TimeVec") << TimeVec;
{//just to be complete transfert RAVec to output
//just to be complete transfert RAVec to output
TArray<r_8> RAVec;
fInPhFreeVisi >> PPFNameTag("RAVec") >> RAVec;
fOut << PPFNameTag("RAVec") << RAVec;
}
double lst;
double ha;
for(sa_size_t i=0; i<TimeVec.Size();i++){
r_8 day0 = iday0 +HdecfrHMS(0,0,0)/24.; //noon
day0 += TimeVec(i)/(24.*3600.);
LocalGeoTime.n_mjd = MJDfrDate(day0,month0,year0);
now_lst(&LocalGeoTime,&lst);
ha = hrrad(lst-raSrc);
if(ha<-M_PI) // to use [-PI; PI] convention
ha += 2*M_PI;
HAVec(i) = ha;
// cout << "t , Ra(rad) " << TimeVec(i) << ", " << HAVec(i) << "\n";
}
//JEC 9/6/18 seems that xephem routine has bugs
TVector<r_4> HAVec(RAVec.Size());
for(sa_size_t i=0; i<RAVec.Size();i++){
HAVec(i) = hrrad(RAVec(i)-raSrc); //hour angle = RA-raSrc hour -> radian
}
//BUG MJDfrDate
// double lst;
// double ha;
// for(sa_size_t i=0; i<TimeVec.Size();i++){
// r_8 day0 = iday0 +HdecfrHMS(0,0,0)/24.; //noon
// day0 += TimeVec(i)/(24.*3600.);
// LocalGeoTime.n_mjd = MJDfrDate(day0,month0,year0);
// now_lst(&LocalGeoTime,&lst);
// ha = hrrad(lst-raSrc);
// if(ha<-M_PI) // to use [-PI; PI] convention
// ha += 2*M_PI;
// HAVec(i) = ha;
// // cout << "t , Ra(rad) " << TimeVec(i) << ", " << HAVec(i) << "\n";
// }
//Selection des Angles Horaires
......@@ -2193,7 +2226,8 @@ void rotateVisi() {
if(iha != -1) {
//cos directeur de la direction de la source dans le repere horizontal
r_8 lat = LocalGeoTime.n_lat;
r_8 ls = cos(decSrc)*cos(ha)*sin(lat)-sin(decSrc)*cos(lat);
//BUG 13/6/18 r_8 ls = cos(decSrc)*cos(ha)*sin(lat)-sin(decSrc)*cos(lat);
r_8 ls = - cos(decSrc)*cos(ha)*sin(lat)+sin(decSrc)*cos(lat);
r_8 ms = cos(decSrc)*sin(ha);
r_8 ns = sin(lat)*sin(decSrc) + cos(lat)*cos(decSrc)*cos(ha);
Vector3d nStar(ls,ms,ns);
......@@ -2227,6 +2261,324 @@ void rotateVisi() {
}//rotateVisi
void BaselinesPhase() {
//Param extractor
int debug = param.debug;
string inputDataDir = param.inputDir;
string inputParamDir = param.inputParamDir;
string inputFile = param.inputFile;
string inputParamFile = param.inputParamFile;
string outputDir = param.outputDir;
string outputFile = param.outputFile;
r_4 FreqMin = param.FreqMin;
r_4 FreqMax = param.FreqMax;
string Type = param.Type; //RE vs IM
string source = param.source;
r_8 raSrcHour = param.raSrcHour;
r_8 raSrcMin = param.raSrcMin;
r_8 raSrcSec = param.raSrcSec;
r_8 raSrc = HdecfrHMS(raSrcHour,raSrcMin,raSrcSec); //J2000
r_8 decSrcDeg = param.decSrcDeg;
r_8 decSrcMin = param.decSrcMin;
r_8 decSrcSec = param.decSrcSec;
r_8 decSrc = degrad(decSrcDeg+decSrcMin/60.+decSrcSec/3600.); //J2000
r_4 hourAngleMax = param.hourAngleMax;
r_4 hourAngleMin = param.hourAngleMin;
int year0 = param.year0;
int month0 = param.month0;
int iday0 = param.iday0;
Now LocalGeoTime = param.LocalGeoTime;
cout << "raSrc (hr): " << raSrc << ", decSrc (rad): " << decSrc << "\n";
TimeStamp meridianT = param.MeridianDate;
size_t Ndish = 4;
size_t Nauto = Ndish; //JEC 3/6/16 NPOLAR dependant
size_t Ncross = Ndish*(Ndish-1)/2; //2*(Ndish*(Ndish-1)/2) and Ndish = Nauto/2
cout << "Npolar, Ndish, Nauto, Ncross: " << TYPE_POLAR << " " << Ndish << " " << Nauto << " " << Ncross << endl;
vector<string> nameCh(Nauto);
#if TYPE_POLAR == 1
nameCh[0] = "1H";
nameCh[1] = "2H";
nameCh[2] = "3H";
nameCh[3] = "4H";
#else
nameCh[0] = "1V";
nameCh[1] = "2V";
nameCh[2] = "3V";
nameCh[3] = "4V";
#endif
//Position par defaut des dipoles.
// ATTENTION !!!!!!!!! Repere Oriente vers North positif X, vers West positif Y, Z vertical
vector<r_8> posXCh(Ndish), posYCh(Ndish), posZCh(Ndish);
//Geometrical initial values
posXCh[0] = 0.;
posYCh[0] = 0.;
posZCh[0] = 0.;
posXCh[1] = 0.;
posYCh[1] = 5.993;
posZCh[1] = 0.;
posXCh[2] = -5.996;
posYCh[2] = -4.380;
posZCh[2] = 0.;
posXCh[3] = 5.995;
posYCh[3] = -4.380;
posZCh[3] = 0.;
vector< pair<int, int> > crossPair; //(i,j) numbering starts at 0
vector<string> crossName;
for(size_t i=0;i<Nauto-1;i++){
for(size_t j=i+1;j<Nauto;j++){
crossPair.push_back(make_pair(i,j));
crossName.push_back(nameCh[i]+nameCh[j]);
}
}
//---------
//File for the output
//---------
string fOutName="";
if(outputFile != ""){
fOutName = outputDir + "/" + outputFile;
} else {
fOutName = outputDir + "/" + "TFM-BaseLine-" + source +".ppf";
}
cout << "Geometrical Phase output in file <"<<fOutName<<">" << endl;
POutPersist fOut(fOutName);
//--------------------
//---------
//Data loading
//---------
//--------------------
//---------
// Original visibilities
//---------
string fName = "";
if(inputFile != ""){
fName = inputDataDir + "/" + inputFile;
} else {
fName = inputDataDir + "/" + source + ".ppf";
}
cout << "Load the raw TFM from <"<<fName<<">\n";
PInPersist fIn(fName);
//---------
//array of frequencies and Time (saved in fNamePhFreeVisi JEC 29/11/17)
//---------
TVector<r_8> FreqVecOrig;
fIn >> PPFNameTag("FreqVec") >> FreqVecOrig;
cout << "FreqVecOrig: nfreq= " << FreqVecOrig.Size() << endl;
//save Freq in output
fOut << PPFNameTag("FreqVec") << FreqVecOrig;
TArray<r_8> TimeVec;
fIn >> PPFNameTag("TimeVec") >> TimeVec;
//save Freq in output
fOut << PPFNameTag("TimeVec") << TimeVec;
//save RAVec in output
TArray<r_8> RAVec;
fIn >> PPFNameTag("RAVec") >> RAVec;
fOut << PPFNameTag("RAVec") << RAVec;
cout << "TimeVec, RAVec n: " << TimeVec.Size() << ", " << RAVec.Size() << endl;
TVector<r_4> HAVec(RAVec.Size());
// double lst;
// double ha;
// for(sa_size_t i=0; i<TimeVec.Size();i++){
// r_8 day0 = iday0 +HdecfrHMS(0,0,0)/24.; //noon
// day0 += TimeVec(i)/(24.*3600.);
// LocalGeoTime.n_mjd = MJDfrDate(day0,month0,year0);
// now_lst(&LocalGeoTime,&lst);
// ha = hrrad(lst-raSrc);
// if(ha<-M_PI) // to use [-PI; PI] convention
// ha += 2*M_PI;
// HAVec(i) = ha;
// // cout << "t , Ra(rad) " << TimeVec(i) << ", " << HAVec(i) << "\n";
// }
//JEC 9/6/18 seems that xephem routine has bugs
for(sa_size_t i=0; i<RAVec.Size();i++){
HAVec(i) = hrrad(RAVec(i)-raSrc); //hour angle = RA-raSrc hour -> radian
}
//---------
//loop on Visibility types and compute Geometrical Baseline Phases
//---------
for(size_t ix=0; ix<Ncross; ix++){ //Loop on cross-corr
// DataTable dt; //pour debugger
// dt.SetShowMinMaxFlag(true);
// dt.AddDoubleColumn("freq"); //0
// dt.AddDoubleColumn("ha"); //1
// dt.AddDoubleColumn("ra"); //2
// dt.AddDoubleColumn("time"); //3
// dt.AddDoubleColumn("blph"); //4
// dt.AddDoubleColumn("rawph"); //5
// DataTableRowPtr rowp=dt.EmptyRowPtr();
int iCh = crossPair[ix].first; //0-indexed
int jCh = crossPair[ix].second;
cout << "Cross["<<ix<<"]: pair ("<< iCh << ", " << jCh
<< "), name " << crossName[ix]
<<endl;
//The values PhaseFree Cross visibility
string objName = "TFM_"+crossName[ix];
TMatrix< complex<r_4> > inVisi; // CARE here we load Matrix to use Range selection afterwards
fIn >> PPFNameTag(objName) >> inVisi;
sa_size_t nrows = inVisi.NRows(); // frequency axis
sa_size_t ncols = inVisi.NCols(); // RA axis
cout << "nrows(nFreq), ncols(nRA)" << nrows << " , " << ncols << endl;
TMatrix<r_8> geomPhase(nrows,ncols);
for (sa_size_t irow=0; irow<nrows; irow++){ //frequence axis
r_8 freq = FreqVecOrig(irow);
if (irow<10)cout << "freq: " << freq << endl;
r_8 waveL = FreqToWave(freq);
//-----------
// baseline parameters
//-----------
r_8 xi = posXCh[iCh]; //O-indexed
r_8 yi = posYCh[iCh];
r_8 zi = posZCh[iCh];
cout << "Chi: " << iCh
<< "x : " << xi << " "
<< "y : " << yi << " "
<< "z : " << zi << " "
<< endl;
r_8 xj = posXCh[jCh]; //O-indexed
r_8 yj = posYCh[jCh];
r_8 zj = posZCh[jCh];
cout << "Chj: " << jCh
<< "x : " << xj << " "
<< "y : " << yj << " "
<< "z : " << zj << " "
<< endl;
//base line parameters
r_8 LSN = xj - xi;
r_8 LEW = yj - yi;
r_8 LVERT = zj - zi;
Vector3d Dx(LSN, LEW, LVERT);
//RA/time/HA servey
for (sa_size_t icol=0; icol<ncols; icol++){
r_8 ha = HAVec(icol);
if (icol<10)cout << "HA (= RA): " << ha << endl;
//cos directeur de la direction de la source dans le repere horizontal
r_8 lat = LocalGeoTime.n_lat;
//BUG 13/6/18 r_8 ls = cos(decSrc)*cos(ha)*sin(lat)-sin(decSrc)*cos(lat);
r_8 ls = - cos(decSrc)*cos(ha)*sin(lat)+sin(decSrc)*cos(lat);
r_8 ms = cos(decSrc)*sin(ha);
r_8 ns = sin(lat)*sin(decSrc) + cos(lat)*cos(decSrc)*cos(ha);
Vector3d nStar(ls,ms,ns);
//Geometrical phase
r_8 phase = -2*M_PI/waveL * nStar.Psc(Dx);
phase = std::remainder(phase,2*M_PI); //modulo 2Pi
// dt.NextRowPtr(rowp);
// rowp(0) = freq;
// rowp(1) = ha; //HA
// rowp(2) = RAVec(icol); //RA
// rowp(3) = TimeVec(icol);
// rowp(4) = phase;
// rowp(5) = (r_8) (std::arg(inVisi(irow,icol)));
geomPhase(irow,icol) = phase;
}//loop on time/RA/HA
}//loop on freq axis