00001
00002
00003 #include "cddefines.h"
00004 #include "cddrive.h"
00005 #include "optimize.h"
00006 #include "grid.h"
00007 #include "save.h"
00008 #include "rfield.h"
00009 #include "prt.h"
00010 #include "input.h"
00011 #include "version.h"
00012 #include "physconst.h"
00013
00014 #define RECORDSIZE 2880
00015 #define LINESIZE 80
00016
00017 #if defined(_BIG_ENDIAN)
00018
00019 # define HtoNL(A) (A)
00020
00021
00022
00023
00024
00025 #else
00026
00027
00028 # define HtoNL(A) ((((A) & 0xff000000) >> 24) | \
00029 (((A) & 0x00ff0000) >> 8) | \
00030 (((A) & 0x0000ff00) << 8) | \
00031 (((A) & 0x000000ff) << 24))
00032
00033
00034
00035
00036
00037
00038
00039 #endif
00040
00041 #define ByteSwap5(x) ByteSwap((unsigned char *) &x,sizeof(x))
00042
00043 #if !defined(_BIG_ENDIAN)
00044 STATIC void ByteSwap(unsigned char * b, int n)
00045 {
00046 register int i = 0;
00047 register int j = n-1;
00048 while (i<j)
00049 {
00050 char temp = b[i];
00051 b[i] = b[j];
00052 b[j] = temp;
00053
00054 i++, j--;
00055 }
00056 return;
00057 }
00058 #endif
00059
00060 static FILE *ioFITS_OUTPUT;
00061 static long bytesAdded = 0;
00062 static long bitpix = 8;
00063 static long pcount = 0;
00064 static long gcount = 1;
00065 static long maxParamValues = 0;
00066 const char ModelUnits[2][17] = {"'dimensionless '", "'photons/cm^2/s'" };
00067
00068 STATIC void punchFITS_PrimaryHeader( bool lgAddModel );
00069 STATIC void punchFITS_ParamHeader( long nintparm, long naddparm );
00070 STATIC void punchFITS_ParamData( char **paramNames, long *paramMethods, realnum **paramRange,
00071 realnum **paramData, long nintparm, long naddparm,
00072 long *numParamValues );
00073 STATIC void punchFITS_EnergyHeader( long numEnergies );
00074 STATIC void punchFITS_EnergyData( const vector<realnum>& Energies, long EnergyOffset );
00075 STATIC void punchFITS_SpectraHeader( bool lgAdditiveModel, long nintparm, long naddparm, long totNumModels, long numEnergies );
00076 STATIC void punchFITS_SpectraData( realnum **interpParameters, multi_arr<realnum,3>& theSpectrum, int option,
00077 long totNumModels, long numEnergies, long nintparm, long naddparm );
00078 STATIC void punchFITS_GenericHeader( long numEnergies );
00079 STATIC void punchFITS_GenericData( long numEnergies, long ipLoEnergy, long ipHiEnergy );
00080 STATIC void writeCloudyDetails( void );
00081 STATIC long addComment( const char *CommentToAdd );
00082 STATIC long addKeyword_txt( const char *theKeyword, const void *theValue, const char *theComment, long Str_Or_Log );
00083 STATIC long addKeyword_num( const char *theKeyword, long theValue, const char *theComment);
00084
00085 void saveFITSfile( FILE* ioPUN, int option )
00086 {
00087
00088 long i;
00089
00090 DEBUG_ENTRY( "saveFITSfile()" );
00091
00092 if( !grid.lgGridDone && option != NUM_OUTPUT_TYPES )
00093 {
00094
00095 return;
00096 }
00097
00098 ioFITS_OUTPUT = ioPUN;
00099
00100 #if 0
00101 {
00102
00103 long i,j;
00104 FILE *asciiDump;
00105
00106 if( (asciiDump = fopen( "gridspectra.con","w" ) ) == NULL )
00107 {
00108 fprintf( ioQQQ, "saveFITSfile could not open save file for writing.\nSorry.\n" );
00109 cdEXIT(EXIT_FAILURE);
00110 }
00111
00112 for( i=0; i<grid.numEnergies-1; i++ )
00113 {
00114 fprintf( asciiDump, "%.12e\t", grid.Energies[i] );
00115 for( j=0; j<grid.totNumModels; j++ )
00116 {
00117 fprintf( asciiDump, "%.12e\t", grid.Spectra[4][j][i] );
00118 }
00119 fprintf( asciiDump, "\n" );
00120 }
00121 fclose( asciiDump );
00122 }
00123 #endif
00124
00125 ASSERT( option >= 0 );
00126
00127
00128 if( option == NUM_OUTPUT_TYPES )
00129 {
00130 punchFITS_PrimaryHeader( false );
00131 punchFITS_GenericHeader( rfield.nflux - 1 );
00132 punchFITS_GenericData( rfield.nflux -1, 0, rfield.nflux -2 );
00133 }
00134
00135 else if( option < NUM_OUTPUT_TYPES )
00136 {
00137 bool lgAdditiveModel;
00138
00139
00140 if( option == 10 )
00141 {
00142
00143 lgAdditiveModel = false;
00144 }
00145 else
00146 {
00147 lgAdditiveModel = true;
00148 }
00149
00150 punchFITS_PrimaryHeader( lgAdditiveModel );
00151
00152 for( i=0; i<grid.nintparm+grid.naddparm; i++ )
00153 {
00154 maxParamValues = MAX2( maxParamValues, grid.numParamValues[i] );
00155 }
00156
00157 ASSERT( maxParamValues >= 2 );
00158
00159 punchFITS_ParamHeader( grid.nintparm, grid.naddparm );
00160 punchFITS_ParamData( grid.paramNames, grid.paramMethods, grid.paramRange, grid.paramData,
00161 grid.nintparm, grid.naddparm, grid.numParamValues );
00162 punchFITS_EnergyHeader( grid.numEnergies );
00163 punchFITS_EnergyData( grid.Energies, grid.ipLoEnergy );
00164 punchFITS_SpectraHeader( lgAdditiveModel, grid.nintparm, grid.naddparm, grid.totNumModels, grid.numEnergies);
00165 punchFITS_SpectraData( grid.interpParameters, grid.Spectra, option,
00166 grid.totNumModels, grid.numEnergies, grid.nintparm, grid.naddparm );
00167 }
00168 else
00169 {
00170 fprintf( ioQQQ, "PROBLEM - undefined option encountered in saveFITSfile. \n" );
00171 cdEXIT( EXIT_FAILURE );
00172 }
00173 return;
00174 }
00175
00176 STATIC void punchFITS_PrimaryHeader( bool lgAddModel )
00177 {
00178 const char *ModelName = "'CLOUDY'";
00179
00180 DEBUG_ENTRY( "punchFITS_PrimaryHeader()" );
00181
00182 bytesAdded = 0;
00183
00184 fixit();
00185
00186 bytesAdded += addKeyword_txt( "SIMPLE" , "T", "file does conform to FITS standard", 1 );
00187 bytesAdded += addKeyword_num( "BITPIX" , bitpix, "number of bits per data pixel" );
00188 bytesAdded += addKeyword_num( "NAXIS" , 0, "number of data axes" );
00189 bytesAdded += addKeyword_txt( "EXTEND" , "T", "FITS dataset may contain extensions", 1 );
00190 bytesAdded += addKeyword_txt( "CONTENT" , "'MODEL '", "spectrum file contains time intervals and event", 0 );
00191 bytesAdded += addKeyword_txt( "MODLNAME", ModelName, "Model name", 0 );
00192 bytesAdded += addKeyword_txt( "MODLUNIT", ModelUnits[lgAddModel], "Model units", 0 );
00193 bytesAdded += addKeyword_txt( "REDSHIFT", "T", "If true then redshift will be included as a par", 1 );
00194 if( lgAddModel == true )
00195 {
00196 bytesAdded += addKeyword_txt( "ADDMODEL", "T", "If true then this is an additive table model", 1 );
00197 }
00198 else
00199 {
00200 bytesAdded += addKeyword_txt( "ADDMODEL", "F", "If true then this is an additive table model", 1 );
00201 }
00202
00203
00204 writeCloudyDetails();
00205
00206 bytesAdded += addKeyword_txt( "HDUCLASS", "'OGIP '", "Format conforms to OGIP/GSFC conventions", 0 );
00207 bytesAdded += addKeyword_txt( "HDUCLAS1", "'XSPEC TABLE MODEL'","Extension contains an image", 0 );
00208 bytesAdded += addKeyword_txt( "HDUVERS" , "'1.0.0 '", "Version of format (OGIP memo OGIP-92-001)", 0 );
00209
00210 bytesAdded += fprintf(ioFITS_OUTPUT, "%-80s", "END" );
00211
00212 ASSERT( bytesAdded%LINESIZE == 0 );
00213
00214
00215 while( bytesAdded%RECORDSIZE > 0 )
00216 {
00217 bytesAdded += fprintf(ioFITS_OUTPUT, "%-1s", " " );
00218 }
00219 return;
00220 }
00221
00222 STATIC void punchFITS_ParamHeader( long nintparm, long naddparm )
00223 {
00224 long numFields = 10;
00225 long naxis, naxis1, naxis2;
00226 char theValue[20];
00227 char theValue_temp[20];
00228
00229 DEBUG_ENTRY( "punchFITS_ParamHeader()" );
00230
00231 ASSERT( nintparm+naddparm <= LIMPAR );
00232
00233
00234 ASSERT( bytesAdded%RECORDSIZE == 0 );
00235
00236 naxis = 2;
00237
00238 naxis1 = 44+4*maxParamValues;
00239 naxis2 = nintparm+naddparm;
00240
00241 bytesAdded += addKeyword_txt( "XTENSION", "'BINTABLE'", "binary table extension", 0 );
00242 bytesAdded += addKeyword_num( "BITPIX" , bitpix, "8-bit bytes" );
00243 bytesAdded += addKeyword_num( "NAXIS" , naxis, "2-dimensional binary table" );
00244 bytesAdded += addKeyword_num( "NAXIS1" , naxis1, "width of table in bytes" );
00245 bytesAdded += addKeyword_num( "NAXIS2" , naxis2, "number of rows in table" );
00246 bytesAdded += addKeyword_num( "PCOUNT" , pcount, "size of special data area" );
00247 bytesAdded += addKeyword_num( "GCOUNT" , gcount, "one data group (required keyword)" );
00248 bytesAdded += addKeyword_num( "TFIELDS" , numFields, "number of fields in each row" );
00249 bytesAdded += addKeyword_txt( "TTYPE1" , "'NAME '", "label for field 1", 0 );
00250 bytesAdded += addKeyword_txt( "TFORM1" , "'12A '", "data format of the field: ASCII Character", 0 );
00251 bytesAdded += addKeyword_txt( "TTYPE2" , "'METHOD '", "label for field 2", 0 );
00252 bytesAdded += addKeyword_txt( "TFORM2" , "'J '", "data format of the field: 4-byte INTEGER", 0 );
00253 bytesAdded += addKeyword_txt( "TTYPE3" , "'INITIAL '", "label for field 3", 0 );
00254 bytesAdded += addKeyword_txt( "TFORM3" , "'E '", "data format of the field: 4-byte REAL", 0 );
00255 bytesAdded += addKeyword_txt( "TTYPE4" , "'DELTA '", "label for field 4", 0 );
00256 bytesAdded += addKeyword_txt( "TFORM4" , "'E '", "data format of the field: 4-byte REAL", 0 );
00257 bytesAdded += addKeyword_txt( "TTYPE5" , "'MINIMUM '", "label for field 5", 0 );
00258 bytesAdded += addKeyword_txt( "TFORM5" , "'E '", "data format of the field: 4-byte REAL", 0 );
00259 bytesAdded += addKeyword_txt( "TTYPE6" , "'BOTTOM '", "label for field 6", 0 );
00260 bytesAdded += addKeyword_txt( "TFORM6" , "'E '", "data format of the field: 4-byte REAL", 0 );
00261 bytesAdded += addKeyword_txt( "TTYPE7" , "'TOP '", "label for field 7", 0 );
00262 bytesAdded += addKeyword_txt( "TFORM7" , "'E '", "data format of the field: 4-byte REAL", 0 );
00263 bytesAdded += addKeyword_txt( "TTYPE8" , "'MAXIMUM '", "label for field 8", 0 );
00264 bytesAdded += addKeyword_txt( "TFORM8" , "'E '", "data format of the field: 4-byte REAL", 0 );
00265 bytesAdded += addKeyword_txt( "TTYPE9" , "'NUMBVALS'", "label for field 9", 0 );
00266 bytesAdded += addKeyword_txt( "TFORM9" , "'J '", "data format of the field: 4-byte INTEGER", 0 );
00267 bytesAdded += addKeyword_txt( "TTYPE10" , "'VALUE '", "label for field 10", 0 );
00268
00269
00270
00271 sprintf( theValue_temp, "%ld%s", maxParamValues, "E" );
00272 sprintf( theValue, "%s%-8s%s", "'",theValue_temp,"'" );
00273 bytesAdded += addKeyword_txt( "TFORM10" , theValue, "data format of the field: 4-byte REAL", 0 );
00274
00275 bytesAdded += addKeyword_txt( "EXTNAME" , "'PARAMETERS'", "name of this binary table extension", 0 );
00276 bytesAdded += addKeyword_txt( "HDUCLASS", "'OGIP '", "Format conforms to OGIP/GSFC conventions", 0 );
00277 bytesAdded += addKeyword_txt( "HDUCLAS1", "'XSPEC TABLE MODEL'","model spectra for XSPEC", 0 );
00278 bytesAdded += addKeyword_txt( "HDUCLAS2", "'PARAMETERS'", "Extension containing paramter info", 0 );
00279 bytesAdded += addKeyword_txt( "HDUVERS" , "'1.0.0 '", "Version of format (OGIP memo OGIP-92-001)", 0 );
00280 bytesAdded += addKeyword_num( "NINTPARM", nintparm, "Number of interpolation parameters" );
00281 bytesAdded += addKeyword_num( "NADDPARM", naddparm, "Number of additional parameters" );
00282
00283 bytesAdded += fprintf(ioFITS_OUTPUT, "%-80s", "END" );
00284
00285 ASSERT( bytesAdded%LINESIZE == 0 );
00286
00287
00288 while( bytesAdded%RECORDSIZE > 0 )
00289 {
00290 bytesAdded += fprintf(ioFITS_OUTPUT, "%-1s", " " );
00291 }
00292 return;
00293 }
00294
00295 STATIC void punchFITS_ParamData( char **paramNames,
00296 long *paramMethods,
00297 realnum **paramRange,
00298 realnum **paramData,
00299 long nintparm,
00300 long naddparm,
00301 long *numParamValues )
00302 {
00303 long i, j;
00304
00305 DEBUG_ENTRY( "punchFITS_ParamData()" );
00306
00307 ASSERT( nintparm+naddparm <= LIMPAR );
00308
00309
00310 for( i=0; i<nintparm+naddparm; i++ )
00311 {
00312 int32 numTemp;
00313
00314 #define LOG2LINEAR 0
00315
00316 paramMethods[i] = HtoNL(paramMethods[i]);
00317
00318 numTemp = HtoNL(numParamValues[i]);
00319
00320 #if LOG2LINEAR
00321
00322 paramRange[i][0] = (realnum)pow( 10., (double)paramRange[i][0] );
00323 paramRange[i][1] = (realnum)pow( 10., (double)paramRange[i][1] );
00324 paramRange[i][2] = (realnum)pow( 10., (double)paramRange[i][2] );
00325 paramRange[i][3] = (realnum)pow( 10., (double)paramRange[i][3] );
00326 paramRange[i][4] = (realnum)pow( 10., (double)paramRange[i][4] );
00327 paramRange[i][5] = (realnum)pow( 10., (double)paramRange[i][5] );
00328 #endif
00329
00330 #if !defined(_BIG_ENDIAN)
00331 ByteSwap5( paramRange[i][0] );
00332 ByteSwap5( paramRange[i][1] );
00333 ByteSwap5( paramRange[i][2] );
00334 ByteSwap5( paramRange[i][3] );
00335 ByteSwap5( paramRange[i][4] );
00336 ByteSwap5( paramRange[i][5] );
00337 #endif
00338
00339
00340 for( j=0; j<numParamValues[i]; j++ )
00341 {
00342
00343 #if LOG2LINEAR
00344 paramData[i][j] = (realnum)pow( 10., (double)paramData[i][j] );
00345 #endif
00346
00347 #if !defined(_BIG_ENDIAN)
00348 ByteSwap5( paramData[i][j] );
00349 #endif
00350 }
00351
00352 bytesAdded += fprintf(ioFITS_OUTPUT, "%-12s", paramNames[i] );
00353 bytesAdded += (long)fwrite( ¶mMethods[i], 1, sizeof(int32), ioFITS_OUTPUT );
00354 bytesAdded += (long)fwrite( paramRange[i], 1, 6*sizeof(realnum), ioFITS_OUTPUT );
00355 bytesAdded += (long)fwrite( &numTemp, 1, sizeof(int32), ioFITS_OUTPUT );
00356
00357 bytesAdded += (long)fwrite( paramData[i], 1, (unsigned)numParamValues[i]*sizeof(realnum), ioFITS_OUTPUT );
00358
00359 for( j=numParamValues[i]+1; j<=maxParamValues; j++ )
00360 {
00361 realnum filler = -10.f;
00362 bytesAdded += (long)fwrite( &filler, 1, sizeof(realnum), ioFITS_OUTPUT );
00363 }
00364 }
00365
00366
00367 for( i=0; i<nintparm+naddparm; i++ )
00368 {
00369 paramMethods[i] = HtoNL(paramMethods[i]);
00370
00371 #if !defined(_BIG_ENDIAN)
00372 ByteSwap5( paramRange[i][0] );
00373 ByteSwap5( paramRange[i][1] );
00374 ByteSwap5( paramRange[i][2] );
00375 ByteSwap5( paramRange[i][3] );
00376 ByteSwap5( paramRange[i][4] );
00377 ByteSwap5( paramRange[i][5] );
00378 #endif
00379
00380
00381 for( j=0; j<numParamValues[i]; j++ )
00382 {
00383 #if !defined(_BIG_ENDIAN)
00384 ByteSwap5( paramData[i][j] );
00385 #endif
00386 }
00387 }
00388
00389 while( bytesAdded%RECORDSIZE > 0 )
00390 {
00391 int tempInt = 0;
00392 bytesAdded += (long)fwrite( &tempInt, 1, 1, ioFITS_OUTPUT );
00393 }
00394 return;
00395 }
00396
00397 STATIC void punchFITS_EnergyHeader( long numEnergies )
00398 {
00399 long numFields = 2;
00400 long naxis, naxis1, naxis2;
00401
00402 DEBUG_ENTRY( "punchFITS_EnergyHeader()" );
00403
00404
00405 ASSERT( bytesAdded%RECORDSIZE == 0 );
00406
00407 naxis = 2;
00408 naxis1 = 2*sizeof(realnum);
00409 naxis2 = numEnergies;
00410
00411 bytesAdded += addKeyword_txt( "XTENSION", "'BINTABLE'", "binary table extension", 0 );
00412 bytesAdded += addKeyword_num( "BITPIX" , bitpix, "8-bit bytes" );
00413 bytesAdded += addKeyword_num( "NAXIS" , naxis, "2-dimensional binary table" );
00414 bytesAdded += addKeyword_num( "NAXIS1" , naxis1, "width of table in bytes" );
00415 bytesAdded += addKeyword_num( "NAXIS2" , naxis2, "number of rows in table" );
00416 bytesAdded += addKeyword_num( "PCOUNT" , pcount, "size of special data area" );
00417 bytesAdded += addKeyword_num( "GCOUNT" , gcount, "one data group (required keyword)" );
00418 bytesAdded += addKeyword_num( "TFIELDS" , numFields, "number of fields in each row" );
00419 bytesAdded += addKeyword_txt( "TTYPE1" , "'ENERG_LO'", "label for field 1", 0 );
00420 bytesAdded += addKeyword_txt( "TFORM1" , "'E '", "data format of the field: 4-byte REAL", 0 );
00421 bytesAdded += addKeyword_txt( "TTYPE2" , "'ENERG_HI'", "label for field 2", 0 );
00422 bytesAdded += addKeyword_txt( "TFORM2" , "'E '", "data format of the field: 4-byte REAL", 0 );
00423 bytesAdded += addKeyword_txt( "EXTNAME" , "'ENERGIES'", "name of this binary table extension", 0 );
00424 bytesAdded += addKeyword_txt( "HDUCLASS", "'OGIP '", "Format conforms to OGIP/GSFC conventions", 0 );
00425 bytesAdded += addKeyword_txt( "HDUCLAS1", "'XSPEC TABLE MODEL'","model spectra for XSPEC", 0 );
00426 bytesAdded += addKeyword_txt( "HDUCLAS2", "'ENERGIES'", "Extension containing energy bin info", 0 );
00427 bytesAdded += addKeyword_txt( "HDUVERS" , "'1.0.0 '", "Version of format (OGIP memo OGIP-92-001)", 0 );
00428
00429 bytesAdded += fprintf(ioFITS_OUTPUT, "%-80s", "END" );
00430
00431 ASSERT( bytesAdded%LINESIZE == 0 );
00432
00433 while( bytesAdded%RECORDSIZE > 0 )
00434 {
00435 bytesAdded += fprintf(ioFITS_OUTPUT, "%-1s", " " );
00436 }
00437 return;
00438 }
00439
00440 STATIC void punchFITS_EnergyData( const vector<realnum>& Energies, long EnergyOffset )
00441 {
00442 DEBUG_ENTRY( "punchFITS_EnergyData()" );
00443
00444
00445 for( unsigned long i=0; i < Energies.size(); i++ )
00446 {
00447 realnum EnergyLow, EnergyHi;
00448 ASSERT( i+EnergyOffset < (unsigned)rfield.nupper );
00449
00450 EnergyLow = 0.001f*(realnum)EVRYD*(Energies[i]-rfield.widflx[i+EnergyOffset]/2.f);
00451
00452 if( i == Energies.size()-1 )
00453 {
00454 EnergyHi = 0.001f*(realnum)EVRYD*(Energies[i] + rfield.widflx[i+EnergyOffset]/2.f);
00455 }
00456 else
00457 {
00458 EnergyHi = 0.001f*(realnum)EVRYD*(Energies[i+1] - rfield.widflx[i+EnergyOffset+1]/2.f);
00459 }
00460
00461 #if !defined(_BIG_ENDIAN)
00462 ByteSwap5(EnergyLow);
00463 ByteSwap5(EnergyHi);
00464 #endif
00465
00466 bytesAdded += (long)fwrite( &EnergyLow , 1, sizeof(realnum), ioFITS_OUTPUT );
00467 bytesAdded += (long)fwrite( &EnergyHi , 1, sizeof(realnum), ioFITS_OUTPUT );
00468 }
00469
00470 while( bytesAdded%RECORDSIZE > 0 )
00471 {
00472 int tempInt = 0;
00473 bytesAdded += (long)fwrite( &tempInt, 1, 1, ioFITS_OUTPUT );
00474 }
00475 return;
00476 }
00477
00478 STATIC void punchFITS_SpectraHeader( bool lgAddModel, long nintparm, long naddparm, long totNumModels, long numEnergies )
00479 {
00480 long i, numFields = 2+naddparm;
00481 long naxis, naxis1, naxis2;
00482 char theKeyword1[8];
00483 char theKeyword2[8];
00484 char theKeyword3[8];
00485 char theValue1[20];
00486 char theValue2[20];
00487 char theValue2temp[20];
00488 char theValue[20];
00489 char theValue_temp[20];
00490 char theComment1[47];
00491
00492 DEBUG_ENTRY( "punchFITS_SpectraHeader()" );
00493
00494 ASSERT( nintparm + naddparm <= LIMPAR );
00495
00496
00497 ASSERT( bytesAdded%RECORDSIZE == 0 );
00498
00499 naxis = 2;
00500 naxis1 = ( numEnergies*(naddparm+1) + nintparm ) * (long)sizeof(realnum);
00501 naxis2 = totNumModels;
00502
00503 bytesAdded += addKeyword_txt( "XTENSION", "'BINTABLE'", "binary table extension", 0 );
00504 bytesAdded += addKeyword_num( "BITPIX" , bitpix, "8-bit bytes" );
00505 bytesAdded += addKeyword_num( "NAXIS" , naxis, "2-dimensional binary table" );
00506 bytesAdded += addKeyword_num( "NAXIS1" , naxis1, "width of table in bytes" );
00507 bytesAdded += addKeyword_num( "NAXIS2" , naxis2, "number of rows in table" );
00508 bytesAdded += addKeyword_num( "PCOUNT" , pcount, "size of special data area" );
00509 bytesAdded += addKeyword_num( "GCOUNT" , gcount, "one data group (required keyword)" );
00510 bytesAdded += addKeyword_num( "TFIELDS" , numFields, "number of fields in each row" );
00511
00512
00513
00514
00515 bytesAdded += addKeyword_txt( "TTYPE1" , "'PARAMVAL'", "label for field 1", 0 );
00516
00517 sprintf( theValue2temp, "%ld%s", nintparm, "E" );
00518 sprintf( theValue2, "%s%-8s%s", "'",theValue2temp,"'" );
00519 bytesAdded += addKeyword_txt( "TFORM1" , theValue2, "data format of the field: 4-byte REAL", 0 );
00520
00521
00522
00523
00524 bytesAdded += addKeyword_txt( "TTYPE2" , "'INTPSPEC'", "label for field 2", 0 );
00525
00526 sprintf( theValue_temp, "%ld%s", numEnergies, "E" );
00527 sprintf( theValue, "%s%-8s%s", "'",theValue_temp,"'" );
00528 bytesAdded += addKeyword_txt( "TFORM2" , theValue, "data format of the field: 4-byte REAL", 0 );
00529 bytesAdded += addKeyword_txt( "TUNIT2" , ModelUnits[lgAddModel], "physical unit of field", 0 );
00530
00531
00532
00533
00534 for( i=1; i<=naddparm; i++ )
00535 {
00536 sprintf( theKeyword1, "%s%ld", "TTYPE", i+2 );
00537 sprintf( theKeyword2, "%s%ld", "TFORM", i+2 );
00538 sprintf( theKeyword3, "%s%ld", "TUNIT", i+2 );
00539
00540 sprintf( theValue1, "%s%2.2ld%s", "'ADDSP", i, "'" );
00541
00542 sprintf( theValue2temp, "%ld%s", numEnergies, "E" );
00543 sprintf( theValue2, "%s%-8s%s", "'",theValue2temp,"'" );
00544
00545 sprintf( theComment1, "%s%ld", "label for field ", i+2 );
00546
00547 bytesAdded += addKeyword_txt( theKeyword1 , theValue1, theComment1, 0 );
00548 bytesAdded += addKeyword_txt( theKeyword2 , theValue2, "data format of the field: 4-byte REAL", 0 );
00549 bytesAdded += addKeyword_txt( theKeyword3 , ModelUnits[lgAddModel], "physical unit of field", 0 );
00550 }
00551
00552 bytesAdded += addKeyword_txt( "EXTNAME" , "'SPECTRA '", "name of this binary table extension", 0 );
00553 bytesAdded += addKeyword_txt( "HDUCLASS", "'OGIP '", "Format conforms to OGIP/GSFC conventions", 0 );
00554 bytesAdded += addKeyword_txt( "HDUCLAS1", "'XSPEC TABLE MODEL'","model spectra for XSPEC", 0 );
00555 bytesAdded += addKeyword_txt( "HDUCLAS2", "'MODEL SPECTRA'", "Extension containing model spectra", 0 );
00556 bytesAdded += addKeyword_txt( "HDUVERS" , "'1.0.0 '", "Version of format (OGIP memo OGIP-92-001)", 0 );
00557
00558 bytesAdded += fprintf(ioFITS_OUTPUT, "%-80s", "END" );
00559
00560 ASSERT( bytesAdded%LINESIZE == 0 );
00561
00562 while( bytesAdded%RECORDSIZE > 0 )
00563 {
00564 bytesAdded += fprintf(ioFITS_OUTPUT, "%-1s", " " );
00565 }
00566 return;
00567 }
00568
00569 STATIC void punchFITS_SpectraData( realnum **interpParameters, multi_arr<realnum,3>& theSpectrum, int option,
00570 long totNumModels, long numEnergies, long nintparm, long naddparm )
00571 {
00572 long i;
00573 long naxis2 = totNumModels;
00574
00575 DEBUG_ENTRY( "punchFITS_SpectraData()" );
00576
00577 ASSERT( nintparm + naddparm <= LIMPAR );
00578
00579
00580 for( i=0; i<naxis2; i++ )
00581 {
00582
00583 #if !defined(_BIG_ENDIAN)
00584 for( long j = 0; j<numEnergies; j++ )
00585 {
00586 ByteSwap5( theSpectrum[option][i][j] );
00587 }
00588
00589 for( long j = 0; j<nintparm; j++ )
00590 {
00591 ByteSwap5( interpParameters[i][j] );
00592 }
00593 #endif
00594
00595
00596 bytesAdded += (long)fwrite( interpParameters[i], 1, (unsigned)nintparm*sizeof(realnum), ioFITS_OUTPUT );
00597
00598 bytesAdded += (long)fwrite( &theSpectrum[option][i][0], 1, (unsigned)numEnergies*sizeof(realnum), ioFITS_OUTPUT );
00599
00600 #if !defined(_BIG_ENDIAN)
00601
00602 for( long j = 0; j<numEnergies; j++ )
00603 {
00604 ByteSwap5( theSpectrum[option][i][j] );
00605 }
00606
00607 for( long j = 0; j<nintparm; j++ )
00608 {
00609 ByteSwap5( interpParameters[i][j] );
00610 }
00611 #endif
00612
00613
00614 if( naddparm > 0 )
00615 {
00616
00617
00618
00619
00620 fprintf( ioQQQ, " Additional parameters not currently supported.\n" );
00621 cdEXIT( EXIT_FAILURE );
00622 }
00623 }
00624
00625 while( bytesAdded%RECORDSIZE > 0 )
00626 {
00627 int tempInt = 0;
00628 bytesAdded += (long)fwrite( &tempInt, 1, 1, ioFITS_OUTPUT );
00629 }
00630 return;
00631 }
00632
00633 STATIC void punchFITS_GenericHeader( long numEnergies )
00634 {
00635 long numFields = 2;
00636 long naxis, naxis1, naxis2;
00637
00638 DEBUG_ENTRY( "punchFITS_GenericHeader()" );
00639
00640
00641 ASSERT( bytesAdded%RECORDSIZE == 0 );
00642
00643 naxis = 2;
00644 naxis1 = numFields*(long)sizeof(realnum);
00645 naxis2 = numEnergies;
00646
00647 bytesAdded += addKeyword_txt( "XTENSION", "'BINTABLE'", "binary table extension", 0 );
00648 bytesAdded += addKeyword_num( "BITPIX" , bitpix, "8-bit bytes" );
00649 bytesAdded += addKeyword_num( "NAXIS" , naxis, "2-dimensional binary table" );
00650 bytesAdded += addKeyword_num( "NAXIS1" , naxis1, "width of table in bytes" );
00651 bytesAdded += addKeyword_num( "NAXIS2" , naxis2, "number of rows in table" );
00652 bytesAdded += addKeyword_num( "PCOUNT" , pcount, "size of special data area" );
00653 bytesAdded += addKeyword_num( "GCOUNT" , gcount, "one data group (required keyword)" );
00654 bytesAdded += addKeyword_num( "TFIELDS" , numFields, "number of fields in each row" );
00655 bytesAdded += addKeyword_txt( "TTYPE1" , "'ENERGY '", "label for field 1", 0 );
00656 bytesAdded += addKeyword_txt( "TFORM1" , "'E '", "data format of the field: 4-byte REAL", 0 );
00657 bytesAdded += addKeyword_txt( "TTYPE2" , "'TRN_SPEC'", "label for field 2", 0 );
00658 bytesAdded += addKeyword_txt( "TFORM2" , "'E '", "data format of the field: 4-byte REAL", 0 );
00659 bytesAdded += addKeyword_txt( "EXTNAME" , "'SPECTRA '", "name of this binary table extension", 0 );
00660 bytesAdded += addKeyword_txt( "HDUCLASS", "'OGIP '", "Format conforms to OGIP/GSFC conventions", 0 );
00661 bytesAdded += addKeyword_txt( "HDUCLAS1", "'XSPEC TABLE MODEL'","model spectra for XSPEC", 0 );
00662 bytesAdded += addKeyword_txt( "HDUCLAS2", "'ENERGIES'", "Extension containing energy bin info", 0 );
00663 bytesAdded += addKeyword_txt( "HDUVERS" , "'1.0.0 '", "Version of format (OGIP memo OGIP-92-001)", 0 );
00664
00665 bytesAdded += fprintf(ioFITS_OUTPUT, "%-80s", "END" );
00666
00667 ASSERT( bytesAdded%LINESIZE == 0 );
00668
00669 while( bytesAdded%RECORDSIZE > 0 )
00670 {
00671 bytesAdded += fprintf(ioFITS_OUTPUT, "%-1s", " " );
00672 }
00673 return;
00674 }
00675
00676 STATIC void punchFITS_GenericData( long numEnergies, long ipLoEnergy, long ipHiEnergy )
00677 {
00678 long i;
00679
00680 DEBUG_ENTRY( "punchFITS_GenericData()" );
00681
00682 realnum *TransmittedSpectrum;
00683
00684 TransmittedSpectrum = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(numEnergies) );
00685
00686 cdSPEC2( 8, numEnergies, ipLoEnergy, ipHiEnergy, TransmittedSpectrum );
00687
00688
00689 for( i=0; i<numEnergies; i++ )
00690 {
00691 realnum Energy;
00692 Energy = rfield.AnuOrg[i];
00693
00694 #if !defined(_BIG_ENDIAN)
00695 ByteSwap5(Energy);
00696 ByteSwap5(TransmittedSpectrum[i]);
00697 #endif
00698
00699 bytesAdded += (long)fwrite( &Energy , 1, sizeof(realnum), ioFITS_OUTPUT );
00700 bytesAdded += (long)fwrite( &TransmittedSpectrum[i],1, sizeof(realnum), ioFITS_OUTPUT );
00701 }
00702
00703 while( bytesAdded%RECORDSIZE > 0 )
00704 {
00705 int tempInt = 0;
00706 bytesAdded += (long)fwrite( &tempInt, 1, 1, ioFITS_OUTPUT );
00707 }
00708
00709 free( TransmittedSpectrum );
00710 return;
00711 }
00712
00713 STATIC void writeCloudyDetails( void )
00714 {
00715 char timeString[30]="";
00716 char tempString[70];
00717 time_t now;
00718 long i, j, k;
00719
00720
00721
00722 now = time(NULL);
00723 if( prt.lgPrintTime )
00724 {
00725
00726
00727 strcpy( timeString , ctime(&now) );
00728 }
00729
00730
00731 for( i=0; i<30; i++ )
00732 {
00733 if( timeString[i] == '\n' )
00734 {
00735 timeString[i] = ' ';
00736 }
00737 }
00738
00739 strcpy( tempString, "Generated by Cloudy " );
00740
00741 strncat( tempString, t_version::Inst().chVersion, sizeof(tempString)-strlen(tempString) );
00742 bytesAdded += addComment( tempString );
00743 bytesAdded += addComment( t_version::Inst().chInfo );
00744 strcpy( tempString, "--- " );
00745 strcat( tempString, timeString );
00746 bytesAdded += addComment( tempString );
00747 bytesAdded += addComment( "Input string was as follows: " );
00748
00749 for( i=0; i<=input.nSave; i++ )
00750 {
00751 char firstLine[70], extraLine[65];
00752
00753 for( j=0; j<INPUT_LINE_LENGTH; j++ )
00754 {
00755 if( input.chCardSav[i][j] =='\0' )
00756 break;
00757 }
00758
00759 ASSERT( j < 200 );
00760 for( k=0; k< MIN2(69, j); k++ )
00761 {
00762 firstLine[k] = input.chCardSav[i][k];
00763 }
00764 firstLine[k] = '\0';
00765 bytesAdded += addComment( firstLine );
00766 if( j >= 69 )
00767 {
00768 for( k=69; k< 133; k++ )
00769 {
00770 extraLine[k-69] = input.chCardSav[i][k];
00771 }
00772
00773 extraLine[64] = '\0';
00774 strcpy( tempString, "more " );
00775 strcat( tempString, extraLine );
00776 bytesAdded += addComment( tempString );
00777 if( j >= 133 )
00778 {
00779 for( k=133; k< 197; k++ )
00780 {
00781 extraLine[k-133] = input.chCardSav[i][k];
00782 }
00783 extraLine[64] = '\0';
00784 strcpy( tempString, "more " );
00785 strcat( tempString, extraLine );
00786 bytesAdded += addComment( tempString );
00787 }
00788 }
00789 }
00790
00791 return;
00792 }
00793
00794 STATIC long addKeyword_txt( const char *theKeyword, const void *theValue, const char *theComment, long Str_Or_Log )
00795 {
00796 long numberOfBytesWritten = 0;
00797
00798 DEBUG_ENTRY( "addKeyword_txt()" );
00799
00800
00801 if( Str_Or_Log == 0 )
00802 {
00803 numberOfBytesWritten = fprintf(ioFITS_OUTPUT, "%-8s%-2s%-20s%3s%-47s",
00804 theKeyword,
00805 "= ",
00806 (char *)theValue,
00807 " / ",
00808 theComment );
00809 }
00810 else
00811 {
00812 ASSERT( Str_Or_Log == 1 );
00813 numberOfBytesWritten = fprintf(ioFITS_OUTPUT, "%-8s%-2s%20s%3s%-47s",
00814 theKeyword,
00815 "= ",
00816 (char *)theValue,
00817 " / ",
00818 theComment );
00819 }
00820
00821 ASSERT( numberOfBytesWritten%LINESIZE == 0 );
00822 return numberOfBytesWritten;
00823 }
00824
00825 STATIC long addKeyword_num( const char *theKeyword, long theValue, const char *theComment)
00826 {
00827 long numberOfBytesWritten = 0;
00828
00829 DEBUG_ENTRY( "addKeyword_num()" );
00830
00831 numberOfBytesWritten = fprintf(ioFITS_OUTPUT, "%-8s%-2s%20ld%3s%-47s",
00832 theKeyword,
00833 "= ",
00834 theValue,
00835 " / ",
00836 theComment );
00837
00838 ASSERT( numberOfBytesWritten%LINESIZE == 0 );
00839 return numberOfBytesWritten;
00840 }
00841
00842 long addComment( const char *CommentToAdd )
00843 {
00844 long i, numberOfBytesWritten = 0;
00845 char tempString[70] = " ";
00846
00847 DEBUG_ENTRY( "addComment()" );
00848
00849 strncpy( &tempString[0], CommentToAdd, 69 );
00850 ASSERT( (int)strlen( tempString ) <= 70 );
00851
00852
00853 for( i=0; i<69; i++ )
00854 {
00855 if( tempString[i] == '\t' )
00856 {
00857 tempString[i] = ' ';
00858 }
00859 }
00860
00861 numberOfBytesWritten = fprintf(ioFITS_OUTPUT, "COMMENT %-70s", tempString );
00862
00863 ASSERT( numberOfBytesWritten%LINESIZE == 0 );
00864 return numberOfBytesWritten;
00865 }