public software.scamp

[/] [trunk/] [src/] [astrefcat.c] - Blame information for rev 334

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

Line No. Rev Author Line
1 247 bertin
/*
2
*                               astrefcat.c
3 2 bertin
*
4 247 bertin
* Manage astrometric reference catalogs (query and load).
5 2 bertin
*
6 247 bertin
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7 2 bertin
*
8 247 bertin
*       This file part of:      SCAMP
9 2 bertin
*
10 331 bertin
*       Copyright:              (C) 2002-2015 Emmanuel Bertin -- IAP/CNRS/UPMC
11 2 bertin
*
12 247 bertin
*       License:                GNU General Public License
13
*
14
*       SCAMP 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
*       SCAMP 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 SCAMP. If not, see <http://www.gnu.org/licenses/>.
24
*
25 334 bertin
*       Last modified:          29/06/2015
26 247 bertin
*
27
*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
28 2 bertin
 
29
#ifdef HAVE_CONFIG_H
30
#include        "config.h"
31
#endif
32
 
33
#include <math.h>
34
#include <netdb.h>
35
#include <stdio.h>
36
#include <stdlib.h>
37
#include <string.h>
38
#include <sys/socket.h>
39
#include <sys/types.h>
40
#include <sys/stat.h>
41
 
42
#include "define.h"
43
#include "globals.h"
44
#include "fits/fitscat.h"
45
#include "fitswcs.h"
46 331 bertin
#include "astrefcat.h"
47 2 bertin
#include "field.h"
48 112 bertin
#include "key.h"
49 2 bertin
#include "prefs.h"
50
#include "samples.h"
51 331 bertin
#include "url.h"
52 2 bertin
 
53
samplestruct    refsample;
54
keystruct       refkey[] = {
55
  {"X_WORLD", "Barycenter position along world x axis",
56
        &refsample.wcspos[0], H_FLOAT, T_DOUBLE, "%15e", "deg"},
57
  {"Y_WORLD", "Barycenter position along world y axis",
58
        &refsample.wcspos[1], H_FLOAT, T_DOUBLE, "%15e", "deg"},
59
  {"ERRA_WORLD", "World RMS position error along major axis",
60
        &refsample.wcsposerr[0], H_FLOAT, T_FLOAT, "%12e", "deg"},
61
  {"ERRB_WORLD", "World RMS position error along minor axis",
62
        &refsample.wcsposerr[1], H_FLOAT, T_FLOAT, "%12e", "deg"},
63
  {"MAG", "Generic magnitude",
64
        &refsample.mag, H_FLOAT, T_FLOAT, "%8.4f", "mag"},
65
  {"MAGERR", "Generic magnitude RMS error",
66
        &refsample.magerr, H_FLOAT, T_FLOAT, "%8.4f", "mag"},
67 305 bertin
  {"OBSDATE", "Observation date",
68
        &refsample.epoch, H_FLOAT, T_DOUBLE, "%13.8f", "yr"},
69 2 bertin
  {""},
70
  };
71
 
72 112 bertin
astrefstruct    astrefcat[] =
73
 {
74 321 bertin
  {"NONE", "", 0, 0, {""}},
75
  {"file", "", 12, 0, {"1","2","3","4","5","6","7","8","9","10","11","12",""}},
76
  {"USNO-A1.0", "pmm1", 2, 0, {"Bj", "Rf",""}, {"B", "R",""}},
77
  {"USNO-A2.0", "pmm2", 2, 0, {"Bj", "Rf",""}, {"B", "R",""}},
78
  {"USNO-B1.0", "usnob1", 3, 0, {"Bj", "Rf", "In",""}, {"B1", "R1", "I",""}},
79
  {"GSC-1.3", "gsc1.3", 1, 0, {"V",""}, {"V",""}},
80
  {"GSC-2.2", "gsc2.2", 4, 0, {"Bj", "V", "Rf", "In",""},
81
                        {"F", "J", "V", "N",""}},
82
  {"GSC-2.3", "gsc2.3", 6, 2, {"U", "B", "Bj", "V", "Rf", "In", ""},
83 233 bertin
                        {"U", "B", "F", "J", "V", "N", ""}},
84 321 bertin
  {"2MASS", "find2m", 3, 0, {"J", "H", "Ks",""}, {"J", "H", "K",""}},
85
  {"DENIS-3", "denis3", 3, 0, {"i", "J", "Ks",""}, {"I", "J", "K",""}},
86
  {"UCAC-1", "ucac1", 1, 0, {"R",""}, {"R",""}},
87
  {"UCAC-2", "ucac2", 1, 0, {"R",""}, {"R",""}},
88
  {"UCAC-3", "ucac3", 1, 0, {"R",""}, {"R",""}},
89
  {"UCAC-4", "ucac4", 1, 0, {"R",""}, {"R",""}},
90 331 bertin
  {"URAT-1", "urat1", 1, 0, {"R",""}, {"R",""}},
91 321 bertin
  {"SDSS-R3", "sdss3", 5, 2, {"u", "g", "r", "i", "z",""},
92
                        {"u", "g", "r", "i", "z",""}},
93
  {"SDSS-R5", "sdss5", 5, 2, {"u", "g", "r", "i", "z",""},
94
                        {"u", "g", "r", "i", "z",""}},
95
  {"SDSS-R6", "sdss6", 5, 2, {"u", "g", "r", "i", "z",""},
96
                        {"u", "g", "r", "i", "z",""}},
97
  {"SDSS-R7", "sdss7", 5, 2, {"u", "g", "r", "i", "z",""},
98
                        {"u", "g", "r", "i", "z",""}},
99
  {"SDSS-R8", "sdss8", 5, 2, {"u", "g", "r", "i", "z",""},
100
                        {"u", "g", "r", "i", "z",""}},
101
  {"SDSS-R9", "sdss9", 5, 2, {"u", "g", "r", "i", "z",""},
102
                        {"u", "g", "r", "i", "z",""}},
103
  {"NOMAD-1.0", "nomad1", 6, 2, {"B", "V", "R", "J", "H", "Ks",""},
104 218 bertin
                        {"B", "V", "R", "J", "H", "K",""}},
105 321 bertin
  {"PPMX", "ppmx", 7, 3, {"B", "V", "R", "Rf", "J", "H", "Ks",""},
106 233 bertin
                        {"B", "V", "R", "C", "J", "H", "K",""}},
107 321 bertin
  {"CMC-14", "cmc14", 4, 0, {"r", "J", "H", "Ks",""}, {"r", "J", "H", "K",""}},
108
  {"TYCHO-2", "tyc2", 2, 1, {"B", "V",""}, {"B", "V",""}},
109 112 bertin
  {""}
110
 };
111 2 bertin
 
112 306 bertin
static char     *astref_strncpy(char *dest, char *src, int n);
113 112 bertin
const char      astref_ctype[NAXIS][8]= {"RA---STG", "DEC--STG"};
114
 
115 2 bertin
const double    astref_crpix[NAXIS] = {8192.5, 8192.5},
116
                astref_cdelt[NAXIS] = {0.2*ARCSEC/DEG, 0.2*ARCSEC/DEG};
117
 
118
const int       astref_naxisn[NAXIS] = {16384, 16384};
119
 
120
/****** get_astreffield *********************************************************
121
PROTO   fieldstruct *get_astreffield(astrefenum refcat, double *wcspos,
122
                                int lng, int lat, int naxis, double maxradius)
123
PURPOSE Download reference catalog.
124
INPUT   Catalog name,
125
        Coordinate vector of the center,
126
        Longitude index,
127
        Latitude index,
128
        Number of axes (dimensions),
129
        Search radius (in degrees).
130
OUTPUT  Pointer to the reference field.
131
NOTES   Global preferences are used.
132
AUTHOR  E. Bertin (IAP)
133 334 bertin
VERSION 29/06/2015
134 2 bertin
*/
135
fieldstruct     *get_astreffield(astrefenum refcat, double *wcspos,
136
                                int lng, int lat, int naxis, double maxradius)
