public software.sextractor

[/] [trunk/] [src/] [astrom.c] - Blame information for rev 235

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 233 bertin 
/*
2 
*                               astrom.c
3 2 bertin 
*
4 233 bertin 
* High level astrometric computations.
5 2 bertin 
*
6 233 bertin 
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 2 bertin 
*
8 233 bertin 
*       This file part of:      SExtractor
9 2 bertin 
*
10 235 bertin 
*       Copyright:              (C) 1993-2010 Emmanuel Bertin -- IAP/CNRS/UPMC
11 2 bertin 
*
12 233 bertin 
*       License:                GNU General Public License
13 
*
14 
*       SExtractor is free software: you can redistribute it and/or modify
15 
*       it under the terms of the GNU General Public License as published by
16 
*       the Free Software Foundation, either version 3 of the License, or
17 
*       (at your option) any later version.
18 
*       SExtractor is distributed in the hope that it will be useful,
19 
*       but WITHOUT ANY WARRANTY; without even the implied warranty of
20 
*       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
21 
*       GNU General Public License for more details.
22 
*       You should have received a copy of the GNU General Public License
23 
*       along with SExtractor. If not, see <http://www.gnu.org/licenses/>.
24 
*
25 
*       Last modified:          11/10/2010
26 
*
27 
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
28 2 bertin 
 
29 
#ifdef HAVE_CONFIG_H
30 
#include        "config.h"
31 
#endif
32 
 
33 
#include        <math.h>
34 
#include        <stdlib.h>
35 
#include        <string.h>
36 
 
37 
#include        "wcs/wcs.h"
38 
#include        "define.h"
39 
#include        "globals.h"
40 
#include        "prefs.h"
41 
#include        "astrom.h"
42 173 bertin 
#include        "fitswcs.h"
43 7 bertin 
#include        "wcs/tnx.h"
44 2 bertin 
 
45 
static obj2struct       *obj2 = &outobj2;
46 
 
47 
/****************************** initastrom **********************************/
48 
/*
49 
Initialize astrometrical structures.
50 
*/
51 
void    initastrom(picstruct *field)
52 
 
53 
  {
54 173 bertin 
   wcsstruct    *wcs;
55 2 bertin 
 
56 173 bertin 
  wcs = field->wcs;
57 2 bertin 
 
58 
/* Test if the WCS is in use */
59 173 bertin 
  if (wcs->lng != wcs->lat)
60 2 bertin 
    {
61 218 bertin 
    if (FLAG(obj2.dtheta2000) || FLAG(obj2.dtheta1950))
62 2 bertin 
      {
63 173 bertin 
      if (fabs(wcs->equinox-2000.0)>0.003)
64 
        precess(wcs->equinox, 0.0, 90.0, 2000.0, &wcs->ap2000, &wcs->dp2000);
65 2 bertin 
      else
66 
        {
67 173 bertin 
        wcs->ap2000 = 0.0;
68 
        wcs->dp2000 = 90.0;
69 2 bertin 
        }
70 
 
71 
      if (FLAG(obj2.theta1950) || FLAG(obj2.poserr_theta1950))
72 173 bertin 
        j2b(wcs->equinox, wcs->ap2000, wcs->dp2000, &wcs->ap1950, &wcs->dp1950);
73 2 bertin 
      }
74 
    }
75 
 
76 
/* Override astrometric definitions only if user supplies a pixel-scale */
77 
  if (prefs.pixel_scale == 0.0)
78 173 bertin 
    field->pixscale = wcs->pixscale*3600.0;     /* in arcsec */
79 2 bertin 
  else
80 173 bertin 
    field->pixscale = prefs.pixel_scale;
81 2 bertin 
 
82 
  return;
83 
  }
84 
 
85 
 
86 199 bertin 
/***************************** astrom_pos **********************************/
87 2 bertin 
/*
88 199 bertin 
Compute real FOCAL and WORLD coordinates according to FITS info.
89 2 bertin 
*/
90 199 bertin 
void    astrom_pos(picstruct *field, objstruct *obj)
91 2 bertin 
 
92 
  {
93 173 bertin 
   wcsstruct    *wcs;
94 
   double       rawpos[NAXIS], wcspos[NAXIS],
95 200 bertin 
                da,dd;
96 173 bertin 
   int          lng,lat;
97 2 bertin 
 
98 173 bertin 
  wcs = field->wcs;
99 
  lng = wcs->lng;
100 
  lat = wcs->lat;
101 2 bertin 
 
102 199 bertin 
/* If working with WCS, compute FOCAL coordinates and local matrix */
103 
  if (FLAG(obj2.mxf))
104 
    {
105 
    rawpos[0] = obj2->posx;
106 
    rawpos[1] = obj2->posy;
107 
    raw_to_red(wcs, rawpos, wcspos);
108 
    obj2->mxf = wcspos[0];
109 
    obj2->myf = wcspos[1];
110 
    }
111 
 
112 2 bertin 
/* If working with WCS, compute WORLD coordinates and local matrix */
113 
  if (FLAG(obj2.mxw))
114 
    {
115 173 bertin 
    rawpos[0] = obj2->posx;
116 
    rawpos[1] = obj2->posy;
117 
    raw_to_wcs(wcs, rawpos, wcspos);
118 
    obj2->mxw = wcspos[0];
119 
    obj2->myw = wcspos[1];
120 
    if (lng != lat)
121 2 bertin 
      {
122 173 bertin 
      obj2->alphas = lng<lat? obj2->mxw : obj2->myw;
123 
      obj2->deltas = lng<lat? obj2->myw : obj2->mxw;
124 2 bertin 
      if (FLAG(obj2.alpha2000))
125 
        {
126 173 bertin 
        if (fabs(wcs->equinox-2000.0)>0.003)
127 
          precess(wcs->equinox, wcspos[lng<lat?0:1], wcspos[lng<lat?1:0],
128 2 bertin 
                2000.0, &obj2->alpha2000, &obj2->delta2000);
129 
        else
130 
          {
131 173 bertin 
          obj2->alpha2000 = lng<lat? obj2->mxw : obj2->myw;
132 
          obj2->delta2000 = lng<lat? obj2->myw : obj2->mxw;
133 2 bertin 
          }
134 173 bertin 
        if (FLAG(obj2.dtheta2000))
135 
          {
136 
          da = wcs->ap2000 - obj2->alpha2000;
137 
          dd = (sin(wcs->dp2000*DEG)
138 
                -sin(obj2->delta2000*DEG)*sin(obj2->deltas*DEG))
139 
                /(cos(obj2->delta2000*DEG)*cos(obj2->deltas*DEG));
140 
          dd = dd<1.0? (dd>-1.0?acos(dd)/DEG:180.0) : 0.0;
141 
          obj2->dtheta2000 = (((da>0.0 && da<180.0) || da<-180.0)?-dd:dd);
142 
          }
143 2 bertin 
        if (FLAG(obj2.alpha1950))
144 173 bertin 
          {
145 
          j2b(wcs->equinox, obj2->alpha2000, obj2->delta2000,
146 2 bertin 
                &obj2->alpha1950, &obj2->delta1950);
147 173 bertin 
          if (FLAG(obj2.dtheta1950))
148 
            {
149 
            da = wcs->ap1950 - obj2->alpha1950;
150 
            dd = (sin(wcs->dp1950*DEG)
151 
                -sin(obj2->delta1950*DEG)*sin(obj2->deltas*DEG))
152 
                /(cos(obj2->delta1950*DEG)*cos(obj2->deltas*DEG));
153 
            dd = dd<1.0? (dd>-1.0?acos(dd)/DEG:180.0) : 0.0;
154 
            obj2->dtheta1950 = (((da>0.0 && da<180.0) || da<-180.0)?-dd:dd);
155 
           }
156 
          }
157 2 bertin 
        }
158 
      }
159 
    }
160 
 
161 199 bertin 
/* Custom coordinate system for the MAMA machine */
162 
  if (FLAG(obj2.mamaposx))
163 
    {
164 
    rawpos[0] = obj2->posx - 0.5;
165 
    rawpos[1] = obj2->posy - 0.5;
166 
    raw_to_wcs(wcs, rawpos, wcspos);
167 
    obj2->mamaposx = wcspos[1]*(MAMA_CORFLEX+1.0);
168 
    obj2->mamaposy = wcspos[0]*(MAMA_CORFLEX+1.0);
169 
    }
170 
 
171 
  return;
172 
  }
