# -*- coding: utf-8 -*-
import os
import time
import math
import numpy as np
import matplotlib.pyplot as plt

# --- celme imports
modulename = 'celme'
if modulename in dir():
    del celme
if modulename not in dir():    
    import celme

from mountaxis import Mountaxis
from mounttools import Mounttools
from mountlog import Mountlog
from mountpad import Mountpad

# #####################################################################
# #####################################################################
# #####################################################################
# Class Mountastro
# #####################################################################
# #####################################################################
# This is an abstract Class
# This class does not communicate to devices
# Use a child class to use the protocol language
# #####################################################################

class Mountastro(Mounttools):
    
    # === Constant for error codes
    NO_ERROR = 0
    ERR_UNKNOWN = 1
    
    ERR_FILE_NOT_EXISTS = 101
    
    # === Constants
    
    # === Private variables
    _last_errno = NO_ERROR
    
    # --- Axis parameters
    axisp = None
    axisb = None

    # === constants for saving coords
    SAVE_NONE = 0
    SAVE_AS_SIMU = 1
    SAVE_AS_REAL = 2
    SAVE_ALL = 3
    
    # === Constants for returning long or short outputs
    OUTPUT_SHORT = 0
    OUTPUT_LONG = 1
    
    # ===
    LX200_0X06 = 0x06
    LX200_HIGH_PRECISION = "HIGH PRECISION"
    LX200_LOW_PRECISION  = "LOW PRECISION"
    
    # === List of Threads for the pad GUI
    _threads = []
    
    _drift_hadec_ha_deg_per_sec = 0
    _drift_hadec_dec_deg_per_sec = 0
    _drift_radec_ra_deg_per_sec = 0
    _drift_radec_dec_deg_per_sec = 0

    _sideral_sec_per_day = 86164.0

    # --- last actions
    _last_goto_drift_deg_per_secs = []
    
# =====================================================================
# =====================================================================
# Private methods
# =====================================================================
# =====================================================================
    
    def _create(self, mount_name, mount_type):
        self._delete()

    def _delete(self):
        self._last_errno = self.NO_ERROR

# =====================================================================
# =====================================================================
# Methods for using a pad GUI
# =====================================================================
# =====================================================================
# to use it:
# mymount = Mountastro("HADEC", name="Test Mount")
# try:
#     mymount.pad_create("pad_dev1")
# except (KeyboardInterrupt, SystemExit):
#     mymount.pad_delete()
# except:
#     raise  
# =====================================================================

    def pad_create(self, pad_type):
        pad = Mountpad(self, pad_type)
        pad.start() # start the thread and continue
        self._pads.append(pad)
        self.pad = pad # to share the last pad

    def pad_delete(self):
        # --- kill gui and associated threads
        for pad in self._pads:
            try:
                self.pad.pad_gui_delete()
            except:
                pass
        self._pads = []
            
# =====================================================================
# =====================================================================
# Methods for debug
# =====================================================================
# =====================================================================

    def disp(self):
        incsimus = ["" for kaxis in range(Mountaxis.AXIS_MAX)]  
        increals, incsimus = self.read_encs(incsimus)
        increalb, rotrealb, celrealb, increalp, rotrealp, celrealp, piersidereal, incsimub, rotsimub, celsimub, incsimup, rotsimup, celsimup, piersidesimu = self.enc2cel(incsimus, output_format=self.OUTPUT_LONG, save=self.SAVE_ALL)
        if self._mount_type=="HADEC":
            hareal, decreal = self.cel2hadec(celrealb, celrealp, "H0.2", "d+090.1")
            hasimu, decsimu = self.cel2hadec(celsimub, celsimup, "H0.2", "d+090.1")
        if self._mount_type=="AZELEV":
            azreal, elevreal = self.cel2azelev(celrealb, celrealp, "D0.2", "d+090.1")
            azsimu, elevsimu = self.cel2azelev(celsimub, celsimup, "D0.2", "d+090.1")
        print("{}".format(20*"-"))
        print("MOUNT name        = {} ".format(self._name))
        print("mount_type        = {} ".format(self._mount_type))
        print("home              = {} ".format(self.site.gps))
        for kaxis in range(Mountaxis.AXIS_MAX):
            current_axis = self.axis[kaxis]
            if current_axis == None:
                continue
            for disp_real in (False,True):
                if disp_real==True and current_axis.real == False:
                    continue
                # --- get inc and compute rot, cel
                if disp_real==True:
                    msg_simu = "REAL"
                    inc = increals[kaxis]
                else:
                    msg_simu = "SIMU"
                    inc = incsimus[kaxis]
                rot, pierside = current_axis.inc2rot(inc, self.SAVE_NONE)
                cel = current_axis.rot2ang(rot, pierside, self.SAVE_NONE)
                # --- additional informations
                msg_coord = ""
                if self._mount_type=="HADEC":
                    if current_axis.axis_type=="HA":
                        if disp_real==True:
                            msg_coord = hareal
                        else:
                            msg_coord = hasimu
                    if current_axis.axis_type=="DEC":
                        if disp_real==True:
                            msg_coord = decreal
                        else:
                            msg_coord = decsimu
                if self._mount_type=="AZELEV":
                    if current_axis.axis_type=="AZ":
                        if disp_real==True:
                            msg_coord = azreal
                        else:
                            msg_coord = azsimu
                    if current_axis.axis_type=="ELEV":
                        if disp_real==True:
                            msg_coord = elevreal
                        else:
                            msg_coord = elevsimu
                # ---
                print("{} {} {} {}".format(20*"-",current_axis.name,msg_simu,current_axis.axis_type))
                print("inc               = {:12.1f} : ".format(inc))
                print("rot               = {:12.7f} : ".format(rot))
                print("pierside          = {:d}            : ".format(pierside))
                print("cel               = {:12.7f} : {} {}".format(cel,current_axis.axis_type,msg_coord))
        
# =====================================================================
# =====================================================================
# Methods to read the encoders
# =====================================================================
# =====================================================================
# Level 1
# =====================================================================
                
    def _my_read_encs(self, incsimus:list)->list:
        """
        Inputs are simulated inc
        Outputs are real inc
        """
        # --- abstract method. 
        # --- Please overload it according your language protocol
        increals = incsimus
        return increals
    
    def read_encs(self, simulation_incs:list)->list:
        """
        Read the raw values of the axes
        
        For the simulation: 
            if simulation_incb=="" the value is calculated from simu_update_inc
            if simulation_incb==34.5 the value is taken equal to 34.5

        For real:
            if the axis is not real then real=simulated value
            else the encoder is read

        Output = Raw calues of encoders (inc)
        """
        # === Simulated values of inc
        incsimus = [0 for kaxis in range(Mountaxis.AXIS_MAX)]
        # --- Loop over all the possible axis types
        for kaxis in range(Mountaxis.AXIS_MAX):
            current_axis = self.axis[kaxis]
            if current_axis == None:
                continue
            # --- This axis is valid. We compute the simulation
            simulation_inc = simulation_incs[kaxis]
            if isinstance(simulation_inc,(int,float))==True:
                incsimus[kaxis] = simulation_inc
            else:
                incsimus[kaxis] = current_axis.simu_update_inc()
        # === Read real values of inc if possible
        increals = self._my_read_encs(incsimus)
        # ---
        if self._record_positions==True:
            with open("positions.txt","a",encoding='utf-8') as fid:
                msg = "{} {} {}\n".format(time.time(), increals[Mountaxis.BASE], incsimus[Mountaxis.BASE])
                fid.write(msg)
        # ---
        return (increals, incsimus)