137
  {
138
   fieldstruct  *field,*tfield;
139
   setstruct    *set;
140
   samplestruct *sample;
141 331 bertin
   URL_FILE     *file;
142 2 bertin
   char         *ctype[NAXIS],
143 331 bertin
                url[MAXCHAR], maglimcmd[128],
144
                col[80], str[MAXCHAR],
145 306 bertin
                smag[MAX_BAND][32],smagerr[MAX_BAND][32],
146
                sprop[NAXIS][32],sproperr[NAXIS][32], sflag[4],
147 218 bertin
                *bandname, *cdsbandname, *catname,
148 287 bertin
                flag1,flag2, smode;
149 2 bertin
   double       poserr[NAXIS],prop[NAXIS],properr[NAXIS],
150 219 bertin
                mag[MAX_BAND],magerr[MAX_BAND], epoch,epocha,epochd,
151 306 bertin
                alpha,delta, dist, poserra,poserrb,poserrtheta, cpt,spt, temp;
152 331 bertin
   int          b,c,d,i,n,s, nsample,nsamplemax, nobs, class, band, nband;
153 2 bertin
 
154
/* One needs 2 angular coordinates here! */
155
  if (naxis<2)
156
    return NULL;
157
 
158
/* If these are not angular coordinates, file mode becomes mandatory */
159
  if (lng == lat && refcat!=ASTREFCAT_FILE)
160
    {
161
    warning("Cartesian coordinates found: ",
162
        "reference catalog switched to FILE mode");
163
    refcat = ASTREFCAT_FILE;
164
    }
165
 
166 112 bertin
  catname = (char *)astrefcat[(int)refcat].name;
167
  nband = astrefcat[(int)refcat].nband;
168
  if (!cistrcmp(prefs.astref_bandname, "DEFAULT", FIND_STRICT))
169
    band = astrefcat[(int)refcat].defband;
170
  else if (!cistrcmp(prefs.astref_bandname, "BLUEST", FIND_STRICT))
171
    band = 0;
172
  else if (!cistrcmp(prefs.astref_bandname, "REDDEST", FIND_STRICT))
173
    band = nband-1;
174
  else if ((band=findkeys(prefs.astref_bandname,
175 113 bertin
                astrefcat[(int)refcat].bandnames, FIND_STRICT))==RETURN_ERROR)
176 2 bertin
    {
177 112 bertin
    sprintf(str, "%s: no such band in astrometric reference catalog %s, use ",
178
        prefs.astref_bandname, catname);
179
    for (b=0; b<nband; b++)
180
      {
181
      if (b)
182
        strcat(str, b==nband-1? " or " : ", ");
183 113 bertin
      strcat(str, astrefcat[(int)refcat].bandnames[b]);
184 112 bertin
      }
185
    error(EXIT_FAILURE, "*Error*: ", str);
186 2 bertin
    }
187 112 bertin
  astrefcat[(int)refcat].band = band;
188 113 bertin
  bandname = astrefcat[(int)refcat].bandname
189
        = astrefcat[(int)refcat].bandnames[band];
190 218 bertin
  cdsbandname = astrefcat[(int)refcat].cdsbandnames[band];
191 2 bertin
 
192
/* Call the right catalog */
193 244 bertin
  if (refcat==ASTREFCAT_FILE)
194 2 bertin
    {
195 244 bertin
    field = NULL;
196
    for (c=0; c<prefs.nastref_name; c++)
197
      {
198
      if ((tfield=load_astreffield(prefs.astref_name[c], wcspos, lng,lat,
199 218 bertin
                naxis, maxradius, band, prefs.astref_maglim)))
200 244 bertin
        {
201 250 bertin
        NFPRINTF(OUTPUT, "");
202
        QPRINTF(OUTPUT, " %d astrometric references loaded from %s\n",
203 2 bertin
                tfield->set[0]->nsample, tfield->rfilename);
204 244 bertin
        if (field)
205
          {
206
          union_samples(tfield->set[0]->sample, field->set[0],
207
                tfield->set[0]->nsample,
208
                ASTREF_ASSOCRADIUS, UNION_WCS);
209
          field->nsample = field->set[0]->nsample;
210
          end_field(tfield);
211
          }
212
        else field=tfield;
213 2 bertin
        }
214 244 bertin
      }
215
    if (!field)
216 2 bertin
        error(EXIT_FAILURE,"*Error*: No appropriate FITS-LDAC astrometric ",
217
                        "reference catalog found");
218 244 bertin
    return field;
219
    }
220 2 bertin
 
221 321 bertin
/* Prepare mag limit section of the command line */
222
  if (prefs.astref_maglim[0]>-99.0 || prefs.astref_maglim[1]<99.0)
223 331 bertin
    sprintf(maglimcmd, "-lm%s&%f,%f&-m&10000000",
224 321 bertin
        cdsbandname, prefs.astref_maglim[0],prefs.astref_maglim[1]);
225
  else
226 331 bertin
    strcpy(maglimcmd, "-m&10000000");
227 2 bertin
 
228 331 bertin
/// Test all provided servers until one replies
229
  for (s=0; s<prefs.nref_server; s++) {
230
    sprintf(url,
231 334 bertin
        "http://%s/viz-bin/aserver.cgi?%s&-c&%12f%+12f&-r&%.8g&%s",
232 331 bertin
        prefs.ref_server[s],
233 321 bertin
        astrefcat[(int)refcat].cdsname,
234
        wcspos[lng], wcspos[lat],
235
        maxradius*DEG/ARCMIN,
236
        maglimcmd);
237
 
238 331 bertin
    sprintf(str,"Querying %s from %s for astrometric reference stars...",
239 244 bertin
        catname,
240 331 bertin
        prefs.ref_server[s]);
241
    NFPRINTF(OUTPUT, str);
242
    file = url_fopen(url, (long)prefs.ref_timeout[s]);
243
    // Test and skip the two first lines
244
    if (url_fgets(str, MAXCHAR, file) && url_fgets(str, MAXCHAR, file))
245
      break;
246
    else {
247
      url_fclose(file);
248
      if (s == prefs.nref_server-1)
249
        error(EXIT_FAILURE,"*Error*: No reply from catalog server ",
250
                prefs.ref_server[s]);
251
      else
252
        warning(prefs.ref_server[s], " connection timed out");
253
    }
254
  }
255 244 bertin
 
256 2 bertin
  set = init_set();
257
  prop[lng] = prop[lat] = properr[lng] = properr[lat] = 0.0;
258
  n = nsample = nsamplemax = 0;
259
 
260
/* Now examine each entry */
261 331 bertin
  while (url_fgets(str, MAXCHAR, file))
262 2 bertin
    if (*str != '#')
263
      {
264
      switch(refcat)
265
        {
266
        case ASTREFCAT_USNOA1:
267
        case ASTREFCAT_USNOA2:
268 306 bertin
          alpha = atof(astref_strncpy(col, str+14, 10));
269
          delta = atof(astref_strncpy(col, str+24, 10));
270
          mag[0] = atof(astref_strncpy(col, str+37, 4));
271
          mag[1] = atof(astref_strncpy(col, str+42, 4));
272
          epoch = atof(astref_strncpy(col, str+47, 8));
273 2 bertin
          poserr[lat] = poserr[lng] = (refcat==ASTREFCAT_USNOA1)?
274
                                USNOA1_POSERR : USNOA2_POSERR;
275 112 bertin
          magerr[0] = magerr[1] = (refcat==ASTREFCAT_USNOA1)?
276 2 bertin
                        USNOA1_BMAGERR : USNOA2_BMAGERR;
277
          break;
278 154 bertin
 
279 2 bertin
        case ASTREFCAT_USNOB1:
280
/*-------- Avoid spikes */
281 306 bertin
          strcpy(sflag, astref_strncpy(col, str+92, 3));
282 2 bertin
          if (sflag[0]=='s' || sflag[1]=='s' || sflag[2]=='s')
283
            continue;
284 306 bertin
          alpha = atof(astref_strncpy(col, str+26, 10));
285
          delta = atof(astref_strncpy(col, str+36, 10));
286
          poserr[lng] = atof(astref_strncpy(col, str+47, 3));
287 2 bertin
          if (poserr[lng] != 999.0)
288
            poserr[lng] *= MAS/DEG;
289
          else
290
            poserr[lng] = USNOB1_POSERR;
291 306 bertin
          poserr[lat] = atof(astref_strncpy(col, str+51, 3));
292 2 bertin
          if (poserr[lat] != 999.0)
293
            poserr[lat] *= MAS/DEG;
294
          else
295
            poserr[lat] = USNOB1_POSERR;
296 306 bertin
          epocha = epochd = atof(astref_strncpy(col, str+55, 6));
297
          prop[lng] = atof(astref_strncpy(col, str+62, 6))*MAS/DEG;
298
          prop[lat] = atof(astref_strncpy(col, str+69, 6))*MAS/DEG;
299
          properr[lng] = atof(astref_strncpy(col, str+78, 3))*MAS/DEG;
300
          properr[lat] = atof(astref_strncpy(col, str+82, 3))*MAS/DEG;
301
          strcpy(smag[0], astref_strncpy(col, str+96, 6));
302
          strcpy(smag[1], astref_strncpy(col, str+127, 6));
303
          strcpy(smag[2], astref_strncpy(col, str+158, 6));
304
          strcpy(smag[3], astref_strncpy(col, str+189, 6));
305
          strcpy(smag[4], astref_strncpy(col, str+220, 6));
306
          epoch = (properr[lng]==0.0 && properr[lat]==0.0)?
307
                0.5*(epocha+epochd) : 2000.0;
308 112 bertin
          for (b=0; b<5; b++)    /* 5, not nband!*/
309 306 bertin
            if (smag[b][4] == '-')
310 112 bertin
              mag[b] = magerr[b] = 99.0;
311 2 bertin
            else
312 112 bertin
              {
313
              mag[b] = atof(smag[b]);
314
              magerr[b] = USNOB1_BMAGERR;
315
              }
316
/*-------- Merge B1 and B2, R1 and R2, shift In */
317
          if (mag[0] > 98.0)
318
            mag[0] = mag[2];
319
          if (mag[2] > 98.0)
320
            mag[2] = mag[0];
321
          mag[0] = (mag[0] + mag[2])*0.5;
322 154 bertin
          if (magerr[0] > 98.0)
323
            magerr[0] = magerr[2];
324
          if (magerr[2] > 98.0)
325
            magerr[2] = magerr[0];
326
          magerr[0] = (magerr[0] + magerr[2])*0.5;
327 112 bertin
          if (mag[1] > 98.0)
328
            mag[1] = mag[3];
329
          if (mag[3] > 98.0)
330
            mag[3] = mag[1];
331
          mag[1] = (mag[1] + mag[3])*0.5;
332 154 bertin
          if (magerr[1] > 98.0)
333
            magerr[1] = magerr[3];
334
          if (magerr[3] > 98.0)
335
            magerr[3] = magerr[1];
336
          magerr[1] = (magerr[1] + magerr[3])*0.5;
337 112 bertin
          mag[2] = mag[4];
338 154 bertin
          magerr[2] = magerr[4];
339 2 bertin
          break;
340 154 bertin
 
341 2 bertin
        case ASTREFCAT_GSC1:
342 306 bertin
          alpha = atof(astref_strncpy(col, str+11, 9));
343
          delta = atof(astref_strncpy(col, str+21, 9));
344
          poserr[lng] = atof(astref_strncpy(col, str+32, 4))*ARCSEC/DEG;
345 2 bertin
          poserr[lat] = poserr[lng];
346 306 bertin
          mag[0] = atof(astref_strncpy(col, str+37, 5));
347
          magerr[0] = atof(astref_strncpy(col, str+43, 4));
348 2 bertin
          epoch = 1960.0;
349
          break;
350
 
351 233 bertin
        case ASTREFCAT_GSC22:
352 306 bertin
          class = atoi(astref_strncpy(col, str+106,4));
353 233 bertin
          if (class==5)
354
            continue;
355 306 bertin
          alpha = atof(astref_strncpy(col, str+15, 10));
356
          delta = atof(astref_strncpy(col, str+26, 10));
357
          epoch = atof(astref_strncpy(col, str+37, 8));
358
          poserr[lng] = atof(astref_strncpy(col, str+46, 5))*ARCSEC/DEG;
359
          poserr[lat] = atof(astref_strncpy(col, str+52, 5))*ARCSEC/DEG;
360
          strcpy(smag[0], astref_strncpy(col, str+58, 5));
361
          strcpy(smagerr[0], astref_strncpy(col, str+65, 4));
362
          strcpy(smag[1], astref_strncpy(col, str+70, 5));
363
          strcpy(smagerr[1], astref_strncpy(col, str+77, 4));
364
          strcpy(smag[2], astref_strncpy(col, str+82, 5));
365
          strcpy(smagerr[2], astref_strncpy(col, str+89, 4));
366
          strcpy(smag[3], astref_strncpy(col, str+94, 5));
367
          strcpy(smagerr[3], astref_strncpy(col, str+101, 4));
368 112 bertin
          for (b=0; b<nband; b++)
369 306 bertin
            if (smag[b][4] == '-')
370 112 bertin
              mag[b] = magerr[b] = 99.0;
371
            else
372
              {
373
              mag[b] = atof(smag[b]);
374
              magerr[b] = atof(smagerr[b]);
375
              }
376 113 bertin
/*-------- Swap Bj with Rf and Rf with V to sort passbands by lambda */
377
          temp=mag[0];mag[0]=mag[1];mag[1]=temp;
378
          temp=magerr[0];magerr[0]=magerr[1];magerr[1]=temp;
379
          temp=mag[1];mag[1]=mag[2];mag[2]=temp;
380
          temp=magerr[1];magerr[1]=magerr[2];magerr[2]=temp;
381 2 bertin
          break;
382
 
383 233 bertin
        case ASTREFCAT_GSC23:
384 306 bertin
          class = atoi(astref_strncpy(col, str+172,2));
385 233 bertin
          if (class==5)
386
            continue;
387 306 bertin
          alpha = atof(astref_strncpy(col, str+33, 10));
388
          delta = atof(astref_strncpy(col, str+43, 10));
389
          poserr[lng] = atof(astref_strncpy(col, str+54, 5))*ARCSEC/DEG;
390
          poserr[lat] = atof(astref_strncpy(col, str+60, 5))*ARCSEC/DEG;
391
          epoch = atof(astref_strncpy(col, str+66, 8));
392
          strcpy(smag[0], astref_strncpy(col, str+76, 5));
393
          strcpy(smagerr[0], astref_strncpy(col, str+83, 4));
394
          strcpy(smag[1], astref_strncpy(col, str+92, 5));
395
          strcpy(smagerr[1], astref_strncpy(col, str+99, 4));
396
          strcpy(smag[2], astref_strncpy(col, str+108, 5));
397
          strcpy(smagerr[2], astref_strncpy(col, str+115, 4));
398
          strcpy(smag[3], astref_strncpy(col, str+124, 5));
399
          strcpy(smagerr[3], astref_strncpy(col, str+131, 4));
400
          strcpy(smag[4], astref_strncpy(col, str+140, 5));
401
          strcpy(smagerr[4], astref_strncpy(col, str+147, 4));
402
          strcpy(smag[5], astref_strncpy(col, str+156, 5));
403
          strcpy(smagerr[5], astref_strncpy(col, str+163, 4));
404 233 bertin
          for (b=0; b<nband; b++)
405 306 bertin
            if (smag[b][4] == '-')
406 233 bertin
              mag[b] = magerr[b] = 99.0;
407
            else
408
              {
409
              mag[b] = atof(smag[b]);
410
              magerr[b] = atof(smagerr[b]);
411
              }
412
/*-------- Swap to sort passbands by lambda */
413
          temp=mag[0];mag[0]=mag[4];mag[4]=temp;
414
          temp=magerr[0];magerr[0]=magerr[4];magerr[4]=temp;
415
          temp=mag[1];mag[1]=mag[5];mag[5]=temp;
416
          temp=magerr[1];magerr[1]=magerr[5];magerr[5]=temp;
417
          temp=mag[3];mag[3]=mag[5];mag[5]=temp;
418
          temp=magerr[3];magerr[3]=magerr[5];magerr[5]=temp;
419
          break;
420
 
421 2 bertin
        case ASTREFCAT_2MASS:
422 287 bertin
/*-------- Avoid contaminated observations */
423 306 bertin
          strcpy(sflag, astref_strncpy(col, str+156, 3));
424 287 bertin
          if (sflag[0]!='0' || sflag[1]!='0' || sflag[2]!='0')
425
            continue;
426 306 bertin
          alpha = atof(astref_strncpy(col, str+0, 10));
427
          delta = atof(astref_strncpy(col, str+11, 10));
428
          poserra = atof(astref_strncpy(col, str+22, 4));
429
          poserrb = atof(astref_strncpy(col, str+27, 4));
430
          poserrtheta = atof(astref_strncpy(col, str+32, 3));
431
          strcpy(smag[0], astref_strncpy(col, str+54, 6));
432
          strcpy(smagerr[0], astref_strncpy(col, str+67, 5));
433
          strcpy(smag[1], astref_strncpy(col, str+84, 6));
434
          strcpy(smagerr[1], astref_strncpy(col, str+97, 5));
435
          strcpy(smag[2], astref_strncpy(col, str+114, 6));
436
          strcpy(smagerr[2], astref_strncpy(col, str+127, 5));
437
          epoch = atof(astref_strncpy(col, str+243, 12));
438 112 bertin
          for (b=0; b<nband; b++)
439 306 bertin
            if (smag[b][4] == '-' || smagerr[b][4] == '-')
440 112 bertin
              mag[b] = magerr[b] = 99.0;
441
            else
442
              {
443
              mag[b] = atof(smag[b]);
444
              magerr[b] = atof(smagerr[b]);
445
              }
446 306 bertin
/*--------- Projet uncertainties on alpha and delta axes */
447
          cpt = cos(poserrtheta*DEG);
448
          spt = sin(poserrtheta*DEG);
449
          poserr[lng] = sqrt(spt*spt*poserra*poserra+cpt*cpt*poserrb*poserrb)
450
                                *ARCSEC/DEG;
451
          poserr[lat] = sqrt(cpt*cpt*poserra*poserra+spt*spt*poserrb*poserrb)
452
                                *ARCSEC/DEG;
453
 
454 305 bertin
/*-------- Convert JDs to epoch */
455 308 bertin
          epoch = 2000.0 + (epoch - JD2000)/365.25;
456 2 bertin
          break;
457
 
458 171 bertin
        case ASTREFCAT_DENIS3:
459 306 bertin
          alpha = atof(astref_strncpy(col, str+31, 10));
460
          delta = atof(astref_strncpy(col, str+42, 10));
461
          strcpy(smag[0], astref_strncpy(col, str+53, 6));
462
          strcpy(smagerr[0], astref_strncpy(col, str+60, 5));
463
          strcpy(smag[1], astref_strncpy(col, str+66, 6));
464
          strcpy(smagerr[1], astref_strncpy(col, str+73, 5));
465
          strcpy(smag[2], astref_strncpy(col, str+79, 6));
466
          strcpy(smagerr[2], astref_strncpy(col, str+86, 5));
467
          epoch = atof(astref_strncpy(col, str+448, 14));
468 171 bertin
          for (b=0; b<nband; b++)
469
            if (smag[b][2] == ' ' || smagerr[b][2] == ' ')
470
              mag[b] = magerr[b] = 99.0;
471
            else
472
              {
473
              mag[b] = atof(smag[b]);
474
              magerr[b] = atof(smagerr[b]);
475
              }
476 306 bertin
          poserr[lat] = poserr[lng] = DENIS3_POSERR;
477 305 bertin
/*-------- Convert JDs to epoch */
478 308 bertin
          epoch = 2000.0 + (epoch - JD2000)/365.25;
479 171 bertin
          break;
480
 
481 2 bertin
        case ASTREFCAT_UCAC1:
482
/*-------- Avoid poor observations */
483 306 bertin
          nobs = atoi(astref_strncpy(col, str+46, 2));
484 2 bertin
          if (nobs<2)
485
            continue;
486 306 bertin
          alpha = atof(astref_strncpy(col, str+9, 11));
487
          delta = atof(astref_strncpy(col, str+20, 11));
488
          poserr[lng] = atof(astref_strncpy(col, str+32, 3))*MAS/DEG;
489
          poserr[lat] = atof(astref_strncpy(col, str+36, 3))*MAS/DEG;
490
          mag[0] = atof(astref_strncpy(col, str+40, 5));
491
          magerr[0] = UCAC_MAGERR;
492
          epoch = atof(astref_strncpy(col, str+49, 8));
493
          prop[lng] = atof(astref_strncpy(col, str+58, 8))*MAS/DEG;
494
          prop[lat] = atof(astref_strncpy(col, str+67, 8))*MAS/DEG;
495
          properr[lng] = atof(astref_strncpy(col, str+76, 4))*MAS/DEG;
496
          properr[lat] = atof(astref_strncpy(col, str+81, 4))*MAS/DEG;
497 2 bertin
          break;
498
 
499
        case ASTREFCAT_UCAC2:
500
/*-------- Avoid poor observations */
501 306 bertin
          nobs = atoi(astref_strncpy(col, str+50, 2));
502 2 bertin
          if (nobs<2)
503
            continue;
504 306 bertin
          alpha = atof(astref_strncpy(col, str+9, 11));
505
          delta = atof(astref_strncpy(col, str+20, 11));
506
          poserr[lng] = atof(astref_strncpy(col, str+32, 3))*MAS/DEG;
507
          poserr[lat] = atof(astref_strncpy(col, str+36, 3))*MAS/DEG;
508
          mag[0] = atof(astref_strncpy(col, str+44, 5));
509
          magerr[0] = UCAC_MAGERR;
510
          epocha = atof(astref_strncpy(col, str+60, 8));
511
          epochd = atof(astref_strncpy(col, str+69, 8));
512
          epoch = 2000.0;
513
          prop[lng] = atof(astref_strncpy(col, str+78, 8))*MAS/DEG;
514
          prop[lat] = atof(astref_strncpy(col, str+87, 8))*MAS/DEG;
515
          properr[lng] = atof(astref_strncpy(col, str+96, 4))*MAS/DEG;
516
          properr[lat] = atof(astref_strncpy(col, str+101, 4))*MAS/DEG;
517 2 bertin
          break;
518
 
519 219 bertin
        case ASTREFCAT_UCAC3:
520
/*-------- Avoid poor observations */
521 306 bertin
          nobs = atoi(astref_strncpy(col, str+88,3));
522 219 bertin
          if (nobs<2)
523
            continue;
524 306 bertin
          alpha = atof(astref_strncpy(col, str+11, 11));
525
          delta = atof(astref_strncpy(col, str+22, 11));
526
          poserr[lng] = poserr[lat]=atof(astref_strncpy(col,str+42,4))*MAS/DEG;
527
          epocha = atof(astref_strncpy(col, str+47, 7));
528
          epochd = atof(astref_strncpy(col, str+55, 7));
529
          mag[0] = atof(astref_strncpy(col, str+63, 6));
530
          strcpy(smagerr[0], astref_strncpy(col, str+77, 5));
531
          strcpy(sprop[lng], astref_strncpy(col, str+104, 8));
532
          strcpy(sprop[lat], astref_strncpy(col, str+113, 8));
533
          strcpy(sproperr[lng], astref_strncpy(col, str+122, 4));
534
          strcpy(sproperr[lat], astref_strncpy(col, str+127, 4));
535
          magerr[0] = (smagerr[0][2]=='-')? 0.9 : atof(smagerr[0]);
536
          if (sprop[lng][6]== ' ' || sprop[lat][6]== ' ')
537
            {
538
            prop[lng] = prop[lat] = properr[lng] = properr[lat] = 0.0;
539
            epoch = 0.5*(epocha+epochd);
540
            }
541
          else
542
            {
543
            prop[lng] = atof(sprop[lng])*MAS/DEG;
544
            prop[lat] = atof(sprop[lat])*MAS/DEG;
545
            properr[lng] = (sproperr[lng][2]=='-'?
546
                                1000.0 : atof(sproperr[lng])) * MAS/DEG;
547
            properr[lat] = (sproperr[lat][2]=='-'?
548
                                1000.0 : atof(sproperr[lat])) * MAS/DEG;
549
            epoch = 2000.0;
550
            }
551 219 bertin
          break;
552
 
553 306 bertin
        case ASTREFCAT_UCAC4:
554
/*-------- Avoid poor observations */
555
          nobs = atoi(astref_strncpy(col, str+89, 3));
556
          if (nobs<2)
557
            continue;
558
          alpha = atof(astref_strncpy(col, str+11, 11));
559
          delta = atof(astref_strncpy(col, str+22, 11));
560
          poserr[lng] = poserr[lat]=atof(astref_strncpy(col,str+42,4))*MAS/DEG;
561
          epocha = atof(astref_strncpy(col, str+47, 7));
562
          epochd = atof(astref_strncpy(col, str+55, 7));
563
          mag[0] = atof(astref_strncpy(col, str+63, 6));
564
          strcpy(smagerr[0], astref_strncpy(col, str+77, 4));
565
          strcpy(sprop[lng], astref_strncpy(col, str+101, 8));
566
          strcpy(sprop[lat], astref_strncpy(col, str+110, 8));
567
          strcpy(sproperr[lng], astref_strncpy(col, str+119, 4));
568
          strcpy(sproperr[lat], astref_strncpy(col, str+124, 4));
569
          magerr[0] = (smagerr[0][2]=='-')? 0.9 : atof(smagerr[0]);
570
          if (sprop[lng][6]== ' ' || sprop[lat][6]== ' ')
571
            {
572
            prop[lng] = prop[lat] = properr[lng] = properr[lat] = 0.0;
573
            epoch = 0.5*(epocha+epochd);
574
            }
575
          else
576
            {
577
            prop[lng] = atof(sprop[lng])*MAS/DEG;
578
            prop[lat] = atof(sprop[lat])*MAS/DEG;
579
            properr[lng] = (sproperr[lng][2]=='-'?
580
                                1000.0 : atof(sproperr[lng])) * MAS/DEG;
581
            properr[lat] = (sproperr[lat][2]=='-'?
582
                                1000.0 : atof(sproperr[lat])) * MAS/DEG;
583
            epoch = 2000.0;
584
            }
585
          break;
586
 
587 331 bertin
        case ASTREFCAT_URAT1:
588
/*-------- Avoid poor observations */
589
          nobs = atoi(astref_strncpy(col, str+42, 2));
590
          if (nobs<2)
591
            continue;
592
          alpha = atof(astref_strncpy(col, str+11, 11));
593
          delta = atof(astref_strncpy(col, str+22, 11));
594
          poserr[lng] = atof(astref_strncpy(col,str+34,3)) * MAS/DEG;
595
          poserr[lat] = atof(astref_strncpy(col,str+38,3)) * MAS/DEG;
596
          epoch = atof(astref_strncpy(col, str+48, 8));
597
          mag[0] = atof(astref_strncpy(col, str+57, 6));
598
          magerr[0] = atof(astref_strncpy(col, str+64, 5));
599
          prop[lng] = atof(astref_strncpy(col, str+91, 8)) * MAS/DEG;
600
          prop[lat] = atof(astref_strncpy(col, str+100, 8)) * MAS/DEG;
601
          properr[lng] = properr[lat] = atof(astref_strncpy(col, str+109, 4))
602
                                        * MAS/DEG;
603
          break;
604
 
605 10 bertin
        case ASTREFCAT_SDSSR3:
606 11 bertin
/*-------- Avoid missing or poor observations */
607 306 bertin
          nobs = atoi(astref_strncpy(col, str+56,1));
608 11 bertin
          if (nobs<2 || nobs>3)
609 10 bertin
            continue;
610 306 bertin
          alpha = atof(astref_strncpy(col, str+25, 10));
611
          delta = atof(astref_strncpy(col, str+35, 10));
612
          epoch = atof(astref_strncpy(col, str+46, 9));
613
          mag[0] = atof(astref_strncpy(col, str+58, 6));
614
          magerr[0] = atof(astref_strncpy(col, str+65, 5));
615
          mag[1] = atof(astref_strncpy(col, str+71, 6));
616
          magerr[1] = atof(astref_strncpy(col, str+78, 5));
617
          mag[2] = atof(astref_strncpy(col, str+84, 6));
618
          magerr[2] = atof(astref_strncpy(col, str+91, 5));
619
          mag[3] = atof(astref_strncpy(col, str+97, 6));
620
          magerr[3] = atof(astref_strncpy(col, str+104, 5));
621
          mag[4] = atof(astref_strncpy(col, str+110, 6));
622
          magerr[4] = atof(astref_strncpy(col, str+117, 5));
623 10 bertin
          poserr[lng] = poserr[lat] = SDSSR3_POSERR;
624
          break;
625
 
626 135 bertin
        case ASTREFCAT_SDSSR5:
627 306 bertin
/*-------- Avoid missing or poor observations, and secondary detections */
628
         smode = str[0];
629
         nobs = atoi(astref_strncpy(col, str+76, 1));
630
         if (nobs<2 || nobs>3 || smode=='2')
631
            continue;
632
          alpha = atof(astref_strncpy(col, str+27, 10));
633
          delta = atof(astref_strncpy(col, str+37, 10));
634
          poserr[lng] = atof(astref_strncpy(col, str+48, 5))*ARCSEC/DEG;
635
          poserr[lat] = atof(astref_strncpy(col, str+54, 5))*ARCSEC/DEG;
636
          epoch = atof(astref_strncpy(col, str+66, 9));
637
          mag[0] = atof(astref_strncpy(col, str+78, 6));
638
          magerr[0] = atof(astref_strncpy(col, str+85, 5));
639
          mag[1] = atof(astref_strncpy(col, str+91, 6));
640
          magerr[1] = atof(astref_strncpy(col, str+98, 5));
641
          mag[2] = atof(astref_strncpy(col, str+104, 6));
642
          magerr[2] = atof(astref_strncpy(col, str+111, 5));
643
          mag[3] = atof(astref_strncpy(col, str+117, 6));
644
          magerr[3] = atof(astref_strncpy(col, str+124, 5));
645
          mag[4] = atof(astref_strncpy(col, str+130, 6));
646
          magerr[4] = atof(astref_strncpy(col, str+137, 5));
647
          break;
648 171 bertin
        case ASTREFCAT_SDSSR6:
649 218 bertin
        case ASTREFCAT_SDSSR7:
650 287 bertin
/*-------- Avoid missing or poor observations, and secondary detections */
651 306 bertin
          smode = str[0];
652
          nobs = atoi(astref_strncpy(col, str+76, 1));
653 287 bertin
          if (nobs<2 || nobs>3 || smode=='2')
654 171 bertin
            continue;
655 306 bertin
          alpha = atof(astref_strncpy(col, str+27, 10));
656
          delta = atof(astref_strncpy(col, str+37, 10));
657
          poserr[lng] = atof(astref_strncpy(col, str+48, 5))*ARCSEC/DEG;
658
          poserr[lat] = atof(astref_strncpy(col, str+54, 5))*ARCSEC/DEG;
659
          epoch = atof(astref_strncpy(col, str+66, 9));
660
          mag[0] = atof(astref_strncpy(col, str+93, 6));
661
          magerr[0] = atof(astref_strncpy(col, str+100, 5));
662
          mag[1] = atof(astref_strncpy(col, str+106, 6));
663
          magerr[1] = atof(astref_strncpy(col, str+113, 5));
664
          mag[2] = atof(astref_strncpy(col, str+119, 6));
665
          magerr[2] = atof(astref_strncpy(col, str+126, 5));
666
          mag[3] = atof(astref_strncpy(col, str+132, 6));
667
          magerr[3] = atof(astref_strncpy(col, str+139, 5));
668
          mag[4] = atof(astref_strncpy(col, str+145, 6));
669
          magerr[4] = atof(astref_strncpy(col, str+152, 5));
670 171 bertin
          break;
671
 
672 286 bertin
        case ASTREFCAT_SDSSR8:
673 287 bertin
/*-------- Avoid missing or poor observations, and secondary detections */
674 306 bertin
          smode = str[0];
675
          nobs = atoi(astref_strncpy(col, str+70, 1));
676 287 bertin
          if (nobs<2 || nobs>3 || smode=='2')
677 286 bertin
            continue;
678 306 bertin
          alpha = atof(astref_strncpy(col, str+27, 10));
679
          delta = atof(astref_strncpy(col, str+37, 10));
680
          poserr[lng] = atof(astref_strncpy(col, str+48, 5))*ARCSEC/DEG;
681
          poserr[lat] = atof(astref_strncpy(col, str+54, 5))*ARCSEC/DEG;
682
          epoch = atof(astref_strncpy(col, str+60, 9));
683
          mag[0] = atof(astref_strncpy(col, str+87, 6));
684
          magerr[0] = atof(astref_strncpy(col, str+94, 5));
685
          mag[1] = atof(astref_strncpy(col, str+100, 6));
686
          magerr[1] = atof(astref_strncpy(col, str+107, 5));
687
          mag[2] = atof(astref_strncpy(col, str+113, 6));
688
          magerr[2] = atof(astref_strncpy(col, str+120, 5));
689
          mag[3] = atof(astref_strncpy(col, str+126, 6));
690
          magerr[3] = atof(astref_strncpy(col, str+133, 5));
691
          mag[4] = atof(astref_strncpy(col, str+139, 6));
692
          magerr[4] = atof(astref_strncpy(col, str+146, 5));
693 286 bertin
          break;
694
 
695 321 bertin
        case ASTREFCAT_SDSSR9:
696
/*-------- Avoid missing or poor observations, and secondary detections */
697
          smode = str[0];
698
          nobs = atoi(astref_strncpy(col, str+69, 1));
699
          if (nobs<2 || nobs>3 || smode=='2')
700
            continue;
701
          alpha = atof(astref_strncpy(col, str+26, 10));
702
          delta = atof(astref_strncpy(col, str+36, 10));
703
          poserr[lng] = atof(astref_strncpy(col, str+47, 5))*ARCSEC/DEG;
704
          poserr[lat] = atof(astref_strncpy(col, str+53, 5))*ARCSEC/DEG;
705
          epoch = atof(astref_strncpy(col, str+59, 9));
706
          mag[0] = atof(astref_strncpy(col, str+71, 6));
707
          magerr[0] = atof(astref_strncpy(col, str+78, 5));
708
          mag[1] = atof(astref_strncpy(col, str+84, 6));
709
          magerr[1] = atof(astref_strncpy(col, str+91, 5));
710
          mag[2] = atof(astref_strncpy(col, str+97, 6));
711
          magerr[2] = atof(astref_strncpy(col, str+104, 5));
712
          mag[3] = atof(astref_strncpy(col, str+110, 6));
713
          magerr[3] = atof(astref_strncpy(col, str+117, 5));
714
          mag[4] = atof(astref_strncpy(col, str+123, 6));
715
          magerr[4] = atof(astref_strncpy(col, str+130, 5));
716
          break;
717
 
718 159 bertin
        case ASTREFCAT_NOMAD1:
719 306 bertin
          alpha = atof(astref_strncpy(col, str+19, 11));
720
          delta = atof(astref_strncpy(col, str+30, 11));
721
          poserr[lng] = atof(astref_strncpy(col, str+44, 4))*ARCSEC/DEG;
722
          poserr[lat] = atof(astref_strncpy(col, str+50, 4))*ARCSEC/DEG;
723
          epocha = atof(astref_strncpy(col, str+54, 6));
724
          epochd = atof(astref_strncpy(col, str+61, 6));
725
          prop[lng] = atof(astref_strncpy(col, str+68, 8))*MAS/DEG;
726
          prop[lat] = atof(astref_strncpy(col, str+77, 8))*MAS/DEG;
727
          properr[lng] = atof(astref_strncpy(col, str+86, 5))*MAS/DEG;
728
          properr[lat] = atof(astref_strncpy(col, str+92, 5))*MAS/DEG;
729
          strcpy(smag[0], astref_strncpy(col, str+98, 6));
730
          strcpy(smag[1], astref_strncpy(col, str+106, 6));
731
          strcpy(smag[2], astref_strncpy(col, str+114, 6));
732
          strcpy(smag[3], astref_strncpy(col, str+122, 6));
733
          strcpy(smag[4], astref_strncpy(col, str+129, 6));
734
          strcpy(smag[5], astref_strncpy(col, str+136, 6));
735
          epoch = (properr[lng]==0.0 && properr[lat]==0.0)?
736
                0.5*(epocha+epochd) : 2000.0;
737 159 bertin
          for (b=0; b<nband; b++)
738 306 bertin
            if (smag[b][2] == '-')
739 159 bertin
              mag[b] = magerr[b] = 99.0;
740
            else
741
              {
742
              mag[b] = atof(smag[b]);
743
              magerr[b] = NOMAD1_MAGERR;
744
              }
745
          break;
746
 
747 233 bertin
        case ASTREFCAT_PPMX:
748 306 bertin
          alpha = atof(astref_strncpy(col, str+17, 10));
749
          delta = atof(astref_strncpy(col, str+27, 10));
750
          prop[lng] = atof(astref_strncpy(col, str+38, 8))*MAS/DEG;
751
          prop[lat] = atof(astref_strncpy(col, str+47, 8))*MAS/DEG;
752
          epocha = atof(astref_strncpy(col, str+56, 7));
753
          epochd = atof(astref_strncpy(col, str+64, 7));
754
          poserr[lng] = atof(astref_strncpy(col, str+72, 3))*MAS/DEG;
755
          poserr[lat] = atof(astref_strncpy(col, str+76, 3))*MAS/DEG;
756
          properr[lng] = atof(astref_strncpy(col, str+80, 4))*MAS/DEG;
757
          properr[lat] = atof(astref_strncpy(col, str+85, 4))*MAS/DEG;
758
          strcpy(smag[3], astref_strncpy(col, str+90, 6));
759
          strcpy(smag[2], astref_strncpy(col, str+97, 6));
760
          strcpy(smag[0], astref_strncpy(col, str+104, 6));
761
          strcpy(smagerr[0], astref_strncpy(col, str+111, 3));
762
          strcpy(smag[1], astref_strncpy(col, str+115, 6));
763
          strcpy(smagerr[1], astref_strncpy(col, str+122, 3));
764
          strcpy(smag[4], astref_strncpy(col, str+126, 6));
765
          strcpy(smagerr[4], astref_strncpy(col, str+133, 3));
766
          strcpy(smag[5], astref_strncpy(col, str+137, 6));
767
          strcpy(smagerr[5], astref_strncpy(col, str+144, 3));
768
          strcpy(smag[6], astref_strncpy(col, str+148, 6));
769
          strcpy(smagerr[6], astref_strncpy(col, str+155, 3));
770
          epoch = 2000.0;
771 233 bertin
          for (b=0; b<nband; b++)
772 306 bertin
            if (smag[b][2] == '-')
773 233 bertin
              mag[b] = magerr[b] = 99.0;
774
            else
775
              {
776
              mag[b] = atof(smag[b]);
777
              magerr[b] = (b!=2 && b!=3)? atof(smagerr[b])*0.001 : GSC_MAGERR;
778
              }
779
          break;
780
 
781 308 bertin
        case ASTREFCAT_CMC14:
782
          alpha = sextodegal(astref_strncpy(col, str+17, 13));
783
          delta = sextodegde(astref_strncpy(col, str+30, 13));
784
          epoch = atof(astref_strncpy(col, str+44, 4));
785
          strcpy(smag[0], astref_strncpy(col, str+49, 6));
786
          poserr[lng] = atof(astref_strncpy(col, str+66, 5))*ARCSEC/DEG;
787
          poserr[lat] = atof(astref_strncpy(col, str+72, 5))*ARCSEC/DEG;
788
          strcpy(smagerr[0], astref_strncpy(col, str+78, 5));
789
          strcpy(smag[1], astref_strncpy(col, str+84, 6));
790
          strcpy(smag[2], astref_strncpy(col, str+91, 6));
791
          strcpy(smag[3], astref_strncpy(col, str+98, 6));
792
/*-------- Convert JDs to epoch */
793
          epoch = 2000.0 + (epoch - (JD2000-2451263.5))/365.25;
794
          for (b=0; b<nband; b++)
795
            if (smag[b][2] == '-')
796
              mag[b] = magerr[b] = 99.0;
797
            else
798
              {
799
              mag[b] = atof(smag[b]);
800 321 bertin
              magerr[b] = (b==0) ? atof(smagerr[b]) : DEFAULT_MAGERR;
801 308 bertin
              }
802
          break;
803 321 bertin
        case ASTREFCAT_TYCHO2:
804 322 bertin
/*-------- Reject (fainter) sources without proper motions */
805
          if (*astref_strncpy(col, str+27, 1) == ' ')
806
            continue;
807 321 bertin
          alpha = atof(astref_strncpy(col, str+24, 12));
808
          delta = atof(astref_strncpy(col, str+37, 12));
809
          prop[lng] = atof(astref_strncpy(col, str+50, 7))*MAS/DEG;
810
          prop[lat] = atof(astref_strncpy(col, str+58, 7))*MAS/DEG;
811
          poserr[lng] = atof(astref_strncpy(col, str+66, 3))*MAS/DEG;
812
          poserr[lat] = atof(astref_strncpy(col, str+70, 3))*MAS/DEG;
813
          properr[lng] = atof(astref_strncpy(col, str+74, 4))*MAS/DEG;
814
          properr[lat] = atof(astref_strncpy(col, str+79, 4))*MAS/DEG;
815
          strcpy(smag[0], astref_strncpy(col, str+123, 6));
816
          strcpy(smagerr[0], astref_strncpy(col, str+130, 6));
817
          strcpy(smag[1], astref_strncpy(col, str+137, 6));
818
          strcpy(smagerr[0], astref_strncpy(col, str+144, 6));
819
          epoch = 2000.0;
820
          for (b=0; b<nband; b++)
821
            if (smag[b][2] == ' ')
822
              mag[b] = magerr[b] = 99.0;
823
            else
824
              {
825
              mag[b] = atof(smag[b]);
826
              magerr[b] = (smagerr[b][2] != ' ') ?
827
                        atof(smagerr[b]) : DEFAULT_MAGERR;
828
              }
829
          break;
830 2 bertin
        case ASTREFCAT_NONE:
831
        default:
832
          break;
833
        }
834
 
835
      if (!(n%10000))
836
        {
837 250 bertin
        sprintf(str,"%-.36s (%s band) : %d / %d references stored",
838
                catname,bandname, nsample,n);
839 2 bertin
        NFPRINTF(OUTPUT, str);
840
        }
841
/*------ Apply some selection over flags, fluxes,.. */
842
      if (1)
843
        {
844
/*------ Allocate memory for the first shipment */
845
        if (!set->sample)
846
          {
847
          nsample = 0;
848
          nsamplemax = LSAMPLE_DEFSIZE;
849
          malloc_samples(set, nsamplemax);
850
          }
851
 
852
/*------ Increase storage space to receive new candidates if needed */
853
        if (nsample>=nsamplemax)
854
          {
855
           int  nadd=(int)(1.62*nsamplemax);
856
          nsamplemax = nadd>nsamplemax?nadd:nsamplemax+1;
857
          realloc_samples(set, nsamplemax);
858
          }
859
 
860
        sample = set->sample + nsample;
861 313 bertin
 
862 112 bertin
        if (mag[band] < 98.0)
863
          {
864
          sample->mag = mag[band];
865
          sample->magerr = magerr[band];
866
          }
867
        else
868
          {
869
          for (b=0; b<nband; b++)
870
            if (mag[b] < 98.0)
871
              {
872
              sample->mag = mag[b];
873
              sample->magerr = magerr[b];
874
              break;
875
              }
876
/*-------- If no magnitude is suitable, drop this source */
877
          if (mag[b] > 98.0)
878
            continue;
879
          }
880 308 bertin
/*
881 112 bertin
        if (nband>1 && mag[0]<98.0 && mag[1]<98.0)
882
          sample->colour = mag[0] - mag[1];
883
        else
884
          sample->colour = 99.0;
885 308 bertin
*/
886 306 bertin
/*
887
printf("%10f %10f +/-%4.0f %4.0f  %4.0f %4.0f +/-%4.0f %4.0f  %7.2f  %5.2f +/- %4.2f\n",
888
alpha,delta,poserr[lng]*DEG/MAS,poserr[lat]*DEG/MAS,
889
prop[lng]*DEG/MAS,prop[lat]*DEG/MAS,properr[lng]*DEG/MAS,properr[lat]*DEG/MAS,
890
epoch, sample->mag, sample->magerr);
891
*/
892 212 bertin
        sample->flux = 0.0;
893 314 bertin
        sample->rawpos[lng] = sample->rawpos[lat] = 0.0;
894
        sample->rawposerr[lng] = sample->rawposerr[lat] = 0.0;
895 2 bertin
        sample->wcspos[lng] = alpha;
896
        sample->wcspos[lat] = delta;
897
        sample->wcsposerr[lng] = poserr[lng];
898
        sample->wcsposerr[lat] = poserr[lat];
899 305 bertin
        sample->epoch = epoch;
900 314 bertin
        sample->spread = sample->spreaderr = 0.0;
901 251 bertin
        sample->sexflags = 0;    /* SEx flags not relevant for ref. sources*/
902 308 bertin
        sample->scampflags = 0;
903 318 bertin
        sample->imaflags = 0;
904 2 bertin
        sample->set = set;
905
        nsample++;
906
        }
907
      n++;
908
      }
909
 
910 331 bertin
  url_fclose(file);
911 2 bertin
 
912
/* Escape here if no source returned */
913
  if (!nsample)
914
    return NULL;
915
 
916
  set->nsample = nsample;
917
 
918
/* Don't waste memory! */
919
  realloc_samples(set, nsample);
920
 
921
/* Now "create" the field */
922
  QCALLOC(field, fieldstruct, 1);
923
  QMALLOC(field->set, setstruct *, 1);
924
  field->set[0] = set;
925
  field->nset = 1;
926
/* This is a reference catalog */
927
  field->astromlabel = field->photomlabel = -1;
928
  set->lng = field->lng = lng;
929
  set->lat = field->lat = lat;
930
  set->naxis = field->naxis = naxis;
931
  field->maxradius = set->radius = maxradius;
932
  set->field = field;
933
  for (d=0; d<naxis; d++)
934
    {
935
    set->wcspos[d] = field->meanwcspos[d] = wcspos[d];
936
    ctype[d] = (char *)astref_ctype[d];
937
    }
938
 
939
/* Create a dummy WCS structure */
940
  set->wcs = create_wcs(ctype, wcspos, (double *)astref_crpix,
941
                        (double *)astref_cdelt,
942
                        (int *)astref_naxisn, naxis);
943
 
944
  field->nsample = set->nsample;
945
 
946
  return field;
947
  }
