horizon.py 12.5 KB
from math import floor, ceil, atan2, asin, sin, cos, tan
import doctest
from .mechanics import Mechanics
from .angles import Angle
from .home import Home

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

class Horizon(Mechanics):
    """ 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, 20), (10,25), (67, 23)
            h = (0,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
        #print("len={}".format(len(args)))
        if (len(args)==0):
            amers = [(0,0), (180,0)]
        # --- first element is the tuple of amers
        amers = args[0]
        #print("amers={}".format(amers))
        # --- verify that the first element of amers is also a tuple
        if isinstance(amers[0], tuple) == False:
            if (len(amers)>=2):
                amers = [amers, ]
            else:
                raise Exception
                return ""
        else:
            amers = list(amers)
        #print("amers={}".format(amers))
        self._amers_raws = amers
        
    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:
            raise Exception
        #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)
    
# ========================================================
# === debug methods
# ========================================================
    
    def infos(self, action) -> None:
        """ To get informations about this class
        
        :param action: A command to run a debug action (see examples).
        :type action: string
        
        :Example:
            
        Horizon().infos("doctest")
        Horizon().infos("doc_methods")
        Horizon().infos("internal_attributes")
        Horizon().infos("public_methods")        
        """
        if (action == "doc_methods"):
            publics = [x for x in dir(self) if x[0]!="_"]
            for public in publics:
                varname = "{}".format(public)
                if (callable(getattr(self,varname))==True):
                    print("\n{:=^40}".format(" method "+varname+" "))
                    t = "Horizon()."+varname+".__doc__"
                    tt =eval(t)
                    print(tt)
        if (action == "doctest"):
            if __name__ == "__main__":
                print("\n{:~^40}".format("doctest"))
                doctest.testmod(verbose=False)
        if (action == "internal_attributes"):
            internals = [x for x in dir(self) if x[0]=="_" and x[1]!="_"]
            for internal in internals:
                varname = "{}".format(internal)
                #if (hasattr(self,varname)==True):
                if (callable(getattr(self,varname))==False):
                    print(varname + "=" + str(getattr(self,varname)))
        if (action == "public_methods"):
            publics = [x for x in dir(self) if x[0]!="_"]
            for public in publics:
                varname = "{}".format(public)
                if (callable(getattr(self,varname))==True):
                    print(varname)

# ========================================================
# === special methods
# ========================================================
        
    def __init__(self, home, *args):
        """ Object initialization
        """
        self._init_horizon(home, *args)
		  # super().__init__()