# =====================================================================
# =====================================================================
# Methods to convert encoders into celestial
# =====================================================================
# =====================================================================
# Level 2
# =====================================================================
        
    def enc2cel(self, simulation_incs:list, output_format=OUTPUT_SHORT, save=SAVE_NONE)->tuple:
        """
        Read encoder values into celestial apparent coordinates
        
        Input  = Raw encoder values (inc)
        Output = HA, Dec, pier side (any celme Angle units)
        """
        increals, incsimus = self.read_encs(simulation_incs)
        # --- shortcuts
        axisb = self.axis[Mountaxis.BASE]
        axisp = self.axis[Mountaxis.POLAR]
        incsimub = incsimus[Mountaxis.BASE]
        incsimup = incsimus[Mountaxis.POLAR]
        increalb = increals[Mountaxis.BASE]
        increalp = increals[Mountaxis.POLAR]
        # --- update the axis parameters (start with polar axis to deduce pierside)
        # --- update for simulations
        if save==self.SAVE_ALL or save==self.SAVE_AS_SIMU:
            savesimu = self.SAVE_AS_SIMU
        else:
            savesimu = self.SAVE_NONE
        rotsimup, piersidesimu = axisp.inc2rot(incsimup, savesimu)
        celsimup = axisp.rot2ang(rotsimup, piersidesimu, savesimu)
        rotsimub, dummy = axisb.inc2rot(incsimub, savesimu)
        celsimub = axisb.rot2ang(rotsimub, piersidesimu, savesimu)
        # --- update for real
        if save==self.SAVE_ALL or save==self.SAVE_AS_REAL:
            savereal = self.SAVE_AS_REAL
        else:
            savereal = self.SAVE_NONE
        rotrealp, piersidereal = axisp.inc2rot(increalp, savereal)
        celrealp = axisp.rot2ang(rotrealp, piersidereal, savereal)
        rotrealb, dummy = axisb.inc2rot(increalb, savereal)
        celrealb = axisb.rot2ang(rotrealb, piersidereal, savereal)
        # --- select results as simulation or real
        if axisp.real == False:
            celp = celsimup
            pierside = piersidesimu
        else:
            celp = celrealp
            pierside = piersidereal
        if axisb.real == False:
            celb = celsimub
        else:
            celb = celrealb
        # ---
        if output_format==self.OUTPUT_SHORT:
            return celb, celp, pierside
        else:
            # all angle output
            return increalb, rotrealb, celrealb, increalp, rotrealp, celrealp, piersidereal, incsimub, rotsimub, celsimub, incsimup, rotsimup, celsimup, piersidesimu

    def cel2enc(self, celb:float, celp:float, pierside:int, output_format=OUTPUT_SHORT, save=SAVE_NONE):
        """
        Convert celestial apparent coordinates into encoder values
        
        Input  = celb, celp, pier side (deg units)
        Output = Raw encoder values (inc)
        """
        # --- shortcuts
        axisb = self.axis[Mountaxis.BASE]
        axisp = self.axis[Mountaxis.POLAR]
        # --- update for simulations
        if save==self.SAVE_ALL or save==self.SAVE_AS_SIMU:
            savesimu = self.SAVE_AS_SIMU
        else:
            savesimu = self.SAVE_NONE
        celsimub = celb
        celsimup = celp
        piersidesimu = pierside
        rotsimup = axisp.ang2rot(celsimup, pierside, savesimu)
        incsimup = axisp.rot2inc(rotsimup, savesimu)
        rotsimub = axisb.ang2rot(celsimub, pierside, savesimu)
        incsimub = axisb.rot2inc(rotsimub, savesimu)
        # --- update for real
        if save==self.SAVE_ALL or save==self.SAVE_AS_REAL:
            savereal = self.SAVE_AS_REAL
        else:
            savereal = self.SAVE_NONE
        celrealb = celb
        celrealp = celp
        piersidereal = pierside
        rotrealp = axisp.ang2rot(celrealp, pierside, savereal)
        increalp = axisp.rot2inc(rotrealp, savereal)
        rotrealb = axisb.ang2rot(celrealb, pierside, savereal)
        increalb = axisb.rot2inc(rotrealb, savereal)
        # --- select results as simulation or real
        if axisp.real == False:
            incp = incsimup
        else:
            incp = increalp
        if axisb.real == False:
            incb = incsimub
        else:
            incb = increalb
        # --- update the output axis parameters
        if output_format==self.OUTPUT_SHORT:
            # short output
            return incb, incp
        else:
            # all angle output
            return increalb, rotrealb, celrealb, increalp, rotrealp, celrealp, piersidereal, incsimub, rotsimub, celsimup, incsimup, rotsimup, celsimup, piersidesimu

    def speedslew(self, *args):
        argc = len(args)
        karg = 0
        speedslew = ()
        for kaxis in range(Mountaxis.AXIS_MAX):
            current_axis = self.axis[kaxis]
            if current_axis == None:
                continue
            # --- read or update the speed
            if karg<argc:
                current_axis.slew_deg_per_sec = abs(args[karg]) # slew speed is always positive
            speedslew += (current_axis.slew_deg_per_sec,)
        return speedslew

# =====================================================================
# =====================================================================
# Methods to convert celestial into astronomical coordinates
# =====================================================================
# =====================================================================
# Level 3
# =====================================================================

    def hadec2cel(self, ha:celme.Angle, dec:celme.Angle)->tuple:
        # ---
        ha = celme.Angle(ha).deg()
        dec = celme.Angle(dec).deg()
        # ---
        meca = celme.Mechanics()
        if self._mount_type.find("HADEC")>=0:
            # --- astro coord are HADEC and cel are HADEC
            celb = ha
            celp = dec
            celr = 0
        if self._mount_type.find("AZELEV")>=0:
            # --- astro coord are HADEC and cel are AZELEV
            ha *= self._d2r
            dec *= self._d2r
            latitude = celme.Angle(self.site.latitude).rad()
            az, elev = meca._mc_hd2ah(ha, dec, latitude)
            rotator = meca._mc_hd2parallactic(ha, dec, latitude)
            celb = az * self._r2d
            celp = elev * self._r2d
            celr = rotator * self._r2d
        # ---
        return celb, celp, celr

    def azelev2cel(self, az:celme.Angle, elev:celme.Angle)->tuple:
        # ---
        az = celme.Angle(az).deg()
        elev = celme.Angle(elev).deg()
        # ---
        meca = celme.Mechanics()
        if self._mount_type.find("HADEC")>=0:
            # --- astro coord are AZELEV and cel are HADEC
            az *= self._d2r
            elev *= self._d2r
            latitude = celme.Angle(self.site.latitude).rad()
            ha, dec= meca._mc_ah2hd(az, elev, latitude)
            rotator = meca._mc_hd2parallactic(ha, dec, latitude)
            celb = ha * self._r2d
            celp = dec * self._r2d
            celr = rotator * self._r2d
        if self._mount_type.find("AZELEV")>=0:
            # --- astro coord are AZELEV and cel are AZELEV
            celb = az
            celp = elev
            celr = 0
        # ---
        return celb, celp, celr

    def astro2cel(self, astro_type:str, base:celme.Angle, polar:celme.Angle, base_deg_per_sec:str="", polar_deg_per_sec:str="")->tuple:
        """        
        drift_type = "HADEC" or "AZELEV"
        base = ha or az (according astro_type)
        polar = dec or elev (according astro_type)
        base_deg_per_sec = dha or daz (according astro_type)
        polar_deg_per_sec = ddec or delev (according astro_type)
        """
        base  = celme.Angle(base).deg()
        polar = celme.Angle(polar).deg()
        # --- special case when the drifts are diurnal
        if base_deg_per_sec=="diurnal" or polar_deg_per_sec=="diurnal":
            # --- compute ha,dec
            if astro_type.find("AZELEV")>=0:
                meca = celme.Mechanics()
                # --- here all angles are in radians
                latitude = celme.Angle(self.site.latitude).rad()
                az = base * self._d2r
                elev = polar * self._d2r
                ha, dec = meca._mc_ah2hd(az, elev, latitude)
                ha *= self._r2d
                dec *= self._r2d
            if astro_type.find("HADEC")>=0:
                ha = base
                dec = polar
            # ---
            if base_deg_per_sec=="diurnal":
                dha = 360.0/self._sideral_sec_per_day
            else:
                dha = base_deg_per_sec
            if base_deg_per_sec=="diurnal":
                ddec = 0
            else:
                ddec = polar_deg_per_sec
            astro_type = "HADEC"
            base = ha
            polar = dec
            base_deg_per_sec = dha
            polar_deg_per_sec = ddec
        # --- compute the derivative
        dt = 1.0 # sec
        base1  = base  - 0.5*dt*base_deg_per_sec
        polar1 = polar - 0.5*dt*polar_deg_per_sec
        base2  = base  + 0.5*dt*base_deg_per_sec
        polar2 = polar + 0.5*dt*polar_deg_per_sec
        if astro_type.find("HADEC")>=0:
            celb1, celp1, celr1 = self.hadec2cel(base1, polar1)
            celb2, celp2, celr2 = self.hadec2cel(base2, polar2)
            celb, celp, celr = self.hadec2cel(base, polar)
        if astro_type.find("AZELEV")>=0:
            celb1, celp1, celr1 = self.azelev2cel(base1, polar1)
            celb2, celp2, celr2 = self.azelev2cel(base2, polar2)
            celb, celp, celr = self.azelev2cel(base, polar)
        # ---
        dcelb = (celb2-celb1)/dt
        dcelp = (celp2-celp1)/dt
        dcelr = (celr2-celr1)/dt
        return celb, celp, celr, dcelb, dcelp, dcelr
    
    def cel2astro(self, celb:float, celp:float, unit_ha:str="", unit_dec:str="", unit_az:str="", unit_elev:str="")->tuple:
        meca = celme.Mechanics()
        # --- here all angles are in radians
        latitude = celme.Angle(self.site.latitude).rad()
        if self._mount_type.find("HADEC")>=0:
            ha = celb * self._d2r
            dec = celp * self._d2r
            az, elev = meca._mc_hd2ah(ha, dec, latitude)
        if self._mount_type.find("AZELEV")>=0:
            az = celb * self._d2r
            elev = celp * self._d2r
            ha, dec = meca._mc_ah2hd(az, elev, latitude)
        # --- rotator computation
        rotator = meca._mc_hd2parallactic(ha, dec, latitude)
        # --- conversion into degrees
        ha *= self._r2d
        dec *= self._r2d
        az *= self._r2d
        elev *= self._r2d
        rotator *= self._r2d
        # --- default units
        if unit_ha=="":
            unit_ha="H0.2"
        if unit_dec=="":
            unit_ha="d+090.1"
        if unit_az=="":
            unit_az="D0.2"
        if unit_elev=="":
            unit_elev="d+090.1"
        # --- conversion into sexagesimal if needed
        if unit_ha!="deg":
            ha = celme.Angle(ha).sexagesimal(unit_ha)
        if unit_dec!="deg":
            dec = celme.Angle(dec).sexagesimal(unit_dec)
        if unit_az!="deg":
            az = celme.Angle(az).sexagesimal(unit_az)
        if unit_elev!="deg":
            elev = celme.Angle(elev).sexagesimal(unit_elev)
        return ha, dec, az, elev, rotator

    def cel2hadec(self, celb:float, celp:float, unit_ha:str="", unit_dec:str="")->tuple:
        # ---
        ha, dec, az, elev, rotator = self.cel2astro(celb, celp, unit_ha, unit_dec, "", "")
        return ha, dec

    def cel2azelev(self, celb:float, celp:float, unit_az:str="", unit_elev:str="")->tuple:
        # ---
        ha, dec, az, elev, rotator = self.cel2astro(celb, celp, "", "", unit_az, unit_elev)
        return az, elev, rotator

# =====================================================================
# =====================================================================
# Methods hadec for users
# =====================================================================
# =====================================================================
# Level 4
# =====================================================================

    def hadec_travel_compute(self, ha_target:celme.Angle, dec_target:celme.Angle, pierside_target:int=Mountaxis.PIERSIDE_AUTO)->tuple:
        # --- shortcuts
        axisb = self.axis[Mountaxis.BASE]
        axisp = self.axis[Mountaxis.POLAR]
        # === Read the current position
        incsimus = ["" for kaxis in range(Mountaxis.AXIS_MAX)]        
        increalb, rotrealb, celrealb, increalp, rotrealp, celrealp, piersidereal, incsimub, rotsimub, celsimup, incsimup, rotsimup, celsimup, piersidesimu = self.enc2cel(incsimus,self.OUTPUT_LONG, self.SAVE_ALL)
        if axisb.real==True:
            incb_start = increalb
        else:
            incb_start = incsimub
        if axisp.real==True:
            incp_start = increalp
            pierside_start = piersidereal
        else:
            incp_start = incsimup
            pierside_start = piersidesimu
        # === Read the target position
        if pierside_target==Mountaxis.PIERSIDE_AUTO:
            pierside_target = pierside_start
        celb, celp, celr, dcelb, dcelp, dcelr = self.astro2cel("HADEC", ha_target, dec_target, self._hadec_speeddrift_ha_deg_per_sec, self._hadec_speeddrift_dec_deg_per_sec)
        increalb, rotrealb, celrealb, increalp, rotrealp, celrealp, piersidereal, incsimub, rotsimub, celsimup, incsimup, rotsimup, celsimup, piersidesimu = self.cel2enc(celb, celp, pierside_target, self.OUTPUT_LONG, self.SAVE_NONE)
        if axisb.real==True:
            incb_target = increalb
        else:
            incb_target = incsimub
        if axisp.real==True:
            incp_target = increalp
        else:
            incp_target = incsimup
        # --- delta incs for the travel (signed incs)
        dincb = incb_target - incb_start
        dincp = incp_target - incp_start
        #print("dincb={} incb_target={} incb_start={} dincp={} incp_target={} incp_start={}".format(dincb, incb_target, incb_start, dincp, incp_target, incp_start))
        # --- slew velocity (positive inc/sec)
        inc_per_secb = axisb.slew_deg_per_sec * axisb.inc_per_deg
        if dincb<0:
            inc_per_secb *= -1
        inc_per_secp = axisp.slew_deg_per_sec * axisb.inc_per_deg
        if dincp<0:
            inc_per_secp *= -1
        #print("inc_per_secb={} inc_per_secp={}".format(inc_per_secb, inc_per_secp))
        # --- delays of slewing (sec)
        delayb = abs(dincb / inc_per_secb)
        delayp = abs(dincp / inc_per_secp)
        if delayb>delayp:
            delay = delayb
        else:
            delay = delayp
        return delayb, delayp
        #print("delayb={} delayp={} delay={}".format(delayb, delayp, delay))
        # --- fraction of a turn of 360 deg (no unit)
        fincb = abs(dincb / axisb._inc_per_sky_rev)
        fincp = abs(dincp / axisp._inc_per_sky_rev)
        if fincb>fincp:
            finc = fincb
        else:
            finc = fincp
        # --- Compute the number of positions during the travel
        ddeg = 1.0 ; # increment in degrees
        npos = math.ceil(finc*360.0/ddeg)
        if npos<3:
            npos = 3
        #print("fincb={} fincp={} finc={}".format(fincb, fincp, finc))
        # --- list of positions
        ts = np.linspace(0,delay,npos) ; # sec
        pincbs = ts * 0
        pincps = ts * 0
        for kpos in range(0,npos):
            t = ts[kpos] ; # sec
            pincbs[kpos] =  incb_start + inc_per_secb * t
            if dincb>=0 and pincbs[kpos]>incb_target:
                pincbs[kpos] = incb_target
            if dincb<0 and pincbs[kpos]<incb_target:
                pincbs[kpos] = incb_target
            pincps[kpos] =  incp_start + inc_per_secp * t
            if dincp>=0 and pincps[kpos]>incp_target:
                pincps[kpos] = incp_target
            if dincp<0 and pincps[kpos]<incp_target:
                pincps[kpos] = incp_target
        
        print("pincbs={} pincps={}".format(pincbs, pincps))        
        # --- utc time
        #jdnow = celme.Date("now").jd()
        #jds = jdnow + ts/86400.
        # --- compute the az,elev
        azims = ts * 0
        elevs = ts * 0
        elevmin = 90
        incposs = ["" for kaxis in range(Mountaxis.AXIS_MAX)]
        for kpos in range(0,npos):
            incposs[Mountaxis.BASE] = pincbs[kpos] ; # inc
            incposs[Mountaxis.POLAR] = pincbs[kpos] ; # inc
            celb, celp, pierside = self.enc2cel(incposs,self.OUTPUT_SHORT, self.SAVE_NONE)
            ha, dec, az, elev, rotator = self.cel2astro(celb, celp, unit_ha="deg", unit_dec="deg", unit_az="deg", unit_elev="deg")
            azims[kpos] = az
            elevs[kpos] = dec
        elevmin = np.amin(elevs)
        # --- return
        return elevmin, ts, elevs

    def _my_hadec_init(self, ha:celme.Angle, dec:celme.Angle, pierside:int="")->tuple:
        # --- abstract method. 
        # --- Please overload it according your language protocol
        err = self.NO_ERROR
        res = 0
        return err, res
    
    def hadec_init(self, ha:celme.Angle, dec:celme.Angle, pierside:int="")->tuple:
        # === Read the current position
        incsimus = ["" for kaxis in range(Mountaxis.AXIS_MAX)] 
        increalb, rotrealb, celrealb, increalp, rotrealp, celrealp, piersidereal, incsimub, rotsimub, celsimub, incsimup, rotsimup, celsimup, piersidesimu = self.enc2cel(incsimus, self.OUTPUT_LONG, self.SAVE_ALL)
        # === Assign the side at the current position if not notified
        if pierside=="":
            pierside = piersidereal
        # === Target celestial position
        celb, celp, celr = self.hadec2cel(ha, dec)
        # --- update inc0 on each axis
        for kaxis in range(Mountaxis.AXIS_MAX):
            current_axis = self.axis[kaxis]
            if current_axis == None:
                continue            
            if current_axis.real==True:
                if kaxis == Mountaxis.BASE:
                    inc = increalb
                    cel = celb
                elif kaxis == Mountaxis.POLAR:
                    inc = increalp
                    cel = celp
            else:
                if kaxis == Mountaxis.BASE:
                    inc = incsimub
                    cel = celb
                elif kaxis == Mountaxis.POLAR:
                    inc = incsimup
                    cel = celp
            print("kaxis={} inc={} cel={}".format(kaxis,inc,cel))
            current_axis.update_inc0(inc, cel, pierside)
        # --- Real inits
        err,res = self._my_hadec_init(ha, dec, pierside)
        # --- return the current new position
        return self.hadec_coord()

    def hadec_coord(self, **kwargs:dict)->tuple:
        # --- Dicos of optional and mandatory parameters
        params_optional = {} 
        params_optional["UNIT_HA"] = (str,'H0.2')
        params_optional["UNIT_DEC"] = (str,'d+090.1')
        params_mandatory = {} 
        # --- Decode parameters
        params = self.decode_kwargs(params_optional, params_mandatory, **kwargs)
        # --- Get the simu and real coordinates of all axis
        incsimus = ["" for kaxis in range(Mountaxis.AXIS_MAX)]        
        celb, celp, pierside = self.enc2cel(incsimus,self.OUTPUT_SHORT, self.SAVE_ALL)
        ha, dec = self.cel2hadec(celb, celp, params["UNIT_HA"], params["UNIT_DEC"])
        return ha, dec, pierside
    
    def _my_hadec_speeddrift(self, deg_per_sec_ha, deg_per_sec_dec):
        # --- abstract method. 
        # --- Please overload it according your language protocol
        return deg_per_sec_ha, deg_per_sec_dec

    def hadec_speeddrift(self, deg_per_sec_ha="", deg_per_sec_dec=""):
        if deg_per_sec_ha!="" and deg_per_sec_dec!="":
            if deg_per_sec_ha=="diurnal":
                deg_per_sec_ha = 360./self._sideral_sec_per_day
            if deg_per_sec_dec=="diurnal":
                deg_per_sec_dec = 0.0
            #print("deg_per_sec_ha={} deg_per_sec_dec={}".format(deg_per_sec_ha,deg_per_sec_dec))
            deg_per_sec_ha, deg_per_sec_dec = self._my_hadec_speeddrift(deg_per_sec_ha, deg_per_sec_dec)
            self._hadec_speeddrift_ha_deg_per_sec = deg_per_sec_ha
            self._hadec_speeddrift_dec_deg_per_sec = deg_per_sec_dec
        return self._hadec_speeddrift_ha_deg_per_sec, self._hadec_speeddrift_dec_deg_per_sec

    def _my_hadec_goto(self, ha_target, dec_target, pierside_target):
        # --- abstract method. 
        # --- Please overload it according your language protocol
        err = self.NO_ERROR
        res = 0
        return err, res
        
    def hadec_goto(self, ha:celme.Angle, dec:celme.Angle, **kwargs):
        err = self.NO_ERROR
        res = 0
        # --- Dicos of optional and mandatory parameters
        params_optional = {} 
        params_optional["BLOCKING"] = (bool,False)
        params_optional["SIDE"] = (int,Mountaxis.PIERSIDE_AUTO)
        params_mandatory = {}
        # --- Decode parameters
        params = self.decode_kwargs(params_optional, params_mandatory, **kwargs)
        # === Read the current position
        incsimus = ["" for kaxis in range(Mountaxis.AXIS_MAX)]        
        celb, celp, pierside = self.enc2cel(incsimus,self.OUTPUT_SHORT, self.SAVE_ALL)
        ha_start, dec_start = self.cel2hadec(celb, celp, "deg", "deg")
        pierside_start = pierside
        # === Target celestial position
        # --- convert angles into deg [-180 ; 180]
        ha_target  = math.fmod(720+celme.Angle(ha).deg(), 360)
        if ha_target>180:
            ha_target -= 360
        dec_target = math.fmod(720+celme.Angle(dec).deg(), 360)
        if dec_target>180:
            dec_target -= 360
        # --- compute the target pierside
        lim_side_east = +30 ; # Tube west = PIERSIDE_POS1 = [-180 : lim_side_east]
        lim_side_west = -30 ; # Tube east = PIERSIDE_POS2 = [lim_side_west : +180]
        if params["SIDE"]==Mountaxis.PIERSIDE_AUTO:
            if ha_target>lim_side_west and ha_target<lim_side_east:
                # --- the target position is in the both possibilitiy range
                pierside_target = pierside_start
            else:
                if ha_target>lim_side_east:
                    # --- the target is after the limit of side=PIERSIDE_POS1
                    pierside_target = Mountaxis.PIERSIDE_POS2
                else:
                    # --- the target is before the limit of side=PIERSIDE_POS2
                    pierside_target = Mountaxis.PIERSIDE_POS1
        else:
            pierside_target = params["SIDE"]
        #print("ha_start={:.4f} ha_target={:.4f} pierside_start={}  params[\"SIDE\"]={} pierside_target={}".format(ha_start, ha_target, pierside_start, params["SIDE"], pierside_target))
        # --- Compute incs of the target and the drifts (for any mount_type)
        celb, celp, celr, dcelb, dcelp, dcelr = self.astro2cel("HADEC", ha_target, dec_target, self._hadec_speeddrift_ha_deg_per_sec, self._hadec_speeddrift_dec_deg_per_sec)
        incb, incp = self.cel2enc(celb, celp, pierside_target, self.OUTPUT_SHORT, self.SAVE_NONE)
        increalb, rotrealb, celrealb, increalp, rotrealp, celrealp, piersidereal, incsimub, rotsimub, celsimup, incsimup, rotsimup, celsimup, piersidesimu = self.cel2enc(celb, celp, pierside_target, self.OUTPUT_LONG, self.SAVE_NONE)
        #print("celrealb={:.4f} rotrealb={:.4f} pierside_target={} -> increalb={} incb={} incp={}".format(celrealb, rotrealb, piersidereal, increalb, incb, incp))
        incr = 0
        # ---
        for kaxis in range(Mountaxis.AXIS_MAX):
            current_axis = self.axis[kaxis]
            if current_axis == None:
                continue
            # === Target position and drift
            if kaxis == Mountaxis.BASE:
                inc = incb
                dcel = dcelb
            elif kaxis == Mountaxis.POLAR:
                inc = incp
                dcel = dcelp
            elif kaxis == Mountaxis.ROTATOR:
                inc = incr
                dcel = dcelr
            dslw = current_axis.slew_deg_per_sec
            # === Slew Velocity inc/sec
            inc_per_sec_slew  = dslw * current_axis.senseinc * current_axis.inc_per_deg
            # === Drift Velocity inc/sec
            inc_per_sec_drift = dcel * current_axis.senseinc * current_axis.inc_per_deg
            if self.site.latitude>=0:
                inc_per_sec_drift *= -1
            # === Simulation. It runs even if there is a real hardware)
            current_axis.simu_motion_start("ABSOLUTE", frame='inc', velocity=inc_per_sec_slew, position=inc, drift=inc_per_sec_drift)
        # === Real hardware
        err, res = self._my_hadec_goto(ha_target, dec_target, pierside_target)
        # === Wait the end of motion if needed
        if params["BLOCKING"]==True and dcelb==0 and dcelb==0 and dcelr==0:
            coord1 = self.hadec_coord(unit_ha="deg", unit_dec="deg")
            t0 = time.time()
            while True:
                time.sleep(0.1)
                dt = time.time()-t0
                if dt>60:
                    break
                coord2 = self.hadec_coord(unit_ha="deg", unit_dec="deg")
                if coord1 == coord2:
                    break
                coord1 = coord2
        # ---            
        return (err, res)

    def _my_hadec_move(self, ha_drift_deg_per_sec, dec_drift_deg_per_sec):
        # --- abstract method. 
        # --- Please overload it according your language protocol
        err = self.NO_ERROR
        res = 0
        return err, res

    def hadec_move(self, ha_drift_deg_per_sec, dec_drift_deg_per_sec, duration_s=0):        
        err = self.NO_ERROR
        res = 0
        # --- shortcuts
        axisb = self.axis[Mountaxis.BASE]
        axisp = self.axis[Mountaxis.POLAR]
        # === Simulation. It runs even if there is a real hardware)
        axisb.simu_motion_start("CONTINUOUS", frame='ang', drift=ha_drift_deg_per_sec)
        axisp.simu_motion_start("CONTINUOUS", frame='ang', drift=dec_drift_deg_per_sec)
        # === Real hardware
        err, res = self._my_hadec_move(ha_drift_deg_per_sec, dec_drift_deg_per_sec)
        return err, res

    def _my_hadec_move_stop(self):
        # --- abstract method. 
        # --- Please overload it according your language protocol
        err = self.NO_ERROR
        res = 0
        return err, res

    def hadec_move_stop(self):        
        err = self.NO_ERROR
        res = 0
        # --- shortcuts
        axisb = self.axis[Mountaxis.BASE]
        axisp = self.axis[Mountaxis.POLAR]
        # === Simulation. It runs even if there is a real hardware)
        ha_drift_deg_per_sec = self._hadec_speeddrift_ha_deg_per_sec
        dec_drift_deg_per_sec = self._hadec_speeddrift_dec_deg_per_sec
        axisb.simu_motion_start("CONTINUOUS", frame='ang', drift=ha_drift_deg_per_sec)
        axisp.simu_motion_start("CONTINUOUS", frame='ang', drift=dec_drift_deg_per_sec)
        # === Real hardware
        err, res = self._my_hadec_move_stop()
        return err, res

    def _my_hadec_stop(self):
        # --- abstract method. 
        # --- Please overload it according your language protocol
        err = self.NO_ERROR
        res = 0
        return err, res

    def hadec_stop(self):
        # === Real hardware
        err, res = self._my_hadec_stop()
        # === Read the current position
        incsimus = ["" for kaxis in range(Mountaxis.AXIS_MAX)]        
        celb, celp, pierside = self.enc2cel(incsimus,self.OUTPUT_SHORT, self.SAVE_ALL)
        # === Stop the simulation and update the simulated position by the real one if exists
        for kaxis in range(Mountaxis.AXIS_MAX):
            current_axis = self.axis[kaxis]
            if current_axis == None:
                continue
            current_axis.simu_motion_stop()
            if current_axis.real == True:
                current_axis.synchro_real2simu()
        return err, res

