public software.psfex

[/] [branches/] [rhl/] [src/] [makeit2.c] - Blame information for rev 191

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 157 bertin
/*
2 191 rhl
*                               makeit2.c
3 2 bertin
*
4 191 rhl
* Main loop, with the I/O parts removed
5 2 bertin
*
6 157 bertin
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 2 bertin
*
8 157 bertin
*       This file part of:      PSFEx
9 2 bertin
*
10 177 bertin
*       Copyright:              (C) 1997-2012 Emmanuel Bertin -- IAP/CNRS/UPMC
11 2 bertin
*
12 157 bertin
*       License:                GNU General Public License
13
*
14
*       PSFEx 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
*       PSFEx 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 PSFEx.  If not, see <http://www.gnu.org/licenses/>.
24
*
25 183 bertin
*       Last modified:          19/07/2012
26 157 bertin
*
27
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
28 2 bertin
 
29
#ifdef HAVE_CONFIG_H
30
#include        "config.h"
31
#endif
32
 
33 183 bertin
#ifdef USE_THREADS
34
 #ifdef HAVE_MKL
35
  #include MKL_H
36
 #endif
37
#endif
38
 
39 8 bertin
#include        <math.h>
40 2 bertin
#include        <stdio.h>
41
#include        <stdlib.h>
42
#include        <string.h>
43
 
44
#include        "define.h"
45
#include        "types.h"
46
#include        "globals.h"
47
#include        "fits/fitscat.h"
48
#include        "check.h"
49 72 bertin
#include        "context.h"
50 77 bertin
#include        "cplot.h"
51 6 bertin
#include        "diagnostic.h"
52 76 bertin
#include        "field.h"
53 61 bertin
#include        "homo.h"
54 57 bertin
#include        "pca.h"
55 2 bertin
#include        "prefs.h"
56
#include        "psf.h"
57
#include        "sample.h"
58 4 bertin
#include        "xml.h"
59 2 bertin
 
60 57 bertin
psfstruct       *make_psf(setstruct *set, float psfstep,
61 77 bertin
                        float *basis, int nbasis, contextstruct *context);
62 4 bertin
 
63 191 rhl
/********************************** makeit_body ******************************/
64 2 bertin
/*
65
*/
66 191 rhl
void
67
makeit_body(
68
   fieldstruct          **fields,
69
   contextstruct        **context_,
70
   contextstruct        **fullcontext_,
71
   int                  free_sets       /* if false, don't free any sets as we don't own them */
72
   )
