cloudy trunk
Loading...
Searching...
No Matches
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
38error 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)
44STATIC 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
60static FILE *ioFITS_OUTPUT;
61static long bytesAdded = 0;
62static long bitpix = 8;
63static long pcount = 0;
64static long gcount = 1;
65static long maxParamValues = 0;
66const char ModelUnits[2][17] = {"'dimensionless '", "'photons/cm^2/s'" };
67
68STATIC void punchFITS_PrimaryHeader( bool lgAddModel );
69STATIC void punchFITS_ParamHeader( /* long *numParamValues, */ long nintparm, long naddparm );
70STATIC void punchFITS_ParamData( char **paramNames, long *paramMethods, realnum **paramRange,
71 realnum **paramData, long nintparm, long naddparm,
72 long *numParamValues );
73STATIC void punchFITS_EnergyHeader( long numEnergies );
74STATIC void punchFITS_EnergyData( const vector<realnum>& Energies, long EnergyOffset );
75STATIC void punchFITS_SpectraHeader( bool lgAdditiveModel, long nintparm, long naddparm, long totNumModels, long numEnergies );
76STATIC void punchFITS_SpectraData( realnum **interpParameters, multi_arr<realnum,3>& theSpectrum, int option,
77 long totNumModels, long numEnergies, long nintparm, long naddparm );
78STATIC void punchFITS_GenericHeader( long numEnergies );
79STATIC void punchFITS_GenericData( long numEnergies, long ipLoEnergy, long ipHiEnergy );
80STATIC void writeCloudyDetails( void );
81STATIC long addComment( const char *CommentToAdd );
82STATIC long addKeyword_txt( const char *theKeyword, const void *theValue, const char *theComment, long Str_Or_Log );
83STATIC long addKeyword_num( const char *theKeyword, long theValue, const char *theComment);
84
85void 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" );
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 {
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" );
172 }
173 return;
174}
175
176STATIC 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 */
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
222STATIC 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 */
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
295STATIC 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
397STATIC 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 */
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
440STATIC 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
478STATIC 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 */
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
569STATIC 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" );
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
633STATIC 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 */
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
676STATIC 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 {
692 Energy = rfield.AnuOrg[i];
693
694#if !defined(_BIG_ENDIAN)
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
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 );
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
794STATIC 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
825STATIC 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
842long 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}
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define MIN2
Definition cddefines.h:761
#define MALLOC(exp)
Definition cddefines.h:501
const int INPUT_LINE_LENGTH
Definition cddefines.h:254
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void fixit(void)
Definition service.cpp:991
void cdSPEC2(int Option, long int nEnergy, long int ipLoEnergy, long int ipHiEnergy, realnum ReturnedSpectrum[])
Definition energy.h:8
static t_version & Inst()
Definition cddefines.h:175
t_grid grid
Definition grid.cpp:5
const int NUM_OUTPUT_TYPES
Definition grid.h:21
t_input input
Definition input.cpp:12
const long LIMPAR
Definition optimize.h:61
UNUSED const double EVRYD
Definition physconst.h:189
t_prt prt
Definition prt.cpp:10
t_rfield rfield
Definition rfield.cpp:8
#define LINESIZE
Definition save_fits.cpp:15
static long gcount
Definition save_fits.cpp:64
static long pcount
Definition save_fits.cpp:63
#define RECORDSIZE
Definition save_fits.cpp:14
STATIC void punchFITS_EnergyHeader(long numEnergies)
STATIC long addComment(const char *CommentToAdd)
STATIC void punchFITS_GenericHeader(long numEnergies)
void saveFITSfile(FILE *ioPUN, int option)
Definition save_fits.cpp:85
const char ModelUnits[2][17]
Definition save_fits.cpp:66
static long bitpix
Definition save_fits.cpp:62
STATIC long addKeyword_num(const char *theKeyword, long theValue, const char *theComment)
static long bytesAdded
Definition save_fits.cpp:61
STATIC void punchFITS_EnergyData(const vector< realnum > &Energies, long EnergyOffset)
STATIC void writeCloudyDetails(void)
STATIC void punchFITS_PrimaryHeader(bool lgAddModel)
STATIC void punchFITS_GenericData(long numEnergies, long ipLoEnergy, long ipHiEnergy)
STATIC void punchFITS_ParamHeader(long nintparm, long naddparm)
STATIC long addKeyword_txt(const char *theKeyword, const void *theValue, const char *theComment, long Str_Or_Log)
static FILE * ioFITS_OUTPUT
Definition save_fits.cpp:60
STATIC void punchFITS_SpectraHeader(bool lgAdditiveModel, long nintparm, long naddparm, long totNumModels, long numEnergies)
#define ByteSwap5(x)
Definition save_fits.cpp:41
STATIC void punchFITS_SpectraData(realnum **interpParameters, multi_arr< realnum, 3 > &theSpectrum, int option, long totNumModels, long numEnergies, long nintparm, long naddparm)
static long maxParamValues
Definition save_fits.cpp:65
STATIC void ByteSwap(unsigned char *b, int n)
Definition save_fits.cpp:44
STATIC void punchFITS_ParamData(char **paramNames, long *paramMethods, realnum **paramRange, realnum **paramData, long nintparm, long naddparm, long *numParamValues)
#define HtoNL(A)
Definition save_fits.cpp:28