948
 
949
 
950
/****** save_astreffield ******************************************************
951
PROTO   void save_astreffield(char *filename, fieldstruct *reffield)
952
PURPOSE Save an astrometric reference catalog in (pseudo-)LDAC format.
953
INPUT   Catalog name,
954
        pointer to the reference field to save.
955
OUTPUT  -.
956
NOTES   Global preferences are used.
957
AUTHOR  E. Bertin (IAP)
958 226 bertin
VERSION 22/10/2009
959 2 bertin
*/
960
void    save_astreffield(char *filename,  fieldstruct *reffield)
961
  {
962
   static char  imtabtemplate[][80] = {
963
"SIMPLE  =                    T / This is a FITS file",
964
"BITPIX  =                    8 / ",
965
"NAXIS   =                    2 / 2D data",
966
"NAXIS1  =                    1 / Number of rows",
967
"NAXIS2  =                    1 / Number of columns",
968
"EXTEND  =                    T / This file may contain FITS extensions",
969
"END                            "};
970
  catstruct     *cat;
971
  tabstruct     *asctab, *imtab, *objtab;
972
  keystruct     *key, *objkeys;
973
  setstruct     *set;
974
  samplestruct  *sample,
975
                objsample;
976
  char          *buf;
977 226 bertin
  long          dptr;
978 2 bertin
  int           i,k,n,s;
979
 
980
/* Create a new output catalog */
981
  cat = new_cat(1);
982
  init_cat(cat);
983
  strcpy(cat->filename, filename);
984
  if (open_cat(cat, WRITE_ONLY) != RETURN_OK)
985
    error(EXIT_FAILURE, "*Error*: cannot open for writing ", filename);
986
 
987
/* Primary header */
988
  save_tab(cat, cat->tab);
989
 
990
  for (s=0; s<reffield->nset; s++)
991
    {
992
    set = reffield->set[s];
993
/*-- We create a dummy table (only used through its header) */
994
    QCALLOC(asctab, tabstruct, 1);
995
    asctab->headnblock = 1 + (sizeof(imtabtemplate)-1)/FBSIZE;
996
    QCALLOC(asctab->headbuf, char, asctab->headnblock*FBSIZE);
997
    memcpy(asctab->headbuf, imtabtemplate, sizeof(imtabtemplate));
998
    for (buf = asctab->headbuf, i=FBSIZE*asctab->headnblock; i--; buf++)
999
      if (!*buf)
1000
        *buf = ' ';
1001
    write_wcs(asctab, set->wcs);
1002
/*-- (dummy) LDAC Image header */
1003
 
1004
    imtab = new_tab("LDAC_IMHEAD");
1005
    key = new_key("Field Header Card");
1006
    key->ptr = asctab->headbuf;
1007
    asctab->headbuf = NULL;
1008
    free_tab(asctab);
1009
    key->naxis = 2;
1010
    QMALLOC(key->naxisn, int, key->naxis);
1011
    key->naxisn[0] = 80;
1012
    key->naxisn[1] = 36;
1013
    key->htype = H_STRING;
1014
    key->ttype = T_STRING;
1015
    key->nobj = 1;
1016
    key->nbytes = 80*(fitsfind(key->ptr, "END     ")+1);
1017
    add_key(key, imtab, 0);
1018
    save_tab(cat, imtab);
1019
    free_tab(imtab);
1020
 
1021
/*-- LDAC Object header */
1022
    objtab = new_tab("LDAC_OBJECTS");
1023
    objtab->cat = cat;
1024
/*-- Set key pointers */
1025
    QCALLOC(objkeys, keystruct, (sizeof(refkey) / sizeof(keystruct)));
1026 226 bertin
    dptr = (long)((char *)&objsample - (char *)&refsample);
1027 2 bertin
    for (k=0; refkey[k].name[0]; k++)
1028
      {
1029
      objkeys[k] = refkey[k];
1030 226 bertin
      objkeys[k].ptr = (void *)((char *)objkeys[k].ptr + dptr); /* a trick */
1031 2 bertin
      add_key(&objkeys[k],objtab, 0);
1032
      }
1033
    init_writeobj(cat, objtab, &buf);
1034
    sample = set->sample;
1035
    for (n=set->nsample; n--;)
1036
      {
1037
      objsample = *(sample++);
1038
      write_obj(objtab, buf);
1039
      }
1040
    end_writeobj(cat, objtab, buf);
1041
    objtab->key = NULL;
1042
    objtab->nkey = 0;
1043
    free_tab(objtab);
1044
    free(objkeys);
1045
    }
1046
 
1047
  free_cat(&cat, 1);
1048
 
1049
  return;
1050
  }
