Add precess_cd function

TODO.md

@@ -78,7 +78,6 @@ Missing in AstroLib.jl
 * `imf`
 * `ismeuv`
 * `planet_coords`
-* `precess_cd`
 * `qdcb_grid`
 * `tdb2tdt`
 * `ticlabels`

docs/src/ref.md

@@ -106,6 +106,7 @@ julia> AstroLib.planets["saturn"].mass
 [`polrec()`](@ref),
 [`posang()`](@ref),
 [`precess()`](@ref),
+[`precess_cd()`](@ref),
 [`precess_xyz()`](@ref),
 [`premat()`](@ref),
 [`radec()`](@ref),

src/precess_cd.jl

@@ -0,0 +1,75 @@
+# This file is a part of AstroLib.jl. License is MIT "Expat".
+
+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)
+
+    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
+    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
+    end
+    pole_ra = zero(T)
+    pole_dec = one(T)*90
+
+    if (epoch1 == 2000 && epoch2 == 1950) || (epoch1 == 1950 && epoch2 == 2000)
+        pole_ra, pole_dec = bprecess(pole_ra, pole_dec)
+    else
+        pole_ra, pole_dec = precess(pole_ra, pole_dec, epoch1, epoch2, FK4=FK4)
+    end
+    sind1, sind2 = sind.([crval_old[2], crval_new[2]])
+    cosd1, cosd2 = cosd.([crval_old[2], crval_new[2]])
+    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]
+    return cd * r
+end
+
+"""
+    precess_cd(cd, epoch1, epoch2, crval_old, crval_new) -> cd
+
+### Purpose ###
+
+Precess the coordinate description matrix.
+
+### Explanation ###
+
+The coordinate matrix is precessed from epoch1 to epoch2.
+
+### Arguments ###
+
+* `cd`: 2 x 2 coordinate description matrix in degrees or radians
+* `epoch1`: original equinox of coordinates, scalar
+* `epoch2`: equinox of precessed coordinates, scalar
+* `crval_old`: 2 element vector containing right ascension and declination
+  in degrees of the reference pixel in the original equinox
+* `crval_new`: 2 element vector giving crval in the new equinox
+* `FK4` (optional boolean keyword): if this keyword is set, then the precession constants
+  are taken in the FK4 reference frame. The default is the FK5 frame
+
+### Output ###
+
+* `cd`: coordinate description containing precessed values
+
+### Example ###
+
+```julia
+julia> precess_cd([20 60; 45 45], 1950, 2000, [34, 58], [12, 83])
+2×2 Array{Float64,2}:
+  48.8944  147.075
+ 110.188   110.365
+```
+
+### Notes ###
+
+Code of this function is based on IDL Astronomy User's Library.
+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(cd::AbstractMatrix{<:Real}, epoch1::Real, epoch2::Real, crval_old::AbstractVector{<:Real},
+           crval_new::AbstractVector{<:Real}, FK4::Bool=false) =
+               _precess_cd(float(cd), promote(float(epoch1), float(epoch2))...,
+                           float(crval_old), float(crval_new), FK4)

src/utils.jl

@@ -121,6 +121,9 @@ export posang
 include("precess.jl")
 export precess
 
+include("precess_cd.jl")
+export precess_cd
+
 include("precess_xyz.jl")
 export precess_xyz
 

test/utils-tests.jl

@@ -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]) ≈
+        [30.919006748724033 62.343071435158905;
+         61.939025084990234 93.56511294650944 ]
+    @test precess_cd([30 60; 60 90], 2000, 1950, [13, 8], [43, 23]) ≈
+        [30.919049149933095 62.34305064076519 ;
+         61.93908876804599 93.56507119523768  ]
+    @test precess_cd([12.45 56.7; 66 89], 2000, 1985, [67.4589455, 0.345345], [37.94291666666666, 89.26405555555556]) ≈
+        [963.4252096895427 4387.890452287613  ;
+         5107.50439824168 6887.559368116342   ]
+    @test precess_cd([30.0 28.967; 60.45 90.65], 2000, 1975, [13, 10.658], [35.54, 67], true) ≈
+        [64.78429401954816 62.637154810889584 ;
+         130.4937981551171 195.9699470010346  ]
+end
+
 # Test precess_xyz
 let
     local x1 ,y1, z1, x2, y2, z2