1 /*
2 psf.c
3
4 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5 *
6 * Part of: SExtractor
7 *
8 * Authors: E.BERTIN (IAP)
9 * P.DELORME (LAOG)
10 *
11 * Contents: Fit the PSF to a detection.
12 *
13 * Last modify: 12/01/2006
14 *
15 *%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
16 */
17
18 #ifdef HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21
22 #include <math.h>
23 #include <stdio.h>
24 #include <stdlib.h>
25 #include <string.h>
26
27 #include "define.h"
28 #include "globals.h"
29 #include "prefs.h"
30 #include "fits/fitscat.h"
31 #include "check.h"
32 #include "filter.h"
33 #include "image.h"
34 #include "poly.h"
35 #include "psf.h"
36
37 /*------------------------------- variables ---------------------------------*/
38
39
40 extern keystruct objkey[];
41 extern objstruct outobj;
42
43 /********************************* psf_init **********************************/
44 /*
45 Allocate memory and stuff for the PSF-fitting.
46 */
psf_init(psfstruct * psf)47 void psf_init(psfstruct *psf)
48 {
49 QMALLOC(thepsfit, psfitstruct, 1);
50 QMALLOC(thepsfit->x, float, prefs.psf_npsfmax);
51 QMALLOC(thepsfit->y, float, prefs.psf_npsfmax);
52 QMALLOC(thepsfit->flux, float, prefs.psf_npsfmax);
53 QMALLOC(ppsfit, psfitstruct, 1); /*?*/
54 QMALLOC(ppsfit->x, float, prefs.psf_npsfmax);
55 QMALLOC(ppsfit->y, float, prefs.psf_npsfmax);
56 QMALLOC(ppsfit->flux, float, prefs.psf_npsfmax);
57
58 return;
59 }
60
61
62 /********************************* psf_end ***********************************/
63 /*
64 Free memory occupied by the PSF-fitting stuff.
65 */
psf_end(psfstruct * psf,psfitstruct * psfit)66 void psf_end(psfstruct *psf, psfitstruct *psfit)
67 {
68 int d, ndim;
69
70 if (psf->pc)
71 pc_end(psf->pc);
72
73 ndim = psf->poly->ndim;
74 for (d=0; d<ndim; d++)
75 free(psf->contextname[d]);
76 free(psf->context);
77 free(psf->contextname);
78 free(psf->contextoffset);
79 free(psf->contextscale);
80 free(psf->contexttyp);
81 poly_end(psf->poly);
82 free(psf->maskcomp);
83 free(psf->maskloc);
84 free(psf->masksize);
85 free(psf);
86
87 free(psfit->x);
88 free(psfit->y);
89 free(psfit->flux);
90 free(psfit);
91
92 return;
93 }
94
95
96 /********************************* psf_load *********************************/
97 /*
98 Read the PSF data from a FITS file.
99 */
psf_load(char * filename)100 psfstruct *psf_load(char *filename)
101 {
102 static objstruct saveobj;
103 static obj2struct saveobj2;
104 psfstruct *psf;
105 catstruct *cat;
106 tabstruct *tab;
107 keystruct *key;
108 char *head, *ci,*co;
109 int deg[POLY_MAXDIM], group[POLY_MAXDIM], ndim, ngroup,
110 i,k;
111
112 /* Open the cat (well it is not a "cat", but simply a FITS file */
113 if (!(cat = read_cat(filename)))
114 error(EXIT_FAILURE, "*Error*: PSF file not found: ", filename);
115
116 /* OK, we now allocate memory for the PSF structure itself */
117 QCALLOC(psf, psfstruct, 1);
118
119 /* Store a short copy of the PSF filename */
120 if ((ci=strrchr(filename, '/')))
121 strcpy(psf->name, ci+1);
122 else
123 strcpy(psf->name, filename);
124
125 if (!(tab = name_to_tab(cat, "PSF_DATA", 0)))
126 error(EXIT_FAILURE, "*Error*: PSF_DATA table not found in catalog ",
127 filename);
128
129 head = tab->headbuf;
130
131 /*-- Dimension of the polynomial */
132 if (fitsread(head, "POLNAXIS", &ndim, H_INT,T_LONG) == RETURN_OK
133 && ndim)
134 {
135 /*-- So we have a polynomial description of the PSF variations */
136 if (ndim > POLY_MAXDIM)
137 {
138 sprintf(gstr, "*Error*: The POLNAXIS parameter in %s exceeds %d",
139 psf->name, POLY_MAXDIM);
140 error(EXIT_FAILURE, gstr, "");
141 }
142
143 QMALLOC(psf->contextname, char *, ndim);
144 QMALLOC(psf->context, double *, ndim);
145 QMALLOC(psf->contexttyp, t_type, ndim);
146 QMALLOC(psf->contextoffset, double, ndim);
147 QMALLOC(psf->contextscale, double, ndim);
148
149 /*-- We will have to use the outobj structs, so we first save their content */
150 saveobj = outobj;
151 saveobj2 = outobj2;
152 /*-- outobj's are used as FLAG arrays, so we initialize them to 0 */
153 memset(&outobj, 0, sizeof(outobj));
154 memset(&outobj2, 0, sizeof(outobj2));
155 for (i=0; i<ndim; i++)
156 {
157 /*---- Polynomial groups */
158 sprintf(gstr, "POLGRP%1d", i+1);
159 if (fitsread(head, gstr, &group[i], H_INT,T_LONG) != RETURN_OK)
160 goto headerror;
161
162 /*---- Contexts */
163 QMALLOC(psf->contextname[i], char, 80);
164 sprintf(gstr, "POLNAME%1d", i+1);
165 if (fitsread(head,gstr,psf->contextname[i],H_STRING,T_STRING)!=RETURN_OK)
166 goto headerror;
167 if (*psf->contextname[i]==(char)':')
168 /*------ It seems we're facing a FITS header parameter */
169 psf->context[i] = NULL; /* This is to tell we'll have to load */
170 /* a FITS header context later on */
171 else
172 /*------ The context element is a dynamic object parameter */
173 {
174 if ((k = findkey(psf->contextname[i], (char *)objkey,
175 sizeof(keystruct)))==RETURN_ERROR)
176 {
177 sprintf(gstr, "*Error*: %s CONTEXT parameter in %s unknown",
178 psf->contextname[i], psf->name);
179 error(EXIT_FAILURE, gstr, "");
180 }
181 key = objkey+k;
182 psf->context[i] = key->ptr;
183 psf->contexttyp[i] = key->ttype;
184 /*------ Declare the parameter "active" to trigger computation by SExtractor */
185 *((char *)key->ptr) = (char)'\1';
186 }
187 /*---- Scaling of the context parameter */
188 sprintf(gstr, "POLZERO%1d", i+1);
189 if (fitsread(head, gstr, &psf->contextoffset[i], H_EXPO, T_DOUBLE)
190 !=RETURN_OK)
191 goto headerror;
192 sprintf(gstr, "POLSCAL%1d", i+1);
193 if (fitsread(head, gstr, &psf->contextscale[i], H_EXPO, T_DOUBLE)
194 !=RETURN_OK)
195 goto headerror;
196 }
197
198 /*-- Number of groups */
199 if (fitsread(head, "POLNGRP ", &ngroup, H_INT, T_LONG) != RETURN_OK)
200 goto headerror;
201
202 for (i=0; i<ngroup; i++)
203 {
204 /*---- Polynomial degree for each group */
205 sprintf(gstr, "POLDEG%1d", i+1);
206 if (fitsread(head, gstr, °[i], H_INT,T_LONG) != RETURN_OK)
207 goto headerror;
208 }
209
210 psf->poly = poly_init(group, ndim, deg, ngroup);
211
212 /*-- Update the permanent FLAG arrays (that is, perform an "OR" on them) */
213 for (ci=(char *)&outobj,co=(char *)&flagobj,i=sizeof(objstruct); i--;)
214 *(co++) |= *(ci++);
215 for (ci=(char *)&outobj2,co=(char *)&flagobj2,i=sizeof(obj2struct); i--;)
216 *(co++) |= *(ci++);
217
218 /*-- Restore previous outobj contents */
219 outobj = saveobj;
220 outobj2 = saveobj2;
221 }
222 else
223 {
224 /*-- This is a simple, constant PSF */
225 psf->poly = poly_init(group, 0, deg, 0);
226 psf->context = NULL;
227 }
228
229 /* Dimensionality of the PSF mask */
230 if (fitsread(head, "PSFNAXIS", &psf->maskdim, H_INT, T_LONG) != RETURN_OK)
231 goto headerror;
232 if (psf->maskdim<2 || psf->maskdim>3)
233 error(EXIT_FAILURE, "*Error*: wrong dimensionality for the PSF "
234 "mask in ", filename);
235 QMALLOC(psf->masksize, int, psf->maskdim);
236 for (i=0; i<psf->maskdim; i++)
237 psf->masksize[i] = 1;
238 psf->masknpix = 1;
239 for (i=0; i<psf->maskdim; i++)
240 {
241 sprintf(gstr, "PSFAXIS%1d", i+1);
242 if (fitsread(head, gstr, &psf->masksize[i], H_INT,T_LONG) != RETURN_OK)
243 goto headerror;
244 psf->masknpix *= psf->masksize[i];
245 }
246
247 /* PSF FWHM: defaulted to 3 pixels */
248 if (fitsread(head, "PSF_FWHM", &psf->fwhm, H_FLOAT,T_DOUBLE) != RETURN_OK)
249 psf->fwhm = 3.0;
250
251 /* PSF oversampling: defaulted to 1 */
252 if (fitsread(head, "PSF_SAMP", &psf->pixstep,H_FLOAT,T_FLOAT) != RETURN_OK)
253 psf->pixstep = 1.0;
254
255 /* Load the PSF mask data */
256 key = read_key(tab, "PSF_MASK");
257 psf->maskcomp = key->ptr;
258
259 psf->pc = pc_load(cat);
260
261 QMALLOC(psf->maskloc, double, psf->masksize[0]*psf->masksize[1]);
262
263 /* But don't touch my arrays!! */
264 blank_keys(tab);
265
266 free_cat(&cat, 1);
267
268 return psf;
269
270 headerror:
271 error(EXIT_FAILURE, "*Error*: Incorrect or obsolete PSF data in ", filename);
272 return NULL;
273 }
274
275
276 /***************************** psf_readcontext *******************************/
277 /*
278 Read the PSF context parameters in the FITS header.
279 */
psf_readcontext(psfstruct * psf,picstruct * field)280 void psf_readcontext(psfstruct *psf, picstruct *field)
281 {
282 static double contextval[POLY_MAXDIM];
283 int i, ndim;
284
285 ndim = psf->poly->ndim;
286 for (i=0; i<ndim; i++)
287 if (!psf->context[i])
288 {
289 psf->context[i] = &contextval[i];
290 psf->contexttyp[i] = T_DOUBLE;
291 if (fitsread(field->fitshead, psf->contextname[i]+1, &contextval[i],
292 H_FLOAT,T_DOUBLE) == RETURN_ERROR)
293 {
294 sprintf(gstr, "*Error*: %s parameter not found in the header of ",
295 psf->contextname[i]+1);
296 error(EXIT_FAILURE, gstr, field->rfilename);
297 }
298 }
299
300 return;
301 }
302
303
304 /******************************** psf_fit ***********************************/
305 /* standart PSF fit for one component */
306 /****************************************************************************/
307
psf_fit(psfstruct * psf,picstruct * field,picstruct * wfield,objstruct * obj)308 void psf_fit(psfstruct *psf, picstruct *field, picstruct *wfield,
309 objstruct *obj)
310 {
311 checkstruct *check;
312 static obj2struct *obj2 = &outobj2;
313 static double x2[PSF_NPSFMAX],y2[PSF_NPSFMAX],xy[PSF_NPSFMAX],
314 deltax[PSF_NPSFMAX],
315 deltay[PSF_NPSFMAX],flux[PSF_NPSFMAX],
316 deltaxb[PSF_NPSFMAX],deltayb[PSF_NPSFMAX],
317 fluxb[PSF_NPSFMAX],
318 sol[PSF_NTOT], covmat[PSF_NTOT*PSF_NTOT],
319 vmat[PSF_NTOT*PSF_NTOT], wmat[PSF_NTOT];
320 double **psfmasks, **psfmaskx,**psfmasky,
321 *data, *data2, *data3, *weight, *mat, *checkdata,
322 *d, *m, *w, *ps, *var,
323 dx,dy,
324 pix,pix2, wthresh,val,
325 backnoise2, gain, radmin2,radmax2,satlevel,chi2,
326 r2, valmax, psf_fwhm;
327 float *dh, *wh, pixstep,fluxerr;
328 PIXTYPE *datah, *weighth, *cpix;
329 int i,j,p, npsf,npsfmax, npix, nppix, ix,iy,niter,
330 width, height, pwidth,pheight, x,y,
331 xmax,ymax, wbad, gainflag, convflag, npsfflag,
332 ival,kill=0;
333
334 checkdata = NULL; /* To avoid gcc -Wall warnings */
335 dx = dy = 0.0;
336 niter = 0;
337 npsfmax = prefs.psf_npsfmax;
338 pixstep = 1.0/psf->pixstep;
339 gain = prefs.gain;
340 backnoise2 = field->backsig*field->backsig;
341 satlevel = prefs.satur_level - obj->bkg;
342 wthresh = wfield?wfield->weight_thresh:BIG;
343 gainflag = prefs.weightgain_flag;
344 psf_fwhm = psf->fwhm*psf->pixstep;
345
346
347 /* Initialize outputs */
348 thepsfit->niter = 0;
349 thepsfit->npsf = 0;
350 for (j=0; j<npsfmax; j++)
351 {
352 thepsfit->x[j] = obj2->posx;
353 thepsfit->y[j] = obj2->posy;
354 thepsfit->flux[j] = 0.0;
355 }
356
357 /* Scale data area with object "size" */
358 ix = (obj->xmax+obj->xmin+1)/2;
359 iy = (obj->ymax+obj->ymin+1)/2;
360 width = obj->xmax-obj->xmin+1+psf_fwhm;
361 if (width < (ival=(int)(psf_fwhm*2)))
362 width = ival;
363 height = obj->ymax-obj->ymin+1+psf_fwhm;
364 if (height < (ival=(int)(psf_fwhm*2)))
365 height = ival;
366 npix = width*height;
367 radmin2 = PSF_MINSHIFT*PSF_MINSHIFT;
368 radmax2 = npix/2.0;
369 fluxerr = 0.0;
370
371 /* Scale total area with PSF FWHM */
372 pwidth = (int)(psf->masksize[0]*psf->pixstep)+width;;
373 pheight = (int)(psf->masksize[1]*psf->pixstep)+height;
374 nppix = pwidth*pheight;
375
376 /* Allocate working space */
377 if (prefs.psf_flag==1)
378 if (prefs.dpsf_flag!=1)
379 if(!FLAG(obj2.fluxerr_psf))
380 QMALLOC(obj2->fluxerr_psf, float, prefs.psf_npsfmax);
381
382 QMALLOC(weighth, PIXTYPE, npix);
383 QMALLOC(weight, double, npix);
384 QMALLOC(datah, PIXTYPE, npix);
385 QMALLOC(data, double, npix);
386 QMALLOC(data2, double, npix);
387 QMALLOC(data3, double, npix);
388 QMALLOC(mat, double, npix*PSF_NTOT);
389 if (prefs.check[CHECK_SUBPSFPROTOS] || prefs.check[CHECK_PSFPROTOS]
390 || prefs.check[CHECK_SUBPCPROTOS] || prefs.check[CHECK_PCPROTOS]
391 || prefs.check[CHECK_PCOPROTOS])
392 {
393 QMALLOC(checkdata, double, nppix);
394 QMALLOC(checkmask, PIXTYPE, nppix);
395 }
396
397 QMALLOC(psfmasks, double *, npsfmax);
398 QMALLOC(psfmaskx, double *, npsfmax);
399 QMALLOC(psfmasky, double *, npsfmax);
400 for (i=0; i<npsfmax; i++)
401 {
402 QMALLOC(psfmasks[i], double, npix);
403 QMALLOC(psfmaskx[i], double, npix);
404 QMALLOC(psfmasky[i], double, npix);
405 }
406
407 copyimage(field, datah, width, height, ix, iy);
408
409 /* Compute weights */
410 wbad = 0;
411 if (wfield)
412 {
413 copyimage(wfield, weighth, width, height, ix, iy);
414 for (wh=weighth, w=weight, dh=datah,p=npix; p--;)
415 if ((pix=*(wh++)) < wthresh && pix>0
416 && (pix2=*(dh++))>-BIG
417 && pix2<satlevel)
418 *(w++) = 1/sqrt(pix+(pix2>0.0?
419 (gainflag? pix2*pix/backnoise2:pix2)/gain
420 :0.0));
421 else
422 {
423 *(w++) = 0.0;
424 wbad++;
425 }
426 }
427 else
428 for (w=weight, dh=datah, p=npix; p--;)
429 if ((pix=*(dh++))>-BIG && pix<satlevel)
430 *(w++) = 1.0/sqrt(backnoise2+(pix>0.0?pix/gain:0.0));
431 else
432 {
433 *(w++) = 0.0;
434 wbad++;
435 }
436
437 /* Special action if most of the weights are zero!! */
438 if (wbad>=npix-3)
439 return;
440
441 /* Weight the data */
442 dh = datah;
443 val = obj->dbkg; /* Take into account a local background change */
444 d = data;
445 w = weight;
446 for (p=npix; p--;)
447 *(d++) = (*(dh++)-val)**(w++);
448
449 /* Get the local PSF */
450 psf_build(psf);
451
452 npsfflag = 1;
453 r2 = psf_fwhm*psf_fwhm/2.0;
454 fluxb[0] = deltaxb[0] = deltayb[0] = 0.0;
455
456 for (npsf=1; npsf<=npsfmax && npsfflag; npsf++)
457 {
458 kill=0;
459 /*-- First compute an optimum initial guess for the positions of components */
460 if (npsf>1)
461 {
462 /*---- Subtract previously fitted components */
463 d = data2;
464 dh = datah;
465 for (p=npix; p--;)
466 *(d++) = (double)*(dh++);
467 for (j=0; j<npsf-1; j++)
468 {
469 d = data2;
470 ps = psfmasks[j];
471 for (p=npix; p--;)
472 *(d++) -= flux[j]**(ps++);
473 }
474 convolve_image(field, data2, data3, width,height);
475 /*---- Ignore regions too close to stellar cores */
476 for (j=0; j<npsf-1; j++)
477 {
478 d = data3;
479 dy = -((double)(height/2)+deltay[j]);
480 for (y=height; y--; dy += 1.0)
481 {
482 dx = -((double)(width/2)+deltax[j]);
483 for (x=width; x--; dx+= 1.0, d++)
484 if (dx*dx+dy*dy<r2)
485 *d = -BIG;
486 }
487 }
488 /*---- Now find the brightest pixel (poor man's guess, to be refined later) */
489 d = data3;
490 valmax = -BIG;
491 xmax = width/2;
492 ymax = height/2;
493 for (y=0; y<height; y++)
494 for (x=0; x<width; x++)
495 {
496 if ((val = *(d++))>valmax)
497 {
498 valmax = val;
499 xmax = x;
500 ymax = y;
501 }
502 }
503 deltax[npsf-1] = (double)(xmax - width/2);
504 deltay[npsf-1] = (double)(ymax - height/2);
505 }
506 else
507 {
508 /*---- Only one component to fit: simply use the barycenter as a guess */
509 deltax[npsf-1] = obj->mx - ix;
510 deltay[npsf-1] = obj->my - iy;
511 }
512
513 niter = 0;
514 convflag = 1;
515 for (i=0; i<PSF_NITER && convflag; i++)
516 {
517 convflag = 0,niter++,m=mat;
518 for (j=0; j<npsf; j++)
519 {
520 /*------ Resample the PSFs here for the 1st iteration */
521 vignet_resample(psf->maskloc, psf->masksize[0], psf->masksize[1],
522 psfmasks[j], width, height,
523 -deltax[j]*pixstep, -deltay[j]*pixstep,
524 pixstep);
525 m=compute_gradient(weight,width,height,
526 psfmasks[j],psfmaskx[j],psfmasky[j],m);
527 }
528
529
530 svdfit(mat, data, npix, npsf*PSF_NA, sol, vmat, wmat);
531
532 compute_pos( &npsf, &convflag, &npsfflag,radmin2,radmax2,
533 r2, sol,flux, deltax, deltay,&dx,&dy);
534 }
535 for (j=0; j<npsf; j++)
536 {
537 /*-- Compute variances and covariances */
538 svdvar(vmat, wmat, npsf*PSF_NA, covmat);
539 var = covmat;
540 /*---- First, the error on the flux estimate */
541 fluxerr = sqrt(*var)>0.0? sqrt(*var):999999.0;
542 //if (flux[j]<12*fluxerr && j>0)
543 // npsfmax--,flux[j]=0;
544 if (flux[j]<12*fluxerr && j>0)
545 {
546 flux[j]=0,kill++,npsfmax--;
547 //if(j==npsfmax-1)
548 // kill++;
549 }
550 }
551 if (npsfflag)
552 {
553 /*--- If we reach this point we know the data are worth backuping */
554 for (j=0; j<npsf; j++)
555 {
556 deltaxb[j] = deltax[j];
557 deltayb[j] = deltay[j];
558 fluxb[j] = flux[j];
559 obj2->fluxerr_psf[j]=fluxerr;
560 }
561 }
562 }
563 npsf=npsf-1-kill;
564
565 /* Now keep only fitted stars that fall within the current detection area */
566 i = 0;
567 for (j=0; j<npsf; j++)
568 {
569 x = (int)(deltaxb[j]+0.4999)+width/2;
570 y = (int)(deltayb[j]+0.4999)+height/2;
571 if (x<0 || x>=width || y<0 || y>=height)
572 continue;
573 if (weight[y*width+x] < 1/BIG)
574 continue;
575 if (10*fluxb[j]<fluxb[0] )
576 continue;
577 if (fluxb[j]<=0 )
578 continue;
579
580 if (FLAG(obj2.poserrmx2_psf))
581 {
582 compute_poserr(j,var,sol,obj2,x2,y2,xy);
583 }
584 else
585 var += 3*PSF_NA+3;
586
587 deltax[i] = deltaxb[j];
588 deltay[i] = deltayb[j];
589 flux[i++] = fluxb[j];
590 }
591
592 npsf = i;
593
594 /* Compute chi2 if asked to
595 if (FLAG(obj2.chi2_psf))
596 {
597 for (j=0; j<npsf; j++)
598 {
599 chi2 = 0.0;
600 for (d=data,w=weight,p=0; p<npix; w++,p++)
601 {
602 pix = *(d++);
603 pix -= psfmasks[j][p]*flux[j]**w;
604 chi2 += pix*pix;
605 if (chi2>1E29) chi2=1E28;
606 }
607 obj2->chi2_psf = obj->sigbkg>0.?
608 chi2/((npix - 3*npsf)*obj->sigbkg*obj->sigbkg):999999;
609
610 }
611
612 }*/
613 /* Compute relative chi2 if asked to */
614 if (FLAG(obj2.chi2_psf))
615 {
616 for (j=0; j<npsf; j++)
617 {
618 chi2 = 0.0;
619 for (d=data,w=weight,p=0; p<npix; w++,p++)
620 {
621 pix = *(d++)/flux[j];
622 pix -= psfmasks[j][p]**w;
623 chi2 += pix*pix;
624 if (chi2>1E29) chi2=1E28;
625 }
626 obj2->chi2_psf = flux[j]>0?
627 chi2/((npix - 3*npsf)*obj->sigbkg*obj->sigbkg):999999;
628
629 }
630
631 }
632 /* CHECK images */
633 if (prefs.check[CHECK_SUBPSFPROTOS] || prefs.check[CHECK_PSFPROTOS])
634 for (j=0; j<npsf; j++)
635 {
636 vignet_resample(psf->maskloc, psf->masksize[0], psf->masksize[1],
637 checkdata, pwidth, pheight,
638 -deltax[j]*pixstep, -deltay[j]*pixstep, pixstep);
639 cpix = checkmask;
640 d = checkdata;
641 for (p=nppix; p--;)
642 *(cpix++) = (PIXTYPE)*(d++);
643 if ((check = prefs.check[CHECK_SUBPSFPROTOS]))
644 addcheck(check, checkmask, pwidth,pheight, ix,iy,-flux[j]);
645 if ((check = prefs.check[CHECK_PSFPROTOS]))
646 addcheck(check, checkmask, pwidth,pheight, ix,iy,flux[j]);
647 }
648
649 thepsfit->niter = niter;
650 thepsfit->npsf = npsf;
651 for (j=0; j<npsf; j++)
652 {
653 thepsfit->x[j] = ix+deltax[j]+1.0;
654 thepsfit->y[j] = iy+deltay[j]+1.0;
655 thepsfit->flux[j] = flux[j];
656 }
657
658
659
660
661 /* Now the morphology stuff */
662 if (prefs.pc_flag)
663 {
664 width = pwidth-1;
665 height = pheight-1;
666 npix = width*height;
667 copyimage(field, datah, width, height, ix, iy);
668
669 /*-- Re-compute weights */
670 if (wfield)
671 {
672 copyimage(wfield, weighth, width, height, ix, iy);
673 for (wh=weighth ,w=weight, p=npix; p--;)
674 *(w++) = (pix=*(wh++))<wthresh? sqrt(pix): 0.0;
675 }
676 else
677 for (w=weight, dh=datah, p=npix; p--;)
678 *(w++) = ((pix = *(dh++))>-BIG && pix<satlevel)?
679 1.0/sqrt(backnoise2+(pix>0.0?pix/gain:0.0))
680 :0.0;
681
682 /*-- Weight the data */
683 dh = datah;
684 d = data;
685 w = weight;
686 for (p=npix; p--;)
687 *(d++) = *(dh++)*(*(w++));
688
689 pc_fit(psf, data, weight, width, height, ix,iy, dx,dy, npix,
690 field->backsig);
691 }
692
693
694 for (i=0; i<prefs.psf_npsfmax; i++)
695 {
696 QFREE(psfmasks[i]);
697 QFREE(psfmaskx[i]);
698 QFREE(psfmasky[i]);
699 }
700
701 QFREE(psfmasks);
702 QFREE(psfmaskx);
703 QFREE(psfmasky);
704 QFREE(datah);
705 QFREE(data);
706 QFREE(data2);
707 QFREE(data3);
708 QFREE(weighth);
709 QFREE(weight);
710 QFREE(data);
711 QFREE(mat);
712
713 if (prefs.check[CHECK_SUBPSFPROTOS] || prefs.check[CHECK_PSFPROTOS]
714 || prefs.check[CHECK_SUBPCPROTOS] || prefs.check[CHECK_PCPROTOS]
715 || prefs.check[CHECK_PCOPROTOS])
716 {
717 QFREE(checkdata);
718 QFREE(checkmask);
719 }
720
721 return;
722 }
723
724
725 /******************************** double_psf_fit *******************************
726 ****/
727 /* double fit to make the psf detection on one image and the photometry on anoth
728 er */
729 /*******************************************************************************
730 ****/
731
double_psf_fit(psfstruct * ppsf,picstruct * pfield,picstruct * pwfield,objstruct * obj,psfstruct * psf,picstruct * field,picstruct * wfield)732 void double_psf_fit(psfstruct *ppsf, picstruct *pfield, picstruct *pwfield,
733 objstruct *obj, psfstruct *psf, picstruct *field,
734 picstruct *wfield)
735 {
736 static double /* sum[PSF_NPSFMAX]*/ pdeltax[PSF_NPSFMAX],
737 pdeltay[PSF_NPSFMAX],psol[PSF_NPSFMAX], pcovmat[PSF_NPSFMAX*PSF_NPSFMAX],
738 pvmat[PSF_NPSFMAX*PSF_NPSFMAX], pwmat[PSF_NPSFMAX],pflux[PSF_NPSFMAX];
739
740 double **ppsfmasks, **ppsfmaskx,**ppsfmasky,
741 *pmat, *checkdata,
742 *pdata, *pdata2, *pdata3,
743 *pm,*pd, *pw, *pps,*pweight,/* *pps, *px, *py,*/
744 dx,dy,pdx,pdy, /* x1,y1, mx,my,mflux, */
745 val, ppix,ppix2, /* dflux, */
746 gain, radmin2,radmax2,satlevel
747 ,chi2,pwthresh,pbacknoise2, /* mr, */
748 r2=0, psf_fwhm,ppsf_fwhm ;
749 float *pdh, *pwh, pixstep,ppixstep;
750 PIXTYPE *pdatah, *pweighth;
751 int i,j,k,p, npsf, npix,ix,iy,
752 width, height, /* hw,hh, */
753 x,y, /* yb, */
754 wbad, gainflag,
755 ival,npsfmax;
756 double *pvar;
757
758 static obj2struct *obj2 = &outobj2;
759 checkdata = NULL; /* To avoid gcc -Wall warnings */
760 pdx = pdy =dx = dy = 0.0;
761 ppixstep = 1.0/ppsf->pixstep;
762 pixstep = 1.0/psf->pixstep;
763 gain = prefs.gain;
764 npsfmax=prefs.psf_npsfmax;
765 pbacknoise2 = pfield->backsig*pfield->backsig;
766 satlevel = prefs.satur_level - obj->bkg;
767 gainflag = prefs.weightgain_flag;
768 psf_fwhm = psf->fwhm*psf->pixstep;
769 ppsf_fwhm = ppsf->fwhm*ppsf->pixstep;
770 pwthresh = pwfield?pwfield->weight_thresh:BIG;
771
772 /* Initialize outputs */
773 ppsfit->niter = 0;
774 ppsfit->npsf = 0;
775 if(!FLAG(obj2.fluxerr_psf))
776 QMALLOC(obj2->fluxerr_psf, float, npsfmax);
777 for (j=0; j<npsfmax; j++)
778 {
779 ppsfit->x[j] = 999999.0;
780 ppsfit->y[j] = 999999.0;
781 ppsfit->flux[j] = 0.0;
782 obj2->fluxerr_psf[j] = 0.0;
783 pdeltax[j]= pdeltay[j]=psol[j]= pwmat[j]=pflux[j]=0.0;
784
785 }
786
787 ix = (obj->xmax+obj->xmin+1)/2;
788 iy = (obj->ymax+obj->ymin+1)/2;
789 width = obj->xmax-obj->xmin+1+psf_fwhm;
790 if (width < (ival=(int)(psf_fwhm*2)))
791 width = ival;
792 height = obj->ymax-obj->ymin+1+psf_fwhm;
793 if (height < (ival=(int)(psf_fwhm*2)))
794 height = ival;
795 npix = width*height;
796 radmin2 = PSF_MINSHIFT*PSF_MINSHIFT;
797 radmax2 = npix/2.0;
798 psf_fit(psf,field, wfield,obj);
799 npsf=thepsfit->npsf;
800
801 QMALLOC(ppsfmasks,double *,npsfmax);
802 QMALLOC(ppsfmaskx,double *,npsfmax);
803 QMALLOC(ppsfmasky,double *,npsfmax);
804
805 for (i=0; i<npsfmax; i++)
806 {
807 QMALLOC(ppsfmasks[i],double,npix);
808 QMALLOC(ppsfmaskx[i],double,npix);
809 QMALLOC(ppsfmasky[i],double,npix);
810 }
811
812 QMALLOC(pweighth, PIXTYPE, npix);
813 QMALLOC(pweight, double, npix);
814 QMALLOC(pdatah, PIXTYPE, npix);
815 QMALLOC(pdata, double, npix);
816 QMALLOC(pdata2, double, npix);
817 QMALLOC(pdata3, double, npix);
818 QMALLOC(pmat, double, npix*npsfmax);
819
820 for (j=0; j<npsf; j++)
821 {
822 pdeltax[j] =thepsfit->x[j]-ix-1 ;
823 pdeltay[j] =thepsfit->y[j]-iy-1 ;
824 ppsfit->flux[j] = 0;
825 }
826
827 /*------------------- Now the photometry fit ---------------------*/
828 copyimage(pfield, pdatah, width, height, ix, iy);
829 /* Compute photometry weights */
830 wbad = 0;
831 if (pwfield)
832 {
833 copyimage(pwfield, pweighth, width, height, ix, iy);
834 for (pwh=pweighth, pw=pweight, pdh=pdatah,p=npix; p--;)
835 {
836 if ((ppix=*(pwh++)) < pwthresh && ppix>0
837 && (ppix2=*(pdh++))>-BIG && ppix2<satlevel)
838 {
839 *(pw++) = 1/sqrt(ppix+(ppix2>0.0?
840 (gainflag? ppix2*ppix/pbacknoise2:ppix2)/gain : 0.0));
841 }
842 else
843 {
844 *(pw++) = 0.0;
845 wbad++;
846 }
847 }
848 }
849 else
850 for (pw=pweight, pdh=pdatah, p=npix; p--;)
851 if ((ppix=*(pdh++))>-BIG && ppix<satlevel)
852 {
853 *(pw++) = 1.0/sqrt(pbacknoise2+(ppix>0.0?ppix/gain:0.0));
854 }
855 else
856 {
857 *(pw++) = 0.0;
858 wbad++;
859 }
860 /* Special action if most of the weights are zero!! */
861 if (wbad>=npix-3)
862 return;
863
864 /* Weight the data */
865 pdh = pdatah;
866 pd = pdata;
867 pw = pweight;
868 val = obj->dbkg;
869 for (p=npix; p--;)
870 *(pd++) = (*(pdh++)-val)**(pw++);
871
872
873 /* Get the photmetry PSF */
874 psf_build(ppsf);
875 for (j=1; j<=npsf; j++)
876 {
877 if (j>1)
878 {
879 /*---- Subtract //previously fitted components in photometry image */
880 pd = pdata2;
881 pdh = pdatah;
882 for (p=npix; p--;)
883 *(pd++) = (double)*(pdh++);
884 for (k=0; k<j-1; k++)
885 {
886 pd = pdata2;
887 pps = ppsfmasks[k];
888 for (p=npix; p--;)
889 *(pd++) -= pflux[k]**(pps++);
890 }
891 convolve_image(pfield, pdata2, pdata3, width,height);
892 /*---- Ignore regions too close to stellar cores */
893 for (k=0; k<j-1; k++)
894 {
895 pd = pdata3;
896 dy = -((double)(height/2)+pdeltay[k]);
897 for (y=height; y--; dy += 1.0)
898 {
899 dx = -((double)(width/2)+pdeltax[k]);
900 for (x=width; x--; dx+= 1.0, pd++)
901 if (dx*dx+dy*dy<r2) /*?*/
902 *pd = -BIG;
903 }
904 }
905 }
906
907 pm=pmat;
908 for (k=0; k<j; k++)
909 {
910 /*------ Resample the PSFs here for the 1st iteration */
911 vignet_resample(ppsf->maskloc,
912 ppsf->masksize[0], ppsf->masksize[1],
913 ppsfmasks[k], width, height,
914 -pdeltax[k]*ppixstep, -pdeltay[k]*ppixstep,
915 ppixstep);
916 pm=compute_gradient_phot(pweight,width,height, ppsfmasks[k],pm);
917 }
918
919 svdfit(pmat, pdata, npix, j, psol, pvmat, pwmat);
920 compute_pos_phot( &j, psol,pflux);
921
922 for (k=0; k<j; k++)
923 {
924 svdvar(pvmat, pwmat, j, pcovmat);
925 pvar = pcovmat;
926 obj2->fluxerr_psf[k]= sqrt(*pvar)>0.0 && sqrt(*pvar)<99?
927 sqrt(*pvar):99;
928 }
929 }
930 /* Compute chi2 if asked to
931 if (FLAG(obj2.chi2_psf))
932 {
933 for (j=0; j<npsf; j++)
934 {
935 chi2 = 0.0;
936 for (pd=pdata,pw=pweight,p=0; p<npix; pw++,p++)
937 {
938 ppix = *(pd++);
939 ppix -= ppsfmasks[j][p]*pflux[j]**pw;
940 chi2 += ppix*ppix;
941 if (chi2>1E29) chi2=1E28;
942 }
943 obj2->chi2_psf = obj->sigbkg>0.?
944 chi2/((npix - 3*npsf)*obj->sigbkg*obj->sigbkg):999999;
945
946 }
947
948 }
949 */
950 /* Compute relative error if asked to */
951 if (FLAG(obj2.chi2_psf))
952 {
953 for (j=0; j<npsf; j++)
954 {
955 chi2 = 0.0;
956 for (pd=pdata,pw=pweight,p=0; p<npix; pw++,p++)
957 {
958 ppix = *(pd++)/pflux[j];
959 ppix -= ppsfmasks[j][p]**pw;
960 chi2 += ppix*ppix;
961 if (chi2>1E29) chi2=1E28;
962 }
963 obj2->chi2_psf = pflux[j]>0?
964 chi2/((npix - 3*npsf)*obj->sigbkg*obj->sigbkg):999999;
965
966 }
967
968 }
969 ppsfit->niter = thepsfit->niter;
970 ppsfit->npsf = npsf;
971
972 for (j=0; j<npsf; j++)
973 {
974 thepsfit->x[j] = ix+pdeltax[j]+1.0;
975 thepsfit->y[j] = iy+pdeltay[j]+1.0;
976 thepsfit->flux[j] = pflux[j];
977 ppsfit->x[j] = ix+pdeltax[j]+1.0;
978 ppsfit->y[j] = iy+pdeltay[j]+1.0;
979 ppsfit->flux[j] = pflux[j];
980 }
981
982
983 for (i=0; i<npsfmax; i++)
984 {
985 QFREE(ppsfmasks[i]);
986 QFREE(ppsfmaskx[i]);
987 QFREE(ppsfmasky[i]);
988 }
989
990
991 QFREE(ppsfmasks);
992 QFREE(ppsfmaskx);
993 QFREE(ppsfmasky);
994 QFREE(pdatah);
995 QFREE(pdata);
996 QFREE(pdata2);
997 QFREE(pdata3);
998 QFREE(pweighth);
999 QFREE(pweight);
1000 QFREE(pdata);
1001 QFREE(pmat);
1002
1003 if (prefs.check[CHECK_SUBPSFPROTOS] || prefs.check[CHECK_PSFPROTOS]
1004 || prefs.check[CHECK_SUBPCPROTOS] || prefs.check[CHECK_PCPROTOS]
1005 || prefs.check[CHECK_PCOPROTOS])
1006 {
1007 QFREE(checkdata);
1008 QFREE(checkmask);
1009 }
1010 return;
1011 }
1012
1013 /******************************* psf_build **********************************/
1014 /*
1015 Build the local PSF (function of "context").
1016 */
psf_build(psfstruct * psf)1017 void psf_build(psfstruct *psf)
1018 {
1019 static double pos[POLY_MAXDIM];
1020 double *pl, *basis, fac;
1021 float *ppc;
1022 int i,n,p, ndim, npix;
1023
1024 npix = psf->masksize[0]*psf->masksize[1];
1025
1026 /* Reset the Local PSF mask */
1027 memset(psf->maskloc, 0, npix*sizeof(double));
1028
1029 /* Grab the context vector */
1030 ndim = psf->poly->ndim;
1031 for (i=0; i<ndim; i++)
1032 {
1033 ttypeconv(psf->context[i], &pos[i], psf->contexttyp[i],T_DOUBLE);
1034 pos[i] = (pos[i] - psf->contextoffset[i]) / psf->contextscale[i];
1035 }
1036 poly_func(psf->poly, pos);
1037
1038 basis = psf->poly->basis;
1039
1040 ppc = psf->maskcomp;
1041 /* Sum each component */
1042 for (n = (psf->maskdim>2?psf->masksize[2]:1); n--;)
1043 {
1044 pl = psf->maskloc;
1045 fac = *(basis++);
1046 for (p=npix; p--;)
1047 *(pl++) += fac**(ppc++);
1048 }
1049
1050 return;
1051 }
1052
1053
1054 /*****************************compute_gradient*********************************/
1055
compute_gradient(double * weight,int width,int height,double * masks,double * maskx,double * masky,double * m)1056 double *compute_gradient(double *weight,int width, int height,
1057 double *masks,double *maskx,double *masky
1058 ,double *m)
1059 {
1060 int x,y;
1061 double *w,*ps,*px,*py;
1062
1063 /*------ copy of the (weighted) PSF, with outer ring set to zero */
1064 ps = masks;
1065 w = weight;
1066 for (y=0; y<height; y++)
1067 for (x=0; x<width; x++, ps++, w++)
1068 *(m++) = y?(y>=(height-1)?0:(x?(x>=(width-1)?0:*ps**w):0)):0;
1069 /*------ (weighted) PSF gradient in x (kernel for first moment in x) */
1070 ps = masks;
1071 px = maskx;
1072 w = weight;
1073 for (y=0; y<height; y++)
1074 for (x=0; x<width; x++, ps++, w++)
1075 *(m++) = ((*px++) = (x?(x>=(width-1)?0:*(ps+1)-*(ps-1)):0))**w/2;
1076 /*------ (weighted) PSF gradient in y (kernel for first moment in y) */
1077 ps = masks;
1078 py = masky;
1079 w = weight;
1080 for (y=0; y<height; y++)
1081 for (x=0; x<width; x++, ps++, w++)
1082 *(m++) = (*(py++)=(y?(y>=(height-1)?0:*(ps+width)-*(ps-width)):0))
1083 **w/2;
1084 return m;
1085 }
1086
1087 /*****************************compute_gradient_phot*****************************
1088 ****/
1089
compute_gradient_phot(double * pweight,int width,int height,double * pmasks,double * pm)1090 double *compute_gradient_phot(double *pweight,int width, int height,
1091 double *pmasks,double *pm)
1092
1093 {
1094 int x,y;
1095 double *pw,*pps;
1096
1097 /*------ copy of the (weighted) PSF, with outer ring set to zero */
1098 pps = pmasks;
1099 pw = pweight;
1100 for (y=0; y<height; y++)
1101 for (x=0; x<width; x++, pps++, pw++)
1102 *(pm++) = y?(y>=(height-1)?0:(x?(x>=(width-1)?0:*pps**pw):0)):0;
1103
1104 return pm;
1105 }
1106
1107 /**************************compute_pos********************************/
1108
compute_pos(int * pnpsf,int * pconvflag,int * pnpsfflag,double radmin2,double radmax2,double r2,double * sol,double * flux,double * deltax,double * deltay,double * pdx,double * pdy)1109 void compute_pos(int *pnpsf,int *pconvflag,int *pnpsfflag,double radmin2,
1110 double radmax2,double r2,double *sol,double *flux
1111 ,double *deltax,double *deltay,double *pdx,double *pdy)
1112 {
1113 int j,k,convflag,npsfflag,npsf;
1114 double dx,dy;
1115
1116 dx=*pdx;
1117 dy=*pdy;
1118 convflag=*pconvflag;
1119 npsfflag=*pnpsfflag;
1120 npsf=*pnpsf;
1121 for (j=0; j<npsf; j++)
1122 {
1123 flux[j] = sol[j*PSF_NA];
1124 /*------ Update the PSF shifts */
1125 if (fabs(flux[j])>0.0)
1126 {
1127 dx = -sol[j*PSF_NA+1]/((npsf>1?2:1)*flux[j]);
1128 dy = -sol[j*PSF_NA+2]/((npsf>1?2:1)*flux[j]);
1129 }
1130
1131 deltax[j] += dx;
1132 deltay[j] += dy;
1133 /*------ Continue until all PSFs have come to a complete stop */
1134 if ((dx*dx+dy*dy) > radmin2)
1135 convflag = 1;
1136 }
1137 for (j=0; j<npsf; j++)
1138 {
1139 /*------ Exit if too much decentering or negative flux */
1140 for (k=j+1; k<npsf; k++)
1141 {
1142 dx = deltax[j]-deltax[k];
1143 dy = deltay[j]-deltay[k];
1144 if (dx*dx+dy*dy<r2/4.0)
1145 {
1146 flux[j] = -BIG;
1147 break;
1148 }
1149 }
1150 if (flux[j]<0.0
1151 || (deltax[j]*deltax[j] + deltay[j]*deltay[j]) > radmax2)
1152 {
1153 npsfflag = 0;
1154 convflag = 0;
1155 npsf--;
1156 break;
1157 }
1158 }
1159 *pdx=dx;
1160 *pdy=dy;
1161 *pconvflag=convflag;
1162 *pnpsfflag= npsfflag;
1163 *pnpsf=npsf;
1164 return;
1165 }
1166
1167 /**************************compute_pos_phot********************************/
1168
compute_pos_phot(int * pnpsf,double * sol,double * flux)1169 void compute_pos_phot(int *pnpsf,double *sol,double *flux)
1170 {
1171 int j,npsf;
1172 npsf=*pnpsf;
1173 for (j=0; j<npsf; j++)
1174 {
1175 flux[j] = sol[j];
1176 }
1177 *pnpsf=npsf;
1178 return;
1179 }
1180
1181
1182 /************************************compute_poserr*****************************
1183 *********/
1184
compute_poserr(int j,double * var,double * sol,obj2struct * obj2,double * x2,double * y2,double * xy)1185 void compute_poserr( int j,double *var,double *sol,obj2struct *obj2,double *x2,
1186 double *y2,double *xy)
1187 {
1188 double vara,covab,varb;
1189
1190 /*------ Variances and covariance along x and y */
1191 vara = *(var += PSF_NA+1);
1192 covab = *(++var);
1193 varb = *(var += PSF_NA);
1194 var += PSF_NA+1;
1195 obj2->poserrmx2_psf = (vara*x2[j]*x2[j]+varb*xy[j]*xy[j]
1196 +2*covab*x2[j]*xy[j])/(sol[0]*sol[0]);
1197 obj2->poserrmy2_psf = (varb*y2[j]*y2[j]+vara*xy[j]*xy[j]
1198 +2*covab*y2[j]*xy[j])/(sol[0]*sol[0]);
1199 obj2->poserrmxy_psf = (vara*x2[j]*xy[j]+varb*y2[j]*xy[j]
1200 +covab*(x2[j]*y2[j]+xy[j]*xy[j]))
1201 /(sol[0]*sol[0]);
1202
1203 /*------ If requested, translate variances to major and minor error axes... */
1204 if (FLAG(obj2.poserra_psf))
1205 {
1206 double pmx2,pmy2,temp,theta;
1207
1208 if (fabs(temp=obj2->poserrmx2_psf-obj2->poserrmy2_psf) > 0.0)
1209 theta = atan2(2.0 * obj2->poserrmxy_psf,temp) / 2.0;
1210 else
1211 theta = PI/4.0;
1212
1213 temp = sqrt(0.25*temp*temp+obj2->poserrmxy_psf*obj2->poserrmxy_psf);
1214 pmy2 = pmx2 = 0.5*(obj2->poserrmx2_psf+obj2->poserrmy2_psf);
1215 pmx2+=temp;
1216 pmy2-=temp;
1217
1218 obj2->poserra_psf = (float)sqrt(pmx2);
1219 obj2->poserrb_psf = (float)sqrt(pmy2);
1220 obj2->poserrtheta_psf = theta*180.0/PI;
1221 }
1222
1223 /*------ ...Or ellipse parameters */
1224 if (FLAG(obj2.poserr_cxx))
1225 {
1226 double xm2,ym2, xym, temp;
1227
1228 xm2 = obj2->poserrmx2_psf;
1229 ym2 = obj2->poserrmy2_psf;
1230 xym = obj2->poserrmxy_psf;
1231 obj2->poserrcxx_psf = (float)(ym2/(temp=xm2*ym2-xym*xym));
1232 obj2->poserrcyy_psf = (float)(xm2/temp);
1233 obj2->poserrcxy_psf = (float)(-2*xym/temp);
1234 }
1235 return;
1236 }
1237
1238
1239 /******************************** svdfit ************************************/
1240 /*
1241 General least-square fit A.x = b, based on Singular Value Decomposition (SVD).
1242 Loosely adapted from Numerical Recipes in C, 2nd Ed. (p. 671).
1243 Note: the a and v matrices are transposed with respect to the N.R. convention.
1244 */
svdfit(double * a,double * b,int m,int n,double * sol,double * vmat,double * wmat)1245 void svdfit(double *a, double *b, int m, int n, double *sol,
1246 double *vmat, double *wmat)
1247 {
1248 #define MAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
1249 (maxarg1) : (maxarg2))
1250 #define PYTHAG(a,b) ((at=fabs(a)) > (bt=fabs(b)) ? \
1251 (ct=bt/at,at*sqrt(1.0+ct*ct)) \
1252 : (bt ? (ct=at/bt,bt*sqrt(1.0+ct*ct)): 0.0))
1253 #define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
1254 #define TOL 1.0e-11
1255
1256 int flag,i,its,j,jj,k,l,nm,mmi,nml;
1257 double c,f,h,s,x,y,z,
1258 anorm, g, scale,
1259 at,bt,ct,maxarg1,maxarg2,
1260 thresh, wmax,
1261 *w,*ap,*ap0,*ap1,*ap10,*rv1p,*vp,*vp0,*vp1,*vp10,
1262 *bp,*tmpp, *rv1,*tmp;
1263
1264 anorm = g = scale = 0.0;
1265 if (m < n)
1266 error(EXIT_FAILURE, "*Error*: Not enough rows for solving the system ",
1267 "in svdfit()");
1268
1269 QMALLOC(rv1, double, n);
1270 QMALLOC(tmp, double, n);
1271 l = nm = nml = 0; /* To avoid gcc -Wall warnings */
1272 for (i=0;i<n;i++)
1273 {
1274 l = i+1;
1275 nml = n-l;
1276 rv1[i] = scale*g;
1277 g = s = scale = 0.0;
1278 if ((mmi = m - i) > 0)
1279 {
1280 ap = ap0 = a+i*(m+1);
1281 for (k=mmi;k--;)
1282 scale += fabs(*(ap++));
1283 if (scale)
1284 {
1285 for (ap=ap0,k=mmi; k--; ap++)
1286 {
1287 *ap /= scale;
1288 s += *ap**ap;
1289 }
1290 f = *ap0;
1291 g = -SIGN(sqrt(s),f);
1292 h = f*g-s;
1293 *ap0 = f-g;
1294 ap10 = a+l*m+i;
1295 for (j=nml;j--; ap10+=m)
1296 {
1297 for (s=0.0,ap=ap0,ap1=ap10,k=mmi; k--;)
1298 s += *(ap1++)**(ap++);
1299 f = s/h;
1300 for (ap=ap0,ap1=ap10,k=mmi; k--;)
1301 *(ap1++) += f**(ap++);
1302 }
1303 for (ap=ap0,k=mmi; k--;)
1304 *(ap++) *= scale;
1305 }
1306 }
1307 wmat[i] = scale*g;
1308 g = s = scale = 0.0;
1309 if (i < m && i+1 != n)
1310 {
1311 ap = ap0 = a+i+m*l;
1312 for (k=nml;k--; ap+=m)
1313 scale += fabs(*ap);
1314 if (scale)
1315 {
1316 for (ap=ap0,k=nml;k--; ap+=m)
1317 {
1318 *ap /= scale;
1319 s += *ap**ap;
1320 }
1321 f=*ap0;
1322 g = -SIGN(sqrt(s),f);
1323 h=f*g-s;
1324 *ap0=f-g;
1325 rv1p = rv1+l;
1326 for (ap=ap0,k=nml;k--; ap+=m)
1327 *(rv1p++) = *ap/h;
1328 ap10 = a+l+m*l;
1329 for (j=m-l; j--; ap10++)
1330 {
1331 for (s=0.0,ap=ap0,ap1=ap10,k=nml; k--; ap+=m,ap1+=m)
1332 s += *ap1**ap;
1333 rv1p = rv1+l;
1334 for (ap1=ap10,k=nml;k--; ap1+=m)
1335 *ap1 += s**(rv1p++);
1336 }
1337 for (ap=ap0,k=nml;k--; ap+=m)
1338 *ap *= scale;
1339 }
1340 }
1341 anorm=MAX(anorm,(fabs(wmat[i])+fabs(rv1[i])));
1342 }
1343
1344 for (i=n-1;i>=0;i--)
1345 {
1346 if (i < n-1)
1347 {
1348 if (g)
1349 {
1350 ap0 = a+l*m+i;
1351 vp0 = vmat+i*n+l;
1352 vp10 = vmat+l*n+l;
1353 g *= *ap0;
1354 for (ap=ap0,vp=vp0,j=nml; j--; ap+=m)
1355 *(vp++) = *ap/g;
1356 for (j=nml; j--; vp10+=n)
1357 {
1358 for (s=0.0,ap=ap0,vp1=vp10,k=nml; k--; ap+=m)
1359 s += *ap**(vp1++);
1360 for (vp=vp0,vp1=vp10,k=nml; k--;)
1361 *(vp1++) += s**(vp++);
1362 }
1363 }
1364 vp = vmat+l*n+i;
1365 vp1 = vmat+i*n+l;
1366 for (j=nml; j--; vp+=n)
1367 *vp = *(vp1++) = 0.0;
1368 }
1369 vmat[i*n+i]=1.0;
1370 g=rv1[i];
1371 l=i;
1372 nml = n-l;
1373 }
1374
1375 for (i=(m<n?m:n); --i>=0;)
1376 {
1377 l=i+1;
1378 nml = n-l;
1379 mmi=m-i;
1380 g=wmat[i];
1381 ap0 = a+i*m+i;
1382 ap10 = ap0 + m;
1383 for (ap=ap10,j=nml;j--;ap+=m)
1384 *ap=0.0;
1385 if (g)
1386 {
1387 g=1.0/g;
1388 for (j=nml;j--; ap10+=m)
1389 {
1390 for (s=0.0,ap=ap0,ap1=ap10,k=mmi; --k;)
1391 s += *(++ap)**(++ap1);
1392 f = (s/(*ap0))*g;
1393 for (ap=ap0,ap1=ap10,k=mmi;k--;)
1394 *(ap1++) += f**(ap++);
1395 }
1396 for (ap=ap0,j=mmi;j--;)
1397 *(ap++) *= g;
1398 }
1399 else
1400 for (ap=ap0,j=mmi;j--;)
1401 *(ap++)=0.0;
1402 ++(*ap0);
1403 }
1404
1405 for (k=n; --k>=0;)
1406 {
1407 for (its=0;its<100;its++)
1408 {
1409 flag=1;
1410 for (l=k;l>=0;l--)
1411 {
1412 nm=l-1;
1413 if (fabs(rv1[l])+anorm == anorm)
1414 {
1415 flag=0;
1416 break;
1417 }
1418 if (fabs(wmat[nm])+anorm == anorm)
1419 break;
1420 }
1421 if (flag)
1422 {
1423 c=0.0;
1424 s=1.0;
1425 ap0 = a+nm*m;
1426 ap10 = a+l*m;
1427 for (i=l; i<=k; i++,ap10+=m)
1428 {
1429 f=s*rv1[i];
1430 if (fabs(f)+anorm == anorm)
1431 break;
1432 g=wmat[i];
1433 h=PYTHAG(f,g);
1434 wmat[i]=h;
1435 h=1.0/h;
1436 c=g*h;
1437 s=(-f*h);
1438 for (ap=ap0,ap1=ap10,j=m; j--;)
1439 {
1440 z = *ap1;
1441 y = *ap;
1442 *(ap++) = y*c+z*s;
1443 *(ap1++) = z*c-y*s;
1444 }
1445 }
1446 }
1447 z=wmat[k];
1448 if (l == k)
1449 {
1450 if (z < 0.0)
1451 {
1452 wmat[k] = -z;
1453 vp = vmat+k*n;
1454 for (j=n; j--; vp++)
1455 *vp = (-*vp);
1456 }
1457 break;
1458 }
1459 if (its == 99)
1460 error(EXIT_FAILURE, "*Error*: No convergence in 100 SVD iterations ",
1461 "in svdfit()");
1462 x=wmat[l];
1463 nm=k-1;
1464 y=wmat[nm];
1465 g=rv1[nm];
1466 h=rv1[k];
1467 f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
1468 g=PYTHAG(f,1.0);
1469 f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x;
1470 c=s=1.0;
1471 ap10 = a+l*m;
1472 vp10 = vmat+l*n;
1473 for (j=l;j<=nm;j++,ap10+=m,vp10+=n)
1474 {
1475 i=j+1;
1476 g=rv1[i];
1477 y=wmat[i];
1478 h=s*g;
1479 g=c*g;
1480 z=PYTHAG(f,h);
1481 rv1[j]=z;
1482 c=f/z;
1483 s=h/z;
1484 f=x*c+g*s;
1485 g=g*c-x*s;
1486 h=y*s;
1487 y=y*c;
1488 for (vp=(vp1=vp10)+n,jj=n; jj--;)
1489 {
1490 z = *vp;
1491 x = *vp1;
1492 *(vp1++) = x*c+z*s;
1493 *(vp++) = z*c-x*s;
1494 }
1495 z=PYTHAG(f,h);
1496 wmat[j]=z;
1497 if (z)
1498 {
1499 z=1.0/z;
1500 c=f*z;
1501 s=h*z;
1502 }
1503 f=c*g+s*y;
1504 x=c*y-s*g;
1505 for (ap=(ap1=ap10)+m,jj=m; jj--;)
1506 {
1507 z = *ap;
1508 y = *ap1;
1509 *(ap1++) = y*c+z*s;
1510 *(ap++) = z*c-y*s;
1511 }
1512 }
1513 rv1[l]=0.0;
1514 rv1[k]=f;
1515 wmat[k]=x;
1516 }
1517 }
1518
1519 wmax=0.0;
1520 w = wmat;
1521 for (j=n;j--; w++)
1522 if (*w > wmax)
1523 wmax=*w;
1524 thresh=TOL*wmax;
1525 w = wmat;
1526 for (j=n;j--; w++)
1527 if (*w < thresh)
1528 *w = 0.0;
1529
1530 w = wmat;
1531 ap = a;
1532 tmpp = tmp;
1533 for (j=n; j--; w++)
1534 {
1535 s=0.0;
1536 if (*w)
1537 {
1538 bp = b;
1539 for (i=m; i--;)
1540 s += *(ap++)**(bp++);
1541 s /= *w;
1542 }
1543 else
1544 ap += m;
1545 *(tmpp++) = s;
1546 }
1547
1548 vp0 = vmat;
1549 for (j=0; j<n; j++,vp0++)
1550 {
1551 s=0.0;
1552 tmpp = tmp;
1553 for (vp=vp0,jj=n; jj--; vp+=n)
1554 s += *vp**(tmpp++);
1555 sol[j]=s;
1556 }
1557
1558 /* Free temporary arrays */
1559 free(tmp);
1560 free(rv1);
1561
1562 return;
1563 }
1564
1565 #undef SIGN
1566 #undef MAX
1567 #undef PYTHAG
1568 #undef TOL
1569
1570 /******************************** svdvar ************************************/
1571 /*
1572 Computation of the covariance matrix from the SVD vmat and wmat matrices.A
1573 dapted from Numerical Recipes in C, 2nd Ed. (p. 679).
1574 */
svdvar(double * v,double * w,int n,double * cov)1575 void svdvar(double *v, double *w, int n, double *cov)
1576 {
1577 static double wti[PSF_NTOT];
1578 double sum;
1579 int i,j,k;
1580
1581 for (i=0; i<n; i++)
1582 wti[i] = w[i]? 1.0/(w[i]*w[i]) : 0.0;
1583
1584 for (i=0; i<n; i++)
1585 for (j=0; j<=i; j++)
1586 {
1587 for (sum=0.0,k=0; k<n; k++)
1588 sum += v[k*n+i]*v[k*n+j]*wti[k];
1589 cov[j*n+i] = cov[i*n+j] = sum;
1590 }
1591
1592 return;
1593 }
1594
1595