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