Commit 386171d2 authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

improvements

parent ca975190
......@@ -20,6 +20,16 @@ setaxelabels "h.a (rad)" "Auto (a.u)" "font=times,bold,20 fixedfontsize"
end
w2eps ${source}-autofit.eps
newwin 2 2 900 900
foreach x ( $vecAUTO )
n/pl dataFit${x}.y%180*x0/M_PI 1 ! "marker=fcircle,5 nsta ntit"
n/pl dataFit_func${x}.y%180*x0/M_PI 1 ! "same cpts red thickdashedline nsta ntit"
settitle "${source} : ${x}"
setaxelabels "h.a (deg)" "Auto (a.u)" "font=times,bold,20 fixedfontsize"
end
set vecCROSS ( 1H2H 1H3H 1H4H 2H3H 2H4H 3H4H )
newwin 3 2 1200 900
......
......@@ -73,6 +73,107 @@ int nloop = 0;
//-------- Utility
//------------------------
//---- MEDIAN : warning the nth_element reorder the input vector => input NO MORE VALID
//------------------------
//Warning; original order is not preserved
class median_of_empty_list_exception:public std::exception{
virtual const char* what() const throw() {
return "Attempt to take the median of an empty list of numbers. "
"The median of an empty list is undefined.";
}
};
template<class RandAccessIter>
double median(RandAccessIter begin, RandAccessIter end) {
if(begin == end){ throw median_of_empty_list_exception(); }
std::size_t size = end - begin;
std::size_t middleIdx = size/2;
RandAccessIter target = begin + middleIdx;
std::nth_element(begin, target, end);
if(size % 2 != 0){ //Odd number of elements
return *target;
}else{ //Even number of elements
double a = *target;
RandAccessIter targetNeighbor= target-1;
std::nth_element(begin, targetNeighbor, end);
return (a+*targetNeighbor)/2.0;
}
}
template<class T>
T median(const vector<T>& vec){
vector<T>tmp(vec);
return median(tmp.begin(),tmp.end());
}
template<class T>
T median(TVector<T>& vec){
vector<T>tmp(vec.ConvertTostdvec());
return median(tmp.begin(),tmp.end());
}
//interquartile range normalized to Gauss(0,1) => sigma estimator
//The sigma on iQRNorm is 1.16639*iQRNorm/sqrt(n) in case of gaussian random variable.
//warning median does not preserve original order
template<class RandAccessIter>
double iQRNorm(RandAccessIter begin, RandAccessIter end) {
if(begin == end){ throw median_of_empty_list_exception(); }
vector<r_4> tmp1(begin,end);
sort(tmp1.begin(),tmp1.end());
vector<r_4> tmp3(tmp1);
std::size_t size = tmp1.end() - tmp1.begin();
std::size_t middleIdx = size/2;
double q1 = median(tmp1.begin(), tmp1.begin()+middleIdx);
double q3 = median(tmp3.begin()+middleIdx,tmp3.end());
return (q3-q1)/1.34898; //normalisation by Q3-Q1 for a Gauss(0,1) value
}
template <class T>
T iQRNorm(TVector<T>& vec){
vector<T>tmp(vec.ConvertTostdvec());
return iQRNorm(tmp.begin(),tmp.end());
}
template <class T>
T iQRNorm(const vector<T>& vec){
vector<T>tmp(vec);
return iQRNorm(tmp.begin(),tmp.end());
}
//input : vecin
//input : nbre de bins sur lesquels est calcule la mediane des vecin[i]
//output: vecout
template<class T>
void medianReduct(int NpackBin, TVector<T>& vecin, vector<T>& vecout){
for (int ibin=0;ibin<vecout.size();ibin++){
sa_size_t ilow = ibin*NpackBin;
sa_size_t ihigh = ilow + NpackBin-1;
TVector<r_8> tmp(vecin(Range(ilow,ihigh)),false); //no data sharing for the time beeing
vecout[ibin] = median(tmp);
}
}
template <class T>
T my_remainder(T x) {
return std::remainder(x,(T)2*M_PI);
}
template <class T>
void Modulo(TMatrix<T> mtx){
MathArray<T>::ApplyFunctionInPlace(mtx,my_remainder);
}
r_8 FreqToWave(r_8 freq){
//21.1061140542 cm <-> 1420405751.7667 Hz
......@@ -2332,9 +2433,9 @@ void BaselinesPhase() {
// 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.;
//Geometrical initial values (UNITS is METER)
posXCh[0] = 0;
posYCh[0] = 0;
posZCh[0] = 0.;
posXCh[1] = 0.;
......@@ -2371,15 +2472,9 @@ void BaselinesPhase() {
cout << "Geometrical Phase output in file <"<<fOutName<<">" << endl;
{//SAve into fOut
POutPersist fOut(fOutName);
//--------------------
//---------
//Data loading
//---------
//--------------------
//---------
// Original visibilities
//---------
......@@ -2423,25 +2518,56 @@ void BaselinesPhase() {
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";
// }
// 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
}
//Selection des Angles Horaires
r_8 haMin = Angle(hourAngleMin, Angle::Degree).ToRadian();
r_8 haMAX = Angle(hourAngleMax, Angle::Degree).ToRadian();
sa_size_t ifirstHA=0;
sa_size_t ilastHA=0;
//Select the first transit of the source during observation
bool first=true;
bool last=true;
for(sa_size_t i=0; i<HAVec.Size();i++){
if(HAVec(i)<haMAX && HAVec(i)>haMin){
if(first){
first = false;
ifirstHA = i;
last = true;
}
if(!first && last)ilastHA=i;
} else {
last = false;
}
}
cout << "Select HA index in ["<<ifirstHA <<","<<ilastHA <<"]\n";
cout << " in ["<<HAVec(ifirstHA) <<","<<HAVec(ilastHA) <<"] rad\n";
if(ifirstHA==ilastHA)
throw RangeCheckError("HA selection not valid");
TVector<r_4> HASelected(HAVec(Range(ifirstHA,ilastHA)));
sa_size_t NHA = HASelected.Size();
//---------
//loop on Visibility types and compute Geometrical Baseline Phases
......@@ -2451,18 +2577,6 @@ void BaselinesPhase() {
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
......@@ -2476,15 +2590,26 @@ void BaselinesPhase() {
sa_size_t nrows = inVisi.NRows(); // frequency axis
sa_size_t ncols = inVisi.NCols(); // RA axis
// sa_size_t ncols = inVisi.NCols(); // RA axis
TMatrix< complex<r_4> > zoomInVisi(inVisi(Range::all(),Range(ifirstHA,ilastHA)));
sa_size_t ncols = NHA; // RA axis
cout << "nrows(nFreq), ncols(nRA)" << nrows << " , " << ncols << endl;
//2D map: freq, RA
TMatrix<r_8> geomPhase(nrows,ncols);
TMatrix<r_8> phaseInstrum(nrows,ncols);
//1D vector: freq
TVector<r_8> vphaseInst(nrows);
for (sa_size_t irow=0; irow<nrows; irow++){ //frequence axis
r_8 freq = FreqVecOrig(irow);
if(freq<FreqMin || freq>FreqMax) continue;
if (irow<10)cout << "freq: " << freq << endl;
r_8 waveL = FreqToWave(freq);
......@@ -2526,7 +2651,8 @@ void BaselinesPhase() {
//RA/time/HA servey
for (sa_size_t icol=0; icol<ncols; icol++){
r_8 ha = HAVec(icol);
// r_8 ha = HAVec(icol);
r_8 ha = HASelected(icol);
if (icol<10)cout << "HA (= RA): " << ha << endl;
//cos directeur de la direction de la source dans le repere horizontal
......@@ -2542,41 +2668,98 @@ void BaselinesPhase() {
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;
phaseInstrum(irow,icol) = std::remainder(std::arg(zoomInVisi(irow,icol))-phase,2*M_PI);
}//loop on time/RA/HA
//take the median on RA-axis of the instrumental phase
{
TVector<r_8>tmp(phaseInstrum(Range(irow),Range::all()));
vphaseInst(irow) = median(tmp);
}
}//loop on freq axis
//---------
//save the zoom original visibilities
//---------
{
string objNameOut = "TFM_RAWZOOM_"+crossName[ix];
fOut << PPFNameTag(objNameOut) << zoomInVisi;
}
//---------
//save the phase-free visibilities
//save the phase-geom visibilities
//---------
{
string objNameOut = "TFM_BASELINE_"+crossName[ix];
fOut << PPFNameTag(objNameOut) << geomPhase;
}
// //save dataTable
// {
// string objNameOut = "DT_"+crossName[ix];
// fOut << PPFNameTag(objNameOut) << dt;
// }
//---------
//save the instrumentaal phase visibilities
//---------
{
string objNameOut = "TFM_INST_"+crossName[ix];
fOut << PPFNameTag(objNameOut) << phaseInstrum;
objNameOut = "TV_INST_"+crossName[ix];
fOut << PPFNameTag(objNameOut) << vphaseInst;
}
}//loop on cross
}//SAve into fOut
TMatrix<r_8>t12;
TMatrix<r_8>t13;
TMatrix<r_8>t14;
TMatrix<r_8>t23;
TMatrix<r_8>t24;
TMatrix<r_8>t34;
{//reload to make closure
PInPersist pIn(fOutName);
#if TYPE_POLAR == 1
pIn >> PPFNameTag(string("TFM_INST_1H2H")) >> t12;
pIn >> PPFNameTag(string("TFM_INST_1H3H")) >> t13;
pIn >> PPFNameTag(string("TFM_INST_1H4H")) >> t14;
pIn >> PPFNameTag(string("TFM_INST_2H3H")) >> t23;
pIn >> PPFNameTag(string("TFM_INST_2H4H")) >> t24;
pIn >> PPFNameTag(string("TFM_INST_3H4H")) >> t34;
#else
pIn >> PPFNameTag(string("TFM_INST_1V2V")) >> t12;
pIn >> PPFNameTag(string("TFM_INST_1V3V")) >> t13;
pIn >> PPFNameTag(string("TFM_INST_1V4V")) >> t14;
pIn >> PPFNameTag(string("TFM_INST_2V3V")) >> t23;
pIn >> PPFNameTag(string("TFM_INST_2V4V")) >> t24;
pIn >> PPFNameTag(string("TFM_INST_3V4V")) >> t34;
#endif
}
{//save closure
string fname = outputDir + "/" + "TFM-ClosureCheck-" + source +".ppf";
POutPersist pOut(fname);
TMatrix<r_8> close1; close1 = t12-t13+t23; Modulo(close1);
TMatrix<r_8> close2; close2 = t12-t14+t24; Modulo(close2);
TMatrix<r_8> close3; close3 = t23-t24+t34; Modulo(close3);
TMatrix<r_8> close4; close4 = t13-t14+t34; Modulo(close4);
pOut << PPFNameTag("CLOS1") << close1;
pOut << PPFNameTag("CLOS2") << close2;
pOut << PPFNameTag("CLOS3") << close3;
pOut << PPFNameTag("CLOS4") << close4;
}
}//BaselinesPhase
//-------------------------------- MAIN --------------------------------------
......
......@@ -292,11 +292,34 @@ void TransformFitData(GeneralFitData& data, std::string tag=""){
idy.erase(idy.begin() + i);
} else {
//Try to see if points before and after a transition point have a difference of phase close to 2pi
r_8 diffPhase = fabs(y[j+1]-y[j-1]);
if(diffPhase < 0.5*(2*M_PI)){
cout << "Transition not validated: diff phase: " << diffPhase << endl;
cout << "y[j-1], y[j], y[j+1]: "
<< y[j-1] << ", " << y[j] << ", " << y[j+1] <<endl;
r_8 diff1= fabs(y[j+1]-y[j-1]);
if( diff1 < 0.5*(2*M_PI)){
cout << "Transition not validated: diff (1) phase: " << diff1 << endl;
idy.erase(idy.begin() + i);
}
}//point not too close to both ends
}//loop on possible transition points
//2nd check
for(int i=idy.size()-1;i>=0;i--){
int j = idy[i];
if(j>2 && j< y.size()-2){
r_8 meanBefore = (y[j-3]+y[j-2])*0.5;
r_8 meanAfter = (y[j+2]+y[j+1])*0.5;
cout << "y[j-2], y[j-1], y[j], y[j+1], y[j+2] "
<< y[j-2] << ", "
<< y[j-1] << ", "
<< y[j] << ", "
<< y[j+1] << ", "
<< y[j+2] <<endl;
cout << "means: " << meanBefore << ", " << meanAfter << endl;
if(fabs(meanAfter-meanBefore)<0.5*(2*M_PI)) {
cout << "Transition not validated: before/after phase: "
<< meanBefore << ", " << meanAfter << endl;
idy.erase(idy.begin() + i);
}
}
}
......@@ -307,6 +330,8 @@ void TransformFitData(GeneralFitData& data, std::string tag=""){
cout << ">" << endl;
{
string outputDir = param.outputDir;
string fNaneDbg = outputDir + "/redpts-" + tag + ".txt";
......
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