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

Satellites/Makefile

+# Makefile for satellite visibility search using SGP4 TLE librairy
+include $(SOPHYABASE)/include/sophyamake.inc
+
+OBJ = ./Objs/
+EXE = ./Objs/
+
+#Include et Librairies pour SGP4
+SGP4INC = ${EXTLIBDIR}/Satellites/sgp4/libsgp4/
+SGP4LIB = ${EXTLIBDIR}/Satellites/sgp4/libsgp4/
+
+#  Define our target list 
+all : predictsatsgp4
+
+clean :
+	rm -f $(EXE)/predictsatsgp4
+	rm -f $(OBJ)/predictsatsgp4.o
+
+### Compilation de .o 
+$(OBJ)/predictsatsgp4.o : predictsatsgp4.cc $(SGP4INC)/Tle.h
+	$(CXXCOMPILE) -o $(OBJ)/predictsatsgp4.o -I$(SGP4INC) predictsatsgp4.cc
+
+######  Compilation et link des executables 
+predictsatsgp4 : $(EXE)/predictsatsgp4
+	echo 'predictsatsgp4 --- made'
+
+$(EXE)/predictsatsgp4 : $(OBJ)/predictsatsgp4.o
+	$(CXXLINK) -o $(EXE)/predictsatsgp4 $(OBJ)/predictsatsgp4.o -L$(SGP4LIB) -lsgp4 -lstdc++ -lc -lm
+

Satellites/README

