|
|
"""
|
|
|
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','<f4')]),480*20)
|
|
|
ha_sky = ha_sky['var'].reshape(20,480).transpose()
|
|
|
|
|
|
dec_sky = numpy.fromfile(f,numpy.dtype([('var','<f4')]),480*20)
|
|
|
dec_sky = dec_sky['var'].reshape((20,480)).transpose()
|
|
|
|
|
|
tmp_sky = numpy.fromfile(f,numpy.dtype([('var','<f4')]),480*20)
|
|
|
tmp_sky = tmp_sky['var'].reshape((20,480)).transpose()
|
|
|
|
|
|
f.close()
|
|
|
|
|
|
nha = 480
|
|
|
tmp_cut = numpy.zeros(nha)
|
|
|
for iha in numpy.arange(nha):
|
|
|
tck = scipy.interpolate.splrep(dec_sky[iha,:],tmp_sky[iha,:],s=0)
|
|
|
tmp_cut[iha] = scipy.interpolate.splev(dec_cut,tck,der=0)
|
|
|
|
|
|
ptr = numpy.nanargmax(tmp_cut)
|
|
|
|
|
|
maxra = ha_sky[ptr,0]
|
|
|
ra = ha_sky[:,0]
|
|
|
|
|
|
return maxra, ra
|
|
|
|
|
|
def skyNoise(self,jd,ut=-5.0,longitude=-76.87,filename='galaxy.txt',filepath=None):
|
|
|
"""
|
|
|
hydrapos returns RA and Dec provided by Bill Coles (Oct 2003).
|
|
|
|
|
|
Parameters
|
|
|
----------
|
|
|
jd = The julian date of the day (and time), scalar or vector.
|
|
|
|
|
|
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.
|
|
|
"""
|
|
|
|
|
|
# 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)) & (lst<numpy.max(otime))); val = val[0]
|
|
|
tck = scipy.interpolate.interp1d(otime,opowr)
|
|
|
ipowr[val] = tck(lst[val])
|
|
|
|
|
|
# Extrapolating above maximum time data (23.75).
|
|
|
uval = numpy.where(lst>numpy.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
|
|
|
|
