The requested changes are too big and content was truncated. Show full diff
This diff has been collapsed as it changes many lines, (1420 lines changed) Show them Hide them | |||
@@ -0,0 +1,1420 | |||
|
1 | """ | |
|
2 | The module ASTRO_COORDS.py gathers classes and functions for coordinates transformation. Additiona- | |
|
3 | lly a class EquatorialCorrections and celestial bodies are defined. The first of these is to correct | |
|
4 | any error in the location of the body and the second to know the location of certain celestial bo- | |
|
5 | dies in the sky. | |
|
6 | ||
|
7 | MODULES CALLED: | |
|
8 | OS, NUMPY, NUMERIC, SCIPY, TIME_CONVERSIONS | |
|
9 | ||
|
10 | MODIFICATION HISTORY: | |
|
11 | Created by Ing. Freddy Galindo (frederickgalindo@gmail.com). ROJ Sep 20, 2009. | |
|
12 | """ | |
|
13 | ||
|
14 | import numpy | |
|
15 | #import Numeric | |
|
16 | import scipy.interpolate | |
|
17 | import os | |
|
18 | import sys | |
|
19 | from utils import TimeTools | |
|
20 | from utils import Misc_Routines | |
|
21 | ||
|
22 | class EquatorialCorrections(): | |
|
23 | def __init__(self): | |
|
24 | """ | |
|
25 | EquatorialCorrections class creates an object to call methods to correct the loca- | |
|
26 | tion of the celestial bodies. | |
|
27 | ||
|
28 | Modification History | |
|
29 | -------------------- | |
|
30 | Converted to Object-oriented Programming by Freddy Galindo, ROJ, 27 September 2009. | |
|
31 | """ | |
|
32 | ||
|
33 | pass | |
|
34 | ||
|
35 | def co_nutate(self,jd,ra,dec): | |
|
36 | """ | |
|
37 | co_nutate calculates changes in RA and Dec due to nutation of the Earth's rotation | |
|
38 | Additionally it returns the obliquity of the ecliptic (eps), nutation in the longi- | |
|
39 | tude of the ecliptic (d_psi) and nutation in the pbliquity of the ecliptic (d_eps). | |
|
40 | ||
|
41 | Parameters | |
|
42 | ---------- | |
|
43 | jd = Julian Date (Scalar or array). | |
|
44 | RA = A scalar o array giving the Right Ascention of interest. | |
|
45 | Dec = A scalar o array giving the Right Ascention of interest. | |
|
46 | ||
|
47 | Return | |
|
48 | ------ | |
|
49 | d_ra = Correction to ra due to nutation. | |
|
50 | d_dec = Correction to dec due to nutation. | |
|
51 | ||
|
52 | Examples | |
|
53 | -------- | |
|
54 | >> Julian = 2462088.7 | |
|
55 | >> Ra = 41.547213 | |
|
56 | >> Dec = 49.348483 | |
|
57 | >> [d_ra,d_dec,eps,d_psi,d_eps] = co_nutate(julian,Ra,Dec) | |
|
58 | >> print d_ra, d_dec, eps, d_psi, d_eps | |
|
59 | [ 15.84276651] [ 6.21641029] [ 0.4090404] [ 14.85990198] [ 2.70408658] | |
|
60 | ||
|
61 | Modification history | |
|
62 | -------------------- | |
|
63 | Written by Chris O'Dell, 2002. | |
|
64 | Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009. | |
|
65 | """ | |
|
66 | ||
|
67 | jd = numpy.atleast_1d(jd) | |
|
68 | ra = numpy.atleast_1d(ra) | |
|
69 | dec = numpy.atleast_1d(dec) | |
|
70 | ||
|
71 | # Useful transformation constants | |
|
72 | d2as = numpy.pi/(180.*3600.) | |
|
73 | ||
|
74 | # Julian centuries from J2000 of jd | |
|
75 | T = (jd - 2451545.0)/36525.0 | |
|
76 | ||
|
77 | # Must calculate obliquity of ecliptic | |
|
78 | [d_psi, d_eps] = self.nutate(jd) | |
|
79 | d_psi = numpy.atleast_1d(d_psi) | |
|
80 | d_eps = numpy.atleast_1d(d_eps) | |
|
81 | ||
|
82 | eps0 = (23.4392911*3600.) - (46.8150*T) - (0.00059*T**2) + (0.001813*T**3) | |
|
83 | # True obliquity of the ecliptic in radians | |
|
84 | eps = (eps0 + d_eps)/3600.*Misc_Routines.CoFactors.d2r | |
|
85 | ||
|
86 | # Useful numbers | |
|
87 | ce = numpy.cos(eps) | |
|
88 | se = numpy.sin(eps) | |
|
89 | ||
|
90 | # Convert Ra-Dec to equatorial rectangular coordinates | |
|
91 | x = numpy.cos(ra*Misc_Routines.CoFactors.d2r)*numpy.cos(dec*Misc_Routines.CoFactors.d2r) | |
|
92 | y = numpy.sin(ra*Misc_Routines.CoFactors.d2r)*numpy.cos(dec*Misc_Routines.CoFactors.d2r) | |
|
93 | z = numpy.sin(dec*Misc_Routines.CoFactors.d2r) | |
|
94 | ||
|
95 | # Apply corrections to each rectangular coordinate | |
|
96 | x2 = x - (y*ce + z*se)*d_psi*Misc_Routines.CoFactors.s2r | |
|
97 | y2 = y + (x*ce*d_psi - z*d_eps)*Misc_Routines.CoFactors.s2r | |
|
98 | z2 = z + (x*se*d_psi + y*d_eps)*Misc_Routines.CoFactors.s2r | |
|
99 | ||
|
100 | # Convert bask to equatorial spherical coordinates | |
|
101 | r = numpy.sqrt(x2**2. + y2**2. + z2**2.) | |
|
102 | xyproj =numpy.sqrt(x2**2. + y2**2.) | |
|
103 | ||
|
104 | ra2 = x2*0.0 | |
|
105 | dec2 = x2*0.0 | |
|
106 | ||
|
107 | xyproj = numpy.atleast_1d(xyproj) | |
|
108 | z = numpy.atleast_1d(z) | |
|
109 | r = numpy.atleast_1d(r) | |
|
110 | x2 = numpy.atleast_1d(x2) | |
|
111 | y2 = numpy.atleast_1d(y2) | |
|
112 | z2 = numpy.atleast_1d(z2) | |
|
113 | ra2 = numpy.atleast_1d(ra2) | |
|
114 | dec2 = numpy.atleast_1d(dec2) | |
|
115 | ||
|
116 | w1 = numpy.where((xyproj==0) & (z!=0)) | |
|
117 | w2 = numpy.where(xyproj!=0) | |
|
118 | ||
|
119 | # Calculate Ra and Dec in radians (later convert to degrees) | |
|
120 | if w1[0].size>0: | |
|
121 | # Places where xyproj=0 (point at NCP or SCP) | |
|
122 | dec2[w1] = numpy.arcsin(z2[w1]/r[w1]) | |
|
123 | ra2[w1] = 0 | |
|
124 | ||
|
125 | if w2[0].size>0: | |
|
126 | # Places other than NCP or SCP | |
|
127 | ra2[w2] = numpy.arctan2(y2[w2],x2[w2]) | |
|
128 | dec2[w2] = numpy.arcsin(z2[w2]/r[w2]) | |
|
129 | ||
|
130 | # Converting to degree | |
|
131 | ra2 = ra2/Misc_Routines.CoFactors.d2r | |
|
132 | dec2 = dec2/Misc_Routines.CoFactors.d2r | |
|
133 | ||
|
134 | w = numpy.where(ra2<0.) | |
|
135 | if w[0].size>0: | |
|
136 | ra2[w] = ra2[w] + 360. | |
|
137 | ||
|
138 | # Return changes in Ra and Dec in arcseconds | |
|
139 | d_ra = (ra2 -ra)*3600. | |
|
140 | d_dec = (dec2 - dec)*3600. | |
|
141 | ||
|
142 | return d_ra, d_dec, eps, d_psi, d_eps | |
|
143 | ||
|
144 | def nutate(self,jd): | |
|
145 | """ | |
|
146 | nutate returns the nutation in longitude and obliquity for a given Julian date. | |
|
147 | ||
|
148 | Parameters | |
|
149 | ---------- | |
|
150 | jd = Julian ephemeris date, scalar or vector. | |
|
151 | ||
|
152 | Return | |
|
153 | ------ | |
|
154 | nut_long = The nutation in longitude. | |
|
155 | nut_obliq = The nutation in latitude. | |
|
156 | ||
|
157 | Example | |
|
158 | ------- | |
|
159 | >> julian = 2446895.5 | |
|
160 | >> [nut_long,nut_obliq] = nutate(julian) | |
|
161 | >> print nut_long, nut_obliq | |
|
162 | -3.78793107711 9.44252069864 | |
|
163 | ||
|
164 | >> julians = 2415020.5 + numpy.arange(50) | |
|
165 | >> [nut_long,nut_obliq] = nutate(julians) | |
|
166 | ||
|
167 | Modification History | |
|
168 | -------------------- | |
|
169 | Written by W.Landsman (Goddard/HSTX), June 1996. | |
|
170 | Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009. | |
|
171 | """ | |
|
172 | ||
|
173 | jd = numpy.atleast_1d(jd) | |
|
174 | ||
|
175 | # Form time in Julian centuries from 1900 | |
|
176 | t = (jd - 2451545.0)/36525.0 | |
|
177 | ||
|
178 | # Mean elongation of the moon | |
|
179 | coeff1 = numpy.array([1/189474.0,-0.0019142,445267.111480,297.85036]) | |
|
180 | d = numpy.poly1d(coeff1) | |
|
181 | d = d(t)*Misc_Routines.CoFactors.d2r | |
|
182 | d = self.cirrange(d,rad=1) | |
|
183 | ||
|
184 | # Sun's mean elongation | |
|
185 | coeff2 = numpy.array([-1./3e5,-0.0001603,35999.050340,357.52772]) | |
|
186 | m = numpy.poly1d(coeff2) | |
|
187 | m = m(t)*Misc_Routines.CoFactors.d2r | |
|
188 | m = self.cirrange(m,rad=1) | |
|
189 | ||
|
190 | # Moon's mean elongation | |
|
191 | coeff3 = numpy.array([1.0/5.625e4,0.0086972,477198.867398,134.96298]) | |
|
192 | mprime = numpy.poly1d(coeff3) | |
|
193 | mprime = mprime(t)*Misc_Routines.CoFactors.d2r | |
|
194 | mprime = self.cirrange(mprime,rad=1) | |
|
195 | ||
|
196 | # Moon's argument of latitude | |
|
197 | coeff4 = numpy.array([-1.0/3.27270e5,-0.0036825,483202.017538,93.27191]) | |
|
198 | f = numpy.poly1d(coeff4) | |
|
199 | f = f(t)*Misc_Routines.CoFactors.d2r | |
|
200 | f = self.cirrange(f,rad=1) | |
|
201 | ||
|
202 | # Longitude fo the ascending node of the Moon's mean orbit on the ecliptic, measu- | |
|
203 | # red from the mean equinox of the date. | |
|
204 | coeff5 = numpy.array([1.0/4.5e5,0.0020708,-1934.136261,125.04452]) | |
|
205 | omega = numpy.poly1d(coeff5) | |
|
206 | omega = omega(t)*Misc_Routines.CoFactors.d2r | |
|
207 | omega = self.cirrange(omega,rad=1) | |
|
208 | ||
|
209 | d_lng = numpy.array([0,-2,0,0,0,0,-2,0,0,-2,-2,-2,0,2,0,2,0,0,-2,0,2,0,0,-2,0,-2,0,0,\ | |
|
210 | 2,-2,0,-2,0,0,2,2,0,-2,0,2,2,-2,-2,2,2,0,-2,-2,0,-2,-2,0,-1,-2,1,0,0,-1,0,\ | |
|
211 | 0,2,0,2]) | |
|
212 | ||
|
213 | m_lng = numpy.array([0,0,0,0,1,0,1,0,0,-1]) | |
|
214 | m_lng = numpy.append(m_lng,numpy.zeros(17)) | |
|
215 | m_lng = numpy.append(m_lng,numpy.array([2,0,2,1,0,-1,0,0,0,1,1,-1,0,0,0,0,0,0,-1,-1,0,0,\ | |
|
216 | 0,1,0,0,1,0,0,0,-1,1,-1,-1,0,-1])) | |
|
217 | ||
|
218 | mp_lng = numpy.array([0,0,0,0,0,1,0,0,1,0,1,0,-1,0,1,-1,-1,1,2,-2,0,2,2,1,0,0, -1, 0,\ | |
|
219 | -1,0,0,1,0,2,-1,1,0,1,0,0,1,2,1,-2,0,1,0,0,2,2,0,1,1,0,0,1,-2,1,1,1,-1,3,0]) | |
|
220 | ||
|
221 | f_lng = numpy.array([0,2,2,0,0,0,2,2,2,2,0,2,2,0,0,2,0,2,0,2,2,2,0,2,2,2,2,0,0,2,0,0,\ | |
|
222 | 0,-2,2,2,2,0,2,2,0,2,2,0,0,0,2,0,2,0,2,-2,0,0,0,2,2,0,0,2,2,2,2]) | |
|
223 | ||
|
224 | om_lng = numpy.array([1,2,2,2,0,0,2,1,2,2,0,1,2,0,1,2,1,1,0,1,2,2,0,2,0,0,1,0,1,2,1, \ | |
|
225 | 1,1,0,1,2,2,0,2,1,0,2,1,1,1,0,1,1,1,1,1,0,0,0,0,0,2,0,0,2,2,2,2]) | |
|
226 | ||
|
227 | sin_lng = numpy.array([-171996,-13187,-2274,2062,1426,712,-517,-386,-301, 217, -158, \ | |
|
228 | 129,123,63,63,-59,-58,-51,48,46,-38,-31,29,29,26,-22,21,17,16,-16,-15,-13,\ | |
|
229 | -12,11,-10,-8,7,-7,-7,-7,6,6,6,-6,-6,5,-5,-5,-5,4,4,4,-4,-4,-4,3,-3,-3,-3,\ | |
|
230 | -3,-3,-3,-3]) | |
|
231 | ||
|
232 | sdelt = numpy.array([-174.2,-1.6,-0.2,0.2,-3.4,0.1,1.2,-0.4,0,-0.5,0, 0.1, 0, 0, 0.1,\ | |
|
233 | 0,-0.1]) | |
|
234 | sdelt = numpy.append(sdelt,numpy.zeros(10)) | |
|
235 | sdelt = numpy.append(sdelt,numpy.array([-0.1, 0, 0.1])) | |
|
236 | sdelt = numpy.append(sdelt,numpy.zeros(33)) | |
|
237 | ||
|
238 | cos_lng = numpy.array([92025,5736,977,-895,54,-7,224,200,129,-95,0,-70,-53,0,-33,26, \ | |
|
239 | 32,27,0,-24,16,13,0,-12,0,0,-10,0,-8,7,9,7,6,0,5,3,-3,0,3,3,0,-3,-3,3,3,0,\ | |
|
240 | 3,3,3]) | |
|
241 | cos_lng = numpy.append(cos_lng,numpy.zeros(14)) | |
|
242 | ||
|
243 | cdelt = numpy.array([8.9,-3.1,-0.5,0.5,-0.1,0.0,-0.6,0.0,-0.1,0.3]) | |
|
244 | cdelt = numpy.append(cdelt,numpy.zeros(53)) | |
|
245 | ||
|
246 | # Sum the periodic terms. | |
|
247 | n = numpy.size(jd) | |
|
248 | nut_long = numpy.zeros(n) | |
|
249 | nut_obliq = numpy.zeros(n) | |
|
250 | ||
|
251 | d_lng = d_lng.reshape(numpy.size(d_lng),1) | |
|
252 | d = d.reshape(numpy.size(d),1) | |
|
253 | matrix_d_lng = numpy.dot(d_lng,d.transpose()) | |
|
254 | ||
|
255 | m_lng = m_lng.reshape(numpy.size(m_lng),1) | |
|
256 | m = m.reshape(numpy.size(m),1) | |
|
257 | matrix_m_lng = numpy.dot(m_lng,m.transpose()) | |
|
258 | ||
|
259 | mp_lng = mp_lng.reshape(numpy.size(mp_lng),1) | |
|
260 | mprime = mprime.reshape(numpy.size(mprime),1) | |
|
261 | matrix_mp_lng = numpy.dot(mp_lng,mprime.transpose()) | |
|
262 | ||
|
263 | f_lng = f_lng.reshape(numpy.size(f_lng),1) | |
|
264 | f = f.reshape(numpy.size(f),1) | |
|
265 | matrix_f_lng = numpy.dot(f_lng,f.transpose()) | |
|
266 | ||
|
267 | om_lng = om_lng.reshape(numpy.size(om_lng),1) | |
|
268 | omega = omega.reshape(numpy.size(omega),1) | |
|
269 | matrix_om_lng = numpy.dot(om_lng,omega.transpose()) | |
|
270 | ||
|
271 | arg = matrix_d_lng + matrix_m_lng + matrix_mp_lng + matrix_f_lng + matrix_om_lng | |
|
272 | ||
|
273 | sarg = numpy.sin(arg) | |
|
274 | carg = numpy.cos(arg) | |
|
275 | ||
|
276 | for ii in numpy.arange(n): | |
|
277 | nut_long[ii] = 0.0001*numpy.sum((sdelt*t[ii] + sin_lng)*sarg[:,ii]) | |
|
278 | nut_obliq[ii] = 0.0001*numpy.sum((cdelt*t[ii] + cos_lng)*carg[:,ii]) | |
|
279 | ||
|
280 | if numpy.size(jd)==1: | |
|
281 | nut_long = nut_long[0] | |
|
282 | nut_obliq = nut_obliq[0] | |
|
283 | ||
|
284 | return nut_long, nut_obliq | |
|
285 | ||
|
286 | def co_aberration(self,jd,ra,dec): | |
|
287 | """ | |
|
288 | co_aberration calculates changes to Ra and Dec due to "the effect of aberration". | |
|
289 | ||
|
290 | Parameters | |
|
291 | ---------- | |
|
292 | jd = Julian Date (Scalar or vector). | |
|
293 | ra = A scalar o vector giving the Right Ascention of interest. | |
|
294 | dec = A scalar o vector giving the Declination of interest. | |
|
295 | ||
|
296 | Return | |
|
297 | ------ | |
|
298 | d_ra = The correction to right ascension due to aberration (must be added to ra to | |
|
299 | get the correct value). | |
|
300 | d_dec = The correction to declination due to aberration (must be added to the dec | |
|
301 | to get the correct value). | |
|
302 | eps = True obliquity of the ecliptic (in radians). | |
|
303 | ||
|
304 | Examples | |
|
305 | -------- | |
|
306 | >> Julian = 2462088.7 | |
|
307 | >> Ra = 41.547213 | |
|
308 | >> Dec = 49.348483 | |
|
309 | >> [d_ra,d_dec,eps] = co_aberration(julian,Ra,Dec) | |
|
310 | >> print d_ra, d_dec, eps | |
|
311 | [ 30.04441796] [ 6.69837858] [ 0.40904059] | |
|
312 | ||
|
313 | Modification history | |
|
314 | -------------------- | |
|
315 | Written by Chris O'Dell , Univ. of Wisconsin, June 2002. | |
|
316 | Converted to Python by Freddy R. Galindo, ROJ, 27 September 2009. | |
|
317 | """ | |
|
318 | ||
|
319 | # Julian centuries from J2000 of jd. | |
|
320 | T = (jd - 2451545.0)/36525.0 | |
|
321 | ||
|
322 | # Getting obliquity of ecliptic | |
|
323 | njd = numpy.size(jd) | |
|
324 | jd = numpy.atleast_1d(jd) | |
|
325 | ra = numpy.atleast_1d(ra) | |
|
326 | dec = numpy.atleast_1d(dec) | |
|
327 | ||
|
328 | d_psi = numpy.zeros(njd) | |
|
329 | d_epsilon = d_psi | |
|
330 | for ii in numpy.arange(njd): | |
|
331 | [dp,de] = self.nutate(jd[ii]) | |
|
332 | d_psi[ii] = dp | |
|
333 | d_epsilon[ii] = de | |
|
334 | ||
|
335 | coeff = 23 + 26/60. + 21.488/3600. | |
|
336 | eps0 = coeff*3600. - 46.8150*T - 0.00059*T**2. + 0.001813*T**3. | |
|
337 | # True obliquity of the ecliptic in radians | |
|
338 | eps = (eps0 + d_epsilon)/3600*Misc_Routines.CoFactors.d2r | |
|
339 | ||
|
340 | celestialbodies = CelestialBodies() | |
|
341 | [sunra,sundec,sunlon,sunobliq] = celestialbodies.sunpos(jd) | |
|
342 | ||
|
343 | # Earth's orbital eccentricity | |
|
344 | e = 0.016708634 - 0.000042037*T - 0.0000001267*T**2. | |
|
345 | ||
|
346 | # longitude of perihelion, in degrees | |
|
347 | pi = 102.93735 + 1.71946*T + 0.00046*T**2 | |
|
348 | ||
|
349 | # Constant of aberration, in arcseconds | |
|
350 | k = 20.49552 | |
|
351 | ||
|
352 | cd = numpy.cos(dec*Misc_Routines.CoFactors.d2r) ; sd = numpy.sin(dec*Misc_Routines.CoFactors.d2r) | |
|
353 | ce = numpy.cos(eps) ; te = numpy.tan(eps) | |
|
354 | cp = numpy.cos(pi*Misc_Routines.CoFactors.d2r) ; sp = numpy.sin(pi*Misc_Routines.CoFactors.d2r) | |
|
355 | cs = numpy.cos(sunlon*Misc_Routines.CoFactors.d2r) ; ss = numpy.sin(sunlon*Misc_Routines.CoFactors.d2r) | |
|
356 | ca = numpy.cos(ra*Misc_Routines.CoFactors.d2r) ; sa = numpy.sin(ra*Misc_Routines.CoFactors.d2r) | |
|
357 | ||
|
358 | term1 = (ca*cs*ce + sa*ss)/cd | |
|
359 | term2 = (ca*cp*ce + sa*sp)/cd | |
|
360 | term3 = (cs*ce*(te*cd - sa*sd) + ca*sd*ss) | |
|
361 | term4 = (cp*ce*(te*cd - sa*sd) + ca*sd*sp) | |
|
362 | ||
|
363 | d_ra = -k*term1 + e*k*term2 | |
|
364 | d_dec = -k*term3 + e*k*term4 | |
|
365 | ||
|
366 | return d_ra, d_dec, eps | |
|
367 | ||
|
368 | def precess(self,ra,dec,equinox1=None,equinox2=None,FK4=0,rad=0): | |
|
369 | """ | |
|
370 | precess coordinates from EQUINOX1 to EQUINOX2 | |
|
371 | ||
|
372 | Parameters | |
|
373 | ----------- | |
|
374 | ra = A scalar o vector giving the Right Ascention of interest. | |
|
375 | dec = A scalar o vector giving the Declination of interest. | |
|
376 | equinox1 = Original equinox of coordinates, numeric scalar. If omitted, the __Pre- | |
|
377 | cess will query for equinox1 and equinox2. | |
|
378 | equinox2 = Original equinox of coordinates. | |
|
379 | FK4 = If this keyword is set and non-zero, the FK4 (B1950) system will be used | |
|
380 | otherwise FK5 (J2000) will be used instead. | |
|
381 | rad = If this keyword is set and non-zero, then the input and output RAD and DEC | |
|
382 | vectors are in radian rather than degree. | |
|
383 | ||
|
384 | Return | |
|
385 | ------ | |
|
386 | ra = Right ascension after precession (scalar or vector) in degrees, unless the rad | |
|
387 | keyword is set. | |
|
388 | dec = Declination after precession (scalar or vector) in degrees, unless the rad | |
|
389 | keyword is set. | |
|
390 | ||
|
391 | Examples | |
|
392 | -------- | |
|
393 | >> Ra = 329.88772 | |
|
394 | >> Dec = -56.992515 | |
|
395 | >> [p_ra,p_dec] = precess(Ra,Dec,1950,1975,FK4=1) | |
|
396 | >> print p_ra, p_dec | |
|
397 | [ 330.31442971] [-56.87186154] | |
|
398 | ||
|
399 | Modification history | |
|
400 | -------------------- | |
|
401 | Written by Wayne Landsman, STI Corporation, August 1986. | |
|
402 | Converted to Python by Freddy R. Galindo, ROJ, 27 September 2009. | |
|
403 | """ | |
|
404 | ||
|
405 | npts = numpy.size(ra) | |
|
406 | ra = numpy.atleast_1d(ra) | |
|
407 | dec = numpy.atleast_1d(dec) | |
|
408 | ||
|
409 | if rad==0: | |
|
410 | ra_rad = ra*Misc_Routines.CoFactors.d2r | |
|
411 | dec_rad = dec*Misc_Routines.CoFactors.d2r | |
|
412 | else: | |
|
413 | ra_rad = ra | |
|
414 | dec_rad = dec | |
|
415 | ||
|
416 | x = numpy.zeros((npts,3)) | |
|
417 | x[:,0] = numpy.cos(dec_rad)*numpy.cos(ra_rad) | |
|
418 | x[:,1] = numpy.cos(dec_rad)*numpy.sin(ra_rad) | |
|
419 | x[:,2] = numpy.sin(dec_rad) | |
|
420 | ||
|
421 | # Use premat function to get precession matrix from equinox1 to equinox2 | |
|
422 | r = self.premat(equinox1,equinox2,FK4) | |
|
423 | ||
|
424 | x2 = numpy.dot(r,x.transpose()) | |
|
425 | ||
|
426 | ra_rad = numpy.arctan2(x2[1,:],x2[0,:]) | |
|
427 | dec_rad = numpy.arcsin(x2[2,:]) | |
|
428 | ||
|
429 | if rad==0: | |
|
430 | ra = ra_rad/Misc_Routines.CoFactors.d2r | |
|
431 | ra = ra + (ra<0)*360. | |
|
432 | dec = dec_rad/Misc_Routines.CoFactors.d2r | |
|
433 | else: | |
|
434 | ra = ra_rad | |
|
435 | ra = ra + (ra<0)*numpy.pi*2. | |
|
436 | dec = dec_rad | |
|
437 | ||
|
438 | return ra, dec | |
|
439 | ||
|
440 | def premat(self,equinox1,equinox2,FK4=0): | |
|
441 | """ | |
|
442 | premat returns the precession matrix needed to go from EQUINOX1 to EQUINOX2. | |
|
443 | ||
|
444 | Parameters | |
|
445 | ---------- | |
|
446 | equinox1 = Original equinox of coordinates, numeric scalar. | |
|
447 | equinox2 = Equinox of precessed coordinates. | |
|
448 | FK4 = If this keyword is set and non-zero, the FK4 (B1950) system precession angles | |
|
449 | are used to compute the precession matrix. The default is to use FK5 (J2000) pre- | |
|
450 | cession angles. | |
|
451 | ||
|
452 | Return | |
|
453 | ------ | |
|
454 | r = Precession matrix, used to precess equatorial rectangular coordinates. | |
|
455 | ||
|
456 | Examples | |
|
457 | -------- | |
|
458 | >> matrix = premat(1950.0,1975.0,FK4=1) | |
|
459 | >> print matrix | |
|
460 | [[ 9.99981438e-01 -5.58774959e-03 -2.42908517e-03] | |
|
461 | [ 5.58774959e-03 9.99984388e-01 -6.78691471e-06] | |
|
462 | [ 2.42908517e-03 -6.78633095e-06 9.99997050e-01]] | |
|
463 | ||
|
464 | Modification history | |
|
465 | -------------------- | |
|
466 | Written by Wayne Landsman, HSTX Corporation, June 1994. | |
|
467 | Converted to Python by Freddy R. Galindo, ROJ, 27 September 2009. | |
|
468 | """ | |
|
469 | ||
|
470 | t = 0.001*(equinox2 - equinox1) | |
|
471 | ||
|
472 | if FK4==0: | |
|
473 | st=0.001*(equinox1 - 2000.) | |
|
474 | # Computing 3 rotation angles. | |
|
475 | A=Misc_Routines.CoFactors.s2r*t*(23062.181+st*(139.656+0.0139*st)+t*(30.188-0.344*st+17.998*t)) | |
|
476 | B=Misc_Routines.CoFactors.s2r*t*t*(79.280+0.410*st+0.205*t)+A | |
|
477 | C=Misc_Routines.CoFactors.s2r*t*(20043.109-st*(85.33+0.217*st)+ t*(-42.665-0.217*st-41.833*t)) | |
|
478 | else: | |
|
479 | st=0.001*(equinox1 - 1900) | |
|
480 | # Computing 3 rotation angles | |
|
481 | A=Misc_Routines.CoFactors.s2r*t*(23042.53+st*(139.75+0.06*st)+t*(30.23-0.27*st+18.0*t)) | |
|
482 | B=Misc_Routines.CoFactors.s2r*t*t*(79.27+0.66*st+0.32*t)+A | |
|
483 | C=Misc_Routines.CoFactors.s2r*t*(20046.85-st*(85.33+0.37*st)+t*(-42.67-0.37*st-41.8*t)) | |
|
484 | ||
|
485 | sina = numpy.sin(A); sinb = numpy.sin(B); sinc = numpy.sin(C) | |
|
486 | cosa = numpy.cos(A); cosb = numpy.cos(B); cosc = numpy.cos(C) | |
|
487 | ||
|
488 | r = numpy.zeros((3,3)) | |
|
489 | r[:,0] = numpy.array([cosa*cosb*cosc-sina*sinb,sina*cosb+cosa*sinb*cosc,cosa*sinc]) | |
|
490 | r[:,1] = numpy.array([-cosa*sinb-sina*cosb*cosc,cosa*cosb-sina*sinb*cosc,-sina*sinc]) | |
|
491 | r[:,2] = numpy.array([-cosb*sinc,-sinb*sinc,cosc]) | |
|
492 | ||
|
493 | return r | |
|
494 | ||
|
495 | def cirrange(self,angle,rad=0): | |
|
496 | """ | |
|
497 | cirrange forces an angle into the range 0<= angle < 360. | |
|
498 | ||
|
499 | Parameters | |
|
500 | ---------- | |
|
501 | angle = The angle to modify, in degrees. Can be scalar or vector. | |
|
502 | rad = Set to 1 if the angle is specified in radians rather than degrees. It is for- | |
|
503 | ced into the range 0 <= angle < 2 PI | |
|
504 | ||
|
505 | Return | |
|
506 | ------ | |
|
507 | angle = The angle after the modification. | |
|
508 | ||
|
509 | Example | |
|
510 | ------- | |
|
511 | >> angle = cirrange(numpy.array([420,400,361])) | |
|
512 | >> print angle | |
|
513 | >> [60, 40, 1] | |
|
514 | ||
|
515 | Modification History | |
|
516 | -------------------- | |
|
517 | Written by Michael R. Greason, Hughes STX, 10 February 1994. | |
|
518 | Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009. | |
|
519 | """ | |
|
520 | ||
|
521 | angle = numpy.atleast_1d(angle) | |
|
522 | ||
|
523 | if rad==1: | |
|
524 | cnst = numpy.pi*2. | |
|
525 | elif rad==0: | |
|
526 | cnst = 360. | |
|
527 | ||
|
528 | # Deal with the lower limit. | |
|
529 | angle = angle % cnst | |
|
530 | ||
|
531 | # Deal with negative values, if way | |
|
532 | neg = numpy.where(angle<0.0) | |
|
533 | if neg[0].size>0: angle[neg] = angle[neg] + cnst | |
|
534 | ||
|
535 | return angle | |
|
536 | ||
|
537 | ||
|
538 | class CelestialBodies(EquatorialCorrections): | |
|
539 | def __init__(self): | |
|
540 | """ | |
|
541 | CelestialBodies class creates a object to call methods of celestial bodies location. | |
|
542 | ||
|
543 | Modification History | |
|
544 | -------------------- | |
|
545 | Converted to Object-oriented Programming by Freddy Galindo, ROJ, 27 September 2009. | |
|
546 | """ | |
|
547 | ||
|
548 | EquatorialCorrections.__init__(self) | |
|
549 | ||
|
550 | def sunpos(self,jd,rad=0): | |
|
551 | """ | |
|
552 | sunpos method computes the RA and Dec of the Sun at a given date. | |
|
553 | ||
|
554 | Parameters | |
|
555 | ---------- | |
|
556 | jd = The julian date of the day (and time), scalar or vector. | |
|
557 | rad = If this keyword is set and non-zero, then the input and output RAD and DEC | |
|
558 | vectors are in radian rather than degree. | |
|
559 | ||
|
560 | Return | |
|
561 | ------ | |
|
562 | ra = The right ascension of the sun at that date in degrees. | |
|
563 | dec = The declination of the sun at that date in degrees. | |
|
564 | elong = Ecliptic longitude of the sun at that date in degrees. | |
|
565 | obliquity = The declination of the sun at that date in degrees. | |
|
566 | ||
|
567 | Examples | |
|
568 | -------- | |
|
569 | >> jd = 2466880 | |
|
570 | >> [ra,dec,elong,obliquity] = sunpos(jd) | |
|
571 | >> print ra, dec, elong, obliquity | |
|
572 | [ 275.53499556] [-23.33840558] [ 275.08917968] [ 23.43596165] | |
|
573 | ||
|
574 | >> [ra,dec,elong,obliquity] = sunpos(jd,rad=1) | |
|
575 | >> print ra, dec, elong, obliquity | |
|
576 | [ 4.80899288] [-0.40733202] [ 4.80121192] [ 0.40903469] | |
|
577 | ||
|
578 | >> jd = 2450449.5 + numpy.arange(365) | |
|
579 | >> [ra,dec,elong,obliquity] = sunpos(jd) | |
|
580 | ||
|
581 | Modification history | |
|
582 | -------------------- | |
|
583 | Written by Micheal R. Greason, STX Corporation, 28 October 1988. | |
|
584 | Converted to Python by Freddy R. Galindo, ROJ, 27 September 2009. | |
|
585 | """ | |
|
586 | ||
|
587 | jd = numpy.atleast_1d(jd) | |
|
588 | ||
|
589 | # Form time in Julian centuries from 1900. | |
|
590 | t = (jd -2415020.0)/36525.0 | |
|
591 | ||
|
592 | # Form sun's mean longitude | |
|
593 | l = (279.696678+((36000.768925*t) % 360.0))*3600.0 | |
|
594 | ||
|
595 | # Allow for ellipticity of the orbit (equation of centre) using the Earth's mean | |
|
596 | # anomoly ME | |
|
597 | me = 358.475844 + ((35999.049750*t) % 360.0) | |
|
598 | ellcor = (6910.1 - 17.2*t)*numpy.sin(me*Misc_Routines.CoFactors.d2r) + 72.3*numpy.sin(2.0*me*Misc_Routines.CoFactors.d2r) | |
|
599 | l = l + ellcor | |
|
600 | ||
|
601 | # Allow for the Venus perturbations using the mean anomaly of Venus MV | |
|
602 | mv = 212.603219 + ((58517.803875*t) % 360.0) | |
|
603 | vencorr = 4.8*numpy.cos((299.1017 + mv - me)*Misc_Routines.CoFactors.d2r) + \ | |
|
604 | 5.5*numpy.cos((148.3133 + 2.0*mv - 2.0*me )*Misc_Routines.CoFactors.d2r) + \ | |
|
605 | 2.5*numpy.cos((315.9433 + 2.0*mv - 3.0*me )*Misc_Routines.CoFactors.d2r) + \ | |
|
606 | 1.6*numpy.cos((345.2533 + 3.0*mv - 4.0*me )*Misc_Routines.CoFactors.d2r) + \ | |
|
607 | 1.0*numpy.cos((318.15 + 3.0*mv - 5.0*me )*Misc_Routines.CoFactors.d2r) | |
|
608 | l = l + vencorr | |
|
609 | ||
|
610 | # Allow for the Mars perturbations using the mean anomaly of Mars MM | |
|
611 | mm = 319.529425 + ((19139.858500*t) % 360.0) | |
|
612 | marscorr = 2.0*numpy.cos((343.8883 - 2.0*mm + 2.0*me)*Misc_Routines.CoFactors.d2r ) + \ | |
|
613 | 1.8*numpy.cos((200.4017 - 2.0*mm + me)*Misc_Routines.CoFactors.d2r) | |
|
614 | l = l + marscorr | |
|
615 | ||
|
616 | # Allow for the Jupiter perturbations using the mean anomaly of Jupiter MJ | |
|
617 | mj = 225.328328 + ((3034.6920239*t) % 360.0) | |
|
618 | jupcorr = 7.2*numpy.cos((179.5317 - mj + me )*Misc_Routines.CoFactors.d2r) + \ | |
|
619 | 2.6*numpy.cos((263.2167 - mj)*Misc_Routines.CoFactors.d2r) + \ | |
|
620 | 2.7*numpy.cos((87.1450 - 2.0*mj + 2.0*me)*Misc_Routines.CoFactors.d2r) + \ | |
|
621 | 1.6*numpy.cos((109.4933 - 2.0*mj + me)*Misc_Routines.CoFactors.d2r) | |
|
622 | l = l + jupcorr | |
|
623 | ||
|
624 | # Allow for Moons perturbations using mean elongation of the Moon from the Sun D | |
|
625 | d = 350.7376814 + ((445267.11422*t) % 360.0) | |
|
626 | mooncorr = 6.5*numpy.sin(d*Misc_Routines.CoFactors.d2r) | |
|
627 | l = l + mooncorr | |
|
628 | ||
|
629 | # Allow for long period terms | |
|
630 | longterm = + 6.4*numpy.sin((231.19 + 20.20*t)*Misc_Routines.CoFactors.d2r) | |
|
631 | l = l + longterm | |
|
632 | l = (l + 2592000.0) % 1296000.0 | |
|
633 | longmed = l/3600.0 | |
|
634 | ||
|
635 | # Allow for Aberration | |
|
636 | l = l - 20.5 | |
|
637 | ||
|
638 | # Allow for Nutation using the longitude of the Moons mean node OMEGA | |
|
639 | omega = 259.183275 - ((1934.142008*t) % 360.0) | |
|
640 | l = l - 17.2*numpy.sin(omega*Misc_Routines.CoFactors.d2r) | |
|
641 | ||
|
642 | # Form the True Obliquity | |
|
643 | oblt = 23.452294 - 0.0130125*t + (9.2*numpy.cos(omega*Misc_Routines.CoFactors.d2r))/3600.0 | |
|
644 | ||
|
645 | # Form Right Ascension and Declination | |
|
646 | l = l/3600.0 | |
|
647 | ra = numpy.arctan2((numpy.sin(l*Misc_Routines.CoFactors.d2r)*numpy.cos(oblt*Misc_Routines.CoFactors.d2r)),numpy.cos(l*Misc_Routines.CoFactors.d2r)) | |
|
648 | ||
|
649 | neg = numpy.where(ra < 0.0) | |
|
650 | if neg[0].size > 0: ra[neg] = ra[neg] + 2.0*numpy.pi | |
|
651 | ||
|
652 | dec = numpy.arcsin(numpy.sin(l*Misc_Routines.CoFactors.d2r)*numpy.sin(oblt*Misc_Routines.CoFactors.d2r)) | |
|
653 | ||
|
654 | if rad==1: | |
|
655 | oblt = oblt*Misc_Routines.CoFactors.d2r | |
|
656 | longmed = longmed*Misc_Routines.CoFactors.d2r | |
|
657 | else: | |
|
658 | ra = ra/Misc_Routines.CoFactors.d2r | |
|
659 | dec = dec/Misc_Routines.CoFactors.d2r | |
|
660 | ||
|
661 | return ra, dec, longmed, oblt | |
|
662 | ||
|
663 | def moonpos(self,jd,rad=0): | |
|
664 | """ | |
|
665 | moonpos method computes the RA and Dec of the Moon at specified Julian date(s). | |
|
666 | ||
|
667 | Parameters | |
|
668 | ---------- | |
|
669 | jd = The julian date of the day (and time), scalar or vector. | |
|
670 | rad = If this keyword is set and non-zero, then the input and output RAD and DEC | |
|
671 | vectors are in radian rather than degree. | |
|
672 | ||
|
673 | Return | |
|
674 | ------ | |
|
675 | ra = The right ascension of the sun at that date in degrees. | |
|
676 | dec = The declination of the sun at that date in degrees. | |
|
677 | dist = The Earth-moon distance in kilometers (between the center of the Earth and | |
|
678 | the center of the moon). | |
|
679 | geolon = Apparent longitude of the moon in degrees, referred to the ecliptic of the | |
|
680 | specified date(s). | |
|
681 | geolat = Apparent latitude the moon in degrees, referred to the ecliptic of the | |
|
682 | specified date(s). | |
|
683 | ||
|
684 | Examples | |
|
685 | -------- | |
|
686 | >> jd = 2448724.5 | |
|
687 | >> [ra,dec,dist,geolon,geolat] = sunpos(jd) | |
|
688 | >> print ra, dec, dist, geolon, geolat | |
|
689 | [ 134.68846855] [ 13.76836663] [ 368409.68481613] [ 133.16726428] [-3.22912642] | |
|
690 | ||
|
691 | >> [ra,dec,dist,geolon, geolat] = sunpos(jd,rad=1) | |
|
692 | >> print ra, dec, dist, geolon, geolat | |
|
693 | [ 2.35075724] [ 0.24030333] [ 368409.68481613] [ 2.32420722] [-0.05635889] | |
|
694 | ||
|
695 | >> jd = 2450449.5 + numpy.arange(365) | |
|
696 | >> [ra,dec,dist,geolon, geolat] = sunpos(jd) | |
|
697 | ||
|
698 | Modification history | |
|
699 | -------------------- | |
|
700 | Written by Micheal R. Greason, STX Corporation, 31 October 1988. | |
|
701 | Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009. | |
|
702 | """ | |
|
703 | ||
|
704 | jd = numpy.atleast_1d(jd) | |
|
705 | ||
|
706 | # Form time in Julian centuries from 1900. | |
|
707 | t = (jd - 2451545.0)/36525.0 | |
|
708 | ||
|
709 | d_lng = numpy.array([0,2,2,0,0,0,2,2,2,2,0,1,0,2,0,0,4,0,4,2,2,1,1,2,2,4,2,0,2,2,1,2,\ | |
|
710 | 0,0,2,2,2,4,0,3,2,4,0,2,2,2,4,0,4,1,2,0,1,3,4,2,0,1,2,2]) | |
|
711 | ||
|
712 | m_lng = numpy.array([0,0,0,0,1,0,0,-1,0,-1,1,0,1,0,0,0,0,0,0,1,1,0,1,-1,0,0,0,1,0,-1,\ | |
|
713 | 0,-2,1,2,-2,0,0,-1,0,0,1,-1,2,2,1,-1,0,0,-1,0,1,0,1,0,0,-1,2,1,0,0]) | |
|
714 | ||
|
715 | mp_lng = numpy.array([1,-1,0,2,0,0,-2,-1,1,0,-1,0,1,0,1,1,-1,3,-2,-1,0,-1,0,1,2,0,-3,\ | |
|
716 | -2,-1,-2,1,0,2,0,-1,1,0,-1,2,-1,1,-2,-1,-1,-2,0,1,4,0,-2,0,2,1,-2,-3,2,1,-1,3,-1]) | |
|
717 | ||
|
718 | f_lng = numpy.array([0,0,0,0,0,2,0,0,0,0,0,0,0,-2,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,\ | |
|
719 | 0,0,0,0,-2,2,0,2,0,0,0,0,0,0,-2,0,0,0,0,-2,-2,0,0,0,0,0,0,0,-2]) | |
|
720 | ||
|
721 | sin_lng = numpy.array([6288774,1274027,658314,213618,-185116,-114332,58793,57066,\ | |
|
722 | 53322,45758,-40923,-34720,-30383,15327,-12528,10980,10675,10034,8548,-7888,\ | |
|
723 | -6766,-5163,4987,4036,3994,3861,3665,-2689,-2602,2390,-2348,2236,-2120,-2069,\ | |
|
724 | 2048,-1773,-1595,1215,-1110,-892,-810,759,-713,-700,691,596,549,537,520,-487,\ | |
|
725 | -399,-381,351,-340,330,327,-323,299,294,0.0]) | |
|
726 | ||
|
727 | cos_lng = numpy.array([-20905355,-3699111,-2955968,-569925,48888,-3149,246158,-152138,\ | |
|
728 | -170733,-204586,-129620,108743,104755,10321,0,79661,-34782,-23210,-21636,24208,\ | |
|
729 | 30824,-8379,-16675,-12831,-10445,-11650,14403,-7003,0,10056,6322, -9884,5751,0,\ | |
|
730 | -4950,4130,0,-3958,0,3258,2616,-1897,-2117,2354,0,0,-1423,-1117,-1571,-1739,0, \ | |
|
731 | -4421,0,0,0,0,1165,0,0,8752.0]) | |
|
732 | ||
|
733 | d_lat = numpy.array([0,0,0,2,2,2,2,0,2,0,2,2,2,2,2,2,2,0,4,0,0,0,1,0,0,0,1,0,4,4,0,4,\ | |
|
734 | 2,2,2,2,0,2,2,2,2,4,2,2,0,2,1,1,0,2,1,2,0,4,4,1,4,1,4,2]) | |
|
735 | ||
|
736 | m_lat = numpy.array([0,0,0,0,0,0,0,0,0,0,-1,0,0,1,-1,-1,-1,1,0,1,0,1,0,1,1,1,0,0,0,0,\ | |
|
737 | 0,0,0,0,-1,0,0,0,0,1,1,0,-1,-2,0,1,1,1,1,1,0,-1,1,0,-1,0,0,0,-1,-2]) | |
|
738 | ||
|
739 | mp_lat = numpy.array([0,1,1,0,-1,-1,0,2,1,2,0,-2,1,0,-1,0,-1,-1,-1,0,0,-1,0,1,1,0,0,\ | |
|
740 | 3,0,-1,1,-2,0,2,1,-2,3,2,-3,-1,0,0,1,0,1,1,0,0,-2,-1,1,-2,2,-2,-1,1,1,-1,0,0]) | |
|
741 | ||
|
742 | f_lat = numpy.array([1,1,-1,-1,1,-1,1,1,-1,-1,-1,-1,1,-1,1,1,-1,-1,-1,1,3,1,1,1,-1,\ | |
|
743 | -1,-1,1,-1,1,-3,1,-3,-1,-1,1,-1,1,-1,1,1,1,1,-1,3,-1,-1,1,-1,-1,1,-1,1,-1,-1, \ | |
|
744 | -1,-1,-1,-1,1]) | |
|
745 | ||
|
746 | sin_lat = numpy.array([5128122,280602,277693,173237,55413,46271, 32573, 17198, 9266, \ | |
|
747 | 8822,8216,4324,4200,-3359,2463,2211,2065,-1870,1828,-1794, -1749, -1565, -1491, \ | |
|
748 | -1475,-1410,-1344,-1335,1107,1021,833,777,671,607,596,491,-451,439,422,421,-366,\ | |
|
749 | -351,331,315,302,-283,-229,223,223,-220,-220,-185,181,-177,176, 166, -164, 132, \ | |
|
750 | -119,115,107.0]) | |
|
751 | ||
|
752 | # Mean longitude of the moon refered to mean equinox of the date. | |
|
753 | coeff0 = numpy.array([-1./6.5194e7,1./538841.,-0.0015786,481267.88123421,218.3164477]) | |
|
754 | lprimed = numpy.poly1d(coeff0) | |
|
755 | lprimed = lprimed(t) | |
|
756 | lprimed = self.cirrange(lprimed,rad=0) | |
|
757 | lprime = lprimed*Misc_Routines.CoFactors.d2r | |
|
758 | ||
|
759 | # Mean elongation of the moon | |
|
760 | coeff1 = numpy.array([-1./1.13065e8,1./545868.,-0.0018819,445267.1114034,297.8501921]) | |
|
761 | d = numpy.poly1d(coeff1) | |
|
762 | d = d(t)*Misc_Routines.CoFactors.d2r | |
|
763 | d = self.cirrange(d,rad=1) | |
|
764 | ||
|
765 | # Sun's mean anomaly | |
|
766 | coeff2 = numpy.array([1.0/2.449e7,-0.0001536,35999.0502909,357.5291092]) | |
|
767 | M = numpy.poly1d(coeff2) | |
|
768 | M = M(t)*Misc_Routines.CoFactors.d2r | |
|
769 | M = self.cirrange(M,rad=1) | |
|
770 | ||
|
771 | # Moon's mean anomaly | |
|
772 | coeff3 = numpy.array([-1.0/1.4712e7,1.0/6.9699e4,0.0087414,477198.8675055,134.9633964]) | |
|
773 | Mprime = numpy.poly1d(coeff3) | |
|
774 | Mprime = Mprime(t)*Misc_Routines.CoFactors.d2r | |
|
775 | Mprime = self.cirrange(Mprime,rad=1) | |
|
776 | ||
|
777 | # Moon's argument of latitude | |
|
778 | coeff4 = numpy.array([1.0/8.6331e8,-1.0/3.526e7,-0.0036539,483202.0175233,93.2720950]) | |
|
779 | F = numpy.poly1d(coeff4) | |
|
780 | F = F(t)*Misc_Routines.CoFactors.d2r | |
|
781 | F = self.cirrange(F,rad=1) | |
|
782 | ||
|
783 | # Eccentricity of Earth's orbit around the sun | |
|
784 | e = 1 - 0.002516*t - 7.4e-6*(t**2.) | |
|
785 | e2 = e**2. | |
|
786 | ||
|
787 | ecorr1 = numpy.where((numpy.abs(m_lng))==1) | |
|
788 | ecorr2 = numpy.where((numpy.abs(m_lat))==1) | |
|
789 | ecorr3 = numpy.where((numpy.abs(m_lng))==2) | |
|
790 | ecorr4 = numpy.where((numpy.abs(m_lat))==2) | |
|
791 | ||
|
792 | # Additional arguments. | |
|
793 | A1 = (119.75 + 131.849*t)*Misc_Routines.CoFactors.d2r | |
|
794 | A2 = (53.09 + 479264.290*t)*Misc_Routines.CoFactors.d2r | |
|
795 | A3 = (313.45 + 481266.484*t)*Misc_Routines.CoFactors.d2r | |
|
796 | suml_add = 3958.*numpy.sin(A1) + 1962.*numpy.sin(lprime - F) + 318*numpy.sin(A2) | |
|
797 | sumb_add = -2235.*numpy.sin(lprime) + 382.*numpy.sin(A3) + 175.*numpy.sin(A1-F) + \ | |
|
798 | 175.*numpy.sin(A1 + F) + 127.*numpy.sin(lprime - Mprime) - 115.*numpy.sin(lprime + Mprime) | |
|
799 | ||
|
800 | # Sum the periodic terms | |
|
801 | geolon = numpy.zeros(jd.size) | |
|
802 | geolat = numpy.zeros(jd.size) | |
|
803 | dist = numpy.zeros(jd.size) | |
|
804 | ||
|
805 | for i in numpy.arange(jd.size): | |
|
806 | sinlng = sin_lng | |
|
807 | coslng = cos_lng | |
|
808 | sinlat = sin_lat | |
|
809 | ||
|
810 | sinlng[ecorr1] = e[i]*sinlng[ecorr1] | |
|
811 | coslng[ecorr1] = e[i]*coslng[ecorr1] | |
|
812 | sinlat[ecorr2] = e[i]*sinlat[ecorr2] | |
|
813 | sinlng[ecorr3] = e2[i]*sinlng[ecorr3] | |
|
814 | coslng[ecorr3] = e2[i]*coslng[ecorr3] | |
|
815 | sinlat[ecorr4] = e2[i]*sinlat[ecorr4] | |
|
816 | ||
|
817 | arg = d_lng*d[i] + m_lng*M[i] + mp_lng*Mprime[i] + f_lng*F[i] | |
|
818 | geolon[i] = lprimed[i] + (numpy.sum(sinlng*numpy.sin(arg)) + suml_add[i] )/1.e6 | |
|
819 | dist[i] = 385000.56 + numpy.sum(coslng*numpy.cos(arg))/1.e3 | |
|
820 | arg = d_lat*d[i] + m_lat*M[i] + mp_lat*Mprime[i] + f_lat*F[i] | |
|
821 | geolat[i] = (numpy.sum(sinlat*numpy.sin(arg)) + sumb_add[i])/1.e6 | |
|
822 | ||
|
823 | [nlon, elon] = self.nutate(jd) | |
|
824 | geolon = geolon + nlon/3.6e3 | |
|
825 | geolon = self.cirrange(geolon,rad=0) | |
|
826 | lamb = geolon*Misc_Routines.CoFactors.d2r | |
|
827 | beta = geolat*Misc_Routines.CoFactors.d2r | |
|
828 | ||
|
829 | # Find mean obliquity and convert lamb, beta to RA, Dec | |
|
830 | c = numpy.array([2.45,5.79,27.87,7.12,-39.05,-249.67,-51.38,1999.25,-1.55,-4680.93, \ | |
|
831 | 21.448]) | |
|
832 | junk = numpy.poly1d(c); | |
|
833 | epsilon = 23. + (26./60.) + (junk(t/1.e2)/3600.) | |
|
834 | # True obliquity in radians | |
|
835 | eps = (epsilon + elon/3600. )*Misc_Routines.CoFactors.d2r | |
|
836 | ||
|
837 | ra = numpy.arctan2(numpy.sin(lamb)*numpy.cos(eps)-numpy.tan(beta)*numpy.sin(eps),numpy.cos(lamb)) | |
|
838 | ra = self.cirrange(ra,rad=1) | |
|
839 | ||
|
840 | dec = numpy.arcsin(numpy.sin(beta)*numpy.cos(eps) + numpy.cos(beta)*numpy.sin(eps)*numpy.sin(lamb)) | |
|
841 | ||
|
842 | if rad==1: | |
|
843 | geolon = lamb | |
|
844 | geolat = beta | |
|
845 | else: | |
|
846 | ra = ra/Misc_Routines.CoFactors.d2r | |
|
847 | dec = dec/Misc_Routines.CoFactors.d2r | |
|
848 | ||
|
849 | return ra, dec, dist, geolon, geolat | |
|
850 | ||
|
851 | def hydrapos(self): | |
|
852 | """ | |
|
853 | hydrapos method returns RA and Dec provided by Bill Coles (Oct 2003). | |
|
854 | ||
|
855 | Parameters | |
|
856 | ---------- | |
|
857 | None | |
|
858 | ||
|
859 | Return | |
|
860 | ------ | |
|
861 | ra = The right ascension of the sun at that date in degrees. | |
|
862 | dec = The declination of the sun at that date in degrees. | |
|
863 | Examples | |
|
864 | -------- | |
|
865 | >> [ra,dec] = hydrapos() | |
|
866 | >> print ra, dec | |
|
867 | 139.45 -12.0833333333 | |
|
868 | ||
|
869 | Modification history | |
|
870 | -------------------- | |
|
871 | Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009. | |
|
872 | """ | |
|
873 | ||
|
874 | ra = (9. + 17.8/60.)*15. | |
|
875 | dec = -(12. + 5./60.) | |
|
876 | ||
|
877 | return ra, dec | |
|
878 | ||
|
879 | ||
|
880 | def skynoise_jro(self,dec_cut=-11.95,filename='skynoise_jro.dat',filepath=None): | |
|
881 | """ | |
|
882 | hydrapos returns RA and Dec provided by Bill Coles (Oct 2003). | |
|
883 | ||
|
884 | Parameters | |
|
885 | ---------- | |
|
886 | dec_cut = A scalar giving the declination to get a cut of the skynoise over Jica- | |
|
887 | marca. The default value is -11.95. | |
|
888 | filename = A string to specify name the skynoise file. The default value is skynoi- | |
|
889 | se_jro.dat | |
|
890 | ||
|
891 | Return | |
|
892 | ------ | |
|
893 | maxra = The maximum right ascension to the declination used to get a cut. | |
|
894 | ra = The right ascension. | |
|
895 | Examples | |
|
896 | -------- | |
|
897 | >> [maxra,ra] = skynoise_jro() | |
|
898 | >> print maxra, ra | |
|
899 | 139.45 -12.0833333333 | |
|
900 | ||
|
901 | Modification history | |
|
902 | -------------------- | |
|
903 | Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009. | |
|
904 | """ | |
|
905 | ||
|
906 | if filepath==None: | |
|
907 | filepath = '/app/utils/' | |
|
908 | ||
|
909 | f = open(os.path.join(filepath,filename),'rb') | |
|
910 | ||
|
911 | # Reading SkyNoise Power (lineal scale) | |
|
912 | ha_sky = numpy.fromfile(f,numpy.dtype([('var','<f4')]),480*20) | |
|
913 | ha_sky = ha_sky['var'].reshape(20,480).transpose() | |
|
914 | ||
|
915 | dec_sky = numpy.fromfile(f,numpy.dtype([('var','<f4')]),480*20) | |
|
916 | dec_sky = dec_sky['var'].reshape((20,480)).transpose() | |
|
917 | ||
|
918 | tmp_sky = numpy.fromfile(f,numpy.dtype([('var','<f4')]),480*20) | |
|
919 | tmp_sky = tmp_sky['var'].reshape((20,480)).transpose() | |
|
920 | ||
|
921 | f.close() | |
|
922 | ||
|
923 | nha = 480 | |
|
924 | tmp_cut = numpy.zeros(nha) | |
|
925 | for iha in numpy.arange(nha): | |
|
926 | tck = scipy.interpolate.splrep(dec_sky[iha,:],tmp_sky[iha,:],s=0) | |
|
927 | tmp_cut[iha] = scipy.interpolate.splev(dec_cut,tck,der=0) | |
|
928 | ||
|
929 | ptr = numpy.nanargmax(tmp_cut) | |
|
930 | ||
|
931 | maxra = ha_sky[ptr,0] | |
|
932 | ra = ha_sky[:,0] | |
|
933 | ||
|
934 | return maxra, ra | |
|
935 | ||
|
936 | def skyNoise(self,jd,ut=-5.0,longitude=-76.87,filename='galaxy.txt',filepath=None): | |
|
937 | """ | |
|
938 | hydrapos returns RA and Dec provided by Bill Coles (Oct 2003). | |
|
939 | ||
|
940 | Parameters | |
|
941 | ---------- | |
|
942 | jd = The julian date of the day (and time), scalar or vector. | |
|
943 | ||
|
944 | dec_cut = A scalar giving the declination to get a cut of the skynoise over Jica- | |
|
945 | marca. The default value is -11.95. | |
|
946 | filename = A string to specify name the skynoise file. The default value is skynoi- | |
|
947 | se_jro.dat | |
|
948 | ||
|
949 | Return | |
|
950 | ------ | |
|
951 | maxra = The maximum right ascension to the declination used to get a cut. | |
|
952 | ra = The right ascension. | |
|
953 | ||
|
954 | Examples | |
|
955 | -------- | |
|
956 | >> [maxra,ra] = skynoise_jro() | |
|
957 | >> print maxra, ra | |
|
958 | 139.45 -12.0833333333 | |
|
959 | ||
|
960 | Modification history | |
|
961 | -------------------- | |
|
962 | Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009. | |
|
963 | """ | |
|
964 | ||
|
965 | # Defining date to compute SkyNoise. | |
|
966 | [year, month, dom, hour, mis, secs] = TimeTools.Julian(jd).change2time() | |
|
967 | is_dom = (month==9) & (dom==21) | |
|
968 | if is_dom: | |
|
969 | tmp = jd | |
|
970 | jd = TimeTools.Time(year,9,22).change2julian() | |
|
971 | dom = 22 | |
|
972 | ||
|
973 | # Reading SkyNoise | |
|
974 | if filepath==None:filepath='./resource' | |
|
975 | f = open(os.path.join(filepath,filename)) | |
|
976 | ||
|
977 | lines = f.read() | |
|
978 | f.close() | |
|
979 | ||
|
980 | nlines = 99 | |
|
981 | lines = lines.split('\n') | |
|
982 | data = numpy.zeros((2,nlines))*numpy.float32(0.) | |
|
983 | for ii in numpy.arange(nlines): | |
|
984 | line = numpy.array([lines[ii][0:6],lines[ii][6:]]) | |
|
985 | data[:,ii] = numpy.float32(line) | |
|
986 | ||
|
987 | # Getting SkyNoise to the date desired. | |
|
988 | otime = data[0,:]*60.0 | |
|
989 | opowr = data[1,:] | |
|
990 | ||
|
991 | hour = numpy.array([0,23]); | |
|
992 | mins = numpy.array([0,59]); | |
|
993 | secs = numpy.array([0,59]); | |
|
994 | LTrange = TimeTools.Time(year,month,dom,hour,mins,secs).change2julday() | |
|
995 | LTtime = LTrange[0] + numpy.arange(1440)*((LTrange[1] - LTrange[0])/(1440.-1)) | |
|
996 | lst = TimeTools.Julian(LTtime + (-3600.*ut/86400.)).change2lst() | |
|
997 | ||
|
998 | ipowr = lst*0.0 | |
|
999 | # Interpolating using scipy (inside max and min "x") | |
|
1000 | otime = otime/3600. | |
|
1001 | val = numpy.where((lst>numpy.min(otime)) & (lst<numpy.max(otime))); val = val[0] | |
|
1002 | tck = scipy.interpolate.interp1d(otime,opowr) | |
|
1003 | ipowr[val] = tck(lst[val]) | |
|
1004 | ||
|
1005 | # Extrapolating above maximum time data (23.75). | |
|
1006 | uval = numpy.where(lst>numpy.max(otime)) | |
|
1007 | if uval[0].size>0: | |
|
1008 | ii = numpy.min(uval[0]) | |
|
1009 | m = (ipowr[ii-1] - ipowr[ii-2])/(lst[ii-1] - lst[ii-2]) | |
|
1010 | b = ipowr[ii-1] - m*lst[ii-1] | |
|
1011 | ipowr[uval] = m*lst[uval] + b | |
|
1012 | ||
|
1013 | if is_dom: | |
|
1014 | lst = numpy.roll(lst,4) | |
|
1015 | ipowr = numpy.roll(ipowr,4) | |
|
1016 | ||
|
1017 | new_lst = numpy.int32(lst*3600.) | |
|
1018 | new_pow = ipowr | |
|
1019 | ||
|
1020 | return ipowr, LTtime, lst | |
|
1021 | ||
|
1022 | ||
|
1023 | class AltAz(EquatorialCorrections): | |
|
1024 | def __init__(self,alt,az,jd,lat=-11.95,lon=-76.8667,WS=0,altitude=500,nutate_=0,precess_=0,\ | |
|
1025 | aberration_=0,B1950=0): | |
|
1026 | """ | |
|
1027 | The AltAz class creates an object which represents the target position in horizontal | |
|
1028 | coordinates (alt-az) and allows to convert (using the methods) from this coordinate | |
|
1029 | system to others (e.g. Equatorial). | |
|
1030 | ||
|
1031 | Parameters | |
|
1032 | ---------- | |
|
1033 | alt = Altitude in degrees. Scalar or vector. | |
|
1034 | az = Azimuth angle in degrees (measured EAST from NORTH, but see keyword WS). Sca- | |
|
1035 | lar or vector. | |
|
1036 | jd = Julian date. Scalar or vector. | |
|
1037 | lat = North geodetic latitude of location in degrees. The default value is -11.95. | |
|
1038 | lon = East longitude of location in degrees. The default value is -76.8667. | |
|
1039 | WS = Set this to 1 to get the azimuth measured westward from south. | |
|
1040 | altitude = The altitude of the observing location, in meters. The default 500. | |
|
1041 | nutate_ = Set this to 1 to force nutation, 0 for no nutation. | |
|
1042 | precess_ = Set this to 1 to force precession, 0 for no precession. | |
|
1043 | aberration_ = Set this to 1 to force aberration correction, 0 for no correction. | |
|
1044 | B1950 = Set this if your RA and DEC are specified in B1950, FK4 coordinates (ins- | |
|
1045 | tead of J2000, FK5) | |
|
1046 | ||
|
1047 | Modification History | |
|
1048 | -------------------- | |
|
1049 | Converted to Object-oriented Programming by Freddy Galindo, ROJ, 26 September 2009. | |
|
1050 | """ | |
|
1051 | ||
|
1052 | EquatorialCorrections.__init__(self) | |
|
1053 | ||
|
1054 | self.alt = numpy.atleast_1d(alt) | |
|
1055 | self.az = numpy.atleast_1d(az) | |
|
1056 | self.jd = numpy.atleast_1d(jd) | |
|
1057 | self.lat = lat | |
|
1058 | self.lon = lon | |
|
1059 | self.WS = WS | |
|
1060 | self.altitude = altitude | |
|
1061 | ||
|
1062 | self.nutate_ = nutate_ | |
|
1063 | self.aberration_ = aberration_ | |
|
1064 | self.precess_ = precess_ | |
|
1065 | self.B1950 = B1950 | |
|
1066 | ||
|
1067 | def change2equatorial(self): | |
|
1068 | """ | |
|
1069 | change2equatorial method converts horizon (Alt-Az) coordinates to equatorial coordi- | |
|
1070 | nates (ra-dec). | |
|
1071 | ||
|
1072 | Return | |
|
1073 | ------ | |
|
1074 | ra = Right ascension of object (J2000) in degrees (FK5). Scalar or vector. | |
|
1075 | dec = Declination of object (J2000), in degrees (FK5). Scalar or vector. | |
|
1076 | ha = Hour angle in degrees. | |
|
1077 | ||
|
1078 | Example | |
|
1079 | ------- | |
|
1080 | >> alt = 88.5401 | |
|
1081 | >> az = -128.990 | |
|
1082 | >> jd = 2452640.5 | |
|
1083 | >> ObjAltAz = AltAz(alt,az,jd) | |
|
1084 | >> [ra, dec, ha] = ObjAltAz.change2equatorial() | |
|
1085 | >> print ra, dec, ha | |
|
1086 | [ 22.20280632] [-12.86610025] [ 1.1638927] | |
|
1087 | ||
|
1088 | Modification History | |
|
1089 | -------------------- | |
|
1090 | Written Chris O'Dell Univ. of Wisconsin-Madison, May 2002. | |
|
1091 | Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009. | |
|
1092 | """ | |
|
1093 | ||
|
1094 | az = self.az | |
|
1095 | alt = self.alt | |
|
1096 | if self.WS>0:az = az -180. | |
|
1097 | ra_tmp = numpy.zeros(numpy.size(self.jd)) + 45. | |
|
1098 | dec_tmp = numpy.zeros(numpy.size(self.jd)) + 45. | |
|
1099 | [dra1,ddec1,eps,d_psi,d_eps] = self.co_nutate(self.jd,ra_tmp, dec_tmp) | |
|
1100 | ||
|
1101 | # Getting local mean sidereal time (lmst) | |
|
1102 | lmst = TimeTools.Julian(self.jd[0]).change2lst() | |
|
1103 | lmst = lmst*Misc_Routines.CoFactors.h2d | |
|
1104 | # Getting local apparent sidereal time (last) | |
|
1105 | last = lmst + d_psi*numpy.cos(eps)/3600. | |
|
1106 | ||
|
1107 | # Now do the spherical trig to get APPARENT hour angle and declination (Degrees). | |
|
1108 | [ha, dec] = self.change2HaDec() | |
|
1109 | ||
|
1110 | # Finding Right Ascension (in degrees, from 0 to 360.) | |
|
1111 | ra = (last - ha + 360.) % 360. | |
|
1112 | ||
|
1113 | # Calculate NUTATION and ABERRATION Correction to Ra-Dec | |
|
1114 | [dra1, ddec1,eps,d_psi,d_eps] = self.co_nutate(self.jd,ra,dec) | |
|
1115 | [dra2,ddec2,eps] = self.co_aberration(self.jd,ra,dec) | |
|
1116 | ||
|
1117 | # Make Nutation and Aberration correction (if wanted) | |
|
1118 | ra = ra - (dra1*self.nutate_ + dra2*self.aberration_)/3600. | |
|
1119 | dec = dec - (ddec1*self.nutate_ + ddec2*self.aberration_)/3600. | |
|
1120 | ||
|
1121 | # Computing current equinox | |
|
1122 | j_now = (self.jd - 2451545.)/365.25 + 2000 | |
|
1123 | ||
|
1124 | # Precess coordinates to current date | |
|
1125 | if self.precess_==1: | |
|
1126 | njd = numpy.size(self.jd) | |
|
1127 | for ii in numpy.arange(njd): | |
|
1128 | ra_i = ra[ii] | |
|
1129 | dec_i = dec[ii] | |
|
1130 | now = j_now[ii] | |
|
1131 | ||
|
1132 | if self.B1950==1: | |
|
1133 | [ra_i,dec_i] = self.precess(ra_i,dec_i,now,1950.,FK4=1) | |
|
1134 | elif self.B1950==0: | |
|
1135 | [ra_i,dec_i] = self.precess(ra_i,dec_i,now,2000.,FK4=0) | |
|
1136 | ||
|
1137 | ra[ii] = ra_i | |
|
1138 | dec[ii] = dec_i | |
|
1139 | ||
|
1140 | return ra, dec, ha | |
|
1141 | ||
|
1142 | def change2HaDec(self): | |
|
1143 | """ | |
|
1144 | change2HaDec method converts from horizon (Alt-Az) coordinates to hour angle and de- | |
|
1145 | clination. | |
|
1146 | ||
|
1147 | Return | |
|
1148 | ------ | |
|
1149 | ha = The local apparent hour angle, in degrees. The hour angle is the time that ri- | |
|
1150 | ght ascension of 0 hours crosses the local meridian. It is unambiguisoly defined. | |
|
1151 | dec = The local apparent declination, in degrees. | |
|
1152 | ||
|
1153 | Example | |
|
1154 | ------- | |
|
1155 | >> alt = 88.5401 | |
|
1156 | >> az = -128.990 | |
|
1157 | >> jd = 2452640.5 | |
|
1158 | >> ObjAltAz = AltAz(alt,az,jd) | |
|
1159 | >> [ha, dec] = ObjAltAz.change2HaDec() | |
|
1160 | >> print ha, dec | |
|
1161 | [ 1.1638927] [-12.86610025] | |
|
1162 | ||
|
1163 | Modification History | |
|
1164 | -------------------- | |
|
1165 | Written Chris O'Dell Univ. of Wisconsin-Madison, May 2002. | |
|
1166 | Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009. | |
|
1167 | """ | |
|
1168 | ||
|
1169 | alt_r = numpy.atleast_1d(self.alt*Misc_Routines.CoFactors.d2r) | |
|
1170 | az_r = numpy.atleast_1d(self.az*Misc_Routines.CoFactors.d2r) | |
|
1171 | lat_r = numpy.atleast_1d(self.lat*Misc_Routines.CoFactors.d2r) | |
|
1172 | ||
|
1173 | # Find local hour angle (in degrees, from 0 to 360.) | |
|
1174 | y_ha = -1*numpy.sin(az_r)*numpy.cos(alt_r) | |
|
1175 | x_ha = -1*numpy.cos(az_r)*numpy.sin(lat_r)*numpy.cos(alt_r) + numpy.sin(alt_r)*numpy.cos(lat_r) | |
|
1176 | ||
|
1177 | ha = numpy.arctan2(y_ha,x_ha) | |
|
1178 | ha = ha/Misc_Routines.CoFactors.d2r | |
|
1179 | ||
|
1180 | w = numpy.where(ha<0.) | |
|
1181 | if w[0].size>0:ha[w] = ha[w] + 360. | |
|
1182 | ha = ha % 360. | |
|
1183 | ||
|
1184 | # Find declination (positive if north of celestial equatorial, negative if south) | |
|
1185 | sindec = numpy.sin(lat_r)*numpy.sin(alt_r) + numpy.cos(lat_r)*numpy.cos(alt_r)*numpy.cos(az_r) | |
|
1186 | dec = numpy.arcsin(sindec)/Misc_Routines.CoFactors.d2r | |
|
1187 | ||
|
1188 | return ha, dec | |
|
1189 | ||
|
1190 | ||
|
1191 | class Equatorial(EquatorialCorrections): | |
|
1192 | def __init__(self,ra,dec,jd,lat=-11.95,lon=-76.8667,WS=0,altitude=500,nutate_=0,precess_=0,\ | |
|
1193 | aberration_=0,B1950=0): | |
|
1194 | """ | |
|
1195 | The Equatorial class creates an object which represents the target position in equa- | |
|
1196 | torial coordinates (ha-dec) and allows to convert (using the class methods) from | |
|
1197 | this coordinate system to others (e.g. AltAz). | |
|
1198 | ||
|
1199 | Parameters | |
|
1200 | ---------- | |
|
1201 | ra = Right ascension of object (J2000) in degrees (FK5). Scalar or vector. | |
|
1202 | dec = Declination of object (J2000), in degrees (FK5). Scalar or vector. | |
|
1203 | jd = Julian date. Scalar or vector. | |
|
1204 | lat = North geodetic latitude of location in degrees. The default value is -11.95. | |
|
1205 | lon = East longitude of location in degrees. The default value is -76.8667. | |
|
1206 | WS = Set this to 1 to get the azimuth measured westward from south. | |
|
1207 | altitude = The altitude of the observing location, in meters. The default 500. | |
|
1208 | nutate = Set this to 1 to force nutation, 0 for no nutation. | |
|
1209 | precess = Set this to 1 to force precession, 0 for no precession. | |
|
1210 | aberration = Set this to 1 to force aberration correction, 0 for no correction. | |
|
1211 | B1950 = Set this if your RA and DEC are specified in B1950, FK4 coordinates (ins- | |
|
1212 | tead of J2000, FK5) | |
|
1213 | ||
|
1214 | Modification History | |
|
1215 | -------------------- | |
|
1216 | Converted to Object-oriented Programming by Freddy Galindo, ROJ, 29 September 2009. | |
|
1217 | """ | |
|
1218 | ||
|
1219 | EquatorialCorrections.__init__(self) | |
|
1220 | ||
|
1221 | self.ra = numpy.atleast_1d(ra) | |
|
1222 | self.dec = numpy.atleast_1d(dec) | |
|
1223 | self.jd = numpy.atleast_1d(jd) | |
|
1224 | self.lat = lat | |
|
1225 | self.lon = lon | |
|
1226 | self.WS = WS | |
|
1227 | self.altitude = altitude | |
|
1228 | ||
|
1229 | self.nutate_ = nutate_ | |
|
1230 | self.aberration_ = aberration_ | |
|
1231 | self.precess_ = precess_ | |
|
1232 | self.B1950 = B1950 | |
|
1233 | ||
|
1234 | def change2AltAz(self): | |
|
1235 | """ | |
|
1236 | change2AltAz method converts from equatorial coordinates (ha-dec) to horizon coordi- | |
|
1237 | nates (alt-az). | |
|
1238 | ||
|
1239 | Return | |
|
1240 | ------ | |
|
1241 | alt = Altitude in degrees. Scalar or vector. | |
|
1242 | az = Azimuth angle in degrees (measured EAST from NORTH, but see keyword WS). Sca- | |
|
1243 | lar or vector. | |
|
1244 | ha = Hour angle in degrees. | |
|
1245 | ||
|
1246 | Example | |
|
1247 | ------- | |
|
1248 | >> ra = 43.370609 | |
|
1249 | >> dec = -28.0000 | |
|
1250 | >> jd = 2452640.5 | |
|
1251 | >> ObjEq = Equatorial(ra,dec,jd) | |
|
1252 | >> [alt, az, ha] = ObjEq.change2AltAz() | |
|
1253 | >> print alt, az, ha | |
|
1254 | [ 65.3546497] [ 133.58753124] [ 339.99609002] | |
|
1255 | ||
|
1256 | Modification History | |
|
1257 | -------------------- | |
|
1258 | Written Chris O'Dell Univ. of Wisconsin-Madison. May 2002 | |
|
1259 | Converted to Python by Freddy R. Galindo, ROJ, 29 September 2009. | |
|
1260 | """ | |
|
1261 | ||
|
1262 | ra = self.ra | |
|
1263 | dec = self.dec | |
|
1264 | ||
|
1265 | # Computing current equinox | |
|
1266 | j_now = (self.jd - 2451545.)/365.25 + 2000 | |
|
1267 | ||
|
1268 | # Precess coordinates to current date | |
|
1269 | if self.precess_==1: | |
|
1270 | njd = numpy.size(self.jd) | |
|
1271 | for ii in numpy.arange(njd): | |
|
1272 | ra_i = ra[ii] | |
|
1273 | dec_i = dec[ii] | |
|
1274 | now = j_now[ii] | |
|
1275 | ||
|
1276 | if self.B1950==1: | |
|
1277 | [ra_i,dec_i] = self.precess(ra_i,dec_i,now,1950.,FK4=1) | |
|
1278 | elif self.B1950==0: | |
|
1279 | [ra_i,dec_i] = self.precess(ra_i,dec_i,now,2000.,FK4=0) | |
|
1280 | ||
|
1281 | ra[ii] = ra_i | |
|
1282 | dec[ii] = dec_i | |
|
1283 | ||
|
1284 | # Calculate NUTATION and ABERRATION Correction to Ra-Dec | |
|
1285 | [dra1, ddec1,eps,d_psi,d_eps] = self.co_nutate(self.jd,ra,dec) | |
|
1286 | [dra2,ddec2,eps] = self.co_aberration(self.jd,ra,dec) | |
|
1287 | ||
|
1288 | # Make Nutation and Aberration correction (if wanted) | |
|
1289 | ra = ra + (dra1*self.nutate_ + dra2*self.aberration_)/3600. | |
|
1290 | dec = dec + (ddec1*self.nutate_ + ddec2*self.aberration_)/3600. | |
|
1291 | ||
|
1292 | # Getting local mean sidereal time (lmst) | |
|
1293 | lmst = TimeTools.Julian(self.jd).change2lst() | |
|
1294 | ||
|
1295 | lmst = lmst*Misc_Routines.CoFactors.h2d | |
|
1296 | # Getting local apparent sidereal time (last) | |
|
1297 | last = lmst + d_psi*numpy.cos(eps)/3600. | |
|
1298 | ||
|
1299 | # Finding Hour Angle (in degrees, from 0 to 360.) | |
|
1300 | ha = last - ra | |
|
1301 | w = numpy.where(ha<0.) | |
|
1302 | if w[0].size>0:ha[w] = ha[w] + 360. | |
|
1303 | ha = ha % 360. | |
|
1304 | ||
|
1305 | # Now do the spherical trig to get APPARENT hour angle and declination (Degrees). | |
|
1306 | [alt, az] = self.HaDec2AltAz(ha,dec) | |
|
1307 | ||
|
1308 | return alt, az, ha | |
|
1309 | ||
|
1310 | def HaDec2AltAz(self,ha,dec): | |
|
1311 | """ | |
|
1312 | HaDec2AltAz convert hour angle and declination (ha-dec) to horizon coords (alt-az). | |
|
1313 | ||
|
1314 | Parameters | |
|
1315 | ---------- | |
|
1316 | ha = The local apparent hour angle, in DEGREES, scalar or vector. | |
|
1317 | dec = The local apparent declination, in DEGREES, scalar or vector. | |
|
1318 | ||
|
1319 | Return | |
|
1320 | ------ | |
|
1321 | alt = Altitude in degrees. Scalar or vector. | |
|
1322 | az = Azimuth angle in degrees (measured EAST from NORTH, but see keyword WS). Sca- | |
|
1323 | lar or vector. | |
|
1324 | ||
|
1325 | Modification History | |
|
1326 | -------------------- | |
|
1327 | Written Chris O'Dell Univ. of Wisconsin-Madison, May 2002. | |
|
1328 | Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009. | |
|
1329 | """ | |
|
1330 | ||
|
1331 | sh = numpy.sin(ha*Misc_Routines.CoFactors.d2r) ; ch = numpy.cos(ha*Misc_Routines.CoFactors.d2r) | |
|
1332 | sd = numpy.sin(dec*Misc_Routines.CoFactors.d2r) ; cd = numpy.cos(dec*Misc_Routines.CoFactors.d2r) | |
|
1333 | sl = numpy.sin(self.lat*Misc_Routines.CoFactors.d2r) ; cl = numpy.cos(self.lat*Misc_Routines.CoFactors.d2r) | |
|
1334 | ||
|
1335 | x = -1*ch*cd*sl + sd*cl | |
|
1336 | y = -1*sh*cd | |
|
1337 | z = ch*cd*cl + sd*sl | |
|
1338 | r = numpy.sqrt(x**2. + y**2.) | |
|
1339 | ||
|
1340 | az = numpy.arctan2(y,x)/Misc_Routines.CoFactors.d2r | |
|
1341 | alt = numpy.arctan2(z,r)/Misc_Routines.CoFactors.d2r | |
|
1342 | ||
|
1343 | # correct for negative az. | |
|
1344 | w = numpy.where(az<0.) | |
|
1345 | if w[0].size>0:az[w] = az[w] + 360. | |
|
1346 | ||
|
1347 | # Convert az to West from South, if desired | |
|
1348 | if self.WS==1: az = (az + 180.) % 360. | |
|
1349 | ||
|
1350 | return alt, az | |
|
1351 | ||
|
1352 | ||
|
1353 | class Geodetic(): | |
|
1354 | def __init__(self,lat=-11.95,alt=0): | |
|
1355 | """ | |
|
1356 | The Geodetic class creates an object which represents the real position on earth of | |
|
1357 | a target (Geodetic Coordinates: lat-alt) and allows to convert (using the class me- | |
|
1358 | thods) from this coordinate system to others (e.g. geocentric). | |
|
1359 | ||
|
1360 | Parameters | |
|
1361 | ---------- | |
|
1362 | lat = Geodetic latitude of location in degrees. The default value is -11.95. | |
|
1363 | ||
|
1364 | alt = Geodetic altitude (km). The default value is 0. | |
|
1365 | ||
|
1366 | Modification History | |
|
1367 | -------------------- | |
|
1368 | Converted to Object-oriented Programming by Freddy R. Galindo, ROJ, 02 October 2009. | |
|
1369 | """ | |
|
1370 | ||
|
1371 | self.lat = numpy.atleast_1d(lat) | |
|
1372 | self.alt = numpy.atleast_1d(alt) | |
|
1373 | ||
|
1374 | self.a = 6378.16 | |
|
1375 | self.ab2 = 1.0067397 | |
|
1376 | self.ep2 = 0.0067397 | |
|
1377 | ||
|
1378 | def change2geocentric(self): | |
|
1379 | """ | |
|
1380 | change2geocentric method converts from Geodetic to Geocentric coordinates. The re- | |
|
1381 | ference geoid is that adopted by the IAU in 1964. | |
|
1382 | ||
|
1383 | Return | |
|
1384 | ------ | |
|
1385 | gclat = Geocentric latitude (in degrees), scalar or vector. | |
|
1386 | gcalt = Geocentric radial distance (km), scalar or vector. | |
|
1387 | ||
|
1388 | Example | |
|
1389 | ------- | |
|
1390 | >> ObjGeoid = Geodetic(lat=-11.95,alt=0) | |
|
1391 | >> [gclat, gcalt] = ObjGeoid.change2geocentric() | |
|
1392 | >> print gclat, gcalt | |
|
1393 | [-11.87227742] [ 6377.25048195] | |
|
1394 | ||
|
1395 | Modification History | |
|
1396 | -------------------- | |
|
1397 | Converted to Python by Freddy R. Galindo, ROJ, 02 October 2009. | |
|
1398 | """ | |
|
1399 | ||
|
1400 | gdl = self.lat*Misc_Routines.CoFactors.d2r | |
|
1401 | slat = numpy.sin(gdl) | |
|
1402 | clat = numpy.cos(gdl) | |
|
1403 | slat2 = slat**2. | |
|
1404 | clat2 = (self.ab2*clat)**2. | |
|
1405 | ||
|
1406 | sbet = slat/numpy.sqrt(slat2 + clat2) | |
|
1407 | sbet2 = (sbet**2.) # < 1 | |
|
1408 | noval = numpy.where(sbet2>1) | |
|
1409 | if noval[0].size>0:sbet2[noval] = 1 | |
|
1410 | cbet = numpy.sqrt(1. - sbet2) | |
|
1411 | ||
|
1412 | rgeoid = self.a/numpy.sqrt(1. + self.ep2*sbet2) | |
|
1413 | ||
|
1414 | x = rgeoid*cbet + self.alt*clat | |
|
1415 | y = rgeoid*sbet + self.alt*slat | |
|
1416 | ||
|
1417 | gcalt = numpy.sqrt(x**2. + y**2.) | |
|
1418 | gclat = numpy.arctan2(y,x)/Misc_Routines.CoFactors.d2r | |
|
1419 | ||
|
1420 | return gclat, gcalt |
@@ -0,0 +1,61 | |||
|
1 | """ | |
|
2 | The module MISC_ROUTINES gathers classes and functions which are useful for daily processing. As an | |
|
3 | example we have conversion factor or universal constants. | |
|
4 | ||
|
5 | MODULES CALLED: | |
|
6 | NUMPY, SYS | |
|
7 | ||
|
8 | MODIFICATION HISTORY: | |
|
9 | Created by Ing. Freddy Galindo (frederickgalindo@gmail.com). ROJ, 21 October 2009. | |
|
10 | """ | |
|
11 | ||
|
12 | import numpy | |
|
13 | import sys | |
|
14 | ||
|
15 | class CoFactors(): | |
|
16 | """ | |
|
17 | CoFactor class used to call pre-defined conversion factor (e.g. degree to radian). The cu- | |
|
18 | The current available factor are: | |
|
19 | ||
|
20 | d2r = degree to radian. | |
|
21 | s2r = seconds to radian?, degree to arcsecond.? | |
|
22 | h2r = hour to radian. | |
|
23 | h2d = hour to degree | |
|
24 | """ | |
|
25 | ||
|
26 | d2r = numpy.pi/180. | |
|
27 | s2r = numpy.pi/(180.*3600.) | |
|
28 | h2r = numpy.pi/12. | |
|
29 | h2d = 15. | |
|
30 | ||
|
31 | ||
|
32 | class Vector: | |
|
33 | """ | |
|
34 | direction = 0 Polar to rectangular; direction=1 rectangular to polar | |
|
35 | """ | |
|
36 | def __init__(self,vect,direction=0): | |
|
37 | nsize = numpy.size(vect) | |
|
38 | if nsize <= 3: | |
|
39 | vect = vect.reshape(1,nsize) | |
|
40 | ||
|
41 | self.vect = vect | |
|
42 | self.dirc = direction | |
|
43 | ||
|
44 | ||
|
45 | ||
|
46 | def Polar2Rect(self): | |
|
47 | if self.dirc == 0: | |
|
48 | jvect = self.vect*numpy.pi/180. | |
|
49 | mmx = numpy.cos(jvect[:,1])*numpy.sin(jvect[:,0]) | |
|
50 | mmy = numpy.cos(jvect[:,1])*numpy.cos(jvect[:,0]) | |
|
51 | mmz = numpy.sin(jvect[:,1]) | |
|
52 | mm = numpy.array([mmx,mmy,mmz]).transpose() | |
|
53 | ||
|
54 | elif self.dirc == 1: | |
|
55 | mm = [numpy.arctan2(self.vect[:,0],self.vect[:,1]),numpy.arcsin(self.vect[:,2])] | |
|
56 | mm = numpy.array(mm)*180./numpy.pi | |
|
57 | ||
|
58 | return mm | |
|
59 | ||
|
60 | ||
|
61 | No newline at end of file |
@@ -0,0 +1,17 | |||
|
1 | attenuation = numpy.array([[[-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], | |
|
2 | [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], | |
|
3 | [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], | |
|
4 | [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], | |
|
5 | [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], | |
|
6 | [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], | |
|
7 | [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], | |
|
8 | [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25]], | |
|
9 | [[21.25,21.25,21.25,21.25,21.25,21.25,21.25,21.25], | |
|
10 | [15.25,15.25,15.25,15.25,15.25,15.25,15.25,15.25], | |
|
11 | [09.25,09.25,09.25,09.25,09.25,09.25,09.25,09.25], | |
|
12 | [03.25,03.25,03.25,03.25,03.25,03.25,03.25,03.25], | |
|
13 | [-03.25,-03.25,-03.25,-03.25,-03.25,-03.25,-03.25,-03.25], | |
|
14 | [-09.25,-09.25,-09.25,-09.25,-09.25,-09.25,-09.25,-09.25], | |
|
15 | [-15.25,-15.25,-15.25,-15.25,-15.25,-15.25,-15.25,-15.25], | |
|
16 | [-21.25,-21.25,-21.25,-21.25,-21.25,-21.25,-21.25,-21.25]]]) | |
|
17 |
|
1 | NO CONTENT: new file 100644, binary diff hidden |
@@ -0,0 +1,34 | |||
|
1 | from numpy import * | |
|
2 | from scipy import optimize | |
|
3 | ||
|
4 | def gaussian(height, center_x, center_y, width_x, width_y): | |
|
5 | """Returns a gaussian function with the given parameters""" | |
|
6 | width_x = float(width_x) | |
|
7 | width_y = float(width_y) | |
|
8 | return lambda x,y: height*exp( | |
|
9 | -(((center_x-x)/width_x)**2+((center_y-y)/width_y)**2)/2) | |
|
10 | ||
|
11 | def moments(data): | |
|
12 | """Returns (height, x, y, width_x, width_y) | |
|
13 | the gaussian parameters of a 2D distribution by calculating its | |
|
14 | moments """ | |
|
15 | total = data.sum() | |
|
16 | X, Y = indices(data.shape) | |
|
17 | x = (X*data).sum()/total | |
|
18 | y = (Y*data).sum()/total | |
|
19 | col = data[:, int(y)] | |
|
20 | width_x = sqrt(abs((arange(col.size)-y)**2*col).sum()/col.sum()) | |
|
21 | row = data[int(x), :] | |
|
22 | width_y = sqrt(abs((arange(row.size)-x)**2*row).sum()/row.sum()) | |
|
23 | height = data.max() | |
|
24 | return height, x, y, width_x, width_y | |
|
25 | ||
|
26 | def fitgaussian(data): | |
|
27 | """Returns (height, x, y, width_x, width_y) | |
|
28 | the gaussian parameters of a 2D distribution found by a fit""" | |
|
29 | params = moments(data) | |
|
30 | errorfunction = lambda p: ravel(gaussian(*p)(*indices(data.shape)) - | |
|
31 | data) | |
|
32 | p, success = optimize.leastsq(errorfunction, params) | |
|
33 | return p | |
|
34 |
|
1 | NO CONTENT: new file 100644, binary diff hidden |
@@ -0,0 +1,196 | |||
|
1 | g/h n m 1900.0 1905.0 1910.0 1915.0 1920.0 1925.0 1930.0 1935.0 1940.0 1945.0 1950.0 1955.0 1960.0 1965.0 1970.0 1975.0 1980.0 1985.0 1990.0 1995.0 2000.0 2005.0 SV | |
|
2 | g 1 0 -31543 -31464 -31354 -31212 -31060 -30926 -30805 -30715 -30654 -30594 -30554 -30500 -30421 -30334 -30220 -30100 -29992 -29873 -29775 -29692 -29619.4 -29556.8 8.8 | |
|
3 | g 1 1 -2298 -2298 -2297 -2306 -2317 -2318 -2316 -2306 -2292 -2285 -2250 -2215 -2169 -2119 -2068 -2013 -1956 -1905 -1848 -1784 -1728.2 -1671.8 10.8 | |
|
4 | h 1 1 5922 5909 5898 5875 5845 5817 5808 5812 5821 5810 5815 5820 5791 5776 5737 5675 5604 5500 5406 5306 5186.1 5080.0 -21.3 | |
|
5 | g 2 0 -677 -728 -769 -802 -839 -893 -951 -1018 -1106 -1244 -1341 -1440 -1555 -1662 -1781 -1902 -1997 -2072 -2131 -2200 -2267.7 -2340.5 -15.0 | |
|
6 | g 2 1 2905 2928 2948 2956 2959 2969 2980 2984 2981 2990 2998 3003 3002 2997 3000 3010 3027 3044 3059 3070 3068.4 3047.0 -6.9 | |
|
7 | h 2 1 -1061 -1086 -1128 -1191 -1259 -1334 -1424 -1520 -1614 -1702 -1810 -1898 -1967 -2016 -2047 -2067 -2129 -2197 -2279 -2366 -2481.6 -2594.9 -23.3 | |
|
8 | g 2 2 924 1041 1176 1309 1407 1471 1517 1550 1566 1578 1576 1581 1590 1594 1611 1632 1663 1687 1686 1681 1670.9 1656.9 -1.0 | |
|
9 | h 2 2 1121 1065 1000 917 823 728 644 586 528 477 381 291 206 114 25 -68 -200 -306 -373 -413 -458.0 -516.7 -14.0 | |
|
10 | g 3 0 1022 1037 1058 1084 1111 1140 1172 1206 1240 1282 1297 1302 1302 1297 1287 1276 1281 1296 1314 1335 1339.6 1335.7 -0.3 | |
|
11 | g 3 1 -1469 -1494 -1524 -1559 -1600 -1645 -1692 -1740 -1790 -1834 -1889 -1944 -1992 -2038 -2091 -2144 -2180 -2208 -2239 -2267 -2288.0 -2305.3 -3.1 | |
|
12 | h 3 1 -330 -357 -389 -421 -445 -462 -480 -494 -499 -499 -476 -462 -414 -404 -366 -333 -336 -310 -284 -262 -227.6 -200.4 5.4 | |
|
13 | g 3 2 1256 1239 1223 1212 1205 1202 1205 1215 1232 1255 1274 1288 1289 1292 1278 1260 1251 1247 1248 1249 1252.1 1246.8 -0.9 | |
|
14 | h 3 2 3 34 62 84 103 119 133 146 163 186 206 216 224 240 251 262 271 284 293 302 293.4 269.3 -6.5 | |
|
15 | g 3 3 572 635 705 778 839 881 907 918 916 913 896 882 878 856 838 830 833 829 802 759 714.5 674.4 -6.8 | |
|
16 | h 3 3 523 480 425 360 293 229 166 101 43 -11 -46 -83 -130 -165 -196 -223 -252 -297 -352 -427 -491.1 -524.5 -2.0 | |
|
17 | g 4 0 876 880 884 887 889 891 896 903 914 944 954 958 957 957 952 946 938 936 939 940 932.3 919.8 -2.5 | |
|
18 | g 4 1 628 643 660 678 695 711 727 744 762 776 792 796 800 804 800 791 782 780 780 780 786.8 798.2 2.8 | |
|
19 | h 4 1 195 203 211 218 220 216 205 188 169 144 136 133 135 148 167 191 212 232 247 262 272.6 281.4 2.0 | |
|
20 | g 4 2 660 653 644 631 616 601 584 565 550 544 528 510 504 479 461 438 398 361 325 290 250.0 211.5 -7.1 | |
|
21 | h 4 2 -69 -77 -90 -109 -134 -163 -195 -226 -252 -276 -278 -274 -278 -269 -266 -265 -257 -249 -240 -236 -231.9 -225.8 1.8 | |
|
22 | g 4 3 -361 -380 -400 -416 -424 -426 -422 -415 -405 -421 -408 -397 -394 -390 -395 -405 -419 -424 -423 -418 -403.0 -379.5 5.9 | |
|
23 | h 4 3 -210 -201 -189 -173 -153 -130 -109 -90 -72 -55 -37 -23 3 13 26 39 53 69 84 97 119.8 145.7 5.6 | |
|
24 | g 4 4 134 146 160 178 199 217 234 249 265 304 303 290 269 252 234 216 199 170 141 122 111.3 100.2 -3.2 | |
|
25 | h 4 4 -75 -65 -55 -51 -57 -70 -90 -114 -141 -178 -210 -230 -255 -269 -279 -288 -297 -297 -299 -306 -303.8 -304.7 0.0 | |
|
26 | g 5 0 -184 -192 -201 -211 -221 -230 -237 -241 -241 -253 -240 -229 -222 -219 -216 -218 -218 -214 -214 -214 -218.8 -227.6 -2.6 | |
|
27 | g 5 1 328 328 327 327 326 326 327 329 334 346 349 360 362 358 359 356 357 355 353 352 351.4 354.4 0.4 | |
|
28 | h 5 1 -210 -193 -172 -148 -122 -96 -72 -51 -33 -12 3 15 16 19 26 31 46 47 46 46 43.8 42.7 0.1 | |
|
29 | g 5 2 264 259 253 245 236 226 218 211 208 194 211 230 242 254 262 264 261 253 245 235 222.3 208.8 -3.0 | |
|
30 | h 5 2 53 56 57 58 58 58 60 64 71 95 103 110 125 128 139 148 150 150 154 165 171.9 179.8 1.8 | |
|
31 | g 5 3 5 -1 -9 -16 -23 -28 -32 -33 -33 -20 -20 -23 -26 -31 -42 -59 -74 -93 -109 -118 -130.4 -136.6 -1.2 | |
|
32 | h 5 3 -33 -32 -33 -34 -38 -44 -53 -64 -75 -67 -87 -98 -117 -126 -139 -152 -151 -154 -153 -143 -133.1 -123.0 2.0 | |
|
33 | g 5 4 -86 -93 -102 -111 -119 -125 -131 -136 -141 -142 -147 -152 -156 -157 -160 -159 -162 -164 -165 -166 -168.6 -168.3 0.2 | |
|
34 | h 5 4 -124 -125 -126 -126 -125 -122 -118 -115 -113 -119 -122 -121 -114 -97 -91 -83 -78 -75 -69 -55 -39.3 -19.5 4.5 | |
|
35 | g 5 5 -16 -26 -38 -51 -62 -69 -74 -76 -76 -82 -76 -69 -63 -62 -56 -49 -48 -46 -36 -17 -12.9 -14.1 -0.6 | |
|
36 | h 5 5 3 11 21 32 43 51 58 64 69 82 80 78 81 81 83 88 92 95 97 107 106.3 103.6 -1.0 | |
|
37 | g 6 0 63 62 62 61 61 61 60 59 57 59 54 47 46 45 43 45 48 53 61 68 72.3 72.9 -0.8 | |
|
38 | g 6 1 61 60 58 57 55 54 53 53 54 57 57 57 58 61 64 66 66 65 65 67 68.2 69.6 0.2 | |
|
39 | h 6 1 -9 -7 -5 -2 0 3 4 4 4 6 -1 -9 -10 -11 -12 -13 -15 -16 -16 -17 -17.4 -20.2 -0.4 | |
|
40 | g 6 2 -11 -11 -11 -10 -10 -9 -9 -8 -7 6 4 3 1 8 15 28 42 51 59 68 74.2 76.6 -0.2 | |
|
41 | h 6 2 83 86 89 93 96 99 102 104 105 100 99 96 99 100 100 99 93 88 82 72 63.7 54.7 -1.9 | |
|
42 | g 6 3 -217 -221 -224 -228 -233 -238 -242 -246 -249 -246 -247 -247 -237 -228 -212 -198 -192 -185 -178 -170 -160.9 -151.1 2.1 | |
|
43 | h 6 3 2 4 5 8 11 14 19 25 33 16 33 48 60 68 72 75 71 69 69 67 65.1 63.7 -0.4 | |
|
44 | g 6 4 -58 -57 -54 -51 -46 -40 -32 -25 -18 -25 -16 -8 -1 4 2 1 4 4 3 -1 -5.9 -15.0 -2.1 | |
|
45 | h 6 4 -35 -32 -29 -26 -22 -18 -16 -15 -15 -9 -12 -16 -20 -32 -37 -41 -43 -48 -52 -58 -61.2 -63.4 -0.4 | |
|
46 | g 6 5 59 57 54 49 44 39 32 25 18 21 12 7 -2 1 3 6 14 16 18 19 16.9 14.7 -0.4 | |
|
47 | h 6 5 36 32 28 23 18 13 8 4 0 -16 -12 -12 -11 -8 -6 -4 -2 -1 1 1 0.7 0.0 -0.2 | |
|
48 | g 6 6 -90 -92 -95 -98 -101 -103 -104 -106 -107 -104 -105 -107 -113 -111 -112 -111 -108 -102 -96 -93 -90.4 -86.4 1.3 | |
|
49 | h 6 6 -69 -67 -65 -62 -57 -52 -46 -40 -33 -39 -30 -24 -17 -7 1 11 17 21 24 36 43.8 50.3 0.9 | |
|
50 | g 7 0 70 70 71 72 73 73 74 74 74 70 65 65 67 75 72 71 72 74 77 77 79.0 79.8 -0.4 | |
|
51 | g 7 1 -55 -54 -54 -54 -54 -54 -54 -53 -53 -40 -55 -56 -56 -57 -57 -56 -59 -62 -64 -72 -74.0 -74.4 0.0 | |
|
52 | h 7 1 -45 -46 -47 -48 -49 -50 -51 -52 -52 -45 -35 -50 -55 -61 -70 -77 -82 -83 -80 -69 -64.6 -61.4 0.8 | |
|
53 | g 7 2 0 0 1 2 2 3 4 4 4 0 2 2 5 4 1 1 2 3 2 1 0.0 -1.4 -0.2 | |
|
54 | h 7 2 -13 -14 -14 -14 -14 -14 -15 -17 -18 -18 -17 -24 -28 -27 -27 -26 -27 -27 -26 -25 -24.2 -22.5 0.4 | |
|
55 | g 7 3 34 33 32 31 29 27 25 23 20 0 1 10 15 13 14 16 21 24 26 28 33.3 38.6 1.1 | |
|
56 | h 7 3 -10 -11 -12 -12 -13 -14 -14 -14 -14 2 0 -4 -6 -2 -4 -5 -5 -2 0 4 6.2 6.9 0.1 | |
|
57 | g 7 4 -41 -41 -40 -38 -37 -35 -34 -33 -31 -29 -40 -32 -32 -26 -22 -14 -12 -6 -1 5 9.1 12.3 0.6 | |
|
58 | h 7 4 -1 0 1 2 4 5 6 7 7 6 10 8 7 6 8 10 16 20 21 24 24.0 25.4 0.2 | |
|
59 | g 7 5 -21 -20 -19 -18 -16 -14 -12 -11 -9 -10 -7 -11 -7 -6 -2 0 1 4 5 4 6.9 9.4 0.4 | |
|
60 | h 7 5 28 28 28 28 28 29 29 29 29 28 36 28 23 26 23 22 18 17 17 17 14.8 10.9 -0.9 | |
|
61 | g 7 6 18 18 18 19 19 19 18 18 17 15 5 9 17 13 13 12 11 10 9 8 7.3 5.5 -0.5 | |
|
62 | h 7 6 -12 -12 -13 -15 -16 -17 -18 -19 -20 -17 -18 -20 -18 -23 -23 -23 -23 -23 -23 -24 -25.4 -26.4 -0.3 | |
|
63 | g 7 7 6 6 6 6 6 6 6 6 5 29 19 18 8 1 -2 -5 -2 0 0 -2 -1.2 2.0 0.9 | |
|
64 | h 7 7 -22 -22 -22 -22 -22 -21 -20 -19 -19 -22 -16 -18 -17 -12 -11 -12 -10 -7 -4 -6 -5.8 -4.8 0.3 | |
|
65 | g 8 0 11 11 11 11 11 11 11 11 11 13 22 11 15 13 14 14 18 21 23 25 24.4 24.8 -0.2 | |
|
66 | g 8 1 8 8 8 8 7 7 7 7 7 7 15 9 6 5 6 6 6 6 5 6 6.6 7.7 0.2 | |
|
67 | h 8 1 8 8 8 8 8 8 8 8 8 12 5 10 11 7 7 6 7 8 10 11 11.9 11.2 -0.2 | |
|
68 | g 8 2 -4 -4 -4 -4 -3 -3 -3 -3 -3 -8 -4 -6 -4 -4 -2 -1 0 0 -1 -6 -9.2 -11.4 -0.2 | |
|
69 | h 8 2 -14 -15 -15 -15 -15 -15 -15 -15 -14 -21 -22 -15 -14 -12 -15 -16 -18 -19 -19 -21 -21.5 -21.0 0.2 | |
|
70 | g 8 3 -9 -9 -9 -9 -9 -9 -9 -9 -10 -5 -1 -14 -11 -14 -13 -12 -11 -11 -10 -9 -7.9 -6.8 0.2 | |
|
71 | h 8 3 7 7 6 6 6 6 5 5 5 -12 0 5 7 9 6 4 4 5 6 8 8.5 9.7 0.2 | |
|
72 | g 8 4 1 1 1 2 2 2 2 1 1 9 11 6 2 0 -3 -8 -7 -9 -12 -14 -16.6 -18.0 -0.2 | |
|
73 | h 8 4 -13 -13 -13 -13 -14 -14 -14 -15 -15 -7 -21 -23 -18 -16 -17 -19 -22 -23 -22 -23 -21.5 -19.8 0.4 | |
|
74 | g 8 5 2 2 2 3 4 4 5 6 6 7 15 10 10 8 5 4 4 4 3 9 9.1 10.0 0.2 | |
|
75 | h 8 5 5 5 5 5 5 5 5 5 5 2 -8 3 4 4 6 6 9 11 12 15 15.5 16.1 0.2 | |
|
76 | g 8 6 -9 -8 -8 -8 -7 -7 -6 -6 -5 -10 -13 -7 -5 -1 0 0 3 4 4 6 7.0 9.4 0.5 | |
|
77 | h 8 6 16 16 16 16 17 17 18 18 19 18 17 23 23 24 21 18 16 14 12 11 8.9 7.7 -0.3 | |
|
78 | g 8 7 5 5 5 6 6 7 8 8 9 7 5 6 10 11 11 10 6 4 2 -5 -7.9 -11.4 -0.7 | |
|
79 | h 8 7 -5 -5 -5 -5 -5 -5 -5 -5 -5 3 -4 -4 1 -3 -6 -10 -13 -15 -16 -16 -14.9 -12.8 0.5 | |
|
80 | g 8 8 8 8 8 8 8 8 8 7 7 2 -1 9 8 4 3 1 -1 -4 -6 -7 -7.0 -5.0 0.5 | |
|
81 | h 8 8 -18 -18 -18 -18 -19 -19 -19 -19 -19 -11 -17 -13 -20 -17 -16 -17 -15 -11 -10 -4 -2.1 -0.1 0.4 | |
|
82 | g 9 0 8 8 8 8 8 8 8 8 8 5 3 4 4 8 8 7 5 5 4 4 5.0 5.6 | |
|
83 | g 9 1 10 10 10 10 10 10 10 10 10 -21 -7 9 6 10 10 10 10 10 9 9 9.4 9.8 | |
|
84 | h 9 1 -20 -20 -20 -20 -20 -20 -20 -20 -21 -27 -24 -11 -18 -22 -21 -21 -21 -21 -20 -20 -19.7 -20.1 | |
|
85 | g 9 2 1 1 1 1 1 1 1 1 1 1 -1 -4 0 2 2 2 1 1 1 3 3.0 3.6 | |
|
86 | h 9 2 14 14 14 14 14 14 14 15 15 17 19 12 12 15 16 16 16 15 15 15 13.4 12.9 | |
|
87 | g 9 3 -11 -11 -11 -11 -11 -11 -12 -12 -12 -11 -25 -5 -9 -13 -12 -12 -12 -12 -12 -10 -8.4 -7.0 | |
|
88 | h 9 3 5 5 5 5 5 5 5 5 5 29 12 7 2 7 6 7 9 9 11 12 12.5 12.7 | |
|
89 | g 9 4 12 12 12 12 12 12 12 11 11 3 10 2 1 10 10 10 9 9 9 8 6.3 5.0 | |
|
90 | h 9 4 -3 -3 -3 -3 -3 -3 -3 -3 -3 -9 2 6 0 -4 -4 -4 -5 -6 -7 -6 -6.2 -6.7 | |
|
91 | g 9 5 1 1 1 1 1 1 1 1 1 16 5 4 4 -1 -1 -1 -3 -3 -4 -8 -8.9 -10.8 | |
|
92 | h 9 5 -2 -2 -2 -2 -2 -2 -2 -3 -3 4 2 -2 -3 -5 -5 -5 -6 -6 -7 -8 -8.4 -8.1 | |
|
93 | g 9 6 -2 -2 -2 -2 -2 -2 -2 -2 -2 -3 -5 1 -1 -1 0 -1 -1 -1 -2 -1 -1.5 -1.3 | |
|
94 | h 9 6 8 8 8 8 9 9 9 9 9 9 8 10 9 10 10 10 9 9 9 8 8.4 8.1 | |
|
95 | g 9 7 2 2 2 2 2 2 3 3 3 -4 -2 2 -2 5 3 4 7 7 7 10 9.3 8.7 | |
|
96 | h 9 7 10 10 10 10 10 10 10 11 11 6 8 7 8 10 11 11 10 9 8 5 3.8 2.9 | |
|
97 | g 9 8 -1 0 0 0 0 0 0 0 1 -3 3 2 3 1 1 1 2 1 1 -2 -4.3 -6.7 | |
|
98 | h 9 8 -2 -2 -2 -2 -2 -2 -2 -2 -2 1 -11 -6 0 -4 -2 -3 -6 -7 -7 -8 -8.2 -7.9 | |
|
99 | g 9 9 -1 -1 -1 -1 -1 -1 -2 -2 -2 -4 8 5 -1 -2 -1 -2 -5 -5 -6 -8 -8.2 -9.2 | |
|
100 | h 9 9 2 2 2 2 2 2 2 2 2 8 -7 5 5 1 1 1 2 2 2 3 4.8 5.9 | |
|
101 | g 10 0 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -8 -3 1 -2 -3 -3 -4 -4 -3 -3 -2.6 -2.2 | |
|
102 | g 10 1 -4 -4 -4 -4 -4 -4 -4 -4 -4 11 4 -5 -3 -3 -3 -3 -4 -4 -4 -6 -6.0 -6.3 | |
|
103 | h 10 1 2 2 2 2 2 2 2 2 2 5 13 -4 4 2 1 1 1 1 2 1 1.7 2.4 | |
|
104 | g 10 2 2 2 2 2 2 2 2 2 2 1 -1 -1 4 2 2 2 2 3 2 2 1.7 1.6 | |
|
105 | h 10 2 1 1 1 1 1 1 1 1 1 1 -2 0 1 1 1 1 0 0 1 0 0.0 0.2 | |
|
106 | g 10 3 -5 -5 -5 -5 -5 -5 -5 -5 -5 2 13 2 0 -5 -5 -5 -5 -5 -5 -4 -3.1 -2.5 | |
|
107 | h 10 3 2 2 2 2 2 2 2 2 2 -20 -10 -8 0 2 3 3 3 3 3 4 4.0 4.4 | |
|
108 | g 10 4 -2 -2 -2 -2 -2 -2 -2 -2 -2 -5 -4 -3 -1 -2 -1 -2 -2 -2 -2 -1 -0.5 -0.1 | |
|
109 | h 10 4 6 6 6 6 6 6 6 6 6 -1 2 -2 2 6 4 4 6 6 6 5 4.9 4.7 | |
|
110 | g 10 5 6 6 6 6 6 6 6 6 6 -1 4 7 4 4 6 5 5 5 4 4 3.7 3.0 | |
|
111 | h 10 5 -4 -4 -4 -4 -4 -4 -4 -4 -4 -6 -3 -4 -5 -4 -4 -4 -4 -4 -4 -5 -5.9 -6.5 | |
|
112 | g 10 6 4 4 4 4 4 4 4 4 4 8 12 4 6 4 4 4 3 3 3 2 1.0 0.3 | |
|
113 | h 10 6 0 0 0 0 0 0 0 0 0 6 6 1 1 0 0 -1 0 0 0 -1 -1.2 -1.0 | |
|
114 | g 10 7 0 0 0 0 0 0 0 0 0 -1 3 -2 1 0 1 1 1 1 1 2 2.0 2.1 | |
|
115 | h 10 7 -2 -2 -2 -2 -2 -2 -2 -1 -1 -4 -3 -3 -1 -2 -1 -1 -1 -1 -2 -2 -2.9 -3.4 | |
|
116 | g 10 8 2 2 2 1 1 1 1 2 2 -3 2 6 -1 2 0 0 2 2 3 5 4.2 3.9 | |
|
117 | h 10 8 4 4 4 4 4 4 4 4 4 -2 6 7 6 3 3 3 4 4 3 1 0.2 -0.9 | |
|
118 | g 10 9 2 2 2 2 3 3 3 3 3 5 10 -2 2 2 3 3 3 3 3 1 0.3 -0.1 | |
|
119 | h 10 9 0 0 0 0 0 0 0 0 0 0 11 -1 0 0 1 1 0 0 -1 -2 -2.2 -2.3 | |
|
120 | g 10 10 0 0 0 0 0 0 0 0 0 -2 3 0 0 0 -1 -1 0 0 0 0 -1.1 -2.2 | |
|
121 | h 10 10 -6 -6 -6 -6 -6 -6 -6 -6 -6 -2 8 -3 -7 -6 -4 -5 -6 -6 -6 -7 -7.4 -8.0 | |
|
122 | g 11 0 2.7 2.9 | |
|
123 | g 11 1 -1.7 -1.6 | |
|
124 | h 11 1 0.1 0.3 | |
|
125 | g 11 2 -1.9 -1.7 | |
|
126 | h 11 2 1.3 1.4 | |
|
127 | g 11 3 1.5 1.5 | |
|
128 | h 11 3 -0.9 -0.7 | |
|
129 | g 11 4 -0.1 -0.2 | |
|
130 | h 11 4 -2.6 -2.4 | |
|
131 | g 11 5 0.1 0.2 | |
|
132 | h 11 5 0.9 0.9 | |
|
133 | g 11 6 -0.7 -0.7 | |
|
134 | h 11 6 -0.7 -0.6 | |
|
135 | g 11 7 0.7 0.5 | |
|
136 | h 11 7 -2.8 -2.7 | |
|
137 | g 11 8 1.7 1.8 | |
|
138 | h 11 8 -0.9 -1.0 | |
|
139 | g 11 9 0.1 0.1 | |
|
140 | h 11 9 -1.2 -1.5 | |
|
141 | g 11 10 1.2 1.0 | |
|
142 | h 11 10 -1.9 -2.0 | |
|
143 | g 11 11 4.0 4.1 | |
|
144 | h 11 11 -0.9 -1.4 | |
|
145 | g 12 0 -2.2 -2.2 | |
|
146 | g 12 1 -0.3 -0.3 | |
|
147 | h 12 1 -0.4 -0.5 | |
|
148 | g 12 2 0.2 0.3 | |
|
149 | h 12 2 0.3 0.3 | |
|
150 | g 12 3 0.9 0.9 | |
|
151 | h 12 3 2.5 2.3 | |
|
152 | g 12 4 -0.2 -0.4 | |
|
153 | h 12 4 -2.6 -2.7 | |
|
154 | g 12 5 0.9 1.0 | |
|
155 | h 12 5 0.7 0.6 | |
|
156 | g 12 6 -0.5 -0.4 | |
|
157 | h 12 6 0.3 0.4 | |
|
158 | g 12 7 0.3 0.5 | |
|
159 | h 12 7 0.0 0.0 | |
|
160 | g 12 8 -0.3 -0.3 | |
|
161 | h 12 8 0.0 0.0 | |
|
162 | g 12 9 -0.4 -0.4 | |
|
163 | h 12 9 0.3 0.3 | |
|
164 | g 12 10 -0.1 0.0 | |
|
165 | h 12 10 -0.9 -0.8 | |
|
166 | g 12 11 -0.2 -0.4 | |
|
167 | h 12 11 -0.4 -0.4 | |
|
168 | g 12 12 -0.4 0.0 | |
|
169 | h 12 12 0.8 1.0 | |
|
170 | g 13 0 -0.2 -0.2 | |
|
171 | g 13 1 -0.9 -0.9 | |
|
172 | h 13 1 -0.9 -0.7 | |
|
173 | g 13 2 0.3 0.3 | |
|
174 | h 13 2 0.2 0.3 | |
|
175 | g 13 3 0.1 0.3 | |
|
176 | h 13 3 1.8 1.7 | |
|
177 | g 13 4 -0.4 -0.4 | |
|
178 | h 13 4 -0.4 -0.5 | |
|
179 | g 13 5 1.3 1.2 | |
|
180 | h 13 5 -1.0 -1.0 | |
|
181 | g 13 6 -0.4 -0.4 | |
|
182 | h 13 6 -0.1 0.0 | |
|
183 | g 13 7 0.7 0.7 | |
|
184 | h 13 7 0.7 0.7 | |
|
185 | g 13 8 -0.4 -0.3 | |
|
186 | h 13 8 0.3 0.2 | |
|
187 | g 13 9 0.3 0.4 | |
|
188 | h 13 9 0.6 0.6 | |
|
189 | g 13 10 -0.1 -0.1 | |
|
190 | h 13 10 0.3 0.4 | |
|
191 | g 13 11 0.4 0.4 | |
|
192 | h 13 11 -0.2 -0.2 | |
|
193 | g 13 12 0.0 -0.1 | |
|
194 | h 13 12 -0.5 -0.5 | |
|
195 | g 13 13 0.1 -0.3 | |
|
196 | h 13 13 -0.9 -1.0 |
@@ -0,0 +1,196 | |||
|
1 | g/h n m 1900.0 1905.0 1910.0 1915.0 1920.0 1925.0 1930.0 1935.0 1940.0 1945.0 1950.0 1955.0 1960.0 1965.0 1970.0 1975.0 1980.0 1985.0 1990.0 1995.0 2000.0 2005.0 2010.0 SV | |
|
2 | g 1 0 -31543 -31464 -31354 -31212 -31060 -30926 -30805 -30715 -30654 -30594 -30554 -30500 -30421 -30334 -30220 -30100 -29992 -29873 -29775 -29692 -29619.4 -29554.63 -29496.5 11.4 | |
|
3 | g 1 1 -2298 -2298 -2297 -2306 -2317 -2318 -2316 -2306 -2292 -2285 -2250 -2215 -2169 -2119 -2068 -2013 -1956 -1905 -1848 -1784 -1728.2 -1669.05 -1585.9 16.7 | |
|
4 | h 1 1 5922 5909 5898 5875 5845 5817 5808 5812 5821 5810 5815 5820 5791 5776 5737 5675 5604 5500 5406 5306 5186.1 5077.99 4945.1 -28.8 | |
|
5 | g 2 0 -677 -728 -769 -802 -839 -893 -951 -1018 -1106 -1244 -1341 -1440 -1555 -1662 -1781 -1902 -1997 -2072 -2131 -2200 -2267.7 -2337.24 -2396.6 -11.3 | |
|
6 | g 2 1 2905 2928 2948 2956 2959 2969 2980 2984 2981 2990 2998 3003 3002 2997 3000 3010 3027 3044 3059 3070 3068.4 3047.69 3026.0 -3.9 | |
|
7 | h 2 1 -1061 -1086 -1128 -1191 -1259 -1334 -1424 -1520 -1614 -1702 -1810 -1898 -1967 -2016 -2047 -2067 -2129 -2197 -2279 -2366 -2481.6 -2594.50 -2707.7 -23.0 | |
|
8 | g 2 2 924 1041 1176 1309 1407 1471 1517 1550 1566 1578 1576 1581 1590 1594 1611 1632 1663 1687 1686 1681 1670.9 1657.76 1668.6 2.7 | |
|
9 | h 2 2 1121 1065 1000 917 823 728 644 586 528 477 381 291 206 114 25 -68 -200 -306 -373 -413 -458.0 -515.43 -575.4 -12.9 | |
|
10 | g 3 0 1022 1037 1058 1084 1111 1140 1172 1206 1240 1282 1297 1302 1302 1297 1287 1276 1281 1296 1314 1335 1339.6 1336.30 1339.7 1.3 | |
|
11 | g 3 1 -1469 -1494 -1524 -1559 -1600 -1645 -1692 -1740 -1790 -1834 -1889 -1944 -1992 -2038 -2091 -2144 -2180 -2208 -2239 -2267 -2288.0 -2305.83 -2326.3 -3.9 | |
|
12 | h 3 1 -330 -357 -389 -421 -445 -462 -480 -494 -499 -499 -476 -462 -414 -404 -366 -333 -336 -310 -284 -262 -227.6 -198.86 -160.5 8.6 | |
|
13 | g 3 2 1256 1239 1223 1212 1205 1202 1205 1215 1232 1255 1274 1288 1289 1292 1278 1260 1251 1247 1248 1249 1252.1 1246.39 1231.7 -2.9 | |
|
14 | h 3 2 3 34 62 84 103 119 133 146 163 186 206 216 224 240 251 262 271 284 293 302 293.4 269.72 251.7 -2.9 | |
|
15 | g 3 3 572 635 705 778 839 881 907 918 916 913 896 882 878 856 838 830 833 829 802 759 714.5 672.51 634.2 -8.1 | |
|
16 | h 3 3 523 480 425 360 293 229 166 101 43 -11 -46 -83 -130 -165 -196 -223 -252 -297 -352 -427 -491.1 -524.72 -536.8 -2.1 | |
|
17 | g 4 0 876 880 884 887 889 891 896 903 914 944 954 958 957 957 952 946 938 936 939 940 932.3 920.55 912.6 -1.4 | |
|
18 | g 4 1 628 643 660 678 695 711 727 744 762 776 792 796 800 804 800 791 782 780 780 780 786.8 797.96 809.0 2.0 | |
|
19 | h 4 1 195 203 211 218 220 216 205 188 169 144 136 133 135 148 167 191 212 232 247 262 272.6 282.07 286.4 0.4 | |
|
20 | g 4 2 660 653 644 631 616 601 584 565 550 544 528 510 504 479 461 438 398 361 325 290 250.0 210.65 166.6 -8.9 | |
|
21 | h 4 2 -69 -77 -90 -109 -134 -163 -195 -226 -252 -276 -278 -274 -278 -269 -266 -265 -257 -249 -240 -236 -231.9 -225.23 -211.2 3.2 | |
|
22 | g 4 3 -361 -380 -400 -416 -424 -426 -422 -415 -405 -421 -408 -397 -394 -390 -395 -405 -419 -424 -423 -418 -403.0 -379.86 -357.1 4.4 | |
|
23 | h 4 3 -210 -201 -189 -173 -153 -130 -109 -90 -72 -55 -37 -23 3 13 26 39 53 69 84 97 119.8 145.15 164.4 3.6 | |
|
24 | g 4 4 134 146 160 178 199 217 234 249 265 304 303 290 269 252 234 216 199 170 141 122 111.3 100.00 89.7 -2.3 | |
|
25 | h 4 4 -75 -65 -55 -51 -57 -70 -90 -114 -141 -178 -210 -230 -255 -269 -279 -288 -297 -297 -299 -306 -303.8 -305.36 -309.2 -0.8 | |
|
26 | g 5 0 -184 -192 -201 -211 -221 -230 -237 -241 -241 -253 -240 -229 -222 -219 -216 -218 -218 -214 -214 -214 -218.8 -227.00 -231.1 -0.5 | |
|
27 | g 5 1 328 328 327 327 326 326 327 329 334 346 349 360 362 358 359 356 357 355 353 352 351.4 354.41 357.2 0.5 | |
|
28 | h 5 1 -210 -193 -172 -148 -122 -96 -72 -51 -33 -12 3 15 16 19 26 31 46 47 46 46 43.8 42.72 44.7 0.5 | |
|
29 | g 5 2 264 259 253 245 236 226 218 211 208 194 211 230 242 254 262 264 261 253 245 235 222.3 208.95 200.3 -1.5 | |
|
30 | h 5 2 53 56 57 58 58 58 60 64 71 95 103 110 125 128 139 148 150 150 154 165 171.9 180.25 188.9 1.5 | |
|
31 | g 5 3 5 -1 -9 -16 -23 -28 -32 -33 -33 -20 -20 -23 -26 -31 -42 -59 -74 -93 -109 -118 -130.4 -136.54 -141.2 -0.7 | |
|
32 | h 5 3 -33 -32 -33 -34 -38 -44 -53 -64 -75 -67 -87 -98 -117 -126 -139 -152 -151 -154 -153 -143 -133.1 -123.45 -118.1 0.9 | |
|
33 | g 5 4 -86 -93 -102 -111 -119 -125 -131 -136 -141 -142 -147 -152 -156 -157 -160 -159 -162 -164 -165 -166 -168.6 -168.05 -163.1 1.3 | |
|
34 | h 5 4 -124 -125 -126 -126 -125 -122 -118 -115 -113 -119 -122 -121 -114 -97 -91 -83 -78 -75 -69 -55 -39.3 -19.57 0.1 3.7 | |
|
35 | g 5 5 -16 -26 -38 -51 -62 -69 -74 -76 -76 -82 -76 -69 -63 -62 -56 -49 -48 -46 -36 -17 -12.9 -13.55 -7.7 1.4 | |
|
36 | h 5 5 3 11 21 32 43 51 58 64 69 82 80 78 81 81 83 88 92 95 97 107 106.3 103.85 100.9 -0.6 | |
|
37 | g 6 0 63 62 62 61 61 61 60 59 57 59 54 47 46 45 43 45 48 53 61 68 72.3 73.60 72.8 -0.3 | |
|
38 | g 6 1 61 60 58 57 55 54 53 53 54 57 57 57 58 61 64 66 66 65 65 67 68.2 69.56 68.6 -0.3 | |
|
39 | h 6 1 -9 -7 -5 -2 0 3 4 4 4 6 -1 -9 -10 -11 -12 -13 -15 -16 -16 -17 -17.4 -20.33 -20.8 -0.1 | |
|
40 | g 6 2 -11 -11 -11 -10 -10 -9 -9 -8 -7 6 4 3 1 8 15 28 42 51 59 68 74.2 76.74 76.0 -0.3 | |
|
41 | h 6 2 83 86 89 93 96 99 102 104 105 100 99 96 99 100 100 99 93 88 82 72 63.7 54.75 44.2 -2.1 | |
|
42 | g 6 3 -217 -221 -224 -228 -233 -238 -242 -246 -249 -246 -247 -247 -237 -228 -212 -198 -192 -185 -178 -170 -160.9 -151.34 -141.4 1.9 | |
|
43 | h 6 3 2 4 5 8 11 14 19 25 33 16 33 48 60 68 72 75 71 69 69 67 65.1 63.63 61.5 -0.4 | |
|
44 | g 6 4 -58 -57 -54 -51 -46 -40 -32 -25 -18 -25 -16 -8 -1 4 2 1 4 4 3 -1 -5.9 -14.58 -22.9 -1.6 | |
|
45 | h 6 4 -35 -32 -29 -26 -22 -18 -16 -15 -15 -9 -12 -16 -20 -32 -37 -41 -43 -48 -52 -58 -61.2 -63.53 -66.3 -0.5 | |
|
46 | g 6 5 59 57 54 49 44 39 32 25 18 21 12 7 -2 1 3 6 14 16 18 19 16.9 14.58 13.1 -0.2 | |
|
47 | h 6 5 36 32 28 23 18 13 8 4 0 -16 -12 -12 -11 -8 -6 -4 -2 -1 1 1 0.7 0.24 3.1 0.8 | |
|
48 | g 6 6 -90 -92 -95 -98 -101 -103 -104 -106 -107 -104 -105 -107 -113 -111 -112 -111 -108 -102 -96 -93 -90.4 -86.36 -77.9 1.8 | |
|
49 | h 6 6 -69 -67 -65 -62 -57 -52 -46 -40 -33 -39 -30 -24 -17 -7 1 11 17 21 24 36 43.8 50.94 54.9 0.5 | |
|
50 | g 7 0 70 70 71 72 73 73 74 74 74 70 65 65 67 75 72 71 72 74 77 77 79.0 79.88 80.4 0.2 | |
|
51 | g 7 1 -55 -54 -54 -54 -54 -54 -54 -53 -53 -40 -55 -56 -56 -57 -57 -56 -59 -62 -64 -72 -74.0 -74.46 -75.0 -0.1 | |
|
52 | h 7 1 -45 -46 -47 -48 -49 -50 -51 -52 -52 -45 -35 -50 -55 -61 -70 -77 -82 -83 -80 -69 -64.6 -61.14 -57.8 0.6 | |
|
53 | g 7 2 0 0 1 2 2 3 4 4 4 0 2 2 5 4 1 1 2 3 2 1 0.0 -1.65 -4.7 -0.6 | |
|
54 | h 7 2 -13 -14 -14 -14 -14 -14 -15 -17 -18 -18 -17 -24 -28 -27 -27 -26 -27 -27 -26 -25 -24.2 -22.57 -21.2 0.3 | |
|
55 | g 7 3 34 33 32 31 29 27 25 23 20 0 1 10 15 13 14 16 21 24 26 28 33.3 38.73 45.3 1.4 | |
|
56 | h 7 3 -10 -11 -12 -12 -13 -14 -14 -14 -14 2 0 -4 -6 -2 -4 -5 -5 -2 0 4 6.2 6.82 6.6 -0.2 | |
|
57 | g 7 4 -41 -41 -40 -38 -37 -35 -34 -33 -31 -29 -40 -32 -32 -26 -22 -14 -12 -6 -1 5 9.1 12.30 14.0 0.3 | |
|
58 | h 7 4 -1 0 1 2 4 5 6 7 7 6 10 8 7 6 8 10 16 20 21 24 24.0 25.35 24.9 -0.1 | |
|
59 | g 7 5 -21 -20 -19 -18 -16 -14 -12 -11 -9 -10 -7 -11 -7 -6 -2 0 1 4 5 4 6.9 9.37 10.4 0.1 | |
|
60 | h 7 5 28 28 28 28 28 29 29 29 29 28 36 28 23 26 23 22 18 17 17 17 14.8 10.93 7.0 -0.8 | |
|
61 | g 7 6 18 18 18 19 19 19 18 18 17 15 5 9 17 13 13 12 11 10 9 8 7.3 5.42 1.6 -0.8 | |
|
62 | h 7 6 -12 -12 -13 -15 -16 -17 -18 -19 -20 -17 -18 -20 -18 -23 -23 -23 -23 -23 -23 -24 -25.4 -26.32 -27.7 -0.3 | |
|
63 | g 7 7 6 6 6 6 6 6 6 6 5 29 19 18 8 1 -2 -5 -2 0 0 -2 -1.2 1.94 4.9 0.4 | |
|
64 | h 7 7 -22 -22 -22 -22 -22 -21 -20 -19 -19 -22 -16 -18 -17 -12 -11 -12 -10 -7 -4 -6 -5.8 -4.64 -3.4 0.2 | |
|
65 | g 8 0 11 11 11 11 11 11 11 11 11 13 22 11 15 13 14 14 18 21 23 25 24.4 24.80 24.3 -0.1 | |
|
66 | g 8 1 8 8 8 8 7 7 7 7 7 7 15 9 6 5 6 6 6 6 5 6 6.6 7.62 8.2 0.1 | |
|
67 | h 8 1 8 8 8 8 8 8 8 8 8 12 5 10 11 7 7 6 7 8 10 11 11.9 11.20 10.9 0.0 | |
|
68 | g 8 2 -4 -4 -4 -4 -3 -3 -3 -3 -3 -8 -4 -6 -4 -4 -2 -1 0 0 -1 -6 -9.2 -11.73 -14.5 -0.5 | |
|
69 | h 8 2 -14 -15 -15 -15 -15 -15 -15 -15 -14 -21 -22 -15 -14 -12 -15 -16 -18 -19 -19 -21 -21.5 -20.88 -20.0 0.2 | |
|
70 | g 8 3 -9 -9 -9 -9 -9 -9 -9 -9 -10 -5 -1 -14 -11 -14 -13 -12 -11 -11 -10 -9 -7.9 -6.88 -5.7 0.3 | |
|
71 | h 8 3 7 7 6 6 6 6 5 5 5 -12 0 5 7 9 6 4 4 5 6 8 8.5 9.83 11.9 0.5 | |
|
72 | g 8 4 1 1 1 2 2 2 2 1 1 9 11 6 2 0 -3 -8 -7 -9 -12 -14 -16.6 -18.11 -19.3 -0.3 | |
|
73 | h 8 4 -13 -13 -13 -13 -14 -14 -14 -15 -15 -7 -21 -23 -18 -16 -17 -19 -22 -23 -22 -23 -21.5 -19.71 -17.4 0.4 | |
|
74 | g 8 5 2 2 2 3 4 4 5 6 6 7 15 10 10 8 5 4 4 4 3 9 9.1 10.17 11.6 0.3 | |
|
75 | h 8 5 5 5 5 5 5 5 5 5 5 2 -8 3 4 4 6 6 9 11 12 15 15.5 16.22 16.7 0.1 | |
|
76 | g 8 6 -9 -8 -8 -8 -7 -7 -6 -6 -5 -10 -13 -7 -5 -1 0 0 3 4 4 6 7.0 9.36 10.9 0.2 | |
|
77 | h 8 6 16 16 16 16 17 17 18 18 19 18 17 23 23 24 21 18 16 14 12 11 8.9 7.61 7.1 -0.1 | |
|
78 | g 8 7 5 5 5 6 6 7 8 8 9 7 5 6 10 11 11 10 6 4 2 -5 -7.9 -11.25 -14.1 -0.5 | |
|
79 | h 8 7 -5 -5 -5 -5 -5 -5 -5 -5 -5 3 -4 -4 1 -3 -6 -10 -13 -15 -16 -16 -14.9 -12.76 -10.8 0.4 | |
|
80 | g 8 8 8 8 8 8 8 8 8 7 7 2 -1 9 8 4 3 1 -1 -4 -6 -7 -7.0 -4.87 -3.7 0.2 | |
|
81 | h 8 8 -18 -18 -18 -18 -19 -19 -19 -19 -19 -11 -17 -13 -20 -17 -16 -17 -15 -11 -10 -4 -2.1 -0.06 1.7 0.4 | |
|
82 | g 9 0 8 8 8 8 8 8 8 8 8 5 3 4 4 8 8 7 5 5 4 4 5.0 5.58 5.4 0.0 | |
|
83 | g 9 1 10 10 10 10 10 10 10 10 10 -21 -7 9 6 10 10 10 10 10 9 9 9.4 9.76 9.4 0.0 | |
|
84 | h 9 1 -20 -20 -20 -20 -20 -20 -20 -20 -21 -27 -24 -11 -18 -22 -21 -21 -21 -21 -20 -20 -19.7 -20.11 -20.5 0.0 | |
|
85 | g 9 2 1 1 1 1 1 1 1 1 1 1 -1 -4 0 2 2 2 1 1 1 3 3.0 3.58 3.4 0.0 | |
|
86 | h 9 2 14 14 14 14 14 14 14 15 15 17 19 12 12 15 16 16 16 15 15 15 13.4 12.69 11.6 0.0 | |
|
87 | g 9 3 -11 -11 -11 -11 -11 -11 -12 -12 -12 -11 -25 -5 -9 -13 -12 -12 -12 -12 -12 -10 -8.4 -6.94 -5.3 0.0 | |
|
88 | h 9 3 5 5 5 5 5 5 5 5 5 29 12 7 2 7 6 7 9 9 11 12 12.5 12.67 12.8 0.0 | |
|
89 | g 9 4 12 12 12 12 12 12 12 11 11 3 10 2 1 10 10 10 9 9 9 8 6.3 5.01 3.1 0.0 | |
|
90 | h 9 4 -3 -3 -3 -3 -3 -3 -3 -3 -3 -9 2 6 0 -4 -4 -4 -5 -6 -7 -6 -6.2 -6.72 -7.2 0.0 | |
|
91 | g 9 5 1 1 1 1 1 1 1 1 1 16 5 4 4 -1 -1 -1 -3 -3 -4 -8 -8.9 -10.76 -12.4 0.0 | |
|
92 | h 9 5 -2 -2 -2 -2 -2 -2 -2 -3 -3 4 2 -2 -3 -5 -5 -5 -6 -6 -7 -8 -8.4 -8.16 -7.4 0.0 | |
|
93 | g 9 6 -2 -2 -2 -2 -2 -2 -2 -2 -2 -3 -5 1 -1 -1 0 -1 -1 -1 -2 -1 -1.5 -1.25 -0.8 0.0 | |
|
94 | h 9 6 8 8 8 8 9 9 9 9 9 9 8 10 9 10 10 10 9 9 9 8 8.4 8.10 8.0 0.0 | |
|
95 | g 9 7 2 2 2 2 2 2 3 3 3 -4 -2 2 -2 5 3 4 7 7 7 10 9.3 8.76 8.4 0.0 | |
|
96 | h 9 7 10 10 10 10 10 10 10 11 11 6 8 7 8 10 11 11 10 9 8 5 3.8 2.92 2.2 0.0 | |
|
97 | g 9 8 -1 0 0 0 0 0 0 0 1 -3 3 2 3 1 1 1 2 1 1 -2 -4.3 -6.66 -8.4 0.0 | |
|
98 | h 9 8 -2 -2 -2 -2 -2 -2 -2 -2 -2 1 -11 -6 0 -4 -2 -3 -6 -7 -7 -8 -8.2 -7.73 -6.1 0.0 | |
|
99 | g 9 9 -1 -1 -1 -1 -1 -1 -2 -2 -2 -4 8 5 -1 -2 -1 -2 -5 -5 -6 -8 -8.2 -9.22 -10.1 0.0 | |
|
100 | h 9 9 2 2 2 2 2 2 2 2 2 8 -7 5 5 1 1 1 2 2 2 3 4.8 6.01 7.0 0.0 | |
|
101 | g 10 0 -3 -3 -3 -3 -3 -3 -3 -3 -3 -3 -8 -3 1 -2 -3 -3 -4 -4 -3 -3 -2.6 -2.17 -2.0 0.0 | |
|
102 | g 10 1 -4 -4 -4 -4 -4 -4 -4 -4 -4 11 4 -5 -3 -3 -3 -3 -4 -4 -4 -6 -6.0 -6.12 -6.3 0.0 | |
|
103 | h 10 1 2 2 2 2 2 2 2 2 2 5 13 -4 4 2 1 1 1 1 2 1 1.7 2.19 2.8 0.0 | |
|
104 | g 10 2 2 2 2 2 2 2 2 2 2 1 -1 -1 4 2 2 2 2 3 2 2 1.7 1.42 0.9 0.0 | |
|
105 | h 10 2 1 1 1 1 1 1 1 1 1 1 -2 0 1 1 1 1 0 0 1 0 0.0 0.10 -0.1 0.0 | |
|
106 | g 10 3 -5 -5 -5 -5 -5 -5 -5 -5 -5 2 13 2 0 -5 -5 -5 -5 -5 -5 -4 -3.1 -2.35 -1.1 0.0 | |
|
107 | h 10 3 2 2 2 2 2 2 2 2 2 -20 -10 -8 0 2 3 3 3 3 3 4 4.0 4.46 4.7 0.0 | |
|
108 | g 10 4 -2 -2 -2 -2 -2 -2 -2 -2 -2 -5 -4 -3 -1 -2 -1 -2 -2 -2 -2 -1 -0.5 -0.15 -0.2 0.0 | |
|
109 | h 10 4 6 6 6 6 6 6 6 6 6 -1 2 -2 2 6 4 4 6 6 6 5 4.9 4.76 4.4 0.0 | |
|
110 | g 10 5 6 6 6 6 6 6 6 6 6 -1 4 7 4 4 6 5 5 5 4 4 3.7 3.06 2.5 0.0 | |
|
111 | h 10 5 -4 -4 -4 -4 -4 -4 -4 -4 -4 -6 -3 -4 -5 -4 -4 -4 -4 -4 -4 -5 -5.9 -6.58 -7.2 0.0 | |
|
112 | g 10 6 4 4 4 4 4 4 4 4 4 8 12 4 6 4 4 4 3 3 3 2 1.0 0.29 -0.3 0.0 | |
|
113 | h 10 6 0 0 0 0 0 0 0 0 0 6 6 1 1 0 0 -1 0 0 0 -1 -1.2 -1.01 -1.0 0.0 | |
|
114 | g 10 7 0 0 0 0 0 0 0 0 0 -1 3 -2 1 0 1 1 1 1 1 2 2.0 2.06 2.2 0.0 | |
|
115 | h 10 7 -2 -2 -2 -2 -2 -2 -2 -1 -1 -4 -3 -3 -1 -2 -1 -1 -1 -1 -2 -2 -2.9 -3.47 -4.0 0.0 | |
|
116 | g 10 8 2 2 2 1 1 1 1 2 2 -3 2 6 -1 2 0 0 2 2 3 5 4.2 3.77 3.1 0.0 | |
|
117 | h 10 8 4 4 4 4 4 4 4 4 4 -2 6 7 6 3 3 3 4 4 3 1 0.2 -0.86 -2.0 0.0 | |
|
118 | g 10 9 2 2 2 2 3 3 3 3 3 5 10 -2 2 2 3 3 3 3 3 1 0.3 -0.21 -1.0 0.0 | |
|
119 | h 10 9 0 0 0 0 0 0 0 0 0 0 11 -1 0 0 1 1 0 0 -1 -2 -2.2 -2.31 -2.0 0.0 | |
|
120 | g 10 10 0 0 0 0 0 0 0 0 0 -2 3 0 0 0 -1 -1 0 0 0 0 -1.1 -2.09 -2.8 0.0 | |
|
121 | h 10 10 -6 -6 -6 -6 -6 -6 -6 -6 -6 -2 8 -3 -7 -6 -4 -5 -6 -6 -6 -7 -7.4 -7.93 -8.3 0.0 | |
|
122 | g 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.7 2.95 3.0 0.0 | |
|
123 | g 11 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1.7 -1.60 -1.5 0.0 | |
|
124 | h 11 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1 0.26 0.1 0.0 | |
|
125 | g 11 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1.9 -1.88 -2.1 0.0 | |
|
126 | h 11 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.3 1.44 1.7 0.0 | |
|
127 | g 11 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.5 1.44 1.6 0.0 | |
|
128 | h 11 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.9 -0.77 -0.6 0.0 | |
|
129 | g 11 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.1 -0.31 -0.5 0.0 | |
|
130 | h 11 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2.6 -2.27 -1.8 0.0 | |
|
131 | g 11 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1 0.29 0.5 0.0 | |
|
132 | h 11 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.9 0.90 0.9 0.0 | |
|
133 | g 11 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.7 -0.79 -0.8 0.0 | |
|
134 | h 11 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.7 -0.58 -0.4 0.0 | |
|
135 | g 11 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.7 0.53 0.4 0.0 | |
|
136 | h 11 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2.8 -2.69 -2.5 0.0 | |
|
137 | g 11 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.7 1.80 1.8 0.0 | |
|
138 | h 11 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.9 -1.08 -1.3 0.0 | |
|
139 | g 11 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1 0.16 0.2 0.0 | |
|
140 | h 11 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1.2 -1.58 -2.1 0.0 | |
|
141 | g 11 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.2 0.96 0.8 0.0 | |
|
142 | h 11 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1.9 -1.90 -1.9 0.0 | |
|
143 | g 11 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4.0 3.99 3.8 0.0 | |
|
144 | h 11 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.9 -1.39 -1.8 0.0 | |
|
145 | g 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2.2 -2.15 -2.1 0.0 | |
|
146 | g 12 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.3 -0.29 -0.2 0.0 | |
|
147 | h 12 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.4 -0.55 -0.8 0.0 | |
|
148 | g 12 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.2 0.21 0.3 0.0 | |
|
149 | h 12 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.3 0.23 0.3 0.0 | |
|
150 | g 12 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.9 0.89 1.0 0.0 | |
|
151 | h 12 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2.5 2.38 2.2 0.0 | |
|
152 | g 12 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.2 -0.38 -0.7 0.0 | |
|
153 | h 12 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2.6 -2.63 -2.5 0.0 | |
|
154 | g 12 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.9 0.96 0.9 0.0 | |
|
155 | h 12 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.7 0.61 0.5 0.0 | |
|
156 | g 12 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.5 -0.30 -0.1 0.0 | |
|
157 | h 12 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.3 0.40 0.6 0.0 | |
|
158 | g 12 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.3 0.46 0.5 0.0 | |
|
159 | h 12 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0 0.01 0.0 0.0 | |
|
160 | g 12 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.3 -0.35 -0.4 0.0 | |
|
161 | h 12 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0 0.02 0.1 0.0 | |
|
162 | g 12 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.4 -0.36 -0.4 0.0 | |
|
163 | h 12 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.3 0.28 0.3 0.0 | |
|
164 | g 12 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.1 0.08 0.2 0.0 | |
|
165 | h 12 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.9 -0.87 -0.9 0.0 | |
|
166 | g 12 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.2 -0.49 -0.8 0.0 | |
|
167 | h 12 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.4 -0.34 -0.2 0.0 | |
|
168 | g 12 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.4 -0.08 0.0 0.0 | |
|
169 | h 12 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.8 0.88 0.8 0.0 | |
|
170 | g 13 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.2 -0.16 -0.2 0.0 | |
|
171 | g 13 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.9 -0.88 -0.9 0.0 | |
|
172 | h 13 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.9 -0.76 -0.8 0.0 | |
|
173 | g 13 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.3 0.30 0.3 0.0 | |
|
174 | h 13 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.2 0.33 0.3 0.0 | |
|
175 | g 13 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1 0.28 0.4 0.0 | |
|
176 | h 13 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.8 1.72 1.7 0.0 | |
|
177 | g 13 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.4 -0.43 -0.4 0.0 | |
|
178 | h 13 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.4 -0.54 -0.6 0.0 | |
|
179 | g 13 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.3 1.18 1.1 0.0 | |
|
180 | h 13 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1.0 -1.07 -1.2 0.0 | |
|
181 | g 13 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.4 -0.37 -0.3 0.0 | |
|
182 | h 13 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.1 -0.04 -0.1 0.0 | |
|
183 | g 13 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.7 0.75 0.8 0.0 | |
|
184 | h 13 7 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.7 0.63 0.5 0.0 | |
|
185 | g 13 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.4 -0.26 -0.2 0.0 | |
|
186 | h 13 8 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.3 0.21 0.1 0.0 | |
|
187 | g 13 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.3 0.35 0.4 0.0 | |
|
188 | h 13 9 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.6 0.53 0.5 0.0 | |
|
189 | g 13 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.1 -0.05 0.0 0.0 | |
|
190 | h 13 10 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.3 0.38 0.4 0.0 | |
|
191 | g 13 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.4 0.41 0.4 0.0 | |
|
192 | h 13 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.2 -0.22 -0.2 0.0 | |
|
193 | g 13 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0 -0.10 -0.3 0.0 | |
|
194 | h 13 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.5 -0.57 -0.5 0.0 | |
|
195 | g 13 13 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.1 -0.18 -0.3 0.0 | |
|
196 | h 13 13 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.9 -0.82 -0.8 0.0 |
|
1 | NO CONTENT: new file 100644 | |
The requested commit or file is too big and content was truncated. Show full diff |
|
1 | NO CONTENT: new file 100644, binary diff hidden |
@@ -148,9 +148,3 h2 { | |||
|
148 | 148 | .container, .container-lg, .container-md, .container-sm, .container-xl { |
|
149 | 149 | max-width: 1020px; |
|
150 | 150 | } |
|
151 | ||
|
152 | @media (min-width: 1380px) { | |
|
153 | .card-columns { | |
|
154 | column-count: 5; | |
|
155 | } | |
|
156 | } |
@@ -46,57 +46,58 | |||
|
46 | 46 | <div class="invalid-tooltip"> |
|
47 | 47 | Please enter a valid date. |
|
48 | 48 | </div> |
|
49 |
<select |
|
|
49 | <select id="overjro-experiment" class="form-control form-control-sm"> | |
|
50 | 50 | <option value="-1">Experiment:</option> |
|
51 | 51 | <option value="-1">------------------</option> |
|
52 |
<option value=" |
|
|
53 |
<option value=" |
|
|
54 |
<option value=" |
|
|
55 |
<option value=" |
|
|
56 |
<option value=" |
|
|
57 |
<option value=" |
|
|
58 |
<option value=" |
|
|
59 |
|
|
|
52 | <option value="50">Vertical Drifts</option> | |
|
53 | <option value="51">East West 1996 (W beam)</option> | |
|
54 | <option value="52">East West 1996 (E beam)</option> | |
|
55 | <option value="61">East West 2003</option> | |
|
56 | <option value="60">Differential Phase 2000</option> | |
|
57 | <option value="63">Differential Phase 2004 High Alt</option> | |
|
58 | <option value="64">Differential Phase 2005 - 2006</option> | |
|
59 | <option value="54">DEWD 2005</option> | |
|
60 | <option value="53">DVD 2006 - 2008</option> | |
|
60 | 61 | <option value="-1">------------------</option> |
|
61 |
<option value=" |
|
|
62 |
<option value=" |
|
|
63 |
<option value=" |
|
|
64 |
<option value=" |
|
|
62 | <option value="4">Oblique ISR On-Axis</option> | |
|
63 | <option value="5">Oblique ISR 4.5</option> | |
|
64 | <option value="6">Oblique ISR 6.0S</option> | |
|
65 | <option value="7">Oblique ISR 3.0N</option> | |
|
65 | 66 | <option value="-1">------------------</option> |
|
66 |
<option value=" |
|
|
67 |
<option value=" |
|
|
68 |
<option value=" |
|
|
69 |
<option value=" |
|
|
70 |
|
|
|
67 | <option value="16">JULIA CP2</option> | |
|
68 | <option value="17">JULIA CP3</option> | |
|
69 | <option value="18">JULIA V (2005-2006)</option> | |
|
70 | <option value="65">JULIA EW 2003</option> | |
|
71 | <option value="19">JULIA EW (2006-2007)</option> | |
|
71 | 72 | <option value="-1">------------------</option> |
|
72 | 73 | <option value="0">Modulo Rx</option> |
|
73 | 74 | <option value="1">1/16 Rx</option> |
|
74 | 75 | <option value="2">1/4 Rx</option> |
|
75 | 76 | <option value="3">All Rx</option> |
|
76 | 77 | <option value="-1">------------------</option> |
|
77 |
<option value=" |
|
|
78 |
<option value=" |
|
|
79 |
|
|
|
78 | <option value="21">EW Imaging 1996</option> | |
|
79 | <option value="22">EW Imaging 2003</option> | |
|
80 | <option value="23">EW Imaging 2006-2008</option> | |
|
80 | 81 | <option value="-1">------------------</option> |
|
81 |
<option value=" |
|
|
82 |
<option value=" |
|
|
83 |
<option value=" |
|
|
84 |
<option value=" |
|
|
82 | <option value="359">MST North (Fritts)</option> | |
|
83 | <option value="360">MST West (Fritts)</option> | |
|
84 | <option value="361">MST South (Fritts)</option> | |
|
85 | <option value="362">MST East (Fritts)</option> | |
|
85 | 86 | <option value="-1">------------------</option> |
|
86 |
<option value=" |
|
|
87 | <option value="67">Vertical (Yellow Cables)</option> | |
|
87 | 88 | </select> |
|
88 | 89 | <br> |
|
89 | 90 | <p class="card-text">Choose object: |
|
90 | 91 | <div class="form-check card-text"> |
|
91 |
<input class="form-check-input" type="checkbox" id="inlineCheckbox1" value=" |
|
|
92 | <input class="form-check-input" type="checkbox" id="inlineCheckbox1" value="bfield" name="celestial"> | |
|
92 | 93 | <label class="form-check-label" for="inlineCheckbox1">B Field</label><br> |
|
93 |
<input class="form-check-input" type="checkbox" id="inlineCheckbox1" value=" |
|
|
94 | <input class="form-check-input" type="checkbox" id="inlineCheckbox1" value="sun" name="celestial"> | |
|
94 | 95 | <label class="form-check-label" for="inlineCheckbox1">Sun</label><br> |
|
95 |
<input class="form-check-input" type="checkbox" id="inlineCheckbox1" value=" |
|
|
96 | <input class="form-check-input" type="checkbox" id="inlineCheckbox1" value="moon" name="celestial"> | |
|
96 | 97 | <label class="form-check-label" for="inlineCheckbox1">Moon</label><br> |
|
97 |
<input class="form-check-input" type="checkbox" id="inlineCheckbox1" value=" |
|
|
98 | <input class="form-check-input" type="checkbox" id="inlineCheckbox1" value="hydra" name="celestial"> | |
|
98 | 99 | <label class="form-check-label" for="inlineCheckbox1">Hydra</label><br> |
|
99 |
<input class="form-check-input" type="checkbox" id="inlineCheckbox1" value=" |
|
|
100 | <input class="form-check-input" type="checkbox" id="inlineCheckbox1" value="galaxy" name="celestial"> | |
|
100 | 101 | <label class="form-check-label" for="inlineCheckbox1">Galaxy Center</label> |
|
101 | 102 | </div> |
|
102 | 103 | <br> |
@@ -106,9 +107,9 | |||
|
106 | 107 | value="5.0" required> |
|
107 | 108 | </div> |
|
108 | 109 | <div class="form-group card-text"> |
|
109 | <label class="form-check-label" for="overjro-height">Height [km]:</label> | |
|
110 | <input type="text" class="form-control form-control-sm tools-date" id="overjro-height" placeholder="Enter Height" | |
|
111 | value="" required> | |
|
110 | <label class="form-check-label" for="overjro-height">Heights [km]:</label> | |
|
111 | <input type="text" class="form-control form-control-sm tools-date" id="overjro-height" placeholder="Enter Heights [km]" | |
|
112 | value="100" required> | |
|
112 | 113 | </div> |
|
113 | 114 | </p> |
|
114 | 115 | </div> |
@@ -142,6 +143,9 | |||
|
142 | 143 | <div class="modal-body text-center"> |
|
143 | 144 | <img class="img-fluid" src=""> |
|
144 | 145 | </div> |
|
146 | <div class="modal-body text-center"> | |
|
147 | <p></p> | |
|
148 | </div> | |
|
145 | 149 | </div> |
|
146 | 150 | </div> |
|
147 | 151 | </div> |
@@ -156,15 +160,39 | |||
|
156 | 160 | //get data attribute of the clicked element |
|
157 | 161 | var title = $(e.relatedTarget).data('title'); |
|
158 | 162 | var image = $(e.relatedTarget).data('image'); |
|
163 | $(e.currentTarget).find('p').text(''); | |
|
164 | $(e.currentTarget).find('img').attr('src', ''); | |
|
159 | 165 | |
|
160 | 166 | if (image.indexOf('skynoise') > 0) { |
|
161 | 167 | var dt = $('#skynoise-date').val(); |
|
162 | 168 | image += '?date=' + dt; |
|
169 | //populate values | |
|
170 | $(e.currentTarget).find('h5').text(title); | |
|
171 | $(e.currentTarget).find('img').attr('src', image); | |
|
163 | 172 | } |
|
164 | 173 | |
|
165 | //populate values | |
|
166 | $(e.currentTarget).find('h5').text(title); | |
|
167 | $(e.currentTarget).find('img').attr('src', image); | |
|
174 | if (image.indexOf('overjro') > 0) { | |
|
175 | $(e.currentTarget).find('h5').text(title); | |
|
176 | ||
|
177 | if ($('#overjro-experiment').val() == '-1'){ | |
|
178 | $(e.currentTarget).find('p').text('Missing Experiment'); | |
|
179 | } else { | |
|
180 | ||
|
181 | var dt = $('#overjro-date').val(); | |
|
182 | var favorite = []; | |
|
183 | $.each($("input[name='celestial']:checked"), function(){ | |
|
184 | favorite.push($(this).val()); | |
|
185 | }); | |
|
186 | ||
|
187 | image += '?date=' + dt; | |
|
188 | image += '&experiment=' + $('#overjro-experiment').val(); | |
|
189 | image += '&angle=' + $('#overjro-angle').val(); | |
|
190 | image += '&height=' + $('#overjro-height').val(); | |
|
191 | image += '&bodys=' + favorite.join(","); | |
|
192 | ||
|
193 | $(e.currentTarget).find('img').attr('src', image); | |
|
194 | } | |
|
195 | } | |
|
168 | 196 | }); |
|
169 | 197 | |
|
170 | 198 | $('#doy-date').change(function() { |
@@ -16,7 +16,7 import mongoengine | |||
|
16 | 16 | |
|
17 | 17 | from plotter.models import Experiment, ExpDetail, PlotMeta, PlotData, JROReport |
|
18 | 18 | |
|
19 | from utils.plots import skynoise_plot | |
|
19 | from utils.plots import skynoise_plot, overjro_plot | |
|
20 | 20 | |
|
21 | 21 | host = os.environ.get('HOST_MONGO', 'localhost') |
|
22 | 22 | mongoengine.connect('dbplots', host=host, port=27017) |
@@ -258,7 +258,12 def plot_overjro(request): | |||
|
258 | 258 | else: |
|
259 | 259 | date = datetime.strptime(date, '%Y-%m-%d') |
|
260 | 260 | |
|
261 | data = skynoise_plot(date.year, date.month, date.day) | |
|
261 | pattern = int(request.GET.get('experiment', '1')) | |
|
262 | angle = float(request.GET.get('angle', '5')) | |
|
263 | height = [float(h) for h in request.GET.get('height', '100').split(',')] | |
|
264 | bodys = (request.GET.get('bodys', '')).split(',') | |
|
265 | ||
|
266 | data = overjro_plot(pattern, date, angle, height, bodys) | |
|
262 | 267 | response = HttpResponse(data.getvalue(), content_type='image/png') |
|
263 | 268 | |
|
264 | 269 | return response |
This diff has been collapsed as it changes many lines, (1661 lines changed) Show them Hide them | |||
@@ -4,11 +4,1637 import time | |||
|
4 | 4 | import numpy |
|
5 | 5 | import scipy |
|
6 | 6 | import base64 |
|
7 | import datetime | |
|
7 | 8 | from io import BytesIO |
|
8 | 9 | import scipy.interpolate |
|
9 | 10 | from matplotlib.figure import Figure |
|
10 | 11 | |
|
11 | 12 | from utils import TimeTools |
|
13 | from utils.patterns import select_pattern | |
|
14 | from utils import Misc_Routines | |
|
15 | from utils import Astro_Coords | |
|
16 | from utils import gaussfit | |
|
17 | ||
|
18 | ||
|
19 | attenuation = numpy.array([[[-21.25,-15.25,-9.25,-3.25,3.25,9.25,15.25,21.25], | |
|
20 | [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], | |
|
21 | [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], | |
|
22 | [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], | |
|
23 | [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], | |
|
24 | [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], | |
|
25 | [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25], | |
|
26 | [-21.25,-15.25,-9.25,-3.25,03.25,09.25,15.25,21.25]], | |
|
27 | [[21.25,21.25,21.25,21.25,21.25,21.25,21.25,21.25], | |
|
28 | [15.25,15.25,15.25,15.25,15.25,15.25,15.25,15.25], | |
|
29 | [09.25,09.25,09.25,09.25,09.25,09.25,09.25,09.25], | |
|
30 | [03.25,03.25,03.25,03.25,03.25,03.25,03.25,03.25], | |
|
31 | [-03.25,-03.25,-03.25,-03.25,-03.25,-03.25,-03.25,-03.25], | |
|
32 | [-09.25,-09.25,-09.25,-09.25,-09.25,-09.25,-09.25,-09.25], | |
|
33 | [-15.25,-15.25,-15.25,-15.25,-15.25,-15.25,-15.25,-15.25], | |
|
34 | [-21.25,-21.25,-21.25,-21.25,-21.25,-21.25,-21.25,-21.25]]]) | |
|
35 | ||
|
36 | class BField(): | |
|
37 | def __init__(self,year=None,doy=None,site=1,heights=None,alpha_i=90): | |
|
38 | """ | |
|
39 | BField class creates an object to get the Magnetic field for a specific date and | |
|
40 | height(s). | |
|
41 | ||
|
42 | Parameters | |
|
43 | ---------- | |
|
44 | year = A scalar giving the desired year. If the value is None (default value) then | |
|
45 | the current year will be used. | |
|
46 | doy = A scalar giving the desired day of the year. If the value is None (default va- | |
|
47 | lue) then the current doy will be used. | |
|
48 | site = An integer to choose the geographic coordinates of the place where the magne- | |
|
49 | tic field will be computed. The default value is over Jicamarca (site=1) | |
|
50 | heights = An array giving the heights (km) where the magnetic field will be modeled By default the magnetic field will be computed at 100, 500 and 1000km. | |
|
51 | alpha_i = Angle to interpolate the magnetic field. | |
|
52 | ||
|
53 | Modification History | |
|
54 | -------------------- | |
|
55 | Converted to Object-oriented Programming by Freddy Galindo, ROJ, 07 October 2009. | |
|
56 | """ | |
|
57 | ||
|
58 | tmp = time.localtime() | |
|
59 | if year==None: year = tmp[0] | |
|
60 | if doy==None: doy = tmp[7] | |
|
61 | self.year = year | |
|
62 | self.doy = doy | |
|
63 | self.site = site | |
|
64 | if heights is None: | |
|
65 | heights = numpy.array([100,500,1000]) | |
|
66 | else: | |
|
67 | heights = numpy.array(heights) | |
|
68 | self.heights = heights | |
|
69 | self.alpha_i = alpha_i | |
|
70 | ||
|
71 | def getBField(self,maglimits=numpy.array([-7,-7,7,7])): | |
|
72 | """ | |
|
73 | getBField models the magnetic field for a different heights in a specific date. | |
|
74 | ||
|
75 | Parameters | |
|
76 | ---------- | |
|
77 | maglimits = An 4-elements array giving ..... The default value is [-7,-7,7,7]. | |
|
78 | ||
|
79 | Return | |
|
80 | ------ | |
|
81 | dcos = An 4-dimensional array giving the directional cosines of the magnetic field | |
|
82 | over the desired place. | |
|
83 | alpha = An 3-dimensional array giving the angle of the magnetic field over the desi- | |
|
84 | red place. | |
|
85 | ||
|
86 | Modification History | |
|
87 | -------------------- | |
|
88 | Converted to Python by Freddy R. Galindo, ROJ, 07 October 2009. | |
|
89 | """ | |
|
90 | ||
|
91 | x_ant = numpy.array([1,0,0]) | |
|
92 | y_ant = numpy.array([0,1,0]) | |
|
93 | z_ant = numpy.array([0,0,1]) | |
|
94 | ||
|
95 | if self.site==0: | |
|
96 | title_site = "Magnetic equator" | |
|
97 | coord_site = numpy.array([-76+52./60.,-11+57/60.,0.5]) | |
|
98 | elif self.site==1: | |
|
99 | title_site = 'Jicamarca' | |
|
100 | coord_site = [-76-52./60.,-11-57/60.,0.5] | |
|
101 | theta = (45+5.35)*numpy.pi/180. # (50.35 and 1.46 from Fleish Thesis) | |
|
102 | delta = -1.46*numpy.pi/180 | |
|
103 | ||
|
104 | x_ant1 = numpy.roll(self.rotvector(self.rotvector(x_ant,1,delta),3,theta),1) | |
|
105 | y_ant1 = numpy.roll(self.rotvector(self.rotvector(y_ant,1,delta),3,theta),1) | |
|
106 | z_ant1 = numpy.roll(self.rotvector(self.rotvector(z_ant,1,delta),3,theta),1) | |
|
107 | ||
|
108 | ang0 = -1*coord_site[0]*numpy.pi/180. | |
|
109 | ang1 = coord_site[1]*numpy.pi/180. | |
|
110 | x_ant = self.rotvector(self.rotvector(x_ant1,2,ang1),3,ang0) | |
|
111 | y_ant = self.rotvector(self.rotvector(y_ant1,2,ang1),3,ang0) | |
|
112 | z_ant = self.rotvector(self.rotvector(z_ant1,2,ang1),3,ang0) | |
|
113 | else: | |
|
114 | # print "No defined Site. Skip..." | |
|
115 | return None | |
|
116 | ||
|
117 | nhei = self.heights.size | |
|
118 | pt_intercep = numpy.zeros((nhei,2)) | |
|
119 | nfields = 1 | |
|
120 | ||
|
121 | grid_res = 0.5 | |
|
122 | nlon = int(numpy.int(maglimits[2] - maglimits[0])/grid_res + 1) | |
|
123 | nlat = int(numpy.int(maglimits[3] - maglimits[1])/grid_res + 1) | |
|
124 | ||
|
125 | location = numpy.zeros((nlon,nlat,2)) | |
|
126 | mlon = numpy.atleast_2d(numpy.arange(nlon)*grid_res + maglimits[0]) | |
|
127 | mrep = numpy.atleast_2d(numpy.zeros(nlat) + 1) | |
|
128 | location0 = numpy.dot(mlon.transpose(),mrep) | |
|
129 | ||
|
130 | mlat = numpy.atleast_2d(numpy.arange(nlat)*grid_res + maglimits[1]) | |
|
131 | mrep = numpy.atleast_2d(numpy.zeros(nlon) + 1) | |
|
132 | location1 = numpy.dot(mrep.transpose(),mlat) | |
|
133 | ||
|
134 | location[:,:,0] = location0 | |
|
135 | location[:,:,1] = location1 | |
|
136 | ||
|
137 | alpha = numpy.zeros((nlon,nlat,nhei)) | |
|
138 | rr = numpy.zeros((nlon,nlat,nhei,3)) | |
|
139 | dcos = numpy.zeros((nlon,nlat,nhei,2)) | |
|
140 | ||
|
141 | global first_time | |
|
142 | ||
|
143 | first_time = None | |
|
144 | for ilon in numpy.arange(nlon): | |
|
145 | for ilat in numpy.arange(nlat): | |
|
146 | outs = self.__bdotk(self.heights, | |
|
147 | self.year + self.doy/366., | |
|
148 | coord_site[1], | |
|
149 | coord_site[0], | |
|
150 | coord_site[2], | |
|
151 | coord_site[1]+location[ilon,ilat,1], | |
|
152 | location[ilon,ilat,0]*720./180.) | |
|
153 | ||
|
154 | alpha[ilon, ilat,:] = outs[1] | |
|
155 | rr[ilon, ilat,:,:] = outs[3] | |
|
156 | ||
|
157 | mrep = numpy.atleast_2d((numpy.zeros(nhei)+1)).transpose() | |
|
158 | tmp = outs[3]*numpy.dot(mrep,numpy.atleast_2d(x_ant)) | |
|
159 | tmp = tmp.sum(axis=1) | |
|
160 | dcos[ilon,ilat,:,0] = tmp/numpy.sqrt((outs[3]**2).sum(axis=1)) | |
|
161 | ||
|
162 | mrep = numpy.atleast_2d((numpy.zeros(nhei)+1)).transpose() | |
|
163 | tmp = outs[3]*numpy.dot(mrep,numpy.atleast_2d(y_ant)) | |
|
164 | tmp = tmp.sum(axis=1) | |
|
165 | dcos[ilon,ilat,:,1] = tmp/numpy.sqrt((outs[3]**2).sum(axis=1)) | |
|
166 | ||
|
167 | return dcos, alpha, nlon, nlat | |
|
168 | ||
|
169 | ||
|
170 | def __bdotk(self,heights,tm,gdlat=-11.95,gdlon=-76.8667,gdalt=0.0,decd=-12.88, ham=-4.61666667): | |
|
171 | ||
|
172 | global first_time | |
|
173 | # Mean Earth radius in Km WGS 84 | |
|
174 | a_igrf = 6371.2 | |
|
175 | ||
|
176 | bk = numpy.zeros(heights.size) | |
|
177 | alpha = numpy.zeros(heights.size) | |
|
178 | bfm = numpy.zeros(heights.size) | |
|
179 | rr = numpy.zeros((heights.size,3)) | |
|
180 | rgc = numpy.zeros((heights.size,3)) | |
|
181 | ||
|
182 | ObjGeodetic = Astro_Coords.Geodetic(gdlat,gdalt) | |
|
183 | [gclat,gcalt] = ObjGeodetic.change2geocentric() | |
|
184 | ||
|
185 | gclat = gclat*numpy.pi/180. | |
|
186 | gclon = gdlon*numpy.pi/180. | |
|
187 | ||
|
188 | # Antenna position from center of Earth | |
|
189 | ca_vector = [numpy.cos(gclat)*numpy.cos(gclon),numpy.cos(gclat)*numpy.sin(gclon),numpy.sin(gclat)] | |
|
190 | ca_vector = gcalt*numpy.array(ca_vector) | |
|
191 | ||
|
192 | dec = decd*numpy.pi/180. | |
|
193 | ||
|
194 | # K vector respect to the center of earth. | |
|
195 | klon = gclon + ham*numpy.pi/720. | |
|
196 | k_vector = [numpy.cos(dec)*numpy.cos(klon),numpy.cos(dec)*numpy.sin(klon),numpy.sin(dec)] | |
|
197 | k_vector = numpy.array(k_vector) | |
|
198 | ||
|
199 | for ih in numpy.arange(heights.size): | |
|
200 | # Vector from Earth's center to volume of interest | |
|
201 | rr[ih,:] = k_vector*heights[ih] | |
|
202 | cv_vector = numpy.squeeze(ca_vector) + rr[ih,:] | |
|
203 | ||
|
204 | cv_gcalt = numpy.sqrt(numpy.sum(cv_vector**2.)) | |
|
205 | cvxy = numpy.sqrt(numpy.sum(cv_vector[0:2]**2.)) | |
|
206 | ||
|
207 | radial = cv_vector/cv_gcalt | |
|
208 | east = numpy.array([-1*cv_vector[1],cv_vector[0],0])/cvxy | |
|
209 | comp1 = east[1]*radial[2] - radial[1]*east[2] | |
|
210 | comp2 = east[2]*radial[0] - radial[2]*east[0] | |
|
211 | comp3 = east[0]*radial[1] - radial[0]*east[1] | |
|
212 | north = -1*numpy.array([comp1, comp2, comp3]) | |
|
213 | ||
|
214 | rr_k = cv_vector - numpy.squeeze(ca_vector) | |
|
215 | u_rr = rr_k/numpy.sqrt(numpy.sum(rr_k**2.)) | |
|
216 | ||
|
217 | cv_gclat = numpy.arctan2(cv_vector[2],cvxy) | |
|
218 | cv_gclon = numpy.arctan2(cv_vector[1],cv_vector[0]) | |
|
219 | ||
|
220 | bhei = cv_gcalt-a_igrf | |
|
221 | blat = cv_gclat*180./numpy.pi | |
|
222 | blon = cv_gclon*180./numpy.pi | |
|
223 | bfield = self.__igrfkudeki(bhei,tm,blat,blon) | |
|
224 | ||
|
225 | B = (bfield[0]*north + bfield[1]*east - bfield[2]*radial)*1.0e-5 | |
|
226 | ||
|
227 | bfm[ih] = numpy.sqrt(numpy.sum(B**2.)) #module | |
|
228 | bk[ih] = numpy.sum(u_rr*B) | |
|
229 | alpha[ih] = numpy.arccos(bk[ih]/bfm[ih])*180/numpy.pi | |
|
230 | rgc[ih,:] = numpy.array([cv_gclon, cv_gclat, cv_gcalt]) | |
|
231 | ||
|
232 | return bk, alpha, bfm, rr, rgc | |
|
233 | ||
|
234 | ||
|
235 | def __igrfkudeki(self,heights,time,latitude,longitude,ae=6371.2): | |
|
236 | """ | |
|
237 | __igrfkudeki calculates the International Geomagnetic Reference Field for given in- | |
|
238 | put conditions based on IGRF2005 coefficients. | |
|
239 | ||
|
240 | Parameters | |
|
241 | ---------- | |
|
242 | heights = Scalar or vector giving the height above the Earth of the point in ques- | |
|
243 | tion in kilometers. | |
|
244 | time = Scalar or vector giving the decimal year of time in question (e.g. 1991.2). | |
|
245 | latitude = Latitude of point in question in decimal degrees. Scalar or vector. | |
|
246 | longitude = Longitude of point in question in decimal degrees. Scalar or vector. | |
|
247 | ae = | |
|
248 | first_time = | |
|
249 | ||
|
250 | Return | |
|
251 | ------ | |
|
252 | bn = | |
|
253 | be = | |
|
254 | bd = | |
|
255 | bmod = | |
|
256 | balpha = | |
|
257 | first_time = | |
|
258 | ||
|
259 | Modification History | |
|
260 | -------------------- | |
|
261 | Converted to Python by Freddy R. Galindo, ROJ, 03 October 2009. | |
|
262 | """ | |
|
263 | ||
|
264 | global first_time | |
|
265 | global gs, hs, nvec, mvec, maxcoef | |
|
266 | ||
|
267 | heights = numpy.atleast_1d(heights) | |
|
268 | time = numpy.atleast_1d(time) | |
|
269 | latitude = numpy.atleast_1d(latitude) | |
|
270 | longitude = numpy.atleast_1d(longitude) | |
|
271 | ||
|
272 | if numpy.max(latitude)==90: | |
|
273 | # print "Field calculations are not supported at geographic poles" | |
|
274 | pass | |
|
275 | ||
|
276 | # output arrays | |
|
277 | bn = numpy.zeros(heights.size) | |
|
278 | be = numpy.zeros(heights.size) | |
|
279 | bd = numpy.zeros(heights.size) | |
|
280 | ||
|
281 | if first_time==None:first_time=0 | |
|
282 | ||
|
283 | time0 = time[0] | |
|
284 | if time!=first_time: | |
|
285 | #print "Getting coefficients for", time0 | |
|
286 | [periods,g,h ] = self.__readIGRFcoeff() | |
|
287 | top_year = numpy.max(periods) | |
|
288 | nperiod = (top_year - 1900)/5 + 1 | |
|
289 | ||
|
290 | maxcoef = 10 | |
|
291 | if time0>=2000:maxcoef = 12 | |
|
292 | ||
|
293 | ||
|
294 | # Normalization array for Schmidt fucntions | |
|
295 | multer = numpy.zeros((2+maxcoef,1+maxcoef)) + 1 | |
|
296 | for cn in (numpy.arange(maxcoef)+1): | |
|
297 | for rm in (numpy.arange(cn)+1): | |
|
298 | tmp = numpy.arange(2*rm) + cn - rm + 1. | |
|
299 | multer[rm+1,cn] = ((-1.)**rm)*numpy.sqrt(2./tmp.prod()) | |
|
300 | ||
|
301 | schmidt = multer[1:,1:].transpose() | |
|
302 | ||
|
303 | # n and m arrays | |
|
304 | nvec = numpy.atleast_2d(numpy.arange(maxcoef)+2) | |
|
305 | mvec = numpy.atleast_2d(numpy.arange(maxcoef+1)).transpose() | |
|
306 | ||
|
307 | # Time adjusted igrf g and h with Schmidt normalization | |
|
308 | # IGRF coefficient arrays: g0(n,m), n=1, maxcoeff,m=0, maxcoeff, ... | |
|
309 | if time0<top_year: | |
|
310 | dtime = (time0 - 1900) % 5 | |
|
311 | ntime = (time0 - 1900 - dtime)/5 | |
|
312 | else: | |
|
313 | # Estimating coefficients for times > top_year | |
|
314 | dtime = (time0 - top_year) + 5 | |
|
315 | ntime = g[:,0,0].size - 2 | |
|
316 | ||
|
317 | g0 = g[ntime,1:maxcoef+1,:maxcoef+1] | |
|
318 | h0 = h[ntime,1:maxcoef+1,:maxcoef+1] | |
|
319 | gdot = g[ntime+1,1:maxcoef+1,:maxcoef+1]-g[ntime,1:maxcoef+1,:maxcoef+1] | |
|
320 | hdot = h[ntime+1,1:maxcoef+1,:maxcoef+1]-h[ntime,1:maxcoef+1,:maxcoef+1] | |
|
321 | gs = (g0 + dtime*(gdot/5.))*schmidt[:maxcoef,0:maxcoef+1] | |
|
322 | hs = (h0 + dtime*(hdot/5.))*schmidt[:maxcoef,0:maxcoef+1] | |
|
323 | ||
|
324 | first_time = time0 | |
|
325 | ||
|
326 | for ii in numpy.arange(heights.size): | |
|
327 | # Height dependence array rad = (ae/(ae+height))**(n+3) | |
|
328 | rad = numpy.atleast_2d((ae/(ae + heights[ii]))**(nvec+1)) | |
|
329 | ||
|
330 | # Sin and Cos of m times longitude phi arrays | |
|
331 | mphi = mvec*longitude[ii]*numpy.pi/180. | |
|
332 | cosmphi = numpy.atleast_2d(numpy.cos(mphi)) | |
|
333 | sinmphi = numpy.atleast_2d(numpy.sin(mphi)) | |
|
334 | ||
|
335 | # Cos of colatitude theta | |
|
336 | c = numpy.cos((90 - latitude[ii])*numpy.pi/180.) | |
|
337 | ||
|
338 | # Legendre functions p(n,m|c) | |
|
339 | [p,dp]= scipy.special.lpmn(maxcoef+1,maxcoef+1,c) | |
|
340 | p = p[:,:-1].transpose() | |
|
341 | s = numpy.sqrt((1. - c)*(1 + c)) | |
|
342 | ||
|
343 | # Generate derivative array dpdtheta = -s*dpdc | |
|
344 | dpdtheta = c*p/s | |
|
345 | for m in numpy.arange(maxcoef+2): dpdtheta[:,m] = m*dpdtheta[:,m] | |
|
346 | dpdtheta = dpdtheta + numpy.roll(p,-1,axis=1) | |
|
347 | ||
|
348 | # Extracting arrays required for field calculations | |
|
349 | p = p[1:maxcoef+1,:maxcoef+1] | |
|
350 | dpdtheta = dpdtheta[1:maxcoef+1,:maxcoef+1] | |
|
351 | ||
|
352 | # Weigh p and dpdtheta with gs and hs coefficients. | |
|
353 | gp = gs*p | |
|
354 | hp = hs*p | |
|
355 | gdpdtheta = gs*dpdtheta | |
|
356 | hdpdtheta = hs*dpdtheta | |
|
357 | # Calcultate field components | |
|
358 | matrix0 = numpy.dot(gdpdtheta,cosmphi) | |
|
359 | matrix1 = numpy.dot(hdpdtheta,sinmphi) | |
|
360 | bn[ii] = numpy.dot(rad,(matrix0 + matrix1)) | |
|
361 | matrix0 = numpy.dot(hp,(mvec*cosmphi)) | |
|
362 | matrix1 = numpy.dot(gp,(mvec*sinmphi)) | |
|
363 | be[ii] = numpy.dot((-1*rad),((matrix0 - matrix1)/s)) | |
|
364 | matrix0 = numpy.dot(gp,cosmphi) | |
|
365 | matrix1 = numpy.dot(hp,sinmphi) | |
|
366 | bd[ii] = numpy.dot((-1*nvec*rad),(matrix0 + matrix1)) | |
|
367 | ||
|
368 | bmod = numpy.sqrt(bn**2. + be**2. + bd**2.) | |
|
369 | btheta = numpy.arctan(bd/numpy.sqrt(be**2. + bn**2.))*180/numpy.pi | |
|
370 | balpha = numpy.arctan(be/bn)*180./numpy.pi | |
|
371 | ||
|
372 | #bn : north | |
|
373 | #be : east | |
|
374 | #bn : radial | |
|
375 | #bmod : module | |
|
376 | ||
|
377 | ||
|
378 | return bn, be, bd, bmod, btheta, balpha | |
|
379 | ||
|
380 | def str2num(self, datum): | |
|
381 | try: | |
|
382 | return int(datum) | |
|
383 | except: | |
|
384 | try: | |
|
385 | return float(datum) | |
|
386 | except: | |
|
387 | return datum | |
|
388 | ||
|
389 | def __readIGRFfile(self, filename): | |
|
390 | list_years=[] | |
|
391 | for i in range(1,24): | |
|
392 | list_years.append(1895.0 + i*5) | |
|
393 | ||
|
394 | epochs=list_years | |
|
395 | epochs.append(epochs[-1]+5) | |
|
396 | nepochs = numpy.shape(epochs) | |
|
397 | ||
|
398 | gg = numpy.zeros((13,14,nepochs[0]),dtype=float) | |
|
399 | hh = numpy.zeros((13,14,nepochs[0]),dtype=float) | |
|
400 | ||
|
401 | coeffs_file=open(filename) | |
|
402 | lines=coeffs_file.readlines() | |
|
403 | ||
|
404 | coeffs_file.close() | |
|
405 | ||
|
406 | for line in lines: | |
|
407 | items = line.split() | |
|
408 | g_h = items[0] | |
|
409 | n = self.str2num(items[1]) | |
|
410 | m = self.str2num(items[2]) | |
|
411 | ||
|
412 | coeffs = items[3:] | |
|
413 | ||
|
414 | for i in range(len(coeffs)-1): | |
|
415 | coeffs[i] = self.str2num(coeffs[i]) | |
|
416 | ||
|
417 | #coeffs = numpy.array(coeffs) | |
|
418 | ncoeffs = numpy.shape(coeffs)[0] | |
|
419 | ||
|
420 | if g_h == 'g': | |
|
421 | # print n," g ",m | |
|
422 | gg[n-1,m,:]=coeffs | |
|
423 | elif g_h=='h': | |
|
424 | # print n," h ",m | |
|
425 | hh[n-1,m,:]=coeffs | |
|
426 | # else : | |
|
427 | # continue | |
|
428 | ||
|
429 | # Ultimo Reordenamiento para almacenar . | |
|
430 | gg[:,:,nepochs[0]-1] = gg[:,:,nepochs[0]-2] + 5*gg[:,:,nepochs[0]-1] | |
|
431 | hh[:,:,nepochs[0]-1] = hh[:,:,nepochs[0]-2] + 5*hh[:,:,nepochs[0]-1] | |
|
432 | ||
|
433 | # return numpy.array([gg,hh]) | |
|
434 | periods = numpy.array(epochs) | |
|
435 | g = gg | |
|
436 | h = hh | |
|
437 | return periods, g, h | |
|
438 | ||
|
439 | ||
|
440 | def __readIGRFcoeff(self,filename="igrf10coeffs.dat"): | |
|
441 | """ | |
|
442 | __readIGRFcoeff reads the coefficients from a binary file which is located in the | |
|
443 | folder "resource." | |
|
444 | ||
|
445 | Parameter | |
|
446 | --------- | |
|
447 | filename = A string to specify the name of the file which contains thec coeffs. The | |
|
448 | default value is "igrf10coeffs.dat" | |
|
449 | ||
|
450 | Return | |
|
451 | ------ | |
|
452 | periods = A lineal array giving... | |
|
453 | g1 = | |
|
454 | h1 = | |
|
455 | ||
|
456 | Modification History | |
|
457 | -------------------- | |
|
458 | Converted to Python by Freddy R. Galindo, ROJ, 03 October 2009. | |
|
459 | """ | |
|
460 | ||
|
461 | # # igrfile = sys.path[-1] + os.sep + "resource" + os.sep + filename | |
|
462 | # igrfile = os.path.join('./resource',filename) | |
|
463 | # f = open(igrfile,'rb') | |
|
464 | # #f = open(os.getcwd() + os.sep + "resource" + os.sep + filename,'rb') | |
|
465 | # | |
|
466 | # # Reading SkyNoise Power (lineal scale) | |
|
467 | # periods = numpy.fromfile(f,numpy.dtype([('var','<f4')]),23) | |
|
468 | # periods = periods['var'] | |
|
469 | # | |
|
470 | # g = numpy.fromfile(f,numpy.dtype([('var','<f8')]),23*14*14) | |
|
471 | # g = g['var'].reshape((14,14,23)).transpose() | |
|
472 | # | |
|
473 | # h = numpy.fromfile(f,numpy.dtype([('var','<f8')]),23*14*14) | |
|
474 | # h = h['var'].reshape((14,14,23)).transpose() | |
|
475 | # | |
|
476 | # f.close() | |
|
477 | base_path = os.path.dirname(os.path.abspath(__file__)) | |
|
478 | filename = os.path.join(base_path, "igrf11coeffs.txt") | |
|
479 | ||
|
480 | period_v, g_v, h_v = self.__readIGRFfile(filename) | |
|
481 | g2 = numpy.zeros((14,14,24)) | |
|
482 | h2 = numpy.zeros((14,14,24)) | |
|
483 | g2[1:14,:,:] = g_v | |
|
484 | h2[1:14,:,:] = h_v | |
|
485 | ||
|
486 | g = numpy.transpose(g2, (2,0,1)) | |
|
487 | h = numpy.transpose(h2, (2,0,1)) | |
|
488 | periods = period_v.copy() | |
|
489 | ||
|
490 | return periods, g, h | |
|
491 | ||
|
492 | def rotvector(self,vector,axis=1,ang=0): | |
|
493 | """ | |
|
494 | rotvector function returns the new vector generated rotating the rectagular coords. | |
|
495 | ||
|
496 | Parameters | |
|
497 | ---------- | |
|
498 | vector = A lineal 3-elements array (x,y,z). | |
|
499 | axis = A integer to specify the axis used to rotate the coord systems. The default | |
|
500 | value is 1. | |
|
501 | axis = 1 -> Around "x" | |
|
502 | axis = 2 -> Around "y" | |
|
503 | axis = 3 -> Around "z" | |
|
504 | ang = Angle of rotation (in radians). The default value is zero. | |
|
505 | ||
|
506 | Return | |
|
507 | ------ | |
|
508 | rotvector = A lineal array of 3 elements giving the new coordinates. | |
|
509 | ||
|
510 | Modification History | |
|
511 | -------------------- | |
|
512 | Converted to Python by Freddy R. Galindo, ROJ, 01 October 2009. | |
|
513 | """ | |
|
514 | ||
|
515 | if axis==1: | |
|
516 | t = [[1,0,0],[0,numpy.cos(ang),numpy.sin(ang)],[0,-numpy.sin(ang),numpy.cos(ang)]] | |
|
517 | elif axis==2: | |
|
518 | t = [[numpy.cos(ang),0,-numpy.sin(ang)],[0,1,0],[numpy.sin(ang),0,numpy.cos(ang)]] | |
|
519 | elif axis==3: | |
|
520 | t = [[numpy.cos(ang),numpy.sin(ang),0],[-numpy.sin(ang),numpy.cos(ang),0],[0,0,1]] | |
|
521 | ||
|
522 | rotvector = numpy.array(numpy.dot(numpy.array(t),numpy.array(vector))) | |
|
523 | ||
|
524 | return rotvector | |
|
525 | ||
|
526 | class AntPatternPlot: | |
|
527 | ||
|
528 | def __init__(self): | |
|
529 | """ | |
|
530 | AntPatternPlot creates an object to call methods to plot the antenna pattern. | |
|
531 | ||
|
532 | Modification History | |
|
533 | -------------------- | |
|
534 | Created by Freddy Galindo, ROJ, 06 October 2009. | |
|
535 | """ | |
|
536 | ||
|
537 | self.fig = Figure(figsize=(8,8), facecolor='white') | |
|
538 | self.ax = self.fig.add_subplot(111) | |
|
539 | ||
|
540 | def contPattern(self,iplot=0,gpath='',filename='',mesg='',amp=None ,x=None ,y=None ,getCut=None,title='', save=False): | |
|
541 | """ | |
|
542 | contPattern plots a contour map of the antenna pattern. | |
|
543 | ||
|
544 | Parameters | |
|
545 | ---------- | |
|
546 | iplot = A integer to specify if the plot is the first, second, ... The default va- | |
|
547 | lue is 0. | |
|
548 | ||
|
549 | Examples | |
|
550 | -------- | |
|
551 | >> Over_Jro.JroPattern(pattern=2).contPattern() | |
|
552 | ||
|
553 | Modification history | |
|
554 | -------------------- | |
|
555 | Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009. | |
|
556 | """ | |
|
557 | ||
|
558 | if getCut == 1: | |
|
559 | return | |
|
560 | ||
|
561 | xmax = numpy.max(x) | |
|
562 | xmin = numpy.min(x) | |
|
563 | ymax = numpy.max(y) | |
|
564 | ymin = numpy.min(y) | |
|
565 | ||
|
566 | levels = numpy.array([1e-3,1e-2,1e-1,0.5,1.0]) | |
|
567 | tmp = numpy.round(10*numpy.log10(levels),decimals=1) | |
|
568 | labels = list(range(5)) | |
|
569 | for i in numpy.arange(5): | |
|
570 | labels[i] = str(numpy.int(tmp[i])) | |
|
571 | ||
|
572 | ||
|
573 | colors = ((0,0,1.),(0,170/255.,0),(127/255.,1.,0),(1.,109/255.,0),(128/255.,0,0)) | |
|
574 | CS = self.ax.contour(x,y,amp.transpose(),levels,colors=colors) | |
|
575 | fmt = {} | |
|
576 | for l,s in zip(CS.levels,labels): | |
|
577 | fmt[l] = s | |
|
578 | ||
|
579 | self.ax.annotate('Ng',xy=(-0.05,1.04),xytext=(0.01,0.962),xycoords='axes fraction',arrowprops=dict(facecolor='black', width=1.,shrink=0.2),fontsize=15.) | |
|
580 | self.ax.annotate(mesg,xy=(0,0),xytext=(0.01,0.01),xycoords='figure fraction') | |
|
581 | self.ax.clabel(CS,CS.levels,inline=True,fmt=fmt,fontsize=10) | |
|
582 | self.ax.set_xlim(xmin,xmax) | |
|
583 | self.ax.set_ylim(ymin,ymax) | |
|
584 | self.ax.set_title("Total Pattern: " + title) | |
|
585 | self.ax.set_xlabel("West to South") | |
|
586 | self.ax.set_ylabel("West to North") | |
|
587 | self.ax.grid(True) | |
|
588 | ||
|
589 | ||
|
590 | def plotRaDec(self,gpath=None,filename=None,jd=2452640.5,ra_obs=None,xg=None,yg=None,x=None,y=None, save=True): | |
|
591 | """ | |
|
592 | plotRaDec draws right ascension and declination lines on a JRO plane. This function | |
|
593 | must call after conPattern. | |
|
594 | ||
|
595 | Parameters | |
|
596 | ---------- | |
|
597 | jd = A scalar giving the Julian date. | |
|
598 | ra_obs = Scalar giving the right ascension of the observatory. | |
|
599 | xg = A 3-element array to specify .. | |
|
600 | yg = A 3-element array to specify .. | |
|
601 | ||
|
602 | Examples | |
|
603 | -------- | |
|
604 | >> Over_Jro.JroPattern(pattern=2).contPattern() | |
|
605 | >> Over_Jro.JroPattern(pattern=2).plotRaDec(jd=jd,ra_obs=ra_obs,xg=xg,yg=yg) | |
|
606 | ||
|
607 | Modification history | |
|
608 | -------------------- | |
|
609 | Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009. | |
|
610 | """ | |
|
611 | ||
|
612 | # Finding RA of observatory for a specific date | |
|
613 | if ra_obs is None:ra_obs = numpy.array([23.37060849]) | |
|
614 | if xg is None:xg = numpy.array([0.62918474,-0.77725579,0.]) | |
|
615 | if yg is None:yg = numpy.array([0.77700346,0.62898048,0.02547905]) | |
|
616 | ||
|
617 | # Getting HA and DEC axes | |
|
618 | mindec = -28; maxdec = 4; incdec = 2. | |
|
619 | ndec = numpy.int((maxdec - mindec)/incdec) + 1 | |
|
620 | ||
|
621 | minha = -20; maxha = 20; incha = 2. | |
|
622 | nha = numpy.int((maxha - minha)/incha) + 1 | |
|
623 | ||
|
624 | #mcosx = numpy.zeros((nha,ndec)) | |
|
625 | #mcosy = numpy.zeros((nha,ndec)) | |
|
626 | ||
|
627 | ha_axes = numpy.reshape(numpy.arange(nha)*incha + minha,(nha,1)) | |
|
628 | ones_dec = numpy.reshape(numpy.zeros(ndec) + 1,(ndec,1)) | |
|
629 | ha_axes = numpy.dot(ha_axes,ones_dec.transpose()) | |
|
630 | ha_axes2 = numpy.array(ra_obs - ha_axes) | |
|
631 | ||
|
632 | dec_axes = numpy.reshape(numpy.arange(ndec)*incdec + mindec,(ndec,1)) | |
|
633 | ones_ra = numpy.reshape(numpy.zeros(nha) + 1,(nha,1)) | |
|
634 | dec_axes = numpy.dot(ones_ra,dec_axes.transpose()) | |
|
635 | dec_axes2 = numpy.array(dec_axes) | |
|
636 | ||
|
637 | ObjHor = Astro_Coords.Equatorial(ha_axes2,dec_axes2,jd) | |
|
638 | [alt,az,ha] = ObjHor.change2AltAz() | |
|
639 | ||
|
640 | z = numpy.transpose(alt)*Misc_Routines.CoFactors.d2r ; z = z.flatten() | |
|
641 | az = numpy.transpose(az)*Misc_Routines.CoFactors.d2r ; az = az.flatten() | |
|
642 | ||
|
643 | vect = numpy.array([numpy.cos(z)*numpy.sin(az),numpy.cos(z)*numpy.cos(az),numpy.sin(z)]) | |
|
644 | ||
|
645 | xg = numpy.atleast_2d(xg) | |
|
646 | dcosx = numpy.array(numpy.dot(xg,vect)) | |
|
647 | yg = numpy.atleast_2d(yg) | |
|
648 | dcosy = numpy.array(numpy.dot(yg,vect)) | |
|
649 | ||
|
650 | mcosx = dcosx.reshape(ndec,nha) | |
|
651 | mcosy = dcosy.reshape(ndec,nha) | |
|
652 | ||
|
653 | # Defining NAN for points outof limits. | |
|
654 | xmax = numpy.max(x) | |
|
655 | xmin = numpy.min(x) | |
|
656 | ymax = numpy.max(y) | |
|
657 | ymin = numpy.min(y) | |
|
658 | ||
|
659 | factor = 1.3 | |
|
660 | noval = numpy.where((mcosx>(xmax*factor)) | (mcosx<(xmin*factor))) | |
|
661 | if noval[0].size>0:mcosx[noval] = numpy.nan | |
|
662 | noval = numpy.where((mcosy>(ymax*factor)) | (mcosy<(ymin*factor))) | |
|
663 | if noval[0].size>0:mcosy[noval] = numpy.nan | |
|
664 | ||
|
665 | # Plotting HA and declination grid. | |
|
666 | iha0 = numpy.int((0 - minha)/incha) | |
|
667 | idec0 = numpy.int((-14 - mindec)/incdec) | |
|
668 | ||
|
669 | colorgrid = (1.,109/255.,0) | |
|
670 | self.ax.plot(mcosx.transpose(),mcosy.transpose(),color=colorgrid,linestyle='--', lw=0.5) | |
|
671 | for idec in numpy.arange(ndec): | |
|
672 | if idec != idec0: | |
|
673 | valx = (mcosx[idec,iha0]<=xmax) & (mcosx[idec,iha0]>=xmin) | |
|
674 | valy = (mcosy[idec,iha0]<=ymax) & (mcosy[idec,iha0]>=ymin) | |
|
675 | if valx & valy: | |
|
676 | text = str(numpy.int(mindec + incdec*idec))+'$^o$' | |
|
677 | self.ax.text(mcosx[idec,iha0],mcosy[idec,iha0],text) | |
|
678 | ||
|
679 | self.ax.plot(mcosx,mcosy,color=colorgrid,linestyle='--',lw=0.5) | |
|
680 | for iha in numpy.arange(nha): | |
|
681 | if iha != iha0: | |
|
682 | valx = (mcosx[idec0,iha]<=xmax) & (mcosx[idec0,iha]>=xmin) | |
|
683 | valy = (mcosy[idec0,iha]<=ymax) & (mcosy[idec0,iha]>=ymin) | |
|
684 | if valx & valy: | |
|
685 | text = str(4*numpy.int(minha + incha*iha))+"'" | |
|
686 | self.ax.text(mcosx[idec0,iha],mcosy[idec0,iha],text) | |
|
687 | ||
|
688 | if save: | |
|
689 | save_fig = os.path.join(gpath,filename) | |
|
690 | self.fig.savefig(save_fig,format='png') | |
|
691 | ||
|
692 | ||
|
693 | def plotBField(self,gpath,filename,dcos,alpha, nlon, nlat, dcosxrange, dcosyrange, heights, alpha_i, save=False): | |
|
694 | """ | |
|
695 | plotBField draws the magnetic field in a directional cosines plot. | |
|
696 | ||
|
697 | Parameters | |
|
698 | ---------- | |
|
699 | dcos = An 4-dimensional array giving the directional cosines of the magnetic field | |
|
700 | over the desired place. | |
|
701 | alpha = An 3-dimensional array giving the angle of the magnetic field over the desi- | |
|
702 | red place. | |
|
703 | nlon = An integer to specify the number of elements per longitude. | |
|
704 | nlat = An integer to specify the number of elements per latitude. | |
|
705 | dcosxrange = A 2-element array giving the range of the directional cosines in the | |
|
706 | "x" axis. | |
|
707 | dcosyrange = A 2-element array giving the range of the directional cosines in the | |
|
708 | "y" axis. | |
|
709 | heights = An array giving the heights (km) where the magnetic field will be modeled By default the magnetic field will be computed at 100, 500 and 1000km. | |
|
710 | alpha_i = Angle to interpolate the magnetic field. | |
|
711 | Modification History | |
|
712 | -------------------- | |
|
713 | Converted to Python by Freddy R. Galindo, ROJ, 07 October 2009. | |
|
714 | """ | |
|
715 | ||
|
716 | handles = [] | |
|
717 | objects = [] | |
|
718 | colors = ['k','m','c','b','g','r','y'] | |
|
719 | marker = ['-+','-*','-D','-x','-s','->','-o','-^'] | |
|
720 | ||
|
721 | alpha_location = numpy.zeros((nlon,2,heights.size)) | |
|
722 | ||
|
723 | for ih in numpy.arange(heights.size): | |
|
724 | alpha_location[:,0,ih] = dcos[:,0,ih,0] | |
|
725 | for ilon in numpy.arange(nlon): | |
|
726 | myx = (alpha[ilon,:,ih])[::-1] | |
|
727 | myy = (dcos[ilon,:,ih,0])[::-1] | |
|
728 | tck = scipy.interpolate.splrep(myx,myy,s=0) | |
|
729 | mydcosx = scipy.interpolate.splev(alpha_i,tck,der=0) | |
|
730 | ||
|
731 | myx = (alpha[ilon,:,ih])[::-1] | |
|
732 | myy = (dcos[ilon,:,ih,1])[::-1] | |
|
733 | tck = scipy.interpolate.splrep(myx,myy,s=0) | |
|
734 | mydcosy = scipy.interpolate.splev(alpha_i,tck,der=0) | |
|
735 | alpha_location[ilon,:,ih] = numpy.array([mydcosx, mydcosy]) | |
|
736 | ||
|
737 | ||
|
738 | ObjFig, = self.ax.plot(alpha_location[:,0,ih],alpha_location[:,1,ih], | |
|
739 | marker[ih % 8],color=colors[numpy.int(ih/8)],ms=4.5,lw=0.5) | |
|
740 | handles.append(ObjFig) | |
|
741 | objects.append(numpy.str(heights[ih]) + ' km') | |
|
742 | ||
|
743 | legend = self.ax.legend(handles, objects,loc="lower right", numpoints=1, handlelength=0.3, | |
|
744 | handletextpad=0.02, borderpad=0.3, labelspacing=0.1) | |
|
745 | self.ax.add_artist(legend) | |
|
746 | ||
|
747 | def plotCelestial(self, jd, main_dec, tod, maxha_min, objects, glat, glon, xg, yg, dcosxrange, dcosyrange): | |
|
748 | ||
|
749 | self.tod = tod | |
|
750 | ||
|
751 | self.dcosx_sun = 1 | |
|
752 | self.dcosy_sun = 1 | |
|
753 | self.ha_sun = 1 | |
|
754 | self.time_sun = 1 | |
|
755 | ||
|
756 | self.dcosx_moon = 1 | |
|
757 | self.dcosy_moon = 1 | |
|
758 | self.ha_moon = 1 | |
|
759 | self.time_moon = 1 | |
|
760 | ||
|
761 | self.dcosx_hydra = 1 | |
|
762 | self.dcosy_hydra = 1 | |
|
763 | self.ha_hydra = 1 | |
|
764 | self.time_hydra = 1 | |
|
765 | ||
|
766 | self.dcosx_galaxy = 1 | |
|
767 | self.dcosy_galaxy = 1 | |
|
768 | self.ha_galaxy = 1 | |
|
769 | self.time_galaxy = 1 | |
|
770 | ||
|
771 | tod = self.tod | |
|
772 | ||
|
773 | maxlev = 24; minlev = 0; maxcol = 39; mincol = 10 | |
|
774 | handles = [] | |
|
775 | titles = ['$Sun$','$Moon$','$Hydra$','$Galaxy$'] | |
|
776 | marker = ['--^','--s','--*','--o'] | |
|
777 | ||
|
778 | # Getting RGB table to plot celestial object over Jicamarca | |
|
779 | colortable = ['olive', 'indigo', 'cyan', 'red'] | |
|
780 | labels = [] | |
|
781 | for io in objects: | |
|
782 | ObjBodies = Astro_Coords.CelestialBodies() | |
|
783 | if io==0: | |
|
784 | [ra,dec,sunlon,sunobliq] = ObjBodies.sunpos(jd) | |
|
785 | elif io==1: | |
|
786 | [ra,dec,dist,moonlon,moonlat] = ObjBodies.moonpos(jd) | |
|
787 | elif io==2: | |
|
788 | [ra,dec] = ObjBodies.hydrapos() | |
|
789 | elif io==3: | |
|
790 | [maxra,ra] = ObjBodies.skynoise_jro(dec_cut=main_dec) | |
|
791 | ra = maxra*15. | |
|
792 | dec = main_dec | |
|
793 | ||
|
794 | ObjEq = Astro_Coords.Equatorial(ra, dec, jd, lat=glat, lon=glon) | |
|
795 | [alt, az, ha] = ObjEq.change2AltAz() | |
|
796 | vect = numpy.array([az,alt]).transpose() | |
|
797 | vect = Misc_Routines.Vector(vect,direction=0).Polar2Rect() | |
|
798 | ||
|
799 | dcosx = numpy.array(numpy.dot(vect,xg)) | |
|
800 | dcosy = numpy.array(numpy.dot(vect,yg)) | |
|
801 | wrap = numpy.where(ha>=180.) | |
|
802 | ||
|
803 | if wrap[0].size>0: | |
|
804 | ha[wrap] = ha[wrap] - 360. | |
|
805 | ||
|
806 | val = numpy.where((numpy.abs(ha))<=(maxha_min*0.25)) | |
|
807 | ||
|
808 | if val[0].size>2: | |
|
809 | tod_1 = tod*1. | |
|
810 | shift_1 = numpy.where(tod>12.) | |
|
811 | tod_1[shift_1] = tod_1[shift_1] - 24. | |
|
812 | tod_2 = tod*1. | |
|
813 | shift_2 = numpy.where(tod<12.) | |
|
814 | tod_2[shift_2] = tod_2[shift_2] + 24. | |
|
815 | ||
|
816 | diff0 = numpy.nanmax(tod[val]) - numpy.nanmin(tod[val]) | |
|
817 | diff1 = numpy.nanmax(tod_1[val]) - numpy.nanmin(tod_1[val]) | |
|
818 | diff2 = numpy.nanmax(tod_2[val]) - numpy.nanmin(tod_2[val]) | |
|
819 | ||
|
820 | if ((diff0<=diff1) & (diff0<=diff2)): | |
|
821 | tod_0 = tod | |
|
822 | elif ((diff1<diff0) & (diff1<diff2)): | |
|
823 | tod_0 = tod_1 | |
|
824 | else: | |
|
825 | tod_0 = tod_2 | |
|
826 | ||
|
827 | if io==0: | |
|
828 | self.dcosx_sun = dcosx[val] | |
|
829 | self.dcosy_sun = dcosy[val] | |
|
830 | self.ha_sun = ha[val] | |
|
831 | self.time_sun = numpy.median(tod_0[val]) | |
|
832 | elif io==1: | |
|
833 | self.dcosx_moon = dcosx[val] | |
|
834 | self.dcosy_moon = dcosy[val] | |
|
835 | self.ha_moon = ha[val] | |
|
836 | self.time_moon = numpy.median(tod_0[val]) | |
|
837 | elif io==2: | |
|
838 | self.dcosx_hydra = dcosx[val] | |
|
839 | self.dcosy_hydra = dcosy[val] | |
|
840 | self.ha_hydra = ha[val] | |
|
841 | self.time_hydra = numpy.mean(tod_0[val]) | |
|
842 | elif io==3: | |
|
843 | self.dcosx_galaxy = dcosx[val] | |
|
844 | self.dcosy_galaxy = dcosy[val] | |
|
845 | self.ha_galaxy = ha[val] | |
|
846 | self.time_galaxy = numpy.mean(tod_0[val]) | |
|
847 | ||
|
848 | index = numpy.mean(tod_0[val]) - minlev | |
|
849 | index = (index*(maxcol - mincol)/(maxlev - minlev)) + mincol | |
|
850 | index = numpy.int(index) | |
|
851 | figobjects, = self.ax.plot(dcosx[val],dcosy[val],marker[io],\ | |
|
852 | lw=1,ms=7,mew=0,color=colortable[io]) | |
|
853 | handles.append(figobjects) | |
|
854 | labels.append(titles[io]) | |
|
855 | ||
|
856 | ||
|
857 | legend = self.ax.legend(handles,labels,loc="lower left", numpoints=1, handlelength=0.5, \ | |
|
858 | borderpad=0.3, handletextpad=0.1,labelspacing=0.2,fontsize='small') | |
|
859 | ||
|
860 | self.ax.add_artist(legend) | |
|
861 | ||
|
862 | ||
|
863 | class JroPattern(): | |
|
864 | def __init__(self,pattern=0,path=None,filename=None,nptsx=101,nptsy=101,maxphi=5,fftopt=0, \ | |
|
865 | getcut=0,dcosx=None,dcosy=None,eomwl=6,airwl=4, **kwargs): | |
|
866 | """ | |
|
867 | JroPattern class creates an object to represent the useful parameters for beam mode- | |
|
868 | lling of the Jicamarca VHF radar. | |
|
869 | ||
|
870 | Parameters | |
|
871 | ---------- | |
|
872 | pattern = An integer (See JroAntSetup to know the available values) to load a prede- | |
|
873 | fined configuration. The default value is 0. To use a user-defined configuration | |
|
874 | pattern must be None. | |
|
875 | path = A string giving the directory that contains the user-configuration file. PATH | |
|
876 | will work if pattern is None. | |
|
877 | filename = A string giving the name of the user-configuration file. FILENAME will | |
|
878 | work if pattern is None. | |
|
879 | nptsx = A scalar to specify the number of points used to define the angular resolu- | |
|
880 | tion in the "x" axis. The default value is 101. | |
|
881 | nptsy = A scalar to specify the number of points used to define the angular resolu- | |
|
882 | tion in the "x" axis. The default value is 101. | |
|
883 | maxphi = A scalar giving the maximum (absolute) angle (in degree) to model the ante- | |
|
884 | nna pattern. The default value is 5 degrees. | |
|
885 | fftopt = Set this input to 1 to model the beam using FFT. To model using antenna | |
|
886 | theory set to 0 (default value). | |
|
887 | getcut = Set to 1 to show an antenna cut instead of a contour plot of itself (set to | |
|
888 | 0). The defautl value is 0. | |
|
889 | dcosx = An array giving the directional cosines for the x-axis. DCOSX will work if | |
|
890 | getcut is actived. | |
|
891 | dcosy = An array giving the directional cosines for the y-axis. DCOSY will work if | |
|
892 | getcut is actived. | |
|
893 | eomwl = A scalar giving the radar wavelength. The default value is 6m (50 MHZ). | |
|
894 | airwl = Set this input to float (or intger) to specify the wavelength (in meters) of | |
|
895 | the transmitted EOM wave in the air. The default value is 4m. | |
|
896 | ||
|
897 | Modification History | |
|
898 | -------------------- | |
|
899 | Converted to Object-oriented Programming by Freddy Galindo, ROJ, 20 September 2009. | |
|
900 | """ | |
|
901 | ||
|
902 | ||
|
903 | ||
|
904 | # Getting antenna configuration. | |
|
905 | if filename: | |
|
906 | setup = JroAntSetup.ReturnSetup(path=path,filename=filename,pattern=pattern) | |
|
907 | ||
|
908 | ues = setup["ues"] | |
|
909 | phase = setup["phase"] | |
|
910 | gaintx = setup["gaintx"] | |
|
911 | gainrx = setup["gainrx"] | |
|
912 | justrx = setup["justrx"] | |
|
913 | self.title = setup["title"] | |
|
914 | else: | |
|
915 | ues = kwargs["ues"] | |
|
916 | phase = kwargs["phases"] | |
|
917 | gaintx = kwargs["gain_tx"] | |
|
918 | gainrx = kwargs["gain_rx"] | |
|
919 | justrx = kwargs["just_rx"] | |
|
920 | self.title = kwargs.get("title", "JRO Pattern") | |
|
921 | ||
|
922 | # Defining attributes for JroPattern class. | |
|
923 | # Antenna configuration | |
|
924 | ||
|
925 | self.uestx = ues | |
|
926 | self.phasetx = phase | |
|
927 | self.gaintx = gaintx | |
|
928 | self.uesrx = ues | |
|
929 | self.phaserx = phase | |
|
930 | self.gainrx = gainrx | |
|
931 | self.justrx = justrx | |
|
932 | ||
|
933 | # Pattern resolution & method to model | |
|
934 | self.maxphi = maxphi | |
|
935 | self.nptsx = nptsx | |
|
936 | self.nptsy = nptsy | |
|
937 | self.fftopt = fftopt | |
|
938 | ||
|
939 | # To get a cut of the pattern. | |
|
940 | self.getcut = getcut | |
|
941 | ||
|
942 | maxdcos = numpy.sin(maxphi*Misc_Routines.CoFactors.d2r) | |
|
943 | if dcosx==None:dcosx = ((numpy.arange(nptsx,dtype=float)/(nptsx-1))-0.5)*2*maxdcos | |
|
944 | if dcosy==None:dcosy = ((numpy.arange(nptsy,dtype=float)/(nptsy-1))-0.5)*2*maxdcos | |
|
945 | self.dcosx = dcosx | |
|
946 | self.dcosy = dcosy | |
|
947 | self.nx = dcosx.size | |
|
948 | self.ny = dcosy.size*(getcut==0) + (getcut==1) | |
|
949 | ||
|
950 | self.eomwl = eomwl | |
|
951 | self.airwl = airwl | |
|
952 | ||
|
953 | self.kk = 2.*numpy.pi/eomwl | |
|
954 | ||
|
955 | self.pattern = None | |
|
956 | self.meanpos = None | |
|
957 | self.norpattern = None | |
|
958 | self.maxpattern = None | |
|
959 | ||
|
960 | ||
|
961 | ||
|
962 | self.getPattern() | |
|
963 | ||
|
964 | def getPattern(self): | |
|
965 | """ | |
|
966 | getpattern method returns the modeled total antenna pattern and its mean position. | |
|
967 | ||
|
968 | Return | |
|
969 | ------ | |
|
970 | pattern = An array giving the Modelled antenna pattern. | |
|
971 | mean_pos = A 2-elements array giving the mean position of the main beam. | |
|
972 | ||
|
973 | Examples | |
|
974 | -------- | |
|
975 | >> [pattern, mean_pos] = JroPattern(pattern=2).getPattern() | |
|
976 | >> print meanpos | |
|
977 | [ 8.08728085e-14 -4.78193873e-14] | |
|
978 | ||
|
979 | Modification history | |
|
980 | -------------------- | |
|
981 | Developed by Jorge L. Chau. | |
|
982 | Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009. | |
|
983 | """ | |
|
984 | ||
|
985 | if (self.fftopt>0) and (self.getcut>0): | |
|
986 | #print "Conflict bewteen fftopt and getcut" | |
|
987 | #print "To get a cut of the antenna pattern uses ffopt=0" | |
|
988 | return None, None | |
|
989 | ||
|
990 | if (self.fftopt==0): | |
|
991 | # Getting antenna pattern using the array method | |
|
992 | self.pattern = self.__usingArray(rx=1) | |
|
993 | if (self.justrx==0):self.pattern = self.pattern*self.__usingArray(rx=0) | |
|
994 | ||
|
995 | elif (self.fftopt>0): | |
|
996 | # Getting antenna pattern using FFT method | |
|
997 | self.pattern = self.__usingFFT(rx=1) | |
|
998 | if (self.justrx==0):self.pattern = self.pattern*self.__usingFFT(rx=0) | |
|
999 | ||
|
1000 | self.maxpattern = numpy.nanmax(self.pattern) | |
|
1001 | self.norpattern = self.pattern/self.maxpattern | |
|
1002 | if self.getcut==0:self.__getBeamPars() | |
|
1003 | ||
|
1004 | def __usingArray(self,rx): | |
|
1005 | """ | |
|
1006 | __usingArray method returns the Jicamarca antenna pattern computed using array model | |
|
1007 | ||
|
1008 | pattern = dipolepattern x modulepattern | |
|
1009 | ||
|
1010 | Parameters | |
|
1011 | ---------- | |
|
1012 | rx = Set to 1 to use the Rx information. Otherwise set to 0 for Tx. | |
|
1013 | ||
|
1014 | Return | |
|
1015 | ------ | |
|
1016 | pattern = An array giving the modelled antenna pattern using the array model. | |
|
1017 | ||
|
1018 | Modification history | |
|
1019 | -------------------- | |
|
1020 | Developed by Jorge L. Chau. | |
|
1021 | Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009. | |
|
1022 | """ | |
|
1023 | ||
|
1024 | if rx==1: | |
|
1025 | ues = self.uesrx | |
|
1026 | phase = self.phaserx | |
|
1027 | gain = self.gainrx | |
|
1028 | elif rx==0: | |
|
1029 | ues = self.uestx | |
|
1030 | phase = self.phasetx | |
|
1031 | gain = self.gaintx | |
|
1032 | ||
|
1033 | ues = ues*360./self.airwl | |
|
1034 | phase = phase*360./self.airwl | |
|
1035 | ||
|
1036 | for ii in range(4): | |
|
1037 | if ii==0:dim = numpy.array([4,0,8,4]) # WEST | |
|
1038 | elif ii==1:dim = numpy.array([0,0,4,4]) # NORTH | |
|
1039 | elif ii==2:dim = numpy.array([0,4,4,8]) # EAST | |
|
1040 | elif ii==3:dim = numpy.array([4,4,8,8]) # SOUTH | |
|
1041 | xi = dim[0]; xf = dim[2]; yi = dim[1]; yf = dim[3] | |
|
1042 | phase[xi:xf,yi:yf] = phase[xi:xf,yi:yf] + ues[ii] | |
|
1043 | ||
|
1044 | phase = -phase | |
|
1045 | ||
|
1046 | ar = self.eomwl*numpy.array([[0.5,6., 24.5],[0.5,6.,24.5]]) | |
|
1047 | nr = numpy.array([[12.,4.,2.],[12.,4.,2.]]) | |
|
1048 | lr = 0.25*self.eomwl*numpy.array([[0,0.,0],[0.,0,0]]) | |
|
1049 | ||
|
1050 | # Computing module and dipole patterns. | |
|
1051 | pattern = (numpy.abs(self.__dipPattern(ar,nr,lr)*self.__modPattern(phase,gain)))**2 | |
|
1052 | ||
|
1053 | return pattern | |
|
1054 | ||
|
1055 | def __usingFFT(self,rx): | |
|
1056 | """ | |
|
1057 | __usingFFT method returns the Jicamarca antenna pattern computed using The Fast Fou- | |
|
1058 | rier Transform. | |
|
1059 | ||
|
1060 | pattern = iFFT(FFT(gain*EXP(j*phase))) | |
|
1061 | ||
|
1062 | Parameters | |
|
1063 | ---------- | |
|
1064 | rx = Set to 1 to use the Rx information. Otherwise set to 0 for Tx. | |
|
1065 | ||
|
1066 | Return | |
|
1067 | ------ | |
|
1068 | pattern = An array giving the modelled antenna pattern using the array model. | |
|
1069 | ||
|
1070 | Modification history | |
|
1071 | -------------------- | |
|
1072 | Developed by Jorge L. Chau. | |
|
1073 | Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009. | |
|
1074 | """ | |
|
1075 | ||
|
1076 | if rx==1: | |
|
1077 | ues = self.uesrx | |
|
1078 | phase = self.phaserx | |
|
1079 | gain = self.gainrx | |
|
1080 | elif rx==0: | |
|
1081 | ues = self.uestx | |
|
1082 | phase = self.phasetx | |
|
1083 | gain = self.gaintx | |
|
1084 | ||
|
1085 | ues = ues*360./self.airwl | |
|
1086 | phase = phase*360./self.airwl | |
|
1087 | ||
|
1088 | for ii in range(4): | |
|
1089 | if ii==0:dim = numpy.array([4,0,8,4]) # WEST | |
|
1090 | elif ii==1:dim = numpy.array([0,0,4,4]) # NORTH | |
|
1091 | elif ii==2:dim = numpy.array([0,4,4,8]) # EAST | |
|
1092 | elif ii==3:dim = numpy.array([4,4,8,8]) # SOUTH | |
|
1093 | xi = dim[0]; xf = dim[2]; yi = dim[1]; yf = dim[3] | |
|
1094 | phase[xi:xf,yi:yf] = phase[xi:xf,yi:yf] + ues[ii] | |
|
1095 | ||
|
1096 | phase = -phase | |
|
1097 | ||
|
1098 | delta_x = self.eomwl/2. | |
|
1099 | delta_y = self.eomwl/2. | |
|
1100 | ||
|
1101 | nxfft = 2048 | |
|
1102 | nyfft = 2048 | |
|
1103 | dcosx = (numpy.arange(nxfft) - (0.5*nxfft))/(nxfft*delta_x)*self.eomwl | |
|
1104 | dcosy = (numpy.arange(nyfft) - (0.5*nyfft))/(nyfft*delta_y)*self.eomwl | |
|
1105 | ||
|
1106 | fft_gain = numpy.zeros((nxfft,nyfft)) | |
|
1107 | fft_phase = numpy.zeros((nxfft,nyfft)) | |
|
1108 | ||
|
1109 | nx = 8 | |
|
1110 | ny = 8 | |
|
1111 | ndx =12 | |
|
1112 | ndy =12 | |
|
1113 | for iy in numpy.arange(ny): | |
|
1114 | for ix in numpy.arange(nx): | |
|
1115 | ix1 = nxfft/2-self.nx/2*ndx+ix*ndx | |
|
1116 | if ix<(nx/2):ix1 = ix1 - 1 | |
|
1117 | if ix>=(nx/2):ix1 = ix1 + 1 | |
|
1118 | ||
|
1119 | iy1 = nyfft/2-ny/2*ndx+iy*ndy | |
|
1120 | if iy<(ny/2):iy1 = iy1 - 1 | |
|
1121 | if iy>=(ny/2):iy1 = iy1 + 1 | |
|
1122 | ||
|
1123 | fft_gain[ix1:ix1+ndx-1,iy1:iy1+ndy-1] = gain[ix,ny-1-iy] | |
|
1124 | fft_phase[ix1:ix1+ndx-1,iy1:iy1+ndy-1] = phase[ix,ny-1-iy] | |
|
1125 | ||
|
1126 | ||
|
1127 | fft_phase = fft_phase*Misc_Routines.CoFactors.d2r | |
|
1128 | ||
|
1129 | pattern = numpy.abs(numpy.fft.fft2(fft_gain*numpy.exp(numpy.complex(0,1)*fft_phase)))**2 | |
|
1130 | pattern = numpy.fft.fftshift(pattern) | |
|
1131 | ||
|
1132 | xvals = numpy.where((dcosx>=(numpy.min(self.dcosx))) & (dcosx<=(numpy.max(self.dcosx)))) | |
|
1133 | yvals = numpy.where((dcosy>=(numpy.min(self.dcosy))) & (dcosy<=(numpy.max(self.dcosy)))) | |
|
1134 | ||
|
1135 | pattern = pattern[xvals[0][0]:xvals[0][-1],yvals[0][0]:yvals[0][-1]] | |
|
1136 | ||
|
1137 | return pattern | |
|
1138 | ||
|
1139 | ||
|
1140 | def __dipPattern(self,ar,nr,lr): | |
|
1141 | """ | |
|
1142 | _dipPattern function computes the dipole's pattern to the Jicamarca radar. The next | |
|
1143 | equation defines the pattern as a function of the mainlobe direction: | |
|
1144 | ||
|
1145 | sincx = SIN(k/2*n0x*(a0x*SIN(phi)*COS(alpha)))/SIN(k/2*(a0x*SIN(phi)*COS(alpha))) | |
|
1146 | sincy = SIN(k/2*n0y*(a0y*SIN(phi)*SIN(alpha)))/SIN(k/2*(a0y*SIN(phi)*SIN(alpha))) | |
|
1147 | A0(phi,alpha) = sincx*sincy | |
|
1148 | Parameters | |
|
1149 | ---------- | |
|
1150 | ar = ? | |
|
1151 | nr = ? | |
|
1152 | lr = ? | |
|
1153 | ||
|
1154 | Return | |
|
1155 | ------ | |
|
1156 | dipole = An array giving antenna pattern from the dipole point of view.. | |
|
1157 | ||
|
1158 | Modification history | |
|
1159 | -------------------- | |
|
1160 | Developed by Jorge L. Chau. | |
|
1161 | Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009. | |
|
1162 | """ | |
|
1163 | ||
|
1164 | dipole = numpy.zeros((self.nx,self.ny),dtype=complex) | |
|
1165 | for iy in range(self.ny): | |
|
1166 | for ix in range(self.nx): | |
|
1167 | yindex = iy*(self.getcut==0) + ix*(self.getcut==1) | |
|
1168 | ||
|
1169 | argx = ar[0,0]*self.dcosx[ix] - lr[0,0] | |
|
1170 | if argx == 0.0: | |
|
1171 | junkx = nr[0,0] | |
|
1172 | else: | |
|
1173 | junkx = numpy.sin(0.5*self.kk*nr[0,0]*argx)/numpy.sin(0.5*self.kk*argx) | |
|
1174 | ||
|
1175 | ||
|
1176 | argy = ar[1,0]*self.dcosy[yindex] - lr[1,0] | |
|
1177 | if argy == 0.0: | |
|
1178 | junky = nr[1,0] | |
|
1179 | else: | |
|
1180 | junky = numpy.sin(0.5*self.kk*nr[1,0]*argy)/numpy.sin(0.5*self.kk*argy) | |
|
1181 | ||
|
1182 | ||
|
1183 | dipole[ix,iy] = junkx*junky | |
|
1184 | ||
|
1185 | return dipole | |
|
1186 | ||
|
1187 | def __modPattern(self,phase,gain): | |
|
1188 | """ | |
|
1189 | ModPattern computes the module's pattern to the Jicamarca radar. The next equation | |
|
1190 | defines the pattern as a function mainlobe direction: | |
|
1191 | ||
|
1192 | phasex = pos(x)*SIN(phi)*COS(alpha) | |
|
1193 | phasey = pos(y)*SIN(phi)*SIN(alpha) | |
|
1194 | ||
|
1195 | A1(phi,alpha) = TOTAL(gain*EXP(COMPLEX(0,k*(phasex+phasey)+phase))) | |
|
1196 | ||
|
1197 | Parameters | |
|
1198 | ---------- | |
|
1199 | phase = Bidimensional array (8x8) giving the phase (in meters) of each module. | |
|
1200 | gain = Bidimensional array (8x8) giving to define modules will be active (ones) | |
|
1201 | and which will not (zeros). | |
|
1202 | ||
|
1203 | Return | |
|
1204 | ------ | |
|
1205 | module = An array giving antenna pattern from the module point of view.. | |
|
1206 | ||
|
1207 | Modification history | |
|
1208 | -------------------- | |
|
1209 | Developed by Jorge L. Chau. | |
|
1210 | Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009. | |
|
1211 | """ | |
|
1212 | ||
|
1213 | pos = self.eomwl*attenuation | |
|
1214 | posx = pos[0,:,:] | |
|
1215 | posy = pos[1,:,:] | |
|
1216 | ||
|
1217 | phase = phase*Misc_Routines.CoFactors.d2r | |
|
1218 | module = numpy.zeros((self.nx,self.ny),dtype=complex) | |
|
1219 | for iy in range(self.ny): | |
|
1220 | for ix in range(self.nx): | |
|
1221 | yindex = iy*(self.getcut==0) + ix*(self.getcut==1) | |
|
1222 | phasex = posx*self.dcosx[ix] | |
|
1223 | phasey = posy*self.dcosy[yindex] | |
|
1224 | tmp = gain*numpy.exp(numpy.complex(0,1.)*(self.kk*(phasex+phasey)+phase)) | |
|
1225 | module[ix,iy] = tmp.sum() | |
|
1226 | ||
|
1227 | return module | |
|
1228 | ||
|
1229 | def __getBeamPars(self): | |
|
1230 | """ | |
|
1231 | _getBeamPars computes the main-beam parameters of the antenna. | |
|
1232 | ||
|
1233 | Modification history | |
|
1234 | -------------------- | |
|
1235 | Developed by Jorge L. Chau. | |
|
1236 | Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009. | |
|
1237 | """ | |
|
1238 | ||
|
1239 | dx = self.dcosx[1] - self.dcosx[0] | |
|
1240 | dy = self.dcosy[1] - self.dcosy[0] | |
|
1241 | ||
|
1242 | amp = self.norpattern | |
|
1243 | ||
|
1244 | xx = numpy.resize(self.dcosx,(self.nx,self.nx)).transpose() | |
|
1245 | yy = numpy.resize(self.dcosy,(self.ny,self.ny)) | |
|
1246 | ||
|
1247 | mm0 = amp[numpy.where(amp > 0.5)] | |
|
1248 | xx0 = xx[numpy.where(amp > 0.5)] | |
|
1249 | yy0 = yy[numpy.where(amp > 0.5)] | |
|
1250 | ||
|
1251 | xc = numpy.sum(mm0*xx0)/numpy.sum(mm0) | |
|
1252 | yc = numpy.sum(mm0*yy0)/numpy.sum(mm0) | |
|
1253 | rc = numpy.sqrt(mm0.size*dx*dy/numpy.pi) | |
|
1254 | ||
|
1255 | nnx = numpy.where(numpy.abs(self.dcosx - xc) < rc) | |
|
1256 | nny = numpy.where(numpy.abs(self.dcosy - yc) < rc) | |
|
1257 | ||
|
1258 | mm1 = amp[numpy.min(nnx):numpy.max(nnx)+1,numpy.min(nny):numpy.max(nny)+1] | |
|
1259 | xx1 = self.dcosx[numpy.min(nnx):numpy.max(nnx)+1] | |
|
1260 | yy1 = self.dcosy[numpy.min(nny):numpy.max(nny)+1] | |
|
1261 | ||
|
1262 | # fitting data into the main beam. | |
|
1263 | ||
|
1264 | params = gaussfit.fitgaussian(mm1) | |
|
1265 | ||
|
1266 | # Tranforming from indexes to axis' values | |
|
1267 | xcenter = xx1[0] + (((xx1[xx1.size-1] - xx1[0])/(xx1.size -1))*(params[1])) | |
|
1268 | ycenter = yy1[0] + (((yy1[yy1.size-1] - yy1[0])/(yy1.size -1))*(params[2])) | |
|
1269 | xwidth = ((xx1[xx1.size-1] - xx1[0])/(xx1.size-1))*(params[3])*(1/Misc_Routines.CoFactors.d2r) | |
|
1270 | ywidth = ((yy1[yy1.size-1] - yy1[0])/(yy1.size-1))*(params[4])*(1/Misc_Routines.CoFactors.d2r) | |
|
1271 | meanwx = (xwidth*ywidth) | |
|
1272 | meanpos = numpy.array([xcenter,ycenter]) | |
|
1273 | ||
|
1274 | #print 'Position: %f %f' %(xcenter,ycenter) | |
|
1275 | #print 'Widths: %f %f' %(xwidth, ywidth) | |
|
1276 | #print 'BWHP: %f' %(2*numpy.sqrt(2*meanwx)*numpy.sqrt(-numpy.log(0.5))) | |
|
1277 | ||
|
1278 | self.meanpos = meanpos | |
|
1279 | ||
|
1280 | ||
|
1281 | class overJroShow: | |
|
1282 | ||
|
1283 | __serverdocspath = '' | |
|
1284 | __tmpDir = '' | |
|
1285 | ||
|
1286 | def __init__(self, title='', heights=None): | |
|
1287 | self.year = None | |
|
1288 | self.month = None | |
|
1289 | self.dom = None | |
|
1290 | self.pattern = None | |
|
1291 | self.maxphi = None | |
|
1292 | self.heights = None | |
|
1293 | self.filename = None | |
|
1294 | self.showType = None | |
|
1295 | self.path = None | |
|
1296 | self.objects = None | |
|
1297 | self.nptsx = 101 | |
|
1298 | self.nptsy = 101 | |
|
1299 | self.fftopt = 0 | |
|
1300 | self.site = 1 | |
|
1301 | self.dcosx = 1 | |
|
1302 | self.dcosy = 1 | |
|
1303 | self.dcosxrange = None | |
|
1304 | self.dcosyrange = None | |
|
1305 | self.maxha_min= 0. | |
|
1306 | self.show_object = None | |
|
1307 | self.dcosx_mag = None | |
|
1308 | self.dcosy_mag = None | |
|
1309 | self.ha_mag = None | |
|
1310 | self.time_mag = None | |
|
1311 | self.main_dec = None | |
|
1312 | self.ObjC = None | |
|
1313 | self.ptitle = title | |
|
1314 | self.path4plotname = None | |
|
1315 | self.plotname0 = None | |
|
1316 | self.plotname1 = None | |
|
1317 | self.plotname2 = None | |
|
1318 | self.scriptHeaders = 0 | |
|
1319 | self.glat = -11.95 | |
|
1320 | self.glon = -76.8667 | |
|
1321 | self.UT = 5 #timezone | |
|
1322 | ||
|
1323 | self.glat = -11.951481 | |
|
1324 | self.glon = -76.874383 | |
|
1325 | if heights is None: | |
|
1326 | self.heights = numpy.array([100.,500.,1000.]) | |
|
1327 | else: | |
|
1328 | self.heights = numpy.array(heights) | |
|
1329 | ||
|
1330 | def initParameters1(self): | |
|
1331 | ||
|
1332 | gui=1 | |
|
1333 | if self.pattern==None: | |
|
1334 | if gui==1: self.filename = self.filename.split(',') | |
|
1335 | ||
|
1336 | pattern = numpy.atleast_1d(self.pattern) | |
|
1337 | filename = numpy.atleast_1d(self.filename) | |
|
1338 | ||
|
1339 | npatterns = numpy.max(numpy.array([pattern.size,filename.size])) | |
|
1340 | ||
|
1341 | self.pattern = numpy.resize(pattern,npatterns) | |
|
1342 | self.filename = numpy.resize(filename,npatterns) | |
|
1343 | ||
|
1344 | self.doy = datetime.datetime(self.year,self.month,self.dom).timetuple().tm_yday | |
|
1345 | ||
|
1346 | ||
|
1347 | if self.objects==None: | |
|
1348 | self.objects=numpy.zeros(5) | |
|
1349 | else: | |
|
1350 | tmp = numpy.atleast_1d(self.objects) | |
|
1351 | self.objects = numpy.zeros(5) | |
|
1352 | self.objects[0:tmp.size] = tmp | |
|
1353 | ||
|
1354 | self.show_object = self.objects | |
|
1355 | ||
|
1356 | self.maxha_min = 4*self.maxphi*numpy.sqrt(2)*1.25 | |
|
1357 | ||
|
1358 | ||
|
1359 | ||
|
1360 | #ROJ geographic coordinates and time zone | |
|
1361 | self.UT = 5 #timezone | |
|
1362 | self.glat = -11.951481 | |
|
1363 | self.glon = -76.874383 | |
|
1364 | ||
|
1365 | ||
|
1366 | self.junkjd = TimeTools.Time(self.year,self.month,self.dom).change2julday() | |
|
1367 | self.junklst = TimeTools.Julian(self.junkjd).change2lst(longitude=self.glon) | |
|
1368 | ||
|
1369 | # Finding RA of observatory for a specific date | |
|
1370 | self.ra_obs = self.junklst*Misc_Routines.CoFactors.h2d | |
|
1371 | ||
|
1372 | def initParameters(self, dt): | |
|
1373 | ||
|
1374 | self.year = dt.year | |
|
1375 | self.month = dt.month | |
|
1376 | self.dom = dt.day | |
|
1377 | # Defining plot filenames | |
|
1378 | self.path4plotname = os.path.join(self.__serverdocspath,self.__tmpDir) | |
|
1379 | self.plotname0 = 'over_jro_0_%i.png'% (time.time()) #plot pattern & objects | |
|
1380 | self.plotname1 = 'over_jro_1_%i.png'% (time.time()) #plot antenna cuts | |
|
1381 | self.plotname2 = 'over_jro_2_%i.png'% (time.time()) #plot sky noise | |
|
1382 | ||
|
1383 | # Defining antenna axes respect to geographic coordinates (See Ochs report). | |
|
1384 | # alfa = 1.46*Misc_Routines.CoFactors.d2r | |
|
1385 | # theta = 51.01*Misc_Routines.CoFactors.d2r | |
|
1386 | ||
|
1387 | alfa = 1.488312*Misc_Routines.CoFactors.d2r | |
|
1388 | th = 6.166710 + 45.0 | |
|
1389 | theta = th*Misc_Routines.CoFactors.d2r | |
|
1390 | ||
|
1391 | self.maxha_min = 4*self.maxphi*numpy.sqrt(2)*1.25 | |
|
1392 | ||
|
1393 | sina = numpy.sin(alfa) | |
|
1394 | cosa = numpy.cos(alfa) | |
|
1395 | MT1 = numpy.array([[1,0,0],[0,cosa,-sina],[0,sina,cosa]]) | |
|
1396 | sinb = numpy.sin(theta) | |
|
1397 | cosb = numpy.cos(theta) | |
|
1398 | MT2 = numpy.array([[cosb,sinb,0],[-sinb,cosb,0],[0,0,1]]) | |
|
1399 | self.MT3 = numpy.array(numpy.dot(MT2, MT1)).transpose() | |
|
1400 | ||
|
1401 | self.xg = numpy.dot(self.MT3.transpose(),numpy.array([1,0,0])) | |
|
1402 | self.yg = numpy.dot(self.MT3.transpose(),numpy.array([0,1,0])) | |
|
1403 | self.zg = numpy.dot(self.MT3.transpose(),numpy.array([0,0,1])) | |
|
1404 | ||
|
1405 | def plotPattern2(self, date, phases, gain_tx, gain_rx, ues, just_rx, angle): | |
|
1406 | # Plotting Antenna patterns. | |
|
1407 | self.maxphi = angle | |
|
1408 | self.initParameters(date) | |
|
1409 | self.doy = datetime.datetime(date.year,date.month,date.day).timetuple().tm_yday | |
|
1410 | self.junkjd = TimeTools.Time(self.year,self.month,self.dom).change2julday() | |
|
1411 | self.junklst = TimeTools.Julian(self.junkjd).change2lst(longitude=self.glon) | |
|
1412 | self.ra_obs = self.junklst*Misc_Routines.CoFactors.h2d | |
|
1413 | ||
|
1414 | date = TimeTools.Time(date.year,date.month,date.day).change2strdate(mode=2) | |
|
1415 | ||
|
1416 | mesg = 'Over Jicamarca: ' + date[0] | |
|
1417 | ||
|
1418 | ObjAnt = JroPattern(pattern=0, | |
|
1419 | filename=None, | |
|
1420 | path=None, | |
|
1421 | nptsx=self.nptsx, | |
|
1422 | nptsy=self.nptsy, | |
|
1423 | maxphi=angle, | |
|
1424 | fftopt=self.fftopt, | |
|
1425 | phases=phases, | |
|
1426 | gain_tx=gain_tx, | |
|
1427 | gain_rx=gain_rx, | |
|
1428 | ues=ues, | |
|
1429 | just_rx=just_rx | |
|
1430 | ) | |
|
1431 | ||
|
1432 | self.pattern_plot = AntPatternPlot() | |
|
1433 | ||
|
1434 | self.pattern_plot.contPattern(iplot=0, | |
|
1435 | gpath=self.path4plotname, | |
|
1436 | filename=self.plotname0, | |
|
1437 | mesg=mesg, | |
|
1438 | amp=ObjAnt.norpattern, | |
|
1439 | x=ObjAnt.dcosx, | |
|
1440 | y=ObjAnt.dcosy, | |
|
1441 | getCut=ObjAnt.getcut, | |
|
1442 | title=self.ptitle, | |
|
1443 | save=False) | |
|
1444 | ||
|
1445 | ||
|
1446 | self.pattern_plot.plotRaDec(gpath=self.path4plotname, | |
|
1447 | filename=self.plotname0, | |
|
1448 | jd=self.junkjd, | |
|
1449 | ra_obs=self.ra_obs, | |
|
1450 | xg=self.xg, | |
|
1451 | yg=self.yg, | |
|
1452 | x=ObjAnt.dcosx, | |
|
1453 | y=ObjAnt.dcosy, | |
|
1454 | save=False) | |
|
1455 | ||
|
1456 | vect_ant = numpy.array([ObjAnt.meanpos[0],ObjAnt.meanpos[1],numpy.sqrt(1-numpy.sum(ObjAnt.meanpos**2.))]) | |
|
1457 | vect_geo = numpy.dot(scipy.linalg.inv(self.MT3),vect_ant) | |
|
1458 | vect_polar = Misc_Routines.Vector(numpy.array(vect_geo),direction=1).Polar2Rect() | |
|
1459 | [ra,dec,ha] = Astro_Coords.AltAz(vect_polar[1],vect_polar[0],self.junkjd).change2equatorial() | |
|
1460 | self.main_dec = dec | |
|
1461 | ||
|
1462 | def plotBfield(self): | |
|
1463 | ||
|
1464 | ObjB = BField(self.year,self.doy,1,self.heights) | |
|
1465 | [dcos, alpha, nlon, nlat] = ObjB.getBField() | |
|
1466 | ||
|
1467 | self.pattern_plot.plotBField('', '',dcos,alpha,nlon,nlat, | |
|
1468 | self.dcosxrange, | |
|
1469 | self.dcosyrange, | |
|
1470 | ObjB.heights, | |
|
1471 | ObjB.alpha_i, | |
|
1472 | save=False) | |
|
1473 | ||
|
1474 | def plotCelestial(self, objects): | |
|
1475 | ||
|
1476 | ntod = 24.*16. | |
|
1477 | ||
|
1478 | tod = numpy.arange(ntod)/ntod*24. | |
|
1479 | ||
|
1480 | [month,dom] = TimeTools.Doy2Date(self.year, self.doy).change2date() | |
|
1481 | ||
|
1482 | jd = TimeTools.Time(self.year, month, dom, tod+self.UT).change2julday() | |
|
1483 | ||
|
1484 | self.pattern_plot.plotCelestial( | |
|
1485 | jd, | |
|
1486 | self.main_dec, | |
|
1487 | tod, | |
|
1488 | self.maxha_min, | |
|
1489 | objects, | |
|
1490 | self.glat, | |
|
1491 | self.glon, | |
|
1492 | self.xg, | |
|
1493 | self.yg, | |
|
1494 | self.dcosxrange, | |
|
1495 | self.dcosyrange | |
|
1496 | ) | |
|
1497 | ||
|
1498 | def plotAntennaCuts(self): | |
|
1499 | # print "Drawing antenna cuts" | |
|
1500 | ||
|
1501 | incha = 0.05 # min | |
|
1502 | nha = numpy.int32(2*self.maxha_min/incha) + 1. | |
|
1503 | newha = numpy.arange(nha)/nha*2.*self.maxha_min - self.maxha_min | |
|
1504 | nha_star = numpy.int32(200./incha) | |
|
1505 | star_ha = (numpy.arange(nha_star) - (nha_star/2))*nha_star | |
|
1506 | ||
|
1507 | #Init ObjCut for PatternCutPlot() | |
|
1508 | view_objects = numpy.where(self.show_object>0) | |
|
1509 | subplots = len(view_objects[0]) | |
|
1510 | ObjCut = Graphics_OverJro.PatternCutPlot(subplots) | |
|
1511 | ||
|
1512 | for io in (numpy.arange(5)): | |
|
1513 | if self.show_object[io]==2: | |
|
1514 | if io==0: | |
|
1515 | if self.dcosx_mag.size!=0: | |
|
1516 | dcosx = self.dcosx_mag | |
|
1517 | dcosy = self.dcosy_mag | |
|
1518 | dcosz = 1 - numpy.sqrt(dcosx**2. + dcosy**2.) | |
|
1519 | ||
|
1520 | # Finding rotation of B respec to antenna coords. | |
|
1521 | [mm,bb] = scipy.polyfit(dcosx,dcosy,1) | |
|
1522 | alfa = 0.0 | |
|
1523 | theta = -1.*numpy.arctan(mm) | |
|
1524 | sina = numpy.sin(alfa); cosa = numpy.cos(alfa) | |
|
1525 | MT1 = [[1,0,0],[0,cosa,-sina],[0,sina,cosa]] | |
|
1526 | MT1 = numpy.array(MT1) | |
|
1527 | sinb = numpy.sin(theta); cosb = numpy.cos(theta) | |
|
1528 | MT2 = [[cosb,sinb,0],[-sinb,cosb,0],[0,0,1]] | |
|
1529 | MT2 = numpy.array(MT2) | |
|
1530 | MT3_mag = numpy.dot(MT2, MT1) | |
|
1531 | MT3_mag = numpy.array(MT3_mag).transpose() | |
|
1532 | # Getting dcos respec to B coords | |
|
1533 | vector = numpy.array([dcosx,dcosy,dcosz]) | |
|
1534 | nvector = numpy.dot(MT3_mag,vector) | |
|
1535 | nvector = numpy.array(nvector).transpose() | |
|
1536 | ||
|
1537 | ## print 'Rotation (deg) %f'%(theta/Misc_Routines.CoFactors.d2r) | |
|
1538 | ||
|
1539 | yoffset = numpy.sum(nvector[:,1])/nvector[:,1].size | |
|
1540 | # print 'Dcosyoffset %f'%(yoffset) | |
|
1541 | ||
|
1542 | ha = self.ha_mag*4. | |
|
1543 | time = self.time_mag | |
|
1544 | width_star = 0.1 # half width in minutes | |
|
1545 | otitle = 'B Perp. cut' | |
|
1546 | # else: | |
|
1547 | # print "No B perp. over Observatory" | |
|
1548 | # | |
|
1549 | # | |
|
1550 | elif io==1: | |
|
1551 | if self.ObjC.dcosx_sun.size!=0: | |
|
1552 | dcosx = self.ObjC.dcosx_sun | |
|
1553 | dcosy = self.ObjC.dcosy_sun | |
|
1554 | ha = self.ObjC.ha_sun*4.0 | |
|
1555 | time = self.ObjC.time_sun | |
|
1556 | width_star = 2. # half width in minutes | |
|
1557 | otitle = 'Sun cut' | |
|
1558 | # else: | |
|
1559 | # print "Sun is not passing over Observatory" | |
|
1560 | ||
|
1561 | elif io==2: | |
|
1562 | if self.ObjC.dcosx_moon.size!=0: | |
|
1563 | dcosx = self.ObjC.dcosx_moon | |
|
1564 | dcosy = self.ObjC.dcosy_moon | |
|
1565 | ha = self.ObjC.ha_moon*4 | |
|
1566 | time = self.ObjC.time_moon | |
|
1567 | m_distance = 404114.6 # distance to the Earth in km | |
|
1568 | m_diameter = 1734.4 # diameter in km. | |
|
1569 | width_star = numpy.arctan(m_distance/m_diameter) | |
|
1570 | width_star = width_star/2./Misc_Routines.CoFactors.d2r*4. | |
|
1571 | otitle = 'Moon cut' | |
|
1572 | # else: | |
|
1573 | # print "Moon is not passing over Observatory" | |
|
1574 | ||
|
1575 | elif io==3: | |
|
1576 | if self.ObjC.dcosx_hydra.size!=0: | |
|
1577 | dcosx = self.ObjC.dcosx_hydra | |
|
1578 | dcosy = self.ObjC.dcosy_hydra | |
|
1579 | ha = self.ObjC.ha_hydra*4. | |
|
1580 | time = self.ObjC.time_hydra | |
|
1581 | width_star = 0.25 # half width in minutes | |
|
1582 | otitle = 'Hydra cut' | |
|
1583 | # else: | |
|
1584 | # print "Hydra is not passing over Observatory" | |
|
1585 | ||
|
1586 | elif io==4: | |
|
1587 | if self.ObjC.dcosx_galaxy.size!=0: | |
|
1588 | dcosx = self.ObjC.dcosx_galaxy | |
|
1589 | dcosy = self.ObjC.dcosy_galaxy | |
|
1590 | ha = self.ObjC.ha_galaxy*4. | |
|
1591 | time = self.ObjC.time_galaxy | |
|
1592 | width_star = 25. # half width in minutes | |
|
1593 | otitle = 'Galaxy cut' | |
|
1594 | # else: | |
|
1595 | # print "Galaxy center is not passing over Jicamarca" | |
|
1596 | # | |
|
1597 | # | |
|
1598 | hour = numpy.int32(time) | |
|
1599 | mins = numpy.int32((time - hour)*60.) | |
|
1600 | secs = numpy.int32(((time - hour)*60. - mins)*60.) | |
|
1601 | ||
|
1602 | ObjT = TimeTools.Time(self.year,self.month,self.dom,hour,mins,secs) | |
|
1603 | subtitle = ObjT.change2strdate() | |
|
1604 | ||
|
1605 | star_cut = numpy.exp(-(star_ha/width_star)**2./2.) | |
|
1606 | ||
|
1607 | pol = scipy.polyfit(ha,dcosx,3.) | |
|
1608 | polx = numpy.poly1d(pol); newdcosx = polx(newha) | |
|
1609 | pol = scipy.polyfit(ha,dcosy,3.) | |
|
1610 | poly = numpy.poly1d(pol);newdcosy = poly(newha) | |
|
1611 | ||
|
1612 | patterns = [] | |
|
1613 | for icut in numpy.arange(self.pattern.size): | |
|
1614 | # Getting Antenna cut. | |
|
1615 | Obj = JroPattern(dcosx=newdcosx, | |
|
1616 | dcosy=newdcosy, | |
|
1617 | getcut=1, | |
|
1618 | pattern=self.pattern[icut], | |
|
1619 | path=self.path, | |
|
1620 | filename=self.filename[icut]) | |
|
1621 | ||
|
1622 | Obj.getPattern() | |
|
1623 | ||
|
1624 | patterns.append(Obj.pattern) | |
|
1625 | ||
|
1626 | ||
|
1627 | ObjCut.drawCut(io, | |
|
1628 | patterns, | |
|
1629 | self.pattern.size, | |
|
1630 | newha, | |
|
1631 | otitle, | |
|
1632 | subtitle, | |
|
1633 | self.ptitle) | |
|
1634 | ||
|
1635 | ObjCut.saveFig(self.path4plotname,self.plotname1) | |
|
1636 | ||
|
1637 | ||
|
12 | 1638 | |
|
13 | 1639 | def skyNoise(jd, ut=-5.0, longitude=-76.87, filename='/app/utils/galaxy.txt'): |
|
14 | 1640 | """ |
@@ -186,3 +1812,38 def skynoise_plot(year, month, day): | |||
|
186 | 1812 | |
|
187 | 1813 | return buf |
|
188 | 1814 | |
|
1815 | def overjro_plot(pattern_id, date, angle, height, bodys): | |
|
1816 | ||
|
1817 | pattern = select_pattern(pattern = pattern_id) | |
|
1818 | ob = overJroShow(pattern['title'], heights=height) | |
|
1819 | ob.plotPattern2( | |
|
1820 | date, | |
|
1821 | pattern['phase'], | |
|
1822 | pattern['gaintx'], | |
|
1823 | pattern['gainrx'], | |
|
1824 | pattern['ues'], | |
|
1825 | pattern['justrx'], | |
|
1826 | angle | |
|
1827 | ) | |
|
1828 | ||
|
1829 | if 'bfield' in bodys: | |
|
1830 | ob.plotBfield() | |
|
1831 | ||
|
1832 | objects = [] | |
|
1833 | for body in bodys: | |
|
1834 | if body=='sun': | |
|
1835 | objects.append(0) | |
|
1836 | elif body=='moon': | |
|
1837 | objects.append(1) | |
|
1838 | elif body=='hydra': | |
|
1839 | objects.append(2) | |
|
1840 | elif body=='galaxy': | |
|
1841 | objects.append(3) | |
|
1842 | ||
|
1843 | if objects: | |
|
1844 | ob.plotCelestial(objects) | |
|
1845 | ||
|
1846 | buf = BytesIO() | |
|
1847 | ob.pattern_plot.fig.savefig(buf, format="png") | |
|
1848 | ||
|
1849 | return buf No newline at end of file |
General Comments 0
You need to be logged in to leave comments.
Login now