# -*- coding: utf-8 -*-
import numpy as np
import os

try:
    from .mountlog import Mountlog
except:
    from mountlog import Mountlog

# #####################################################################
# #####################################################################
# #####################################################################
# Class Mountmodpoi
# #####################################################################
# #####################################################################
# This class computes pointing model
# #####################################################################

class Mountmodpoi:

    _symbols = None
    _coefs = None
    _datas = None
    _ddatas = None
    _latitude = 45
    _name = "Not referenced"
    _mount_type = "HADEC"

    def load_datas(self, filename):
        with open(filename, 'r', errors='surrogateescape', newline='') as fid:
            lignes = fid.read()
        lignes = lignes.split("\n")
        ddatas = []
        for ligne in lignes:
            if len(ligne)<1:
                continue
            if ligne[0]=="#":
                key = "# --- Input parameters of the pointing model:"
                k = ligne.find(key)
                if k==0:
                    k = len(key)
                    self.name = ligne[k:].strip()
                key = "# --- mount_type:"
                k = ligne.find(key)
                if k==0:
                    k = len(key)
                    self.mount_type = ligne[k:].strip()
                key = "# --- latitude:"
                k = ligne.find(key)
                if k==0:
                    k = len(key)
                    self.latitude = ligne[k:].strip()
                continue
            if ligne[0]=="-":
                break
            self.log.print("{}".format(ligne))
            mots= ligne.split()
            if len(mots)<4:
                continue
            data = []
            #self.log.print("mots={}".format(mots))
            data.append(float(mots[0]))
            data.append(float(mots[1]))
            data.append(float(mots[2]))
            data.append(float(mots[3]))
            #self.log.print("{}".format(data))
            ddatas.append(data)
        self._ddatas = np.array(ddatas)
        return self._ddatas

    def save_datas(self, filename):
        if isinstance(self._ddatas,type(None)) == True:
            return
        ddatas = self._ddatas
        nl, nc = np.shape(ddatas)
        with open(filename, 'wt', errors='surrogateescape', newline='') as fid:
            # --- data
            texte = self.print_datas(False)
            fid.write(texte)
            # --- model coefs
            texte = self.print_coefs(False)
            fid.write(texte)

    def print_datas(self, toprint:bool=True):
        if isinstance(self._ddatas,type(None)) == True:
            return
        ddatas = self._ddatas
        nl, nc = np.shape(ddatas)
        texte = ""
        texte += f"# --- Input parameters of the pointing model: {self.name}\n"
        texte += f"# --- mount_type: {self.mount_type}\n"
        texte += f"# --- latitude: {self.latitude}\n"
        texte += "# " + "-"*80 + "\n"
        if self.mount_type.find("AZ")>=0:
            texte += "#    Az_app   Elev_app     dAz_obs dElev_obs   dAz_mod dElev_mod   dAz_res dElev_res\n"
        else:
            texte += "#    HA_app    Dec_app     dHA_obs  dDec_obs   dHA_mod  dDec_mod   dHA_res  dDec_res\n"
        texte += "#       deg        deg      arcmin    arcmin    arcmin    arcmin    arcmin    arcmin\n"
        for kl in range(nl):
            base  = ddatas[kl,0]
            pole = ddatas[kl,1]
            texte += f"{base:+11.6f} {pole:+10.6f}  "
            for kc in range(2,nc):
                d  = ddatas[kl,kc]
                texte += f" {d:+9.4f}"
            texte += "\n"
        texte += "-"*84 + "\n"
        means = np.mean(ddatas,0)
        stds = np.std(ddatas,0)
        maxi = np.max(ddatas,0)
        mini = np.min(ddatas,0)
        drange = maxi-mini
        texte += f"{means[0]:+11.4f} {means[1]:+10.4f}  "
        for kc in range(2,nc):
            texte += f" {means[kc]:+9.4f}"
        texte += " mean values\n"
        texte += f"{drange[0]:+11.4f} {drange[1]:+10.4f}  "
        for kc in range(2,nc):
            texte += f" {drange[kc]:+9.4f}"
        texte += " range values\n"
        texte += f"{stds[0]:+11.4f} {stds[1]:+10.4f}  "
        for kc in range(2,nc):
            texte += f" {stds[kc]:+9.4f}"
        texte += " std values\n"
        texte += "-"*84 + "\n"
        if toprint == True:
            self.log.print(texte, end ="")
        return texte

    def print_coefs(self, toprint:bool=True):
        if isinstance(self._coefs,type(None)) == True:
            return ""
        if isinstance(self._symbols,type(None)) == True:
            return ""
        texte = ""
        texte += f"--- Coefficients of the pointing model: {self.name}\n"
        coefs = self._coefs
        symbols = self._symbols
        nc, = np.shape(coefs)
        for kc in range(nc):
            coef = coefs[kc]
            symbol = symbols[kc]
            descr = ""
            if self.mount_type.find("AZ")>=0:
                if symbol == "IA":
                    descr = "Index error Az"
                elif symbol == "IE":
                    descr = "Index error Elev"
                elif symbol == "NPAE":
                    descr = "Non perpendicularity of Az and Elev axes"
                elif symbol == "CA":
                    descr = "Collimation error in Az"
                elif symbol == "AN":
                    descr = "North-South misalignment of Az axis"
                elif symbol == "AW":
                    descr = "East-West misalignment of Az axis"
                elif symbol == "ACEC":
                    descr = "Az centering error by cos(Az)"
                elif symbol == "ECEC":
                    descr = "Elev centering error by cos(Elev)"
                elif symbol == "ACES":
                    descr = "Az centering error by sin(Az)"
                elif symbol == "ECES":
                    descr = "Elev centering error by sin(Elev)"
                elif symbol == "NRX":
                    descr = "Nasmyth rotator displacement (vertical)"
                elif symbol == "NRY":
                    descr = "Nasmyth rotator displacement (horizontal)"
                for k in range(2,7):
                    sym = "ACEC"+str(k)
                    if symbol == sym:
                        descr = f"Az centering error by cos({k}*Az)"
                    sym = "ACES"+str(k)
                    if symbol == sym:
                        descr = f"Az centering error by sin({k}*Az)"
                    sym = "AN"+str(k)
                    if symbol == sym:
                        descr = f"North-South misalignment of Az axis by {k}*Az"
                    sym = "AW"+str(k)
                    if symbol == sym:
                        descr = f"East-West misalignment of Az axis by {k}*Az"
            else:
                if symbol == "IH":
                    descr = "Index error HA"
                elif symbol == "ID":
                    descr = "Index error Dec"
                elif symbol == "NP":
                    descr = "Non perpendicularity HA/Dec"
                elif symbol == "CH":
                    descr = "Non perpendicularity Dec/Opt"
                elif symbol == "ME":
                    descr = "Polar axis error in elevation"
                elif symbol == "MA":
                    descr = "Polar axis error E-W"
                elif symbol == "TF":
                    descr = "Tube flexure in sin(z)"
                elif symbol == "FO":
                    descr = "Fork Flexure"
                elif symbol == "DAF":
                    descr = "Flexure of cantilevered Dec axis"  # cantilevered = port-a-faux
                elif symbol == "HF":
                    descr = "Horseshoe Flexion"
                elif symbol == "TX":
                    descr = "Tube flexure in sin(z)"
                elif symbol == "DNP":
                    descr = "Dynamic non perpendicularity"
                elif symbol == "EHS":
                    descr = "sin(HA) effect in polar axis elevation"
                elif symbol == "EHC":
                    descr = "cos(HA) effect in polar axis elevation"
                elif symbol == "HCEC":
                    descr = "HA centering error by cos(HA)"
                elif symbol == "DCEC":
                    descr = "Dec centering error by cos(Dec)"
                elif symbol == "HCES":
                    descr = "HA centering error by sin(HA)"
                elif symbol == "DCES":
                    descr = "Dec centering error by sin(Dec)"
                elif symbol == "RWHS":
                    descr = "Ratio of wheels HA axis in sin(ratio*HA)"
                elif symbol == "RWHC":
                    descr = "Ratio of wheels HA axis in cos(ratio*HA)"
                elif symbol == "RWDS":
                    descr = "Ratio of wheels Dec axis in sin(ratio*Dec)"
                elif symbol == "RWDC":
                    descr = "Ratio of wheels Dec axis in cos(ratio*Dec)"
                else:
                    for k in range(2,5):
                        sym = "D"+str(k)+"HS"
                        if symbol == sym:
                            descr = f"Effect sin({k}*Dec) in Dec"
                        sym = "D"+str(k)+"HC"
                        if symbol == sym:
                            descr = f"Effect cos({k}*Dec) in Dec"
                        sym = "X"+str(k)+"HS"
                        if symbol == sym:
                            descr = f"Effect sin({k}*HA) in HA"
                        sym = "X"+str(k)+"HC"
                        if symbol == sym:
                            descr = f"Effect cos({k}*HA) in HA"
            if descr == "":
                descr = "Coefficient not described"
            texte += f"{symbol:4s} = {coef:+10.4f} arcmin ({descr})\n"
        texte += "-"*84 + "\n"
        if toprint == True:
            self.log.print(texte)
        return texte

    def _twolines(self, base, pole, latitude, symbols):
        xbases = []
        xpoles = []
        if self.mount_type.find("AZ")>=0:
            azr = np.radians(base)
            elevr = np.radians(pole)
            tane = np.tan(elevr)
            sine = np.sin(elevr)
            cose = np.cos(elevr)
            sina = np.sin(azr)
            cosa = np.cos(azr)
        else:
            har = np.radians(base)
            decr = np.radians(pole)
            tand = np.tan(decr)
            sind = np.sin(decr)
            cosd = np.cos(decr)
            sinh = np.sin(har)
            cosh = np.cos(har)
        latr = np.radians(latitude)
        cosl = np.cos(latr)
        sinl = np.sin(latr)
        #cosz = sinl*sind+cosl*cosh*cosd # = sinelev
        #sinz = np.sqrt(1-cosz*cosz)
        ks = -1
        for symbol in symbols:
            ks += 1
            if self.mount_type.find("AZ")>=0:
                xaz = None
                xelev = None
                if symbol == "IA":
                    xaz = 1.0
                    xelev = 0.0
                if symbol == "IE":
                    xaz = 0.0
                    xelev = 1.0
                if symbol == "NPAE":
                    xaz = tane
                    xelev = 0.0
                if symbol == "CA":
                    xaz = 1./cose
                    xelev = 0.0
                if symbol == "AN":
                    xaz = sina*tane
                    xelev = cosa
                if symbol == "AW":
                    xaz = -cosa*tane
                    xelev = sina
                if symbol == "ACEC":
                    xaz = cosa
                    xelev = 0.0
                if symbol == "ECEC":
                    xaz = 0.
                    xelev = cose
                if symbol == "ACES":
                    xaz = sina
                    xelev = 0.0
                if symbol == "ECES":
                    xaz = 0.
                    xelev = sine
                if symbol == "NRX":
                    xaz = 1.0
                    xelev = -sine
                if symbol == "NRY":
                    xaz = tane
                    xelev = cose
                for k in range(2,7):
                    sym = "ACEC"+str(k)
                    if symbol == sym:
                        xaz = np.cos(k*azr)
                        xelev = 0.0
                    sym = "ACES"+str(k)
                    if symbol == sym:
                        xaz = np.sin(k*azr)
                        xelev = 0.0
                    sym = "AN"+str(k)
                    if symbol == sym:
                        xaz = np.sin(k*azr)*tane
                        xelev = np.cos(k*azr)
                    sym = "AW"+str(k)
                    if symbol == sym:
                        xaz = -np.cos(k*azr)*tane
                        xelev = np.sin(k*azr)
                if xaz == None or xelev == None:
                    continue
                xbases.append(xaz)
                xpoles.append(xelev)
            else:
                xha = None
                xdec = None
                if symbol == "IH":
                    xha  = 1
                    xdec = 0
                if symbol == "ID":
                    xha  = 0
                    xdec = 1
                if symbol == "NP":
                    xha  = tand
                    xdec = 0
                if symbol == "CH":
                    xha  = 1./cosd
                    xdec = 0
                if symbol == "ME":
                    xha  = sinh*tand
                    xdec = cosh
                if symbol == "MA":
                    xha  = -cosh*tand
                    xdec = sinh
                if symbol == "TF":
                    xha  = cosl*sinh/cosd
                    xdec = cosl*cosh*sind-sinl*cosd
                if symbol == "FO":
                    xha  = 0
                    xdec = cosh
                if symbol == "DAF":
                    xha  = -cosl*cosh-sinl*tand
                    xdec = 0
                if symbol == "HF":
                    xha  = -sinh/cosd
                    xdec = 0
                if symbol == "TX":
                    xha  = cosl*sinh*cosd / (sind*sinl+cosd*cosh*cosl)
                    xdec = (cosl*cosh*sind-sinl*cosd) / (sind*sinl+cosd*cosh*cosl)
                if symbol == "DNP":
                    xha  = sinh*tand
                    xdec = 0
                if symbol == "EHS":
                    xha  = sinh*sinh*tand
                    xdec = sinh*cosh
                if symbol == "EHC":
                    xha  = sinh*cosh*tand
                    xdec = cosh*cosh
                if symbol == "HCEC":
                    xha  = cosh
                    xdec = 0
                if symbol == "DCEC":
                    xha  = 0
                    xdec = cosd
                if symbol == "HCES":
                    xha  = sinh
                    xdec = 0
                if symbol == "DCES":
                    xha  = 0
                    xdec = sind
                for k in range(2,5):
                    sym = "D"+str(k)+"HS"
                    if symbol == sym:
                        xha  = 0
                        xdec = np.sin(k*decr)
                    sym = "D"+str(k)+"HC"
                    if symbol == sym:
                        xha  = 0
                        xdec = np.cos(k*decr)
                    sym = "X"+str(k)+"HS"
                    if symbol == sym:
                        xha  = np.sin(k*har)/cosd
                        xdec = 0
                    sym = "X"+str(k)+"HC"
                    if symbol == sym:
                        xha  = np.cos(k*har)/cosd
                        xdec = 0
                if symbol == "RWHS":
                    xha  = np.sin(self._ratio_basis_wheel_puley*har)
                    xdec = 0
                if symbol == "RWHC":
                    xha  = np.cos(self._ratio_basis_wheel_puley*har)
                    xdec = 0
                if symbol == "RWDS":
                    xha  = 0
                    xdec = np.sin(self._ratio_polar_wheel_puley*decr)
                if symbol == "RWDC":
                    xha  = 0
                    xdec = np.cos(self._ratio_polar_wheel_puley*decr)
                if xha == None or xdec == None:
                    continue
                xbases.append(xha)
                xpoles.append(xdec)
        return xbases, xpoles

    def ddata_sigmaclipping(self, kappa):
        if isinstance(self._ddatas,type(None)) == True:
            return
        ddatas = self._ddatas
        means = np.mean(ddatas,0)
        stds = np.std(ddatas,0)
        threshold_minis = means - kappa*stds
        threshold_maxis = means + kappa*stds
        ddatanews = []
        nl, nc = np.shape(ddatas)
        for kl in range(nl):
            valid = True
            for kc in list([2, 3]):
                val = ddatas[kl,kc]
                threshold_mini = threshold_minis[kc]
                threshold_maxi = threshold_maxis[kc]
                if val<threshold_mini or val>threshold_maxi:
                    valid = False
            if valid == True:
                ddatanew = list(ddatas[kl,:])
                ddatanews.append(ddatanew)
            else:
                self.log.print(f"Elimine {ddatas[kl,0:4]} ...")
        ddatanews = np.array(ddatanews)
        self._ddatas = ddatanews

    def ddata2coefs(self):
        if isinstance(self._symbols,type(None)) == True:
            return
        if isinstance(self._ddatas,type(None)) == True:
            return
        if isinstance(self._latitude,type(None)) == True:
            return
        ddatas = self._ddatas
        symbols = self._symbols
        latitude = self._latitude
        # === conversion array ha,dec,dha,ddec -> matrixes
        nl, nc = np.shape(ddatas)
        ns = len(symbols)
        x = np.zeros((2*nl,ns))
        y = np.zeros(2*nl)
        # --- Loop over pointings
        for kl in range(nl):
            kl_base = 2*kl
            kl_pole = kl_base+1
            ddata = ddatas[kl,:]
            base, pole , dbase, dpole = ddata
            y[kl_base]  = dbase
            y[kl_pole] = dpole
            xbases, xpoles = self._twolines(base, pole, latitude, symbols)
            # --- Loop over symbols
            for ks in range(ns):
                x[kl_base, ks] = xbases[ks]
                x[kl_pole, ks] = xpoles[ks]

        # === Solve: y = x*a
        a, residuals, rank, singular_values = np.linalg.lstsq(x, y, rcond=None)

        # === Another way to solve: x = x*a
        # res = sm.OLS(y,x).fit()
        # a = res.params
        coefs = a

        # === Compute yp = x*a
        #invx = np.linalg.inv(x)
        yp = np.dot(x,a)
        dy = yp-y

        # === Add new columns in the inputs
        ddata2s = np.zeros((nl,8))
        for ky in range(nl):
            # --- copy only the first four columns
            for kc in range(4):
                ddata2s[ky,kc] = ddatas[ky,kc]
            # --- add or modify the new columns
            ddata2s[ky,4] = yp[2*ky]
            ddata2s[ky,5] = yp[2*ky+1]
            ddata2s[ky,6] = dy[2*ky]
            ddata2s[ky,7] = dy[2*ky+1]

        # ---
        self.coefs = coefs
        self.ddatas = ddata2s
        return coefs, ddata2s

    def coefs2ddata(self, base, pole):
        """
        Compute the dbase and dpole from a couple of base,pole and a pointing model
        """
        if isinstance(self._symbols,type(None)) == True or isinstance(self._coefs,type(None)) == True:
            dbase = 0
            dpole = 0
        else:
            latitude = self._latitude
            coefs = self._coefs
            symbols = self._symbols
            ns = len(symbols)
            x = np.zeros((2,ns))
            xbases, xpoles = self._twolines(base, pole, latitude, symbols)
            # --- Loop over symbols
            for ks in range(ns):
                x[0, ks] = xbases[ks]
                x[1, ks] = xpoles[ks]
            # === Compute yp = x*a
            dbase, dpole = np.dot(x,coefs)
        return dbase, dpole

    def load_coefs(self, filename):
        with open(filename, 'r', errors='surrogateescape', newline='') as fid:
            lignes = fid.read()
        lignes = lignes.split("\n")
        symbs = []
        coefs = []
        found = False
        key = "--- Coefficients of the pointing model:"
        for ligne in lignes:
            if ligne.find(key)==0:
                found = True
                k = len(key)
                self.name = ligne[k:].strip()
                continue
            elif found == False:
                continue
            if len(ligne)<1:
                continue
            if ligne[0]=="-":
                break
            #self.log.print("{}".format(ligne))
            mots= ligne.split()
            if len(mots)<4:
                continue
            #self.log.print(f"mots = {mots}")
            symbs.append(mots[0].strip())
            coefs.append(float(mots[2].strip()))
            #self.log.print(f"mot0 = {mots[0].strip()}")
            #self.log.print(f"mot2 = {mots[2].strip()}")
        if len(symbs)>0:
            self.symbols = symbs
            self.coefs = np.array(coefs)
        return symbs

    def _set_symbols(self, symbols):
        if symbols=="" or symbols=="default":
            if self.mount_type.find("AZ")>=0:
                self._symbols = ["IA", "IE", "AN", "AW", "NPAE", "CA", "ACEC", "ECEC", "ACES", "ECES"]
            else:
                self._symbols = ["IH", "ID", "ME", "MA", "DNP",  "CH", "HCEC", "DCEC", "HCES", "DCES"]
        else:
            self._symbols = symbols

    def _get_symbols(self):
        return self._symbols

    def _set_coefs(self, coefs):
        self._coefs = coefs

    def _get_coefs(self):
        return self._coefs

    def _set_datas(self, datas):
        self._datas = datas

    def _get_datas(self):
        return self._datas

    def _set_ddatas(self, ddatas):
        self._ddatas = ddatas

    def _get_ddatas(self):
        return self._ddatas

    def _set_name(self, name):
        self._name = name

    def _get_name(self):
        return self._name

    def _set_ratio_basis_wheel_puley(self, ratio):
        self._ratio_basis_wheel_puley = ratio

    def _get_ratio_basis_wheel_puley(self):
        return self._ratio_basis_wheel_puley

    def _set_ratio_polar_wheel_puley(self, ratio):
        self._ratio_polar_wheel_puley = ratio

    def _get_ratio_polar_wheel_puley(self):
        return self._ratio_polar_wheel_puley

    def _set_mount_type(self, mount_type):
        if self._mount_type != mount_type:
            self._mount_type = mount_type
            self._symbols = None
            self._coefs = None
            self._datas = None
            self._ddatas = None

    def _get_mount_type(self):
        return self._mount_type

    symbols = property(_get_symbols, _set_symbols)
    coefs = property(_get_coefs, _set_coefs)
    datas = property(_get_datas, _set_datas)
    ddatas = property(_get_ddatas, _set_ddatas)
    name = property(_get_name, _set_name)
    ratio_polar_wheel_puley = property(_get_ratio_polar_wheel_puley, _set_ratio_polar_wheel_puley)
    ratio_basis_wheel_puley = property(_get_ratio_basis_wheel_puley, _set_ratio_basis_wheel_puley)
    mount_type = property(_get_mount_type, _set_mount_type)

    def __init__(self, *args, **kwars):
        self._symbols = None
        self._coefs = None
        self._datas = None
        self._ddatas = None
        self._latitude = 45
        self._ratio_polar_wheel_puley = 1
        self._ratio_basis_wheel_puley = 1
        self._name = "Not referenced"
        self._mount_type = "HADEC"
        # === Log
        self.log = Mountlog(os.path.basename(__file__))


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

