AzEl issue

[twillis449]

As Oleg points out, Wim Brouw is infallible but that doesn't prevent the rest of us from going astray. I've concluded that the AzEl and AzElRaDec nodes in MeqTrees, while correct as they stand, do not provide the right AzEl for the working astronomer. Wim's 'default' AzEl ( MDirection::AZEL) gives the AzEl of a radio source as seen from the centre of the Earth using the geocentric latitude. This is useless for practical astronomy as the AzEl as seen at the site of an observatory is given by MDirection::AZELGEO the AzEl of a source as seen from a position on the actual oblate Earth using the geodetic latitude and not the geocentric latitude. I suggest creating two new MeqTrees nodes AzElGeo and AzElGeoRaDec. All that has to be done is replace instances of MDirection::AZEL by MDirection::AZELGEO in the *.cc files and references to topocentric by geodetic in the corresponding header files'

Since no-one has ever noticed this error, I suspect that conversions between ra, dec and az, el haven't been done too often !!

[twillis449]

PS I'd do this myself, but I remember that there's some obscure 'special' incantations necessary when adding nodes. An alternative would be to just update the existing nodes with the changes mentioned above. I don't see any practical use for nodes giving geocentric Az and El.

[twillis449]

BTW on the (really old) astron wiiki there used to be a page that contained small test scripts for most MeqTrees nodes. Is this page still hidden in some alternative universe where it can be accessed in some way?

[o-smirnov]

Ahhh, well spotted Tony. The sacred tools must be wielded well or they can give incorrect results. I had completely overlooked this -- in fact the intention of MeqAzEl was to give the geodetic az/el. So I wonder if we shouldn't just change MeqAzEl itself (or perhaps provide an option in the init state that selects geodetic versus geocentric)?

Do you have a feeling for how large the difference can be? I'm now wondering if this also affects MeqParAngle -- I've had some puzzling results with the polarization of 3C286 around transit which I could never explain...

[sjperkins]

Supporting argument. @bmerry changed AZEL to AZELGEO when comparing meqtrees and astropy parallactic angles. To get them to agree more closely I think.

[twillis449]

The differences can be up to 12 arcmin or so, and depends on latitude and peaks at about +/-49 deg (guess what the latitude of the DRAO telescope is?) which is why I started to notice the problem - we were comparing CASA Measures output vs astropy vs JPL HORIZON vs PyEphem for the Sun. And yes, it affects parallactic angle. There the problem is that MeqParAngle calls the CASA ParAngleMachine (@tammojan) which has MDirection::AZEL hard wired in the code, so that's an issue for CASA to look at. My personal suggestion would be for options to be provided to select between AZELGEO and AZEL, but if thats hard to do, then change AZEL to AZELGEO as that's where Az and El are most likely to be used.

@o-smirnov If you were observing 3C286 with the WSRT, which is at 55 deg latitude or thereabouts, yes, you should be seeing some discrepancies due to ParAngleMachine using a geocentric definition for AzEl.

[twillis449]

Here's a plot of the difference in geocentric and geodetic latitude (0 at the poles and equator - peak at somewhere close to 49 deg)

[twillis449]

@o-smirnov Of course the parallactic angle difference at transit is zero, but increases as hour angle increases. (Further investigation shows that the above statement is not necessarily true - see the output from my next message).

Here's the output of a run I did for 3C286 as seen from DRAO over a 24 hr period

astropy 3C286 J2000 position RA 13h31m8.287896s Dec = 30d30m32.958s starting 3C286 test at 13jul2017 0 hrs columns are: utc, ha( deg), elev (deg), par_angle_difference (arcmin) 2017-07-13 00:00:00, -31.5470796713, 59.6606265431, 20.2431429995 2017-07-13 01:00:00, -16.5060246335, 67.3621324507, 19.0001225853 2017-07-13 02:00:00, -1.46497090709, 71.073633595, 2.41623046051 2017-07-13 03:00:00, 13.5760827059, 68.4995286627, -17.3287871344 2017-07-13 04:00:00, 28.6171374012, 61.3368250813, -20.5619110418 2017-07-13 05:00:00, 43.6581943053, 52.2495930473, -18.1543667764 2017-07-13 06:00:00, 58.6992543989, 42.5358774264, -15.4965697269 2017-07-13 07:00:00, 73.7403184335, 32.7795959076, -13.3678270527 2017-07-13 08:00:00, 88.7813868989, 23.327245502, -11.6695155851 2017-07-13 09:00:00, 103.822459978, 14.4619174011, -10.1932699575 2017-07-13 10:00:00, 118.863537536, 6.47054068704, -8.73234462784 2017-07-13 11:00:00, 133.904619139, -0.33001787016, -7.09492501601 2017-07-13 12:00:00, 148.945704067, -5.5954470066, -5.12973164661 2017-07-13 13:00:00, 163.986791376, -8.98948433497, -2.78573514753 2017-07-13 14:00:00, 179.027879959, -10.2522664803, -0.172629040492 2017-07-13 15:00:00, -165.931031376, -9.27671593682, 2.45882773172 2017-07-13 16:00:00, -150.889943834, -6.14608521304, 4.84748700835 2017-07-13 17:00:00, -135.848858533, -1.10445891552, 6.86173459936 2017-07-13 18:00:00, -120.80777644, 5.51737813945, 8.5344124133 2017-07-13 19:00:00, -105.766698312, 13.3733913114, 10.0078198155 2017-07-13 20:00:00, -90.7256246201, 22.1423955269, 11.470680677 2017-07-13 21:00:00, -75.6845555388, 31.5344473046, 13.1283194415 2017-07-13 22:00:00, -60.6434909251, 41.268560906, 15.1909246937 2017-07-13 23:00:00, -45.6024303324, 51.0127661141, 17.7890534714

[twillis449]

When I make the same run at the VLA site (latitude ~ 34 deg) on the same day I get (just for hour angles in the range -90 deg to 90 deg) I get:

astropy 3C286 J2000 position RA 13h31m8.287896s Dec = 30d30m32.958s starting 3C286 test at 13jul2017 0 hrs utc,ha(deg),el(deg),par_angle_diff(arcmin) 2017-07-13 00:00:00, -19.5346278106, 73.1084458029, 36.5865480434 2017-07-13 01:00:00, -4.49357856546, 84.729070724, 87.6388644618 2017-07-13 02:00:00, 10.5474702201, 80.3667872774, -60.6740682777 2017-07-13 03:00:00, 25.5885200826, 68.1122773707, -28.6831032444 2017-07-13 04:00:00, 40.6295724746, 55.6692030008, -18.880641269 2017-07-13 05:00:00, 55.6706286806, 43.3563539577, -14.4032344934 2017-07-13 06:00:00, 70.711689729, 31.3103842237, -11.9233838637 2017-07-13 07:00:00, 85.7527563061, 19.6728358025, -10.3722509609 2017-07-13 20:00:00, -78.7131496006, 25.0583140884, 11.0199052754 2017-07-13 21:00:00, -63.6720856493, 36.9066455609, 12.9253148947 2017-07-13 22:00:00, -48.6310269459, 49.0936754264, 16.1386811822 2017-07-13 23:00:00, -33.5899726392, 61.4866830311, 22.3950745577 2017-07-14 00:00:00, -18.5489215685, 73.917785376, 38.3031954131