← Back to branch summary

~ubuntu-branches/ubuntu/karmic/kst/karmic

« back to all changes in this revision

Viewing changes to kst/kst/datasources/healpix/healpix_tools_fits.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T Chen
  • Date: 2006-06-30 19:11:30 UTC
  • mfrom: (1.2.3 upstream)
  • Revision ID: james.westby@ubuntu.com-20060630191130-acumuar75bz4puty
Tags: 1.2.1-1ubuntu1
Merge from debian unstable.

Show diffs side-by-side

added added

removed removed

Lines of Context:
kst/kst/datasources/healpix/healpix_tools_fits.cpp
 
1
/***************************************************************************
 
2
                  healpix_tools.cpp  -  tools for healpix datasource
 
3
                             -------------------
 
4
    begin                : Wed June 01 2005
 
5
    copyright            : (C) 2005 Ted Kisner
 
6
    email                : tskisner.public@gmail.com
 
7
 ***************************************************************************/
 
8
 
 
9
/***************************************************************************
 
10
 *                                                                         *
 
11
 *   This program is free software; you can redistribute it and/or modify  *
 
12
 *   it under the terms of the GNU General Public License as published by  *
 
13
 *   the Free Software Foundation; either version 2 of the License, or     *
 
14
 *   (at your option) any later version.                                   *
 
15
 *                                                                         *
 
16
 ***************************************************************************/
 
17
 
 
18
#include "healpix_tools.h"
 
19
 
 
20
/* FITS keys */
 
21
 
 
22
healpix_keys *healpix_keys_alloc()
 
23
{
 
24
  healpix_keys *keys;
 
25
  keys = (healpix_keys *) calloc(1, sizeof(healpix_keys));
 
26
  keys->nskeys = 0;
 
27
  keys->nikeys = 0;
 
28
  keys->nfkeys = 0;
 
29
  keys->skeynames = NULL;
 
30
  keys->skeyvals = NULL;
 
31
  keys->skeycoms = NULL;
 
32
  keys->ikeynames = NULL;
 
33
  keys->ikeyvals = NULL;
 
34
  keys->ikeycoms = NULL;
 
35
  keys->fkeynames = NULL;
 
36
  keys->fkeyvals = NULL;
 
37
  keys->fkeycoms = NULL;
 
38
  return keys;
 
39
}
 
40
 
 
41
int healpix_keys_free(healpix_keys * keys)
 
42
{
 
43
  size_t i;
 
44
  if (keys) {
 
45
    for (i = 0; i < (keys->nskeys); i++) {
 
46
      free(keys->skeynames[i]);
 
47
      free(keys->skeyvals[i]);
 
48
      free(keys->skeycoms[i]);
 
49
    }
 
50
    if (keys->nskeys != 0) {
 
51
      free(keys->skeynames);
 
52
      free(keys->skeyvals);
 
53
      free(keys->skeycoms);
 
54
    }
 
55
    for (i = 0; i < (keys->nikeys); i++) {
 
56
      free(keys->ikeynames[i]);
 
57
      free(keys->ikeycoms[i]);
 
58
    }
 
59
    if (keys->nikeys != 0) {
 
60
      free(keys->ikeynames);
 
61
      free(keys->ikeyvals);
 
62
      free(keys->ikeycoms);
 
63
    }
 
64
    for (i = 0; i < (keys->nfkeys); i++) {
 
65
      free(keys->fkeynames[i]);
 
66
      free(keys->fkeycoms[i]);
 
67
    }
 
68
    if (keys->nfkeys != 0) {
 
69
      free(keys->fkeynames);
 
70
      free(keys->fkeyvals);
 
71
      free(keys->fkeycoms);
 
72
    }
 
73
    free(keys);
 
74
  }
 
75
  return 0;
 
76
}
 
77
 
 
78
int healpix_keys_clear(healpix_keys * keys)
 
