horizon.py 12.7 KB
from math import floor, ceil, atan2, asin, sin, cos, tan

try:
    from .mechanics import Mechanics
except:
    from mechanics import Mechanics

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

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

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

# ========================================================
# ========================================================
# === HORIZON
# ========================================================
# ========================================================

class HorizonException(GuitastroException):
    """Exception raised for errors in the Horizon class.
    """

    BAD_INPUT_FORMAT = 0
    ELEMENTS_HAVE_NOT_SAME_LEN = 1
    ELEMENTS_MUST_BE_CONTAINER = 2

    errors = [""]*3
    errors[BAD_INPUT_FORMAT] = "Bad input format"
    errors[ELEMENTS_HAVE_NOT_SAME_LEN] = "All elements have not the same length"
    errors[ELEMENTS_MUST_BE_CONTAINER] = "Elements must be containers (list, tuple)"


class Horizon(Mechanics, GuitastroTools):
    """ Class to describe an Horizon
    """
# ========================================================
# === attributs
# ========================================================

    _home = Home()
    _TYPE_HORIZON_ALTAZ = 0
    _TYPE_HORIZON_HADEC = 1
    _horizon_type = _TYPE_HORIZON_ALTAZ
    _amers_raws = [] ; # raw values
    _amers_altaz = [] ; # raw values in degrees
    _computed_interpolation_altaz = 0
    _horizon_az = [] ; # resampled values
    _horizon_elev = [] ; # resampled values
    _computed_interpolation_hadec = 0
    _horizon_dec = [] ; # resampled values
    _horizon_harise = [] ; # resampled values
    _horizon_haset = [] ; # resampled values

# ========================================================
# === internal methods : Generals
# ========================================================

    def _init_horizon(self, home, *args):
        """ Object initialization
        """
        # --- home
        if isinstance(home, Home) == True:
            self._home = home
        else:
            self._home = Home(home)
        #print("self._home={}".format(self._home.gps))
        # --- *args
        #print("args={}".format(args))
        self._amers_raws = []
        self._computed_interpolation_altaz = 0
        self._computed_interpolation_hadec = 0
        self._horizon_type = self._TYPE_HORIZON_ALTAZ
        if (len(args)==0):
            h = ((0,0), (360,0))
            self._horizon_args2amers(h)
        else:
            self._horizon_args2amers(*args)

