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 [ ]: