On 23 November 2012 09:57, Christoph Deil [email protected] wrote:

I would like to add tests for conversion to the HorizonSystem.

@dsberry https://github.com/dsberry Are there any existing reference
results?

If not we could simply use a few observatory locations and a grid of dates?

I'm not aware of any reference results myself (that's not to say they do
not exist though). Probably best to do as you suggest - use the various
packages to convert the table of standard FK5 J2000 positions to (az,el)
values assuming some arbitrary geographical position and time, and compare
the results. Maybe do this for several different geographical positions
and times.

David

Since this is not implemented in astropy.coordinates yet, I would say this is lower priority, but it would indeed be good.

I need this test anyways for another non-public code, so I want to do it now. It won't hurt to have reference results available when someone implements the transformation to HorizontalCoordinates in the future.

I'm thinking of taking the existing input.txt as (lon, lat) positions on Earth (maybe a subset) and for each (lon, lat) checking a small grid of altitudes and observation times, e.g.:

altitudes = [1e3, 6e3, 10e3, 20e3] # unit: km
times = [1900, 1950, 1990, 2000, 2010, 2020]

@dsberry Is there anything I have to be aware of, e.g. that certain coordinates are only defined for certain time ranges, or that certain Earth ellipsoid models are not defined for very low or high altitudes or ... or is this straightforward?

On 23 November 2012 11:31, Christoph Deil [email protected] wrote:

I need this test anyways for another non-public code, so I want to do it
now. It won't hurt to have reference results available when someone
implements the transformation to HorizontalCoordinates in the future.

I'm thinking of taking the existing input.txt as (lon, lat) positions on
Earth (maybe a subset) and for each (lon, lat) checking a small grid of
altitudes and observation times, e.g.:

altitudes = [1e3, 6e3, 10e3, 20e3] # unit: km
times = [1900, 1950, 1990, 2000, 2010, 2020]

@dsberry https://github.com/dsberry Is there anything I have to be
aware of, e.g. that certain coordinates are only defined for certain time
ranges, or that certain Earth ellipsoid models are not defined for very low
or high altitudes or ... or is this straightforward?

We've only used horizon coords for ground-based work, so I'm not sure about
the altitude limits. The notion of an horizon becomes a bit abstract when
you are 20,000 km high!

Make sure you specify the time correctly - for instance pyast expects the
epoch to be given in the TDB timescale, although using the TT timescale
(=TAI + 32.184) instead should introduce errors of no more than 0.03
arc-sec. And pyast also needs the dut1 value at the epoch to be supplied.
This is the UT1-UTC correction, in seconds. See
http://www.starlink.rl.ac.uk/docs/sun211.htx/node486.html.

For each earth position, how many sky positions are you going to test?

Information about AST's (az,el) coordinate system is available at
http://www.starlink.rl.ac.uk/docs/sun211.htx/node612.html#System

(search for "AZEL").

David

@cdeil - FYI, you now have push access to this repository, since you've been doing so much work on it!

In the attached generate_observers.py I take 100 random (lon, lat) and for each one these altitude and time values:

# Sample of altitudes (distances from the Earth's center) in km
altitude = np.array([6e3, 10e3], dtype='int')

# Sample of UTC times
time = ['1900-01-01', '1950-01-01', '1990-01-01', '2000-01-01', '2010-01-01', '2020-01-01']

When running the script I get this warning:

WARNING: iauUtctai: dubious year for UTC (before 1960.0 or 5 years beyond last known leap second) [astropy.time.core]

So there might be problems for pre-1960 or post-2012 times.
Should we simply not test those?

On 23 November 2012 13:00, Christoph Deil [email protected] wrote:

In the attached generate_observers.py I take 100 random (lon, lat) and
for each one these altitude and time values:

