Molonglo coordinate transforms

Useful coordinate transforms for the molonglo radio telescope


In [ ]:
import numpy as np
import ephem as e
from scipy.optimize import minimize
import matplotlib.pyplot as plt
np.set_printoptions(precision=5,suppress =True)

Below we define the rotation and reflection matrices


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

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 np.array(m[d])

Define a position vectors


In [ ]:
def pos_vector(a,b):
    return np.array([[np.cos(b)*np.cos(a)],
                     [np.cos(b)*np.sin(a)],
                     [np.sin(b)]])

def pos_from_vector(vec):
    a,b,c = vec
    a_ = np.arctan2(b,a)
    c_ = np.arcsin(c)   
    return a_,c_

Generic transform


In [ ]:
def transform(a,b,R,inverse=True):
    P = pos_vector(a,b)
    if inverse:
        R = R.T
    V = np.dot(R,P).ravel()
    a,b = pos_from_vector(V)
    a = 0 if np.isnan(a) else a
    b = 0 if np.isnan(a) else b
    return a,b

Reference conversion formula from Duncan's old TCC


In [ ]:
def hadec_to_nsew(ha,dec):
    ew = np.arcsin((0.9999940546 * np.cos(dec) * np.sin(ha))
                    - (0.0029798011806 * np.cos(dec) * np.cos(ha))
                    + (0.002015514993 * np.sin(dec)))
    ns = np.arcsin(((-0.0000237558704 * np.cos(dec) * np.sin(ha))
                     + (0.578881847 * np.cos(dec) * np.cos(ha))
                     + (0.8154114339 * np.sin(dec)))
                     / np.cos(ew))
    return ns,ew

New conversion formula using rotation matrices

What do we think we should have:

\begin{equation} \begin{bmatrix} \cos(\rm EW)\cos(\rm NS) \\ \cos(\rm EW)\sin(\rm NS) \\ \sin(\rm EW) \end{bmatrix} = \mathbf{R} \begin{bmatrix} \cos(\delta)\cos(\rm HA) \\ \cos(\delta)\sin(\rm HA) \\ \sin(\delta) \end{bmatrix} \end{equation}

Where $\mathbf{R}$ is a composite rotation matrix.

We need a rotations in axis of array plus orthogonal rotation w.r.t. to array centre. Note that the NS convention is flipped so HA and NS go clockwise and anti-clockwise respectively when viewed from the north pole in both coordinate systems.

\begin{equation} \mathbf{R}_x = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos(\theta) & -\sin(\theta) \\ 0 & \sin(\theta) & \cos(\theta) \end{bmatrix} \end{equation}\begin{equation} \mathbf{R}_y = \begin{bmatrix} \cos(\phi) & 0 & \sin(\phi) \\ 0 & 1 & 0 \\ -\sin(\phi) & 0 & \cos(\phi) \end{bmatrix} \end{equation}\begin{equation} \mathbf{R}_z = \begin{bmatrix} \cos(\eta) & -\sin(\eta) & 0\\ \sin(\eta) & \cos(\eta) & 0\\ 0 & 0 & 1 \end{bmatrix} \end{equation}\begin{equation} \mathbf{R} = \mathbf{R}_x \mathbf{R}_y \mathbf{R}_z \end{equation}

Here I think $\theta$ is a $3\pi/2$ rotation to put the telescope pole (west) at the telescope zenith and $\phi$ is also $\pi/2$ to rotate the telescope meridian (which is lengthwise on the array, what we traditionally think of as the meridian is actually the equator of the telescope) into the position of $Az=0$.

However rotation of NS and HA are opposite, so a reflection is needed. For example reflection around a plane in along which the $z$ axis lies:

\begin{equation} \mathbf{\bar{R}}_z = \begin{bmatrix} 1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & -1 \end{bmatrix} \end{equation}

Conversion to azimuth and elevations should therefore require $\theta=-\pi/2$ and $\phi=\pi/2$ with a reflection about $x$.

Taking into account the EW skew and slope of the telescope:

