public software.sextractor

[/] [branches/] [multi/] [src/] [profit.c] - Blame information for rev 285

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 233 bertin
/*
2
*                               profit.c
3 42 bertin
*
4 233 bertin
* Fit a range of galaxy models to an image.
5 42 bertin
*
6 233 bertin
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 42 bertin
*
8 233 bertin
*       This file part of:      SExtractor
9 42 bertin
*
10 278 bertin
*       Copyright:              (C) 2006-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
11 42 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 285 bertin
*       Last modified:          06/05/2012
26 233 bertin
*
27
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
28 42 bertin
 
29
#ifdef HAVE_CONFIG_H
30
#include        "config.h"
31
#endif
32
 
33 206 bertin
#ifndef HAVE_MATHIMF_H
34
#define _GNU_SOURCE
35
#endif
36
 
37 201 bertin
#include        <math.h>
38 42 bertin
#include        <stdio.h>
39
#include        <stdlib.h>
40
#include        <string.h>
41
 
42
#include        "define.h"
43
#include        "globals.h"
44
#include        "prefs.h"
45
#include        "fits/fitscat.h"
46 228 bertin
#include        "levmar/levmar.h"
47 50 bertin
#include        "fft.h"
48 123 bertin
#include        "fitswcs.h"
49 42 bertin
#include        "check.h"
50 274 bertin
#include        "graph.h"
51 229 bertin
#include        "image.h"
52 274 bertin
#include        "misc.h"
53 141 bertin
#include        "pattern.h"
54 50 bertin
#include        "psf.h"
55 42 bertin
#include        "profit.h"
56 283 bertin
#include        "subimage.h"
57 42 bertin
 
58 207 bertin
static double   prof_gammainc(double x, double a),
59
                prof_gamma(double x);
60 201 bertin
static float    prof_interpolate(profstruct *prof, float *posin);
61
static float    interpolate_pix(float *posin, float *pix, int *naxisn,
62 52 bertin
                interpenum interptype);
63 50 bertin
 
64 201 bertin
static void     make_kernel(float pos, float *kernel, interpenum interptype);
65 50 bertin
 
66 42 bertin
/*------------------------------- variables ---------------------------------*/
67
 
68 209 bertin
const int       interp_kernwidth[5]={1,2,4,6,8};
69
 
70 233 bertin
const int       flux_flag[PARAM_NPARAM] = {0,
71 285 bertin
                                        1,2,3,4,5,0,0,
72
                                        1,2,3,4,5,0,0,0,0,
73
                                        1,2,3,4,5,0,0,0,
74 209 bertin
                                        1,0,0,0,0,0,0,0,
75
                                        1,0,0,
76
                                        1,0,0,
77
                                        1,0,0
78
                                        };
79
 
80 267 bertin
int     the_gal;
81
 
82 46 bertin
/****** profit_init ***********************************************************
83 285 bertin
PROTO   profitstruct profit_init(obj2struct *obj2, unsigned int modeltype)
84 46 bertin
PURPOSE Allocate and initialize a new profile-fitting structure.
85 285 bertin
INPUT   pointer to the obj2,
86 258 bertin
        model type.
87 46 bertin
OUTPUT  A pointer to an allocated profit structure.
88 117 bertin
NOTES   -.
89 42 bertin
AUTHOR  E. Bertin (IAP)
90 285 bertin
VERSION 06/05/2012
91 42 bertin
 ***/
92 285 bertin
profitstruct    *profit_init(obj2struct *obj2, unsigned int modeltype)
93 42 bertin
  {
94 46 bertin
   profitstruct         *profit;
95 285 bertin
   subprofitstruct      *subprofit;
96
   subimagestruct       *subimage;
97 267 bertin
   double               psf_fwhm;
98 285 bertin
   int                  s,t, nmodels, npix, nsub;
99 42 bertin
 
100 49 bertin
  QCALLOC(profit, profitstruct, 1);
101 285 bertin
 
102 258 bertin
  profit->obj2 = obj2;
103
 
104 285 bertin
  QCALLOC(profit->subprofit, subprofitstruct, obj2->nsubimage);
105
  subprofit = profit->subprofit;
106
  subimage = obj2->subimage;
107
  nsub = 0;
108
  for (s=0; s<obj2->nsubimage; s++, subimage++)
109
    {
110
    subprofit->field = subimage->field;
111
    subprofit->wfield = subimage->wfield;
112
    subprofit->psf = subimage->field->psf;
113
    subprofit->pixstep = subprofit->psf->pixstep;
114
    subprofit->sigma = obj2->sigbkg[s];
115
/*-- Set initial guesses and boundaries */
116
    if ((subprofit->guessradius = 0.5*subprofit->psf->fwhm) < obj2->hl_radius)
117
      subprofit->guessradius = obj2->hl_radius;
118
    if ((subprofit->guessflux = obj2->flux_auto[s]) <= 0.0)
119
      subprofit->guessflux = 0.0;
120
    if ((subprofit->guessfluxmax = 10.0*obj2->fluxerr_auto[s])
121
        <= subprofit->guessflux)
122
      subprofit->guessfluxmax = subprofit->guessflux;
123
    if (subprofit->guessfluxmax <= 0.0)
124
      subprofit->guessfluxmax = 1.0;
125
    subprofit->fluxfac = 1.0;   /* Default */
126
/*-- Create pixmaps at image resolution */
127
    subprofit->ix = (int)(obj2->mx + 0.49999); /* 1st pix=0 */
128
    subprofit->iy = (int)(obj2->my + 0.49999); /* 1st pix=0 */
129
    psf_fwhm = subprofit->psf->masksize[0]*subprofit->psf->pixstep;
130
    subprofit->objnaxisn[0] = ((int)((subimage->imsize[0] + psf_fwhm + 0.499)
131 258 bertin
                *1.2)/2)*2 + 1;
132 285 bertin
    subprofit->objnaxisn[1] = ((int)((subimage->imsize[1] + psf_fwhm + 0.499)
133 258 bertin
                *1.2)/2)*2 + 1;
134 285 bertin
    if (subprofit->objnaxisn[1] < subprofit->objnaxisn[0])
135
      subprofit->objnaxisn[1] = subprofit->objnaxisn[0];
136
    else
137
      subprofit->objnaxisn[0] = subprofit->objnaxisn[1];
138 258 bertin
 
139 285 bertin
    if (subprofit->objnaxisn[0]>PROFIT_MAXOBJSIZE)
140
      {
141
      subprofit->subsamp
142
                = ceil((float)subprofit->objnaxisn[0]/PROFIT_MAXOBJSIZE);
143
      subprofit->objnaxisn[1]
144
                = (subprofit->objnaxisn[0] /= (int)subprofit->subsamp);
145
      profit->flag |= PROFLAG_OBJSUB;   /* !CHECK */
146
      }
147
    else
148
      subprofit->subsamp = 1.0;
149
    subprofit->nobjpix = subprofit->objnaxisn[0]*subprofit->objnaxisn[1];
150 258 bertin
 
151
/* Allocate memory for the data and the resampled model */
152 285 bertin
    QMALLOC16(subprofit->objpix, PIXTYPE, subprofit->nobjpix);
153
    QMALLOC16(subprofit->objweight, PIXTYPE, subprofit->nobjpix);
154
    QMALLOC16(subprofit->lmodpix, PIXTYPE, subprofit->nobjpix);
155
    QMALLOC16(subprofit->lmodpix2, PIXTYPE, subprofit->nobjpix);
156 258 bertin
 
157 285 bertin
/*-- Create pixmap at PSF resolution */
158
    subprofit->modnaxisn[0] =
159
        ((int)(subprofit->objnaxisn[0]*subprofit->subsamp
160
                /subprofit->pixstep+0.4999)/2+1)*2;
161
    subprofit->modnaxisn[1] =
162
        ((int)(subprofit->objnaxisn[1]*subprofit->subsamp
163
                /subprofit->pixstep+0.4999)/2+1)*2;
164
    if (subprofit->modnaxisn[1] < subprofit->modnaxisn[0])
165
      subprofit->modnaxisn[1] = subprofit->modnaxisn[0];
166
    else
167
      subprofit->modnaxisn[0] = subprofit->modnaxisn[1];
168
    if (subprofit->modnaxisn[0]>PROFIT_MAXMODSIZE)
169
      {
170
      subprofit->pixstep = (double)subprofit->modnaxisn[0] / PROFIT_MAXMODSIZE;
171
      subprofit->modnaxisn[0] = subprofit->modnaxisn[1] = PROFIT_MAXMODSIZE;
172
      profit->flag |= PROFLAG_MODSUB;
173
      }
174
    subprofit->nmodpix = subprofit->modnaxisn[0]*subprofit->modnaxisn[1];
175 258 bertin
 
176 285 bertin
/* Allocate memory for the complete model */
177
    QMALLOC16(subprofit->modpix, float, subprofit->nmodpix);
178
    QMALLOC16(subprofit->modpix2, float, subprofit->nmodpix);
179
    QMALLOC16(subprofit->cmodpix, float, subprofit->nmodpix);
180
    QMALLOC16(subprofit->psfpix, float, subprofit->nmodpix);
181
    if ((npix = subprofit_copyobjpix(subprofit, subimage)))
182
      {
183
      profit->nresi += npix;
184
/*---- Compute the local PSF */
185
      subprofit_psf(subprofit, obj2);
186
      subprofit->index = nsub++;
187
      subprofit++;
188
      }
189
    else
190
      subprofit_end(subprofit);
191 258 bertin
    }
192
 
193 285 bertin
  profit->nsubprofit = nsub;
194 258 bertin
 
195 285 bertin
  QMALLOC16(profit->resi, float, profit->nresi);
196 258 bertin
 
197 245 bertin
  QMALLOC(profit->prof, profstruct *, MODEL_NMAX);
198
  nmodels = 0;
199
  for (t=1; t<(1<<MODEL_NMAX); t<<=1)
200
    if (modeltype&t)
201
      profit->prof[nmodels++] = prof_init(profit, t);
202
  profit->nprof = nmodels;
203 50 bertin
 
204 258 bertin
  QMALLOC16(profit->covar, float, profit->nparam*profit->nparam);
205
 
206
  profit_resetparams(profit);
207
 
208 46 bertin
  return profit;
209 42 bertin
  }
210
 
211
 
212 46 bertin
/****** profit_end ************************************************************
213 285 bertin
PROTO   void profit_end(profitstruct *profit)
214 46 bertin
PURPOSE End (deallocate) a profile-fitting structure.
215 42 bertin
INPUT   Prof structure.
216
OUTPUT  -.
217 117 bertin
NOTES   -.
218 42 bertin
AUTHOR  E. Bertin (IAP)
219 285 bertin
VERSION 06/05/2012
220 42 bertin
 ***/
221 46 bertin
void    profit_end(profitstruct *profit)
222 42 bertin
  {
223 285 bertin
   subprofitstruct      *subprofit;
224
   int                  p,s;
225 49 bertin
 
226
  for (p=0; p<profit->nprof; p++)
227
    prof_end(profit->prof[p]);
228 270 bertin
  free(profit->prof);
229 285 bertin
  subprofit = profit->subprofit;
230
  for (s=profit->nsubprofit; s--;)
231
    subprofit_end(subprofit++);
232
 
233
  free(profit->subprofit);
234 221 bertin
  free(profit->resi);
235 123 bertin
  free(profit->covar);
236 285 bertin
 
237 46 bertin
  free(profit);
238 42 bertin
 
239
  return;
240
  }
241
 
242
 
243 285 bertin
/****** subprofit_end *********************************************************
244
PROTO   void subprofit_end(subprofitstruct *profit)
245
PURPOSE End (free content) a subprofile-fitting structure.
246
INPUT   SubProfit structure.
247
OUTPUT  -.
248
NOTES   -.
249
AUTHOR  E. Bertin (IAP)
250
VERSION 03/05/2012
251
 ***/
252
void    subprofit_end(subprofitstruct *subprofit)
253
  {
254
  free(subprofit->objpix);
255
  free(subprofit->objweight);
256
  free(subprofit->lmodpix);
257
  free(subprofit->lmodpix2);
258
  free(subprofit->modpix);
259
  free(subprofit->modpix2);
260
  free(subprofit->cmodpix);
261
  free(subprofit->psfpix);
262
  free(subprofit->psfdft);
263
 
264
  return;
265
  }
266
 
267
 
268 65 bertin
/****** profit_fit ************************************************************
269 259 bertin
PROTO   int profit_fit(profitstruct *profit, obj2struct *obj2)
270 46 bertin
PURPOSE Fit profile(s) convolved with the PSF to a detected object.
271 259 bertin
INPUT   Pointer to the model structure,
272
        pointer to the obj2.
273
OUTPUT  Number of minimization iterations (0 in case of error).
274
NOTES   -.
275 46 bertin
AUTHOR  E. Bertin (IAP)
276 267 bertin
VERSION 06/10/2011
277 46 bertin
 ***/
278 259 bertin
int     profit_fit(profitstruct *profit, obj2struct *obj2)
279 46 bertin
  {
280 259 bertin
   int                  p, nparam, nparam2;
281 46 bertin
 
282 128 bertin
  nparam = profit->nparam;
283 208 bertin
  nparam2 = nparam*nparam;
284 52 bertin
 
285 258 bertin
/* reset all flags except those set in profit_init() */
286 259 bertin
  profit->flag &= PROFLAG_OBJSUB|PROFLAG_MODSUB;
287 52 bertin
 
288 211 bertin
/* Check if the number of constraints exceeds the number of free parameters */
289 128 bertin
  if (profit->nresi < nparam)
290 49 bertin
    {
291 176 bertin
    if (FLAG(obj2.prof_vector))
292
      for (p=0; p<nparam; p++)
293
        obj2->prof_vector[p] = 0.0;
294 208 bertin
    if (FLAG(obj2.prof_errvector))
295
      for (p=0; p<nparam; p++)
296
        obj2->prof_errvector[p] = 0.0;
297
    if (FLAG(obj2.prof_errmatrix))
298
      for (p=0; p<nparam2; p++)
299
        obj2->prof_errmatrix[p] = 0.0;
300 50 bertin
    obj2->prof_niter = 0;
301 259 bertin
    profit->flag |= PROFLAG_NOTCONST;
302
    return 0;
303 49 bertin
    }
304 48 bertin
 
305 123 bertin
/* Actual minimisation */
306 209 bertin
  fft_reset();
307 259 bertin
  return (profit->niter = profit_minimize(profit, PROFIT_MAXITER));
308
  }
309 235 bertin
 
310 259 bertin
 
311
/****** profit_measure ********************************************************
312
PROTO   void profit_measure(profitstruct *profit, obj2struct *obj2)
313
PURPOSE Perform measurements on the fitted model.
314
INPUT   Pointer to the model structure,
315
        pointer to the obj2.
316
OUTPUT  -.
317
NOTES   -.
318
AUTHOR  E. Bertin (IAP)
319 285 bertin
VERSION 06/05/2012
320 259 bertin
 ***/
