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

modif transit approx. usage de TimeVec pour obtenir une premiere plage des...

modif transit approx. usage de TimeVec pour obtenir une premiere plage des indices de RAVec pour la selection des transits quand il y a plusieurs jours et plusieurs sources par jour.
parent 6fae0b71
......@@ -72,7 +72,7 @@ int nloop = 0;
#define POSXYZ_FIXED 1
//JEC 20/7/18 auto-corr azimut contraint >0 (en general non => 0)
#define AZIMUT_FIXED 1
#define AZIMUT_FIXED 0
......@@ -295,6 +295,7 @@ struct PARAM {
r_8 decSrcSec;
TimeStamp MeridianDate; //passage au meridien de la source
TimeStamp ReferenceDate; //JEC 18/2/19 reference des temps pour TimeVec
Now LocalGeoTime;
int year0;
......@@ -1068,6 +1069,9 @@ void fitAutoCrossCorr(){
TimeStamp meridianT = param.MeridianDate;
//JEC 18/2/19 new input the reference date HH::MM::SS
TimeStamp referenceDate = param.ReferenceDate;
r_4 hourAngleMax = param.hourAngleMax;
r_4 hourAngleMin = param.hourAngleMin;
......@@ -1107,14 +1111,14 @@ void fitAutoCrossCorr(){
// 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
//Geometrical initial values (units = meters)
posXCh[0] = 0.;
posYCh[0] = 0.;
posZCh[0] = 0.;
posXCh[1] = 0.;
posYCh[1] = 5.993;
posZCh[1] = 0.;
posZCh[1] = 0.; //was 0.100 pour tester => phase se translate
posXCh[2] = -5.996;
posYCh[2] = -4.380;
......@@ -1272,18 +1276,53 @@ void fitAutoCrossCorr(){
//RAVec
//JEC 18/2/19 : attention we have during long data acquisition many [0h,24h] a cyclic RAVec
//2: select indices according to TimeVec: #seconds since a 0h0 of a given day
TArray<r_8> RAVec;
fIn >> PPFNameTag("RAVec") >> RAVec;
TArray<r_8> TimeVec;
fIn >> PPFNameTag("TimeVec") >> TimeVec;
//number of second since reference date of source meridian transit
r_8 tOfMeridianSinceRef = TimeStamp::TimeDifferenceSeconds(meridianT,referenceDate);
//select +/- 5h arround
r_8 rangeOfTime = 5.0*3600.; //sec
sa_size_t ifirstTime=0;
sa_size_t ilastTime=0;
{
bool first=true;
bool last=true;
for(sa_size_t i=0; i<TimeVec.Size();i++){
if(fabs(TimeVec(i)-tOfMeridianSinceRef)<=rangeOfTime){
if(first){
first = false;
ifirstTime = i;
last = true;
}
if(!first && last)ilastTime=i;
} else {
last = false;
}
}
}//find ifirstTime/ilastTime
if(ifirstTime==ilastTime)
throw RangeCheckError("Time selection not valid");
cout << "Select Time/RA range in index [" << ifirstTime << ", " << ilastTime << "]" << endl;
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
}
//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;
......@@ -1292,7 +1331,10 @@ void fitAutoCrossCorr(){
//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++){
//JEC 18/2/19 restriction of the good HA range according to the ifirstTime/ilastTime range
//was for(sa_size_t i=0; i<HAVec.Size();i++){
for(sa_size_t i=ifirstTime; i<ilastTime;i++){
if(HAVec(i)<haMAX && HAVec(i)>haMin){
if(first){
first = false;
......@@ -2889,25 +2931,38 @@ int main(int narg, char* arg[]) {
string source="CasA31oct17";
string meridianDate="2017-10-31T20:00:00"; //passage au meridien format YYYY-MM-DDTHH:MM:SS
//JEC 18/2/19 reference date
string referenceDate="2019-01-14T00:00:00"; //YYYY-MM-DDTHH:MM:SS
string inputFile="";
string inputParamFile="";
string outputFile="";
//CasA RA-DEC
r_8 raSrcHour = 23.; //RA H:M:S.S
r_8 raSrcMin = 23.;
r_8 raSrcSec = 27.9;
r_8 decSrcDeg = 58; //DEC D:M:S.S
r_8 decSrcMin = 48;
r_8 decSrcSec = 42;
r_4 raSrcHour = 23.; //RA H:M:S.S
r_4 raSrcMin = 23.;
r_4 raSrcSec = 27.9;
r_4 decSrcDeg = 58; //DEC D:M:S.S
r_4 decSrcMin = 48;
r_4 decSrcSec = 42;
//Crabe RA-DEC
// r_4 raSrcHour = 5.; //RA H:M:S.S
// r_4 raSrcMin = 34.;
// r_4 raSrcSec = 31.94;
// r_4 decSrcDeg = 22; //DEC D:M:S.S
// r_4 decSrcMin = 0;
// r_4 decSrcSec = 52.2;
//CygA RA-DEC
// r_8 raSrcHour = 19.; //RA H:M:S.S
// r_8 raSrcMin = 59.;
// r_8 raSrcSec = 28.3;
// r_8 decSrcDeg = 40.; //DEC D:M:S.S
// r_8 decSrcMin = 44.;
// r_8 decSrcSec = 2.0;
// r_4 raSrcHour = 19.; //RA H:M:S.S
// r_4 raSrcMin = 59.;
// r_4 raSrcSec = 28.3;
// r_4 decSrcDeg = 40.; //DEC D:M:S.S
// r_4 decSrcMin = 44.;
// r_4 decSrcSec = 2.0;
r_4 FreqMin = 1250.0; //MHz //Freq minimal, maxi, step OR Freq. centrale + bandwidth
r_4 FreqMax = 1500.0; //MHz
......@@ -2929,7 +2984,9 @@ int main(int narg, char* arg[]) {
int ka=1;
while (ka<narg) {
cout << "option " << arg[ka] << " process...." << endl;
if (strcmp(arg[ka],"-h")==0) {
cout << "usage ex: ./bin/fringeV2 -hasel -10,10 >& tmp\n";
cout << "options: \n";
......@@ -2944,6 +3001,7 @@ int main(int narg, char* arg[]) {
cout << "-ra: Right Ascension (J2000 H:M:S.S) (ex. 23:23:27.9)\n";
cout << "-dec: Declinaison (J2000 D:M:S.S) (ex. 58:48:42.0)\n";
cout << "-meridian: approximate meridian date (ex. 2017-10-18T21:00:00)\n";
cout << "-refdate: date de reference (ex. 2019-01-14T00:00:00)\n";
cout << "-f: freq. range MHz freqMin,freqMax,freqBand (ex. 1300,1450,1): notice freqBand not used\n";
cout << "-hasel: hour angle selected in [hamin,hamax] (Degree): (ex. -10,10)\n";
cout << "-type: type de cross REal or IMag [RE] (WARNING PARTIE_REELLE define sould be tuned accordingly for the moment)\n";
......@@ -2986,11 +3044,14 @@ int main(int narg, char* arg[]) {
ka+=2;
}
else if (strcmp(arg[ka],"-ra")==0) {
scanf(arg[ka+1],"%f:%f:%f",&raSrcHour,&raSrcMin,&raSrcSec);
cout << "option -ra <" << arg[ka+1] << ">" << endl;
sscanf(arg[ka+1],"%f:%f:%f",&raSrcHour,&raSrcMin,&raSrcSec);
cout << raSrcHour << "," << raSrcMin<< "," << raSrcSec << endl;
ka+=2;
}
else if (strcmp(arg[ka],"-dec")==0) {
scanf(arg[ka+1],"%f:%f:%f",&decSrcDeg,&decSrcMin,&decSrcSec);
cout << "option -dec <" << arg[ka+1] << ">" << endl;
sscanf(arg[ka+1],"%f:%f:%f",&decSrcDeg,&decSrcMin,&decSrcSec);
ka+=2;
}
else if (strcmp(arg[ka],"-f")==0){
......@@ -3002,6 +3063,10 @@ int main(int narg, char* arg[]) {
meridianDate=arg[ka+1];
ka+=2;
}
else if (strcmp(arg[ka],"-refdate")==0){
referenceDate=arg[ka+1];
ka+=2;
}
else if (strcmp(arg[ka],"-hasel")==0){
sscanf(arg[ka+1],"%f,%f",&hourAngleMin,&hourAngleMax);
ka+=2;
......@@ -3046,7 +3111,7 @@ int main(int narg, char* arg[]) {
}
//Computation of the exact median passage
r_8 raSrc = HdecfrHMS(raSrcHour,raSrcMin,raSrcSec); //J2000
r_8 raSrc = HdecfrHMS((int)raSrcHour,(int)raSrcMin,(r_8)raSrcSec); //J2000
Now lGeoTime;
//Old value
......@@ -3093,76 +3158,80 @@ int main(int narg, char* arg[]) {
//calcul du passage au meridien
r_8 day0 = iday0 +HdecfrHMS(0,0,0)/24.; //noon
r_8 day0 = iday0 +HdecfrHMS(hr0,0,0)/24.; //noon
// { //OLD
// double dt = 0;
// double lst;
// #define MAXLOOPS 10
// #define MAXERR (1./3600.) /* hours */
// int i =0;
// do {
// day0 += dt/24.;
// lGeoTime.n_mjd = MJDfrDate(day0,month0,year0);
// // calcul passage meridien
// now_lst(&lGeoTime,&lst);
// dt = raSrc-lst;
// cout << "METHOD OLD: meridien dt = " << dt << " " << ToStringHdec(dt) << endl;
// } while (++i < MAXLOOPS && fabs(dt) > MAXERR);
// if (i == MAXLOOPS) throw NotFoundExc("Meridian date not found");
{ //OLD
double dt = 0;
double lst;
#define MAXLOOPS 10
#define MAXERR (1./3600.) /* hours */
int i =0;
do {
day0 += dt/24.;
lGeoTime.n_mjd = MJDfrDate(day0,month0,year0);
// calcul passage meridien
now_lst(&lGeoTime,&lst);
dt = raSrc-lst;
cout << "METHOD OLD: meridien dt = " << dt << " " << ToStringHdec(dt) << endl;
} while (++i < MAXLOOPS && fabs(dt) > MAXERR);
if (i == MAXLOOPS) throw NotFoundExc("Meridian date not found");
// #undef MAXLOOPS
// #undef MAXERR
#undef MAXLOOPS
#undef MAXERR
// meridianT.SetHour(ToStringHdec((day0-iday0)*24.));
meridianT.SetHour(ToStringHdec((day0-iday0)*24.));
// }//OLD
}//OLD
// cout << "OLD passage ay meridien " << meridianT << endl;
{//NEW
// {//NEW
double dt = 0;
#define MAXLOOPS 10
#define MAXERR (1./3600.) /* hours */
int i =0;
// double dt = 0;
// #define MAXLOOPS 10
// #define MAXERR (1./3600.) /* hours */
// int i =0;
r_8 hr = 0;
// r_8 hr = 0;
JMA day;
day.Annee = year0;
day.Mois = month0;
day.Jour = iday0;
do {
hr += dt;
HMS hour = HtoHMS(hr);
HMS tsid;
tsid = TUtoTSid(day, hour);
r_8 ra = HMStoH(tsid); // le ra si lobservatoire etait a Greenwich (longitude 0 deg)
ra += lGeoTime.n_lng/degrad(15.);
// JMA day;
// day.Annee = year0;
// day.Mois = month0;
// day.Jour = iday0;
// do {
// hr += dt;
// HMS hour = HtoHMS(hr);
// HMS tsid;
// tsid = TUtoTSid(day, hour);
// r_8 ra = HMStoH(tsid); // le ra si lobservatoire etait a Greenwich (longitude 0 deg)
// ra += lGeoTime.n_lng/degrad(15.);
dt = raSrc-ra;
cout << "METHOD NEW: meridien dt = " << dt << " " << ToStringHdec(dt) << endl;
// dt = raSrc-ra;
// cout << "METHOD NEW: meridien dt = " << dt << " " << ToStringHdec(dt) << endl;
} while (++i < MAXLOOPS && fabs(dt) > MAXERR);
if (i == MAXLOOPS) throw NotFoundExc("Meridian date not found");
// } while (++i < MAXLOOPS && fabs(dt) > MAXERR);
// if (i == MAXLOOPS) throw NotFoundExc("Meridian date not found");
#undef MAXLOOPS
#undef MAXERR
// #undef MAXLOOPS
// #undef MAXERR
meridianT.SetHour(ToStringHdec(hr));
// meridianT.SetHour(ToStringHdec(hr));
}//NEW
// }//NEW
param.MeridianDate = meridianT;
//JEC 18/2/19
TimeStamp referenceT(referenceDate);
param.ReferenceDate = referenceT;
cout << "APana initial parameter: \n"
<< "Input dir: " << param.inputDir << "\n"
......@@ -3184,7 +3253,7 @@ int main(int narg, char* arg[]) {
<< param.raSrcHour << ":" << param.raSrcMin << ":" << param.raSrcSec << "), DEC d:m:s ("
<< param.decSrcDeg << ":" << param.decSrcMin << ":" << param.decSrcSec << ")\n"
<< "Passage au meridien Nancay a: " << param.MeridianDate << endl;
cout << "Reference date a: " << param.ReferenceDate << endl; //JEC 18/2/19
cout << " Action = " << action
<< " Type of cross REal or IMag: " << param.Type
<< endl;
......
......@@ -3,8 +3,8 @@
include $(SOPHYABASE)/include/sophyamake.inc
### CCIN2P3 MINUITDIR = /sps/baoradio/JEC/PAON4/Minuit2-5.34.14/
MINUITDIR = /Users/campagne/Travail/kits/Minuit2-5.34.14/
MINUITDIR = /sps/baoradio/JEC/PAON4/Minuit2-5.34.14/
### MacJEC MINUITDIR = /Users/campagne/Travail/kits/Minuit2-5.34.14/
MINUITLIB = -L$(MINUITDIR)/lib
MINUITINC = -I$(MINUITDIR)/include
......
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