79
{
 
80
  size_t i;
 
81
  if (keys) {
 
82
    for (i = 0; i < (keys->nskeys); i++) {
 
83
      free(keys->skeynames[i]);
 
84
      free(keys->skeyvals[i]);
 
85
      free(keys->skeycoms[i]);
 
86
    }
 
87
    if (keys->nskeys != 0) {
 
88
      free(keys->skeynames);
 
89
      free(keys->skeyvals);
 
90
      free(keys->skeycoms);
 
91
    }
 
92
    for (i = 0; i < (keys->nikeys); i++) {
 
93
      free(keys->ikeynames[i]);
 
94
      free(keys->ikeycoms[i]);
 
95
    }
 
96
    if (keys->nikeys != 0) {
 
97
      free(keys->ikeynames);
 
98
      free(keys->ikeyvals);
 
99
      free(keys->ikeycoms);
 
100
    }
 
101
    for (i = 0; i < (keys->nfkeys); i++) {
 
102
      free(keys->fkeynames[i]);
 
103
      free(keys->fkeycoms[i]);
 
104
    }
 
105
    if (keys->nfkeys != 0) {
 
106
      free(keys->fkeynames);
 
107
      free(keys->fkeyvals);
 
108
      free(keys->fkeycoms);
 
109
    }
 
110
    keys->nskeys = 0;
 
111
    keys->nikeys = 0;
 
112
    keys->nfkeys = 0;
 
113
    keys->skeynames = NULL;
 
114
    keys->skeyvals = NULL;
 
115
    keys->skeycoms = NULL;
 
116
    keys->ikeynames = NULL;
 
117
    keys->ikeyvals = NULL;
 
118
    keys->ikeycoms = NULL;
 
119
    keys->fkeynames = NULL;
 
120
    keys->fkeyvals = NULL;
 
121
    keys->fkeycoms = NULL;
 
122
  }
 
123
  return 0;
 
124
}
 
125
 
 
126
/* add keys */
 
127
 
 
128
int healpix_keys_sadd(healpix_keys * keys, char *keyname, char *keyval, char *keycom)
 
129
{
 
130
  if (keys) {
 
131
    
 
132
    (keys->skeynames) =
 
133
        (char **)realloc(keys->skeynames,
 
134
    ((keys->nskeys) + 1) * sizeof(char *));
 
135
    (keys->skeynames[keys->nskeys]) =
 
136
        (char *)calloc(HEALPIX_STRNL, sizeof(char));
 
137
    
 
138
    (keys->skeyvals) = (char **)realloc(keys->skeyvals, ((keys->nskeys) + 1) * sizeof(char *));
 
139
    (keys->skeyvals[keys->nskeys]) = (char *)calloc(HEALPIX_STRNL, sizeof(char));
 
140
    
 
141
    (keys->skeycoms) = (char **)realloc(keys->skeycoms, ((keys->nskeys) + 1) * sizeof(char *));
 
142
    (keys->skeycoms[keys->nskeys]) = (char *)calloc(HEALPIX_STRNL, sizeof(char));
 
143
    
 
144
    strncpy(keys->skeynames[keys->nskeys], keyname, HEALPIX_STRNL);
 
145
    strncpy(keys->skeyvals[keys->nskeys], keyval, HEALPIX_STRNL);
 
146
    strncpy(keys->skeycoms[keys->nskeys], keycom, HEALPIX_STRNL);
 
147
    (keys->nskeys)++;
 
148
  }
 
149
  return 0;
 
150
}
 
151
 
 
152
int healpix_keys_iadd(healpix_keys * keys, char *keyname, int keyval, char *keycom)
 
153
{
 
154
  if (keys) {
 
155
    (keys->ikeynames) =
 
156
        (char **)realloc(keys->ikeynames,
 
157
    ((keys->nikeys) + 1) * sizeof(char *));
 
158
    (keys->ikeynames[keys->nikeys]) =
 
159
        (char *)calloc(HEALPIX_STRNL, sizeof(char));
 
160
    
 
161
    (keys->ikeyvals) =
 
162
        (int *)realloc(keys->ikeyvals, ((keys->nikeys) + 1) * sizeof(int));
 
163
    
 
164
    (keys->ikeycoms) =
 
165
        (char **)realloc(keys->ikeycoms, ((keys->nikeys) + 1) * sizeof(char *));
 
166
    (keys->ikeycoms[keys->nikeys]) = (char *)calloc(HEALPIX_STRNL, sizeof(char));
 
167
    
 
168
    strncpy(keys->ikeynames[keys->nikeys], keyname, HEALPIX_STRNL);
 
169
    keys->ikeyvals[keys->nikeys] = keyval;
 
170
    strncpy(keys->ikeycoms[keys->nikeys], keycom, HEALPIX_STRNL);
 
171
    (keys->nikeys)++;
 
172
  }
 
173
  return 0;
 
174
}
 