173 
 
174 
 
175 
/***************************** astrom_peakpos *******************************/
176 
/*
177 
Compute real FOCAL and WORLD peak coordinates according to FITS info.
178 
*/
179 
void    astrom_peakpos(picstruct *field, objstruct *obj)
180 
 
181 
  {
182 
   wcsstruct    *wcs;
183 
   double       rawpos[NAXIS], wcspos[NAXIS];
184 
   int          lng,lat;
185 
 
186 
  wcs = field->wcs;
187 
  lng = wcs->lng;
188 
  lat = wcs->lat;
189 
 
190 
  if (FLAG(obj2.peakxf))
191 
    {
192 
    rawpos[0] = obj->peakx;
193 
    rawpos[1] = obj->peaky;
194 
    raw_to_red(wcs, rawpos, wcspos);
195 
    obj2->peakxf = wcspos[0];
196 
    obj2->peakyf = wcspos[1];
197 
    }
198 2 bertin 
  if (FLAG(obj2.peakxw))
199 
    {
200 173 bertin 
    rawpos[0] = obj->peakx;
201 
    rawpos[1] = obj->peaky;
202 
    raw_to_wcs(wcs, rawpos, wcspos);
203 
    obj2->peakxw = wcspos[0];
204 
    obj2->peakyw = wcspos[1];
205 
    if (lng != lat)
206 2 bertin 
      {
207 173 bertin 
      obj2->peakalphas = lng<lat? obj2->peakxw : obj2->peakyw;
208 
      obj2->peakdeltas = lng<lat? obj2->peakyw : obj2->peakxw;
209 2 bertin 
      if (FLAG(obj2.peakalpha2000))
210 
        {
211 173 bertin 
        if (fabs(wcs->equinox-2000.0)>0.003)
212 
          precess(wcs->equinox, wcspos[lng<lat?0:1], wcspos[lng<lat?1:0],
213 2 bertin 
                2000.0, &obj2->peakalpha2000, &obj2->peakdelta2000);
214 
        else
215 
          {
216 173 bertin 
          obj2->peakalpha2000 = lng<lat? obj2->peakxw : obj2->peakyw;
217 
          obj2->peakdelta2000 = lng<lat? obj2->peakyw : obj2->peakxw;
218 2 bertin 
          }
219 
        if (FLAG(obj2.peakalpha1950))
220 173 bertin 
          j2b(wcs->equinox, obj2->peakalpha2000, obj2->peakdelta2000,
221 2 bertin 
                &obj2->peakalpha1950, &obj2->peakdelta1950);
222 
        }
223 
      }
224 
    }
225 
 
226 199 bertin 
  return;
227 
  }
228 
 
229 
 
230 
/****************************** astrom_winpos *******************************/
231 
/*
232 
Compute real FOCAL and WORLD windowed coordinates according to FITS info.
233 
*/
234 
void    astrom_winpos(picstruct *field, objstruct *obj)
235 
 
236 
  {
237 
   wcsstruct    *wcs;
238 
   double       rawpos[NAXIS], wcspos[NAXIS];
239 
   int          lng,lat;
240 
 
241 
  wcs = field->wcs;
242 
  lng = wcs->lng;
243 
  lat = wcs->lat;
244 
 
245 
  if (FLAG(obj2.winpos_xf))
246 
    {
247 
    rawpos[0] = obj2->winpos_x;
248 
    rawpos[1] = obj2->winpos_y;
249 
    raw_to_red(wcs, rawpos, wcspos);
250 
    obj2->winpos_xf = wcspos[0];
251 
    obj2->winpos_yf = wcspos[1];
252 
    }
253 2 bertin 
  if (FLAG(obj2.winpos_xw))
254 
    {
255 173 bertin 
    rawpos[0] = obj2->winpos_x;
256 
    rawpos[1] = obj2->winpos_y;
257 
    raw_to_wcs(wcs, rawpos, wcspos);
258 
    obj2->winpos_xw = wcspos[0];
259 
    obj2->winpos_yw = wcspos[1];
260 
    if (lng != lat)
261 2 bertin 
      {
262 173 bertin 
      obj2->winpos_alphas = lng<lat? obj2->winpos_xw : obj2->winpos_yw;
263 
      obj2->winpos_deltas = lng<lat? obj2->winpos_yw : obj2->winpos_xw;
264 2 bertin 
      if (FLAG(obj2.winpos_alpha2000))
265 
        {
266 173 bertin 
        if (fabs(wcs->equinox-2000.0)>0.003)
267 
          precess(wcs->equinox, wcspos[0], wcspos[1],
268 2 bertin 
                2000.0, &obj2->winpos_alpha2000, &obj2->winpos_delta2000);
269 
        else
270 
          {
271 173 bertin 
          obj2->winpos_alpha2000 = lng<lat? obj2->winpos_xw : obj2->winpos_yw;
272 
          obj2->winpos_delta2000 = lng<lat? obj2->winpos_yw : obj2->winpos_xw;
273 2 bertin 
          }
274 
        if (FLAG(obj2.winpos_alpha1950))
275 173 bertin 
          j2b(wcs->equinox, obj2->winpos_alpha2000, obj2->winpos_delta2000,
276 2 bertin 
                &obj2->winpos_alpha1950, &obj2->winpos_delta1950);
277 
        }
278 
      }
279 173 bertin 
    }
280 2 bertin 
 
281 199 bertin 
  return;
282 
  }
283 
 
284 
 
285 218 bertin 
/****************************** astrom_psfpos *******************************/
286 
/*
287 
Compute real FOCAL and WORLD PSF coordinates according to FITS info.
288 
*/
289 
void    astrom_psfpos(picstruct *field, objstruct *obj)
290 
 
291 
  {
292 
   wcsstruct    *wcs;
293 
   double       rawpos[NAXIS], wcspos[NAXIS];
294 
   int          lng,lat;
295 
 
296 
  wcs = field->wcs;
297 
  lng = wcs->lng;
298 
  lat = wcs->lat;
299 
 
300 
  if (FLAG(obj2.xf_psf))
301 
    {
302 
    rawpos[0] = obj2->x_psf;
303 
    rawpos[1] = obj2->y_psf;
304 
    raw_to_red(wcs, rawpos, wcspos);
305 
    obj2->xf_psf = wcspos[0];
306 
    obj2->yf_psf = wcspos[1];
307 
    }
308 
  if (FLAG(obj2.xw_psf))
309 
    {
310 
    rawpos[0] = obj2->x_psf;
311 
    rawpos[1] = obj2->y_psf;
312 
    raw_to_wcs(wcs, rawpos, wcspos);
313 
    obj2->xw_psf = wcspos[0];
314 
    obj2->yw_psf = wcspos[1];
315 
    if (lng != lat)
316 
      {
317 
      obj2->alphas_psf = lng<lat? obj2->xw_psf : obj2->yw_psf;
318 
      obj2->deltas_psf = lng<lat? obj2->yw_psf : obj2->xw_psf;
319 
      if (FLAG(obj2.alpha2000_psf))
320 
        {
321 
        if (fabs(wcs->equinox-2000.0)>0.003)
322 
          precess(wcs->equinox, wcspos[0], wcspos[1],
323 
                2000.0, &obj2->alpha2000_psf, &obj2->delta2000_psf);
324 
        else
325 
          {
326 
          obj2->alpha2000_psf = lng<lat? obj2->xw_psf : obj2->yw_psf;
327 
          obj2->delta2000_psf = lng<lat? obj2->yw_psf : obj2->xw_psf;
328 
          }
329 
        if (FLAG(obj2.alpha1950_psf))
330 
          j2b(wcs->equinox, obj2->alpha2000_psf, obj2->delta2000_psf,
331 
                &obj2->alpha1950_psf, &obj2->delta1950_psf);
332 
        }
333 
      }
334 
    }
335 
 
336 
  return;
337 
  }
338 
 
339 
 
340 199 bertin 
/****************************** astrom_profpos *******************************/
341 
/*
342 
Compute real FOCAL and WORLD profit coordinates according to FITS info.
343 
*/
344 
void    astrom_profpos(picstruct *field, objstruct *obj)
345 
 
