In [500]:

    

import sympy as sym
import numpy as np
sym.init_printing(fontsize='5pt',use_latex='mathjax', wrap_line=True)
sym.init_session()
from sympy.physics.vector import vlatex
from sympy.printing import latex


    

    




IPython console for SymPy 1.0 (Python 2.7.6-64-bit) (ground types: python)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at http://docs.sympy.org/1.0/

In [501]:

    

#symbols
skew, slope, lat, lon = sym.symbols('skew slope lat lon')

#values
consts = dict(skew = 2.3755870374367263e-05,
              slope = 0.0034494653328734047,                                                          
              lat = np.radians(-35.370707),
              lon = np.radians(149.424702))


    

In [502]:

    

def rotation_matrix(angle,d):
    directions = {
            "x":[1,0,0],
            "y":[0,1,0],
            "z":[0,0,1]
            }
    direction = np.array(directions[d])
    sina = sym.sin(angle)
    cosa = sym.cos(angle)
    # rotation matrix around unit vector                                                              
    R = sym.Matrix([[cosa,0,0],[0,cosa,0],[0,0,cosa]])
    R += sym.Matrix(np.outer(direction, direction)) * (1 - cosa)
    direction = sym.Matrix(direction)
    direction *= sina
    R += sym.Matrix([[ 0,           -direction[2],  direction[1]],
                    [ direction[2], 0,             -direction[0]],
                    [-direction[1], direction[0],   0]])
    return R


    

In [503]:

    

#reflections
def reflection_matrix(d):
    m = {
    "x":[[-1,0,0],[0, 1,0],[0,0, 1]],
    "y":[[1, 0,0],[0,-1,0],[0,0, 1]],
    "z":[[1, 0,0],[0, 1,0],[1,0,-1]]
    }
    return sym.Matrix(m[d])


    

In [504]:

    

def position_vector(a,b):
    return sym.Matrix([[sym.cos(b)*sym.cos(a)],
                      [sym.cos(b)*sym.sin(a)],
                      [sym.sin(b)]])


    

In [525]:

    

# z-axis lies along telescope east to west
# using right-hand rule
# thumb is z
# index is x
# middle is y
telescope_to_nsew = (rotation_matrix(-skew,"x") * rotation_matrix(slope,"y"))
telescope_to_nsew


    

    
Out[525]:





$$\left[\begin{matrix}\cos{\left (slope \right )} & 0 & \sin{\left (slope \right )}\\- \sin{\left (skew \right )} \sin{\left (slope \right )} & \cos{\left (skew \right )} & \sin{\left (skew \right )} \cos{\left (slope \right )}\\- \sin{\left (slope \right )} \cos{\left (skew \right )} & - \sin{\left (skew \right )} & \cos{\left (skew \right )} \cos{\left (slope \right )}\end{matrix}\right]$$

In [506]:

    

#use reflecion to convert to left-hand rule
nsew_to_azel = rotation_matrix(-sym.pi/2,"x") * rotation_matrix(sym.pi/2,"y") * reflection_matrix("x")    
nsew_to_azel


    

    
Out[506]:





$$\left[\begin{matrix}0 & 0 & 1\\1 & 0 & 0\\0 & -1 & 0\end{matrix}\right]$$

In [507]:

    

telescope_to_azel = (telescope_to_nsew * nsew_to_azel)
telescope_to_azel


    

    
Out[507]:





$$\left[\begin{matrix}0 & - \sin{\left (slope \right )} & \cos{\left (slope \right )}\\\cos{\left (skew \right )} & - \sin{\left (skew \right )} \cos{\left (slope \right )} & - \sin{\left (skew \right )} \sin{\left (slope \right )}\\- \sin{\left (skew \right )} & - \cos{\left (skew \right )} \cos{\left (slope \right )} & - \sin{\left (slope \right )} \cos{\left (skew \right )}\end{matrix}\right]$$

In [508]:

    

azel_to_hadec = rotation_matrix(sym.pi/2-lat,"y") * rotation_matrix(sym.pi,"z")
azel_to_hadec


    

    
Out[508]:





$$\left[\begin{matrix}- \sin{\left (lat \right )} & 0 & \cos{\left (lat \right )}\\0 & -1 & 0\\\cos{\left (lat \right )} & 0 & \sin{\left (lat \right )}\end{matrix}\right]$$

In [550]:

    

telescope_to_hadec = telescope_to_azel * azel_to_hadec
telescope_to_hadec


    

    
Out[550]:





