Skip to content

GitLab

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

Merge branch 'master' of gitlab.in2p3.fr:baoradio/AnaPAON4

parents d34bbc52 080707b4
Hide whitespace changes
Inline Side-by-side
Showing with 2164 additions and 45 deletions
Python/viewtfmap.py
View file @ d193defb
......@@ -10,23 +10,29 @@ from skimage import exposure
import mynormalize as mynorm
from matplotlib.widgets import Slider, Button, RadioButtons
import time
from colorsys import hsv_to_rgb, hls_to_rgb
# Choices: ['auto', 'gtk', 'gtk3', 'inline', 'nbagg', 'notebook', 'osx', 'qt', 'qt4', 'qt5', 'tk', 'wx']
rfreq = [1370.,1390.,1420.5,1430.]
rfreq = [1381.,1390.,1420.5,1430.]
normf=[1325.,1465.]
names = ['1H','2H','3H','4H','1V','2V','3V','4V','1Hx2H','1Hx3H','1Hx4H','2Hx3H','2Hx4H','3Hx4H',
'1Vx2V','1Vx3V','1Vx4V','2Vx3V','2Vx4V','3Vx4V','1Hx1V','1Hx2V','1Hx3V','1Hx4V','2Hx1V','2Hx2V','2Hx3V','2Hx4V',
'3Hx1V','3Hx2V','3Hx3V','3Hx4V','4Hx1V','4Hx2V','4Hx3V','4Hx4V']
def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=False,mtitle="",clean=True,mod=False,arg=False,noplo=False,verb=False,yrcm=False,cor=[],cimag=False,creal=False,hasvar=True):
print len(cor)
def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=False,mtitle="",clean=True,mod=False,
arg=False,cmplx=False,noplo=False,verb=False,yrcm=False,cor=[],cimag=False,creal=False,hasvar=True,
ramod=False,plras=[],nmras=[],mask=False):
if verb :
print len(cor)
filename = folder
#print i
#print num
if noplo :
clean=False
if clean:
plt.close("all")
hdulist = fits.open(folder+'.fits')
winwid = (18,8)
......@@ -34,7 +40,10 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
times=hdulist[len(hdulist)-1].data
ras=hdulist[len(hdulist)-2].data #41 in 43
freqs= hdulist[len(hdulist)-3].data
if ramod :
ras=hdulist[len(hdulist)-1].data #41 in 43
freqs= hdulist[len(hdulist)-2].data
i=2*num
if num>7 :
i= 16 + 4*(num-8)
......@@ -51,13 +60,19 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
print hdulist[i].header
bfreq = [ (np.abs(freqs-rf)).argmin() for rf in rfreq ]
print ' number of hds ' + str(len(hdulist))
bfnor = [(np.abs(freqs-rf)).argmin() for rf in normf ]
if (verb) :
print ' number of hds ' + str(len(hdulist))
nfre = np.size(freqs)
ntim = np.size(times)
timin = min(times)/3600.
if timin ==0. :
timin=times[0]/3600.
timax = max(times)/3600.
if ramod :
timin = min(ras)
timax= max(ras)
frmin = min(freqs)
frmax=max(freqs)
if (not noplo):
......@@ -66,7 +81,9 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
img = hdulist[ i].data
if(cimag):
img = hdulist[ i+1].data
print np.shape(img)
if (verb) :
print 'Image shape :'
print np.shape(img)
if mod :
if num>7 :
......@@ -84,14 +101,47 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
img = np.arctan2(imgi,img) * 180. / np.pi
mvmin=-180.
mvmax=180.
mtitle= names[num]+" arg(deg)"
mtitle= names[num]+" arg.(deg)"
else :
print "arg demande mais auto selectionne : ERREUR "
exit
vmin = np.asarray(img[300:900,]).min()
vmax = np.asarray(img[300:900,]).max()
mskdat=img*0.
if cmplx :
imgi= hdulist[ i+1].data
zmod = np.sqrt(img**2+imgi**2)
arg = np.arctan2(imgi,img)
aa = (arg+np.pi) /(2.*np.pi)
aa = (aa+0.5)%1
bb = 1./(abs(zmod)**0.3 +1.)
bb = -bb+1.
print np.shape(aa)
cc = np.ndarray(np.shape(aa))
print np.shape(cc)
cc.fill(1.)
iimg =[]
for ii in np.arange(np.shape(aa)[0]):
for jj in np.arange(np.shape(aa)[1]):
iimg.append( hls_to_rgb(aa[ii,jj], bb[ii,jj],cc[ii,jj] ))
img = np.reshape( iimg,(np.shape(aa)[0],np.shape(aa)[1],3))
if mask :
if cmplx :
print "mask & cmplx not yet compatible "
exit
msknam = filename +'_masks.fits'
print msknam+ " " + str(num)
hdmskl = fits.open(msknam)
mskdat = np.fabs(hdmskl[num+1].data-1.)
img = img*mskdat
filename=filename+"_masked"
if not cmplx :
vmin = np.asarray(img[bfnor[0]:bfnor[1],]).min()
vmax = np.asarray(img[bfnor[0]:bfnor[1],]).max()
else :
vmin=0
vmax=1
if mvmin !=0. :
vmin = mvmin
......@@ -118,24 +168,40 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
plt.title(filename + ' ' +mtitle,y=1.08)
plt.ylabel("Frequency (Mhz)")
plt.xlabel("Temps (h)")
if ramod:
plt.xlabel("RA (h)")
else:
plt.xlabel("Temps (h)")
if ramod :
for rap in plras:
ax1.plot([rap,rap],[frmax,frmin],color='red',linestyle='--')
#ax1.annotate(nmras[plras.index(rap)],(rap,frmin),horizontalalignment='center',verticalalignment='bottom')
ax1.annotate(nmras[plras.index(rap)],(rap,frmin),horizontalalignment='left',verticalalignment='bottom',rotation=30)
ax1.yaxis.grid(which="major", color='gray', linestyle='-', linewidth=.5)
ax1.xaxis.grid(which="major", color='gray', linestyle='-', linewidth=.5)
ax2=ax1.twinx()
ax2.yaxis.set_view_interval(len(freqs),0.,ignore=True)
ax2.set_ylabel("Freq. bin",rotation=270.,y=0.9)
ax3=ax1.twiny()
ax3.xaxis.set_view_interval(0.,len(times))
ax3.set_xlabel("time bin",x=0.85)
if not ramod :
ax3=ax1.twiny()
ax3.xaxis.set_view_interval(0.,len(times))
ax3.set_xlabel("time bin",x=0.85)
if save:
ext=""
if creal :
ext = ext+"_real"
plt.savefig(filename+"_"+ names[num]+"_rawTFM.png")
if mask :
ext=ext+"_mask"
if ramod :
if len(plras) >0 :
plt.savefig(filename+"_"+ names[num]+"_RAM"+ext+"_srcs.png")
else :
plt.savefig(filename+"_"+ names[num]+"_RAM"+ext+".png")
else :
plt.savefig(filename+"_"+ names[num]+"_TFM"+ext+".png")
axcolor = 'lightgoldenrodyellow'
axmin = fig.add_axes([0.075, 0.05, 0.35, 0.03], facecolor = axcolor )
......@@ -240,14 +306,18 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
leg = plt.legend()
for legobj in leg.legendHandles:
legobj.set_linewidth(2.0)
plt.ylim((mvmin,mvmax))
plt.xlim((0.,24.))
if mtitle =="":
plt.title(filename + ' ' + names[num],y=1.08)
else :
plt.title(filename + ' ' +mtitle,y=1.08)
if save:
plt.savefig(filename+"_"+ names[num]+"_rawTFM_freq.png")
if mask :
plt.savefig(filename+"_"+ names[num]+"_mask_TFM_ra.png")
else:
plt.savefig(filename+"_"+ names[num]+"_TFM_ra.png")
if timeplot :
......@@ -259,6 +329,7 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
leg = plt.legend()
for legobj in leg.legendHandles:
legobj.set_linewidth(2.0)
plt.ylim((mvmin,mvmax))
if mtitle =="":
......@@ -266,15 +337,20 @@ def viewtfmap(folder,num,mvmin=0.,mvmax=0.,save=False, raplot=False,timeplot=F
else :
plt.title(filename + ' ' +mtitle,y=1.08)
if save:
plt.savefig(filename+"_"+ names[num]+"_rawTFM_freq.png")
if mask :
plt.savefig(filename+"_"+ names[num]+"_TFM_mask_freq.png")
else:
plt.savefig(filename+"_"+ names[num]+"_TFM_freq.png")
if (verb) :
print times.size
if (verb) :
print np.shape(img)
if noplo :
return times/3600.,img,freqs,ras,[ img[b,] for b in bfreq]
return times/3600.,img,freqs,ras,[ img[b,] for b in bfreq],mskdat
elif cmplx:
return times/3600.,img,freqs,ras
else :
return times/3600.,img,freqs,ras,[ img[b,] for b in bfreq],ax1,bnext,bprev,smin,smax
return times/3600.,img,freqs,ras,[ img[b,] for b in bfreq],ax1,bnext,bprev,smin,smax,mskdat
......
Satellites/Makefile
View file @ d193defb
......@@ -8,16 +8,26 @@ 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
rm -f $(OBJ)/predictsatsgp4.o
rm -f $(EXE)/predictsatsgp4 $(OBJ)/predictsatsgp4.o
rm -f $(EXE)/trk2dt $(OBJ)/trk2dt.o
### 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
$(OBJ)/trk_ac_fit.o : trk_ac_fit.cc
$(CXXCOMPILE) -o $(OBJ)/trk_ac_fit.o trk_ac_fit.cc
###### Compilation et link des executables
predictsatsgp4 : $(EXE)/predictsatsgp4
......@@ -26,3 +36,11 @@ 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 (programme track file to PPF/fits )
trk2dt : $(EXE)/trk2dt
echo 'trk2dt --- made'
$(EXE)/trk2dt : $(OBJ)/trk2dt.o
$(CXXLINK) -o $(EXE)/trk2dt $(OBJ)/trk2dt.o $(SOPHYAEXTSLBLIST)
Satellites/README
View file @ d193defb
......@@ -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/
......
Satellites/predictsatsgp4.cc
View file @ d193defb
......@@ -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
......@@ -247,6 +277,9 @@ int main(int narg, char *arg[])
cout <<" velocity: |V|="<<velsat<<" km/s, angular |Vangle|="<<vangmod<<", Vaz="<<vangaz<<", Valt="<<vangalt<<" deg/min"<<endl;
cout <<" closest-los: sep="<<sepmin<<" "<<datemin<<" "<<topomin<<" "<<vsatfilename[isat]<<endl;
if(sepmin>sep) cout<<"WARNING: sep="<<sep<<"< sepmin="<<sepmin<<endl;
if (sepmin<2.) cout << "GOOD SAT " <<satname << " " << vsatfilename[isat]<<" " << " closest-los: sep="<<sepmin<<" "<<datemin<<" "<<topomin << endl;
//
} catch (SatelliteException& e) {
if(lp>0) cerr<<"Sat="<<satname<<" skipped! "<<e.what()<<" "<<vsatfilename[isat]<<endl;
......@@ -357,3 +390,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%02d%02d.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;
}
Satellites/trk2dt.cc 0 → 100644
View file @ d193defb
#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[