Sample of altitudes (distances from the Earth's center) in km

altitude = np.array([6e3, 10e3], dtype='int')

Sample of UTC times

time = ['1900-01-01', '1950-01-01', '1990-01-01', '2000-01-01', '2010-01-01', '2020-01-01']

When running the script I get this warning:

WARNING: iauUtctai: dubious year for UTC (before 1960.0 or 5 years beyond last known leap second) [astropy.time.core]

So there might be problems for pre-1960 or post-2012 times.
Should we simply not test those?

I think so. The current formulation in AST is based on UTC and leap
seconds, so it's problematic to use times before UTC began (in 1960) or too
far into the future (when other leap seconds may have been introduced).

David

I've added a few horizontal coordinates computed with pyephem to get started.
I'll look tomorrow how to compute them with AST and TPM.

On 23 November 2012 19:46, Christoph Deil [email protected] wrote:

I've added a few horizontal coordinates computed with pyephem to get
started.
I'll look tomorrow how to compute them with AST and TPM

For AST, it's just a case of setting the ObsLon, ObsLat, ObsAlt and DUT1
attributes in the SkyFrame (see
http://starlink.jach.hawaii.edu/docs/sun211.htx/node651.html), and then
setting the SkyFrame's System attribute to "AZEL" . The rest is just the
same as for the other systems.

David

@dsberry Is there some python package that let's me compute DUT1 for a given UTC time?

Note to self: Before doing this full test, first try to get approximately equal results from pyast, pytpm and pyephem for this one case:
http://phn.github.com/pytpm/conversions.html?highlight=azimuth#fk5-equinox-and-epoch-j2000-0-to-topocentric-observed

@cdeil - about DUT1, if astropy.time (http://docs.astropy.org/en/latest/time/index.html) doesn't do it, should it be implemented?

On 26 November 2012 13:13, Christoph Deil [email protected] wrote:

@dsberry https://github.com/dsberry Is there some python package that
let's me compute DUT1 for a given UTC time?

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/10#issuecomment-10714782.

DUT1 is determined experimentally. Historical values are available at

http://data.iers.org/products/214/14443/orig/eopc04_08_IAU2000.62-now

AST does not incorporate this table since it gets extended (in principle)
every day. Instead, AST expects the DUT1 value to be supplied with the
data, in the same way that epoch, equinox, etc, is supplied with the data.

David

If I understand correctly, DUT1 = UT1 - UTC is not implemented in astropy.time at the moment.
If it's needed to compute Horizontal coordinates, we will want it sooner or later in astropy, right?
@taldcroft and others: should I open an issue?

I see this: astropy/astropy#352

@cdeil - yes, you should go ahead and open an issue. IIRC this is nominally assigned to @iguananaut since he had someone working on an implementation that would handle getting updated tables automatically. I'm not sure about the status of that work.

I'm not so sure either--it was an intern who I think may be slowly, painfully flaking out. I should have assigned an easier task but it's hard for me to judge what "easy" is--I thought it would be. They did come up with some working regular expressions that I might be able to salvage, but otherwise if you want to move ahead with that we can do it.

I believe astropy/astropy#351 is the issue for this (I re-assigned from astropy/astropy#502, as I believe this is a duplicate).

I rebased this PR for sky to horizon coordinate checks and plan to work on this as soon as sky to horizon coordinate transformations become available in Astropy.

@eteq is there a branch or issue I can watch for sky to horizon coordinate transformations in Astropy?

@cdeil not yet because it requires the ERFA vectorizing PR to be merged. But I will be sure to @ mention you as soon as there is something to test.

@eteq I had a quick look at the ERFA vectorizing PR ... it's huge and almost impossible to review on Github. Thanks for all that work!

Yeah, it's a bit annoying, but we had to check some auto-generated code because it wasn't clear how to integrate it into the setup. The hope is to get that in for 1.1 or something so that future reviews were be easier. I can't take all the credit, though - @jwoillez did a lot of the nitty-gritty templating and numpy stuff!

@cdeil - it would be great if you could manage to get pyast and pyephem (or any other two codes) to at least roughly agree before we do it in astropy. That is, we want to make sure we have independent reference values that are validated by at least two packages. So if possible, we should try and get this PR to work and merge it before @eteq finishes implementing it in astropy.

@cdeil - one thing that will be in astropy (possibly in v1.0, possibly later, depending on time constraints) will be at least an option for atmospheric refraction corrections. So that requires a few more observer parameters:

1. pressure
2. temperature
3. humidity
4. wavelength of observation

Of course if pyast and pyephem don't support these it's not that critical to get them in here for comparison, but it might be useful to at least have an option for accuracy-checking purposes.

pyast corrects for diurnal aberattion but not for atmospheric refraction or polar motion.

@eteq In this PR I only want to do obs.pressure = 0 ... as a reminder I've opened a new issue for the extensions you suggest: #40

@cdeil - it looks like the values from the two tools here are still very different? Do you have any ideas what might be going on?

No ... I didn't have time to look at this and figure it out yet.
I want to look at composite models now and review that, but I might have time to come back to this in a few hours.

@astrofrog - Working on this now

Astropy now gives roughly consistent results with PyEphem.

I did not update the pyast tests, because I'd first like to re-structure coordinates-benchmarks into a Python package, so that we can have a utils.py and import shared functions / constants from there.

@astrofrog - OK to merge this now and continue in new PRs? Do you want to review this?

@cdeil - thanks! Would it be possible to commit the changes to the summary files before we merge this PR?

@astrofrog - The summary and ./run_benchmark.py --tasks horizontal isn't implemented yet.
Fo now you have to do this:

$ cd tools/astropy
$ python convert_horizontal.py
$ cd tools/pyephem
$ python convert_horizontal.py

I will continue working on coordinated-benchmarks this week, but I'd like to do it in separate PRs.

Ok, then feel free to merge this for now - thanks!