321 267 bertin
void    profit_measure(profitstruct *profit, obj2struct *obj2)
322 259 bertin
  {
323
    profitstruct        *pprofit, *qprofit;
324 285 bertin
    subprofitstruct     *subprofit;
325 259 bertin
    patternstruct       *pattern;
326
    checkstruct         *check;
327
    double              emx2,emy2,emxy, a , cp,sp, cn, bn, n;
328
    float               param0[PARAM_NPARAM], param1[PARAM_NPARAM],
329
                        param[PARAM_NPARAM],
330
                        **list,
331
                        *cov,
332 285 bertin
                        psf_fwhm, err, aspect, chi2, flux;
333 259 bertin
    int                 *index,
334 285 bertin
                        i,j,p,s, nparam, nparam2, ncomp, nsub;
335 269 bertin
 
336
  nparam = profit->nparam;
337
  nparam2 = nparam*nparam;
338 285 bertin
  nsub = profit->nsubprofit;
339
  subprofit = profit->subprofit;
340 269 bertin
 
341 247 bertin
  if (profit->nlimmin)
342 259 bertin
    profit->flag |= PROFLAG_MINLIM;
343 247 bertin
  if (profit->nlimmax)
344 259 bertin
    profit->flag |= PROFLAG_MAXLIM;
345 247 bertin
 
346 128 bertin
  for (p=0; p<nparam; p++)
347
    profit->paramerr[p]= sqrt(profit->covar[p*(nparam+1)]);
348 126 bertin
 
349 268 bertin
/* Clean up FFTW plans */
350
  fft_reset();
351
 
352 50 bertin
/* CHECK-Images */
353 221 bertin
  if ((check = prefs.check[CHECK_PROFILES]))
354
    {
355 258 bertin
    profit_residuals(profit, 0.0, profit->paraminit, NULL);
356 285 bertin
    check_add(check, subprofit->lmodpix,
357
                subprofit->objnaxisn[0],subprofit->objnaxisn[1],
358
                subprofit->ix,subprofit->iy, 1.0);
359 221 bertin
    }
360 235 bertin
 
361 50 bertin
  if ((check = prefs.check[CHECK_SUBPROFILES]))
362 51 bertin
    {
363 258 bertin
    profit_residuals(profit, 0.0, profit->paraminit, NULL);
364 285 bertin
    check_add(check, subprofit->lmodpix,
365
                subprofit->objnaxisn[0],subprofit->objnaxisn[1],
366
                subprofit->ix,subprofit->iy, -1.0);
367 51 bertin
    }
368 221 bertin
  if ((check = prefs.check[CHECK_SPHEROIDS]))
369 51 bertin
    {
370 221 bertin
/*-- Set to 0 flux components that do not belong to spheroids */
371
    for (p=0; p<profit->nparam; p++)
372
      param[p] = profit->paraminit[p];
373
    list = profit->paramlist;
374
    index = profit->paramindex;
375
    for (i=0; i<PARAM_NPARAM; i++)
376
      if (list[i] && flux_flag[i] && i!= PARAM_SPHEROID_FLUX)
377
        param[index[i]] = 0.0;
378 258 bertin
    profit_residuals(profit, 0.0, param, NULL);
379 285 bertin
    check_add(check, subprofit->lmodpix,
380
                subprofit->objnaxisn[0],subprofit->objnaxisn[1],
381
                subprofit->ix,subprofit->iy, 1.0);
382 51 bertin
    }
383 221 bertin
  if ((check = prefs.check[CHECK_SUBSPHEROIDS]))
384
    {
385
/*-- Set to 0 flux components that do not belong to spheroids */
386
    for (p=0; p<profit->nparam; p++)
387
      param[p] = profit->paraminit[p];
388
    list = profit->paramlist;
389
    index = profit->paramindex;
390
    for (i=0; i<PARAM_NPARAM; i++)
391
      if (list[i] && flux_flag[i] && i!= PARAM_SPHEROID_FLUX)
392
        param[index[i]] = 0.0;
393 258 bertin
    profit_residuals(profit, 0.0, param, NULL);
394 285 bertin
    check_add(check, subprofit->lmodpix,
395
                subprofit->objnaxisn[0],subprofit->objnaxisn[1],
396
                subprofit->ix,subprofit->iy, -1.0);
397 221 bertin
    }
398
  if ((check = prefs.check[CHECK_DISKS]))
399
    {
400
/*-- Set to 0 flux components that do not belong to disks */
401
    for (p=0; p<profit->nparam; p++)
402
      param[p] = profit->paraminit[p];
403
    list = profit->paramlist;
404
    index = profit->paramindex;
405
    for (i=0; i<PARAM_NPARAM; i++)
406
      if (list[i] && flux_flag[i] && i!= PARAM_DISK_FLUX)
407
        param[index[i]] = 0.0;
408 258 bertin
    profit_residuals(profit, 0.0, param, NULL);
409 285 bertin
    check_add(check, subprofit->lmodpix,
410
                subprofit->objnaxisn[0],subprofit->objnaxisn[1],
411
                subprofit->ix,subprofit->iy, 1.0);
412 221 bertin
    }
413
  if ((check = prefs.check[CHECK_SUBDISKS]))
414
    {
415
/*-- Set to 0 flux components that do not belong to disks */
416
    for (p=0; p<profit->nparam; p++)
417
      param[p] = profit->paraminit[p];
418
    list = profit->paramlist;
419
    index = profit->paramindex;
420
    for (i=0; i<PARAM_NPARAM; i++)
421
      if (list[i] && flux_flag[i] && i!= PARAM_DISK_FLUX)
422
        param[index[i]] = 0.0;
423 258 bertin
    profit_residuals(profit, 0.0, param, NULL);
424 285 bertin
    check_add(check, subprofit->lmodpix,
425
                subprofit->objnaxisn[0],subprofit->objnaxisn[1],
426
                subprofit->ix,subprofit->iy, -1.0);
427 221 bertin
    }
428 235 bertin
 
429 221 bertin
/* Compute compressed residuals */
430 258 bertin
  profit_residuals(profit, 10.0, profit->paraminit,profit->resi);
431 221 bertin
 
432 73 bertin
/* Fill measurement parameters */
433 69 bertin
  if (FLAG(obj2.prof_vector))
434
    {
435 128 bertin
    for (p=0; p<nparam; p++)
436 221 bertin
      obj2->prof_vector[p]= profit->paraminit[p];
437 69 bertin
    }
438 127 bertin
  if (FLAG(obj2.prof_errvector))
439
    {
440 128 bertin
    for (p=0; p<nparam; p++)
441
      obj2->prof_errvector[p]= profit->paramerr[p];
442 127 bertin
    }
443 208 bertin
  if (FLAG(obj2.prof_errmatrix))
444
    {
445
    for (p=0; p<nparam2; p++)
446
      obj2->prof_errmatrix[p]= profit->covar[p];
447
    }
448 89 bertin
 
449 79 bertin
  obj2->prof_niter = profit->niter;
450 285 bertin
  subprofit = profit->subprofit;
451
  for (s=0; s<nsub; s++, subprofit++)
452 194 bertin
    {
453 285 bertin
    if (FLAG(obj2.flux_prof))
454
      obj2->flux_prof[s] = subprofit->flux;
455
    if (FLAG(obj2.mag_prof))
456
      obj2->mag_prof[s] = subprofit->flux>0.0?
457
                -FDMAG * log(subprofit->flux) + prefs.mag_zeropoint[s] : 99.0f;
458
    if (FLAG(obj2.fluxerr_prof) || FLAG(obj2.magerr_prof))
459
      {
460
      err = 0.0;
461
      cov = profit->covar;
462
      index = profit->paramindex;
463
      list = profit->paramlist;
464
      for (i=0; i<PARAM_NPARAM; i++)
465
        if (flux_flag[i]==s+1 && list[i])
466
          {
467
          cov = profit->covar + nparam*index[i];
468
          for (j=0; j<PARAM_NPARAM; j++)
469
            if (flux_flag[j]==s+1 && list[j])
470
              err += cov[index[j]];
471
          }
472
      if (err<0.0)
473
        err = 0.0;
474
      if (FLAG(obj2.fluxerr_prof))
475
        obj2->fluxerr_prof[s] = sqrt(err);
476
      if (FLAG(obj2.magerr_prof))
477
        obj2->magerr_prof[s] = subprofit->flux>0.0?
478
                FDMAG * sqrt(err) / subprofit->flux : 99.0f;
479
      }
480 194 bertin
    }
481
 
482 180 bertin
  obj2->prof_chi2 = (profit->nresi > profit->nparam)?
483
                profit->chi2 / (profit->nresi - profit->nparam) : 0.0;
484
 
485 245 bertin
/* Position */
486 131 bertin
  if (FLAG(obj2.x_prof))
487
    {
488
    i = profit->paramindex[PARAM_X];
489
    j = profit->paramindex[PARAM_Y];
490 145 bertin
/*-- Model coordinates follow the FITS convention (first pixel at 1,1) */
491 179 bertin
    if (profit->paramlist[PARAM_X])
492
      {
493 285 bertin
      obj2->x_prof = (double)profit->subprofit->ix + *profit->paramlist[PARAM_X]
494
                        + 1.0;  /* !CHECK */    /* FITS convention */
495 179 bertin
      obj2->poserrmx2_prof = emx2 = profit->covar[i*(nparam+1)];
496
      }
497
    else
498
      emx2 = 0.0;
499
    if (profit->paramlist[PARAM_Y])
500
      {
501 285 bertin
      obj2->y_prof = (double)profit->subprofit->iy + *profit->paramlist[PARAM_Y]
502
                        + 1.0;  /* !CHECK */    /* FITS convention */
503 179 bertin
      obj2->poserrmy2_prof = emy2 = profit->covar[j*(nparam+1)];
504
      }
505
    else
506
      emy2 = 0.0;
507
    if (profit->paramlist[PARAM_X] && profit->paramlist[PARAM_Y])
508
      obj2->poserrmxy_prof = emxy = profit->covar[i+j*nparam];
509
    else
510
      emxy = 0.0;
511 131 bertin
 
512 135 bertin
/*-- Error ellipse parameters */
513
    if (FLAG(obj2.poserra_prof))
514
      {
515
       double   pmx2,pmy2,temp,theta;
516 131 bertin
 
517 135 bertin
      if (fabs(temp=emx2-emy2) > 0.0)
518
        theta = atan2(2.0 * emxy,temp) / 2.0;
519
      else
520
        theta = PI/4.0;
521 131 bertin
 
522 135 bertin
      temp = sqrt(0.25*temp*temp+ emxy*emxy);
523
      pmy2 = pmx2 = 0.5*(emx2+emy2);
524
      pmx2+=temp;
525
      pmy2-=temp;
526 131 bertin
 
527 135 bertin
      obj2->poserra_prof = (float)sqrt(pmx2);
528
      obj2->poserrb_prof = (float)sqrt(pmy2);
529 201 bertin
      obj2->poserrtheta_prof = (float)(theta/DEG);
530 135 bertin
      }
531 131 bertin
 
532 135 bertin
    if (FLAG(obj2.poserrcxx_prof))
533
      {
534
       double   temp;
535 131 bertin
 
536 135 bertin
      obj2->poserrcxx_prof = (float)(emy2/(temp=emx2*emy2-emxy*emxy));
537
      obj2->poserrcyy_prof = (float)(emx2/temp);
538
      obj2->poserrcxy_prof = (float)(-2*emxy/temp);
539
      }
540 131 bertin
    }
541
 
542 245 bertin
/* Equivalent noise area */
543
  if (FLAG(obj2.prof_noisearea))
544 285 bertin
    obj2->prof_noisearea = subprofit_noisearea(profit->subprofit); /* !CHECK */
545 245 bertin
 
546
/* Second order moments and ellipticities */
547
  if (FLAG(obj2.prof_mx2))
548
    profit_moments(profit, obj2);
549
 
550
/* Second order moments of the convolved model (used by other parameters) */
551
  if (FLAG(obj2.prof_convmx2))
552 285 bertin
    subprofit_convmoments(profit->subprofit, obj2); /* !CHECK */
553 245 bertin
 
554
/* "Hybrid" magnitudes */
555 235 bertin
  if (FLAG(obj2.fluxcor_prof))
556
    {
557 258 bertin
    profit_residuals(profit, 0.0, profit->paraminit, NULL);
558 285 bertin
    subprofit_fluxcor(profit->subprofit, obj2); /* !CHECK */
559 235 bertin
    }
560
 
561 221 bertin
/* Do measurements on the rasterised model (surface brightnesses) */
562 266 bertin
  if (FLAG(obj2.fluxeff_prof))
563 285 bertin
    subprofit_surface(profit, profit->subprofit, obj2);  /* !CHECK */
564 207 bertin
 
565 233 bertin
/* Background offset */
566
  if (FLAG(obj2.prof_offset_flux))
567
    {
568
    obj2->prof_offset_flux = *profit->paramlist[PARAM_BACK];
569
    obj2->prof_offset_fluxerr=profit->paramerr[profit->paramindex[PARAM_BACK]];
570
    }
571
 
572
/* Point source */
573
  if (FLAG(obj2.prof_dirac_flux))
574
    {
575 285 bertin
    for (s=0; s<nsub; s++)
576
      {
577
      flux = obj2->prof_dirac_flux[s] = *profit->paramlist[PARAM_DIRAC_FLUX+s];
578
      if (FLAG(obj2.prof_dirac_mag))
579
        obj2->prof_dirac_mag[s] = flux>0.0f?
580
                 -FDMAG * log(flux) + prefs.mag_zeropoint[s] : 99.0f;
581
      if (FLAG(obj2.prof_dirac_fluxerr))
582
        obj2->prof_dirac_fluxerr[s] =
583
                profit->paramerr[profit->paramindex[PARAM_DIRAC_FLUX+s]];
584
      if (FLAG(obj2.prof_dirac_magerr))
585
        obj2->prof_dirac_magerr[s] = flux>0.0f?
586
                 FDMAG*profit->paramerr[profit->paramindex[PARAM_DIRAC_FLUX+s]]
587
                / flux : 99.0f;
588
      }
589 233 bertin
    }
590
 
591 207 bertin
/* Spheroid */
592 115 bertin
  if (FLAG(obj2.prof_spheroid_flux))
593
    {
594 218 bertin
    if ((aspect = *profit->paramlist[PARAM_SPHEROID_ASPECT]) > 1.0)
595
      {
596
      *profit->paramlist[PARAM_SPHEROID_REFF] *= aspect;
597
      profit->paramerr[profit->paramindex[PARAM_SPHEROID_REFF]] *= aspect;
598
      profit->paramerr[profit->paramindex[PARAM_SPHEROID_ASPECT]]
599
                        /= (aspect*aspect);
600
      *profit->paramlist[PARAM_SPHEROID_ASPECT] = 1.0 / aspect;
601
      *profit->paramlist[PARAM_SPHEROID_POSANG] += 90.0;
602
      }
603 285 bertin
    for (s=0; s<nsub; s++)
604
      {
605
      flux = obj2->prof_spheroid_flux[s]
606
                = *profit->paramlist[PARAM_SPHEROID_FLUX+s];
607
      if (FLAG(obj2.prof_spheroid_mag))
608
        obj2->prof_spheroid_mag[s] = flux>0.0f?
609
                 -FDMAG * log(flux) + prefs.mag_zeropoint[s] : 99.0f;
610
      if (FLAG(obj2.prof_spheroid_fluxerr))
611
        obj2->prof_spheroid_fluxerr[s] =
612
                profit->paramerr[profit->paramindex[PARAM_SPHEROID_FLUX+s]];
613
      if (FLAG(obj2.prof_spheroid_magerr))
614
        obj2->prof_spheroid_magerr[s] = flux>0.0f?
615
              FDMAG*profit->paramerr[profit->paramindex[PARAM_SPHEROID_FLUX+s]]
616
                / flux : 99.0f;
617
      }
618 115 bertin
    obj2->prof_spheroid_reff = *profit->paramlist[PARAM_SPHEROID_REFF];
619 128 bertin
    obj2->prof_spheroid_refferr =
620
                profit->paramerr[profit->paramindex[PARAM_SPHEROID_REFF]];
621 115 bertin
    obj2->prof_spheroid_aspect = *profit->paramlist[PARAM_SPHEROID_ASPECT];
622 128 bertin
    obj2->prof_spheroid_aspecterr =
623
                profit->paramerr[profit->paramindex[PARAM_SPHEROID_ASPECT]];
624 123 bertin
    obj2->prof_spheroid_theta =
625
                        fmod_m90_p90(*profit->paramlist[PARAM_SPHEROID_POSANG]);
626 128 bertin
    obj2->prof_spheroid_thetaerr =
627
                profit->paramerr[profit->paramindex[PARAM_SPHEROID_POSANG]];
628 115 bertin
    if (FLAG(obj2.prof_spheroid_sersicn))
629 128 bertin
      {
630 116 bertin
      obj2->prof_spheroid_sersicn = *profit->paramlist[PARAM_SPHEROID_SERSICN];
631 128 bertin
      obj2->prof_spheroid_sersicnerr =
632
                profit->paramerr[profit->paramindex[PARAM_SPHEROID_SERSICN]];
633
      }
634 206 bertin
    else
635
      obj2->prof_spheroid_sersicn = 4.0;
636
    if (FLAG(obj2.prof_spheroid_peak))
637
      {
638
      n = obj2->prof_spheroid_sersicn;
639 207 bertin
      bn = 2.0*n - 1.0/3.0 + 4.0/(405.0*n) + 46.0/(25515.0*n*n)
640 206 bertin
                + 131.0/(1148175*n*n*n);        /* Ciotti & Bertin 1999 */
641 207 bertin
      cn = n * prof_gamma(2.0*n) * pow(bn, -2.0*n);
642 206 bertin
      obj2->prof_spheroid_peak = obj2->prof_spheroid_reff>0.0?
643 285 bertin
        obj2->prof_spheroid_flux[0]
644 207 bertin
                / (2.0 * PI * cn
645 206 bertin
                * obj2->prof_spheroid_reff*obj2->prof_spheroid_reff
646 207 bertin
                * obj2->prof_spheroid_aspect)
647 206 bertin
        : 0.0;
648
      if (FLAG(obj2.prof_spheroid_fluxeff))
649 207 bertin
        obj2->prof_spheroid_fluxeff = obj2->prof_spheroid_peak * exp(-bn);
650
      if (FLAG(obj2.prof_spheroid_fluxmean))
651
        obj2->prof_spheroid_fluxmean = obj2->prof_spheroid_peak * cn;
652 206 bertin
      }
653 115 bertin
    }
654 50 bertin
 
655 132 bertin
/* Disk */
656 115 bertin
  if (FLAG(obj2.prof_disk_flux))
657
    {
658 218 bertin
    if ((aspect = *profit->paramlist[PARAM_DISK_ASPECT]) > 1.0)
659
      {
660
      *profit->paramlist[PARAM_DISK_SCALE] *= aspect;
661
      profit->paramerr[profit->paramindex[PARAM_DISK_SCALE]] *= aspect;
662
      profit->paramerr[profit->paramindex[PARAM_DISK_ASPECT]]
663
                        /= (aspect*aspect);
664
      *profit->paramlist[PARAM_DISK_ASPECT] = 1.0 / aspect;
665
      *profit->paramlist[PARAM_DISK_POSANG] += 90.0;
666
      }
667 285 bertin
    for (s=0; s<nsub; s++)
668
      {
669
      flux = obj2->prof_disk_flux[s] = *profit->paramlist[PARAM_DISK_FLUX+s];
670
      if (FLAG(obj2.prof_disk_mag))
671
        obj2->prof_disk_mag[s] = flux>0.0f?
672
                 -FDMAG * log(flux) + prefs.mag_zeropoint[s] : 99.0f;
673
      if (FLAG(obj2.prof_disk_fluxerr))
674
        obj2->prof_disk_fluxerr[s] =
675
                profit->paramerr[profit->paramindex[PARAM_DISK_FLUX+s]];
676
      if (FLAG(obj2.prof_disk_magerr))
677
        obj2->prof_disk_magerr[s] = flux>0.0f?
678
                 FDMAG*profit->paramerr[profit->paramindex[PARAM_DISK_FLUX+s]]
679
                / flux : 99.0f;
680
      }
681 115 bertin
    obj2->prof_disk_scale = *profit->paramlist[PARAM_DISK_SCALE];
682 128 bertin
    obj2->prof_disk_scaleerr =
683
                profit->paramerr[profit->paramindex[PARAM_DISK_SCALE]];
684 115 bertin
    obj2->prof_disk_aspect = *profit->paramlist[PARAM_DISK_ASPECT];
685 128 bertin
    obj2->prof_disk_aspecterr =
686
                profit->paramerr[profit->paramindex[PARAM_DISK_ASPECT]];
687 123 bertin
    obj2->prof_disk_theta = fmod_m90_p90(*profit->paramlist[PARAM_DISK_POSANG]);
688 128 bertin
    obj2->prof_disk_thetaerr =
689
                profit->paramerr[profit->paramindex[PARAM_DISK_POSANG]];
690 132 bertin
    if (FLAG(obj2.prof_disk_inclination))
691
      {
692
      obj2->prof_disk_inclination = acos(obj2->prof_disk_aspect) / DEG;
693
      if (FLAG(obj2.prof_disk_inclinationerr))
694
        {
695
        a = sqrt(1.0-obj2->prof_disk_aspect*obj2->prof_disk_aspect);
696
        obj2->prof_disk_inclinationerr = obj2->prof_disk_aspecterr
697
                                        /(a>0.1? a : 0.1)/DEG;
698
        }
699
      }
700 142 bertin
 
701 206 bertin
    if (FLAG(obj2.prof_disk_peak))
702
      {
703
      obj2->prof_disk_peak = obj2->prof_disk_scale>0.0?
704 285 bertin
        obj2->prof_disk_flux[0]
705 206 bertin
        / (2.0 * PI * obj2->prof_disk_scale*obj2->prof_disk_scale
706
                * obj2->prof_disk_aspect)
707
        : 0.0;
708
      if (FLAG(obj2.prof_disk_fluxeff))
709 207 bertin
        obj2->prof_disk_fluxeff = obj2->prof_disk_peak * 0.186682; /* e^-(b_n)*/
710
      if (FLAG(obj2.prof_disk_fluxmean))
711
        obj2->prof_disk_fluxmean = obj2->prof_disk_peak * 0.355007;/* b_n^(-2)*/
712 206 bertin
      }
713
 
714 143 bertin
/* Disk pattern */
715 147 bertin
    if (prefs.pattern_flag)
716 143 bertin
      {
717 258 bertin
      profit_residuals(profit, PROFIT_DYNPARAM, profit->paraminit,profit->resi);
718 171 bertin
      pattern = pattern_init(profit, prefs.pattern_type,
719 147 bertin
                prefs.prof_disk_patternncomp);
720 143 bertin
      pattern_fit(pattern, profit);
721 152 bertin
      if (FLAG(obj2.prof_disk_patternspiral))
722
        obj2->prof_disk_patternspiral = pattern_spiral(pattern);
723 147 bertin
      if (FLAG(obj2.prof_disk_patternvector))
724
        {
725
        ncomp = pattern->size[2];
726
        for (p=0; p<ncomp; p++)
727
          obj2->prof_disk_patternvector[p] = (float)pattern->coeff[p];
728
        }
729
      if (FLAG(obj2.prof_disk_patternmodvector))
730
        {
731
        ncomp = pattern->ncomp*pattern->nfreq;
732
        for (p=0; p<ncomp; p++)
733
          obj2->prof_disk_patternmodvector[p] = (float)pattern->mcoeff[p];
734
        }
735
      if (FLAG(obj2.prof_disk_patternargvector))
736
        {
737
        ncomp = pattern->ncomp*pattern->nfreq;
738
        for (p=0; p<ncomp; p++)
739
          obj2->prof_disk_patternargvector[p] = (float)pattern->acoeff[p];
740
        }
741 143 bertin
      pattern_end(pattern);
742
      }
743 142 bertin
 
744 132 bertin
/* Bar */
745 116 bertin
    if (FLAG(obj2.prof_bar_flux))
746
      {
747
      obj2->prof_bar_flux = *profit->paramlist[PARAM_BAR_FLUX];
748 130 bertin
      obj2->prof_bar_fluxerr =
749
                profit->paramerr[profit->paramindex[PARAM_BAR_FLUX]];
750 116 bertin
      obj2->prof_bar_length = *profit->paramlist[PARAM_ARMS_START]
751
                                **profit->paramlist[PARAM_DISK_SCALE];
752 130 bertin
      obj2->prof_bar_lengtherr = *profit->paramlist[PARAM_ARMS_START]
753
                  * profit->paramerr[profit->paramindex[PARAM_DISK_SCALE]]
754
                + *profit->paramlist[PARAM_DISK_SCALE]
755
                  * profit->paramerr[profit->paramindex[PARAM_ARMS_START]];
756 116 bertin
      obj2->prof_bar_aspect = *profit->paramlist[PARAM_BAR_ASPECT];
757 130 bertin
      obj2->prof_bar_aspecterr =
758
                profit->paramerr[profit->paramindex[PARAM_BAR_ASPECT]];
759 123 bertin
      obj2->prof_bar_posang =
760
                        fmod_m90_p90(*profit->paramlist[PARAM_ARMS_POSANG]);
761 130 bertin
      obj2->prof_bar_posangerr =
762
                profit->paramerr[profit->paramindex[PARAM_ARMS_POSANG]];
763 123 bertin
      if (FLAG(obj2.prof_bar_theta))
764 130 bertin
        {
765
        cp = cos(obj2->prof_bar_posang*DEG);
766
        sp = sin(obj2->prof_bar_posang*DEG);
767
        a = obj2->prof_disk_aspect;
768
        obj2->prof_bar_theta = fmod_m90_p90(atan2(a*sp,cp)/DEG
769
                                + obj2->prof_disk_theta);
770
        obj2->prof_bar_thetaerr = obj2->prof_bar_posangerr*a/(cp*cp+a*a*sp*sp);
771
        }
772 132 bertin
 
773
/* Arms */
774 116 bertin
      if (FLAG(obj2.prof_arms_flux))
775
        {
776
        obj2->prof_arms_flux = *profit->paramlist[PARAM_ARMS_FLUX];
777 130 bertin
        obj2->prof_arms_fluxerr =
778
                profit->paramerr[profit->paramindex[PARAM_ARMS_FLUX]];
779 132 bertin
        obj2->prof_arms_pitch =
780
                fmod_m90_p90(*profit->paramlist[PARAM_ARMS_PITCH]);
781 130 bertin
        obj2->prof_arms_pitcherr =
782
                profit->paramerr[profit->paramindex[PARAM_ARMS_PITCH]];
783 116 bertin
        obj2->prof_arms_start = *profit->paramlist[PARAM_ARMS_START]
784
                                **profit->paramlist[PARAM_DISK_SCALE];
785 130 bertin
        obj2->prof_arms_starterr = *profit->paramlist[PARAM_ARMS_START]
786
                  * profit->paramerr[profit->paramindex[PARAM_DISK_SCALE]]
787
                + *profit->paramlist[PARAM_DISK_SCALE]
788
                  * profit->paramerr[profit->paramindex[PARAM_ARMS_START]];
789 116 bertin
        obj2->prof_arms_quadfrac = *profit->paramlist[PARAM_ARMS_QUADFRAC];
790 130 bertin
        obj2->prof_arms_quadfracerr =
791
                profit->paramerr[profit->paramindex[PARAM_ARMS_QUADFRAC]];
792 123 bertin
        obj2->prof_arms_posang =
793
                        fmod_m90_p90(*profit->paramlist[PARAM_ARMS_POSANG]);
794 130 bertin
        obj2->prof_arms_posangerr =
795
                profit->paramerr[profit->paramindex[PARAM_ARMS_POSANG]];
796 116 bertin
        }
797
      }
798 115 bertin
    }
799
 
800 259 bertin
  obj2->prof_flag = profit->flag;
801
 
802 258 bertin
  return;
803
  }