1051
 
1052
/****** load_astreffield ******************************************************
1053
PROTO   fieldstruct *load_astreffield(char *filename, double *wcspos,
1054 218 bertin
                        int lng, int lat, int naxis, double maxradius, int band,
1055
                        double *maglim)
1056 2 bertin
PURPOSE Load a reference catalog in (pseudo-)LDAC format.
1057
INPUT   Catalog name,
1058 211 bertin
        pointer to the field center coordinates,
1059
        longitude index,
1060
        latitude index,
1061
        number of axes (dimensions),
1062
        search radius (in degrees),
1063 218 bertin
        band index,
1064
        magnitude range.
1065 2 bertin
OUTPUT  Pointer to the reference catalog field structure.
1066
NOTES   Global preferences are used.
1067
AUTHOR  E. Bertin (IAP)
1068 218 bertin
VERSION 02/10/2009
1069 2 bertin
*/
1070
fieldstruct     *load_astreffield(char *filename, double *wcspos,
1071
                                int lng, int lat,
1072 218 bertin
                                int naxis, double maxradius, int band,
1073
                                double *maglim)
1074 2 bertin
  {
1075
   catstruct    *cat;
1076
   tabstruct    *tab;
1077
   fieldstruct  *field;
1078
   setstruct    *set;
1079
   char         str[MAXCHAR],
1080
                *rfilename, *pstr, *pspath;
1081
   int          i, n, nsample;
1082
 
1083
 
1084
/* A short, "relative" version of the filename */
1085
  if (!(rfilename = strrchr(filename, '/')))
1086
    rfilename = filename;
1087
  else
1088
    rfilename++;
1089
 
1090 211 bertin
  sprintf(str,"Examining Catalog %s...", rfilename);
1091 2 bertin
  NFPRINTF(OUTPUT, str);
1092
 
1093
/*-- Read input catalog */
1094
  if (!(cat = read_cat(filename)))
1095
    error(EXIT_FAILURE, "*Error*: No such catalog: ", filename);
1096
  QCALLOC(field, fieldstruct, 1);
1097
  QMALLOC(field->set, setstruct *, 1);
1098
  strcpy (field->filename, filename);
1099
  field->rfilename = rfilename;
1100
 
1101
/* Create a file name with a "header" extension */
1102
  strcpy(field->hfilename, filename);
1103
  if (!(pstr = strrchr(field->hfilename, '.')))
1104
    pstr = field->hfilename+strlen(field->hfilename);
1105
  sprintf(pstr, "%s", prefs.ahead_suffix);
1106
 
1107
/* Extract the path from the filename */
1108
#ifdef HAVE_GETENV
1109
  pspath = getenv("PWD");
1110
#else
1111
  pspath = NULL;
1112
#endif
1113
  if (*field->filename == '/')
1114
    strcpy(field->path, field->filename);
1115
  else
1116
    {
1117
    strcpy(field->path, pspath? pspath: ".");
1118
    if (*field->filename != '.' && (pstr = strchr(field->filename, '/')))
1119
      {
1120
      strcat(field->path, "/");
1121
      strcat(field->path, pstr+1);
1122
      }
1123
    }
1124
  if ((pstr = strrchr(field->path, '/')))
1125
    *pstr = '\0';
1126
 
1127
/* Identify image headers in catalog  */
1128
  tab = cat->tab;
1129
 
1130
  n = 0;
1131
 
1132
/* Find the object table */
1133 211 bertin
  sprintf(str,"Loading Catalog %s...", rfilename);
1134
  NFPRINTF(OUTPUT, str);
1135 2 bertin
  tab = cat->tab;
1136
  set = NULL;
1137
  nsample = 0;
1138
  n = 0;
1139
  for (i=cat->ntab; i--; tab=tab->nexttab)
1140
    if (!strcmp("LDAC_OBJECTS", tab->extname)
1141
        || !strcmp("OBJECTS", tab->extname))
1142
    {
1143
    if (field->nset>1)
1144
      sprintf(str, "%s [%d]", rfilename, n+1);
1145
    else
1146
      strcpy(str, rfilename);
1147 218 bertin
    set = read_astrefsamples(set,tab, str, wcspos,lng,lat,naxis,maxradius,band,
1148
                        maglim);
1149 2 bertin
    nsample += set->nsample;
1150
    n++;
1151
    }
1152
 
1153
  field->nsample = nsample;
1154
  free_cat(&cat, 1);
1155
 
1156
  if (!n)
1157
    {
1158
/*-- No source found: return a NULL pointer */
1159
    end_field(field);
1160
    return NULL;
1161
    }
1162
 
1163
/* This is a reference catalog */
1164
  field->astromlabel = field->photomlabel = -1;
1165
  field->set[0] = set;
1166
  field->nset = 1;
1167
  set->lng = field->lng = lng;
1168
  set->lat = field->lat = lat;
1169
  set->naxis = field->naxis = naxis;
1170
  set->radius = maxradius;
1171
  set->field = field;
1172
 
1173
  field->nsample = set->nsample;
1174
 
1175
  return field;
1176
  }