if __name__ == "__main__":

    example = 1
    print("Example       = {}".format(example))

    if example == 1:
        modpoi = Mountmodpoi()
        # --- Initiate a pointing model
        modpoi.latitude = 45 # deg
        modpoi.symbols = ["IH", "ID", "ME", "MA"]
        modpoi.coefs   = [  20,  -10,   -5,   15] # arcmin
        # --- Some pointings
        hadecs = []
        hadecs.extend([[-30, -10], [-10, -10], [10, -10], [30, -10]])
        hadecs.extend([[-30, -10], [-10, -10], [10, -10], [30, -10]])
        # --- Compute the dhadecs as a numpy array
        dhadecs = []
        nstar = len(hadecs)
        for hadec in hadecs:
            ha, dec = hadec
            dha, ddec = modpoi.coefs2ddata(ha, dec)
            dhadec = [ha, dec, dha, ddec]
            dhadecs.append(dhadec)
        modpoi.ddatas = np.array(dhadecs)
        # --- Compute the pointing model
        modpoi.ddata2coefs()
        # --- print results
        modpoi.print_datas()
        modpoi.print_coefs()

    if example == 2:
        modpoi = Mountmodpoi()
        # --- load a data file: HA(deg) Dec(deg) dHA(arcmin) dDec(arcmin)
        filename = "C:/d/t1m_pic_du_midi/mission_20210329/20210331/pointages_zwo.txt"
        modpoi.load_datas(filename)
        modpoi.print_datas()
        # --- Initiate a pointing model
        modpoi.latitude = 42 # deg
        #modpoi.symbols = ["IH", "ID", "ME", "MA", "NP", "CH", "TF", "FO", "DNP", "EHS", "EHC", "HCEC", "DCEC", "HCES", "DCES", "D2HS", "D2HC", "X2HS", "X2HC"]
        modpoi.symbols = ["IH", "ID", "ME", "MA", "NP", "CH", "TF", "FO", "DNP", "EHS", "EHC", "HCEC", "DCEC", "HCES", "DCES", "D2HS", "D2HC", "D3HS", "D3HC", "D4HS", "D4HC", "X2HS", "X2HC", "X3HS", "X3HC", "X4HS", "X4HC"]
        modpoi.symbols = ["IH", "ID", "ME", "MA", "DNP", "HCEC", "DCEC", "HCES", "DCES", "D2HS", "D2HC", "D4HS", "D4HC"]
        # --- Compute the pointing model
        modpoi.name = "{} coefs".format(len(modpoi.symbols))
        modpoi.ddata2coefs()
        # --- print results
        modpoi.print_datas()
        modpoi.print_coefs()
        filename = "C:/d/t1m_pic_du_midi/mission_20210329/20210331/modpoi_zwo.txt"
        modpoi.save_datas(filename)