175
 
 
176
int healpix_keys_fadd(healpix_keys * keys, char *keyname, float keyval, char *keycom)
 
177
{
 
178
  if (keys) {
 
179
    (keys->fkeynames) =
 
180
        (char **)realloc(keys->fkeynames,
 
181
    ((keys->nfkeys) + 1) * sizeof(char *));
 
182
    (keys->fkeynames[keys->nfkeys]) =
 
183
        (char *)calloc(HEALPIX_STRNL, sizeof(char));
 
184
    
 
185
    (keys->fkeyvals) =
 
186
        (float *)realloc(keys->fkeyvals, ((keys->nfkeys) + 1) * sizeof(float));
 
187
    
 
188
    (keys->fkeycoms) =
 
189
        (char **)realloc(keys->fkeycoms, ((keys->nfkeys) + 1) * sizeof(char *));
 
190
    (keys->fkeycoms[keys->nfkeys]) = (char *)calloc(HEALPIX_STRNL, sizeof(char));
 
191
    
 
192
    strncpy(keys->fkeynames[keys->nfkeys], keyname, HEALPIX_STRNL);
 
193
    keys->fkeyvals[keys->nfkeys] = keyval;
 
194
    strncpy(keys->fkeycoms[keys->nfkeys], keycom, HEALPIX_STRNL);
 
195
    (keys->nfkeys)++;
 
196
  }
 
197
  return 0;  
 
198
}
 
199
 
 
200
/* read keys */
 
201
 
 
202
int healpix_keys_read(healpix_keys * keys, fitsfile * fp, int *ret)
 
203
{
 
204
  int nread = 0;
 
205
  int rec = 0;
 
206
  size_t nexc = 21;
 
207
  char **exclist;
 
208
  char **inclist;
 
209
  char card[HEALPIX_STRNL];
 
210
  char keyval[HEALPIX_STRNL];
 
211
  char keycom[HEALPIX_STRNL];
 
212
  char keytype;
 
213
  char keyname[HEALPIX_STRNL];
 
214
  int keylen;
 
215
 
 
216
  exclist = healpix_strarr_alloc(nexc);
 
217
  inclist = healpix_strarr_alloc(1);
 
218
 
 
219
  /*exclude required keywords */
 
220
  strcpy(exclist[0], "XTENSION");
 
221
  strcpy(exclist[1], "BITPIX");
 
222
  strcpy(exclist[2], "NAXIS*");
 
223
  strcpy(exclist[3], "PCOUNT");
 
224
  strcpy(exclist[4], "GCOUNT");
 
225
  strcpy(exclist[5], "TFIELDS");
 
226
  strcpy(exclist[6], "TTYPE*");
 
227
  strcpy(exclist[7], "TFORM*");
 
228
  strcpy(exclist[8], "TUNIT*");
 
229
  strcpy(exclist[9], "EXTNAME");
 
230
  strcpy(exclist[10], "PIXTYPE");
 
231
  strcpy(exclist[11], "ORDERING");
 
232
  strcpy(exclist[12], "NSIDE");
 
233
  strcpy(exclist[13], "COORDSYS");
 
234
  strcpy(exclist[14], "INDXSCHM");
 
235
  strcpy(exclist[15], "GRAIN");
 
236
  strcpy(exclist[16], "COMMENT");
 
237
  strcpy(exclist[17], "TBCOL*");
 
238
  strcpy(exclist[18], "SIMPLE");
 
239
  strcpy(exclist[19], "EXTEND");
 
240
  strcpy(exclist[19], "CREATOR");
 
241
  strcpy(inclist[0], "*");
 
242
 
 
243
  (*ret) = 0;
 
244
  if (fits_read_record(fp, rec, card, ret)) {
 
245
    return 0;
 
246
  }
 
247
  while (!fits_find_nextkey(fp, inclist, 1, exclist, (int)nexc, card, ret)) {
 
248
    fits_get_keyname(card, keyname, &keylen, ret);
 
249
    fits_parse_value(card, keyval, keycom, ret);
 
250
    fits_get_keytype(keyval, &keytype, ret);
 
251
    switch (keytype) {
 
252
      case 'I':
 
253
        healpix_keys_iadd(keys, keyname, atoi(keyval), keycom);
 
254
        break;
 
255
      case 'F':
 
256
        healpix_keys_fadd(keys, keyname, (float)atof(keyval), keycom);
 
257
        break;
 
258
      default:
 
259
        healpix_keys_sadd(keys, keyname, keyval, keycom);
 
260
        break;
 
261
    }
 
262
    nread++;
 
263
  }
 
264
  (*ret = 0);
 
265
 
 
266
  healpix_strarr_free(exclist, nexc);
 
267
  healpix_strarr_free(inclist, 1);
 
268
 
 
269
  return nread;
 
270
}
 
