# -*- coding: utf-8 -*-
"""
@author: aklotz@irap.omp.eu
"""

import os
import time
import skyfield.api
from skyfield.api import wgs84
from astroquery.simbad import Simbad
from astroquery.mpc import MPC
import astropy.units as u
from astropy.coordinates import ICRS, AltAz, HADec, EarthLocation, SkyCoord
from astropy.time import Time
import numpy as np
import requests

try:
    from .home import Home
except:
    from home import Home

try:
    from .dates import Date
except:
    from dates import Date

try:
    from .durations import Duration
except:
    from durations import Duration

try:
    from .angles import Angle
except:
    from angles import Angle

try:
    from .coords import Coords
except:
    from coords import Coords

try:
    from .filenames import FileNames
except:
    from filenames import FileNames

try:
    from .guitastrotools import GuitastroTools, GuitastroException
except:
    from guitastrotools import GuitastroTools, GuitastroException

# #####################################################################
# #####################################################################
# #####################################################################
# Class Ephemeris
# #####################################################################
# #####################################################################
#
#
# #####################################################################

class EphemerisException(GuitastroException):
    """Exception raised for errors in the Ephemeris class.
    """

    TARGET_NOT_FOUND = 0
    DATE_OUTSIDE_THE_NIGHT = 1

    errors = [""]*2
    errors[TARGET_NOT_FOUND] = "Target not found"
    errors[DATE_OUTSIDE_THE_NIGHT] = "The date is outside the night limits"