+--------------------------------------------------------------------------
+------ Predict satellite visibility using SGP4 library based on TLE ------
+--------------------------------------------------------------------------
+
+-----
+1-/ Step 1 : SGP4 library
+-----
+.Download and install the SPG4 library (https://www.danrw.com/sgp4/)
+ > set HERESGP4 = "...where I want to install SGP4..."
+ > cd $HERESGP4
+ > git clone https://github.com/dnwrnr/sgp4.git
+ 
+ .compil
+ > cd $HERESGP4/sgp4/
+ > cmake $HERESGP4/sgp4/
+ > make
+  -> library:  $HERESGP4/sgp4/libsgp4/libsgp4.a
+  -> includes: $HERESGP4/sgp4/libsgp4/*.h
+
+ .tests
+ > cd $EXTLIBDIR/Satellites/sgp4/runtest/
+ > ln -s ../SGP4-VER.TLE .
+ > ./runtest
+
+ > cd $EXTLIBDIR/Satellites/sgp4/passpredict/
+ > ./passpredict
+
+ > cd $EXTLIBDIR/Satellites/sgp4/sattrack/
+ > ./sattrack
+
+ .Documentation
+ https://www.danrw.com/sgp4/
+ telecharger:
+   SpaceTrack Report No. 3  -> spacetrk.pdf
+   Revisiting Spacetrack Report #3 -> AIAA-2006-6753-Rev2.pdf
+
+
+-----
+2-/ Step 2 : telecharger les TLE (two-line element sets)
+-----
+
+Utiliser le script bash :
+> getTLE.sh
+Le script va chercher les TLE pour un certain nomre de satellites.
+Il cree un repertoire "YYYMMDD" dans lequel sont telecharges les fichers de TLE.
+Le telechargement doit etre fait regulierement car les TLE changent au cours du temps pour un satellite.
+Il est conseille de garder les repertoires (aux divers dates).
+Les satellites dont on telecharge les TLE sont listes dans les variables: navigation communic rescue.
+  On peut y ajouter/enlever des satellites.
+
+-----
+3-/ Step 3 : compiler le programme de prediction
+-----
+> cd .../AnaPAON4/Satellites/
+> make SGP4INC=${HERESGP4}/sgp4/libsgp4/ SGP4LIB=${HERESGP4}/sgp4/libsgp4/
+.par defaut les objets et executables sont mis dans "./Objs/" (faire mkdir Objs)
+.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"
+
+-----
+4-/ Step 4 : executer le programme de prediction
+-----
+.help
+> ./Objs/predictsatsgp4 -h
+
+.exemples
+> ./Objs/predictsatsgp4 YYYMMDD/*.tle
+> ./Objs/predictsatsgp4 YYYMMDD/gps-ops.txt YYYMMDD/galileo.txt
+  debug for satellite GSAT0102 (PRN E12) i.e. "37847"
+> ./Objs/predictsatsgp4 -D 37847 satname YYYMMDD/galileo.txt [...]
+
+
+-----
+5-/ Remarques
+-----
+.Est aussi fourni un petit job de test (testpredsatsgp4.py) python (2.7) utilisant le module "sgp4"
+    https://pypi.python.org/pypi/sgp4/
+    OU
+    anaconda search -t conda sgp4
+ (partiellement teste, TBC)
+
+.Gpredict est un logiciel graphique s'installant tres facilement sous Linux (testé Ubuntu >=16.04)
+ code source:
+    http://gpredict.oz9aec.net/
+    https://github.com/csete/gpredict/releases
+ ou installation directe par paquet dans les depots
+    sudo apt-get install gpredict
+

Satellites/getTLE.sh

+#!/bin/bash
+# $Id: getTLE.sh
+# $auth: Steve Torchinsky <steve.torchinsky@obspm.fr>
+# $created: Mon 21 Nov 2016 11:37:49 CET
+# Modified by cmv (January 2018)
+#--------------------------------------------------------------------------------
+# get Two Line Element data (TLE), and saved it in a file with the date in the name
+#      http://www.celestrak.com/NORAD/elements/
+#--------------------------------------------------------------------------------
+#>>> L Band part 1.2 - 1.8 GHz
+#This frequency range includes a very diverse range of satellites and encompasses many sub-allocations.
+#This range includes the GPS and other GNSS (Global Navigation Satellite Systems - Russian Glonass, European Galileo, Chinese Beidou).
+#It also hosts SARSAT/COSPAS search and rescue satellites which are carried on board US and Russian meteorological satellites.
+#It also includes a mobile satellite communication band.
+#--------------------------------------------------------------------------------
+
+http='http://www.celestrak.com/NORAD/elements'
+
+#Navigation Satellites
+navigation="galileo gps-ops glo-ops beidou nnss sbas"
+#Communication Satellites
+communic="intelsat ses iridium iridium-NEXT orbcomm globalstar amateur x-comm other-comm gorizont raduga molniya stations"
+#Rescue Satellites
+rescue="sarsat "
+
+#
+todo=( ${navigation[*]} ${communic[*]} ${rescue[*]} )
+echo 'number of satellites requested:' ${#todo[*]}
+REP_NAME=`date +"%Y%m%d"`
+mkdir ${REP_NAME}
+cd ${REP_NAME}
+pwd
+for F in ${todo[*]} ; do
+    wget ${http}/${F}.txt
+    #FILE_DATE=`date -r ${F}.txt +"%Y%m%d-%H%M"`
+    #FILE_DATE=`date -r ${F}.txt +"%Y%m%d"`
+    #mv -fv ${F}.txt ${F}_${FILE_DATE}.txt
+done
+
+echo 'number of satellites requested:' ${#todo[*]}
+echo 'number of satellites found:' `ls -1 *.txt | wc -l`
+
+exit 0

Satellites/predictsatsgp4.cc

+#include <stdlib.h>
+#include <unistd.h>
+
+#include <Tle.h>
+#include <SGP4.h>
+#include <Observer.h>
+#include <CoordGeodetic.h>
+#include <CoordTopocentric.h>
+#include <DateTime.h>
+#include <SolarPosition.h>
+
+#include <list>
+#include <string>
+#include <iomanip>
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <cstdlib>
+
+using namespace std;
+using namespace Util;
+void ReadFile(string infile,vector<string>& vline1,vector<string>& vline2,vector<string>& vsatfilename);
+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 usage(void);
+void usage(void)
+{
+cout<<"predictsatsgp4 [options] filetle1.txt [filetle2.txt ...] : find satellites closed to a az-alt pointing"<<endl
+    <<" -T \"yyyy/mm/dd hh:mm:ss\" : UTC date and time of observation (def=now)"<<endl
+    <<" -O latitude(dd.ddd,N+),longitude(dd.ddd,E+),height(km) : observer location (def=paon4)"<<endl
+    <<" -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
+    <<" -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
+    <<" -v dang4vel : angular increment (deg) to compute the angular velocity of the satellite (def=1.)"<<endl
+    <<" -p lp : print debug"<<endl
+    <<endl;
+}
+
+//
+int main(int narg, char *arg[])
+{
+  //--- Default values for arguments
+  //Observateur: longitude (deg, E+), latitude (deg, N+), altitude (km)
+  double obslong = +2.1995750, obslat = +47.3820972, obsaltkm =0.155;  //Paon4 antenne centrale
+  //Heure de l'observation UTC
+  string datestr;
+  DateTime dateobs = DateTime::Now();
+  //Pointage az,alt telescope (los)
+  // azimuth N=0, E->W , altitude Z+
+  double pointaz=180., pointalt=83.40; //degres (CygA a Paon4 au transit sud du 5/10/2017 ~18:52:45)
+  //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;
+  //Time increment (min) to compute the angular velocity of the satellite
+  double dang4vel = 1.;
+  //Satellite name to debug
+  string debugsatname;
+  //print level
+  int lp=0;
+  
+  //--- Decodage des arguments
+  char c;
+  while((c = getopt(narg,arg,"hT:O:H:A:s:D:v:p:")) != -1) {
+    switch (c) {
+    case 'T' :
+      datestr = optarg;
+      break;
+    case 'O' :
+      sscanf(optarg,"%lf,%lf,%lf",&obslong,&obslat,&obsaltkm);
+      break;
+    case 'H' :
+      sscanf(optarg,"%lf,%lf",&pointaz,&pointalt);
+      break;
+    case 'A' :
+      scanlos = atof(optarg);
+      break;
+    case 's' :
+      sscanf(optarg,"%u,%u,%u",&dtspanhour,&dtspanmin,&satincmin);
+      break;
+    case 'D' :
+      debugsatname = optarg;
+      break;
+    case 'v' :
+      dang4vel = atof(optarg);
+      break;
+    case 'p' :
+      lp = atoi(optarg);
+      break;
+    case 'h' :
+    default :
+      usage(); return -1;
+    }
+  }
+  if(optind>=narg) {
+    cout<<"Give TLE files: paonsat4 [options] filetle1.txt [filetle2.txt ...]"<<endl;
+    usage();
+    return -1;
+  }
+ 
+  //--- Observateur
+  Observer obspaon4(obslat,obslong,obsaltkm);
+  cout<<"Observer: "<<obspaon4.GetLocation()<<endl;
+  
+  //--- Observation date and time (UTC)
+  if(datestr.size() > 2) {
+    int year=0, month=0, day=0, hour=0, minute=0, second=0;
+    sscanf(datestr.c_str(),"%d/%d/%d %d:%d:%d",&year,&month,&day,&hour,&minute,&second);
+    dateobs.Initialise(year,month,day,hour,minute,second,0.);
+  }
+  cout<<"DateObs: "<<dateobs<<", jd="<<setprecision(16)<<dateobs.ToJulian()<<setprecision(8)<<endl;
+  cout<<"Sideral time: Greenwich="<<dateobs.ToGreenwichSiderealTime()/M_PI*12.<<" h, "
+      <<"local="<<dateobs.ToLocalMeanSiderealTime(obspaon4.GetLocation().longitude/*rd*/)/M_PI*12.<<" h"<<endl;
+   
+  //--- Pointage az,alt telescope
+  pointaz = Wrap360(pointaz);
+  if(pointalt<-90. || pointalt>90.) {cout<<"Bad altitude="<<pointalt<<endl; return -3;}
+  Vector Vpoint = AzAlt2Vec(pointaz,pointalt);
+  cout<<"Pointage: az="<<pointaz<<", alt="<<pointalt<<" deg : "<<Vpoint<<endl;
+
+  //--- Scan autour de +/- l'heure d'observation
+  TimeSpan tspan(dtspanhour,dtspanmin,0);
+  if(satincmin<=0) satincmin = 1;
+  TimeSpan tspaninc(0,satincmin,0); 
+  cout<<"TSpan: "<<tspan<<" / inc="<<tspaninc<<endl;
+  DateTime datestart = dateobs - tspan;
+  DateTime dateend = dateobs + tspan;
+  cout<<"DateScan: "<<datestart<<" -> "<<dateend<<endl;
+  
+  //--- Angular increment to compute the angular velocity of the satellite
+  if(dang4vel<=0.) dang4vel = 1.;  //deg
+  cout<<"Angular increment for velocity computation: "<<dang4vel<<" deg"<<endl;
+
+  //--- Debug satellite:
+  if(debugsatname.size()>0) cout<<"Debug satellite: "<<debugsatname<<endl;
+    
+  //--- Sun (ephemerides de la journee d'observation)
+  cout<<endl;
+  DateTime dt0(dateobs.Year(),dateobs.Month(),dateobs.Day(),0,0,0);
+  for(int i=0;i<=24;i+=1) {
+    SolarPosition sunpos;
+    DateTime dt = dt0 + TimeSpan(i,0,0);  //inc=1h
+    Eci ecisun = sunpos.FindPosition(dt);
+    CoordTopocentric toposun = obspaon4.GetLookAngle(ecisun);
+    CoordGeodetic geosun = ecisun.ToGeodetic();
+    if(RadiansToDegrees(toposun.elevation) < -18.0) continue;
+    cout <<"Sun: "<<dt<<" "<<toposun<<" "<<geosun<<endl;
+  }
+
+  //--- Read TLE parameters for all satellite files and fill TLE into vector<>
+  cout<<endl;
+  vector<string> vline1, vline2, vsatfilename;
+  int nsat= 0;
+  for(int i=optind;i<narg;i++) {
+    ReadFile(string(arg[i]),vline1,vline2,vsatfilename);
+    if(lp>0) cout<<"FILE: "<<arg[i]<<" : nsat="<<vline1.size()-nsat<<endl;
+    nsat = vline1.size();
+  }
+  nsat = vline1.size();
+  cout<<"Number of Satellites found: "<<nsat<<endl;
+
+  //--- Loop on satellites
+  if(nsat<1) return -2;
+  for(int isat=0;isat<nsat;isat++) {
+    int idum; char strsatname[64];
+    sscanf(vline2[isat].c_str(),"%d %s",&idum,strsatname);
+    string satname(strsatname);
+    Tle tle = Tle(satname,vline1[isat],vline2[isat]);
+    long ifounddbg = (debugsatname.size()>0) ? satname.find(debugsatname,0): -1;
+    //
+    try {
+      SGP4 sgp4(tle);
+      //...Calculate satellite position and its velocity
+      Eci eci = sgp4.FindPosition(dateobs);
+      Vector Velsat = eci.Velocity();
+      double velsat = Velsat.Magnitude(); //km/s
+      //...Get look angle for observer to satellite
+      CoordTopocentric topo = obspaon4.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));
+      //
+      //...Separation / telescope pointing
+      double sep = SepAngle(Vpoint,Vsat);
+      if(sep>scanlos && ifounddbg!=0) continue;
+      //
+      //...Compute satellite angular speed
+      double vangmod = -1., vangaz=0., vangalt=0.;
+      {
+      //satellite velocity and altitude -> angular speed -> choose the right time increment
+      //  approx Vlos << Vtranserse + neglect observer velocity due to earth rotation
+      double vangsat = velsat / geo.altitude;  //angular velocity rd/s
+      double tsec = 10.;  // if velocity==0, give any value
+      if(vangsat>0.) tsec = DegreesToRadians(dang4vel) / vangsat;  //sec
+      TimeSpan dt4vel(int64_t(tsec*TicksPerSecond));
+      //compute new satellite position after the small time increment
+      Eci eci1 = sgp4.FindPosition(dateobs+dt4vel);
+      //compute the satellite speed in azimuth and altitude
+      //  approximation des petites variations (ATTENTION pres des poles, indetermination!)
+      CoordTopocentric topo1 = obspaon4.GetLookAngle(eci1);
+      Vector Vsat1 = AzAlt2Vec(RadiansToDegrees(topo1.azimuth),RadiansToDegrees(topo1.elevation));
+      /*
+      cout<<"--DBG-- vangsat="<<RadiansToDegrees(vangsat)*60.<<" deg/min, tinc="<<dt4vel
+	  <<", az: "<<RadiansToDegrees(topo.azimuth)<<" -> "<<RadiansToDegrees(topo1.azimuth)
+	  <<" , alt: "<<RadiansToDegrees(topo.elevation)<<" -> "<<RadiansToDegrees(topo1.elevation)<<" deg"<<endl;
+      */
+      vangmod = SepAngle(Vsat,Vsat1) / dt4vel.TotalMinutes();  //en deg/min
+      DelAzAltSmall(topo,topo1,vangaz,vangalt);
+      vangaz /= dt4vel.TotalMinutes(); vangalt /= dt4vel.TotalMinutes();  //en deg/min
+      }
+      //...Search the closest approch to los in the specified time range
+      //   (initialize at observing date)
+      CoordTopocentric topomin = topo;
+      DateTime datemin = dateobs;
+      double sepmin = sep;
+      DateTime datecur = datestart;
+      while(datecur <= dateend+tspaninc) {
+        Eci ecisep = sgp4.FindPosition(datecur);
+        CoordTopocentric toposep = obspaon4.GetLookAngle(ecisep);
+        Vector Vsatsep = AzAlt2Vec(RadiansToDegrees(toposep.azimuth),RadiansToDegrees(toposep.elevation));
+        double s = SepAngle(Vpoint,Vsatsep);
+        ////cout <<"--DBG-- sep="<<s<<" datecur="<<datecur<<" topo="<<toposep<<endl;
+	if(s<sepmin) {sepmin = s; datemin = datecur; topomin = toposep;}
+	datecur = datecur + tspaninc;
+      }
+      //...Print results
+      //INFO: "Rng" is the "slant range" (distance observer-satellite) in km
+      //  and "Rng Rt" is the "slant range rate" in km/s
+      if(ifounddbg==0) {
+	cout <<"***DBG-SAT: "<<endl;
+	if(lp>1) {
+	  cout<<"Eci: date="<<eci.GetDateTime()<<endl
+	      <<"     pos="<<eci.Position()<<" , |pos|="<<eci.Position().Magnitude()<<endl
+	      <<"     vel="<<eci.Velocity()<<" , |vel|="<<eci.Velocity().Magnitude()<<endl;
+	  cout<<"CoordTopocentric: "<<topo<<" , Vunit: "<<Vsat<<endl;
+	  cout<<"CoordGeodetic: "<<geo<<" , Vunit: "<<Vungeo<<endl;
+	}
+      }
+      cout <<"sat="<<satname<<" sep="<<sep<<" "<<dateobs<<" "<<topo<<" "<<geo<<endl;
+      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;
+    //
+    } catch (SatelliteException& e) {
+      if(lp>0) cerr<<"Sat="<<satname<<" skipped! "<<e.what()<<" "<<vsatfilename[isat]<<endl;
+      continue;
+    //
+    } catch (DecayedException& e) {
+      if(lp>0) cerr<<"Sat="<<satname<<" skipped! "<<e.what()<<" "<<vsatfilename[isat]<<endl;
+      continue;
+    }
+    //
+  }
+  
+  return 0;
+}
+
+//-----------------------------------------------------------------------
+void ReadFile(string infile,vector<string>& vline1,vector<string>& vline2,vector<string>& vsatfilename)
+{
+  ifstream file;
+  file.open(infile.c_str());
+  if(!file.is_open()) {cerr << "Error opening file" << endl; return;}
+  //get file name without directory
+  {
+  unsigned long ip=infile.find_last_of('/') + 1;
+  if(ip<infile.size()) infile = infile.substr(ip);
+  }
+
+  bool got_first_line = false;
+  string line1, line2, parameters;
+  while (!file.eof()) {  //Begin while
+    string line;
+    getline(file,line);
+    Util::Trim(line);
+    // skip blank lines or lines starting with #
+    if(line.length() == 0 || line[0] == '#') { got_first_line = false; continue;}
+    // find first line
+    if(!got_first_line) {  //111111111111111111111111111111111111111111111111111
+      try {
+        if(line.length() >= Tle::LineLength()) {
+          //Tle::IsValidLine(line.substr(0, Tle::LineLength()), 1);
+          // store line and now read in second line
+          got_first_line = true;
+          line1 = line;
+	}
+      } catch (TleException& e) {
+	cerr << "Error: " << e.what() << endl;
+	cerr << infile <<"\n"<< line << endl;
+      }
+    } else {  //111111111111111111111111111111111111111111111111111
+      // no second chances, second line should follow the first
+      got_first_line = false;
+      // split line, first 69 is the second line of the tle
+      // the rest is the test parameters, if there is any
+      line2 = line.substr(0, Tle::LineLength());
+      try { // following line must be the second line
+	if(line.length() >= Tle::LineLength()) {
+          //Tle::IsValidLine(line.substr(0, Tle::LineLength()), 2);
+	  vline1.push_back(line1);
+	  vline2.push_back(line2);
+	  vsatfilename.push_back(infile);
+	}
+      } catch (TleException& e) {
+	cerr << "Error: " << e.what() << endl;
+	cerr << infile <<"\n"<< line << endl;
+      }
+    }  //111111111111111111111111111111111111111111111111111
+  }  //End while
+
+  // close file
+  file.close();
+  return;
+}
+
+//---------------------------------------------------------
+Vector AzAlt2Vec(double az,double alt)
+// convert azimuth and altitude (deg) into cartesian unitary vector
+//WARNING: coordinate system azimuth (N=0,E->W), altitude from horizon
+//  x'=N  |  y'=E |  z=zenith ==> left handed frame (indirect)
+// Return vector with: x=N | y=W | z=zenith ==> right handed frame (direct)
+{
+  az = DegreesToRadians(az);
+  alt = DegreesToRadians(alt);
+  double z = sin(alt);
+  double x = cos(az)*cos(alt);
+  double y = sin(az)*cos(alt); y = -y; // E->W
+  return Vector(x,y,z);
+}
+
+//---------------------------------------------------------
+double SepAngle(Vector& v1, Vector& v2)
+// return separation angle between 2 vector in degree
+{
+  double dot = v1.Dot(v2);
+  dot /= v1.Magnitude() * v2.Magnitude();
+  dot = acos(dot);
+  return RadiansToDegrees(dot);
+}
+
+//---------------------------------------------------------
+void DelAzAltSmall(const CoordTopocentric& topo1,const CoordTopocentric& topo2,double &daz,double &dalt)
+//Return daz=az2-az1 and dalt=alt2-alt1 for "SMALL" drift (all angles in degrees)
+//Warning: we take the smallest possibility, be carefull for drift closed to the poles even if small
+{
+  double az1 = RadiansToDegrees(topo1.azimuth), alt1 = RadiansToDegrees(topo1.elevation);
+  double az2 = RadiansToDegrees(topo2.azimuth), alt2 = RadiansToDegrees(topo2.elevation);
+  dalt = alt2 - alt1;  //no problem, alt in [-90,90]
+  daz = az2 - az1;  //deal with az in [0,360[
+  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
+}

