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 precess_cd function #16

Merged
merged 8 commits into from Jun 9, 2017

Conversation

Projects
None yet
4 participants
@TestSubjector
Contributor

TestSubjector commented Jun 3, 2017

*TODO.md: Remove 'precess_cd' from list.
*docs/src/ref.md: Add "precess_cd" entry to the manual.
*src/utils.jl: Add precess_cd entry to "utils".
*src/precess_cd.jl: Contains code of precess_cd function.
*test/util-tests.jl: Include tests for precess_cd.

TestSubjector added some commits Jun 3, 2017

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

This comment has been minimized.

codecov-io commented Jun 3, 2017

Codecov Report

Merging #16 into master will decrease coverage by 0.17%.
The diff coverage is 100%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master      #16      +/-   ##
==========================================
- Coverage   99.71%   99.54%   -0.18%     
==========================================
  Files          61       62       +1     
  Lines        1055     1095      +40     
==========================================
+ Hits         1052     1090      +38     
- Misses          3        5       +2
Impacted Files Coverage Δ
src/precess_cd.jl 100% <100%> (ø)
src/ten.jl 85.71% <0%> (-14.29%) ⬇️
src/types.jl 82.35% <0%> (-5.15%) ⬇️
src/get_date.jl 100% <0%> (ø) ⬆️
src/helio_jd.jl 100% <0%> (ø) ⬆️
src/adstring.jl 100% <0%> (ø) ⬆️
src/gal_uvw.jl 100% <0%> (ø) ⬆️
src/helio_rv.jl 100% <0%> (ø) ⬆️
src/precess.jl 100% <0%> (ø) ⬆️
src/precess_xyz.jl 100% <0%> (ø) ⬆️
... and 6 more

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 9deab8a...1c65feb. Read the comment docs.

function _precess_cd{T<:AbstractFloat}(cd::Matrix{T}, epoch1::T, epoch2::T, crval_old::Vector{T},
crval_new::Vector{T}, FK4::Bool)
crvalold = deg2rad.(crval_old)
crvalnew = deg2rad.(crval_new)

This comment has been minimized.

@giordano

giordano Jun 3, 2017

Member

Julia sind and cosd functions that allows you to avoid defining crvalold and crvalnew.

c = t*(20043.109 - st*(85.33 + st*0.217) + t*(-42.665 - st*0.217 - t*41.8))/3600
end
pole_ra = 0
pole_dec = 90

This comment has been minimized.

@giordano

giordano Jun 3, 2017

Member

pole_ra and pole_dec are type-unstable! Use

    pole_ra = zero(T)
    pole_dec = one(T) * 90

instead ;-)

pole_ra, pole_dec = precess(pole_ra, pole_dec, epoch1, epoch2, FK4=FK4)
end
sind1, sind2 = sin.([crvalold[2], crvalnew[2]])
cosd1, cosd2 = cos.([crvalold[2], crvalnew[2]])

This comment has been minimized.

@giordano

giordano Jun 3, 2017

Member

This writing is concise, but you're also allocating useless arrays. You can simply do

sind1 = sind(crval_old[2])
sind2 = ...

(remember the comment above about sind and cosd)

sinra = sin(crvalnew[1] - deg2rad.(pole_ra))
cosfi = (cos(c) - sind1*sind2)/(cosd1*cosd2)
sinfi = (abs(sin(c))* sinra)/cosd1
r = [[cosfi -sinfi]; [sinfi cosfi]]

This comment has been minimized.

@giordano

giordano Jun 3, 2017

Member

This can be simplified to

r = [cosfi -sinfi; sinfi cosfi]

which is also much faster and allocates a single array, instead of 3.

@@ -0,0 +1,78 @@
# This file is a part of AstroLib.jl. License is MIT "Expat".
function _precess_cd{T<:AbstractFloat}(cd::Matrix{T}, epoch1::T, epoch2::T, crval_old::Vector{T},

This comment has been minimized.

@giordano

giordano Jun 3, 2017

Member

For more generality, use AbstractMatrix and AbstractVector in place of Matrix and Vector ;-)