class Ephemeris(EphemerisException, GuitastroTools):
    """Compute coordinates of an object from various input formats

    :Usage:

    First, instanciate an object from the class:

    ::

        eph = Ephemeris()

    Second, compute coordinates of the Sun at Guitalens at the current time:

    ::

        target  = "sun"
        eph.set_home("guitalens")
        date = "Now"
        ephem = eph.radec(target, unit_ra="H0.2", unit_dec="d+090.1")
        ra, dec, equinox, epoch = ephem['ra_equinox'], ephem['dec_equinox'], ephem['header']['equinox'], ephem['jd']
        print(f"{name} ra={ra} dec={dec}")

    To get the drift:

        target  = "sun"
        eph.set_home("guitalens")
        date = "Now"
        ephem = eph.radec_speed(target, unit_ra="H0.2", unit_dec="d+090.1")
        ra, dec, equinox, epoch, dra, ddec = ephem['ra_equinox'], ephem['dec_equinox'], ephem['header']['equinox'], ephem['jd'], ephem['dra_equinox'], ephem['ddec_equinox']
        print(f"{name} ra={ra} dec={dec} dra={dra} ddec={ddec}")

    Note the output of the method radec is ephemeris dictionary.

    Other examples of inputs:

    ::

        eph.set_home("GPS 2.567 E -21.456 1205")
        date = "2022-01-31T22:34:15"
        ephem = eph.radec(target, unit_ra="H0.2", unit_dec="d+090.1")
        ra, dec, equinox, epoch = ephem['ra_equinox'], ephem['dec_equinox'], ephem['header']['equinox'], ephem['jd']

    The target is a major planet of the solar system (+ Moon and Sun):

    ::

        target = "neptune"

    The target is a name which can be solved by Simbad:

    ::

        target = "UGC 3462"

    The target is a name which can be solved by MPCEphem:

    ::

        target = "1994 PC1345"

    The target is a RA, Dec coordinates:

    ::

        target = "RADEC 13h29m40s +47d04m09s"

    The target is a RA, Dec coordinates followed by the drift dRA, dDec (drift in deg/s):

    ::

        target = "RADECDRIFT 13h29m40s +47d04m09s 0.0123 -0.0452"

    The target is a list of coordinates at given times.
    The output is an interpolation:

    ::

        target  = "DATERADECS "
        target += "2022-01-31T22:34:00 : 12 34 32.12 -01 45 46.2\\n"
        target += "2022-01-31T22:35:00 : 12 34 33.23 -01 56 04.5\\n"

    The target is designed in GCN Circulars:

    ::

        target = "GCNC 220408A"

    The target is a name which can be solved by SGP4 and the Celestrack database.
    This example shows how to compute the derivatives of the coordinates:

    ::

        target = "CELESTRACK ISS"
        ephem = eph.radec_speed(target, date="now")
        ra, dec, equinox, epoch, dra, ddec = ephem['ra_equinox'], ephem['dec_equinox'], ephem['header']['equinox'], ephem['jd'], ephem['dra_equinox'], ephem['ddec_equinox']
        print(f"{name} ra={ra} dec={dec} dra={dra*3600:.3f} ddec={ddec*3600:.3f}")


    The target is a name which can be solved by SGP4 and its TLE.

    ::

        name = "GSAT0203 (PRN E26)"
        tle  = f"TLE {name}\\n"
        tle += "1 40544U 15017A   22029.59832018 -.00000064  00000+0  00000+0 0  9991\\n"
        tle += "2 40544  56.7391  25.1641 0005059 269.6475  90.3121  1.70475734 42591\\n"
        ephem = eph.radec_speed(tle)
        ra, dec, equinox, epoch, dra, ddec = ephem['ra_equinox'], ephem['dec_equinox'], ephem['header']['equinox'], ephem['jd'], ephem['dra_equinox'], ephem['ddec_equinox']
        print(f"{name} ra={ra:.6f} dec={dec:.6f} dra={dra*3600:.5f} ddec={ddec*3600:.5f}")

    """

    TARGET_TYPE_NAME = 1
    TARGET_TYPE_RADEC = 2
    TARGET_TYPE_DATERADECS = 3
    TARGET_TYPE_TLEFILES = 4
    TARGET_TYPE_TLE = 5
    TARGET_TYPE_GCNC = 6
    TARGET_TYPE_MPC = 7
    TARGET_TYPE_RADECDRIFT = 8
    TARGET_TYPE_HADEC = 9
    TARGET_TYPE_AZELEV = 9

    def __init__(self):
        self._observatory = None
        self.home = None
        self._obscodes = None
        self._earthsatellites = None
        # ---
        self._planets = skyfield.api.load('de421.bsp')
        self._earth = self._planets['earth']
        # ---
        self.home = Home("GPS 0.1423 E 42.936639 2880")
        self.set_home(self.home)
        # ---
        self._ts = skyfield.api.load.timescale()
        self._earthsatellites = None
        # ---
        self._computed_ephem = {}

    def set_home(self, home):
        """Set the home position (on Earth)

        home can be "148", "guitalens"
        """
        found = False
        if type(home)==Home:
            self.home = home
            found = True
        else:
            try:
                self.home = Home(home)
                found = True
            except:
                home = str(home).upper()
                #print(f"home={home}")
                if isinstance(self._obscodes, type(None)) == True:
                    self._obscodes = MPC.get_observatory_codes()
                # astropy.table.table.Table
                # self._obscodes[0].colnames
                # ['Code', 'Longitude', 'cos', 'sin', 'Name']
                # self._obscodes[lig][col]
                for klig in range(len(self._obscodes)):
                    lig = self._obscodes[klig]
                    code, longitude, cos, sin, name = lig
                    if home==code:
                        found = True
                        break
                    if name.upper().find(home) >= 0:
                        found = True
                        #print(f"code={code}")
                        break
                if found==True:
                    home = f"MPC {longitude} {cos} {sin}"
                    self.home = Home(home)
        if found==False:
            if self.home==None:
                self.home = Home("GPS 0.1423 E 42.936639 2880")
        latitude = self.home.latitude
        longitude = self.home.longitude
        self._observatory = self._earth + wgs84.latlon(latitude, longitude)

    def _date2ts(self, date):
        date = Date(date)
        y, m, d, hh, mm, ss = date.ymdhms()
        self._t = self._ts.utc(y, m, d, hh, mm, ss)

    def name2solar_system_planet(self, name=""):
        """ Return the identificator of skyfield planet from a human name
        """
        lps = []
        lplanets = list(self._planets.names().values())
        for lplanet in lplanets:
            for planet in lplanet:
                lps.append(planet)
        res = None
        if name=="":
            return lps
        else:
            name = name.upper()
            if name in lps:
                return name
            name_bar = name + "_BARYCENTER"
            if name_bar in lps:
                return name_bar
            if len(name)>=3:
                name = name[:3]
                for lp in lps:
                    if name == lp[:3]:
                        res = lp
            return res

    def tle_files(self):
        # ---
        tlefiles = []
        # --- Weather & Earth Resources Satellites
        tlefiles.append("weather.txt")
        tlefiles.append("noaa.txt")
        tlefiles.append("goes.txt")
        tlefiles.append("resource.txt")
        tlefiles.append("sarsat.txt")
        tlefiles.append("dmc.txt")
        tlefiles.append("tdrss.txt")
        tlefiles.append("argos.txt")
        tlefiles.append("planet.txt")
        tlefiles.append("spire.txt")
        # --- Communications Satellites
        tlefiles.append("geo.txt")
        tlefiles.append("gpz.txt")
        tlefiles.append("gpz-plus.txt")
        tlefiles.append("intelsat.txt")
        tlefiles.append("ses.txt")
        tlefiles.append("iridium.txt")
        tlefiles.append("iridium-NEXT.txt")
        tlefiles.append("starlink.txt")
        tlefiles.append("orbcomm.txt")
        tlefiles.append("globalstar.txt")
        tlefiles.append("amateur.txt")
        tlefiles.append("x-comm.txt")
        tlefiles.append("other-comm.txt")
        tlefiles.append("satnogs.txt")
        tlefiles.append("gorizont.txt")
        tlefiles.append("raduga.txt")
        tlefiles.append("molniya.txt")
        # --- Navigation Satellites
        tlefiles.append("gps-ops.txt")
        tlefiles.append("glo-ops.txt")
        tlefiles.append("galileo.txt")
        tlefiles.append("beidou.txt")
        tlefiles.append("sbas.txt")
        tlefiles.append("nnss.txt")
        tlefiles.append("musson.txt")
        # --- Scientific Satellites
        tlefiles.append("science.txt")
        tlefiles.append("geodetic.txt")
        tlefiles.append("engineering.txt")
        tlefiles.append("education.txt")
        # --- Miscellaneous Satellites
        tlefiles.append("military.txt")
        tlefiles.append("radar.txt")
        tlefiles.append("cubesat.txt")
        tlefiles.append("other.txt")
        return tlefiles

    def tle_delete(self):
        tlefiles = self.tle_files()
        for tlefile in tlefiles:
            if os.path.exists(tlefile)==True:
                os.remove(tlefile)

    def tle_download(self, reload=False):
        """Download the Celestrack TLE files if needed.

        The TLE files are considered obsolete after 3 days.
        """
        dtmax = 86400*3 # s
        if self._earthsatellites != None:
            dt = time.time() - self._earthsatellite_time
            if dt < dtmax:
                return
        # ---
        tlefiles = self.tle_files()
        stations_url0 = 'http://celestrak.com/NORAD/elements/'
        kf = 0
        satellites = []
        for tlefile in tlefiles:
            reload = False
            if os.path.exists(tlefile)==True:
                dt = time.time() - os.path.getmtime(tlefile)
                # print(f"tlefile={tlefile} dt={dt}")
                if dt > dtmax:
                    reload = True
            # lire la date du fichier et forcer ou non la mise à jour
            # 'http://celestrak.com/NORAD/elements/stations.txt'
            stations_url = stations_url0 + tlefile
            try:
                sats = skyfield.api.load.tle_file(stations_url, reload = reload)
                if kf==0:
                    satellites = sats
                else:
                    satellites.extend(sats)
            except:
                pass
            kf += 1
        self._earthsatellites = satellites
        self._earthsatellite_time = time.time()
        return len(satellites)

    def gcnc_download(self, name):
        """Download the gcnc files if needed.
        From https://gcn.gsfc.nasa.gov/selected.html
        Valid only from 020305 to 230414B
        TODO https://github.com/nasa-gcn/gcn-kafka-python for earliests.
        or https://gcn.nasa.gov/docs/contributing
        """
        equinox = "J2000"
        sra = "0"
        sdec = "0"
        gcncfile = name.upper()+".gcn3"
        yy = gcncfile[0:2]
        mm = gcncfile[2:4]
        dd = gcncfile[4:6]
        t0 = f"20{yy}-{mm}-{dd}T00:00:00"
        gcnc_url0 = 'https://gcn.gsfc.nasa.gov/other/'
        ident = gcncfile.replace(" ","+")
        url = f"{gcnc_url0}{ident}"
        res = requests.get(url, timeout = 10)
        if res.ok==False:
            # GCNC not found
            return sra, sdec, equinox, t0
        texte = res.text
        lignes = texte.split("\n")
        st0 = ""
        # --- Swift
        # At 11:54:30 UT, the Swift Burst Alert Telescope (BAT) triggered and
        if sra=="0" or sdec=="0":
            for ligne in lignes:
                key = "the Swift Burst Alert Telescope (BAT) triggered"
                k1 = ligne.find(key)
                if k1 >= 0:
                    st0 = ligne.split()[1]
                    t0 = t0[0:10]+"T"+st0
                ligne = ligne.upper()
                keys = ["RA(J2000):","RA (J2000):","RA(J2000)  =", "RA(J2000) ="]
                for key in keys:
                    k1 = ligne.find(key)
                    if k1 >= 0:
                        k1 += len(key)
                        sra = ligne[k1:]
                keys = ["DEC(J2000):","DEC (J2000):","DEC(J2000) ="]
                for key in keys:
                    k1 = ligne.find(key)
                    if k1 >= 0:
                        k1 += len(key)
                        sdec = ligne[k1:]
        # --- Fermi
        # At 22:23:51 UT on 10 Mar 2022, the Fermi Gamma-ray Burst Monitor (GBM) triggered and located
        if sra=="0" or sdec=="0":
            for ligne in lignes:
                key = "the Fermi Gamma-ray Burst Monitor (GBM) triggered"
                k1 = ligne.find(key)
                if k1 >= 0:
                    st0 = ligne.split()[1]
                    t0 = t0[0:10]+"T"+st0
                ligne = ligne.upper()
                ligne = ligne.replace(",","")
                key = "IS RA = "
                k1 = ligne.find(key)
                if k1 >= 0:
                    k1 += len(key)
                    lig = ligne[k1:]
                    lig = lig.replace("DEC =","")
                    ligs = lig.split()
                    sra = ligs[0]
                    sdec = ligs[1]
        elif st0=="":
            for ligne in lignes:
                key = "the Fermi Gamma-ray Burst Monitor (GBM) triggered"
                k1 = ligne.find(key)
                if k1 >= 0:
                    st0 = ligne.split()[1]
                    t0 = t0[0:10]+"T"+st0
                    #print(f"ligne={ligne}")
        return sra, sdec, equinox, t0

    def radec(self, target: str, **kwargs)-> tuple:
        """Return ra, dec, equinox, epoch of a target.

        Args:

            target: A string that starts with a keyword:

                * 'RADEC': Followed by an equatorial Right Ascension, Declination position.
                * 'RADECDRIFT': Followed by an equatorial Right Ascension, Declination position and the drift.
                * 'CELESTRACK' or 'TLEFILES': Followed by a satellite name in the Celestrack TLE files.
                * 'DATERADECS': Followed by a list of equatorial Right Ascension, Declination positions.
                * 'HADEC': Followed by a list of equatorial true Hour Angle, Declination positions.
                * 'AZELEV': Followed by a list of equatorial true Azimut, Elevation positions.
                * 'TLE': Followed by a satellite defined by its TLE (Two Line Elements).
                * 'GCNC': Followed by a number which is a GRB name (e.g. GCNC 990123) to search information in GCN circulars.
                * 'MPC': Followed by a name which is a Solar System Body.

            **kwargs: A dictionary of options, keys can be:

                * 'date': To compute ephemeris for a given date ("now" for now).
                * 'unit_ra': To indicate the output format (see Angle)
                * 'unit_dec': To indicate the output format (see Angle)
                * 'target_type': To indicate the input format if the keyword is not indicated in the target string.
                * 'target_type_only': To return the identified input format.

        Returns:

            ra: Target Right Ascention position at the given date for a given equinox
            dec: Target Declination at the given date for a given equinox
            equinox: Equinox of the position
            epoch: Date of the position when the object is moving

            or

            target_type: The identified input format.

        """
        equinox = "J2000"
        date = "now"
        unit_ra = "deg" # "H0.2"
        unit_dec = "deg" # "d+090.1"
        target_type = self.TARGET_TYPE_NAME
        target_type_only = False
        if len(kwargs) > 0:
            keys = kwargs.keys()
            if "date" in keys:
                date = kwargs["date"]
            if "unit_ra" in keys:
                unit_ra = kwargs["unit_ra"]
            if "unit_dec" in keys:
                unit_dec = kwargs["unit_dec"]
            if "target_type_only" in keys:
                target_type_only = kwargs["target_type_only"]
            if "target_type" in keys:
                target_t = kwargs["target_type"].upper()
                if target_t=="TLEFILES" or target_t=="CELESTRACK":
                    target_type = self.TARGET_TYPE_TLEFILES
                elif target_t=="TLE":
                    target_type = self.TARGET_TYPE_TLE
                elif target_t=="DATERADECS":
                    target_type = self.TARGET_TYPE_EPHEMRADEC
                elif target_t=="RADEC":
                    target_type = self.TARGET_TYPE_RADEC
                elif target_t=="RADECDRIFT":
                    target_type = self.TARGET_TYPE_RADECDRIFT
                elif target_t=="HADEC":
                    target_type = self.TARGET_TYPE_HADEC
                elif target_t=="AZELEV" or target_t=="ALTAZ":
                    target_type = self.TARGET_TYPE_AZELEV
                elif target_t=="GCNC":
                    target_type = self.TARGET_TYPE_GCNC
                elif target_t=="MPC":
                    target_type = self.TARGET_TYPE_MPC
                else:
                    target_type = self.TARGET_TYPE_NAME
        epoch = Date(date).iso(3)
        # ---
        target_init = target
        target = str(target)
        res = target.split()
        if len(res) > 0:
            target_t = res[0].upper()
            if target_t=="TLEFILES" or target_t=="CELESTRACK":
                target_type = self.TARGET_TYPE_TLEFILES
                target = target[len(res[0])+1:]
            elif target_t=="DATERADECS":
                target_type = self.TARGET_TYPE_DATERADECS
                target = target[len(res[0])+1:]
            elif target_t=="RADEC":
                target_type = self.TARGET_TYPE_RADEC
                target = target[len(res[0])+1:]
            elif target_t=="RADECDRIFT":
                target_type = self.TARGET_TYPE_RADECDRIFT
                target = target[len(res[0])+1:]
            elif target_t=="HADEC":
                target_type = self.TARGET_TYPE_HADEC
                target = target[len(res[0])+1:]
            elif target_t=="AZELEV" or target_t=="ALTAZ":
                target_type = self.TARGET_TYPE_AZELEV
                target = target[len(res[0])+1:]
            elif target_t=="GCNC":
                target_type = self.TARGET_TYPE_GCNC
                target = target[len(res[0])+1:]
            elif target_t=="MPC":
                target_type = self.TARGET_TYPE_MPC
                target = target[len(res[0])+1:]
            elif target_t=="TLE":
                target_type = self.TARGET_TYPE_TLE
                target = target[len(res[0])+1:]
        # ---
        if target_type_only:
            if target_type == self.TARGET_TYPE_NAME:
                target_t = "NAME"
            return target_t
        # ---
        self._date2ts(date)
        # ---
        if target_type == self.TARGET_TYPE_MPC:
            start = Date(date).iso()
            location = ( f"{self.home.longitude}d", f"{self.home.latitude}d", f"{self.home.altitude:.1f}m")
            try:
                table = MPC.get_ephemeris(target, location=location, number=1, start=start)
            except:
                table = None
            if table != None:
                colnames = table.colnames
                data= table.as_array()[0]
                index = colnames.index('RA')
                ra = data[index]
                index = colnames.index('Dec')
                dec = data[index]
                found = True
        # ---
        if target_type == self.TARGET_TYPE_RADEC:
            c = SkyCoord(target.lower(), unit=(u.hourangle, u.deg))
            ra, dec = c.to_string("hmsdms").split()
        # ---
        if target_type == self.TARGET_TYPE_HADEC:
            time = Time(epoch)
            location = EarthLocation(lat=self.home.latitude*u.deg, lon=self.home.longitude*u.deg, height=self.home.altitude*u.m)
            c = SkyCoord(target.lower(), unit=(u.hourangle, u.deg), frame="hadec", obstime = time, location=location)
            ra, dec = c.icrs.to_string("hmsdms").split()
        # ---
        if target_type == self.TARGET_TYPE_AZELEV:
            time = Time(epoch)
            location = EarthLocation(lat=self.home.latitude*u.deg, lon=self.home.longitude*u.deg, height=self.home.altitude*u.m)
            c = SkyCoord(target.lower(), unit=(u.deg, u.deg), frame="altaz", obstime = time, location=location)
            ra, dec = c.icrs.to_string("hmsdms").split()
        # ---
        if target_type == self.TARGET_TYPE_RADECDRIFT:
            res = target.split()
            # last two elements are dra, ddec. Not used here.
            target = ""
            for k in range(len(res)-2):
                target += res[k] + " "
            target = target.strip()
            c = SkyCoord(target.lower(), unit=(u.hourangle, u.deg))
            ra, dec = c.to_string("hmsdms").split()
        # ---
        if target_type == self.TARGET_TYPE_GCNC:
            res = self.gcnc_download(target)
            ra, dec, equinox, epoch = res
        # ---
        if target_type == self.TARGET_TYPE_DATERADECS:
            targets = target.split("\n")
            # you must separate the date to coordinates by an isolated ' : '
            # 2022-01-31T22:34:00 : 12 34 32.12 -01 45 46.2
            # 2022-01-31T22:35:00 : 12 34 33.23 -01 56 04.5
            jd0 = Date(date).jd()
            jds = []
            radecs = []
            for target in targets:
                k = target.find(" : ")
                if k == -1:
                    continue
                date, radec = target.split(" : ")
                print(f"date={date} radec={radec}")
                jd = Date(date).jd()
                jds.append(jd)
                c = SkyCoord(radec.lower(), unit=(u.hourangle, u.deg))
                radec = c.to_string("hmsdms")
                radecs.append(radec)
            # sort the vector radecs with jds increasing
            inds = np.argsort(jds)
            sorted_jds = []
            sorted_radecs = []
            for i in inds:
                sorted_jds.append(jds[i])
                sorted_radecs.append(radecs[i])
            # place jd inside the range of jds if needed
            jd1 = sorted_jds[0]
            jd2 = sorted_jds[-1]
            if jd<jd1: jd = jd1
            if jd>jd2: jd = jd2
            # search the two jds values that are just before and just after jd0
            n = len(sorted_jds)
            if jd0 < sorted_jds[0]:
                # extrapol before
                k1 = 0
                k2 = 1
            elif jd0 > sorted_jds[n-1]:
                # extrapol after
                k1 = n-2
                k2 = n-1
            else:
                # interpol
                jd1 = sorted_jds[0]
                for k in range(1,n):
                    jd2 = sorted_jds[k]
                    if jd0>=jd1 and jd0<=jd2:
                        k1 = k-1
                        k2 = k
                        break
                    jd1 = jd2
            jd1 = sorted_jds[k1]
            jd2 = sorted_jds[k2]
            frac = (jd0 - jd1) / (jd2 - jd1)
            radec1 = sorted_radecs[k1]
            radec2 = sorted_radecs[k2]
            ra1, dec1 = radec1.split()
            ra2, dec2 = radec2.split()
            c1 = Coords((1, Angle(ra1), Angle(dec1)))
            c2 = Coords((1, Angle(ra2), Angle(dec2)))
            dra, ddec = c2.difference_angles_with_reference(c1)
            ra = Angle(ra1) + frac*dra
            dec = Angle(dec1) + frac*ddec
            ra = ra.deg()
            dec = dec.deg()
        # ---
        if target_type == self.TARGET_TYPE_NAME:
            found = False
            # --- Try a solar system planet
            if found == False:
                #print("==> Try solar system")
                ident = self.name2solar_system_planet(target)
                if ident != None:
                    planet = self._planets[ident]
                    astrometric = self._observatory.at(self._t).observe(planet)
                    ra, dec, distance = astrometric.radec()
                    #difference = planet - self._observatory
                    #topocentric = difference.at(self._t)
                    #elev, az, distance = topocentric.altaz()
                    ra = ra._degrees
                    dec = dec._degrees
                    found = True
            # --- Try a Simbad name
            if found == False:
                #print("==> Try Simbad")
                try:
                    table = Simbad.query_object(target)
                except:
                    pass
                if table != None:
                    colnames = table.colnames
                    data= table.as_array()[0]
                    index = colnames.index('RA')
                    ra = data[index]
                    index = colnames.index('DEC')
                    dec = data[index]
                    ra = Angle(ra).deg()*15
                    dec = Angle(dec).deg()
                    found = True
            # --- Try a MPC object
            if found == False:
                #print("==> Try MPC")
                start = Date(date).iso()
                location = ( f"{self.home.longitude}d", f"{self.home.latitude}d", f"{self.home.altitude:.1f}m")
                try:
                    table = MPC.get_ephemeris(target, location=location, number=1, start=start)
                except:
                    table = None
                if table != None:
                    colnames = table.colnames
                    data= table.as_array()[0]
                    index = colnames.index('RA')
                    ra = data[index]
                    index = colnames.index('Dec')
                    dec = data[index]
                    found = True
            # ---
            if found == False:
                msg = f"The target {target} was not found in Skyfield, Simbad and MPC databases"
                raise EphemerisException(EphemerisException.TARGET_NOT_FOUND, msg)
        # ---
        if target_type == self.TARGET_TYPE_TLEFILES:
            found = False
            # --- Try an Earth satellite
            target = str(target).upper()
            self.tle_download()
            by_names = {sat.name: sat for sat in self._earthsatellites}
            for name, satellite in by_names.items():
                # https://rhodesmill.org/skyfield/api-satellites.html
                norad = satellite.model.satnum
                international_designator = satellite.model.intldesg
                if name.find(target) >= 0:
                    found = True
                if target == str(norad):
                    found = True
                elif target == international_designator:
                    found = True
                if found == True:
                    break
            if found == False:
                msg = f"The target {target} was not found in Celestrack database"
                raise EphemerisException(EphemerisException.TARGET_NOT_FOUND, msg)
            # --- Compute ephemeris
            #print(f"Satel found = {name}")
            latitude = self.home.latitude
            longitude = self.home.longitude
            bluffton = wgs84.latlon(latitude, longitude)
            ts = skyfield.api.load.timescale()
            t = ts.now()
            difference = satellite - bluffton
            topocentric = difference.at(t)
            ra, dec, distance = topocentric.radec()  # ICRF ("J2000")
            ra = ra._degrees
            dec = dec._degrees
            #res = topocentric.speed()
        if target_type == self.TARGET_TYPE_TLE:
            found = False
            target = str(target).upper()
            lines = target.split("\n")
            latitude = self.home.latitude
            longitude = self.home.longitude
            bluffton = wgs84.latlon(latitude, longitude)
            ts = skyfield.api.load.timescale()
            t = ts.now()
            satellite = skyfield.api.EarthSatellite(lines[1], lines[2], lines[0], ts)
            difference = satellite - bluffton
            topocentric = difference.at(t)
            ra, dec, distance = topocentric.radec()  # ICRF ("J2000")
            ra = ra._degrees
            dec = dec._degrees
        # ---
        if unit_ra == "deg":
            ra = Angle(ra).deg()
        else:
            ra = Angle(ra).sexagesimal(unit_ra)
        if unit_dec == "deg":
            dec = Angle(dec).deg()
        else:
            dec = Angle(dec).sexagesimal(unit_dec)
        ra_equinox = ra
        dec_equinox = dec
         # --- dictionary
        eph = {}
        eph['header'] = {}
        eph['header']['home'] = self.home.gps
        eph['header']['target'] = target_init
        eph['header']['equinox'] = equinox
        eph['jd'] = Date(epoch).jd()
        eph['ra_equinox'] = ra_equinox
        eph['dec_equinox'] = dec_equinox
        return eph

    def radec_speed(self, target, **kwargs):
        """Return ra, dec J2000 of a target as method radec and add the speed (deg/s)

        Args:

            target: See radec for explanations
            **kwargs: See radec for explanations

        Returns:

            ra: Target Right Ascention position at the given date for a given equinox
            dec: Target Declination at the given date for a given equinox
            equinox: Equinox of the position
            epoch: Date of the position when the object is moving
            dra: Velocity on Right Ascension axis (deg/s)
            ddec: Velocity on Declination axis (deg/s)

        """
        ephem = self.radec(target, **kwargs)
        ra = ephem['ra_equinox']
        dec = ephem['dec_equinox']
        epoch1 = ephem['jd']
        target = str(target)
        res = target.split()
        target_t = res[0].upper()
        if target_t=="RADECDRIFT":
            target = target[len(res[0])+1:]
            res = target.split()
            ddec = float(res[-1])
            dra = float(res[-2])
        else:
            kwarg0s = kwargs.copy()
            if len(kwarg0s) == 0:
                kwarg0s = { "date":"now" }
            keys = kwarg0s.keys()
            if "date" in keys:
                date1 = kwarg0s["date"]
            else:
                date1 = "now"
            date_1 = Date(date1)
            dt = Duration("60s")
            date_2 = date_1 + dt.day()
            kwarg0s["unit_ra"] = "deg"
            kwarg0s["unit_dec"] = "deg"
            kwarg0s["date"] = date_1.iso(3)
            eph = self.radec(target, **kwarg0s)
            ra1 = eph['ra_equinox']
            dec1 = eph['dec_equinox']
            epoch1 = eph['jd']
            kwarg0s["date"] = date_2.iso(3)
            eph = self.radec(target, **kwarg0s)
            ra2 = eph['ra_equinox']
            dec2 = eph['dec_equinox']
            c1 = Coords((1, Angle(ra1), Angle(dec1)))
            c2 = Coords((1, Angle(ra2), Angle(dec2)))
            dra, ddec = c2.difference_angles_with_reference(c1)
            dra /= (86400*dt.day())
            ddec /= (86400*dt.day())
        eph = {}
        eph['header'] = {}
        eph['header'] = ephem['header']
        eph['jd'] = Date(epoch1).jd()
        eph['ra_equinox'] = Date(epoch1).jd()
        eph['ra_equinox'] = ra
        eph['dec_equinox'] = dec
        eph['dra_equinox'] = dra
        eph['ddec_equinox'] = ddec
        return eph

    def altitude2tp(self, alti:float, p0m:float=101325):
        """Compute the theoretical pressure and temperature in the Earth atmosphere given an altitude

        Args:

            alti: The altitude of the observation site in meters
            p0m: The pressure at the sea level. Default is 101325 Pascal.

        Returns:

            pressure: The pressure in Pascal
            temperature: The temperature in Kelvin
        """
        tk0m=273.15+15
        if alti<11000:
            tk=tk0m-0.0065*alti
            p=p0m*pow(tk/tk0m,5.255)
        elif alti<15000:
            tk0m=273.15+15
            tk=tk0m-0.0065*11000
            p=p0m*pow((tk0m-0.0065*alti)/tk0m,5.255)
        elif (alti>=15000) and (alti<20000):
            h1=15000; p1=p0m*pow((tk0m-0.0065*h1)/tk0m,5.255); t1=tk0m-0.0065*11000
            h2=20000; p2=5500; t2=273.15-46
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        elif (alti>=20000) and (alti<30000):
            h1=20000; p1=5500; t1=273.15-46
            h2=30000; p2=1100; t2=273.15-38
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        elif (alti>=30000) and (alti<40000):
            h1=30000; p1=1100; t1=273.15-38
            h2=40000; p2=300; t2=273.15-5
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        elif (alti>=40000) and (alti<50000):
            h1=40000; p1=300; t1=273.15-5
            h2=50000; p2=90; t2=273.15+1
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        elif (alti>=50000) and (alti<60000):
            h1=50000; p1=90; t1=273.15+1
            h2=60000; p2=25; t2=273.15-20
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        elif (alti>=60000) and (alti<100000):
            h1=60000; p1=25; t1=273.15-20
            h2=100000; p2=0.04; t2=273.15-64
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        elif (alti>=100000) and (alti<200000):
            h1=100000; p1=0.04; t1=273.15-64
            h2=200000; p2=1.3e-4; t2=273.15-82.2
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        elif (alti>=200000) and (alti<400000):
            h1=200000; p1=1.3e-4; t1=273.15-82.2
            h2=400000; p2=4.4e-6; t2=273.15-97.3
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        elif (alti>=200000) and (alti<400000):
            h1=400000; p1=4.4e-6; t1=273.15-97.3
            h2=500000; p2=0; t2=273.15-97.7
            frac=(alti-h1)/(h2-h1)
            p=p1+frac*(p2-p1)
            tk=t1+frac*(t2-t1)
        else:
            p=0
            tk=273.15-97.7
        return p, tk

    def date_ephem(self, ephem_night:dict, date:Date="now")->dict:
        """Extract the ephemeris for a given date asked from a night ephemeris.

        Arg:

            ephem_night: A night ephemeris returned by the method night_ephem
            date: The date of the calculation

        Returns:

            eph: The ephemeris at the date
            deph: The differential ephemeris in units of inverse of seconds (deg/s for angles)

        """
        # jd = jd0 + k*djd
        d = Date(date)
        jd = d.jd()
        jjds = ephem_night['jd']
        njd = len(jjds)
        djd = jjds[1]-jjds[0]
        n = round((jd-jjds[0])/djd)
        if n<0 or n>njd:
            msg = f"The date {d.iso(0)} is not inside the night {ephem_night['header']['night']}"
            raise EphemerisException(EphemerisException.DATE_OUTSIDE_THE_NIGHT, msg)
        eph = {}
        for key, val in ephem_night.items():
            if isinstance(val, np.ndarray):
                eph[key] = val[n] # deg
            else:
                eph[key] = val # str
        return eph

    def night_ephem(self, target, night:str, ephem_sun:dict=None, ephem_moon:dict=None, **kwargs)->dict:
        """Same as target2night but avoids to recompute many times the same night ephemeris if target is the same

        All args and returns are the same than target2night.
        """
        #
        # _computed_ephem["night"]
        # _computed_ephem["night"]["targets"] = List of targets (0..n-1)
        # _computed_ephem["night"][0] = targets index 0
        # _computed_ephem["night"][...]
        # _computed_ephem["night"][n-1] = targets index n-1
        if night in self._computed_ephem.keys():
            # --- this night is known
            targets = self._computed_ephem[night]["targets"]
            try:
                # --- case target is ever computer. Return the history instead to recompute
                indx = targets.index(target)
                ephem = self._computed_ephem[night][indx]
                return ephem
            except:
                # --- new target to be computed
                indx = len(targets)
        else:
            # --- this night is not known, clear all the history
            del self._computed_ephem
            self._computed_ephem = {}
            self._computed_ephem[night] = {}
            indx = 0
            self._computed_ephem[night]["targets"] = []
        # --- Compute the night ephemeris
        ephem = self.target2night(target, night, ephem_sun, ephem_moon, **kwargs)
        # --- Add the ephem into the history
        self._computed_ephem[night]["targets"].append(target)
        self._computed_ephem[night][indx] = ephem
        return ephem

    def target2night(self, target, night:str, ephem_sun:dict=None, ephem_moon:dict=None, **kwargs)->dict:
        """Compute the ephemeris at every second of local coodinates for a night

        Two computations are importants:

            * 'visibility': An integer defining why the target is visible or not.
            * 'observability': A float value from 0 (not observable) to 100 (best conditions to observe)

        Args:

            target: A target in the formalism of the method radec
            night: A night symbol (e.g. 20230330)
            ephem_sun: A dictionary result of the method called target2night("sun", night)
            ephem_moon: A dictionary result of the method called target2night("moon", night)
            kwargs: A dictionary of options:

                * horizon: An object of the class Horizon of Guitastro that defines the horizon line.
                * preference: A string "bestelev" or "immediate" to compute the observability
                * duskelev: A float, the elevation of the Sun defining the start and end of the night.
                * wavelength_nm: A float, the observation wavelength (in nanometers)
                * humidity: A float, the relative humidiy (between 0 and 1)
                * distmin_sun: A float to define the minimum angular distance from the Sun (degrees)
                * distmin_moon: A float to define the minimum angular distance from the Moon (degrees)
                * speed: A boolean, True to compute the derivative of coordinates
                * nsec: An integer, default is 86400 to compute for every seconds of the night. Else, compute only centered on the Date indicated in the night. Speed can be computed only for nsec >= 3.

        Returns:

            A dictionary:

                * 'night': String of the night
                * 'home': String of the GPS position
                * 'target': String of the input target
                * 'ndate': Number of dates used to compute amers before interpolations
                * 'jd': Numpy array of Julian days
                * 'alt': Numpy array of the elevation (deg)
                * 'az': Numpy array of the azimut (deg)
                * 'ha': Numpy array of the apparent hour angle (deg)
                * 'parallactic': Numpy array of the parallactic angle (deg)
                * 'dec': Numpy array of the apparent declination (deg)
                * 'ra_equinox': Numpy array of the Right Ascension at the input equinox (deg)
                * 'dec_equinox': Numpy array of the Declination at the input equinox (deg)
                * 'cosphi': Numpy array of cos(ra_equinox)
                * 'sinphi': Numpy array of sin(ra_equinox)
                * 'costheta': Numpy array of cos(dec_equinox)
                * 'sintheta': Numpy array of sin(dec_equinox)
                * 'distsun': Numpy array of the distance from the Sun (deg)
                * 'distmoon': Numpy array of the distance from the Moon (deg)
                * 'visibility': Numpy array of the visibility
                * 'observability': Numpy array of the observability (>0 means observable)
                * 'horizon': Numpy array of horizon elevations (deg)

        The Numpy arrays are 1D of 86400 elements.
        if ephem_sun==None or ephem_moon==None, the distances between the target and the Sun and Moon will not be calculated.

        """
        if ephem_sun==None or ephem_moon==None:
            sunmoon = False
        else:
            sunmoon = True
        if 'horizon' in kwargs.keys():
            hor = kwargs['horizon']
            hor_az, hor_elev = hor.horizon_altaz
        elif 'siteobs' in kwargs.keys():
            siteobs = kwargs['siteobs']
            hor_az, hor_elev = siteobs.horizon_altaz
        else:
            hor_az = np.arange(0, 361)
            hor_elev = np.zeros(len(hor_az))
        if 'preference' in kwargs.keys():
            preference = kwargs['preference'].lower()
        else:
            preference="bestelev"
        if 'duskelev' in kwargs.keys():
            duskelev = kwargs['duskelev']
        else:
            duskelev = -7
        if 'humidity' in kwargs.keys():
            rel_humidity = kwargs['humidity']
        else:
            rel_humidity = 0.6
        if 'wavelength_nm' in kwargs.keys():
            wavelength_nm = kwargs['wavelength_nm']
        else:
            wavelength_nm = 600
        if 'distmin_sun' in kwargs.keys():
            distmin_sun = float(kwargs['distmin_sun'])
        else:
            distmin_sun = 30
        if 'distmin_moon' in kwargs.keys():
            distmin_moon = float(kwargs['distmin_moon'])
        else:
            distmin_moon = 0
        if 'speed' in kwargs.keys():
            speed = kwargs['speed']
        else:
            speed = False
        if 'nsec' in kwargs.keys():
            nsec = int(kwargs['nsec'])
        else:
            nsec = 86400
        location = EarthLocation(lat=self.home.latitude*u.deg, lon=self.home.longitude*u.deg, height=self.home.altitude*u.m)
        tan_lat = np.tan(np.radians(self.home.latitude))
        temp_k, pres_pa = self.altitude2tp(self.home.altitude)
        temperature = (temp_k + 273.15) * u.deg_C
        pressure = pres_pa * u.pascal
        relative_humidity = rel_humidity
        obswl = wavelength_nm * u.nm
        fn = FileNames()
        fn.longitude(self.home.longitude)
        # --- compute jd1, jd2 the limits of dates to compute the ephemeris
        if nsec==86400:
            # - we compute ephemeris for all the night
            jd1, jd2 = fn.night2date(night)
        else:
            # - we compute ephemeris only for a duration centered on the date indicated by night
            jd = Date(night).jd()
            night = fn.date2night(jd)
            #nsec = round(nsec)
            if nsec <= 1:
                nsec = 1
            djd = (nsec-1)/86400.
            jd1, jd2 = jd-djd, jd+djd
        # --- identify the type of target
        target_type = self.radec(target, target_type_only=True)
        # --- compute the drift
        date = Date((jd1+jd2)/2).iso()
        speedephem = self.radec_speed(target, date=date, unit_ra="deg", unit_dec="deg")
        ra_equinox = speedephem['ra_equinox']
        dec_equinox = speedephem['dec_equinox']
        epoch = speedephem['jd']
        dra = speedephem['dra_equinox']
        ddec = speedephem['ddec_equinox']
        ddrift = np.sqrt(dra*dra+ddec*ddec) # deg/s
        # --- adapt ndate according the drift
        ndate = 0
        if dra==0 and ddec==0:
            drift = False
            epoch = Date((jd1+jd2)/2).iso()
            targ = SkyCoord(frame=ICRS, ra=ra_equinox*u.deg, dec=dec_equinox*u.deg, obstime=epoch)
        else:
            ndate = int(np.floor(86400*abs(ddrift)/10.0))
            drift = True
        #print(f"ndate={ndate} ddrift={ddrift}")
        # --- one computation every 30 min at minimum
        lim = 1440/30
        if ndate < lim:
            ndate = round(lim)
        # --- one computation every 5 min at maximum
        lim = 1440/5
        if ndate > lim:
            ndate = round(lim)
        # --- one computation every second at maximum
        if ndate > nsec:
            ndate = nsec
        # --- prepare angles
        jds = np.linspace(jd1, jd2, ndate)
        nangle = 13
        angles = np.zeros(nangle*ndate).reshape((nangle,ndate))
        angle_offsets = np.zeros(nangle)
        angle_prevs = np.zeros(nangle)
        angle_curs = np.zeros(nangle)
        # --- compute angles for a restricted number of dates
        for k in range(ndate):
            # --- compute celestial local angles
            jd = jds[k]
            obstime = Time(jd, format="jd")
            if drift == True:
                # --- recompute ra,dec for each date
                date = jd
                speedephem = self.radec_speed(target, date=date, unit_ra="deg", unit_dec="deg")
                ra_equinox = speedephem['ra_equinox']
                dec_equinox = speedephem['dec_equinox']
                epoch = Time(speedephem['jd'], format="jd")
                dra = speedephem['dra_equinox']
                ddec = speedephem['ddec_equinox']
                targ = SkyCoord(frame=ICRS, ra=ra_equinox*u.deg, dec=dec_equinox*u.deg, obstime=epoch)
            if target_type == "HADEC":
                target = str(target)
                res = target.split()
                mtarget = target[len(res[0])+1:]
                targ = SkyCoord(mtarget.lower(), unit=(u.hourangle, u.deg), frame="hadec", obstime = obstime, location=location)
                hadec= targ
            else:
                hadec = targ.transform_to(HADec(obstime=obstime, location=location, pressure=pressure, temperature=temperature, relative_humidity=relative_humidity, obswl=obswl))
            ha = hadec.ha.deg
            dec = hadec.dec.deg
            altaz = targ.transform_to(AltAz(obstime=obstime, location=location, pressure=pressure, temperature=temperature, relative_humidity=relative_humidity, obswl=obswl))
            alt = altaz.alt.deg
            az = altaz.az.deg
            az -= 180 # astro azimut instead geo
            ha_rad = np.radians(ha)
            dec_rad = np.radians(dec)
            y = np.sin(ha_rad)
            x = tan_lat * np.cos(dec_rad) - np.sin(dec_rad) * np.cos(ha_rad)
            parallactic = np.degrees(np.arctan2(y,x))
            ra_rad = np.radians(ra_equinox)
            dec_rad = np.radians(dec_equinox)
            cos_phi = np.cos(ra_rad)
            sin_phi = np.sin(ra_rad)
            cos_theta = np.cos(dec_rad)
            sin_theta = np.sin(dec_rad)
            if sunmoon==True:
                kk = int(np.round(k/ndate*nsec))
                # --- distance from the Sun
                _cos_phi = ephem_sun['cosphi'][kk]
                _sin_phi = ephem_sun['sinphi'][kk]
                _cos_theta = ephem_sun['costheta'][kk]
                _sin_theta = ephem_sun['sintheta'][kk]
                dx = _cos_theta * _cos_phi - cos_theta * cos_phi
                dy = _cos_theta * _sin_phi - cos_theta * sin_phi
                dz = _sin_theta - _sin_theta
                c = np.sqrt(dx*dx + dy*dy + dz*dz)
                dist_sun = np.degrees(2*np.arcsin(c/2))
                # --- distance from the Moon
                _cos_phi = ephem_moon['cosphi'][kk]
                _sin_phi = ephem_moon['sinphi'][kk]
                _cos_theta = ephem_moon['costheta'][kk]
                _sin_theta = ephem_moon['sintheta'][kk]
                dx = _cos_theta * _cos_phi - cos_theta * cos_phi
                dy = _cos_theta * _sin_phi - cos_theta * sin_phi
                dz = _sin_theta - _sin_theta
                c = np.sqrt(dx*dx + dy*dy + dz*dz)
                dist_moon = np.degrees(2*np.arcsin(c/2))
            else:
                dist_sun = 180
                dist_moon = 180
            # --- ensure continue angles
            angle_curs[0], angle_curs[1] = alt, az
            angle_curs[2], angle_curs[3] = ha, dec
            angle_curs[4], angle_curs[5] = ra_equinox, dec_equinox
            angle_curs[6], angle_curs[7] = cos_phi, sin_phi
            angle_curs[8], angle_curs[9] = cos_theta, sin_theta
            angle_curs[10], angle_curs[11] = dist_sun, dist_moon
            angle_curs[12] = parallactic
            # --- ensure continue angles
            if k>0:
                dif = angle_curs - angle_prevs
                for kk in range(nangle):
                    if dif[kk]>180:
                        angle_offsets[kk] -= 360
                    elif dif[kk]<-180:
                        angle_offsets[kk] += 360
            for kk in range(nangle):
                angles[kk][k] = angle_curs[kk] + angle_offsets[kk]
            angle_prevs = angle_curs.copy()
        # --- interpolate angles for a full number of dates
        jjds = np.linspace(jd1, jd2, nsec)
        aangles = np.zeros(nangle*nsec).reshape((nangle,nsec))
        if nsec == ndate:
            for kk in range(nangle):
                aangles[kk] = angles[kk]
        else:
            for kk in range(nangle):
                aangles[kk] = np.interp(jjds, jds, angles[kk])
        # --- interpolated angles
        alts = aangles[0]
        azs = aangles[1]
        has = aangles[2]
        decs = aangles[3]
        raequinoxs = aangles[4]
        decequinoxs = aangles[5]
        cosphis = aangles[6]
        sinphis = aangles[7]
        costhetas = aangles[8]
        sinthetas = aangles[9]
        distsuns = aangles[10]
        distmoons = aangles[11]
        parallactics = aangles[12]
        # --- compute speed angles if needed
        if speed and len(jjds)>=3:
            djd = jjds[1]-jjds[0]
            dt = djd*86400
            dt2 = dt*2
            dangles = np.zeros(nangle*nsec).reshape((nangle,nsec))
            for kk in range(nangle):
                y = aangles[kk]
                dangles[kk, 1:nsec-1] = (y[2:] - y[:-2])/dt2
                dangles[kk, 0] = (y[1]-y[0])/dt
                dangles[kk, nsec-1] = (y[nsec-1]-y[nsec-2])/dt
            dalts = dangles[0]
            dazs = dangles[1]
            dhas = dangles[2]
            ddecs = dangles[3]
            draequinoxs = dangles[4]
            ddecequinoxs = dangles[5]
            dparallactics = dangles[12]
        # --- observability and visibility
        visibilitys = np.zeros(nsec)
        observabilitys = np.zeros(nsec)
        horizons = np.zeros(nsec)
        # --- visibility
        kk = 0
        for alt, az in zip(alts, azs):
            k = int(np.floor(az))%360
            altmini = hor_elev[k]
            horizons[kk] = altmini
            if alt < altmini:
                visibilitys[kk] += 1
            if sunmoon==True:
                if ephem_sun['alt'][kk] > duskelev:
                    visibilitys[kk] += 2
            kk += 1
        visibilitys = np.where(distsuns > distmin_sun, visibilitys, visibilitys+4)
        visibilitys = np.where(distmoons > distmin_moon, visibilitys, visibilitys+8)
        # --- observability
        nvis = 0
        altmaxi = 0
        for kk in range(nsec):
            if visibilitys[kk] == 0:
                nvis += 1
                if alts[kk] > altmaxi:
                    altmaxi = alts[kk]
        if preference=="immediate":
            kvis = nvis
            for kk in range(nsec):
                if visibilitys[kk] == 0:
                    observabilitys[kk] = kvis/nvis*100
                    kvis -= 1
        if preference=="bestelev":
            for kk in range(nsec):
                if visibilitys[kk] == 0:
                    observabilitys[kk] = alts[kk]/altmaxi*100
         # --- dictionary
        eph = {}
        eph['header'] = {}
        eph['header']['night'] = night
        eph['header']['home'] = self.home.gps
        eph['header']['target'] = target
        eph['header']['ndate'] = ndate
        eph['header']['duskelev'] = duskelev
        eph['header']['preference'] = preference
        eph['header']['temperature_k'] = temp_k
        eph['header']['pressure_pa'] = pres_pa
        eph['header']['humidity_rel'] = rel_humidity
        eph['header']['wavelength_nm'] = wavelength_nm
        eph['jd'] = jjds
        eph['alt'] = alts
        eph['az'] = azs
        eph['parallactic'] = parallactics
        eph['ha'] = has
        eph['dec'] = decs
        eph['ra_equinox'] = raequinoxs
        eph['dec_equinox'] = decequinoxs
        eph['cosphi'] = cosphis
        eph['sinphi'] = sinphis
        eph['costheta'] = costhetas
        eph['sintheta'] = sinthetas
        eph['distsun'] = distsuns
        eph['distmoon'] = distmoons
        eph['visibility'] = visibilitys
        eph['observability'] = observabilitys
        eph['horizon'] = horizons
        if speed and len(jjds)>=3:
            eph['dalt'] = dalts
            eph['daz'] = dazs
            eph['dparallactic'] = dparallactics
            eph['dha'] = dhas
            eph['ddec'] = ddecs
            eph['dra_equinox'] = draequinoxs
            eph['ddec_equinox'] = ddecequinoxs
        return eph

    def hadec2ephem(self, ha, dec, date, **kwargs):
        if 'humidity' in kwargs.keys():
            rel_humidity = kwargs['humidity']
        else:
            rel_humidity = 0.6
        if 'wavelength_nm' in kwargs.keys():
            wavelength_nm = kwargs['wavelength_nm']
        else:
            wavelength_nm = 600
        # ---
        jd = Date(date).jd()
        obstime = Time(jd, format="jd")
        location = EarthLocation(lat=self.home.latitude*u.deg, lon=self.home.longitude*u.deg, height=self.home.altitude*u.m)
        temp_k, pres_pa = self.altitude2tp(self.home.altitude)
        temperature = (temp_k + 273.15) * u.deg_C
        pressure = pres_pa * u.pascal
        relative_humidity = rel_humidity
        obswl = wavelength_nm * u.nm
        # ---
        c = SkyCoord(frame="hadec", ha=ha*u.deg, dec=dec*u.deg, obstime=obstime, location=location)
        # ---
        altaz = c.transform_to(AltAz(obstime=obstime, location=location, pressure=pressure, temperature=temperature, relative_humidity=relative_humidity, obswl=obswl))
        alt = altaz.alt.deg
        az = altaz.az.deg
        az -= 180
        # ---
        radec = c.icrs
        ra_equinox, dec_equinox = radec.ra.deg, radec.dec.deg
        # ---
        ephem = {}
        ephem['header'] = {}
        ephem['header']['home'] = self.home.gps
        ephem['header']['temperature_k'] = temp_k
        ephem['header']['pressure_pa'] = pres_pa
        ephem['header']['humidity_rel'] = rel_humidity
        ephem['header']['wavelength_nm'] = wavelength_nm
        ephem['jd'] = jd
        ephem['ra_equinox'] = ra_equinox
        ephem['dec_equinox'] = dec_equinox
        ephem['ha'] = ha
        ephem['dec'] = dec
        ephem['az'] = az
        ephem['alt'] = alt
        return ephem

    def altaz2ephem(self, az, alt, date, **kwargs):
        if 'humidity' in kwargs.keys():
            rel_humidity = kwargs['humidity']
        else:
            rel_humidity = 0.6
        if 'wavelength_nm' in kwargs.keys():
            wavelength_nm = kwargs['wavelength_nm']
        else:
            wavelength_nm = 600
        # ---
        jd = Date(date).jd()
        obstime = Time(jd, format="jd")
        location = EarthLocation(lat=self.home.latitude*u.deg, lon=self.home.longitude*u.deg, height=self.home.altitude*u.m)
        temp_k, pres_pa = self.altitude2tp(self.home.altitude)
        temperature = (temp_k + 273.15) * u.deg_C
        pressure = pres_pa * u.pascal
        relative_humidity = rel_humidity
        obswl = wavelength_nm * u.nm
        # ---
        azg = az + 180
        c = SkyCoord(frame="altaz", az=azg*u.deg, alt=alt*u.deg, obstime=obstime, location=location)
        # ---
        hadec = c.transform_to(HADec(obstime=obstime, location=location, pressure=pressure, temperature=temperature, relative_humidity=relative_humidity, obswl=obswl))
        ha = hadec.ha.deg
        dec = hadec.dec.deg
        # ---
        radec = c.icrs
        ra_equinox, dec_equinox = radec.ra.deg, radec.dec.deg
        # ---
        ephem = {}
        ephem['header'] = {}
        ephem['header']['home'] = self.home.gps
        ephem['header']['temperature_k'] = temp_k
        ephem['header']['pressure_pa'] = pres_pa
        ephem['header']['humidity_rel'] = rel_humidity
        ephem['header']['wavelength_nm'] = wavelength_nm
        ephem['jd'] = jd
        ephem['ra_equinox'] = ra_equinox
        ephem['dec_equinox'] = dec_equinox
        ephem['ha'] = ha
        ephem['dec'] = dec
        ephem['az'] = az
        ephem['alt'] = alt
        return ephem

