cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
punch_fits.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 #include "cddefines.h"
4 #include "cddrive.h"
5 #include "optimize.h"
6 #include "grid.h"
7 #include "punch.h"
8 #include "rfield.h"
9 #include "prt.h"
10 #include "input.h"
11 #include "version.h"
12 #include "physconst.h"
13 
14 #define RECORDSIZE 2880
15 #define LINESIZE 80
16 
17 #if defined(_BIG_ENDIAN)
18  /* the value of A will not be manipulated */
19 # define HtoNL(A) (A)
20 /*
21 # define HtoNS(A) (A)
22 # define NtoHS(A) (A)
23 # define NtoHL(A) (A)
24 */
25 #else
26 /* defined(_LITTLE_ENDIAN) */
27 /* the value of A will be byte swapped */
28 # define HtoNL(A) ((((A) & 0xff000000) >> 24) | \
29  (((A) & 0x00ff0000) >> 8) | \
30  (((A) & 0x0000ff00) << 8) | \
31  (((A) & 0x000000ff) << 24))
32 /*
33 # define HtoNS(A) ((((A) & 0xff00) >> 8) | (((A) & 0x00ff) << 8))
34 # define NtoHS HtoNS
35 # define NtoHL HtoNL
36 */
37 /*#else
38 error One of BIG_ENDIAN or LITTLE_ENDIAN must be #defined.*/
39 #endif
40 
41 #define ByteSwap5(x) ByteSwap((unsigned char *) &x,sizeof(x))
42 
43 #if !defined(_BIG_ENDIAN)
44 STATIC void ByteSwap(unsigned char * b, int n)
45 {
46  register int i = 0;
47  register int j = n-1;
48  while (i<j)
49  {
50  char temp = b[i];
51  b[i] = b[j];
52  b[j] = temp;
53  /* std::swap(b[i], b[j]); */
54  i++, j--;
55  }
56  return;
57 }
58 #endif
59 
60 static FILE *ioFITS_OUTPUT;
61 static long bytesAdded = 0;
62 static long bitpix = 8;
63 static long pcount = 0;
64 static long gcount = 1;
65 static long maxParamValues = 0;
66 const char ModelUnits[2][17] = {"'dimensionless '", "'photons/cm^2/s'" };
67 
68 STATIC void punchFITS_PrimaryHeader( bool lgAddModel );
69 STATIC void punchFITS_ParamHeader( /* long *numParamValues, */ long nintparm, long naddparm );
70 STATIC void punchFITS_ParamData( char **paramNames, long *paramMethods, realnum **paramRange,
71  realnum **paramData, long nintparm, long naddparm,
72  long *numParamValues );
73 STATIC void punchFITS_EnergyHeader( long numEnergies );
74 STATIC void punchFITS_EnergyData( realnum *Energies, long numEnergies );
75 STATIC void punchFITS_SpectraHeader( bool lgAdditiveModel, long nintparm, long naddparm, long totNumModels, long numEnergies );
76 STATIC void punchFITS_SpectraData( realnum **interpParameters, realnum **theSpectrum,
77  long totNumModels, long numEnergies, long nintparm, long naddparm );
78 STATIC void punchFITS_GenericHeader( long numEnergies );
79 STATIC void punchFITS_GenericData( long numEnergies );
80 STATIC void writeCloudyDetails( void );
81 STATIC long addComment( const char *CommentToAdd );
82 STATIC long addKeyword_txt( const char *theKeyword, const void *theValue, const char *theComment, long Str_Or_Log );
83 STATIC long addKeyword_num( const char *theKeyword, long theValue, const char *theComment);
84 
85 void punchFITSfile( FILE* ioPUN, int option )
86 {
87 
88  long i;
89 
90  DEBUG_ENTRY( "punchFITSfile()" );
91 
92  if( !grid.lgGridDone && option != NUM_OUTPUT_TYPES )
93  {
94  /* don't do any of this until flag set true. */
95  return;
96  }
97 
98  ioFITS_OUTPUT = ioPUN;
99 
100 #if 0
101  {
102 
103  long i,j;
104  FILE *asciiDump;
105 
106  if( (asciiDump = fopen( "gridspectra.con","w" ) ) == NULL )
107  {
108  fprintf( ioQQQ, "punchFITSfile could not open punch file for writing.\nSorry.\n" );
109  cdEXIT(EXIT_FAILURE);
110  }
111 
112  for( i=0; i<grid.numEnergies-1; i++ )
113  {
114  fprintf( asciiDump, "%.12e\t", grid.Energies[i] );
115  for( j=0; j<grid.totNumModels; j++ )
116  {
117  fprintf( asciiDump, "%.12e\t", grid.Spectra[10][j][i] );
118  }
119  fprintf( asciiDump, "\n" );
120  }
121  fclose( asciiDump );
122  }
123 #endif
124 
125  /* This is generic FITS option */
126  if( option == NUM_OUTPUT_TYPES )
127  {
128  punchFITS_PrimaryHeader( false );
131  }
132  /* These are specially designed XSPEC outputs. */
133  else if( option >= 1 && option < NUM_OUTPUT_TYPES )
134  {
135  bool lgAdditiveModel;
136 
137  /* option 10 is exp(-tau). */
138  if( option == 10 )
139  {
140  /* false says not an additive model */
141  lgAdditiveModel = false;
142  }
143  else
144  {
145  lgAdditiveModel = true;
146  }
147 
148  punchFITS_PrimaryHeader( lgAdditiveModel );
149 
150  for( i=0; i<grid.nintparm+grid.naddparm; i++ )
151  {
153  }
154 
155  ASSERT( maxParamValues >= 2 );
156 
157  punchFITS_ParamHeader( /* grid.numParamValues, */ grid.nintparm, grid.naddparm );
165  }
166  else
167  {
168  fprintf( ioQQQ, "PROBLEM - undefined option encountered in punchFITSfile. \n" );
169  cdEXIT( EXIT_FAILURE );
170  }
171  return;
172 }
173 
174 STATIC void punchFITS_PrimaryHeader( bool lgAddModel )
175 {
176  const char *ModelName = "'CLOUDY'";
177 
178  DEBUG_ENTRY( "punchFITS_PrimaryHeader()" );
179 
180  bytesAdded = 0;
181 
182  fixit(); /* bitpix is wrong when realnum is double? */
183 
184  bytesAdded += addKeyword_txt( "SIMPLE" , "T", "file does conform to FITS standard", 1 );
185  bytesAdded += addKeyword_num( "BITPIX" , bitpix, "number of bits per data pixel" );
186  bytesAdded += addKeyword_num( "NAXIS" , 0, "number of data axes" );
187  bytesAdded += addKeyword_txt( "EXTEND" , "T", "FITS dataset may contain extensions", 1 );
188  bytesAdded += addKeyword_txt( "CONTENT" , "'MODEL '", "spectrum file contains time intervals and event", 0 );
189  bytesAdded += addKeyword_txt( "MODLNAME", ModelName, "Model name", 0 );
190  bytesAdded += addKeyword_txt( "MODLUNIT", ModelUnits[lgAddModel], "Model units", 0 );
191  bytesAdded += addKeyword_txt( "REDSHIFT", "T", "If true then redshift will be included as a par", 1 );
192  if( lgAddModel == true )
193  {
194  bytesAdded += addKeyword_txt( "ADDMODEL", "T", "If true then this is an additive table model", 1 );
195  }
196  else
197  {
198  bytesAdded += addKeyword_txt( "ADDMODEL", "F", "If true then this is an additive table model", 1 );
199  }
200 
201  /* bytes are added here as well */
203 
204  bytesAdded += addKeyword_txt( "HDUCLASS", "'OGIP '", "Format conforms to OGIP/GSFC conventions", 0 );
205  bytesAdded += addKeyword_txt( "HDUCLAS1", "'XSPEC TABLE MODEL'","Extension contains an image", 0 );
206  bytesAdded += addKeyword_txt( "HDUVERS" , "'1.0.0 '", "Version of format (OGIP memo OGIP-92-001)", 0 );
207  /* After everything else */
208  bytesAdded += fprintf(ioFITS_OUTPUT, "%-80s", "END" );
209 
210  ASSERT( bytesAdded%LINESIZE == 0 );
211 
212  /* Now add blanks */
213  while( bytesAdded%RECORDSIZE > 0 )
214  {
215  bytesAdded += fprintf(ioFITS_OUTPUT, "%-1s", " " );
216  }
217  return;
218 }
219 
220 STATIC void punchFITS_ParamHeader( /* long *numParamValues, */ long nintparm, long naddparm )
221 {
222  long numFields = 10;
223  long naxis, naxis1, naxis2;
224  char theValue[20];
225  char theValue_temp[20];
226 
227  DEBUG_ENTRY( "punchFITS_ParamHeader()" );
228 
229  ASSERT( nintparm+naddparm <= LIMPAR );
230 
231  /* Make sure the previous blocks are the right size */
232  ASSERT( bytesAdded%RECORDSIZE == 0 );
233 
234  naxis = 2;
235  /* >>chng 06 aug 23, change to maximum number of parameter values. */
236  naxis1 = 44+4*maxParamValues;
237  naxis2 = nintparm+naddparm;
238 
239  bytesAdded += addKeyword_txt( "XTENSION", "'BINTABLE'", "binary table extension", 0 );
240  bytesAdded += addKeyword_num( "BITPIX" , bitpix, "8-bit bytes" );
241  bytesAdded += addKeyword_num( "NAXIS" , naxis, "2-dimensional binary table" );
242  bytesAdded += addKeyword_num( "NAXIS1" , naxis1, "width of table in bytes" );
243  bytesAdded += addKeyword_num( "NAXIS2" , naxis2, "number of rows in table" );
244  bytesAdded += addKeyword_num( "PCOUNT" , pcount, "size of special data area" );
245  bytesAdded += addKeyword_num( "GCOUNT" , gcount, "one data group (required keyword)" );
246  bytesAdded += addKeyword_num( "TFIELDS" , numFields, "number of fields in each row" );
247  bytesAdded += addKeyword_txt( "TTYPE1" , "'NAME '", "label for field 1", 0 );
248  bytesAdded += addKeyword_txt( "TFORM1" , "'12A '", "data format of the field: ASCII Character", 0 );
249  bytesAdded += addKeyword_txt( "TTYPE2" , "'METHOD '", "label for field 2", 0 );
250  bytesAdded += addKeyword_txt( "TFORM2" , "'J '", "data format of the field: 4-byte INTEGER", 0 );
251  bytesAdded += addKeyword_txt( "TTYPE3" , "'INITIAL '", "label for field 3", 0 );
252  bytesAdded += addKeyword_txt( "TFORM3" , "'E '", "data format of the field: 4-byte REAL", 0 );
253  bytesAdded += addKeyword_txt( "TTYPE4" , "'DELTA '", "label for field 4", 0 );
254  bytesAdded += addKeyword_txt( "TFORM4" , "'E '", "data format of the field: 4-byte REAL", 0 );
255  bytesAdded += addKeyword_txt( "TTYPE5" , "'MINIMUM '", "label for field 5", 0 );
256  bytesAdded += addKeyword_txt( "TFORM5" , "'E '", "data format of the field: 4-byte REAL", 0 );
257  bytesAdded += addKeyword_txt( "TTYPE6" , "'BOTTOM '", "label for field 6", 0 );
258  bytesAdded += addKeyword_txt( "TFORM6" , "'E '", "data format of the field: 4-byte REAL", 0 );
259  bytesAdded += addKeyword_txt( "TTYPE7" , "'TOP '", "label for field 7", 0 );
260  bytesAdded += addKeyword_txt( "TFORM7" , "'E '", "data format of the field: 4-byte REAL", 0 );
261  bytesAdded += addKeyword_txt( "TTYPE8" , "'MAXIMUM '", "label for field 8", 0 );
262  bytesAdded += addKeyword_txt( "TFORM8" , "'E '", "data format of the field: 4-byte REAL", 0 );
263  bytesAdded += addKeyword_txt( "TTYPE9" , "'NUMBVALS'", "label for field 9", 0 );
264  bytesAdded += addKeyword_txt( "TFORM9" , "'J '", "data format of the field: 4-byte INTEGER", 0 );
265  bytesAdded += addKeyword_txt( "TTYPE10" , "'VALUE '", "label for field 10", 0 );
266 
267  /* >>chng 06 aug 23, use maxParamValues instead of numParamValues */
268  /* The size of this array is dynamic, set to size of the maximum of the numParamValues vector */
269  sprintf( theValue_temp, "%ld%s", maxParamValues, "E" );
270  sprintf( theValue, "%s%-8s%s", "'",theValue_temp,"'" );
271  bytesAdded += addKeyword_txt( "TFORM10" , theValue, "data format of the field: 4-byte REAL", 0 );
272 
273  bytesAdded += addKeyword_txt( "EXTNAME" , "'PARAMETERS'", "name of this binary table extension", 0 );
274  bytesAdded += addKeyword_txt( "HDUCLASS", "'OGIP '", "Format conforms to OGIP/GSFC conventions", 0 );
275  bytesAdded += addKeyword_txt( "HDUCLAS1", "'XSPEC TABLE MODEL'","model spectra for XSPEC", 0 );
276  bytesAdded += addKeyword_txt( "HDUCLAS2", "'PARAMETERS'", "Extension containing paramter info", 0 );
277  bytesAdded += addKeyword_txt( "HDUVERS" , "'1.0.0 '", "Version of format (OGIP memo OGIP-92-001)", 0 );
278  bytesAdded += addKeyword_num( "NINTPARM", nintparm, "Number of interpolation parameters" );
279  bytesAdded += addKeyword_num( "NADDPARM", naddparm, "Number of additional parameters" );
280  /* After everything else */
281  bytesAdded += fprintf(ioFITS_OUTPUT, "%-80s", "END" );
282 
283  ASSERT( bytesAdded%LINESIZE == 0 );
284 
285  /* Now add blanks */
286  while( bytesAdded%RECORDSIZE > 0 )
287  {
288  bytesAdded += fprintf(ioFITS_OUTPUT, "%-1s", " " );
289  }
290  return;
291 }
292 
293 STATIC void punchFITS_ParamData( char **paramNames,
294  long *paramMethods,
295  realnum **paramRange,
296  realnum **paramData,
297  long nintparm,
298  long naddparm,
299  long *numParamValues )
300 {
301  long i, j;
302 
303  DEBUG_ENTRY( "punchFITS_ParamData()" );
304 
305  ASSERT( nintparm+naddparm <= LIMPAR );
306 
307  /* Now add the parameters data */
308  for( i=0; i<nintparm+naddparm; i++ )
309  {
310  int32 numTemp;
311 
312 #define LOG2LINEAR 0
313 
314  paramMethods[i] = HtoNL(paramMethods[i]);
315  /* >>chng 06 aug 23, numParamValues is now an array. */
316  numTemp = HtoNL(numParamValues[i]);
317 
318 #if LOG2LINEAR
319  /* change to linear */
320  paramRange[i][0] = (realnum)pow( 10., (double)paramRange[i][0] );
321  paramRange[i][1] = (realnum)pow( 10., (double)paramRange[i][1] );
322  paramRange[i][2] = (realnum)pow( 10., (double)paramRange[i][2] );
323  paramRange[i][3] = (realnum)pow( 10., (double)paramRange[i][3] );
324  paramRange[i][4] = (realnum)pow( 10., (double)paramRange[i][4] );
325  paramRange[i][5] = (realnum)pow( 10., (double)paramRange[i][5] );
326 #endif
327 
328 #if !defined(_BIG_ENDIAN)
329  ByteSwap5( paramRange[i][0] );
330  ByteSwap5( paramRange[i][1] );
331  ByteSwap5( paramRange[i][2] );
332  ByteSwap5( paramRange[i][3] );
333  ByteSwap5( paramRange[i][4] );
334  ByteSwap5( paramRange[i][5] );
335 #endif
336 
337  /* >>chng 06 aug 23, numParamValues is now an array. */
338  for( j=0; j<numParamValues[i]; j++ )
339  {
340 
341 #if LOG2LINEAR
342  paramData[i][j] = (realnum)pow( 10., (double)paramData[i][j] );
343 #endif
344 
345 #if !defined(_BIG_ENDIAN)
346  ByteSwap5( paramData[i][j] );
347 #endif
348  }
349 
350  bytesAdded += fprintf(ioFITS_OUTPUT, "%-12s", paramNames[i] );
351  bytesAdded += (long)fwrite( &paramMethods[i], 1, sizeof(int32), ioFITS_OUTPUT );
352  bytesAdded += (long)fwrite( paramRange[i], 1, 6*sizeof(realnum), ioFITS_OUTPUT );
353  bytesAdded += (long)fwrite( &numTemp, 1, sizeof(int32), ioFITS_OUTPUT );
354  /* >>chng 06 aug 23, numParamValues is now an array. */
355  bytesAdded += (long)fwrite( paramData[i], 1, (unsigned)numParamValues[i]*sizeof(realnum), ioFITS_OUTPUT );
356 
357  for( j=numParamValues[i]+1; j<=maxParamValues; j++ )
358  {
359  realnum filler = -10.f;
360  bytesAdded += (long)fwrite( &filler, 1, sizeof(realnum), ioFITS_OUTPUT );
361  }
362  }
363 
364  /* Switch the endianness again */
365  for( i=0; i<nintparm+naddparm; i++ )
366  {
367  paramMethods[i] = HtoNL(paramMethods[i]);
368 
369 #if !defined(_BIG_ENDIAN)
370  ByteSwap5( paramRange[i][0] );
371  ByteSwap5( paramRange[i][1] );
372  ByteSwap5( paramRange[i][2] );
373  ByteSwap5( paramRange[i][3] );
374  ByteSwap5( paramRange[i][4] );
375  ByteSwap5( paramRange[i][5] );
376 #endif
377 
378  /* >>chng 06 aug 23, numParamValues is now an array. */
379  for( j=0; j<numParamValues[i]; j++ )
380  {
381 #if !defined(_BIG_ENDIAN)
382  ByteSwap5( paramData[i][j] );
383 #endif
384  }
385  }
386 
387  while( bytesAdded%RECORDSIZE > 0 )
388  {
389  int tempInt = 0;
390  bytesAdded += (long)fwrite( &tempInt, 1, 1, ioFITS_OUTPUT );
391  }
392  return;
393 }
394 
395 STATIC void punchFITS_EnergyHeader( long numEnergies )
396 {
397  long numFields = 2;
398  long naxis, naxis1, naxis2;
399 
400  DEBUG_ENTRY( "punchFITS_EnergyHeader()" );
401 
402  /* Make sure the previous blocks are the right size */
403  ASSERT( bytesAdded%RECORDSIZE == 0 );
404 
405  naxis = 2;
406  naxis1 = 2*sizeof(realnum);
407  naxis2 = numEnergies;
408 
409  bytesAdded += addKeyword_txt( "XTENSION", "'BINTABLE'", "binary table extension", 0 );
410  bytesAdded += addKeyword_num( "BITPIX" , bitpix, "8-bit bytes" );
411  bytesAdded += addKeyword_num( "NAXIS" , naxis, "2-dimensional binary table" );
412  bytesAdded += addKeyword_num( "NAXIS1" , naxis1, "width of table in bytes" );
413  bytesAdded += addKeyword_num( "NAXIS2" , naxis2, "number of rows in table" );
414  bytesAdded += addKeyword_num( "PCOUNT" , pcount, "size of special data area" );
415  bytesAdded += addKeyword_num( "GCOUNT" , gcount, "one data group (required keyword)" );
416  bytesAdded += addKeyword_num( "TFIELDS" , numFields, "number of fields in each row" );
417  bytesAdded += addKeyword_txt( "TTYPE1" , "'ENERG_LO'", "label for field 1", 0 );
418  bytesAdded += addKeyword_txt( "TFORM1" , "'E '", "data format of the field: 4-byte REAL", 0 );
419  bytesAdded += addKeyword_txt( "TTYPE2" , "'ENERG_HI'", "label for field 2", 0 );
420  bytesAdded += addKeyword_txt( "TFORM2" , "'E '", "data format of the field: 4-byte REAL", 0 );
421  bytesAdded += addKeyword_txt( "EXTNAME" , "'ENERGIES'", "name of this binary table extension", 0 );
422  bytesAdded += addKeyword_txt( "HDUCLASS", "'OGIP '", "Format conforms to OGIP/GSFC conventions", 0 );
423  bytesAdded += addKeyword_txt( "HDUCLAS1", "'XSPEC TABLE MODEL'","model spectra for XSPEC", 0 );
424  bytesAdded += addKeyword_txt( "HDUCLAS2", "'ENERGIES'", "Extension containing energy bin info", 0 );
425  bytesAdded += addKeyword_txt( "HDUVERS" , "'1.0.0 '", "Version of format (OGIP memo OGIP-92-001)", 0 );
426  /* After everything else */
427  bytesAdded += fprintf(ioFITS_OUTPUT, "%-80s", "END" );
428 
429  ASSERT( bytesAdded%LINESIZE == 0 );
430 
431  while( bytesAdded%RECORDSIZE > 0 )
432  {
433  bytesAdded += fprintf(ioFITS_OUTPUT, "%-1s", " " );
434  }
435  return;
436 }
437 
438 STATIC void punchFITS_EnergyData( realnum *Energies, long numEnergies )
439 {
440  long i;
441 
442  DEBUG_ENTRY( "punchFITS_EnergyData()" );
443 
444  /* Now add the energies data */
445  for( i=0; i<numEnergies; i++ )
446  {
447  realnum EnergyLow, EnergyHi;
448  /* Convert all of these to kev */
449  EnergyLow = 0.001f*(realnum)EVRYD*(Energies[i]-rfield.widflx[i]/2.f);
450 
451  if( i == numEnergies-1 )
452  {
453  EnergyHi = 0.001f*(realnum)EVRYD*(Energies[i] + rfield.widflx[i]/2.f);
454  }
455  else
456  {
457  EnergyHi = 0.001f*(realnum)EVRYD*(Energies[i+1] - rfield.widflx[i+1]/2.f);
458  }
459 
460 #if !defined(_BIG_ENDIAN)
461  ByteSwap5(EnergyLow);
462  ByteSwap5(EnergyHi);
463 #endif
464 
465  bytesAdded += (long)fwrite( &EnergyLow , 1, sizeof(realnum), ioFITS_OUTPUT );
466  bytesAdded += (long)fwrite( &EnergyHi , 1, sizeof(realnum), ioFITS_OUTPUT );
467  }
468 
469  while( bytesAdded%RECORDSIZE > 0 )
470  {
471  int tempInt = 0;
472  bytesAdded += (long)fwrite( &tempInt, 1, 1, ioFITS_OUTPUT );
473  }
474  return;
475 }
476 
477 STATIC void punchFITS_SpectraHeader( bool lgAddModel, long nintparm, long naddparm, long totNumModels, long numEnergies )
478 {
479  long i, numFields = 2+naddparm;
480  long naxis, naxis1, naxis2;
481  char theKeyword1[8];
482  char theKeyword2[8];
483  char theKeyword3[8];
484  char theValue1[20];
485  char theValue2[20];
486  char theValue2temp[20];
487  char theValue[20];
488  char theValue_temp[20];
489  char theComment1[47];
490 
491  DEBUG_ENTRY( "punchFITS_SpectraHeader()" );
492 
493  ASSERT( nintparm + naddparm <= LIMPAR );
494 
495  /* Make sure the previous blocks are the right size */
496  ASSERT( bytesAdded%RECORDSIZE == 0 );
497 
498  naxis = 2;
499  naxis1 = ( numEnergies*(naddparm+1) + nintparm ) * (long)sizeof(realnum);
500  naxis2 = totNumModels;
501 
502  bytesAdded += addKeyword_txt( "XTENSION", "'BINTABLE'", "binary table extension", 0 );
503  bytesAdded += addKeyword_num( "BITPIX" , bitpix, "8-bit bytes" );
504  bytesAdded += addKeyword_num( "NAXIS" , naxis, "2-dimensional binary table" );
505  bytesAdded += addKeyword_num( "NAXIS1" , naxis1, "width of table in bytes" );
506  bytesAdded += addKeyword_num( "NAXIS2" , naxis2, "number of rows in table" );
507  bytesAdded += addKeyword_num( "PCOUNT" , pcount, "size of special data area" );
508  bytesAdded += addKeyword_num( "GCOUNT" , gcount, "one data group (required keyword)" );
509  bytesAdded += addKeyword_num( "TFIELDS" , numFields, "number of fields in each row" );
510 
511  /******************************************/
512  /* These are the interpolation parameters */
513  /******************************************/
514  bytesAdded += addKeyword_txt( "TTYPE1" , "'PARAMVAL'", "label for field 1", 0 );
515  /* The size of this array is dynamic, set to size of nintparm */
516  sprintf( theValue2temp, "%ld%s", nintparm, "E" );
517  sprintf( theValue2, "%s%-8s%s", "'",theValue2temp,"'" );
518  bytesAdded += addKeyword_txt( "TFORM1" , theValue2, "data format of the field: 4-byte REAL", 0 );
519 
520  /******************************************/
521  /* This is the interpolated spectrum */
522  /******************************************/
523  bytesAdded += addKeyword_txt( "TTYPE2" , "'INTPSPEC'", "label for field 2", 0 );
524  /* The size of this array is dynamic, set to size of numEnergies */
525  sprintf( theValue_temp, "%ld%s", numEnergies, "E" );
526  sprintf( theValue, "%s%-8s%s", "'",theValue_temp,"'" );
527  bytesAdded += addKeyword_txt( "TFORM2" , theValue, "data format of the field: 4-byte REAL", 0 );
528  bytesAdded += addKeyword_txt( "TUNIT2" , ModelUnits[lgAddModel], "physical unit of field", 0 );
529 
530  /******************************************/
531  /* These are the additional parameters */
532  /******************************************/
533  for( i=1; i<=naddparm; i++ )
534  {
535  sprintf( theKeyword1, "%s%ld", "TTYPE", i+2 );
536  sprintf( theKeyword2, "%s%ld", "TFORM", i+2 );
537  sprintf( theKeyword3, "%s%ld", "TUNIT", i+2 );
538 
539  sprintf( theValue1, "%s%2.2ld%s", "'ADDSP", i, "'" );
540  /* The size of this array is dynamic, set to size of numEnergies */
541  sprintf( theValue2temp, "%ld%s", numEnergies, "E" );
542  sprintf( theValue2, "%s%-8s%s", "'",theValue2temp,"'" );
543 
544  sprintf( theComment1, "%s%ld", "label for field ", i+2 );
545 
546  bytesAdded += addKeyword_txt( theKeyword1 , theValue1, theComment1, 0 );
547  bytesAdded += addKeyword_txt( theKeyword2 , theValue2, "data format of the field: 4-byte REAL", 0 );
548  bytesAdded += addKeyword_txt( theKeyword3 , ModelUnits[lgAddModel], "physical unit of field", 0 );
549  }
550 
551  bytesAdded += addKeyword_txt( "EXTNAME" , "'SPECTRA '", "name of this binary table extension", 0 );
552  bytesAdded += addKeyword_txt( "HDUCLASS", "'OGIP '", "Format conforms to OGIP/GSFC conventions", 0 );
553  bytesAdded += addKeyword_txt( "HDUCLAS1", "'XSPEC TABLE MODEL'","model spectra for XSPEC", 0 );
554  bytesAdded += addKeyword_txt( "HDUCLAS2", "'MODEL SPECTRA'", "Extension containing model spectra", 0 );
555  bytesAdded += addKeyword_txt( "HDUVERS" , "'1.0.0 '", "Version of format (OGIP memo OGIP-92-001)", 0 );
556  /* After everything else */
557  bytesAdded += fprintf(ioFITS_OUTPUT, "%-80s", "END" );
558 
559  ASSERT( bytesAdded%LINESIZE == 0 );
560 
561  while( bytesAdded%RECORDSIZE > 0 )
562  {
563  bytesAdded += fprintf(ioFITS_OUTPUT, "%-1s", " " );
564  }
565  return;
566 }
567 
568 STATIC void punchFITS_SpectraData( realnum **interpParameters, realnum **theSpectrum,
569  long totNumModels, long numEnergies, long nintparm, long naddparm )
570 {
571  long i;
572  long naxis2 = totNumModels;
573 
574  DEBUG_ENTRY( "punchFITS_SpectraData()" );
575 
576  ASSERT( nintparm + naddparm <= LIMPAR );
577 
578  /* Now add the spectra data */
579  for( i=0; i<naxis2; i++ )
580  {
581 
582 #if !defined(_BIG_ENDIAN)
583  for( long j = 0; j<numEnergies; j++ )
584  {
585  ByteSwap5( theSpectrum[i][j] );
586  }
587 
588  for( long j = 0; j<nintparm; j++ )
589  {
590  ByteSwap5( interpParameters[i][j] );
591  }
592 #endif
593 
594  /* The interpolation parameters vector */
595  bytesAdded += (long)fwrite( interpParameters[i], 1, (unsigned)nintparm*sizeof(realnum), ioFITS_OUTPUT );
596  /* The interpolated spectrum */
597  bytesAdded += (long)fwrite( theSpectrum[i], 1, (unsigned)numEnergies*sizeof(realnum), ioFITS_OUTPUT );
598 
599 #if !defined(_BIG_ENDIAN)
600  /* Switch the endianness back to native. */
601  for( long j = 0; j<numEnergies; j++ )
602  {
603  ByteSwap5( theSpectrum[i][j] );
604  }
605 
606  for( long j = 0; j<nintparm; j++ )
607  {
608  ByteSwap5( interpParameters[i][j] );
609  }
610 #endif
611 
612  /* >>chng 06 aug 23, disable additional parameters for now */
613  if( naddparm > 0 )
614  {
615  /* The additional parameters */
616 
617  /* bytesAdded += (long)fwrite( theSpectrum[i], 1, (unsigned)numEnergies*sizeof(realnum), ioFITS_OUTPUT ); */
618  /* \todo 2 must create another array if we are to punch additional parameter information. */
619  fprintf( ioQQQ, " Additional parameters not currently supported.\n" );
620  cdEXIT( EXIT_FAILURE );
621  }
622  }
623 
624  while( bytesAdded%RECORDSIZE > 0 )
625  {
626  int tempInt = 0;
627  bytesAdded += (long)fwrite( &tempInt, 1, 1, ioFITS_OUTPUT );
628  }
629  return;
630 }
631 
632 STATIC void punchFITS_GenericHeader( long numEnergies )
633 {
634  long numFields = 2;
635  long naxis, naxis1, naxis2;
636 
637  DEBUG_ENTRY( "punchFITS_GenericHeader()" );
638 
639  /* Make sure the previous blocks are the right size */
640  ASSERT( bytesAdded%RECORDSIZE == 0 );
641 
642  naxis = 2;
643  naxis1 = numFields*(long)sizeof(realnum);
644  naxis2 = numEnergies;
645 
646  bytesAdded += addKeyword_txt( "XTENSION", "'BINTABLE'", "binary table extension", 0 );
647  bytesAdded += addKeyword_num( "BITPIX" , bitpix, "8-bit bytes" );
648  bytesAdded += addKeyword_num( "NAXIS" , naxis, "2-dimensional binary table" );
649  bytesAdded += addKeyword_num( "NAXIS1" , naxis1, "width of table in bytes" );
650  bytesAdded += addKeyword_num( "NAXIS2" , naxis2, "number of rows in table" );
651  bytesAdded += addKeyword_num( "PCOUNT" , pcount, "size of special data area" );
652  bytesAdded += addKeyword_num( "GCOUNT" , gcount, "one data group (required keyword)" );
653  bytesAdded += addKeyword_num( "TFIELDS" , numFields, "number of fields in each row" );
654  bytesAdded += addKeyword_txt( "TTYPE1" , "'ENERGY '", "label for field 1", 0 );
655  bytesAdded += addKeyword_txt( "TFORM1" , "'E '", "data format of the field: 4-byte REAL", 0 );
656  bytesAdded += addKeyword_txt( "TTYPE2" , "'TRN_SPEC'", "label for field 2", 0 );
657  bytesAdded += addKeyword_txt( "TFORM2" , "'E '", "data format of the field: 4-byte REAL", 0 );
658  bytesAdded += addKeyword_txt( "EXTNAME" , "'SPECTRA '", "name of this binary table extension", 0 );
659  bytesAdded += addKeyword_txt( "HDUCLASS", "'OGIP '", "Format conforms to OGIP/GSFC conventions", 0 );
660  bytesAdded += addKeyword_txt( "HDUCLAS1", "'XSPEC TABLE MODEL'","model spectra for XSPEC", 0 );
661  bytesAdded += addKeyword_txt( "HDUCLAS2", "'ENERGIES'", "Extension containing energy bin info", 0 );
662  bytesAdded += addKeyword_txt( "HDUVERS" , "'1.0.0 '", "Version of format (OGIP memo OGIP-92-001)", 0 );
663  /* After everything else */
664  bytesAdded += fprintf(ioFITS_OUTPUT, "%-80s", "END" );
665 
666  ASSERT( bytesAdded%LINESIZE == 0 );
667 
668  while( bytesAdded%RECORDSIZE > 0 )
669  {
670  bytesAdded += fprintf(ioFITS_OUTPUT, "%-1s", " " );
671  }
672  return;
673 }
674 
675 STATIC void punchFITS_GenericData( long numEnergies )
676 {
677  long i;
678 
679  DEBUG_ENTRY( "punchFITS_GenericData()" );
680 
681  realnum *TransmittedSpectrum;
682 
683  TransmittedSpectrum = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(numEnergies) );
684 
685  cdSPEC2( 8, numEnergies, TransmittedSpectrum );
686 
687  /* Now add the energies data */
688  for( i=0; i<numEnergies; i++ )
689  {
690  realnum Energy;
691  Energy = rfield.AnuOrg[i];
692 
693 #if !defined(_BIG_ENDIAN)
694  ByteSwap5(Energy);
695  ByteSwap5(TransmittedSpectrum[i]);
696 #endif
697 
698  bytesAdded += (long)fwrite( &Energy , 1, sizeof(realnum), ioFITS_OUTPUT );
699  bytesAdded += (long)fwrite( &TransmittedSpectrum[i],1, sizeof(realnum), ioFITS_OUTPUT );
700  }
701 
702  while( bytesAdded%RECORDSIZE > 0 )
703  {
704  int tempInt = 0;
705  bytesAdded += (long)fwrite( &tempInt, 1, 1, ioFITS_OUTPUT );
706  }
707 
708  free( TransmittedSpectrum );
709  return;
710 }
711 
713 {
714  char timeString[30]="";
715  char tempString[70];
716  time_t now;
717  long i, j, k;
718 
719  /* usually print date and time info - do not if "no times" command entered,
720  * which set this flag false */
721  now = time(NULL);
722  if( prt.lgPrintTime )
723  {
724  /* now add date of this run */
725  /* now print this time at the end of the string. the system put cr at the end of the string */
726  strcpy( timeString , ctime(&now) );
727  }
728  /* ctime puts a carriage return at the end, but we can't have that in a fits file.
729  * remove the carriage return here. */
730  for( i=0; i<30; i++ )
731  {
732  if( timeString[i] == '\n' )
733  {
734  timeString[i] = ' ';
735  }
736  }
737 
738  strcpy( tempString, "Generated by Cloudy " );
739  // strncat guarantees that terminating 0 byte will always be written...
740  strncat( tempString, t_version::Inst().chVersion, sizeof(tempString)-strlen(tempString) );
741  bytesAdded += addComment( tempString );
742  bytesAdded += addComment( t_version::Inst().chInfo );
743  strcpy( tempString, "--- " );
744  strcat( tempString, timeString );
745  bytesAdded += addComment( tempString );
746  bytesAdded += addComment( "Input string was as follows: " );
747  /* >>chng 05 nov 24, from <nSave to <=nSave bug caught by PvH */
748  for( i=0; i<=input.nSave; i++ )
749  {
750  char firstLine[70], extraLine[65];
751 
752  for( j=0; j<INPUT_LINE_LENGTH; j++ )
753  {
754  if( input.chCardSav[i][j] =='\0' )
755  break;
756  }
757 
758  ASSERT( j < 200 );
759  for( k=0; k< MIN2(69, j); k++ )
760  {
761  firstLine[k] = input.chCardSav[i][k];
762  }
763  firstLine[k] = '\0';
764  bytesAdded += addComment( firstLine );
765  if( j >= 69 )
766  {
767  for( k=69; k< 133; k++ )
768  {
769  extraLine[k-69] = input.chCardSav[i][k];
770  }
771  /* >> chng 06 jan 05, this was exceeding array bounds. */
772  extraLine[64] = '\0';
773  strcpy( tempString, "more " );
774  strcat( tempString, extraLine );
775  bytesAdded += addComment( tempString );
776  if( j >= 133 )
777  {
778  for( k=133; k< 197; k++ )
779  {
780  extraLine[k-133] = input.chCardSav[i][k];
781  }
782  extraLine[64] = '\0';
783  strcpy( tempString, "more " );
784  strcat( tempString, extraLine );
785  bytesAdded += addComment( tempString );
786  }
787  }
788  }
789 
790  return;
791 }
792 
793 STATIC long addKeyword_txt( const char *theKeyword, const void *theValue, const char *theComment, long Str_Or_Log )
794 {
795  long numberOfBytesWritten = 0;
796 
797  DEBUG_ENTRY( "addKeyword_txt()" );
798 
799  /* False means string, true means logical */
800  if( Str_Or_Log == 0 )
801  {
802  numberOfBytesWritten = fprintf(ioFITS_OUTPUT, "%-8s%-2s%-20s%3s%-47s",
803  theKeyword,
804  "= ",
805  (char *)theValue,
806  " / ",
807  theComment );
808  }
809  else
810  {
811  ASSERT( Str_Or_Log == 1 );
812  numberOfBytesWritten = fprintf(ioFITS_OUTPUT, "%-8s%-2s%20s%3s%-47s",
813  theKeyword,
814  "= ",
815  (char *)theValue,
816  " / ",
817  theComment );
818  }
819 
820  ASSERT( numberOfBytesWritten%LINESIZE == 0 );
821  return numberOfBytesWritten;
822 }
823 
824 STATIC long addKeyword_num( const char *theKeyword, long theValue, const char *theComment)
825 {
826  long numberOfBytesWritten = 0;
827 
828  DEBUG_ENTRY( "addKeyword_num()" );
829 
830  numberOfBytesWritten = fprintf(ioFITS_OUTPUT, "%-8s%-2s%20ld%3s%-47s",
831  theKeyword,
832  "= ",
833  theValue,
834  " / ",
835  theComment );
836 
837  ASSERT( numberOfBytesWritten%LINESIZE == 0 );
838  return numberOfBytesWritten;
839 }
840 
841 long addComment( const char *CommentToAdd )
842 {
843  long i, numberOfBytesWritten = 0;
844  char tempString[70] = " ";
845 
846  DEBUG_ENTRY( "addComment()" );
847 
848  strncpy( &tempString[0], CommentToAdd, 69 );
849  ASSERT( (int)strlen( tempString ) <= 70 );
850 
851  /* tabs violate FITS standard, replace them with spaces. */
852  for( i=0; i<69; i++ )
853  {
854  if( tempString[i] == '\t' )
855  {
856  tempString[i] = ' ';
857  }
858  }
859 
860  numberOfBytesWritten = fprintf(ioFITS_OUTPUT, "COMMENT %-70s", tempString );
861 
862  ASSERT( numberOfBytesWritten%LINESIZE == 0 );
863  return numberOfBytesWritten;
864 }

Generated for cloudy by doxygen 1.8.3.1