Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add hor2eq function #31

Merged
merged 8 commits into from Jul 6, 2017
Merged

Add hor2eq function #31

merged 8 commits into from Jul 6, 2017

Conversation

TestSubjector
Copy link
Contributor

@TestSubjector TestSubjector commented Jul 1, 2017

The function was a bit tricky to do (the Real type keyword arguments, if placed before the Bool type keyword arguments lead to type-instability).

Additionally, this is an important function because of the many functions it is dependent upon (good chance of this function's tests breaking if some other function is modified in the future).

* TODO.md: Remove hor2eq from list.
* docs/src/ref.md: Add hor2eq entry to the manual.
* src/utils.jl: Add hor2eq entry to "utils".
* src/hor2eq.jl: Contains code of hor2eq function.
* test/util-tests.jl: Include tests for hor2eq.
@codecov-io
Copy link

codecov-io commented Jul 1, 2017

Codecov Report

Merging #31 into master will increase coverage by 0.01%.
The diff coverage is 100%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master      #31      +/-   ##
==========================================
+ Coverage   99.53%   99.55%   +0.01%     
==========================================
  Files          71       73       +2     
  Lines        1292     1344      +52     
==========================================
+ Hits         1286     1338      +52     
  Misses          6        6
Impacted Files Coverage Δ
src/hor2eq.jl 100% <100%> (ø)
src/planet_coords.jl 100% <0%> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 41086a5...6f66a9b. Read the comment docs.


Code of this function is based on IDL Astronomy User's Library.
"""
hor2eq(alt::Real, az::Real, jd::Real; ws::Bool=false, B1950::Bool=false,
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In _hor2eq, put all floating point arguments one after the other, so that you can promote all of them together:

_hor2eq(promote(float(alt), float(az), float(jd), float(lat), float(lon), float(altitude),
                float(temperature), float(pressure))..., ws,
        B1950, precession, nutate, aberration, refract,
        obsname)

* `alt`: altitude of horizon coords, in degrees
* `az`: azimuth angle measured East from North (unless ws is `true`), in degrees
* `jd`: julian date
* `ws` (optional boolean keyword): set this to `true` to get the azimuth measured
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Document the default value

src/hor2eq.jl Outdated
if ws
az -= 180
end
dra1, ddec1, eps, d_psi, _ = co_nutate(jd, 45, 45)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

_ is useful to discard leading or intermediate elements of an iterable (_, b, _, d = 1, 2, 3, 4), but for trailing elements there is no need to use _,

_, _, eps, d_psi = co_nutate(jd, 45, 45)

should be sufficient ;-)

src/hor2eq.jl Outdated
ha, dec = altaz2hadec(alt_b, az, lat)
ra = mod(last - ha, 360)
ha /= 15
dra1, ddec1, eps, _, _ = co_nutate(jd, ra, dec)
Copy link
Member

@giordano giordano Jul 3, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Like above,

dra1, ddec1, eps = co_nutate(jd, ra, dec)

In addition, can this line be moved inside the if nutate?

src/hor2eq.jl Outdated
return ra, dec, ha
end
"""

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Something missing here? Please, leave a blank line before the docstring

ra_o, dec_o, ha_o = hor2eq(43.6879, 56.684, AstroLib.J2000, B1950=true)
@test ra_o ≈ 259.20005705918317
@test dec_o ≈ 49.674706171288655
@test ha_o ≈ 19.40545272279752
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

GDL> hor2eq, 43.6879d0, 56.684d0, 2451545d0, ra_o, dec_o, ha_o, /B1950
GDL> print, ra_o, dec_o, ha_o
       259.20239       49.674907       291.08179

Do you have a clue about why ra_o and dec_o differ already on the 4th-5th most significant digit?

ha_o is completely different: you converted to hours, but docstring says it's in degrees, but when I convert back to degrees the value of this quantity exactly matches the value in IDL AstroLib, so there is no problem with it, only the RA and Dec coordinates.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the difference stems from the use of the updated values in true_obliquity and mean_obliquity; in co_aberration and co_nutate respectively ( #24 )

src/hor2eq.jl Outdated
last = 15 * ct2lst(lon, jd) + d_psi * cos(eps) / 3600
ha, dec = altaz2hadec(alt_b, az, lat)
ra = mod(last - ha, 360)
ha /= 15
Copy link
Member

@giordano giordano Jul 3, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

According to docstring, ha is returned in degrees, not hours, this instruction should be removed (see comment about the result of the test below)

src/hor2eq.jl Outdated
* `lon` (optional keyword): AST longitude of location, in degrees. You can specify west
longitude with a negative sign. Dafault is NaN.
* `altitude` (optional keyword): the altitude of the observing location, in meters.
It it zero by default
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

"It is"

src/hor2eq.jl Outdated
Default is NaN
* `obsname` (optional keyword): set this to a valid observatory name to
be used by the observatory type in [types](@ref), which will return
the latitude and longitude to be used by this program. Default is NaN.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should be defined better. The default is the empty string "", not NaN. In addition, it should be explained that when it's an empty string, lat and lon defaults to the coordinates of Pine Bluff Observatory.

src/hor2eq.jl Outdated
ra = mod(last - ha, 360)
ha /= 15
dra1, ddec1, eps, _, _ = co_nutate(jd, ra, dec)
dra2, ddec2 = co_aberration(jd, ra, dec, eps)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can this be moved inside the if aberration?

src/hor2eq.jl Outdated
dra2, ddec2 = co_aberration(jd, ra, dec, eps)

if nutate
ra -= dra1
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Shouldn't dra1 and ddec1 in following line be divided by 3600?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

co_aberration.pro returned the ra and dec in arcseconds (that's why the division by 3600 to make it into degrees), however co_aberration.jl returns the ra and dec in degrees itself.

src/hor2eq.jl Outdated
altitude = observatories[obsname].altitude
end

if refract
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this if can be simplified to

    if refract
        alt = co_refract(alt, altitude, pressure, temperature)
    end

and use alt instead of alt_b below

src/hor2eq.jl Outdated
lon = T(-89.865)
end
else
lat = observatories[obsname].latitude
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Like above (good!), use

lat = T(...)
lon = T(...)
altitude = T(...)
...

because in observatories those quantities are stored as Float64.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks. I had done that till f818e8e actually, changed that in attempt for type stability and didn't revert it back later

@giordano
Copy link
Member

giordano commented Jul 4, 2017

As you noted, type-stability of this function is a bit tricky because of the keyword arguments. Keyword arguments in Julia are handy, but also often a pain when you seek performance, because they aren't part of dispatch.

The inner function _hor2eq instead is completely inferrable, for example try

@code_warntype AstroLib._hor2eq(43.6879, 56.684, Float64(AstroLib.J2000), NaN, NaN, 0.0, NaN, NaN, false, false, true, true, true, true, "")

@code_warntype AstroLib._hor2eq(big(43.6879), big(56.684), BigFloat(AstroLib.J2000), big(NaN), big(NaN), big(0.0), big(NaN), big(NaN), false, false, true, true, true, true, "")

Since at least the inner function looks stable (which is the most important thing), I'd probably merge the PR after the latest requested fixes, but please add at the end of the file drop a TODO note saying that hor2eq isn't stable because of the keyword arguments.

src/hor2eq.jl Outdated

### Output ###

* `ra`: right ascension of object, in degrees (FK5)
* `dec`: declination of the object, in degrees (FK5)
* `ha`: hour angle, in degrees
* `ha`: hour angle, in hours
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why this change? I'd have kept the unit in degrees, for consistency with IDL AstroLib

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I assumed that you wanted to change the docstring to hours, my bad. Updated.

src/hor2eq.jl Outdated

julia> adstring(ra_o, dec_o)
" 00 13 17.3 +15 11 26"
" 00 13 17.3 +15 11 26
Copy link
Member

@giordano giordano Jul 4, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed the closing quote by error?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes

* Clean code
* Add TODO to make hor2eq type-stable in the future
* First attempt at using @inferred
@TestSubjector
Copy link
Contributor Author

Made a small attempt at using @inferred in tests, as it gives the side knowledge that _hor2eq has a inferable return type.

@giordano giordano merged commit 771b977 into JuliaAstro:master Jul 6, 2017
@TestSubjector TestSubjector deleted the hor2eq branch July 6, 2017 16:09
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

3 participants