$$\left[\begin{matrix}\cos{\left (lat \right )} \cos{\left (slope \right )} & \sin{\left (slope \right )} & \sin{\left (lat \right )} \cos{\left (slope \right )}\\- \sin{\left (lat \right )} \cos{\left (skew \right )} - \sin{\left (skew \right )} \sin{\left (slope \right )} \cos{\left (lat \right )} & \sin{\left (skew \right )} \cos{\left (slope \right )} & - \sin{\left (lat \right )} \sin{\left (skew \right )} \sin{\left (slope \right )} + \cos{\left (lat \right )} \cos{\left (skew \right )}\\\sin{\left (lat \right )} \sin{\left (skew \right )} - \sin{\left (slope \right )} \cos{\left (lat \right )} \cos{\left (skew \right )} & \cos{\left (skew \right )} \cos{\left (slope \right )} & - \sin{\left (lat \right )} \sin{\left (slope \right )} \cos{\left (skew \right )} - \sin{\left (skew \right )} \cos{\left (lat \right )}\end{matrix}\right]$$

In [510]:

    

telescope_to_hadec.T


    

    
Out[510]:





$$\left[\begin{matrix}\cos{\left (lat \right )} \cos{\left (slope \right )} & - \sin{\left (lat \right )} \cos{\left (skew \right )} - \sin{\left (skew \right )} \sin{\left (slope \right )} \cos{\left (lat \right )} & \sin{\left (lat \right )} \sin{\left (skew \right )} - \sin{\left (slope \right )} \cos{\left (lat \right )} \cos{\left (skew \right )}\\\sin{\left (slope \right )} & \sin{\left (skew \right )} \cos{\left (slope \right )} & \cos{\left (skew \right )} \cos{\left (slope \right )}\\\sin{\left (lat \right )} \cos{\left (slope \right )} & - \sin{\left (lat \right )} \sin{\left (skew \right )} \sin{\left (slope \right )} + \cos{\left (lat \right )} \cos{\left (skew \right )} & - \sin{\left (lat \right )} \sin{\left (slope \right )} \cos{\left (skew \right )} - \sin{\left (skew \right )} \cos{\left (lat \right )}\end{matrix}\right]$$

In [511]:

    

HA,Dec = sym.symbols('HA Dec')
hadec_vec = position_vector(HA,Dec)
hadec_vec


    

    
Out[511]:





$$\left[\begin{matrix}\cos{\left (Dec \right )} \cos{\left (HA \right )}\\\sin{\left (HA \right )} \cos{\left (Dec \right )}\\\sin{\left (Dec \right )}\end{matrix}\right]$$

In [512]:

    

ns,ew = sym.symbols('NS EW')
nsew_vec = position_vector(ns,ew)
nsew_vec


    

    
Out[512]:





$$\left[\begin{matrix}\cos{\left (EW \right )} \cos{\left (NS \right )}\\\sin{\left (NS \right )} \cos{\left (EW \right )}\\\sin{\left (EW \right )}\end{matrix}\right]$$

In [537]:

    

az,el = sym.symbols('Az El')
azel_vec = position_vector(az,el)
azel_vec


    

    
Out[537]:





$$\left[\begin{matrix}\cos{\left (Az \right )} \cos{\left (El \right )}\\\sin{\left (Az \right )} \cos{\left (El \right )}\\\sin{\left (El \right )}\end{matrix}\right]$$

In [513]:

    

equations = sym.Eq(nsew_vec,(telescope_to_hadec * hadec_vec))
equations


    

    
Out[513]:





$$\left[\begin{matrix}\cos{\left (EW \right )} \cos{\left (NS \right )}\\\sin{\left (NS \right )} \cos{\left (EW \right )}\\\sin{\left (EW \right )}\end{matrix}\right] = \left[\begin{matrix}\sin{\left (Dec \right )} \sin{\left (lat \right )} \cos{\left (slope \right )} + \sin{\left (HA \right )} \sin{\left (slope \right )} \cos{\left (Dec \right )} + \cos{\left (Dec \right )} \cos{\left (HA \right )} \cos{\left (lat \right )} \cos{\left (slope \right )}\\\left(- \sin{\left (lat \right )} \cos{\left (skew \right )} - \sin{\left (skew \right )} \sin{\left (slope \right )} \cos{\left (lat \right )}\right) \cos{\left (Dec \right )} \cos{\left (HA \right )} + \left(- \sin{\left (lat \right )} \sin{\left (skew \right )} \sin{\left (slope \right )} + \cos{\left (lat \right )} \cos{\left (skew \right )}\right) \sin{\left (Dec \right )} + \sin{\left (HA \right )} \sin{\left (skew \right )} \cos{\left (Dec \right )} \cos{\left (slope \right )}\\\left(\sin{\left (lat \right )} \sin{\left (skew \right )} - \sin{\left (slope \right )} \cos{\left (lat \right )} \cos{\left (skew \right )}\right) \cos{\left (Dec \right )} \cos{\left (HA \right )} + \left(- \sin{\left (lat \right )} \sin{\left (slope \right )} \cos{\left (skew \right )} - \sin{\left (skew \right )} \cos{\left (lat \right )}\right) \sin{\left (Dec \right )} + \sin{\left (HA \right )} \cos{\left (Dec \right )} \cos{\left (skew \right )} \cos{\left (slope \right )}\end{matrix}\right]$$

In [514]:

    

values = consts.copy()
values.update({HA:1.0,Dec:0.8})
e = equations.subs(values)
e


    

    
Out[514]:





$$\left[\begin{matrix}\cos{\left (EW \right )} \cos{\left (NS \right )}\\\sin{\left (NS \right )} \cos{\left (EW \right )}\\\sin{\left (EW \right )}\end{matrix}\right] = \left[\begin{matrix}-0.106277123772499\\0.802866409308139\\0.586609496826723\end{matrix}\right]$$

In [515]:

    

# caluclate ew using full matrix
a,b,c = np.array((telescope_to_hadec * hadec_vec).subs(values)).astype('float64').squeeze()
calc_ew = np.arcsin(c)
calc_ns = np.arctan2(b,a)
print "(EB) NS, EW (radians) = ",calc_ns,", ",calc_ew


    

    




(EB) NS, EW (radians) =  1.70240331085 ,  0.626865982991

In [516]:

    

def jer_delay(ha_source,dec_source, tilt=slope, skew=skew, latitude=lat):
    from numpy import sin,cos
    return sin(ha_source)*cos(dec_source)*cos(tilt) + \
            cos(ha_source)*cos(dec_source)* (skew*sin(latitude)-sin(tilt)*cos(latitude))- \
            sin(dec_source)*(skew*cos(latitude)+sin(tilt)*sin(latitude))


    

In [517]:

    

# calculate ew using jer delay
calc_ew_jer = jer_delay(values[HA],values[Dec],tilt=values['slope'],skew=values['skew'],latitude=values['lat'])
print "(JER) EW (radians) = ",np.arcsin(calc_ew_jer)


    

    




(JER) EW (radians) =  0.626865983195

In [518]:

    

#difference between JER and full matrix
print "Offset between EB and JER of:",np.degrees(calc_ew - np.arcsin(calc_ew_jer))*3600.0,"arcseconds"


    

    




Offset between EB and JER of: -4.2158041114e-05 arcseconds

In [519]:

    

#JER equation
sym.Eq(sin(ew),sin(HA)*cos(Dec)*cos(slope) + \
cos(HA)*cos(Dec)* (skew*sin(lat)-sin(slope)*cos(lat))- \
sin(Dec)*(skew*cos(lat)+sin(slope)*sin(lat)))


    

    
Out[519]:





$$\sin{\left (EW \right )} = \left(skew \sin{\left (lat \right )} - \sin{\left (slope \right )} \cos{\left (lat \right )}\right) \cos{\left (Dec \right )} \cos{\left (HA \right )} - \left(skew \cos{\left (lat \right )} + \sin{\left (lat \right )} \sin{\left (slope \right )}\right) \sin{\left (Dec \right )} + \sin{\left (HA \right )} \cos{\left (Dec \right )} \cos{\left (slope \right )}$$

In [520]:

    

new_eq = sym.Eq(nsew_vec[2],(telescope_to_hadec * hadec_vec)[2])
new_eq


    

    
Out[520]:





$$\sin{\left (EW \right )} = \left(\sin{\left (lat \right )} \sin{\left (skew \right )} - \sin{\left (slope \right )} \cos{\left (lat \right )} \cos{\left (skew \right )}\right) \cos{\left (Dec \right )} \cos{\left (HA \right )} + \left(- \sin{\left (lat \right )} \sin{\left (slope \right )} \cos{\left (skew \right )} - \sin{\left (skew \right )} \cos{\left (lat \right )}\right) \sin{\left (Dec \right )} + \sin{\left (HA \right )} \cos{\left (Dec \right )} \cos{\left (skew \right )} \cos{\left (slope \right )}$$

In [521]:

    

#might not be obvious, but the skew terms now become fractions or go to 1.
simplified = nsimplify(new_eq.subs({skew:values['skew']}),tolerance=0.000000001)
simplified


    

    
Out[521]:





$$\sin{\left (EW \right )} = \left(- \sin{\left (lat \right )} \sin{\left (slope \right )} - \frac{19325}{813483139} \cos{\left (lat \right )}\right) \sin{\left (Dec \right )} + \left(\frac{19325}{813483139} \sin{\left (lat \right )} - \sin{\left (slope \right )} \cos{\left (lat \right )}\right) \cos{\left (Dec \right )} \cos{\left (HA \right )} + \sin{\left (HA \right )} \cos{\left (Dec \right )} \cos{\left (slope \right )}$$