\begin{equation} \begin{bmatrix} \cos(\rm EW)\cos(\rm NS) \\ \cos(\rm EW)\sin(\rm NS) \\ \sin(\rm EW) \end{bmatrix} = \begin{bmatrix} \cos(\alpha) & -\sin(\alpha) & 0\\ \sin(\alpha) & \cos(\alpha) & 0\\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \cos(\beta) & 0 & \sin(\beta) \\ 0 & 1 & 0 \\ -\sin(\beta) & 0 & \cos(\beta) \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & -1 & 0 \end{bmatrix} \begin{bmatrix} 0 & 0 & -1 \\ 0 & 1 & 0 \\ 1 & 0 & 0 \end{bmatrix} \begin{bmatrix} -1 & 0 & 0\\ 0 & 1 & 0\\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \cos(\delta)\cos(\rm HA) \\ \cos(\delta)\sin(\rm HA) \\ \sin(\delta) \end{bmatrix} \end{equation}

So the correction matrix to take telescope coordinates to ns,ew

\begin{bmatrix} \cos(\alpha)\sin(\beta) & -\sin(\beta) & \cos(\alpha)\sin(\beta) \\ \sin(\alpha)\cos(\beta) & \cos(\alpha) & \sin(\alpha)\sin(\beta) \\ -\sin(\beta) & 0 & \cos(\beta) \end{bmatrix}

and to Az Elv

\begin{bmatrix} \sin(\alpha) & -\cos(\alpha)\sin(\beta) & -\cos(\alpha)\cos(\beta) \\ \cos(\alpha) & -\sin(\alpha)\sin(\beta) & -\sin(\alpha)\cos(\beta) \\ -\cos(\beta) & 0 & \sin(\beta) \end{bmatrix}

In [ ]:
# There should be a slope and tilt conversion to get accurate change
#skew = 4.363323129985824e-05
#slope = 0.0034602076124567475

#skew =  0.00004
#slope = 0.00346

skew = 0.01297 # <- this is the skew I get if I optimize for the same results as duncan's system
slope= 0.00343

def telescope_to_nsew_matrix(skew,slope):
    R = rotation_matrix(skew,"z")
    R = np.dot(R,rotation_matrix(slope,"y"))
    return R

def nsew_to_azel_matrix(skew,slope):
    pre_R = telescope_to_nsew_matrix(skew,slope)
    x_rot = rotation_matrix(-np.pi/2,"x")
    y_rot = rotation_matrix(np.pi/2,"y")
    R = np.dot(x_rot,y_rot)
    R = np.dot(pre_R,R)
    R_bar = reflection_matrix("x")
    R = np.dot(R,R_bar)
    return R

def nsew_to_azel(ns, ew):    
    az,el = transform(ns,ew,nsew_to_azel_matrix(skew,slope))
    return az,el

print nsew_to_azel(0,np.pi/2) # should be -pi/2 and 0
print nsew_to_azel(-np.pi/2,0)# should be -pi and 0
print nsew_to_azel(0.0,.5)    # should be pi/2 and something near pi/2
print nsew_to_azel(-.5,.5)    # less than pi/2 and less than pi/2
print nsew_to_azel(.5,-.5)    
print nsew_to_azel(.5,.5)

The inverse of this is:


In [ ]:
def azel_to_nsew(az, el):  
    ns,ew = transform(az,el,nsew_to_azel_matrix(skew,slope).T)
    return ns,ew

Extending this to HA Dec


In [ ]:
mol_lat = -0.6043881274183919 # in radians

def azel_to_hadec_matrix(lat):
    rot_y = rotation_matrix(np.pi/2-lat,"y")
    rot_z = rotation_matrix(np.pi,"z")
    R = np.dot(rot_y,rot_z)
    return R

def azel_to_hadec(az,el,lat):
    ha,dec = transform(az,el,azel_to_hadec_matrix(lat))
    return ha,dec

