Repository navigation
Conversation
See the comments added to the code for detail about why this approach improves things. As a side benefit, the special case handling for the poles can be removed. This will change the results at the poles, in that orientation will be determined by the given RA. The golden RA/DEC in the unit tests changed. These are just a paste of what is coming out, rather than independently verified results.
|
@mauch FYI. |
ludwigschwardt
left a comment
ludwigschwardt
left a comment
There was a problem hiding this comment.
Great stuff!
It's interesting that you only need to fix the m direction / v vector via an offset target. The distortion you mention therefore probably keeps l and m pretty much orthogonal locally.
DiFX also introduces an l offset target to fix the u vector independently of v (instead of insisting that they are orthogonal by using a cross product). I wonder how big that effect is...
| # conversion doesn't simply rotate the celestial sphere, but also | ||
| # distorts it, and so that introduced errors. | ||
| # | ||
| # To avoid issues close to the poles, we always offset towards the |
| ra, dec = self.radec(None, antenna) | ||
| else: | ||
| ra, dec = self.radec(timestamp, antenna) | ||
| offset_sign = -1 if dec > 0 else -1 |
There was a problem hiding this comment.
You got -1 twice...
There was a problem hiding this comment.
Oops, thanks! Will fix.
| # To avoid issues close to the poles, we always offset towards the | ||
| # equator. We also can't offset by too little, as ephem uses only | ||
| # single precision and this method suffers from loss of precision. | ||
| # 0.03 was found by experimentation (albeit on a single data set) to |
There was a problem hiding this comment.
0.03 radians, I assume. That's about 2 degrees... Interesting that one can go that big. I guess it is a mild distortion to begin with.
There was a problem hiding this comment.
I suspect we'd get slightly better results if we could go smaller, but then the floating point errors climb rapidly.
| # u axis is orthogonal to z and w, and row_stack makes it 2-D array of column vectors | ||
| u = np.row_stack(np.cross(z, w, axis=0)) * offset_sign | ||
| u_norm = np.sqrt(np.sum(u ** 2, axis=0)) | ||
| # If the target is a celestial pole (so that w equals z or -z), u and v become degenerate |
That's a good point. A related issue is scaling - we just normalise the vector we get back, but a unit vector in geocentric coordinates won't be a unit long in topocentric coordinates. But I think that locally the effects scale with β² rather than β (v/c ~= 1e-4) so are too small to worry about. |
It doesn't affect the sign of the result, just leads to stability problems within a few degrees of the south pole.
The previous bug-fix on whether to displace north or south caused the results to change by enough to throw off the unit tests (a few ppm).
See the comments added to the code for detail about why this approach
improves things. As a side benefit, the special case handling for the
poles can be removed. This will change the results at the poles, in that
orientation will be determined by the given RA.
The golden RA/DEC in the unit tests changed. These are just a paste of
what is coming out, rather than independently verified results.