In [522]:

    

simple_ew = float(solve(simplified.subs(values),ew)[0])
print "Offset between EB (simple) and JER of:",np.degrees(simple_ew - np.arcsin(calc_ew_jer))*3600.0,"arcseconds"


    

    




Offset between EB (simple) and JER of: 3.89299893011e-10 arcseconds

In [527]:

    

new_eq.subs(consts)


    

    
Out[527]:





$$\sin{\left (EW \right )} = 0.00197739746120948 \sin{\left (Dec \right )} + 0.999994050318189 \sin{\left (HA \right )} \cos{\left (Dec \right )} - 0.00282652215698604 \cos{\left (Dec \right )} \cos{\left (HA \right )}$$

In [528]:

    

ns_eq = sym.Eq(nsew_vec[0],(telescope_to_hadec * hadec_vec)[0])
ns_eq = ns_eq.subs(consts)
ns_eq


    

    
Out[528]:





$$\cos{\left (EW \right )} \cos{\left (NS \right )} = - 0.578860911092984 \sin{\left (Dec \right )} + 0.00344945849212142 \sin{\left (HA \right )} \cos{\left (Dec \right )} + 0.815419000787148 \cos{\left (Dec \right )} \cos{\left (HA \right )}$$

In [535]:

    

inv_dec = sym.Eq(hadec_vec[2],(telescope_to_hadec.T * nsew_vec)[2])
inv_dec.subs(consts)


    

    
Out[535]:





$$\sin{\left (Dec \right )} = 0.00197739746120948 \sin{\left (EW \right )} + 0.815423899274408 \sin{\left (NS \right )} \cos{\left (EW \right )} - 0.578860911092984 \cos{\left (EW \right )} \cos{\left (NS \right )}$$

In [536]:

    

inv_ha = sym.Eq(hadec_vec[1],(telescope_to_hadec.T * nsew_vec)[1])
inv_ha.subs(consts)


    

    
Out[536]:





$$\sin{\left (HA \right )} \cos{\left (Dec \right )} = 0.999994050318189 \sin{\left (EW \right )} + 2.37557290389662 \cdot 10^{-5} \sin{\left (NS \right )} \cos{\left (EW \right )} + 0.00344945849212142 \cos{\left (EW \right )} \cos{\left (NS \right )}$$

In [545]:

    

ew_from_azel = sym.Eq(nsew_vec[2],(telescope_to_azel * azel_vec)[2])
ew_from_azel.subs(consts)


    

    
Out[545]:





$$\sin{\left (EW \right )} = - 0.999994050318189 \sin{\left (Az \right )} \cos{\left (El \right )} - 0.00344945849114808 \sin{\left (El \right )} - 2.37558703721329 \cdot 10^{-5} \cos{\left (Az \right )} \cos{\left (El \right )}$$

In [544]:

    

ns_from_azel = sym.Eq(nsew_vec[1],(telescope_to_azel * azel_vec)[1])
ns_from_azel.subs(consts)


    

    
Out[544]:





$$\sin{\left (NS \right )} \cos{\left (EW \right )} = - 2.37557290389662 \cdot 10^{-5} \sin{\left (Az \right )} \cos{\left (El \right )} - 8.19448887928893 \cdot 10^{-8} \sin{\left (El \right )} + 0.999999999717829 \cos{\left (Az \right )} \cos{\left (El \right )}$$

In [548]:

    

el_from_nsew = sym.Eq(azel_vec[2],(telescope_to_azel.T * nsew_vec)[2])
el_from_nsew.subs(consts)


    

    
Out[548]:





$$\sin{\left (El \right )} = - 0.00344945849114808 \sin{\left (EW \right )} - 8.19448887928893 \cdot 10^{-8} \sin{\left (NS \right )} \cos{\left (EW \right )} + 0.999994050600358 \cos{\left (EW \right )} \cos{\left (NS \right )}$$

In [549]:

    

az_from_nsew = sym.Eq(azel_vec[1],(telescope_to_azel.T * nsew_vec)[1])
az_from_nsew.subs(consts)


    

    
Out[549]:





$$\sin{\left (Az \right )} \cos{\left (El \right )} = - 0.999994050318189 \sin{\left (EW \right )} - 2.37557290389662 \cdot 10^{-5} \sin{\left (NS \right )} \cos{\left (EW \right )} - 0.00344945849212142 \cos{\left (EW \right )} \cos{\left (NS \right )}$$

In [574]:

    

#func_ = (asin((telescope_to_azel.T * nsew_vec)[2])).subs(consts)
#z = sym.lambdify((ns,ew),func_,'numpy',dummify=False)
#z(0.1,0.2)


    

In [573]:

    




    

In [ ]: