From 90f8767c066409ca5174fa93a1da7596624dc801 Mon Sep 17 00:00:00 2001 From: TestSubjector Date: Sun, 4 Jun 2017 00:30:43 +0530 Subject: [PATCH 1/8] 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. --- TODO.md | 1 - docs/src/ref.md | 1 + src/precess_cd.jl | 78 +++++++++++++++++++++++++++++++++++++++++++++ src/utils.jl | 3 ++ test/utils-tests.jl | 16 ++++++++++ 5 files changed, 98 insertions(+), 1 deletion(-) create mode 100644 src/precess_cd.jl diff --git a/TODO.md b/TODO.md index 6d5c8c7..d73b72f 100644 --- a/TODO.md +++ b/TODO.md @@ -78,7 +78,6 @@ Missing in AstroLib.jl * `imf` * `ismeuv` * `planet_coords` -* `precess_cd` * `qdcb_grid` * `tdb2tdt` * `ticlabels` diff --git a/docs/src/ref.md b/docs/src/ref.md index 9663375..d169caa 100644 --- a/docs/src/ref.md +++ b/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), diff --git a/src/precess_cd.jl b/src/precess_cd.jl new file mode 100644 index 0000000..ff471d1 --- /dev/null +++ b/src/precess_cd.jl @@ -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}, + crval_new::Vector{T}, FK4::Bool) + crvalold = deg2rad.(crval_old) + crvalnew = deg2rad.(crval_new) + 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 = 0 + pole_dec = 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 = sin.([crvalold[2], crvalnew[2]]) + cosd1, cosd2 = cos.([crvalold[2], crvalnew[2]]) + 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]] + cd = cd * r + return cd +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}: + 49.1292 146.996 + 110.365 110.188 +``` + +### 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{R<:Real}(cd::Matrix{R}, epoch1::Real, epoch2::Real, crval_old::Vector{R}, + crval_new::Vector{R}, FK4::Bool=false) = + _precess_cd(float(cd), promote(float(epoch1), float(epoch2))..., + float(crval_old), float(crval_new), FK4) diff --git a/src/utils.jl b/src/utils.jl index 3cc1321..ef8a290 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -121,6 +121,9 @@ export posang include("precess.jl") export precess +include("precess_cd") +export precess_cd + include("precess_xyz.jl") export precess_xyz diff --git a/test/utils-tests.jl b/test/utils-tests.jl index 5d67afb..57cde8e 100644 --- a/test/utils-tests.jl +++ b/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]) ≈ + [31.323053098892704 62.14104826007457; + 62.54509461024324 93.16106659634077 ] + @test precess_cd([[30 60]; [60 90]], 2000, 1950, [13, 8], [43, 23]) ≈ + [31.323011022652295 62.14106970440559; + 62.54503157712479 93.16110932251847 ] + @test precess_cd([[12.45 56.7]; [66 89]], 2000, 1985, [67.4589455, 0.345345], [37.94291666666666, 89.26405555555556]) ≈ + [963.5270623437744 4387.868087815929 ; + 5107.6642727783155 6887.440809471227] + @test precess_cd([[12.45 56.7]; [66 89]], 2000, 1985, [67.4589455, 0.345345], [37.94291666666666, 89.26405555555556]) ≈ + [963.5270506245276 4387.868092741147 ; + 5107.664256221472 6887.440827208123 ] +end + # Test precess_xyz let local x1 ,y1, z1, x2, y2, z2 From 138addab2d7f9bf58dee453be8c4cf8f61bf1a3c Mon Sep 17 00:00:00 2001 From: TestSubjector Date: Sun, 4 Jun 2017 00:48:38 +0530 Subject: [PATCH 2/8] Quick file name fix --- src/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index ef8a290..ca4a2a6 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -121,7 +121,7 @@ export posang include("precess.jl") export precess -include("precess_cd") +include("precess_cd.jl") export precess_cd include("precess_xyz.jl") From a095ad3208a8a4e37973855fd36a3fe47e656433 Mon Sep 17 00:00:00 2001 From: TestSubjector Date: Sun, 4 Jun 2017 01:15:42 +0530 Subject: [PATCH 3/8] Fix test for FK4 frame --- test/utils-tests.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/utils-tests.jl b/test/utils-tests.jl index 57cde8e..a7c3df0 100644 --- a/test/utils-tests.jl +++ b/test/utils-tests.jl @@ -452,9 +452,9 @@ end @test precess_cd([[12.45 56.7]; [66 89]], 2000, 1985, [67.4589455, 0.345345], [37.94291666666666, 89.26405555555556]) ≈ [963.5270623437744 4387.868087815929 ; 5107.6642727783155 6887.440809471227] - @test precess_cd([[12.45 56.7]; [66 89]], 2000, 1985, [67.4589455, 0.345345], [37.94291666666666, 89.26405555555556]) ≈ - [963.5270506245276 4387.868092741147 ; - 5107.664256221472 6887.440827208123 ] + @test precess_cd([[30.0 28.967]; [60.45 90.65]], 2000, 1975, [13, 10.658], [35.54, 67], true) ≈ + [64.86784270539363 62.55062667282337 ; + 130.7552573456406 195.79559280283118] end # Test precess_xyz From 221b889d2c966f7bed0b2569e51dbe37a5f0dccb Mon Sep 17 00:00:00 2001 From: TestSubjector Date: Mon, 5 Jun 2017 06:16:16 +0530 Subject: [PATCH 4/8] Improvements and Correcting matrix r --- src/precess_cd.jl | 35 ++++++++++++++++------------------- test/utils-tests.jl | 24 ++++++++++++------------ 2 files changed, 28 insertions(+), 31 deletions(-) diff --git a/src/precess_cd.jl b/src/precess_cd.jl index ff471d1..d641407 100644 --- a/src/precess_cd.jl +++ b/src/precess_cd.jl @@ -1,9 +1,7 @@ # 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}, - crval_new::Vector{T}, FK4::Bool) - crvalold = deg2rad.(crval_old) - crvalnew = deg2rad.(crval_new) +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 @@ -13,22 +11,21 @@ function _precess_cd{T<:AbstractFloat}(cd::Matrix{T}, epoch1::T, epoch2::T, crva 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 = 0 - pole_dec = 90 + 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 = sin.([crvalold[2], crvalnew[2]]) - cosd1, cosd2 = cos.([crvalold[2], crvalnew[2]]) - sinra = sin(crvalnew[1] - deg2rad.(pole_ra)) + 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]] - cd = cd * r - return cd + r = [cosfi sinfi; -sinfi cosfi] + return cd * r end """ @@ -60,10 +57,10 @@ The coordinate matrix is precessed from epoch1 to epoch2. ### Example ### ```julia -julia> precess_cd([[20 60]; [45 45]], 1950, 2000, [34, 58], [12, 83]) +julia> precess_cd([20 60; 45 45], 1950, 2000, [34, 58], [12, 83]) 2×2 Array{Float64,2}: - 49.1292 146.996 - 110.365 110.188 + 48.8944 147.075 + 110.188 110.365 ``` ### Notes ### @@ -72,7 +69,7 @@ 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{R<:Real}(cd::Matrix{R}, epoch1::Real, epoch2::Real, crval_old::Vector{R}, - crval_new::Vector{R}, FK4::Bool=false) = - _precess_cd(float(cd), promote(float(epoch1), float(epoch2))..., - float(crval_old), float(crval_new), FK4) +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) diff --git a/test/utils-tests.jl b/test/utils-tests.jl index a7c3df0..3826d9e 100644 --- a/test/utils-tests.jl +++ b/test/utils-tests.jl @@ -443,18 +443,18 @@ end # Test precess_cd @testset "precess_cd" begin - @test precess_cd([[30 60]; [60 90]], 1950, 2000, [13, 8], [43, 23]) ≈ - [31.323053098892704 62.14104826007457; - 62.54509461024324 93.16106659634077 ] - @test precess_cd([[30 60]; [60 90]], 2000, 1950, [13, 8], [43, 23]) ≈ - [31.323011022652295 62.14106970440559; - 62.54503157712479 93.16110932251847 ] - @test precess_cd([[12.45 56.7]; [66 89]], 2000, 1985, [67.4589455, 0.345345], [37.94291666666666, 89.26405555555556]) ≈ - [963.5270623437744 4387.868087815929 ; - 5107.6642727783155 6887.440809471227] - @test precess_cd([[30.0 28.967]; [60.45 90.65]], 2000, 1975, [13, 10.658], [35.54, 67], true) ≈ - [64.86784270539363 62.55062667282337 ; - 130.7552573456406 195.79559280283118] + @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 From bfae37c893383465a0971bc384d95d8ac5f505df Mon Sep 17 00:00:00 2001 From: TestSubjector Date: Tue, 6 Jun 2017 20:18:23 +0530 Subject: [PATCH 5/8] 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]) --- src/precess_cd.jl | 20 +++++++++++--------- test/utils-tests.jl | 16 ++++++++-------- 2 files changed, 19 insertions(+), 17 deletions(-) diff --git a/src/precess_cd.jl b/src/precess_cd.jl index d641407..dbb330f 100644 --- a/src/precess_cd.jl +++ b/src/precess_cd.jl @@ -2,25 +2,27 @@ 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) + t = (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 + c = sec2rad(t*(20046.85 - st*(85.33 + st*0.37) + t*(-42.67 - st*0.37 - t*41.8))) 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 + 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(T)*90 + pole_dec = one(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]]) + sind1 = sind(crval_old[2]) + sind2 = sind(crval_new[2]) + cosd1 = cosd(crval_old[2]) + cosd2 = cosd(crval_new[2]) sinra = sind(crval_new[1] - pole_ra) cosfi = (cos(c) - sind1*sind2)/(cosd1*cosd2) sinfi = (abs(sin(c))* sinra)/cosd1 @@ -41,7 +43,7 @@ The coordinate matrix is precessed from epoch1 to epoch2. ### Arguments ### -* `cd`: 2 x 2 coordinate description matrix in degrees or radians +* `cd`: 2 x 2 coordinate description matrix in degrees * `epoch1`: original equinox of coordinates, scalar * `epoch2`: equinox of precessed coordinates, scalar * `crval_old`: 2 element vector containing right ascension and declination @@ -59,8 +61,8 @@ The coordinate matrix is precessed from epoch1 to epoch2. ```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 + 48.8914 147.076 + 110.186 110.367 ``` ### Notes ### diff --git a/test/utils-tests.jl b/test/utils-tests.jl index 3826d9e..60842a6 100644 --- a/test/utils-tests.jl +++ b/test/utils-tests.jl @@ -444,17 +444,17 @@ 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 ] + [30.917848605585117 62.343650719942836; + 61.93728791292476 93.56627143079152 ] @test precess_cd([30 60; 60 90], 2000, 1950, [13, 8], [43, 23]) ≈ - [30.919049149933095 62.34305064076519 ; - 61.93908876804599 93.56507119523768 ] + [30.917848605585117 62.343650719942836; + 61.93728791292476 93.56627143079152 ] @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 ] + [963.4250989585853 4387.890476297739 ; + 5107.5042241937 6887.559496480489 ] @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 ] + [64.78412957049132 62.63732507732793 ; + 130.49328355065168 195.97029006259285] end # Test precess_xyz From 520dd90ce452769dda2f3be50df43f0c12d61e8c Mon Sep 17 00:00:00 2001 From: TestSubjector Date: Wed, 7 Jun 2017 16:51:38 +0530 Subject: [PATCH 6/8] Replaced one(90) with T(90) * Fixed silly mistake * Updated example and tests --- src/precess_cd.jl | 6 +++--- test/utils-tests.jl | 16 ++++++++-------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/src/precess_cd.jl b/src/precess_cd.jl index dbb330f..9561de5 100644 --- a/src/precess_cd.jl +++ b/src/precess_cd.jl @@ -12,7 +12,7 @@ function _precess_cd{T<:AbstractFloat}(cd::AbstractMatrix{T}, epoch1::T, epoch2: 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) + pole_dec = T(90) if (epoch1 == 2000 && epoch2 == 1950) || (epoch1 == 1950 && epoch2 == 2000) pole_ra, pole_dec = bprecess(pole_ra, pole_dec) @@ -61,8 +61,8 @@ The coordinate matrix is precessed from epoch1 to epoch2. ```julia julia> precess_cd([20 60; 45 45], 1950, 2000, [34, 58], [12, 83]) 2×2 Array{Float64,2}: - 48.8914 147.076 - 110.186 110.367 + 48.8944 147.075 + 110.188 110.365 ``` ### Notes ### diff --git a/test/utils-tests.jl b/test/utils-tests.jl index 60842a6..00faeb4 100644 --- a/test/utils-tests.jl +++ b/test/utils-tests.jl @@ -444,17 +444,17 @@ end # Test precess_cd @testset "precess_cd" begin @test precess_cd([30 60; 60 90], 1950, 2000, [13, 8], [43, 23]) ≈ - [30.917848605585117 62.343650719942836; - 61.93728791292476 93.56627143079152 ] + [30.919029003435927 62.343060521017435; + 61.93905850970097 93.56509103294071 ] @test precess_cd([30 60; 60 90], 2000, 1950, [13, 8], [43, 23]) ≈ - [30.917848605585117 62.343650719942836; - 61.93728791292476 93.56627143079152 ] + [30.919029003435927 62.343060521017435; + 61.93905850970097 93.56509103294071 ] @test precess_cd([12.45 56.7; 66 89], 2000, 1985, [67.4589455, 0.345345], [37.94291666666666, 89.26405555555556]) ≈ - [963.4250989585853 4387.890476297739 ; - 5107.5042241937 6887.559496480489 ] + [963.4252080520984 4387.890452343343 ; + 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.78412957049132 62.63732507732793 ; - 130.49328355065168 195.97029006259285] + [64.78429186351575 62.637156996728194 ; + 130.49379143419722 195.9699513801844 ] end # Test precess_xyz From 2e611f9cb89c499794bef70b2173cddfb44f82aa Mon Sep 17 00:00:00 2001 From: TestSubjector Date: Wed, 7 Jun 2017 19:02:25 +0530 Subject: [PATCH 7/8] Add reference to sec2rad --- src/precess_cd.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/precess_cd.jl b/src/precess_cd.jl index 9561de5..9576e95 100644 --- a/src/precess_cd.jl +++ b/src/precess_cd.jl @@ -69,7 +69,7 @@ julia> precess_cd([20 60; 45 45], 1950, 2000, [34, 58], [12, 83]) 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). +This function calls [sec2rad](@ref), [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) = From 1c65feb204c0a79f63e2aba54c6982b02b1c9dc5 Mon Sep 17 00:00:00 2001 From: TestSubjector Date: Fri, 9 Jun 2017 14:34:11 +0530 Subject: [PATCH 8/8] Comment information about test values --- test/utils-tests.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/utils-tests.jl b/test/utils-tests.jl index 00faeb4..bd59272 100644 --- a/test/utils-tests.jl +++ b/test/utils-tests.jl @@ -442,6 +442,9 @@ let end # Test precess_cd +# The values used for the testset are from running the code. However they have been +# correlated with the output from precess_cd routine of IDL AstroLib, with +# differences only in the least significant digits. @testset "precess_cd" begin @test precess_cd([30 60; 60 90], 1950, 2000, [13, 8], [43, 23]) ≈ [30.919029003435927 62.343060521017435;