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