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