##// END OF EJS Templates
schain
RSS
changes from v4
joabAM -
r1664:a9d92c6e95c5
parent child
Show More
Side by Side Unified
Context file:
r1664:a9d92c6e95c5 Loading diff...:
Show file before | Show file after | Hide comments
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 schainpy.model.utils import TimeTools
20 from schainpy.model.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
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (507 lines changed) Show them Hide them
@@ -0,0 +1,507
1 """
2
3 Converted to Object-oriented Programming by Freddy Galindo, ROJ, 07 October 2009.
4 Added to signal Chain by Joab Apaza, ROJ, Jun 2023.
5 """
6
7 import numpy
8 import time
9 import os
10 from scipy.special import lpmn
11 from schainpy.model.utils import Astro_Coords
12
13 class BField():
14 def __init__(self,year=None,doy=None,site=1,heights=None,alpha_i=90):
15 """
16 BField class creates an object to get the Magnetic field for a specific date and
17 height(s).
18
19 Parameters
20 ----------
21 year = A scalar giving the desired year. If the value is None (default value) then
22 the current year will be used.
23 doy = A scalar giving the desired day of the year. If the value is None (default va-
24 lue) then the current doy will be used.
25 site = An integer to choose the geographic coordinates of the place where the magne-
26 tic field will be computed. The default value is over Jicamarca (site=1)
27 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.
28 alpha_i = Angle to interpolate the magnetic field.
29
30 Modification History
31 --------------------
32 Converted to Object-oriented Programming by Freddy Galindo, ROJ, 07 October 2009.
33 Added to signal Chain by Joab Apaza, ROJ, Jun 2023.
34 """
35
36 tmp = time.localtime()
37 if year==None: year = tmp[0]
38 if doy==None: doy = tmp[7]
39 self.year = year
40 self.doy = doy
41 self.site = site
42 if heights is None:
43 heights = numpy.array([100,500,1000])
44 else:
45 heights = numpy.array(heights)
46 self.heights = heights
47 self.alpha_i = alpha_i
48
49 def getBField(self,maglimits=numpy.array([-37,-37,37,37])):
50 """
51 getBField models the magnetic field for a different heights in a specific date.
52
53 Parameters
54 ----------
55 maglimits = An 4-elements array giving ..... The default value is [-7,-7,7,7].
56
57 Return
58 ------
59 dcos = An 4-dimensional array giving the directional cosines of the magnetic field
60 over the desired place.
61 alpha = An 3-dimensional array giving the angle of the magnetic field over the desi-
62 red place.
63
64 Modification History
65 --------------------
66 Converted to Python by Freddy R. Galindo, ROJ, 07 October 2009.
67 """
68
69 x_ant = numpy.array([1,0,0])
70 y_ant = numpy.array([0,1,0])
71 z_ant = numpy.array([0,0,1])
72
73 if self.site==0:
74 title_site = "Magnetic equator"
75 coord_site = numpy.array([-76+52./60.,-11+57/60.,0.5])
76 elif self.site==1:
77 title_site = 'Jicamarca'
78 coord_site = [-76-52./60.,-11-57/60.,0.5]
79 heta = (45+5.35)*numpy.pi/180. # (50.35 and 1.46 from Fleish Thesis)
80 delta = -1.46*numpy.pi/180
81
82
83 x_ant1 = numpy.roll(self.rotvector(self.rotvector(x_ant,1,delta),3,theta),1)
84 y_ant1 = numpy.roll(self.rotvector(self.rotvector(y_ant,1,delta),3,theta),1)
85 z_ant1 = numpy.roll(self.rotvector(self.rotvector(z_ant,1,delta),3,theta),1)
86
87 ang0 = -1*coord_site[0]*numpy.pi/180.
88 ang1 = coord_site[1]*numpy.pi/180.
89 x_ant = self.rotvector(self.rotvector(x_ant1,2,ang1),3,ang0)
90 y_ant = self.rotvector(self.rotvector(y_ant1,2,ang1),3,ang0)
91 z_ant = self.rotvector(self.rotvector(z_ant1,2,ang1),3,ang0)
92
93 elif self.site==2: #AMISR
94 title_site = 'AMISR 14'
95
96 coord_site = [-76.874913, -11.953371, 0.52984]
97
98 theta = (0.0977)*numpy.pi/180. # 0.0977
99 delta = 0.110*numpy.pi/180 # 0.11
100
101 x_ant1 = numpy.roll(self.rotvector(self.rotvector(x_ant,1,delta),3,theta),1)
102 y_ant1 = numpy.roll(self.rotvector(self.rotvector(y_ant,1,delta),3,theta),1)
103 z_ant1 = numpy.roll(self.rotvector(self.rotvector(z_ant,1,delta),3,theta),1)
104
105 ang0 = -1*coord_site[0]*numpy.pi/180.
106 ang1 = coord_site[1]*numpy.pi/180.
107 x_ant = self.rotvector(self.rotvector(x_ant1,2,ang1),3,ang0)
108 y_ant = self.rotvector(self.rotvector(y_ant1,2,ang1),3,ang0)
109 z_ant = self.rotvector(self.rotvector(z_ant1,2,ang1),3,ang0)
110 else:
111 # print "No defined Site. Skip..."
112 return None
113
114 nhei = self.heights.size
115 pt_intercep = numpy.zeros((nhei,2))
116 nfields = 1
117
118 grid_res = 2.5
119 nlon = int(int(maglimits[2] - maglimits[0])/grid_res + 1)
120 nlat = int(int(maglimits[3] - maglimits[1])/grid_res + 1)
121
122 location = numpy.zeros((nlon,nlat,2))
123 mlon = numpy.atleast_2d(numpy.arange(nlon)*grid_res + maglimits[0])
124 mrep = numpy.atleast_2d(numpy.zeros(nlat) + 1)
125 location0 = numpy.dot(mlon.transpose(),mrep)
126
127 mlat = numpy.atleast_2d(numpy.arange(nlat)*grid_res + maglimits[1])
128 mrep = numpy.atleast_2d(numpy.zeros(nlon) + 1)
129 location1 = numpy.dot(mrep.transpose(),mlat)
130
131 location[:,:,0] = location0
132 location[:,:,1] = location1
133
134 alpha = numpy.zeros((nlon,nlat,nhei))
135 rr = numpy.zeros((nlon,nlat,nhei,3))
136 dcos = numpy.zeros((nlon,nlat,nhei,2))
137
138 global first_time
139
140 first_time = None
141 for ilon in numpy.arange(nlon):
142 for ilat in numpy.arange(nlat):
143 outs = self.__bdotk(self.heights,
144 self.year + self.doy/366.,
145 coord_site[1],
146 coord_site[0],
147 coord_site[2],
148 coord_site[1]+location[ilon,ilat,1],
149 location[ilon,ilat,0]*720./180.)
150
151 alpha[ilon, ilat,:] = outs[1]
152 rr[ilon, ilat,:,:] = outs[3]
153
154 mrep = numpy.atleast_2d((numpy.zeros(nhei)+1)).transpose()
155 tmp = outs[3]*numpy.dot(mrep,numpy.atleast_2d(x_ant))
156 tmp = tmp.sum(axis=1)
157 dcos[ilon,ilat,:,0] = tmp/numpy.sqrt((outs[3]**2).sum(axis=1))
158
159 mrep = numpy.atleast_2d((numpy.zeros(nhei)+1)).transpose()
160 tmp = outs[3]*numpy.dot(mrep,numpy.atleast_2d(y_ant))
161 tmp = tmp.sum(axis=1)
162 dcos[ilon,ilat,:,1] = tmp/numpy.sqrt((outs[3]**2).sum(axis=1))
163
164 return dcos, alpha, nlon, nlat
165
166
167 def __bdotk(self,heights,tm,gdlat=-11.95,gdlon=-76.8667,gdalt=0.0,decd=-12.88, ham=-4.61666667):
168
169 global first_time
170 # Mean Earth radius in Km WGS 84
171 a_igrf = 6371.2
172
173 bk = numpy.zeros(heights.size)
174 alpha = numpy.zeros(heights.size)
175 bfm = numpy.zeros(heights.size)
176 rr = numpy.zeros((heights.size,3))
177 rgc = numpy.zeros((heights.size,3))
178
179 ObjGeodetic = Astro_Coords.Geodetic(gdlat,gdalt)
180 [gclat,gcalt] = ObjGeodetic.change2geocentric()
181
182 gclat = gclat*numpy.pi/180.
183 gclon = gdlon*numpy.pi/180.
184
185 # Antenna position from center of Earth
186 ca_vector = [numpy.cos(gclat)*numpy.cos(gclon),numpy.cos(gclat)*numpy.sin(gclon),numpy.sin(gclat)]
187 ca_vector = gcalt*numpy.array(ca_vector)
188
189 dec = decd*numpy.pi/180.
190
191 # K vector respect to the center of earth.
192 klon = gclon + ham*numpy.pi/720.
193 k_vector = [numpy.cos(dec)*numpy.cos(klon),numpy.cos(dec)*numpy.sin(klon),numpy.sin(dec)]
194 k_vector = numpy.array(k_vector)
195
196 for ih in numpy.arange(heights.size):
197 # Vector from Earth's center to volume of interest
198 rr[ih,:] = k_vector*heights[ih]
199 cv_vector = numpy.squeeze(ca_vector) + rr[ih,:]
200
201 cv_gcalt = numpy.sqrt(numpy.sum(cv_vector**2.))
202 cvxy = numpy.sqrt(numpy.sum(cv_vector[0:2]**2.))
203
204 radial = cv_vector/cv_gcalt
205 east = numpy.array([-1*cv_vector[1],cv_vector[0],0])/cvxy
206 comp1 = east[1]*radial[2] - radial[1]*east[2]
207 comp2 = east[2]*radial[0] - radial[2]*east[0]
208 comp3 = east[0]*radial[1] - radial[0]*east[1]
209 north = -1*numpy.array([comp1, comp2, comp3])
210
211 rr_k = cv_vector - numpy.squeeze(ca_vector)
212 u_rr = rr_k/numpy.sqrt(numpy.sum(rr_k**2.))
213
214 cv_gclat = numpy.arctan2(cv_vector[2],cvxy)
215 cv_gclon = numpy.arctan2(cv_vector[1],cv_vector[0])
216
217 bhei = cv_gcalt-a_igrf
218 blat = cv_gclat*180./numpy.pi
219 blon = cv_gclon*180./numpy.pi
220 bfield = self.__igrfkudeki(bhei,tm,blat,blon)
221
222 B = (bfield[0]*north + bfield[1]*east - bfield[2]*radial)*1.0e-5
223
224 bfm[ih] = numpy.sqrt(numpy.sum(B**2.)) #module
225 bk[ih] = numpy.sum(u_rr*B)
226 alpha[ih] = numpy.arccos(bk[ih]/bfm[ih])*180/numpy.pi
227 rgc[ih,:] = numpy.array([cv_gclon, cv_gclat, cv_gcalt])
228
229 return bk, alpha, bfm, rr, rgc
230
231
232 def __igrfkudeki(self,heights,time,latitude,longitude,ae=6371.2):
233 """
234 __igrfkudeki calculates the International Geomagnetic Reference Field for given in-
235 put conditions based on IGRF2005 coefficients.
236
237 Parameters
238 ----------
239 heights = Scalar or vector giving the height above the Earth of the point in ques-
240 tion in kilometers.
241 time = Scalar or vector giving the decimal year of time in question (e.g. 1991.2).
242 latitude = Latitude of point in question in decimal degrees. Scalar or vector.
243 longitude = Longitude of point in question in decimal degrees. Scalar or vector.
244 ae =
245 first_time =
246
247 Return
248 ------
249 bn =
250 be =
251 bd =
252 bmod =
253 balpha =
254 first_time =
255
256 Modification History
257 --------------------
258 Converted to Python by Freddy R. Galindo, ROJ, 03 October 2009.
259 """
260
261 global first_time
262 global gs, hs, nvec, mvec, maxcoef
263
264 heights = numpy.atleast_1d(heights)
265 time = numpy.atleast_1d(time)
266 latitude = numpy.atleast_1d(latitude)
267 longitude = numpy.atleast_1d(longitude)
268
269 if numpy.max(latitude)==90:
270 # print "Field calculations are not supported at geographic poles"
271 pass
272
273 # output arrays
274 bn = numpy.zeros(heights.size)
275 be = numpy.zeros(heights.size)
276 bd = numpy.zeros(heights.size)
277
278 if first_time==None:first_time=0
279
280 time0 = time[0]
281 if time!=first_time:
282 #print "Getting coefficients for", time0
283 [periods,g,h ] = self.__readIGRFcoeff()
284 top_year = numpy.max(periods)
285 nperiod = (top_year - 1900)/5 + 1
286
287 maxcoef = 10
288
289 if time0>=2000:maxcoef = 12
290
291
292 # Normalization array for Schmidt fucntions
293 multer = numpy.zeros((2+maxcoef,1+maxcoef)) + 1
294 for cn in (numpy.arange(maxcoef)+1):
295 for rm in (numpy.arange(cn)+1):
296 tmp = numpy.arange(2*rm) + cn - rm + 1.
297 multer[rm+1,cn] = ((-1.)**rm)*numpy.sqrt(2./tmp.prod())
298
299 schmidt = multer[1:,1:].transpose()
300
301 # n and m arrays
302 nvec = numpy.atleast_2d(numpy.arange(maxcoef)+2)
303 mvec = numpy.atleast_2d(numpy.arange(maxcoef+1)).transpose()
304
305 # Time adjusted igrf g and h with Schmidt normalization
306 # IGRF coefficient arrays: g0(n,m), n=1, maxcoeff,m=0, maxcoeff, ...
307 if time0<top_year:
308 dtime = (time0 - 1900) % 5
309 ntime = int((time0 - 1900 - dtime)/5)
310 else:
311 # Estimating coefficients for times > top_year
312 dtime = (time0 - top_year) + 5
313 ntime = int(g[:,0,0].size - 2)
314
315
316 g0 = g[ntime,1:maxcoef+1,:maxcoef+1]
317 h0 = h[ntime,1:maxcoef+1,:maxcoef+1]
318 gdot = g[ntime+1,1:maxcoef+1,:maxcoef+1]-g[ntime,1:maxcoef+1,:maxcoef+1]
319 hdot = h[ntime+1,1:maxcoef+1,:maxcoef+1]-h[ntime,1:maxcoef+1,:maxcoef+1]
320 gs = (g0 + dtime*(gdot/5.))*schmidt[:maxcoef,0:maxcoef+1]
321 hs = (h0 + dtime*(hdot/5.))*schmidt[:maxcoef,0:maxcoef+1]
322
323 first_time = time0
324
325 for ii in numpy.arange(heights.size):
326 # Height dependence array rad = (ae/(ae+height))**(n+3)
327 rad = numpy.atleast_2d((ae/(ae + heights[ii]))**(nvec+1))
328
329 # Sin and Cos of m times longitude phi arrays
330 mphi = mvec*longitude[ii]*numpy.pi/180.
331 cosmphi = numpy.atleast_2d(numpy.cos(mphi))
332 sinmphi = numpy.atleast_2d(numpy.sin(mphi))
333
334 # Cos of colatitude theta
335 c = numpy.cos((90 - latitude[ii])*numpy.pi/180.)
336
337 # Legendre functions p(n,m|c)
338 [p,dp]= lpmn(maxcoef+1,maxcoef+1,c)
339 p = p[:,:-1].transpose()
340 s = numpy.sqrt((1. - c)*(1 + c))
341
342 # Generate derivative array dpdtheta = -s*dpdc
343 dpdtheta = c*p/s
344 for m in numpy.arange(maxcoef+2): dpdtheta[:,m] = m*dpdtheta[:,m]
345 dpdtheta = dpdtheta + numpy.roll(p,-1,axis=1)
346
347 # Extracting arrays required for field calculations
348 p = p[1:maxcoef+1,:maxcoef+1]
349 dpdtheta = dpdtheta[1:maxcoef+1,:maxcoef+1]
350
351 # Weigh p and dpdtheta with gs and hs coefficients.
352 gp = gs*p
353 hp = hs*p
354 gdpdtheta = gs*dpdtheta
355 hdpdtheta = hs*dpdtheta
356 # Calcultate field components
357 matrix0 = numpy.dot(gdpdtheta,cosmphi)
358 matrix1 = numpy.dot(hdpdtheta,sinmphi)
359 bn[ii] = numpy.dot(rad,(matrix0 + matrix1))
360 matrix0 = numpy.dot(hp,(mvec*cosmphi))
361 matrix1 = numpy.dot(gp,(mvec*sinmphi))
362 be[ii] = numpy.dot((-1*rad),((matrix0 - matrix1)/s))
363 matrix0 = numpy.dot(gp,cosmphi)
364 matrix1 = numpy.dot(hp,sinmphi)
365 bd[ii] = numpy.dot((-1*nvec*rad),(matrix0 + matrix1))
366
367 bmod = numpy.sqrt(bn**2. + be**2. + bd**2.)
368 btheta = numpy.arctan(bd/numpy.sqrt(be**2. + bn**2.))*180/numpy.pi
369 balpha = numpy.arctan(be/bn)*180./numpy.pi
370
371 #bn : north
372 #be : east
373 #bn : radial
374 #bmod : module
375
376
377 return bn, be, bd, bmod, btheta, balpha
378
379 def str2num(self, datum):
380 try:
381 return int(datum)
382 except:
383 try:
384 return float(datum)
385 except:
386 return datum
387
388 def __readIGRFfile(self, filename):
389 list_years=[]
390 for i in range(1,26):
391 list_years.append(1895.0 + i*5)
392
393 epochs=list_years
394 epochs.append(epochs[-1]+5)
395 nepochs = numpy.shape(epochs)
396
397 gg = numpy.zeros((13,14,nepochs[0]),dtype=float)
398 hh = numpy.zeros((13,14,nepochs[0]),dtype=float)
399
400 coeffs_file=open(filename)
401 lines=coeffs_file.readlines()
402
403 coeffs_file.close()
404
405 for line in lines:
406 items = line.split()
407 g_h = items[0]
408 n = self.str2num(items[1])
409 m = self.str2num(items[2])
410
411 coeffs = items[3:]
412
413 for i in range(len(coeffs)-1):
414 coeffs[i] = self.str2num(coeffs[i])
415
416 #coeffs = numpy.array(coeffs)
417 ncoeffs = numpy.shape(coeffs)[0]
418
419 if g_h == 'g':
420 # print n," g ",m
421 gg[n-1,m,:]=coeffs
422 elif g_h=='h':
423 # print n," h ",m
424 hh[n-1,m,:]=coeffs
425 # else :
426 # continue
427
428 # Ultimo Reordenamiento para almacenar .
429 gg[:,:,nepochs[0]-1] = gg[:,:,nepochs[0]-2] + 5*gg[:,:,nepochs[0]-1]
430 hh[:,:,nepochs[0]-1] = hh[:,:,nepochs[0]-2] + 5*hh[:,:,nepochs[0]-1]
431
432 # return numpy.array([gg,hh])
433 periods = numpy.array(epochs)
434 g = gg
435 h = hh
436 return periods, g, h
437
438
439 def __readIGRFcoeff(self,filename="igrf10coeffs.dat"):
440 """
441 __readIGRFcoeff reads the coefficients from a binary file which is located in the
442 folder "resource."
443
444 Parameter
445 ---------
446 filename = A string to specify the name of the file which contains thec coeffs. The
447 default value is "igrf10coeffs.dat"
448
449 Return
450 ------
451 periods = A lineal array giving...
452 g1 =
453 h1 =
454
455 Modification History
456 --------------------
457 Converted to Python by Freddy R. Galindo, ROJ, 03 October 2009.
458 """
459
460 base_path = os.path.dirname(os.path.abspath(__file__))
461 filename = os.path.join(base_path, "igrf13coeffs.txt")
462
463 period_v, g_v, h_v = self.__readIGRFfile(filename)
464 g2 = numpy.zeros((14,14,26))
465 h2 = numpy.zeros((14,14,26))
466 g2[1:14,:,:] = g_v
467 h2[1:14,:,:] = h_v
468
469 g = numpy.transpose(g2, (2,0,1))
470 h = numpy.transpose(h2, (2,0,1))
471 periods = period_v.copy()
472
473 return periods, g, h
474
475 def rotvector(self,vector,axis=1,ang=0):
476 """
477 rotvector function returns the new vector generated rotating the rectagular coords.
478
479 Parameters
480 ----------
481 vector = A lineal 3-elements array (x,y,z).
482 axis = A integer to specify the axis used to rotate the coord systems. The default
483 value is 1.
484 axis = 1 -> Around "x"
485 axis = 2 -> Around "y"
486 axis = 3 -> Around "z"
487 ang = Angle of rotation (in radians). The default value is zero.
488
489 Return
490 ------
491 rotvector = A lineal array of 3 elements giving the new coordinates.
492
493 Modification History
494 --------------------
495 Converted to Python by Freddy R. Galindo, ROJ, 01 October 2009.
496 """
497
498 if axis==1:
499 t = [[1,0,0],[0,numpy.cos(ang),numpy.sin(ang)],[0,-numpy.sin(ang),numpy.cos(ang)]]
500 elif axis==2:
501 t = [[numpy.cos(ang),0,-numpy.sin(ang)],[0,1,0],[numpy.sin(ang),0,numpy.cos(ang)]]
502 elif axis==3:
503 t = [[numpy.cos(ang),numpy.sin(ang),0],[-numpy.sin(ang),numpy.cos(ang),0],[0,0,1]]
504
505 rotvector = numpy.array(numpy.dot(numpy.array(t),numpy.array(vector)))
506
507 return rotvector
Show file before | Show file after | Hide comments
@@ -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
Show file before | Show file after | Hide comments
@@ -0,0 +1,430
1 """
2 The TIME_CONVERSIONS.py module gathers classes and functions for time system transformations
3 (e.g. between seconds from 1970 to datetime format).
4
5 MODULES CALLED:
6 NUMPY, TIME, DATETIME, CALENDAR
7
8 MODIFICATION HISTORY:
9 Created by Ing. Freddy Galindo (frederickgalindo@gmail.com). ROJ Aug 13, 2009.
10 """
11
12 import numpy
13 import time
14 from datetime import datetime as dt
15 import calendar
16
17 class Time:
18 """
19 time(year,month,dom,hour,min,secs)
20
21 An object represents a date and time of certain event..
22
23 Parameters
24 ----------
25 YEAR = Number of the desired year. Year must be valid values from the civil calendar.
26 Years B.C.E must be represented as negative integers. Years in the common era are repre-
27 sented as positive integers. In particular, note that there is no year 0 in the civil
28 calendar. 1 B.C.E. (-1) is followed by 1 C.E. (1).
29
30 MONTH = Number of desired month (1=Jan, ..., 12=December).
31
32 DOM = Number of day of the month.
33
34 HOUR = Number of the hour of the day. By default hour=0
35
36 MINS = Number of the minute of the hour. By default min=0
37
38 SECS = Number of the second of the minute. By default secs=0.
39
40 Examples
41 --------
42 time_info = time(2008,9,30,12,30,00)
43
44 time_info = time(2008,9,30)
45 """
46
47 def __init__(self, year=None, month=None, dom=None, hour=0, mins=0, secs=0):
48 # If one the first three inputs are not defined, it takes the current date.
49 date = time.localtime()
50 if year==None:year=date[0]
51 if month==None:month=date[1]
52 if dom==None:dom=date[2]
53
54 # Converting to arrays
55 year = numpy.array([year]); month = numpy.array([month]); dom = numpy.array([dom])
56 hour = numpy.array([hour]); mins = numpy.array([mins]); secs = numpy.array([secs])
57
58 # Defining time information object.
59 self.year = numpy.atleast_1d(year)
60 self.month = numpy.atleast_1d(month)
61 self.dom = numpy.atleast_1d(dom)
62 self.hour = numpy.atleast_1d(hour)
63 self.mins = numpy.atleast_1d(mins)
64 self.secs = numpy.atleast_1d(secs)
65
66 def change2julday(self):
67 """
68 Converts a datetime to Julian days.
69 """
70
71 # Defining constants
72 greg = 2299171 # incorrect Julian day for Oct, 25, 1582.
73 min_calendar = -4716
74 max_calendar = 5000000
75
76 min_year = numpy.nanmin(self.year)
77 max_year = numpy.nanmax(self.year)
78 if (min_year<min_calendar) or (max_year>max_calendar):
79 print ("Value of Julian date is out of allowed range")
80 return -1
81
82 noyear = numpy.sum(self.year==0)
83 if noyear>0:
84 print ("There is no year zero in the civil calendar")
85 return -1
86
87 # Knowing if the year is less than 0.
88 bc = self.year<0
89
90 # Knowing if the month is less than March.
91 inJanFeb = self.month<=2
92
93 jy = self.year + bc - inJanFeb
94 jm = self.month + (1 + 12*inJanFeb)
95
96 # Computing Julian days.
97 jul= numpy.floor(365.25*jy) + numpy.floor(30.6001*jm) + (self.dom+1720995.0)
98
99 # Test whether to change to Gregorian Calendar
100 if numpy.min(jul) >= greg:
101 ja = numpy.int32(0.01*jy)
102 jul = jul + 2 - ja + numpy.int32(0.25*ja)
103 else:
104 gregchange = numpy.where(jul >= greg)
105 if gregchange[0].size>0:
106 ja = numpy.int32(0.01 + jy[gregchange])
107 jy[gregchange] = jy[gregchange] + 2 - ja + numpy.int32(0.25*ja)
108
109 # Determining machine-specific parameters affecting floating-point.
110 eps = 0.0 # Replace this line for a function to get precision.
111 eps = abs(jul)*0.0 > eps
112
113 jul = jul + (self.hour/24. -0.5) + (self.mins/1440.) + (self.secs/86400.) + eps
114
115 return jul[0]
116
117 def change2secs(self):
118 """
119 Converts datetime to number of seconds respect to 1970.
120 """
121
122 dtime = dt(self.year[0], self.month[0], self.dom[0])
123 return (dtime-dt(1970, 1, 1)).total_seconds()
124
125 year = self.year
126 if year.size>1: year = year[0]
127
128 month = self.month
129 if month.size>1: month = month[0]
130
131 dom = self.dom
132 if dom.size>1: dom = dom[0]
133
134 # Resizing hour, mins and secs if it was necessary.
135 hour = self.hour
136 if hour.size>1:hour = hour[0]
137 if hour.size==1:hour = numpy.resize(hour,year.size)
138
139 mins = self.mins
140 if mins.size>1:mins = mins[0]
141 if mins.size==1:mins = numpy.resize(mins,year.size)
142
143 secs = self.secs
144 if secs.size>1:secs = secs[0]
145 if secs.size==1:secs = numpy.resize(secs,year.size)
146
147 # Using time.mktime to compute seconds respect to 1970.
148 secs1970 = numpy.zeros(year.size)
149 for ii in numpy.arange(year.size):
150 secs1970[ii] = time.mktime((int(year[ii]),int(month[ii]),int(dom[ii]),\
151 int(hour[ii]),int(mins[ii]),int(secs[ii]),0,0,0))
152
153 secs1970 = numpy.int32(secs1970 - time.timezone)
154
155 return secs1970
156
157 def change2strdate(self,mode=1):
158 """
159 change2strdate method converts a date and time of certain event to date string. The
160 string format is like localtime (e.g. Fri Oct 9 15:00:19 2009).
161
162 Parameters
163 ----------
164 None.
165
166 Return
167 ------
168
169 Modification History
170 --------------------
171 Created by Freddy R. Galindo, ROJ, 09 October 2009.
172
173 """
174
175 secs = numpy.atleast_1d(self.change2secs())
176 strdate = []
177 for ii in numpy.arange(numpy.size(secs)):
178 secs_tmp = time.localtime(secs[ii] + time.timezone)
179 if mode==1:
180 strdate.append(time.strftime("%d-%b-%Y (%j) %H:%M:%S",secs_tmp))
181 elif mode==2:
182 strdate.append(time.strftime("%d-%b-%Y (%j)",secs_tmp))
183
184 strdate = numpy.array(strdate)
185
186 return strdate
187
188
189 class Secs:
190 """
191 secs(secs):
192
193 An object represents the number of seconds respect to 1970.
194
195 Parameters
196 ----------
197
198 SECS = A scalar or array giving the number of seconds respect to 1970.
199
200 Example:
201 --------
202 secs_info = secs(1251241373)
203
204 secs_info = secs([1251241373,1251241383,1251241393])
205 """
206 def __init__(self,secs):
207 self.secs = secs
208
209 def change2julday(self):
210 """
211 Convert seconds from 1970 to Julian days.
212 """
213
214 secs_1970 = time(1970,1,1,0,0,0).change2julday()
215
216 julian = self.secs/86400.0 + secs_1970
217
218 return julian
219
220 def change2time(self):
221 """
222 Converts seconds from 1970 to datetime.
223 """
224
225 secs1970 = numpy.atleast_1d(self.secs)
226
227 datetime = numpy.zeros((9,secs1970.size))
228 for ii in numpy.arange(secs1970.size):
229 tuple = time.gmtime(secs1970[ii])
230 datetime[0,ii] = tuple[0]
231 datetime[1,ii] = tuple[1]
232 datetime[2,ii] = tuple[2]
233 datetime[3,ii] = tuple[3]
234 datetime[4,ii] = tuple[4]
235 datetime[5,ii] = tuple[5]
236 datetime[6,ii] = tuple[6]
237 datetime[7,ii] = tuple[7]
238 datetime[8,ii] = tuple[8]
239
240 datetime = numpy.int32(datetime)
241
242 return datetime
243
244
245 class Julian:
246 """
247 julian(julian):
248
249 An object represents julian days.
250
251 Parameters
252 ----------
253
254 JULIAN = A scalar or array giving the julina days.
255
256 Example:
257 --------
258 julian_info = julian(2454740)
259
260 julian_info = julian([2454740,2454760,2454780])
261 """
262 def __init__(self,julian):
263 self.julian = numpy.atleast_1d(julian)
264
265 def change2time(self):
266 """
267 change2time method converts from julian day to calendar date and time.
268
269 Return
270 ------
271 year = An array giving the year of the desired julian day.
272 month = An array giving the month of the desired julian day.
273 dom = An array giving the day of the desired julian day.
274 hour = An array giving the hour of the desired julian day.
275 mins = An array giving the minute of the desired julian day.
276 secs = An array giving the second of the desired julian day.
277
278 Examples
279 --------
280 >> jd = 2455119.0
281 >> [yy,mo,dd,hh,mi,ss] = TimeTools.julian(jd).change2time()
282 >> print [yy,mo,dd,hh,mi,ss]
283 [2009] [10] [ 14.] [ 12.] [ 0.] [ 0.]
284
285 Modification history
286 --------------------
287 Translated from "Numerical Recipies in C", by William H. Press, Brian P. Flannery,
288 Saul A. Teukolsky, and William T. Vetterling. Cambridge University Press, 1988.
289 Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009.
290 """
291
292 min_julian = -1095
293 max_julian = 1827933925
294 if (numpy.min(self.julian) < min_julian) or (numpy.max(self.julian) > max_julian):
295 print ('Value of Julian date is out of allowed range.')
296 return None
297
298 # Beginning of Gregorian calendar
299 igreg = 2299161
300 julLong = numpy.floor(self.julian + 0.5)
301 minJul = numpy.min(julLong)
302
303 if (minJul >= igreg):
304 # All are Gregorian
305 jalpha = numpy.int32(((julLong - 1867216) - 0.25)/36524.25)
306 ja = julLong + 1 + jalpha - numpy.int32(0.25*jalpha)
307 else:
308 ja = julLong
309 gregChange = numpy.where(julLong >= igreg)
310 if gregChange[0].size>0:
311 jalpha = numpy.int32(((julLong[gregChange]-1867216) - 0.25)/36524.25)
312 ja[gregChange] = julLong[gregChange]+1+jalpha-numpy.int32(0.25*jalpha)
313
314 # clear memory.
315 jalpha = -1
316
317 jb = ja + 1524
318 jc = numpy.int32(6680. + ((jb-2439870)-122.1)/365.25)
319 jd = numpy.int32(365.*jc + (0.25*jc))
320 je = numpy.int32((jb - jd)/30.6001)
321
322 dom = jb - jd - numpy.int32(30.6001*je)
323 month = je - 1
324 month = ((month - 1) % 12) + 1
325 month = numpy.atleast_1d(month)
326 year = jc - 4715
327 year = year - (month > 2)*1
328 year = year - (year <= 0)*1
329 year = numpy.atleast_1d(year)
330
331 # Getting hours, minutes, seconds
332 fraction = self.julian + 0.5 - julLong
333 eps_0 = dom*0.0 + 1.0e-12
334 eps_1 = 1.0e-12*numpy.abs(julLong)
335 eps = (eps_0>eps_1)*eps_0 + (eps_0<=eps_1)*eps_1
336
337 hour_0 = dom*0 + 23
338 hour_2 = dom*0 + 0
339 hour_1 = numpy.floor(fraction*24.0 + eps)
340 hour = ((hour_1>hour_0)*23) + ((hour_1<=hour_0)*hour_1)
341 hour = ((hour_1<hour_2)*0) + ((hour_1>=hour_2)*hour_1)
342
343 fraction = fraction - (hour/24.0)
344 mins_0 = dom*0 + 59
345 mins_2 = dom*0 + 0
346 mins_1 = numpy.floor(fraction*1440.0 + eps)
347 mins = ((mins_1>mins_0)*59) + ((mins_1<=mins_0)*mins_1)
348 mins = ((mins_1<mins_2)*0) + ((mins_1>=mins_2)*mins_1)
349
350 secs_2 = dom*0 + 0
351 secs_1 = (fraction - mins/1440.0)*86400.0
352 secs = ((secs_1<secs_2)*0) + ((secs_1>=secs_2)*secs_1)
353
354 return year,month,dom,hour,mins,secs
355
356 def change2secs(self):
357 """
358 Converts from Julian days to seconds from 1970.
359 """
360
361 jul_1970 = Time(1970,1,1,0,0,0).change2julday()
362
363 secs = numpy.int32((self.julian - jul_1970)*86400)
364
365 return secs
366
367 def change2lst(self,longitude=-76.874369):
368 """
369 CT2LST converts from local civil time to local mean sideral time
370
371 longitude = The longitude in degrees (east of Greenwich) of the place for which
372 the local sideral time is desired, scalar. The Greenwich mean sideral time (GMST)
373 can be found by setting longitude=0.
374 """
375
376 # Useful constants, see Meus, p. 84
377 c = numpy.array([280.46061837, 360.98564736629, 0.000387933, 38710000.0])
378 jd2000 = 2451545.0
379 t0 = self.julian - jd2000
380 t = t0/36525.
381
382 # Computing GST in seconds
383 theta = c[0] + (c[1]*t0) + (t**2)*(c[2]-t/c[3])
384
385 # Computing LST in hours
386 lst = (theta + longitude)/15.0
387 neg = numpy.where(lst < 0.0)
388 if neg[0].size>0:lst[neg] = 24.0 + (lst[neg] % 24)
389 lst = lst % 24.0
390
391 return lst
392
393
394 class date2doy:
395 def __init__(self,year,month,day):
396 self.year = year
397 self.month = month
398 self.day = day
399
400 def change2doy(self):
401 if calendar.isleap(self.year) == True:
402 tfactor = 1
403 else:
404 tfactor = 2
405
406 day = self.day
407 month = self.month
408
409 doy = numpy.floor((275*month)/9.0) - (tfactor*numpy.floor((month+9)/12.0)) + day - 30
410
411 return numpy.int32(doy)
412
413
414 class Doy2Date:
415 def __init__(self,year,doy):
416 self.year = year
417 self.doy = doy
418
419 def change2date(self):
420 months = numpy.arange(12) + 1
421
422 first_dem = date2doy(self.year,months,1)
423 first_dem = first_dem.change2doy()
424
425 imm = numpy.where((self.doy - first_dem) > 0)
426
427 month = imm[0].size
428 dom = self.doy -first_dem[month - 1] + 1
429
430 return month, dom
Show file before | Show file after | Hide comments
@@ -258,6 +258,13 class Plot(Operation):
258 258 self.tmin = kwargs.get('tmin', None)
259 259 self.t_units = kwargs.get('t_units', "h_m")
260 260 self.selectedHeightsList = kwargs.get('selectedHeightsList', [])
261 self.extFile = kwargs.get('filename', None)
262 self.bFieldList = kwargs.get('bField', [])
263 self.celestialList = kwargs.get('celestial', [])
264
265 if isinstance(self.bFieldList, int):
266 self.bFieldList = [self.bFieldList]
267
261 268 if isinstance(self.selectedHeightsList, int):
262 269 self.selectedHeightsList = [self.selectedHeightsList]
263 270
Show file before | Show file after | Hide comments
@@ -13,6 +13,10 from schainpy.model.graphics.jroplot_base import Plot, plt, log
13 13 from itertools import combinations
14 14 from matplotlib.ticker import LinearLocator
15 15
16 from schainpy.model.utils.BField import BField
17 from scipy.interpolate import splrep
18 from scipy.interpolate import splev
19
16 20 from matplotlib import __version__ as plt_version
17 21
18 22 if plt_version >='3.3.4':
@@ -67,6 +71,11 class SpectraPlot(Plot):
67 71 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
68 72 noise = 10*numpy.log10(dataOut.getNoise()/norm)
69 73
74 <<<<<<< HEAD
75 =======
76
77
78 >>>>>>> 37cccf17c7b80521b59b978cb30e4ab2e6f37fce
70 79 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
71 80 for ch in range(dataOut.nChannels):
72 81 if hasattr(dataOut.normFactor,'ndim'):
@@ -1466,4 +1475,251 class GeneralProfilePlot(Plot):
1466 1475 #self.xmax = max(self.z)
1467 1476 ax.plt_r = ax.plot(self.z[i], self.y)[0]
1468 1477 else:
1469 ax.plt_r.set_data(self.z[i], self.y) No newline at end of file
1478 ax.plt_r.set_data(self.z[i], self.y)
1479
1480
1481 ##########################################################################################################
1482 ########################################## AMISR_V4 ######################################################
1483
1484 class RTIMapPlot(Plot):
1485 '''
1486 Plot for RTI data
1487
1488 Example:
1489
1490 controllerObj = Project()
1491 controllerObj.setup(id = '11', name='eej_proc', description=desc)
1492 ##.......................................................................................
1493 ##.......................................................................................
1494 readUnitConfObj = controllerObj.addReadUnit(datatype='AMISRReader', path=inPath, startDate='2023/05/24',endDate='2023/05/24',
1495 startTime='12:00:00',endTime='12:45:59',walk=1,timezone='lt',margin_days=1,code = code,nCode = nCode,
1496 nBaud = nBaud,nOsamp = nosamp,nChannels=nChannels,nFFT=NFFT,
1497 syncronization=False,shiftChannels=0)
1498
1499 volts_proc = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
1500
1501 opObj01 = volts_proc.addOperation(name='Decoder', optype='other')
1502 opObj01.addParameter(name='code', value=code, format='floatlist')
1503 opObj01.addParameter(name='nCode', value=1, format='int')
1504 opObj01.addParameter(name='nBaud', value=nBaud, format='int')
1505 opObj01.addParameter(name='osamp', value=nosamp, format='int')
1506
1507 opObj12 = volts_proc.addOperation(name='selectHeights', optype='self')
1508 opObj12.addParameter(name='minHei', value='90', format='float')
1509 opObj12.addParameter(name='maxHei', value='150', format='float')
1510
1511 proc_spc = controllerObj.addProcUnit(datatype='SpectraProc', inputId=volts_proc.getId())
1512 proc_spc.addParameter(name='nFFTPoints', value='8', format='int')
1513
1514 opObj11 = proc_spc.addOperation(name='IncohInt', optype='other')
1515 opObj11.addParameter(name='n', value='1', format='int')
1516
1517 beamMapFile = "/home/japaza/Documents/AMISR_sky_mapper/UMET_beamcodes.csv"
1518
1519 opObj12 = proc_spc.addOperation(name='RTIMapPlot', optype='external')
1520 opObj12.addParameter(name='selectedHeightsList', value='95, 100, 105, 110 ', format='int')
1521 opObj12.addParameter(name='bField', value='100', format='int')
1522 opObj12.addParameter(name='filename', value=beamMapFile, format='str')
1523
1524 '''
1525
1526 CODE = 'rti_skymap'
1527
1528 plot_type = 'scatter'
1529 titles = None
1530 colormap = 'jet'
1531 channelList = []
1532 elevationList = []
1533 azimuthList = []
1534 last_noise = None
1535 flag_setIndex = False
1536 heights = []
1537 dcosx = []
1538 dcosy = []
1539 fullDcosy = None
1540 fullDcosy = None
1541 hindex = []
1542 mapFile = False
1543 ##### BField ####
1544 flagBField = False
1545 dcosxB = []
1546 dcosyB = []
1547 Bmarker = ['+','*','D','x','s','>','o','^']
1548
1549
1550
1551 def setup(self):
1552
1553 self.xaxis = 'Range (Km)'
1554 if len(self.selectedHeightsList) > 0:
1555 self.nplots = len(self.selectedHeightsList)
1556 else:
1557 self.nplots = 4
1558 self.ncols = int(numpy.ceil(self.nplots/2))
1559 self.nrows = int(numpy.ceil(self.nplots/self.ncols))
1560 self.ylabel = 'dcosy'
1561 self.xlabel = 'dcosx'
1562 self.colorbar = True
1563 self.width = 6 + 4.1*self.nrows
1564 self.height = 3 + 3.5*self.ncols
1565
1566
1567 if self.extFile!=None:
1568 try:
1569 pointings = numpy.genfromtxt(self.extFile, delimiter=',')
1570 full_azi = pointings[:,1]
1571 full_elev = pointings[:,2]
1572 self.fullDcosx = numpy.cos(numpy.radians(full_elev))*numpy.sin(numpy.radians(full_azi))
1573 self.fullDcosy = numpy.cos(numpy.radians(full_elev))*numpy.cos(numpy.radians(full_azi))
1574 mapFile = True
1575 except Exception as e:
1576 self.extFile = None
1577 print(e)
1578
1579
1580
1581
1582 def update_list(self,dataOut):
1583 if len(self.channelList) == 0:
1584 self.channelList = dataOut.channelList
1585 if len(self.elevationList) == 0:
1586 self.elevationList = dataOut.elevationList
1587 if len(self.azimuthList) == 0:
1588 self.azimuthList = dataOut.azimuthList
1589 a = numpy.radians(numpy.asarray(self.azimuthList))
1590 e = numpy.radians(numpy.asarray(self.elevationList))
1591 self.heights = dataOut.heightList
1592 self.dcosx = numpy.cos(e)*numpy.sin(a)
1593 self.dcosy = numpy.cos(e)*numpy.cos(a)
1594
1595 if len(self.bFieldList)>0:
1596 datetObj = datetime.datetime.fromtimestamp(dataOut.utctime)
1597 doy = datetObj.timetuple().tm_yday
1598 year = datetObj.year
1599 # self.dcosxB, self.dcosyB
1600 ObjB = BField(year=year,doy=doy,site=2,heights=self.bFieldList)
1601 [dcos, alpha, nlon, nlat] = ObjB.getBField()
1602
1603 alpha_location = numpy.zeros((nlon,2,len(self.bFieldList)))
1604 for ih in range(len(self.bFieldList)):
1605 alpha_location[:,0,ih] = dcos[:,0,ih,0]
1606 for ilon in numpy.arange(nlon):
1607 myx = (alpha[ilon,:,ih])[::-1]
1608 myy = (dcos[ilon,:,ih,0])[::-1]
1609 tck = splrep(myx,myy,s=0)
1610 mydcosx = splev(ObjB.alpha_i,tck,der=0)
1611
1612 myx = (alpha[ilon,:,ih])[::-1]
1613 myy = (dcos[ilon,:,ih,1])[::-1]
1614 tck = splrep(myx,myy,s=0)
1615 mydcosy = splev(ObjB.alpha_i,tck,der=0)
1616 alpha_location[ilon,:,ih] = numpy.array([mydcosx, mydcosy])
1617 self.dcosxB.append(alpha_location[:,0,ih])
1618 self.dcosyB.append(alpha_location[:,1,ih])
1619 self.flagBField = True
1620
1621 if len(self.celestialList)>0:
1622 #getBField(self.bFieldList, date)
1623 #pass = kwargs.get('celestial', [])
1624 pass
1625
1626
1627
1628 def update(self, dataOut):
1629
1630 if len(self.channelList) == 0:
1631 self.update_list(dataOut)
1632
1633 if not self.flag_setIndex:
1634 if len(self.selectedHeightsList)>0:
1635 for sel_height in self.selectedHeightsList:
1636 index_list = numpy.where(self.heights >= sel_height)
1637 index_list = index_list[0]
1638 self.hindex.append(index_list[0])
1639 # else:
1640 # k = len(self.heights)
1641 # self.hindex.append(int(k/2))
1642 self.flag_setIndex = True
1643
1644 data = {}
1645 meta = {}
1646
1647 data['rti_skymap'] = dataOut.getPower()
1648 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1649 noise = 10*numpy.log10(dataOut.getNoise()/norm)
1650 data['noise'] = noise
1651
1652 return data, meta
1653
1654 def plot(self):
1655
1656 ######
1657 self.x = self.dcosx
1658 self.y = self.dcosy
1659 self.z = self.data[-1]['rti_skymap']
1660 self.z = numpy.array(self.z, dtype=float)
1661
1662 #print("inde x1 ", self.height_index)
1663 if len(self.hindex) > 0:
1664 index = self.hindex
1665 else:
1666 index = numpy.arange(0, len(self.heights), int((len(self.heights))/4.2))
1667
1668 #print(index)
1669 self.titles = ['Height {:.2f} km '.format(self.heights[i])+" " for i in index]
1670 for n, ax in enumerate(self.axes):
1671
1672 if ax.firsttime:
1673
1674
1675 self.xmax = self.xmax if self.xmax else numpy.nanmax(self.x)
1676 self.xmin = self.xmin if self.xmin else numpy.nanmin(self.x)
1677
1678 self.ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
1679 self.ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
1680
1681 self.zmax = self.zmax if self.zmax else numpy.nanmax(self.z)
1682 self.zmin = self.zmin if self.zmin else numpy.nanmin(self.z)
1683
1684
1685 if self.extFile!=None:
1686 ax.scatter(self.fullDcosx, self.fullDcosy, marker="+", s=20)
1687 #print(self.fullDcosx)
1688 pass
1689
1690
1691 ax.plt = ax.scatter(self.x, self.y, c=self.z[:,index[n]], cmap = 'jet',vmin = self.zmin,
1692 s=60, marker="s", vmax = self.zmax)
1693
1694
1695 ax.minorticks_on()
1696 ax.grid(which='major', axis='both')
1697 ax.grid(which='minor', axis='x')
1698
1699 if self.flagBField :
1700
1701 for ih in range(len(self.bFieldList)):
1702 label = str(self.bFieldList[ih]) + ' km'
1703 ax.plot(self.dcosxB[ih], self.dcosyB[ih], color='k', marker=self.Bmarker[ih % 8],
1704 label=label, linestyle='--', ms=4.0,lw=0.5)
1705 handles, labels = ax.get_legend_handles_labels()
1706 a = -0.05
1707 b = 1.15 - 1.19*(self.nrows)
1708 self.axes[0].legend(handles,labels, bbox_to_anchor=(a,b), prop={'size': (5.8+ 1.1*self.nplots)}, title='B Field ⊥')
1709
1710 else:
1711
1712 ax.plt = ax.scatter(self.x, self.y, c=self.z[:,index[n]], cmap = 'jet',vmin = self.zmin,
1713 s=80, marker="s", vmax = self.zmax)
1714
1715 if self.flagBField :
1716 for ih in range(len(self.bFieldList)):
1717 ax.plot (self.dcosxB[ih], self.dcosyB[ih], color='k', marker=self.Bmarker[ih % 8],
1718 linestyle='--', ms=4.0,lw=0.5)
1719
1720 # handles, labels = ax.get_legend_handles_labels()
1721 # a = -0.05
1722 # b = 1.15 - 1.19*(self.nrows)
1723 # self.axes[0].legend(handles,labels, bbox_to_anchor=(a,b), prop={'size': (5.8+ 1.1*self.nplots)}, title='B Field ⊥')
1724
1725
Show file before | Show file after | Hide comments
@@ -305,6 +305,15 class PulsepairSignalPlot(ScopePlot):
305 305 class Spectra2DPlot(Plot):
306 306 '''
307 307 Plot for 2D Spectra data
308 Necessary data as Block
309 you could use profiles2Block Operation
310
311 Example:
312 # opObj11 = volts_proc.addOperation(name='profiles2Block', optype='other')
313 # # opObj11.addParameter(name='n', value=10, format='int')
314 # opObj11.addParameter(name='timeInterval', value='2', format='int')
315
316 # opObj12 = volts_proc.addOperation(name='Spectra2DPlot', optype='external')
308 317 '''
309 318
310 319 CODE = 'spc'
Show file before | Show file after | Hide comments
@@ -61,8 +61,9 class AMISRReader(ProcessingUnit):
61 61 self.elevationList = []
62 62 self.dataShape = None
63 63 self.flag_old_beams = False
64
65
64
65 self.flagAsync = False #Use when the experiment has no syncronization
66 self.shiftChannels = 0
66 67 self.profileIndex = 0
67 68
68 69
@@ -109,6 +110,8 class AMISRReader(ProcessingUnit):
109 110 ignEndDate=None,
110 111 ignStartTime=None,
111 112 ignEndTime=None,
113 syncronization=True,
114 shiftChannels=0
112 115 ):
113 116
114 117
@@ -123,7 +126,8 class AMISRReader(ProcessingUnit):
123 126 self.nOsamp = int(nOsamp)
124 127 self.margin_days = margin_days
125 128 self.__sampleRate = None
126
129 self.flagAsync = not syncronization
130 self.shiftChannels = shiftChannels
127 131 self.nFFT = nFFT
128 132 self.nChannels = nChannels
129 133 if ignStartTime!=None and ignEndTime!=None:
@@ -178,7 +182,14 class AMISRReader(ProcessingUnit):
178 182 a = [line for line in linesExp if "nbeamcodes" in line]
179 183 self.nChannels = int(a[0][11:])
180 184
185 if not self.flagAsync: #for experiments with no syncronization
186 self.shiftChannels = 0
187
188
189
181 190 self.beamCodeByPulse = fp.get(header+'/BeamCode') # LIST OF BEAMS PER PROFILE, TO BE USED ON REARRANGE
191
192
182 193 if (self.startDate > datetime.date(2021, 7, 15)) or self.flag_old_beams: #Se cambió la forma de extracción de Apuntes el 17 o forzar con flag de reorganización
183 194 self.beamcodeFile = fp['Setup/Beamcodefile'][()].decode()
184 195 self.trueBeams = self.beamcodeFile.split("\n")
@@ -189,10 +200,17 class AMISRReader(ProcessingUnit):
189 200 beams = [self.trueBeams[b] for b in beams_idx]
190 201 self.beamCode = [int(x, 16) for x in beams]
191 202
203 if(self.flagAsync and self.shiftChannels == 0):
204 initBeam = self.beamCodeByPulse[0, 0]
205 self.shiftChannels = numpy.argwhere(self.beamCode ==initBeam)[0,0]
206
192 207 else:
193 208 _beamCode= fp.get('Raw11/Data/Beamcodes') #se usa la manera previa al cambio de apuntes
194 209 self.beamCode = _beamCode[0,:]
195 210
211
212
213
196 214 if self.beamCodeMap == None:
197 215 self.beamCodeMap = fp['Setup/BeamcodeMap']
198 216 for beam in self.beamCode:
@@ -219,7 +237,8 class AMISRReader(ProcessingUnit):
219 237 self.nblocks = self.pulseCount.shape[0] #nblocks
220 238 self.profPerBlockRAW = self.pulseCount.shape[1] #profiles per block in raw data
221 239 self.nprofiles = self.pulseCount.shape[1] #nprofile
222 self.nsa = self.nsamplesPulse[0,0] #ngates
240 #self.nsa = self.nsamplesPulse[0,0] #ngates
241 self.nsa = len(self.rangeFromFile[0])
223 242 self.nchannels = len(self.beamCode)
224 243 self.ippSeconds = (self.radacTime[0][1] -self.radacTime[0][0]) #Ipp in seconds
225 244 #print("IPPS secs: ",self.ippSeconds)
@@ -521,8 +540,9 class AMISRReader(ProcessingUnit):
521 540 #profPerCH = int(self.profPerBlockRAW / self.nChannels)
522 541 for thisChannel in range(nchan):
523 542
543 ich = thisChannel
524 544
525 idx_ch = [self.nFFT*(thisChannel + nchan*k) for k in range(profPerCH)]
545 idx_ch = [self.nFFT*(ich + nchan*k) for k in range(profPerCH)]
526 546 #print(idx_ch)
527 547 if self.nFFT > 1:
528 548 aux = [numpy.arange(i, i+self.nFFT) for i in idx_ch]
@@ -532,14 +552,15 class AMISRReader(ProcessingUnit):
532 552 else:
533 553 idx_ch = numpy.array(idx_ch, dtype=int)
534 554
535 #print(thisChannel,profPerCH,idx_ch)
536 #print(numpy.where(channels==self.beamCode[thisChannel])[0])
537 #new_block[:,thisChannel,:,:] = self.dataset[:,numpy.where(channels==self.beamCode[thisChannel])[0],:]
538 new_block[:,thisChannel,:,:] = self.dataset[:,idx_ch,:]
555 #print(ich,profPerCH,idx_ch)
556 #print(numpy.where(channels==self.beamCode[ich])[0])
557 #new_block[:,ich,:,:] = self.dataset[:,numpy.where(channels==self.beamCode[ich])[0],:]
558 new_block[:,ich,:,:] = self.dataset[:,idx_ch,:]
539 559
540 560 new_block = numpy.transpose(new_block, (1,0,2,3))
541 561 new_block = numpy.reshape(new_block, (nchan,-1, nsamples))
542
562 if self.flagAsync:
563 new_block = numpy.roll(new_block, self.shiftChannels, axis=0)
543 564 return new_block
544 565
545 566 def updateIndexes(self):
@@ -575,10 +596,13 class AMISRReader(ProcessingUnit):
575 596 # self.dataOut.channelIndexList = None
576 597
577 598
578 self.dataOut.azimuthList = numpy.array(self.azimuthList)
579 self.dataOut.elevationList = numpy.array(self.elevationList)
580 self.dataOut.codeList = numpy.array(self.beamCode)
581
599 #self.dataOut.azimuthList = numpy.roll( numpy.array(self.azimuthList) ,self.shiftChannels)
600 #self.dataOut.elevationList = numpy.roll(numpy.array(self.elevationList) ,self.shiftChannels)
601 #self.dataOut.codeList = numpy.roll(numpy.array(self.beamCode), self.shiftChannels)
602
603 self.dataOut.azimuthList = self.azimuthList
604 self.dataOut.elevationList = self.elevationList
605 self.dataOut.codeList = self.beamCode
582 606
583 607
584 608
@@ -637,7 +661,11 class AMISRReader(ProcessingUnit):
637 661 self.dataOut.radarControllerHeaderObj.elevationList = self.elevationList
638 662 self.dataOut.radarControllerHeaderObj.dtype = "Voltage"
639 663 self.dataOut.ippSeconds = self.ippSeconds
664 <<<<<<< HEAD
640 665 self.dataOut.ippFactor = self.nchannels*self.nFFT
666 =======
667 self.dataOut.ippFactor = self.nFFT
668 >>>>>>> 37cccf17c7b80521b59b978cb30e4ab2e6f37fce
641 669 pass
642 670
643 671 def readNextFile(self,online=False):
General Comments 0
You need to be logged in to leave comments. Login now