346 
  {
347 
   wcsstruct    *wcs;
348 
   double       rawpos[NAXIS], wcspos[NAXIS];
349 
   int          lng,lat;
350 
 
351 
  wcs = field->wcs;
352 
  lng = wcs->lng;
353 
  lat = wcs->lat;
354 
 
355 
  if (FLAG(obj2.xf_prof))
356 
    {
357 
    rawpos[0] = obj2->x_prof;
358 
    rawpos[1] = obj2->y_prof;
359 
    raw_to_red(wcs, rawpos, wcspos);
360 
    obj2->xf_prof = wcspos[0];
361 
    obj2->yf_prof = wcspos[1];
362 
    }
363 173 bertin 
  if (FLAG(obj2.xw_prof))
364 
    {
365 
    rawpos[0] = obj2->x_prof;
366 
    rawpos[1] = obj2->y_prof;
367 
    raw_to_wcs(wcs, rawpos, wcspos);
368 
    obj2->xw_prof = wcspos[0];
369 
    obj2->yw_prof = wcspos[1];
370 
    if (lng != lat)
371 
      {
372 
      obj2->alphas_prof = lng<lat? obj2->xw_prof : obj2->yw_prof;
373 
      obj2->deltas_prof = lng<lat? obj2->yw_prof : obj2->xw_prof;
374 
      if (FLAG(obj2.alpha2000_prof))
375 
        {
376 
        if (fabs(wcs->equinox-2000.0)>0.003)
377 
          precess(wcs->equinox, wcspos[0], wcspos[1],
378 
                2000.0, &obj2->alpha2000_prof, &obj2->delta2000_prof);
379 
        else
380 
          {
381 
          obj2->alpha2000_prof = lng<lat? obj2->xw_prof : obj2->yw_prof;
382 
          obj2->delta2000_prof = lng<lat? obj2->yw_prof : obj2->xw_prof;
383 
          }
384 
        if (FLAG(obj2.alpha1950_prof))
385 
          j2b(wcs->equinox, obj2->alpha2000_prof, obj2->delta2000_prof,
386 
                &obj2->alpha1950_prof, &obj2->delta1950_prof);
387 
        }
388 2 bertin 
      }
389 
    }
390 
 
391 
  return;
392 
  }
393 
 
394 
 
395 
/****************************** astrom_shapeparam ****************************/
396 
/*
397 
Compute shape parameters in WORLD and SKY coordinates.
398 
*/
399 
void    astrom_shapeparam(picstruct *field, objstruct *obj)
400 
  {
401 173 bertin 
   wcsstruct    *wcs;
402 
   double       dx2,dy2,dxy, xm2,ym2,xym, temp,pm2, lm0,lm1,lm2,lm3;
403 7 bertin 
   int          lng,lat, naxis;
404 2 bertin 
 
405 173 bertin 
  wcs = field->wcs;
406 
  naxis = wcs->naxis;
407 
  lng = wcs->lng;
408 
  lat = wcs->lat;
409 7 bertin 
  if (lng == lat)
410 
    {
411 
    lng = 0;
412 
    lat = 1;
413 
    }
414 173 bertin 
  lm0 = obj2->jacob[lng+naxis*lng];
415 
  lm1 = obj2->jacob[lat+naxis*lng];
416 
  lm2 = obj2->jacob[lng+naxis*lat];
417 
  lm3 = obj2->jacob[lat+naxis*lat];
418 7 bertin 
 
419 
 
420 2 bertin 
/* All WORLD params based on 2nd order moments have to pass through here */
421 
  dx2 = obj->mx2;
422 
  dy2 = obj->my2;
423 
  dxy = obj->mxy;
424 7 bertin 
  obj2->mx2w = xm2 = lm0*lm0*dx2 + lm1*lm1*dy2 + lm0*lm1*dxy;
425 
  obj2->my2w = ym2 = lm2*lm2*dx2 + lm3*lm3*dy2 + lm2*lm3*dxy;
426 
  obj2->mxyw = xym = lm0*lm2*dx2 + lm1*lm3*dy2 + (lm0*lm3+lm1*lm2)*dxy;
427 2 bertin 
  temp=xm2-ym2;
428 
  if (FLAG(obj2.thetaw))
429 
    {
430 173 bertin 
    obj2->thetaw = fmod_m90_p90((temp == 0.0)?
431 
                                (45.0) : (0.5*atan2(2.0 * xym,temp)/DEG));
432 2 bertin 
 
433 
/*-- Compute position angles in J2000 or B1950 reference frame */
434 173 bertin 
    if (wcs->lng != wcs->lat)
435 2 bertin 
      {
436 173 bertin 
      if (FLAG(obj2.thetas))
437 
        obj2->thetas = fmod_m90_p90(lng<lat?
438 
                                ((obj2->thetaw>0.0?90:-90.0) - obj2->thetaw)
439 
                                : obj2->thetaw);
440 2 bertin 
      if (FLAG(obj2.theta2000))
441 173 bertin 
        obj2->theta2000 = fmod_m90_p90(obj2->thetas + obj2->dtheta2000);
442 2 bertin 
      if (FLAG(obj2.theta1950))
443 173 bertin 
        obj2->theta1950 = fmod_m90_p90(obj2->thetas + obj2->dtheta1950);
444 2 bertin 
      }
445 
    }
446 
 
447 
  if (FLAG(obj2.aw))
448 
    {
449 
    temp = sqrt(0.25*temp*temp+xym*xym);
450 
    pm2 = 0.5*(xm2+ym2);
451 
    obj2->aw = (float)sqrt(pm2+temp);
452 
    obj2->bw = (float)sqrt(pm2-temp);
453 
    obj2->polarw = temp / pm2;
454 
    }
455 
 
456 
  if (FLAG(obj2.cxxw))
457 
    {
458 
/*-- Handle large, fully correlated profiles (can cause a singularity...) */
459 
    if ((temp=xm2*ym2-xym*xym)<1e-6)
460 
      {
461 
      temp = 1e-6;
462 
      xym *= 0.99999;
463 
      }
464 
    obj2->cxxw = (float)(ym2/temp);
465 
    obj2->cyyw = (float)(xm2/temp);
466 
    obj2->cxyw = (float)(-2*xym/temp);
467 
    }
468 
 
469 
  return;
470 
  }
471 
 
472 
 
