Commit 0a520dda authored by Jean-Eric Campagne's avatar Jean-Eric Campagne
Browse files

Deux bugs fix: MJDfrDate a un bug, donc le calcul de angle horaire se fait a...

Deux bugs fix: MJDfrDate a un bug, donc le calcul de angle horaire se fait a partir de RAVec; le calcul du vecteur ns pointant vers la source a un faute de signe.
parent 812b86e3
......@@ -237,7 +237,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 +261,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 +479,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 +926,7 @@ void fitAutoCrossCorr(){
r_8 decSrc = degrad(decSrcDeg+decSrcMin/60.+decSrcSec/3600.); //J2000
cout << "raSrc (hr): " << raSrc << ", decSrc (rad): " << decSrc << "\n";
......@@ -936,6 +940,9 @@ void fitAutoCrossCorr(){
if(fabs(hourAngleMin)<HATHRESHOLD)
throw ParmError("fitAllCross |hourAngleMin|<HATHRESHOLD");
Now LocalGeoTime = param.LocalGeoTime;
//--------------------
// Prepare numbering & naming stuff
//--------------------
......@@ -965,10 +972,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 +1123,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;
//RAVec
TArray<r_8> RAVec;
fIn >> 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";
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 +2063,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
TArray<r_8> RAVec;
fInPhFreeVisi >> PPFNameTag("RAVec") >> RAVec;
fOut << PPFNameTag("RAVec") << RAVec;
}
//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 +2225,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 +2260,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
//---------
//save the phase-free visibilities
//---------
{
string objNameOut = "TFM_BASELINE_"+crossName[ix];
fOut << PPFNameTag(objNameOut) << geomPhase;
}
// //save dataTable
// {
// string objNameOut = "DT_"+crossName[ix];
// fOut << PPFNameTag(objNameOut) << dt;
// }
}//loop on cross
}//BaselinesPhase
//-------------------------------- MAIN --------------------------------------
//-------------------------------- MAIN --------------------------------------
//-------------------------------- MAIN --------------------------------------
......@@ -2270,10 +2621,9 @@ int main(int narg, char* arg[]) {
// r_8 decSrcMin = 44.;
// r_8 decSrcSec = 2.0;
r_4 FreqMin = 1300.0; //MHz //Freq minimal, maxi, step OR Freq. centrale + bandwidth
r_4 FreqMax = 1450.0; //MHz
// r_4 FreqMin = 1349.0; //MHz //Freq minimal, maxi, step OR Freq. centrale + bandwidth
// r_4 FreqMax = 1349.15; //MHz
r_4 FreqMin = 1255.0; //MHz //Freq minimal, maxi, step OR Freq. centrale + bandwidth
r_4 FreqMax = 1485.0; //MHz
#if DEBUG_THRESH > 200
cout << "WARNING !!! High DEBUG LEVEL: " << DEBUG_THRESH << endl;
......@@ -2514,6 +2864,8 @@ int main(int narg, char* arg[]) {
fitAutoCrossCorr();
} else if (action == "rotate") { //remove NS ( & Vertical) phases
rotateVisi();
} else if (action == "basephase") { // compute baseline-only phases TFM map (8/6/18 JEC)
BaselinesPhase();
} else {
throw ParmError("Unknown action. Bye");
}
......
......@@ -20,21 +20,9 @@
#include "poly.h"
// //Minuit
// #include "Math/MinimizerOptions.h"
#include "Minuit2/FCNGradientBase.h"
// #include "Minuit2/FunctionMinimum.h"
// #include "Minuit2/MnUserParameterState.h"
// #include "Minuit2/MnPrint.h"
// #include "Minuit2/MnMigrad.h"
// #include "Minuit2/MnMinimize.h"
// #include "Minuit2/MnMinos.h"
// #include "Minuit2/MnContours.h"
// #include "Minuit2/MnPlot.h"
// #include "Minuit2/MinosError.h"
// #include "Minuit2/ContoursError.h"
using namespace ROOT::Minuit2;
// using namespace ROOT::Math;
//nbre of polar per dish: 1 => default H-type, 2 => V-type
#define TYPE_POLAR 1
......@@ -86,25 +74,6 @@ void split(const string& str, const string& delimiters , vector<string>& tokens)
}
}
//-------
// char *regexp (const char *string, const char *patrn, int *begin, int *end) {
// int i, w=0, len;
// char *word = NULL;
// regex_t rgT;
// regmatch_t match;
// regcomp(&rgT,patrn,REG_EXTENDED);
// if ((regexec(&rgT,string,1,&match,0)) == 0) {
// *begin = (int)match.rm_so;
// *end = (int)match.rm_eo;
// len = *end-*begin;
// word=(char*)malloc(len+1);
// for (i=*begin; i<*end; i++) {
// word[w] = string[i];
// w++; }
// word[w]=0;
// }
// regfree(&rgT);
// return word;
// }
sa_size_t round_sa(r_4 r) {
return static_cast<sa_size_t>((r > 0.0) ? (r + 0.5) : (r - 0.5));
}
......@@ -364,47 +333,6 @@ void TransformFitData(GeneralFitData& data, std::string tag=""){
//----
//--------------------------------- MINUIT func CLASS --------------------
// void DumpMinuitParam(const MnUserParameters& upar) {
// int pr = cout.precision();
// #define PRECISION 13
// #define WIDTH 12
// for(std::vector<MinuitParameter>::const_iterator ipar = upar.Parameters().begin();
// ipar != upar.Parameters().end(); ipar++) {
// cout << std::setw(4) << (*ipar).Number() << std::setw(5) << "||";
// cout << std::setw(10) << (*ipar).Name() << std::setw(3) << "||";
// if((*ipar).IsConst()) {
// cout << " const ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value() << " ||" << std::endl;
// } else if((*ipar).IsFixed()) {
// cout << " fixed ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value() << " ||" << std::endl;
// } else if((*ipar).HasLimits()) {
// if((*ipar).Error() > 0.) {
// cout << " limited ||" << std::setprecision(PRECISION) << std::setw(WIDTH) << (*ipar).Value();
// cout << " || " << std::setw(12) << (*ipar).Error() << " ||";