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