804
 
805
 
806
/****** profit_spread ********************************************************
807 275 bertin
PROTO   void profit_spread(profitstruct *profit, fieldstruct *field,
808
                fieldstruct *wfield, obj2struct *obj2)
809 258 bertin
PURPOSE Perform star/galaxy separation by comparing the best-fitting local PSF
810
        and a compact exponential profile to the actual data using linear
811
        discriminant analysis.
812
INPUT   Pointer to the profile-fitting structure,
813
        pointer to the field,
814
        pointer to the field weight,
815
        pointer to the obj2.
816
OUTPUT  -.
817
NOTES   -.
818
AUTHOR  E. Bertin (IAP)
819 285 bertin
VERSION 06/05/2012
820 258 bertin
 ***/
821 275 bertin
void    profit_spread(profitstruct *profit,  fieldstruct *field,
822
                fieldstruct *wfield, obj2struct *obj2)
823 258 bertin
  {
824 285 bertin
   profitstruct         *pprofit,*qprofit;
825
   subprofitstruct      *subprofit, *psubprofit,*qsubprofit;
826
   PIXTYPE              valp,valq, sig2;
827
   double               sump,sumq, sumpw2,sumqw2,sumpqw, sump0,sumq0;
828
   float                dchi2;
829
   int                  p,s, nsub;
830 258 bertin
 
831 285 bertin
  nsub = profit->nsubprofit;
832 258 bertin
 
833 285 bertin
  pprofit = profit_init(obj2, MODEL_DIRAC);
834
  qprofit = profit_init(obj2, MODEL_EXPONENTIAL);
835
 
836 258 bertin
  profit_residuals(profit, PROFIT_DYNPARAM, profit->paraminit, profit->resi);
837
  profit_resetparams(pprofit);
838
  if (profit->paramlist[PARAM_X] && profit->paramlist[PARAM_Y])
839 179 bertin
    {
840 258 bertin
    pprofit->paraminit[pprofit->paramindex[PARAM_X]]
841
                = *profit->paramlist[PARAM_X];
842
    pprofit->paraminit[pprofit->paramindex[PARAM_Y]]
843
                = *profit->paramlist[PARAM_Y];
844
    }
845
 
846 285 bertin
  subprofit = profit->subprofit;
847
  for (s=0; s<nsub; s++)
848
    pprofit->paraminit[pprofit->paramindex[PARAM_DIRAC_FLUX+s]]
849
        = (subprofit++)->flux;
850
 
851 258 bertin
  pprofit->niter = profit_minimize(pprofit, PROFIT_MAXITER);
852
 
853
  profit_residuals(pprofit, PROFIT_DYNPARAM,pprofit->paraminit,pprofit->resi);
854
 
855
  qprofit->paraminit[qprofit->paramindex[PARAM_X]]
856 247 bertin
                = pprofit->paraminit[pprofit->paramindex[PARAM_X]];
857 258 bertin
  qprofit->paraminit[qprofit->paramindex[PARAM_Y]]
858 247 bertin
                = pprofit->paraminit[pprofit->paramindex[PARAM_Y]];
859 285 bertin
  for (s=0; s<nsub; s++)
860
    qprofit->paraminit[qprofit->paramindex[PARAM_DISK_FLUX+s]]
861
                = pprofit->paraminit[pprofit->paramindex[PARAM_DIRAC_FLUX+s]];
862 258 bertin
  qprofit->paraminit[qprofit->paramindex[PARAM_DISK_SCALE]]
863 285 bertin
        = profit->subprofit->psf->fwhm/16.0;    /* !CHECK */
864 258 bertin
  qprofit->paraminit[qprofit->paramindex[PARAM_DISK_ASPECT]] = 1.0;
865
  qprofit->paraminit[qprofit->paramindex[PARAM_DISK_POSANG]] = 0.0;
866 245 bertin
 
867 258 bertin
  profit_residuals(qprofit,PROFIT_DYNPARAM, qprofit->paraminit,qprofit->resi);
868
 
869
  if (FLAG(obj2.prof_class_star))
870
    {
871
    dchi2 = 0.5*(pprofit->chi2 - profit->chi2);
872
    obj2->prof_class_star = dchi2 < 50.0?
873
        (dchi2 > -50.0? 2.0/(1.0+expf(dchi2)) : 2.0) : 0.0;
874
    }
875 285 bertin
 
876
  psubprofit = pprofit->subprofit;
877
  qsubprofit = qprofit->subprofit;
878
  for (s=0; s<profit->nsubprofit; s++)
879 258 bertin
    {
880 285 bertin
    sump = sumq = sumpw2 = sumqw2 = sumpqw = sump0 = sumq0 = 0.0;
881
    for (p=0; p<psubprofit->nobjpix; p++)
882
      if (psubprofit->objweight[p]>0 && psubprofit->objpix[p]>-BIG)
883
        {
884
        valp = psubprofit->lmodpix[p];
885
        sump += (double)(valp*psubprofit->objpix[p]);
886
        valq = qsubprofit->lmodpix[p];
887
        sumq += (double)(valq*psubprofit->objpix[p]);
888
        sump0 += (double)(valp*valp);
889
        sumq0 += (double)(valp*valq);
890
        sig2 = 1.0f/(psubprofit->objweight[p]*psubprofit->objweight[p]);
891
        sumpw2 += valp*valp*sig2;
892
        sumqw2 += valq*valq*sig2;
893
        sumpqw += valp*valq*sig2;
894
        }
895
 
896
    if (FLAG(obj2.prof_concentration))
897
      {
898
      obj2->prof_concentration[s] = sump>0.0? (sumq/sump - sumq0/sump0) : 1.0;
899
      if (FLAG(obj2.prof_concentrationerr))
900
        obj2->prof_concentrationerr[s] = sump>0.0?
901 245 bertin
                sqrt(sumqw2*sump*sump+sumpw2*sumq*sumq-2.0*sumpqw*sump*sumq)
902
                        / (sump*sump) : 0.0;
903 285 bertin
      }
904 179 bertin
    }
905
 
906 258 bertin
  profit_end(pprofit);
907
  profit_end(qprofit);
908 208 bertin
 
909 65 bertin
  return;
910
  }
911
 
912 285 bertin
/****** subprofit_noisearea ***************************************************
913
PROTO   float subprofit_noisearea(profitstruct *profit)
914 245 bertin
PURPOSE Return the equivalent noise area (see King 1983) of a model.
915 285 bertin
INPUT   Sub-profile-fitting structure,
916 245 bertin
OUTPUT  Equivalent noise area, in pixels.
917
NOTES   -.
918
AUTHOR  E. Bertin (IAP)
919 285 bertin
VERSION 06/05/2012
920 245 bertin
 ***/