# =====================================================================
# =====================================================================
# Methods radec for users
# =====================================================================
# =====================================================================
# Level 4
# =====================================================================

    def _my_radec_speeddrift(self, deg_per_sec_ra, deg_per_sec_dec):
        # --- abstract method. 
        # --- Please overload it according your language protocol
        return deg_per_sec_ra, deg_per_sec_dec

    def radec_speeddrift(self, deg_per_sec_ra="", deg_per_sec_dec=""):
        if deg_per_sec_ra!="" and deg_per_sec_dec!="":
            if deg_per_sec_ra=="diurnal":
                deg_per_sec_ra = 0
            if deg_per_sec_dec=="diurnal":
                deg_per_sec_dec = 0.0
            deg_per_sec_ra, deg_per_sec_dec = self._my_hadec_speeddrift(deg_per_sec_ra, deg_per_sec_dec)
            self._drift_radec_ra_deg_per_sec = deg_per_sec_ra
            self._drift_radec_dec_deg_per_sec = deg_per_sec_dec
        return self._drift_radec_ra_deg_per_sec, self._drift_radec_dec_deg_per_sec

    def radec_app2cat(self, jd, ha, dec, params):
        longuai = self.site.longitude*self._d2r
        ha *= self._d2r
        dec *= self._d2r
        meca = celme.Mechanics()
        ra = meca._mc_hd2ad(jd, longuai, ha)
        equinox = celme.Date(params["EQUINOX"]).jd()
        # --- calcul de la precession ---*/
        ra_equinox, dec_equinox = meca._mc_precad(jd,ra,dec,equinox)
        #ra_equinox, dec_equinox = (ra,dec)
        ra_equinox *= self._r2d
        dec_equinox *= self._r2d
        return ra_equinox, dec_equinox
    
    def radec_coord(self, **kwargs):
        # --- Dicos of optional and mandatory parameters
        params_optional = {} 
        params_optional["UNIT_RA"] = (str,'H0.2')
        params_optional["UNIT_DEC"] = (str,'d+090.1')
        params_optional["EQUINOX"] = (str,'J2000')
        params_mandatory = {}
        # --- Decode parameters
        params = self.decode_kwargs(params_optional, params_mandatory, **kwargs)
        # --- Get HA,Dec
        ha, dec, pierside = self.hadec_coord(UNIT_HA="deg",UNIT_DEC="deg")
        jd = celme.Date("now").jd()
        # --- app2cat
        ra_equinox, dec_equinox = self.radec_app2cat(jd, ha, dec, params)
        # --- Output unit conversion
        if params["UNIT_RA"]!="deg":
            ra_equinox = celme.Angle(ra_equinox).sexagesimal(params["UNIT_RA"])
        if params["UNIT_DEC"]!="deg":
            dec_equinox = celme.Angle(dec_equinox).sexagesimal(params["UNIT_DEC"])       
        return ra_equinox, dec_equinox, pierside

    def radec_cat2app(self, jd, ra_equinox, dec_equinox, params):
        equinox = celme.Date(params["EQUINOX"]).jd()
        meca = celme.Mechanics()
        ra_equinox *= self._d2r
        dec_equinox *= self._d2r
        # --- calcul de la precession ---*/
        ra, dec = meca._mc_precad(equinox,ra_equinox,dec_equinox,jd)
        #print("ra({})={} dec({})={}".format(params["EQUINOX"], ra_equinox*self._r2d, params["EQUINOX"], dec_equinox*self._r2d))
        #print("ra(date)={} dec(date)={}".format(ra*self._r2d, dec*self._r2d))
        #ra, dec = (ra_equinox,dec_equinox)
        longuai = celme.Angle(self.site.longitude).rad()
        ha = meca._mc_ad2hd(jd, longuai, ra)
        ha *= self._r2d
        dec *= self._r2d
        return ha, dec

    def radec_init(self, ra_angle:celme.Angle, dec_angle:celme.Angle, pierside:int="", **kwargs)->tuple:
        # --- Dicos of optional and mandatory parameters
        params_optional = {} 
        params_optional["EQUINOX"] = (str,'J2000')
        params_optional["BLOCKING"] = (bool,False)
        params_mandatory = {} 
        # --- Decode parameters
        params = self.decode_kwargs(params_optional, params_mandatory, **kwargs)
        # --- Get HA,Dec
        ra_equinox = celme.Angle(ra_angle).deg()
        dec_equinox = celme.Angle(dec_angle).deg()
        jd = celme.Date("now").jd()
        # --- cat2app
        ha_target, dec_target = self.radec_cat2app(jd, ra_equinox, dec_equinox,params)
        self.hadec_init(ha_target, dec_target, pierside)
        return self.radec_coord()
        
    def radec_goto(self, ra_angle:celme.Angle, dec_angle:celme.Angle, **kwargs):
        # --- Dicos of optional and mandatory parameters
        params_optional = {} 
        params_optional["EQUINOX"] = (str,'J2000')
        params_optional["BLOCKING"] = (bool,False)
        params_optional["SIDE"] = (int,Mountaxis.PIERSIDE_AUTO)
        params_mandatory = {} 
        # --- Decode parameters
        params = self.decode_kwargs(params_optional, params_mandatory, **kwargs)
        pierside_target = params["SIDE"]
        # --- Get HA,Dec
        ra_equinox = celme.Angle(ra_angle).deg()
        dec_equinox = celme.Angle(dec_angle).deg()
        jd = celme.Date("now").jd()
        # --- cat2app
        ha_target, dec_target = self.radec_cat2app(jd, ra_equinox, dec_equinox, params)
        # --- compute the time of slewing
        delayb, delayp = self.hadec_travel_compute(ha_target, dec_target, pierside_target)
        # --- Temporal offset to anticipate the delay of slewing
        dt = 2.4 + delayb ; # (sec)
        print("delayb={} dt={}".format(delayb,dt))
        jd += dt/86400
        # --- cat2app
        ha_target, dec_target = self.radec_cat2app(jd, ra_equinox, dec_equinox, params)
        # --- store the current HADEC drift values
        drift_hadec_ha_deg_per_sec = self._hadec_speeddrift_ha_deg_per_sec
        drift_hadec_dec_deg_per_sec = self._hadec_speeddrift_ha_deg_per_sec
        # --- change the HADEC drift values by the RADEC ones
        self._hadec_speeddrift_ha_deg_per_sec  = 360./self._sideral_sec_per_day - self._drift_radec_ra_deg_per_sec
        self._hadec_speeddrift_dec_deg_per_sec = self._drift_radec_dec_deg_per_sec
        # --- cal the HADEC GOTO
        #print("=== _drift_hadec_ha_deg_per_sec={}".format(self._drift_hadec_ha_deg_per_sec))
        err, res = self.hadec_goto(ha_target, dec_target, **kwargs)
        # --- restore the current HADEC drift values
        self._hadec_speeddrift_ha_deg_per_sec = drift_hadec_ha_deg_per_sec
        self._hadec_speeddrift_dec_deg_per_sec = drift_hadec_dec_deg_per_sec        
        # ---
        return (err, res)