1177
 
1178
 
1179
/****** read_astrefsamples ****************************************************
1180
PROTO   setstruct *read_astrefsamples(setstruct *set, tabstruct *tab,
1181
                                char *rfilename,
1182
                                double *wcspos, int lng, int lat, int naxis,
1183 218 bertin
                                double maxradius, int band, double *maglim)
1184 2 bertin
PURPOSE Read a set of astrometric reference samples.
1185
INPUT   Set structure pointer,
1186
        Pointer to the tab that contains the catalog,
1187
        Reduced filename.
1188
        Coordinate vector of the center,
1189
        Longitude index,
1190
        Latitude index,
1191
        Number of axes (dimensions),
1192 211 bertin
        Search radius (in degrees),
1193 218 bertin
        band index,
1194
        magnitude range.
1195 2 bertin
OUTPUT  setstruct pointer (allocated if the input setstruct pointer is NULL).
1196
NOTES   The filename is used for error messages only. Global preferences are
1197
        used.
1198 92 bertin
AUTHOR  E. Bertin (IAP)
1199 318 bertin
VERSION 12/11/2013
1200 2 bertin
*/
1201
setstruct *read_astrefsamples(setstruct *set, tabstruct *tab, char *rfilename,
1202
                                double *wcspos, int lng, int lat, int naxis,
1203 218 bertin
                                double maxradius, int band, double *maglim)
