##// END OF EJS Templates
Add over JRO tool
jespinoza -
r49:78d40e1a09aa
parent child
Show More

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 name="experiment" class="form-control form-control-sm">
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="20">Vertical Drifts</option>
53 <option value="[21,22]">East West 1996</option>
54 <option value="[25,26]">East West 2003</option>
55 <option value="23">Differential Phase 2000</option>
56 <option value="24">Differential Phase 2004 High Alt</option>
57 <option value="27">Differential Phase 2005 - 2006</option>
58 <option value="[28,29]">DEWD 2005</option>
59 <option value="2710">DVD 2006 - 2008</option>
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="10">Oblique ISR On-Axis</option>
62 <option value="11">Oblique ISR 4.5</option>
63 <option value="12">Oblique ISR 6.0S</option>
64 <option value="13">Oblique ISR 3.0N</option>
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="[30,31]">JULIA CP2</option>
67 <option value="32">JULIA CP3</option>
68 <option value="35">JULIA V (2005-2006)</option>
69 <option value="[33,34]">JULIA EW 2003</option>
70 <option value="[35,36]">JULIA EW (2006-2007)</option>
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="40">EW Imaging 1996</option>
78 <option value="41">EW Imaging 2003</option>
79 <option value="43">EW Imaging 2006-2008</option>
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="50">MST North (Fritts)</option>
82 <option value="51">MST West (Fritts)</option>
83 <option value="52">MST South (Fritts)</option>
84 <option value="53">MST East (Fritts)</option>
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="54">Vertical (Yellow Cables)</option>
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="option1">
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="option1">
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="option1">
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="option1">
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="option1">
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