/* wfits.c Copyright (c) Kapteyn Laboratorium Groningen 1993 All Rights Reserved. #> wfits.dc1 Program: WFITS Purpose: Writes disk set and subsets onto magnetic tape (or disk) in FITS format. Category: FITS, TAPES, UTILITY File: wfits.c Author: K.G. Begeman Keywords: INSET= Set (and subset(s)) where to copy data from. Maximum number of subsets is 2048. BOX= Frame of input subsets [entire subset]. FITSFILE= Name of file to write FITS data to [write to tape]. If more than one subset is to be written to FITS file, a generic filename should be entered. A generic filename contains fields of consecutive hash (#) symbols. The filenames of the FITS files will be generated by replacing each field with the subset sequence number. The number of hash symbols in a field denotes the field width. Examples: FITSFILE=FILE#.FTS generates FILE1.FTS, FILE2.FTS ... FITSFILE=FILE##.FTS generates FILE01.FTS, FILE02.FTS ... FITSFILE=FILE###d.FTS generates FILE001.FTS, FILE002.FTS ... If no hash symbols are present in a generic name, one hash symbol will be appended. If data is to be written to file, the following two keywords will be skipped. OUTTAPE= Name of output tape device [list of all tape devices]. SKIP= Number of files to skip on OUTTAPE [SKIP to End of Tape]. The first input file will be stored at this location, the second at the next location etc. Note that for new (empty) tapes you should give SKIP=0. BITPIX= Specify number of bits to represent pixel value [16]. This number can be 8, 16, 32, or -32 for IEEE floating point. For antenna patterns it is advised to take BITPIX=32 (or -32). ** BLOCKED= Specify the number of FITS records / tape record [10]. The maximum blocking factor is 20. ** MNMX= Always calculate min, max Y/[N]? If the minimum and maximum of the subsets are not found in the header, they are automatically calculated. If you are not sure that these header values are correct, give MNMX=Y. Description: FITS: A Flexible Image Transport System, is a format for the interchange of astronomical images and other digital arrays on magnetic tape. This format provides a simple but powerful mechanism for the unambiguous transmission of n-dimensional, regularly spaced data arrays. It also provides a method for the transmission of a virtually unlimited number of auxiliary parameters that may be associated with the image. The parameters are written in a form which is easily interpreted by both humans and computers (D. C. Wells, Astron. Astrophys. Suppl. Ser. 44, (1981) 363-370 ). WFITS creates a FITS header record from the GIPSY header items and converts the image data (INSET=) to 8, 16, 32 bits integers or to 32 bits IEEE floating point (BITPIX=). Programs to read FITS files like RFITS, use this header item to convert tape data to image data. Updates: Oct 21, 1991: KGB, Modified code for new mtopen routine. Nov 12, 1991: VOG, 'LONG_MIN' casted twice before assigning to float because of bug in 'gcc' Dec 9, 1991: VOG, Put OBJECT in fits header on users request. May 14, 1992: VOG, New defaults for MINMAX= Jul 20, 1993: KGB, Completely rewritten, replaced fts_ and ftsd_ routines. Apr 20, 1994: KGB, Minor repairs. Feb 7, 1996: KGB, Keyword FITSFILE= implemented. Okt 8, 1999: VOG, Removed non printable characters in record in gdsd_rfits. #< */ #include "ctype.h" /* */ #include "float.h" /* */ #include "limits.h" /* */ #include "math.h" /* */ #include "stdio.h" /* */ #include "stdlib.h" /* */ #include "string.h" /* */ #include "time.h" /* */ #include "gipsyc.h" /* GIPSY definitions */ #include "cmain.h" /* C program */ #include "cancel.h" /* cancel_c */ #include "finis.h" /* finis_c */ #include "ftsd_mkhead.h" /* ftsd_mkhead_c */ #include "gdsbox.h" /* gdsbox_c */ #include "gdsinp.h" /* gdsinp_c */ #include "gdsa_istable.h" /* gdsa_istable_c */ #include "gdsc_fill.h" /* gdsc_fill_c */ #include "gdsc_grid.h" /* gdsc_grid_c */ #include "gdsc_ndims.h" /* gdsc_ndims_c */ #include "gdsc_range.h" /* gdsc_range_c */ #include "gdsd_find.h" /* gdsd_find_c */ #include "gdsd_length.h" /* gdsd_length_c */ #include "gdsd_rfits.h" /* gdsd_rfits_c */ #include "gdsd_rvar.h" /* gdsd_rvar_c */ #include "gdsd_rewind.h" /* gdsd_rewind_c */ #include "gdsd_rreal.h" /* gdsd_rreal_c */ #include "gdsi_read.h" /* gdsi_read_c */ #include "init.h" /* init_c */ #include "minmax1.h" /* minmax_c */ #include "mtbsf.h" /* mtbsf_c */ #include "mtclose.h" /* mtclose_c */ #include "mtfsf.h" /* mrfsf_c */ #include "mtname.h" /* mtname_c */ #include "mtopen.h" /* mtopen_c */ #include "mtread.h" /* mtread_c */ #include "mtweof.h" /* mtweof_c */ #include "mtwrite.h" /* mtwrite_c */ #include "myname.h" /* myname_c */ #include "nelc.h" /* nelc_c */ #include "reject.h" /* reject_c */ #include "scaler.h" /* scaler_c */ #include "setfblank.h" /* setfblank_c */ #include "spfplf.h" /* spfplf_c */ #include "stabar.h" /* stabar_c */ #include "userfio.h" /* errorf etc. */ #define FITSBLOCKLEN 2880 /* Length FITS block */ #define FITSRECORDLEN 80 /* Length FITS header record */ #define GDSDESCRLEN 20 /* Length GDS descriptor name */ #define GDSVARLEN 132 /* Length GDS var. record */ #define MAXBLOCK 20 /* Max. blocking factor */ #define MAXDATA 4096 /* Buffer size for I/O */ #define MAXAXES 10 /* Max. number of axes */ #define MAXSUBSETS 2048 /* Max. number of subsets */ #define TEXTLEN 80 /* Length text strings */ #define VERSION "1.3" /* Version number */ #define NINT(r) ( r > 0.0 ? (int) ( r + 0.5 ) : (int) ( r - 0.5 ) ) #define MKCHARREC( name, value, comment ) \ { \ int k = 0; \ int l = sprintf( record.a, "%-8.8s= '%-8.68s'", name, value ); \ if (l < FITSRECORDLEN) record.a[l++] = ' '; \ if (l < FITSRECORDLEN) record.a[l++] = '/'; \ if (l < FITSRECORDLEN) record.a[l++] = ' '; \ while (l < FITSRECORDLEN && comment[k]) record.a[l++] = comment[k++]; \ while (l < FITSRECORDLEN) record.a[l++] = ' '; \ } #define MKDOUBLEREC( name, value, comment ) \ { \ int l = sprintf( record.a, "%-8.8s= %20.14E / %-47.47s", name, value, comment ); \ if ( l != FITSRECORDLEN ) { \ errorf( FATAL, "MKDOUBLEREC: Wrong record size (%d)", l ); \ } \ } #define MKFLOATREC( name, value, comment ) \ { \ int l = sprintf( record.a, "%-8.8s= %#20.6E / %-47.47s", name, value, comment ); \ if ( l != FITSRECORDLEN ) { \ errorf( FATAL, "MKFLOATREC: Wrong record size (%d)", l ); \ } \ } #define MKINTREC( name, value, comment ) \ { \ int l = sprintf( record.a, "%-8.8s= %20d / %-47.47s", name, value, comment ); \ if ( l != FITSRECORDLEN ) { \ errorf( FATAL, "MKINTREC: Wrong record size (%d)", l ); \ } \ } #define MKLOGREC( name, value, comment ) \ { \ int l; \ if (value) { \ l = sprintf( record.a, "%-8.8s= %20.1s / %-47.47s", name, "T", comment ); \ } else { \ l = sprintf( record.a, "%-8.8s= %20.1s / %-47.47s", name, "F", comment ); \ } \ if ( l != FITSRECORDLEN ) { \ errorf( FATAL, "MKLOGREC: Wrong record size (%d)", l ); \ } \ } #define PUTOUT( f, m, p, n ) \ { \ fint size = ( blocked * FITSBLOCKLEN ); \ fint plen = (n); \ byte *b = (byte *)p; \ if ((b == NULL) || (plen == 0)) { \ while (nfts % FITSBLOCKLEN) ftsdata[nfts++] = 0; \ if (nfts) { \ if (f != NULL) { \ if ( fwrite( ftsdata, sizeof( char ), nfts, f ) != nfts ) { \ errorf( FATAL, "Error writing to FITS file" ); \ } \ } else { \ fataltape( mtwrite_c( &m, (fint *) ftsdata, &nfts ) ); \ } \ nfts = 0; \ } \ } else { \ fint l = size - nfts; \ while ( plen > l ) { \ if (l) { \ memmove( &ftsdata[nfts], b, l ); \ nfts += l; \ plen -= l; \ b += l; \ } \ if (f != NULL) { \ if ( fwrite( ftsdata, sizeof( char ), nfts, f ) != nfts ) { \ errorf( FATAL, "Error writing to FITS file" ); \ } \ } else { \ fataltape( mtwrite_c( &m, (fint *) ftsdata, &nfts ) ); \ } \ nfts = 0; \ l = size - nfts; \ } \ if ( plen ) { \ memmove( &ftsdata[nfts], b, plen ); \ nfts += plen; \ } \ } \ } typedef unsigned char byte; typedef struct { char name[GDSDESCRLEN+1]; fint level; } desc; /* * Input of set, subsets: */ static char setb[TEXTLEN]; /* Buffer for name of set */ static fchar set = { setb, TEXTLEN }; /* Name of set */ static fint subsets[MAXSUBSETS]; /* The subsets */ static fint nsubs; /* Number of input subsets */ static fint subdim; /* Dimension of subset */ static fint setdim; /* Dimension of set */ static fint fgridlo[MAXAXES]; /* Grids of entire frame */ static fint fgridhi[MAXAXES]; static fint bgridlo[MAXAXES]; /* Grids of part of frame */ static fint bgridhi[MAXAXES]; static float blank; /* the blank value */ /* * Tape related: */ static byte ftsdata[MAXBLOCK*FITSBLOCKLEN]; /* The fits data buffer */ static char recordb[FITSRECORDLEN+1]; /* Fits header record buffer */ static fchar record = { recordb, FITSRECORDLEN }; /* Header record */ static fint nfts; /* Bytes in fits data buffer */ static char tnameb[TEXTLEN]; /* Buffer for Ftname */ static fchar tname = { tnameb, TEXTLEN }; /* Name of tape */ /* * Generate file name from a generic name. */ static int genname( char name[], char gname[], int n, int generate ) { char b[256]; int h = 0; int i = 0; int j = 0; if (!generate) { strcpy( name, gname ); return( strlen( name ) ); } while (gname[i]) { int k = 0; while (gname[i] == '#') { i++; k++; }; if (k) { int l = 0; h++; (void) sprintf( b, "%*.*d", k, k, n ); while (b[l]) name[j++] = b[l++]; } if (gname[i]) name[j++] = gname[i++]; } if (!h) { int l = 0; (void) sprintf( b, "%d", n ); while (b[l]) name[j++] = b[l++]; } name[j++] = 0; return( j ); } /* * Calculate minimum and maximum values in given box in given * subset. */ static void getmnmx( fint cwlo, fint cwhi, float mnmx[] ) { fint maxdata = MAXDATA, nread, tid = 0; float datamin, datamax; float data[MAXDATA]; mnmx[0] = mnmx[1] = blank; do { gdsi_read_c( set, &cwlo, &cwhi, data, &maxdata, &nread, &tid ); minmax1_c( data, &nread, &datamin, &datamax ); if (datamin != blank) { if (mnmx[0] == blank) { mnmx[0] = datamin; } else if (datamin < mnmx[0]) { mnmx[0] = datamin; } } if (datamax != blank) { if (mnmx[1] == blank) { mnmx[1] = datamax; } else if (datamax > mnmx[1]) { mnmx[1] = datamax; } } } while ( tid > 0 ); } /* * Tape routines return error code, use code to give error message. */ void fataltape( fint mterr ) { if (mterr >= 0) return; switch( mterr ) { case -1: { errorf( FATAL, "tape error: (-1) operation error" ); break; } case -2: { errorf( FATAL, "tape error: (-2) end of tape detected" ); break; } case -3: { errorf( FATAL, "tape error: (-3) device already opened" ); break; } case -4: { errorf( FATAL, "tape error: (-4) device not opened" ); break; } case -5: { errorf( FATAL, "tape error: (-5) tape mark encountered" ); break; } case -6: { errorf( FATAL, "tape error: (-6) not available" ); break; } case -7: { errorf( FATAL, "tape error: (-7) tape at BOT" ); break; } case -8: { errorf( FATAL, "tape error: (-8) record too long" ); break; } case -9: { errorf( FATAL, "tape error: (-9) error in call" ); break; } case -10: { errorf( FATAL, "tape error: (-10) write protected" ); break; } default: { errorf( FATAL, "tape error: (%d) Unknown error", mterr ); break; } } } /* * The main program. */ MAIN_PROGRAM_ENTRY { FILE *fid = NULL; /* id of FITS file */ bool domnmx; /* Calculate min and max */ char gname[FILENAME_MAX+1]; /* generic fits file name */ fint bitpix; /* Bits/pixel */ fint blocked; /* Blocking factor */ fint mtid; /* Tape transfer id */ float stabarv[3]; /* Values for stabar */ int isub; /* Subset counter */ /* * Initialize some things. */ init_c(); /* Contact Hermes */ setfblank_c( &blank ); /* Get system blank */ if (CHAR_BIT != 8) { errorf( FATAL, "Program expects 8 bit bytes!" ); } if (sizeof( char ) != 1) { errorf( FATAL, "Program expects 1 byte chars!" ); } if (sizeof( short ) != 2) { errorf( FATAL, "Program expects 2 byte shorts!" ); } if (sizeof( int ) != 4) { errorf( FATAL, "Program expects 4 byte ints!" ); } if (sizeof( float ) != 4) { errorf( FATAL, "Program expects 4 byte floats!" ); } /* * Task identification */ { char b[20]; /* Buffer for task name */ fchar task; /* Name of current task */ task.a = b; task.l = sizeof(b) - 1; /* Initialize */ myname_c( task ); /* Get task name */ task.a[nelc_c(task)] = '\0'; /* Terminate name with null */ IDENTIFICATION( task.a, VERSION ); /* Show task and version */ } /* * Get the set and subsets to write to tape. */ { fchar keyword; fchar message; fint axcount[MAXAXES]; fint axnum[MAXAXES]; fint class; fint cwhi; fint cwlo; fint input; fint error; fint maxaxes; fint maxsubs; fint output; fint wholeset; int m; keyword = tofchar("INSET="); message = tofchar("Give set (, subsets) to transfer to tape" ); class = 1; input = 0; maxaxes = MAXAXES; maxsubs = MAXSUBSETS; wholeset = 0; output = 3; subdim = 0; nsubs = gdsinp_c( set, subsets, &maxsubs, &input, keyword, message, &output, axnum, axcount, &maxaxes, &class, &subdim ); wholeset= 0; setdim = gdsc_ndims_c( set, &wholeset ); error = 0; gdsc_range_c( set, &wholeset, &cwlo, &cwhi, &error ); for (m = 0; m < setdim; m++) { fgridlo[m] = gdsc_grid_c( set, &axnum[m], &cwlo, &error ); fgridhi[m] = gdsc_grid_c( set, &axnum[m], &cwhi, &error ); } } /* * Prepare a box for INSET. Default is entire subset */ { fchar keyword; fchar message; fint boxopt; fint input; fint output; int m; int totpixels; keyword = tofchar("BOX="); message = tofchar(" "); boxopt = 0; input = 1; output = 3; gdsbox_c( bgridlo, bgridhi, set, subsets, &input, keyword, message, &output, &boxopt ); totpixels = 1; for ( m = 0; m < subdim; m++) { totpixels *= (bgridhi[m] - bgridlo[m] + 1); } totpixels *= nsubs; stabarv[0] = 0.0; stabarv[1] = (float) totpixels; stabarv[2] = 0.0; } /* * Get name of FITS file to write. */ { fchar Name; memset( gname, ' ', sizeof( gname ) ); Name.a = gname; Name.l = sizeof( gname ) - 1; (void) userfchar( Name, 1, 1, "FITSFILE=", "%sname for output FITS file%s[write to tape]", (nsubs > 1) ? "generic " : " ", (nsubs > 1) ? "s " : " " ); gname[nelc_c(Name)] = 0; } /* * Open a tape device and position it. */ if (!gname[0]) { fchar keymes; fint count; fint mterr; fint ndone; fint nskip; keymes = tofchar( "?OUTTAPE=Name of output tape device [list of all tape devices]" ); mtid = mtopen_c( keymes ); fataltape( mtid ); mterr = mtname_c( &mtid, tname ); fataltape( mterr ); count = 0; nskip = 0; if (!userfint( &nskip, 1, 1, "SKIP=", "Number of files to skip on tape [skip to EOT]" )) { fint data; fint ndat; nskip = 1; mterr = mtbsf_c( &mtid, &nskip ); fataltape( mterr ); if (mterr == 0) { anyoutf( 3, "Tape at BOT" ); } else { count -= 1; } ndat = 1; ndone = 1; while ( ( (mterr = mtfsf_c( &mtid, &nskip ) ) == 1 ) && ( ( ndone = mtread_c( &mtid, &data, &ndat ) ) > 0 ) ) { count++; } if ( ndone == -2 ) ndone = 0; fataltape( mterr ); fataltape( ndone ); if (( mterr == 1 ) && ( ndone == 0 )) { anyoutf( 3, "Skipped %d files forward", ++count ); fataltape( mtbsf_c( &mtid, &nskip )); } } else if (nskip > 0) { mterr = mtfsf_c( &mtid, &nskip ); fataltape( mterr ); anyoutf( 3, "Skipped %d files forward", mterr ); } else if (nskip < 0) { count = nskip = 1 - nskip; mterr = mtbsf_c( &mtid, &nskip ); fataltape( mterr ); if (mterr == nskip) { nskip = 1; mterr = mtfsf_c( &mtid, &nskip ); fataltape( mterr ); anyoutf( 3, "Skipped %d files backward", --count ); } else if (mterr >= 0) { anyoutf( 3, "Tape at BOT" ); anyoutf( 3, "Skipped %d files backward", mterr ); } } } /* * Get BITPIX= and BLOCKING factor. */ { int agreed; bitpix = 16; do { (void) userfint( &bitpix, 1, 1, "BITPIX=", "Specify number of bits to represent pix. val. [16]" ); agreed = ( (bitpix == 8) || (bitpix == 16) || (bitpix == 32) || (bitpix == -32) ); if (!agreed) { reject_c( tofchar( "BITPIX=" ), tofchar("Only 8,16,32 or -32") ); bitpix = 16; } } while (!agreed); blocked = 10; do { (void) userfint( &blocked, 1, 2, "BLOCKED=", "Specify number FITS BLOCKS/ tape block (1..20) [10]" ); agreed = ( blocked > 0 && blocked <= MAXBLOCK ); if (!agreed) { reject_c( tofchar( "BLOCKED=" ), tofchar("Blocking factor out of range") ); blocked = 10; } } while (!agreed); } /* * Do we calculate the min and the max? */ { domnmx = FALSE; (void) userflog( &domnmx, 1, 2, "MNMX=", "Always calculate min, max? Y/[N]" ); } /* * Loop over all subsets, create header, read data and write to tape. */ for ( isub = 0; isub < nsubs; isub++ ) { fint cwlo; fint cwhi; fint fitsblank; float mnmx[2]; double bscale; double bzero; if (gname[0]) { char name[FILENAME_MAX+1]; (void) genname( name, gname, isub+1, nsubs > 1 ); fid = NULL; do { fid = fopen( name, "r+b" ); if (fid != NULL) { bool ok = TRUE; (void) userflog( &ok, 1, 1, "OKAY=", "File \"%s\" exists, overwrite [Y]/N", name ); cancel( "OKAY=" ); if (!tobool( ok )) { (void) fclose( fid ); fid = NULL; } } else { fid = fopen( name, "wb" ); } if (fid == NULL) { fchar Name; Name.a = name; Name.l = sizeof( name ) - 1; (void) userfchar( Name, 1, 0, "FILENAME=", "New name of file" ); cancel( "FILENAME=" ); name[nelc_c(Name)] = 0; } } while (fid == NULL); anyoutf( 0, "Writing to FITS file \"%s\"", name ); } nfts = 0; /* reset ftsdata pointer */ /* * Items SIMPLE and BITPIX are created first */ { MKLOGREC( "SIMPLE", 1, "SIMPLE FITS FORMAT" ); PUTOUT( fid, mtid, record.a, FITSRECORDLEN ); MKINTREC( "BITPIX", bitpix, "NUMBER OF BITS PER PIXEL" ); PUTOUT( fid, mtid, record.a, FITSRECORDLEN ); } /* * Create coordinate words for use in 'mkhead'. The routine 'mkhead' * creates all axes related keywords and the keyword 'INSTRUME' */ { bool blckd; fchar head; fint maxrec; fint ret; if (blocked != 1) blckd = TRUE; else blckd = FALSE; cwlo = gdsc_fill_c( set, &subsets[isub], bgridlo ); cwhi = gdsc_fill_c( set, &subsets[isub], bgridhi ); head.a = (char *) &ftsdata[nfts]; head.l = sizeof( ftsdata ) - nfts; maxrec = head.l / FITSRECORDLEN; ret = ftsd_mkhead_c( set, &cwlo, &cwhi, &blckd, head, &maxrec ); switch( ret ) { case -1: { errorf( FATAL, "Set does not exist" ); break; } case -2: { errorf( FATAL, "Memory allocation error" ); break; } case -3: { errorf( FATAL, "Not enough space in header" ); break; } default: { break; } } PUTOUT( fid, mtid, head.a, ret * FITSRECORDLEN ); } /* * Determine min. and max. image value and calculate BSCALE and BZERO * ( image = tape*BSCALE + BZERO ). Write these parameters to the FITS * header. */ { double bitmax; double bitmin; fint ret1; fint ret2; ret1 = ret2 = 0; gdsd_rreal_c( set, tofchar("DATAMIN"), &subsets[isub], &mnmx[0], &ret1 ); gdsd_rreal_c( set, tofchar("DATAMAX"), &subsets[isub], &mnmx[1], &ret2 ); if ( (ret1 != subsets[isub]) || (ret2 != subsets[isub]) || tobool(domnmx) ) { getmnmx( cwlo, cwhi, mnmx ); } if (mnmx[0] == blank) { anyoutf( 1, "No defined values" ); mnmx[0] = -1.0; mnmx[1] = 1.0; } else if (mnmx[0] == mnmx[1]) { anyoutf( 1, "Values are equal" ); } else { anyoutf( 1, "Minimum and maximum in box : [%g, %g]", mnmx[0], mnmx[1] ); } switch( bitpix ) { case 8: { fitsblank = 0; bitmin = 2.0; bitmax = (float)(int)UCHAR_MAX; break; } case 16: { fitsblank = SHRT_MIN; bitmin = (float)(int)SHRT_MIN + 16.0; bitmax = (float)(int)SHRT_MAX - 16.0; break; } case 32: { fitsblank = INT_MIN; bitmin = (float)(int)INT_MIN + 4096.0; bitmax = (float)(int)INT_MAX - 4096.0; break; } case -32: { fint type = 0; fint items = 1; spfplf_c( &type, &blank, (float *)&fitsblank, &items ); bitmin = mnmx[0]; bitmax = mnmx[1]; break; } default: { bitmax = bitmin = 0.0; /* ha ha */ break; } } if ( bitpix > 0 ) { bscale = ( mnmx[0] - mnmx[1] ) / ( bitmin - bitmax ); if (bscale == 0.0) { bscale = mnmx[0] / bitmax; bzero = 0.0; } else { bzero = ( mnmx[1] * bitmin - mnmx[0] * bitmax ) / ( bitmin - bitmax ) ; } } else { bscale = 1.0; bzero = 0.0; } MKINTREC( "BLANK", fitsblank, "TAPE VALUE FOR UNDEFINED PIXELS" ); PUTOUT( fid, mtid, record.a, FITSRECORDLEN ); MKFLOATREC( "BSCALE", bscale, "TRUE = TAPE * BSCALE + BZERO" ); PUTOUT( fid, mtid, record.a, FITSRECORDLEN ); MKFLOATREC( "BZERO", bzero, "BIAS FOR THE REAL DATA" ); PUTOUT( fid, mtid, record.a, FITSRECORDLEN ); MKFLOATREC( "DATAMAX", mnmx[1], "MAXIMUM DATA VALUE" ); PUTOUT( fid, mtid, record.a, FITSRECORDLEN ); MKFLOATREC( "DATAMIN", mnmx[0], "MINIMUM DATA VALUE" ); PUTOUT( fid, mtid, record.a, FITSRECORDLEN ); } /* * Our ID. */ { char b[FITSRECORDLEN-12+1]; sprintf( b, "WFITS VERSION %s", VERSION ); MKCHARREC( "ORIGIN", b, "VERSION OF THE GIPSY PROGRAM" ); PUTOUT( fid, mtid, record.a, FITSRECORDLEN ); } /* * Now get descriptor items. */ { char descrb[GDSDESCRLEN]; desc *dsc = NULL; fchar descr; fint recordnum; int first_table = 1; int m, n; int ndsc = 0; recordnum = 0; /* Position in header file */ descr.a = descrb; descr.l = sizeof( descrb ); do { fint level; fint ret; int table_status; /* * The main purpose of the following is to generate all available * keywords with the function 'gdsd_find'. The HISTORY and * COMMENT fields are read separately as variable length records. */ level = subsets[isub]; ret = 0; gdsd_find_c( descr, set, &level, &recordnum, &ret ); if ((recordnum != 0) && (ret < 0)) { errorf( WARNING, "Header error: Found entry but no item" ); } if ((recordnum != 0) && (ret >= 0)) { table_status = gdsa_istable_c( descr ); if ( (table_status) && (first_table) ) { first_table = 0; anyoutf( 1, "WFITS: Descriptor contains tables which are NOT included in FITS header" ); } if ( !table_status && strncmp( descr.a, "BITPIX", 6 ) && strncmp( descr.a, "BLANK", 5 ) && strncmp( descr.a, "BLOCKED", 7 ) && strncmp( descr.a, "BSCALE", 6 ) && strncmp( descr.a, "BZERO", 5 ) && strncmp( descr.a, "CDELT", 5 ) && strncmp( descr.a, "CRPIX", 5 ) && strncmp( descr.a, "CROTA", 5 ) && strncmp( descr.a, "CTYPE", 5 ) && strncmp( descr.a, "CUNIT", 5 ) && strncmp( descr.a, "CRVAL", 5 ) && strncmp( descr.a, "DATAMAX", 7 ) && strncmp( descr.a, "DATAMIN", 7 ) && strncmp( descr.a, "DDELT", 5 ) && strncmp( descr.a, "DRPIX", 5 ) && strncmp( descr.a, "DROTA", 5 ) && strncmp( descr.a, "DTYPE", 5 ) && strncmp( descr.a, "DUNIT", 5 ) && strncmp( descr.a, "DRVAL", 5 ) && strncmp( descr.a, "EPOCH", 5 ) && strncmp( descr.a, "FREQ0", 5 ) && strncmp( descr.a, "INSTRUME", 8 ) && strncmp( descr.a, "NAXIS", 5 ) && strncmp( descr.a, "NBLANK", 6 ) && strncmp( descr.a, "ORIGIN", 6 ) && strncmp( descr.a, "SIMPLE", 6) ) { for ( n = 0; n < ndsc && strncmp( descr.a, dsc[n].name, GDSDESCRLEN ); n++); if (n == ndsc) { dsc = realloc( dsc, sizeof( desc ) * (++ndsc) ); if (dsc == NULL) { errorf( FATAL, "Memory allocation problems" ); } memmove( dsc[n].name, descr.a, GDSDESCRLEN ); dsc[n].level = ret; } else if (ret == subsets[isub]) { dsc[n].level = ret; } } } } while (recordnum != 0); for ( n = 0; n < ndsc; n++) { fint ret = 0; if (strncmp( dsc[n].name, "HISTORY", 7 ) && strncmp( dsc[n].name, "COMMENT", 7 )) { memmove( descr.a, dsc[n].name, GDSDESCRLEN ); gdsd_rfits_c( set, descr, &subsets[isub], record, &ret ); /* Remove zero's from record. */ { int i, l1, l2, len; l1 = nelc_c(record); l2 = strlen(record.a); len = l1; if (l2 < len) len = l2; for (i = 0; i < len; i++) { if (!isprint(record.a[i])) record.a[i] = ' ' ; } for (i = len; i < FITSRECORDLEN; i++) record.a[i] = ' '; } if (ret >= 0) { PUTOUT( fid, mtid, record.a, FITSRECORDLEN ); } } } /* * History and comments on selected or higher level. * There are separate loops for both HISTORY and COMMENT. * For each item we have to reset the current read position at * the beginning of the descriptor item by applying 'gdsd_rewind'. */ for (m = 0; m < 2; m++ ) { fint l, n; if (m == 0) descr = tofchar( "COMMENT" ); if (m == 1) descr = tofchar( "HISTORY" ); n = 0; for (n = 0; n < ndsc && strncmp( dsc[n].name, descr.a, descr.l ); n++); if (n < ndsc) { fint level; fint ret = 0; level = dsc[n].level; l = gdsd_length_c( set, descr, &level, &ret ); if (ret >= 0 && l > 0) { gdsd_rewind_c( set, descr, &level, &ret ); while ( ret >= 0 ) { char varrecb[GDSVARLEN]; fchar varrec; varrec.a = varrecb; varrec.l = sizeof( varrecb ); ret = 0; gdsd_rvar_c( set, descr, &level, varrec, &ret ); if (ret >= 0) { int j = 0; int l = nelc_c( varrec ); while (j < l) { int k = 0; memmove( record.a, descr.a, descr.l ); k = k + descr.l; record.a[k++] = ' '; while (k mnmx[1]) { d[n] = fitsblank; nblank++; } else { double t = ( datum - bzero ) / bscale; d[n] = NINT( t ); } } PUTOUT( fid, mtid, d, nread ); break; } case 16: { short d[MAXDATA]; int n; for ( n = 0; n < nread; n++ ) { datum = image[n]; if (datum == blank) { d[n] = fitsblank; } else if (datum < mnmx[0]) { d[n] = fitsblank; nblank++; } else if (datum > mnmx[1]) { d[n] = fitsblank; nblank++; } else { double t = ( datum - bzero ) / bscale; d[n] = NINT( t ); } } #if ( OS_INTEGER_TYPE == 1 ) { union { short s; byte b[sizeof(short)]; } ui, uo; for ( n = 0; n < nread; n++ ) { ui.s = d[n]; uo.b[0] = ui.b[1]; uo.b[1] = ui.b[0]; d[n] = uo.s; } } #endif PUTOUT( fid, mtid, d, nread * 2 ); break; } case 32: { int d[MAXDATA]; int n; for ( n = 0; n < nread; n++ ) { datum = image[n]; if (datum == blank) { d[n] = fitsblank; } else if (datum < mnmx[0]) { d[n] = fitsblank; nblank++; } else if (datum > mnmx[1]) { d[n] = fitsblank; nblank++; } else { double t = ( datum - bzero ) / bscale; d[n] = NINT( t ); } } #if ( OS_INTEGER_TYPE == 1 ) { union { int l; byte b[sizeof(int)]; } ui, uo; for ( n = 0; n < nread; n++ ) { ui.l = d[n]; uo.b[0] = ui.b[3]; uo.b[1] = ui.b[2]; uo.b[2] = ui.b[1]; uo.b[3] = ui.b[0]; d[n] = uo.l; } } #endif PUTOUT( fid, mtid, d, nread * 4 ); break; } case -32: { fint type = 0; spfplf_c( &type, image, image, &nread ); PUTOUT( fid, mtid, image, nread * 4 ); break; } default: { break; } } stabarv[2] += (float) nread; (void) stabar_c( &stabarv[0], &stabarv[1], &stabarv[2] ); } while (tid != 0); if (nblank) { anyoutf( 1, "Introduced %d Undefined value%s", nblank, ( ( nblank == 1 ) ? "!" : "s!") ); } } /* * Fix up last FITS block. */ { fint tm = 1; if (nfts) { PUTOUT( fid, mtid, NULL, 0 ); } if (fid == NULL) { fataltape( mtweof_c( &mtid, &tm ) ); } else { (void) fclose( fid ); } } } /* end of subset loop */ /* * Write final Tape Mark. */ if (fid == NULL) { fint tm = 1; fint skip = 1; fataltape( mtweof_c( &mtid, &tm ) ); /* write final tape mark */ fataltape( mtbsf_c( &mtid, &skip ) ); /* skip back a file */ fataltape( mtclose_c( &mtid ) ); /* close tape device */ } finis_c( ); /* last routine */ return( EXIT_SUCCESS ); /* last statement */ }