1 /*
2 astrom.c
3
4 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 *
6 * Part of: SExtractor
7 *
8 * Author: E.BERTIN (IAP)
9 *
10 * Contents: Astrometrical computations.
11 *
12 * Last modify: 13/07/2006
13 *
14 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15 */
16
17 #ifdef HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20
21 #include <math.h>
22 #include <stdlib.h>
23 #include <string.h>
24
25 #include "wcs/wcs.h"
26 #include "define.h"
27 #include "globals.h"
28 #include "prefs.h"
29 #include "astrom.h"
30 #include "wcs/tnx.h"
31
32 static obj2struct *obj2 = &outobj2;
33
34 /****************************** initastrom **********************************/
35 /*
36 Initialize astrometrical structures.
37 */
initastrom(picstruct * field)38 void initastrom(picstruct *field)
39
40 {
41 astromstruct *as;
42 double *lm;
43 int l,n, lng,lat, naxis;
44
45 as = field->astrom;
46
47 naxis = as->naxis;
48
49 /* Test if the WCS is in use */
50 if (as->wcs_flag)
51 /*-- ...Yes: launch the WCS stuff! */
52 {
53 QCALLOC(as->lin, struct linprm, 1);
54 QMALLOC(as->cel, struct celprm, 1);
55 QMALLOC(as->prj, struct prjprm, 1);
56 QMALLOC(as->lin->cdelt, double, naxis);
57 QMALLOC(as->lin->crpix, double, naxis);
58 QMALLOC(as->lin->pc, double, naxis*naxis);
59 /* Set WCS flags to 0: structures will be reinitialized by the WCS library */
60 as->wcs->flag = as->lin->flag = as->cel->flag = as->prj->flag = 0;
61 as->lin->naxis = as->naxis;
62
63 /* wcsprm structure */
64 lng = as->lng = as->wcs->lng;
65 lat = as->lat = as->wcs->lat;
66
67 /*-- linprm structure */
68 for (l=0; l<naxis; l++)
69 {
70 as->lin->crpix[l] = as->crpix[l];
71 as->lin->cdelt[l] = as->cdelt[l];
72 as->cel->ref[l] = as->crval[l];
73 }
74 for (l=0; l<naxis*naxis; l++)
75 as->lin->pc[l] = as->pc[l];
76
77 /*-- celprm structure */
78 if (lng>=0)
79 {
80 as->cel->ref[0] = as->crval[lng];
81 as->cel->ref[1] = as->crval[lat];
82 }
83 else
84 {
85 as->cel->ref[0] = as->crval[0];
86 as->cel->ref[1] = as->crval[1];
87 }
88 as->cel->ref[2] = as->longpole;
89 as->cel->ref[3] = as->latpole;
90
91 /* prjprm structure */
92 as->prj->r0 = as->r0;
93 as->prj->tnx_lngcor = as->tnx_lngcor;
94 as->prj->tnx_latcor = as->tnx_latcor;
95 if (lng>=0)
96 {
97 n = 0;
98 for (l=100; l--;)
99 {
100 as->prj->p[l] = as->projp[l+lng*100];
101 as->prj->p[l+100] = as->projp[l+lat*100];
102 if (!n && (as->prj->p[l] || as->prj->p[l+100]))
103 n = l+1;
104 }
105 }
106
107 /*-- Compute an "average linear matrix" (at field center) */
108 compute_wcs(field, (field->width+1)/2.0, (field->height+1)/2.0);
109
110 /*---- Compute Pole coordinates in J2000 and/or B1950 for THETAs */
111 if (FLAG(obj2.theta2000) || FLAG(obj2.theta1950)
112 || FLAG(obj2.poserr_theta2000) || FLAG(obj2.poserr_theta1950)
113 || FLAG(obj2.win_theta2000) || FLAG(obj2.win_theta1950)
114 || FLAG(obj2.winposerr_theta2000) || FLAG(obj2.winposerr_theta1950))
115 {
116 if (fabs(as->equinox-2000.0)>0.003)
117 precess(as->equinox, 0.0, 90.0, 2000.0, &as->ap2000, &as->dp2000);
118 else
119 {
120 as->ap2000 = 0.0;
121 as->dp2000 = 90.0;
122 }
123
124 if (FLAG(obj2.theta1950) || FLAG(obj2.poserr_theta1950))
125 j2b(as->equinox, as->ap2000, as->dp2000, &as->ap1950, &as->dp1950);
126 }
127 }
128 else
129 /*-- ...No: compute only the determinant */
130 {
131 /*-- Simplify the original FITS PC matrix */
132 lm = as->linmat;
133 for (l=0; l<naxis*naxis; l++)
134 lm[l] = as->pc[l]*as->cdelt[l/naxis];
135 /*-- Check valid only in 2D */
136 if ((as->lindet = lm[0]*lm[3] - lm[1]*lm[2]) == 0.0)
137 warning ("Null determinant in the global distortion matrix:\n",
138 " Some WORLD-parameters will be incorrect");
139 }
140
141 /* Override astrometric definitions only if user supplies a pixel-scale */
142 if (prefs.pixel_scale == 0.0)
143 {
144 as->pixscale = sqrt(fabs(as->lindet));
145 field->pixscale = 3600.0*as->pixscale; /* in arcsec2 */
146 }
147 else
148 as->pixscale = (field->pixscale=prefs.pixel_scale)/3600.0;
149
150 return;
151 }
152
153
154 /**************************** computeastrom *********************************/
155 /*
156 Compute real WORLD coordinates and dimensions according to FITS info.
157 */
computeastrom(picstruct * field,objstruct * obj)158 void computeastrom(picstruct *field, objstruct *obj)
159
160 {
161 astromstruct *as;
162 double *lm, *wcspos;
163
164 as = field->astrom;
165 lm = as->linmat;
166
167 /* If working with WCS, compute WORLD coordinates and local matrix */
168 if (FLAG(obj2.mxw))
169 {
170 if (as->wcs_flag)
171 {
172 wcspos = compute_wcs(field, obj2->posx, obj2->posy);
173 obj2->alphas = obj2->mxw = wcspos[0];
174 obj2->deltas = obj2->myw = wcspos[1];
175 if (FLAG(obj2.alpha2000))
176 {
177 if (fabs(as->equinox-2000.0)>0.003)
178 precess(as->equinox, wcspos[0], wcspos[1],
179 2000.0, &obj2->alpha2000, &obj2->delta2000);
180 else
181 {
182 obj2->alpha2000 = obj2->mxw;
183 obj2->delta2000 = obj2->myw;
184 }
185 if (FLAG(obj2.alpha1950))
186 j2b(as->equinox, obj2->alpha2000, obj2->delta2000,
187 &obj2->alpha1950, &obj2->delta1950);
188 }
189 }
190 else
191 {
192 double dx,dy;
193
194 dx = obj2->posx - as->crpix[0];
195 dy = obj2->posy - as->crpix[1];
196 obj2->mxw = as->crval[0]+ lm[0]*dx + lm[1]*dy; /* CDELT included! */
197 obj2->myw = as->crval[1]+ lm[2]*dx + lm[3]*dy; /* CDELT included! */
198 }
199 }
200
201 /* Idem for peak-flux positions */
202 if (FLAG(obj2.peakxw))
203 {
204 if (as->wcs_flag)
205 {
206 wcspos = compute_wcs(field, (double)obj->peakx, (double)obj->peaky);
207 obj2->peakalphas = obj2->peakxw = wcspos[0];
208 obj2->peakdeltas = obj2->peakyw = wcspos[1];
209 if (FLAG(obj2.peakalpha2000))
210 {
211 if (fabs(as->equinox-2000.0)>0.003)
212 precess(as->equinox, wcspos[0], wcspos[1],
213 2000.0, &obj2->peakalpha2000, &obj2->peakdelta2000);
214 else
215 {
216 obj2->peakalpha2000 = obj2->peakxw;
217 obj2->peakdelta2000 = obj2->peakyw;
218 }
219 if (FLAG(obj2.peakalpha1950))
220 j2b(as->equinox, obj2->peakalpha2000, obj2->peakdelta2000,
221 &obj2->peakalpha1950, &obj2->peakdelta1950);
222 }
223 }
224 else
225 {
226 double dx,dy;
227
228 dx = obj->peakx - as->crpix[0];
229 dy = obj->peaky - as->crpix[1];
230 obj2->peakxw = as->crval[0]+ lm[0]*dx + lm[1]*dy; /* CDELT included! */
231 obj2->peakyw = as->crval[1]+ lm[2]*dx + lm[3]*dy; /* CDELT included! */
232 }
233 }
234
235 /* Idem for Windowed positions */
236 if (FLAG(obj2.winpos_xw))
237 {
238 if (as->wcs_flag)
239 {
240 wcspos = compute_wcs(field, obj2->winpos_x, obj2->winpos_y);
241 obj2->winpos_alphas = obj2->winpos_xw = wcspos[0];
242 obj2->winpos_deltas = obj2->winpos_yw = wcspos[1];
243 if (FLAG(obj2.winpos_alpha2000))
244 {
245 if (fabs(as->equinox-2000.0)>0.003)
246 precess(as->equinox, wcspos[0], wcspos[1],
247 2000.0, &obj2->winpos_alpha2000, &obj2->winpos_delta2000);
248 else
249 {
250 obj2->winpos_alpha2000 = obj2->winpos_xw;
251 obj2->winpos_delta2000 = obj2->winpos_yw;
252 }
253 if (FLAG(obj2.winpos_alpha1950))
254 j2b(as->equinox, obj2->winpos_alpha2000, obj2->winpos_delta2000,
255 &obj2->winpos_alpha1950, &obj2->winpos_delta1950);
256 }
257 }
258 else
259 {
260 double dx,dy;
261
262 dx = obj2->winpos_x - as->crpix[0];
263 dy = obj2->winpos_y - as->crpix[1];
264 obj2->winpos_xw = as->crval[0]+ lm[0]*dx + lm[1]*dy;/* CDELT included! */
265 obj2->winpos_yw = as->crval[1]+ lm[2]*dx + lm[3]*dy;/* CDELT included! */
266 }
267 }
268
269 /* Custom coordinate system for the MAMA machine */
270 if (FLAG(obj2.mamaposx))
271 {
272 double dx,dy;
273
274 dx = obj2->posx - 0.5;
275 dy = obj2->posy - 0.5;
276 obj2->mamaposx = (as->crval[1]+lm[2]*dx+lm[3]*dy)
277 *(MAMA_CORFLEX+1.0); /* CDELT included! */
278 obj2->mamaposy = (as->crval[0]+lm[0]*dx+lm[1]*dy); /* CDELT included! */
279 }
280
281 /* Express shape parameters in WORLD frame */
282 if (FLAG(obj2.mx2w))
283 astrom_shapeparam(field, obj);
284 if (FLAG(obj2.win_mx2w))
285 astrom_winshapeparam(field, obj);
286
287 /* Express position error parameters in WORLD frame */
288 if (FLAG(obj2.poserr_mx2w))
289 astrom_errparam(field, obj);
290 if (FLAG(obj2.winposerr_mx2w))
291 astrom_winerrparam(field, obj);
292
293 if (FLAG(obj2.npixw))
294 obj2->npixw = obj->npix*as->pixscale*as->pixscale;
295 if (FLAG(obj2.fdnpixw))
296 obj2->fdnpixw = obj->fdnpix*as->pixscale*as->pixscale;
297
298 if (FLAG(obj2.fwhmw))
299 obj2->fwhmw = obj->fwhm*as->pixscale;
300
301 return;
302 }
303
304
305 /****************************** compute_wcs *********************************/
306 /*
307 Compute real WORLD coordinates and local distortion matrix according to the
308 WCS info.
309 */
compute_wcs(picstruct * field,double mx,double my)310 double *compute_wcs(picstruct *field, double mx, double my)
311
312 {
313 astromstruct *as;
314 static double pixpos[NAXIS], wcspos[NAXIS],wcspos0[2], imgcrd[NAXIS],
315 phi,theta;
316 double *lm, al,da,de,cde;
317 int rcode, lng,lat, naxis;
318
319 as = field->astrom;
320 lm = as->linmat;
321
322 naxis = as->naxis;
323 lng = as->lng;
324 lat = as->lat;
325 if (lng == lat)
326 {
327 lng = 0;
328 lat = 1;
329 }
330
331 pixpos[lng] = mx;
332 pixpos[lat] = my;
333
334 if ((rcode=wcsrev((const char(*)[9])as->ctype, as->wcs, pixpos, as->lin,
335 imgcrd, as->prj, &phi, &theta, as->crval, as->cel, wcspos)))
336 error(EXIT_FAILURE, "*Error* in WCSlib: ", (char *)wcsrev_errmsg[rcode]);
337
338 /* Compute the local distortion matrix */
339 al = wcspos0[lng<lat?0:1] = wcspos[lng];
340 de = wcspos0[lng<lat?1:0] = wcspos[lat];
341
342 /* Get world coordinates for vector 1,0 */
343 pixpos[lng] = mx + 1.0;
344 pixpos[lat] = my;
345 if ((rcode=wcsrev((const char(*)[9])as->ctype, as->wcs, pixpos, as->lin,
346 imgcrd, as->prj, &phi, &theta, as->crval, as->cel, wcspos)))
347 error(EXIT_FAILURE, "*Error* in WCSlib: ", (char *)wcsrev_errmsg[rcode]);
348
349 da = wcspos[lng]-al;
350 if (da>180.0)
351 da -= 360.0;
352 else if (da<-180.0)
353 da += 360.0;
354
355 lm[lng] = da*(cde=cos(de*DEG));
356 lm[lat] = wcspos[lat] - de;
357
358 /* Get world coordinates for vector 0,1 */
359 /* Second one */
360 pixpos[lng] = mx;
361 pixpos[lat] = my + 1.0;
362 if ((rcode=wcsrev((const char(*)[9])as->ctype, as->wcs, pixpos, as->lin,
363 imgcrd, as->prj, &phi, &theta, as->crval, as->cel, wcspos)))
364 error(EXIT_FAILURE, "*Error* in WCSlib: ", (char *)wcsrev_errmsg[rcode]);
365
366 da = wcspos[lng]-al;
367 if (da>180.0)
368 da -= 360.0;
369 else if (da<-180.0)
370 da += 360.0;
371
372 lm[2] = da*cde;
373 lm[3] = wcspos[lat] - de;
374
375 as->lindet = lm[lng+lng*naxis]*lm[lat+lat*naxis]
376 - lm[lng+lat*naxis]*lm[lat+lng*naxis];
377 if (as->lindet == 0.0)
378 warning ("Null determinant in the local distortion matrix:\n",
379 " Some WORLD-parameters will be incorrect");
380
381 if (prefs.pixel_scale == 0.0)
382 as->pixscale = sqrt(fabs(as->lindet));
383
384 return wcspos0;
385 }
386
387
388 /****************************** astrom_shapeparam ****************************/
389 /*
390 Compute shape parameters in WORLD and SKY coordinates.
391 */
astrom_shapeparam(picstruct * field,objstruct * obj)392 void astrom_shapeparam(picstruct *field, objstruct *obj)
393 {
394 astromstruct *as;
395 double *lm,
396 dx2,dy2,dxy, xm2,ym2,xym, temp,pm2, lm0,lm1,lm2,lm3;
397 int lng,lat, naxis;
398
399 as = field->astrom;
400 lm = as->linmat;
401
402 naxis = as->naxis;
403 lng = as->lng;
404 lat = as->lat;
405 if (lng == lat)
406 {
407 lng = 0;
408 lat = 1;
409 }
410 lm0 = lm[lng+naxis*lng];
411 lm1 = lm[lat+naxis*lng];
412 lm2 = lm[lng+naxis*lat];
413 lm3 = lm[lat+naxis*lat];
414
415
416 /* All WORLD params based on 2nd order moments have to pass through here */
417 dx2 = obj->mx2;
418 dy2 = obj->my2;
419 dxy = obj->mxy;
420 obj2->mx2w = xm2 = lm0*lm0*dx2 + lm1*lm1*dy2 + lm0*lm1*dxy;
421 obj2->my2w = ym2 = lm2*lm2*dx2 + lm3*lm3*dy2 + lm2*lm3*dxy;
422 obj2->mxyw = xym = lm0*lm2*dx2 + lm1*lm3*dy2 + (lm0*lm3+lm1*lm2)*dxy;
423 temp=xm2-ym2;
424 if (FLAG(obj2.thetaw))
425 {
426 obj2->thetaw = (temp == 0.0)? (45.0) : (0.5*atan2(2.0 * xym,temp)/DEG);
427 if (as->wcs_flag && FLAG(obj2.thetas))
428 obj2->thetas = obj2->thetaw + (obj2->thetaw>0.0?-90:90.0);
429
430 /*-- Compute position angles in J2000 or B1950 reference frame */
431 if (as->wcs_flag)
432 {
433 double da,dd;
434
435 if (FLAG(obj2.theta2000))
436 {
437 da = as->ap2000 - obj2->alpha2000;
438 dd = (sin(as->dp2000*DEG)
439 -sin(obj2->delta2000*DEG)*sin(obj2->deltas*DEG))
440 /(cos(obj2->delta2000*DEG)*cos(obj2->deltas*DEG));
441 dd = dd<1.0? (dd>-1.0?acos(dd)/DEG:180.0) : 0.0;
442 obj2->theta2000 = obj2->thetas
443 + (((da>0.0 && da<180.0) || da<-180.0)?-dd:dd);
444 }
445
446 if (FLAG(obj2.theta1950))
447 {
448 da = as->ap1950 - obj2->alpha1950;
449 dd = (sin(as->dp1950*DEG)
450 -sin(obj2->delta1950*DEG)*sin(obj2->deltas*DEG))
451 /(cos(obj2->delta1950*DEG)*cos(obj2->deltas*DEG));
452 dd = dd<1.0? (dd>-1.0?acos(dd)/DEG:180.0) : 0.0;
453 obj2->theta1950 = obj2->thetas
454 + (((da>0.0 && da<180.0) || da<-180.0)?-dd:dd);
455 }
456 }
457 }
458
459 if (FLAG(obj2.aw))
460 {
461 temp = sqrt(0.25*temp*temp+xym*xym);
462 pm2 = 0.5*(xm2+ym2);
463 obj2->aw = (float)sqrt(pm2+temp);
464 obj2->bw = (float)sqrt(pm2-temp);
465 obj2->polarw = temp / pm2;
466 }
467
468 if (FLAG(obj2.cxxw))
469 {
470 /*-- Handle large, fully correlated profiles (can cause a singularity...) */
471 if ((temp=xm2*ym2-xym*xym)<1e-6)
472 {
473 temp = 1e-6;
474 xym *= 0.99999;
475 }
476 obj2->cxxw = (float)(ym2/temp);
477 obj2->cyyw = (float)(xm2/temp);
478 obj2->cxyw = (float)(-2*xym/temp);
479 }
480
481 return;
482 }
483
484
485 /**************************** astrom_winshapeparam ***************************/
486 /*
487 Compute shape parameters in WORLD and SKY coordinates.
488 */
astrom_winshapeparam(picstruct * field,objstruct * obj)489 void astrom_winshapeparam(picstruct *field, objstruct *obj)
490 {
491 astromstruct *as;
492 double *lm,
493 dx2,dy2,dxy, xm2,ym2,xym, temp,pm2, lm0,lm1,lm2,lm3;
494 int lng,lat, naxis;
495
496 as = field->astrom;
497 lm = as->linmat;
498
499 naxis = as->naxis;
500 lng = as->lng;
501 lat = as->lat;
502 if (lng == lat)
503 {
504 lng = 0;
505 lat = 1;
506 }
507 lm0 = lm[lng+naxis*lng];
508 lm1 = lm[lat+naxis*lng];
509 lm2 = lm[lng+naxis*lat];
510 lm3 = lm[lat+naxis*lat];
511
512 /* All WORLD params based on 2nd order moments have to pass through here */
513 dx2 = obj2->win_mx2;
514 dy2 = obj2->win_my2;
515 dxy = obj2->win_mxy;
516 obj2->win_mx2w = xm2 = lm0*lm0*dx2 + lm1*lm1*dy2 + lm0*lm1*dxy;
517 obj2->win_my2w = ym2 = lm2*lm2*dx2 + lm3*lm3*dy2 + lm2*lm3*dxy;
518 obj2->win_mxyw = xym = lm0*lm2*dx2 + lm1*lm3*dy2 + (lm0*lm3+lm1*lm2)*dxy;
519 temp=xm2-ym2;
520 if (FLAG(obj2.win_thetaw))
521 {
522 obj2->win_thetaw = (temp == 0.0)? (45.0) : (0.5*atan2(2.0*xym,temp)/DEG);
523 if (as->wcs_flag && FLAG(obj2.win_thetas))
524 obj2->win_thetas = obj2->win_thetaw +
525 (obj2->win_thetaw>0.0?-90:90.0);
526
527 /*-- Compute position angles in J2000 or B1950 reference frame */
528 if (as->wcs_flag)
529 {
530 double da,dd;
531
532 if (FLAG(obj2.win_theta2000))
533 {
534 da = as->ap2000 - obj2->winpos_alpha2000;
535 dd = (sin(as->dp2000*DEG)
536 -sin(obj2->winpos_delta2000*DEG)*sin(obj2->winpos_deltas*DEG))
537 /(cos(obj2->winpos_delta2000*DEG)*cos(obj2->winpos_deltas*DEG));
538 dd = dd<1.0? (dd>-1.0?acos(dd)/DEG:180.0) : 0.0;
539 obj2->win_theta2000 = obj2->win_thetas
540 + (((da>0.0 && da<180.0) || da<-180.0)?-dd:dd);
541 }
542
543 if (FLAG(obj2.win_theta1950))
544 {
545 da = as->ap1950 - obj2->winpos_alpha1950;
546 dd = (sin(as->dp1950*DEG)
547 -sin(obj2->winpos_delta1950*DEG)*sin(obj2->winpos_deltas*DEG))
548 /(cos(obj2->winpos_delta1950*DEG)*cos(obj2->winpos_deltas*DEG));
549 dd = dd<1.0? (dd>-1.0?acos(dd)/DEG:180.0) : 0.0;
550 obj2->win_theta1950 = obj2->win_thetas
551 + (((da>0.0 && da<180.0) || da<-180.0)?-dd:dd);
552 }
553 }
554 }
555
556 if (FLAG(obj2.win_aw))
557 {
558 temp = sqrt(0.25*temp*temp+xym*xym);
559 pm2 = 0.5*(xm2+ym2);
560 obj2->win_aw = (float)sqrt(pm2+temp);
561 obj2->win_bw = (float)sqrt(pm2-temp);
562 obj2->win_polarw = temp / pm2;
563 }
564
565 if (FLAG(obj2.win_cxxw))
566 {
567 /*-- Handle large, fully correlated profiles (can cause a singularity...) */
568 if ((temp=xm2*ym2-xym*xym)<1e-6)
569 {
570 temp = 1e-6;
571 xym *= 0.99999;
572 }
573 obj2->win_cxxw = (float)(ym2/temp);
574 obj2->win_cyyw = (float)(xm2/temp);
575 obj2->win_cxyw = (float)(-2*xym/temp);
576 }
577
578 return;
579 }
580
581
582 /******************************* astrom_errparam *****************************/
583 /*
584 Compute error ellipse parameters in WORLD and SKY coordinates.
585 */
astrom_errparam(picstruct * field,objstruct * obj)586 void astrom_errparam(picstruct *field, objstruct *obj)
587 {
588 astromstruct *as;
589 double *lm,
590 dx2,dy2,dxy, xm2,ym2,xym, temp,pm2, lm0,lm1,lm2,lm3;
591 int lng,lat, naxis;
592
593 as = field->astrom;
594 lm = as->linmat;
595
596 naxis = as->naxis;
597 lng = as->lng;
598 lat = as->lat;
599 if (lng == lat)
600 {
601 lng = 0;
602 lat = 1;
603 }
604 lm0 = lm[lng+naxis*lng];
605 lm1 = lm[lat+naxis*lng];
606 lm2 = lm[lng+naxis*lat];
607 lm3 = lm[lat+naxis*lat];
608
609 /* All WORLD params based on 2nd order moments have to pass through here */
610 dx2 = obj->poserr_mx2;
611 dy2 = obj->poserr_my2;
612 dxy = obj->poserr_mxy;
613 obj2->poserr_mx2w = xm2 = lm0*lm0*dx2 + lm1*lm1*dy2 + lm0*lm1*dxy;
614 obj2->poserr_my2w = ym2 = lm2*lm2*dx2 + lm3*lm3*dy2 + lm2*lm3*dxy;
615 obj2->poserr_mxyw = xym = lm0*lm2*dx2 + lm1*lm3*dy2 + (lm0*lm3+lm1*lm2)*dxy;
616 temp=xm2-ym2;
617 if (FLAG(obj2.poserr_thetaw))
618 {
619 obj2->poserr_thetaw = (temp==0.0)? (45.0):(0.5*atan2(2.0*xym,temp)/DEG);
620 if (as->wcs_flag && FLAG(obj2.poserr_thetas))
621 obj2->poserr_thetas = obj2->poserr_thetaw
622 + (obj2->poserr_thetaw>0.0? -90:90.0);
623
624 /*-- Compute position angles in J2000 or B1950 reference frame */
625 if (as->wcs_flag)
626 {
627 double da,dd;
628
629 if (FLAG(obj2.poserr_theta2000))
630 {
631 da = as->ap2000 - obj2->alpha2000;
632 dd = (sin(as->dp2000*DEG)
633 -sin(obj2->delta2000*DEG)*sin(obj2->deltas*DEG))
634 /(cos(obj2->delta2000*DEG)*cos(obj2->deltas*DEG));
635 dd = dd<1.0? (dd>-1.0?acos(dd)/DEG:180.0) : 0.0;
636 obj2->poserr_theta2000 = obj2->poserr_thetas
637 + (((da>0.0 && da<180.0) || da<-180.0)?-dd:dd);
638 }
639
640 if (FLAG(obj2.poserr_theta1950))
641 {
642 da = as->ap1950 - obj2->alpha1950;
643 dd = (sin(as->dp1950*DEG)
644 -sin(obj2->delta1950*DEG)*sin(obj2->deltas*DEG))
645 /(cos(obj2->delta1950*DEG)*cos(obj2->deltas*DEG));
646 dd = dd<1.0? (dd>-1.0?acos(dd)/DEG:180.0) : 0.0;
647 obj2->poserr_theta1950 = obj2->poserr_thetas
648 + (((da>0.0 && da<180.0) || da<-180.0)?-dd:dd);
649 }
650 }
651 }
652
653 if (FLAG(obj2.poserr_aw))
654 {
655 temp = sqrt(0.25*temp*temp+xym*xym);
656 pm2 = 0.5*(xm2+ym2);
657 obj2->poserr_aw = (float)sqrt(pm2+temp);
658 obj2->poserr_bw = (float)sqrt(pm2-temp);
659 }
660
661 if (FLAG(obj2.poserr_cxxw))
662 {
663 /*-- Handle large, fully correlated profiles (can cause a singularity...) */
664 if ((temp=xm2*ym2-xym*xym)<1e-6)
665 {
666 temp = 1e-6;
667 xym *= 0.99999;
668 }
669 obj2->poserr_cxxw = (float)(ym2/temp);
670 obj2->poserr_cyyw = (float)(xm2/temp);
671 obj2->poserr_cxyw = (float)(-2*xym/temp);
672 }
673
674 return;
675 }
676
677
678 /***************************** astrom_winerrparam ***************************/
679 /*
680 Compute error ellipse parameters in WORLD and SKY coordinates.
681 */
astrom_winerrparam(picstruct * field,objstruct * obj)682 void astrom_winerrparam(picstruct *field, objstruct *obj)
683 {
684 astromstruct *as;
685 double *lm,
686 dx2,dy2,dxy, xm2,ym2,xym, temp,pm2, lm0,lm1,lm2,lm3;
687 int lng,lat, naxis;
688
689 as = field->astrom;
690 lm = as->linmat;
691
692 naxis = as->naxis;
693 lng = as->lng;
694 lat = as->lat;
695 if (lng == lat)
696 {
697 lng = 0;
698 lat = 1;
699 }
700 lm0 = lm[lng+naxis*lng];
701 lm1 = lm[lat+naxis*lng];
702 lm2 = lm[lng+naxis*lat];
703 lm3 = lm[lat+naxis*lat];
704
705 /* All WORLD params based on 2nd order moments have to pass through here */
706 dx2 = obj2->winposerr_mx2;
707 dy2 = obj2->winposerr_my2;
708 dxy = obj2->winposerr_mxy;
709 obj2->winposerr_mx2w = xm2 = lm0*lm0*dx2 + lm1*lm1*dy2 + lm0*lm1*dxy;
710 obj2->winposerr_my2w = ym2 = lm2*lm2*dx2 + lm3*lm3*dy2 + lm2*lm3*dxy;
711 obj2->winposerr_mxyw = xym = lm0*lm2*dx2 + lm1*lm3*dy2 + (lm0*lm3+lm1*lm2)*dxy;
712 temp=xm2-ym2;
713 if (FLAG(obj2.winposerr_thetaw))
714 {
715 obj2->winposerr_thetaw = (temp==0.0)? (45.0):(0.5*atan2(2.0*xym,temp)/DEG);
716 if (as->wcs_flag && FLAG(obj2.winposerr_thetas))
717 obj2->winposerr_thetas = obj2->winposerr_thetaw
718 + (obj2->winposerr_thetaw>0.0? -90:90.0);
719
720 /*-- Compute position angles in J2000 or B1950 reference frame */
721 if (as->wcs_flag)
722 {
723 double da,dd;
724
725 if (FLAG(obj2.winposerr_theta2000))
726 {
727 da = as->ap2000 - obj2->winpos_alpha2000;
728 dd = (sin(as->dp2000*DEG)
729 -sin(obj2->winpos_delta2000*DEG)*sin(obj2->winpos_deltas*DEG))
730 /(cos(obj2->winpos_delta2000*DEG)*cos(obj2->winpos_deltas*DEG));
731 dd = dd<1.0? (dd>-1.0?acos(dd)/DEG:180.0) : 0.0;
732 obj2->winposerr_theta2000 = obj2->winposerr_thetas
733 + (((da>0.0 && da<180.0) || da<-180.0)?-dd:dd);
734 }
735
736 if (FLAG(obj2.winposerr_theta1950))
737 {
738 da = as->ap1950 - obj2->winpos_alpha1950;
739 dd = (sin(as->dp1950*DEG)
740 -sin(obj2->winpos_delta1950*DEG)*sin(obj2->winpos_deltas*DEG))
741 /(cos(obj2->winpos_delta1950*DEG)*cos(obj2->winpos_deltas*DEG));
742 dd = dd<1.0? (dd>-1.0?acos(dd)/DEG:180.0) : 0.0;
743 obj2->winposerr_theta1950 = obj2->winposerr_thetas
744 + (((da>0.0 && da<180.0) || da<-180.0)?-dd:dd);
745 }
746 }
747 }
748
749 if (FLAG(obj2.winposerr_aw))
750 {
751 temp = sqrt(0.25*temp*temp+xym*xym);
752 pm2 = 0.5*(xm2+ym2);
753 obj2->winposerr_aw = (float)sqrt(pm2+temp);
754 obj2->winposerr_bw = (float)sqrt(pm2-temp);
755 }
756
757 if (FLAG(obj2.winposerr_cxxw))
758 {
759 /*-- Handle large, fully correlated profiles (can cause a singularity...) */
760 if ((temp=xm2*ym2-xym*xym)<1e-6)
761 {
762 temp = 1e-6;
763 xym *= 0.99999;
764 }
765 obj2->winposerr_cxxw = (float)(ym2/temp);
766 obj2->winposerr_cyyw = (float)(xm2/temp);
767 obj2->winposerr_cxyw = (float)(-2*xym/temp);
768 }
769
770 return;
771 }
772
773
774 /******************************* copyastrom *********************************/
775 /*
776 Copy astrometrical structures.
777 */
copyastrom(picstruct * infield,picstruct * outfield)778 void copyastrom(picstruct *infield, picstruct *outfield)
779
780 {
781 astromstruct *inas, *outas;
782 int naxis;
783
784 if (infield->astrom)
785 {
786 QMEMCPY(infield->astrom, outfield->astrom, astromstruct, 1);
787 inas = infield->astrom;
788 outas = outfield->astrom;
789 naxis = inas->naxis;
790 if (inas->wcs_flag)
791 {
792 QMEMCPY(inas->wcs, outas->wcs, struct wcsprm, 1);
793 QMEMCPY(inas->lin, outas->lin, struct linprm, 1);
794 QMEMCPY(inas->cel, outas->cel, struct celprm, 1);
795 QMEMCPY(inas->prj, outas->prj, struct prjprm, 1);
796 QMEMCPY(inas->lin->cdelt, outas->lin->cdelt, double, naxis);
797 QMEMCPY(inas->lin->crpix, outas->lin->crpix, double, naxis);
798 QMEMCPY(inas->lin->pc, outas->lin->pc, double, naxis*naxis);
799 outas->tnx_lngcor = copy_tnxaxis(inas->tnx_lngcor);
800 outas->tnx_latcor = copy_tnxaxis(inas->tnx_latcor);
801 }
802 }
803
804 return;
805 }
806
807
808 /******************************* endastrom ***********************************/
809 /*
810 Free astrometrical structures.
811 */
endastrom(picstruct * field)812 void endastrom(picstruct *field)
813
814 {
815 astromstruct *as;
816
817 as = field->astrom;
818 if (as->wcs_flag)
819 {
820 free(as->lin->cdelt);
821 free(as->lin->crpix);
822 free(as->lin->pc);
823 free(as->wcs);
824 free(as->lin);
825 free(as->cel);
826 free(as->prj);
827 free_tnxaxis(as->tnx_lngcor);
828 free_tnxaxis(as->tnx_latcor);
829 }
830
831 free(as);
832
833 return;
834 }
835
836
837 /********************************* precess ***********************************/
838 /*
839 precess equatorial coordinates according to the equinox (from Ephemerides du
840 Bureau des Longitudes 1992). Epoch for coordinates should be J2000
841 (FK5 system).
842 */
precess(double yearin,double alphain,double deltain,double yearout,double * alphaout,double * deltaout)843 void precess(double yearin, double alphain, double deltain,
844 double yearout, double *alphaout, double *deltaout)
845
846 {
847 double dzeta,theta,z, t1,t1t1, t2,t2t2,t2t2t2,
848 cddsadz, cddcadz, cdd, sdd, adz, cdin,sdin,ct,st,caindz;
849
850 alphain *= DEG;
851 deltain *= DEG;
852
853 t1 = (yearin - 2000.0)/1000.0;
854 t2 = (yearout - yearin)/1000.0;
855 t1t1 = t1*t1;
856 t2t2t2 = (t2t2 = t2*t2)*t2;
857 theta = (97171.735e-06 - 413.691e-06*t1 - 1.052e-06 * t1t1) * t2
858 + (-206.846e-06 - 1.052e-06*t1) * t2t2 - 202.812e-06 * t2t2t2;
859 dzeta = (111808.609e-06 + 677.071e-06*t1 - 0.674e-06 * t1t1) * t2
860 + (146.356e-06 - 1.673e-06*t1) * t2t2 + 87.257e-06 * t2t2t2;
861 z = (111808.609e-06 +677.071e-06*t1 - 0.674e-06 * t1t1) * t2
862 + (530.716e-06 + 0.320e-06*t1) * t2t2 + 88.251e-06 * t2t2t2;
863 cddsadz = (cdin=cos(deltain)) * sin(alphain+dzeta);
864 cddcadz = -(sdin=sin(deltain))*(st=sin(theta))
865 +cdin*(ct=cos(theta))*(caindz=cos(alphain+dzeta));
866 sdd = sdin*ct + cdin*st*caindz;
867 cdd = cos(*deltaout = asin(sdd));
868 adz = asin(cddsadz/cdd);
869 if (cddcadz<0.0)
870 adz = PI - adz;
871 if (adz<0.0)
872 adz += 2.0*PI;
873 adz += z;
874 *alphaout = adz/DEG;
875 *deltaout /= DEG;
876
877 return;
878 }
879
880
881 /*********************************** j2b *************************************/
882 /*
883 conver equatorial coordinates from equinox and epoch J2000 to equinox and
884 epoch B1950 for extragalactic sources (from Aoki et al. 1983, after
885 inversion of their matrix and some custom arrangements).
886 */
j2b(double yearobs,double alphain,double deltain,double * alphaout,double * deltaout)887 void j2b(double yearobs, double alphain, double deltain,
888 double *alphaout, double *deltaout)
889 {
890 int i,j;
891 static double a[3] = {-1.62557e-6, -0.31919e-6, -0.13843e-6},
892 ap[3] = {1.245e-3, -1.580e-3, -0.659e-3},
893 m[6][6] = {
894 { 0.9999256794678425, 0.01118148281196562, 0.004859003848996022,
895 -2.423898417033081e-06,-2.710547600126671e-08,-1.177738063266745e-08},
896 {-0.01118148272969232, 0.9999374849247641, -2.717708936468247e-05,
897 2.710547578707874e-08,-2.423927042585208e-06, 6.588254898401055e-11},
898 {-0.00485900399622881, -2.715579322970546e-05, 0.999988194643078,
899 1.177738102358923e-08, 6.582788892816657e-11,-2.424049920613325e-06},
900 {-0.0005508458576414713, 0.2384844384742432, -0.4356144527773499,
901 0.9999043171308133, 0.01118145410120206, 0.004858518651645554},
902 {-0.2385354433560954, -0.002664266996872802, 0.01225282765749546,
903 -0.01118145417187502, 0.9999161290795875, -2.717034576263522e-05},
904 { 0.4357269351676567, -0.008536768476441086, 0.002113420799663768,
905 -0.004858518477064975, -2.715994547222661e-05, 0.9999668385070383}},
906 a1[3], r[3], ro[3], r1[3], r2[3], v1[3], v[3];
907 static double cai, sai, cdi, sdi, dotp, rmod, alpha, delta, t1;
908
909 /* Convert Julian years from J2000.0 to tropic centuries from B1950.0 */
910 t1 = ((yearobs - 2000.0) + (MJD2000 - MJD1950)/365.25)*JU2TROP/100.0;
911 alphain *= DEG;
912 deltain *= DEG;
913 cai = cos(alphain);
914 sai = sin(alphain);
915 cdi = cos(deltain);
916 sdi = sin(deltain);
917 r[0] = cdi*cai;
918 r[1] = cdi*sai;
919 r[2] = sdi;
920 for (i=0; i<3; i++)
921 v[i] = r2[i] = v1[i] = 0.0;
922 for (j=0; j<6; j++)
923 for (i=0; i<6; i++)
924 if (j<3)
925 r2[j] += m[j][i]*(i<3?r[i]:v[i-3]);
926 else
927 v1[j-3] += m[j][i]*(i<3?r[i]:v[i-3]);
928
929 for (i=0; i<3; i++)
930 r1[i] = r2[i]+v1[i]*ARCSEC*t1;
931
932 dotp = 0.0;
933 for (i=0; i<3; i++)
934 {
935 a1[i] = a[i]+ap[i]*ARCSEC*t1;
936 dotp += a1[i]*(r1[i]+a1[i]);
937 }
938 dotp = 2.0/(sqrt(1+4.0*dotp)+1.0);
939 rmod = 0.0;
940 for (i=0; i<3; i++)
941 {
942 ro[i] = dotp*(r1[i]+a1[i]);
943 rmod += ro[i]*ro[i];
944 }
945 rmod = sqrt(rmod);
946 delta = asin(ro[2]/rmod);
947 alpha = acos(ro[0]/cos(delta)/rmod);
948 if (ro[1]<0)
949 alpha = 2.0*PI - alpha;
950 *alphaout = alpha/DEG;
951 *deltaout = delta/DEG;
952
953 return;
954 }
955
956