921 285 bertin
float   subprofit_noisearea(subprofitstruct *subprofit)
922 245 bertin
  {
923
   double       dval, flux,flux2;
924
   PIXTYPE      *pix;
925
   int          p;
926
 
927
  flux = flux2 = 0.0;
928 285 bertin
  pix = subprofit->lmodpix;
929
  for (p=subprofit->nobjpix; p--;)
930 245 bertin
    {
931
    dval = (double)*(pix++);
932
    flux += dval;
933
    flux2 += dval*dval;
934
    }
935
 
936
  return (float)(flux2>0.0? flux*flux / flux2 : 0.0);
937
  }
938
 
939
 
940 285 bertin
/****** subprofit_fluxcor ****************************************************
941
PROTO   void subprofit_fluxcor(subprofitstruct *subprofit, obj2struct *obj2)
942 235 bertin
PURPOSE Integrate the flux within an ellipse and complete it with the wings of
943 285 bertin
        the fitted model.
944
INPUT   Sub-profile-fitting structure,
945
        pointer to the obj2 structure.
946
OUTPUT  -.
947 235 bertin
NOTES   -.
948
AUTHOR  E. Bertin (IAP)
949 282 bertin
VERSION 16/03/2012
950 235 bertin
 ***/
951 285 bertin
void    subprofit_fluxcor(subprofitstruct *subprofit, obj2struct *obj2)
952 235 bertin
  {
953 245 bertin
    checkstruct         *check;
954
    double              mx,my, dx,dy, cx2,cy2,cxy, klim,klim2, tvobj,sigtvobj,
955
                        tvm,tvmin,tvmout, r1,v1;
956
    PIXTYPE             *objpix,*objpixt,*objweight,*objweightt, *lmodpix,
957 235 bertin
                        pix, weight,var;
958 245 bertin
    int                 x,y, x2,y2, pos, w,h, area, corrflag;
959 235 bertin
 
960 245 bertin
 
961 235 bertin
  corrflag = (prefs.mask_type==MASK_CORRECT);
962 285 bertin
  w = subprofit->objnaxisn[0];
963
  h = subprofit->objnaxisn[1];
964 235 bertin
  mx = (float)(w/2);
965
  my = (float)(h/2);
966 285 bertin
/*
967 235 bertin
  if (FLAG(obj2.x_prof))
968
    {
969
    if (profit->paramlist[PARAM_X])
970
      mx += *profit->paramlist[PARAM_X];
971
    if (profit->paramlist[PARAM_Y])
972
      my += *profit->paramlist[PARAM_Y];
973
    }
974 285 bertin
*/
975 282 bertin
  if (obj2->auto_kronfactor>0.0)
976 235 bertin
    {
977 258 bertin
    cx2 = obj2->cxx;
978
    cy2 = obj2->cyy;
979
    cxy = obj2->cxy;
980 282 bertin
    klim2 = 0.64*obj2->auto_kronfactor*obj2->auto_kronfactor;
981 235 bertin
    }
982
  else
983
/*-- ...if not, use the circular aperture provided by the user */
984
    {
985
    cx2 = cy2 = 1.0;
986
    cxy = 0.0;
987
    klim2 = (prefs.autoaper[1]/2.0)*(prefs.autoaper[1]/2.0);
988
    }
989 245 bertin
/*
990
  cx2 = obj2->prof_convcxx;
991
  cy2 = obj2->prof_convcyy;
992
  cxy = obj2->prof_convcxy;
993 235 bertin
 
994 245 bertin
  lmodpix = profit->lmodpix;
995
  r1 = v1 = 0.0;
996
  for (y=0; y<h; y++)
997
    {
998
    dy = y - my;
999
    for (x=0; x<w; x++)
1000
      {
1001
      dx = x - mx;
1002
      pix = *(lmodpix++);
1003
      r1 += sqrt(cx2*dx*dx + cy2*dy*dy + cxy*dx*dy)*pix;
1004
      v1 += pix;
1005
      }
1006
    }
1007
 
1008
  klim = r1/v1*2.0;
1009
  klim2 = klim*klim;
1010
 
1011
if ((check = prefs.check[CHECK_APERTURES]))
1012 274 bertin
sexellipse(check->pix, check->width, check->height,
1013 245 bertin
obj2->x_prof-1.0, obj2->y_prof-1.0, klim*obj2->prof_conva,klim*obj2->prof_convb,
1014
obj2->prof_convtheta, check->overlay, 0);
1015
*/
1016
 
1017 235 bertin
  area = 0;
1018
  tvmin = tvmout = tvobj = sigtvobj = 0.0;
1019 285 bertin
  lmodpix = subprofit->lmodpix;
1020
  objpixt = objpix = subprofit->objpix;
1021
  objweightt = objweight = subprofit->objweight;
1022 235 bertin
  for (y=0; y<h; y++)
1023
    {
1024
    for (x=0; x<w; x++, objpixt++,objweightt++)
1025
      {
1026
      dx = x - mx;
1027
      dy = y - my;
1028
      if ((cx2*dx*dx + cy2*dy*dy + cxy*dx*dy) <= klim2)
1029
        {
1030
        area++;
1031
/*------ Here begin tests for pixel and/or weight overflows. Things are a */
1032
/*------ bit intricated to have it running as fast as possible in the most */
1033
/*------ common cases */
1034
        if ((weight=*objweightt)<=0.0)
1035
          {
1036
          if (corrflag
1037
                && (x2=(int)(2*mx+0.49999-x))>=0 && x2<w
1038
                && (y2=(int)(2*my+0.49999-y))>=0 && y2<h
1039
                && (weight=objweight[pos = y2*w + x2])>0.0)
1040
            {
1041
            pix = objpix[pos];
1042
            var = 1.0/(weight*weight);
1043
            }
1044
          else
1045
            pix = var = 0.0;
1046
          }
1047
        else
1048
          {
1049
          pix = *objpixt;
1050
          var = 1.0/(weight*weight);
1051
          }
1052
        tvobj += pix;
1053
        sigtvobj += var;
1054 245 bertin
        tvmin += *(lmodpix++);
1055
//        *(lmodpix++) = pix;
1056 235 bertin
        }
1057
      else
1058
        tvmout += *(lmodpix++);
1059
      }
1060
    }
1061
 
1062
//  tv -= area*bkg;
1063
 
1064
  tvm = tvmin + tvmout;
1065
  if (tvm != 0.0)
1066
    {
1067 285 bertin
    obj2->fluxcor_prof = tvobj+obj2->flux_prof[0]*tvmout/tvm;
1068 235 bertin
    obj2->fluxcorerr_prof = sqrt(sigtvobj
1069 285 bertin
                +obj2->fluxerr_prof[0]*obj2->fluxerr_prof[0]*tvmout/tvm);
1070 235 bertin
    }
1071
  else
1072
    {
1073
    obj2->fluxcor_prof = tvobj;
1074
    obj2->fluxcorerr_prof = sqrt(sigtvobj);
1075
    }
1076
 
1077 245 bertin
/*
1078
  if ((check = prefs.check[CHECK_OTHER]))
1079 275 bertin
    check_add(check, profit->lmodpix, w, h, profit->ix,profit->iy, 1.0);
1080 245 bertin
*/
1081 235 bertin
  return;
1082
  }
1083
 
1084
 
1085 207 bertin
/****i* prof_gammainc *********************************************************
1086
PROTO   double prof_gammainc(double x, double a)
1087 233 bertin
PURPOSE Returns the incomplete Gamma function (based on algorithm described in
1088
        Numerical Recipes in C, chap. 6.1).
1089 207 bertin
INPUT   A double,
1090
        upper integration limit.
1091
OUTPUT  Incomplete Gamma function.
1092
NOTES   -.
1093
AUTHOR  E. Bertin (IAP)
1094 233 bertin
VERSION 08/10/2010
1095 207 bertin
*/
1096
static double   prof_gammainc (double x, double a)
1097
 
1098
  {
1099
   double       b,c,d,h, xn,xp, del,sum;
1100
   int          i;
1101
 
1102
  if (a < 0.0 || x <= 0.0)
1103
    return 0.0;
1104
 
1105
  if (a < (x+1.0))
1106
    {
1107
/*-- Use the series representation */
1108
    xp = x;
1109
    del = sum = 1.0/x;
1110
    for (i=100;i--;)    /* Iterate to convergence */
1111
      {
1112
      sum += (del *= a/(++xp));
1113
      if (fabs(del) < fabs(sum)*3e-7)
1114
        return sum*exp(-a+x*log(a)) / prof_gamma(x);
1115
      }
1116
    }
1117
  else
1118
    {
1119
/*-- Use the continued fraction representation and take its complement */
1120
    b = a + 1.0 - x;
1121
    c = 1e30;
1122
    h = d = 1.0/b;
1123
    for (i=1; i<=100; i++)      /* Iterate to convergence */
1124
      {
1125
      xn = -i*(i-x);
1126
      b += 2.0;
1127
      if (fabs(d=xn*d+b) < 1e-30)
1128
        d = 1e-30;
1129
      if (fabs(c=b+xn/c) < 1e-30)
1130
        c = 1e-30;
1131
      del= c * (d = 1.0/d);
1132
      h *= del;
1133
      if (fabs(del-1.0) < 3e-7)
1134
        return 1.0 - exp(-a+x*log(a))*h / prof_gamma(x);
1135
      }
1136
    }
1137
  error(EXIT_FAILURE, "*Error*: out of bounds in ",
1138
                "prof_gammainc()");
1139
  return 0.0;
1140
  }
1141
 
1142
 
1143 206 bertin
/****i* prof_gamma ************************************************************
1144 207 bertin
PROTO   double prof_gamma(double xx)
1145 233 bertin
PURPOSE Returns the Gamma function (based on algorithm described in Numerical
1146
        Recipes in C, chap 6.1).
1147 206 bertin
INPUT   A double.
1148
OUTPUT  Gamma function.
1149
NOTES   -.
1150
AUTHOR  E. Bertin (IAP)
1151 233 bertin
VERSION 11/09/2009
1152 206 bertin
*/
1153
static double   prof_gamma(double xx)
1154 65 bertin
 
1155 206 bertin
  {
1156
   double               x,tmp,ser;
1157
   static double        cof[6]={76.18009173,-86.50532033,24.01409822,
1158
                        -1.231739516,0.120858003e-2,-0.536382e-5};
1159
   int                  j;
1160
 
1161
  tmp=(x=xx-1.0)+5.5;
1162
  tmp -= (x+0.5)*log(tmp);
1163
  ser=1.0;
1164
  for (j=0;j<6;j++)
1165
    ser += cof[j]/(x+=1.0);
1166
 
1167
  return 2.50662827465*ser*exp(-tmp);
1168
  }
1169
 
1170
 
1171 207 bertin
/****** profit_minradius ******************************************************
1172
PROTO   float profit_minradius(profitstruct *profit, float refffac)
1173
PURPOSE Returns the minimum disk radius that guarantees that each and
1174
        every model component fits within some margin in that disk.
1175
INPUT   Profit structure pointer,
1176
        margin in units of (r/r_eff)^(1/n)).
1177
OUTPUT  Radius (in pixels).
1178
NOTES   -.
1179
AUTHOR  E. Bertin (IAP)
1180 233 bertin
VERSION 08/10/2010
1181 207 bertin
*/
1182
float   profit_minradius(profitstruct *profit, float refffac)
1183
 
1184
  {
1185
   double       r,reff,rmax;
1186
   int          p;
1187
 
1188
  rmax = reff = 0.0;
1189
  for (p=0; p<profit->nprof; p++)
1190
    {
1191
    switch (profit->prof[p]->code)
1192
      {
1193 245 bertin
      case MODEL_BACK:
1194
      case MODEL_DIRAC:
1195 233 bertin
        reff = 0.0;
1196
      break;
1197 245 bertin
      case MODEL_SERSIC:
1198 207 bertin
        reff = *profit->paramlist[PARAM_SPHEROID_REFF];
1199 233 bertin
        break;
1200 245 bertin
      case MODEL_DEVAUCOULEURS:
1201 207 bertin
        reff = *profit->paramlist[PARAM_SPHEROID_REFF];
1202 233 bertin
        break;
1203 245 bertin
      case MODEL_EXPONENTIAL:
1204 207 bertin
        reff = *profit->paramlist[PARAM_DISK_SCALE]*1.67835;
1205 233 bertin
        break;
1206 207 bertin
      default:
1207
        error(EXIT_FAILURE, "*Internal Error*: Unknown profile parameter in ",
1208
                "profit_minradius()");
1209 233 bertin
        break;
1210 207 bertin
      }
1211
    r = reff*(double)refffac;
1212
    if (r>rmax)
1213
      rmax = r;
1214
    }
1215
 
1216
  return (float)rmax;
1217
  }
1218
 
1219
 
1220 285 bertin
/****** subprofit_psf *********************************************************
1221
PROTO   void    subprofit_psf(subprofitstruct *subprofit, obj2struct *obj2)
1222 118 bertin
PURPOSE Build the local PSF at a given resolution.
1223 285 bertin
INPUT   Pointer to sub-profile-fitting structure,
1224
        pointer to obj2.
1225 118 bertin
OUTPUT  -.
1226
NOTES   -.
1227
AUTHOR  E. Bertin (IAP)
1228 285 bertin
VERSION 06/05/2012
1229 118 bertin
 ***/
1230 285 bertin
void    subprofit_psf(subprofitstruct *subprofit, obj2struct *obj2)
1231 118 bertin
  {
1232 285 bertin
   psfstruct    *psf;
1233 201 bertin
   double       flux;
1234
   float        posin[2], posout[2], dnaxisn[2],
1235 118 bertin
                *pixout,
1236 201 bertin
                xcout,ycout, xcin,ycin, invpixstep, norm;
1237 118 bertin
   int          d,i;
1238
 
1239 285 bertin
  psf = subprofit->psf;
1240
  psf_build(psf, obj2);
1241 118 bertin
 
1242 285 bertin
  xcout = (float)(subprofit->modnaxisn[0]/2) + 1.0;      /* FITS convention */
1243
  ycout = (float)(subprofit->modnaxisn[1]/2) + 1.0;     /* FITS convention */
1244 118 bertin
  xcin = (psf->masksize[0]/2) + 1.0;                     /* FITS convention */
1245
  ycin = (psf->masksize[1]/2) + 1.0;                    /* FITS convention */
1246 285 bertin
  invpixstep = subprofit->pixstep / psf->pixstep;
1247 118 bertin
 
1248
/* Initialize multi-dimensional counters */
1249
  for (d=0; d<2; d++)
1250
    {
1251
    posout[d] = 1.0;                                    /* FITS convention */
1252 285 bertin
    dnaxisn[d] = subprofit->modnaxisn[d]+0.5;
1253 118 bertin
    }
1254
 
1255
/* Remap each pixel */
1256 285 bertin
  pixout = subprofit->psfpix;
1257 118 bertin
  flux = 0.0;
1258 285 bertin
  for (i=subprofit->modnaxisn[0]*subprofit->modnaxisn[1]; i--;)
1259 118 bertin
    {
1260
    posin[0] = (posout[0] - xcout)*invpixstep + xcin;
1261
    posin[1] = (posout[1] - ycout)*invpixstep + ycin;
1262
    flux += ((*(pixout++) = interpolate_pix(posin, psf->maskloc,
1263
                psf->masksize, INTERP_LANCZOS3)));
1264
    for (d=0; d<2; d++)
1265
      if ((posout[d]+=1.0) < dnaxisn[d])
1266
        break;
1267
      else
1268
        posout[d] = 1.0;
1269
    }
1270
 
1271 164 bertin
/* Normalize PSF flux (just in case...) */
1272 285 bertin
  flux *= subprofit->pixstep*subprofit->pixstep;
1273 221 bertin
  if (fabs(flux) <= 0.0)
1274
    error(EXIT_FAILURE, "*Error*: PSF model is empty or negative: ", psf->name);
1275 164 bertin
 
1276 221 bertin
  norm = 1.0/flux;
1277 285 bertin
  pixout = subprofit->psfpix;
1278
  for (i=subprofit->modnaxisn[0]*subprofit->modnaxisn[1]; i--;)
1279 221 bertin
    *(pixout++) *= norm;
1280
 
1281 118 bertin
  return;
1282
  }
1283
 
1284
 
1285 65 bertin
/****** profit_minimize *******************************************************
1286
PROTO   void profit_minimize(profitstruct *profit)
1287
PURPOSE Provide a function returning residuals to lmfit.
1288
INPUT   Pointer to the profit structure involved in the fit,
1289
        maximum number of iterations.
1290
OUTPUT  Number of iterations used.
1291
NOTES   -.
1292
AUTHOR  E. Bertin (IAP)
1293 247 bertin
VERSION 20/05/2011
1294 65 bertin
 ***/
