1 /*
2 # This file is part of the Astrometry.net suite.
3 # Licensed under a 3-clause BSD style license - see LICENSE
4 */
5
6 #include <math.h>
7 #include <stdio.h>
8 #include <stdlib.h>
9 #include <assert.h>
10 #include <sys/types.h>
11 #include <regex.h>
12
13 #include "os-features.h"
14 #include "keywords.h"
15 #include "mathutil.h"
16 #include "starutil.h"
17 #include "errors.h"
18
19 #define POGSON 2.51188643150958
20 #define LOGP 0.92103403719762
21
22 #define InlineDefine InlineDefineC
23 #include "starutil.inc"
24 #undef InlineDefine
25
radec_derivatives(double ra,double dec,double * dra,double * ddec)26 void radec_derivatives(double ra, double dec, double* dra, double* ddec) {
27 double cosd = cos(deg2rad(dec));
28 double cosra = cos(deg2rad(ra));
29 double sinra = sin(deg2rad(ra));
30 if (dra) {
31 dra[0] = cosd * -sinra;
32 dra[1] = cosd * cosra;
33 dra[2] = 0.0;
34 normalize_3(dra);
35 }
36 if (ddec) {
37 double sind = sin(deg2rad(dec));
38 ddec[0] = -sind * cosra;
39 ddec[1] = -sind * sinra;
40 ddec[2] = cosd;
41 normalize_3(ddec);
42 }
43 }
44
radecrange2xyzrange(double ralo,double declo,double rahi,double dechi,double * minxyz,double * maxxyz)45 void radecrange2xyzrange(double ralo, double declo, double rahi, double dechi,
46 double* minxyz, double* maxxyz) {
47 double minmult, maxmult;
48 double uxlo, uxhi, uylo, uyhi;
49 // Dec only affects z, and is monotonic (z = sin(dec))
50 minxyz[2] = radec2z(0, declo);
51 maxxyz[2] = radec2z(0, dechi);
52
53 // min,max of cos(dec). cos(dec) is concave down.
54 minmult = MIN(cos(deg2rad(declo)), cos(deg2rad(dechi)));
55 maxmult = MAX(cos(deg2rad(declo)), cos(deg2rad(dechi)));
56 if (declo < 0 && dechi > 0)
57 maxmult = 1.0;
58 // unscaled x (ie, cos(ra))
59 uxlo = MIN(cos(deg2rad(ralo)), cos(deg2rad(rahi)));
60 if (ralo < 180 && rahi > 180)
61 uxlo = -1.0;
62 uxhi = MAX(cos(deg2rad(ralo)), cos(deg2rad(rahi)));
63 minxyz[0] = MIN(uxlo * minmult, uxlo * maxmult);
64 maxxyz[0] = MAX(uxhi * minmult, uxhi * maxmult);
65 // unscaled y (ie, sin(ra))
66 uylo = MIN(sin(deg2rad(ralo)), sin(deg2rad(rahi)));
67 if (ralo < 270 && rahi > 270)
68 uylo = -1.0;
69 uyhi = MAX(sin(deg2rad(ralo)), sin(deg2rad(rahi)));
70 if (ralo < 90 && rahi > 90)
71 uyhi = -1.0;
72 minxyz[1] = MIN(uylo * minmult, uylo * maxmult);
73 maxxyz[1] = MAX(uyhi * minmult, uyhi * maxmult);
74 }
75
76
parse_hms_string(const char * str,int * sign,int * term1,int * term2,double * term3)77 static int parse_hms_string(const char* str,
78 int* sign, int* term1, int* term2, double* term3) {
79 anbool matched;
80 regmatch_t matches[6];
81 int nmatches = 6;
82 regex_t re;
83 regmatch_t* m;
84 const char* s;
85
86 const char* restr =
87 "^([+-])?([[:digit:]]{1,2}):"
88 "([[:digit:]]{1,2}):"
89 "([[:digit:]]*(\\.[[:digit:]]*)?)$";
90
91 if (!str)
92 return 1;
93
94 if (regcomp(&re, restr, REG_EXTENDED)) {
95 ERROR("Failed to compile H:M:S regex \"%s\"", restr);
96 return -1;
97 }
98 matched = (regexec(&re, str, nmatches, matches, 0) == 0);
99 regfree(&re);
100 if (!matched)
101 return 1;
102
103 // sign
104 m = matches + 1;
105 s = str + m->rm_so;
106 if (m->rm_so == -1 || s[0] == '+')
107 *sign = 1;
108 else
109 *sign = -1;
110 // hrs / deg
111 m = matches + 2;
112 s = str + m->rm_so;
113 if (s[0] == '0')
114 s++;
115 *term1 = atoi(s);
116 // Min
117 m = matches + 3;
118 s = str + m->rm_so;
119 if (s[0] == '0')
120 s++;
121 *term2 = atoi(s);
122 // Sec
123 m = matches + 4;
124 s = str + m->rm_so;
125 *term3 = atof(s);
126 return 0;
127 }
128
atora(const char * str)129 double atora(const char* str) {
130 char* eptr;
131 double ra;
132 int sgn, hr, min;
133 double sec;
134 int rtn;
135
136 if (!str) {
137 //ERROR("Null string to atora()");
138 return HUGE_VAL;
139 }
140 rtn = parse_hms_string(str, &sgn, &hr, &min, &sec);
141 if (rtn == -1) {
142 ERROR("Failed to run regex");
143 return HUGE_VAL;
144 }
145 if (rtn == 0)
146 return sgn * hms2ra(hr, min, sec);
147
148 ra = strtod(str, &eptr);
149 if (eptr == str)
150 // no conversion
151 return HUGE_VAL;
152 return ra;
153 }
154
atodec(const char * str)155 double atodec(const char* str) {
156 char* eptr;
157 double dec;
158 int sgn, deg, min;
159 double sec;
160 int rtn;
161
162 rtn = parse_hms_string(str, &sgn, °, &min, &sec);
163 if (rtn == -1) {
164 ERROR("Failed to run regex");
165 return HUGE_VAL;
166 }
167 if (rtn == 0)
168 return dms2dec(sgn, deg, min, sec);
169
170 dec = strtod(str, &eptr);
171 if (eptr == str)
172 // no conversion
173 return HUGE_VAL;
174 return dec;
175 }
176
mag2flux(double mag)177 double mag2flux(double mag) {
178 return pow(POGSON, -mag);
179 }
180
flux2mag(double flux)181 double flux2mag(double flux) {
182 return -log(flux) * LOGP;
183 }
184
xyzarr2radecarr(const double * xyz,double * radec)185 inline void xyzarr2radecarr(const double* xyz, double *radec) {
186 xyz2radec(xyz[0], xyz[1], xyz[2], radec+0, radec+1);
187 }
188
distsq_between_radecdeg(double ra1,double dec1,double ra2,double dec2)189 double distsq_between_radecdeg(double ra1, double dec1,
190 double ra2, double dec2) {
191 double xyz1[3];
192 double xyz2[3];
193 radecdeg2xyzarr(ra1, dec1, xyz1);
194 radecdeg2xyzarr(ra2, dec2, xyz2);
195 return distsq(xyz1, xyz2, 3);
196 }
197
arcsec_between_radecdeg(double ra1,double dec1,double ra2,double dec2)198 double arcsec_between_radecdeg(double ra1, double dec1,
199 double ra2, double dec2) {
200 return distsq2arcsec(distsq_between_radecdeg(ra1, dec1, ra2, dec2));
201 }
202
deg_between_radecdeg(double ra1,double dec1,double ra2,double dec2)203 double deg_between_radecdeg(double ra1, double dec1, double ra2, double dec2) {
204 return arcsec2deg(arcsec_between_radecdeg(ra1, dec1, ra2, dec2));
205 }
206
project_equal_area(double x,double y,double z,double * projx,double * projy)207 void project_equal_area(double x, double y, double z, double* projx, double* projy) {
208 double Xp = x*sqrt(1./(1. + z));
209 double Yp = y*sqrt(1./(1. + z));
210 Xp = 0.5 * (1.0 + Xp);
211 Yp = 0.5 * (1.0 + Yp);
212 assert(Xp >= 0.0 && Xp <= 1.0);
213 assert(Yp >= 0.0 && Yp <= 1.0);
214 *projx = Xp;
215 *projy = Yp;
216 }
217
project_hammer_aitoff_x(double x,double y,double z,double * projx,double * projy)218 void project_hammer_aitoff_x(double x, double y, double z, double* projx, double* projy) {
219 double theta = atan(x/z);
220 double r = sqrt(x*x+z*z);
221 double zp, xp;
222 /* Hammer-Aitoff projection with x-z plane compressed to purely +z side
223 * of xz plane */
224 if (z < 0) {
225 if (x < 0) {
226 theta = theta - M_PI;
227 } else {
228 theta = M_PI + theta;
229 }
230 }
231 theta /= 2.0;
232 zp = r*cos(theta);
233 xp = r*sin(theta);
234 assert(zp >= -0.01);
235 project_equal_area(xp, y, zp, projx, projy);
236 }
237
238 /* makes a star object located uniformly at random within the limits given
239 on the sphere */
make_rand_star(double * star,double ramin,double ramax,double decmin,double decmax)240 void make_rand_star(double* star, double ramin, double ramax,
241 double decmin, double decmax)
242 {
243 double decval, raval;
244 if (ramin < 0.0)
245 ramin = 0.0;
246 if (ramax > (2*M_PI))
247 ramax = 2 * M_PI;
248 if (decmin < -M_PI / 2.0)
249 decmin = -M_PI / 2.0;
250 if (decmax > M_PI / 2.0)
251 decmax = M_PI / 2.0;
252
253 decval = asin(uniform_sample(sin(decmin), sin(decmax)));
254 raval = uniform_sample(ramin, ramax);
255 star[0] = radec2x(raval, decval);
256 star[1] = radec2y(raval, decval);
257 star[2] = radec2z(raval, decval);
258 }
259
260
261 // RA in degrees to Mercator X coordinate [0, 1).
ra2mercx(double ra)262 inline double ra2mercx(double ra) {
263 double mx = ra / 360.0;
264 if (mx < 0.0 || mx > 1.0) {
265 mx = fmod(mx, 1.0);
266 if (mx < 0.0)
267 mx += 1.0;
268 }
269 return mx;
270 }
271
272 // Dec in degrees to Mercator Y coordinate [0, 1).
dec2mercy(double dec)273 inline double dec2mercy(double dec) {
274 return 0.5 + (asinh(tan(deg2rad(dec))) / (2.0 * M_PI));
275 }
276
hms2ra(int h,int m,double s)277 double hms2ra(int h, int m, double s) {
278 return 15.0 * (h + ((m + (s / 60.0)) / 60.0));
279 }
280
dms2dec(int sgn,int d,int m,double s)281 double dms2dec(int sgn, int d, int m, double s) {
282 return sgn * (d + ((m + (s / 60.0)) / 60.0));
283 }
284
285 // RA in degrees to H:M:S
ra2hms(double ra,int * h,int * m,double * s)286 inline void ra2hms(double ra, int* h, int* m, double* s) {
287 double rem;
288 ra = fmod(ra, 360.0);
289 if (ra < 0.0)
290 ra += 360.0;
291 rem = ra / 15.0;
292 *h = floor(rem);
293 // remaining (fractional) hours
294 rem -= *h;
295 // -> minutes
296 rem *= 60.0;
297 *m = floor(rem);
298 // remaining (fractional) minutes
299 rem -= *m;
300 // -> seconds
301 rem *= 60.0;
302 *s = rem;
303 }
304
305 // Dec in degrees to D:M:S
dec2dms(double dec,int * sign,int * d,int * m,double * s)306 void dec2dms(double dec, int* sign, int* d, int* m, double* s) {
307 double rem;
308 double flr;
309 *sign = (dec >= 0.0) ? 1 : -1;
310 dec *= (*sign);
311 flr = floor(dec);
312 *d = flr;
313 // remaining degrees:
314 rem = dec - flr;
315 // -> minutes
316 rem *= 60.0;
317 *m = floor(rem);
318 // remaining (fractional) minutes
319 rem -= *m;
320 // -> seconds
321 rem *= 60.0;
322 *s = rem;
323 }
324
ra2hmsstring(double ra,char * str)325 void ra2hmsstring(double ra, char* str) {
326 int h, m;
327 double s;
328 int ss;
329 int ds;
330 ra2hms(ra, &h, &m, &s);
331
332 // round to display to 3 decimal places
333 ss = (int)floor(s);
334 ds = (int)round((s - ss) * 1000.0);
335 if (ds >= 1000) {
336 ss++;
337 ds -= 1000;
338 }
339 if (ss >= 60) {
340 ss -= 60;
341 m += 1;
342 }
343 if (m >= 60) {
344 m -= 60;
345 h += 1;
346 }
347 sprintf(str, "%02i:%02i:%02i.%03i", h, m, ss, ds);
348 }
349
dec2dmsstring(double dec,char * str)350 void dec2dmsstring(double dec, char* str) {
351 int sign, d, m;
352 double s;
353 int ss, ds;
354 dec2dms(dec, &sign, &d, &m, &s);
355 ss = (int)floor(s);
356 ds = (int)round((s - ss) * 1000.0);
357 if (ds >= 1000) {
358 ss++;
359 ds -= 1000;
360 }
361 if (ss >= 60) {
362 ss -= 60;
363 m += 1;
364 }
365 if (m >= 60) {
366 m -= 60;
367 d += 1;
368 }
369 sprintf(str, "%c%02i:%02i:%02i.%03i", (sign==1 ? '+':'-'), d, m, ss, ds);
370 }
371
372