# =====================================================================
# =====================================================================
# Methods inc for users
# =====================================================================
# =====================================================================

    def _my_inc_goto(self, axis_id:int, inc:float, inc_per_sec_slew:float):
        # --- abstract method. 
        # --- Please overload it according your language protocol
        err = self.NO_ERROR
        res = 0
        return err, res
        
    def inc_goto(self, axis_id:int, inc:float, inc_per_sec_slew:float):
        # === Simulation. It runs even if there is a real hardware)
        err = self.NO_ERROR
        res = 0
        kaxis = axis_id
        current_axis = self.axis[kaxis]
        current_axis.simu_motion_start("ABSOLUTE", frame='inc', velocity=inc_per_sec_slew, position=inc, drift=0)
        # === Real
        err, res = self._my_inc_goto(axis_id, inc, inc_per_sec_slew)
        return err, res

    def _my_goto_park(self, incb:float, incp:float, incr:float):
        # --- abstract method. 
        # --- Please overload it according your language protocol
        err = self.NO_ERROR
        res = 0
        return err, res

    def goto_park(self, incb="", incp="", incr=""):
        err, res = self._my_goto_park(incb, incp, incr)
        return err, res
        
    def plot_rot(self, lati, azim, elev, rotb, rotp, outfile=""):
        """
        Vizualize the rotation angles of the mount according the local coordinates
        # --- Site latitude
        lati (deg) : Site latitude
        # --- Observer view
        elev = 15 # turn around the X axis
        azim = 140 # turn around the Z axis (azim=0 means W foreground, azim=90 means N foreground) 
        # --- rob, rotp
        """
        if lati=="":
            lati = self.site.latitude
        if lati>=0:
            latitude = lati
            cards="SNEW"
        else:
            latitude = -lati
            cards="NSWE"
        toplots = []
        # --- cartesian frame
        options = {"linewisth":0.5}
        toplots.append(["text",[0, 0, 0],"o",'k',{}])
        toplots.append(["line",[0, 0, 0],[1, 0, 0],'k',options])
        toplots.append(["line",[0, 0, 0],[-1, 0, 0],'k',options])
        toplots.append(["text",[1, 0, 0],cards[0],'k',{}])
        toplots.append(["text",[-1, 0, 0],cards[1],'k',{}])
        toplots.append(["text",[0, 1, 0],cards[2],'k',{}])
        toplots.append(["text",[0, -1, 0],cards[3],'k',{}])
        toplots.append(["line",[0, 0, 0],[0, 1, 0],'k',options])
        toplots.append(["line",[0, 0, 0],[0, -1, 0],'k',options])
        toplots.append(["line",[0, 0, 0],[0, 0, 1],'k',options])
        toplots.append(["line",[0, 0, 0],[0, 0, -1],'k',options])
        toplots.append(["text",[0, 0, 1.1],"zenith",'k',options])
        # --- great circle tangeant to the projection plane
        toplots.append(["circle",1,[[90,0,0], [-elev,0,0], [0,0,-azim]],0,360,'k',options])
        #toplots.append(["circle",1,[[90,-elev,-azim]],0,360,'k',{"linewisth":1}])
        # --- local horizon
        toplots.append(["circle",1,[[0,0,0]],0,360,'k',{"linewisth":0.5}])
        # --- local meridian
        toplots.append(["circle",1,[[90,0,0]],0,360,'k',{"linewisth":0.5}])
        # --- local equator
        rotxyzs = [[0, latitude, 0]]
        options = {"linewisth":2}
        color = 'r'
        toplots.append(["circle",1,rotxyzs,0,360,color,options])
        toplots.append(["linefrom0",1,rotxyzs,color,options])
        toplots.append(["linefrom0",-1,rotxyzs,color,options])
        toplots.append(["textfrom0",1.05,rotxyzs,"x",color,options])
        rotxyzs = [[0, 0, 90]]
        toplots.append(["linefrom0",1,rotxyzs,color,options])
        toplots.append(["linefrom0",-1,rotxyzs,color,options])
        toplots.append(["textfrom0",1.15,rotxyzs,"y",color,options])
        rotxyzs = [[0, latitude+90, 0]]
        toplots.append(["linefrom0",1,rotxyzs,color,options])
        toplots.append(["linefrom0",-1,rotxyzs,color,options])
        toplots.append(["textfrom0",1.05,rotxyzs,"z = visible pole ({})".format(cards[1]),color,options])
        # --- pointing Rotp
        rotxyzs = [[90,0,0], [0,0,rotb], [0, latitude, 0]]
        options = {"linewisth":0.5}
        color = 'r'
        toplots.append(["circle",1,rotxyzs,0,360,color,options])
        options = {"linewisth":2}
        color = 'g'
        toplots.append(["circle",1,rotxyzs,90,90-rotp,color,options])
        rotxyzs = [[0,90-rotp,0], [0,0,rotb], [0, latitude, 0]]
        options = {"linewisth":1}
        toplots.append(["linefrom0",1,rotxyzs,color,options])
        toplots.append(["textfrom0",1.1,rotxyzs,"rotp",color,options])
        rotxyzs = [[0,90,0], [0,0,rotb], [0, latitude, 0]]
        toplots.append(["linefrom0",1,rotxyzs,color,options])
        # --- pointing Rotb
        rotxyzs = [[0,0,rotb], [0, latitude, 0]]
        options = {"linewisth":2}
        color = 'b'
        toplots.append(["circle",1,rotxyzs,0,-rotb,color,options])
        options = {"linewisth":1}
        toplots.append(["linefrom0",1,rotxyzs,color,options])
        toplots.append(["textfrom0",1.3,rotxyzs,"rotb",color,options])
        rotxyzs = [[0, latitude, 0]]
        options = {"linewisth":1}
        toplots.append(["linefrom0",1,rotxyzs,color,options])
        # ===
        fig = plt.figure()
        #ax = fig.add_subplot(1,1,1,projection='3d')
        #ax.view_init(elev=elev, azim=azim)
        ax = fig.add_subplot(1,1,1)
        for toplot in toplots:
            ptype = toplot[0]
            if ptype=="line":
                dummy, xyz1, xyz2, color, options = toplot
                xyz0s = []
                xyz0 = np.array(xyz1)
                xyz0s.append(xyz0)
                xyz0 = np.array(xyz2)
                xyz0s.append(xyz0)
                na = len(xyz0s)
            elif ptype=="text":
                dummy, xyz, text, color, options = toplot
                xyz0s = []
                xyz0 = np.array(xyz)
                xyz0s.append(xyz0)
                na = len(xyz0s)
            elif ptype=="linefrom0":
                dummy, length, rotxyzs, color, options = toplot
                xyz0s = []
                xyz0 = np.array([0, 0, 0])
                xyz0s.append(xyz0)
                xyz0 = np.array([length, 0 ,0])
                xyz0s.append(xyz0)
                na = len(xyz0s)
            elif ptype=="textfrom0":
                dummy, length, rotxyzs, text, color, options = toplot
                xyz0s = []
                xyz0 = np.array([length, 0 ,0])
                xyz0s.append(xyz0)
                na = len(xyz0s)
            elif ptype=="circle":
                dummy, radius, rotxyzs, ang1, ang2, color, options = toplot
                na = 50
                alphas =np.linspace(ang1,ang2,na)
                # --- cercle dans le plan (x,y)
                xyz0s = []
                for alpha in alphas:
                    alpha = math.radians(alpha)
                    x = radius*math.cos(alpha)
                    y = radius*math.sin(alpha)
                    z = 0
                    xyz0 = np.array([x, y, z])
                    xyz0s.append(xyz0)
                na = len(xyz0s)
            else:
                continue
            # --- rotations
            if ptype=="linefrom0" or ptype=="textfrom0" or ptype=="circle":
                for rotxyz in rotxyzs:
                    rotx, roty, rotz = rotxyz
                    cosrx = math.cos(math.radians(rotx))
                    sinrx = math.sin(math.radians(rotx))
                    cosry = math.cos(math.radians(roty))
                    sinry = math.sin(math.radians(roty))
                    cosrz = math.cos(math.radians(rotz))
                    sinrz = math.sin(math.radians(rotz))
                    rotrx = np.array([ [1, 0, 0], [0, cosrx, -sinrx], [0, sinrx, cosrx] ])
                    rotry = np.array([ [cosry, 0, -sinry], [0, 1, 0], [sinry, 0, cosry] ])
                    rotrz = np.array([ [cosrz, -sinrz, 0], [sinrz, cosrz, 0] , [0, 0, 1] ])
                    xyz1s = []
                    for xyz in xyz0s:
                        xyz = np.dot(rotrx, xyz)
                        xyz = np.dot(rotry, xyz)
                        xyz = np.dot(rotrz, xyz)
                        xyz1s.append(xyz)
                    xyz0s = xyz1s # ready for a second rotation
            else:
                xyz1s = xyz0s
            # --- projections
            cosaz = math.cos(math.radians(azim))
            sinaz = math.sin(math.radians(azim))
            cosel = math.cos(math.radians(elev))
            sinel = math.sin(math.radians(elev))
            rotaz = np.array([ [cosaz, -sinaz, 0], [sinaz, cosaz, 0] , [0, 0, 1] ])
            #rotel = np.array([ [cosel, 0, -sinel], [0, 1, 0], [sinel, 0, cosel] ])
            rotel = np.array([ [1, 0, 0], [0, cosel, -sinel], [0, sinel, cosel] ])
            xyz2s = []
            for xyz in xyz1s:
                xyz = np.dot(rotaz, xyz)
                xyz = np.dot(rotel, xyz)
                xyz2s.append(xyz)
            # --- plot
            if ptype=="textfrom0" or ptype=="text":
                x,y,z = xyz2s[0]
                h = ax.text(x,z,text)
                h.set_color(color)
            else:
                for ka in range(0,na-1):
                    xyz1 = xyz2s[ka]
                    xyz2 = xyz2s[ka+1]
                    x = [xyz1[0], xyz2[0]]
                    y = [xyz1[1], xyz2[1]]
                    z = [xyz1[2], xyz2[2]]
                    if y[0]>1e-3 or y[1]>1e-3:
                        symbol = color+':'
                    else:
                        symbol = color+'-'
                    h = ax.plot(x,z,symbol,'linewidth',1.0)
                    for option in options.items():
                        key = option[0]
                        val = option[1]
                        if key=="linewisth":
                            h[0].set_linewidth(val)
        # ---
        ax.axis('equal')
        #fig.patch.set_visible(False)
        ax.axis('off')        
        plt.title("Rotation angles for latitude {} deg".format(lati))
        plt.show()
        if outfile!="":
            plt.savefig(outfile, facecolor='w', edgecolor='w')
        
    def pad(self, pad_type="dev"):
        self.pad = Mountpad(self,pad_type)
        self.pad.start()

    def remote_command_protocol(self, remote_command_protocol="LX200"):
        remote_command_protocol = remote_command_protocol.upper()
        if remote_command_protocol=="LX200":
            self._remote_command_protocol="LX200"
            self._lx200_precision = self.LX200_HIGH_PRECISION
            self._lx200_setra_deg = 0.0
            self._lx200_setdec_deg = 0.0
            self._lx200_24_12_hour_time_format = 24

    def _my_remote_command_processing(self, command):
        # --- abstract method. 
        # --- Please overload it according your language protocol
        res = ""
        return res

    def remote_command_processing(self, command):
        res = ""
        if self._remote_command_protocol=="LX200":
            if command.startswith(chr(6))==True:
                res = "P"
            elif command.startswith(":Gg")==True:
                deg = self.site.longitude
                res = celme.Angle(deg).sexagesimal("d +0180")
                res = res[0:4] + "*"  + res[5:7]
            elif command.startswith(":Gm")==True:
                ha, dec, side = self.hadec_coord()
                if side==1:
                    res = "W"
                else:
                    res = "E"
            elif command.startswith(":Gt")==True:
                deg = self.site.latitude
                res = celme.Angle(deg).sexagesimal("d +090")
                res = res[0:3] + "*" + res[4:6]
            elif command.startswith(":GD")==True:
                ra, dec, side = self.radec_coord(UNIT_RA="H:0.0",UNIT_DEC="d:+090.0")
                dec = dec[0:3] + "*" + dec[4:6] + "'" + dec[7:9]
                if self._lx200_precision == self.LX200_LOW_PRECISION:
                    dec = dec[0:6]
                res = dec
            elif command.startswith(":GR")==True:
                ra, dec, side = self.radec_coord(UNIT_RA="H:0.0",UNIT_DEC="d:+090.0")
                if self._lx200_precision == self.LX200_LOW_PRECISION:
                    sec = int(ra[6:8])
                    sec = "{:.1f}".format(sec/60.0)
                    ra = ra[0:5]+sec[1:]
                res = ra
            elif command.startswith(":H")==True:                
                if self._lx200_24_12_hour_time_format == 24:
                    self._lx200_24_12_hour_time_format = 12
                else:
                    self._lx200_24_12_hour_time_format = 24
            elif command.startswith(":Gc")==True:
                res = str(self._lx200_24_12_hour_time_format)
            elif command.startswith(":GC")==True:
                iso = celme.Date("now").iso()
                res = iso[5:7] + "/" + iso[8:10] + "/" + iso[2:4]                
            elif command.startswith(":GL")==True:
                iso = celme.Date("now").iso()
                res = iso[11:19]
            elif command.startswith(":Q")==True:
                self.hadec_stop()
            elif command.startswith(":P")==True:
                res = self._lx200_precision
            elif command.startswith(":U")==True:
                if self._lx200_precision == self.LX200_HIGH_PRECISION:
                    self._lx200_precision = self.LX200_LOW_PRECISION
                else:
                    self._lx200_precision = self.LX200_HIGH_PRECISION
            elif command.startswith(":Sr")==True:
                angle = command[3:]
                self._lx200_setra_deg = celme.Angle(angle).deg() * 15.0
                res = str(1)
            elif command.startswith(":Sd")==True:
                angle = command[3:]
                self._lx200_setdec_deg = celme.Angle(angle).deg()
                res = str(1)
            elif command.startswith(":Sg")==True:
                angle = command[3:]
                self.site.longitude = celme.Angle(angle).deg()
                res = str(1)
            elif command.startswith(":Sts")==True:
                angle = command[3:]
                self.site.latitude = celme.Angle(angle).deg()
                res = str(1)                
            elif command.startswith(":CM")==True:
                self.radec_init(self._lx200_setra_deg, self._lx200_setdec_deg)
            elif command.startswith(":MS")==True:
                self.radec_goto(self._lx200_setra_deg, self._lx200_setdec_deg)
        if self._remote_command_protocol=="astromecca":
            if command.startswith(chr(6))==True:
                res = "P"
            elif command.startswith(":Gg")==True:
                deg = self.site.longitude
                res = celme.Angle(deg).sexagesimal("d +0180")
                res = res[0:4] + "*"  + res[5:7]       
                """
        while True:
            # --- read the ASCOM commands
            lignes = ""
            try:
                err, lignes = ascom_chan.read_chan()
                #print("lignes = {}".format(lignes))
                lignes = str(lignes[0])
            except:
                pass
            if err==ascom_chan.NO_ERROR and lignes != "":
                mount_scx11.log.print("===\nlignes recu ={}".format(lignes))
                mots = lignes.split()
                mount_scx11.log.print("mots recus ={}".format(mots))
                msg = ""
                # --- traiter la ligne recue de ASCOM
                if lignes.startswith("READ_RA")==True:
                    # --- on traite un goto
                    ra, dec, side = mount_scx11.radec_coord(UNIT_RA="deg")
                    mount_scx11.log.print("{} ra={}".format(mots[0],ra))
                    msg = str( float(ra)/15.)
                if lignes.startswith("READ_DEC")==True:
                    # --- on traite un goto
                    ra, dec, side = mount_scx11.radec_coord(UNIT_DEC="deg")
                    mount_scx11.log.print("{} dec={}".format(mots[0],dec))
                    msg = str(dec)
                # --- traiter la ligne recue de ASCOM
                if lignes.startswith("GOTO_RA_DEC")==True:
                    ra = float(mots[1])*15.
                    dec = float(mots[2])
                    mount_scx11.log.print("{} ra={} dec={}".format(mots[0],ra,dec))
                    mount_scx11.radec_goto(ra, dec, EQUINOX="NOW")
                if lignes.startswith("PARK")==True:
                    ha = 270.0
                    dec = 90.0
                    side = 1
                    mount_scx11.log.print("{} ha={} dec={} side={}".format(mots[0],ra,dec,side))
                    mount_scx11.hadec_goto(ha, dec, side=side)
                # --- send the EQMOD answer to ASCOM
                mount_scx11.log.print("msg renvoyé ={}".format(msg))
                ascom_chan.put_chan(msg)
                 """
        else:
            pass
        if res=="":
            res = self._my_remote_command_processing(command)
        return res
    