473 
/**************************** astrom_winshapeparam ***************************/
474 
/*
475 
Compute shape parameters in WORLD and SKY coordinates.
476 
*/
477 
void    astrom_winshapeparam(picstruct *field, objstruct *obj)
478 
  {
479 173 bertin 
   wcsstruct    *wcs;
480 
   double       dx2,dy2,dxy, xm2,ym2,xym, temp,pm2, lm0,lm1,lm2,lm3;
481 7 bertin 
   int          lng,lat, naxis;
482 2 bertin 
 
483 173 bertin 
  wcs = field->wcs;
484 
  naxis = wcs->naxis;
485 
  lng = wcs->lng;
486 
  lat = wcs->lat;
487 7 bertin 
  if (lng == lat)
488 
    {
489 
    lng = 0;
490 
    lat = 1;
491 
    }
492 173 bertin 
  lm0 = obj2->jacob[lng+naxis*lng];
493 
  lm1 = obj2->jacob[lat+naxis*lng];
494 
  lm2 = obj2->jacob[lng+naxis*lat];
495 
  lm3 = obj2->jacob[lat+naxis*lat];
496 7 bertin 
 
497 2 bertin 
/* All WORLD params based on 2nd order moments have to pass through here */
498 
  dx2 = obj2->win_mx2;
499 
  dy2 = obj2->win_my2;
500 
  dxy = obj2->win_mxy;
501 7 bertin 
  obj2->win_mx2w = xm2 = lm0*lm0*dx2 + lm1*lm1*dy2 + lm0*lm1*dxy;
502 
  obj2->win_my2w = ym2 = lm2*lm2*dx2 + lm3*lm3*dy2 + lm2*lm3*dxy;
503 
  obj2->win_mxyw = xym = lm0*lm2*dx2 + lm1*lm3*dy2 + (lm0*lm3+lm1*lm2)*dxy;
504 2 bertin 
  temp=xm2-ym2;
505 
  if (FLAG(obj2.win_thetaw))
506 
    {
507 173 bertin 
    obj2->win_thetaw = fmod_m90_p90((temp == 0.0)?
508 
                        (45.0) : (0.5*atan2(2.0*xym,temp)/DEG));
509 2 bertin 
 
510 
/*-- Compute position angles in J2000 or B1950 reference frame */
511 173 bertin 
    if (wcs->lng != wcs->lat)
512 2 bertin 
      {
513 173 bertin 
      if (FLAG(obj2.win_thetas))
514 
        obj2->win_thetas = fmod_m90_p90(lng<lat?
515 
                        ((obj2->win_thetaw>0.0?90:-90.0) - obj2->win_thetaw)
516 
                        : obj2->win_thetaw);
517 2 bertin 
      if (FLAG(obj2.win_theta2000))
518 173 bertin 
        obj2->win_theta2000 = fmod_m90_p90(obj2->win_thetas + obj2->dtheta2000);
519 2 bertin 
      if (FLAG(obj2.win_theta1950))
520 173 bertin 
        obj2->win_theta1950 = fmod_m90_p90(obj2->win_thetas + obj2->dtheta1950);
521 2 bertin 
      }
522 
    }
523 
 
524 
  if (FLAG(obj2.win_aw))
525 
    {
526 
    temp = sqrt(0.25*temp*temp+xym*xym);
527 
    pm2 = 0.5*(xm2+ym2);
528 3 bertin 
    obj2->win_aw = (float)sqrt(pm2+temp);
529 
    obj2->win_bw = (float)sqrt(pm2-temp);
530 2 bertin 
    obj2->win_polarw = temp / pm2;
531 
    }
532 
 
533 
  if (FLAG(obj2.win_cxxw))
534 
    {
535 
/*-- Handle large, fully correlated profiles (can cause a singularity...) */
536 
    if ((temp=xm2*ym2-xym*xym)<1e-6)
537 
      {
538 
      temp = 1e-6;
539 
      xym *= 0.99999;
540 
      }
541 3 bertin 
    obj2->win_cxxw = (float)(ym2/temp);
542 
    obj2->win_cyyw = (float)(xm2/temp);
543 
    obj2->win_cxyw = (float)(-2*xym/temp);
544 2 bertin 
    }
545 
 
546 
  return;
547 
  }
548 
 
549 
 
550 
/******************************* astrom_errparam *****************************/
551 
/*
552 
Compute error ellipse parameters in WORLD and SKY coordinates.
553 
*/
554 
void    astrom_errparam(picstruct *field, objstruct *obj)
555 
  {
556 173 bertin 
   wcsstruct    *wcs;
557 
   double       dx2,dy2,dxy, xm2,ym2,xym, temp,pm2, lm0,lm1,lm2,lm3;
558 7 bertin 
   int          lng,lat, naxis;
559 2 bertin 
 
560 173 bertin 
  wcs = field->wcs;
561 
  naxis = wcs->naxis;
562 
  lng = wcs->lng;
563 
  lat = wcs->lat;
564 7 bertin 
  if (lng == lat)
565 
    {
566 
    lng = 0;
567 
    lat = 1;
568 
    }
569 173 bertin 
  lm0 = obj2->jacob[lng+naxis*lng];
570 
  lm1 = obj2->jacob[lat+naxis*lng];
571 
  lm2 = obj2->jacob[lng+naxis*lat];
572 
  lm3 = obj2->jacob[lat+naxis*lat];
573 7 bertin 
 
574 2 bertin 
/* All WORLD params based on 2nd order moments have to pass through here */
575 
  dx2 = obj->poserr_mx2;
576 
  dy2 = obj->poserr_my2;
577 
  dxy = obj->poserr_mxy;
578 7 bertin 
  obj2->poserr_mx2w = xm2 = lm0*lm0*dx2 + lm1*lm1*dy2 + lm0*lm1*dxy;
579 
  obj2->poserr_my2w = ym2 = lm2*lm2*dx2 + lm3*lm3*dy2 + lm2*lm3*dxy;
580 
  obj2->poserr_mxyw = xym = lm0*lm2*dx2 + lm1*lm3*dy2 + (lm0*lm3+lm1*lm2)*dxy;
581 2 bertin 
  temp=xm2-ym2;
582 
  if (FLAG(obj2.poserr_thetaw))
583 
    {
584 173 bertin 
    obj2->poserr_thetaw = fmod_m90_p90((temp==0.0)?
585 
                                (45.0):(0.5*atan2(2.0*xym,temp)/DEG));
586 2 bertin 
 
587 
/*-- Compute position angles in J2000 or B1950 reference frame */
588 173 bertin 
    if (wcs->lng != wcs->lat)
589 2 bertin 
      {
590 173 bertin 
      if (FLAG(obj2.poserr_thetas))
591 
        obj2->poserr_thetas = fmod_m90_p90(lng<lat?
592 
                ((obj2->poserr_thetaw>0.0?90:-90.0) - obj2->poserr_thetaw)
593 
                : obj2->poserr_thetaw);
594 2 bertin 
      if (FLAG(obj2.poserr_theta2000))
595 173 bertin 
        obj2->poserr_theta2000 = fmod_m90_p90(obj2->poserr_thetas
596 
                                + obj2->dtheta2000);
597 2 bertin 
      if (FLAG(obj2.poserr_theta1950))
598 173 bertin 
        obj2->poserr_theta1950 = fmod_m90_p90(obj2->poserr_thetas
599 
                                + obj2->dtheta1950);
600 2 bertin 
      }
601 
    }
602 
 
603 
  if (FLAG(obj2.poserr_aw))
604 
    {
605 
    temp = sqrt(0.25*temp*temp+xym*xym);
606 
    pm2 = 0.5*(xm2+ym2);
607 
    obj2->poserr_aw = (float)sqrt(pm2+temp);
608 
    obj2->poserr_bw = (float)sqrt(pm2-temp);
609 
    }
610 
 
611 
  if (FLAG(obj2.poserr_cxxw))
612 
    {
613 
/*-- Handle large, fully correlated profiles (can cause a singularity...) */
614 
    if ((temp=xm2*ym2-xym*xym)<1e-6)
615 
      {
616 
      temp = 1e-6;
617 
      xym *= 0.99999;
618 
      }
619 
    obj2->poserr_cxxw = (float)(ym2/temp);
620 
    obj2->poserr_cyyw = (float)(xm2/temp);
621 
    obj2->poserr_cxyw = (float)(-2*xym/temp);
622 
    }
623 
 
624 
  return;
625 
  }
626 
 
627 
 
