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

nouveau calcul du temps de passage au meridien

parent 0a520dda
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;
......@@ -1125,9 +1126,9 @@ void fitAutoCrossCorr(){
//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);
// TArray<r_8> TimeVec;
// fIn >> PPFNameTag("TimeVec") >> TimeVec;
// TimeVec.Show(cout);
//RAVec
......@@ -2621,8 +2622,8 @@ int main(int narg, char* arg[]) {
// r_8 decSrcMin = 44.;
// r_8 decSrcSec = 2.0;
r_4 FreqMin = 1255.0; //MHz //Freq minimal, maxi, step OR Freq. centrale + bandwidth
r_4 FreqMax = 1485.0; //MHz
r_4 FreqMin = 1260.0; //MHz //Freq minimal, maxi, step OR Freq. centrale + bandwidth
r_4 FreqMax = 1480.0; //MHz
#if DEBUG_THRESH > 200
......@@ -2761,8 +2762,13 @@ int main(int narg, char* arg[]) {
r_8 raSrc = HdecfrHMS(raSrcHour,raSrcMin,raSrcSec); //J2000
Now lGeoTime;
lGeoTime.n_lat = degrad(47.+21./60.+2./3600.);
lGeoTime.n_lng = degrad(2.+11./60.+43./3600.);
//Old value
// lGeoTime.n_lat = degrad(47.+21./60.+2./3600.);
// lGeoTime.n_lng = degrad(2.+11./60.+43./3600.);
//NEW value JEC 14/6/18 (change very minor)
lGeoTime.n_lat = degrad(47.+22./60.+55.1/3600.);
lGeoTime.n_lng = degrad(2.+11./60.+58.7/3600.);
lGeoTime.n_tz = 0.0; // temps en UTC apres
lGeoTime.n_temp = 0.; //Unused
lGeoTime.n_pressure = 0.; //Unused
......@@ -2803,26 +2809,69 @@ int main(int narg, char* arg[]) {
r_8 day0 = iday0 +HdecfrHMS(0,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");
// #undef MAXLOOPS
// #undef MAXERR
// meridianT.SetHour(ToStringHdec((day0-iday0)*24.));
double dt = 0;
double lst;
// }//OLD
{//NEW
double dt = 0;
#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 << "meridien dt = " << dt << " " << ToStringHdec(dt) << endl;
} while (++i < MAXLOOPS && fabs(dt) > MAXERR);
if (i == MAXLOOPS) throw NotFoundExc("Meridian date not found");
int i =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.);
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");
#undef MAXLOOPS
#undef MAXERR
meridianT.SetHour(ToStringHdec((day0-iday0)*24.));
meridianT.SetHour(ToStringHdec(hr));
}//NEW
param.MeridianDate = meridianT;
......
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