# =====================================================================
# =====================================================================
# Special methods
# =====================================================================
# =====================================================================
        
    def __init__(self, *args, **kwargs): 
        """
        Conversion from Uniform Python object into protocol language
        Usage : Mountastro("HADEC", name="SCX11")
        """
        # --- Dico of optional parameters for all axis_types
        param_optionals = {}
        param_optionals["MODEL"] = (str, "")
        param_optionals["MANUFACTURER"] = (str, "")
        param_optionals["SERIAL_NUMBER"] = (str, "")
        param_optionals["REAL"] = (bool, False)
        param_optionals["DESCRIPTION"] = (str, "No description.")
        param_optionals["SITE"] = (celme.Site,"GPS 0 E 45 100")
        # --- Dico of axis_types and their parameters
        mount_types = {}
        mount_types["HADEC"]= {"MANDATORY" : {"NAME":[str,"Unknown"]}, "OPTIONAL" : {"LABEL":[str,"Equatorial"]} }
        mount_types["HADECROT"]= {"MANDATORY" : {"NAME":[str,"Unknown"]}, "OPTIONAL" : {"LABEL":[str,"Equatorial"]} }
        mount_types["AZELEV"]= {"MANDATORY" : {"NAME":[str,"Unknown"]}, "OPTIONAL" : {"LABEL":[str,"Altazimutal"]} }
        mount_types["AZELEVROT"]= {"MANDATORY" : {"NAME":[str,"Unknown"]}, "OPTIONAL" : {"LABEL":[str,"Altazimutal"]} }
        # --- Decode args and kwargs parameters
        self._mount_params = self.decode_args_kwargs(0, mount_types, param_optionals, *args, **kwargs)
        # ===
        self._mount_type = self._mount_params["SELECTED_ARG"]
        self._name = self._mount_params["NAME"]
        self._description = self._mount_params["DESCRIPTION"]
        self._model = self._mount_params["MODEL"]
        self._manufacturer = self._mount_params["MANUFACTURER"]
        self._serial_number= self._mount_params["SERIAL_NUMBER"]
        # === local
        if type(self._mount_params["SITE"])==celme.site.Site:
            self.site = self._mount_params["SITE"]
        else:
            self.site = celme.Home(self._mount_params["SITE"])
        # === Initial state real or simulation
        real = self._mount_params["REAL"]
        # === Initialisation of axis list according the mount_type
        self.axis = [None for kaxis in range(Mountaxis.AXIS_MAX)] # or  Mountaxis.AXIS_NOT_DEFINED            
        if self._mount_type.find("HA")>=0:
            current_axis = Mountaxis("HA", name = "Hour angle")
            current_axis.update_inc0(0,0,current_axis.PIERSIDE_POS1)
            self.axis[Mountaxis.BASE] = current_axis
        if self._mount_type.find("DEC")>=0: 
            current_axis = Mountaxis("DEC", name = "Declination")
            current_axis.update_inc0(0,90,current_axis.PIERSIDE_POS1)
            self.axis[Mountaxis.POLAR] = current_axis
        if self._mount_type.find("AZ")>=0:
            current_axis = Mountaxis("AZ", name = "Azimuth")
            current_axis.update_inc0(0,0,current_axis.PIERSIDE_POS1)
            self.axis[Mountaxis.BASE] = current_axis
        if self._mount_type.find("ELEV")>=0: 
            current_axis = Mountaxis("ELEV", name = "Elevation")
            current_axis.update_inc0(0,90,current_axis.PIERSIDE_POS1)
            self.axis[Mountaxis.POLAR] = current_axis
        if self._mount_type.find("ROT")>=0: 
            current_axis = Mountaxis("ROT", name = "Rotator")
            current_axis.update_inc0(0,0,current_axis.PIERSIDE_POS1)
            self.axis[Mountaxis.ROT] = current_axis
        if self._mount_type.find("ROLL")>=0:
            current_axis = Mountaxis("ROLL", name = "Roll")
            current_axis.update_inc0(0,0,current_axis.PIERSIDE_POS1)
        if self._mount_type.find("PITCH")>=0: 
            current_axis = Mountaxis("PITCH", name = "Pitch")
            current_axis.update_inc0(0,90,current_axis.PIERSIDE_POS1)
            self.axis[Mountaxis.POLAR] = current_axis
        if self._mount_type.find("YAW")>=0: 
            current_axis = Mountaxis("YAW", name = "Yaw")
            current_axis.update_inc0(0,0,current_axis.PIERSIDE_POS1)
            self.axis[Mountaxis.YAW] = current_axis
        # === Default Setup
        ratio_wheel_puley = 1 ; # 5.25
        ratio_puley_motor = 100.0 ; # harmonic reducer
        inc_per_motor_rev = 1000.0 ; # IMC parameter. System Confg -> System Parameters - Distance/Revolution
        for kaxis in range(Mountaxis.AXIS_MAX):
            current_axis = self.axis[kaxis]
            if current_axis == None:
                continue
            current_axis.slewmax_deg_per_sec = 30    
            current_axis.slew_deg_per_sec = 30    
            current_axis.real = real
            current_axis.latitude = self.site.latitude
            current_axis.ratio_wheel_puley = ratio_wheel_puley    
            current_axis.ratio_puley_motor = ratio_puley_motor    
            current_axis.inc_per_motor_rev = inc_per_motor_rev
        # ---
        self.hadec_speeddrift(0,0)
        # === Pads
        self._pads = []
        # === Remote commands        
        self._remote_command_protocol=""
        self.remote_command_protocol("LX200")
        # === Log positions
        self._record_positions = True        
        # === Log
        path_data = os.getcwd()
        self.log = Mountlog(self._name,self.site.gps,path_data)
        self.log.print("Launch log")

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