1295
int     profit_minimize(profitstruct *profit, int niter)
1296
  {
1297 247 bertin
   double       lm_opts[5], info[LM_INFO_SZ],
1298
                dcovar[PARAM_NPARAM*PARAM_NPARAM], dparam[PARAM_NPARAM];
1299
   int          nfree;
1300 65 bertin
 
1301 221 bertin
  profit->iter = 0;
1302
  memset(dcovar, 0, profit->nparam*profit->nparam*sizeof(double));
1303 65 bertin
 
1304 72 bertin
/* Perform fit */
1305 245 bertin
  lm_opts[0] = 1.0e-3;           /* Initial mu */
1306 269 bertin
  lm_opts[1] = 1.0e-8;          /* ||J^T e||_inf stopping factor */
1307
  lm_opts[2] = 1.0e-8;          /* |Dp||_2 stopping factor */
1308
  lm_opts[3] = 1.0e-8;          /* ||e||_2 stopping factor */
1309 245 bertin
  lm_opts[4] = 1.0e-4;          /* Jacobian step */
1310 65 bertin
 
1311 247 bertin
  nfree = profit_boundtounbound(profit, profit->paraminit, dparam,
1312
                                PARAM_ALLPARAMS);
1313 65 bertin
 
1314 247 bertin
  niter = dlevmar_dif(profit_evaluate, dparam, NULL, nfree, profit->nresi,
1315
                        niter, lm_opts, info, NULL, dcovar, profit);
1316 65 bertin
 
1317 221 bertin
  profit_unboundtobound(profit, dparam, profit->paraminit, PARAM_ALLPARAMS);
1318 52 bertin
 
1319 221 bertin
/* Convert covariance matrix to bounded space */
1320
  profit_covarunboundtobound(profit, dcovar, profit->covar);
1321
 
1322 72 bertin
  return niter;
1323 46 bertin
  }
1324
 
1325 48 bertin
 
1326
/****** profit_printout *******************************************************
1327 201 bertin
PROTO   void profit_printout(int n_par, float* par, int m_dat, float* fvec,
1328 48 bertin
                void *data, int iflag, int iter, int nfev )
1329
PURPOSE Provide a function to print out results to lmfit.
1330
INPUT   Number of fitted parameters,
1331
        pointer to the vector of parameters,
1332
        number of data points,
1333
        pointer to the vector of residuals (output),
1334
        pointer to the data structure (unused),
1335
 
1336
        outer loop counter,
1337
        number of calls to evaluate().
1338
OUTPUT  -.
1339
NOTES   Input arguments are there only for compatibility purposes (unused)
1340
AUTHOR  E. Bertin (IAP)
1341 285 bertin
VERSION 06/05/2012
1342 48 bertin
 ***/
1343 201 bertin
void    profit_printout(int n_par, float* par, int m_dat, float* fvec,
1344 48 bertin
                void *data, int iflag, int iter, int nfev )
1345
  {
1346 275 bertin
   fieldstruct  *field;
1347 77 bertin
   checkstruct  *check;
1348
   profitstruct *profit;
1349
   char         filename[256];
1350
   static int   itero;
1351
 
1352
  profit = (profitstruct *)data;
1353
 
1354 116 bertin
  if (0 && (iter!=itero || iter<0))
1355 77 bertin
    {
1356
    if (iter<0)
1357
      itero++;
1358
    else
1359
      itero = iter;
1360 258 bertin
    field = NULL;       /*! Should be replaced with a global variable to work!*/
1361 77 bertin
    sprintf(filename, "check_%d_%04d.fits", the_gal, itero);
1362 278 bertin
    check=check_init(filename, CHECK_PROFILES, 0, 1);
1363 275 bertin
    check_reinit(field, check);
1364 285 bertin
    check_add(check, profit->subprofit->lmodpix,
1365
                profit->subprofit->objnaxisn[0],profit->subprofit->objnaxisn[1],
1366
                profit->subprofit->ix,profit->subprofit->iy, 1.0);
1367 66 bertin
 
1368 275 bertin
    check_reend(field, check);
1369
    check_end(check);
1370 77 bertin
    }
1371
 
1372 48 bertin
  return;
1373
  }
1374
 
1375
 
1376 102 bertin
/****** profit_evaluate ******************************************************
1377 221 bertin
PROTO   void profit_evaluate(double *par, double *fvec, int m, int n,
1378
                        void *adata)
1379 73 bertin
PURPOSE Provide a function returning residuals to levmar.
1380 57 bertin
INPUT   Pointer to the vector of parameters,
1381 73 bertin
        pointer to the vector of residuals (output),
1382
        number of model parameters,
1383 57 bertin
        number of data points,
1384 221 bertin
        pointer to a data structure (we use it for the profit structure here).
1385 57 bertin
OUTPUT  -.
1386 221 bertin
NOTES   -.
1387 57 bertin
AUTHOR  E. Bertin (IAP)
1388 258 bertin
VERSION 26/07/2011
1389 57 bertin
 ***/
1390 221 bertin
void    profit_evaluate(double *dpar, double *fvec, int m, int n, void *adata)
1391 57 bertin
  {
1392 221 bertin
   profitstruct         *profit;
1393
   profstruct           **prof;
1394
   double               *dpar0, *dresi;
1395
   float                *modpixt, *profpixt, *resi,
1396
                        tparam, val;
1397
   PIXTYPE              *lmodpixt,*lmodpix2t, *objpix,*weight,
1398
                        wval;
1399
   int                  c,f,i,p,q, fd,pd, jflag,sflag, nprof;
1400 57 bertin
 
1401
  profit = (profitstruct *)adata;
1402 221 bertin
 
1403
/* Detect "Jacobian-related" calls */
1404
  jflag = pd = fd = 0;
1405
  dpar0 = profit->dparam;
1406
  if (profit->iter)
1407
    {
1408
    f = q = 0;
1409
    for (p=0; p<profit->nparam; p++)
1410
      {
1411
      if (dpar[f] - dpar0[f] != 0.0)
1412
        {
1413
        pd = p;
1414
        fd = f;
1415
        q++;
1416
        }
1417 247 bertin
      if (profit->parfittype[p]!=PARFIT_FIXED)
1418 221 bertin
        f++;
1419
      }
1420
    if (f>0 && q==1)
1421
      jflag = 1;
1422
    }
1423 245 bertin
jflag = 0;       /* Temporarily deactivated (until problems are fixed) */
1424 285 bertin
/*
1425 245 bertin
  if (jflag && !(profit->nprof==1 && profit->prof[0]->code == MODEL_DIRAC))
1426 221 bertin
    {
1427
    prof = profit->prof;
1428
    nprof = profit->nprof;
1429
 
1430 285 bertin
*-- "Jacobian call" *
1431 221 bertin
    tparam = profit->param[pd];
1432
    profit_unboundtobound(profit, &dpar[fd], &profit->param[pd], pd);
1433
    sflag = 1;
1434
    switch(profit->paramrevindex[pd])
1435
      {
1436 233 bertin
      case PARAM_BACK:
1437
        lmodpixt = profit->lmodpix;
1438
        lmodpix2t = profit->lmodpix2;
1439
        val = (profit->param[pd] - tparam);
1440
        for (i=profit->nobjpix;i--;)
1441
          *(lmodpix2t++) = val;
1442
        break;
1443 221 bertin
      case PARAM_X:
1444
      case PARAM_Y:
1445
        profit_resample(profit, profit->cmodpix, profit->lmodpix2, 1.0);
1446
        lmodpixt = profit->lmodpix;
1447
        lmodpix2t = profit->lmodpix2;
1448
        for (i=profit->nobjpix;i--;)
1449
          *(lmodpix2t++) -= *(lmodpixt++);
1450
        break;
1451 233 bertin
      case PARAM_DIRAC_FLUX:
1452 221 bertin
      case PARAM_SPHEROID_FLUX:
1453
      case PARAM_DISK_FLUX:
1454
      case PARAM_ARMS_FLUX:
1455
      case PARAM_BAR_FLUX:
1456
        if (nprof==1 && tparam != 0.0)
1457
          {
1458
          lmodpixt = profit->lmodpix;
1459
          lmodpix2t = profit->lmodpix2;
1460
          val = (profit->param[pd] - tparam) / tparam;
1461
          for (i=profit->nobjpix;i--;)
1462
            *(lmodpix2t++) = val**(lmodpixt++);
1463
          }
1464
        else
1465
          {
1466
          for (c=0; c<nprof; c++)
1467 285 bertin
            if (prof[c]->flux[0] == &profit->param[pd])
1468 221 bertin
              break;
1469
          memcpy(profit->modpix, prof[c]->pix, profit->nmodpix*sizeof(float));
1470
          profit_convolve(profit, profit->modpix);
1471
          profit_resample(profit, profit->modpix, profit->lmodpix2,
1472
                profit->param[pd] - tparam);
1473
          }
1474
        break;
1475
      case PARAM_SPHEROID_REFF:
1476
      case PARAM_SPHEROID_ASPECT:
1477
      case PARAM_SPHEROID_POSANG:
1478
      case PARAM_SPHEROID_SERSICN:
1479 285 bertin
        sflag = 0;                      * We are in the same switch *
1480 221 bertin
        for (c=0; c<nprof; c++)
1481 245 bertin
          if (prof[c]->code == MODEL_SERSIC
1482
                || prof[c]->code == MODEL_DEVAUCOULEURS)
1483 221 bertin
            break;
1484
      case PARAM_DISK_SCALE:
1485
      case PARAM_DISK_ASPECT:
1486
      case PARAM_DISK_POSANG:
1487
        if (sflag)
1488
          for (c=0; c<nprof; c++)
1489 245 bertin
            if (prof[c]->code == MODEL_EXPONENTIAL)
1490 221 bertin
              break;
1491
        sflag = 0;
1492
      case PARAM_ARMS_QUADFRAC:
1493
      case PARAM_ARMS_SCALE:
1494
      case PARAM_ARMS_START:
1495
      case PARAM_ARMS_POSANG:
1496
      case PARAM_ARMS_PITCH:
1497
      case PARAM_ARMS_PITCHVAR:
1498
      case PARAM_ARMS_WIDTH:
1499
        if (sflag)
1500
          for (c=0; c<nprof; c++)
1501 245 bertin
            if (prof[c]->code == MODEL_ARMS)
1502 221 bertin
              break;
1503
        sflag = 0;
1504
      case PARAM_BAR_ASPECT:
1505
      case PARAM_BAR_POSANG:
1506
        if (sflag)
1507
          for (c=0; c<nprof; c++)
1508 245 bertin
            if (prof[c]->code == MODEL_ARMS)
1509 221 bertin
              break;
1510
        modpixt = profit->modpix;
1511
        profpixt = prof[c]->pix;
1512 285 bertin
        val = -*prof[c]->flux[0];
1513 221 bertin
        for (i=profit->nmodpix;i--;)
1514
          *(modpixt++) = val**(profpixt++);
1515
        memcpy(profit->modpix2, prof[c]->pix, profit->nmodpix*sizeof(float));
1516
        prof_add(profit, prof[c], 0);
1517
        memcpy(prof[c]->pix, profit->modpix2, profit->nmodpix*sizeof(float));
1518
        profit_convolve(profit, profit->modpix);
1519
        profit_resample(profit, profit->modpix, profit->lmodpix2, 1.0);
1520
        break;
1521
      default:
1522
        error(EXIT_FAILURE, "*Internal Error*: ",
1523
                        "unknown parameter index in profit_jacobian()");
1524
        break;
1525
      }
1526
    objpix = profit->objpix;
1527
    weight = profit->objweight;
1528
    lmodpixt = profit->lmodpix;
1529
    lmodpix2t = profit->lmodpix2;
1530
    resi = profit->resi;
1531
    dresi = fvec;
1532
    if (PROFIT_DYNPARAM > 0.0)
1533
      for (i=profit->nobjpix;i--; lmodpixt++, lmodpix2t++)
1534
        {
1535
        val = *(objpix++);
1536
        if ((wval=*(weight++))>0.0)
1537
          *(dresi++) = *(resi++) + *lmodpix2t
1538
                * wval/(1.0+wval*fabs(*lmodpixt - val)/PROFIT_DYNPARAM);
1539
        }
1540
    else
1541
      for (i=profit->nobjpix;i--; lmodpix2t++)
1542
        if ((wval=*(weight++))>0.0)
1543
          *(dresi++) = *(resi++) + *lmodpix2t * wval;
1544
    }
1545
  else
1546 285 bertin
*/
1547 221 bertin
    {
1548
/*-- "Regular call" */
1549
    for (p=0; p<profit->nparam; p++)
1550
      dpar0[p] = dpar[p];
1551
    profit_unboundtobound(profit, dpar, profit->param, PARAM_ALLPARAMS);
1552
 
1553 258 bertin
    profit_residuals(profit, PROFIT_DYNPARAM, profit->param, profit->resi);
1554 221 bertin
 
1555
    for (p=0; p<profit->nresi; p++)
1556
      fvec[p] = profit->resi[p];
1557
    }
1558
 
1559
//  profit_printout(m, par, n, fvec, adata, 0, -1, 0 );
1560
  profit->iter++;
1561
 
1562 57 bertin
  return;
1563
  }
1564
 
1565
 
1566 46 bertin
/****** profit_residuals ******************************************************
1567 258 bertin
PROTO   float *profit_residuals(profitstruct *profit, float dynparam,
1568
                        float *param, float *resi)
1569 46 bertin
PURPOSE Compute the vector of residuals between the data and the galaxy
1570
        profile model.
1571
INPUT   Profile-fitting structure,
1572 180 bertin
        dynamic compression parameter (0=no compression),
1573 51 bertin
        pointer to the model parameters (output),
1574
        pointer to the computed residuals (output).
1575 42 bertin
OUTPUT  Vector of residuals.
1576 117 bertin
NOTES   -.
1577 42 bertin
AUTHOR  E. Bertin (IAP)
1578 285 bertin
VERSION 06/05/2012
1579 42 bertin
 ***/
1580 258 bertin
float   *profit_residuals(profitstruct *profit, float dynparam, float *param,
1581
                float *resi)
1582 42 bertin
  {
1583 285 bertin
   subprofitstruct      *subprofit;
1584
   float                *pixin, *pixout,
1585
                        pflux;
1586
   int                  i,p,s;
1587 42 bertin
 
1588 285 bertin
  subprofit = profit->subprofit;
1589
  for (s=0; s<profit->nsubprofit; s++, subprofit++)
1590
    memset(subprofit->modpix, 0, subprofit->nmodpix*sizeof(float));
1591 51 bertin
  for (p=0; p<profit->nparam; p++)
1592
    profit->param[p] = param[p];
1593 179 bertin
/* Simple PSF shortcut */
1594 245 bertin
  if (profit->nprof == 1 && profit->prof[0]->code == MODEL_DIRAC)
1595 208 bertin
    {
1596 285 bertin
    subprofit = profit->subprofit;
1597
    for (s=0; s<profit->nsubprofit; s++, subprofit++)
1598
      {
1599
      profit_resample(profit, subprofit, subprofit->psfpix,
1600
                subprofit->lmodpix, *profit->prof[0]->flux[s]);
1601
      subprofit->flux = *profit->prof[0]->flux[s];
1602
      }
1603 208 bertin
    }
1604 179 bertin
  else
1605
    {
1606 285 bertin
    subprofit = profit->subprofit;
1607
    for (s=0; s<profit->nsubprofit; s++, subprofit++)
1608
      {
1609
      subprofit->flux = 0.0;
1610
      for (p=0; p<profit->nprof; p++)
1611
if(!s)
1612
        subprofit->flux += prof_add(subprofit, profit->prof[p], 0);
1613
else
1614
{
1615
pflux = *profit->prof[p]->flux[s];
1616
pixin = profit->prof[p]->pix;
1617
pixout = subprofit->modpix;
1618
for (i=subprofit->nmodpix; i--;)
1619
  *(pixout++) += pflux**(pixin++);
1620
subprofit->flux += pflux;
1621
        }
1622
      memcpy(subprofit->cmodpix, subprofit->modpix,
1623
                subprofit->nmodpix*sizeof(float));
1624
      subprofit_convolve(subprofit, subprofit->cmodpix);
1625
      profit_resample(profit, subprofit,
1626
        subprofit->cmodpix, subprofit->lmodpix, 1.0);
1627
      }
1628 179 bertin
    }
1629 42 bertin
 
1630 221 bertin
  if (resi)
1631
    profit_compresi(profit, dynparam, resi);
1632
 
1633 51 bertin
  return resi;
1634 42 bertin
  }
1635
 
1636
 