Satellites/quick_compil.csh

+#!/bin/csh
+
+if( "$1" == "-c" ) then
+  rm -f predictsatsgp4 predictsatsgp4.o
+  exit 0
+endif
+
+set SGP4INC = $EXTLIBDIR/Satellites/sgp4/libsgp4/
+set SGP4LIB = $EXTLIBDIR/Satellites/sgp4/libsgp4/
+
+gcc -c -o predictsatsgp4.o -I$SGP4INC predictsatsgp4.cc
+gcc -o predictsatsgp4 predictsatsgp4.o -L$SGP4LIB -lsgp4 -lstdc++ -lc -lm
+

Satellites/testpredsatsgp4.py

+#!/usr/bin/env python
+#http://www.astropy.org/
+#aussi la version C++: https://www.danrw.com/sgp4/
+
+import sys, os, io
+import numpy as np
+from datetime import datetime
+import matplotlib.pyplot as plt
+
+from astropy.time import Time
+from astropy import units as u
+from astropy.coordinates import SkyCoord, EarthLocation, AltAz, Angle, CartesianRepresentation
+from astropy.coordinates import GCRS, PrecessedGeocentric, ICRS
+from astropy.coordinates import solar_system_ephemeris, get_body, get_moon
+from astroplan import Observer , FixedTarget
+
+from sgp4.earth_gravity import wgs72
+from sgp4.io import twoline2rv
+
+#---observations
+locpaon4 = EarthLocation.from_geodetic(lon=+2.1995750*u.degree,lat=+47.3820972*u.degree,height=155.*u.m) #Nancay-Paon4
+###locpaon4 = EarthLocation.from_geodetic(lon=2.1333*u.degree,lat=48.6833*u.degree,height=155.*u.m) #Gif-sur-Yvette
+print "EarthLocation paon4:",locpaon4
+obspaon4 = Observer(name='Paon4',location=locpaon4,pressure=1.013*u.bar,relative_humidity=0.30, \
+                   temperature=0*u.deg_C,timezone="Europe/Paris", \
+                   description="4 Parabolas Interferometer")
+print "Observer paon4:",obspaon4
+tuobs = Time(datetime.utcnow(),format="datetime",scale="utc",location=locpaon4)
+equinobs = tuobs.decimalyear
+print "TU date:",tuobs,":",equinobs
+obslst = obspaon4.local_sidereal_time(tuobs)  #ou tuobs.sidereal_time('apparent')
+print "Local sideral time:",obslst
+print "Moon: alt: {0.alt}, az: {0.az}".format(obspaon4.moon_altaz(tuobs))
+print "Sun: rise=",obspaon4.sun_rise_time(tuobs).iso,"set=",obspaon4.sun_set_time(tuobs).iso
+
+pointpaon4 = AltAz(az=180*u.deg,alt=83.40*u.deg,obstime=tuobs,location=locpaon4)
+print "Pointage telescope: az=",pointpaon4.az.degree,"d (N=0,E+), alt=",pointpaon4.alt.degree
+
+#---Sun, Moon and Jupiter
+print ""
+solar_system_ephemeris.set('builtin')
+
+suncoord = get_body('sun',time=tuobs,location=locpaon4)
+sunaltaz = suncoord.transform_to(frame=obspaon4.altaz(tuobs))
+hasun = obspaon4.target_hour_angle(time=tuobs,target=suncoord)
+print "Sun: up?",obspaon4.target_is_up(time=tuobs,target=suncoord), \
+      "az=",sunaltaz.az.degree,"d, alt=",sunaltaz.alt.degree,"d", \
+      "ha=",hasun.hour,"h"
+
+mooncoord = get_body('moon',time=tuobs,location=locpaon4)
+moonaltaz = mooncoord.transform_to(frame=obspaon4.altaz(tuobs))
+hamoon = obspaon4.target_hour_angle(time=tuobs,target=mooncoord)
+print "Moon: up?",obspaon4.target_is_up(time=tuobs,target=mooncoord), \
+      "az=",moonaltaz.az.degree,"d, alt=",moonaltaz.alt.degree,"d", \
+      "ha=",hamoon.hour,"h"
+
+jupcoord = get_body('jupiter',time=tuobs,location=locpaon4)
+jupaltaz = jupcoord.transform_to(frame=obspaon4.altaz(tuobs))
+hajup = obspaon4.target_hour_angle(time=tuobs,target=jupcoord)
+print "Jupiter up?:",obspaon4.target_is_up(time=tuobs,target=jupcoord), \
+      "az=",jupaltaz.az.degree,"d, alt=",jupaltaz.alt.degree,"d", \
+      "ha=",hajup.hour,"h"