73 2 bertin
  {
74 177 bertin
   wcsstruct            *wcs;
75 191 rhl
   fieldstruct          *field;
76 73 bertin
   psfstruct            **cpsf,
77
                        *psf;
78 79 bertin
   setstruct            *set, *set2;
79 73 bertin
   contextstruct        *context, *fullcontext;
80 112 bertin
   char                 str[MAXCHAR];
81 191 rhl
   char                 **incatnames;
82 112 bertin
   float                **psfbasiss,
83
                        *psfsteps, *psfbasis, *basis,
84
                        psfstep, step;
85 109 bertin
   int                  c,i,p, ncat, ext, next, nmed, nbasis;
86 2 bertin
 
87 73 bertin
  incatnames = prefs.incat_name;
88
  ncat = prefs.ncat;
89 76 bertin
  next = fields[0]->next;
90 75 bertin
 
91 59 bertin
  psfstep = prefs.psf_step;
92
  psfsteps = NULL;
93
  nbasis = 0;
94 112 bertin
  psfbasis = NULL;
95
  psfbasiss = NULL;
96 59 bertin
 
97 72 bertin
/* Initialize context */
98 111 bertin
  NFPRINTF(OUTPUT, "Initializing contexts...");
99 72 bertin
  context = context_init(prefs.context_name, prefs.context_group,
100
                prefs.ncontext_group, prefs.group_deg, prefs.ngroup_deg,
101 110 bertin
                CONTEXT_REMOVEHIDDEN);
102 73 bertin
  fullcontext = context->npc?
103
                  context_init(prefs.context_name, prefs.context_group,
104
                    prefs.ncontext_group, prefs.group_deg, prefs.ngroup_deg,
105 110 bertin
                    CONTEXT_KEEPHIDDEN)
106 73 bertin
                : context;
107 191 rhl
  *context_ = context;
108
  *fullcontext_ = fullcontext;
109 72 bertin
 
110 114 bertin
  if (context->npc && ncat<2)
111
    warning("Hidden dependencies cannot be derived from",
112
        " a single catalog");
113
  else if (context->npc && prefs.stability_type == STABILITY_EXPOSURE)
114
    warning("Hidden dependencies have no effect",
115
        " in STABILITY_TYPE EXPOSURE mode");
116
 
117 111 bertin
/* Compute PSF steps */
118
  if (!prefs.psf_step)
119 56 bertin
    {
120 111 bertin
    NFPRINTF(OUTPUT, "Computing optimum PSF sampling steps...");
121 112 bertin
    if (prefs.newbasis_type==NEWBASIS_PCACOMMON
122
        || (prefs.stability_type == STABILITY_SEQUENCE
123
                && prefs.psf_mef_type == PSF_MEF_COMMON))
124
      {
125 119 bertin
      set = load_samples(incatnames, 0, ncat, ALL_EXTENSIONS, next, context);
126 112 bertin
      psfstep = (float)((set->fwhm/2.35)*0.5);
127 191 rhl
      if (free_sets) end_set(set);
128 112 bertin
      }
129 111 bertin
/*-- Need to derive a common pixel step for each ext */
130 112 bertin
    else if (prefs.newbasis_type == NEWBASIS_PCAINDEPENDENT
131
        || context->npc
132 111 bertin
        || (prefs.stability_type == STABILITY_SEQUENCE
133
                && prefs.psf_mef_type == PSF_MEF_INDEPENDENT))
134 57 bertin
      {
135 111 bertin
/*-- Run through all samples to derive a different pixel step for each extension */
136 72 bertin
      QMALLOC(psfsteps, float, next);
137
      for (ext=0 ; ext<next; ext++)
138
        {
139 119 bertin
        set = load_samples(incatnames, 0, ncat, ext, next, context);
140 111 bertin
        psfsteps[ext] = (float)(psfstep? psfstep : (set->fwhm/2.35)*0.5);
141 191 rhl
        if (free_sets) end_set(set);
142 72 bertin
        }
143
      }
144 111 bertin
    }
145
 
146
/* Derive a new common PCA basis for all extensions */
147
  if (prefs.newbasis_type==NEWBASIS_PCACOMMON)
148
    {
149 72 bertin
    QMALLOC(cpsf, psfstruct *, ncat*next);
150
    for (ext=0 ; ext<next; ext++)
151
      for (c=0; c<ncat; c++)
152
        {
153 112 bertin
        sprintf(str, "Computing new PCA image basis from %s...",
154
                fields[c]->rtcatname);
155
        NFPRINTF(OUTPUT, str);
156 119 bertin
        set = load_samples(incatnames, c, 1, ext, next, context);
157 112 bertin
        step = psfstep;
158 77 bertin
        cpsf[c+ext*ncat] = make_psf(set, psfstep, NULL, 0, context);
159 191 rhl
        if (free_sets) end_set(set);
160 72 bertin
        }
161
    nbasis = prefs.newbasis_number;
162 112 bertin
    psfbasis = pca_onsnaps(cpsf, ncat*next, nbasis);
163 72 bertin
    for (i=0 ; i<ncat*next; i++)
164
      psf_end(cpsf[i]);
165
    free(cpsf);
166
    }
167 111 bertin
/* Derive a new PCA basis for each extension */
168 112 bertin
  else if (prefs.newbasis_type == NEWBASIS_PCAINDEPENDENT)
169 111 bertin
    {
170
    nbasis = prefs.newbasis_number;
171 112 bertin
    QMALLOC(psfbasiss, float *, next);
172 111 bertin
    for (ext=0; ext<next; ext++)
173
      {
174 112 bertin
      if (psfsteps)
175
        step = psfsteps[ext];
176
      else
177
        step = psfstep;
178 111 bertin
      QMALLOC(cpsf, psfstruct *, ncat);
179
      for (c=0; c<ncat; c++)
180
        {
181 112 bertin
        sprintf(str, "Computing new PCA image basis from %s...",
182
                fields[c]->rtcatname);
183
        NFPRINTF(OUTPUT, str);
184 119 bertin
        set = load_samples(incatnames, c, 1, ext, next, context);
185 112 bertin
        cpsf[c] = make_psf(set, step, NULL, 0, context);
186 191 rhl
        if (free_sets) end_set(set);
187 111 bertin
        }
188 112 bertin
      psfbasiss[ext] = pca_onsnaps(cpsf, ncat, nbasis);
189 111 bertin
      for (c=0 ; c<ncat; c++)
190
        psf_end(cpsf[c]);
191
      free(cpsf);
192
      }
193
    }
194
 
195
  if (context->npc && prefs.hidden_mef_type == HIDDEN_MEF_COMMON)
196 109 bertin
/*-- Derive principal components of PSF variation from the whole mosaic */
197 57 bertin
    {
198 110 bertin
    p = 0;
199
    QMALLOC(cpsf, psfstruct *, ncat*next);
200
    for (c=0; c<ncat; c++)
201 2 bertin
      {
202 112 bertin
      sprintf(str, "Computing hidden dependency parameter(s) from %s...",
203
                fields[c]->rtcatname);
204
      NFPRINTF(OUTPUT, str);
205 110 bertin
      for (ext=0 ; ext<next; ext++)
206 57 bertin
        {
207 119 bertin
        set = load_samples(incatnames, c, 1, ext, next, context);
208 112 bertin
        if (psfsteps)
209
          step = psfsteps[ext];
210
        else
211
          step = psfstep;
212
        basis = psfbasiss? psfbasiss[ext] : psfbasis;
213
        cpsf[p++] = make_psf(set, step, basis, nbasis, context);
214 191 rhl
        if (free_sets) end_set(set);
215 57 bertin
        }
216
      }
217 110 bertin
    free(fullcontext->pc);
218
    fullcontext->pc = pca_oncomps(cpsf, next, ncat, context->npc);
219
    for (c=0 ; c<ncat*next; c++)
220
      psf_end(cpsf[c]);
221
    free(cpsf);
222
    }
223 57 bertin
 
224 119 bertin
/* Compute "final" PSF models */
225 110 bertin
  if (prefs.psf_mef_type == PSF_MEF_COMMON)
226
    {
227 107 bertin
    if (prefs.stability_type == STABILITY_SEQUENCE)
228
      {
229
/*---- Load all the samples at once */
230 119 bertin
      set = load_samples(incatnames, 0, ncat, ALL_EXTENSIONS, next, context);
231 112 bertin
      step = psfstep;
232
      basis = psfbasis;
233 119 bertin
      field_count(fields, set, COUNT_LOADED);
234 112 bertin
      psf = make_psf(set, step, basis, nbasis, context);
235 119 bertin
      field_count(fields, set, COUNT_ACCEPTED);
236 191 rhl
      if (free_sets) end_set(set);
237 112 bertin
      NFPRINTF(OUTPUT, "Computing final PSF model...");
238 110 bertin
      context_apply(context, psf, fields, ALL_EXTENSIONS, 0, ncat);
239 107 bertin
      psf_end(psf);
240
      }
241 110 bertin
    else
242
      for (c=0; c<ncat; c++)
243 107 bertin
        {
244
/*------ Load the samples for current exposure */
245 112 bertin
        sprintf(str, "Computing final PSF model from %s...",
246
                fields[c]->rtcatname);
247
        NFPRINTF(OUTPUT, str);
248 119 bertin
        set = load_samples(incatnames, c, 1, ALL_EXTENSIONS, next, context);
249 112 bertin
        if (psfstep)
250
          step = psfstep;
251
        else
252
          step = (float)((set->fwhm/2.35)*0.5);
253
        basis = psfbasis;
254 119 bertin
        field_count(fields, set, COUNT_LOADED);
255 112 bertin
        psf = make_psf(set, step, basis, nbasis, fullcontext);
256 119 bertin
        field_count(fields, set, COUNT_ACCEPTED);
257 191 rhl
        if (free_sets) end_set(set);
258 110 bertin
        context_apply(fullcontext, psf, fields, ALL_EXTENSIONS, c, 1);
259 107 bertin
        psf_end(psf);
260
        }
261 110 bertin
    }
262
  else
263
    for (ext=0 ; ext<next; ext++)
264
      {
265 112 bertin
      basis = psfbasiss? psfbasiss[ext] : psfbasis;
266
      if (context->npc && prefs.hidden_mef_type == HIDDEN_MEF_INDEPENDENT)
267 110 bertin
/*------ Derive principal components of PSF components */
268
        {
269
        QMALLOC(cpsf, psfstruct *, ncat);
270 112 bertin
        if (psfsteps)
271
          step = psfsteps[ext];
272
        else
273
          step = psfstep;
274 110 bertin
        for (c=0; c<ncat; c++)
275
          {
276 112 bertin
          if (next>1)
277
            sprintf(str,
278
                "Computing hidden dependency parameter(s) from %s[%d/%d]...",
279
                fields[c]->rtcatname, ext+1, next);
280
          else
281
            sprintf(str, "Computing hidden dependency parameter(s) from %s...",
282
                fields[c]->rtcatname);
283
          NFPRINTF(OUTPUT, str);
284 119 bertin
          set = load_samples(incatnames, c, 1, ext, next, context);
285
          field_count(fields, set, COUNT_LOADED);
286 112 bertin
          cpsf[c] = make_psf(set, step, basis, nbasis, context);
287 119 bertin
          field_count(fields, set, COUNT_ACCEPTED);
288 191 rhl
          if (free_sets) end_set(set);
289 110 bertin
          }
290
        free(fullcontext->pc);
291
        fullcontext->pc = pca_oncomps(cpsf, 1, ncat, context->npc);
292
        for (c=0 ; c<ncat; c++)
293
          psf_end(cpsf[c]);
294
        free(cpsf);
295
        }
296
 
297
      if (prefs.stability_type == STABILITY_SEQUENCE)
298
        {
299
/*------ Load all the samples at once */
300 112 bertin
        if (next>1)
301
          {
302
          sprintf(str, "Computing final PSF model for extension [%d/%d]...",
303
                ext+1, next);
304
          NFPRINTF(OUTPUT, str);
305
          }
306
        else
307
          NFPRINTF(OUTPUT, "Computing final PSF model...");
308 119 bertin
        set = load_samples(incatnames, 0, ncat, ext, next, fullcontext);
309 112 bertin
        if (psfstep)
310
          step = psfstep;
311
        else if (psfsteps)
312
          step = psfsteps[ext];
313
        else
314
          step = (float)((set->fwhm/2.35)*0.5);
315 119 bertin
        field_count(fields, set, COUNT_LOADED);
316 112 bertin
        psf = make_psf(set, step, basis, nbasis, fullcontext);
317 119 bertin
        field_count(fields, set, COUNT_ACCEPTED);
318 191 rhl
        if (free_sets) end_set(set);
319 110 bertin
        context_apply(fullcontext, psf, fields, ext, 0, ncat);
320
        psf_end(psf);
321
        }
322
      else
323
        for (c=0; c<ncat; c++)
324
          {
325
/*-------- Load the samples for current exposure */
326 112 bertin
          if (next>1)
327 177 bertin
            sprintf(str, "Reading data from %s[%d/%d]...",
328 112 bertin
                fields[c]->rtcatname, ext+1, next);
329
          else
330 177 bertin
            sprintf(str, "Reading data from %s...",
331 112 bertin
                fields[c]->rtcatname);
332
          NFPRINTF(OUTPUT, str);
333 119 bertin
          set = load_samples(incatnames, c, 1, ext, next, context);
334 112 bertin
          if (psfstep)
335
            step = psfstep;
336
          else if (psfsteps)
337
            step = psfsteps[ext];
338
          else
339
            step = (float)((set->fwhm/2.35)*0.5);
340 177 bertin
          if (next>1)
341
            sprintf(str, "Computing final PSF model for %s[%d/%d]...",
342
                fields[c]->rtcatname, ext+1, next);
343
          else
344
            sprintf(str, "Computing final PSF model for %s...",
345
                fields[c]->rtcatname);
346
          NFPRINTF(OUTPUT, str);
347 119 bertin
          field_count(fields, set, COUNT_LOADED);
348 112 bertin
          psf = make_psf(set, step, basis, nbasis, context);
349 119 bertin
          field_count(fields, set, COUNT_ACCEPTED);
350 191 rhl
          if (free_sets) end_set(set);
351 112 bertin
          context_apply(context, psf, fields, ext, c, 1);
352 110 bertin
          psf_end(psf);
353
          }
354
      }
355
 
356
  free(psfsteps);
357 112 bertin
  if (psfbasiss)
358 110 bertin
    {
359
    for (ext=0; ext<next; ext++)
360 112 bertin
      free(psfbasiss[ext]);
361
    free(psfbasiss);
362 110 bertin
    }
363 112 bertin
  else if (psfbasis)
364
    free(psfbasis);
365 110 bertin
 
366
/* Compute diagnostics and check-images */
367 183 bertin
#ifdef USE_THREADS
368
/* Force MKL using a single thread as diagnostic code is already multithreaded*/
369
 #ifdef HAVE_MKL
370
  if (prefs.context_nsnap>2)
371
    mkl_set_num_threads(1);
372
 #endif
373
#endif
374 110 bertin
  QIPRINTF(OUTPUT,
375 111 bertin
        "   filename      [ext] accepted/total samp. chi2/dof FWHM ellip."
376
        " resi. asym.");
377 110 bertin
  for (c=0; c<ncat; c++)
378
    {
379 176 bertin
    field = fields[c];
380 110 bertin
    for (ext=0 ; ext<next; ext++)
381
      {
382 176 bertin
      psf = field->psf[ext];
383 177 bertin
      wcs = field->wcs[ext];
384 114 bertin
      if (next>1)
385
        sprintf(str, "Computing diagnostics for %s[%d/%d]...",
386 176 bertin
                field->rtcatname, ext+1, next);
387 114 bertin
      else
388
        sprintf(str, "Computing diagnostics for %s...",
389 176 bertin
                field->rtcatname);
390 114 bertin
      NFPRINTF(OUTPUT, str);
391 79 bertin
/*---- Check PSF with individual datasets */
392 119 bertin
      set2 = load_samples(incatnames, c, 1, ext, next, context);
393 79 bertin
      psf->samples_loaded = set2->nsample;
394
      if (set2->nsample>1)
395
        {
396
/*------ Remove bad PSF candidates */
397 80 bertin
        psf_clean(psf, set2, prefs.prof_accuracy);
398 79 bertin
        psf->chi2 = set2->nsample? psf_chi2(psf, set2) : 0.0;
399
        }
400
      psf->samples_accepted = set2->nsample;
401 150 bertin
/*---- Compute diagnostics and field statistics */
402 77 bertin
      psf_diagnostic(psf);
403 177 bertin
      psf_wcsdiagnostic(psf, wcs);
404 110 bertin
      nmed = psf->nmed;
405 150 bertin
      field_stats(fields, set2);
406 110 bertin
/*---- Display stats for current catalog/extension */
407 111 bertin
      if (next>1)
408
        sprintf(str, "[%d/%d]", ext+1, next);
409
      else
410
        str[0] = '\0';
411
      QPRINTF(OUTPUT, "%-17.17s%-7.7s %5d/%-5d %6.2f %6.2f %6.2f  %4.2f"
412
        " %5.2f %5.2f\n",
413 176 bertin
        ext==0? field->rtcatname : "",
414 111 bertin
        str,
415 73 bertin
        psf->samples_accepted, psf->samples_loaded,
416
        psf->pixstep,
417
        psf->chi2,
418 142 bertin
        psf->moffat_fwhm,
419
        psf->moffat_ellipticity,
420 145 bertin
        psf->pfmoffat_residuals,
421 142 bertin
        psf->sym_residuals);
422 87 bertin
/*---- Save "Check-images" */
423
      for (i=0; i<prefs.ncheck_type; i++)
424
        if (prefs.check_type[i])
425
          {
426
          sprintf(str, "Saving CHECK-image #%d...", i+1);
427
          NFPRINTF(OUTPUT, str);
428 176 bertin
          check_write(field, set2, prefs.check_name[i], prefs.check_type[i],
429 87 bertin
                ext, next, prefs.check_cubeflag);
430
          }
431
/*---- Free memory */
432 191 rhl
      if (free_sets) end_set(set2);
433 77 bertin
      }
434 56 bertin
    }
435 2 bertin
 
436 183 bertin
#ifdef USE_THREADS
437
/* Back to multithreaded MKL */
438
 #ifdef HAVE_MKL
439
   mkl_set_num_threads(prefs.nthreads);
440
 #endif
441
#endif
442
 
443 2 bertin
  return;
444
  }
