Commit 189ea7e8 authored by perdereau's avatar perdereau
Browse files
parents 3b848e3f 6903cf30
......@@ -8,8 +8,11 @@ EXE = ./Objs/
SGP4INC = ${EXTLIBDIR}/Satellites/sgp4/libsgp4/
SGP4LIB = ${EXTLIBDIR}/Satellites/sgp4/libsgp4/
#Flag de compilation optionnel
SGP4CCFLG =
# Define our target list
all : predictsatsgp4
all : predictsatsgp4 trk2dt
clean :
rm -f $(EXE)/predictsatsgp4
......@@ -17,7 +20,11 @@ clean :
### Compilation de .o
$(OBJ)/predictsatsgp4.o : predictsatsgp4.cc $(SGP4INC)/Tle.h
$(CXXCOMPILE) -o $(OBJ)/predictsatsgp4.o -I$(SGP4INC) predictsatsgp4.cc
$(CXXCOMPILE) $(SGP4CCFLG) -o $(OBJ)/predictsatsgp4.o -I$(SGP4INC) predictsatsgp4.cc
$(OBJ)/trk2dt.o : trk2dt.cc
$(CXXCOMPILE) -o $(OBJ)/trk2dt.o trk2dt.cc
###### Compilation et link des executables
predictsatsgp4 : $(EXE)/predictsatsgp4
......@@ -26,3 +33,10 @@ predictsatsgp4 : $(EXE)/predictsatsgp4
$(EXE)/predictsatsgp4 : $(OBJ)/predictsatsgp4.o
$(CXXLINK) -o $(EXE)/predictsatsgp4 $(OBJ)/predictsatsgp4.o -L$(SGP4LIB) -lsgp4 -lstdc++ -lc -lm
###### Compilation et link des executables
trk2dt : $(EXE)/trk2dt
echo 'trk2dt --- made'
$(EXE)/trk2dt : $(OBJ)/trk2dt.o
$(CXXLINK) -o $(EXE)/trk2dt $(OBJ)/trk2dt.o $(SOPHYAEXTSLBLIST)
......@@ -34,6 +34,21 @@
SpaceTrack Report No. 3 -> spacetrk.pdf
Revisiting Spacetrack Report #3 -> AIAA-2006-6753-Rev2.pdf
.Notes Reza 25 Nov 2018
o Voir aussi le site web http://celestrak.com
et les liens sur ce site, en particulier le livre http://celestrak.com/software/vallado-sw.php
o Pour faire le build sur mon Mac , aprs avoir tlcharger sgp4.git ,
j'ai fait :
> cd $HERESGP4/
> cmake sgp4
-> library: $HERESGP4/libsgp4/libsgp4.a
-> includes: $HERESGP4/sgp4/libsgp4/*.h
A noter aussi que j'ai du modifier les trois fichiers suivants pour enlever la rfrence a librt.a (supprimer -lrt)
sattrack/CMakeFiles/sattrack.dir/link.txt
runtest/CMakeFiles/runtest.dir/link.txt
passpredict/CMakeFiles/passpredict.dir/link.txt
-----
2-/ Step 2 : telecharger les TLE (two-line element sets)
......@@ -57,6 +72,13 @@ Les satellites dont on telecharge les TLE sont listes dans les variables: naviga
.pour changer l'emplacement, ajouter a l'ordre "make ..." les arguments:
OBJ="la ou je veux mettre mes objets" EXE="la ou je veux mettre mes executables"
.Reza 25/11/2018 : Avec le build comme indique ci-dessus, j'ai du faire
> make SGP4INC=${HERESGP4}/sgp4/libsgp4/ SGP4LIB=${HERESGP4}/libsgp4/
.Reza 25/11/2018 : Au CC-IN2P3, il faut ajouter l'option de compilation -std=c++11
J'ai donc modifie le Makefile en ajoutant un flag optionnel de compilation - pour faire make au CC
> set HERESGP4 = /pbs/throng/baoradio/Library/SGP4_Code/
> make SGP4INC=${HERESGP4}/sgp4/libsgp4/ SGP4LIB=${HERESGP4}/libsgp4/ SGP4CCFLG="-std=c++11"
-----
4-/ Step 4 : executer le programme de prediction
-----
......@@ -69,9 +91,17 @@ Les satellites dont on telecharge les TLE sont listes dans les variables: naviga
debug for satellite GSAT0102 (PRN E12) i.e. "37847"
> ./Objs/predictsatsgp4 -D 37847 satname YYYMMDD/galileo.txt [...]
.Reza 25/11/2018 : j'ai ajoute l'option -K au programme, pour creer un fichier avec la trace (track) du satellite
..../AnaPAON4/Satellites/Objs/predictsatsgp4 -T "2018/10/20 15:16:00" -H 0.,81. -K 41859 -p 2 galileo.txt
-----
5-/ Fichiers trace
-----
.le programme predictsatsgp4 peut crer des fichiers de trace (track) pour satellite slectionn et pour le solei.
On peut utiliser le programme trk2dt.cc pour lire le fichier trace et le convertir en DataTable (fichier PPF, ou FITS)
-----
5-/ Remarques
6-/ Remarques
-----
.Est aussi fourni un petit job de test (testpredsatsgp4.py) python (2.7) utilisant le module "sgp4"
https://pypi.python.org/pypi/sgp4/
......
......@@ -23,6 +23,9 @@ void ReadFile(string infile,vector<string>& vline1,vector<string>& vline2,vector
Vector AzAlt2Vec(double az,double alt);
double SepAngle(Vector& v1, Vector& v2);
void DelAzAltSmall(const CoordTopocentric& topo1,const CoordTopocentric& topo2,double &daz,double &dalt);
void SaveTrack(DateTime & dateobs, string & satname, SGP4 & sgp4, Observer & obs, DateTime& datestart,
DateTime& dateend, TimeSpan& tspaninc, int prtlev=0);
void SaveSunTrack(DateTime & dateobs, Observer & obs, int sincmin, int prtlev=0);
void usage(void);
void usage(void)
......@@ -33,10 +36,12 @@ cout<<"predictsatsgp4 [options] filetle1.txt [filetle2.txt ...] : find satellite
<<" -H azimuth(dd.ddd,N=0,E->W),altitude(dd.ddd,Z+) : (def=180d,83.40d)"<<endl
<<" telescope horizontal coordinates pointing"<<endl
<<" -D satname : print debug for this satellite (even if not in the search area)"<<endl
<<" -K satname : Create a track file for this satellite (same time span as the one defined for scan)"<<endl
<<" -S sunincmin : Create the sun track file over 24 hours for obs. date, sun pos every sunincmin minutes"<<endl
<<" -A scanlos(dd.ddd) : search area around los (def=10d)"<<endl
<<" find satellites closer than \"scanlos\" degrees of los"<<endl
<<" -s spanhour(int),spanmin(min),spaninc(min) :"<<endl
<<" for selected satellites (-A,-D) accuratly scan at +/-(hh:mm) with step \"spaninc\" (def=1h,30mn,1mn)"<<endl
<<" -s spanhour(int),spanmin(min),spanincm(min),spanics(sec) :"<<endl
<<" for selected satellites (-A,-D) accuratly scan at +/-(hh:mm) with step \"spaninc\" (def=1h,30mn,1mn,0s)"<<endl
<<" -v dang4vel : angular increment (deg) to compute the angular velocity of the satellite (def=1.)"<<endl
<<" -p lp : print debug"<<endl
<<endl;
......@@ -57,17 +62,21 @@ int main(int narg, char *arg[])
//Scan area around los (deg)
double scanlos=10.0; //degres (CygA a Paon4 au transit sud du 5/10/2017 ~18:52:45)
//Scan autour de +/- l'heure d'observation et increment de scan (en minute)
int dtspanhour=1, dtspanmin = 30, satincmin = 1;
int dtspanhour=1, dtspanmin = 30, satincmin = 1, satincsec = 0;
//Time increment (min) to compute the angular velocity of the satellite
double dang4vel = 1.;
//Satellite name to debug
string debugsatname;
//Satellite name for which to compute and save the track
string tracksatname;
//time increment in minutes for creating sun track file (create track file if sunincmin>0)
int sunincmin = 0;
//print level
int lp=0;
//--- Decodage des arguments
char c;
while((c = getopt(narg,arg,"hT:O:H:A:s:D:v:p:")) != -1) {
while((c = getopt(narg,arg,"hT:O:H:A:s:D:K:S:v:p:")) != -1) {
switch (c) {
case 'T' :
datestr = optarg;
......@@ -82,11 +91,17 @@ int main(int narg, char *arg[])
scanlos = atof(optarg);
break;
case 's' :
sscanf(optarg,"%u,%u,%u",&dtspanhour,&dtspanmin,&satincmin);
sscanf(optarg,"%u,%u,%u,%u",&dtspanhour,&dtspanmin,&satincmin,&satincsec);
break;
case 'D' :
debugsatname = optarg;
break;
case 'K' :
tracksatname = optarg;
break;
case 'S' :
sunincmin = atoi(optarg);
break;
case 'v' :
dang4vel = atof(optarg);
break;
......@@ -126,8 +141,13 @@ int main(int narg, char *arg[])
//--- Scan autour de +/- l'heure d'observation
TimeSpan tspan(dtspanhour,dtspanmin,0);
if(satincmin<=0) satincmin = 1;
TimeSpan tspaninc(0,satincmin,0);
if(satincmin<=0) {
satincmin = 0;
if(satincsec<=0) satincmin=1;
}
if(satincsec<=0) satincsec=0;
TimeSpan tspaninc(0,satincmin,satincsec);
cout<<"TSpan: "<<tspan<<" / inc="<<tspaninc<<endl;
DateTime datestart = dateobs - tspan;
DateTime dateend = dateobs + tspan;
......@@ -139,7 +159,9 @@ int main(int narg, char *arg[])
//--- Debug satellite:
if(debugsatname.size()>0) cout<<"Debug satellite: "<<debugsatname<<endl;
//--- Debug satellite:
if(tracksatname.size()>0) cout<<"Creating Track file for satellite: "<<tracksatname<<endl;
//--- Sun (ephemerides de la journee d'observation)
cout<<endl;
DateTime dt0(dateobs.Year(),dateobs.Month(),dateobs.Day(),0,0,0);
......@@ -152,7 +174,10 @@ int main(int narg, char *arg[])
if(RadiansToDegrees(toposun.elevation) < -18.0) continue;
cout <<"Sun: "<<dt<<" "<<toposun<<" "<<geosun<<endl;
}
if (sunincmin>0) {
cout<<"***SUN-Creating Track file for SUN over 24 hours with every "<<sunincmin<<" minutes..."<<endl;
SaveSunTrack(dateobs, obspaon4, sunincmin, lp);
}
//--- Read TLE parameters for all satellite files and fill TLE into vector<>
cout<<endl;
vector<string> vline1, vline2, vsatfilename;
......@@ -173,6 +198,7 @@ int main(int narg, char *arg[])
string satname(strsatname);
Tle tle = Tle(satname,vline1[isat],vline2[isat]);
long ifounddbg = (debugsatname.size()>0) ? satname.find(debugsatname,0): -1;
long ifoundtrk = (tracksatname.size()>0) ? satname.find(tracksatname,0): -1;
//
try {
SGP4 sgp4(tle);
......@@ -230,6 +256,10 @@ int main(int narg, char *arg[])
if(s<sepmin) {sepmin = s; datemin = datecur; topomin = toposep;}
datecur = datecur + tspaninc;
}
if(ifoundtrk==0) { // Create track file for this satellite
cout<<"***TRK-Creating Track file for SAT="<<satname<<endl;
SaveTrack(dateobs, satname, sgp4, obspaon4, datestart, dateend, tspaninc, lp);
}
//...Print results
//INFO: "Rng" is the "slant range" (distance observer-satellite) in km
// and "Rng Rt" is the "slant range rate" in km/s
......@@ -357,3 +387,90 @@ void DelAzAltSmall(const CoordTopocentric& topo1,const CoordTopocentric& topo2,d
if(daz > +180.) daz -= 360.; //ex: az1=10d, az2=350d, az2-az1=350-10=+340d -> daz=+340-360=-20d
if(daz <= -180.) daz += 360.; //ex: az1=350d, az2=10d, az2-az1=10-350=-340d -> daz=-340+360=+20d
}
//---------------------------------------------------------
void SaveTrack(DateTime & dateobs, string & satname, SGP4 & sgp4, Observer & obs,
DateTime& datestart, DateTime& dateend, TimeSpan& tspaninc, int prtlev)
{
char flnm[128];
sprintf(flnm,"trk_%s_%4d%2d%2d.txt",satname.c_str(),(int)dateobs.Year(), (int)dateobs.Month(), (int)dateobs.Day());
ofstream ofs(flnm);
CoordGeodetic obsgeo=obs.GetLocation();
ofs << "#### Track for satellite : " << satname << " Date (yyyy/mm/dd): "
<<dateobs.Year()<<"/"<<dateobs.Month()<<"/"<<dateobs.Day()<<endl;
ofs << "#### Observations position, latitude,longitude (degrees) : "
<< RadiansToDegrees(obsgeo.latitude)<<" , "<<RadiansToDegrees(obsgeo.longitude) << endl;
ofs << "## Date Time UTC : yyyy-mm-dd hh:mm:ss.ssss "<<endl;
ofs << "## Satellite Azimuth and Elevation , corresponding Latitude and Longitude (in degrees)"<<endl;
ofs << "## Date Time UTC Azimuth Elevation Range Latitude Longitude Altitude (Angles in degree, Range,Altitude in km)" << endl;
if (prtlev>1) {
cout << "#### Track for satellite : " << satname << " Date (yyyy/mm/dd): "
<< dateobs.Year()<<"/"<<dateobs.Month()<<"/"<<dateobs.Day()<<endl;
cout << "## Date Time UTC Azimuth Elevation Range Latitude Longitude Altitude (Angles in degree, Range,Altitude in km)" << endl;
}
DateTime datecur = datestart;
while(datecur <= dateend+tspaninc) {
Eci eci = sgp4.FindPosition(datecur);
//-- Vector Velsat = eci.Velocity();
//-- double velsat = Velsat.Magnitude(); //km/s
//...Get look angle for observer to satellite
CoordTopocentric topo = obs.GetLookAngle(eci);
// Vector Vsat = AzAlt2Vec(RadiansToDegrees(topo.azimuth),RadiansToDegrees(topo.elevation));
//...Convert satellite position to geodetic coordinates
CoordGeodetic geo = eci.ToGeodetic();
// Vector Vungeo = AzAlt2Vec(RadiansToDegrees(geo.longitude),RadiansToDegrees(geo.latitude));
ofs <<datecur<<" "<<setw(10)<<RadiansToDegrees(topo.azimuth)
<<" "<<setw(10)<<RadiansToDegrees(topo.elevation)<<" "<<setw(10)<<topo.range
<<" "<<setw(10)<<RadiansToDegrees(geo.latitude)<<" "<<setw(10)<<RadiansToDegrees(geo.longitude)
<<" "<<setw(10)<<geo.altitude<<endl;
if (prtlev>1) {
cout <<datecur<<" "<<setw(10)<<RadiansToDegrees(topo.azimuth)
<<" "<<setw(10)<<RadiansToDegrees(topo.elevation)<<" "<<setw(10)<<topo.range
<<" "<<setw(10)<<RadiansToDegrees(geo.latitude)<<" "<<setw(10)<<RadiansToDegrees(geo.longitude)
<<" "<<setw(10)<<geo.altitude<<endl;
}
datecur = datecur + tspaninc;
}
return;
}
void SaveSunTrack(DateTime & dateobs, Observer & obs, int sincmin, int prtlev)
{
char flnm[128];
sprintf(flnm,"suntrk_%4d%2d%2d.txt",(int)dateobs.Year(), (int)dateobs.Month(), (int)dateobs.Day());
ofstream ofs(flnm);
CoordGeodetic obsgeo=obs.GetLocation();
ofs << "#### Sun track for Date (yyyy/mm/dd): "
<< dateobs.Year()<<"/"<<dateobs.Month()<<"/"<<dateobs.Day()<<endl;
ofs << "#### Observations position, latitude,longitude (degrees) : "
<< RadiansToDegrees(obsgeo.latitude)<<" , "<<RadiansToDegrees(obsgeo.longitude) << endl;
ofs << "## Date Time UTC Azimuth Elevation Range Latitude Longitude Altitude (Angles in degree, Range,Altitude in km)" << endl;
DateTime dt0(dateobs.Year(),dateobs.Month(),dateobs.Day(),0,0,0);
DateTime dt = dt0;
TimeSpan tsp(0,sincmin,0);
int jmax=(24*60)/sincmin;
SolarPosition sunpos;
//DBG cout << "SUN*DBG* tspan="<<tsp<<" sincmin="<<sincmin<<endl;
for(int j=0;j<=jmax;j++) {
Eci ecisun = sunpos.FindPosition(dt);
CoordTopocentric topo = obs.GetLookAngle(ecisun);
CoordGeodetic geo = ecisun.ToGeodetic();
//DBG cout << "SUN*DBG*"<<dt<<" topo.elevation="<< RadiansToDegrees(topo.elevation)<<endl;
if(RadiansToDegrees(topo.elevation) < -18.0) {
dt = dt + tsp; //inc=sincmin minutes
continue;
}
ofs <<dt<<" "<<setw(10)<<RadiansToDegrees(topo.azimuth)
<<" "<<setw(10)<<RadiansToDegrees(topo.elevation)<<" "<<setw(10)<<topo.range
<<" "<<setw(10)<<RadiansToDegrees(geo.latitude)<<" "<<setw(10)<<RadiansToDegrees(geo.longitude)
<<" "<<setw(10)<<geo.altitude<<endl;
dt = dt + tsp; //inc=sincmin minutes
}
return;
}
#include "machdefs.h"
/* ----------------------------------------------------------
Projet BAORadio/PAON4 - (C) LAL/IRFU 2008-2015
Programme de conversion fichiers ascii track (satellites, soleil)
en DataTable, sauve en PPF (ou fits) Reza, Mai 2018
---------------------------------------------------------- */
// include standard c/c++
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <iostream>
#include <string>
// include SOPHYA
#include "pexceptions.h"
#include "timestamp.h"
#include "histats.h"
#include "datatable.h"
#include "histinit.h"
#include "fitsioserver.h"
using namespace std;
using namespace SOPHYA;
int Usage(void)
{
cout<<"--- trk2dt.cc : Converting track file to DataTable \n"<<endl;
cout<<"Usage: trk2dt track_file.txt outFile [StartDate yyyy/mm/dd] \n";
cout<<" - outFile : PPF or FITS format based on the extension (.ppf or .fits) \n";
cout<<" - StartDate : is used as offset (at 00h00) to compute the timesec value \n";
cout<<" timesec = (DataTime-StarDate_00h00).ToSeconds \n";
cout<<" Default StartDate (if not specified) is extracted from the first DateTime in track file \n";
cout<<endl;
return 1;
}
//----------------------------------------------------
int main(int narg, const char* arg[])
{
int rc=0;
// --- Decoding parameters
if ((narg<3)||((narg>1)&&(strcmp(arg[1],"-h")==0))) return Usage();
string infile = arg[1];
string outfile = arg[2];
cout << " ==== trk2dt reading input track file "<<infile<<" to create DataTable"<<endl
<< " ... saved to output file"<<outfile<<endl;
TimeStamp date0(2018, 1, 1, 0, 0, 0.); // first date
bool fgdate0set=false;
if (narg>3) {
int yy=2018 , mm=1, dd=1;
sscanf(arg[3],"%d/%d/%d",&yy,&mm,&dd);
date0.SetDate(yy,mm,dd);
fgdate0set=true;
cout << "StartDate (Offset for timesec column) set from command line argument ->"<<date0<<endl;
}
try {
DataTable dt(128);
dt.AddDateTimeColumn("datetime");
dt.AddDoubleColumn("timesec");
dt.AddDoubleColumn("azimuth");
dt.AddDoubleColumn("elevation");
dt.AddDoubleColumn("range");
dt.AddDoubleColumn("latitude");
dt.AddDoubleColumn("longitude");
dt.AddDoubleColumn("altitude");
dt.Info().Comment()="Satellite track DataTable from predictsatsgp4.cc";
dt.Info()["datetime"] = "Date/Time";
dt.Info()["timesec"] = "time of day, in seconds";
dt.Info()["azimuth"] = "Azimuth angle , in degree";
dt.Info()["elevation"] = "Elevation angle , in degree";
dt.Info()["range"] = "Range (distance ?) in km";
dt.Info()["latitude"] = "Latitude , in degree";
dt.Info()["longitude"] = "Longitude , in degree";
dt.Info()["altitude"] = "Altitude in km";
int nlines = 0;
int nok = 0;
DataTableRow row = dt.EmptyRow();
ifstream is(infile.c_str());
string sline;
while (!is.eof()) {
sline = "";
getline(is, sline); nlines++;
if (!is.good() && !is.eof()) continue;
if ((sline.length()<69) || (sline[0] == '#')) continue;
int yy=2018 , mm=1, dd=1;
sscanf(sline.c_str(), "%d-%d-%d",&yy,&mm,&dd);
int hh=0 , mn=0;
double ss=0.;
sscanf(sline.c_str()+11, "%d:%d:%lg",&hh,&mn,&ss);
TimeStamp datetm(yy, mm, dd, hh, mn, ss);
if (!fgdate0set) {
date0.SetDate(yy,mm,dd);
fgdate0set=true;
cout << "StartDate (Offset for timesec column) set from first date in track file ->"<<date0<<endl;
}
double azim, elev, range, latitude, longitude, alt;
sscanf(sline.c_str()+30, "%lg %lg %lg %lg %lg %lg", &azim, &elev, &range, &latitude, &longitude, &alt);
double tmsec = TimeStamp::TimeDifferenceSeconds(datetm, date0);
row[0] = datetm;
row[1] = tmsec;
row[2] = azim;
row[3] = elev;
row[4] = range;
row[5] = latitude;
row[6] = longitude;
row[7] = alt;
dt.AddRow(row); nok++;
}
dt.SetShowMinMaxFlag(true);
cout << dt ;
size_t ll = outfile.length();
size_t pp = outfile.rfind('.');
if ((pp<ll)&&(outfile.substr(pp)==".fits")) {
cout << "----- Saving track DataTable to FITS file "<<outfile<<endl;
FitsInOutFile fos(outfile,FitsInOutFile::Fits_Create);
fos << dt;
}
else {
cout << "----- Saving track DataTable to PPF file "<<outfile<<endl;
POutPersist pof(outfile);
pof << dt;
}
}
catch (PException& exc) {
cerr << " rdvisip4.cc catched PException " << exc.Msg() << endl;
rc = 77;
}
catch (std::exception& sex) {
cerr << "\n trk2dt.cc std::exception :"
<< (string)typeid(sex).name() << "\n msg= " << sex.what() << endl;
rc = 78;
}
catch (...) {
cerr << " trk2dt.cc catched unknown (...) exception " << endl;
rc = 79;
}
cout << ">>>> trk2dt.cc ------- END ----------- RC=" << rc << endl;
return rc;
}
/*
Analyse PAON4 , Premier essais sur les signaux des satellites
R. Ansari, Novembre 2018
Code destine a etre utilise avec anasat.pic
Il faut definir les 4 arguments suivants au niveau de l'interpreteur par c++args
avant l'appel a anasat.cc
FREQ_MIN , FREQ_MAX : definition de la bande de frequence en MHz , 1275-1282 pour Galileo
TIM_START , TIM_END : bornes en temps (en minutes) pour les calculs de normalisations
c++args $FREQ_MIN $FREQ_MAX $TIM_START $TIM_END
c++execfrf anasat.cc
*/
if (args.size() < 4) {
cout << " anasat.cc / Missing arguments FREQ_MIN , FREQ_MAX , TIM_START , TIM_END"<<endl
<< " Define FREQ_MIN , FREQ_MAX , TIM_START , TIM_END "<<endl
<< " Use c++args $FREQ_MIN $FREQ_MAX $TIM_START $TIM_END"<<endl
<< " before calling anasat.cc through c++exec " << endl;
return 9;
}
string FREQ_MIN = args[0];
string FREQ_MAX = args[1];
string TIM_START = args[2];
string TIM_END = args[3];
// Frequency Binning
// int JMIN=51, JMAX=65;
// int JMIN=102, JMAX=132;
double delta_freq=250./(double)TFM_1H.SizeY();
sa_size_t JMIN = (atof(FREQ_MIN.c_str())-1250.+0.5)/delta_freq;
sa_size_t JMAX = (atof(FREQ_MAX.c_str())-1250.+0.5)/delta_freq;
cout << "anasat.cc/Info: FREQ_MIN,MAX="<<FREQ_MIN<<","<<FREQ_MAX<<" -> JMIN,MAX="<<JMIN<<","<<JMAX<<endl;
// Normalisation Time-range
// sa_size_t ITMA=180, ITMB=220;
// sa_size_t ITMA=360, ITMB=440; /* B01 41859 */
// sa_size_t ITMA=1010, ITMB=1070; /* B04 40889 */
double tm_start=TimeVec(0);
double tm_end=TimeVec(TimeVec.Size()-1);
double delta_temps=(tm_end-tm_start)/(double)TimeVec.Size();
sa_size_t ITMA=(atof(TIM_START.c_str())*60.-tm_start+0.5)/delta_temps;
sa_size_t ITMB=(atof(TIM_END.c_str())*60.-tm_start+0.5)/delta_temps;
cout << "anasat.cc/Info: TIM_START,END="<<TIM_START<<","<<TIM_END<<" -> ITMA,B="<<ITMA<<","<<ITMB<<endl;
sa_size_t SZX=TFM_1H.SizeX();
Vector V11(SZX); V11=0.; Vector V22(SZX); V22=0.;
Vector V33(SZX); V33=0.; Vector V44(SZX); V44=0.;
TVector< complex<double> > V12(SZX); V12= complex<double>(0.,0.);
TVector< complex<double> > V13(SZX); V13= complex<double>(0.,0.);
TVector< complex<double> > V14(SZX); V14= complex<double>(0.,0.);
TVector< complex<double> > V23(SZX); V23= complex<double>(0.,0.);
TVector< complex<double> > V24(SZX); V24= complex<double>(0.,0.);
TVector< complex<double> > V34(SZX); V34= complex<double>(0.,0.);
for(sa_size_t i=0; i<V11.Size(); i++) {
for(sa_size_t j=JMIN; j<=JMAX; j++) {
V11(i)+=TFM_1H(i,j); V22(i)+=TFM_2H(i,j);
V33(i)+=TFM_3H(i,j); V44(i)+=TFM_4H(i,j);
V12(i)+=TFM_1H2H(i,j);
V13(i)+=TFM_1H3H(i,j);
V14(i)+=TFM_1H4H(i,j);
V23(i)+=TFM_2H3H(i,j);
V24(i)+=TFM_2H4H(i,j);
V34(i)+=TFM_3H4H(i,j);
}
}
double fnorm=1./(double)(JMAX-JMIN+1);
complex<double> znorm=complex<double>(fnorm,0.);
V11 /= fnorm; V22 /= fnorm;
V33 /= fnorm; V44 /= fnorm;
V12 /= znorm; V13 /= fnorm; V14 /= znorm;
V23 /= znorm; V24 /= fnorm; V34 /= znorm;
double min,max;
vector<double> vfnorm;
V11(Range(ITMA,ITMB)).MinMax(min,max); vfnorm.push_back(max);
cout << " V11->Max="<<max<<endl; V11 /= max;
V22(Range(ITMA,ITMB)).MinMax(min,max); vfnorm.push_back(max);
cout << " V22->Max="<<max<<endl; V22 /= max;
V33(Range(ITMA,ITMB)).MinMax(min,max); vfnorm.push_back(max);
cout << " V33->Max="<<max<<endl; V33 /= max;
V44(Range(ITMA,ITMB)).MinMax(min,max); vfnorm.push_back(max);
cout << " V44->Max="<<max<<endl; V44 /= max;
V12 /= complex<double>(sqrt(vfnorm[0]*vfnorm[1]));
V13 /= complex<double>(sqrt(vfnorm[0]*vfnorm[2]));
V14 /= complex<double>(sqrt(vfnorm[0]*vfnorm[3]));
V23 /= complex<double>(sqrt(vfnorm[1]*vfnorm[2]));
V24 /= complex<double>(sqrt(vfnorm[1]*vfnorm[3]));
V34 /= complex<double>(sqrt(vfnorm[2]*vfnorm[3]));
cout << " Finished normalisation --- Creating DataTable ..."<<endl;
DataTable dtvis(128);
dtvis.AddDoubleColumn("timesec");
dtvis.AddDoubleColumn("ra_hour");
dtvis.AddDoubleColumn("V11");
dtvis.AddDoubleColumn("V22");
dtvis.AddDoubleColumn("V33");
dtvis.AddDoubleColumn("V44");
dtvis.AddDoubleComplexColumn("V12");
dtvis.AddDoubleComplexColumn("V13");
dtvis.AddDoubleComplexColumn("V14");
dtvis.AddDoubleComplexColumn("V23");
dtvis.AddDoubleComplexColumn("V24");
dtvis.AddDoubleComplexColumn("V34");
DataTableRow row = dtvis.EmptyRow();
for(size_t i=0; i<(size_t)SZX; i++) {
row[0] = TimeVec(i);
row[1] = RAVec(i);
row[2] = V11(i); row[3] = V22(i);
row[4] = V33(i); row[5] = V44(i);
row[6] = V12(i); row[7] = V13(i);
row[8] = V14(i); row[9] = V23(i);
row[10] = V24(i); row[11] = V34(i);
dtvis.AddRow(row);
}
dtvis.SetShowMinMaxFlag(true);
cout<<dtvis;
KeepObj(dtvis);
##########################################################################
####### Analyse PAON4 , Premier essais sur les signaux des satellites
####### R. Ansari, Novembre 2018
## Analyse / affichage a partir des cartes TFM Auto et Cross-Correlation
## et DataTable visibilites fait par le programme JSkyMap/trk2vis.cc
## Ce script utilise anasat.cc
##########################################################################
setaxesatt 'font=helvetica,bold,16 fixedfontsize minorticks'
###########################################
###### Part 1 : Analyse B01