horizon.py
12.5 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
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__()