1637 47 bertin
/****** profit_compresi ******************************************************
1638 285 bertin
PROTO   float *profit_compresi(profitstruct *profit,float dynparam,float *resi)
1639 47 bertin
PURPOSE Compute the vector of residuals between the data and the galaxy
1640
        profile model.
1641
INPUT   Profile-fitting structure,
1642 180 bertin
        dynamic-compression parameter (0=no compression),
1643 51 bertin
        vector of residuals (output).
1644 47 bertin
OUTPUT  Vector of residuals.
1645 117 bertin
NOTES   -.
1646 47 bertin
AUTHOR  E. Bertin (IAP)
1647 285 bertin
VERSION 06/05/2012
1648 47 bertin
 ***/
1649 201 bertin
float   *profit_compresi(profitstruct *profit, float dynparam, float *resi)
1650 47 bertin
  {
1651 285 bertin
   subprofitstruct      *subprofit;
1652
   double               error;
1653
   float                *resit;
1654
   PIXTYPE              *objpix, *objweight, *lmodpix,
1655
                        val,val2,wval, invsig;
1656
   int                  i,s;
1657 50 bertin
 
1658
/* Compute vector of residuals */
1659 51 bertin
  resit = resi;
1660 79 bertin
  error = 0.0;
1661 285 bertin
  subprofit = profit->subprofit;
1662 180 bertin
  if (dynparam > 0.0)
1663 88 bertin
    {
1664 180 bertin
    invsig = (PIXTYPE)(1.0/dynparam);
1665 285 bertin
    for (s=profit->nsubprofit; s--; subprofit++)
1666 117 bertin
      {
1667 285 bertin
      objpix = subprofit->objpix;
1668
      objweight = subprofit->objweight;
1669
      lmodpix = subprofit->lmodpix;
1670
      for (i=subprofit->objnaxisn[0]*subprofit->objnaxisn[1]; i--; lmodpix++)
1671 88 bertin
        {
1672 285 bertin
        val = *(objpix++);
1673
        if ((wval=*(objweight++))>0.0)
1674
          {
1675
          val2 = (*lmodpix - val)*wval*invsig;
1676
          val2 = val2>0.0? logf(1.0+val2) : -logf(1.0-val2);
1677
          *(resit++) = val2*dynparam;
1678
          error += val2*val2;
1679
          }
1680 88 bertin
        }
1681 117 bertin
      }
1682 180 bertin
    profit->chi2 = dynparam*dynparam*error;
1683 88 bertin
    }
1684 180 bertin
  else
1685
    {
1686 285 bertin
    for (s=profit->nsubprofit; s--; subprofit++)
1687 180 bertin
      {
1688 285 bertin
      objpix = subprofit->objpix;
1689
      objweight = subprofit->objweight;
1690
      lmodpix = subprofit->lmodpix;
1691
      for (i=subprofit->objnaxisn[0]*subprofit->objnaxisn[1]; i--; lmodpix++)
1692 180 bertin
        {
1693 285 bertin
        val = *(objpix++);
1694
        if ((wval=*(objweight++))>0.0)
1695
          {
1696
          val2 = (*lmodpix - val)*wval;
1697
          *(resit++) = val2;
1698
          error += val2*val2;
1699
          }
1700 180 bertin
        }
1701
      }
1702
    profit->chi2 = error;
1703
    }
1704 50 bertin
 
1705 51 bertin
  return resi;
1706 50 bertin
  }
1707
 
1708
 
1709
/****** profit_resample ******************************************************
1710 285 bertin
PROTO   int     prof_resample(profitstruct *profit, subprofitstruct *subprofit,
1711
                float *inpix, PIXTYPE *outpix, float factor)
1712 50 bertin
PURPOSE Resample the current full resolution model to image resolution.
1713 221 bertin
INPUT   Profile-fitting structure,
1714 285 bertin
        sub-profile-fitting structure,
1715 221 bertin
        pointer to input raster,
1716
        pointer to output raster,
1717
        multiplicating factor.
1718
OUTPUT  RETURN_ERROR if the rasters don't overlap, RETURN_OK otherwise.
1719 117 bertin
NOTES   -.
1720 50 bertin
AUTHOR  E. Bertin (IAP)
1721 285 bertin
VERSION 05/05/2012
1722 50 bertin
 ***/
1723 285 bertin
int     profit_resample(profitstruct *profit, subprofitstruct *subprofit,
1724
                float *inpix, PIXTYPE *outpix, float factor)
1725 50 bertin
  {
1726 221 bertin
   PIXTYPE      *pixout,*pixout0;
1727
   float        *pixin,*pixin0, *mask,*maskt, *pixinout, *dpixin,*dpixin0,
1728
                *dpixout,*dpixout0, *dx,*dy,
1729
                xcin,xcout,ycin,ycout, xsin,ysin, xin,yin, x,y, dxm,dym, val,
1730 229 bertin
                invpixstep, norm;
1731 221 bertin
   int          *start,*startt, *nmask,*nmaskt,
1732
                i,j,k,n,t,
1733
                ixsout,iysout, ixout,iyout, dixout,diyout, nxout,nyout,
1734
                iysina, nyin, hmw,hmh, ix,iy, ixin,iyin;
1735 51 bertin
 
1736 285 bertin
  invpixstep = subprofit->subsamp/subprofit->pixstep;
1737 221 bertin
 
1738 285 bertin
  xcin = (subprofit->modnaxisn[0]/2);
1739
  xcout = (float)(subprofit->objnaxisn[0]/2);
1740 71 bertin
  if ((dx=(profit->paramlist[PARAM_X])))
1741 74 bertin
    xcout += *dx;
1742 221 bertin
 
1743
  xsin = xcin - xcout*invpixstep;                       /* Input start x-coord*/
1744 285 bertin
  if ((int)xsin >= subprofit->modnaxisn[0])
1745 221 bertin
    return RETURN_ERROR;
1746
  ixsout = 0;                            /* Int. part of output start x-coord */
1747
  if (xsin<0.0)
1748
    {
1749
    dixout = (int)(1.0-xsin/invpixstep);
1750
/*-- Simply leave here if the images do not overlap in x */
1751 285 bertin
    if (dixout >= subprofit->objnaxisn[0])
1752 221 bertin
      return RETURN_ERROR;
1753
    ixsout += dixout;
1754
    xsin += dixout*invpixstep;
1755
    }
1756 285 bertin
  nxout = (int)((subprofit->modnaxisn[0]-xsin)/invpixstep);/* nb of interpolated
1757 221 bertin
                                                        input pixels along x */
1758 285 bertin
  if (nxout>(ixout=subprofit->objnaxisn[0]-ixsout))
1759 221 bertin
    nxout = ixout;
1760
  if (!nxout)
1761
    return RETURN_ERROR;
1762
 
1763 285 bertin
  ycin = (subprofit->modnaxisn[1]/2);
1764
  ycout = (float)(subprofit->objnaxisn[1]/2);
1765 71 bertin
  if ((dy=(profit->paramlist[PARAM_Y])))
1766 74 bertin
    ycout += *dy;
1767 57 bertin
 
1768 221 bertin
  ysin = ycin - ycout*invpixstep;               /* Input start y-coord*/
1769 285 bertin
  if ((int)ysin >= subprofit->modnaxisn[1])
1770 221 bertin
    return RETURN_ERROR;
1771
  iysout = 0;                            /* Int. part of output start y-coord */
1772
  if (ysin<0.0)
1773 47 bertin
    {
1774 221 bertin
    diyout = (int)(1.0-ysin/invpixstep);
1775
/*-- Simply leave here if the images do not overlap in y */
1776 285 bertin
    if (diyout >= subprofit->objnaxisn[1])
1777 221 bertin
      return RETURN_ERROR;
1778
    iysout += diyout;
1779
    ysin += diyout*invpixstep;
1780 47 bertin
    }
1781 285 bertin
  nyout = (int)((subprofit->modnaxisn[1]-ysin)/invpixstep);/* nb of interpolated
1782 221 bertin
                                                        input pixels along y */
1783 285 bertin
  if (nyout>(iyout=subprofit->objnaxisn[1]-iysout))
1784 221 bertin
    nyout = iyout;
1785
  if (!nyout)
1786
    return RETURN_ERROR;
1787 47 bertin
 
1788 221 bertin
/* Set the yrange for the x-resampling with some margin for interpolation */
1789
  iysina = (int)ysin;   /* Int. part of Input start y-coord with margin */
1790 229 bertin
  hmh = INTERPW/2 - 1;  /* Interpolant start */
1791 221 bertin
  if (iysina<0 || ((iysina -= hmh)< 0))
1792
    iysina = 0;
1793 229 bertin
  nyin = (int)(ysin+nyout*invpixstep)+INTERPW-hmh;/* Interpolated Input y size*/
1794 285 bertin
  if (nyin>subprofit->modnaxisn[1])                                             /* with margin */
1795
    nyin = subprofit->modnaxisn[1];
1796 221 bertin
/* Express everything relative to the effective Input start (with margin) */
1797
  nyin -= iysina;
1798
  ysin -= (float)iysina;
1799
 
1800
/* Allocate interpolant stuff for the x direction */
1801
  QMALLOC(mask, float, nxout*INTERPW);  /* Interpolation masks */
1802
  QMALLOC(nmask, int, nxout);           /* Interpolation mask sizes */
1803
  QMALLOC(start, int, nxout);           /* Int. part of Input conv starts */
1804
/* Compute the local interpolant and data starting points in x */
1805
  hmw = INTERPW/2 - 1;
1806
  xin = xsin;
1807
  maskt = mask;
1808
  nmaskt = nmask;
1809
  startt = start;
1810
  for (j=nxout; j--; xin+=invpixstep)
1811 47 bertin
    {
1812 221 bertin
    ix = (ixin=(int)xin) - hmw;
1813
    dxm = ixin - xin - hmw;     /* starting point in the interpolation func */
1814
    if (ix < 0)
1815 212 bertin
      {
1816 221 bertin
      n = INTERPW+ix;
1817
      dxm -= (float)ix;
1818
      ix = 0;
1819 212 bertin
      }
1820 221 bertin
    else
1821
      n = INTERPW;
1822 285 bertin
    if (n>(t=subprofit->modnaxisn[0]-ix))
1823 221 bertin
      n=t;
1824
    *(startt++) = ix;
1825
    *(nmaskt++) = n;
1826 229 bertin
    norm = 0.0;
1827 221 bertin
    for (x=dxm, i=n; i--; x+=1.0)
1828 229 bertin
      norm += (*(maskt++) = INTERPF(x));
1829
    norm = norm>0.0? 1.0/norm : 1.0;
1830
    maskt -= n;
1831
    for (i=n; i--;)
1832
      *(maskt++) *= norm;
1833 212 bertin
    }
1834 221 bertin
 
1835
  QCALLOC(pixinout, float, nxout*nyin); /* Intermediary frame-buffer */
1836
 
1837
/* Make the interpolation in x (this includes transposition) */
1838 285 bertin
  pixin0 = inpix + iysina*subprofit->modnaxisn[0];
1839 221 bertin
  dpixout0 = pixinout;
1840 285 bertin
  for (k=nyin; k--; pixin0+=subprofit->modnaxisn[0], dpixout0++)
1841 221 bertin
    {
1842
    maskt = mask;
1843
    nmaskt = nmask;
1844
    startt = start;
1845
    dpixout = dpixout0;
1846
    for (j=nxout; j--; dpixout+=nyin)
1847 212 bertin
      {
1848 221 bertin
      pixin = pixin0+*(startt++);
1849
      val = 0.0;
1850
      for (i=*(nmaskt++); i--;)
1851
        val += *(maskt++)**(pixin++);
1852
      *dpixout = val;
1853 212 bertin
      }
1854 221 bertin
    }
1855 47 bertin
 
1856 221 bertin
/* Reallocate interpolant stuff for the y direction */
1857 229 bertin
  QREALLOC(mask, float, nyout*INTERPW); /* Interpolation masks */
1858 221 bertin
  QREALLOC(nmask, int, nyout);                  /* Interpolation mask sizes */
1859
  QREALLOC(start, int, nyout);          /* Int. part of Input conv starts */
1860
 
1861
/* Compute the local interpolant and data starting points in y */
1862 229 bertin
  hmh = INTERPW/2 - 1;
1863 221 bertin
  yin = ysin;
1864
  maskt = mask;
1865
  nmaskt = nmask;
1866
  startt = start;
1867
  for (j=nyout; j--; yin+=invpixstep)
1868
    {
1869
    iy = (iyin=(int)yin) - hmh;
1870
    dym = iyin - yin - hmh;     /* starting point in the interpolation func */
1871
    if (iy < 0)
1872
      {
1873 229 bertin
      n = INTERPW+iy;
1874 221 bertin
      dym -= (float)iy;
1875
      iy = 0;
1876
      }
1877
    else
1878 229 bertin
      n = INTERPW;
1879 221 bertin
    if (n>(t=nyin-iy))
1880
      n=t;
1881
    *(startt++) = iy;
1882
    *(nmaskt++) = n;
1883 229 bertin
    norm = 0.0;
1884 221 bertin
    for (y=dym, i=n; i--; y+=1.0)
1885 229 bertin
      norm += (*(maskt++) = INTERPF(y));
1886
    norm = norm>0.0? 1.0/norm : 1.0;
1887
    maskt -= n;
1888
    for (i=n; i--;)
1889
      *(maskt++) *= norm;
1890 221 bertin
    }
1891
 
1892
/* Initialize destination buffer to zero */
1893 285 bertin
  memset(outpix, 0, (size_t)subprofit->nobjpix*sizeof(PIXTYPE));
1894 221 bertin
 
1895
/* Make the interpolation in y and transpose once again */
1896
  dpixin0 = pixinout;
1897 285 bertin
  pixout0 = outpix+ixsout+iysout*subprofit->objnaxisn[0];
1898 221 bertin
  for (k=nxout; k--; dpixin0+=nyin, pixout0++)
1899
    {
1900
    maskt = mask;
1901
    nmaskt = nmask;
1902
    startt = start;
1903
    pixout = pixout0;
1904 285 bertin
    for (j=nyout; j--; pixout+=subprofit->objnaxisn[0])
1905 221 bertin
      {
1906
      dpixin = dpixin0+*(startt++);
1907
      val = 0.0;
1908
      for (i=*(nmaskt++); i--;)
1909
        val += *(maskt++)**(dpixin++);
1910
       *pixout = (PIXTYPE)(factor*val);
1911
      }
1912
    }
1913
 
1914
/* Free memory */
1915
  free(pixinout);
1916
  free(mask);
1917
  free(nmask);
1918
  free(start);
1919
 
1920
  return RETURN_OK;
1921 47 bertin
  }
1922
 
1923
 
1924 285 bertin
/****** subprofit_convolve ****************************************************
1925
PROTO   void subprofit_convolve(subprofitstruct *subprofit, float *modpix)
1926 142 bertin
PURPOSE Convolve a model image with the local PSF.
1927 285 bertin
INPUT   Pointer to the subprofit structure,
1928 142 bertin
        Pointer to the image raster.
1929 46 bertin
OUTPUT  -.
1930 117 bertin
NOTES   -.
1931 46 bertin
AUTHOR  E. Bertin (IAP)
1932 285 bertin
VERSION 05/05/2012
1933 46 bertin
 ***/
1934 285 bertin
void    subprofit_convolve(subprofitstruct *subprofit, float *modpix)
1935 46 bertin
  {
1936 285 bertin
  if (!subprofit->psfdft)
1937
    subprofit_makedft(subprofit);
1938 46 bertin
 
1939 285 bertin
  fft_conv(modpix, subprofit->psfdft, subprofit->modnaxisn);
1940 52 bertin
 
1941
  return;
1942
  }
1943
 
1944
 
1945 285 bertin
/****** subprofit_makedft ****************************************************
1946
PROTO   void subprofit_makedft(subprofitstruct *subprofit)
1947 52 bertin
PURPOSE Create the Fourier transform of the descrambled PSF component.
1948 285 bertin
INPUT   Pointer to the subprofit structure.
1949 52 bertin
OUTPUT  -.
1950 117 bertin
NOTES   -.
1951 52 bertin
AUTHOR  E. Bertin (IAP)
1952 285 bertin
VERSION 06/05/2012
1953 52 bertin
 ***/
