Commit 35796736 authored by Magneville's avatar Magneville
Browse files

Prediction des passages de satellites (cmv 28/01/2018)

parent 2564a2b4
# 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
--------------------------------------------------------------------------
------ 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
#!/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
#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,N+) : (def=180d,83.40d)"<<endl
<<" telescope horizontal coordinates pointing"<<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 (-S) finely scan at +/-(hh:mm) with step \"spaninc\" (def=1h,30mn,1mn)"<<endl
<<" -D satname : print debug for this satellite (even if not in the search area)"<<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 N+
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
CoordTopocentric topomin = topo;
DateTime datemin = datestart, datecur = datestart;
double sepmin = sep;
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
}
#!/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
#!/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"
#---CygA et CasA
print ""
cygcoord = SkyCoord.from_name("cyga",frame=u'gcrs')
cygaltaz = cygcoord.transform_to(frame=obspaon4.altaz(tuobs))
hacyg = obspaon4.target_hour_angle(time=tuobs,target=cygcoord)
print "CygA up?:",obspaon4.target_is_up(time=tuobs,target=cygcoord), \
"az=",cygaltaz.az.degree,"d, alt=",cygaltaz.alt.degree,"d", \