if __name__ == "__main__":
    cwd = os.getcwd()

    example = 4
    print("Example       = {}".format(example))
            
    if example == 1:
        # === SCX11
        home = celme.Home("GPS 2.25 E 43.567 148")
        site = celme.Site(home)
        mount = Mountastro("HADEC", name="Test Mount",site=site)

    if example == 2:
        home = celme.Home("GPS 2.25 E 43.567 148")
        #horizon = [(0,0), (360,0)]
        site = celme.Site(home)
        mount_scx11 = Mountastro("HADEC", name="Guitalens Mount", manufacturer="Astro MECCA", model="TM350", site=site)
        mount_eqmod = Mountastro("HADEC", name="Guitalens Mount", manufacturer="Astro MECCA", model="EQ 6",  site=site)
        # --- shortcuts
        mount_scx11_axisb = mount_scx11.axis[Mountaxis.BASE]
        mount_scx11_axisp = mount_scx11.axis[Mountaxis.POLAR]
        mount_eqmod_axisb = mount_eqmod.axis[Mountaxis.BASE]
        mount_eqmod_axisp = mount_eqmod.axis[Mountaxis.POLAR]
        # --- simulation or not
        mount_scx11_axisb.real = False
        mount_scx11_axisp.real = False
        mount_eqmod_axisb.real = False
        mount_eqmod_axisp.real = False
        # ======= SCX11
        incsimus = [0 for kaxis in range(Mountaxis.AXIS_MAX)]
        increals, incsimus = mount_scx11.read_encs(incsimus)
        # ---
        incb, rotb, celb, incp, rotp, celp, pierside, incsimub, rotsimub, celsimub, incsimup, rotsimup, celsimup, piersidesimu = mount_scx11.enc2cel(incsimus, mount_scx11.OUTPUT_LONG, mount_scx11.SAVE_ALL)        
        print("SCX11 incb={:.0f} rotb={:.3f} celb={:.3f} incp={:.0f} rotp={:.3f} celp={:.3f} pierside={:.0f}".format(incsimub, rotsimub, celsimub, incsimup, rotsimup, celsimup, piersidesimu))
        incb, rotb, celb, incp, rotp, celp, pierside, incsimub, rotsimub, celsimub, incsimup, rotsimup, celsimup, piersidesimu = mount_scx11.enc2cel(incsimus, mount_scx11.OUTPUT_LONG, mount_scx11.SAVE_ALL)        
        print("SCX11 incb={:.0f} rotb={:.3f} celb={:.3f} incp={:.0f} rotp={:.3f} celp={:.3f} pierside={:.0f}".format(incsimub, rotsimub, celsimub, incsimup, rotsimup, celsimup, piersidesimu))
        # --- update EQMOD parameters
        incb, rotb, celb, incp, rotp, celp, pierside, incsimub, rotsimub, celsimub, incsimup, rotsimup, celsimup, piersidesimu = mount_eqmod.enc2cel(incsimus, mount_scx11.OUTPUT_LONG, mount_scx11.SAVE_ALL)        
        print("EQMOD incb={:.0f} rotb={:.3f} celb={:.3f} incp={:.0f} rotp={:.3f} celp={:.3f} pierside={:.0f}".format(incsimub, rotsimub, celsimub, incsimup, rotsimup, celsimup, piersidesimu))
        incb, rotb, celb, incp, rotp, celp, pierside, incsimub, rotsimub, celsimub, incsimup, rotsimup, celsimup, piersidesimu = mount_eqmod.enc2cel(incsimus, mount_scx11.OUTPUT_LONG, mount_scx11.SAVE_ALL)        
        print("EQMOD incb={:.0f} rotb={:.3f} celb={:.3f} incp={:.0f} rotp={:.3f} celp={:.3f} pierside={:.0f}".format(incsimub, rotsimub, celsimub, incsimup, rotsimup, celsimup, piersidesimu))
               
    if example == 3:
        if True:
            # --- Site latitude
            latitude = 30
            # --- observer view
            elev = 15 # tourne autour de l'axe x
            azim = 140 # tourne autour de de l'axe z (azim=0 on regarde W devant, azim=90 N devant) 
            outfile = cwd+"/rotbp_n.png"
        else:
            # --- Site latitude
            latitude = -30
            # --- observer view
            elev = 15 # tourne autour de l'axe x
            azim = 140 # tourne autour de de l'axe z (azim=0 on regarde W devant, azim=90 N devant) 
            outfile = cwd+"/rotbp_s.png"
        # --- rob, rotp
        rotb = 60
        rotp = 30
        home = celme.Home("GPS 2.25 E 43.567 148")
        site = celme.Site(home)
        mount = Mountastro("HADEC", name="Example Mount", site=site)
        mount.plot_rot(latitude, azim, elev, rotb, rotp, outfile)

    if example==4:
        home = celme.Home("GPS 2.25 E 43.567 148")
        #horizon = [(0,0), (360,0)]
        site = celme.Site(home)
        mount_scx11 = Mountastro("HADEC", name="Guitalens Mount", manufacturer="Astro MECCA", model="TM350", site=site)
        # --- shortcuts
        mount_scx11_axisb = mount_scx11.axis[Mountaxis.BASE]
        mount_scx11_axisp = mount_scx11.axis[Mountaxis.POLAR]
        # --- simulation or not
        mount_scx11_axisb.real = False
        mount_scx11_axisp.real = False
        # ======= SCX11
        mount_scx11.speedslew(10,10)
        ha_start = 70
        dec_start = 60
        ha_target = 10
        dec_target = -25
        pierside_target = 1
        res = mount_scx11.hadec_travel_compute(ha_target, dec_target, pierside_target)