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

Je ne sais pas ce que j'ai modifie sur mon mac de bureau - Reza

Merge branch 'master' of gitlab.in2p3.fr:baoradio/AnaPAON4
parents df765ff6 d193defb
......@@ -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