628 
/***************************** astrom_winerrparam ***************************/
629 
/*
630 
Compute error ellipse parameters in WORLD and SKY coordinates.
631 
*/
632 
void    astrom_winerrparam(picstruct *field, objstruct *obj)
633 
  {
634 173 bertin 
   wcsstruct    *wcs;
635 
   double       dx2,dy2,dxy, xm2,ym2,xym, temp,pm2, lm0,lm1,lm2,lm3;
636 7 bertin 
   int          lng,lat, naxis;
637 2 bertin 
 
638 173 bertin 
  wcs = field->wcs;
639 
  naxis = wcs->naxis;
640 
  lng = wcs->lng;
641 
  lat = wcs->lat;
642 7 bertin 
  if (lng == lat)
643 
    {
644 
    lng = 0;
645 
    lat = 1;
646 
    }
647 173 bertin 
  lm0 = obj2->jacob[lng+naxis*lng];
648 
  lm1 = obj2->jacob[lat+naxis*lng];
649 
  lm2 = obj2->jacob[lng+naxis*lat];
650 
  lm3 = obj2->jacob[lat+naxis*lat];
651 7 bertin 
 
652 2 bertin 
/* All WORLD params based on 2nd order moments have to pass through here */
653 
  dx2 = obj2->winposerr_mx2;
654 
  dy2 = obj2->winposerr_my2;
655 
  dxy = obj2->winposerr_mxy;
656 7 bertin 
  obj2->winposerr_mx2w = xm2 = lm0*lm0*dx2 + lm1*lm1*dy2 + lm0*lm1*dxy;
657 
  obj2->winposerr_my2w = ym2 = lm2*lm2*dx2 + lm3*lm3*dy2 + lm2*lm3*dxy;
658 
  obj2->winposerr_mxyw = xym = lm0*lm2*dx2 + lm1*lm3*dy2 + (lm0*lm3+lm1*lm2)*dxy;
659 2 bertin 
  temp=xm2-ym2;
660 
  if (FLAG(obj2.winposerr_thetaw))
661 
    {
662 173 bertin 
    obj2->winposerr_thetaw = (fmod_m90_p90(temp==0.0)?
663 
                                (45.0):(0.5*atan2(2.0*xym,temp)/DEG));
664 2 bertin 
 
665 
/*-- Compute position angles in J2000 or B1950 reference frame */
666 173 bertin 
    if (wcs->lng != wcs->lat)
667 2 bertin 
      {
668 173 bertin 
      if (FLAG(obj2.winposerr_thetas))
669 
        obj2->winposerr_thetas = fmod_m90_p90(lng<lat?
670 
                ((obj2->winposerr_thetaw>0.0?90:-90.0) - obj2->winposerr_thetaw)
671 
                : obj2->winposerr_thetaw);
672 2 bertin 
      if (FLAG(obj2.winposerr_theta2000))
673 173 bertin 
        obj2->winposerr_theta2000 = fmod_m90_p90(obj2->winposerr_thetas
674 
                                + obj2->dtheta2000);
675 2 bertin 
      if (FLAG(obj2.winposerr_theta1950))
676 173 bertin 
        obj2->winposerr_theta1950 = fmod_m90_p90(obj2->winposerr_thetas
677 
                                + obj2->dtheta1950);
678 2 bertin 
      }
679 
    }
680 
 
681 
  if (FLAG(obj2.winposerr_aw))
682 
    {
683 
    temp = sqrt(0.25*temp*temp+xym*xym);
684 
    pm2 = 0.5*(xm2+ym2);
685 
    obj2->winposerr_aw = (float)sqrt(pm2+temp);
686 
    obj2->winposerr_bw = (float)sqrt(pm2-temp);
687 
    }
688 
 
689 
  if (FLAG(obj2.winposerr_cxxw))
690 
    {
691 
/*-- Handle large, fully correlated profiles (can cause a singularity...) */
692 
    if ((temp=xm2*ym2-xym*xym)<1e-6)
693 
      {
694 
      temp = 1e-6;
695 
      xym *= 0.99999;
696 
      }
697 
    obj2->winposerr_cxxw = (float)(ym2/temp);
698 
    obj2->winposerr_cyyw = (float)(xm2/temp);
699 
    obj2->winposerr_cxyw = (float)(-2*xym/temp);
700 
    }
701 
 
702 
  return;
703 
  }
704 
 
705 
 
706 218 bertin 
/***************************** astrom_psferrparam ***************************/
707 
/*
708 
Compute error ellipse parameters in WORLD and SKY coordinates.
709 
*/
710 
void    astrom_psferrparam(picstruct *field, objstruct *obj)
711 
  {
712 
   wcsstruct    *wcs;
713 
   double       dx2,dy2,dxy, xm2,ym2,xym, temp,pm2, lm0,lm1,lm2,lm3;
714 
   int          lng,lat, naxis;
715 
 
716 
  wcs = field->wcs;
717 
  naxis = wcs->naxis;
718 
  lng = wcs->lng;
719 
  lat = wcs->lat;
720 
  if (lng == lat)
721 
    {
722 
    lng = 0;
723 
    lat = 1;
724 
    }
725 
  lm0 = obj2->jacob[lng+naxis*lng];
726 
  lm1 = obj2->jacob[lat+naxis*lng];
727 
  lm2 = obj2->jacob[lng+naxis*lat];
728 
  lm3 = obj2->jacob[lat+naxis*lat];
729 
 
730 
/* All WORLD params based on 2nd order moments have to pass through here */
731 
  dx2 = obj2->poserrmx2_psf;
732 
  dy2 = obj2->poserrmy2_psf;
733 
  dxy = obj2->poserrmxy_psf;
734 
  obj2->poserrmx2w_psf = xm2 = lm0*lm0*dx2+lm1*lm1*dy2+lm0*lm1*dxy;
735 
  obj2->poserrmy2w_psf = ym2 = lm2*lm2*dx2+lm3*lm3*dy2+lm2*lm3*dxy;
736 
  obj2->poserrmxyw_psf = xym = lm0*lm2*dx2+lm1*lm3*dy2+(lm0*lm3+lm1*lm2)*dxy;
737 
  temp=xm2-ym2;
738 
  if (FLAG(obj2.poserrthetaw_psf))
739 
    {
740 
    obj2->poserrthetaw_psf = (fmod_m90_p90(temp==0.0)?
741 
                                (45.0):(0.5*atan2(2.0*xym,temp)/DEG));
742 
 
743 
/*-- Compute position angles in J2000 or B1950 reference frame */
744 
    if (wcs->lng != wcs->lat)
745 
      {
746 
      if (FLAG(obj2.poserrthetas_psf))
747 
        obj2->poserrthetas_psf = fmod_m90_p90(lng<lat?
748 
                ((obj2->poserrthetaw_psf>0.0?90:-90.0)-obj2->poserrthetaw_psf)
749 
                : obj2->poserrthetaw_psf);
750 
      if (FLAG(obj2.poserrtheta2000_psf))
751 
        obj2->poserrtheta2000_psf = fmod_m90_p90(obj2->poserrthetas_psf
752 
                                + obj2->dtheta2000);
753 
      if (FLAG(obj2.poserrtheta1950_psf))
754 
        obj2->poserrtheta1950_psf = fmod_m90_p90(obj2->poserrthetas_psf
755 
                                + obj2->dtheta1950);
756 
      }
757 
    }
758 
 
759 
  if (FLAG(obj2.poserraw_psf))
760 
    {
761 
    temp = sqrt(0.25*temp*temp+xym*xym);
762 
    pm2 = 0.5*(xm2+ym2);
763 
    obj2->poserraw_psf = (float)sqrt(pm2+temp);
764 
    obj2->poserrbw_psf = (float)sqrt(pm2-temp);
765 
    }
766 
 
767 
  if (FLAG(obj2.poserrcxxw_psf))
768 
    {
769 
/*-- Handle large, fully correlated profiles (can cause a singularity...) */
770 
    if ((temp=xm2*ym2-xym*xym)<1e-6)
771 
      {
772 
      temp = 1e-6;
773 
      xym *= 0.99999;
774 
      }
775 
    obj2->poserrcxxw_psf = (float)(ym2/temp);
776 
    obj2->poserrcyyw_psf = (float)(xm2/temp);
777 
    obj2->poserrcxyw_psf = (float)(-2*xym/temp);
778 
    }
779 
 
780 
  return;
781 
  }
782 
 
783 
 