1954 285 bertin
void    subprofit_makedft(subprofitstruct *subprofit)
1955 52 bertin
  {
1956
   psfstruct    *psf;
1957 201 bertin
   float      *mask,*maskt, *ppix;
1958
   float       dx,dy, r,r2,rmin,rmin2,rmax,rmax2,rsig,invrsig2;
1959 52 bertin
   int          width,height,npix,offset, psfwidth,psfheight,psfnpix,
1960
                cpwidth, cpheight,hcpwidth,hcpheight, i,j,x,y;
1961
 
1962 285 bertin
  if (!(psf=subprofit->psf))
1963 52 bertin
    return;
1964
 
1965 285 bertin
  psfwidth = subprofit->modnaxisn[0];
1966
  psfheight = subprofit->modnaxisn[1];
1967 52 bertin
  psfnpix = psfwidth*psfheight;
1968 285 bertin
  width = subprofit->modnaxisn[0];
1969
  height = subprofit->modnaxisn[1];
1970 52 bertin
  npix = width*height;
1971 201 bertin
  QCALLOC(mask, float, npix);
1972 52 bertin
  cpwidth = (width>psfwidth)?psfwidth:width;
1973
  hcpwidth = cpwidth>>1;
1974
  cpwidth = hcpwidth<<1;
1975
  offset = width - cpwidth;
1976
  cpheight = (height>psfheight)?psfheight:height;
1977
  hcpheight = cpheight>>1;
1978
  cpheight = hcpheight<<1;
1979
 
1980
/* Frame and descramble the PSF data */
1981 285 bertin
  ppix = subprofit->psfpix + (psfheight/2)*psfwidth + psfwidth/2;
1982 52 bertin
  maskt = mask;
1983
  for (j=hcpheight; j--; ppix+=psfwidth)
1984 46 bertin
    {
1985 52 bertin
    for (i=hcpwidth; i--;)
1986
      *(maskt++) = *(ppix++);
1987
    ppix -= cpwidth;
1988
    maskt += offset;
1989
    for (i=hcpwidth; i--;)
1990
      *(maskt++) = *(ppix++);
1991 46 bertin
    }
1992
 
1993 285 bertin
  ppix = subprofit->psfpix + ((psfheight/2)-hcpheight)*psfwidth + psfwidth/2;
1994 52 bertin
  maskt += width*(height-cpheight);
1995
  for (j=hcpheight; j--; ppix+=psfwidth)
1996
    {
1997
    for (i=hcpwidth; i--;)
1998
      *(maskt++) = *(ppix++);
1999
    ppix -= cpwidth;
2000
    maskt += offset;
2001
    for (i=hcpwidth; i--;)
2002
      *(maskt++) = *(ppix++);
2003
    }
2004 46 bertin
 
2005 52 bertin
/* Truncate to a disk that has diameter = (box width) */
2006
  rmax = cpwidth - 1.0 - hcpwidth;
2007
  if (rmax > (r=hcpwidth))
2008
    rmax = r;
2009
  if (rmax > (r=cpheight-1.0-hcpheight))
2010
    rmax = r;
2011
  if (rmax > (r=hcpheight))
2012
    rmax = r;
2013
  if (rmax<1.0)
2014
    rmax = 1.0;
2015
  rmax2 = rmax*rmax;
2016 285 bertin
  rsig = psf->fwhm/subprofit->pixstep;
2017 52 bertin
  invrsig2 = 1/(2*rsig*rsig);
2018
  rmin = rmax - (3*rsig);     /* 3 sigma annulus (almost no aliasing) */
2019
  rmin2 = rmin*rmin;
2020
 
2021
  maskt = mask;
2022
  dy = 0.0;
2023
  for (y=hcpheight; y--; dy+=1.0)
2024
    {
2025
    dx = 0.0;
2026
    for (x=hcpwidth; x--; dx+=1.0, maskt++)
2027
      if ((r2=dx*dx+dy*dy)>rmin2)
2028 201 bertin
        *maskt *= (r2>rmax2)?0.0:expf((2*rmin*sqrtf(r2)-r2-rmin2)*invrsig2);
2029 52 bertin
    dx = -hcpwidth;
2030
    maskt += offset;
2031
    for (x=hcpwidth; x--; dx+=1.0, maskt++)
2032
      if ((r2=dx*dx+dy*dy)>rmin2)
2033 201 bertin
        *maskt *= (r2>rmax2)?0.0:expf((2*rmin*sqrtf(r2)-r2-rmin2)*invrsig2);
2034 52 bertin
    }
2035
  dy = -hcpheight;
2036
  maskt += width*(height-cpheight);
2037
  for (y=hcpheight; y--; dy+=1.0)
2038
    {
2039
    dx = 0.0;
2040
    for (x=hcpwidth; x--; dx+=1.0, maskt++)
2041
      if ((r2=dx*dx+dy*dy)>rmin2)
2042 201 bertin
        *maskt *= (r2>rmax2)?0.0:expf((2*rmin*sqrtf(r2)-r2-rmin2)*invrsig2);
2043 52 bertin
    dx = -hcpwidth;
2044
    maskt += offset;
2045
    for (x=hcpwidth; x--; dx+=1.0, maskt++)
2046
      if ((r2=dx*dx+dy*dy)>rmin2)
2047 201 bertin
        *maskt *= (r2>rmax2)?0.0:expf((2*rmin*sqrtf(r2)-r2-rmin2)*invrsig2);
2048 52 bertin
    }
2049
 
2050
/* Finally move to Fourier space */
2051 285 bertin
  subprofit->psfdft = fft_rtf(mask, subprofit->modnaxisn);
2052 52 bertin
 
2053
  free(mask);
2054
 
2055 46 bertin
  return;
2056
  }
2057
 
2058
 
2059 285 bertin
/****** subprofit_copyobjpix **************************************************
2060
PROTO   int subprofit_copyobjpix(subprofitstruct *subprofit,
2061
                                subimagestruct *subimage)
2062
PURPOSE Copy a piece of the input object image to a subprofit structure.
2063
INPUT   Pointer to the subprofit structure,
2064
        subimagestruct *subimage.
2065 49 bertin
OUTPUT  The number of valid pixels copied.
2066 285 bertin
NOTES   -.
2067 49 bertin
AUTHOR  E. Bertin (IAP)
2068 285 bertin
VERSION 06/05/2012
2069 49 bertin
 ***/
2070 285 bertin
int     subprofit_copyobjpix(subprofitstruct *subprofit,
2071
                                subimagestruct *subimage)
2072 49 bertin
  {
2073 285 bertin
   fieldstruct          *field, *wfield;
2074 283 bertin
   float                dx, dy2, dr2, rad2;
2075
   PIXTYPE              *pixin,*spixin, *wpixin,*swpixin, *pixout,*wpixout,
2076
                        backnoise2, invgain, satlevel, wthresh, pix,spix,
2077
                        wpix,swpix;
2078 285 bertin
   int                  i,x,y, xmin,xmax,ymin,ymax, w,h,dw, npix, off,
2079
                        gainflag,badflag, sflag, sx,sy,sn, ix,iy, win,hin;
2080 49 bertin
 
2081 285 bertin
  field = subimage->field;
2082
  wfield = subimage->wfield;
2083
 
2084 49 bertin
/* First put the image background to -BIG */
2085 285 bertin
  pixout = subprofit->objpix;
2086
  wpixout = subprofit->objweight;
2087
  for (i=subprofit->objnaxisn[0]*subprofit->objnaxisn[1]; i--;)
2088 117 bertin
    {
2089 49 bertin
    *(pixout++) = -BIG;
2090 117 bertin
    *(wpixout++) = 0.0;
2091
    }
2092 285 bertin
 
2093
  ix = subprofit->ix - subimage->immin[0];
2094
  iy = subprofit->iy - subimage->immin[1];
2095 283 bertin
  win = subimage->imsize[0];
2096
  hin = subimage->imsize[1];
2097 285 bertin
 
2098 117 bertin
  backnoise2 = field->backsig*field->backsig;
2099 285 bertin
  sn = (int)subprofit->subsamp;
2100 211 bertin
  sflag = (sn>1);
2101 285 bertin
  w = subprofit->objnaxisn[0]*sn;
2102
  h = subprofit->objnaxisn[1]*sn;
2103 211 bertin
  if (sflag)
2104
    backnoise2 *= (PIXTYPE)sn;
2105 117 bertin
  invgain = (field->gain > 0.0) ? 1.0/field->gain : 0.0;
2106 285 bertin
  satlevel = field->satur_level - subimage->bkg;
2107 117 bertin
  rad2 = h/2.0;
2108
  if (rad2 > w/2.0)
2109
    rad2 = w/2.0;
2110
  rad2 *= rad2;
2111 285 bertin
 
2112 49 bertin
/* Set the image boundaries */
2113 285 bertin
  pixout = subprofit->objpix;
2114
  wpixout = subprofit->objweight;
2115 49 bertin
  ymin = iy-h/2;
2116
  ymax = ymin + h;
2117 248 bertin
  if (ymin<0)
2118 49 bertin
    {
2119 248 bertin
    off = (-ymin-1)/sn + 1;
2120 285 bertin
    pixout += off*subprofit->objnaxisn[0];
2121
    wpixout += off*subprofit->objnaxisn[0];
2122 211 bertin
    ymin += off*sn;
2123 49 bertin
    }
2124 248 bertin
  if (ymax>hin)
2125
    ymax -= ((ymax-hin-1)/sn + 1)*sn;
2126 285 bertin
 
2127 49 bertin
  xmin = ix-w/2;
2128
  xmax = xmin + w;
2129 211 bertin
  dw = 0;
2130 248 bertin
  if (xmax>win)
2131 49 bertin
    {
2132 248 bertin
    off = (xmax-win-1)/sn + 1;
2133 211 bertin
    dw += off;
2134
    xmax -= off*sn;
2135 49 bertin
    }
2136
  if (xmin<0)
2137
    {
2138 211 bertin
    off = (-xmin-1)/sn + 1;
2139
    pixout += off;
2140
    wpixout += off;
2141
    dw += off;
2142
    xmin += off*sn;
2143 49 bertin
    }
2144 285 bertin
 
2145 49 bertin
/* Copy the right pixels to the destination */
2146
  npix = 0;
2147 117 bertin
  if (wfield)
2148 49 bertin
    {
2149 117 bertin
    wthresh = wfield->weight_thresh;
2150 285 bertin
    gainflag = field->weightgain_flag;
2151 211 bertin
    if (sflag)
2152 117 bertin
      {
2153 211 bertin
/*---- Sub-sampling case */
2154
      for (y=ymin; y<ymax; y+=sn, pixout+=dw,wpixout+=dw)
2155 117 bertin
        {
2156 211 bertin
        for (x=xmin; x<xmax; x+=sn)
2157 117 bertin
          {
2158 211 bertin
          pix = wpix = 0.0;
2159
          badflag = 0;
2160
          for (sy=0; sy<sn; sy++)
2161
            {
2162
            dy2 = (y+sy-iy);
2163
            dy2 *= dy2;
2164
            dx = (x-ix);
2165 248 bertin
            off = x + (y+sy)*win;
2166 283 bertin
            spixin = subimage->image + off;
2167
            swpixin = subimage->weight + off;
2168 211 bertin
            for (sx=sn; sx--;)
2169
              {
2170
              dr2 = dy2 + dx*dx;
2171
              dx++;
2172
              spix = *(spixin++);
2173
              swpix = *(swpixin++);
2174
              if (dr2<rad2 && spix>-BIG && spix<satlevel && swpix<wthresh)
2175
                {
2176
                pix += spix;
2177
                wpix += swpix;
2178
                }
2179
              else
2180
                badflag=1;
2181
              }
2182
            }
2183
          *(pixout++) = pix;
2184 248 bertin
          if (!badflag) /* A single bad pixel ruins is all (saturation, etc.)*/
2185 211 bertin
            {
2186
            *(wpixout++) = 1.0 / sqrt(wpix+(pix>0.0?
2187 117 bertin
                (gainflag? pix*wpix/backnoise2:pix)*invgain : 0.0));
2188 211 bertin
            npix++;
2189
            }
2190
          else
2191
            *(wpixout++) = 0.0;
2192 117 bertin
          }
2193
        }
2194
      }
2195 211 bertin
    else
2196
      for (y=ymin; y<ymax; y++, pixout+=dw,wpixout+=dw)
2197
        {
2198
        dy2 = y-iy;
2199
        dy2 *= dy2;
2200 283 bertin
        pixin = subimage->image + xmin + y*win;
2201
        wpixin = subimage->weight + xmin + y*win;
2202 211 bertin
        for (x=xmin; x<xmax; x++)
2203
          {
2204
          dx = x-ix;
2205
          dr2 = dy2 + dx*dx;
2206
          pix = *(pixin++);
2207
          wpix = *(wpixin++);
2208
          if (dr2<rad2 && pix>-BIG && pix<satlevel && wpix<wthresh)
2209
            {
2210
            *(pixout++) = pix;
2211
            *(wpixout++) = 1.0 / sqrt(wpix+(pix>0.0?
2212
                (gainflag? pix*wpix/backnoise2:pix)*invgain : 0.0));
2213
            npix++;
2214
            }
2215
          else
2216
            *(pixout++) = *(wpixout++) = 0.0;
2217
          }
2218
        }
2219 49 bertin
    }
2220 117 bertin
  else
2221 211 bertin
    {
2222
    if (sflag)
2223 117 bertin
      {
2224 211 bertin
/*---- Sub-sampling case */
2225
      for (y=ymin; y<ymax; y+=sn, pixout+=dw, wpixout+=dw)
2226 117 bertin
        {
2227 211 bertin
        for (x=xmin; x<xmax; x+=sn)
2228 117 bertin
          {
2229 211 bertin
          pix = 0.0;
2230
          badflag = 0;
2231
          for (sy=0; sy<sn; sy++)
2232
            {
2233
            dy2 = y+sy-iy;
2234
            dy2 *= dy2;
2235
            dx = x-ix;
2236 283 bertin
            spixin = subimage->image + x + (y+sy)*win;
2237 211 bertin
            for (sx=sn; sx--;)
2238
              {
2239
              dr2 = dy2 + dx*dx;
2240
              dx++;
2241
              spix = *(spixin++);
2242
              if (dr2<rad2 && spix>-BIG && spix<satlevel)
2243
                pix += spix;
2244
              else
2245
                badflag=1;
2246
              }
2247
            }
2248
          *(pixout++) = pix;
2249 248 bertin
          if (!badflag) /* A single bad pixel ruins is all (saturation, etc.)*/
2250 211 bertin
            {
2251
            *(wpixout++) = 1.0 / sqrt(backnoise2 + (pix>0.0?pix*invgain:0.0));
2252
            npix++;
2253
            }
2254
          else
2255
            *(wpixout++) = 0.0;
2256 117 bertin
          }
2257
        }
2258
      }
2259 211 bertin
    else
2260
      for (y=ymin; y<ymax; y++, pixout+=dw,wpixout+=dw)
2261
        {
2262
        dy2 = y-iy;
2263
        dy2 *= dy2;
2264 285 bertin
        pixin = subimage->image + xmin + y*win;
2265 211 bertin
        for (x=xmin; x<xmax; x++)
2266
          {
2267
          dx = x-ix;
2268
          dr2 = dy2 + dx*dx;
2269
          pix = *(pixin++);
2270
          if (dr2<rad2 && pix>-BIG && pix<satlevel)
2271
            {
2272
            *(pixout++) = pix;
2273
            *(wpixout++) = 1.0 / sqrt(backnoise2 + (pix>0.0?pix*invgain : 0.0));
2274
            npix++;
2275
            }
2276
          else
2277
            *(pixout++) = *(wpixout++) = 0.0;
2278
          }
2279
        }
2280
    }
2281 117 bertin
 
2282 49 bertin
  return npix;
2283
  }
2284 285 bertin
 
2285 49 bertin
 
2286 285 bertin
/****** subprofit_submodpix *****************************************************
2287
PROTO   void    subprofit_submodpix(subprofitstruct *subprofitobj,
2288
                        subprofitstruct *subprofitmod, float fac)
2289 259 bertin
PURPOSE Subtract a rasterized model from the objpix of another model structure.
2290 285 bertin
INPUT   Pointer to the sub-profile structure containing the object pixels,
2291
        Pointer to the sub-profile structure containing the model raster,
2292 259 bertin
        multiplicative factor to apply to the model pixels.
2293
OUTPUT  -.
2294
NOTES   -.
2295
AUTHOR  E. Bertin (IAP)
2296 285 bertin
VERSION 05/05/2012
2297 259 bertin
 ***/
2298 285 bertin
void    subprofit_submodpix(subprofitstruct *subprofitobj,
2299
                        subprofitstruct *subprofitmod, float fac)