1204 2 bertin
 
1205
 
1206
  {
1207
   tabstruct            *keytab;
1208
   keystruct            *key;
1209
   samplestruct         *sample;
1210
   char                 str[MAXCHAR];
1211
   char                 *buf;
1212
   unsigned short       *flags;
1213 305 bertin
   float                *xm,*ym, *mag, *magerr, *obsdate, *erra,*errb;
1214
   double               *dxm, *dym, *dmag, *dmagerr, *dobsdate, *derra, *derrb,
1215 218 bertin
                        x,y, dx,dy,dfac, ea,eb, maxradius2, mmag;
1216
   int                  n, nsample,nsamplemax, nobj, objflags, maglimflag;
1217 2 bertin
 
1218
/* One needs 2 angular coordinates here! */
1219
  dxm = dym = dmag = derra = derrb = NULL;
1220
  xm = ym = mag = erra = errb = NULL;
1221
  maxradius2 = maxradius*maxradius;
1222
  dfac = (lng!=lat)? cos(wcspos[lat]*DEG) : 1.0;
1223
 
1224
/* If a NULL pointer is provided, we allocate a new set */
1225
  if (!set)
1226
    {
1227
    set = init_set();
1228
    nsample = nsamplemax = 0;
1229
    }
1230
  else
1231
    nsample = nsamplemax = set->nsample;
1232
 
1233
/* Init the single-row tab */
1234
  keytab = init_readobj(tab, &buf);
1235
 
1236
  if (!(key = name_to_key(keytab, prefs.astrefcent_key[0])))
1237
    {
1238
    sprintf(str, "*Error*: %s parameter not found in catalog ",
1239
                prefs.astrefcent_key[0]);
1240
    error(EXIT_FAILURE, str, rfilename);
1241
    }
1242
  if (key->ttype == T_DOUBLE)
1243
    dxm = (double *)key->ptr;
1244
  else
1245
    xm = (float *)key->ptr;
1246
 
1247
  nobj = key->nobj;
1248
 
1249
  if (!(key = name_to_key(keytab, prefs.astrefcent_key[1])))
1250
    {
1251
    sprintf(str, "*Error*: %s parameter not found in catalog ",
1252
                prefs.astrefcent_key[1]);
1253
    error(EXIT_FAILURE, str, rfilename);
1254
    }
1255
  if (key->ttype == T_DOUBLE)
1256
    dym = (double *)key->ptr;
1257
  else
1258
    ym = (float *)key->ptr;
1259
 
1260
  if (!(key = name_to_key(keytab, prefs.astreferr_key[0])))
1261
    {
1262
    sprintf(str, "*Error*: %s parameter not found in catalog ",
1263
                prefs.astreferr_key[0]);
1264
    error(EXIT_FAILURE, str, rfilename);
1265
    }
1266
  if (key->ttype == T_DOUBLE)
1267
    derra = (double *)key->ptr;
1268
  else
1269
    erra = (float *)key->ptr;
1270
 
1271
  if (!(key = name_to_key(keytab, prefs.astreferr_key[1])))
1272
    {
1273
    sprintf(str, "*Error*: %s parameter not found in catalog ",
1274
                prefs.astreferr_key[1]);
1275
    error(EXIT_FAILURE, str, rfilename);
1276
    }
1277
  if (key->ttype == T_DOUBLE)
1278
    derrb = (double *)key->ptr;
1279
  else
1280
    errb = (float *)key->ptr;
1281
 
1282
  if (!(key = name_to_key(keytab, prefs.astrefmag_key)))
1283
    {
1284
    sprintf(str, "*Error*: %s parameter not found in catalog ",
1285
                prefs.astrefmag_key);
1286
    error(EXIT_FAILURE, str, rfilename);
1287
    }
1288
  if (key->ttype == T_DOUBLE)
1289
    dmag = (double *)key->ptr;
1290
  else
1291
    mag = (float *)key->ptr;
1292
 
1293 305 bertin
  if (!(key = name_to_key(keytab, prefs.astrefmagerr_key)))
1294
    {
1295
    sprintf(str, "%s parameter not found in catalog ", prefs.astrefmagerr_key);
1296
    warning(str, rfilename);
1297
    dmagerr = NULL;
1298
    magerr = NULL;
1299
    }
1300
  else
1301
    {
1302
    if (key->ttype == T_DOUBLE)
1303
      dmagerr = (double *)key->ptr;
1304
    else
1305
      magerr = (float *)key->ptr;
1306
    }
1307
 
1308
  if (!(key = name_to_key(keytab, prefs.astrefobsdate_key)))
1309
    {
1310
    sprintf(str, "%s parameter not found in catalog ", prefs.astrefobsdate_key);
1311
    warning(str, rfilename);
1312
    dobsdate = NULL;
1313
    obsdate = NULL;
1314
    }
1315
  else
1316
    {
1317
    if (key->ttype == T_DOUBLE)
1318
      dobsdate = (double *)key->ptr;
1319
    else
1320
      obsdate = (float *)key->ptr;
1321
    }
1322
 
1323 211 bertin
/* Check that catalog contains enough bands if needed */
1324
  if (band && (!key->naxis || band>=*key->naxisn))
1325
    {
1326
    sprintf(str, "*Error*: band #%d not found in catalog ", band+1);
1327
    error(EXIT_FAILURE, str, rfilename);
1328
    }
1329
 
1330 2 bertin
  if (!(key = name_to_key(keytab, "FLAGS")))
1331
    warning("FLAGS parameter not found in catalog ", rfilename);
1332
  flags = key? (unsigned short *)key->ptr : NULL;
1333
 
1334 218 bertin
/* Check magnitude limits only if needed */
1335
  maglimflag = (maglim[0]>-99.0 || maglim[1]<99.0)? 1 : 0;
1336
 
1337 2 bertin
/* Now examine each vector of the shipment */
1338
  for (n=0; nobj--; n++)
1339
    {
1340 92 bertin
    objflags = 0;
1341 2 bertin
    read_obj(keytab,tab, buf);
1342
    if (!(n%10000))
1343
      {
1344 250 bertin
      sprintf(str,"%-.36s: %d / %d references stored",
1345
                rfilename,nsample,n);
1346 2 bertin
      NFPRINTF(OUTPUT, str);
1347
      }
1348
/*---- Apply some selection over flags, fluxes... */
1349 305 bertin
    mmag = mag? mag[band] : dmag[band];
1350 218 bertin
    if (maglimflag && (mmag<maglim[0] || mmag>maglim[1]))
1351
      continue;
1352 2 bertin
    if (flags)
1353
      {
1354 318 bertin
      if (*flags & prefs.sexflags_mask)
1355 92 bertin
        continue;
1356
/*---- Mapping from SExtractor flags is straightforward */
1357 306 bertin
      objflags = *flags;
1358
      if (objflags & OBJ_SATUR)         /* A saturated object */
1359 2 bertin
        set->nsaturated++;
1360
      }
1361 96 bertin
    x = fmod((xm? *xm : *dxm) +360.0,360.0);
1362 2 bertin
    y = ym? *ym : *dym;
1363 158 bertin
    dx = x-wcspos[lng];
1364
    if (dx>180.0)
1365
      dx -= 360.0;
1366
    else if (dx<-180.0)
1367
      dx += 360.0;
1368
    dx *= dfac;
1369 2 bertin
    dy = y-wcspos[lat];
1370
    if (dx*dx+dy*dy > maxradius2)
1371
      continue;
1372
/*-- ... and check the integrity of the sample */
1373
/*-- Allocate memory for the first shipment */
1374
    if (!set->sample)
1375
      {
1376
      nsample = 0;
1377
      nsamplemax = LSAMPLE_DEFSIZE;
1378
      malloc_samples(set, nsamplemax);
1379
      }
1380
 
1381
/*-- Increase storage space to receive new candidates if needed */
1382
    if (nsample>=nsamplemax)
1383
      {
1384
       int      nadd=(int)(1.62*nsamplemax);
1385
      nsamplemax = nadd>nsamplemax?nadd:nsamplemax+1;
1386
      realloc_samples(set, nsamplemax);
1387
      }
1388
 
1389
    sample = set->sample + nsample;
1390
    sample->set = set;
1391 251 bertin
    sample->sexflags = objflags;
1392 211 bertin
    sample->mag = mag? mag[band] : dmag[band];
1393 305 bertin
    sample->magerr = magerr? magerr[band] : (dmagerr? dmagerr[band] : 0.0);
1394
    sample->epoch = dobsdate? *dobsdate : (obsdate? *obsdate : 0.0);
1395 212 bertin
    sample->flux = 0.0;
1396 2 bertin
    sample->wcspos[lng] = x;
1397
    sample->wcspos[lat] = y;
1398
    ea = erra? *erra : *derra;
1399
    eb = errb? *errb : *derrb;
1400
    sample->wcsposerr[lng] = sample->wcsposerr[lat] = sqrt(ea*ea+eb*eb);
1401
/*-- In case of a contamination, position errors are easily doubled */
1402
    if (flags && *flags>0)
1403
      sample->wcsposerr[lng] = (sample->wcsposerr[lat] *= 2.0);
1404
    nsample++;
1405
    }
1406
 
1407
  end_readobj(keytab, tab, buf);
1408
 
1409
  set->nsample = nsample;
1410
 
1411
/* Don't waste memory! */
1412
  if (nsample)
1413
    realloc_samples(set, nsample);
1414
 
1415
  return set;
1416
  }
1417
 
1418 306 bertin
/***i** astref_strncpy ****************************************************
1419
PROTO   char *astref_strncpy(char *dest, char *src, int n)
1420
PURPOSE Copy a piece of string with max length n, and end it with '\0'.
1421
INPUT   Destination string,
1422
        input string,
1423
        max number of characters.
1424
OUTPUT  Pointer to dest.
1425
NOTES   -.
1426
AUTHOR  E. Bertin (IAP)
1427
VERSION 10/09/2012
1428
*/
1429
static char     *astref_strncpy(char *dest, char *src, int n)
1430
 
1431
  {
1432
   char *destt;
1433
   int  i;
1434
 
1435
  destt = dest;
1436
  for (i=n; i-- && *src != '\0';)
1437
    *(destt++) = *(src++);
1438
  *destt = '\0';
1439
 
1440
  return dest;
1441
  }
1442
 
1443