784 173 bertin 
/***************************** astrom_proferrparam ***************************/
785 2 bertin 
/*
786 173 bertin 
Compute error ellipse parameters in WORLD and SKY coordinates.
787 2 bertin 
*/
788 173 bertin 
void    astrom_proferrparam(picstruct *field, objstruct *obj)
789 2 bertin 
  {
790 173 bertin 
   wcsstruct    *wcs;
791 
   double       dx2,dy2,dxy, xm2,ym2,xym, temp,pm2, lm0,lm1,lm2,lm3;
792 
   int          lng,lat, naxis;
793 2 bertin 
 
794 173 bertin 
  wcs = field->wcs;
795 
  naxis = wcs->naxis;
796 
  lng = wcs->lng;
797 
  lat = wcs->lat;
798 
  if (lng == lat)
799 2 bertin 
    {
800 173 bertin 
    lng = 0;
801 
    lat = 1;
802 
    }
803 
  lm0 = obj2->jacob[lng+naxis*lng];
804 
  lm1 = obj2->jacob[lat+naxis*lng];
805 
  lm2 = obj2->jacob[lng+naxis*lat];
806 
  lm3 = obj2->jacob[lat+naxis*lat];
807 
 
808 
/* All WORLD params based on 2nd order moments have to pass through here */
809 
  dx2 = obj2->poserrmx2_prof;
810 
  dy2 = obj2->poserrmy2_prof;
811 
  dxy = obj2->poserrmxy_prof;
812 
  obj2->poserrmx2w_prof = xm2 = lm0*lm0*dx2+lm1*lm1*dy2+lm0*lm1*dxy;
813 
  obj2->poserrmy2w_prof = ym2 = lm2*lm2*dx2+lm3*lm3*dy2+lm2*lm3*dxy;
814 
  obj2->poserrmxyw_prof = xym = lm0*lm2*dx2+lm1*lm3*dy2+(lm0*lm3+lm1*lm2)*dxy;
815 
  temp=xm2-ym2;
816 
  if (FLAG(obj2.poserrthetaw_prof))
817 
    {
818 
    obj2->poserrthetaw_prof = (fmod_m90_p90(temp==0.0)?
819 
                                (45.0):(0.5*atan2(2.0*xym,temp)/DEG));
820 
 
821 
/*-- Compute position angles in J2000 or B1950 reference frame */
822 
    if (wcs->lng != wcs->lat)
823 2 bertin 
      {
824 173 bertin 
      if (FLAG(obj2.poserrthetas_prof))
825 
        obj2->poserrthetas_prof = fmod_m90_p90(lng<lat?
826 
                ((obj2->poserrthetaw_prof>0.0?90:-90.0)-obj2->poserrthetaw_prof)
827 
                : obj2->poserrthetaw_prof);
828 
      if (FLAG(obj2.poserrtheta2000_prof))
829 
        obj2->poserrtheta2000_prof = fmod_m90_p90(obj2->poserrthetas_prof
830 
                                + obj2->dtheta2000);
831 
      if (FLAG(obj2.poserrtheta1950_prof))
832 
        obj2->poserrtheta1950_prof = fmod_m90_p90(obj2->poserrthetas_prof
833 
                                + obj2->dtheta1950);
834 2 bertin 
      }
835 
    }
836 
 
837 173 bertin 
  if (FLAG(obj2.poserraw_prof))
838 
    {
839 
    temp = sqrt(0.25*temp*temp+xym*xym);
840 
    pm2 = 0.5*(xm2+ym2);
841 
    obj2->poserraw_prof = (float)sqrt(pm2+temp);
842 
    obj2->poserrbw_prof = (float)sqrt(pm2-temp);
843 
    }
844 2 bertin 
 
845 173 bertin 
  if (FLAG(obj2.poserrcxxw_prof))
846 2 bertin 
    {
847 173 bertin 
/*-- Handle large, fully correlated profiles (can cause a singularity...) */
848 
    if ((temp=xm2*ym2-xym*xym)<1e-6)
849 
      {
850 
      temp = 1e-6;
851 
      xym *= 0.99999;
852 
      }
853 
    obj2->poserrcxxw_prof = (float)(ym2/temp);
854 
    obj2->poserrcyyw_prof = (float)(xm2/temp);
855 
    obj2->poserrcxyw_prof = (float)(-2*xym/temp);
856 2 bertin 
    }
857 
 
858 
  return;
859 
  }
860 
 
861 
 