### Example ###
```julia
julia> precess_cd([[20 60]; [45 45]], 1950, 2000, [34, 58], [12, 83])

This comment has been minimized.

@giordano

giordano Jun 3, 2017

Member

[[20 60]; [45 45]] can be simplified to [20 60; 45 45]. Also in tests.

This function should not be used for values more than 2.5 centuries from the year 1900.
This function calls [precess](@ref) and [bprecess](@ref).
"""
precess_cd{R<:Real}(cd::Matrix{R}, epoch1::Real, epoch2::Real, crval_old::Vector{R},

This comment has been minimized.

@giordano

giordano Jun 3, 2017

Member

Like above, for more generality use AbstractMatrix and AbstractVector in place of Matrix and Vector. In addition, here we're not interested in annotating the R type so exploiting new syntax you can simply write

precess_cd(cd::AbstractMatrix{<:Real}, epoch1::Real, epoch2::Real, crval_old::AbstractVector{<:Real},
           crval_new::AbstractVector{<:Real}, FK4::Bool=false) =

This couldn't be done in previous Julia <= 0.5. Have a look at the third point of "New type system capabilities" at https://github.com/JuliaLang/julia/blob/release-0.6/NEWS.md#new-language-features

cosfi = (cos(c) - sind1*sind2)/(cosd1*cosd2)
sinfi = (abs(sin(c))* sinra)/cosd1
r = [[cosfi -sinfi]; [sinfi cosfi]]
cd = cd * r

This comment has been minimized.

@giordano

giordano Jun 3, 2017

Member

Pay attention here! You're mutating the input cd array. This may a good idea, but if you do so, you should rename the function to precess_cd! (and document that the input array is mutated), because the exclamation mark appended to function name is used to denote a mutating function. Furthermore, with improved loop fusion this could be better written as cd .= cd * r

This comment has been minimized.

@giordano

giordano Jun 4, 2017

Member

On a second thought, I think that's probably better not to mutate the input array. Thus, I suggest to delete the assignment cd = ... and simply return cd * r (but also look at the order of operands, I'm not sure it's correct).

This comment has been minimized.

@kbarbary

kbarbary Jun 4, 2017

Member

For the record, this assignment does not mutate the input array. It first creates a new array cd * r then assigns the value to the identifier cd. The identifier cd no longer points to the input array, which is unmodified.

But yeah, simply return cd * r would be good here.

This comment has been minimized.

@giordano

giordano Jun 4, 2017

Member

Uhm, you're right! I was doing tests and noted that the input cd matrix was changing, but probably I had replaced the line with cd .= cd * r, which does mutate.

This comment has been minimized.

@TestSubjector

TestSubjector Jun 5, 2017

Contributor

In precess_cd.pro, the order of operands is r # cd (r * cd). However, from the IDL documentation, the matrix multiplication in IDL is done by multiplying the columns of the first array by the rows of the second array. So I reversed the the operands to make it a normal matrix multiplication.

@@ -441,6 +441,22 @@ let
@test dec2 ≈ [-56.87186126487889]
end
# Test precess_cd
@testset "precess_cd" begin
@test precess_cd([[30 60]; [60 90]], 1950, 2000, [13, 8], [43, 23]) ≈

This comment has been minimized.

@giordano

giordano Jun 3, 2017

Member

Did you compare the results of these tests with IDL AstroLib?

This comment has been minimized.

@TestSubjector

TestSubjector Jun 5, 2017

Contributor

Nope, haven't been able to run IDL on my distro till now.

This comment has been minimized.

@giordano

giordano Jun 5, 2017

Member

Check out http://gnudatalanguage.sourceforge.net/ then, most AstroLib procedures work in GDL ;-)

sinra = sind(crval_new[1] - pole_ra)
cosfi = (cos(c) - sind1*sind2)/(cosd1*cosd2)
sinfi = (abs(sin(c))* sinra)/cosd1
r = [cosfi sinfi; -sinfi cosfi]

This comment has been minimized.

@TestSubjector

TestSubjector Jun 5, 2017

Contributor

I had transposed this r = [cosfi -sinfi; sinfi cosfi] because I initially wanted to keep operands sequence as r * cd, but then I changed my mind about the sequence and forgot to rectify the transposing. Corrected.

### Arguments ###
* `cd`: 2 x 2 coordinate description matrix in degrees or radians

This comment has been minimized.

@giordano

giordano Jun 5, 2017

Member

It seems that the matrix has to be in degrees. I think that the original description in inaccurate, they also say the matrix is altered to match the units of input matrix, but they don't do so.

else
pole_ra, pole_dec = precess(pole_ra, pole_dec, epoch1, epoch2, FK4=FK4)
end
sind1, sind2 = sind.([crval_old[2], crval_new[2]])

This comment has been minimized.

@giordano

giordano Jun 5, 2017

Member

Don't use the arrays here ;-) If you still want to have two calculation in one line you can do

sind1, sind2 = sind(crval_old[2]), sind(crval_new[2])
...

but I don't think you can gain much with respect to

sind1 = sind(crval_old[2])
sind2 = sind(crval_new[2])
...

That's up to you

function _precess_cd{T<:AbstractFloat}(cd::AbstractMatrix{T}, epoch1::T, epoch2::T, crval_old::AbstractVector{T},
crval_new::AbstractVector{T}, FK4::Bool)
t = deg2rad((epoch2 - epoch1)*0.001)

This comment has been minimized.

@giordano

giordano Jun 5, 2017

Member

I think that you shouldn't call deg2rad here, t isn't simply a global multiplicative factor in definition of c, which is rather a third-order polynomial in t.

if FK4
st = (epoch1 - 1900)*0.001
c = t*(20046.85 - st*(85.33 + st*0.37) + t*(-42.67 - st*0.37 - t*41.8))/3600

This comment has been minimized.

@giordano

giordano Jun 5, 2017

Member

Remove division by 3600 and wrap the whole right-hand-side with sec2rad (this is an AstroLib.jl's function).

c = t*(20043.109 - st*(85.33 + st*0.217) + t*(-42.665 - st*0.217 - t*41.8))/3600
end
pole_ra = zero(T)
pole_dec = one(T)*90

This comment has been minimized.

@giordano

giordano Jun 5, 2017

Member

I told you to use one(T) * 90, but T(90) is actually better. Sorry for the wrong suggestion before.

c = t*(20046.85 - st*(85.33 + st*0.37) + t*(-42.67 - st*0.37 - t*41.8))/3600
else
st = (epoch1 - 2000)*0.001
c = t*(20043.109 - st*(85.33 + st*0.217) + t*(-42.665 - st*0.217 - t*41.8))/3600

This comment has been minimized.

@giordano

giordano Jun 5, 2017

Member

Same as at line 9.

Final fixes
Test cases values found to be in an acceptable range with those received from precess_cd.pro (using GDL - GNU Data Language [Version 0.9.4])
c = sec2rad(t*(20043.109 - st*(85.33 + st*0.217) + t*(-42.665 - st*0.217 - t*41.8)))
end
pole_ra = zero(T)
pole_dec = one(90)

This comment has been minimized.

@giordano

giordano Jun 7, 2017

Member

T(90), not one(90) :-) (you have to adjust the example and tests)

This comment has been minimized.

@TestSubjector

TestSubjector Jun 7, 2017

Contributor

Whoops! Fixed now

TestSubjector added some commits Jun 7, 2017

Replaced one(90) with T(90)
* Fixed silly mistake
* Updated example and tests
5107.504395433958 6887.55936949333 ]
@test precess_cd([30.0 28.967; 60.45 90.65], 2000, 1975, [13, 10.658], [35.54, 67], true) ≈
[64.78429186351575 62.637156996728194 ;
130.49379143419722 195.9699513801844 ]

This comment has been minimized.

@kbarbary

kbarbary Jun 7, 2017

Member

Were these values generated by running the Julia code or by the IDL version? I think we're really going to want to have tests against the IDL version output. But at the very least, there should be a comment here explaining where the "truth" values come from.

This comment has been minimized.

@giordano

giordano Jun 7, 2017

Member

I compared some of the results with IDL AstroLib. There's a difference on the leas significant digits, but it happens the same with bprecess.

I don't think it's a good idea to always compare with results from IDL AstroLib, there are often 32- and 64-bit floating point numbers mixed together, so that they're often numerically unreliable if one expects full 64-bit accuracy

This comment has been minimized.

@kbarbary

kbarbary Jun 7, 2017

Member

One could then just decrease the required precision for isapprox, right? I think a test against an independent implementation is more valuable than a test of this implementation against itself.

Or we could have both tests, or maybe just put the IDL results in a comment.

This comment has been minimized.

@giordano

giordano Jun 7, 2017

Member

One could then just decrease the required precision for isapprox, right?

Yes, and for the specific case probably this isn't even needed, because rtol defaults to sqrt(eps(T)) (pretty high, IMHO) and atol to 0, I think the results will be well consistent within 1 part in ~10^8.

I think a test against an independent implementation is more valuable than a test of this implementation against itself.

If you don't know the "truth", testing an implementation against itself doesn't ensure the code is doing the Right Thing, but at least checks that nothing will be broken by future changes.

This comment has been minimized.

@TestSubjector

TestSubjector Jun 8, 2017

Contributor

Or we could have both tests, or maybe just put the IDL results in a comment.

Commenting the IDL results seems doable. I could work on it during and after the missing astronomical routines are done.

@giordano giordano merged commit c85a71a into JuliaAstro:master Jun 9, 2017

5 checks passed

codecov/patch 100% of diff hit (target 99.71%)
Details
codecov/project Absolute coverage decreased by -0.17% but relative coverage increased by +0.28% compared to c645430
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details
coverage/coveralls First build on master at 99.543%
Details
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment