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__()