/
precess_cd.jl
78 lines (61 loc) · 2.71 KB
/
precess_cd.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
# This file is a part of AstroLib.jl. License is MIT "Expat".
function precess_cd(cd::AbstractMatrix{T}, epoch1::T, epoch2::T, crval_old::AbstractVector{T},
crval_new::AbstractVector{T}, FK4::Bool) where {T<:AbstractFloat}
t = (epoch2 - epoch1) / 1000
if FK4
st = (epoch1 - 1900) / 1000
c = sec2rad(t*(20046.85 - st*(85.33 + st*0.37) + t*(-42.67 - st*0.37 - t*41.8)))
else
st = (epoch1 - 2000) / 1000
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 = 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, cosd1 = sincos(deg2rad(crval_old[2]))
sind2, cosd2 = sincos(deg2rad(crval_new[2]))
sinra = sind(crval_new[1] - pole_ra)
sin_c, cos_c = sincos(c)
cosfi = (cos_c - sind1 * sind2) / (cosd1 * cosd2)
sinfi = (abs(sin_c) * sinra) / cosd1
return cd * SMatrix{2,2}(cosfi, -sinfi, sinfi, cosfi)
end
"""
precess_cd(cd, epoch1, epoch2, crval_old, crval_new[, FK4=true]) -> 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
* `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 to `true`,
then the precession constants are taken in the FK4 reference frame. When it
is `false`, the default is the FK5 frame
### Output ###
* `cd`: coordinate description containing precessed values
### Example ###
```jldoctest
julia> using AstroLib
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 [`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) =
precess_cd(float(cd), promote(float(epoch1), float(epoch2))...,
float(crval_old), float(crval_new), FK4)