# ========================================================
# === internal methods :
# ========================================================

    def _horizon_args2amers(self, *args):
        """ Record the horizon definition as amers
        """
        # --- case of void args
        if len(args)==0:
            amers = [(0,0), (360,0)]
        else:
            amers = args[0]
        # --- Verify that all the elements have the same length
        na = len(amers)
        lamers = []
        for ka in range(na):
            try:
                n = len(amers[ka])
            except:
                msg = f"The {ka+1}th element {amers[ka]} is not compatible"
                raise HorizonException("ELEMENTS_MUST_BE_CONTAINER", msg)
            if ka == 0:
                n0 = n
            elif n != n0:
                msg = f"The {ka+1}th element {amers[ka]} does not have the same length than the first one {amers[0]}"
                raise HorizonException("ELEMENTS_HAVE_NOT_SAME_LEN", msg)
            lamers.append(amers[ka])
        self._amers_raws = lamers

    def _amers_decode(self):
        """ Decodes the amers

        _TYPE_HORIZON_ALTAZ : (az1,elev1) (az2,elev2) etc.
        _TYPE_HORIZON_HADEC : (dec1,harise1,haset1) (dec2,harise2,haset2) etc.

        """
        amers = (self._amers_raws).copy()
        # --- identify the type of coordinates.
        # --- how many coordinates in each element ? (2=ALTAZ 3=HADEC)
        n2 = 0
        n3 = 0
        print("amers={}".format(amers))
        for amer in amers:
            #print("amer={}".format(amer))
            n = len(amer)
            if (n==2):
                n2 += 1
            if (n==3):
                n3 += 1
        #print("n2={} n3={}".format(n2,n3))
        if (n2>0) and (n3==0):
            self._horizon_type = self._TYPE_HORIZON_ALTAZ
        elif (n2==0) and (n3>0):
            self._horizon_type = self._TYPE_HORIZON_HADEC
        else:
            msg = f"n2={n2} n3={n3}"
            raise HorizonException("BAD_INPUT_FORMAT", msg)
        #print("self._horizon_type={}".format(self._horizon_type))
        # --- convert coordinates into degrees
        angle = Angle()
        deg_amers = []
        for amer in amers:
            degs = []
            for angle_raw in amer:
                angle.angle(angle_raw)
                angle_deg = (angle%360).deg()
                degs.append(angle_deg)
            deg_amers.append(degs)
        # --- if horizontype = HADEC, convert each amer into altaz
        if self._horizon_type == self._TYPE_HORIZON_HADEC:
            latitude = (self._home).latitude * self._DR
            deg_amers = []
            for amer in amers:
                dec = amer[0] * self._DR
                ha1 = amer[1] * self._DR
                ha2 = amer[2] * self._DR
                az1, elev1 = self._mc_hd2ah(ha1, dec, latitude)
                degs = [ az1/self._DR, elev1/self._DR]
                deg_amers.append(degs)
                az2, elev2 = self._mc_hd2ah(ha2, dec, latitude)
                degs = [ az2/self._DR, elev2/self._DR]
                deg_amers.append(degs)
        # --- sort in azimuts before interpolation
        amers = list(sorted(deg_amers, key=lambda item: item[0]))
        #print("amers={}".format(amers))
        # --- add two extra items to be able to interpolate between [0-360]
        first = (amers[0]).copy()
        last = (amers[-1]).copy()
        last[0] -= 360
        amers.insert(0,last)
        first[0] += 360
        amers.append(first)
        #print("amers={}".format(amers))
        self._amers_altaz = amers
        self._computed_interpolation_altaz = 0
        # --- we can make linear interpolations for a resampling

    def _amers_interpolation_altaz(self,az_sampling_deg=1):
        """
        """
        amers = self._amers_altaz
        if amers == [] or self._computed_interpolation_altaz == 1:
            return
        self._horizon_az = []
        self._horizon_elev = []
        az = 0
        k1 = -1
        az2 = -1000 ; # < -360 to oblige the first if az>az2
        while (az<=360):
            if az > az2:
                k1 += 1
                az1, elev1 = amers[k1]
                az2, elev2 = amers[k1+1]
                daz = az2-az1
                delev = elev2-elev1
                continue
            # --- we can interpolate
            elev = (az-az1)/daz*delev+elev1
            self._horizon_az.append(az)
            self._horizon_elev.append(elev)
            az += az_sampling_deg
        self._computed_interpolation_altaz = 1

    def _amers_interpolation_hadec(self,dec_sampling_deg=1):
        """
        """
        amers = self._amers_altaz
        if amers == [] or self._computed_interpolation_hadec == 1:
            return
        if self._computed_interpolation_altaz == 0:
            self._amers_interpolation_altaz(1)
        # ---
        az_sampling = self._horizon_az[1] - self._horizon_az[0]
        # --- declination limits for visibility
        latitude_deg = (self._home).latitude
        latitude_rad = latitude_deg * self._DR
        coslatitude = cos(latitude_rad)
        sinlatitude = sin(latitude_rad)
        declim_inf = -90
        declim_sup = 90
        if (latitude_deg>=0):
            declim_inf = latitude_deg - 90
        else:
            declim_sup = latitude_deg + 90
        dec1 = ceil(declim_inf)
        dec2 = floor(declim_sup)
        # ---
        self._horizon_dec = []
        self._horizon_harise = []
        self._horizon_haset = []
        dec_deg = -90
        while (dec_deg<=90):
            if (dec_deg <= dec1) or (dec_deg >= dec2):
                # --- never visible
                ha_set = 0
                ha_rise = 0
                continue
            else:
                ha_rise_computed = 0
                ha_set_computed = 0
                ha_rise = -180
                ha_set = 180
                dec = dec_deg * self._DR
                cosdec = cos(dec)
                sindec = sin(dec)
                tandec = tan(dec)
                sinlatitude_sindec = sinlatitude*sindec
                coslatitude_cosdec = coslatitude*cosdec
                tandec_coslatitude = tandec*coslatitude
                ha_deg = -180
                while (ha_deg <= 180):
                    ha = ha_deg * self._PI
                    sinha = sin(ha)
                    cosha = cos(ha)
                    az=atan2(sinha,cosha*sinlatitude-tandec_coslatitude)
                    el=asin(sinlatitude_sindec+coslatitude_cosdec*cosha)
                    az /= self._DR
                    el /= self._DR
                    if (az<0):
                        az += 360
                    kaltaz = int ((az - self._horizon_az[0])/360.*az_sampling)
                    horizon_elev = self._horizon_elev[kaltaz]
                    if (ha_rise_computed == 0):
                        if (el >= horizon_elev):
                            ha_rise = ha_deg
                            ha_rise_computed = 1
                    else:
                        if (ha_set_computed == 0) and (el <= horizon_elev):
                            ha_set = ha_deg
                            ha_set_computed = 1
                            break
            self._horizon_dec.append(dec_deg)
            self._horizon_harise.append(ha_rise)
            self._horizon_haset.append(ha_set)
        self._computed_interpolation_hadec = 1

# ========================================================
# === Horizon methods
# ========================================================

    def horizon(self, home, *args):
        """ Object initialization

        """
        self._init_horizon(home, *args)
        return args

# ========================================================
# === get/set methods
# ========================================================

    def _get_horizon_altaz(self):
        if (self._computed_interpolation_altaz == 0):
            self._amers_decode()
            az_sampling_deg=1
            self._amers_interpolation_altaz(az_sampling_deg)
        return self._horizon_az , self._horizon_elev

    def _set_horizon_altaz(self, *args):
        self._horizon_args2amers(*args)

    def _get_horizon_hadec(self):
        if (self._computed_interpolation_hadec == 0):
            self._amers_decode()
            hadec_sampling_deg=1
            self._amers_interpolation_hadec(hadec_sampling_deg)
        return self._horizon_dec , self._horizon_harise , self._horizon_haset

    def _set_horizon_hadec(self, *args):
        self._horizon_args2amers(*args)

    def get_horizon_raw_amers(self):
        return self._amers_raws

    horizon_hadec = property(_get_horizon_hadec, _set_horizon_hadec)
    horizon_altaz = property(_get_horizon_altaz, _set_horizon_altaz)

# ========================================================
# === special methods
# ========================================================

    def __init__(self, home, *args):
        """ Object initialization
        """
        self._init_horizon(home, *args)



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

if __name__ == "__main__":

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

    if example == 1:
        """
        Difference coord angles
        """
        home = Home("GPS 2.345 E +56.3456 2300")
        hor1 = Horizon(home)
        #hor1 = Horizon(home, ( (0, 10), (60, 12), (180, 20), (270, 22), (360, 10)))
        #hor1 = Horizon(home, [ [0, 10], [60, 12], [180, 20], [270, 22], [360, 10]])
        print(f"{home.gps} = {home.mpc}")
        print(f"horizon = {hor1.horizon_altaz}")