445
 
446 56 bertin
/****** make_psf *************************************************************
447 77 bertin
PROTO   psfstruct *make_psf(setstruct *set, float psfstep,
448
                        float *basis, int nbasis, contextstruct *context)
449 56 bertin
PURPOSE Make PSFs from a set of FITS binary catalogs.
450
INPUT   Pointer to a sample set,
451 57 bertin
        PSF sampling step,
452
        Pointer to basis image vectors,
453
        Number of basis vectors,
454 72 bertin
        Pointer to context structure.
455 56 bertin
OUTPUT  Pointer to the PSF structure.
456 57 bertin
NOTES   Diagnostics are computed only if diagflag != 0.
457 56 bertin
AUTHOR  E. Bertin (IAP)
458 136 bertin
VERSION 03/09/2009
459 56 bertin
 ***/
460 57 bertin
psfstruct       *make_psf(setstruct *set, float psfstep,
461 77 bertin
                        float *basis, int nbasis, contextstruct *context)
462 56 bertin
  {
463
   psfstruct            *psf;
464 65 bertin
   basistypenum         basistype;
465 136 bertin
   float                pixsize[2];
466 56 bertin
 
467 136 bertin
  pixsize[0] = (float)prefs.psf_pixsize[0];
468
  pixsize[1] = (float)prefs.psf_pixsize[1];
469 112 bertin
//  NFPRINTF(OUTPUT,"Initializing PSF modules...");
470 136 bertin
  psf = psf_init(context, prefs.psf_size, psfstep, pixsize, set->nsample);
471 56 bertin
 
472
  psf->samples_loaded = set->nsample;
473
  psf->fwhm = set->fwhm;
474 57 bertin
 
475 60 bertin
/* Make the basic PSF-model (1st pass) */
476 112 bertin
//  NFPRINTF(OUTPUT,"Modeling the PSF (1/3)...");
477 80 bertin
  psf_make(psf, set, 0.2);
478 57 bertin
  if (basis && nbasis)
479
    {
480
    QMEMCPY(basis, psf->basis, float, nbasis*psf->size[0]*psf->size[1]);
481
    psf->nbasis = nbasis;
482
    }
483 60 bertin
  else
484
    {
485 112 bertin
//    NFPRINTF(OUTPUT,"Generating the PSF model...");
486 65 bertin
    basistype = prefs.basis_type;
487
    if (basistype==BASIS_PIXEL_AUTO)
488
      basistype = (psf->fwhm < PSF_AUTO_FWHM)? BASIS_PIXEL : BASIS_NONE;
489 66 bertin
    psf_makebasis(psf, set, basistype, prefs.basis_number);
490 60 bertin
    }
491
  psf_refine(psf, set);
492 56 bertin
 
493
/* Remove bad PSF candidates */
494
  if (set->nsample>1)
495
    {
496 80 bertin
    psf_clean(psf, set, 0.2);
497 150 bertin
 
498 56 bertin
/*-- Make the basic PSF-model (2nd pass) */
499 112 bertin
//    NFPRINTF(OUTPUT,"Modeling the PSF (2/3)...");
500 80 bertin
    psf_make(psf, set, 0.1);
501 60 bertin
    psf_refine(psf, set);
502 56 bertin
    }
503 150 bertin
 
504 56 bertin
/* Remove bad PSF candidates */
505
  if (set->nsample>1)
506 80 bertin
    {
507
    psf_clean(psf, set, 0.1);
508 150 bertin
 
509 80 bertin
/*-- Make the basic PSF-model (3rd pass) */
510 112 bertin
//    NFPRINTF(OUTPUT,"Modeling the PSF (3/3)...");
511 80 bertin
    psf_make(psf, set, 0.05);
512
    psf_refine(psf, set);
513
    }
514 150 bertin
 
515 80 bertin
/* Remove bad PSF candidates */
516
  if (set->nsample>1)
517
    psf_clean(psf, set, 0.05);
518 150 bertin
 
519 56 bertin
  psf->samples_accepted = set->nsample;
520 150 bertin
 
521 56 bertin
/* Refine the PSF-model */
522 80 bertin
  psf_make(psf, set, prefs.prof_accuracy);
523 60 bertin
  psf_refine(psf, set);
524 150 bertin
 
525 56 bertin
/* Clip the PSF-model */
526
  psf_clip(psf);
527
 
528
/*-- Just check the Chi2 */
529
  psf->chi2 = set->nsample? psf_chi2(psf, set) : 0.0;
530
 
531
  return psf;
532
  }