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