-
Notifications
You must be signed in to change notification settings - Fork 21
/
bprecess.jl
217 lines (181 loc) · 8.24 KB
/
bprecess.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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.
const Mbprec =
reshape([+0.9999256795, -0.0111814828, -0.0048590040,
-0.000551, -0.238560, +0.435730,
+0.0111814828, +0.9999374849, -0.0000271557,
+0.238509, -0.002667, -0.008541,
+0.0048590039, -0.0000271771, +0.9999881946,
-0.435614, +0.012254, +0.002117,
-0.00000242389840, +0.00000002710544, +0.00000001177742,
+0.99990432, -0.01118145, -0.00485852,
-0.00000002710544, -0.00000242392702, +0.00000000006585,
+0.01118145, +0.99991613, -0.00002716,
-0.00000001177742, +0.00000000006585, -0.00000242404995,
+0.00485852, -0.00002717, +0.99996684], 6, 6)
# Note: IDL version of `bprecess' changes in-place "muradec", "parallax" and
# "radvel". We don't do anything like this, but calculations are below,
# commented, in case someone is interested.
function _bprecess{T<:AbstractFloat}(ra::T, dec::T, parallax::T,
radvel::T, epoch::T, muradec::Vector{T})
@assert length(muradec) == 2
cosra = cosd(ra)
sinra = sind(ra)
cosdec = cosd(dec)
sindec = sind(dec)
ra1950 = dec1950 = zero(T)
A = copy(A_precess)
r0 = [cosra*cosdec, sinra*cosdec, sindec]
r0_dot = [-muradec[1]*sinra*cosdec - muradec[2]*cosra*sindec,
muradec[1]*cosra*cosdec - muradec[2]*sinra*sindec,
muradec[2]*cosdec] + 21.095*radvel*parallax*r0
R_0 = vcat(r0, r0_dot)
R_1 = Mbprec*R_0
r1 = R_1[1:3]
r1_dot = R_1[4:6]
if isfinite(epoch)
r1 += deg2rad.(r1_dot*(epoch - 1950.0)*inv(360000.0))
A += deg2rad.(A_dot_precess*(epoch - 1950.0)*inv(360000.0))
end
rmag = vecnorm(r1)
s1 = r1*inv(rmag)
s1_dot = r1_dot*inv(rmag)
r = Array{T}(3)
s = s1
for j = 0:2
r = s1 + A - sum(s.*A)*s
s = r*inv(rmag)
end
rmag = vecnorm(r)
# r_dot = s1_dot + A_dot_precess - sum(s.*A_dot_precess)*s
x = r[1]
y = r[2]
z = r[3]
# x_dot = r_dot[1]
# y_dot = r_dot[2]
# z_dot = r_dot[3]
# muradec[1] = (x*y_dot - y*x_dot)/(x*x + y*y)
# muradec[2] = (z_dot*(x*x + y*y) - z*(x*x_dot + y*y_dot))/(rmag*rmag*hypot(x, y))
dec1950 = asin(z*inv(rmag))
ra1950 = atan2(y, x)
# if parallax > 0
# radvel = (x*x_dot + y*y_dot + z*z_dot)/(21.095*parallax*rmag)
# parallax = parallax*inv(rmag)
# end
if ra1950 < 0
ra1950 += tau
end
ra1950 = rad2deg(ra1950)
dec1950 = rad2deg(dec1950)
return ra1950, dec1950
end
# Main interface.
bprecess{R<:Real}(ra::Real, dec::Real, muradec::Vector{R};
parallax::Real=0.0, radvel::Real=0.0) =
_bprecess(promote(float(ra), float(dec), float(parallax),
float(radvel), NaN)...,
float(muradec))
bprecess(ra::Real, dec::Real, epoch::Real=2000.0) =
_bprecess(promote(float(ra), float(dec), 0.0, 0.0, float(epoch))...,
zeros(typeof(float(ra)), 2))
# Tuple arguments.
bprecess{R1<:Real,R2<:Real,R3<:Real}(radec::Tuple{R1,R2}, muradec::Vector{R3};
parallax::Real=0.0, radvel::Real=0.0) =
bprecess(radec..., muradec,
parallax=parallax,
radvel=radvel)
bprecess{R1<:Real,R2<:Real}(radec::Tuple{R1,R2}, epoch::Real=2000.0) =
bprecess(radec..., epoch)
# Vectorial arguments.
function bprecess{R<:Real,D<:Real,M<:Real,P<:Real,V<:Real}(ra::AbstractArray{R},
dec::AbstractArray{D},
muradec::AbstractArray{M};
parallax::AbstractArray{P}=zeros(R, length(ra)),
radvel::AbstractArray{V}=zeros(R, length(ra)))
@assert length(ra) == length(dec) == size(muradec)[2] == length(parallax) == length(radvel)
typer = typeof(float(one(R)))
ra1950 = similar(ra, typer)
dec1950 = similar(dec, typer)
for i in eachindex(ra)
ra1950[i], dec1950[i] = bprecess(ra[i], dec[i], muradec[:,i],
parallax=parallax[i], radvel=radvel[i])
end
return ra1950, dec1950
end
function bprecess{R<:Real,D<:Real}(ra::AbstractArray{R},
dec::AbstractArray{D},
epoch::Real=2000.0)
@assert length(ra) == length(dec)
typer = float(R)
ra1950 = similar(ra, typer)
dec1950 = similar(dec, typer)
for i in eachindex(ra)
ra1950[i], dec1950[i] = bprecess(ra[i], dec[i], epoch[i])
end
return ra1950, dec1950
end
"""
bprecess(ra, dec[, epoch]) -> ra1950, dec1950
bprecess(ra, dec, muradec[, parallax=parallax, radvel=radvel]) -> ra1950, dec1950
### Purpose ###
Precess positions from J2000.0 (FK5) to B1950.0 (FK4).
### Explanation ###
Calculates the mean place of a star at B1950.0 on the FK4 system from the mean
place at J2000.0 on the FK5 system.
`bprecess` function has two methods, one for each of the following cases:
* the proper motion is known and non-zero
* the proper motion is unknown or known to be exactly zero (i.e. extragalactic
radio sources). Better precision can be achieved in this case by inputting
the epoch of the original observations.
### Arguments ###
The function has 2 methods. The common mandatory arguments are:
* `ra`: input J2000 right ascension, in degrees.
* `dec`: input J2000 declination, in degrees.
The two methods have a different third argument (see "Explanation" section for
more details). It can be one of the following:
* `muradec`: 2-element vector containing the proper motion in seconds of arc
per tropical *century* in right ascension and declination.
* `epoch`: scalar giving epoch of original observations.
If none of these two arguments is provided (so `bprecess` is fed only with right
ascension and declination), it is assumed that proper motion is exactly zero and
`epoch = 2000`.
If it is used the method involving `muradec` argument, the following keywords
are available:
* `parallax` (optional numerical keyword): stellar parallax, in seconds of arc.
* `radvel` (optional numerical keyword): radial velocity in km/s.
Right ascension and declination can be passed as the 2-tuple `(ra, dec)`. You
can also pass `ra`, `dec`, `parallax`, and `radvel` as arrays, all of the same
length N. In that case, `muradec` should be a matrix 2×N.
### Output ###
The 2-tuple of right ascension and declination in 1950, in degrees, of input
coordinates is returned. If `ra` and `dec` (and other possible optional
arguments) are arrays, the 2-tuple of arrays `(ra1950, dec1950)` of the same
length as the input coordinates is returned.
### Method ###
The algorithm is taken from the Explanatory Supplement to the Astronomical
Almanac 1992, page 186. See also Aoki et al (1983), A&A, 128, 263. URL:
http://adsabs.harvard.edu/abs/1983A%26A...128..263A.
### Example ###
The SAO2000 catalogue gives the J2000 position and proper motion for the star HD
119288. Find the B1950 position.
* RA(2000) = 13h 42m 12.740s
* Dec(2000) = 8d 23' 17.69''
* Mu(RA) = -.0257 s/yr
* Mu(Dec) = -.090 ''/yr
``` julia
muradec = 100*[-15*0.0257, -0.090]; # convert to century proper motion
ra = ten(13, 42, 12.74)*15;
decl = ten(8, 23, 17.69);
adstring(bprecess(ra, decl, muradec), precision=2)
# => " 13 39 44.526 +08 38 28.63"
```
### Notes ###
"When transferring individual observations, as opposed to catalog mean place,
the safest method is to transform the observations back to the epoch of the
observation, on the FK4 system (or in the system that was used to to produce the
observed mean place), convert to the FK5 system, and transform to the the epoch
and equinox of J2000.0" -- from the Explanatory Supplement (1992), p. 180
`jprecess` performs the precession to J2000 coordinates.
Code of this function is based on IDL Astronomy User's Library.
"""
bprecess