public software.psfex

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

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

Line No. Rev Author Line
1 157 bertin
/*
2
*                               makeit.c
3 2 bertin
*
4 157 bertin
* Main loop.
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 4 bertin
time_t          thetime, thetime2;
61 191 rhl
void            write_error(char *msg1, char *msg2);
62 4 bertin
 
63 2 bertin
/********************************** makeit ***********************************/
64
/*
65
*/
66 73 bertin
void    makeit(void)
67 2 bertin
 
68
  {
69 177 bertin
   wcsstruct            *wcs;
70 176 bertin
   fieldstruct          **fields,
71
                        *field;
72 73 bertin
   psfstruct            **cpsf,
73
                        *psf;
74 79 bertin
   setstruct            *set, *set2;
75 73 bertin
   contextstruct        *context, *fullcontext;
76 56 bertin
   struct tm            *tm;
77 112 bertin
   char                 str[MAXCHAR];
78 73 bertin
   char                 **incatnames,
79
                        *pstr;
80 112 bertin
   float                **psfbasiss,
81
                        *psfsteps, *psfbasis, *basis,
82
                        psfstep, step;
83 109 bertin
   int                  c,i,p, ncat, ext, next, nmed, nbasis;
84 2 bertin
 
85 4 bertin
/* Install error logging */
86
  error_installfunc(write_error);
87
 
88 73 bertin
  incatnames = prefs.incat_name;
89
  ncat = prefs.ncat;
90
 
91 4 bertin
/* Processing start date and time */
92
  thetime = time(NULL);
93
  tm = localtime(&thetime);
94
  sprintf(prefs.sdate_start,"%04d-%02d-%02d",
95
        tm->tm_year+1900, tm->tm_mon+1, tm->tm_mday);
96
  sprintf(prefs.stime_start,"%02d:%02d:%02d",
97
        tm->tm_hour, tm->tm_min, tm->tm_sec);
98
 
99
  NFPRINTF(OUTPUT, "");
100
  QPRINTF(OUTPUT,
101
        "----- %s %s started on %s at %s with %d thread%s\n\n",
102
                BANNER,
103
                MYVERSION,
104
                prefs.sdate_start,
105
                prefs.stime_start,
106
                prefs.nthreads,
107
                prefs.nthreads>1? "s":"");
108
 
109 73 bertin
 
110
/* End here if no filename has been provided */
111
  if (!ncat)
112
    {
113
/*-- Processing end date and time */
114
    thetime2 = time(NULL);
115
    tm = localtime(&thetime2);
116
    sprintf(prefs.sdate_end,"%04d-%02d-%02d",
117
        tm->tm_year+1900, tm->tm_mon+1, tm->tm_mday);
118
    sprintf(prefs.stime_end,"%02d:%02d:%02d",
119
        tm->tm_hour, tm->tm_min, tm->tm_sec);
120
    prefs.time_diff = difftime(thetime2, thetime);
121
 
122
/*-- Write XML */
123
    if (prefs.xml_flag)
124
      {
125
      init_xml(0);
126
      write_xml(prefs.xml_name);
127
      end_xml();
128
      }
129
    return;
130
    }
131
 
132 75 bertin
/* Create an array of PSFs (one PSF for each extension) */
133 76 bertin
  QMALLOC(fields, fieldstruct *, ncat);
134 111 bertin
 
135
  NFPRINTF(OUTPUT, "");
136
  QPRINTF(OUTPUT, "----- %d input catalogues:\n", ncat);
137 75 bertin
  for (c=0; c<ncat; c++)
138 111 bertin
    {
139 76 bertin
    fields[c] = field_init(incatnames[c]);
140 111 bertin
    QPRINTF(OUTPUT, "%-20.20s:  \"%-16.16s\"  %3d extension%s %7d detection%s\n",
141
        fields[c]->rcatname, fields[c]->ident,
142
        fields[c]->next, fields[c]->next>1 ? "s":"",
143
        fields[c]->ndet, fields[c]->ndet>1 ? "s":"");
144
    }
145
  QPRINTF(OUTPUT, "\n");
146 75 bertin
 
147 8 bertin
  if (prefs.xml_flag)
148 84 bertin
    init_xml(ncat);
149 191 rhl
 
150
  makeit_body(fields, &context, &fullcontext, 1);
151
  next = fields[0]->next;
152 8 bertin
 
153 191 rhl
/* Write XML */
154
  if (prefs.xml_flag)
155 56 bertin
    {
156 191 rhl
    NFPRINTF(OUTPUT, "Writing XML file...");
157
    write_xml(prefs.xml_name);
158
    end_xml();
159 111 bertin
    }
160
 
161 73 bertin
/* Save result */
162
  for (c=0; c<ncat; c++)
163
    {
164 114 bertin
    sprintf(str, "Saving PSF model and metadata for %s...",
165
        fields[c]->rtcatname);
166 112 bertin
    NFPRINTF(OUTPUT, str);
167 73 bertin
/*-- Create a file name with a "PSF" extension */
168 128 bertin
    if (*prefs.psf_dir)
169
      {
170
      if ((pstr = strrchr(incatnames[c], '/')))
171
        pstr++;
172
      else
173
        pstr = incatnames[c];
174
      sprintf(str, "%s/%s", prefs.psf_dir, pstr);
175
      }
176
    else
177
      strcpy(str, incatnames[c]);
178 73 bertin
    if (!(pstr = strrchr(str, '.')))
179
      pstr = str+strlen(str);
180
    sprintf(pstr, "%s", prefs.psf_suffix);
181 76 bertin
    field_psfsave(fields[c], str);
182 73 bertin
/* Create homogenisation kernels */
183
    if (prefs.homobasis_type != HOMOBASIS_NONE)
184
      {
185
      for (ext=0; ext<next; ext++)
186 112 bertin
        {
187
        if (next>1)
188
          sprintf(str, "Computing homogenisation kernel for %s[%d/%d]...",
189
                fields[c]->rtcatname, ext+1, next);
190
        else
191
          sprintf(str, "Computing homogenisation kernel for %s...",
192
                fields[c]->rtcatname);
193
        NFPRINTF(OUTPUT, str);
194 128 bertin
        if (*prefs.homokernel_dir)
195
          {
196
          if ((pstr = strrchr(incatnames[c], '/')))
197
            pstr++;
198
          else
199
            pstr = incatnames[c];
200
          sprintf(str, "%s/%s", prefs.homokernel_dir, pstr);
201
          }
202
        else
203
          strcpy(str, incatnames[c]);
204 113 bertin
        if (!(pstr = strrchr(str, '.')))
205
          pstr = str+strlen(str);
206
        sprintf(pstr, "%s", prefs.homokernel_suffix);
207 76 bertin
        psf_homo(fields[c]->psf[ext], str, prefs.homopsf_params,
208 73 bertin
                prefs.homobasis_number, prefs.homobasis_scale, ext, next);
209 112 bertin
        }
210 73 bertin
      }
211 75 bertin
#ifdef HAVE_PLPLOT
212 78 bertin
/* Plot diagnostic maps for all catalogs */
213 119 bertin
    cplot_ellipticity(fields[c]);
214 76 bertin
    cplot_fwhm(fields[c]);
215 89 bertin
    cplot_moffatresi(fields[c]);
216
    cplot_asymresi(fields[c]);
217 119 bertin
    cplot_counts(fields[c]);
218
    cplot_countfrac(fields[c]);
219 150 bertin
    cplot_modchi2(fields[c]);
220
    cplot_modresi(fields[c]);
221 75 bertin
#endif
222 191 rhl
 
223 84 bertin
/*-- Update XML */
224
    if (prefs.xml_flag)
225
      update_xml(fields[c]);
226 73 bertin
    }
227
 
228 4 bertin
/* Processing end date and time */
229
  thetime2 = time(NULL);
230
  tm = localtime(&thetime2);
231
  sprintf(prefs.sdate_end,"%04d-%02d-%02d",
232
        tm->tm_year+1900, tm->tm_mon+1, tm->tm_mday);
233
  sprintf(prefs.stime_end,"%02d:%02d:%02d",
234
        tm->tm_hour, tm->tm_min, tm->tm_sec);
235
  prefs.time_diff = difftime(thetime2, thetime);
236
 
237 8 bertin
/* Free memory */
238 73 bertin
  for (c=0; c<ncat; c++)
239 76 bertin
    field_end(fields[c]);
240 100 bertin
  free(fields);
241 2 bertin
 
242 73 bertin
  if (context->npc)
243
    context_end(fullcontext);
244
  context_end(context);
245
 
246 2 bertin
  return;
247
  }
248
 
249 4 bertin
 
250 56 bertin
 
251
 
252 4 bertin
/****** write_error ********************************************************
253
PROTO   void    write_error(char *msg1, char *msg2)
254
PURPOSE Manage files in case of a catched error
255
INPUT   a character string,
256
        another character string
257
OUTPUT  RETURN_OK if everything went fine, RETURN_ERROR otherwise.
258
NOTES   -.
259
AUTHOR  E. Bertin (IAP)
260
VERSION 23/02/2007
261
 ***/
262
void    write_error(char *msg1, char *msg2)
263
  {
264
   char error[MAXCHAR];
265
 
266
  sprintf(error, "%s%s", msg1,msg2);
267
  if (prefs.xml_flag)
268
    write_xmlerror(prefs.xml_name, error);
269
  end_xml();
270
 
271
  return;
272
  }