271
 
 
272
/* FITS file tests */
 
273
 
 
274
int healpix_fits_map_test(char *filename, size_t * nside, int *order, int *coord, int *type, size_t * nmaps)
 
275
{
 
276
  fitsfile *fp;
 
277
  int ret = 0;
 
278
  int ttype;
 
279
  int inside;
 
280
  long nrows;
 
281
  long pcount;
 
282
  int tfields;
 
283
  int grain;
 
284
  char pixtype[HEALPIX_STRNL];
 
285
  char coordstr[HEALPIX_STRNL];
 
286
  char orderstr[HEALPIX_STRNL];
 
287
  char indxstr[HEALPIX_STRNL];
 
288
  char comment[HEALPIX_STRNL];
 
289
  char extname[HEALPIX_STRNL];
 
290
  float nullval = HEALPIX_NULL;
 
291
  int nnull;
 
292
  float testval;
 
293
  long keynpix;
 
294
  long keyfirst;
 
295
  int ischunk = 0;
 
296
  int lastcol;
 
297
 
 
298
  /* open file */
 
299
  if (fits_open_file(&fp, filename, READONLY, &ret)) {
 
300
    return 0;
 
301
  }
 
302
 
 
303
  /* make sure extension is a binary table */
 
304
  if (fits_movabs_hdu(fp, 2, &ttype, &ret)) {
 
305
    ret = 0;
 
306
    fits_close_file(fp, &ret);
 
307
    return 0;
 
308
  }
 
309
  if (ttype != BINARY_TBL) {
 
310
    ret = 0;
 
311
    fits_close_file(fp, &ret);
 
312
    return 0;
 
313
  }
 
314
 
 
315
  /* determine the number of maps */
 
316
  if (fits_read_btblhdr(fp, HEALPIX_FITS_MAXCOL, &nrows, &tfields, NULL, NULL, NULL, extname, &pcount, &ret)) {
 
317
    ret = 0;
 
318
    fits_close_file(fp, &ret);
 
319
    return 0;
 
320
  }
 
321
 
 
322
  /* make sure this is a HEALPIX file */
 
323
  if (fits_read_key(fp, TSTRING, "PIXTYPE", pixtype, comment, &ret)) {
 
324
    ret = 0;
 
325
    fits_close_file(fp, &ret);
 
326
    return 0;
 
327
  }
 
328
  if (strncmp(pixtype, "HEALPIX", HEALPIX_STRNL) != 0) {
 
329
    ret = 0;
 
330
    fits_close_file(fp, &ret);
 
331
    return 0;
 
332
  }
 
333
 
 
334
  /* check required keywords (NSIDE, ORDERING, COORDSYS) */
 
335
  if (fits_read_key(fp, TINT, "NSIDE", &inside, comment, &ret)) {
 
336
    ret = 0;
 
337
    fits_close_file(fp, &ret);
 
338
    return 0;
 
339
  }
 
340
  (*nside) = (size_t) inside;
 
341
  if (healpix_nsidecheck(*nside)) {
 
342
    ret = 0;
 
343
    fits_close_file(fp, &ret);
 
344
    return 0;
 
345
  }
 
346
  if (fits_read_key(fp, TSTRING, "ORDERING", orderstr, comment, &ret)) {
 
347
    ret = 0;
 
348
    fits_close_file(fp, &ret);
 
349
    return 0;
 
350
  }
 
351
  if (strncmp(orderstr, "RING", HEALPIX_STRNL) == 0) {
 
352
    (*order) = HEALPIX_RING;
 
353
  } else if (strncmp(orderstr, "NESTED", HEALPIX_STRNL) == 0) {
 
354
    (*order) = HEALPIX_NEST;
 
355
  } else {
 
356
    ret = 0;
 
357
    fits_close_file(fp, &ret);
 
358
    return 0;
 
359
  }
 
360
  if (fits_read_key(fp, TSTRING, "COORDSYS", coordstr, comment, &ret)) {
 
361
    ret = 0;
 
362
    (*coord) = HEALPIX_COORD_C;
 
363
  } else {
 
364
    if (strncmp(coordstr, "C", HEALPIX_STRNL) == 0) {
 
365
      (*coord) = HEALPIX_COORD_C;
 
366
    } else if (strncmp(coordstr, "G", HEALPIX_STRNL) == 0) {
 
367
      (*coord) = HEALPIX_COORD_G;
 
368
    } else if (strncmp(coordstr, "E", HEALPIX_STRNL) == 0) {
 
369
      (*coord) = HEALPIX_COORD_E;
 
370
    } else {
 
371
      (*coord) = HEALPIX_COORD_O;
 
372
    }
 
373
  }
 
374
 
 
375
  /* check for INDXSCHM and GRAIN.  if not found, assume implicit */
 
376
  
 
377
  strcpy(indxstr,"");
 
378
  if (fits_read_key(fp, TSTRING, "OBJECT", indxstr, comment, &ret)) {
 
379
    ret = 0;
 
380
    if (fits_read_key(fp, TSTRING, "INDXSCHM", indxstr, comment, &ret)) {
 
381
      ret = 0;
 
382
      (*type) = HEALPIX_FITS_FULL;
 
383
    } else {
 
384
      if (strncmp(indxstr, "EXPLICIT", HEALPIX_STRNL) == 0) {
 
385
        (*type) = HEALPIX_FITS_CUT;
 
386
      } else {
 
387
        (*type) = HEALPIX_FITS_FULL;
 
388
      }
 
389
    }
 
390
    if (fits_read_key(fp, TINT, "GRAIN", &grain, comment, &ret)) {
 
391
      ret = 0;
 
392
      grain = 0;
 
393
    }
 
394
    if (((grain == 0) && ((*type) == HEALPIX_FITS_CUT)) ||
 
395
        ((grain != 0) && ((*type) == HEALPIX_FITS_FULL))) {
 
396
      ret = 0;
 
397
      fits_close_file(fp, &ret);
 
398
      return 0;
 
399
    }
 
400
  } else {   
 
401
    if (strncmp(indxstr, "PARTIAL", HEALPIX_STRNL) == 0) {
 
402
      (*type) = HEALPIX_FITS_CUT;
 
403
    } else {
 
404
      if (strncmp(indxstr, "FULLSKY", HEALPIX_STRNL) == 0) {
 
405
        (*type) = HEALPIX_FITS_FULL;
 
406
      } else {
 
407
        if (fits_read_key(fp, TSTRING, "INDXSCHM", indxstr, comment, &ret)) {
 
408
          ret = 0;
 
409
          (*type) = HEALPIX_FITS_FULL;
 
410
        } else {
 
411
          if (strncmp(indxstr, "EXPLICIT", HEALPIX_STRNL) == 0) {
 
412
            (*type) = HEALPIX_FITS_CUT;
 
413
          } else {
 
414
            (*type) = HEALPIX_FITS_FULL;
 
415
          }
 
416
        }
 
417
        if (fits_read_key(fp, TINT, "GRAIN", &grain, comment, &ret)) {
 
418
          ret = 0;
 
419
          grain = 0;
 
420
        }
 
421
        if (((grain == 0) && ((*type) == HEALPIX_FITS_CUT)) ||
 
422
            ((grain != 0) && ((*type) == HEALPIX_FITS_FULL))) {
 
423
          ret = 0;
 
424
          fits_close_file(fp, &ret);
 
425
          return 0;
 
426
        } 
 
427
      } 
 
428
    }
 
429
  }
 
430
  
 
431
  /* check for chunk file and truncation */
 
432
  
 
433
  if ((*type) == HEALPIX_FITS_FULL) {
 
434
    (*nmaps) = tfields;
 
435
    if ((nrows != (long)(12*inside*inside))&&(1024*nrows != (long)(12*inside*inside))) {
 
436
      /*is this a chunk file?*/
 
437
      if (fits_read_key(fp, TLONG, "FIRSTPIX", &keyfirst, comment, &ret)) {
 
438
        /*must at least have FIRSTPIX key*/
 
439
        ret = 0;
 
440
        ischunk = 0;
 
441
      } else {
 
442
        if (fits_read_key(fp, TLONG, "NPIX", &keynpix, comment, &ret)) {
 
443
          ret = 0;
 
444
          /*might be using LASTPIX instead*/
 
445
          if (fits_read_key(fp, TLONG, "LASTPIX", &keynpix, comment, &ret)) {
 
446
            ret = 0;
 
447
            ischunk = 0;
 
448
          } else {
 
449
            keynpix = keynpix - keyfirst + 1;
 
450
            if ((keyfirst < 0)||(keynpix < 0)||(keynpix+keyfirst > (long)(12*inside*inside))) {
 
451
              ischunk = 0;
 
452
            } else {
 
453
              ischunk = 1;
 
454
            }
 
455
          }
 
456
        } else {
 
457
          if ((keyfirst < 0)||(keynpix < 0)||(keynpix+keyfirst > (long)(12*inside*inside))) {
 
458
            ischunk = 0;
 
459
          } else {
 
460
            ischunk = 1;
 
461
          }
 
462
        }
 
463
      }
 
464
    } else {
 
465
      /*header doesn't matter, since we have entire map*/
 
466
      if (nrows == (long)(12*inside*inside)) {
 
467
        lastcol = 1;
 
468
      } else {
 
469
        lastcol = 1024;
 
470
      }
 
471
    }
 
472
    if (ischunk) {
 
473
      if (nrows == keynpix) {
 
474
        lastcol = 1;
 
475
      } else {
 
476
        lastcol = keynpix % 1024;
 
477
        if (lastcol == 0) {
 
478
          lastcol = 1024;
 
479
        }
 
480
      }
 
481
    }
 
482
    if (fits_read_col(fp, TFLOAT, 1, nrows, lastcol, 1, &nullval, &testval, &nnull, &ret)) {
 
483
      ret = 0;
 
484
      fits_close_file(fp, &ret);
 
485
      return 0;
 
486
    }
 
487
  } else {
 
488
    (*nmaps) = tfields - 3;
 
489
    if (fits_read_col(fp, TFLOAT, 2, nrows, 1, 1, &nullval, &testval, &nnull, &ret)) {
 
490
      ret = 0;
 
491
      fits_close_file(fp, &ret);
 
492
      return 0;
 
493
    } 
 
494
  }
 
495
  
 
496
  fits_close_file(fp, &ret);
 
497
  return 1;
 
498
}
 