862 173 bertin 
/*************************** astrom_profshapeparam ***************************/
863 2 bertin 
/*
864 173 bertin 
Compute profile-fitting shape parameters in WORLD and SKY coordinates.
865 2 bertin 
*/
866 173 bertin 
void    astrom_profshapeparam(picstruct *field, objstruct *obj)
867 2 bertin 
  {
868 173 bertin 
   wcsstruct    *wcs;
869 226 bertin 
   double       mat[9], tempmat[9], mx2wcov[9], dpdmx2[6], cov[4],
870 
                dx2,dy2,dxy, xm2,ym2,xym, pm2, lm0,lm1,lm2,lm3, ct,st,
871 
                temp, invstemp, den, invden, dval;
872 173 bertin 
   int          lng,lat, naxis;
873 2 bertin 
 
874 173 bertin 
  wcs = field->wcs;
875 
  naxis = wcs->naxis;
876 
  lng = wcs->lng;
877 
  lat = wcs->lat;
878 
  if (lng == lat)
879 
    {
880 
    lng = 0;
881 
    lat = 1;
882 
    }
883 
  lm0 = obj2->jacob[lng+naxis*lng];
884 
  lm1 = obj2->jacob[lat+naxis*lng];
885 
  lm2 = obj2->jacob[lng+naxis*lat];
886 
  lm3 = obj2->jacob[lat+naxis*lat];
887 2 bertin 
 
888 173 bertin 
/* Spheroid World coordinates */
889 
  obj2->prof_spheroid_thetaerrw=obj2->prof_spheroid_thetaerr; /* quick & dirty*/
890 
  if (FLAG(obj2.prof_spheroid_reffw))
891 
    {
892 
    ct = cos(obj2->prof_spheroid_theta*DEG);
893 
    st = sin(obj2->prof_spheroid_theta*DEG);
894 
    dx2 = obj2->prof_spheroid_reff*obj2->prof_spheroid_reff * (ct*ct
895 
        + st*st * obj2->prof_spheroid_aspect*obj2->prof_spheroid_aspect);
896 
    dy2 = obj2->prof_spheroid_reff*obj2->prof_spheroid_reff * (st*st
897 
        + ct*ct * obj2->prof_spheroid_aspect*obj2->prof_spheroid_aspect);
898 
    dxy = ct*st * obj2->prof_spheroid_reff*obj2->prof_spheroid_reff
899 
        *(1.0 - obj2->prof_spheroid_aspect*obj2->prof_spheroid_aspect);
900 
    xm2 = lm0*lm0*dx2 + lm1*lm1*dy2 + lm0*lm1*dxy;
901 
    ym2 = lm2*lm2*dx2 + lm3*lm3*dy2 + lm2*lm3*dxy;
902 
    xym = lm0*lm2*dx2 + lm1*lm3*dy2 + (lm0*lm3+lm1*lm2)*dxy;
903 
    temp=xm2-ym2;
904 
    if (FLAG(obj2.prof_spheroid_thetaw))
905 
      {
906 
      obj2->prof_spheroid_thetaw = fmod_m90_p90((temp == 0.0)?
907 
                (45.0) : (0.5*atan2(2.0 * xym,temp)/DEG));
908 2 bertin 
 
909 173 bertin 
      if (wcs->lng != wcs->lat)
910 
        {
911 
        if (FLAG(obj2.prof_spheroid_thetas))
912 
          obj2->prof_spheroid_thetas = fmod_m90_p90(lng<lat?
913 
                        ((obj2->prof_spheroid_thetaw>0.0?90:-90.0)
914 
                                - obj2->prof_spheroid_thetaw)
915 
                        : obj2->prof_spheroid_thetaw);
916 
        if (FLAG(obj2.prof_spheroid_theta2000))
917 
          obj2->prof_spheroid_theta2000=fmod_m90_p90(obj2->prof_spheroid_thetas
918 
                                        + obj2->dtheta2000);
919 
        if (FLAG(obj2.prof_spheroid_theta1950))
920 
          obj2->prof_spheroid_theta1950=fmod_m90_p90(obj2->prof_spheroid_thetas
921 
                                        + obj2->dtheta1950);
922 
        }
923 
      }
924 
    temp = sqrt(0.25*temp*temp+xym*xym);
925 
    pm2 = 0.5*(xm2+ym2);
926 
    obj2->prof_spheroid_reffw = sqrt(pm2+temp);
927 
    obj2->prof_spheroid_refferrw = obj2->prof_spheroid_reff > 0.0?
928 
        obj2->prof_spheroid_refferr/obj2->prof_spheroid_reff
929 
                                *obj2->prof_spheroid_reffw
930 
        : 0.0;
931 233 bertin 
    obj2->prof_spheroid_aspectw = (obj2->prof_spheroid_reffw>0.0 && pm2>temp)?
932 173 bertin 
                                  sqrt(pm2-temp) / obj2->prof_spheroid_reffw
933 
                                : obj2->prof_spheroid_aspect;
934 
    obj2->prof_spheroid_aspecterrw = obj2->prof_spheroid_aspect > 0.0?
935 
        obj2->prof_spheroid_aspecterr/obj2->prof_spheroid_aspect
936 
                                *obj2->prof_spheroid_aspectw
937 
        : 0.0;
938 
    }
939 2 bertin 
 
940 173 bertin 
/* Disk World coordinates */
941 
  obj2->prof_disk_thetaerrw = obj2->prof_disk_thetaerr; /* quick & dirty*/
942 
  if (FLAG(obj2.prof_disk_scalew))
943 
    {
944 
    ct = cos(obj2->prof_disk_theta*DEG);
945 
    st = sin(obj2->prof_disk_theta*DEG);
946 
    dx2 = obj2->prof_disk_scale*obj2->prof_disk_scale * (ct*ct
947 
        + st*st * obj2->prof_disk_aspect*obj2->prof_disk_aspect);
948 
    dy2 = obj2->prof_disk_scale*obj2->prof_disk_scale * (st*st
949 
        + ct*ct * obj2->prof_disk_aspect*obj2->prof_disk_aspect);
950 
    dxy = ct*st * obj2->prof_disk_scale*obj2->prof_disk_scale
951 
        *(1.0 - obj2->prof_disk_aspect*obj2->prof_disk_aspect);
952 
    xm2 = lm0*lm0*dx2 + lm1*lm1*dy2 + lm0*lm1*dxy;
953 
    ym2 = lm2*lm2*dx2 + lm3*lm3*dy2 + lm2*lm3*dxy;
954 
    xym = lm0*lm2*dx2 + lm1*lm3*dy2 + (lm0*lm3+lm1*lm2)*dxy;
955 
    temp=xm2-ym2;
956 
    if (FLAG(obj2.prof_disk_thetaw))
957 
      {
958 
      obj2->prof_disk_thetaw = fmod_m90_p90((temp == 0.0)?
959 
                (45.0) : (0.5*atan2(2.0 * xym,temp)/DEG));
960 2 bertin 
 
961 173 bertin 
/*---- Compute position angles in J2000 or B1950 reference frame */
962 
      if (wcs->lng != wcs->lat)
963 
        {
964 
        if (FLAG(obj2.prof_disk_thetas))
965 
          obj2->prof_disk_thetas = fmod_m90_p90(lng<lat?
966 
                        ((obj2->prof_disk_thetaw>0.0?90:-90.0)
967 
                                - obj2->prof_disk_thetaw)
968 
                        : obj2->prof_disk_thetaw);
969 
        if (FLAG(obj2.prof_disk_theta2000))
970 
          obj2->prof_disk_theta2000 = fmod_m90_p90(obj2->prof_disk_thetas
971 
                                        + obj2->dtheta2000);
972 
        if (FLAG(obj2.prof_disk_theta1950))
973 
          obj2->prof_disk_theta1950 = fmod_m90_p90(obj2->prof_disk_thetas
974 
                                        + obj2->dtheta1950);
975 
        }
976 
      }
977 
    temp = sqrt(0.25*temp*temp+xym*xym);
978 
    pm2 = 0.5*(xm2+ym2);
979 
    obj2->prof_disk_scalew = sqrt(pm2+temp);
980 
    obj2->prof_disk_scaleerrw = obj2->prof_disk_scale > 0.0?
981 
        obj2->prof_disk_scaleerr/obj2->prof_disk_scale*obj2->prof_disk_scalew
982 
        : 0.0;
983 233 bertin 
    obj2->prof_disk_aspectw = (obj2->prof_disk_scalew>0.0 && pm2>temp)?
984 173 bertin 
                                  sqrt(pm2-temp) / obj2->prof_disk_scalew
985 
                                : obj2->prof_disk_aspect;
986 
    obj2->prof_disk_aspecterrw = obj2->prof_disk_aspect > 0.0?
987 
        obj2->prof_disk_aspecterr/obj2->prof_disk_aspect*obj2->prof_disk_aspectw
988 
        : 0.0;
989 
/*-- Arms World coordinates */
990 
    if (FLAG(obj2.prof_arms_scalew))
991 
      {
992 
      obj2->prof_arms_scalew = (obj2->prof_disk_scale > 0.0) ?
993 
          obj2->prof_arms_scale*obj2->prof_disk_scalew/obj2->prof_disk_scale
994 
        : 0.0;
995 
      obj2->prof_arms_scaleerrw = (obj2->prof_arms_scale > 0.0) ?
996 
        obj2->prof_arms_scaleerr/obj2->prof_arms_scale*obj2->prof_arms_scalew
997 
        : 0.0;
998 
      obj2->prof_arms_startw = (obj2->prof_disk_scale > 0.0) ?
999 
          obj2->prof_arms_start*obj2->prof_disk_scalew/obj2->prof_disk_scale
1000 
        : 0.0;
1001 
      obj2->prof_arms_starterrw = (obj2->prof_arms_start > 0.0) ?
1002 
        obj2->prof_arms_starterr/obj2->prof_arms_start*obj2->prof_arms_startw
1003 
        : 0.0;
1004 
      }
1005 
    }
1006 2 bertin 
 
1007 173 bertin 
/* Bar World coordinates */
1008 
  obj2->prof_bar_thetaerrw = obj2->prof_bar_thetaerr;
1009 
  if (FLAG(obj2.prof_bar_lengthw))
1010 2 bertin 
    {
1011 173 bertin 
    ct = cos(obj2->prof_bar_theta*DEG);
1012 
    st = sin(obj2->prof_bar_theta*DEG);
1013 
    dx2 = obj2->prof_bar_length*obj2->prof_bar_length * (ct*ct
1014 
        + st*st * obj2->prof_bar_aspect*obj2->prof_bar_aspect);
1015 
    dy2 = obj2->prof_bar_length*obj2->prof_bar_length * (st*st
1016 
        + ct*ct * obj2->prof_bar_aspect*obj2->prof_bar_aspect);
1017 
    dxy = ct*st * obj2->prof_bar_length*obj2->prof_bar_length
1018 
        *(1.0 - obj2->prof_bar_aspect*obj2->prof_bar_aspect);
1019 
    xm2 = lm0*lm0*dx2 + lm1*lm1*dy2 + lm0*lm1*dxy;
1020 
    ym2 = lm2*lm2*dx2 + lm3*lm3*dy2 + lm2*lm3*dxy;
1021 
    xym = lm0*lm2*dx2 + lm1*lm3*dy2 + (lm0*lm3+lm1*lm2)*dxy;
1022 
    temp=xm2-ym2;
1023 
    if (FLAG(obj2.prof_bar_thetaw))
1024 
      {
1025 
      obj2->prof_bar_thetaw = fmod_m90_p90((temp == 0.0)?
1026 
                (45.0) : (0.5*atan2(2.0 * xym,temp)/DEG));
1027 
 
1028 
/*---- Compute position angles in J2000 or B1950 reference frame */
1029 
      if (wcs->lng != wcs->lat)
1030 
        {
1031 
        if (FLAG(obj2.prof_bar_thetas))
1032 
          obj2->prof_bar_thetas = fmod_m90_p90(lng<lat?
1033 
                        ((obj2->prof_bar_thetaw>0.0?90:-90.0)
1034 
                                - obj2->prof_bar_thetaw)
1035 
                        : obj2->prof_bar_thetaw);
1036 
        if (FLAG(obj2.prof_bar_theta2000))
1037 
          obj2->prof_bar_theta2000 = fmod_m90_p90(obj2->prof_bar_thetas
1038 
                                        + obj2->dtheta2000);
1039 
        if (FLAG(obj2.prof_bar_theta1950))
1040 
          obj2->prof_bar_theta1950 = fmod_m90_p90(obj2->prof_bar_thetas
1041 
                                        + obj2->dtheta1950);
1042 
        }
1043 
      }
1044 
    temp = sqrt(0.25*temp*temp+xym*xym);
1045 
    pm2 = 0.5*(xm2+ym2);
1046 
    obj2->prof_bar_lengthw = sqrt(pm2+temp);
1047 
    obj2->prof_bar_lengtherrw = obj2->prof_bar_length > 0.0?
1048 
        obj2->prof_bar_lengtherr/obj2->prof_bar_length*obj2->prof_bar_lengthw
1049 
        : 0.0;
1050 
    obj2->prof_bar_aspectw = obj2->prof_bar_lengthw>0.0?
1051 
                                  sqrt(pm2-temp) / obj2->prof_bar_lengthw
1052 
                                : obj2->prof_bar_aspect;
1053 
    obj2->prof_bar_aspecterrw = obj2->prof_bar_aspect > 0.0?
1054 
        obj2->prof_bar_aspecterr/obj2->prof_bar_aspect*obj2->prof_bar_aspectw
1055 
        : 0.0;
1056 2 bertin 
    }
1057 
 
1058 206 bertin 
/* Global 2nd-order moments */
1059 
  if (FLAG(obj2.prof_mx2w))
1060 
    {
1061 
    dx2 = obj2->prof_mx2;
1062 
    dy2 = obj2->prof_my2;
1063 
    dxy = obj2->prof_mxy;
1064 226 bertin 
    mat[0] = lm0*lm0;
1065 
    mat[3] = lm1*lm1;
1066 
    mat[6] = lm0*lm1;
1067 
    mat[1] = lm2*lm2;
1068 
    mat[4] = lm3*lm3;
1069 
    mat[7] = lm2*lm3;
1070 
    mat[2] = lm0*lm2;
1071 
    mat[5] = lm1*lm3;
1072 
    mat[8] = lm0*lm3+lm1*lm2;
1073 
    obj2->prof_mx2w = xm2 = mat[0]*dx2 + mat[3]*dy2 + mat[6]*dxy;
1074 
    obj2->prof_my2w = ym2 = mat[1]*dx2 + mat[4]*dy2 + mat[7]*dxy;
1075 
    obj2->prof_mxyw = xym = mat[2]*dx2 + mat[5]*dy2 + mat[8]*dxy;
1076 206 bertin 
    temp=xm2-ym2;
1077 
    if (FLAG(obj2.prof_thetaw))
1078 
      {
1079 
      obj2->prof_thetaw = fmod_m90_p90((temp == 0.0)?
1080 
                        (45.0) : (0.5*atan2(2.0*xym,temp)/DEG));
1081 
 
1082 
/*---- Compute position angles in J2000 or B1950 reference frame */
1083 
      if (wcs->lng != wcs->lat)
1084 
        {
1085 
        if (FLAG(obj2.prof_thetas))
1086 
          obj2->prof_thetas = fmod_m90_p90(lng<lat?
1087 
                        ((obj2->prof_thetaw>0.0?90:-90.0) - obj2->prof_thetaw)
1088 
                        : obj2->prof_thetaw);
1089 
        if (FLAG(obj2.prof_theta2000))
1090 
          obj2->prof_theta2000 = fmod_m90_p90(obj2->prof_thetas
1091 
                + obj2->dtheta2000);
1092 
        if (FLAG(obj2.prof_theta1950))
1093 
          obj2->prof_theta1950 = fmod_m90_p90(obj2->prof_thetas
1094 
                + obj2->dtheta1950);
1095 
        }
1096 
      }
1097 
 
1098 
    if (FLAG(obj2.prof_aw))
1099 
      {
1100 
      temp = sqrt(0.25*temp*temp+xym*xym);
1101 
      pm2 = 0.5*(xm2+ym2);
1102 
      obj2->prof_aw = (float)sqrt(pm2+temp);
1103 
      obj2->prof_bw = (float)sqrt(pm2-temp);
1104 
      }
1105 
 
1106 
    if (FLAG(obj2.prof_cxxw))
1107 
      {
1108 
/*---- Handle large, fully correlated profiles (can cause a singularity...) */
1109 
      if ((temp=xm2*ym2-xym*xym)<1e-6)
1110 
        {
1111 
        temp = 1e-6;
1112 
        xym *= 0.99999;
1113 
        }
1114 
      obj2->prof_cxxw = (float)(ym2/temp);
1115 
      obj2->prof_cyyw = (float)(xm2/temp);
1116 
      obj2->prof_cxyw = (float)(-2*xym/temp);
1117 
      }
1118 
 
1119 226 bertin 
/*-- Use the Jacobians to compute the moment covariance matrix */
1120 
    if (FLAG(obj2.prof_pol1errw) || FLAG(obj2.prof_e1errw))
1121 
      propagate_covar(obj2->prof_mx2cov, mat, mx2wcov, 3, 3, tempmat);
1122 
 
1123 225 bertin 
    if (FLAG(obj2.prof_pol1w))
1124 206 bertin 
      {
1125 
      if (xm2+ym2 > 1.0/BIG)
1126 
        {
1127 
        obj2->prof_pol1w = (xm2 - ym2) / (xm2+ym2);
1128 
        obj2->prof_pol2w = 2.0*xym / (xm2 + ym2);
1129 226 bertin 
        if (FLAG(obj2.prof_pol1errw))
1130 
          {
1131 
/*-------- Compute the Jacobian of polarisation */
1132 
          invden = 1.0/(xm2+ym2);
1133 
          dpdmx2[0] =  2.0*ym2*invden*invden;
1134 
          dpdmx2[1] = -2.0*xm2*invden*invden;
1135 
          dpdmx2[2] =  0.0;
1136 
          dpdmx2[3] = -2.0*xym*invden*invden;
1137 
          dpdmx2[4] = -2.0*xym*invden*invden;
1138 
          dpdmx2[5] =  2.0*invden;
1139 
          propagate_covar(mx2wcov, dpdmx2, cov, 3, 2, tempmat);
1140 
          obj2->prof_pol1errw = (float)sqrt(cov[0]<0.0? 0.0: cov[0]);
1141 
          obj2->prof_pol2errw = (float)sqrt(cov[3]<0.0? 0.0: cov[3]);
1142 
          obj2->prof_pol12corrw = (dval=cov[0]*cov[3]) > 0.0?
1143 
                                        (float)(cov[1]/sqrt(dval)) : 0.0;
1144 
          }
1145 225 bertin 
        }
1146 
      else
1147 226 bertin 
        obj2->prof_pol1w = obj2->prof_pol2w = obj2->prof_pol1errw
1148 
                = obj2->prof_pol2errw = obj2->prof_pol12corrw = 0.0;
1149 225 bertin 
      }
1150 
 
1151 
    if (FLAG(obj2.prof_e1w))
1152 
      {
1153 
      if (xm2+ym2 > 1.0/BIG)
1154 
        {
1155 226 bertin 
        temp = xm2*ym2 - xym*xym;
1156 
        den = (temp>=0.0) ? xm2+ym2+2.0*sqrt(temp) : xm2+ym2;
1157 
        invden = 1.0/den;
1158 
        obj2->prof_e1w = (float)(invden*(xm2 - ym2));
1159 
        obj2->prof_e2w = (float)(2.0 * invden * xym);
1160 
        if (FLAG(obj2.prof_e1errw))
1161 
        {
1162 
/*------ Compute the Jacobian of ellipticity */
1163 
        invstemp = (temp>=0.0) ? 1.0/sqrt(temp) : 0.0;
1164 
        dpdmx2[0] = ( den - (1.0+ym2*invstemp)*(xm2-ym2))*invden*invden;
1165 
        dpdmx2[1] = (-den - (1.0+xm2*invstemp)*(xm2-ym2))*invden*invden;
1166 
        dpdmx2[2] = 2.0*xym*invstemp*(xm2-ym2)*invden*invden;
1167 
        dpdmx2[3] = -2.0*xym*(1.0+ym2*invstemp)*invden*invden;
1168 
        dpdmx2[4] = -2.0*xym*(1.0+xm2*invstemp)*invden*invden;
1169 
        dpdmx2[5] =  (2.0*den+4.0*xym*xym*invstemp)*invden*invden;
1170 
 
1171 
/*------ Use the Jacobian to compute the ellipticity covariance matrix */
1172 
        propagate_covar(mx2wcov, dpdmx2, cov, 3, 2, tempmat);
1173 
        obj2->prof_e1errw = (float)sqrt(cov[0]<0.0? 0.0: cov[0]);
1174 
        obj2->prof_e2errw = (float)sqrt(cov[3]<0.0? 0.0: cov[3]);
1175 
        obj2->prof_e12corrw = (dval=cov[0]*cov[3]) > 0.0?
1176 
                                        (float)(cov[1]/sqrt(dval)) : 0.0;
1177 206 bertin 
        }
1178 226 bertin 
        }
1179 206 bertin 
      else
1180 226 bertin 
        obj2->prof_e1w = obj2->prof_e2w = obj2->prof_e1errw
1181 
                = obj2->prof_e2errw = obj2->prof_e12corrw = 0.0;
1182 206 bertin 
      }
1183 
    }
1184 
 
1185 2 bertin 
  return;
1186 
  }
1187