home *** CD-ROM | disk | FTP | other *** search
/ PC Pro 2002 April / pcpro0402.iso / essentials / graphics / Gimp / gimp-src-20001226.exe / src / gimp / plug-ins / fits / fitsrw.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-05-02  |  79.9 KB  |  2,097 lines
  1. /******************************************************************************/
  2. /*                      Peter Kirchgessner                                    */
  3. /*                      e-mail: pkirchg@aol.com                               */
  4. /*                      WWW   : http://members.aol.com/pkirchg                */
  5. /******************************************************************************/
  6. /*  #BEG-HDR                                                                  */
  7. /*                                                                            */
  8. /*  Package       : FITS reading/writing library                              */
  9. /*  Modul-Name    : fitsrw.c                                                  */
  10. /*  Description   : Support of reading/writing FITS-files                     */
  11. /*  Function(s)   : fits_new_filestruct    - (local) initialize file structure*/
  12. /*                  fits_new_hdulist       - (local) initialize hdulist struct*/
  13. /*                  fits_delete_filestruct - (local) delete file structure    */
  14. /*                  fits_delete_recordlist - (local) delete record list       */
  15. /*                  fits_delete_hdulist    - (local) delete hdu list          */
  16. /*                  fits_nan_32            - (local) check IEEE NaN values    */
  17. /*                  fits_nan_64            - (local) check IEEE NaN values    */
  18. /*                  fits_get_error         - get error message                */
  19. /*                  fits_set_error         - (local) set error message        */
  20. /*                  fits_drop_error        - (local) remove an error message  */
  21. /*                  fits_open              - open a FITS file                 */
  22. /*                  fits_close             - close a FITS file                */
  23. /*                  fits_add_hdu           - add a HDU to a FITS file         */
  24. /*                  fits_add_card          - add a card to the HDU            */
  25. /*                  fits_print_header      - print a single FITS header       */
  26. /*                  fits_read_header       - (local) read in FITS header      */
  27. /*                  fits_write_header      - write a FITS header              */
  28. /*                  fits_decode_header     - (local) decode a header          */
  29. /*                  fits_eval_pixrange     - (local) evaluate range of pixels */
  30. /*                  fits_decode_card       - decode a card                    */
  31. /*                  fits_search_card       - search a card in a record list   */
  32. /*                  fits_image_info        - get information about image      */
  33. /*                  fits_seek_image        - position to an image             */
  34. /*                  fits_read_pixel        - read pixel values from file      */
  35. /*                  fits_to_pgmraw         - convert FITS-file to PGM-file    */
  36. /*                  pgmraw_to_fits         - convert PGM-file to FITS-file    */
  37. /*                                                                            */
  38. /*  Author        : P. Kirchgessner                                           */
  39. /*  Date of Gen.  : 12-Apr-97                                                 */
  40. /*  Last modified : 20-Dec-97                                                 */
  41. /*  Version       : 0.11                                                      */
  42. /*  Compiler Opt. :                                                           */
  43. /*  Changes       :                                                           */
  44. /*  #MOD-0001, nn, 20-Dec-97, Initialize some variables                       */
  45. /*                                                                            */
  46. /*  #END-HDR                                                                  */
  47. /******************************************************************************/
  48. /*  References:                                                               */
  49. /*  - NOST, Definition of the Flexible Image Transport System (FITS),         */
  50. /*    September 29, 1995, Standard, NOST 100-1.1                              */
  51. /*    (ftp://nssdc.gsfc.nasa.gov/pub/fits/fits_standard_ps.Z)                 */
  52. /*  - The FITS IMAGE Extension. A Proposal. J.D. Ponz, R.W. Thompson,         */
  53. /*    J.R. Munoz, Feb. 7, 1992                                                */
  54. /*    (ftp://www.cv.nrao.edu/fits/documents/standards/image.ps.gz)            */
  55. /*                                                                            */
  56. /******************************************************************************/
  57.  
  58. #define                                    VERSIO  0.11
  59. /* Identifikation:    "@(#) <product>              <ver> <dd-mmm-yy>" */
  60. static char ident[] = "@(#) libfits.c              0.11  20-Dec-97  (%I%)";
  61.  
  62. /******************************************************************************/
  63. /* FITS reading/writing library                                               */
  64. /* Copyright (C) 1997 Peter Kirchgessner                                      */
  65. /* (email: pkirchg@aol.com, WWW: http://members.aol.com/pkirchg)              */
  66. /* The library was developed for a FITS-plug-in to GIMP, the GNU Image        */
  67. /* Manipulation Program. But it is completely independant to that. If someone */
  68. /* finds it useful for other purposes, try to keep it independant from your   */
  69. /* application.                                                               */
  70. /******************************************************************************/
  71. /*                                                                            */
  72. /* This program is free software; you can redistribute it and/or modify       */
  73. /* it under the terms of the GNU General Public License as published by       */
  74. /* the Free Software Foundation; either version 2 of the License, or          */
  75. /* (at your option) any later version.                                        */
  76. /*                                                                            */
  77. /* This program is distributed in the hope that it will be useful,            */
  78. /* but WITHOUT ANY WARRANTY; without even the implied warranty of             */
  79. /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              */
  80. /* GNU General Public License for more details.                               */
  81. /*                                                                            */
  82. /* You should have received a copy of the GNU General Public License          */
  83. /* along with this program; if not, write to the Free Software                */
  84. /* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
  85. /******************************************************************************/
  86.  
  87. #include "config.h"
  88.  
  89. #include <stdio.h>
  90. #include <stdlib.h>
  91. #include <string.h>
  92. #ifdef HAVE_UNISTD_H
  93. #include <unistd.h>
  94. #endif
  95.  
  96. #include "fitsrw.h"
  97.  
  98.  
  99. /* Declaration of local funtions */
  100.  
  101. static FITS_FILE *fits_new_filestruct (void);
  102. static FITS_HDU_LIST *fits_new_hdulist (void);
  103. static void fits_delete_filestruct (FITS_FILE *ff);
  104. static void fits_delete_recordlist (FITS_RECORD_LIST *rl);
  105. static void fits_delete_hdulist (FITS_HDU_LIST *hl);
  106. static int  fits_nan_32 (unsigned char *value);
  107. static int  fits_nan_64 (unsigned char *value);
  108. static void fits_set_error (char *errmsg);
  109. static void fits_drop_error (void);
  110. static FITS_RECORD_LIST *fits_read_header (FILE *fp, int *nrec);
  111. static FITS_HDU_LIST *fits_decode_header (FITS_RECORD_LIST *hdr,
  112.                         long hdr_offset, long dat_offset);
  113. static int  fits_eval_pixrange (FILE *fp, FITS_HDU_LIST *hdu);
  114.  
  115.  
  116. /* Error handling like a FIFO */
  117. #define FITS_MAX_ERROR      16
  118. #define FITS_ERROR_LENGTH  256
  119. static int fits_n_error = 0;
  120. static char fits_error[FITS_MAX_ERROR][FITS_ERROR_LENGTH];
  121.  
  122. /* What byte ordering for IEEE-format we are running on ? */
  123. static int fits_ieee32_intel = 0;
  124. static int fits_ieee32_motorola = 0;
  125. static int fits_ieee64_intel = 0;
  126. static int fits_ieee64_motorola = 0;
  127.  
  128. /* Macros */
  129. #define FITS_RETURN(msg, val) { fits_set_error (msg); return (val); }
  130. #define FITS_VRETURN(msg) { fits_set_error (msg); return; }
  131.  
  132. /* Get pixel values from memory. p must be an (unsigned char *) */
  133. #define FITS_GETBITPIX16(p,val) val = ((p[0] << 8) | (p[1]))
  134. #define FITS_GETBITPIX32(p,val) val = \
  135.           ((p[0] << 24) | (p[1] << 16) | (p[2] << 8) | p[3])
  136.  
  137. /* Get floating point values from memory. p must be an (unsigned char *). */
  138. /* The floating point values must directly correspond */
  139. /* to machine representation. Otherwise it does not work. */
  140. #define FITS_GETBITPIXM32(p,val) \
  141.  { if (fits_ieee32_intel) {unsigned char uc[4]; \
  142.    uc[0] = p[3]; uc[1] = p[2]; uc[2] = p[1]; uc[3] = p[0]; \
  143.    val = *(FITS_BITPIXM32 *)uc; } \
  144.    else if (fits_ieee32_motorola) { val = *(FITS_BITPIXM32 *)p; } \
  145.    else if (fits_ieee64_motorola) {FITS_BITPIXM64 m64; \
  146.    unsigned char *uc= (unsigned char *)&m64; \
  147.    uc[0]=p[0]; uc[1]=p[1]; uc[2]=p[2]; uc[3]=p[3]; uc[4]=uc[5]=uc[6]=uc[7]=0; \
  148.    val = (FITS_BITPIXM32)m64; } \
  149.    else if (fits_ieee64_intel) {FITS_BITPIXM64 i64; \
  150.    unsigned char *uc= (unsigned char *)&i64; \
  151.    uc[0]=uc[1]=uc[2]=uc[3]=0; uc[7]=p[0]; uc[6]=p[1]; uc[5]=p[2]; uc[4]=p[3]; \
  152.    val = (FITS_BITPIXM32)i64;}\
  153. }
  154.  
  155. #define FITS_GETBITPIXM64(p,val) \
  156.  { if (fits_ieee64_intel) {unsigned char uc[8]; \
  157.    uc[0] = p[7]; uc[1] = p[6]; uc[2] = p[5]; uc[3] = p[4]; \
  158.    uc[4] = p[3]; uc[5] = p[2]; uc[6] = p[1]; uc[7] = p[0]; \
  159.    val = *(FITS_BITPIXM64 *)uc; } else val = *(FITS_BITPIXM64 *)p; }
  160.  
  161. #define FITS_WRITE_BOOLCARD(fp,key,value) \
  162. {char card[81]; \
  163.  sprintf (card, "%-8.8s= %20s%50s", key, value ? "T" : "F", " "); \
  164.  fwrite (card, 1, 80, fp); }
  165.  
  166. #define FITS_WRITE_LONGCARD(fp,key,value) \
  167. {char card[81]; \
  168.  sprintf (card, "%-8.8s= %20ld%50s", key, (long)value, " "); \
  169.  fwrite (card, 1, 80, fp); }
  170.  
  171. #define FITS_WRITE_DOUBLECARD(fp,key,value) \
  172. {char card[81], dbl[21], *istr; \
  173.  sprintf (dbl, "%20f", (double)value); istr = strstr (dbl, "e"); \
  174.  if (istr) *istr = 'E'; \
  175.  sprintf (card, "%-8.8s= %20.20s%50s", key, dbl, " "); \
  176.  fwrite (card, 1, 80, fp); }
  177.  
  178. #define FITS_WRITE_STRINGCARD(fp,key,value) \
  179. {char card[81]; int k;\
  180.  sprintf (card, "%-8.8s= \'%s", key, value); \
  181.  for (k = strlen (card); k < 81; k++) card[k] = ' '; \
  182.  k = strlen (key); if (k < 8) card[19] = '\''; else card[11+k] = '\''; \
  183.  fwrite (card, 1, 80, fp); }
  184.  
  185. #define FITS_WRITE_CARD(fp,value) \
  186. {char card[81]; \
  187.  sprintf (card, "%-80.80s", value); \
  188.  fwrite (card, 1, 80, fp); }
  189.  
  190.  
  191. /*****************************************************************************/
  192. /* #BEG-PAR                                                                  */
  193. /*                                                                           */
  194. /* Function  : fits_new_filestruct - (local) initialize file structure       */
  195. /*                                                                           */
  196. /* Parameters:                                                               */
  197. /* -none-                                                                    */
  198. /*                                                                           */
  199. /* Returns a pointer to an initialized fits file structure.                  */
  200. /* On failure, a NULL-pointer is returned.                                   */
  201. /*                                                                           */
  202. /* #END-PAR                                                                  */
  203. /*****************************************************************************/
  204.  
  205. static FITS_FILE *fits_new_filestruct (void)
  206.  
  207. {FITS_FILE *ff;
  208.  
  209.  ff = (FITS_FILE *)malloc (sizeof (FITS_FILE));
  210.  if (ff == NULL) return (NULL);
  211.  
  212.  memset ((char *)ff, 0, sizeof (*ff));
  213.  return (ff);
  214. }
  215.  
  216.  
  217. /*****************************************************************************/
  218. /* #BEG-PAR                                                                  */
  219. /*                                                                           */
  220. /* Function  : fits_new_hdulist    - (local) initialize hdulist structure    */
  221. /*                                                                           */
  222. /* Parameters:                                                               */
  223. /* -none-                                                                    */
  224. /*                                                                           */
  225. /* Returns a pointer to an initialized hdulist structure.                    */
  226. /* On failure, a NULL-pointer is returned.                                   */
  227. /*                                                                           */
  228. /* #END-PAR                                                                  */
  229. /*****************************************************************************/
  230.  
  231. static FITS_HDU_LIST *fits_new_hdulist (void)
  232.  
  233. {FITS_HDU_LIST *hdl;
  234.  
  235.  hdl = (FITS_HDU_LIST *)malloc (sizeof (FITS_HDU_LIST));
  236.  if (hdl == NULL) return (NULL);
  237.  
  238.  memset ((char *)hdl, 0, sizeof (*hdl));
  239.  hdl->pixmin = hdl->pixmax = hdl->datamin = hdl->datamax = 0.0;
  240.  hdl->bzero = hdl->bscale = 0.0;
  241.  
  242.  return (hdl);
  243. }
  244.  
  245.  
  246. /*****************************************************************************/
  247. /* #BEG-PAR                                                                  */
  248. /*                                                                           */
  249. /* Function  : fits_delete_filestruct - (local) delete file structure        */
  250. /*                                                                           */
  251. /* Parameters:                                                               */
  252. /* FITS_FILE *ff   [I] : pointer to fits file structure                      */
  253. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  254. /*                                                                           */
  255. /* Frees all memory allocated by the file structure.                         */
  256. /*                                                                           */
  257. /* #END-PAR                                                                  */
  258. /*****************************************************************************/
  259.  
  260. static void fits_delete_filestruct (FITS_FILE *ff)
  261.  
  262. {
  263.  if (ff == NULL) return;
  264.  
  265.  fits_delete_hdulist (ff->hdu_list);
  266.  ff->hdu_list = NULL;
  267.  
  268.  ff->fp = NULL;
  269.  free ((char *)ff);
  270. }
  271.  
  272.  
  273. /*****************************************************************************/
  274. /* #BEG-PAR                                                                  */
  275. /*                                                                           */
  276. /* Function  : fits_delete_recordlist - (local) delete record list           */
  277. /*                                                                           */
  278. /* Parameters:                                                               */
  279. /* FITS_RECORD_LIST *rl  [I] : record list to delete                         */
  280. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  281. /*                                                                           */
  282. /* #END-PAR                                                                  */
  283. /*****************************************************************************/
  284.  
  285. static void fits_delete_recordlist (FITS_RECORD_LIST *rl)
  286.  
  287. {FITS_RECORD_LIST *next;
  288.  
  289.  while (rl != NULL)
  290.  {
  291.    next = rl->next_record;
  292.    rl->next_record = NULL;
  293.    free ((char *)rl);
  294.    rl = next;
  295.  }
  296. }
  297.  
  298.  
  299. /*****************************************************************************/
  300. /* #BEG-PAR                                                                  */
  301. /*                                                                           */
  302. /* Function  : fits_delete_hdulist - (local) delete hdu list                 */
  303. /*                                                                           */
  304. /* Parameters:                                                               */
  305. /* FITS_HDU_LIST *hl  [I] : hdu list to delete                               */
  306. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  307. /*                                                                           */
  308. /* #END-PAR                                                                  */
  309. /*****************************************************************************/
  310.  
  311. static void fits_delete_hdulist (FITS_HDU_LIST *hl)
  312.  
  313. {FITS_HDU_LIST *next;
  314.  
  315.  while (hl != NULL)
  316.  {
  317.    fits_delete_recordlist (hl->header_record_list);
  318.    next = hl->next_hdu;
  319.    hl->next_hdu = NULL;
  320.    free ((char *)hl);
  321.    hl = next;
  322.  }
  323. }
  324.  
  325.  
  326. /*****************************************************************************/
  327. /* #BEG-PAR                                                                  */
  328. /*                                                                           */
  329. /* Function  : fits_nan_32 - (local) check for IEEE NaN values (32 bit)      */
  330. /*                                                                           */
  331. /* Parameters:                                                               */
  332. /* unsigned char *v     [I] : value to check                                 */
  333. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  334. /*                                                                           */
  335. /* The function returns 1 if the value is a NaN. Otherwise 0 is returned.    */
  336. /* The byte sequence at v must start with the sign/eponent byte.             */
  337. /*                                                                           */
  338. /* #END-PAR                                                                  */
  339. /*****************************************************************************/
  340.  
  341. static int fits_nan_32 (unsigned char *v)
  342.  
  343. {register unsigned long k;
  344.  
  345.  k = (v[0] << 24) | (v[1] << 16) | (v[2] << 8) | v[3];
  346.  k &= 0x7fffffff;  /* Dont care about the sign bit */
  347.  
  348.  /* See NOST Definition of the Flexible Image Transport System (FITS), */
  349.  /* Appendix F, IEEE special formats. */
  350.  return (   ((k >= 0x7f7fffff) && (k <= 0x7fffffff))
  351.          || ((k >= 0x00000001) && (k <= 0x00800000)));
  352. }
  353.      
  354.  
  355. /*****************************************************************************/
  356. /* #BEG-PAR                                                                  */
  357. /*                                                                           */
  358. /* Function  : fits_nan_64 - (local) check for IEEE NaN values (64 bit)      */
  359. /*                                                                           */
  360. /* Parameters:                                                               */
  361. /* unsigned char *v     [I] : value to check                                 */
  362. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  363. /*                                                                           */
  364. /* The function returns 1 if the value is a NaN. Otherwise 0 is returned.    */
  365. /* The byte sequence at v must start with the sign/eponent byte.             */
  366. /* (currently we ignore the low order 4 bytes of the mantissa. Therefore     */
  367. /* this function is the same as for 32 bits).                                */
  368. /*                                                                           */
  369. /* #END-PAR                                                                  */
  370. /*****************************************************************************/
  371.  
  372. static int fits_nan_64 (unsigned char *v)
  373.  
  374. {register unsigned long k;
  375.  
  376.  k = (v[0] << 24) | (v[1] << 16) | (v[2] << 8) | v[3];
  377.  k &= 0x7fffffff;  /* Dont care about the sign bit */
  378.  
  379.  /* See NOST Definition of the Flexible Image Transport System (FITS), */
  380.  /* Appendix F, IEEE special formats. */
  381.  return (   ((k >= 0x7f7fffff) && (k <= 0x7fffffff))
  382.          || ((k >= 0x00000001) && (k <= 0x00800000)));
  383. }
  384.  
  385.  
  386. /*****************************************************************************/
  387. /* #BEG-PAR                                                                  */
  388. /*                                                                           */
  389. /* Function  : fits_get_error - get an error message                         */
  390. /*                                                                           */
  391. /* Parameters:                                                               */
  392. /* -none-                                                                    */
  393. /*                                                                           */
  394. /* If an error message has been set, a pointer to the message is returned.   */
  395. /* Otherwise a NULL pointer is returned.                                     */
  396. /* An inquired error message is removed from the error FIFO.                 */
  397. /*                                                                           */
  398. /* #END-PAR                                                                  */
  399. /*****************************************************************************/
  400.  
  401. char *fits_get_error (void)
  402.  
  403. {static char errmsg[FITS_ERROR_LENGTH];
  404.  int k;
  405.  
  406.  if (fits_n_error <= 0) return (NULL);
  407.  strcpy (errmsg, fits_error[0]);
  408.  
  409.  for (k = 1; k < fits_n_error; k++)
  410.    strcpy (fits_error[k-1], fits_error[k]);
  411.  
  412.  fits_n_error--;
  413.  
  414.  return (errmsg);
  415. }
  416.  
  417.  
  418. /*****************************************************************************/
  419. /* #BEG-PAR                                                                  */
  420. /*                                                                           */
  421. /* Function  : fits_set_error - (local) set an error message                 */
  422. /*                                                                           */
  423. /* Parameters:                                                               */
  424. /* char *errmsg   [I] : Error message to set                                 */
  425. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  426. /*                                                                           */
  427. /* Places the error message in the FIFO. If the FIFO is full,                */
  428. /* the message is discarded.                                                 */
  429. /*                                                                           */
  430. /* #END-PAR                                                                  */
  431. /*****************************************************************************/
  432.  
  433. static void fits_set_error (char *errmsg)
  434.  
  435. {
  436.   if (fits_n_error < FITS_MAX_ERROR)
  437.   {
  438.     strncpy (fits_error[fits_n_error], errmsg, FITS_ERROR_LENGTH);
  439.     fits_error[fits_n_error++][FITS_ERROR_LENGTH-1] = '\0';
  440.   }
  441. }
  442.  
  443.  
  444. /*****************************************************************************/
  445. /* #BEG-PAR                                                                  */
  446. /*                                                                           */
  447. /* Function  : fits_drop_error - (local) remove an error message             */
  448. /*                                                                           */
  449. /* Parameters:                                                               */
  450. /* -none-                                                                    */
  451. /*                                                                           */
  452. /* Removes the last error message from the error message FIFO                */
  453. /*                                                                           */
  454. /* #END-PAR                                                                  */
  455. /*****************************************************************************/
  456.  
  457. static void fits_drop_error (void)
  458.  
  459. {
  460.  if (fits_n_error > 0) fits_n_error--;
  461. }
  462.  
  463.  
  464. /*****************************************************************************/
  465. /* #BEG-PAR                                                                  */
  466. /*                                                                           */
  467. /* Function  : fits_open - open a FITS file                                  */
  468. /*                                                                           */
  469. /* Parameters:                                                               */
  470. /* char *filename   [I] : name of file to open                               */
  471. /* char *openmode   [I] : mode to open the file ("r", "w")                   */
  472. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  473. /*                                                                           */
  474. /* On success, a FITS_FILE-pointer is returned. On failure, a NULL-          */
  475. /* pointer is returned.                                                      */
  476. /* The functions scans through the file loading each header and analyzing    */
  477. /* them.                                                                     */
  478. /*                                                                           */
  479. /* #END-PAR                                                                  */
  480. /*****************************************************************************/
  481.  
  482. FITS_FILE *fits_open (char *filename, char *openmode)
  483.  
  484. {int reading, writing, n_rec, n_hdr;
  485.  long fpos_header, fpos_data;
  486.  FILE *fp;
  487.  FITS_FILE *ff;
  488.  FITS_RECORD_LIST *hdrlist;
  489.  FITS_HDU_LIST *hdulist, *last_hdulist;
  490.  
  491.  /* initialize */
  492.  
  493.  hdulist = NULL;
  494.  last_hdulist = NULL;
  495.  
  496.  /* Check the IEEE-format we are running on */
  497.  {float one32 = 1.0;
  498.   double one64 = 1.0;
  499.   unsigned char *op32 = (unsigned char *)&one32;
  500.   unsigned char *op64 = (unsigned char *)&one64;
  501.  
  502.   if (sizeof (float) == 4)
  503.   {
  504.     fits_ieee32_intel = (op32[3] == 0x3f);
  505.     fits_ieee32_motorola = (op32[0] == 0x3f);
  506.   }
  507.   if (sizeof (double) == 8)
  508.   {
  509.     fits_ieee64_intel = (op64[7] == 0x3f);
  510.     fits_ieee64_motorola = (op64[0] == 0x3f);
  511.   }
  512.  }
  513.  
  514.  if ((filename == NULL) || (*filename == '\0') || (openmode == NULL))
  515.    FITS_RETURN ("fits_open: Invalid parameters", NULL);
  516.  
  517.  reading = (strcmp (openmode, "r") == 0);
  518.  writing = (strcmp (openmode, "w") == 0);
  519.  if ((!reading) && (!writing))
  520.    FITS_RETURN ("fits_open: Invalid openmode", NULL);
  521.  
  522.  fp = fopen (filename, reading ? "rb" : "wb");
  523.  if (fp == NULL) FITS_RETURN ("fits_open: fopen() failed", NULL);
  524.  
  525.  ff = fits_new_filestruct ();
  526.  if (ff == NULL)
  527.  {
  528.    fclose (fp);
  529.    FITS_RETURN ("fits_open: No more memory", NULL);
  530.  }
  531.  
  532.  ff->fp = fp;
  533.  ff->openmode = *openmode;
  534.  
  535.  if (writing) return (ff);
  536.  
  537.  for (n_hdr = 0; ; n_hdr++)   /* Read through all HDUs */
  538.  {
  539.    fpos_header = ftell (fp);    /* Save file position of header */
  540.    hdrlist = fits_read_header (fp, &n_rec);
  541.  
  542.    if (hdrlist == NULL)
  543.    {
  544.      if (n_hdr > 0)        /* At least one header must be present. */
  545.        fits_drop_error (); /* If we got a header already, drop the error */
  546.      break;
  547.    }
  548.    fpos_data = ftell (fp);      /* Save file position of data */
  549.  
  550.                            /* Decode the header */
  551.    hdulist = fits_decode_header (hdrlist, fpos_header, fpos_data);
  552.    if (hdulist == NULL)
  553.    {
  554.      fits_delete_recordlist (hdrlist);
  555.      break;
  556.    }
  557.    ff->n_hdu++;
  558.    ff->n_pic += hdulist->numpic;
  559.  
  560.    if (hdulist->used.blank_value) ff->blank_used = 1;
  561.    if (hdulist->used.nan_value) ff->nan_used = 1;
  562.  
  563.    if (n_hdr == 0)   
  564.      ff->hdu_list = hdulist;
  565.    else
  566.      last_hdulist->next_hdu = hdulist;
  567.    last_hdulist = hdulist;
  568.                            /* Evaluate the range of pixel data */
  569.    fits_eval_pixrange (fp, hdulist);
  570.  
  571.    /* Reposition to start of next header */
  572.    if (fseek (fp, hdulist->data_offset+hdulist->data_size, SEEK_SET) < 0)
  573.      break;
  574.  }
  575.  
  576.  return (ff);
  577. }
  578.  
  579.  
  580. /*****************************************************************************/
  581. /* #BEG-PAR                                                                  */
  582. /*                                                                           */
  583. /* Function  : fits_close - close a FITS file                                */
  584. /*                                                                           */
  585. /* Parameters:                                                               */
  586. /* FITS_FILE *ff  [I] : FITS file pointer                                    */
  587. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  588. /*                                                                           */
  589. /* #END-PAR                                                                  */
  590. /*****************************************************************************/
  591.  
  592. void fits_close (FITS_FILE *ff)
  593.  
  594. {
  595.  if (ff == NULL) FITS_VRETURN ("fits_close: Invalid parameter");
  596.  
  597.  fclose (ff->fp);
  598.  
  599.  fits_delete_filestruct (ff);
  600. }
  601.  
  602.  
  603. /*****************************************************************************/
  604. /* #BEG-PAR                                                                  */
  605. /*                                                                           */
  606. /* Function  : fits_add_hdu - add a HDU to the file                          */
  607. /*                                                                           */
  608. /* Parameters:                                                               */
  609. /* FITS_FILE *ff  [I] : FITS file pointer                                    */
  610. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  611. /*                                                                           */
  612. /* Adds a new HDU to the list kept in ff. A pointer to the new HDU is        */
  613. /* returned. On failure, a NULL-pointer is returned.                         */
  614. /*                                                                           */
  615. /* #END-PAR                                                                  */
  616. /*****************************************************************************/
  617.  
  618. FITS_HDU_LIST *fits_add_hdu (FITS_FILE *ff)
  619.  
  620. {FITS_HDU_LIST *newhdu, *hdu;
  621.  
  622.  if (ff->openmode != 'w')
  623.    FITS_RETURN ("fits_add_hdu: file not open for writing", NULL);
  624.  
  625.  newhdu = fits_new_hdulist ();
  626.  if (newhdu == NULL) return (NULL);
  627.  
  628.  if (ff->hdu_list == NULL)
  629.  {
  630.    ff->hdu_list = newhdu;
  631.  }
  632.  else
  633.  {
  634.    hdu = ff->hdu_list;
  635.    while (hdu->next_hdu != NULL)
  636.      hdu = hdu->next_hdu;
  637.    hdu->next_hdu = newhdu;
  638.  }
  639.  
  640.  return (newhdu);
  641. }
  642.  
  643.  
  644. /*****************************************************************************/
  645. /* #BEG-PAR                                                                  */
  646. /*                                                                           */
  647. /* Function  : fits_add_card - add a card to the HDU                         */
  648. /*                                                                           */
  649. /* Parameters:                                                               */
  650. /* FITS_HDU_LIST *hdulist [I] : HDU listr                                    */
  651. /* char *card             [I] : card to add                                  */
  652. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  653. /*                                                                           */
  654. /* The card must follow the standards of FITS. The card must not use a       */
  655. /* keyword that is written using *hdulist itself. On success 0 is returned.  */
  656. /* On failure -1 is returned.                                                */
  657. /*                                                                           */
  658. /* #END-PAR                                                                  */
  659. /*****************************************************************************/
  660.  
  661. int fits_add_card (FITS_HDU_LIST *hdulist, char *card)
  662.  
  663. {int k;
  664.  
  665.  if (hdulist->naddcards >= FITS_NADD_CARDS) return (-1);
  666.  
  667.  k = strlen (card);
  668.  if (k < FITS_CARD_SIZE)
  669.  {
  670.    memset (&(hdulist->addcards[hdulist->naddcards][k]), ' ', FITS_CARD_SIZE-k);
  671.    memcpy (hdulist->addcards[(hdulist->naddcards)++], card, k);
  672.  }
  673.  else
  674.  {
  675.    memcpy (hdulist->addcards[(hdulist->naddcards)++], card, FITS_CARD_SIZE);
  676.  }
  677.  return (0);
  678. }
  679.  
  680.  
  681. /*****************************************************************************/
  682. /* #BEG-PAR                                                                  */
  683. /*                                                                           */
  684. /* Function  : fits_print_header - print the internal representation         */
  685. /*                                 of a single header                        */
  686. /* Parameters:                                                               */
  687. /* FITS_HDU_LIST *hdr  [I] : pointer to the header                           */
  688. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  689. /*                                                                           */
  690. /* #END-PAR                                                                  */
  691. /*****************************************************************************/
  692.  
  693. void fits_print_header (FITS_HDU_LIST *hdr)
  694.  
  695. {int k;
  696.  
  697.  if (hdr->used.simple)
  698.    printf ("Content of SIMPLE-header:\n");
  699.  else
  700.    printf ("Content of XTENSION-header %s:\n", hdr->xtension);
  701.  printf ("header_offset : %ld\n", hdr->header_offset);
  702.  printf ("data_offset   : %ld\n", hdr->data_offset);
  703.  printf ("data_size     : %ld\n", hdr->data_size);
  704.  printf ("used data_size: %ld\n", hdr->udata_size);
  705.  printf ("bytes p.pixel : %d\n", hdr->bpp);
  706.  printf ("pixmin        : %f\n", hdr->pixmin);
  707.  printf ("pixmax        : %f\n", hdr->pixmax);
  708.  
  709.  printf ("naxis         : %d\n", hdr->naxis);
  710.  for (k = 1; k <= hdr->naxis; k++)
  711.    printf ("naxis%-3d      : %d\n", k, hdr->naxisn[k-1]);
  712.  
  713.  printf ("bitpix        : %d\n", hdr->bitpix);
  714.  
  715.  if (hdr->used.blank)
  716.    printf ("blank         : %ld\n", hdr->blank);
  717.  else
  718.    printf ("blank         : not used\n");
  719.  
  720.  if (hdr->used.datamin)
  721.    printf ("datamin       : %f\n", hdr->datamin);
  722.  else
  723.    printf ("datamin       : not used\n");
  724.  if (hdr->used.datamax)
  725.    printf ("datamax       : %f\n", hdr->datamax);
  726.  else
  727.    printf ("datamax       : not used\n");
  728.  
  729.  if (hdr->used.gcount)
  730.    printf ("gcount        : %ld\n", hdr->gcount);
  731.  else
  732.    printf ("gcount        : not used\n");
  733.  if (hdr->used.pcount)
  734.    printf ("pcount        : %ld\n", hdr->pcount);
  735.  else
  736.    printf ("pcount        : not used\n");
  737.  
  738.  if (hdr->used.bscale)
  739.    printf ("bscale        : %f\n", hdr->bscale);
  740.  else
  741.    printf ("bscale        : not used\n");
  742.  if (hdr->used.bzero)
  743.    printf ("bzero         : %f\n", hdr->bzero);
  744.  else
  745.    printf ("bzero         : not used\n");
  746. }
  747.  
  748.  
  749. /*****************************************************************************/
  750. /* #BEG-PAR                                                                  */
  751. /*                                                                           */
  752. /* Function  : fits_read_header - (local) read FITS header                   */
  753. /*                                                                           */
  754. /* Parameters:                                                               */
  755. /* FILE *fp   [I] : file pointer                                             */
  756. /* int  *nrec [O] : number of records read                                   */
  757. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  758. /*                                                                           */
  759. /* Reads in all header records up to the record that keeps the END-card.     */
  760. /* A pointer to the record list is returned on success.                      */
  761. /* On failure, a NULL-pointer is returned.                                   */
  762. /*                                                                           */
  763. /* #END-PAR                                                                  */
  764. /*****************************************************************************/
  765.  
  766. static FITS_RECORD_LIST *fits_read_header (FILE *fp, int *nrec)
  767.  
  768. {char record[FITS_RECORD_SIZE];
  769.  FITS_RECORD_LIST *start_list = NULL, *cu_record = NULL, *new_record;
  770.  FITS_DATA *fdat;
  771.  int k, simple, xtension;
  772.  
  773.  *nrec = 0;
  774.  
  775.  k = fread (record, 1, FITS_RECORD_SIZE, fp);
  776.  if (k != FITS_RECORD_SIZE)
  777.    FITS_RETURN ("fits_read_header: Error in read of first record", NULL);
  778.  
  779.  simple = (strncmp (record, "SIMPLE  ", 8) == 0);
  780.  xtension = (strncmp (record, "XTENSION", 8) == 0);
  781.  if ((!simple) && (!xtension))
  782.    FITS_RETURN ("fits_read_header: Missing keyword SIMPLE or XTENSION", NULL);
  783.  
  784.  if (simple)
  785.  {
  786.    fdat = fits_decode_card (record, typ_fbool);
  787.    if (fdat && !fdat->fbool)
  788.      fits_set_error ("fits_read_header (warning): keyword SIMPLE does not have\
  789.  value T");
  790.  }
  791.  
  792.  for (;;)   /* Process all header records */
  793.  {
  794.    new_record = (FITS_RECORD_LIST *)malloc (sizeof (FITS_RECORD_LIST));
  795.    if (new_record == NULL)
  796.    {
  797.      fits_delete_recordlist (start_list);
  798.      FITS_RETURN ("fits_read_header: Not enough memory", NULL);
  799.    }
  800.    memcpy (new_record->data, record, FITS_RECORD_SIZE);
  801.    new_record->next_record = NULL;
  802.    (*nrec)++;
  803.  
  804.    if (start_list == NULL)      /* Add new record to the list */
  805.      start_list = new_record;
  806.    else
  807.      cu_record->next_record = new_record;
  808.  
  809.    cu_record = new_record;
  810.                                 /* Was this the last record ? */
  811.    if (fits_search_card (cu_record, "END") != NULL) break;
  812.  
  813.     k = fread (record, 1, FITS_RECORD_SIZE, fp);
  814.     if (k != FITS_RECORD_SIZE)
  815.       FITS_RETURN ("fits_read_header: Error in read of record", NULL);
  816.  }
  817.  return (start_list);
  818. }
  819.  
  820.  
  821. /*****************************************************************************/
  822. /* #BEG-PAR                                                                  */
  823. /*                                                                           */
  824. /* Function  : fits_write_header - write a FITS header                       */
  825. /*                                                                           */
  826. /* Parameters:                                                               */
  827. /* FITS_FILE *ff [I] : FITS-file pointer                                     */
  828. /* FITS_HDU_LIST [I] : pointer to header                                     */
  829. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  830. /*                                                                           */
  831. /* Writes a header to the file. On success, 0 is returned. On failure,       */
  832. /* -1 is returned.                                                           */
  833. /*                                                                           */
  834. /* #END-PAR                                                                  */
  835. /*****************************************************************************/
  836.  
  837. int fits_write_header (FITS_FILE *ff, FITS_HDU_LIST *hdulist)
  838.  
  839. {int numcards;
  840.  int k;
  841.  
  842.  if (ff->openmode != 'w')
  843.    FITS_RETURN ("fits_write_header: file not open for writing", -1);
  844.  
  845.  numcards = 0;
  846.  
  847.  if (hdulist->used.simple)
  848.  {
  849.    FITS_WRITE_BOOLCARD (ff->fp, "SIMPLE", 1);
  850.    numcards++;
  851.  }
  852.  else if (hdulist->used.xtension)
  853.  {
  854.    FITS_WRITE_STRINGCARD (ff->fp, "XTENSION", hdulist->xtension);
  855.    numcards++;
  856.  }
  857.  
  858.  FITS_WRITE_LONGCARD (ff->fp, "BITPIX", hdulist->bitpix);
  859.  numcards++;
  860.  
  861.  FITS_WRITE_LONGCARD (ff->fp, "NAXIS", hdulist->naxis);
  862.  numcards++;
  863.  
  864.  for (k = 0; k < hdulist->naxis; k++)
  865.  {char naxisn[10];
  866.    sprintf (naxisn, "NAXIS%d", k+1);
  867.    FITS_WRITE_LONGCARD (ff->fp, naxisn, hdulist->naxisn[k]);
  868.    numcards++;
  869.  }
  870.  
  871.  if (hdulist->used.extend)
  872.  {
  873.    FITS_WRITE_BOOLCARD (ff->fp, "EXTEND", hdulist->extend);
  874.    numcards++;
  875.  }
  876.  
  877.  if (hdulist->used.groups)
  878.  {
  879.    FITS_WRITE_BOOLCARD (ff->fp, "GROUPS", hdulist->groups);
  880.    numcards++;
  881.  }
  882.  
  883.  if (hdulist->used.pcount)
  884.  {
  885.    FITS_WRITE_LONGCARD (ff->fp, "PCOUNT", hdulist->pcount);
  886.    numcards++;
  887.  }
  888.  if (hdulist->used.gcount)
  889.  {
  890.    FITS_WRITE_LONGCARD (ff->fp, "GCOUNT", hdulist->gcount);
  891.    numcards++;
  892.  }
  893.  
  894.  if (hdulist->used.bzero)
  895.  {
  896.    FITS_WRITE_DOUBLECARD (ff->fp, "BZERO", hdulist->bzero);
  897.    numcards++;
  898.  }
  899.  if (hdulist->used.bscale)
  900.  {
  901.    FITS_WRITE_DOUBLECARD (ff->fp, "BSCALE", hdulist->bscale);
  902.    numcards++;
  903.  }
  904.  
  905.  if (hdulist->used.datamin)
  906.  {
  907.    FITS_WRITE_DOUBLECARD (ff->fp, "DATAMIN", hdulist->datamin);
  908.    numcards++;
  909.  }
  910.  if (hdulist->used.datamax)
  911.  {
  912.    FITS_WRITE_DOUBLECARD (ff->fp, "DATAMAX", hdulist->datamax);
  913.    numcards++;
  914.  }
  915.  
  916.  if (hdulist->used.blank)
  917.  {
  918.    FITS_WRITE_LONGCARD (ff->fp, "BLANK", hdulist->blank);
  919.    numcards++;
  920.  }
  921.  
  922.  /* Write additional cards */
  923.  if (hdulist->naddcards > 0)
  924.  {
  925.    fwrite (hdulist->addcards, FITS_CARD_SIZE, hdulist->naddcards, ff->fp);
  926.    numcards += hdulist->naddcards;
  927.  }
  928.  
  929.  FITS_WRITE_CARD (ff->fp, "END");
  930.  numcards++;
  931.  
  932.  k = (numcards*FITS_CARD_SIZE) % FITS_RECORD_SIZE;
  933.  if (k)  /* Must the record be filled up ? */
  934.  {
  935.    while (k++ < FITS_RECORD_SIZE)
  936.      putc (' ', ff->fp);
  937.  }
  938.  
  939.  return (ferror (ff->fp) ? -1 : 0);
  940. }
  941.  
  942.  
  943. /*****************************************************************************/
  944. /* #BEG-PAR                                                                  */
  945. /*                                                                           */
  946. /* Function  : fits_decode_header - (local) decode a header                  */
  947. /*                                                                           */
  948. /* Parameters:                                                               */
  949. /* FITS_RECORD_LIST *hdr  [I] : the header record list                       */
  950. /* long hdr_offset        [I] : fileposition of header                       */
  951. /* long dat_offset        [I] : fileposition of data (end of header)         */
  952. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  953. /*                                                                           */
  954. /* The function decodes the mostly used data within the header and generates */
  955. /* a FITS_HDU_LIST-entry. On failure, a NULL-pointer is returned.            */
  956. /*                                                                           */
  957. /* #END-PAR                                                                  */
  958. /*****************************************************************************/
  959.  
  960. static FITS_HDU_LIST *fits_decode_header (FITS_RECORD_LIST *hdr,
  961.                         long hdr_offset, long dat_offset)
  962.  
  963. {FITS_HDU_LIST *hdulist;
  964.  FITS_DATA *fdat;
  965.  char errmsg[80], key[9];
  966.  int k, bpp, random_groups;
  967.  long mul_axis, data_size, bitpix_supported;
  968.  
  969. #define FITS_DECODE_CARD(mhdr,mkey,mfdat,mtyp) \
  970.  {strcpy (key, mkey); \
  971.   mfdat = fits_decode_card (fits_search_card (mhdr, mkey), mtyp); \
  972.   if (mfdat == NULL) goto err_missing; }
  973.  
  974. #define FITS_TRY_CARD(mhdr,mhdu,mkey,mvar,mtyp,unionvar) \
  975.  {FITS_DATA *mfdat = fits_decode_card (fits_search_card (mhdr,mkey), mtyp); \
  976.   mhdu->used.mvar = (mfdat != NULL); \
  977.   if (mhdu->used.mvar) mhdu->mvar = mfdat->unionvar; }
  978.  
  979.  hdulist = fits_new_hdulist ();
  980.  if (hdulist == NULL)
  981.    FITS_RETURN ("fits_decode_header: Not enough memory", NULL);
  982.  
  983.  /* Initialize the header data */
  984.  hdulist->header_offset = hdr_offset;
  985.  hdulist->data_offset = dat_offset;
  986.  
  987.  hdulist->used.simple = (strncmp (hdr->data, "SIMPLE  ", 8) == 0);
  988.  hdulist->used.xtension = (strncmp (hdr->data, "XTENSION", 8) == 0);
  989.  if (hdulist->used.xtension)
  990.  {
  991.    fdat = fits_decode_card (fits_search_card (hdr, "XTENSION"), typ_fstring);
  992.    strcpy (hdulist->xtension, fdat->fstring);
  993.  }
  994.  
  995.  FITS_DECODE_CARD (hdr, "NAXIS", fdat, typ_flong);
  996.  hdulist->naxis = fdat->flong;
  997.  
  998.  FITS_DECODE_CARD (hdr, "BITPIX", fdat, typ_flong);
  999.  bpp = hdulist->bitpix = (int)fdat->flong;
  1000.  if (   (bpp != 8) && (bpp != 16) && (bpp != 32)
  1001.      && (bpp != -32) && (bpp != -64))
  1002.  {
  1003.    strcpy (errmsg, "fits_decode_header: Invalid BITPIX-value");
  1004.    goto err_return;
  1005.  }
  1006.  if (bpp < 0) bpp = -bpp;
  1007.  bpp /= 8;
  1008.  hdulist->bpp = bpp;
  1009.  
  1010.  FITS_TRY_CARD (hdr, hdulist, "GCOUNT", gcount, typ_flong, flong);
  1011.  FITS_TRY_CARD (hdr, hdulist, "PCOUNT", pcount, typ_flong, flong);
  1012.  
  1013.  FITS_TRY_CARD (hdr, hdulist, "GROUPS", groups, typ_fbool, fbool);
  1014.  random_groups = hdulist->used.groups && hdulist->groups;
  1015.  
  1016.  FITS_TRY_CARD (hdr, hdulist, "EXTEND", extend, typ_fbool, fbool);
  1017.  
  1018.  if (hdulist->used.xtension)  /* Extension requires GCOUNT and PCOUNT */
  1019.  {
  1020.    if ((!hdulist->used.gcount) || (!hdulist->used.pcount))
  1021.    {
  1022.      strcpy (errmsg, "fits_decode_header: Missing GCOUNT/PCOUNT for XTENSION");
  1023.      goto err_return;
  1024.    }
  1025.  }
  1026.  
  1027.  mul_axis = 1;
  1028.  
  1029.  /* Find all NAXISx-cards */
  1030.  for (k = 1; k <= FITS_MAX_AXIS; k++)
  1031.  {char naxisn[9];
  1032.  
  1033.    sprintf (naxisn, "NAXIS%-3d", k);
  1034.    fdat = fits_decode_card (fits_search_card (hdr, naxisn), typ_flong);
  1035.    if (fdat == NULL)
  1036.    {
  1037.      k--;   /* Save the last NAXISk read */
  1038.      break;
  1039.    }
  1040.    hdulist->naxisn[k-1] = (int)fdat->flong;
  1041.    if (hdulist->naxisn[k-1] < 0)
  1042.    {
  1043.      strcpy (errmsg, "fits_decode_header: Negative value in NAXISn");
  1044.      goto err_return;
  1045.    }
  1046.    if ((k == 1) && (random_groups))
  1047.    {
  1048.      if (hdulist->naxisn[0] != 0)
  1049.      {
  1050.        strcpy (errmsg, "fits_decode_header: Random groups with NAXIS1 != 0");
  1051.        goto err_return;
  1052.      }
  1053.    }
  1054.    else
  1055.      mul_axis *= hdulist->naxisn[k-1];
  1056.  }
  1057.  
  1058.  if ((hdulist->naxis > 0) && (k < hdulist->naxis))
  1059.  {
  1060.    strcpy (errmsg, "fits_decode_card: Not enough NAXISn-cards");
  1061.    goto err_return;
  1062.  }
  1063.  
  1064.  /* If we have only one dimension, just set the second to size one. */
  1065.  /* So we dont have to check for naxis < 2 in some places. */
  1066.  if (hdulist->naxis < 2)
  1067.    hdulist->naxisn[1] = 1;
  1068.  if (hdulist->naxis < 1)
  1069.  {
  1070.    mul_axis = 0;
  1071.    hdulist->naxisn[0] = 1;
  1072.  }
  1073.  
  1074.  if (hdulist->used.xtension)
  1075.    data_size = bpp*hdulist->gcount*(hdulist->pcount + mul_axis);
  1076.  else
  1077.    data_size = bpp*mul_axis;
  1078.  hdulist->udata_size = data_size;  /* Used data size without padding */
  1079.  
  1080.  /* Datasize must be a multiple of the FITS logical record size */
  1081.  data_size = (data_size + FITS_RECORD_SIZE - 1) / FITS_RECORD_SIZE;
  1082.  data_size *= FITS_RECORD_SIZE;
  1083.  hdulist->data_size = data_size;
  1084.  
  1085.  
  1086.  FITS_TRY_CARD (hdr, hdulist, "BLANK", blank, typ_flong, flong);
  1087.  
  1088.  FITS_TRY_CARD (hdr, hdulist, "DATAMIN", datamin, typ_fdouble, fdouble);
  1089.  FITS_TRY_CARD (hdr, hdulist, "DATAMAX", datamax, typ_fdouble, fdouble);
  1090.  
  1091.  FITS_TRY_CARD (hdr, hdulist, "BZERO", bzero, typ_fdouble, fdouble);
  1092.  FITS_TRY_CARD (hdr, hdulist, "BSCALE", bscale, typ_fdouble, fdouble);
  1093.  
  1094.  /* Evaluate number of interpretable images for this HDU */
  1095.  hdulist->numpic = 0;
  1096.  
  1097.  /* We must support this format */
  1098.  bitpix_supported =    (hdulist->bitpix > 0)
  1099.                     || (   (hdulist->bitpix == -64)
  1100.                         && (fits_ieee64_intel || fits_ieee64_motorola))
  1101.                     || (   (hdulist->bitpix == -32)
  1102.                         && (   fits_ieee32_intel || fits_ieee32_motorola
  1103.                             || fits_ieee64_intel || fits_ieee64_motorola));
  1104.    
  1105.  if (bitpix_supported)
  1106.  {
  1107.    if (hdulist->used.simple)
  1108.    {
  1109.      if (hdulist->naxis > 0)
  1110.      {
  1111.        hdulist->numpic = 1;
  1112.        for (k = 3; k <= hdulist->naxis; k++)
  1113.          hdulist->numpic *= hdulist->naxisn[k-1];
  1114.      }
  1115.    }
  1116.    else if (   hdulist->used.xtension
  1117.             && (strncmp (hdulist->xtension, "IMAGE", 5) == 0))
  1118.    {
  1119.      if (hdulist->naxis > 0)
  1120.      {
  1121.        hdulist->numpic = 1;
  1122.        for (k = 3; k <= hdulist->naxis; k++)
  1123.          hdulist->numpic *= hdulist->naxisn[k-1];
  1124.      }
  1125.    }
  1126.  }
  1127.  else
  1128.  {char msg[160];
  1129.    sprintf (msg, "fits_decode_header: IEEE floating point format required for\
  1130.  BITPIX=%d\nis not supported on this machine", hdulist->bitpix);
  1131.    fits_set_error (msg);
  1132.  }
  1133.  
  1134.  hdulist->header_record_list = hdr;  /* Add header records to the list */
  1135.  return (hdulist);
  1136.  
  1137. err_missing:
  1138.  sprintf (errmsg, "fits_decode_header: missing/invalid %s card", key);
  1139.  
  1140. err_return:
  1141.  fits_delete_hdulist (hdulist);
  1142.  fits_set_error (errmsg);
  1143.  return (NULL);
  1144.  
  1145. #undef FITS_DECODE_CARD
  1146. }
  1147.  
  1148.  
  1149. /*****************************************************************************/
  1150. /* #BEG-PAR                                                                  */
  1151. /*                                                                           */
  1152. /* Function  : fits_eval_pixrange - (local) evaluate range of pixel data     */
  1153. /*                                                                           */
  1154. /* Parameters:                                                               */
  1155. /* FILE *fp               [I] : file pointer                                 */
  1156. /* FITS_HDU_LIST *hdu     [I] : pointer to header                            */
  1157. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  1158. /*                                                                           */
  1159. /* The Function sets the values hdu->pixmin and hdu->pixmax. On success 0    */
  1160. /* is returned. On failure, -1 is returned.                                  */
  1161. /*                                                                           */
  1162. /* #END-PAR                                                                  */
  1163. /*****************************************************************************/
  1164.  
  1165. static int fits_eval_pixrange (FILE *fp, FITS_HDU_LIST *hdu)
  1166.  
  1167. {register int maxelem;
  1168. #define FITSNPIX 4096 
  1169.  unsigned char pixdat[FITSNPIX];
  1170.  int nelem, bpp;
  1171.  int blank_found = 0, nan_found = 0;
  1172.    
  1173.  if (fseek (fp, hdu->data_offset, SEEK_SET) < 0)
  1174.    FITS_RETURN ("fits_eval_pixrange: cant position file", -1);
  1175.  
  1176.  bpp = hdu->bpp;                  /* Number of bytes per pixel */
  1177.  nelem = hdu->udata_size / bpp;   /* Number of data elements */
  1178.  
  1179.  switch (hdu->bitpix)
  1180.  {
  1181.    case 8: {
  1182.      register FITS_BITPIX8 pixval;
  1183.      register unsigned char *ptr;
  1184.      FITS_BITPIX8 minval = 255, maxval = 0;
  1185.      FITS_BITPIX8 blankval;
  1186.  
  1187.      while (nelem > 0)
  1188.      {
  1189.        maxelem = sizeof (pixdat)/bpp;
  1190.        if (nelem < maxelem) maxelem = nelem;
  1191.        nelem -= maxelem;
  1192.        if (fread ((char *)pixdat, bpp, maxelem, fp) != maxelem)
  1193.          FITS_RETURN ("fits_eval_pixrange: error on read bitpix 8 data", -1);
  1194.   
  1195.        ptr = pixdat;
  1196.        if (hdu->used.blank)
  1197.        {
  1198.          blankval = (FITS_BITPIX8)hdu->blank;
  1199.          while (maxelem-- > 0)
  1200.          {
  1201.            pixval = (FITS_BITPIX8)*(ptr++);
  1202.            if (pixval != blankval)
  1203.            {
  1204.              if (pixval < minval) minval = pixval;
  1205.              else if (pixval > maxval) maxval = pixval;
  1206.            }
  1207.            else blank_found = 1;
  1208.          }
  1209.        }
  1210.        else
  1211.        {
  1212.          while (maxelem-- > 0)
  1213.          {
  1214.            pixval = (FITS_BITPIX8)*(ptr++);
  1215.            if (pixval < minval) minval = pixval;
  1216.            else if (pixval > maxval) maxval = pixval;
  1217.          }
  1218.        }
  1219.      }
  1220.      hdu->pixmin = minval;
  1221.      hdu->pixmax = maxval;
  1222.      break; }
  1223.  
  1224.    case 16: {
  1225.      register FITS_BITPIX16 pixval;
  1226.      register unsigned char *ptr;
  1227.      FITS_BITPIX16 minval = 0x7fff, maxval = ~0x7fff;
  1228.  
  1229.      while (nelem > 0)
  1230.      {
  1231.        maxelem = sizeof (pixdat)/bpp;
  1232.        if (nelem < maxelem) maxelem = nelem;
  1233.        nelem -= maxelem;
  1234.        if (fread ((char *)pixdat, bpp, maxelem, fp) != maxelem)
  1235.          FITS_RETURN ("fits_eval_pixrange: error on read bitpix 16 data", -1);
  1236.   
  1237.        ptr = pixdat;
  1238.        if (hdu->used.blank)
  1239.        {FITS_BITPIX16 blankval = (FITS_BITPIX16)hdu->blank;
  1240.  
  1241.          while (maxelem-- > 0)
  1242.          {
  1243.            FITS_GETBITPIX16 (ptr, pixval);
  1244.            ptr += 2;
  1245.            if (pixval != blankval)
  1246.            {
  1247.              if (pixval < minval) minval = pixval;
  1248.              else if (pixval > maxval) maxval = pixval;
  1249.            }
  1250.            else blank_found = 1;
  1251.          }
  1252.        }
  1253.        else
  1254.        {
  1255.          while (maxelem-- > 0)
  1256.          {
  1257.            FITS_GETBITPIX16 (ptr, pixval);
  1258.            ptr += 2;
  1259.            if (pixval < minval) minval = pixval;
  1260.            else if (pixval > maxval) maxval = pixval;
  1261.          }
  1262.        }
  1263.      }
  1264.      hdu->pixmin = minval;
  1265.      hdu->pixmax = maxval;
  1266.      break; }
  1267.  
  1268.  
  1269.    case 32: {
  1270.      register FITS_BITPIX32 pixval;
  1271.      register unsigned char *ptr;
  1272.      FITS_BITPIX32 minval = 0x7fffffff, maxval = ~0x7fffffff;
  1273.  
  1274.      while (nelem > 0)
  1275.      {
  1276.        maxelem = sizeof (pixdat)/bpp;
  1277.        if (nelem < maxelem) maxelem = nelem;
  1278.        nelem -= maxelem;
  1279.        if (fread ((char *)pixdat, bpp, maxelem, fp) != maxelem)
  1280.          FITS_RETURN ("fits_eval_pixrange: error on read bitpix 32 data", -1);
  1281.   
  1282.        ptr = pixdat;
  1283.        if (hdu->used.blank)
  1284.        {FITS_BITPIX32 blankval = (FITS_BITPIX32)hdu->blank;
  1285.  
  1286.          while (maxelem-- > 0)
  1287.          {
  1288.            FITS_GETBITPIX32 (ptr, pixval);
  1289.            ptr += 4;
  1290.            if (pixval != blankval)
  1291.            {
  1292.              if (pixval < minval) minval = pixval;
  1293.              else if (pixval > maxval) maxval = pixval;
  1294.            }
  1295.            else blank_found = 1;
  1296.          }
  1297.        }
  1298.        else
  1299.        {
  1300.          while (maxelem-- > 0)
  1301.          {
  1302.            FITS_GETBITPIX32 (ptr, pixval);
  1303.            ptr += 4;
  1304.            if (pixval < minval) minval = pixval;
  1305.            else if (pixval > maxval) maxval = pixval;
  1306.          }
  1307.        }
  1308.      }
  1309.      hdu->pixmin = minval;
  1310.      hdu->pixmax = maxval;
  1311.      break; }
  1312.  
  1313.    case -32: {
  1314.      register FITS_BITPIXM32 pixval;
  1315.      register unsigned char *ptr;
  1316.      FITS_BITPIXM32 minval, maxval;
  1317.      int first = 1;
  1318.  
  1319.      /* initialize */
  1320.  
  1321.      pixval = 0;
  1322.      minval = 0;
  1323.      maxval = 0;
  1324.  
  1325.      while (nelem > 0)
  1326.      {
  1327.        maxelem = sizeof (pixdat)/bpp;
  1328.        if (nelem < maxelem) maxelem = nelem;
  1329.        nelem -= maxelem;
  1330.        if (fread ((char *)pixdat, bpp, maxelem, fp) != maxelem)
  1331.          FITS_RETURN ("fits_eval_pixrange: error on read bitpix -32 data", -1);
  1332.   
  1333.        ptr = pixdat;
  1334.        while (maxelem-- > 0)
  1335.        {
  1336.          if (!fits_nan_32 (ptr))
  1337.          {
  1338.            FITS_GETBITPIXM32 (ptr, pixval);
  1339.            ptr += 4;
  1340.            if (first)
  1341.            {
  1342.              first = 0;
  1343.              minval = maxval = pixval;
  1344.            }
  1345.            else if (pixval < minval) { minval = pixval; }
  1346.            else if (pixval > maxval) { maxval = pixval; }
  1347.          }
  1348.          else nan_found = 1;
  1349.        }
  1350.      }
  1351.      hdu->pixmin = minval;
  1352.      hdu->pixmax = maxval;
  1353.      break; }
  1354.  
  1355.    case -64: {
  1356.      register FITS_BITPIXM64 pixval;
  1357.      register unsigned char *ptr;
  1358.      FITS_BITPIXM64 minval, maxval;
  1359.      int first = 1;
  1360.  
  1361.      /* initialize */
  1362.  
  1363.      minval = 0;
  1364.      maxval = 0;
  1365.  
  1366.      while (nelem > 0)
  1367.      {
  1368.        maxelem = sizeof (pixdat)/bpp;
  1369.        if (nelem < maxelem) maxelem = nelem;
  1370.        nelem -= maxelem;
  1371.        if (fread ((char *)pixdat, bpp, maxelem, fp) != maxelem)
  1372.          FITS_RETURN ("fits_eval_pixrange: error on read bitpix -64 data", -1);
  1373.   
  1374.        ptr = pixdat;
  1375.        while (maxelem-- > 0)
  1376.        {
  1377.          if (!fits_nan_64 (ptr))
  1378.          {
  1379.            FITS_GETBITPIXM64 (ptr, pixval);
  1380.            ptr += 8;
  1381.            if (first)
  1382.            {
  1383.              first = 0;
  1384.              minval = maxval = pixval;
  1385.            }
  1386.            else if (pixval < minval) { minval = pixval; }
  1387.            else if (pixval > maxval) { maxval = pixval; }
  1388.          }
  1389.          else nan_found = 1;
  1390.        }
  1391.      }
  1392.      hdu->pixmin = minval;
  1393.      hdu->pixmax = maxval;
  1394.      break; }
  1395.  }
  1396.  if (nan_found) hdu->used.nan_value = 1;
  1397.  if (blank_found) hdu->used.blank_value = 1;
  1398.  
  1399.  return (0);
  1400. }
  1401.  
  1402.  
  1403. /*****************************************************************************/
  1404. /* #BEG-PAR                                                                  */
  1405. /*                                                                           */
  1406. /* Function  : fits_decode_card - decode a card                              */
  1407. /*                                                                           */
  1408. /* Parameters:                                                               */
  1409. /* const char *card   [I] : pointer to card image                            */
  1410. /* FITS_DATA_TYPES data_type [I] : datatype to decode                        */
  1411. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  1412. /*                                                                           */
  1413. /* Decodes a card and returns a pointer to the union, keeping the data.      */
  1414. /* If card is NULL or on failure, a NULL-pointer is returned.                */
  1415. /* If the card does not have the value indicator, an error is generated,     */
  1416. /* but its tried to decode the card. The data is only valid up to the next   */
  1417. /* call of the function.                                                     */
  1418. /*                                                                           */
  1419. /* #END-PAR                                                                  */
  1420. /*****************************************************************************/
  1421.  
  1422. FITS_DATA *fits_decode_card (const char *card, FITS_DATA_TYPES data_type)
  1423.  
  1424. {static FITS_DATA data;
  1425.  long l_long;
  1426.  double l_double;
  1427.  char l_card[FITS_CARD_SIZE+1], msg[256];
  1428.  char *cp = ident, *dst, *end;
  1429.  
  1430.  if (card == NULL) return (NULL);
  1431.  
  1432.  memcpy (l_card, card, FITS_CARD_SIZE);
  1433.  l_card[FITS_CARD_SIZE] = '\0';
  1434.  
  1435.  if (strncmp (card+8, "= ", 2) != 0)
  1436.  {
  1437.    sprintf (msg, "fits_decode_card (warning): Missing value indicator\
  1438.  '= ' for %8.8s", l_card);
  1439.    fits_set_error (msg);
  1440.  }
  1441.  
  1442.  switch (data_type)
  1443.  {
  1444.    case typ_bitpix8:
  1445.      data.bitpix8 = (FITS_BITPIX8)(l_card[10]);
  1446.      break;
  1447.  
  1448.    case typ_bitpix16:
  1449.      if (sscanf (l_card+10, "%ld", &l_long) != 1)
  1450.        FITS_RETURN ("fits_decode_card: error decoding typ_bitpix16", NULL);
  1451.      data.bitpix16 = (FITS_BITPIX16)l_long;
  1452.      break;
  1453.  
  1454.    case typ_bitpix32:
  1455.      if (sscanf (l_card+10, "%ld", &l_long) != 1)
  1456.        FITS_RETURN ("fits_decode_card: error decoding typ_bitpix32", NULL);
  1457.      data.bitpix32 = (FITS_BITPIX32)l_long;
  1458.      break;
  1459.  
  1460.    case typ_bitpixm32:
  1461.      if (sscanf (l_card+10, "%lf", &l_double) != 1)
  1462.        FITS_RETURN ("fits_decode_card: error decoding typ_bitpixm32", NULL);
  1463.      data.bitpixm32 = (FITS_BITPIXM32)l_double;
  1464.      break;
  1465.  
  1466.    case typ_bitpixm64:
  1467.      if (sscanf (l_card+10, "%lf", &l_double) != 1)
  1468.        FITS_RETURN ("fits_decode_card: error decoding typ_bitpixm64", NULL);
  1469.      data.bitpixm64 = (FITS_BITPIXM64)l_double;
  1470.      break;
  1471.      
  1472.    case typ_fbool:
  1473.      cp = l_card+10;
  1474.      while (*cp == ' ') cp++;
  1475.      if (*cp == 'T') data.fbool = 1;
  1476.      else if (*cp == 'F') data.fbool = 0;
  1477.      else FITS_RETURN ("fits_decode_card: error decoding typ_fbool", NULL);
  1478.      break;
  1479.  
  1480.    case typ_flong:
  1481.      if (sscanf (l_card+10, "%ld", &l_long) != 1)
  1482.        FITS_RETURN ("fits_decode_card: error decoding typ_flong", NULL);
  1483.      data.flong = (FITS_BITPIX32)l_long;
  1484.      break;
  1485.  
  1486.    case typ_fdouble:
  1487.      if (sscanf (l_card+10, "%lf", &l_double) != 1)
  1488.        FITS_RETURN ("fits_decode_card: error decoding typ_fdouble", NULL);
  1489.      data.fdouble = (FITS_BITPIXM32)l_double;
  1490.      break;
  1491.  
  1492.    case typ_fstring:
  1493.      cp = l_card+10;
  1494.      if (*cp != '\'')
  1495.        FITS_RETURN ("fits_decode_card: missing \' decoding typ_fstring", NULL);
  1496.  
  1497.      dst = data.fstring;
  1498.      cp++;
  1499.      end = l_card+FITS_CARD_SIZE-1;
  1500.      for (;;)   /* Search for trailing quote */
  1501.      {
  1502.        if (*cp != '\'')    /* All characters but quote are used. */
  1503.        {
  1504.          *(dst++) = *cp;
  1505.        }
  1506.        else                /* Maybe there is a quote in the string */
  1507.        {
  1508.          if (cp >= end) break;  /* End of card ? finished */
  1509.          if (*(cp+1) != '\'') break;
  1510.          *(dst++) = *(cp++);
  1511.        }
  1512.        if (cp >= end) break;
  1513.        cp++;
  1514.      }
  1515.      *dst = '\0';
  1516.      break;
  1517.  }
  1518.  return (&data);
  1519. }
  1520.  
  1521.  
  1522. /*****************************************************************************/
  1523. /* #BEG-PAR                                                                  */
  1524. /*                                                                           */
  1525. /* Function  : fits_search_card - search a card in the record list           */
  1526. /*                                                                           */
  1527. /* Parameters:                                                               */
  1528. /* FITS_RECORD_LIST *rl  [I] : record list to search                         */
  1529. /* char *keyword         [I] : keyword identifying the card                  */
  1530. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  1531. /*                                                                           */
  1532. /* A card is searched in the reord list. Only the first eight characters of  */
  1533. /* keyword are significant. If keyword is less than 8 characters, its filled */
  1534. /* with blanks.                                                              */
  1535. /* If the card is found, a pointer to the card is returned.                  */
  1536. /* The pointer does not point to a null-terminated string. Only the next     */
  1537. /* 80 bytes are allowed to be read.                                          */
  1538. /* On failure a NULL-pointer is returned.                                    */
  1539. /*                                                                           */
  1540. /* #END-PAR                                                                  */
  1541. /*****************************************************************************/
  1542.  
  1543. char *fits_search_card (FITS_RECORD_LIST *rl, char *keyword)
  1544.  
  1545. {int key_len, k;
  1546.  char *card;
  1547.  char key[9];
  1548.  
  1549.  key_len = strlen (keyword);
  1550.  if (key_len > 8) key_len = 8;
  1551.  if (key_len == 0)
  1552.    FITS_RETURN ("fits_search_card: Invalid parameter", NULL);
  1553.  
  1554.  strcpy (key, "        ");
  1555.  memcpy (key, keyword, key_len);
  1556.  
  1557.  while (rl != NULL)
  1558.  {
  1559.    card = (char *)rl->data;
  1560.    for (k = 0; k < FITS_RECORD_SIZE / FITS_CARD_SIZE; k++)
  1561.    {
  1562.      if (strncmp (card, key, 8) == 0) return (card);
  1563.      card += FITS_CARD_SIZE;
  1564.    }
  1565.    rl = rl->next_record;
  1566.  }
  1567.  return (NULL);
  1568. }
  1569.  
  1570.  
  1571. /*****************************************************************************/
  1572. /* #BEG-PAR                                                                  */
  1573. /*                                                                           */
  1574. /* Function  : fits_image_info - get information about an image              */
  1575. /*                                                                           */
  1576. /* Parameters:                                                               */
  1577. /* FITS_FILE *ff         [I] : FITS file structure                           */
  1578. /* int picind            [I] : Index of picture in file (1,2,...)            */
  1579. /* int *hdupicind        [O] : Index of picture in HDU (1,2,...)             */
  1580. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  1581. /*                                                                           */
  1582. /* The function returns on success a pointer to a FITS_HDU_LIST. hdupicind   */
  1583. /* then gives the index of the image within the HDU.                         */
  1584. /* On failure, NULL is returned.                                             */
  1585. /*                                                                           */
  1586. /* #END-PAR                                                                  */
  1587. /*****************************************************************************/
  1588.  
  1589. FITS_HDU_LIST *fits_image_info (FITS_FILE *ff, int picind, int *hdupicind)
  1590.  
  1591. {FITS_HDU_LIST *hdulist;
  1592.  int firstpic, lastpic;
  1593.  
  1594.  if (ff == NULL)
  1595.    FITS_RETURN ("fits_image_info: ff is NULL", NULL);
  1596.  
  1597.  if (ff->openmode != 'r')
  1598.    FITS_RETURN ("fits_image_info: file not open for reading", NULL);
  1599.    
  1600.  if ((picind < 1) || (picind > ff->n_pic))
  1601.    FITS_RETURN ("fits_image_info: picind out of range", NULL);
  1602.  
  1603.  firstpic = 1;
  1604.  for (hdulist = ff->hdu_list; hdulist != NULL; hdulist = hdulist->next_hdu)
  1605.  {
  1606.    if (hdulist->numpic <= 0) continue;
  1607.    lastpic = firstpic+hdulist->numpic-1;
  1608.    if (picind <= lastpic)  /* Found image in current HDU ? */
  1609.      break;
  1610.  
  1611.    firstpic = lastpic+1;
  1612.  }
  1613.  *hdupicind = picind - firstpic + 1;
  1614.  return (hdulist);
  1615. }
  1616.  
  1617.  
  1618. /*****************************************************************************/
  1619. /* #BEG-PAR                                                                  */
  1620. /*                                                                           */
  1621. /* Function  : fits_seek_image - position to a specific image                */
  1622. /*                                                                           */
  1623. /* Parameters:                                                               */
  1624. /* FITS_FILE *ff         [I] : FITS file structure                           */
  1625. /* int picind            [I] : Index of picture to seek (1,2,...)            */
  1626. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  1627. /*                                                                           */
  1628. /* The function positions the file pointer to a specified image.             */
  1629. /* The function returns on success a pointer to a FITS_HDU_LIST. This pointer*/
  1630. /* must also be used when reading data from the image.                       */
  1631. /* On failure, NULL is returned.                                             */
  1632. /*                                                                           */
  1633. /* #END-PAR                                                                  */
  1634. /*****************************************************************************/
  1635.  
  1636. FITS_HDU_LIST *fits_seek_image (FITS_FILE *ff, int picind)
  1637.  
  1638. {FITS_HDU_LIST *hdulist;
  1639.  int hdupicind;
  1640.  long offset, pic_size;
  1641.  
  1642.  hdulist = fits_image_info (ff, picind, &hdupicind);
  1643.  if (hdulist == NULL) return (NULL);
  1644.  
  1645.  pic_size = hdulist->bpp * hdulist->naxisn[0] * hdulist->naxisn[1];
  1646.  offset = hdulist->data_offset + (hdupicind-1)*pic_size;
  1647.  if (fseek (ff->fp, offset, SEEK_SET) < 0)
  1648.    FITS_RETURN ("fits_seek_image: Unable to position to image", NULL);
  1649.  
  1650.  return (hdulist);
  1651. }
  1652.  
  1653.  
  1654. /*****************************************************************************/
  1655. /* #BEG-PAR                                                                  */
  1656. /*                                                                           */
  1657. /* Function  : fits_read_pixel - read pixel values from a file               */
  1658. /*                                                                           */
  1659. /* Parameters:                                                               */
  1660. /* FITS_FILE *ff           [I] : FITS file structure                         */
  1661. /* FITS_HDU_LIST *hdulist  [I] : pointer to hdulist that describes image     */
  1662. /* int npix                [I] : number of pixel values to read              */
  1663. /* FITS_PIX_TRANSFORM *trans [I]: pixel transformation                       */
  1664. /* void *buf               [O] : buffer where to place transformed pixels    */
  1665. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  1666. /*                                                                           */
  1667. /* The function reads npix pixel values from the file, transforms them       */
  1668. /* checking for blank/NaN pixels and stores the transformed values in buf.   */
  1669. /* The number of transformed pixels is returned. If the returned value is    */
  1670. /* less than npix (or even -1), an error has occured.                        */
  1671. /* hdulist must be a pointer returned by fits_seek_image(). Before starting  */
  1672. /* to read an image, fits_seek_image() must be called. Even for successive   */
  1673. /* images.                                                                   */
  1674. /*                                                                           */
  1675. /* #END-PAR                                                                  */
  1676. /*****************************************************************************/
  1677.  
  1678. int fits_read_pixel (FITS_FILE *ff, FITS_HDU_LIST *hdulist, int npix,
  1679.                      FITS_PIX_TRANSFORM *trans, void *buf)
  1680.  
  1681. {double offs, scale;
  1682.  double datadiff, pixdiff;
  1683.  unsigned char pixbuffer[4096], *pix, *cdata;
  1684.  unsigned char creplace;
  1685.  int transcount = 0;
  1686.  long tdata, tmin, tmax;
  1687.  int maxelem;
  1688.  FITS_BITPIX8 bp8, bp8blank;
  1689.  FITS_BITPIX16 bp16, bp16blank;
  1690.  FITS_BITPIX32 bp32, bp32blank;
  1691.  FITS_BITPIXM32 bpm32;
  1692.  FITS_BITPIXM64 bpm64;
  1693.  
  1694.  /* initialize */
  1695.  
  1696.  bpm32 = 0;
  1697.  
  1698.  if (ff->openmode != 'r') return (-1);   /* Not open for reading */
  1699.  if (trans->dsttyp != 'c') return (-1);  /* Currently we only return chars */
  1700.  if (npix <= 0) return (npix);
  1701.  
  1702.  datadiff = trans->datamax - trans->datamin;
  1703.  pixdiff = trans->pixmax - trans->pixmin;
  1704.  
  1705.  offs = trans->datamin - trans->pixmin*datadiff/pixdiff;
  1706.  scale = datadiff / pixdiff;
  1707.  
  1708.  tmin = (long)trans->datamin;
  1709.  tmax = (long)trans->datamax;
  1710.  if (tmin < 0) tmin = 0; else if (tmin > 255) tmin = 255;
  1711.  if (tmax < 0) tmax = 0; else if (tmax > 255) tmax = 255;
  1712.  
  1713.  cdata = (unsigned char *)buf;
  1714.  creplace = (unsigned char)trans->replacement;
  1715.  
  1716.  switch (hdulist->bitpix)
  1717.  {
  1718.    case 8:
  1719.      while (npix > 0)  /* For all pixels to read */
  1720.      {
  1721.        maxelem = sizeof (pixbuffer) / hdulist->bpp;
  1722.        if (maxelem > npix) maxelem = npix;
  1723.        if (fread ((char *)pixbuffer, hdulist->bpp, maxelem, ff->fp) != maxelem)
  1724.          return (-1);
  1725.        npix -= maxelem;
  1726.  
  1727.        pix = pixbuffer;
  1728.        if (hdulist->used.blank)
  1729.        {
  1730.          bp8blank = (FITS_BITPIX8)hdulist->blank;
  1731.          while (maxelem--)
  1732.          {
  1733.            bp8 = (FITS_BITPIX8)*(pix++);
  1734.            if (bp8 == bp8blank)      /* Is it a blank pixel ? */
  1735.              *(cdata++) = creplace;
  1736.            else                      /* Do transform */
  1737.            {
  1738.              tdata = (long)(bp8 * scale + offs);
  1739.              if (tdata < tmin) tdata = tmin;
  1740.              else if (tdata > tmax) tdata = tmax;
  1741.              *(cdata++) = (unsigned char)tdata;
  1742.            }
  1743.            transcount++;
  1744.          }
  1745.        }
  1746.        else
  1747.        {
  1748.          while (maxelem--)
  1749.          {
  1750.            bp8 = (FITS_BITPIX8)*(pix++);
  1751.            tdata = (long)(bp8 * scale + offs);
  1752.            if (tdata < tmin) tdata = tmin;
  1753.            else if (tdata > tmax) tdata = tmax;
  1754.            *(cdata++) = (unsigned char)tdata;
  1755.            transcount++;
  1756.          }
  1757.        }
  1758.      }
  1759.      break;
  1760.  
  1761.    case 16:
  1762.      while (npix > 0)  /* For all pixels to read */
  1763.      {
  1764.        maxelem = sizeof (pixbuffer) / hdulist->bpp;
  1765.        if (maxelem > npix) maxelem = npix;
  1766.        if (fread ((char *)pixbuffer, hdulist->bpp, maxelem, ff->fp) != maxelem)
  1767.          return (-1);
  1768.        npix -= maxelem;
  1769.  
  1770.        pix = pixbuffer;
  1771.        if (hdulist->used.blank)
  1772.        {
  1773.          bp16blank = (FITS_BITPIX16)hdulist->blank;
  1774.          while (maxelem--)
  1775.          {
  1776.            FITS_GETBITPIX16 (pix, bp16);
  1777.            if (bp16 == bp16blank)
  1778.              *(cdata++) = creplace;
  1779.            else
  1780.            {
  1781.              tdata = (long)(bp16 * scale + offs);
  1782.              if (tdata < tmin) tdata = tmin;
  1783.              else if (tdata > tmax) tdata = tmax;
  1784.              *(cdata++) = (unsigned char)tdata;
  1785.            }
  1786.            transcount++;
  1787.            pix += 2;
  1788.          }
  1789.        }
  1790.        else
  1791.        {
  1792.          while (maxelem--)
  1793.          {
  1794.            FITS_GETBITPIX16 (pix, bp16);
  1795.            tdata = (long)(bp16 * scale + offs);
  1796.            if (tdata < tmin) tdata = tmin;
  1797.            else if (tdata > tmax) tdata = tmax;
  1798.            *(cdata++) = (unsigned char)tdata;
  1799.            transcount++;
  1800.            pix += 2;
  1801.          }
  1802.        }
  1803.      }
  1804.      break;
  1805.  
  1806.    case 32:
  1807.      while (npix > 0)  /* For all pixels to read */
  1808.      {
  1809.        maxelem = sizeof (pixbuffer) / hdulist->bpp;
  1810.        if (maxelem > npix) maxelem = npix;
  1811.        if (fread ((char *)pixbuffer, hdulist->bpp, maxelem, ff->fp) != maxelem)
  1812.          return (-1);
  1813.        npix -= maxelem;
  1814.  
  1815.        pix = pixbuffer;
  1816.        if (hdulist->used.blank)
  1817.        {
  1818.          bp32blank = (FITS_BITPIX32)hdulist->blank;
  1819.          while (maxelem--)
  1820.          {
  1821.            FITS_GETBITPIX32 (pix, bp32);
  1822.            if (bp32 == bp32blank)
  1823.              *(cdata++) = creplace;
  1824.            else
  1825.            {
  1826.              tdata = (long)(bp32 * scale + offs);
  1827.              if (tdata < tmin) tdata = tmin;
  1828.              else if (tdata > tmax) tdata = tmax;
  1829.              *(cdata++) = (unsigned char)tdata;
  1830.            }
  1831.            transcount++;
  1832.            pix += 4;
  1833.          }
  1834.        }
  1835.        else
  1836.        {
  1837.          while (maxelem--)
  1838.          {
  1839.            FITS_GETBITPIX32 (pix, bp32);
  1840.            tdata = (long)(bp32 * scale + offs);
  1841.            if (tdata < tmin) tdata = tmin;
  1842.            else if (tdata > tmax) tdata = tmax;
  1843.            *(cdata++) = (unsigned char)tdata;
  1844.            transcount++;
  1845.            pix += 4;
  1846.          }
  1847.        }
  1848.      }
  1849.      break;
  1850.  
  1851.    case -32:
  1852.      while (npix > 0)  /* For all pixels to read */
  1853.      {
  1854.        maxelem = sizeof (pixbuffer) / hdulist->bpp;
  1855.        if (maxelem > npix) maxelem = npix;
  1856.        if (fread ((char *)pixbuffer, hdulist->bpp, maxelem, ff->fp) != maxelem)
  1857.          return (-1);
  1858.        npix -= maxelem;
  1859.  
  1860.        pix = pixbuffer;
  1861.        while (maxelem--)
  1862.        {
  1863.          if (fits_nan_32 (pix))    /* An IEEE special value ? */
  1864.            *(cdata++) = creplace;
  1865.          else                      /* Do transform */
  1866.          {
  1867.            FITS_GETBITPIXM32 (pix, bpm32);
  1868.            tdata = (long)(bpm32 * scale + offs);
  1869.            if (tdata < tmin) tdata = tmin;
  1870.            else if (tdata > tmax) tdata = tmax;
  1871.            *(cdata++) = (unsigned char)tdata;
  1872.          }
  1873.          transcount++;
  1874.          pix += 4;
  1875.        }
  1876.      }
  1877.      break;
  1878.  
  1879.    case -64:
  1880.      while (npix > 0)  /* For all pixels to read */
  1881.      {
  1882.        maxelem = sizeof (pixbuffer) / hdulist->bpp;
  1883.        if (maxelem > npix) maxelem = npix;
  1884.        if (fread ((char *)pixbuffer, hdulist->bpp, maxelem, ff->fp) != maxelem)
  1885.          return (-1);
  1886.        npix -= maxelem;
  1887.  
  1888.        pix = pixbuffer;
  1889.        while (maxelem--)
  1890.        {
  1891.          if (fits_nan_64 (pix))
  1892.            *(cdata++) = creplace;
  1893.          else
  1894.          {
  1895.            FITS_GETBITPIXM64 (pix, bpm64);
  1896.            tdata = (long)(bpm64 * scale + offs);
  1897.            if (tdata < tmin) tdata = tmin;
  1898.            else if (tdata > tmax) tdata = tmax;
  1899.            *(cdata++) = (unsigned char)tdata;
  1900.          }
  1901.          transcount++;
  1902.          pix += 8;
  1903.        }
  1904.      }
  1905.      break;
  1906.  }
  1907.  return (transcount);
  1908. }
  1909.  
  1910.  
  1911. #ifndef FITS_NO_DEMO
  1912. /*****************************************************************************/
  1913. /* #BEG-PAR                                                                  */
  1914. /*                                                                           */
  1915. /* Function  : fits_to_pgmraw - convert FITS-file to raw PGM-file            */
  1916. /*                                                                           */
  1917. /* Parameters:                                                               */
  1918. /* char *fitsfile          [I] : name of fitsfile                            */
  1919. /* char *pgmfile           [I] : name of pgmfile                             */
  1920. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  1921. /*                                                                           */
  1922. /* The function converts a FITS-file to a raw PGM-file. The PGM-file will    */
  1923. /* be upside down, because the orientation for storing the image is          */
  1924. /* different. On success, 0 is returned. On failure, -1 is returned.         */
  1925. /*                                                                           */
  1926. /* #END-PAR                                                                  */
  1927. /*****************************************************************************/
  1928.  
  1929. int fits_to_pgmraw (char *fitsfile, char *pgmfile)
  1930.  
  1931. {FITS_FILE *fitsin = NULL;
  1932.  FILE *pgmout = NULL;
  1933.  FITS_HDU_LIST *hdu;
  1934.  FITS_PIX_TRANSFORM trans;
  1935.  int retval = -1, nbytes, maxbytes;
  1936.  char buffer[1024];
  1937.  
  1938.  fitsin = fits_open (fitsfile, "r");  /* Open FITS-file for reading */
  1939.  if (fitsin == NULL) goto err_return;
  1940.  
  1941.  if (fitsin->n_pic < 1) goto err_return;  /* Any picture in it ? */
  1942.  
  1943.  hdu = fits_seek_image (fitsin, 1);       /* Position to the first image */
  1944.  if (hdu == NULL) goto err_return;
  1945.  if (hdu->naxis < 2) goto err_return;     /* Enough dimensions ? */
  1946.  
  1947.  pgmout = fopen (pgmfile, "wb");
  1948.  if (pgmout == NULL) goto err_return;
  1949.  
  1950.                                    /* Write PGM header with width/height */
  1951.  fprintf (pgmout, "P5\n%d %d\n255\n", hdu->naxisn[0], hdu->naxisn[1]);
  1952.  
  1953.  /* Set up transformation for FITS pixel values to 0...255 */
  1954.  /* It maps trans.pixmin to trans.datamin and trans.pixmax to trans.datamax. */
  1955.  /* Values out of range [datamin, datamax] are clamped */
  1956.  trans.pixmin = hdu->pixmin;
  1957.  trans.pixmax = hdu->pixmax;
  1958.  trans.datamin = 0.0;
  1959.  trans.datamax = 255.0;
  1960.  trans.replacement = 0.0;  /* Blank/NaN replacement value */
  1961.  trans.dsttyp = 'c';       /* Output type is character */
  1962.  
  1963.  nbytes = hdu->naxisn[0]*hdu->naxisn[1];
  1964.  while (nbytes > 0)
  1965.  {
  1966.    maxbytes = sizeof (buffer);
  1967.    if (maxbytes > nbytes) maxbytes = nbytes;
  1968.  
  1969.    /* Read pixels and transform them */
  1970.    if (fits_read_pixel (fitsin, hdu, maxbytes, &trans, buffer) != maxbytes)
  1971.      goto err_return;
  1972.  
  1973.    if (fwrite (buffer, 1, maxbytes, pgmout) != maxbytes)
  1974.      goto err_return;
  1975.  
  1976.    nbytes -= maxbytes;
  1977.  }
  1978.  retval = 0;
  1979.  
  1980. err_return:
  1981.  
  1982.  if (fitsin) fits_close (fitsin);
  1983.  if (pgmout) fclose (pgmout);
  1984.  
  1985.  return (retval);
  1986. }
  1987.  
  1988.  
  1989. /*****************************************************************************/
  1990. /* #BEG-PAR                                                                  */
  1991. /*                                                                           */
  1992. /* Function  : pgmraw to fits - convert raw PGM-file to FITS-file            */
  1993. /*                                                                           */
  1994. /* Parameters:                                                               */
  1995. /* char *pgmfile           [I] : name of pgmfile                             */
  1996. /* char *fitsfile          [I] : name of fitsfile                            */
  1997. /*                  ( mode : I=input, O=output, I/O=input/output )           */
  1998. /*                                                                           */
  1999. /* The function converts a raw PGM-file to a FITS-file. The FITS-file will   */
  2000. /* be upside down, because the orientation for storing the image is          */
  2001. /* different. On success, 0 is returned. On failure, -1 is returned.         */
  2002. /*                                                                           */
  2003. /* #END-PAR                                                                  */
  2004. /*****************************************************************************/
  2005.  
  2006. int pgmraw_to_fits (char *pgmfile, char *fitsfile)
  2007.  
  2008. {FITS_FILE *fitsout = NULL;
  2009.  FILE *pgmin = NULL;
  2010.  FITS_HDU_LIST *hdu;
  2011.  char buffer[1024];
  2012.  int width, height, numbytes, maxbytes;
  2013.  int retval = -1;
  2014.  
  2015.  fitsout = fits_open (fitsfile, "w");
  2016.  if (fitsout == NULL) goto err_return;
  2017.  
  2018.  pgmin = fopen (pgmfile, "r");
  2019.  if (pgmin == NULL) goto err_return;
  2020.  
  2021.  /* Read signature of PGM file */
  2022.  if (fgets (buffer, sizeof (buffer), pgmin) == NULL) goto err_return;
  2023.  if ((buffer[0] != 'P') || (buffer[1] != '5')) goto err_return;
  2024.  
  2025.  /* Skip comments upto width/height */
  2026.  do
  2027.  {
  2028.    if (fgets (buffer, sizeof (buffer), pgmin) == NULL) goto err_return;
  2029.  } while (buffer[0] == '#');
  2030.  
  2031.  if (sscanf (buffer, "%d%d", &width, &height) != 2) goto err_return;
  2032.  if ((width < 1) || (height < 1)) goto err_return;
  2033.  
  2034.  /* Skip comments upto maxval */
  2035.  do
  2036.  {
  2037.    if (fgets (buffer, sizeof (buffer), pgmin) == NULL) goto err_return;
  2038.  } while (buffer[0] == '#');
  2039.  /* Ignore maxval */
  2040.  
  2041.  hdu = fits_add_hdu (fitsout);     /* Create a HDU for the FITS file */
  2042.  if (hdu == NULL) goto err_return;
  2043.  
  2044.  hdu->used.simple = 1;         /* Set proper values */
  2045.  hdu->bitpix = 8;
  2046.  hdu->naxis = 2;
  2047.  hdu->naxisn[0] = width;
  2048.  hdu->naxisn[1] = height;
  2049.  hdu->used.datamin = 1;
  2050.  hdu->datamin = 0.0;
  2051.  hdu->used.datamax = 1;
  2052.  hdu->datamax = 255.0;
  2053.  hdu->used.bzero = 1;
  2054.  hdu->bzero = 0.0;
  2055.  hdu->used.bscale = 1;
  2056.  hdu->bscale = 1.0;
  2057.  
  2058.  fits_add_card (hdu, "");
  2059.  fits_add_card (hdu, "HISTORY THIS FITS FILE WAS GENERATED BY FITSRW\
  2060.  USING PGMRAW_TO_FITS"); 
  2061.  
  2062.  /* Write the header. Blocking is done automatically */
  2063.  if (fits_write_header (fitsout, hdu) < 0) goto err_return;
  2064.  
  2065.  /* The primary array plus blocking must be written manually */
  2066.  numbytes = width * height;
  2067.  
  2068.  while (numbytes > 0)
  2069.  {
  2070.    maxbytes = sizeof (buffer);
  2071.    if (maxbytes > numbytes) maxbytes = numbytes;
  2072.  
  2073.    if (fread (buffer, 1, maxbytes, pgmin) != maxbytes) goto err_return;
  2074.    if (fwrite (buffer, 1, maxbytes, fitsout->fp) != maxbytes) goto err_return;
  2075.  
  2076.    numbytes -= maxbytes;
  2077.  }
  2078.  
  2079.  /* Do blocking */
  2080.  numbytes = (width * height) % FITS_RECORD_SIZE;
  2081.  if (numbytes)
  2082.  {
  2083.    while (numbytes++ < FITS_RECORD_SIZE)
  2084.      if (putc (0, fitsout->fp) == EOF) goto err_return;
  2085.  }
  2086.  retval = 0;
  2087.  
  2088. err_return:
  2089.  
  2090.  if (fitsout) fits_close (fitsout);
  2091.  if (pgmin) fclose (pgmin);
  2092.  
  2093.  return (retval);
  2094. }
  2095.  
  2096. #endif
  2097.