499
 
 
500
int healpix_fits_map_info(char *filename, size_t * nside, int *order, int *coord, int *type, size_t * nmaps, char *creator, char *extname, char **names, char **units, healpix_keys *keys)
 
501
{
 
502
  fitsfile *fp;
 
503
  int ret = 0;
 
504
  int ttype;
 
505
  int inside;
 
506
  long nrows;
 
507
  long pcount;
 
508
  int tfields;
 
509
  int grain;
 
510
  char pixtype[HEALPIX_STRNL];
 
511
  char coordstr[HEALPIX_STRNL];
 
512
  char orderstr[HEALPIX_STRNL];
 
513
  char indxstr[HEALPIX_STRNL];
 
514
  char comment[HEALPIX_STRNL];
 
515
  float nullval = HEALPIX_NULL;
 
516
  int nnull;
 
517
  float testval;
 
518
  long keynpix;
 
519
  long keyfirst;
 
520
  int ischunk = 0;
 
521
  int lastcol;
 
522
 
 
523
  /* open file */
 
524
  if (fits_open_file(&fp, filename, READONLY, &ret)) {
 
525
    return 0;
 
526
  }
 
527
  
 
528
  /* read header info */
 
529
  fits_read_key(fp, TSTRING, "CREATOR", creator, comment, &ret);
 
530
  ret = 0;
 
531
 
 
532
  /* make sure extension is a binary table */
 
533
  if (fits_movabs_hdu(fp, 2, &ttype, &ret)) {
 
534
    ret = 0;
 
535
    fits_close_file(fp, &ret);
 
536
    return 0;
 
537
  }
 
538
  if (ttype != BINARY_TBL) {
 
539
    ret = 0;
 
540
    fits_close_file(fp, &ret);
 
541
    return 0;
 
542
  }
 
543
 
 
544
  /* determine the number of maps */
 
545
  if (fits_read_btblhdr(fp, HEALPIX_FITS_MAXCOL, &nrows, &tfields, names, NULL, units, extname, &pcount, &ret)) {
 
546
    ret = 0;
 
547
    fits_close_file(fp, &ret);
 
548
    return 0;
 
549
  }
 
550
  
 
551
  /* make sure this is a HEALPIX file */
 
552
  if (fits_read_key(fp, TSTRING, "PIXTYPE", pixtype, comment, &ret)) {
 
553
    ret = 0;
 
554
    fits_close_file(fp, &ret);
 
555
    return 0;
 
556
  }
 
557
  if (strncmp(pixtype, "HEALPIX", HEALPIX_STRNL) != 0) {
 
558
    ret = 0;
 
559
    fits_close_file(fp, &ret);
 
560
    return 0;
 
561
  }
 
562
 
 
563
  /* check required keywords (NSIDE, ORDERING, COORDSYS) */
 
564
  if (fits_read_key(fp, TINT, "NSIDE", &inside, comment, &ret)) {
 
565
    ret = 0;
 
566
    fits_close_file(fp, &ret);
 
567
    return 0;
 
568
  }
 
569
  (*nside) = (size_t) inside;
 
570
  if (healpix_nsidecheck(*nside)) {
 
571
    ret = 0;
 
572
    fits_close_file(fp, &ret);
 
573
    return 0;
 
574
  }
 
575
  if (fits_read_key(fp, TSTRING, "ORDERING", orderstr, comment, &ret)) {
 
576
    ret = 0;
 
577
    fits_close_file(fp, &ret);
 
578
    return 0;
 
579
  }
 
580
  if (strncmp(orderstr, "RING", HEALPIX_STRNL) == 0) {
 
581
    (*order) = HEALPIX_RING;
 
582
  } else if (strncmp(orderstr, "NESTED", HEALPIX_STRNL) == 0) {
 
583
    (*order) = HEALPIX_NEST;
 
584
  } else {
 
585
    ret = 0;
 
586
    fits_close_file(fp, &ret);
 
587
    return 0;
 
588
  }
 
589
  if (fits_read_key(fp, TSTRING, "COORDSYS", coordstr, comment, &ret)) {
 
590
    ret = 0;
 
591
    (*coord) = HEALPIX_COORD_C;
 
592
  } else {
 
593
    if (strncmp(coordstr, "C", HEALPIX_STRNL) == 0) {
 
594
      (*coord) = HEALPIX_COORD_C;
 
595
    } else if (strncmp(coordstr, "G", HEALPIX_STRNL) == 0) {
 
596
      (*coord) = HEALPIX_COORD_G;
 
597
    } else if (strncmp(coordstr, "E", HEALPIX_STRNL) == 0) {
 
598
      (*coord) = HEALPIX_COORD_E;
 
599
    } else {
 
600
      (*coord) = HEALPIX_COORD_O;
 
601
    }
 
602
  }
 
603
 
 
604
  /* check for INDXSCHM and GRAIN.  if not found, assume implicit */
 
605
  
 
606
  strcpy(indxstr,"");
 
607
  if (fits_read_key(fp, TSTRING, "OBJECT", indxstr, comment, &ret)) {
 
608
    ret = 0;
 
609
    if (fits_read_key(fp, TSTRING, "INDXSCHM", indxstr, comment, &ret)) {
 
610
      ret = 0;
 
611
      (*type) = HEALPIX_FITS_FULL;
 
612
    } else {
 
613
      if (strncmp(indxstr, "EXPLICIT", HEALPIX_STRNL) == 0) {
 
614
        (*type) = HEALPIX_FITS_CUT;
 
615
      } else {
 
616
        (*type) = HEALPIX_FITS_FULL;
 
617
      }
 
618
    }
 
619
    if (fits_read_key(fp, TINT, "GRAIN", &grain, comment, &ret)) {
 
620
      ret = 0;
 
621
      grain = 0;
 
622
    }
 
623
    if (((grain == 0) && ((*type) == HEALPIX_FITS_CUT)) ||
 
624
        ((grain != 0) && ((*type) == HEALPIX_FITS_FULL))) {
 
625
      ret = 0;
 
626
      fits_close_file(fp, &ret);
 
627
      return 0;
 
628
    }
 
629
  } else {   
 
630
    if (strncmp(indxstr, "PARTIAL", HEALPIX_STRNL) == 0) {
 
631
      (*type) = HEALPIX_FITS_CUT;
 
632
    } else {
 
633
      if (strncmp(indxstr, "FULLSKY", HEALPIX_STRNL) == 0) {
 
634
        (*type) = HEALPIX_FITS_FULL;
 
635
      } else {
 
636
        if (fits_read_key(fp, TSTRING, "INDXSCHM", indxstr, comment, &ret)) {
 
637
          ret = 0;
 
638
          (*type) = HEALPIX_FITS_FULL;
 
639
        } else {
 
640
          if (strncmp(indxstr, "EXPLICIT", HEALPIX_STRNL) == 0) {
 
641
            (*type) = HEALPIX_FITS_CUT;
 
642
          } else {
 
643
            (*type) = HEALPIX_FITS_FULL;
 
644
          }
 
645
        }
 
646
        if (fits_read_key(fp, TINT, "GRAIN", &grain, comment, &ret)) {
 
647
          ret = 0;
 
648
          grain = 0;
 
649
        }
 
650
        if (((grain == 0) && ((*type) == HEALPIX_FITS_CUT)) ||
 
651
            ((grain != 0) && ((*type) == HEALPIX_FITS_FULL))) {
 
652
          ret = 0;
 
653
          fits_close_file(fp, &ret);
 
654
          return 0;
 
655
        } 
 
656
      } 
 
657
    }
 
658
  }
 
659
      
 
660
  /* read extension keys */
 
661
  healpix_keys_read(keys, fp, &ret);
 
662
        
 
663
  /* check for chunk file and truncation */
 
664
  
 
665
  if ((*type) == HEALPIX_FITS_FULL) {
 
666
    (*nmaps) = tfields;
 
667
    if ((nrows != (long)(12*inside*inside))&&(1024*nrows != (long)(12*inside*inside))) {
 
668
      /*is this a chunk file?*/
 
669
      if (fits_read_key(fp, TLONG, "FIRSTPIX", &keyfirst, comment, &ret)) {
 
670
        /*must at least have FIRSTPIX key*/
 
671
        ret = 0;
 
672
        ischunk = 0;
 
673
      } else {
 
674
        if (fits_read_key(fp, TLONG, "NPIX", &keynpix, comment, &ret)) {
 
675
          ret = 0;
 
676
          /*might be using LASTPIX instead*/
 
677
          if (fits_read_key(fp, TLONG, "LASTPIX", &keynpix, comment, &ret)) {
 
678
            ret = 0;
 
679
            ischunk = 0;
 
680
          } else {
 
681
            keynpix = keynpix - keyfirst + 1;
 
682
            if ((keyfirst < 0)||(keynpix < 0)||(keynpix+keyfirst > (long)(12*inside*inside))) {
 
683
              ischunk = 0;
 
684
            } else {
 
685
              ischunk = 1;
 
686
            }
 
687
          }
 
688
        } else {
 
689
          if ((keyfirst < 0)||(keynpix < 0)||(keynpix+keyfirst > (long)(12*inside*inside))) {
 
690
            ischunk = 0;
 
691
          } else {
 
692
            ischunk = 1;
 
693
          }
 
694
        }
 
695
      }
 
696
    } else {
 
697
      /*header doesn't matter, since we have entire map*/
 
698
      if (nrows == (long)(12*inside*inside)) {
 
699
        lastcol = 1;
 
700
      } else {
 
701
        lastcol = 1024;
 
702
      }
 
703
    }
 
704
    if (ischunk) {
 
705
      if (nrows == keynpix) {
 
706
        lastcol = 1;
 
707
      } else {
 
708
        lastcol = keynpix % 1024;
 
709
        if (lastcol == 0) {
 
710
          lastcol = 1024;
 
711
        }
 
712
      }
 
713
    }
 
714
    if (fits_read_col(fp, TFLOAT, 1, nrows, lastcol, 1, &nullval, &testval, &nnull, &ret)) {
 
715
      ret = 0;
 
716
      fits_close_file(fp, &ret);
 
717
      return 0;
 
718
    }
 
719
  } else {
 
720
    (*nmaps) = tfields - 3;
 
721
    if (fits_read_col(fp, TFLOAT, 2, nrows, 1, 1, &nullval, &testval, &nnull, &ret)) {
 
722
      ret = 0;
 
723
      fits_close_file(fp, &ret);
 
724
      return 0;
 
725
    } 
 
726
  }
 
727
  
 
728
  fits_close_file(fp, &ret);
 
729
  return 1;
 
730
}
 
731