2300 259 bertin
  {
2301
   float        dx, dy2, dr2, rad2;
2302
   PIXTYPE      *pixin,*spixin, *pixout,
2303
                pix,spix;
2304
   int          i,x,y, xmin,xmax,ymin,ymax, w,h,dw, npix, off, gainflag,
2305
                badflag, sflag, sx,sy,sn, ix,iy, win,hin;
2306
 
2307
/* Don't go further if out of frame!! */
2308 285 bertin
  win = subprofitmod->objnaxisn[0];
2309
  hin = subprofitmod->objnaxisn[1];
2310
  ix = subprofitobj->ix - (subprofitmod->ix - win/2);
2311
  iy = subprofitobj->iy - (subprofitmod->iy - hin/2);
2312 259 bertin
 
2313 285 bertin
  sn = (int)subprofitobj->subsamp;
2314 259 bertin
  sflag = (sn>1);
2315 285 bertin
  w = subprofitobj->objnaxisn[0]*sn;
2316
  h = subprofitobj->objnaxisn[1]*sn;
2317 259 bertin
  rad2 = h/2.0;
2318
  if (rad2 > w/2.0)
2319
    rad2 = w/2.0;
2320
  rad2 *= rad2;
2321
 
2322
/* Set the image boundaries */
2323 285 bertin
  pixout = subprofitobj->objpix;
2324 259 bertin
  ymin = iy-h/2;
2325
  ymax = ymin + h;
2326
  if (ymin<0)
2327
    {
2328
    off = (-ymin-1)/sn + 1;
2329 285 bertin
    pixout += off*subprofitobj->objnaxisn[0];
2330 259 bertin
    ymin += off*sn;
2331
    }
2332
  if (ymax>hin)
2333
    ymax -= ((ymax-hin-1)/sn + 1)*sn;
2334
 
2335
  xmin = ix-w/2;
2336
  xmax = xmin + w;
2337
  dw = 0;
2338
  if (xmax>win)
2339
    {
2340
    off = (xmax-win-1)/sn + 1;
2341
    dw += off;
2342
    xmax -= off*sn;
2343
    }
2344
  if (xmin<0)
2345
    {
2346
    off = (-xmin-1)/sn + 1;
2347
    pixout += off;
2348
    dw += off;
2349
    xmin += off*sn;
2350
    }
2351
 
2352
/* Copy the right pixels to the destination */
2353
  npix = 0;
2354
  if (sflag)
2355
    {
2356
/*-- Sub-sampling case */
2357
    for (y=ymin; y<ymax; y+=sn, pixout+=dw)
2358
      {
2359
      for (x=xmin; x<xmax; x+=sn)
2360
        {
2361
        pix = 0.0;
2362
        badflag = 0;
2363
        for (sy=0; sy<sn; sy++)
2364
          {
2365
          dy2 = y+sy-iy;
2366
          dy2 *= dy2;
2367
          dx = x-ix;
2368 285 bertin
          spixin = subprofitmod->lmodpix + x + (y+sy)*win;
2369 259 bertin
          for (sx=sn; sx--;)
2370
            {
2371
            dr2 = dy2 + dx*dx;
2372
            dx++;
2373
            spix = *(spixin++);
2374
            pix += spix;
2375
            }
2376
          }
2377
        *(pixout++) -= fac*pix;
2378
        }
2379
      }
2380 267 bertin
    }
2381
  else
2382
    for (y=ymin; y<ymax; y++, pixout+=dw)
2383
      {
2384
      dy2 = y-iy;
2385
      dy2 *= dy2;
2386 285 bertin
      pixin = subprofitmod->lmodpix + xmin + y*win;
2387 267 bertin
      for (x=xmin; x<xmax; x++)
2388 259 bertin
        {
2389 267 bertin
        dx = x-ix;
2390
        dr2 = dy2 + dx*dx;
2391
        *(pixout++) -= fac**(pixin++);
2392 259 bertin
        }
2393 267 bertin
      }
2394
 
2395 259 bertin
  return;
2396
  }
2397
 
2398
 
2399 285 bertin
/****** subprofit_spiralindex *************************************************
2400
PROTO   float subprofit_spiralindex(subprofitstruct *subprofit)
2401 77 bertin
PURPOSE Compute the spiral index of a galaxy image (positive for arms
2402
        extending counter-clockwise and negative for arms extending CW, 0 for
2403
        no spiral pattern).
2404 146 bertin
INPUT   Profile-fitting structure.
2405 77 bertin
OUTPUT  Vector of residuals.
2406 117 bertin
NOTES   -.
2407 77 bertin
AUTHOR  E. Bertin (IAP)
2408 285 bertin
VERSION 06/05/2012
2409 77 bertin
 ***/
2410 285 bertin
float subprofit_spiralindex(subprofitstruct *subprofit, obj2struct *obj2)
2411 77 bertin
  {
2412 201 bertin
   float        *dx,*dy, *fdx,*fdy, *gdx,*gdy, *gdxt,*gdyt, *pix,
2413 77 bertin
                fwhm, invtwosigma2, hw,hh, ohw,ohh, x,y,xstart, tx,ty,txstart,
2414 78 bertin
                gx,gy, r2, spirindex, invsig, val, sep;
2415
   PIXTYPE      *fpix;
2416 77 bertin
   int          i,j, npix;
2417 78 bertin
 
2418 285 bertin
  npix = subprofit->objnaxisn[0]*subprofit->objnaxisn[1];
2419 77 bertin
 
2420
/* Compute simple derivative vectors at a fraction of the object scale */
2421 285 bertin
  fwhm = subprofit->guessradius * 2.0 / 4.0;
2422 77 bertin
  if (fwhm < 2.0)
2423
    fwhm = 2.0;
2424 78 bertin
  sep = 2.0;
2425 77 bertin
 
2426
  invtwosigma2 = -(2.35*2.35/(2.0*fwhm*fwhm));
2427 285 bertin
  hw = (float)(subprofit->objnaxisn[0]/2);
2428
  ohw = subprofit->objnaxisn[0] - hw;
2429
  hh = (float)(subprofit->objnaxisn[1]/2);
2430
  ohh = subprofit->objnaxisn[1] - hh;
2431 77 bertin
  txstart = -hw;
2432
  ty = -hh;
2433 201 bertin
  QMALLOC(dx, float, npix);
2434 77 bertin
  pix = dx;
2435 285 bertin
  for (j=subprofit->objnaxisn[1]; j--; ty+=1.0)
2436 77 bertin
    {
2437
    tx = txstart;
2438
    y = ty < -0.5? ty + hh : ty - ohh;
2439 285 bertin
    for (i=subprofit->objnaxisn[0]; i--; tx+=1.0)
2440 77 bertin
      {
2441
      x = tx < -0.5? tx + hw : tx - ohw;
2442 78 bertin
      *(pix++) = exp(invtwosigma2*((x+sep)*(x+sep)+y*y))
2443
                - exp(invtwosigma2*((x-sep)*(x-sep)+y*y));
2444 77 bertin
      }
2445
    }
2446 201 bertin
  QMALLOC(dy, float, npix);
2447 77 bertin
  pix = dy;
2448
  ty = -hh;
2449 285 bertin
  for (j=subprofit->objnaxisn[1]; j--; ty+=1.0)
2450 77 bertin
    {
2451
    tx = txstart;
2452
    y = ty < -0.5? ty + hh : ty - ohh;
2453 285 bertin
    for (i=subprofit->objnaxisn[0]; i--; tx+=1.0)
2454 77 bertin
      {
2455
      x = tx < -0.5? tx + hw : tx - ohw;
2456 78 bertin
      *(pix++) = exp(invtwosigma2*(x*x+(y+sep)*(y+sep)))
2457
                - exp(invtwosigma2*(x*x+(y-sep)*(y-sep)));
2458 77 bertin
      }
2459
    }
2460
 
2461 201 bertin
  QMALLOC(gdx, float, npix);
2462 77 bertin
  gdxt = gdx;
2463 285 bertin
  fpix = subprofit->objpix;
2464
  invsig = npix/subprofit->sigma;
2465 78 bertin
  for (i=npix; i--; fpix++)
2466 77 bertin
    {
2467 78 bertin
    val = *fpix > -1e29? *fpix*invsig : 0.0;
2468
    *(gdxt++) = (val>0.0? log(1.0+val) : -log(1.0-val));
2469 77 bertin
    }
2470 123 bertin
  gdy = NULL;                   /* to avoid gcc -Wall warnings */
2471 201 bertin
  QMEMCPY(gdx, gdy, float, npix);
2472 285 bertin
  fdx = fft_rtf(dx, subprofit->objnaxisn);
2473
  fft_conv(gdx, fdx, subprofit->objnaxisn);
2474
  fdy = fft_rtf(dy, subprofit->objnaxisn);
2475
  fft_conv(gdy, fdy, subprofit->objnaxisn);
2476 77 bertin
 
2477
/* Compute estimator */
2478 285 bertin
  invtwosigma2 = -1.18*1.18/(2.0*subprofit->guessradius*subprofit->guessradius);
2479 258 bertin
  xstart = -hw - obj2->mx + (int)(obj2->mx+0.49999);
2480
  y = -hh -  obj2->my + (int)(obj2->my+0.49999);;
2481 77 bertin
  spirindex = 0.0;
2482
  gdxt = gdx;
2483
  gdyt = gdy;
2484 285 bertin
  for (j=subprofit->objnaxisn[1]; j--; y+=1.0)
2485 77 bertin
    {
2486
    x = xstart;
2487 285 bertin
    for (i=subprofit->objnaxisn[0]; i--; x+=1.0)
2488 77 bertin
      {
2489
      gx = *(gdxt++);
2490
      gy = *(gdyt++);
2491
      if ((r2=x*x+y*y)>0.0)
2492 78 bertin
        spirindex += (x*y*(gx*gx-gy*gy)+gx*gy*(y*y-x*x))/r2
2493 77 bertin
                        * exp(invtwosigma2*r2);
2494
      }
2495
    }
2496
 
2497
  free(dx);
2498
  free(dy);
2499
  free(fdx);
2500
  free(fdy);
2501
  free(gdx);
2502
  free(gdy);
2503
 
2504
  return spirindex;
2505
  }
2506
 
2507
 
2508 285 bertin
/****** profit_moments *******************************************************
2509 207 bertin
PROTO   void profit_moments(profitstruct *profit, obj2struct *obj2)
2510 89 bertin
PURPOSE Compute the 2nd order moments from the unconvolved object model.
2511 207 bertin
INPUT   Profile-fitting structure,
2512
        Pointer to obj2 structure.
2513 89 bertin
OUTPUT  -.
2514 117 bertin
NOTES   -.
2515 89 bertin
AUTHOR  E. Bertin (IAP)
2516 245 bertin
VERSION 22/04/2011
2517 89 bertin
 ***/
2518 207 bertin
void     profit_moments(profitstruct *profit, obj2struct *obj2)
2519 89 bertin
  {
2520 221 bertin
   profstruct   *prof;
2521 226 bertin
   double       dpdmx2[6], cov[4],
2522
                *jac,*jact, *pjac,*pjact, *dcovar,*dcovart,
2523
                *dmx2,*dmy2,*dmxy,
2524
                m0,invm0, mx2,my2,mxy, den,invden,
2525
                temp, temp2,invtemp2,invstemp2,
2526
                pmx2,theta, flux, dval;
2527 225 bertin
   float         *covart;
2528 245 bertin
   int          findex[MODEL_NMAX],
2529 225 bertin
                i,j,p, nparam;
2530 89 bertin
 
2531 221 bertin
/*  hw = (float)(profit->modnaxisn[0]/2);*/
2532
/*  hh = (float)(profit->modnaxisn[1]/2);*/
2533
/*  r2max = hw<hh? hw*hw : hh*hh;*/
2534
/*  xstart = -hw;*/
2535
/*  y = -hh;*/
2536
/*  pix = profit->modpix;*/
2537
/*  mx2 = my2 = mxy = mx = my = sum = 0.0;*/
2538
/*  for (iy=profit->modnaxisn[1]; iy--; y+=1.0)*/
2539
/*    {*/
2540
/*    x = xstart;*/
2541
/*    for (ix=profit->modnaxisn[0]; ix--; x+=1.0)*/
2542
/*      if (y*y+x*x <= r2max)*/
2543
/*        {*/
2544
/*        val = *(pix++);*/
2545
/*        sum += val;*/
2546
/*        mx  += val*x;*/
2547
/*        my  += val*y;*/
2548
/*        mx2 += val*x*x;*/
2549
/*        mxy += val*x*y;*/
2550
/*        my2 += val*y*y;*/
2551
/*        }*/
2552
/*      else*/
2553
/*        pix++;*/
2554
/*    }*/
2555
 
2556
/*  if (sum <= 1.0/BIG)*/
2557
/*    sum = 1.0;*/
2558
/*  mx /= sum;*/
2559
/*  my /= sum;*/
2560
/*  obj2->prof_mx2 = mx2 = mx2/sum - mx*mx;*/
2561
/*  obj2->prof_my2 = my2 = my2/sum - my*my;*/
2562
/*  obj2->prof_mxy = mxy = mxy/sum - mx*my;*/
2563
 
2564 225 bertin
  nparam = profit->nparam;
2565 226 bertin
  if (FLAG(obj2.prof_e1err) || FLAG(obj2.prof_pol1err))
2566 225 bertin
    {
2567
/*-- Set up Jacobian matrices */
2568
    QCALLOC(jac, double, nparam*3);
2569 226 bertin
    QMALLOC(pjac, double, (nparam<2? 6 : nparam*3));
2570
    QMALLOC(dcovar, double, nparam*nparam);
2571
    dcovart = dcovar;
2572
    covart = profit->covar;
2573
    for (i=nparam*nparam; i--;)
2574
      *(dcovart++) = (double)(*(covart++));
2575 225 bertin
    dmx2 = jac;
2576
    dmy2 = jac+nparam;
2577
    dmxy = jac+2*nparam;
2578
    }
2579
  else
2580 226 bertin
    jac = pjac = dcovar = dmx2 = dmy2 = dmxy = NULL;
2581 225 bertin
 
2582 221 bertin
  m0 = mx2 = my2 = mxy = 0.0;
2583
  for (p=0; p<profit->nprof; p++)
2584 89 bertin
    {
2585 221 bertin
    prof = profit->prof[p];
2586 225 bertin
    findex[p] = prof_moments(profit, prof, pjac);
2587 285 bertin
    flux = *prof->flux[0];
2588 225 bertin
    m0 += flux;
2589
    mx2 += prof->mx2*flux;
2590
    my2 += prof->my2*flux;
2591
    mxy += prof->mxy*flux;
2592
    if (jac)
2593
      {
2594
      jact = jac;
2595
      pjact = pjac;
2596
      for (j=nparam*3; j--;)
2597
        *(jact++) += flux * *(pjact++);
2598
      }
2599 89 bertin
    }
2600 225 bertin
  invm0 = 1.0 / m0;
2601
  obj2->prof_mx2 = (mx2 *= invm0);
2602
  obj2->prof_my2 = (my2 *= invm0);
2603
  obj2->prof_mxy = (mxy *= invm0);
2604
/* Complete the flux derivative of moments */
2605
  if (jac)
2606
    {
2607
    for (p=0; p<profit->nprof; p++)
2608
      {
2609 226 bertin
      prof = profit->prof[p];
2610
      dmx2[findex[p]] = prof->mx2 - mx2;
2611
      dmy2[findex[p]] = prof->my2 - my2;
2612
      dmxy[findex[p]] = prof->mxy - mxy;
2613 225 bertin
      }
2614
    jact = jac;
2615
    for (j=nparam*3; j--;)
2616
      *(jact++) *= invm0;
2617
    }
2618 89 bertin
 
2619 206 bertin
/* Handle fully correlated profiles (which cause a singularity...) */
2620
  if ((temp2=mx2*my2-mxy*mxy)<0.00694)
2621 89 bertin
    {
2622 206 bertin
    mx2 += 0.0833333;
2623
    my2 += 0.0833333;
2624
    temp2 = mx2*my2-mxy*mxy;
2625
    }
2626
 
2627 226 bertin
/* Use the Jacobians to compute the moment covariance matrix */
2628
  if (jac)
2629
    propagate_covar(dcovar, jac, obj2->prof_mx2cov, nparam, 3,
2630
                                                pjac);  /* We re-use pjac */
2631
 
2632 225 bertin
  if (FLAG(obj2.prof_pol1))
2633 206 bertin
    {
2634 226 bertin
/*--- "Polarisation", i.e. module = (a^2-b^2)/(a^2+b^2) */
2635 206 bertin
    if (mx2+my2 > 1.0/BIG)
2636
      {
2637
      obj2->prof_pol1 = (mx2 - my2) / (mx2+my2);
2638
      obj2->prof_pol2 = 2.0*mxy / (mx2 + my2);
2639 226 bertin
      if (FLAG(obj2.prof_pol1err))
2640 225 bertin
        {
2641 226 bertin
/*------ Compute the Jacobian of polarisation */
2642
        invden = 1.0/(mx2+my2);
2643
        dpdmx2[0] =  2.0*my2*invden*invden;
2644
        dpdmx2