+
+
+{% endblock %}
diff --git a/apps/abs/utils/Astro_Coords.py b/apps/abs/utils/Astro_Coords.py
new file mode 100644
index 0000000..73ca019
--- /dev/null
+++ b/apps/abs/utils/Astro_Coords.py
@@ -0,0 +1,1419 @@
+"""
+The module ASTRO_COORDS.py gathers classes and functions for coordinates transformation. Additiona-
+lly a class EquatorialCorrections and celestial bodies are defined. The first of these is to correct
+any error in the location of the body and the second to know the location of certain celestial bo-
+dies in the sky.
+
+MODULES CALLED:
+OS, NUMPY, NUMERIC, SCIPY, TIME_CONVERSIONS
+
+MODIFICATION HISTORY:
+Created by Ing. Freddy Galindo (frederickgalindo@gmail.com). ROJ Sep 20, 2009.
+"""
+
+import numpy
+#import Numeric
+import scipy.interpolate
+import os
+import sys
+import TimeTools
+import Misc_Routines
+
+class EquatorialCorrections():
+ def __init__(self):
+ """
+ EquatorialCorrections class creates an object to call methods to correct the loca-
+ tion of the celestial bodies.
+
+ Modification History
+ --------------------
+ Converted to Object-oriented Programming by Freddy Galindo, ROJ, 27 September 2009.
+ """
+
+ pass
+
+ def co_nutate(self,jd,ra,dec):
+ """
+ co_nutate calculates changes in RA and Dec due to nutation of the Earth's rotation
+ Additionally it returns the obliquity of the ecliptic (eps), nutation in the longi-
+ tude of the ecliptic (d_psi) and nutation in the pbliquity of the ecliptic (d_eps).
+
+ Parameters
+ ----------
+ jd = Julian Date (Scalar or array).
+ RA = A scalar o array giving the Right Ascention of interest.
+ Dec = A scalar o array giving the Right Ascention of interest.
+
+ Return
+ ------
+ d_ra = Correction to ra due to nutation.
+ d_dec = Correction to dec due to nutation.
+
+ Examples
+ --------
+ >> Julian = 2462088.7
+ >> Ra = 41.547213
+ >> Dec = 49.348483
+ >> [d_ra,d_dec,eps,d_psi,d_eps] = co_nutate(julian,Ra,Dec)
+ >> print d_ra, d_dec, eps, d_psi, d_eps
+ [ 15.84276651] [ 6.21641029] [ 0.4090404] [ 14.85990198] [ 2.70408658]
+
+ Modification history
+ --------------------
+ Written by Chris O'Dell, 2002.
+ Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009.
+ """
+
+ jd = numpy.atleast_1d(jd)
+ ra = numpy.atleast_1d(ra)
+ dec = numpy.atleast_1d(dec)
+
+ # Useful transformation constants
+ d2as = numpy.pi/(180.*3600.)
+
+ # Julian centuries from J2000 of jd
+ T = (jd - 2451545.0)/36525.0
+
+ # Must calculate obliquity of ecliptic
+ [d_psi, d_eps] = self.nutate(jd)
+ d_psi = numpy.atleast_1d(d_psi)
+ d_eps = numpy.atleast_1d(d_eps)
+
+ eps0 = (23.4392911*3600.) - (46.8150*T) - (0.00059*T**2) + (0.001813*T**3)
+ # True obliquity of the ecliptic in radians
+ eps = (eps0 + d_eps)/3600.*Misc_Routines.CoFactors.d2r
+
+ # Useful numbers
+ ce = numpy.cos(eps)
+ se = numpy.sin(eps)
+
+ # Convert Ra-Dec to equatorial rectangular coordinates
+ x = numpy.cos(ra*Misc_Routines.CoFactors.d2r)*numpy.cos(dec*Misc_Routines.CoFactors.d2r)
+ y = numpy.sin(ra*Misc_Routines.CoFactors.d2r)*numpy.cos(dec*Misc_Routines.CoFactors.d2r)
+ z = numpy.sin(dec*Misc_Routines.CoFactors.d2r)
+
+ # Apply corrections to each rectangular coordinate
+ x2 = x - (y*ce + z*se)*d_psi*Misc_Routines.CoFactors.s2r
+ y2 = y + (x*ce*d_psi - z*d_eps)*Misc_Routines.CoFactors.s2r
+ z2 = z + (x*se*d_psi + y*d_eps)*Misc_Routines.CoFactors.s2r
+
+ # Convert bask to equatorial spherical coordinates
+ r = numpy.sqrt(x2**2. + y2**2. + z2**2.)
+ xyproj =numpy.sqrt(x2**2. + y2**2.)
+
+ ra2 = x2*0.0
+ dec2 = x2*0.0
+
+ xyproj = numpy.atleast_1d(xyproj)
+ z = numpy.atleast_1d(z)
+ r = numpy.atleast_1d(r)
+ x2 = numpy.atleast_1d(x2)
+ y2 = numpy.atleast_1d(y2)
+ z2 = numpy.atleast_1d(z2)
+ ra2 = numpy.atleast_1d(ra2)
+ dec2 = numpy.atleast_1d(dec2)
+
+ w1 = numpy.where((xyproj==0) & (z!=0))
+ w2 = numpy.where(xyproj!=0)
+
+ # Calculate Ra and Dec in radians (later convert to degrees)
+ if w1[0].size>0:
+ # Places where xyproj=0 (point at NCP or SCP)
+ dec2[w1] = numpy.arcsin(z2[w1]/r[w1])
+ ra2[w1] = 0
+
+ if w2[0].size>0:
+ # Places other than NCP or SCP
+ ra2[w2] = numpy.arctan2(y2[w2],x2[w2])
+ dec2[w2] = numpy.arcsin(z2[w2]/r[w2])
+
+ # Converting to degree
+ ra2 = ra2/Misc_Routines.CoFactors.d2r
+ dec2 = dec2/Misc_Routines.CoFactors.d2r
+
+ w = numpy.where(ra2<0.)
+ if w[0].size>0:
+ ra2[w] = ra2[w] + 360.
+
+ # Return changes in Ra and Dec in arcseconds
+ d_ra = (ra2 -ra)*3600.
+ d_dec = (dec2 - dec)*3600.
+
+ return d_ra, d_dec, eps, d_psi, d_eps
+
+ def nutate(self,jd):
+ """
+ nutate returns the nutation in longitude and obliquity for a given Julian date.
+
+ Parameters
+ ----------
+ jd = Julian ephemeris date, scalar or vector.
+
+ Return
+ ------
+ nut_long = The nutation in longitude.
+ nut_obliq = The nutation in latitude.
+
+ Example
+ -------
+ >> julian = 2446895.5
+ >> [nut_long,nut_obliq] = nutate(julian)
+ >> print nut_long, nut_obliq
+ -3.78793107711 9.44252069864
+
+ >> julians = 2415020.5 + numpy.arange(50)
+ >> [nut_long,nut_obliq] = nutate(julians)
+
+ Modification History
+ --------------------
+ Written by W.Landsman (Goddard/HSTX), June 1996.
+ Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009.
+ """
+
+ jd = numpy.atleast_1d(jd)
+
+ # Form time in Julian centuries from 1900
+ t = (jd - 2451545.0)/36525.0
+
+ # Mean elongation of the moon
+ coeff1 = numpy.array([1/189474.0,-0.0019142,445267.111480,297.85036])
+ d = numpy.poly1d(coeff1)
+ d = d(t)*Misc_Routines.CoFactors.d2r
+ d = self.cirrange(d,rad=1)
+
+ # Sun's mean elongation
+ coeff2 = numpy.array([-1./3e5,-0.0001603,35999.050340,357.52772])
+ m = numpy.poly1d(coeff2)
+ m = m(t)*Misc_Routines.CoFactors.d2r
+ m = self.cirrange(m,rad=1)
+
+ # Moon's mean elongation
+ coeff3 = numpy.array([1.0/5.625e4,0.0086972,477198.867398,134.96298])
+ mprime = numpy.poly1d(coeff3)
+ mprime = mprime(t)*Misc_Routines.CoFactors.d2r
+ mprime = self.cirrange(mprime,rad=1)
+
+ # Moon's argument of latitude
+ coeff4 = numpy.array([-1.0/3.27270e5,-0.0036825,483202.017538,93.27191])
+ f = numpy.poly1d(coeff4)
+ f = f(t)*Misc_Routines.CoFactors.d2r
+ f = self.cirrange(f,rad=1)
+
+ # Longitude fo the ascending node of the Moon's mean orbit on the ecliptic, measu-
+ # red from the mean equinox of the date.
+ coeff5 = numpy.array([1.0/4.5e5,0.0020708,-1934.136261,125.04452])
+ omega = numpy.poly1d(coeff5)
+ omega = omega(t)*Misc_Routines.CoFactors.d2r
+ omega = self.cirrange(omega,rad=1)
+
+ 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,\
+ 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,\
+ 0,2,0,2])
+
+ m_lng = numpy.array([0,0,0,0,1,0,1,0,0,-1])
+ m_lng = numpy.append(m_lng,numpy.zeros(17))
+ 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,\
+ 0,1,0,0,1,0,0,0,-1,1,-1,-1,0,-1]))
+
+ 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,\
+ -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])
+
+ 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,\
+ 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])
+
+ 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, \
+ 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])
+
+ sin_lng = numpy.array([-171996,-13187,-2274,2062,1426,712,-517,-386,-301, 217, -158, \
+ 129,123,63,63,-59,-58,-51,48,46,-38,-31,29,29,26,-22,21,17,16,-16,-15,-13,\
+ -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,\
+ -3,-3,-3,-3])
+
+ 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,\
+ 0,-0.1])
+ sdelt = numpy.append(sdelt,numpy.zeros(10))
+ sdelt = numpy.append(sdelt,numpy.array([-0.1, 0, 0.1]))
+ sdelt = numpy.append(sdelt,numpy.zeros(33))
+
+ cos_lng = numpy.array([92025,5736,977,-895,54,-7,224,200,129,-95,0,-70,-53,0,-33,26, \
+ 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,\
+ 3,3,3])
+ cos_lng = numpy.append(cos_lng,numpy.zeros(14))
+
+ cdelt = numpy.array([8.9,-3.1,-0.5,0.5,-0.1,0.0,-0.6,0.0,-0.1,0.3])
+ cdelt = numpy.append(cdelt,numpy.zeros(53))
+
+ # Sum the periodic terms.
+ n = numpy.size(jd)
+ nut_long = numpy.zeros(n)
+ nut_obliq = numpy.zeros(n)
+
+ d_lng = d_lng.reshape(numpy.size(d_lng),1)
+ d = d.reshape(numpy.size(d),1)
+ matrix_d_lng = numpy.dot(d_lng,d.transpose())
+
+ m_lng = m_lng.reshape(numpy.size(m_lng),1)
+ m = m.reshape(numpy.size(m),1)
+ matrix_m_lng = numpy.dot(m_lng,m.transpose())
+
+ mp_lng = mp_lng.reshape(numpy.size(mp_lng),1)
+ mprime = mprime.reshape(numpy.size(mprime),1)
+ matrix_mp_lng = numpy.dot(mp_lng,mprime.transpose())
+
+ f_lng = f_lng.reshape(numpy.size(f_lng),1)
+ f = f.reshape(numpy.size(f),1)
+ matrix_f_lng = numpy.dot(f_lng,f.transpose())
+
+ om_lng = om_lng.reshape(numpy.size(om_lng),1)
+ omega = omega.reshape(numpy.size(omega),1)
+ matrix_om_lng = numpy.dot(om_lng,omega.transpose())
+
+ arg = matrix_d_lng + matrix_m_lng + matrix_mp_lng + matrix_f_lng + matrix_om_lng
+
+ sarg = numpy.sin(arg)
+ carg = numpy.cos(arg)
+
+ for ii in numpy.arange(n):
+ nut_long[ii] = 0.0001*numpy.sum((sdelt*t[ii] + sin_lng)*sarg[:,ii])
+ nut_obliq[ii] = 0.0001*numpy.sum((cdelt*t[ii] + cos_lng)*carg[:,ii])
+
+ if numpy.size(jd)==1:
+ nut_long = nut_long[0]
+ nut_obliq = nut_obliq[0]
+
+ return nut_long, nut_obliq
+
+ def co_aberration(self,jd,ra,dec):
+ """
+ co_aberration calculates changes to Ra and Dec due to "the effect of aberration".
+
+ Parameters
+ ----------
+ jd = Julian Date (Scalar or vector).
+ ra = A scalar o vector giving the Right Ascention of interest.
+ dec = A scalar o vector giving the Declination of interest.
+
+ Return
+ ------
+ d_ra = The correction to right ascension due to aberration (must be added to ra to
+ get the correct value).
+ d_dec = The correction to declination due to aberration (must be added to the dec
+ to get the correct value).
+ eps = True obliquity of the ecliptic (in radians).
+
+ Examples
+ --------
+ >> Julian = 2462088.7
+ >> Ra = 41.547213
+ >> Dec = 49.348483
+ >> [d_ra,d_dec,eps] = co_aberration(julian,Ra,Dec)
+ >> print d_ra, d_dec, eps
+ [ 30.04441796] [ 6.69837858] [ 0.40904059]
+
+ Modification history
+ --------------------
+ Written by Chris O'Dell , Univ. of Wisconsin, June 2002.
+ Converted to Python by Freddy R. Galindo, ROJ, 27 September 2009.
+ """
+
+ # Julian centuries from J2000 of jd.
+ T = (jd - 2451545.0)/36525.0
+
+ # Getting obliquity of ecliptic
+ njd = numpy.size(jd)
+ jd = numpy.atleast_1d(jd)
+ ra = numpy.atleast_1d(ra)
+ dec = numpy.atleast_1d(dec)
+
+ d_psi = numpy.zeros(njd)
+ d_epsilon = d_psi
+ for ii in numpy.arange(njd):
+ [dp,de] = self.nutate(jd[ii])
+ d_psi[ii] = dp
+ d_epsilon[ii] = de
+
+ coeff = 23 + 26/60. + 21.488/3600.
+ eps0 = coeff*3600. - 46.8150*T - 0.00059*T**2. + 0.001813*T**3.
+ # True obliquity of the ecliptic in radians
+ eps = (eps0 + d_epsilon)/3600*Misc_Routines.CoFactors.d2r
+
+ celestialbodies = CelestialBodies()
+ [sunra,sundec,sunlon,sunobliq] = celestialbodies.sunpos(jd)
+
+ # Earth's orbital eccentricity
+ e = 0.016708634 - 0.000042037*T - 0.0000001267*T**2.
+
+ # longitude of perihelion, in degrees
+ pi = 102.93735 + 1.71946*T + 0.00046*T**2
+
+ # Constant of aberration, in arcseconds
+ k = 20.49552
+
+ cd = numpy.cos(dec*Misc_Routines.CoFactors.d2r) ; sd = numpy.sin(dec*Misc_Routines.CoFactors.d2r)
+ ce = numpy.cos(eps) ; te = numpy.tan(eps)
+ cp = numpy.cos(pi*Misc_Routines.CoFactors.d2r) ; sp = numpy.sin(pi*Misc_Routines.CoFactors.d2r)
+ cs = numpy.cos(sunlon*Misc_Routines.CoFactors.d2r) ; ss = numpy.sin(sunlon*Misc_Routines.CoFactors.d2r)
+ ca = numpy.cos(ra*Misc_Routines.CoFactors.d2r) ; sa = numpy.sin(ra*Misc_Routines.CoFactors.d2r)
+
+ term1 = (ca*cs*ce + sa*ss)/cd
+ term2 = (ca*cp*ce + sa*sp)/cd
+ term3 = (cs*ce*(te*cd - sa*sd) + ca*sd*ss)
+ term4 = (cp*ce*(te*cd - sa*sd) + ca*sd*sp)
+
+ d_ra = -k*term1 + e*k*term2
+ d_dec = -k*term3 + e*k*term4
+
+ return d_ra, d_dec, eps
+
+ def precess(self,ra,dec,equinox1=None,equinox2=None,FK4=0,rad=0):
+ """
+ precess coordinates from EQUINOX1 to EQUINOX2
+
+ Parameters
+ -----------
+ ra = A scalar o vector giving the Right Ascention of interest.
+ dec = A scalar o vector giving the Declination of interest.
+ equinox1 = Original equinox of coordinates, numeric scalar. If omitted, the __Pre-
+ cess will query for equinox1 and equinox2.
+ equinox2 = Original equinox of coordinates.
+ FK4 = If this keyword is set and non-zero, the FK4 (B1950) system will be used
+ otherwise FK5 (J2000) will be used instead.
+ rad = If this keyword is set and non-zero, then the input and output RAD and DEC
+ vectors are in radian rather than degree.
+
+ Return
+ ------
+ ra = Right ascension after precession (scalar or vector) in degrees, unless the rad
+ keyword is set.
+ dec = Declination after precession (scalar or vector) in degrees, unless the rad
+ keyword is set.
+
+ Examples
+ --------
+ >> Ra = 329.88772
+ >> Dec = -56.992515
+ >> [p_ra,p_dec] = precess(Ra,Dec,1950,1975,FK4=1)
+ >> print p_ra, p_dec
+ [ 330.31442971] [-56.87186154]
+
+ Modification history
+ --------------------
+ Written by Wayne Landsman, STI Corporation, August 1986.
+ Converted to Python by Freddy R. Galindo, ROJ, 27 September 2009.
+ """
+
+ npts = numpy.size(ra)
+ ra = numpy.atleast_1d(ra)
+ dec = numpy.atleast_1d(dec)
+
+ if rad==0:
+ ra_rad = ra*Misc_Routines.CoFactors.d2r
+ dec_rad = dec*Misc_Routines.CoFactors.d2r
+ else:
+ ra_rad = ra
+ dec_rad = dec
+
+ x = numpy.zeros((npts,3))
+ x[:,0] = numpy.cos(dec_rad)*numpy.cos(ra_rad)
+ x[:,1] = numpy.cos(dec_rad)*numpy.sin(ra_rad)
+ x[:,2] = numpy.sin(dec_rad)
+
+ # Use premat function to get precession matrix from equinox1 to equinox2
+ r = self.premat(equinox1,equinox2,FK4)
+
+ x2 = numpy.dot(r,x.transpose())
+
+ ra_rad = numpy.arctan2(x2[1,:],x2[0,:])
+ dec_rad = numpy.arcsin(x2[2,:])
+
+ if rad==0:
+ ra = ra_rad/Misc_Routines.CoFactors.d2r
+ ra = ra + (ra<0)*360.
+ dec = dec_rad/Misc_Routines.CoFactors.d2r
+ else:
+ ra = ra_rad
+ ra = ra + (ra<0)*numpy.pi*2.
+ dec = dec_rad
+
+ return ra, dec
+
+ def premat(self,equinox1,equinox2,FK4=0):
+ """
+ premat returns the precession matrix needed to go from EQUINOX1 to EQUINOX2.
+
+ Parameters
+ ----------
+ equinox1 = Original equinox of coordinates, numeric scalar.
+ equinox2 = Equinox of precessed coordinates.
+ FK4 = If this keyword is set and non-zero, the FK4 (B1950) system precession angles
+ are used to compute the precession matrix. The default is to use FK5 (J2000) pre-
+ cession angles.
+
+ Return
+ ------
+ r = Precession matrix, used to precess equatorial rectangular coordinates.
+
+ Examples
+ --------
+ >> matrix = premat(1950.0,1975.0,FK4=1)
+ >> print matrix
+ [[ 9.99981438e-01 -5.58774959e-03 -2.42908517e-03]
+ [ 5.58774959e-03 9.99984388e-01 -6.78691471e-06]
+ [ 2.42908517e-03 -6.78633095e-06 9.99997050e-01]]
+
+ Modification history
+ --------------------
+ Written by Wayne Landsman, HSTX Corporation, June 1994.
+ Converted to Python by Freddy R. Galindo, ROJ, 27 September 2009.
+ """
+
+ t = 0.001*(equinox2 - equinox1)
+
+ if FK4==0:
+ st=0.001*(equinox1 - 2000.)
+ # Computing 3 rotation angles.
+ A=Misc_Routines.CoFactors.s2r*t*(23062.181+st*(139.656+0.0139*st)+t*(30.188-0.344*st+17.998*t))
+ B=Misc_Routines.CoFactors.s2r*t*t*(79.280+0.410*st+0.205*t)+A
+ C=Misc_Routines.CoFactors.s2r*t*(20043.109-st*(85.33+0.217*st)+ t*(-42.665-0.217*st-41.833*t))
+ else:
+ st=0.001*(equinox1 - 1900)
+ # Computing 3 rotation angles
+ A=Misc_Routines.CoFactors.s2r*t*(23042.53+st*(139.75+0.06*st)+t*(30.23-0.27*st+18.0*t))
+ B=Misc_Routines.CoFactors.s2r*t*t*(79.27+0.66*st+0.32*t)+A
+ C=Misc_Routines.CoFactors.s2r*t*(20046.85-st*(85.33+0.37*st)+t*(-42.67-0.37*st-41.8*t))
+
+ sina = numpy.sin(A); sinb = numpy.sin(B); sinc = numpy.sin(C)
+ cosa = numpy.cos(A); cosb = numpy.cos(B); cosc = numpy.cos(C)
+
+ r = numpy.zeros((3,3))
+ r[:,0] = numpy.array([cosa*cosb*cosc-sina*sinb,sina*cosb+cosa*sinb*cosc,cosa*sinc])
+ r[:,1] = numpy.array([-cosa*sinb-sina*cosb*cosc,cosa*cosb-sina*sinb*cosc,-sina*sinc])
+ r[:,2] = numpy.array([-cosb*sinc,-sinb*sinc,cosc])
+
+ return r
+
+ def cirrange(self,angle,rad=0):
+ """
+ cirrange forces an angle into the range 0<= angle < 360.
+
+ Parameters
+ ----------
+ angle = The angle to modify, in degrees. Can be scalar or vector.
+ rad = Set to 1 if the angle is specified in radians rather than degrees. It is for-
+ ced into the range 0 <= angle < 2 PI
+
+ Return
+ ------
+ angle = The angle after the modification.
+
+ Example
+ -------
+ >> angle = cirrange(numpy.array([420,400,361]))
+ >> print angle
+ >> [60, 40, 1]
+
+ Modification History
+ --------------------
+ Written by Michael R. Greason, Hughes STX, 10 February 1994.
+ Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009.
+ """
+
+ angle = numpy.atleast_1d(angle)
+
+ if rad==1:
+ cnst = numpy.pi*2.
+ elif rad==0:
+ cnst = 360.
+
+ # Deal with the lower limit.
+ angle = angle % cnst
+
+ # Deal with negative values, if way
+ neg = numpy.where(angle<0.0)
+ if neg[0].size>0: angle[neg] = angle[neg] + cnst
+
+ return angle
+
+
+class CelestialBodies(EquatorialCorrections):
+ def __init__(self):
+ """
+ CelestialBodies class creates a object to call methods of celestial bodies location.
+
+ Modification History
+ --------------------
+ Converted to Object-oriented Programming by Freddy Galindo, ROJ, 27 September 2009.
+ """
+
+ EquatorialCorrections.__init__(self)
+
+ def sunpos(self,jd,rad=0):
+ """
+ sunpos method computes the RA and Dec of the Sun at a given date.
+
+ Parameters
+ ----------
+ jd = The julian date of the day (and time), scalar or vector.
+ rad = If this keyword is set and non-zero, then the input and output RAD and DEC
+ vectors are in radian rather than degree.
+
+ Return
+ ------
+ ra = The right ascension of the sun at that date in degrees.
+ dec = The declination of the sun at that date in degrees.
+ elong = Ecliptic longitude of the sun at that date in degrees.
+ obliquity = The declination of the sun at that date in degrees.
+
+ Examples
+ --------
+ >> jd = 2466880
+ >> [ra,dec,elong,obliquity] = sunpos(jd)
+ >> print ra, dec, elong, obliquity
+ [ 275.53499556] [-23.33840558] [ 275.08917968] [ 23.43596165]
+
+ >> [ra,dec,elong,obliquity] = sunpos(jd,rad=1)
+ >> print ra, dec, elong, obliquity
+ [ 4.80899288] [-0.40733202] [ 4.80121192] [ 0.40903469]
+
+ >> jd = 2450449.5 + numpy.arange(365)
+ >> [ra,dec,elong,obliquity] = sunpos(jd)
+
+ Modification history
+ --------------------
+ Written by Micheal R. Greason, STX Corporation, 28 October 1988.
+ Converted to Python by Freddy R. Galindo, ROJ, 27 September 2009.
+ """
+
+ jd = numpy.atleast_1d(jd)
+
+ # Form time in Julian centuries from 1900.
+ t = (jd -2415020.0)/36525.0
+
+ # Form sun's mean longitude
+ l = (279.696678+((36000.768925*t) % 360.0))*3600.0
+
+ # Allow for ellipticity of the orbit (equation of centre) using the Earth's mean
+ # anomoly ME
+ me = 358.475844 + ((35999.049750*t) % 360.0)
+ ellcor = (6910.1 - 17.2*t)*numpy.sin(me*Misc_Routines.CoFactors.d2r) + 72.3*numpy.sin(2.0*me*Misc_Routines.CoFactors.d2r)
+ l = l + ellcor
+
+ # Allow for the Venus perturbations using the mean anomaly of Venus MV
+ mv = 212.603219 + ((58517.803875*t) % 360.0)
+ vencorr = 4.8*numpy.cos((299.1017 + mv - me)*Misc_Routines.CoFactors.d2r) + \
+ 5.5*numpy.cos((148.3133 + 2.0*mv - 2.0*me )*Misc_Routines.CoFactors.d2r) + \
+ 2.5*numpy.cos((315.9433 + 2.0*mv - 3.0*me )*Misc_Routines.CoFactors.d2r) + \
+ 1.6*numpy.cos((345.2533 + 3.0*mv - 4.0*me )*Misc_Routines.CoFactors.d2r) + \
+ 1.0*numpy.cos((318.15 + 3.0*mv - 5.0*me )*Misc_Routines.CoFactors.d2r)
+ l = l + vencorr
+
+ # Allow for the Mars perturbations using the mean anomaly of Mars MM
+ mm = 319.529425 + ((19139.858500*t) % 360.0)
+ marscorr = 2.0*numpy.cos((343.8883 - 2.0*mm + 2.0*me)*Misc_Routines.CoFactors.d2r ) + \
+ 1.8*numpy.cos((200.4017 - 2.0*mm + me)*Misc_Routines.CoFactors.d2r)
+ l = l + marscorr
+
+ # Allow for the Jupiter perturbations using the mean anomaly of Jupiter MJ
+ mj = 225.328328 + ((3034.6920239*t) % 360.0)
+ jupcorr = 7.2*numpy.cos((179.5317 - mj + me )*Misc_Routines.CoFactors.d2r) + \
+ 2.6*numpy.cos((263.2167 - mj)*Misc_Routines.CoFactors.d2r) + \
+ 2.7*numpy.cos((87.1450 - 2.0*mj + 2.0*me)*Misc_Routines.CoFactors.d2r) + \
+ 1.6*numpy.cos((109.4933 - 2.0*mj + me)*Misc_Routines.CoFactors.d2r)
+ l = l + jupcorr
+
+ # Allow for Moons perturbations using mean elongation of the Moon from the Sun D
+ d = 350.7376814 + ((445267.11422*t) % 360.0)
+ mooncorr = 6.5*numpy.sin(d*Misc_Routines.CoFactors.d2r)
+ l = l + mooncorr
+
+ # Allow for long period terms
+ longterm = + 6.4*numpy.sin((231.19 + 20.20*t)*Misc_Routines.CoFactors.d2r)
+ l = l + longterm
+ l = (l + 2592000.0) % 1296000.0
+ longmed = l/3600.0
+
+ # Allow for Aberration
+ l = l - 20.5
+
+ # Allow for Nutation using the longitude of the Moons mean node OMEGA
+ omega = 259.183275 - ((1934.142008*t) % 360.0)
+ l = l - 17.2*numpy.sin(omega*Misc_Routines.CoFactors.d2r)
+
+ # Form the True Obliquity
+ oblt = 23.452294 - 0.0130125*t + (9.2*numpy.cos(omega*Misc_Routines.CoFactors.d2r))/3600.0
+
+ # Form Right Ascension and Declination
+ l = l/3600.0
+ ra = numpy.arctan2((numpy.sin(l*Misc_Routines.CoFactors.d2r)*numpy.cos(oblt*Misc_Routines.CoFactors.d2r)),numpy.cos(l*Misc_Routines.CoFactors.d2r))
+
+ neg = numpy.where(ra < 0.0)
+ if neg[0].size > 0: ra[neg] = ra[neg] + 2.0*numpy.pi
+
+ dec = numpy.arcsin(numpy.sin(l*Misc_Routines.CoFactors.d2r)*numpy.sin(oblt*Misc_Routines.CoFactors.d2r))
+
+ if rad==1:
+ oblt = oblt*Misc_Routines.CoFactors.d2r
+ longmed = longmed*Misc_Routines.CoFactors.d2r
+ else:
+ ra = ra/Misc_Routines.CoFactors.d2r
+ dec = dec/Misc_Routines.CoFactors.d2r
+
+ return ra, dec, longmed, oblt
+
+ def moonpos(self,jd,rad=0):
+ """
+ moonpos method computes the RA and Dec of the Moon at specified Julian date(s).
+
+ Parameters
+ ----------
+ jd = The julian date of the day (and time), scalar or vector.
+ rad = If this keyword is set and non-zero, then the input and output RAD and DEC
+ vectors are in radian rather than degree.
+
+ Return
+ ------
+ ra = The right ascension of the sun at that date in degrees.
+ dec = The declination of the sun at that date in degrees.
+ dist = The Earth-moon distance in kilometers (between the center of the Earth and
+ the center of the moon).
+ geolon = Apparent longitude of the moon in degrees, referred to the ecliptic of the
+ specified date(s).
+ geolat = Apparent latitude the moon in degrees, referred to the ecliptic of the
+ specified date(s).
+
+ Examples
+ --------
+ >> jd = 2448724.5
+ >> [ra,dec,dist,geolon,geolat] = sunpos(jd)
+ >> print ra, dec, dist, geolon, geolat
+ [ 134.68846855] [ 13.76836663] [ 368409.68481613] [ 133.16726428] [-3.22912642]
+
+ >> [ra,dec,dist,geolon, geolat] = sunpos(jd,rad=1)
+ >> print ra, dec, dist, geolon, geolat
+ [ 2.35075724] [ 0.24030333] [ 368409.68481613] [ 2.32420722] [-0.05635889]
+
+ >> jd = 2450449.5 + numpy.arange(365)
+ >> [ra,dec,dist,geolon, geolat] = sunpos(jd)
+
+ Modification history
+ --------------------
+ Written by Micheal R. Greason, STX Corporation, 31 October 1988.
+ Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009.
+ """
+
+ jd = numpy.atleast_1d(jd)
+
+ # Form time in Julian centuries from 1900.
+ t = (jd - 2451545.0)/36525.0
+
+ 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,\
+ 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])
+
+ 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,\
+ 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])
+
+ 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,\
+ -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])
+
+ 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,\
+ 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])
+
+ sin_lng = numpy.array([6288774,1274027,658314,213618,-185116,-114332,58793,57066,\
+ 53322,45758,-40923,-34720,-30383,15327,-12528,10980,10675,10034,8548,-7888,\
+ -6766,-5163,4987,4036,3994,3861,3665,-2689,-2602,2390,-2348,2236,-2120,-2069,\
+ 2048,-1773,-1595,1215,-1110,-892,-810,759,-713,-700,691,596,549,537,520,-487,\
+ -399,-381,351,-340,330,327,-323,299,294,0.0])
+
+ cos_lng = numpy.array([-20905355,-3699111,-2955968,-569925,48888,-3149,246158,-152138,\
+ -170733,-204586,-129620,108743,104755,10321,0,79661,-34782,-23210,-21636,24208,\
+ 30824,-8379,-16675,-12831,-10445,-11650,14403,-7003,0,10056,6322, -9884,5751,0,\
+ -4950,4130,0,-3958,0,3258,2616,-1897,-2117,2354,0,0,-1423,-1117,-1571,-1739,0, \
+ -4421,0,0,0,0,1165,0,0,8752.0])
+
+ 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,\
+ 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])
+
+ 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,\
+ 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])
+
+ 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,\
+ 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])
+
+ 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,\
+ -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, \
+ -1,-1,-1,-1,1])
+
+ sin_lat = numpy.array([5128122,280602,277693,173237,55413,46271, 32573, 17198, 9266, \
+ 8822,8216,4324,4200,-3359,2463,2211,2065,-1870,1828,-1794, -1749, -1565, -1491, \
+ -1475,-1410,-1344,-1335,1107,1021,833,777,671,607,596,491,-451,439,422,421,-366,\
+ -351,331,315,302,-283,-229,223,223,-220,-220,-185,181,-177,176, 166, -164, 132, \
+ -119,115,107.0])
+
+ # Mean longitude of the moon refered to mean equinox of the date.
+ coeff0 = numpy.array([-1./6.5194e7,1./538841.,-0.0015786,481267.88123421,218.3164477])
+ lprimed = numpy.poly1d(coeff0)
+ lprimed = lprimed(t)
+ lprimed = self.cirrange(lprimed,rad=0)
+ lprime = lprimed*Misc_Routines.CoFactors.d2r
+
+ # Mean elongation of the moon
+ coeff1 = numpy.array([-1./1.13065e8,1./545868.,-0.0018819,445267.1114034,297.8501921])
+ d = numpy.poly1d(coeff1)
+ d = d(t)*Misc_Routines.CoFactors.d2r
+ d = self.cirrange(d,rad=1)
+
+ # Sun's mean anomaly
+ coeff2 = numpy.array([1.0/2.449e7,-0.0001536,35999.0502909,357.5291092])
+ M = numpy.poly1d(coeff2)
+ M = M(t)*Misc_Routines.CoFactors.d2r
+ M = self.cirrange(M,rad=1)
+
+ # Moon's mean anomaly
+ coeff3 = numpy.array([-1.0/1.4712e7,1.0/6.9699e4,0.0087414,477198.8675055,134.9633964])
+ Mprime = numpy.poly1d(coeff3)
+ Mprime = Mprime(t)*Misc_Routines.CoFactors.d2r
+ Mprime = self.cirrange(Mprime,rad=1)
+
+ # Moon's argument of latitude
+ coeff4 = numpy.array([1.0/8.6331e8,-1.0/3.526e7,-0.0036539,483202.0175233,93.2720950])
+ F = numpy.poly1d(coeff4)
+ F = F(t)*Misc_Routines.CoFactors.d2r
+ F = self.cirrange(F,rad=1)
+
+ # Eccentricity of Earth's orbit around the sun
+ e = 1 - 0.002516*t - 7.4e-6*(t**2.)
+ e2 = e**2.
+
+ ecorr1 = numpy.where((numpy.abs(m_lng))==1)
+ ecorr2 = numpy.where((numpy.abs(m_lat))==1)
+ ecorr3 = numpy.where((numpy.abs(m_lng))==2)
+ ecorr4 = numpy.where((numpy.abs(m_lat))==2)
+
+ # Additional arguments.
+ A1 = (119.75 + 131.849*t)*Misc_Routines.CoFactors.d2r
+ A2 = (53.09 + 479264.290*t)*Misc_Routines.CoFactors.d2r
+ A3 = (313.45 + 481266.484*t)*Misc_Routines.CoFactors.d2r
+ suml_add = 3958.*numpy.sin(A1) + 1962.*numpy.sin(lprime - F) + 318*numpy.sin(A2)
+ sumb_add = -2235.*numpy.sin(lprime) + 382.*numpy.sin(A3) + 175.*numpy.sin(A1-F) + \
+ 175.*numpy.sin(A1 + F) + 127.*numpy.sin(lprime - Mprime) - 115.*numpy.sin(lprime + Mprime)
+
+ # Sum the periodic terms
+ geolon = numpy.zeros(jd.size)
+ geolat = numpy.zeros(jd.size)
+ dist = numpy.zeros(jd.size)
+
+ for i in numpy.arange(jd.size):
+ sinlng = sin_lng
+ coslng = cos_lng
+ sinlat = sin_lat
+
+ sinlng[ecorr1] = e[i]*sinlng[ecorr1]
+ coslng[ecorr1] = e[i]*coslng[ecorr1]
+ sinlat[ecorr2] = e[i]*sinlat[ecorr2]
+ sinlng[ecorr3] = e2[i]*sinlng[ecorr3]
+ coslng[ecorr3] = e2[i]*coslng[ecorr3]
+ sinlat[ecorr4] = e2[i]*sinlat[ecorr4]
+
+ arg = d_lng*d[i] + m_lng*M[i] + mp_lng*Mprime[i] + f_lng*F[i]
+ geolon[i] = lprimed[i] + (numpy.sum(sinlng*numpy.sin(arg)) + suml_add[i] )/1.e6
+ dist[i] = 385000.56 + numpy.sum(coslng*numpy.cos(arg))/1.e3
+ arg = d_lat*d[i] + m_lat*M[i] + mp_lat*Mprime[i] + f_lat*F[i]
+ geolat[i] = (numpy.sum(sinlat*numpy.sin(arg)) + sumb_add[i])/1.e6
+
+ [nlon, elon] = self.nutate(jd)
+ geolon = geolon + nlon/3.6e3
+ geolon = self.cirrange(geolon,rad=0)
+ lamb = geolon*Misc_Routines.CoFactors.d2r
+ beta = geolat*Misc_Routines.CoFactors.d2r
+
+ # Find mean obliquity and convert lamb, beta to RA, Dec
+ c = numpy.array([2.45,5.79,27.87,7.12,-39.05,-249.67,-51.38,1999.25,-1.55,-4680.93, \
+ 21.448])
+ junk = numpy.poly1d(c);
+ epsilon = 23. + (26./60.) + (junk(t/1.e2)/3600.)
+ # True obliquity in radians
+ eps = (epsilon + elon/3600. )*Misc_Routines.CoFactors.d2r
+
+ ra = numpy.arctan2(numpy.sin(lamb)*numpy.cos(eps)-numpy.tan(beta)*numpy.sin(eps),numpy.cos(lamb))
+ ra = self.cirrange(ra,rad=1)
+
+ dec = numpy.arcsin(numpy.sin(beta)*numpy.cos(eps) + numpy.cos(beta)*numpy.sin(eps)*numpy.sin(lamb))
+
+ if rad==1:
+ geolon = lamb
+ geolat = beta
+ else:
+ ra = ra/Misc_Routines.CoFactors.d2r
+ dec = dec/Misc_Routines.CoFactors.d2r
+
+ return ra, dec, dist, geolon, geolat
+
+ def hydrapos(self):
+ """
+ hydrapos method returns RA and Dec provided by Bill Coles (Oct 2003).
+
+ Parameters
+ ----------
+ None
+
+ Return
+ ------
+ ra = The right ascension of the sun at that date in degrees.
+ dec = The declination of the sun at that date in degrees.
+ Examples
+ --------
+ >> [ra,dec] = hydrapos()
+ >> print ra, dec
+ 139.45 -12.0833333333
+
+ Modification history
+ --------------------
+ Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009.
+ """
+
+ ra = (9. + 17.8/60.)*15.
+ dec = -(12. + 5./60.)
+
+ return ra, dec
+
+
+ def skynoise_jro(self,dec_cut=-11.95,filename='skynoise_jro.dat',filepath=None):
+ """
+ hydrapos returns RA and Dec provided by Bill Coles (Oct 2003).
+
+ Parameters
+ ----------
+ dec_cut = A scalar giving the declination to get a cut of the skynoise over Jica-
+ marca. The default value is -11.95.
+ filename = A string to specify name the skynoise file. The default value is skynoi-
+ se_jro.dat
+
+ Return
+ ------
+ maxra = The maximum right ascension to the declination used to get a cut.
+ ra = The right ascension.
+ Examples
+ --------
+ >> [maxra,ra] = skynoise_jro()
+ >> print maxra, ra
+ 139.45 -12.0833333333
+
+ Modification history
+ --------------------
+ Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009.
+ """
+
+ if filepath==None:filepath = './resource'
+
+ f = open(os.path.join(filepath,filename),'rb')
+
+ # Reading SkyNoise Power (lineal scale)
+ ha_sky = numpy.fromfile(f,numpy.dtype([('var','> [maxra,ra] = skynoise_jro()
+ >> print maxra, ra
+ 139.45 -12.0833333333
+
+ Modification history
+ --------------------
+ Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009.
+ """
+
+ # Defining date to compute SkyNoise.
+ [year, month, dom, hour, mis, secs] = TimeTools.Julian(jd).change2time()
+ is_dom = (month==9) & (dom==21)
+ if is_dom:
+ tmp = jd
+ jd = TimeTools.Time(year,9,22).change2julian()
+ dom = 22
+
+ # Reading SkyNoise
+ if filepath==None:filepath='./resource'
+ f = open(os.path.join(filepath,filename))
+
+ lines = f.read()
+ f.close()
+
+ nlines = 99
+ lines = lines.split('\n')
+ data = numpy.zeros((2,nlines))*numpy.float32(0.)
+ for ii in numpy.arange(nlines):
+ line = numpy.array([lines[ii][0:6],lines[ii][6:]])
+ data[:,ii] = numpy.float32(line)
+
+ # Getting SkyNoise to the date desired.
+ otime = data[0,:]*60.0
+ opowr = data[1,:]
+
+ hour = numpy.array([0,23]);
+ mins = numpy.array([0,59]);
+ secs = numpy.array([0,59]);
+ LTrange = TimeTools.Time(year,month,dom,hour,mins,secs).change2julday()
+ LTtime = LTrange[0] + numpy.arange(1440)*((LTrange[1] - LTrange[0])/(1440.-1))
+ lst = TimeTools.Julian(LTtime + (-3600.*ut/86400.)).change2lst()
+
+ ipowr = lst*0.0
+ # Interpolating using scipy (inside max and min "x")
+ otime = otime/3600.
+ val = numpy.where((lst>numpy.min(otime)) & (lstnumpy.max(otime))
+ if uval[0].size>0:
+ ii = numpy.min(uval[0])
+ m = (ipowr[ii-1] - ipowr[ii-2])/(lst[ii-1] - lst[ii-2])
+ b = ipowr[ii-1] - m*lst[ii-1]
+ ipowr[uval] = m*lst[uval] + b
+
+ if is_dom:
+ lst = numpy.roll(lst,4)
+ ipowr = numpy.roll(ipowr,4)
+
+ new_lst = numpy.int32(lst*3600.)
+ new_pow = ipowr
+
+ return ipowr, LTtime, lst
+
+
+class AltAz(EquatorialCorrections):
+ def __init__(self,alt,az,jd,lat=-11.95,lon=-76.8667,WS=0,altitude=500,nutate_=0,precess_=0,\
+ aberration_=0,B1950=0):
+ """
+ The AltAz class creates an object which represents the target position in horizontal
+ coordinates (alt-az) and allows to convert (using the methods) from this coordinate
+ system to others (e.g. Equatorial).
+
+ Parameters
+ ----------
+ alt = Altitude in degrees. Scalar or vector.
+ az = Azimuth angle in degrees (measured EAST from NORTH, but see keyword WS). Sca-
+ lar or vector.
+ jd = Julian date. Scalar or vector.
+ lat = North geodetic latitude of location in degrees. The default value is -11.95.
+ lon = East longitude of location in degrees. The default value is -76.8667.
+ WS = Set this to 1 to get the azimuth measured westward from south.
+ altitude = The altitude of the observing location, in meters. The default 500.
+ nutate_ = Set this to 1 to force nutation, 0 for no nutation.
+ precess_ = Set this to 1 to force precession, 0 for no precession.
+ aberration_ = Set this to 1 to force aberration correction, 0 for no correction.
+ B1950 = Set this if your RA and DEC are specified in B1950, FK4 coordinates (ins-
+ tead of J2000, FK5)
+
+ Modification History
+ --------------------
+ Converted to Object-oriented Programming by Freddy Galindo, ROJ, 26 September 2009.
+ """
+
+ EquatorialCorrections.__init__(self)
+
+ self.alt = numpy.atleast_1d(alt)
+ self.az = numpy.atleast_1d(az)
+ self.jd = numpy.atleast_1d(jd)
+ self.lat = lat
+ self.lon = lon
+ self.WS = WS
+ self.altitude = altitude
+
+ self.nutate_ = nutate_
+ self.aberration_ = aberration_
+ self.precess_ = precess_
+ self.B1950 = B1950
+
+ def change2equatorial(self):
+ """
+ change2equatorial method converts horizon (Alt-Az) coordinates to equatorial coordi-
+ nates (ra-dec).
+
+ Return
+ ------
+ ra = Right ascension of object (J2000) in degrees (FK5). Scalar or vector.
+ dec = Declination of object (J2000), in degrees (FK5). Scalar or vector.
+ ha = Hour angle in degrees.
+
+ Example
+ -------
+ >> alt = 88.5401
+ >> az = -128.990
+ >> jd = 2452640.5
+ >> ObjAltAz = AltAz(alt,az,jd)
+ >> [ra, dec, ha] = ObjAltAz.change2equatorial()
+ >> print ra, dec, ha
+ [ 22.20280632] [-12.86610025] [ 1.1638927]
+
+ Modification History
+ --------------------
+ Written Chris O'Dell Univ. of Wisconsin-Madison, May 2002.
+ Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009.
+ """
+
+ az = self.az
+ alt = self.alt
+ if self.WS>0:az = az -180.
+ ra_tmp = numpy.zeros(numpy.size(self.jd)) + 45.
+ dec_tmp = numpy.zeros(numpy.size(self.jd)) + 45.
+ [dra1,ddec1,eps,d_psi,d_eps] = self.co_nutate(self.jd,ra_tmp, dec_tmp)
+
+ # Getting local mean sidereal time (lmst)
+ lmst = TimeTools.Julian(self.jd[0]).change2lst()
+ lmst = lmst*Misc_Routines.CoFactors.h2d
+ # Getting local apparent sidereal time (last)
+ last = lmst + d_psi*numpy.cos(eps)/3600.
+
+ # Now do the spherical trig to get APPARENT hour angle and declination (Degrees).
+ [ha, dec] = self.change2HaDec()
+
+ # Finding Right Ascension (in degrees, from 0 to 360.)
+ ra = (last - ha + 360.) % 360.
+
+ # Calculate NUTATION and ABERRATION Correction to Ra-Dec
+ [dra1, ddec1,eps,d_psi,d_eps] = self.co_nutate(self.jd,ra,dec)
+ [dra2,ddec2,eps] = self.co_aberration(self.jd,ra,dec)
+
+ # Make Nutation and Aberration correction (if wanted)
+ ra = ra - (dra1*self.nutate_ + dra2*self.aberration_)/3600.
+ dec = dec - (ddec1*self.nutate_ + ddec2*self.aberration_)/3600.
+
+ # Computing current equinox
+ j_now = (self.jd - 2451545.)/365.25 + 2000
+
+ # Precess coordinates to current date
+ if self.precess_==1:
+ njd = numpy.size(self.jd)
+ for ii in numpy.arange(njd):
+ ra_i = ra[ii]
+ dec_i = dec[ii]
+ now = j_now[ii]
+
+ if self.B1950==1:
+ [ra_i,dec_i] = self.precess(ra_i,dec_i,now,1950.,FK4=1)
+ elif self.B1950==0:
+ [ra_i,dec_i] = self.precess(ra_i,dec_i,now,2000.,FK4=0)
+
+ ra[ii] = ra_i
+ dec[ii] = dec_i
+
+ return ra, dec, ha
+
+ def change2HaDec(self):
+ """
+ change2HaDec method converts from horizon (Alt-Az) coordinates to hour angle and de-
+ clination.
+
+ Return
+ ------
+ ha = The local apparent hour angle, in degrees. The hour angle is the time that ri-
+ ght ascension of 0 hours crosses the local meridian. It is unambiguisoly defined.
+ dec = The local apparent declination, in degrees.
+
+ Example
+ -------
+ >> alt = 88.5401
+ >> az = -128.990
+ >> jd = 2452640.5
+ >> ObjAltAz = AltAz(alt,az,jd)
+ >> [ha, dec] = ObjAltAz.change2HaDec()
+ >> print ha, dec
+ [ 1.1638927] [-12.86610025]
+
+ Modification History
+ --------------------
+ Written Chris O'Dell Univ. of Wisconsin-Madison, May 2002.
+ Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009.
+ """
+
+ alt_r = numpy.atleast_1d(self.alt*Misc_Routines.CoFactors.d2r)
+ az_r = numpy.atleast_1d(self.az*Misc_Routines.CoFactors.d2r)
+ lat_r = numpy.atleast_1d(self.lat*Misc_Routines.CoFactors.d2r)
+
+ # Find local hour angle (in degrees, from 0 to 360.)
+ y_ha = -1*numpy.sin(az_r)*numpy.cos(alt_r)
+ x_ha = -1*numpy.cos(az_r)*numpy.sin(lat_r)*numpy.cos(alt_r) + numpy.sin(alt_r)*numpy.cos(lat_r)
+
+ ha = numpy.arctan2(y_ha,x_ha)
+ ha = ha/Misc_Routines.CoFactors.d2r
+
+ w = numpy.where(ha<0.)
+ if w[0].size>0:ha[w] = ha[w] + 360.
+ ha = ha % 360.
+
+ # Find declination (positive if north of celestial equatorial, negative if south)
+ sindec = numpy.sin(lat_r)*numpy.sin(alt_r) + numpy.cos(lat_r)*numpy.cos(alt_r)*numpy.cos(az_r)
+ dec = numpy.arcsin(sindec)/Misc_Routines.CoFactors.d2r
+
+ return ha, dec
+
+
+class Equatorial(EquatorialCorrections):
+ def __init__(self,ra,dec,jd,lat=-11.95,lon=-76.8667,WS=0,altitude=500,nutate_=0,precess_=0,\
+ aberration_=0,B1950=0):
+ """
+ The Equatorial class creates an object which represents the target position in equa-
+ torial coordinates (ha-dec) and allows to convert (using the class methods) from
+ this coordinate system to others (e.g. AltAz).
+
+ Parameters
+ ----------
+ ra = Right ascension of object (J2000) in degrees (FK5). Scalar or vector.
+ dec = Declination of object (J2000), in degrees (FK5). Scalar or vector.
+ jd = Julian date. Scalar or vector.
+ lat = North geodetic latitude of location in degrees. The default value is -11.95.
+ lon = East longitude of location in degrees. The default value is -76.8667.
+ WS = Set this to 1 to get the azimuth measured westward from south.
+ altitude = The altitude of the observing location, in meters. The default 500.
+ nutate = Set this to 1 to force nutation, 0 for no nutation.
+ precess = Set this to 1 to force precession, 0 for no precession.
+ aberration = Set this to 1 to force aberration correction, 0 for no correction.
+ B1950 = Set this if your RA and DEC are specified in B1950, FK4 coordinates (ins-
+ tead of J2000, FK5)
+
+ Modification History
+ --------------------
+ Converted to Object-oriented Programming by Freddy Galindo, ROJ, 29 September 2009.
+ """
+
+ EquatorialCorrections.__init__(self)
+
+ self.ra = numpy.atleast_1d(ra)
+ self.dec = numpy.atleast_1d(dec)
+ self.jd = numpy.atleast_1d(jd)
+ self.lat = lat
+ self.lon = lon
+ self.WS = WS
+ self.altitude = altitude
+
+ self.nutate_ = nutate_
+ self.aberration_ = aberration_
+ self.precess_ = precess_
+ self.B1950 = B1950
+
+ def change2AltAz(self):
+ """
+ change2AltAz method converts from equatorial coordinates (ha-dec) to horizon coordi-
+ nates (alt-az).
+
+ Return
+ ------
+ alt = Altitude in degrees. Scalar or vector.
+ az = Azimuth angle in degrees (measured EAST from NORTH, but see keyword WS). Sca-
+ lar or vector.
+ ha = Hour angle in degrees.
+
+ Example
+ -------
+ >> ra = 43.370609
+ >> dec = -28.0000
+ >> jd = 2452640.5
+ >> ObjEq = Equatorial(ra,dec,jd)
+ >> [alt, az, ha] = ObjEq.change2AltAz()
+ >> print alt, az, ha
+ [ 65.3546497] [ 133.58753124] [ 339.99609002]
+
+ Modification History
+ --------------------
+ Written Chris O'Dell Univ. of Wisconsin-Madison. May 2002
+ Converted to Python by Freddy R. Galindo, ROJ, 29 September 2009.
+ """
+
+ ra = self.ra
+ dec = self.dec
+
+ # Computing current equinox
+ j_now = (self.jd - 2451545.)/365.25 + 2000
+
+ # Precess coordinates to current date
+ if self.precess_==1:
+ njd = numpy.size(self.jd)
+ for ii in numpy.arange(njd):
+ ra_i = ra[ii]
+ dec_i = dec[ii]
+ now = j_now[ii]
+
+ if self.B1950==1:
+ [ra_i,dec_i] = self.precess(ra_i,dec_i,now,1950.,FK4=1)
+ elif self.B1950==0:
+ [ra_i,dec_i] = self.precess(ra_i,dec_i,now,2000.,FK4=0)
+
+ ra[ii] = ra_i
+ dec[ii] = dec_i
+
+ # Calculate NUTATION and ABERRATION Correction to Ra-Dec
+ [dra1, ddec1,eps,d_psi,d_eps] = self.co_nutate(self.jd,ra,dec)
+ [dra2,ddec2,eps] = self.co_aberration(self.jd,ra,dec)
+
+ # Make Nutation and Aberration correction (if wanted)
+ ra = ra + (dra1*self.nutate_ + dra2*self.aberration_)/3600.
+ dec = dec + (ddec1*self.nutate_ + ddec2*self.aberration_)/3600.
+
+ # Getting local mean sidereal time (lmst)
+ lmst = TimeTools.Julian(self.jd).change2lst()
+
+ lmst = lmst*Misc_Routines.CoFactors.h2d
+ # Getting local apparent sidereal time (last)
+ last = lmst + d_psi*numpy.cos(eps)/3600.
+
+ # Finding Hour Angle (in degrees, from 0 to 360.)
+ ha = last - ra
+ w = numpy.where(ha<0.)
+ if w[0].size>0:ha[w] = ha[w] + 360.
+ ha = ha % 360.
+
+ # Now do the spherical trig to get APPARENT hour angle and declination (Degrees).
+ [alt, az] = self.HaDec2AltAz(ha,dec)
+
+ return alt, az, ha
+
+ def HaDec2AltAz(self,ha,dec):
+ """
+ HaDec2AltAz convert hour angle and declination (ha-dec) to horizon coords (alt-az).
+
+ Parameters
+ ----------
+ ha = The local apparent hour angle, in DEGREES, scalar or vector.
+ dec = The local apparent declination, in DEGREES, scalar or vector.
+
+ Return
+ ------
+ alt = Altitude in degrees. Scalar or vector.
+ az = Azimuth angle in degrees (measured EAST from NORTH, but see keyword WS). Sca-
+ lar or vector.
+
+ Modification History
+ --------------------
+ Written Chris O'Dell Univ. of Wisconsin-Madison, May 2002.
+ Converted to Python by Freddy R. Galindo, ROJ, 26 September 2009.
+ """
+
+ sh = numpy.sin(ha*Misc_Routines.CoFactors.d2r) ; ch = numpy.cos(ha*Misc_Routines.CoFactors.d2r)
+ sd = numpy.sin(dec*Misc_Routines.CoFactors.d2r) ; cd = numpy.cos(dec*Misc_Routines.CoFactors.d2r)
+ sl = numpy.sin(self.lat*Misc_Routines.CoFactors.d2r) ; cl = numpy.cos(self.lat*Misc_Routines.CoFactors.d2r)
+
+ x = -1*ch*cd*sl + sd*cl
+ y = -1*sh*cd
+ z = ch*cd*cl + sd*sl
+ r = numpy.sqrt(x**2. + y**2.)
+
+ az = numpy.arctan2(y,x)/Misc_Routines.CoFactors.d2r
+ alt = numpy.arctan2(z,r)/Misc_Routines.CoFactors.d2r
+
+ # correct for negative az.
+ w = numpy.where(az<0.)
+ if w[0].size>0:az[w] = az[w] + 360.
+
+ # Convert az to West from South, if desired
+ if self.WS==1: az = (az + 180.) % 360.
+
+ return alt, az
+
+
+class Geodetic():
+ def __init__(self,lat=-11.95,alt=0):
+ """
+ The Geodetic class creates an object which represents the real position on earth of
+ a target (Geodetic Coordinates: lat-alt) and allows to convert (using the class me-
+ thods) from this coordinate system to others (e.g. geocentric).
+
+ Parameters
+ ----------
+ lat = Geodetic latitude of location in degrees. The default value is -11.95.
+
+ alt = Geodetic altitude (km). The default value is 0.
+
+ Modification History
+ --------------------
+ Converted to Object-oriented Programming by Freddy R. Galindo, ROJ, 02 October 2009.
+ """
+
+ self.lat = numpy.atleast_1d(lat)
+ self.alt = numpy.atleast_1d(alt)
+
+ self.a = 6378.16
+ self.ab2 = 1.0067397
+ self.ep2 = 0.0067397
+
+ def change2geocentric(self):
+ """
+ change2geocentric method converts from Geodetic to Geocentric coordinates. The re-
+ ference geoid is that adopted by the IAU in 1964.
+
+ Return
+ ------
+ gclat = Geocentric latitude (in degrees), scalar or vector.
+ gcalt = Geocentric radial distance (km), scalar or vector.
+
+ Example
+ -------
+ >> ObjGeoid = Geodetic(lat=-11.95,alt=0)
+ >> [gclat, gcalt] = ObjGeoid.change2geocentric()
+ >> print gclat, gcalt
+ [-11.87227742] [ 6377.25048195]
+
+ Modification History
+ --------------------
+ Converted to Python by Freddy R. Galindo, ROJ, 02 October 2009.
+ """
+
+ gdl = self.lat*Misc_Routines.CoFactors.d2r
+ slat = numpy.sin(gdl)
+ clat = numpy.cos(gdl)
+ slat2 = slat**2.
+ clat2 = (self.ab2*clat)**2.
+
+ sbet = slat/numpy.sqrt(slat2 + clat2)
+ sbet2 = (sbet**2.) # < 1
+ noval = numpy.where(sbet2>1)
+ if noval[0].size>0:sbet2[noval] = 1
+ cbet = numpy.sqrt(1. - sbet2)
+
+ rgeoid = self.a/numpy.sqrt(1. + self.ep2*sbet2)
+
+ x = rgeoid*cbet + self.alt*clat
+ y = rgeoid*sbet + self.alt*slat
+
+ gcalt = numpy.sqrt(x**2. + y**2.)
+ gclat = numpy.arctan2(y,x)/Misc_Routines.CoFactors.d2r
+
+ return gclat, gcalt
diff --git a/apps/abs/utils/Files.py b/apps/abs/utils/Files.py
new file mode 100644
index 0000000..f392d95
--- /dev/null
+++ b/apps/abs/utils/Files.py
@@ -0,0 +1,18 @@
+'''
+Created on Jun 19, 2013
+
+@author: Jose Antonio Sal y Rosas Celi
+@contact: jose.salyrosas@jro.igp.gob.pe
+'''
+
+from datetime import datetime
+
+class Files(object):
+
+ def setFilename(self):
+ return datetime.today().strftime("%Y%m%d%H%M%S%f")
+
+ def save(self, filename, contentFile):
+ f = open(filename, 'a+')
+ f.write(contentFile)
+ f.close()
diff --git a/apps/abs/utils/Graphics_Miscens.py b/apps/abs/utils/Graphics_Miscens.py
new file mode 100644
index 0000000..ae73466
--- /dev/null
+++ b/apps/abs/utils/Graphics_Miscens.py
@@ -0,0 +1,51 @@
+"""
+The GRAPHICS_MISC.py module gathers classes and/or functions useful for generation of plots.
+
+MODULES CALLED:
+NUMPY, OS
+
+MODIFICATION HISTORY:
+Created by Ing. Freddy Galindo (frederickgalindo@gmail.com). ROJ, 13 August 2009.
+"""
+
+import os
+import numpy
+import sys
+
+
+class ColorTable:
+ def __init__(self,table=1,filepath=None):
+ self.table = table
+ #set to path for data folder, file: col_koki.dat
+ if filepath==None:
+ filepath= './data/'
+ self.filepath = filepath
+
+ def readTable(self):
+ if self.table>0:
+ if self.table==1:
+ f = open(os.path.join(self.filepath,'col_koki.dat'),'rb')
+
+ #f = open('./col_koki.dat','rb')
+
+ # Reading SkyNoise Power (lineal scale)
+ blue = numpy.fromfile(f,numpy.dtype([('var','b')]),256)
+ blue = numpy.int32(blue['var'])
+ val = numpy.where(blue<0)
+ if val[0].size:blue[val] = blue[val] + numpy.int32(256)
+
+ green = numpy.fromfile(f,numpy.dtype([('var','b')]),256)
+ green = numpy.int32(green['var'])
+ val = numpy.where(green<0)
+ if val[0].size:green[val] = green[val] + numpy.int32(256)
+
+ red = numpy.fromfile(f,numpy.dtype([('var','b')]),256)
+ red = numpy.int32(red['var'])
+ val = numpy.where(red<0)
+ if val[0].size:red[val] = red[val] + numpy.int32(256)
+
+ f.close()
+
+ colortable = numpy.array([red/255.,green/255.,blue/255.])
+
+ return colortable
diff --git a/apps/abs/utils/Graphics_OverJro.py b/apps/abs/utils/Graphics_OverJro.py
new file mode 100644
index 0000000..6012cd2
--- /dev/null
+++ b/apps/abs/utils/Graphics_OverJro.py
@@ -0,0 +1,563 @@
+"""
+The module GRAPHICS_OVERJRO.py gathers classes or/and functions to create graphics from OVER-JRO
+project (e.g. antenna patterns, skynoise, ...).
+
+MODULES CALLED:
+TIME, NUMPY, MATPLOTLIB, TIMETOOLS
+
+MODIFICATION HISTORY:
+Created by Ing. Freddy Galindo (frederickgalindo@gmail.com). ROJ Oct 18, 2009.
+"""
+
+import time
+import numpy
+import sys
+import os
+
+# set HOME environment variable to a directory the httpd server can write to
+#os.environ[ 'HOME' ] = '/usr/local/www/htdocs/overJro/tempReports'
+#os.environ[ 'HOME' ] = '/home/dsuarez/Pictures'
+#os.environ[ 'HOME' ] = '/tmp/'
+import matplotlib
+#if ide==1:
+# matplotlib.use('Qt4Agg')
+#elif ide==2:
+# matplotlib.use("Agg")
+#else:
+# matplotlib.use('TKAgg')
+#matplotlib.use("Agg")
+#matplotlib.interactive(1)
+import matplotlib.pyplot
+#import Numeric
+#import scipy
+import scipy.interpolate
+
+import Astro_Coords
+import TimeTools
+import Graphics_Miscens
+
+import Misc_Routines
+
+class AntPatternPlot:
+ def __init__(self):
+ """
+ AntPatternPlot creates an object to call methods to plot the antenna pattern.
+
+ Modification History
+ --------------------
+ Created by Freddy Galindo, ROJ, 06 October 2009.
+ """
+ self.figure = None
+ pass
+
+ def contPattern(self,iplot=0,gpath='',filename='',mesg='',amp=None ,x=None ,y=None ,getCut=None,title=''):
+ """
+ contPattern plots a contour map of the antenna pattern.
+
+ Parameters
+ ----------
+ iplot = A integer to specify if the plot is the first, second, ... The default va-
+ lue is 0.
+
+ Examples
+ --------
+ >> Over_Jro.JroPattern(pattern=2).contPattern()
+
+ Modification history
+ --------------------
+ Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009.
+ """
+
+ if getCut == 1:
+ return
+
+ xmax = numpy.max(x)
+ xmin = numpy.min(x)
+ ymax = numpy.max(y)
+ ymin = numpy.min(y)
+
+ levels = numpy.array([1e-3,1e-2,1e-1,0.5,1.0])
+ tmp = numpy.round(10*numpy.log10(levels),decimals=1)
+ labels = range(5)
+ for i in numpy.arange(5):labels[i] = str(numpy.int(tmp[i]))
+
+ if iplot==0:
+ xsize = 8.0
+ if matplotlib.get_backend()=='QT4Agg':xsize = 6.0
+ ysize = 8.0
+ self.figure = matplotlib.pyplot.figure(num=2,figsize=(xsize,ysize))
+ matplotlib.pyplot.clf()
+
+ colors = ((0,0,1.),(0,170/255.,0),(127/255.,1.,0),(1.,109/255.,0),(128/255.,0,0))
+ CS = matplotlib.pyplot.contour(x,y,amp.transpose(),levels,colors=colors)
+ fmt = {}
+ for l,s in zip(CS.levels,labels):fmt[l] = s
+
+ matplotlib.pyplot.annotate('Ng',xy=(-0.05,1.04),xytext=(0.01,0.962),xycoords='axes fraction',arrowprops=dict(facecolor='black', width=1.,shrink=0.2),fontsize=15.)
+ matplotlib.pyplot.annotate(mesg,xy=(0,0),xytext=(0.01,0.01),xycoords='figure fraction')
+ matplotlib.pyplot.clabel(CS,CS.levels,inline=True,fmt=fmt,fontsize=10)
+ matplotlib.pyplot.xlim(xmin,xmax)
+ matplotlib.pyplot.ylim(ymin,ymax)
+ matplotlib.pyplot.title("Total Pattern" + title)
+ matplotlib.pyplot.xlabel("West to South")
+ matplotlib.pyplot.ylabel("West to North")
+ matplotlib.pyplot.grid(True)
+ print "SAVE_FIG"
+ print gpath
+ print filename
+ save_fig = os.path.join(gpath,filename)
+ matplotlib.pyplot.savefig(save_fig,format='png')
+
+ def plotRaDec(self,gpath=None,filename=None,jd=2452640.5,ra_obs=None,xg=None,yg=None,x=None,y=None):
+ """
+ plotRaDec draws right ascension and declination lines on a JRO plane. This function
+ must call after conPattern.
+
+ Parameters
+ ----------
+ jd = A scalar giving the Julian date.
+ ra_obs = Scalar giving the right ascension of the observatory.
+ xg = A 3-element array to specify ..
+ yg = A 3-element array to specify ..
+
+ Examples
+ --------
+ >> Over_Jro.JroPattern(pattern=2).contPattern()
+ >> Over_Jro.JroPattern(pattern=2).plotRaDec(jd=jd,ra_obs=ra_obs,xg=xg,yg=yg)
+
+ Modification history
+ --------------------
+ Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009.
+ """
+
+ # Finding RA of observatory for a specific date
+ if ra_obs==None:ra_obs = numpy.array([23.37060849])
+ if xg==None:xg = numpy.array([0.62918474,-0.77725579,0.])
+ if yg==None:yg = numpy.array([0.77700346,0.62898048,0.02547905])
+
+ # Getting HA and DEC axes
+ mindec = -28; maxdec = 4; incdec = 2.
+ ndec = numpy.int((maxdec - mindec)/incdec) + 1
+
+ minha = -20; maxha = 20; incha = 2.
+ nha = numpy.int((maxha - minha)/incha) + 1
+
+ mcosx = numpy.zeros((nha,ndec))
+ mcosy = numpy.zeros((nha,ndec))
+
+ ha_axes = numpy.reshape(numpy.arange(nha)*incha + minha,(nha,1))
+ ones_dec = numpy.reshape(numpy.zeros(ndec) + 1,(ndec,1))
+ ha_axes = numpy.dot(ha_axes,ones_dec.transpose())
+ ha_axes2 = numpy.array(ra_obs - ha_axes)
+
+ dec_axes = numpy.reshape(numpy.arange(ndec)*incdec + mindec,(ndec,1))
+ ones_ra = numpy.reshape(numpy.zeros(nha) + 1,(nha,1))
+ dec_axes = numpy.dot(ones_ra,dec_axes.transpose())
+ dec_axes2 = numpy.array(dec_axes)
+
+ ObjHor = Astro_Coords.Equatorial(ha_axes2,dec_axes2,jd)
+ [alt,az,ha] = ObjHor.change2AltAz()
+
+ z = numpy.transpose(alt)*Misc_Routines.CoFactors.d2r ; z = z.flatten()
+ az = numpy.transpose(az)*Misc_Routines.CoFactors.d2r ; az = az.flatten()
+
+ vect = numpy.array([numpy.cos(z)*numpy.sin(az),numpy.cos(z)*numpy.cos(az),numpy.sin(z)])
+
+ xg = numpy.atleast_2d(xg)
+ dcosx = numpy.array(numpy.dot(xg,vect))
+ yg = numpy.atleast_2d(yg)
+ dcosy = numpy.array(numpy.dot(yg,vect))
+
+ mcosx = dcosx.reshape(ndec,nha)
+ mcosy = dcosy.reshape(ndec,nha)
+
+ # Defining NAN for points outof limits.
+ xmax = numpy.max(x)
+ xmin = numpy.min(x)
+ ymax = numpy.max(y)
+ ymin = numpy.min(y)
+
+ factor = 1.3
+ noval = numpy.where((mcosx>(xmax*factor)) | (mcosx<(xmin*factor)))
+ if noval[0].size>0:mcosx[noval] = numpy.nan
+ noval = numpy.where((mcosy>(ymax*factor)) | (mcosy<(ymin*factor)))
+ if noval[0].size>0:mcosy[noval] = numpy.nan
+
+ # Plotting HA and declination grid.
+ iha0 = numpy.int((0 - minha)/incha)
+ idec0 = numpy.int((-14 - mindec)/incdec)
+
+ colorgrid = (1.,109/255.,0)
+ matplotlib.pyplot.plot(mcosx.transpose(),mcosy.transpose(),color=colorgrid,linestyle='--')
+ for idec in numpy.arange(ndec):
+ if idec != idec0:
+ valx = (mcosx[idec,iha0]<=xmax) & (mcosx[idec,iha0]>=xmin)
+ valy = (mcosy[idec,iha0]<=ymax) & (mcosy[idec,iha0]>=ymin)
+ if valx & valy:
+ text = str(numpy.int(mindec + incdec*idec))+'$^o$'
+ matplotlib.pyplot.text(mcosx[idec,iha0],mcosy[idec,iha0],text)
+
+ matplotlib.pyplot.plot(mcosx,mcosy,color=colorgrid,linestyle='--')
+ for iha in numpy.arange(nha):
+ if iha != iha0:
+ valx = (mcosx[idec0,iha]<=xmax) & (mcosx[idec0,iha]>=xmin)
+ valy = (mcosy[idec0,iha]<=ymax) & (mcosy[idec0,iha]>=ymin)
+ if valx & valy:
+ text = str(4*numpy.int(minha + incha*iha))+"'"
+ matplotlib.pyplot.text(mcosx[idec0,iha],mcosy[idec0,iha],text)
+
+ matplotlib.pyplot.xlim(xmin,xmax)
+ matplotlib.pyplot.ylim(ymin,ymax)
+
+ save_fig = os.path.join(gpath,filename)
+ matplotlib.pyplot.savefig(save_fig,format='png')
+
+
+class BFieldPlot:
+ def __init__(self):
+ """
+ BFieldPlot creates an object for drawing magnetic Field lines over Jicamarca.
+
+ Modification History
+ --------------------
+ Created by Freddy Galindo, ROJ, 07 October 2009.
+ """
+
+ self.alpha_location = 1
+# pass
+
+ def plotBField(self,gpath,filename,dcos,alpha, nlon, nlat, dcosxrange, dcosyrange, heights, alpha_i):
+ """
+ plotBField draws the magnetic field in a directional cosines plot.
+
+ Parameters
+ ----------
+ dcos = An 4-dimensional array giving the directional cosines of the magnetic field
+ over the desired place.
+ alpha = An 3-dimensional array giving the angle of the magnetic field over the desi-
+ red place.
+ nlon = An integer to specify the number of elements per longitude.
+ nlat = An integer to specify the number of elements per latitude.
+ dcosxrange = A 2-element array giving the range of the directional cosines in the
+ "x" axis.
+ dcosyrange = A 2-element array giving the range of the directional cosines in the
+ "y" axis.
+ 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.
+ alpha_i = Angle to interpolate the magnetic field.
+ Modification History
+ --------------------
+ Converted to Python by Freddy R. Galindo, ROJ, 07 October 2009.
+ """
+
+ handles = []
+ objects = []
+ colors = ['k','m','c','b','g','r','y']
+ marker = ['-+','-*','-D','-x','-s','->','-o','-^']
+
+ alpha_location = numpy.zeros((nlon,2,heights.size))
+
+ for ih in numpy.arange(heights.size):
+ alpha_location[:,0,ih] = dcos[:,0,ih,0]
+ for ilon in numpy.arange(nlon):
+ myx = (alpha[ilon,:,ih])[::-1]
+ myy = (dcos[ilon,:,ih,0])[::-1]
+ tck = scipy.interpolate.splrep(myx,myy,s=0)
+ mydcosx = scipy.interpolate.splev(alpha_i,tck,der=0)
+
+ myx = (alpha[ilon,:,ih])[::-1]
+ myy = (dcos[ilon,:,ih,1])[::-1]
+ tck = scipy.interpolate.splrep(myx,myy,s=0)
+ mydcosy = scipy.interpolate.splev(alpha_i,tck,der=0)
+ alpha_location[ilon,:,ih] = numpy.array([mydcosx, mydcosy])
+
+
+ ObjFig, = matplotlib.pyplot.plot(alpha_location[:,0,ih],alpha_location[:,1,ih], \
+ marker[ih % 8],color=colors[numpy.int(ih/8)],ms=4.5,lw=0.5)
+ handles.append(ObjFig)
+ objects.append(numpy.str(heights[ih]) + ' km')
+
+ matplotlib.pyplot.xlim(dcosxrange[0],dcosxrange[1])
+ matplotlib.pyplot.ylim(dcosyrange[0],dcosyrange[1])
+
+ try:
+ ObjlegB = matplotlib.pyplot.legend(handles,objects,loc="lower right", numpoints=1, handlelength=0.3, \
+ handletextpad=0.02, borderpad=0.3, labelspacing=0.1)
+ except:
+ ObjlegB = matplotlib.pyplot.legend(handles,objects,loc=[0.01,0.75], numpoints=1, handlelength=0, \
+ pad=0.015, handletextsep=0.02,labelsep=0.01)
+
+ matplotlib.pyplot.setp(ObjlegB.get_texts(),fontsize='small')
+ matplotlib.pyplot.gca().add_artist(ObjlegB)
+
+ save_fig = os.path.join(gpath,filename)
+ matplotlib.pyplot.savefig(save_fig,format='png')
+ self.alpha_location = alpha_location
+
+
+class CelestialObjectsPlot:
+ def __init__(self,jd,dec,tod,maxha_min,show_object=None):
+
+ self.jd = jd
+ self.dec = dec
+ self.tod = tod
+ self.maxha_min = maxha_min
+
+ if show_object==None:show_object=numpy.zeros(4)+2
+ self.show_object = show_object
+
+ self.dcosx_sun = 1
+ self.dcosy_sun = 1
+ self.ha_sun = 1
+ self.time_sun = 1
+
+ self.dcosx_moon = 1
+ self.dcosy_moon = 1
+ self.ha_moon = 1
+ self.time_moon = 1
+
+ self.dcosx_hydra = 1
+ self.dcosy_hydra = 1
+ self.ha_hydra = 1
+ self.time_hydra = 1
+
+ self.dcosx_galaxy = 1
+ self.dcosy_galaxy = 1
+ self.ha_galaxy = 1
+ self.time_galaxy = 1
+
+ def drawObject(self,glat,glon,xg,yg,dcosxrange,dcosyrange,gpath='',filename=''):
+
+ jd = self.jd
+ main_dec = self.dec
+ tod = self.tod
+ maxha_min = self.maxha_min
+
+ mesg = "Drawing celestial objects over Observatory"
+# print mesg
+# if textid!=None:textid.append(mesg)
+
+ maxlev = 24; minlev = 0; maxcol = 39; mincol = 10
+ handles = []
+ objects = ['$Sun$','$Moon$','$Hydra$','$Galaxy$']
+ marker = ['--^','--s','--*','--o']
+
+ # Getting RGB table to plot celestial object over Jicamarca
+ colortable = Graphics_Miscens.ColorTable(table=1).readTable()
+
+ for io in (numpy.arange(4)+1):
+ if self.show_object[io]!=0:
+ ObjBodies = Astro_Coords.CelestialBodies()
+ if io==1:
+ [ra,dec,sunlon,sunobliq] = ObjBodies.sunpos(jd)
+ elif io==2:
+ [ra,dec,dist,moonlon,moonlat] = ObjBodies.moonpos(jd)
+ elif io==3:
+ [ra,dec] = ObjBodies.hydrapos()
+ elif io==4:
+ [maxra,ra] = ObjBodies.skynoise_jro(dec_cut=main_dec)
+ ra = maxra*15.
+ dec = main_dec
+
+ ObjEq = Astro_Coords.Equatorial(ra,dec,jd,lat=glat,lon=glon)
+ [alt, az, ha] = ObjEq.change2AltAz()
+ vect = numpy.array([az,alt]).transpose()
+ vect = Misc_Routines.Vector(vect,direction=0).Polar2Rect()
+
+ dcosx = numpy.array(numpy.dot(vect,xg))
+ dcosy = numpy.array(numpy.dot(vect,yg))
+ wrap = numpy.where(ha>=180.)
+ if wrap[0].size>0:ha[wrap] = ha[wrap] - 360.
+
+ val = numpy.where((numpy.abs(ha))<=(maxha_min*0.25))
+ if val[0].size>2:
+ tod_1 = tod*1.
+ shift_1 = numpy.where(tod>12.)
+ tod_1[shift_1] = tod_1[shift_1] - 24.
+ tod_2 = tod*1.
+ shift_2 = numpy.where(tod<12.)
+ tod_2[shift_2] = tod_2[shift_2] + 24.
+
+ diff0 = numpy.nanmax(tod[val]) - numpy.nanmin(tod[val])
+ diff1 = numpy.nanmax(tod_1[val]) - numpy.nanmin(tod_1[val])
+ diff2 = numpy.nanmax(tod_2[val]) - numpy.nanmin(tod_2[val])
+
+ if ((diff0<=diff1) & (diff0<=diff2)):
+ tod_0 = tod
+ elif ((diff10)
+ objects = numpy.array(objects)
+ objects = list(objects[val])
+ try:
+ ObjlegC = matplotlib.pyplot.legend(handles,objects,loc="lower left", numpoints=1, handlelength=0.3, \
+ borderpad=0.3, handletextpad=0.02,labelspacing=0.1)
+ except:
+ ObjlegC = matplotlib.pyplot.legend(handles,objects,loc=[0.01,0.75], numpoints=1, handlelength=0, \
+ pad=0.015, handletextsep=0.02,labelsep=0.01)
+
+ matplotlib.pyplot.setp(ObjlegC.get_texts(),fontsize='small')
+ ObjlegC.isaxes = False
+ save_fig = os.path.join(gpath,filename)
+ matplotlib.pyplot.savefig(save_fig,format='png')
+
+
+class PatternCutPlot:
+ def __init__(self,nsubplots):
+ self.nsubplots = nsubplots
+
+ self.fig = None
+
+ self.__plot_width = 8
+
+ if self.nsubplots == 5:
+ self.__plot_height = 11
+
+ if self.nsubplots == 4:
+ self.__plot_height = 9
+
+ if self.nsubplots == 3:
+ self.__plot_height = 7
+
+ if self.nsubplots == 2:
+ self.__plot_height = 5
+
+ if self.nsubplots == 1:
+ self.__plot_height = 3
+
+ self.fig = matplotlib.pyplot.figure(num = 4,figsize = (self.__plot_width, self.__plot_height))
+
+ if self.nsubplots < 5:
+ self.__height_inch = 1.1 #altura de los subplots (pulgadas)
+ top_inch = 1.5/2.7 #espacio entre el primer subplot y el limite superior del plot
+ self.__vspace_plot_inch = 1.0#1.5/2 # espacio vertical entre subplots
+ self.__left = 0.1
+ else:
+ self.__height_inch = 1.1 #altura de los subplots (pulgadas)
+ top_inch = 1.5/2.7 #espacio entre el primer subplot y el limite superior del plot
+ self.__vspace_plot_inch = 1.0 # espacio vertical entre subplots
+ self.__left = 0.1
+
+ self.__bottom_inch = self.__plot_height - (self.__height_inch + top_inch)
+ self.__height = self.__height_inch/self.__plot_height
+
+ self.__width = 0.8
+
+
+ def drawCut(self,io,patterns,npatterns,ha,otitle,subtitle,ptitle):
+
+ t_cuts = ['B','Sun','Moon','Hydra','Galaxy']
+ self.__bottom = self.__bottom_inch/self.__plot_height
+
+
+ subp = self.fig.add_axes([self.__left,self.__bottom,self.__width,self.__height])
+
+ on_axis_angle = -4.65562
+ for icut in numpy.arange(npatterns):
+ # Getting Antenna cut.
+ pattern = patterns[icut]
+ power = numpy.abs(pattern/numpy.nanmax(pattern))
+ max_power_db = numpy.round(10.*numpy.log10(numpy.nanmax(pattern)),2)
+
+ bval = numpy.where(power[:,0]==numpy.nanmax(power))
+ beta = -0.25*(ha[bval[0]] + on_axis_angle)
+# print 'Angle (deg): '+"%f"%(beta)
+
+ subp.plot(ha,power)
+
+
+ xmax = numpy.max(numpy.nanmin(ha))
+ xmin = numpy.min(numpy.nanmax(ha))
+ ymax = numpy.max(1)
+ ymin = numpy.min(0)
+
+
+ subp.set_xlim(xmin, xmax)
+
+ subp.set_ylim(ymin, ymax)
+
+ subp.set_title(otitle + ' ' + ptitle,size="medium")
+
+ subp.text(0.5, 1.26,subtitle[0],
+ horizontalalignment='center',
+ verticalalignment='center',
+ transform = subp.transAxes)
+
+ xlabels = subp.get_xticks()
+
+ subp.set_xticklabels(xlabels,size="small")
+
+ ylabels = subp.get_yticks()
+
+ subp.set_yticklabels(ylabels,size="small")
+
+ subp.set_xlabel('Hour angle (min) (+ve to West)',size="small")
+
+ subp.set_ylabel("Power [Max: " + str(max_power_db) + ' dB]',size="small")
+
+ subp.grid()
+
+
+ self.__bottom_inch = self.__bottom_inch - (self.__height_inch + self.__vspace_plot_inch)
+
+
+class SkyNoisePlot:
+ def __init__(self,date,powr,time,time_lst):
+ """
+ SkyNoisePlot class creates an object which represents the SkyNoise Object to genera-
+ te a SkyNoise map.
+
+ Parameters
+ ----------
+ date = A List of 3 elements to define the desired date ([year, month, day]).
+ powr = An array giving the SkyNoise power for the desired time.
+ time = An array giving the number of seconds since 1970 to the desired time.
+ time_lst = Set this input to an array to define the Local Sidereal Time of the desi-
+ red time.
+
+ Modification History
+ --------------------
+ Created by Freddy Galindo, ROJ, 18 October 2009.
+ """
+
+ self.date = date
+ self.powr = powr
+ self.time = time
+ self.time_lst = time_lst
diff --git a/apps/abs/utils/JroAntSetup.py b/apps/abs/utils/JroAntSetup.py
new file mode 100644
index 0000000..1c98353
--- /dev/null
+++ b/apps/abs/utils/JroAntSetup.py
@@ -0,0 +1,1035 @@
+'''
+The module JroAntSetup contains the pre-defined parameters for beam modelling of the Jicamarca ante-
+nna. Any new configuration must be added in this module (if the user decides that) using a specific
+ID (pattern value) or it would be read from a file using pattern=None.
+
+MODULES CALLED:
+OS, NUMPY
+
+MODIFICATION HISTORY:
+Created by Ing. Freddy Galindo (frederickgalindo@gmail.com). ROJ Sep 20, 2009.
+'''
+
+import os
+import numpy
+
+def ReturnSetup(path=None,filename=None,pattern=0):
+ """
+ ReturnSetup is a pre-defined list of Jicamarca antenna configurations which returns a dic-
+ tionary giving the configuration parameters (e.g. transmitted phases). To choose one, the
+ user must define the input "pattern" (See valid values below).
+
+ Parameters:
+ -----------
+
+ pattern = A integer (>=0) to specify the setup to choose. The default value is zero. If the
+ antenna configuration is user-defined pattern must be None.
+
+ path = Set this input a string to specifiy the folder path where the user-defined configu-
+ ration file is placed. If this value is not defined ReturnSetup will return None.
+
+ file = Set this input a string to specifiy the name of the user-defined configuration file
+ (*.txt). if this value is not defined ReturnSEtup will return None.
+
+ Examples:
+ ---------
+
+ Choosing a pre-defined antenna configuration
+ setup = ReturnSetup(pattern=1)
+
+ Reading a user-defined antenna configuration
+ setup = ReturnSetup(path="/users/users/Progs/Patterns/",file="ExpSep232009.txt")
+ """
+
+
+ if pattern == 0:
+ title = "for module (rx)"
+
+ ues = numpy.array([1.,2.,2.,1.])
+ phase = numpy.zeros([8,8])
+ phase[0:4,:] = 4
+ phase[4:8,:] = 5
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[0,0] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0,0] = 1
+
+ justrx = 1
+
+ elif pattern==1:
+ # Configuration 1/16 on-axis (rx)
+ title = " for 1/16 on-axis (rx)"
+
+ ues = numpy.array([1.,2.,2.,1.])
+ phase = numpy.zeros([8,8])
+ phase[0:4,:] = 4
+ phase[4:8,:] = 5
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[0:2,0:2] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0:2,0:2] = 1
+
+ justrx = 1
+
+ elif pattern == 2:
+ # Configuration for On-Axis
+ title = " for 1/4 on-axis (rx)"
+
+ ues = numpy.array([1.,2.,2.,1.])
+ phase = numpy.zeros([8,8])
+ phase[0:4,:] = 4
+ phase[4:8,:] = 5
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[0:4,0:4] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0:4,0:4] = 1
+
+ justrx = 1
+
+ elif pattern == 3:
+ # Configuration for On-Axis
+ title = " for all on-axis (rx)"
+
+ ues = numpy.array([1.,2.,2.,1.])
+ phase = numpy.zeros([8,8])
+ phase[0:4,:] = 4
+ phase[4:8,:] = 5
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[:,:] = 1
+
+ justrx = 0
+
+ elif pattern == 4:
+ # Configuration for oblique ISR On-Axis
+ title = " for Oblique ISR On-axis"
+
+ ues = numpy.array([1.,2.,2.,1.])
+ phase = numpy.zeros([8,8])
+ phase[0:4,:] = 4
+ phase[4:8,:] = 5
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[:,:] = 1
+
+ justrx = 0
+
+ elif pattern == 5:
+ # Configuration for oblique ISR "4.5"
+ title = " for Oblique ISR '4.5'"
+
+ ues = numpy.array([1.,2.,2.,1.])
+ phase = numpy.array([[4,4,5,5,2,2,3,3],
+ [4,5,5,2,2,3,3,4],
+ [5,5,2,2,3,3,4,4],
+ [5,2,2,3,3,4,4,5],
+ [3,3,4,4,5,5,2,2],
+ [3,4,4,5,5,2,2,3],
+ [4,4,5,5,2,2,3,3],
+ [4,5,5,2,2,3,3,4]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[:,:] = 1
+
+ justrx = 0
+
+ elif pattern == 6:
+ # Configuration for oblique ISR "6.0S"
+ title = " for Oblique ISR '6.0S'"
+
+ ues = numpy.array([1.,2.,2.,1.])
+ phase = numpy.array([[4,5,2,3,4,5,2,3],
+ [5,2,3,4,5,2,3,4],
+ [2,3,4,5,2,3,4,5],
+ [3,4,5,2,3,4,5,2],
+ [5,2,3,4,5,2,3,4],
+ [2,3,4,5,2,3,4,5],
+ [3,4,5,2,3,4,5,2],
+ [4,5,2,3,4,5,2,3]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[:,:] = 1
+
+ justrx = 0
+
+ elif pattern == 7:
+ # Configuration for oblique ISR "3.0N"
+ title = " for Oblique ISR '3.0N'"
+
+ ues = numpy.array([1.,2.,2.,1.])
+ phase = numpy.array([[4,3,2,5,4,3,2,5],
+ [3,2,5,4,3,2,5,4],
+ [2,5,4,3,2,5,4,3],
+ [5,4,3,2,5,4,3,2],
+ [5,4,3,2,5,4,3,2],
+ [4,3,2,5,4,3,2,5],
+ [3,2,5,4,3,2,5,4],
+ [2,5,4,3,2,5,4,3]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[:,:] = 1
+
+ justrx = 0
+
+ elif pattern == 8:
+ # Configuration for North Fritts"
+ title = " for North (Fritts)"
+
+ ues = numpy.array([2.513, 1.0, 3.0, 0.413])
+ phase = numpy.array([[4.29, 3.55, 2.82, 2.08, 4.20, 3.47, 2.73, 2.00],
+ [2.94, 2.20, 5.44, 4.70, 4.32, 3.59, 2.85, 2.12],
+ [5.56, 4.82, 4.09, 3.35, 4.44, 3.71, 2.97, 2.24],
+ [4.20, 3.47, 2.73, 2.00, 4.56, 3.82, 3.09, 2.35],
+ [4.20, 3.47, 2.73, 2.00, 4.56, 3.82, 3.09, 2.35],
+ [4.32, 3.59, 2.85, 2.12, 2.94, 2.20, 5.44, 4.70],
+ [4.44, 3.71, 2.97, 2.24, 5.56, 4.82, 4.09, 3.35],
+ [4.56, 3.82, 3.09, 2.35, 4.20, 3.47, 2.73, 2.00]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[0:4,0:4] = 1
+ gaintx[4:8,4:8] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0:4,0:4] = 1
+ gainrx[4:8,4:8] = 1
+
+ justrx = 0
+
+ elif pattern == 9:
+ # Configuration for West Fritts"
+ title = " for West (Fritts)"
+
+ ues = numpy.array([2.513, 1.0, 3.0, 0.413])
+ phase = numpy.array([[4.29, 3.55, 2.82, 2.08, 4.20, 3.47, 2.73, 2.00],
+ [2.94, 2.20, 5.44, 4.70, 4.32, 3.59, 2.85, 2.12],
+ [5.56, 4.82, 4.09, 3.35, 4.44, 3.71, 2.97, 2.24],
+ [4.20, 3.47, 2.73, 2.00, 4.56, 3.82, 3.09, 2.35],
+ [4.20, 3.47, 2.73, 2.00, 4.56, 3.82, 3.09, 2.35],
+ [4.32, 3.59, 2.85, 2.12, 2.94, 2.20, 5.44, 4.70],
+ [4.44, 3.71, 2.97, 2.24, 5.56, 4.82, 4.09, 3.35],
+ [4.56, 3.82, 3.09, 2.35, 4.20, 3.47, 2.73, 2.00]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[4:8,0:4] = 1
+ gaintx[4:8,0:4] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[4:8,0:4] = 1
+ gainrx[4:8,0:4] = 1
+
+ justrx = 0
+
+ elif pattern == 10:
+ # Configuration for South Fritts"
+ title = " for South (Fritts)"
+
+ ues = numpy.array([0.413, 2.0, 1.0, 1.513])
+ phase = numpy.array([[2.0 , 2.73, 3.47, 4.2 , 2.08, 2.82, 3.55, 4.29],
+ [2.12, 2.85, 3.59, 4.32, 4.7 , 5.44, 2.20, 2.94],
+ [2.24, 2.97, 3.71, 4.44, 3.35, 4.09, 4.82, 5.56],
+ [2.35, 3.09, 3.82, 4.56, 2.0 , 2.73, 3.47, 4.20],
+ [2.08, 2.82, 3.55, 4.29, 2.0 , 2.73, 3.47, 4.20],
+ [4.70, 5.44, 2.20, 2.94, 2.12, 2.85, 3.59, 4.32],
+ [3.35, 4.09, 4.82, 5.56, 2.24, 2.97, 3.71, 4.44],
+ [2.00, 2.73, 3.47, 4.20, 2.35, 3.09, 3.82, 4.56]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[0:4,0:4] = 1
+ gaintx[4:8,4:8] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0:4,0:4] = 1
+ gainrx[4:8,4:8] = 1
+
+ justrx = 0
+
+ elif pattern == 11:
+ # Configuration for East Fritts"
+ title = " for East (Fritts)"
+
+ ues = numpy.array([0.413, 2.0, 1.0, 1.513])
+ phase = numpy.array([[2.0 , 2.73, 3.47, 4.2 , 2.08, 2.82, 3.55, 4.29],
+ [2.12, 2.85, 3.59, 4.32, 4.7 , 5.44, 2.20, 2.94],
+ [2.24, 2.97, 3.71, 4.44, 3.35, 4.09, 4.82, 5.56],
+ [2.35, 3.09, 3.82, 4.56, 2.0 , 2.73, 3.47, 4.20],
+ [2.08, 2.82, 3.55, 4.29, 2.0 , 2.73, 3.47, 4.20],
+ [4.70, 5.44, 2.20, 2.94, 2.12, 2.85, 3.59, 4.32],
+ [3.35, 4.09, 4.82, 5.56, 2.24, 2.97, 3.71, 4.44],
+ [2.00, 2.73, 3.47, 4.20, 2.35, 3.09, 3.82, 4.56]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[4:8,0:4] = 1
+ gaintx[4:8,0:4] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[4:8,0:4] = 1
+ gainrx[4:8,0:4] = 1
+
+ justrx = 0
+
+ elif pattern == 12:
+ # Configuration for DEWD position (2009)
+ title = " for DEWD position (2009) East Beam"
+
+ ues = numpy.array([0.,0.,0.75,0.75])
+ phase = numpy.array([[2,3,3,3,3,4,4,4],
+ [5,2,2,2,2,3,3,3],
+ [3,4,4,4,4,5,5,5],
+ [2,3,3,3,3,4,4,4],
+ [4,5,5,5,5,2,2,2],
+ [3,4,4,4,4,5,5,5],
+ [5,2,2,2,2,3,3,3],
+ [4,5,5,5,5,2,2,2]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[:,0:4] = 1
+
+ justrx = 0
+
+ elif pattern == 13:
+ # Configuration for DEWD position (2009)
+ title = " for DEWD position (2009) West Beam"
+
+ ues = numpy.array([1.0,0.5,1.5,2.0])
+ phase = numpy.array([[5,4,2,5,3,2,4,3],
+ [2,5,3,2,4,3,5,4],
+ [2,5,3,2,4,3,5,4],
+ [3,2,4,3,5,4,2,5],
+ [3,2,4,3,5,4,2,5],
+ [4,3,5,4,2,5,3,2],
+ [4,3,5,4,2,5,3,2],
+ [5,4,2,5,3,2,4,3]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[4:8,:] = 1
+
+ justrx = 0
+
+ elif pattern == 14:
+ # Configuration for DVD position (2009)
+ title = " for DVD position (2009)"
+
+ ues = numpy.array([1.0,2.0,2.0,1.25])
+ phase = numpy.array([[2,2,5,5,4,4,3,3],
+ [2,5,5,4,4,3,3,2],
+ [5,5,4,4,3,3,2,2],
+ [5,4,4,3,3,2,2,5],
+ [5,5,4,4,3,3,2,2],
+ [5,4,4,3,3,2,2,5],
+ [4,4,3,3,2,2,5,5],
+ [4,3,3,2,2,5,5,4]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0:4,0:4] = 1
+
+ justrx = 0
+
+ elif pattern == 15:
+ # Configuration for Julia CP2
+ title = " for Julia CP2 Ew"
+
+ ues = numpy.array([0.0,1.0,1.0,0.0])
+ phase = numpy.array([[2,2,5,4,3,3,2,5],
+ [2,5,4,4,3,2,5,5],
+ [5,4,3,3,2,5,4,4],
+ [4,4,3,2,5,5,4,3],
+ [4,4,3,2,5,5,4,3],
+ [4,3,2,2,5,4,3,3],
+ [3,2,5,5,4,3,2,2],
+ [2,2,5,4,3,3,2,5]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[0:4,4:8] = 1
+ gaintx[4:8,0:4] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0,0] = 1
+
+ justrx = 0
+
+ elif pattern == 16:
+ # Configuration for Julia CP2
+ title = " for Julia CP2 NS"
+
+ ues = numpy.array([1.0,2.0,2.0,1.0])
+ phase = numpy.array([[4,4,3,2,5,5,4,3],
+ [4,3,2,2,5,4,3,3],
+ [3,2,5,5,4,3,2,2],
+ [2,2,5,4,3,3,2,5],
+ [2,2,5,4,3,3,2,5],
+ [2,5,4,4,3,2,5,5],
+ [5,4,3,3,2,5,4,4],
+ [4,4,3,2,5,5,4,3]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[0:4,0:4] = 1
+ gaintx[4:8,4:8] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0:4,0:4] = 1
+
+ justrx = 0
+
+ elif pattern == 17:
+ # Configuration for Julia CP3
+ title = " for Julia CP3 NS"
+
+ ues = numpy.array([1.0,1.0,1.0,1.0])
+ phase = numpy.array([[4,4,3,2,5,5,4,3],
+ [4,3,2,2,5,4,3,3],
+ [3,2,5,5,4,3,2,2],
+ [2,2,5,4,3,3,2,5],
+ [2,2,5,4,3,3,2,5],
+ [2,5,4,4,3,2,5,5],
+ [5,4,3,3,2,5,4,4],
+ [4,4,3,2,5,5,4,3]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[0:4,0:4] = 1
+ gaintx[4:8,4:8] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0:4,0:4] = 1
+
+ justrx = 0
+
+ elif pattern == 18:
+ # Configuration for Julia V
+ title = " for Julia V"
+
+ ues = (2/3.)*numpy.array([1.5,3.0+0.75,3.0,1.5-0.75])
+ phase = numpy.array([[4,4,3,3,2,2,5,5],
+ [4,3,3,2,2,5,5,4],
+ [3,3,2,2,5,5,4,4],
+ [3,2,2,5,5,4,4,3],
+ [3,3,2,2,5,5,4,4],
+ [3,2,2,5,5,4,4,3],
+ [2,2,5,5,4,4,3,3],
+ [2,5,5,4,4,3,3,2]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[0:4,0:4] = 1
+ gaintx[4:8,4:8] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0:4,0:4] = 1
+
+ justrx = 0
+
+ elif pattern == 19:
+ # Configuration for Julia V
+ title = " for Julia EW 2006-2007 (W)"
+
+ ues = numpy.array([1.0+0.66,2.0+0.66,2.0,1.0])
+ phase = numpy.array([[4,3,2,5,4,3,2,5],
+ [4,3,2,5,4,3,2,5],
+ [4,3,2,5,4,3,2,5],
+ [4,3,2,5,4,3,2,5],
+ [5,4,3,2,5,4,3,2],
+ [5,4,3,2,5,4,3,2],
+ [5,4,3,2,5,4,3,2],
+ [5,4,3,2,5,4,3,2]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[:,:] = 1
+
+ justrx = 0
+
+ elif pattern == 20:
+ # Configuration for Julia V
+ title = " for Julia EW 2006-2007 (E)"
+
+ ues = numpy.array([1.0,1.0,1.0,1.0])
+ phase = numpy.array([[4,4,4,4,5,5,5,5],
+ [3,3,3,3,4,4,4,4],
+ [5,5,5,5,2,2,2,2],
+ [4,4,4,4,5,5,5,5],
+ [2,2,2,2,3,3,3,3],
+ [5,5,5,5,2,2,2,2],
+ [3,3,3,3,4,4,4,4],
+ [2,2,2,2,3,3,3,3]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0:4,0:4] = 1
+ gainrx[4:8,4:8] = 1
+
+ justrx = 0
+
+ elif pattern == 21:
+ # Configuration for EW Imaging 1996
+ title = " for EW Imaging 1996"
+
+ ues = numpy.array([1.0,2.0,2.0,1.0])
+ phase = numpy.array([[4,4,3,2,5,5,4,3],
+ [4,3,2,2,5,4,3,3],
+ [3,2,5,5,4,3,2,2],
+ [2,2,5,4,3,3,2,5],
+ [2,2,5,4,3,3,2,5],
+ [2,5,4,4,3,2,5,5],
+ [5,4,3,3,2,5,4,4],
+ [4,4,3,2,5,5,4,3]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[0:4,0:4] = 1
+ gaintx[4:8,4:8] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0,0] = 1
+
+ justrx = 0
+
+ elif pattern == 22:
+ # Configuration for EW Imaging 2003
+ title = " for EW Imaging 2003"
+
+ ues = numpy.array([1.0,1.0,1.0,1.0])
+ phase = numpy.array([[4,4,3,2,0,0,0,0],
+ [2,3,2,2,0,0,0,0],
+ [5,0,2,5,0,0,0,0],
+ [2,4,3,4,0,0,0,0],
+ [0,0,0,0,3,3,2,5],
+ [0,0,0,0,2,2,5,5],
+ [0,0,0,0,4,3,5,4],
+ [0,0,0,0,5,3,2,3]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[0:4,0:4] = 1
+ gaintx[4:8,4:8] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0,0] = 1
+
+ justrx = 0
+
+ elif pattern == 23:
+ # Configuration for EW Imaging 2003
+ title = " for EW Imaging 2006-2008"
+
+ ues = numpy.array([1.0,1.0,1.0,2.0])
+ phase = numpy.array([[4,4,3,2,0,0,0,0],
+ [2,3,2,2,0,0,0,0],
+ [5,0,2,5,0,0,0,0],
+ [2,4,3,4,0,0,0,0],
+ [0,0,0,0,3,3,2,5],
+ [0,0,0,0,2,2,5,5],
+ [0,0,0,0,4,3,5,4],
+ [0,0,0,0,5,3,2,3]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[0:4,0:4] = 1
+ gaintx[4:8,4:8] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0,0] = 1
+
+ justrx = 0
+
+ elif pattern == 50:
+ # Configuration for vertical drift 1996
+ title = " for Vertical drift 1996"
+
+ ues = (2/3.)*numpy.array([0.,1.5,1.5,0.])
+ phase = numpy.array([[4,4,3,2,5,5,4,3],
+ [4,3,2,2,5,4,3,3],
+ [3,2,5,5,4,3,2,2],
+ [2,2,5,4,3,3,2,5],
+ [2,2,5,4,3,3,2,5],
+ [2,5,4,4,3,2,5,5],
+ [5,4,3,3,2,5,4,4],
+ [4,4,3,2,5,5,4,3]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[:,:] = 1
+
+ justrx = 0
+
+ elif pattern == 51:
+ # Configuration for vertical drift 1996
+ title = " for East-West Drifts 1996 (W beam)"
+
+ ues = numpy.array([0.0,1.0,2.0,1.0])
+ phase = numpy.array([[4,3,5,4,2,5,3,2],
+ [4,3,5,4,2,5,3,2],
+ [4,3,5,4,2,5,3,2],
+ [4,3,5,4,2,5,3,2],
+ [5,4,2,5,3,2,4,3],
+ [5,4,2,5,3,2,4,3],
+ [5,4,2,5,3,2,4,3],
+ [5,4,2,5,3,2,4,3]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[4:8,:] = 1
+
+ justrx = 0
+
+ elif pattern == 52:
+ # Configuration for vertical drift 1996
+ title = " for East-West Drifts 1996 (E Beam)"
+
+ ues = numpy.array([1.0,1.0,0.0,0.0])
+ phase = numpy.array([[4,4,4,4,5,5,5,5],
+ [3,3,3,3,4,4,4,4],
+ [5,5,5,5,2,2,2,2],
+ [4,4,4,4,5,5,5,5],
+ [2,2,2,2,3,3,3,3],
+ [5,5,5,5,2,2,2,2],
+ [3,3,3,3,4,4,4,4],
+ [2,2,2,2,3,3,3,3]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[:,0:4] = 1
+
+ justrx = 0
+
+ elif pattern == 53:
+ # Configuration for vertical drift 1996
+ title = " for DVD position 3 (2006-2008)"
+
+ ues = numpy.array([1.,2,2,1])
+ phase = numpy.array([[4,4,3,3,2,2,5,5],
+ [4,3,3,2,2,5,5,4],
+ [3,3,2,2,5,5,4,4],
+ [3,2,2,5,5,4,4,3],
+ [3,3,2,2,5,5,4,4],
+ [3,2,2,5,5,4,4,3],
+ [2,2,5,5,4,4,3,3],
+ [2,5,5,4,4,3,3,2]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0:4,4:8] = 1
+
+ justrx = 0
+
+ elif pattern == 54:
+ # Configuration for vertical drift 1996
+ title = " for DEWD (Mar 2005)"
+
+ ues = numpy.array([0.,1.,1/3.,1])
+ phase = numpy.array([[4,3,2,5,3,3,3,3],
+ [4,3,2,5,2,2,2,2],
+ [4,3,2,4,5,5,5,5],
+ [4,3,2,4,4,4,3,3],
+ [5,4,3,2,2,2,2,2],
+ [5,4,3,2,5,5,5,5],
+ [5,4,3,5,4,4,4,4],
+ [5,4,3,5,3,3,2,2]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[:,0:4] = 1
+
+ justrx = 0
+
+ elif pattern == 55:
+ # Configuration for vertical drift 1996
+ title = " for DEWD (Mar 2005)"
+
+ ues = numpy.array([0.,1.,1/3.,1])
+ phase = numpy.array([[4,3,2,5,3,3,3,3],
+ [4,3,2,5,2,2,2,2],
+ [4,3,2,4,5,5,5,5],
+ [4,3,2,4,4,4,3,3],
+ [5,4,3,2,2,2,2,2],
+ [5,4,3,2,5,5,5,5],
+ [5,4,3,5,4,4,4,4],
+ [5,4,3,5,3,3,2,2]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0:4:,4:8] = 1
+
+ justrx = 0
+
+ elif pattern ==56:
+ # Configuration using antenna compression
+ title = " for antenna compression AA*"
+
+ ues = numpy.array([0.0,0.0,0.0,0.0])
+ phase = numpy.array([[4,4,4,2,4,4,2,4],
+ [4,4,4,2,4,4,2,4],
+ [2,2,2,4,2,2,4,2],
+ [4,4,4,2,4,4,2,4],
+ [2,2,2,4,2,2,4,2],
+ [2,2,2,4,2,2,4,2],
+ [4,4,4,2,4,4,2,4],
+ [2,2,2,4,2,2,4,2]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0:4,0:4] = 1
+
+ justrx = 0
+
+ elif pattern ==57:
+ # Configuration using antenna compression
+ title = " for antenna compression AB*"
+
+ ues = numpy.array([0.0,0.0,0.0,0.0])
+ phase = numpy.array([[4,4,2,4,2,2,4,2],
+ [4,4,2,4,2,2,4,2],
+ [2,2,4,2,4,4,2,4],
+ [4,4,2,4,2,2,4,2],
+ [2,2,4,2,4,4,2,4],
+ [2,2,4,2,4,4,2,4],
+ [4,4,2,4,2,2,4,2],
+ [2,2,4,2,4,4,2,4]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0:4,0:4] = 1
+
+ justrx = 0
+
+ elif pattern ==58:
+ # Configuration using in Oblique ISR 4.5
+ title = " for Oblique ISR 4.5"
+
+ ues = numpy.array([1.0,2.0,2.0,1.0])
+ phase = numpy.array([[4,4,5,5,2,2,3,3],
+ [4,5,5,2,2,3,3,4],
+ [5,5,2,2,3,3,4,4],
+ [5,2,2,3,3,4,4,5],
+ [3,3,4,4,5,5,2,2],
+ [3,4,4,5,5,2,2,3],
+ [4,4,5,5,2,2,3,3],
+ [4,5,5,2,2,3,3,4]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[:,:] = 1
+
+ justrx = 1
+
+ elif pattern == 60:
+ title=" for Differential phase 2000"
+ ues = (2/3.)*numpy.array([0.,1.5-0.5,1.5,0.+0.5])
+
+ phase = numpy.array([[4,4,3,2,5,5,4,3],
+ [4,3,2,2,5,4,3,3],
+ [3,2,5,5,4,3,2,2],
+ [2,2,5,4,3,3,2,5],
+ [2,2,5,4,3,3,2,5],
+ [2,5,4,4,3,2,5,5],
+ [5,4,3,3,2,5,4,4],
+ [4,4,3,2,5,5,4,3]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0:4,0:4] = 1
+
+ justrx = 0
+
+ elif pattern == 61:
+ #for East-West 2003 W
+ title=" for East-West 2003"
+
+ ues = numpy.array([1.+0.66,2.+0.66,2.,1.])
+
+ phase = numpy.array([[4,3,2,5,4,3,2,5],
+ [4,3,2,5,4,3,2,5],
+ [4,3,2,5,4,3,2,5],
+ [4,3,2,5,4,3,2,5],
+ [5,4,3,2,5,4,3,2],
+ [5,4,3,2,5,4,3,2],
+ [5,4,3,2,5,4,3,2],
+ [5,4,3,2,5,4,3,2]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[4:8,:] = 1
+
+ justrx = 0
+
+ elif pattern == 62:
+ #for East-West 2003 E
+ title=" for East-West 2003"
+
+ ues = numpy.array([1.,1.,0.+1.0,0.+1.0])
+
+ phase = numpy.array([[4,4,4,4,5,5,5,5],
+ [3,3,3,3,4,4,4,4],
+ [5,5,5,5,2,2,2,2],
+ [4,4,4,4,5,5,5,5],
+ [2,2,2,2,3,3,3,3],
+ [5,5,5,5,2,2,2,2],
+ [3,3,3,3,4,4,4,4],
+ [2,2,2,2,3,3,3,3]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[:,0:4] = 1
+
+ justrx = 0
+
+ elif pattern == 63:
+
+ title=" for Differential phase 2004 High Alt."
+
+ ues = (2/3.)*numpy.array([0.,1.5-1.0,1.5,0.+1.0])
+
+ phase = numpy.array([[4,4,3,2,5,5,4,3],
+ [4,3,2,2,5,4,3,3],
+ [3,2,5,5,4,3,2,2],
+ [2,2,5,4,3,3,2,5],
+ [2,2,5,4,3,3,2,5],
+ [2,5,4,4,3,2,5,5],
+ [5,4,3,3,2,5,4,4],
+ [4,4,3,2,5,5,4,3]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[:,:] = 1
+
+ justrx = 0
+
+ elif pattern == 64:
+
+ title=" for Differential Phase Perp to B 2005-2006"
+
+ ues = (2/3.)*numpy.array([1.5,3.0+0.75,3.0,1.5-0.75])
+
+ phase = numpy.array([[4,4,3,3,2,2,5,5],
+ [4,3,3,2,2,5,5,4],
+ [3,3,2,2,5,5,4,4],
+ [3,2,2,5,5,4,4,3],
+ [3,3,2,2,5,5,4,4],
+ [3,2,2,5,5,4,4,3],
+ [2,2,5,5,4,4,3,3],
+ [2,5,5,4,4,3,3,2]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[0:4,4:8] = 1
+
+ justrx = 0
+
+ elif pattern == 65:
+ #for JULIA EW 2003 W
+ title=" for JULIA EW 2003"
+
+ ues = numpy.array([1+0.66,2+0.66,2.,1.])
+
+ phase = numpy.array([[4,3,2,5,4,3,2,5],
+ [4,3,2,5,4,3,2,5],
+ [4,3,2,5,4,3,2,5],
+ [4,3,2,5,4,3,2,5],
+ [5,4,3,2,5,4,3,2],
+ [5,4,3,2,5,4,3,2],
+ [5,4,3,2,5,4,3,2],
+ [5,4,3,2,5,4,3,2]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[4:8,:] = 1
+
+ justrx = 0
+
+ elif pattern == 66:
+ #for JULIA EW 2003 E
+ title=" for JULIA EW 2003"
+
+ ues = numpy.array([1.,1.,0.,0.])
+
+ phase = numpy.array([[4,4,4,4,5,5,5,5],
+ [3,3,3,3,4,4,4,4],
+ [5,5,5,5,2,2,2,2],
+ [4,4,4,4,5,5,5,5],
+ [2,2,2,2,3,3,3,3],
+ [5,5,5,5,2,2,2,2],
+ [3,3,3,3,4,4,4,4],
+ [2,2,2,2,3,3,3,3]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[:,0:4] = 1
+
+ justrx = 0
+
+ elif pattern == 67:
+
+ title=" for Vertical (Yellow Cables)"
+
+ ues = numpy.array([0.25, 0.25, 0.25, 0.25])
+
+ phase = numpy.array([[3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41],
+ [2.78, 2.78, 2.78, 2.78, 2.78, 2.78, 2.78, 2.78],
+ [2.15, 2.15, 2.15, 2.15, 2.15, 2.15, 2.15, 2.15],
+ [5.52, 5.52, 5.52, 5.52, 5.52, 5.52, 5.52, 5.52],
+ [4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89, 4.89],
+ [4.26, 4.26, 4.26, 4.26, 4.26, 4.26, 4.26, 4.26],
+ [3.63, 3.63, 3.63, 3.63, 3.63, 3.63, 3.63, 3.63],
+ [3.00, 3.00, 3.00, 3.00, 3.00, 3.00, 3.00, 3.00]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[:,:] = 1
+
+ justrx = 0
+
+ elif pattern == 100:
+
+ title=" for High Altitude Drift"
+
+ ues = numpy.array([ 2.0,0.8,1.0,2.2])
+
+ phase = numpy.array([[5,5,4,4, 4,4,3,3],
+ [5,4,4,4, 4,3,3,3],
+ [4,4,4,4, 3,3,3,3],
+ [4,4,4,3, 3,3,3,2],
+ [3,3,2,2, 2,2,5,5],
+ [3,2,2,2, 2,5,5,5],
+ [2,2,2,2, 5,5,5,5],
+ [2,2,2,5, 5,5,5,4]],dtype=float)
+
+ gaintx = numpy.zeros([8,8])
+ gaintx[:,:] = 1
+
+ gainrx = numpy.zeros([8,8])
+ gainrx[:,:] = 1
+
+ justrx = 0
+
+
+ elif pattern==None:
+
+ inputs = numpy.array(["title","ues_tx","phase_tx","gain_tx","gain_rx","just_rx"])
+
+ # Reading user-defined configuration.
+ if path==None:path = os.getcwd() + os.sep + "patterns" + os.sep
+ if filename==None:filename = "jropattern.txt"
+
+ ff = open(os.path.join(path,filename),'r')
+
+ while 1:
+ # Checking EOF.
+ init = ff.tell()
+ if not ff.readline():break
+ else:ff.seek(init)
+
+ line = ff.readline().lstrip()
+ if line.__len__()!=0:
+ if line[0]!='#':
+ keys = line.split("=")
+ key = keys[0].lstrip().rstrip().lower()
+ vv = numpy.where(inputs==key)
+ if vv[0][0]==0:
+ title = keys[1].lstrip().rstrip()
+ elif vv[0][0]==1:
+ ues = (keys[1].lstrip().rstrip())
+ ues = numpy.float32(ues[1:-1].split(","))
+ elif vv[0][0]==2:
+ phase = numpy.zeros([8,8])
+ tx = (keys[1].lstrip().rstrip())
+ tx = numpy.float32(tx[2:-3].split(","))
+ phase[0,:] = tx
+ for ii in numpy.arange(7):
+ tx = ff.readline().lstrip().rstrip()
+ tx = numpy.float32(tx[1:-3+(ii==6)].split(","))
+ phase[ii+1,:] = tx
+ elif vv[0][0]==3:
+ gaintx = numpy.zeros([8,8])
+ gg = (keys[1].lstrip().rstrip())
+ gg = numpy.float32(gg[2:-3].split(","))
+ gaintx[0,:] = gg
+ for ii in numpy.arange(7):
+ gg = ff.readline().lstrip().rstrip()
+ gg = numpy.float32(gg[1:-3+(ii==6)].split(","))
+ gaintx[ii+1,:] = gg
+ elif vv[0][0]==4:
+ gainrx = numpy.zeros([8,8])
+ gg = (keys[1].lstrip().rstrip())
+ gg = numpy.float32(gg[2:-3].split(","))
+ gainrx[0,:] = gg
+ for ii in numpy.arange(7):
+ gg = ff.readline().lstrip().rstrip()
+ gg = numpy.float32(gg[1:-3+(ii==6)].split(","))
+ gainrx[ii+1,:] = gg
+ elif vv[0][0]==5:
+ justrx = numpy.float(keys[1].lstrip().rstrip())
+
+ ff.close()
+
+
+
+ setup = {"ues":ues, "phase":phase, "gaintx":gaintx, "gainrx":gainrx, "justrx":justrx, \
+ "title":title}
+
+ return setup
+
+
+
+
\ No newline at end of file
diff --git a/apps/abs/utils/Misc_Routines.py b/apps/abs/utils/Misc_Routines.py
new file mode 100644
index 0000000..886cbb1
--- /dev/null
+++ b/apps/abs/utils/Misc_Routines.py
@@ -0,0 +1,81 @@
+"""
+The module MISC_ROUTINES gathers classes and functions which are useful for daily processing. As an
+example we have conversion factor or universal constants.
+
+MODULES CALLED:
+NUMPY, SYS
+
+MODIFICATION HISTORY:
+Created by Ing. Freddy Galindo (frederickgalindo@gmail.com). ROJ, 21 October 2009.
+"""
+
+import numpy
+import sys
+
+class CoFactors():
+ """
+ CoFactor class used to call pre-defined conversion factor (e.g. degree to radian). The cu-
+ The current available factor are:
+
+ d2r = degree to radian.
+ s2r = seconds to radian?, degree to arcsecond.?
+ h2r = hour to radian.
+ h2d = hour to degree
+ """
+
+ d2r = numpy.pi/180.
+ s2r = numpy.pi/(180.*3600.)
+ h2r = numpy.pi/12.
+ h2d = 15.
+
+class Redirect:
+ def __init__(self,stdout):
+ self.stdout = stdout
+
+ def write(self,message):
+ self.stdout.insertPlainText(message)
+
+class WidgetPrint:
+ """
+ WidgetPrint class allows to define the standard output.
+ """
+ def __init__(self,textid=None):
+ self.__stdout = sys.stdout
+ self.textid = textid
+ self.wPrint()
+
+ def wPrint(self):
+ if self.textid == None: sys.stdout = self.__stdout
+ if self.textid != None: sys.stdout = Redirect(self.textid)
+ print ("")
+
+class Vector:
+ """
+ direction = 0 Polar to rectangular; direction=1 rectangular to polar
+ """
+ def __init__(self,vect,direction=0):
+ nsize = numpy.size(vect)
+ if nsize <= 3:
+ vect = vect.reshape(1,nsize)
+
+ self.vect = vect
+ self.dirc = direction
+
+
+
+ def Polar2Rect(self):
+ if self.dirc == 0:
+ jvect = self.vect*numpy.pi/180.
+ mmx = numpy.cos(jvect[:,1])*numpy.sin(jvect[:,0])
+ mmy = numpy.cos(jvect[:,1])*numpy.cos(jvect[:,0])
+ mmz = numpy.sin(jvect[:,1])
+ mm = numpy.array([mmx,mmy,mmz]).transpose()
+
+ elif self.dirc == 1:
+ mm = [numpy.arctan2(self.vect[:,0],self.vect[:,1]),numpy.arcsin(self.vect[:,2])]
+ mm = numpy.array(mm)*180./numpy.pi
+
+ return mm
+
+
+
diff --git a/apps/abs/utils/OverJRO.py b/apps/abs/utils/OverJRO.py
new file mode 100644
index 0000000..77c3506
--- /dev/null
+++ b/apps/abs/utils/OverJRO.py
@@ -0,0 +1,97 @@
+'''
+Created on May 8, 2013
+
+@author: Jose Antonio Sal y Rosas Celi
+@contact: jose.salyrosas@jro.igp.gob.pe
+'''
+
+import os
+from Files import Files
+
+class OverJRO(Files):
+
+ __scriptName = "OverJRO.py"
+
+ def __init__(self):
+ pass
+
+ def setParameters(self, path, exp_name, phase_tx, gain_tx, gain_rx, ues_tx, just_rx):
+ self.path = path
+ self.exp_name = exp_name
+ self.phase_tx = phase_tx
+ self.gain_tx = gain_tx
+ self.gain_rx = gain_rx
+ self.ues_tx = ues_tx
+ self.just_rx = just_rx
+
+ def saveFile(self, contentFile):
+ filename = self.setFilename()
+ finalpath = os.path.join(self.path, self.setFileExtension(filename))
+ print "HAHAH"
+ finalpath = "apps/abs/static/data/"+finalpath
+ self.save(finalpath, contentFile)
+ return finalpath
+
+ def setTextContent(self):
+ title = "title ='%s'" % self.exp_name
+ ues_tx = "ues_tx = %s" % self.ues_tx
+ phase_tx = "phase_tx = %s" % (self.convertValue(self.phase_tx))
+ gain_tx = "gain_tx = %s" % (self.convertValue(self.gain_tx))
+ gain_rx = "gain_rx = %s" % (self.convertValue(self.gain_rx))
+ just_rx = "just_rx = %d" % self.just_rx
+ content = " %s\r\n\n %s\r\n\n %s\r\n %s\r\n %s\r\n %s\r\n" % (title, ues_tx, phase_tx, gain_tx, gain_rx, just_rx)
+ return content
+
+ def setFileExtension(self, filename):
+ txtFile = filename + ".txt"
+
+ return txtFile
+
+ def convertValue(self, strAntenna):
+ value = ""
+ strAntenna = strAntenna.replace("],[","]+[")
+ lsAntenna = strAntenna.split("+")
+ for i,element in enumerate(lsAntenna):
+ if i == 0:
+ value += "%s,$\n" % element
+ elif i == 7:
+ value += " %s\n" % element
+ else:
+ value += " %s,$\n" % element
+
+ return value
+
+
+if __name__ == '__main__':
+ path = "/home/fquino/workspace/radarsys/webapp/apps/abs/static/data"
+ exp_name = "MST-ISR 2009 (NS-Up)"
+ phase_tx = "[[0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5]," \
+ "[1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0]," \
+ "[0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5]," \
+ "[0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5]," \
+ "[1.0,1.0,1.0,1.0,1.5,1.5,1.5,1.5]," \
+ "[0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5]," \
+ "[1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0]," \
+ "[0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5]]"
+ gain_tx = "[[1,1,1,1,1,1,1,1]," \
+ "[1,1,1,1,1,1,1,1]," \
+ "[1,1,1,1,1,1,1,1]," \
+ "[1,1,1,1,1,1,1,1]," \
+ "[1,1,1,1,1,1,1,1]," \
+ "[1,1,1,1,1,1,1,1]," \
+ "[1,1,1,1,1,1,1,1]," \
+ "[1,1,1,1,1,1,1,1]]"
+ gain_rx = "[[1,1,1,1,0,0,0,0]," \
+ "[1,1,1,1,0,0,0,0]," \
+ "[1,1,1,1,0,0,0,0]," \
+ "[1,1,1,1,0,0,0,0]," \
+ "[0,0,0,0,1,1,1,1]," \
+ "[0,0,0,0,1,1,1,1]," \
+ "[0,0,0,0,1,1,1,1]," \
+ "[0,0,0,0,1,1,1,1]]"
+ ues_tx = "[0.533333,0.00000,1.06667,0.00000]"
+ just_rx = 0
+ data = OverJRO()
+ data.setParameters(path, exp_name, phase_tx, gain_tx, gain_rx, ues_tx, just_rx)
+ contentFile = data.setTextContent()
+ data.saveFile(contentFile)
diff --git a/apps/abs/utils/TimeTools.py b/apps/abs/utils/TimeTools.py
new file mode 100644
index 0000000..19f8f9c
--- /dev/null
+++ b/apps/abs/utils/TimeTools.py
@@ -0,0 +1,427 @@
+"""
+The TIME_CONVERSIONS.py module gathers classes and functions for time system transformations
+(e.g. between seconds from 1970 to datetime format).
+
+MODULES CALLED:
+NUMPY, TIME, DATETIME, CALENDAR
+
+MODIFICATION HISTORY:
+Created by Ing. Freddy Galindo (frederickgalindo@gmail.com). ROJ Aug 13, 2009.
+"""
+
+import numpy
+import time
+#import datetime as dt
+import calendar
+
+class Time:
+ """
+ time(year,month,dom,hour,min,secs)
+
+ An object represents a date and time of certain event..
+
+ Parameters
+ ----------
+ YEAR = Number of the desired year. Year must be valid values from the civil calendar.
+ Years B.C.E must be represented as negative integers. Years in the common era are repre-
+ sented as positive integers. In particular, note that there is no year 0 in the civil
+ calendar. 1 B.C.E. (-1) is followed by 1 C.E. (1).
+
+ MONTH = Number of desired month (1=Jan, ..., 12=December).
+
+ DOM = Number of day of the month.
+
+ HOUR = Number of the hour of the day. By default hour=0
+
+ MINS = Number of the minute of the hour. By default min=0
+
+ SECS = Number of the second of the minute. By default secs=0.
+
+ Examples
+ --------
+ time_info = time(2008,9,30,12,30,00)
+
+ time_info = time(2008,9,30)
+ """
+
+ def __init__(self,year=None,month=None,dom=None,hour=0,mins=0,secs=0):
+ # If one the first three inputs are not defined, it takes the current date.
+ date = time.localtime()
+ if year==None:year=date[0]
+ if month==None:month=date[1]
+ if dom==None:dom=date[2]
+
+ # Converting to arrays
+ year = numpy.array([year]); month = numpy.array([month]); dom = numpy.array([dom])
+ hour = numpy.array([hour]); mins = numpy.array([mins]); secs = numpy.array([secs])
+
+ # Defining time information object.
+ self.year = numpy.atleast_1d(year)
+ self.month = numpy.atleast_1d(month)
+ self.dom = numpy.atleast_1d(dom)
+ self.hour = numpy.atleast_1d(hour)
+ self.mins = numpy.atleast_1d(mins)
+ self.secs = numpy.atleast_1d(secs)
+
+ def change2julday(self):
+ """
+ Converts a datetime to Julian days.
+ """
+
+ # Defining constants
+ greg = 2299171 # incorrect Julian day for Oct, 25, 1582.
+ min_calendar = -4716
+ max_calendar = 5000000
+
+ min_year = numpy.nanmin(self.year)
+ max_year = numpy.nanmax(self.year)
+ if (min_yearmax_calendar):
+ print ("Value of Julian date is out of allowed range")
+ return -1
+
+ noyear = numpy.sum(self.year==0)
+ if noyear>0:
+ print ("There is no year zero in the civil calendar")
+ return -1
+
+ # Knowing if the year is less than 0.
+ bc = self.year<0
+
+ # Knowing if the month is less than March.
+ inJanFeb = self.month<=2
+
+ jy = self.year + bc - inJanFeb
+ jm = self.month + (1 + 12*inJanFeb)
+
+ # Computing Julian days.
+ jul= numpy.floor(365.25*jy) + numpy.floor(30.6001*jm) + (self.dom+1720995.0)
+
+ # Test whether to change to Gregorian Calendar
+ if numpy.min(jul) >= greg:
+ ja = numpy.int32(0.01*jy)
+ jul = jul + 2 - ja + numpy.int32(0.25*ja)
+ else:
+ gregchange = numpy.where(jul >= greg)
+ if gregchange[0].size>0:
+ ja = numpy.int32(0.01 + jy[gregchange])
+ jy[gregchange] = jy[gregchange] + 2 - ja + numpy.int32(0.25*ja)
+
+ # Determining machine-specific parameters affecting floating-point.
+ eps = 0.0 # Replace this line for a function to get precision.
+ eps = abs(jul)*0.0 > eps
+
+ jul = jul + (self.hour/24. -0.5) + (self.mins/1440.) + (self.secs/86400.) + eps
+
+ return jul[0]
+
+ def change2secs(self):
+ """
+ Converts datetime to number of seconds respect to 1970.
+ """
+
+ year = self.year
+ if year.size>1: year = year[0]
+
+ month = self.month
+ if month.size>1: month = month[0]
+
+ dom = self.dom
+ if dom.size>1: dom = dom[0]
+
+ # Resizing hour, mins and secs if it was necessary.
+ hour = self.hour
+ if hour.size>1:hour = hour[0]
+ if hour.size==1:hour = numpy.resize(hour,year.size)
+
+ mins = self.mins
+ if mins.size>1:mins = mins[0]
+ if mins.size==1:mins = numpy.resize(mins,year.size)
+
+ secs = self.secs
+ if secs.size>1:secs = secs[0]
+ if secs.size==1:secs = numpy.resize(secs,year.size)
+
+ # Using time.mktime to compute seconds respect to 1970.
+ secs1970 = numpy.zeros(year.size)
+ for ii in numpy.arange(year.size):
+ secs1970[ii] = time.mktime((int(year[ii]),int(month[ii]),int(dom[ii]),\
+ int(hour[ii]),int(mins[ii]),int(secs[ii]),0,0,0))
+
+ secs1970 = numpy.int32(secs1970 - time.timezone)
+
+ return secs1970
+
+ def change2strdate(self,mode=1):
+ """
+ change2strdate method converts a date and time of certain event to date string. The
+ string format is like localtime (e.g. Fri Oct 9 15:00:19 2009).
+
+ Parameters
+ ----------
+ None.
+
+ Return
+ ------
+
+ Modification History
+ --------------------
+ Created by Freddy R. Galindo, ROJ, 09 October 2009.
+
+ """
+
+ secs = numpy.atleast_1d(self.change2secs())
+ strdate = []
+ for ii in numpy.arange(numpy.size(secs)):
+ secs_tmp = time.localtime(secs[ii] + time.timezone)
+ if mode==1:
+ strdate.append(time.strftime("%d-%b-%Y (%j) %H:%M:%S",secs_tmp))
+ elif mode==2:
+ strdate.append(time.strftime("%d-%b-%Y (%j)",secs_tmp))
+
+ strdate = numpy.array(strdate)
+
+ return strdate
+
+
+class Secs:
+ """
+ secs(secs):
+
+ An object represents the number of seconds respect to 1970.
+
+ Parameters
+ ----------
+
+ SECS = A scalar or array giving the number of seconds respect to 1970.
+
+ Example:
+ --------
+ secs_info = secs(1251241373)
+
+ secs_info = secs([1251241373,1251241383,1251241393])
+ """
+ def __init__(self,secs):
+ self.secs = secs
+
+ def change2julday(self):
+ """
+ Convert seconds from 1970 to Julian days.
+ """
+
+ secs_1970 = time(1970,1,1,0,0,0).change2julday()
+
+ julian = self.secs/86400.0 + secs_1970
+
+ return julian
+
+ def change2time(self):
+ """
+ Converts seconds from 1970 to datetime.
+ """
+
+ secs1970 = numpy.atleast_1d(self.secs)
+
+ datetime = numpy.zeros((9,secs1970.size))
+ for ii in numpy.arange(secs1970.size):
+ tuple = time.gmtime(secs1970[ii])
+ datetime[0,ii] = tuple[0]
+ datetime[1,ii] = tuple[1]
+ datetime[2,ii] = tuple[2]
+ datetime[3,ii] = tuple[3]
+ datetime[4,ii] = tuple[4]
+ datetime[5,ii] = tuple[5]
+ datetime[6,ii] = tuple[6]
+ datetime[7,ii] = tuple[7]
+ datetime[8,ii] = tuple[8]
+
+ datetime = numpy.int32(datetime)
+
+ return datetime
+
+
+class Julian:
+ """
+ julian(julian):
+
+ An object represents julian days.
+
+ Parameters
+ ----------
+
+ JULIAN = A scalar or array giving the julina days.
+
+ Example:
+ --------
+ julian_info = julian(2454740)
+
+ julian_info = julian([2454740,2454760,2454780])
+ """
+ def __init__(self,julian):
+ self.julian = numpy.atleast_1d(julian)
+
+ def change2time(self):
+ """
+ change2time method converts from julian day to calendar date and time.
+
+ Return
+ ------
+ year = An array giving the year of the desired julian day.
+ month = An array giving the month of the desired julian day.
+ dom = An array giving the day of the desired julian day.
+ hour = An array giving the hour of the desired julian day.
+ mins = An array giving the minute of the desired julian day.
+ secs = An array giving the second of the desired julian day.
+
+ Examples
+ --------
+ >> jd = 2455119.0
+ >> [yy,mo,dd,hh,mi,ss] = TimeTools.julian(jd).change2time()
+ >> print ([yy,mo,dd,hh,mi,ss])
+ [2009] [10] [ 14.] [ 12.] [ 0.] [ 0.]
+
+ Modification history
+ --------------------
+ Translated from "Numerical Recipies in C", by William H. Press, Brian P. Flannery,
+ Saul A. Teukolsky, and William T. Vetterling. Cambridge University Press, 1988.
+ Converted to Python by Freddy R. Galindo, ROJ, 06 October 2009.
+ """
+
+ min_julian = -1095
+ max_julian = 1827933925
+ if (numpy.min(self.julian) < min_julian) or (numpy.max(self.julian) > max_julian):
+ print ('Value of Julian date is out of allowed range.')
+ return None
+
+ # Beginning of Gregorian calendar
+ igreg = 2299161
+ julLong = numpy.floor(self.julian + 0.5)
+ minJul = numpy.min(julLong)
+
+ if (minJul >= igreg):
+ # All are Gregorian
+ jalpha = numpy.int32(((julLong - 1867216) - 0.25)/36524.25)
+ ja = julLong + 1 + jalpha - numpy.int32(0.25*jalpha)
+ else:
+ ja = julLong
+ gregChange = numpy.where(julLong >= igreg)
+ if gregChange[0].size>0:
+ jalpha = numpy.int32(((julLong[gregChange]-1867216) - 0.25)/36524.25)
+ ja[gregChange] = julLong[gregChange]+1+jalpha-numpy.int32(0.25*jalpha)
+
+ # clear memory.
+ jalpha = -1
+
+ jb = ja + 1524
+ jc = numpy.int32(6680. + ((jb-2439870)-122.1)/365.25)
+ jd = numpy.int32(365.*jc + (0.25*jc))
+ je = numpy.int32((jb - jd)/30.6001)
+
+ dom = jb - jd - numpy.int32(30.6001*je)
+ month = je - 1
+ month = ((month - 1) % 12) + 1
+ month = numpy.atleast_1d(month)
+ year = jc - 4715
+ year = year - (month > 2)*1
+ year = year - (year <= 0)*1
+ year = numpy.atleast_1d(year)
+
+ # Getting hours, minutes, seconds
+ fraction = self.julian + 0.5 - julLong
+ eps_0 = dom*0.0 + 1.0e-12
+ eps_1 = 1.0e-12*numpy.abs(julLong)
+ eps = (eps_0>eps_1)*eps_0 + (eps_0<=eps_1)*eps_1
+
+ hour_0 = dom*0 + 23
+ hour_2 = dom*0 + 0
+ hour_1 = numpy.floor(fraction*24.0 + eps)
+ hour = ((hour_1>hour_0)*23) + ((hour_1<=hour_0)*hour_1)
+ hour = ((hour_1=hour_2)*hour_1)
+
+ fraction = fraction - (hour/24.0)
+ mins_0 = dom*0 + 59
+ mins_2 = dom*0 + 0
+ mins_1 = numpy.floor(fraction*1440.0 + eps)
+ mins = ((mins_1>mins_0)*59) + ((mins_1<=mins_0)*mins_1)
+ mins = ((mins_1=mins_2)*mins_1)
+
+ secs_2 = dom*0 + 0
+ secs_1 = (fraction - mins/1440.0)*86400.0
+ secs = ((secs_1=secs_2)*secs_1)
+
+ return year,month,dom,hour,mins,secs
+
+ def change2secs(self):
+ """
+ Converts from Julian days to seconds from 1970.
+ """
+
+ jul_1970 = Time(1970,1,1,0,0,0).change2julday()
+
+ secs = numpy.int32((self.julian - jul_1970)*86400)
+
+ return secs
+
+ def change2lst(self,longitude=-76.8667):
+ """
+ CT2LST converts from local civil time to local mean sideral time
+
+ longitude = The longitude in degrees (east of Greenwich) of the place for which
+ the local sideral time is desired, scalar. The Greenwich mean sideral time (GMST)
+ can be found by setting longitude=0.
+ """
+
+ # Useful constants, see Meus, p. 84
+ c = numpy.array([280.46061837, 360.98564736629, 0.000387933, 38710000.0])
+ jd2000 = 2451545.0
+ t0 = self.julian - jd2000
+ t = t0/36525.
+
+ # Computing GST in seconds
+ theta = c[0] + (c[1]*t0) + (t**2)*(c[2]-t/c[3])
+
+ # Computing LST in hours
+ lst = (theta + longitude)/15.0
+ neg = numpy.where(lst < 0.0)
+ if neg[0].size>0:lst[neg] = 24.0 + (lst[neg] % 24)
+ lst = lst % 24.0
+
+ return lst
+
+
+class date2doy:
+ def __init__(self,year,month,day):
+ self.year = year
+ self.month = month
+ self.day = day
+
+ def change2doy(self):
+ if calendar.isleap(self.year) == True:
+ tfactor = 1
+ else:
+ tfactor = 2
+
+ day = self.day
+ month = self.month
+
+ doy = numpy.floor((275*month)/9.0) - (tfactor*numpy.floor((month+9)/12.0)) + day - 30
+
+ return numpy.int32(doy)
+
+
+class Doy2Date:
+ def __init__(self,year,doy):
+ self.year = year
+ self.doy = doy
+
+ def change2date(self):
+ months = numpy.arange(12) + 1
+
+ first_dem = date2doy(self.year,months,1)
+ first_dem = first_dem.change2doy()
+
+ imm = numpy.where((self.doy - first_dem) > 0)
+
+ month = imm[0].size
+ dom = self.doy -first_dem[month - 1] + 1
+
+ return month, dom
diff --git a/apps/abs/utils/__init__.py b/apps/abs/utils/__init__.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/apps/abs/utils/__init__.py
diff --git a/apps/abs/utils/col_koki.dat b/apps/abs/utils/col_koki.dat
new file mode 100755
index 0000000000000000000000000000000000000000..1b22c9d39f1893aef283867e94b9a5ab7ee78390
GIT binary patch
literal 768
zc${NkW?=Zw@Q>mD&i`Przuuk!8~6_g4Kp{Lc?w~Ee0=r9(xx~YNfCZdR;=0}K*l1p
zuy@t5NB{qSd;8?}ZDtY1EFYEN@@QD#E8zniUzwz8BU2SZP9U;l)OlO|7@I&J!l
znX_iknLBU(f`y9~FIl>5`HGdRRKIRx^4T8ox67L*}HH5frEz*A31vL
z_=%IJPM=S6+fWcx&6lTn`mg|3kcFoHp|e
znEUqj$?eNplARRMJ(R{y1R%{{{p&{Qn<vw=Za&C(~KxNYXzExQk&x^(m5ix1zC)C`dSGl2lR{{c04
B)}jCa
literal 0
Hc$@mD&i`Przuuk!8~6_g4Kp{Lc?w~Ee0=r9(xx~YNfCZdR;=0}K*l1p
zuy@t5NB{qSd;8?}ZDtY1EFYEN@@QD#E8zniUzwz8BU2SZP9U;l)OlO|7@I&J!l
znX_iknLBU(f`y9~FIl>5`HGdRRKIRx^4T8ox67L*}HH5frEz*A31vL
z_=%IJPM=S6+fWcx&6lTn`mg|3kcFoHp|e
znEUqj$?eNplARRMJ(R{y1R%{{{p&{Qn<vw=Za&C(~KxNYXzExQk&x^(m5ix1zC)C`dSGl2lR{{c04
B)}jCa
literal 0
Hc$@> [pattern, mean_pos] = JroPattern(pattern=2).getPattern()
+ >> print meanpos
+ [ 8.08728085e-14 -4.78193873e-14]
+
+ Modification history
+ --------------------
+ Developed by Jorge L. Chau.
+ Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009.
+ """
+
+ if (self.fftopt>0) and (self.getcut>0):
+ #print "Conflict bewteen fftopt and getcut"
+ #print "To get a cut of the antenna pattern uses ffopt=0"
+ return None, None
+
+ if (self.fftopt==0):
+ # Getting antenna pattern using the array method
+ self.pattern = self.__usingArray(rx=1)
+ if (self.justrx==0):self.pattern = self.pattern*self.__usingArray(rx=0)
+
+ elif (self.fftopt>0):
+ # Getting antenna pattern using FFT method
+ self.pattern = self.__usingFFT(rx=1)
+ if (self.justrx==0):self.pattern = self.pattern*self.__usingFFT(rx=0)
+
+ self.maxpattern = numpy.nanmax(self.pattern)
+ self.norpattern = self.pattern/self.maxpattern
+ if self.getcut==0:self.__getBeamPars()
+
+ def __usingArray(self,rx):
+ """
+ __usingArray method returns the Jicamarca antenna pattern computed using array model
+
+ pattern = dipolepattern x modulepattern
+
+ Parameters
+ ----------
+ rx = Set to 1 to use the Rx information. Otherwise set to 0 for Tx.
+
+ Return
+ ------
+ pattern = An array giving the modelled antenna pattern using the array model.
+
+ Modification history
+ --------------------
+ Developed by Jorge L. Chau.
+ Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009.
+ """
+
+ if rx==1:
+ ues = self.uesrx
+ phase = self.phaserx
+ gain = self.gainrx
+ elif rx==0:
+ ues = self.uestx
+ phase = self.phasetx
+ gain = self.gaintx
+
+ ues = ues*360./self.airwl
+ phase = phase*360./self.airwl
+
+ for ii in range(4):
+ if ii==0:dim = numpy.array([4,0,8,4]) # WEST
+ elif ii==1:dim = numpy.array([0,0,4,4]) # NORTH
+ elif ii==2:dim = numpy.array([0,4,4,8]) # EAST
+ elif ii==3:dim = numpy.array([4,4,8,8]) # SOUTH
+ xi = dim[0]; xf = dim[2]; yi = dim[1]; yf = dim[3]
+ phase[xi:xf,yi:yf] = phase[xi:xf,yi:yf] + ues[ii]
+
+ phase = -phase
+
+ ar = self.eomwl*numpy.array([[0.5,6., 24.5],[0.5,6.,24.5]])
+ nr = numpy.array([[12.,4.,2.],[12.,4.,2.]])
+ lr = 0.25*self.eomwl*numpy.array([[0,0.,0],[0.,0,0]])
+
+ # Computing module and dipole patterns.
+ pattern = (numpy.abs(self.__dipPattern(ar,nr,lr)*self.__modPattern(phase,gain)))**2
+
+ return pattern
+
+ def __usingFFT(self,rx):
+ """
+ __usingFFT method returns the Jicamarca antenna pattern computed using The Fast Fou-
+ rier Transform.
+
+ pattern = iFFT(FFT(gain*EXP(j*phase)))
+
+ Parameters
+ ----------
+ rx = Set to 1 to use the Rx information. Otherwise set to 0 for Tx.
+
+ Return
+ ------
+ pattern = An array giving the modelled antenna pattern using the array model.
+
+ Modification history
+ --------------------
+ Developed by Jorge L. Chau.
+ Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009.
+ """
+
+ if rx==1:
+ ues = self.uesrx
+ phase = self.phaserx
+ gain = self.gainrx
+ elif rx==0:
+ ues = self.uestx
+ phase = self.phasetx
+ gain = self.gaintx
+
+ ues = ues*360./self.airwl
+ phase = phase*360./self.airwl
+
+ for ii in range(4):
+ if ii==0:dim = numpy.array([4,0,8,4]) # WEST
+ elif ii==1:dim = numpy.array([0,0,4,4]) # NORTH
+ elif ii==2:dim = numpy.array([0,4,4,8]) # EAST
+ elif ii==3:dim = numpy.array([4,4,8,8]) # SOUTH
+ xi = dim[0]; xf = dim[2]; yi = dim[1]; yf = dim[3]
+ phase[xi:xf,yi:yf] = phase[xi:xf,yi:yf] + ues[ii]
+
+ phase = -phase
+
+ delta_x = self.eomwl/2.
+ delta_y = self.eomwl/2.
+
+ nxfft = 2048
+ nyfft = 2048
+ dcosx = (numpy.arange(nxfft) - (0.5*nxfft))/(nxfft*delta_x)*self.eomwl
+ dcosy = (numpy.arange(nyfft) - (0.5*nyfft))/(nyfft*delta_y)*self.eomwl
+
+ fft_gain = numpy.zeros((nxfft,nyfft))
+ fft_phase = numpy.zeros((nxfft,nyfft))
+
+ nx = 8
+ ny = 8
+ ndx =12
+ ndy =12
+ for iy in numpy.arange(ny):
+ for ix in numpy.arange(nx):
+ ix1 = nxfft/2-self.nx/2*ndx+ix*ndx
+ if ix<(nx/2):ix1 = ix1 - 1
+ if ix>=(nx/2):ix1 = ix1 + 1
+
+ iy1 = nyfft/2-ny/2*ndx+iy*ndy
+ if iy<(ny/2):iy1 = iy1 - 1
+ if iy>=(ny/2):iy1 = iy1 + 1
+
+ fft_gain[ix1:ix1+ndx-1,iy1:iy1+ndy-1] = gain[ix,ny-1-iy]
+ fft_phase[ix1:ix1+ndx-1,iy1:iy1+ndy-1] = phase[ix,ny-1-iy]
+
+
+ fft_phase = fft_phase*Misc_Routines.CoFactors.d2r
+
+ pattern = numpy.abs(numpy.fft.fft2(fft_gain*numpy.exp(numpy.complex(0,1)*fft_phase)))**2
+ pattern = numpy.fft.fftshift(pattern)
+
+ xvals = numpy.where((dcosx>=(numpy.min(self.dcosx))) & (dcosx<=(numpy.max(self.dcosx))))
+ yvals = numpy.where((dcosy>=(numpy.min(self.dcosy))) & (dcosy<=(numpy.max(self.dcosy))))
+
+ pattern = pattern[xvals[0][0]:xvals[0][-1],yvals[0][0]:yvals[0][-1]]
+
+ return pattern
+
+ def __readAttenuation(self):
+ """
+ _readAttenuation reads the attenuations' file and returns an array giving these va-
+ lues (dB). The ext file must be in the directory "resource".
+
+ Return
+ ------
+ attenuation = An array giving attenuation values read from the text file.
+
+ Modification history
+ --------------------
+ Developed by Jorge L. Chau.
+ Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009.
+ """
+
+ attenuation = None
+# foldr = sys.path[-1] + os.sep + "resource" + os.sep
+ base_path = os.path.dirname(os.path.abspath(__file__))
+ #foldr = './resource'
+ #filen = "attenuation.txt"
+ attenuationFile = os.path.join(base_path,"resource","attenuation.txt")
+ #ff = open(os.path.join(foldr,filen),'r')
+ ff = open(attenuationFile,'r')
+ exec(ff.read())
+ ff.close()
+
+ return attenuation
+
+ def __dipPattern(self,ar,nr,lr):
+ """
+ _dipPattern function computes the dipole's pattern to the Jicamarca radar. The next
+ equation defines the pattern as a function of the mainlobe direction:
+
+ sincx = SIN(k/2*n0x*(a0x*SIN(phi)*COS(alpha)))/SIN(k/2*(a0x*SIN(phi)*COS(alpha)))
+ sincy = SIN(k/2*n0y*(a0y*SIN(phi)*SIN(alpha)))/SIN(k/2*(a0y*SIN(phi)*SIN(alpha)))
+ A0(phi,alpha) = sincx*sincy
+ Parameters
+ ----------
+ ar = ?
+ nr = ?
+ lr = ?
+
+ Return
+ ------
+ dipole = An array giving antenna pattern from the dipole point of view..
+
+ Modification history
+ --------------------
+ Developed by Jorge L. Chau.
+ Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009.
+ """
+
+ dipole = numpy.zeros((self.nx,self.ny),dtype=complex)
+ for iy in range(self.ny):
+ for ix in range(self.nx):
+ yindex = iy*(self.getcut==0) + ix*(self.getcut==1)
+
+ argx = ar[0,0]*self.dcosx[ix] - lr[0,0]
+ junkx = numpy.sin(0.5*self.kk*nr[0,0]*argx)/numpy.sin(0.5*self.kk*argx)
+ if argx == 0.0: junkx = nr[0,0]
+
+ argy = ar[1,0]*self.dcosy[yindex] - lr[1,0]
+ junky = numpy.sin(0.5*self.kk*nr[1,0]*argy)/numpy.sin(0.5*self.kk*argy)
+ if argy == 0.0: junky = nr[1,0]
+
+ dipole[ix,iy] = junkx*junky
+
+ return dipole
+
+ def __modPattern(self,phase,gain):
+ """
+ ModPattern computes the module's pattern to the Jicamarca radar. The next equation
+ defines the pattern as a function mainlobe direction:
+
+ phasex = pos(x)*SIN(phi)*COS(alpha)
+ phasey = pos(y)*SIN(phi)*SIN(alpha)
+
+ A1(phi,alpha) = TOTAL(gain*EXP(COMPLEX(0,k*(phasex+phasey)+phase)))
+
+ Parameters
+ ----------
+ phase = Bidimensional array (8x8) giving the phase (in meters) of each module.
+ gain = Bidimensional array (8x8) giving to define modules will be active (ones)
+ and which will not (zeros).
+
+ Return
+ ------
+ module = An array giving antenna pattern from the module point of view..
+
+ Modification history
+ --------------------
+ Developed by Jorge L. Chau.
+ Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009.
+ """
+
+ pos = self.eomwl*self.__readAttenuation()
+ posx = pos[0,:,:]
+ posy = pos[1,:,:]
+
+ phase = phase*Misc_Routines.CoFactors.d2r
+ module = numpy.zeros((self.nx,self.ny),dtype=complex)
+ for iy in range(self.ny):
+ for ix in range(self.nx):
+ yindex = iy*(self.getcut==0) + ix*(self.getcut==1)
+ phasex = posx*self.dcosx[ix]
+ phasey = posy*self.dcosy[yindex]
+ tmp = gain*numpy.exp(numpy.complex(0,1.)*(self.kk*(phasex+phasey)+phase))
+ module[ix,iy] = tmp.sum()
+
+ return module
+
+ def __getBeamPars(self):
+ """
+ _getBeamPars computes the main-beam parameters of the antenna.
+
+ Modification history
+ --------------------
+ Developed by Jorge L. Chau.
+ Converted to Python by Freddy R. Galindo, ROJ, 20 September 2009.
+ """
+
+ dx = self.dcosx[1] - self.dcosx[0]
+ dy = self.dcosy[1] - self.dcosy[0]
+
+ amp = self.norpattern
+
+ xx = numpy.resize(self.dcosx,(self.nx,self.nx)).transpose()
+ yy = numpy.resize(self.dcosy,(self.ny,self.ny))
+
+ mm0 = amp[numpy.where(amp > 0.5)]
+ xx0 = xx[numpy.where(amp > 0.5)]
+ yy0 = yy[numpy.where(amp > 0.5)]
+
+ xc = numpy.sum(mm0*xx0)/numpy.sum(mm0)
+ yc = numpy.sum(mm0*yy0)/numpy.sum(mm0)
+ rc = numpy.sqrt(mm0.size*dx*dy/numpy.pi)
+
+ nnx = numpy.where(numpy.abs(self.dcosx - xc) < rc)
+ nny = numpy.where(numpy.abs(self.dcosy - yc) < rc)
+
+ mm1 = amp[numpy.min(nnx):numpy.max(nnx)+1,numpy.min(nny):numpy.max(nny)+1]
+ xx1 = self.dcosx[numpy.min(nnx):numpy.max(nnx)+1]
+ yy1 = self.dcosy[numpy.min(nny):numpy.max(nny)+1]
+
+ # fitting data into the main beam.
+ import gaussfit
+ params = gaussfit.fitgaussian(mm1)
+
+ # Tranforming from indexes to axis' values
+ xcenter = xx1[0] + (((xx1[xx1.size-1] - xx1[0])/(xx1.size -1))*(params[1]))
+ ycenter = yy1[0] + (((yy1[yy1.size-1] - yy1[0])/(yy1.size -1))*(params[2]))
+ xwidth = ((xx1[xx1.size-1] - xx1[0])/(xx1.size-1))*(params[3])*(1/Misc_Routines.CoFactors.d2r)
+ ywidth = ((yy1[yy1.size-1] - yy1[0])/(yy1.size-1))*(params[4])*(1/Misc_Routines.CoFactors.d2r)
+ meanwx = (xwidth*ywidth)
+ meanpos = numpy.array([xcenter,ycenter])
+
+ #print 'Position: %f %f' %(xcenter,ycenter)
+ #print 'Widths: %f %f' %(xwidth, ywidth)
+ #print 'BWHP: %f' %(2*numpy.sqrt(2*meanwx)*numpy.sqrt(-numpy.log(0.5)))
+
+ self.meanpos = meanpos
+
+
+class BField():
+ def __init__(self,year=None,doy=None,site=1,heights=None,alpha_i=90):
+ """
+ BField class creates an object to get the Magnetic field for a specific date and
+ height(s).
+
+ Parameters
+ ----------
+ year = A scalar giving the desired year. If the value is None (default value) then
+ the current year will be used.
+ doy = A scalar giving the desired day of the year. If the value is None (default va-
+ lue) then the current doy will be used.
+ site = An integer to choose the geographic coordinates of the place where the magne-
+ tic field will be computed. The default value is over Jicamarca (site=1)
+ 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.
+ alpha_i = Angle to interpolate the magnetic field.
+
+ Modification History
+ --------------------
+ Converted to Object-oriented Programming by Freddy Galindo, ROJ, 07 October 2009.
+ """
+
+ tmp = time.localtime()
+ if year==None: year = tmp[0]
+ if doy==None: doy = tmp[7]
+ self.year = year
+ self.doy = doy
+ self.site = site
+ if heights==None:heights = numpy.array([100,500,1000])
+ self.heights = heights
+ self.alpha_i = alpha_i
+
+ def getBField(self,maglimits=numpy.array([-7,-7,7,7])):
+ """
+ getBField models the magnetic field for a different heights in a specific date.
+
+ Parameters
+ ----------
+ maglimits = An 4-elements array giving ..... The default value is [-7,-7,7,7].
+
+ Return
+ ------
+ dcos = An 4-dimensional array giving the directional cosines of the magnetic field
+ over the desired place.
+ alpha = An 3-dimensional array giving the angle of the magnetic field over the desi-
+ red place.
+
+ Modification History
+ --------------------
+ Converted to Python by Freddy R. Galindo, ROJ, 07 October 2009.
+ """
+
+ x_ant = numpy.array([1,0,0])
+ y_ant = numpy.array([0,1,0])
+ z_ant = numpy.array([0,0,1])
+
+ if self.site==0:
+ title_site = "Magnetic equator"
+ coord_site = numpy.array([-76+52./60.,-11+57/60.,0.5])
+ elif self.site==1:
+ title_site = 'Jicamarca'
+ coord_site = [-76-52./60.,-11-57/60.,0.5]
+ theta = (45+5.35)*numpy.pi/180. # (50.35 and 1.46 from Fleish Thesis)
+ delta = -1.46*numpy.pi/180
+
+ x_ant1 = numpy.roll(self.rotvector(self.rotvector(x_ant,1,delta),3,theta),1)
+ y_ant1 = numpy.roll(self.rotvector(self.rotvector(y_ant,1,delta),3,theta),1)
+ z_ant1 = numpy.roll(self.rotvector(self.rotvector(z_ant,1,delta),3,theta),1)
+
+ ang0 = -1*coord_site[0]*numpy.pi/180.
+ ang1 = coord_site[1]*numpy.pi/180.
+ x_ant = self.rotvector(self.rotvector(x_ant1,2,ang1),3,ang0)
+ y_ant = self.rotvector(self.rotvector(y_ant1,2,ang1),3,ang0)
+ z_ant = self.rotvector(self.rotvector(z_ant1,2,ang1),3,ang0)
+ else:
+# print "No defined Site. Skip..."
+ return None
+
+ nhei = self.heights.size
+ pt_intercep = numpy.zeros((nhei,2))
+ nfields = 1
+
+ grid_res = 0.5
+ nlon = numpy.int(maglimits[2] - maglimits[0])/grid_res + 1
+ nlat = numpy.int(maglimits[3] - maglimits[1])/grid_res + 1
+
+ location = numpy.zeros((nlon,nlat,2))
+ mlon = numpy.atleast_2d(numpy.arange(nlon)*grid_res + maglimits[0])
+ mrep = numpy.atleast_2d(numpy.zeros(nlat) + 1)
+ location0 = numpy.dot(mlon.transpose(),mrep)
+
+ mlat = numpy.atleast_2d(numpy.arange(nlat)*grid_res + maglimits[1])
+ mrep = numpy.atleast_2d(numpy.zeros(nlon) + 1)
+ location1 = numpy.dot(mrep.transpose(),mlat)
+
+ location[:,:,0] = location0
+ location[:,:,1] = location1
+
+ alpha = numpy.zeros((nlon,nlat,nhei))
+ rr = numpy.zeros((nlon,nlat,nhei,3))
+ dcos = numpy.zeros((nlon,nlat,nhei,2))
+
+ global first_time
+
+ first_time = None
+ for ilon in numpy.arange(nlon):
+ for ilat in numpy.arange(nlat):
+ outs = self.__bdotk(self.heights,
+ self.year + self.doy/366.,
+ coord_site[1],
+ coord_site[0],
+ coord_site[2],
+ coord_site[1]+location[ilon,ilat,1],
+ location[ilon,ilat,0]*720./180.)
+
+ alpha[ilon, ilat,:] = outs[1]
+ rr[ilon, ilat,:,:] = outs[3]
+
+ mrep = numpy.atleast_2d((numpy.zeros(nhei)+1)).transpose()
+ tmp = outs[3]*numpy.dot(mrep,numpy.atleast_2d(x_ant))
+ tmp = tmp.sum(axis=1)
+ dcos[ilon,ilat,:,0] = tmp/numpy.sqrt((outs[3]**2).sum(axis=1))
+
+ mrep = numpy.atleast_2d((numpy.zeros(nhei)+1)).transpose()
+ tmp = outs[3]*numpy.dot(mrep,numpy.atleast_2d(y_ant))
+ tmp = tmp.sum(axis=1)
+ dcos[ilon,ilat,:,1] = tmp/numpy.sqrt((outs[3]**2).sum(axis=1))
+
+ return dcos, alpha, nlon, nlat
+
+
+ def __bdotk(self,heights,tm,gdlat=-11.95,gdlon=-76.8667,gdalt=0.0,decd=-12.88, ham=-4.61666667):
+
+ global first_time
+ # Mean Earth radius in Km WGS 84
+ a_igrf = 6371.2
+
+ bk = numpy.zeros(heights.size)
+ alpha = numpy.zeros(heights.size)
+ bfm = numpy.zeros(heights.size)
+ rr = numpy.zeros((heights.size,3))
+ rgc = numpy.zeros((heights.size,3))
+
+ ObjGeodetic = Astro_Coords.Geodetic(gdlat,gdalt)
+ [gclat,gcalt] = ObjGeodetic.change2geocentric()
+
+ gclat = gclat*numpy.pi/180.
+ gclon = gdlon*numpy.pi/180.
+
+ # Antenna position from center of Earth
+ ca_vector = [numpy.cos(gclat)*numpy.cos(gclon),numpy.cos(gclat)*numpy.sin(gclon),numpy.sin(gclat)]
+ ca_vector = gcalt*numpy.array(ca_vector)
+
+ dec = decd*numpy.pi/180.
+
+ # K vector respect to the center of earth.
+ klon = gclon + ham*numpy.pi/720.
+ k_vector = [numpy.cos(dec)*numpy.cos(klon),numpy.cos(dec)*numpy.sin(klon),numpy.sin(dec)]
+ k_vector = numpy.array(k_vector)
+
+ for ih in numpy.arange(heights.size):
+ # Vector from Earth's center to volume of interest
+ rr[ih,:] = k_vector*heights[ih]
+ cv_vector = numpy.squeeze(ca_vector) + rr[ih,:]
+
+ cv_gcalt = numpy.sqrt(numpy.sum(cv_vector**2.))
+ cvxy = numpy.sqrt(numpy.sum(cv_vector[0:2]**2.))
+
+ radial = cv_vector/cv_gcalt
+ east = numpy.array([-1*cv_vector[1],cv_vector[0],0])/cvxy
+ comp1 = east[1]*radial[2] - radial[1]*east[2]
+ comp2 = east[2]*radial[0] - radial[2]*east[0]
+ comp3 = east[0]*radial[1] - radial[0]*east[1]
+ north = -1*numpy.array([comp1, comp2, comp3])
+
+ rr_k = cv_vector - numpy.squeeze(ca_vector)
+ u_rr = rr_k/numpy.sqrt(numpy.sum(rr_k**2.))
+
+ cv_gclat = numpy.arctan2(cv_vector[2],cvxy)
+ cv_gclon = numpy.arctan2(cv_vector[1],cv_vector[0])
+
+ bhei = cv_gcalt-a_igrf
+ blat = cv_gclat*180./numpy.pi
+ blon = cv_gclon*180./numpy.pi
+ bfield = self.__igrfkudeki(bhei,tm,blat,blon)
+
+ B = (bfield[0]*north + bfield[1]*east - bfield[2]*radial)*1.0e-5
+
+ bfm[ih] = numpy.sqrt(numpy.sum(B**2.)) #module
+ bk[ih] = numpy.sum(u_rr*B)
+ alpha[ih] = numpy.arccos(bk[ih]/bfm[ih])*180/numpy.pi
+ rgc[ih,:] = numpy.array([cv_gclon, cv_gclat, cv_gcalt])
+
+ return bk, alpha, bfm, rr, rgc
+
+
+ def __igrfkudeki(self,heights,time,latitude,longitude,ae=6371.2):
+ """
+ __igrfkudeki calculates the International Geomagnetic Reference Field for given in-
+ put conditions based on IGRF2005 coefficients.
+
+ Parameters
+ ----------
+ heights = Scalar or vector giving the height above the Earth of the point in ques-
+ tion in kilometers.
+ time = Scalar or vector giving the decimal year of time in question (e.g. 1991.2).
+ latitude = Latitude of point in question in decimal degrees. Scalar or vector.
+ longitude = Longitude of point in question in decimal degrees. Scalar or vector.
+ ae =
+ first_time =
+
+ Return
+ ------
+ bn =
+ be =
+ bd =
+ bmod =
+ balpha =
+ first_time =
+
+ Modification History
+ --------------------
+ Converted to Python by Freddy R. Galindo, ROJ, 03 October 2009.
+ """
+
+ global first_time
+ global gs, hs, nvec, mvec, maxcoef
+
+ heights = numpy.atleast_1d(heights)
+ time = numpy.atleast_1d(time)
+ latitude = numpy.atleast_1d(latitude)
+ longitude = numpy.atleast_1d(longitude)
+
+ if numpy.max(latitude)==90:
+# print "Field calculations are not supported at geographic poles"
+ pass
+
+ # output arrays
+ bn = numpy.zeros(heights.size)
+ be = numpy.zeros(heights.size)
+ bd = numpy.zeros(heights.size)
+
+ if first_time==None:first_time=0
+
+ time0 = time[0]
+ if time!=first_time:
+ #print "Getting coefficients for", time0
+ [periods,g,h ] = self.__readIGRFcoeff()
+ top_year = numpy.max(periods)
+ nperiod = (top_year - 1900)/5 + 1
+
+ maxcoef = 10
+ if time0>=2000:maxcoef = 12
+
+
+ # Normalization array for Schmidt fucntions
+ multer = numpy.zeros((2+maxcoef,1+maxcoef)) + 1
+ for cn in (numpy.arange(maxcoef)+1):
+ for rm in (numpy.arange(cn)+1):
+ tmp = numpy.arange(2*rm) + cn - rm + 1.
+ multer[rm+1,cn] = ((-1.)**rm)*numpy.sqrt(2./tmp.prod())
+
+ schmidt = multer[1:,1:].transpose()
+
+ # n and m arrays
+ nvec = numpy.atleast_2d(numpy.arange(maxcoef)+2)
+ mvec = numpy.atleast_2d(numpy.arange(maxcoef+1)).transpose()
+
+ # Time adjusted igrf g and h with Schmidt normalization
+ # IGRF coefficient arrays: g0(n,m), n=1, maxcoeff,m=0, maxcoeff, ...
+ if time0 top_year
+ dtime = (time0 - top_year) + 5
+ ntime = g[:,0,0].size - 2
+
+ g0 = g[ntime,1:maxcoef+1,:maxcoef+1]
+ h0 = h[ntime,1:maxcoef+1,:maxcoef+1]
+ gdot = g[ntime+1,1:maxcoef+1,:maxcoef+1]-g[ntime,1:maxcoef+1,:maxcoef+1]
+ hdot = h[ntime+1,1:maxcoef+1,:maxcoef+1]-h[ntime,1:maxcoef+1,:maxcoef+1]
+ gs = (g0 + dtime*(gdot/5.))*schmidt[:maxcoef,0:maxcoef+1]
+ hs = (h0 + dtime*(hdot/5.))*schmidt[:maxcoef,0:maxcoef+1]
+
+ first_time = time0
+
+ for ii in numpy.arange(heights.size):
+ # Height dependence array rad = (ae/(ae+height))**(n+3)
+ rad = numpy.atleast_2d((ae/(ae + heights[ii]))**(nvec+1))
+
+ # Sin and Cos of m times longitude phi arrays
+ mphi = mvec*longitude[ii]*numpy.pi/180.
+ cosmphi = numpy.atleast_2d(numpy.cos(mphi))
+ sinmphi = numpy.atleast_2d(numpy.sin(mphi))
+
+ # Cos of colatitude theta
+ c = numpy.cos((90 - latitude[ii])*numpy.pi/180.)
+
+ # Legendre functions p(n,m|c)
+ [p,dp]= scipy.special.lpmn(maxcoef+1,maxcoef+1,c)
+ p = p[:,:-1].transpose()
+ s = numpy.sqrt((1. - c)*(1 + c))
+
+ # Generate derivative array dpdtheta = -s*dpdc
+ dpdtheta = c*p/s
+ for m in numpy.arange(maxcoef+2): dpdtheta[:,m] = m*dpdtheta[:,m]
+ dpdtheta = dpdtheta + numpy.roll(p,-1,axis=1)
+
+ # Extracting arrays required for field calculations
+ p = p[1:maxcoef+1,:maxcoef+1]
+ dpdtheta = dpdtheta[1:maxcoef+1,:maxcoef+1]
+
+ # Weigh p and dpdtheta with gs and hs coefficients.
+ gp = gs*p
+ hp = hs*p
+ gdpdtheta = gs*dpdtheta
+ hdpdtheta = hs*dpdtheta
+ # Calcultate field components
+ matrix0 = numpy.dot(gdpdtheta,cosmphi)
+ matrix1 = numpy.dot(hdpdtheta,sinmphi)
+ bn[ii] = numpy.dot(rad,(matrix0 + matrix1))
+ matrix0 = numpy.dot(hp,(mvec*cosmphi))
+ matrix1 = numpy.dot(gp,(mvec*sinmphi))
+ be[ii] = numpy.dot((-1*rad),((matrix0 - matrix1)/s))
+ matrix0 = numpy.dot(gp,cosmphi)
+ matrix1 = numpy.dot(hp,sinmphi)
+ bd[ii] = numpy.dot((-1*nvec*rad),(matrix0 + matrix1))
+
+ bmod = numpy.sqrt(bn**2. + be**2. + bd**2.)
+ btheta = numpy.arctan(bd/numpy.sqrt(be**2. + bn**2.))*180/numpy.pi
+ balpha = numpy.arctan(be/bn)*180./numpy.pi
+
+ #bn : north
+ #be : east
+ #bn : radial
+ #bmod : module
+
+
+ return bn, be, bd, bmod, btheta, balpha
+
+ def str2num(self, datum):
+ try:
+ return int(datum)
+ except:
+ try:
+ return float(datum)
+ except:
+ return datum
+
+ def __readIGRFfile(self, filename):
+ list_years=[]
+ for i in range(1,24):
+ list_years.append(1895.0 + i*5)
+
+ epochs=list_years
+ epochs.append(epochs[-1]+5)
+ nepochs = numpy.shape(epochs)
+
+ gg = numpy.zeros((13,14,nepochs[0]),dtype=float)
+ hh = numpy.zeros((13,14,nepochs[0]),dtype=float)
+
+ coeffs_file=open(filename)
+ lines=coeffs_file.readlines()
+
+ coeffs_file.close()
+
+ for line in lines:
+ items = line.split()
+ g_h = items[0]
+ n = self.str2num(items[1])
+ m = self.str2num(items[2])
+
+ coeffs = items[3:]
+
+ for i in range(len(coeffs)-1):
+ coeffs[i] = self.str2num(coeffs[i])
+
+ #coeffs = numpy.array(coeffs)
+ ncoeffs = numpy.shape(coeffs)[0]
+
+ if g_h == 'g':
+ # print n," g ",m
+ gg[n-1,m,:]=coeffs
+ elif g_h=='h':
+ # print n," h ",m
+ hh[n-1,m,:]=coeffs
+ # else :
+ # continue
+
+ # Ultimo Reordenamiento para almacenar .
+ gg[:,:,nepochs[0]-1] = gg[:,:,nepochs[0]-2] + 5*gg[:,:,nepochs[0]-1]
+ hh[:,:,nepochs[0]-1] = hh[:,:,nepochs[0]-2] + 5*hh[:,:,nepochs[0]-1]
+
+# return numpy.array([gg,hh])
+ periods = numpy.array(epochs)
+ g = gg
+ h = hh
+ return periods, g, h
+
+
+ def __readIGRFcoeff(self,filename="igrf10coeffs.dat"):
+ """
+ __readIGRFcoeff reads the coefficients from a binary file which is located in the
+ folder "resource."
+
+ Parameter
+ ---------
+ filename = A string to specify the name of the file which contains thec coeffs. The
+ default value is "igrf10coeffs.dat"
+
+ Return
+ ------
+ periods = A lineal array giving...
+ g1 =
+ h1 =
+
+ Modification History
+ --------------------
+ Converted to Python by Freddy R. Galindo, ROJ, 03 October 2009.
+ """
+
+# # igrfile = sys.path[-1] + os.sep + "resource" + os.sep + filename
+# igrfile = os.path.join('./resource',filename)
+# f = open(igrfile,'rb')
+# #f = open(os.getcwd() + os.sep + "resource" + os.sep + filename,'rb')
+#
+# # Reading SkyNoise Power (lineal scale)
+# periods = numpy.fromfile(f,numpy.dtype([('var',' Around "x"
+ axis = 2 -> Around "y"
+ axis = 3 -> Around "z"
+ ang = Angle of rotation (in radians). The default value is zero.
+
+ Return
+ ------
+ rotvector = A lineal array of 3 elements giving the new coordinates.
+
+ Modification History
+ --------------------
+ Converted to Python by Freddy R. Galindo, ROJ, 01 October 2009.
+ """
+
+ if axis==1:
+ t = [[1,0,0],[0,numpy.cos(ang),numpy.sin(ang)],[0,-numpy.sin(ang),numpy.cos(ang)]]
+ elif axis==2:
+ t = [[numpy.cos(ang),0,-numpy.sin(ang)],[0,1,0],[numpy.sin(ang),0,numpy.cos(ang)]]
+ elif axis==3:
+ t = [[numpy.cos(ang),numpy.sin(ang),0],[-numpy.sin(ang),numpy.cos(ang),0],[0,0,1]]
+
+ rotvector = numpy.array(numpy.dot(numpy.array(t),numpy.array(vector)))
+
+ return rotvector
+
+
+class overJroShow:
+
+# __serverdocspath = '/usr/local/www/htdocs'
+# __tmpDir = 'overJro/tempReports'
+# __serverdocspath = '/Users/dsuarez/Pictures'
+# __tmpDir = 'overjro'
+ __serverdocspath = None
+ __tmpDir = None
+
+ def __init__(self):
+ self.year = None
+ self.month = None
+ self.dom = None
+ self.pattern = None
+ self.maxphi = None
+ self.heights = None
+ self.filename = None
+ self.showType = None
+ self.path = None
+ self.objects = None
+ self.nptsx = 101
+ self.nptsy = 101
+ self.fftopt = 0
+ self.site = 1
+ self.dcosx = 1
+ self.dcosy = 1
+ self.dcosxrange = None
+ self.dcosyrange = None
+ self.maxha_min= 0.
+ self.show_object = None
+ self.dcosx_mag = None
+ self.dcosy_mag = None
+ self.ha_mag = None
+ self.time_mag = None
+ self.main_dec = None
+ self.ObjC = None
+ self.ptitle = ''
+ self.path4plotname = None
+ self.plotname0 = None
+ self.plotname1 = None
+ self.plotname2 = None
+ self.scriptHeaders = 0
+# self.outputHead('Show Plot')
+# self.printBody()
+
+ def setScriptState(self):
+ self.madForm = cgi.FieldStorage()
+
+ if self.madForm.has_key('serverdocspath'):
+ self.__serverdocspath = self.madForm.getvalue('serverdocspath')#'/usr/local/www/htdocs'
+
+ if self.madForm.has_key('tmpdir'):
+ self.__tmpDir = self.madForm.getvalue('tmpdir')#'overJro/tempReports'
+
+ if self.madForm.has_key('showType'):
+ self.showType = int(self.madForm.getvalue('showType'))
+
+ if self.showType == 0 or self.showType == 1:
+
+# if self.madForm.has_key('year') and \
+# self.madForm.has_key('month') and \
+# self.madForm.has_key('dom') and \
+# self.madForm.has_key('pattern') and \
+# self.madForm.has_key('maxphi') and \
+# self.madForm.has_key('objects') and \
+# self.madForm.has_key('heights'):
+
+ if self.madForm.has_key('year') and \
+ self.madForm.has_key('month') and \
+ self.madForm.has_key('dom') and \
+ self.madForm.has_key('maxphi') and \
+ self.madForm.has_key('objects') and \
+ self.madForm.has_key('heights'):
+
+ self.year = int(self.madForm.getvalue('year'))
+ self.month = int(self.madForm.getvalue('month'))
+ self.dom = int(self.madForm.getvalue('dom'))
+ self.maxphi = float(self.madForm.getvalue('maxphi'))
+
+ if self.madForm.has_key('pattern'):
+
+ tmp_pattern = self.madForm.getvalue('pattern') #pattern es predifinido en listado o definido por el usuario
+ self.pattern=[]
+ if tmp_pattern[0] == '[':
+ tmp_pattern=tmp_pattern[1:]
+
+ if tmp_pattern[-1] == ']':
+ tmp_pattern=tmp_pattern[0:len(tmp_pattern)-1]
+
+ for s in tmp_pattern.split(','):
+ self.pattern.append(float(s))
+ elif self.madForm.has_key('filename'):
+ if self.madForm.has_key('filename'):
+ self.filename = self.madForm.getvalue('filename') # nombre de archivo: patron de radiacion definido por el usuario
+
+ if self.madForm.has_key('path'):
+ self.path = self.madForm.getvalue('path') #path donde se encuentra el archivo: patron de radiacion del usuario
+
+ else:
+ print "Content-Type: text/html\n"
+ print '
This cgi plot script was called without the proper arguments.
'
+ print '
This is a script used to plot Antenna Cuts over Jicamarca Antenna
'
+ print '
Required arguments:
'
+ print '
pattern - chekbox indicating objects over jicamarca antenna
'
+ print '
or'
+ print '
filename - The pattern defined by users is a file text'
+ print '
path - folder with pattern files'
+ sys.exit(0)
+
+
+ tmp_heights = self.madForm.getvalue('heights')
+ self.heights=[]
+ if tmp_heights[0] == '[':
+ tmp_heights=tmp_heights[1:]
+
+ if tmp_heights[-1] == ']':
+ tmp_heights=tmp_heights[0:len(tmp_heights)-1]
+
+ for s in tmp_heights.split(','):
+ self.heights.append(float(s))
+ self.heights = numpy.array(self.heights)
+
+ tmp_objects = self.madForm.getvalue('objects') #lista con los objetos a graficar en el patron de radiacion
+ self.objects=[]
+ if tmp_objects[0] == '[':
+ tmp_objects=tmp_objects[1:]
+
+ if tmp_objects[-1] == ']':
+ tmp_objects=tmp_objects[0:len(tmp_objects)-1]
+
+ for s in tmp_objects.split(','):
+ self.objects.append(int(s))
+
+ if self.showType == 1:
+ if numpy.sum(self.objects) == 0:
+ if self.scriptHeaders == 0:
+ print "Content-Type: text/html\n"
+ print '
This cgi plot script was called without the proper arguments.
'
+ print '
This is a script used to plot Antenna Cuts over Jicamarca Antenna
'
+ print '
Required arguments:
'
+ print '
objects - chekbox indicating objects over jicamarca antenna
'
+ print '
Please, options in "Select Object" must be checked'
+ sys.exit(0)
+
+ #considerar para futura implementacion
+ if self.madForm.has_key('filename'):
+ self.filename = self.madForm.getvalue('filename') # nombre de archivo: patron de radiacion definido por el usuario
+
+ if self.madForm.has_key('path'):
+ self.path = self.madForm.getvalue('path') #path donde se encuentra el archivo: patron de radiacion del usuario
+
+
+ else:
+ if self.scriptHeaders == 0:
+ print "Content-Type: text/html\n"
+
+ print '
This cgi plot script was called without the proper arguments.
'
+ print '
This is a script used to plot Pattern Field and Celestial Objects over Jicamarca Antenna
'
+ print '
Required arguments:
'
+ print '
year - year of event
'
+ print '
month - month of event
'
+ print '
dom - day of month
'
+ print '
pattern - pattern is defined by "Select an Experiment" list box
'
+ print '
maxphi - maxphi is defined by "Max Angle" text box
'
+ print '
objects - objects is a list defined by checkbox in "Select Object"
'
+ print '
heights - heights is defined by "Heights" text box, for default heights=[100,500,1000]
'
+ print '
showType - showType is a hidden element for show plot of Pattern&Object or Antenna Cuts or Sky Noise
'
+
+ sys.exit(0)
+
+ if self.showType == 2:
+ if self.madForm.has_key('year') and \
+ self.madForm.has_key('month') and \
+ self.madForm.has_key('dom'):
+
+ self.year = int(self.madForm.getvalue('year'))
+ self.month = int(self.madForm.getvalue('month'))
+ self.dom = int(self.madForm.getvalue('dom'))
+
+ else:
+ if self.scriptHeaders == 0:
+ print "Content-Type: text/html\n"
+ print '
This cgi plot script was called without the proper arguments.
'
+ print '
This is a script used to plot Sky Noise over Jicamarca Antenna
+
+