# #####################################################################
# #####################################################################
# #####################################################################
# Main
# #####################################################################
# #####################################################################
# #####################################################################

if __name__ == "__main__":

    default = 14
    example = input(f"Select the example (0 to 14) ({default}) ")
    try:
        example = int(example)
    except:
        example = default
    print("Example       = {}".format(example))

    if example == 1:
        """
        Simple ephemeris of a planet
        """
        eph = Ephemeris()
        name = "neptune"
        ephem = eph.radec(name, date="now", unit_ra="H0.2", unit_dec="d+090.1")
        print(f"{name} ra={ephem['ra_equinox']} dec={ephem['dec_equinox']} equinox={ephem['header']['equinox']} epoch={ephem['jd']}")

    if example == 2:
        """
        Simple ask to Simbad
        """
        eph = Ephemeris()
        name = "UGC 3462"
        ephem = eph.radec(name)
        print(f"{name} ra={ephem['ra_equinox']} dec={ephem['dec_equinox']} equinox={ephem['header']['equinox']} epoch={ephem['jd']}")

    if example == 3:
        """
        Simple ask to MPC
        """
        eph = Ephemeris()
        eph.set_home("guitalens")
        name = "2022 AB"
        ephem = eph.radec(name, date="now")
        print(f"{name} ra={ephem['ra_equinox']} dec={ephem['dec_equinox']} equinox={ephem['header']['equinox']} epoch={ephem['jd']}")

    if example == 4:
        """
        Simple ask from CELESTRACK
        """
        eph = Ephemeris()
        eph.set_home("guitalens")
        nsat = eph.tle_download()
        print(f"TLE files contain {nsat} satellites.")
        # --- First method
        name = "ISS"
        ephem = eph.radec(name, date="now", target_type="celestrack")
        print(f"1st method: ra={ephem['ra_equinox']} dec={ephem['dec_equinox']} equinox={ephem['header']['equinox']} epoch={ephem['jd']}")
        # --- Second method
        name = "CELESTRACK ISS"
        ephem = eph.radec(name, date="now")
        print(f"2nd method: ra={ephem['ra_equinox']} dec={ephem['dec_equinox']} equinox={ephem['header']['equinox']} epoch={ephem['jd']}")

    if example == 5:
        """
        Simple ask ISS with speed
        """
        eph = Ephemeris()
        name = "CELESTRACK ISS"
        ephem = eph.radec_speed(name)
        print(f"{name} ra={ephem['ra_equinox']} dec={ephem['dec_equinox']} dra={ephem['dra_equinox']*3600:.3f} ddec={ephem['ddec_equinox']*3600:.3f} equinox={ephem['header']['equinox']} epoch={ephem['jd']}")

    if example == 6:
        """
        Simple ask from TLE
        """
        eph = Ephemeris()
        eph.set_home("guitalens")
        name = "GSAT0203 (PRN E26)"
        tle  = f"TLE {name}\n"
        tle += "1 40544U 15017A   22029.59832018 -.00000064  00000+0  00000+0 0  9991\n"
        tle += "2 40544  56.7391  25.1641 0005059 269.6475  90.3121  1.70475734 42591_n"
        ephem = eph.radec_speed(tle)
        print(f"{name} ra={ephem['ra_equinox']} dec={ephem['dec_equinox']} dra={ephem['dra_equinox']*3600:.3f} ddec={ephem['ddec_equinox']*3600:.3f} equinox={ephem['header']['equinox']} epoch={ephem['jd']}")

    if example == 7:
        """
        Simple ask to DATERADECS
        """
        eph = Ephemeris()
        eph.set_home("guitalens")
        name = "My object"
        target  = "DATERADECS "
        target += "2022-01-31T22:34:00 : 12 34 32.12 -01 45 46.2\n"
        target += "2022-01-31T22:35:00 : 12 34 33.23 -01 56 04.5\n"
        ephem = eph.radec_speed(target, date="2022-01-31T22:34:15", unit_ra="H0.2", unit_dec="d+090.1")
        print(f"{name} ra={ephem['ra_equinox']} dec={ephem['dec_equinox']} dra={ephem['dra_equinox']*3600:.3f} ddec={ephem['ddec_equinox']*3600:.3f} equinox={ephem['header']['equinox']} epoch={ephem['jd']}")

    if example == 8:
        """
        Simple ask to GCNC
        """
        eph = Ephemeris()
        name = "220412A"
        target = f"GCNC {name}"
        ephem = eph.radec(target, unit_ra="H0.2", unit_dec="d+090.1")
        print(f"{name} ra={ephem['ra_equinox']} dec={ephem['dec_equinox']} equinox={ephem['header']['equinox']} epoch={ephem['jd']}")

    if example == 9:
        """
        Simple ask to Simbad
        """
        eph = Ephemeris()
        #
        #names = ["AF Vir", "DR And", "HH Aqr", "LN Boo", "TU Com", "TU Per", "V348 Vir"]
        names = ["OX Aqr"]
        for name in names:
            ephem = eph.radec(name)
            print(f"{name} ra={ephem['ra_equinox']} dec={ephem['dec_equinox']} equinox={ephem['header']['equinox']} epoch={ephem['jd']}")

    if example == 10:
        """
        Simple ask to DATERADECDRIFT
        """
        eph = Ephemeris()
        eph.set_home("guitalens")
        name = "My object"
        target  = "RADECDRIFT 12 34 32.12 -01 45 46.2 1 2"
        #target  = "RADEC 0H10M -16d "
        ephem = eph.radec_speed(target)
        print(f"{name} ra={ephem['ra_equinox']:.6f} dec={ephem['dec_equinox']:.6f} dra={ephem['dra_equinox']*3600:.5f} ddec={ephem['ddec_equinox']*3600:.5f} equinox={ephem['header']['equinox']} epoch={ephem['jd']}")

    if example == 11:
        """
        Simple ask to HADEC
        """
        eph = Ephemeris()
        eph.set_home("guitalens")
        name = "My object"
        target  = "HADEC 12 34 32.12 -01 45 46.2"
        #target_type = eph.radec(target, target_type_only=True)
        #ra, dec, equinox, epoch= eph.radec(target)
        #print(f"{name} ra={ra:.6f} dec={dec:.6f} equinox={equinox} epoch={epoch}")
        ephem = eph.radec_speed(target)
        print(f"{name} ra={ephem['ra_equinox']:.6f} dec={ephem['dec_equinox']:.6f} dra={ephem['dra_equinox']*3600:.5f} ddec={ephem['ddec_equinox']*3600:.5f} equinox={ephem['header']['equinox']} epoch={ephem['jd']}")

    if example == 12:
        """
        Compute the ephemeris of a target along all a night
        """
        import matplotlib.pyplot as plt
        def compute_hours(ephem):
            jd0 = ephem['jd'][0]
            frac = jd0 - np.floor(jd0)
            offset = frac - 0.5
            hours = (ephem['jd'] - jd0 + offset)*24
            date = Date(jd0).iso(0)
            return hours, date
        # ---
        eph = Ephemeris()
        eph.set_home("guitalens")
        nsec = 86400
        night = "20230320"
        preference = "bestelev"
        duskelev = -7
        #preference = "immediate"
        # --- horizon
        from horizon import Horizon
        hor = Horizon(eph.home)
        hor.horizon_altaz = [(0,40), (180,0), (360,40)]
        hor_az, hor_elev = hor.horizon_altaz
        # --- Test nsec
        if False:
            nsec = 3
            night = "2023-03-20T12:00:00"
        # --- sun
        target = "sun"
        t0 = time.time()
        ephem_sun = eph.target2night(target, night, None, None, nsec=nsec)
        dt = time.time()-t0
        print(f"SUN dt={dt}")
        # --- moon
        target = "moon"
        t0 = time.time()
        ephem_moon = eph.target2night(target, night, None, None, nsec=nsec)
        dt = time.time()-t0
        print(f"MOON dt={dt}")
        # --- target
        target = "RADEC 4h56m -12d23m"
        #target = "HADEC 2h +16d"
        t0 = time.time()
        speed = True
        ephem = eph.target2night(target, night, ephem_sun, ephem_moon, horizon=hor, preference=preference, duskelev=duskelev, speed=speed, nsec=nsec)
        dt = time.time()-t0
        print(f"TARGET dt={dt}")
        hours, date = compute_hours(ephem)
        plt.plot(hours, ephem_sun['alt'], "y-")
        plt.plot(hours, ephem['alt'], "b-")
        plt.plot(hours, ephem['horizon'], "g-")
        plt.plot(hours, ephem['visibility'], "k:")
        plt.plot(hours, ephem['observability'], "k-")
        plt.grid()
        plt.ylabel('Degrees')
        plt.xlabel(f'Hours since {date}')

    if example == 13:
        """
        Transform HADEC or ALTAZ to ephem (a dict of all coordinates)
        """
        eph = Ephemeris()
        eph.set_home("guitalens")
        ephem = eph.altaz2ephem(-90, 0, "now")
        print(ephem)

    if example == 14:
        """
        Compute the ephemeris for PYROS test
        """
        import matplotlib.pyplot as plt
        def compute_hours(ephem):
            jd0 = ephem['jd'][0]
            frac = jd0 - np.floor(jd0)
            offset = frac - 0.5
            hours = (ephem['jd'] - jd0 + offset)*24
            date = Date(jd0).iso(0)
            return hours, date
        # ---
        eph = Ephemeris()
        eph.set_home("GPS 2.0375 E 43.6443484725 136.9")
        nsec = 86400
        night = "now"
        preference = "bestelev"
        duskelev = -7
        #preference = "immediate"
        # --- horizon
        from horizon import Horizon
        hor = Horizon(eph.home)
        hor.horizon_altaz = [ [0,0], [360,0] ]
        hor_az, hor_elev = hor.horizon_altaz
        # --- sun
        target = "sun"
        t0 = time.time()
        ephem_sun = eph.target2night(target, night, None, None, nsec=nsec)
        dt = time.time()-t0
        print(f"SUN dt={dt}")
        # --- moon
        target = "moon"
        t0 = time.time()
        ephem_moon = eph.target2night(target, night, None, None, nsec=nsec)
        dt = time.time()-t0
        print(f"MOON dt={dt}")
        # --- target
        target = "RADEC 18h -15d"
        t0 = time.time()
        speed = True
        ephem = eph.target2night(target, night, ephem_sun, ephem_moon, horizon=hor, preference=preference, duskelev=duskelev, speed=speed, nsec=nsec)
        #ephem = eph.target2night(target, night, None, None, preference=preference, duskelev=duskelev, speed=speed, nsec=nsec)
        dt = time.time()-t0
        print(f"TARGET dt={dt}")
        hours, date = compute_hours(ephem)
        plt.plot(hours, ephem_sun['alt'], "y-")
        plt.plot(hours, ephem['distmoon'], "c-")
        plt.plot(hours, ephem['alt'], "b-")
        plt.plot(hours, ephem['horizon'], "g-")
        plt.plot(hours, ephem['visibility'], "k:")
        plt.plot(hours, ephem['observability'], "k-")
        plt.grid()
        plt.ylabel('Degrees')
        plt.xlabel(f'Hours since {date}')