1 | ;
|
2 |
|
3 | Object.defineProperty(exports, "__esModule", {
|
4 | value: true
|
5 | });
|
6 |
|
7 | var _slicedToArray = function () { function sliceIterator(arr, i) { var _arr = []; var _n = true; var _d = false; var _e = undefined; try { for (var _i = arr[Symbol.iterator](), _s; !(_n = (_s = _i.next()).done); _n = true) { _arr.push(_s.value); if (i && _arr.length === i) break; } } catch (err) { _d = true; _e = err; } finally { try { if (!_n && _i["return"]) _i["return"](); } finally { if (_d) throw _e; } } return _arr; } return function (arr, i) { if (Array.isArray(arr)) { return arr; } else if (Symbol.iterator in Object(arr)) { return sliceIterator(arr, i); } else { throw new TypeError("Invalid attempt to destructure non-iterable instance"); } }; }(); /**
|
8 | * @copyright 2013 Sonia Keys
|
9 | * @copyright 2016 commenthol
|
10 | * @license MIT
|
11 | * @module angle
|
12 | */
|
13 | /**
|
14 | * Angle: Chapter 17: Angular Separation.
|
15 | *
|
16 | * Functions in this package are useful for Ecliptic, Equatorial, or any
|
17 | * similar coordinate frame. To avoid suggestion of a particular frame,
|
18 | * function parameters are specified simply as "r1, d1" to correspond to a
|
19 | * right ascenscion, declination pair or to a longitude, latitude pair.
|
20 | *
|
21 | * In function Sep, Meeus recommends 10 arc min as a threshold. This
|
22 | * value is in package base as base.SmallAngle because it has general utility.
|
23 | *
|
24 | * All angles are in radians.
|
25 | */
|
26 |
|
27 | exports.sep = sep;
|
28 | exports.minSep = minSep;
|
29 | exports.minSepRect = minSepRect;
|
30 | exports.hav = hav;
|
31 | exports.sepHav = sepHav;
|
32 | exports.minSepHav = minSepHav;
|
33 | exports.sepPauwels = sepPauwels;
|
34 | exports.minSepPauwels = minSepPauwels;
|
35 | exports.relativePosition = relativePosition;
|
36 |
|
37 | var _base = require('./base');
|
38 |
|
39 | var _base2 = _interopRequireDefault(_base);
|
40 |
|
41 | var _interpolation = require('./interpolation');
|
42 |
|
43 | var _interpolation2 = _interopRequireDefault(_interpolation);
|
44 |
|
45 | function _interopRequireDefault(obj) { return obj && obj.__esModule ? obj : { default: obj }; }
|
46 |
|
47 | var abs = Math.abs,
|
48 | acos = Math.acos,
|
49 | asin = Math.asin,
|
50 | atan2 = Math.atan2,
|
51 | cos = Math.cos,
|
52 | hypot = Math.hypot,
|
53 | sin = Math.sin,
|
54 | sqrt = Math.sqrt,
|
55 | tan = Math.tan;
|
56 |
|
57 | /**
|
58 | * `sep` returns the angular separation between two celestial bodies.
|
59 | *
|
60 | * The algorithm is numerically naïve, and while patched up a bit for
|
61 | * small separations, remains unstable for separations near π.
|
62 | *
|
63 | * @param {base.Coord} c1 - coordinate of celestial body 1
|
64 | * @param {base.Coord} c2 - coordinate of celestial body 2
|
65 | * @return {Number} angular separation between two celestial bodies
|
66 | */
|
67 |
|
68 | function sep(c1, c2) {
|
69 | var _base$sincos = _base2.default.sincos(c1.dec),
|
70 | _base$sincos2 = _slicedToArray(_base$sincos, 2),
|
71 | sind1 = _base$sincos2[0],
|
72 | cosd1 = _base$sincos2[1];
|
73 |
|
74 | var _base$sincos3 = _base2.default.sincos(c2.dec),
|
75 | _base$sincos4 = _slicedToArray(_base$sincos3, 2),
|
76 | sind2 = _base$sincos4[0],
|
77 | cosd2 = _base$sincos4[1];
|
78 |
|
79 | var cd = sind1 * sind2 + cosd1 * cosd2 * cos(c1.ra - c2.ra); // (17.1) p. 109
|
80 | if (cd < _base2.default.CosSmallAngle) {
|
81 | return acos(cd);
|
82 | }
|
83 | return hypot((c2.ra - c1.ra) * cosd1, c2.dec - c1.dec); // (17.2) p. 109
|
84 | }
|
85 |
|
86 | /**
|
87 | * `minSep` returns the minimum separation between two moving objects.
|
88 | *
|
89 | * The motion is represented as an ephemeris of three rows, equally spaced
|
90 | * in time. Jd1, jd3 are julian day times of the first and last rows.
|
91 | * R1, d1, r2, d2 are coordinates at the three times. They must each be
|
92 | * slices of length 3.0
|
93 | *
|
94 | * Result is obtained by computing separation at each of the three times
|
95 | * and interpolating a minimum. This may be invalid for sufficiently close
|
96 | * approaches.
|
97 | *
|
98 | * @throws Error
|
99 | * @param {Number} jd1 - Julian day - time at cs1[0], cs2[0]
|
100 | * @param {Number} jd3 - Julian day - time at cs1[2], cs2[2]
|
101 | * @param {base.Coord[]} cs1 - 3 coordinates of moving object 1
|
102 | * @param {base.Coord[]} cs2 - 3 coordinates of moving object 2
|
103 | * @param {function} [fnSep] - alternative `sep` function e.g. `angle.sepPauwels`, `angle.sepHav`
|
104 | * @return {Number} angular separation between two celestial bodies
|
105 | */
|
106 | function minSep(jd1, jd3, cs1, cs2, fnSep) {
|
107 | fnSep = fnSep || sep;
|
108 | if (cs1.length !== 3 || cs2.length !== 3) {
|
109 | throw _interpolation2.default.errorNot3;
|
110 | }
|
111 | var y = new Array(3);
|
112 | cs1.forEach(function (c, x) {
|
113 | y[x] = sep(cs1[x], cs2[x]);
|
114 | });
|
115 | var d3 = new _interpolation2.default.Len3(jd1, jd3, y);
|
116 | var dMin = d3.extremum()[1];
|
117 | return dMin;
|
118 | }
|
119 |
|
120 | /**
|
121 | * `minSepRect` returns the minimum separation between two moving objects.
|
122 | *
|
123 | * Like `minSep`, but using a method of rectangular coordinates that gives
|
124 | * accurate results even for close approaches.
|
125 | *
|
126 | * @throws Error
|
127 | * @param {Number} jd1 - Julian day - time at cs1[0], cs2[0]
|
128 | * @param {Number} jd3 - Julian day - time at cs1[2], cs2[2]
|
129 | * @param {base.Coord[]} cs1 - 3 coordinates of moving object 1
|
130 | * @param {base.Coord[]} cs2 - 3 coordinates of moving object 2
|
131 | * @return {Number} angular separation between two celestial bodies
|
132 | */
|
133 | function minSepRect(jd1, jd3, cs1, cs2) {
|
134 | if (cs1.length !== 3 || cs2.length !== 3) {
|
135 | throw _interpolation2.default.ErrorNot3;
|
136 | }
|
137 | var uv = function uv(c1, c2) {
|
138 | var _base$sincos5 = _base2.default.sincos(c1.dec),
|
139 | _base$sincos6 = _slicedToArray(_base$sincos5, 2),
|
140 | sind1 = _base$sincos6[0],
|
141 | cosd1 = _base$sincos6[1];
|
142 |
|
143 | var Δr = c2.ra - c1.ra;
|
144 | var tanΔr = tan(Δr);
|
145 | var tanhΔr = tan(Δr / 2);
|
146 | var K = 1 / (1 + sind1 * sind1 * tanΔr * tanhΔr);
|
147 | var sinΔd = sin(c2.dec - c1.dec);
|
148 | var u = -K * (1 - sind1 / cosd1 * sinΔd) * cosd1 * tanΔr;
|
149 | var v = K * (sinΔd + sind1 * cosd1 * tanΔr * tanhΔr);
|
150 | return [u, v];
|
151 | };
|
152 | var us = new Array(3).fill(0);
|
153 | var vs = new Array(3).fill(0);
|
154 | cs1.forEach(function (c, x) {
|
155 | var _uv = uv(cs1[x], cs2[x]);
|
156 |
|
157 | var _uv2 = _slicedToArray(_uv, 2);
|
158 |
|
159 | us[x] = _uv2[0];
|
160 | vs[x] = _uv2[1];
|
161 | });
|
162 | var u3 = new _interpolation2.default.Len3(-1, 1, us); // if line throws then bug not caller's fault.
|
163 | var v3 = new _interpolation2.default.Len3(-1, 1, vs); // if line throws then bug not caller's fault.
|
164 | var up0 = (us[2] - us[0]) / 2;
|
165 | var vp0 = (vs[2] - vs[0]) / 2;
|
166 | var up1 = us[0] + us[2] - 2 * us[1];
|
167 | var vp1 = vs[0] + vs[2] - 2 * vs[1];
|
168 | var up = up0;
|
169 | var vp = vp0;
|
170 | var dn = -(us[1] * up + vs[1] * vp) / (up * up + vp * vp);
|
171 | var n = dn;
|
172 | var u = void 0;
|
173 | var v = void 0;
|
174 | for (var limit = 0; limit < 10; limit++) {
|
175 | u = u3.interpolateN(n);
|
176 | v = v3.interpolateN(n);
|
177 | if (abs(dn) < 1e-5) {
|
178 | return hypot(u, v); // success
|
179 | }
|
180 | var _up = up0 + n * up1;
|
181 | var _vp = vp0 + n * vp1;
|
182 | dn = -(u * _up + v * _vp) / (_up * _up + _vp * _vp);
|
183 | n += dn;
|
184 | }
|
185 | throw new Error('minSepRect: failure to converge');
|
186 | }
|
187 |
|
188 | /**
|
189 | * haversine function (17.5) p. 115
|
190 | */
|
191 | function hav(a) {
|
192 | return 0.5 * (1 - Math.cos(a));
|
193 | }
|
194 |
|
195 | /**
|
196 | * `sepHav` returns the angular separation between two celestial bodies.
|
197 | *
|
198 | * The algorithm uses the haversine function and is superior to the naïve
|
199 | * algorithm of the Sep function.
|
200 | *
|
201 | * @param {base.Coord} c1 - coordinate of celestial body 1
|
202 | * @param {base.Coord} c2 - coordinate of celestial body 2
|
203 | * @return {Number} angular separation between two celestial bodies
|
204 | */
|
205 | function sepHav(c1, c2) {
|
206 | // using (17.5) p. 115
|
207 | return 2 * asin(sqrt(hav(c2.dec - c1.dec) + cos(c1.dec) * cos(c2.dec) * hav(c2.ra - c1.ra)));
|
208 | }
|
209 |
|
210 | /**
|
211 | * Same as `minSep` but uses function `sepHav` to return the minimum separation
|
212 | * between two moving objects.
|
213 | */
|
214 | function minSepHav(jd1, jd3, cs1, cs2) {
|
215 | return minSep(jd1, jd3, cs1, cs2, sepHav);
|
216 | }
|
217 |
|
218 | /**
|
219 | * `sepPauwels` returns the angular separation between two celestial bodies.
|
220 | *
|
221 | * The algorithm is a numerically stable form of that used in `sep`.
|
222 | *
|
223 | * @param {base.Coord} c1 - coordinate of celestial body 1
|
224 | * @param {base.Coord} c2 - coordinate of celestial body 2
|
225 | * @return {Number} angular separation between two celestial bodies
|
226 | */
|
227 | function sepPauwels(c1, c2) {
|
228 | var _base$sincos7 = _base2.default.sincos(c1.dec),
|
229 | _base$sincos8 = _slicedToArray(_base$sincos7, 2),
|
230 | sind1 = _base$sincos8[0],
|
231 | cosd1 = _base$sincos8[1];
|
232 |
|
233 | var _base$sincos9 = _base2.default.sincos(c2.dec),
|
234 | _base$sincos10 = _slicedToArray(_base$sincos9, 2),
|
235 | sind2 = _base$sincos10[0],
|
236 | cosd2 = _base$sincos10[1];
|
237 |
|
238 | var cosdr = cos(c2.ra - c1.ra);
|
239 | var x = cosd1 * sind2 - sind1 * cosd2 * cosdr;
|
240 | var y = cosd2 * sin(c2.ra - c1.ra);
|
241 | var z = sind1 * sind2 + cosd1 * cosd2 * cosdr;
|
242 | return atan2(hypot(x, y), z);
|
243 | }
|
244 |
|
245 | /**
|
246 | * Same as `minSep` but uses function `sepPauwels` to return the minimum
|
247 | * separation between two moving objects.
|
248 | */
|
249 | function minSepPauwels(jd1, jd3, cs1, cs2) {
|
250 | return minSep(jd1, jd3, cs1, cs2, sepPauwels);
|
251 | }
|
252 |
|
253 | /**
|
254 | * RelativePosition returns the position angle of one body with respect to
|
255 | * another.
|
256 | *
|
257 | * The position angle result `p` is measured counter-clockwise from North.
|
258 | * If negative then `p` is in the range of 90° ... 270°
|
259 | *
|
260 | * ````
|
261 | * North
|
262 | * |
|
263 | * (p) ..|
|
264 | * . |
|
265 | * V |
|
266 | * c1 x------------x c2
|
267 | * |
|
268 | * ````
|
269 | *
|
270 | * @param {base.Coord} c1 - coordinate of celestial body 1
|
271 | * @param {base.Coord} c2 - coordinate of celestial body 2
|
272 | * @return {Number} position angle (p)
|
273 | */
|
274 | function relativePosition(c1, c2) {
|
275 | var _base$sincos11 = _base2.default.sincos(c2.ra - c1.ra),
|
276 | _base$sincos12 = _slicedToArray(_base$sincos11, 2),
|
277 | sinΔr = _base$sincos12[0],
|
278 | cosΔr = _base$sincos12[1];
|
279 |
|
280 | var _base$sincos13 = _base2.default.sincos(c2.dec),
|
281 | _base$sincos14 = _slicedToArray(_base$sincos13, 2),
|
282 | sind2 = _base$sincos14[0],
|
283 | cosd2 = _base$sincos14[1];
|
284 |
|
285 | var p = atan2(sinΔr, cosd2 * tan(c1.dec) - sind2 * cosΔr);
|
286 | return p;
|
287 | }
|
288 |
|
289 | exports.default = {
|
290 | sep: sep,
|
291 | minSep: minSep,
|
292 | minSepRect: minSepRect,
|
293 | hav: hav,
|
294 | sepHav: sepHav,
|
295 | minSepHav: minSepHav,
|
296 | sepPauwels: sepPauwels,
|
297 | minSepPauwels: minSepPauwels,
|
298 | relativePosition: relativePosition
|
299 | }; |
\ | No newline at end of file |