def nsew_to_hadec(ns,ew,lat,skew=skew,slope=slope):
    R = np.dot(nsew_to_azel_matrix(skew,slope),azel_to_hadec_matrix(lat))
    ha,dec = transform(ns,ew,R)
    return ha,dec

ns,ew = 0.8,0.8
az,el = nsew_to_azel(ns,ew)
print "AzEl:",az,el
ha,dec = azel_to_hadec(az,el,mol_lat)
print "HADec:",ha,dec
ha,dec = nsew_to_hadec(ns,ew,mol_lat)
print "HADec2:",ha,dec

# This is Duncan's version
ns_,ew_ = hadec_to_nsew(ha,dec)
print "NSEW Duncan:",ns_,ew_
print "NS offset:",ns_-ns,"   EW offset:",ew_-ew

In [ ]:
def test(ns,ew,skew,slope):
    ha,dec = nsew_to_hadec(ns,ew,mol_lat,skew,slope)
    ns_,ew_ = hadec_to_nsew(ha,dec)
    no,eo = ns-ns_,ew-ew_
    no = 0 if np.isnan(no) else no
    eo = 0 if np.isnan(eo) else eo
    return no,eo

ns = np.linspace(-np.pi/2+0.1,np.pi/2-0.1,10)
ew = np.linspace(-np.pi/2+0.1,np.pi/2-0.1,10)

def test2(a):
    skew,slope = a
    out_ns = np.empty([10,10])
    out_ew = np.empty([10,10])
    for ii,n in enumerate(ns):
        for jj,k in enumerate(ew):
            a,b = test(n,k,skew,slope)
            out_ns[ii,jj] = a
            out_ew[ii,jj] = b
    a = abs(out_ns).sum()#abs(np.median(out_ns))
    b = abs(out_ew).sum()#abs(np.median(out_ew))
    print a,b
    print max(a,b)
    return max(a,b)   

#minimize(test2,[skew,slope])

In [ ]:
# Plotting out the conversion error as a function of HA and Dec. 
# Colour scale is log of the absolute difference between original system and new system


ns = np.linspace(-np.pi/2,np.pi/2,10)
ew = np.linspace(-np.pi/2,np.pi/2,10)
out_ns = np.empty([10,10])
out_ew = np.empty([10,10])
for ii,n in enumerate(ns):
    for jj,k in enumerate(ew):
        print jj
        a,b = test(n,k,skew,slope)
        out_ns[ii,jj] = a
        out_ew[ii,jj] = b
plt.figure()
plt.subplot(121)
plt.imshow(abs(out_ns),aspect="auto")
plt.colorbar()

plt.subplot(122)
plt.imshow(abs(out_ew),aspect="auto")
plt.colorbar()
plt.show()

In [ ]:
from mpl_toolkits.mplot3d import Axes3D
from itertools import product, combinations
from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d import proj3d
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.set_aspect("equal")

#draw sphere
u, v = np.mgrid[0:2*np.pi:20j, 0:np.pi:10j]
x=np.cos(u)*np.sin(v)
y=np.sin(u)*np.sin(v)
z=np.cos(v)
ax.plot_wireframe(x, y, z, color="r",lw=1)

R = rotation_matrix(np.pi/2,"x")
pos_v = np.array([[x],[y],[z]])
p = pos_v.T
for i in p:
    for j in i:
        j[0] = np.dot(R,j[0])

        
class Arrow3D(FancyArrowPatch):
    def __init__(self, xs, ys, zs, *args, **kwargs):
        FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
        self._verts3d = xs, ys, zs

    def draw(self, renderer):
        xs3d, ys3d, zs3d = self._verts3d
        xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
        self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
        FancyArrowPatch.draw(self, renderer)

a = Arrow3D([0,1],[0,0.1],[0,.10], mutation_scale=20, lw=1, arrowstyle="-|>", color="k")
ax.add_artist(a)        
        
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
        
x=p.T[0,0]
y=p.T[1,0]
z=p.T[2,0]
ax.plot_wireframe(x, y, z, color="b",lw=1)
plt.show()

In [ ]:


In [ ]:


In [ ]: