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 /*ParseTable parse the table read command */ 00004 /*lines_table invoked by table lines command, check if we can find all lines in a given list */ 00005 /*read_hm05 read in the data files and interpolate the continuum to 00006 * the correct redshift */ 00007 #include "cddefines.h" 00008 #include "cddrive.h" 00009 #include "physconst.h" 00010 #include "optimize.h" 00011 #include "rfield.h" 00012 #include "trace.h" 00013 #include "radius.h" 00014 #include "input.h" 00015 #include "stars.h" 00016 #include "lines.h" 00017 #include "prt.h" 00018 #include "parse.h" 00019 /* HP cc cannot compile this routine with any optimization */ 00020 #if defined(__HP_aCC) 00021 # pragma OPT_LEVEL 1 00022 #endif 00023 00024 /* these will become the label and wl for a possible list of lines, 00025 * obtained when tables lines used */ 00026 static char **chLabel; 00027 static realnum *wl; 00028 static long int nLINE_TABLE = 0; 00029 static char chLINE_LIST[FILENAME_PATH_LENGTH]; 00030 00031 /* tables of various built in continue */ 00032 /* crab nebula */ 00033 #define NCRAB 10 00034 static double tnucrb[NCRAB], 00035 fnucrb[NCRAB]; 00036 00037 /* Bob Rubin's corrected theta 1 Ori C continuum */ 00038 #define NRUBIN 56 00039 double tnurbn[NRUBIN] = {1.05E-08,1.05E-07,1.05E-06,1.04E-05,1.00E-04,1.00E-03,1.00E-02,3.01E-02,1.00E-01, 00040 1.50E-01,2.50E-01,4.01E-01,6.01E-01,9.8E-01,9.96E-01,1.00E+00,1.02445,1.07266,1.12563,1.18411,1.23881, 00041 1.29328,1.35881,1.42463,1.48981,1.55326,1.6166,1.68845,1.76698,1.8019,1.808,1.84567,1.9317,2.04891,2.14533, 00042 2.19702,2.27941,2.37438,2.43137,2.49798,2.56113,2.59762,2.66597,2.80543,2.95069,3.02911,3.11182,3.22178, 00043 3.3155,3.42768,3.50678,3.56157,3.61811,3.75211,3.89643,4.}, 00044 fnurbn[NRUBIN] = {1.56E-20,1.55E-17,1.54E-14,1.53E-11,1.35E-08,1.34E-05,1.35E-02,3.62E-01,1.27E+01, 00045 3.90E+01,1.48E+02,4.52E+02,1.02E+03,2.27E+03,8.69E+02,8.04E+02,6.58E+02,6.13E+02,6.49E+02,6.06E+02, 00046 6.30E+02,5.53E+02,5.55E+02,5.24E+02,4.86E+02,4.49E+02,4.42E+02,3.82E+02,3.91E+02,2.91E+02,2.61E+02, 00047 1.32E+02,1.20E+02,1.16E+02,9.56E+01,9.94E+01,9.10E+01,4.85E+01,7.53E+01,9.53E+01,8.52E+01,5.76E+01, 00048 6.72E+01,5.20E+01,8.09E+00,3.75E+00,5.57E+00,3.80E+00,2.73E+00,2.22E+00,3.16E-01,1.26E-01,1.39E-01, 00049 6.15E-02,3.22E-02,7.98E-03}; 00050 00051 /* stores numbers for table cooling flow */ 00052 #define NCFL 40 00053 static double cfle[NCFL], 00054 cflf[NCFL]; 00055 00056 /* Brad Peterson's AKN 120 continuum */ 00057 #define NKN120 11 00058 static double tnuakn[NKN120], 00059 fnuakn[NKN120]; 00060 00061 /* Black's ISM continuum, with He hole filled in */ 00062 #define NISM 23 00063 static double tnuism[NISM], 00064 fnuism[NISM]; 00065 00066 /* z=2 background, 00067 * >>refer continuum background Haardt, Francesco, & Madau, Piero, 1996, 00068 * >>refercon ApJ, 461, 20 */ 00069 #define NHM96 14 00070 /* log energy in Ryd */ 00071 static double tnuHM96[NHM96]={-8,-1.722735683,-0.351545683,-0.222905683,-0.133385683, 00072 /* changeg these two energies to prevent degeneracy */ 00073 -0.127655683,-0.004575683,0.297544317,0.476753,0.476756,0.588704317, 00074 0.661374317,1.500814317,2.245164317}; 00075 /*-0.127655683,-0.004575683,0.297544317,0.476754317,0.476754317,0.588704317,*/ 00076 /*log J in the units of (erg cm^{-2} s^{-1} Hz^{-1} sr^{-1})*/ 00077 static double fnuHM96[NHM96]={-32.53342863,-19.9789,-20.4204,-20.4443,-20.5756,-20.7546, 00078 -21.2796,-21.6256,-21.8404,-21.4823,-22.2102,-22.9263,-23.32,-24.2865}; 00079 00080 /* Mathews and Ferland generic AGN continuum */ 00081 #define NAGN 8 00082 static double tnuagn[NAGN], 00083 tslagn[NAGN]; 00084 00085 /* table Draine ISM continuum */ 00086 #define NDRAINE 15 00087 static double tnudrn[NDRAINE] , tsldrn[NDRAINE]; 00088 00089 /* routine that stores values for above vectors */ 00090 STATIC void ZeroContin(void); 00091 00092 /*read_hm05 read in the data files and interpolate the Haardt & Madau continuum to 00093 * the correct redshift */ 00094 STATIC void read_hm05( const char chFile[] , double **tnuHM , double **fnuHM , 00095 long int *nhm , double redshift ) 00096 { 00097 00098 FILE *ioFILE; 00099 double *table_redshifts = NULL, 00100 *table_wl = NULL , 00101 **table_fn=NULL, 00102 frac_hi; 00103 char chLine[1000]; 00104 long int nRedshift , i , n , nz , ipLo , ipHi; 00105 bool lgEOL; 00106 00107 DEBUG_ENTRY( "read_hm05()" ); 00108 00109 ioFILE = open_data( chFile, "r", AS_LOCAL_DATA ); 00110 00111 /* this will be the number of continuum points in their table */ 00112 *nhm = 0; 00113 /* the number of redshifts in their table */ 00114 nRedshift = -1; 00115 /* first read past comments and get magic number */ 00116 while( read_whole_line( chLine , (int)sizeof(chLine) , ioFILE ) != NULL ) 00117 { 00118 /* we want to count the lines that do not start with # 00119 * since these contain data */ 00120 if( chLine[0] != '#') 00121 { 00122 ++*nhm; 00123 /* check magic number if first valid line */ 00124 if( *nhm==1 ) 00125 { 00126 long mag_read; 00127 /* this is the magic number - read it in */ 00128 i = 1; 00129 mag_read = (long)FFmtRead(chLine,&i,(int)sizeof(chLine),&lgEOL); 00130 if( mag_read != 50808 ) 00131 { 00132 fprintf(ioQQQ, 00133 " Magic number in Haardt & Madau file is not correct, I expected 50808 and found %li\n", 00134 mag_read ); 00135 cdEXIT(EXIT_FAILURE); 00136 } 00137 } 00138 /* second valid line count number of redshifts */ 00139 else if( *nhm==2 ) 00140 { 00141 double a; 00142 i = 2; 00143 lgEOL = false; 00144 nRedshift = 0; 00145 while( !lgEOL ) 00146 { 00147 ++nRedshift; 00148 a = FFmtRead(chLine,&i,(int)sizeof(chLine),&lgEOL); 00149 ASSERT( a >= 0. ); 00150 } 00151 /* have over counted by one - back up one */ 00152 --nRedshift; 00153 /* highest redshift data missing in file i received */ 00154 --nRedshift; 00155 /*fprintf(ioQQQ," number of z is %li\n", nRedshift);*/ 00156 /* malloc some space */ 00157 table_redshifts = (double*)MALLOC( (size_t)nRedshift*sizeof(double) ); 00158 /* now read in the redshifts */ 00159 i = 2; 00160 for( n=0; n<nRedshift; ++n ) 00161 { 00162 table_redshifts[n] = FFmtRead(chLine,&i,(int)sizeof(chLine),&lgEOL); 00163 /*fprintf(ioQQQ,"DEBUG %li z %.3f\n", n , table_redshifts[n] );*/ 00164 } 00165 if( lgEOL ) 00166 TotalInsanity(); 00167 } 00168 } 00169 } 00170 00171 /* malloc space for the arrays first wavelength array */ 00172 table_wl = (double*)MALLOC( (size_t)*nhm*sizeof(double) ); 00173 /* the spectrum array table_fn[z][nu] */ 00174 table_fn = (double**)MALLOC( (size_t)nRedshift*sizeof(double*) ); 00175 for(n=0; n<nRedshift; ++n ) 00176 { 00177 table_fn[n] = (double*)MALLOC( (size_t)(*nhm)*sizeof(double) ); 00178 } 00179 00180 /* rewind the file so we can read it a second time*/ 00181 if( fseek( ioFILE , 0 , SEEK_SET ) != 0 ) 00182 { 00183 fprintf( ioQQQ, " read_hm05 could not rewind HM05 date file.\n"); 00184 cdEXIT(EXIT_FAILURE); 00185 } 00186 n = 0; 00187 *nhm = 0; 00188 while( read_whole_line( chLine , (int)sizeof(chLine) , ioFILE ) != NULL ) 00189 { 00190 /* we want to count the lines that do not start with # 00191 * since these contain data */ 00192 if( chLine[0] != '#') 00193 { 00194 /* this is number of non-comment lines - will skip magic number 00195 * and line giving redshift */ 00196 ++n; 00197 if( n>2 ) 00198 { 00199 i = 1; 00200 table_wl[*nhm] = FFmtRead(chLine,&i,(int)sizeof(chLine),&lgEOL); 00201 if( lgEOL ) 00202 TotalInsanity(); 00203 for( nz=0; nz<nRedshift; ++nz ) 00204 { 00205 /* >>chng 07 feb 18, PvH change from last branck to first */ 00206 table_fn[nz][*nhm] = FFmtRead(chLine,&i,(int)sizeof(chLine),&lgEOL); 00207 } 00208 ++*nhm; 00209 } 00210 } 00211 } 00212 00213 fclose( ioFILE ); 00214 00215 { 00216 /* change following to true to print their original table */ 00217 enum {DEBUG_LOC=false}; 00218 if( DEBUG_LOC) 00219 { 00220 fprintf(ioQQQ,"wavelength/z"); 00221 for(nz=0; nz<nRedshift; ++nz ) 00222 fprintf(ioQQQ,"\t%.3f", table_redshifts[nz] ); 00223 fprintf(ioQQQ,"\n"); 00224 for( i=0; i<*nhm; ++i ) 00225 { 00226 fprintf(ioQQQ,"%.3e", table_wl[i] ); 00227 for( nz=0; nz<nRedshift; ++nz ) 00228 fprintf(ioQQQ,"\t%.3e", table_fn[nz][i] ); 00229 fprintf(ioQQQ,"\n"); 00230 } 00231 } 00232 } 00233 00234 /* this is just to shut up the lint and must not be ASSERT */ 00235 assert( table_redshifts!=NULL ); 00236 00237 /* first check that desired redshift is within range of table */ 00238 if( redshift < table_redshifts[0] || 00239 redshift > table_redshifts[nRedshift-1] ) 00240 { 00241 fprintf(ioQQQ," The redshift requested on table HM05 is %g but is not within the range of the table, which goes from %g to %g.\n", 00242 redshift, table_redshifts[0] , table_redshifts[nRedshift-1] ); 00243 fprintf(ioQQQ," Sorry.\n"); 00244 cdEXIT(EXIT_FAILURE); 00245 } 00246 00247 ipLo = -1; 00248 ipHi = -1; 00249 /* find which redshift bin we need */ 00250 for( nz=0; nz<nRedshift-1; ++nz ) 00251 { 00252 if( redshift >= table_redshifts[nz] && 00253 redshift <= table_redshifts[nz+1] ) 00254 { 00255 ipLo = nz; 00256 ipHi = nz+1; 00257 break; 00258 } 00259 } 00260 ASSERT( ipLo>=0 && ipHi >=0 ); 00261 00262 /* make sure that the wavelengths are in increasing order - they were not in the 00263 * original data table, but had repeated values near important ionization edges */ 00264 for( i=0; i<*nhm-1; ++i ) 00265 { 00266 if( table_wl[i]>=table_wl[i+1] ) 00267 { 00268 fprintf(ioQQQ," The wavelengths in the table HM05 data table are not in increasing order. This is required.\n"); 00269 fprintf(ioQQQ," Sorry.\n"); 00270 cdEXIT(EXIT_FAILURE); 00271 } 00272 } 00273 00274 /* malloc the space for the returned arrays */ 00275 *fnuHM = (double *)MALLOC( (size_t)(*nhm)*sizeof(double ) ); 00276 *tnuHM = (double *)MALLOC( (size_t)(*nhm)*sizeof(double ) ); 00277 00278 /* now fill in the arrays with the interpolated table 00279 * we will interpolate linearly in redshift 00280 * in general the redshift will be between the tabulated redshift */ 00281 frac_hi = (redshift-table_redshifts[ipLo]) / (table_redshifts[ipHi]-table_redshifts[ipLo]); 00282 for( i=0; i<*nhm; ++i ) 00283 { 00284 /* we have wavelengths in angstroms and want log Ryd 00285 * original table was in decreasing energy order, must also 00286 * flip it around since need increasing energy order */ 00287 (*tnuHM)[*nhm-1-i] = RYDLAM / table_wl[i]; 00288 /* original table in correct units so no conversion needed */ 00289 (*fnuHM)[*nhm-1-i] = table_fn[ipLo][i]*(1.-frac_hi) + table_fn[ipHi][i]*frac_hi; 00290 /*fprintf(ioQQQ,"DEBUG1\t%.3e\t%.3e\n",(*tnuHM)[*nhm-1-i] , (*fnuHM)[*nhm-1-i] );*/ 00291 } 00292 00293 for( n=0; n<nRedshift; ++n ) 00294 free( table_fn[n] ); 00295 free( table_fn ); 00296 free( table_wl ); 00297 free( table_redshifts ); 00298 return; 00299 } 00300 00301 void ParseTable(long int *nqh, 00302 char *chCard, 00303 realnum *ar1) 00304 { 00305 char chFile[FILENAME_PATH_LENGTH_2]; /*file name for table read */ 00306 00307 IntMode imode = IM_ILLEGAL_MODE; 00308 bool lgEOL, 00309 lgHit, 00310 lgLogSet; 00311 static bool lgCalled=false; 00312 00313 long int i, 00314 j, 00315 k, 00316 ncont, 00317 nstar, 00318 ipNORM; 00319 00320 double alpha, 00321 brakmm, 00322 brakxr, 00323 ConBreak, 00324 fac, 00325 scale, 00326 slopir, 00327 slopxr; 00328 00329 bool lgNoContinuum = false, 00330 lgQuoteFound; 00331 00332 DEBUG_ENTRY( "ParseTable()" ); 00333 00334 /* if first call then set up values for table */ 00335 if( !lgCalled ) 00336 { 00337 ZeroContin(); 00338 lgCalled = true; 00339 } 00340 00341 if( rfield.nspec >= LIMSPC ) 00342 { 00343 fprintf( ioQQQ, " %ld is too many spectra entered. Increase LIMSPC\n Sorry.\n", 00344 rfield.nspec ); 00345 cdEXIT(EXIT_FAILURE); 00346 } 00347 00348 /* three commands, tables line, read, and star, have quotes on the 00349 * lines giving file names. must get quotes first so that filename 00350 * does not confuse parser */ 00351 lgQuoteFound = true; 00352 if( GetQuote( chFile , chCard , false ) ) 00353 lgQuoteFound = false; 00354 00355 /* set flag telling interpolate */ 00356 strcpy( rfield.chSpType[rfield.nspec], "INTER" ); 00357 /* >>chng 06 jul 16, fix memory leak when optimizing, PvH */ 00358 if( !rfield.lgContMalloc[rfield.nspec] ) 00359 { 00360 rfield.tNuRyd[rfield.nspec] = (realnum*)CALLOC( (size_t)NCELL , sizeof(realnum) ); 00361 rfield.tslop[rfield.nspec] = (realnum*)CALLOC( (size_t)NCELL , sizeof(realnum) ); 00362 rfield.tFluxLog[rfield.nspec] = (realnum*)CALLOC( (size_t)NCELL , sizeof(realnum) ); 00363 rfield.lgContMalloc[rfield.nspec] = true; 00364 } 00365 00366 /* NB when adding more keys also change the comment at the end */ 00367 if( nMatch(" AGN",chCard) ) 00368 { 00369 /* do Mathews and Ferland generic AGN continuum */ 00370 ASSERT( NAGN < NCELL); 00371 for( i=0; i < NAGN; i++ ) 00372 { 00373 rfield.tNuRyd[rfield.nspec][i] = (realnum)tnuagn[i]; 00374 rfield.tslop[rfield.nspec][i] = (realnum)tslagn[i]; 00375 } 00376 ncont = NAGN - 1; 00377 00378 /* optional keyword break, to adjust IR cutoff */ 00379 if( nMatch("BREA",chCard) ) 00380 { 00381 i = 4; 00382 ConBreak = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00383 00384 if( lgEOL ) 00385 { 00386 /* no break, set to low energy limit of code */ 00387 if( nMatch(" NO ",chCard) ) 00388 { 00389 ConBreak = rfield.emm*1.01f; 00390 } 00391 else 00392 { 00393 fprintf( ioQQQ, " There must be a number for the break.\n Sorry.\n" ); 00394 cdEXIT(EXIT_FAILURE); 00395 } 00396 } 00397 00398 if( ConBreak == 0. ) 00399 { 00400 fprintf( ioQQQ, " The break must be greater than 0.2 Ryd.\n Sorry.\n" ); 00401 cdEXIT(EXIT_FAILURE); 00402 } 00403 00404 if( nMatch("MICR",chCard) ) 00405 { 00406 /* optional keyword, ``microns'', convert to Rydbergs */ 00407 ConBreak = 0.0912/ConBreak; 00408 } 00409 00410 if( ConBreak < 0. ) 00411 { 00412 /* option to enter break as LOG10 */ 00413 ConBreak = pow(10.,ConBreak); 00414 } 00415 00416 else if( ConBreak == 0. ) 00417 { 00418 fprintf( ioQQQ, " An energy of 0 is not allowed.\n Sorry.\n" ); 00419 cdEXIT(EXIT_FAILURE); 00420 } 00421 00422 if( ConBreak >= rfield.tNuRyd[rfield.nspec][2] ) 00423 { 00424 fprintf( ioQQQ, " The energy of the break cannot be greater than%10.2e Ryd.\n Sorry.\n", 00425 rfield.tNuRyd[rfield.nspec][2] ); 00426 cdEXIT(EXIT_FAILURE); 00427 } 00428 00429 else if( ConBreak <= rfield.tNuRyd[rfield.nspec][0] ) 00430 { 00431 fprintf( ioQQQ, " The energy of the break cannot be less than%10.2e Ryd.\n Sorry.\n", 00432 rfield.tNuRyd[rfield.nspec][0] ); 00433 cdEXIT(EXIT_FAILURE); 00434 } 00435 00436 rfield.tNuRyd[rfield.nspec][1] = (realnum)ConBreak; 00437 00438 rfield.tslop[rfield.nspec][1] = 00439 (realnum)(rfield.tslop[rfield.nspec][2] + 00440 log10(rfield.tNuRyd[rfield.nspec][2]/rfield.tNuRyd[rfield.nspec][1])); 00441 00442 rfield.tslop[rfield.nspec][0] = 00443 (realnum)(rfield.tslop[rfield.nspec][1] - 00444 2.5*log10(rfield.tNuRyd[rfield.nspec][1]/rfield.tNuRyd[rfield.nspec][0])); 00445 } 00446 } 00447 00448 else if( nMatch("AKN1",chCard) ) 00449 { 00450 /* AKN120 continuum from Brad Peterson */ 00451 ASSERT( NKN120 < NCELL ); 00452 for( i=0; i < NKN120; i++ ) 00453 { 00454 rfield.tNuRyd[rfield.nspec][i] = (realnum)tnuakn[i]; 00455 rfield.tslop[rfield.nspec][i] = (realnum)log10(fnuakn[i]); 00456 } 00457 ncont = NKN120 - 1; 00458 } 00459 00460 else if( nMatch("CRAB",chCard) ) 00461 { 00462 ASSERT( NCRAB < NCELL ); 00463 /* Crab Nebula continuum from Davidson and Fesen 1985 Ann Rev article */ 00464 for( i=0; i < NCRAB; i++ ) 00465 { 00466 rfield.tNuRyd[rfield.nspec][i] = (realnum)tnucrb[i]; 00467 rfield.tslop[rfield.nspec][i] = (realnum)log10(fnucrb[i]); 00468 } 00469 ncont = NCRAB - 1; 00470 } 00471 00472 else if( nMatch("COOL",chCard) ) 00473 { 00474 ASSERT( NCFL < NCELL ); 00475 /* cooling flow from Johnstone et al. 1992, MN 255, 431. */ 00476 for( i=0; i < NCFL; i++ ) 00477 { 00478 rfield.tNuRyd[rfield.nspec][i] = (realnum)cfle[i]; 00479 rfield.tslop[rfield.nspec][i] = (realnum)cflf[i]; 00480 } 00481 ncont = NCFL - 1; 00482 } 00483 00484 else if( (i=nMatch("HM96",chCard))>0 ) 00485 { 00486 /* this is the old Haardt & Madau continuum, one set of points 00487 * with only the quasars 00488 * this command does not include the CMB - do that separately with the CMB command */ 00489 /* set flag telling interpolate */ 00490 strcpy( rfield.chSpType[rfield.nspec], "INTER" ); 00491 00492 ASSERT( NHM96 < NCELL ); 00493 /* z=2 background, 00494 * >>refer continuum background Haardt, Francesco, & Madau, Piero, 1996, ApJ, 461, 20 */ 00495 for( j=0; j < NHM96; j++ ) 00496 { 00497 /* frequency was stored as log of ryd */ 00498 rfield.tNuRyd[rfield.nspec][j] = (realnum)pow( 10. , tnuHM96[j] ); 00499 rfield.tslop[rfield.nspec][j] = (realnum)fnuHM96[j]; 00500 } 00501 ncont = NHM96 - 1; 00502 00503 /* optional scale factor to change default intensity from their value 00504 * assumed to be log if negative, and linear otherwise 00505 * increment i to move past the 96 in the keyword */ 00506 i += 4; 00507 scale = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00508 if( scale > 0. ) 00509 scale = log10(scale); 00510 00511 /* this also sets continuum intensity*/ 00512 if( *nqh >= LIMSPC ) 00513 { 00514 fprintf( ioQQQ, " %ld is too many continua entered. Increase LIMSPC\n Sorry.\n", 00515 *nqh ); 00516 cdEXIT(EXIT_FAILURE); 00517 } 00518 00519 /* check that stack of shape and luminosity specifications 00520 * is parallel, stop if not - this happens is background comes 00521 * BETWEEN another set of shape and luminosity commands */ 00522 if( rfield.nspec != *nqh ) 00523 { 00524 fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" ); 00525 fprintf( ioQQQ, " Sorry.\n" ); 00526 cdEXIT(EXIT_FAILURE); 00527 } 00528 00529 strcpy( rfield.chRSpec[*nqh], "SQCM" ); 00530 strcpy( rfield.chSpNorm[*nqh], "FLUX" ); 00531 00532 /* this is an isotropic radiation field */ 00533 rfield.lgBeamed[*nqh] = false; 00534 00535 /* this will be flux density at some frequency on the table. the numbers 00536 * are per Hz and sr so must multiply by 4 pi 00537 * [2] is not special, could have been any within array*/ 00538 rfield.range[*nqh][0] = pow(10. , tnuHM96[2] )*1.0001; 00539 00540 /* convert intensity HM96 give to current units of code */ 00541 rfield.totpow[*nqh] = (fnuHM96[2] + log10(PI4) + scale); 00542 00543 /* set radius to very large value if not already set */ 00544 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */ 00545 if( !radius.lgRadiusKnown ) 00546 { 00547 *ar1 = (realnum)radius.rdfalt; 00548 radius.Radius = pow(10.,radius.rdfalt); 00549 } 00550 ++*nqh; 00551 00552 /* set radius to very large value if not already set */ 00553 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */ 00554 if( !radius.lgRadiusKnown ) 00555 { 00556 *ar1 = (realnum)radius.rdfalt; 00557 radius.Radius = pow(10.,radius.rdfalt); 00558 } 00559 } 00560 00561 else if( nMatch("HM05",chCard) ) 00562 { 00563 double *tnuHM , *fnuHM; 00564 double redshift; 00565 long int nhm; 00566 /* the Haardt & Madau 2005 continuum, read in a table and interpolate 00567 * for any redshift, background includes both quasars and galaxies 00568 * >>refer continuum background Haardt, Francesco, & Madau, Piero, 2005, in preparation 00569 * this command does not include the CMB - do that separately with the CMB command 00570 * set flag telling interpolate */ 00571 strcpy( rfield.chSpType[rfield.nspec], "INTER" ); 00572 if( nMatch("QUAS",chCard) ) 00573 { 00574 /* quasar-only continuum */ 00575 strcpy( chFile , "haardt_madau_quasar.dat" ); 00576 } 00577 else 00578 { 00579 /* the default, quasar plus galaxy continuum */ 00580 strcpy( chFile , "haardt_madau_galaxy.dat" ); 00581 } 00582 00583 /* find the redshift - must be within bounds of table, which we 00584 * do not know yet */ 00585 i = 4; 00586 /* this must pick up the 05 in the key word */ 00587 k = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00588 if( lgEOL || k!=5 ) 00589 TotalInsanity(); 00590 redshift = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00591 if( lgEOL ) 00592 { 00593 fprintf(ioQQQ," The redshift MUST be specified on the table HM05 command. I did not find one.\n Sorry.\n"); 00594 cdEXIT(EXIT_FAILURE); 00595 } 00596 00597 /* read in the data files and interpolate the continuum to 00598 * the correct redshift */ 00599 read_hm05( chFile , &tnuHM , &fnuHM , &nhm , redshift ); 00600 /* now find a point near 1 Ryd where continuum may not be far too faint */ 00601 ipNORM = -1; 00602 for( j=0; j < nhm-1; j++ ) 00603 { 00604 /* looking for continuum frequency near 1 Ryd 00605 * does not need to be precise */ 00606 if( tnuHM[j]<=1. && 1. <= tnuHM[j+1] ) 00607 { 00608 ipNORM = j; 00609 break; 00610 } 00611 } 00612 ASSERT( ipNORM>=0 ); 00613 00614 /* optional scale factor to change default intensity from their value 00615 * assumed to be log if negative, and linear otherwise 00616 * increment i to move past the 96 in the keyword */ 00617 scale = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00618 if( scale > 0. ) 00619 scale = log10(scale); 00620 00621 /* this will be flux density at some frequency on the table. the numbers 00622 * are per Hz and sr so must multiply by 4 pi */ 00623 rfield.range[*nqh][0] = tnuHM[ipNORM]; 00624 00625 /* convert intensity HM96 give to current units of code */ 00626 rfield.totpow[*nqh] = log10(fnuHM[ipNORM]) + log10(PI4) + scale; 00627 00628 /* this also sets continuum intensity*/ 00629 if( *nqh >= LIMSPC ) 00630 { 00631 fprintf( ioQQQ, " %ld is too many continua entered. Increase LIMSPC\n Sorry.\n", 00632 *nqh ); 00633 cdEXIT(EXIT_FAILURE); 00634 } 00635 00636 /* check that stack of shape and luminosity specifications 00637 * is parallel, stop if not - this happens is background comes 00638 * BETWEEN another set of shape and luminosity commands */ 00639 if( rfield.nspec != *nqh ) 00640 { 00641 fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" ); 00642 fprintf( ioQQQ, " Sorry.\n" ); 00643 cdEXIT(EXIT_FAILURE); 00644 } 00645 00646 strcpy( rfield.chRSpec[*nqh], "SQCM" ); 00647 strcpy( rfield.chSpNorm[*nqh], "FLUX" ); 00648 00649 /* this is an isotropic radiation field */ 00650 rfield.lgBeamed[*nqh] = false; 00651 00652 /* many of their values are outside the range of a 32 bit cpu and we don't 00653 * want to fill array with zeros - so rescale to cont near 1 ryd */ 00654 scale = SDIV( fnuHM[ipNORM] ); 00655 /* their table does not extend to the low-energy limit of the code 00656 * usually does not matter since CMB will take over, but there is 00657 * an obvious gap at low, z = 0, redshift. assume Rayleigh-Jeans tail 00658 * and add a first continuum point */ 00659 rfield.tNuRyd[rfield.nspec][0] = (realnum)rfield.emm; 00660 rfield.tslop[rfield.nspec][0] = 00661 (realnum)log10(fnuHM[0]*pow((double)(rfield.emm/tnuHM[0]),2.)/scale); 00662 /*fprintf(ioQQQ,"DEBUG2\t%.3e\t%.3e\n", 00663 rfield.tNuRyd[rfield.nspec][0] , 00664 rfield.tslop[rfield.nspec][0] );*/ 00665 for( j=0; j < nhm; j++ ) 00666 { 00667 /* frequency is already Ryd */ 00668 rfield.tNuRyd[rfield.nspec][j+1] = (realnum)tnuHM[j]; 00669 rfield.tslop[rfield.nspec][j+1] = (realnum)log10(fnuHM[j]/scale); 00670 /*fprintf(ioQQQ,"DEBUG2\t%.3e\t%.3e\n", 00671 rfield.tNuRyd[rfield.nspec][j+1] , 00672 rfield.tslop[rfield.nspec][j+1] );*/ 00673 } 00674 ++nhm; 00675 ncont = nhm - 1; 00676 00677 /* now make sure they are in increasing order */ 00678 for( j=0; j < nhm-1; j++ ) 00679 ASSERT( rfield.tNuRyd[rfield.nspec][j] < rfield.tNuRyd[rfield.nspec][j+1] ); 00680 00681 /* set radius to very large value if not already set */ 00682 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */ 00683 if( !radius.lgRadiusKnown ) 00684 { 00685 *ar1 = (realnum)radius.rdfalt; 00686 radius.Radius = pow(10.,radius.rdfalt); 00687 } 00688 ++*nqh; 00689 00690 /* set radius to very large value if not already set */ 00691 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */ 00692 if( !radius.lgRadiusKnown ) 00693 { 00694 *ar1 = (realnum)radius.rdfalt; 00695 radius.Radius = pow(10.,radius.rdfalt); 00696 } 00697 free( tnuHM ); 00698 free( fnuHM ); 00699 00700 } 00701 else if( nMatch(" ISM",chCard) ) 00702 { 00703 ASSERT( NISM < NCELL ); 00704 /* local ISM radiation field from Black 1987, Interstellar Processes */ 00705 /* >>chng 04 mar 16, rm CMB from field so that it can be used at 00706 * any redshift */ 00707 rfield.tNuRyd[rfield.nspec][0] = 6.; 00708 rfield.tslop[rfield.nspec][0] = -21.21f - 6.f; 00709 for( i=6; i < NISM; i++ ) 00710 { 00711 rfield.tNuRyd[rfield.nspec][i-5] = (realnum)tnuism[i]; 00712 rfield.tslop[rfield.nspec][i-5] = (realnum)(fnuism[i] - 00713 tnuism[i]); 00714 } 00715 00716 ncont = NISM - 1 -5; 00717 i = 5; 00718 00719 /* optional scale factor to change default luminosity 00720 * from observed value 00721 * want final number to be log 00722 * assumed to be log if negative, and linear otherwise unless log option is present */ 00723 scale = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00724 if( scale > 0. && !nMatch(" LOG",chCard) ) 00725 scale = log10(scale); 00726 00727 /* this also sets continuum intensity*/ 00728 if( *nqh >= LIMSPC ) 00729 { 00730 fprintf( ioQQQ, " %4ld is too many continua entered. Increase LIMSPC\n Sorry.\n", 00731 *nqh ); 00732 cdEXIT(EXIT_FAILURE); 00733 } 00734 00735 /* check that stack of shape and luminosity specifications 00736 * is parallel, stop if not - this happens is background comes 00737 * BETWEEN another set of shape and luminosity commands */ 00738 if( rfield.nspec != *nqh ) 00739 { 00740 fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" ); 00741 fprintf( ioQQQ, " Sorry.\n" ); 00742 cdEXIT(EXIT_FAILURE); 00743 } 00744 00745 strcpy( rfield.chRSpec[*nqh], "SQCM" ); 00746 strcpy( rfield.chSpNorm[*nqh], "FLUX" ); 00747 00748 /* this is an isotropic radiation field */ 00749 rfield.lgBeamed[*nqh] = false; 00750 00751 /* this will be flux density at 1 Ryd 00752 * >>chng 96 dec 18, from 1 Ryd to H mass Rydberg 00753 * >>chng 97 jan 10, had HLevNIonRyd but not defined yet */ 00754 rfield.range[*nqh][0] = HIONPOT; 00755 00756 /* interpolated from Black 1987 */ 00757 rfield.totpow[*nqh] = (-18.517 + scale); 00758 00759 /* set radius to very large value if not already set */ 00760 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */ 00761 if( !radius.lgRadiusKnown ) 00762 { 00763 *ar1 = (realnum)radius.rdfalt; 00764 radius.Radius = pow(10.,radius.rdfalt); 00765 } 00766 ++*nqh; 00767 00768 if( optimize.lgVarOn ) 00769 { 00770 optimize.nvarxt[optimize.nparm] = 1; 00771 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE ISM LOG %f"); 00772 /* pointer to where to write */ 00773 optimize.nvfpnt[optimize.nparm] = input.nRead; 00774 /* the scale factor */ 00775 optimize.vparm[0][optimize.nparm] = (realnum)scale; 00776 optimize.vincr[optimize.nparm] = 0.2f; 00777 ++optimize.nparm; 00778 } 00779 } 00780 else if( nMatch("DRAI",chCard) ) 00781 { 00782 ASSERT( NDRAINE < NCELL ); 00783 rfield.lgMustBlockHIon = true; 00784 /* local ISM radiation field from equation 23 00785 *>>refer ISM continuum Draine & Bertoldi 1996 */ 00786 for( i=0; i < NDRAINE; i++ ) 00787 { 00788 rfield.tNuRyd[rfield.nspec][i] = (realnum)tnudrn[i]; 00789 rfield.tslop[rfield.nspec][i] = (realnum)tsldrn[i]; 00790 } 00791 00792 ncont = NDRAINE - 1; 00793 i = 5; 00794 00795 /* optional scale factor to change default luminosity 00796 * from observed value 00797 * assumed to be log if negative, and linear otherwise unless log option is present */ 00798 scale = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00799 if( scale > 0. && !nMatch(" LOG",chCard) ) 00800 scale = log10(scale); 00801 00802 /* this also sets continuum intensity*/ 00803 if( *nqh >= LIMSPC ) 00804 { 00805 fprintf( ioQQQ, " %4ld is too many continua entered. Increase LIMSPC\n Sorry.\n", 00806 *nqh ); 00807 cdEXIT(EXIT_FAILURE); 00808 } 00809 00810 /* check that stack of shape and luminosity specifications 00811 * is parallel, stop if not - this happens is background comes 00812 * BETWEEN another set of shape and luminosity commands */ 00813 if( rfield.nspec != *nqh ) 00814 { 00815 fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" ); 00816 fprintf( ioQQQ, " Sorry.\n" ); 00817 cdEXIT(EXIT_FAILURE); 00818 } 00819 00820 strcpy( rfield.chRSpec[*nqh], "SQCM" ); 00821 strcpy( rfield.chSpNorm[*nqh], "FLUX" ); 00822 00823 /* this is an isotropic radiation field */ 00824 rfield.lgBeamed[*nqh] = false; 00825 00826 /* continuum normalization given by flux density at first point, 00827 * must set energy a bit higher to make sure it is within energy bounds 00828 * that results from float arithmetic */ 00829 rfield.range[*nqh][0] = tnudrn[0]*1.01; 00830 00831 /* this is f_nu at this first point */ 00832 rfield.totpow[*nqh] = tsldrn[0] + scale; 00833 00834 if( optimize.lgVarOn ) 00835 { 00836 optimize.nvarxt[optimize.nparm] = 1; 00837 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE DRAINE LOG %f"); 00838 /* pointer to where to write */ 00839 optimize.nvfpnt[optimize.nparm] = input.nRead; 00840 /* the scale factor */ 00841 optimize.vparm[0][optimize.nparm] = (realnum)scale; 00842 optimize.vincr[optimize.nparm] = 0.2f; 00843 ++optimize.nparm; 00844 } 00845 00846 /* set radius to very large value if not already set */ 00847 /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */ 00848 if( !radius.lgRadiusKnown ) 00849 { 00850 *ar1 = (realnum)radius.rdfalt; 00851 radius.Radius = pow(10.,radius.rdfalt); 00852 } 00853 ++*nqh; 00854 } 00855 00856 else if( nMatch("LINE",chCard) ) 00857 { 00858 /* table lines command - way to check that lines within a data 00859 * file are still valid */ 00860 static bool lgLines_Malloc = false; 00861 00862 /* say that this is not a continuum command, so don't try to work with unallocated space */ 00863 lgNoContinuum = true; 00864 00865 /* this is not a continuum source - it is to read a table of lines */ 00866 if( lgLines_Malloc ) 00867 { 00868 fprintf(ioQQQ," sorry, only one table line per input stream\n"); 00869 cdEXIT(EXIT_FAILURE); 00870 } 00871 00872 /* get file name within double quotes, if not present will use default 00873 * return value of 1 indicates did not find double quotes on line */ 00874 if( lgQuoteFound ) 00875 strcpy( chLINE_LIST , chFile ); 00876 else 00877 strcpy( chLINE_LIST , "LineList_BLR.dat" ); 00878 00879 lgLines_Malloc = true; 00880 /* this flag says this is not a continuum source - nearly all table XXX 00881 * commands deal with continuum sources */ 00882 ncont = -10; 00883 00884 /* actually get the lines, and malloc the space in the arrays */ 00885 if( (nLINE_TABLE = cdGetLineList( chLINE_LIST , &chLabel , &wl) )<0 ) 00886 { 00887 /* did not find file, abort */ 00888 fprintf(ioQQQ,"\n DISASTER PROBLEM ParseTable could not find " 00889 "line list file %s\n", chLINE_LIST ); 00890 fprintf(ioQQQ," Please check the spelling of the file name and that it " 00891 "is in either the local or data directory.\n\n"); 00892 cdEXIT(EXIT_FAILURE); 00893 } 00894 } 00895 00896 else if( nMatch("POWE",chCard) ) 00897 { 00898 /* simple power law continuum between 10 micron and 50 keV 00899 * option to read in any slope for the intermediate continuum */ 00900 i = 5; 00901 alpha = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00902 00903 /* default (no number on line) is f_nu proportional nu^-1 */ 00904 if( lgEOL ) 00905 alpha = -1.; 00906 00907 /* this is low energy for code */ 00908 rfield.tNuRyd[rfield.nspec][0] =(realnum) rfield.emm; 00909 /* and the value of the flux at this point (f_nu units)*/ 00910 rfield.tslop[rfield.nspec][0] = -5.f; 00911 00912 lgLogSet = false; 00913 00914 /* option to adjust sub-millimeter break */ 00915 brakmm = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00916 00917 /* default is 10 microns */ 00918 if( lgEOL ) 00919 { 00920 lgLogSet = true; 00921 brakmm = 9.115e-3; 00922 } 00923 00924 else if( brakmm == 0. ) 00925 { 00926 /* if second number on line is zero then set lower limit to 00927 * low-energy limit of the code. Also set linear mode, 00928 * so that last number will also be linear. */ 00929 lgLogSet = false; 00930 brakmm = rfield.tNuRyd[rfield.nspec][0]*1.001; 00931 } 00932 00933 else if( brakmm < 0. ) 00934 { 00935 /* if number is negative then this and next are logs */ 00936 lgLogSet = true; 00937 brakmm = pow(10.,brakmm); 00938 } 00939 00940 /* optional microns keyword - convert to Rydbergs */ 00941 if( nMatch("MICR",chCard) ) 00942 brakmm = RYDLAM / (1e4*brakmm); 00943 00944 rfield.tNuRyd[rfield.nspec][1] = (realnum)brakmm; 00945 00946 /* check whether this is a reasonable mm break */ 00947 if( brakmm > 1. ) 00948 fprintf(ioQQQ, 00949 " Check the order of parameters on this table power law command - the low-energy break of %f Ryd seems high to me.\n", 00950 brakmm ); 00951 00952 /* this is spectral slope, in F_nu units, between the low energy limit 00953 * and the break that may have been adjusted above 00954 * this is the slope appropriate for self-absorbed synchrotron, see eq 6.54, p.190 00955 *>>refer continuum synchrotron Rybicki, G. B., & Lightman, A.P. 1979, 00956 *>>refercon Radiative Processes in Astrophysics (New York: Wiley)*/ 00957 slopir = 5./2.; 00958 00959 /* now extrapolate a flux at this energy using the flux entered for 00960 * the first point, and this slope */ 00961 rfield.tslop[rfield.nspec][1] = 00962 (realnum)(rfield.tslop[rfield.nspec][0] + 00963 slopir*log10(rfield.tNuRyd[rfield.nspec][1]/rfield.tNuRyd[rfield.nspec][0])); 00964 00965 /* this is energy of the high-energy limit to code */ 00966 rfield.tNuRyd[rfield.nspec][3] = (realnum)rfield.egamry; 00967 00968 /* option to adjust hard X-ray break */ 00969 brakxr = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00970 00971 /* default is 50 keV */ 00972 if( lgEOL ) 00973 { 00974 brakxr = 3676.; 00975 } 00976 00977 else if( brakxr == 0. ) 00978 { 00979 brakxr = rfield.tNuRyd[rfield.nspec][3]/1.001; 00980 } 00981 00982 else if( lgLogSet ) 00983 { 00984 /* first number was negative this is a logs */ 00985 brakxr = pow(10.,brakxr); 00986 } 00987 00988 /* note that this second cutoff does not have the micron keyword */ 00989 rfield.tNuRyd[rfield.nspec][2] = (realnum)brakxr; 00990 00991 /* >>chng 03 jul 19, check that upper energy is greater than lower energy, 00992 * quit if this is not the case */ 00993 if( brakmm >= brakxr ) 00994 { 00995 fprintf( ioQQQ, " HELP!! The lower energy for the power law, %f, is greater than the upper energy, %f. This is not possible.\n Sorry.\n", 00996 brakmm , brakxr ); 00997 cdEXIT(EXIT_FAILURE); 00998 } 00999 01000 /* alpha was first option on line, is slope of mid-range */ 01001 rfield.tslop[rfield.nspec][2] = 01002 (realnum)(rfield.tslop[rfield.nspec][1] + 01003 alpha*log10(rfield.tNuRyd[rfield.nspec][2]/rfield.tNuRyd[rfield.nspec][1])); 01004 01005 /* high energy range is nu^-2 */ 01006 slopxr = -2.; 01007 01008 rfield.tslop[rfield.nspec][3] = 01009 (realnum)(rfield.tslop[rfield.nspec][2] + 01010 slopxr*log10(rfield.tNuRyd[rfield.nspec][3]/rfield.tNuRyd[rfield.nspec][2])); 01011 01012 /* following is number of portions of continuum */ 01013 ncont = 3; 01014 } 01015 01016 else if( nMatch("READ",chCard) ) 01017 { 01018 /* make sure not more than one table read command */ 01019 if( rfield.lgTableRead ) 01020 { 01021 fprintf(ioQQQ," Oops, there are more than one TABLE READ command in this input stream. I can only deal with a single TABLE READ.\n Sorry.\n"); 01022 cdEXIT(EXIT_FAILURE); 01023 } 01024 rfield.lgTableRead = true; 01025 01026 /* set up eventual read of table of points previously punched by code 01027 * get file name within double quotes, return as null terminated string 01028 * also blank out original, chCard version of name so do not trigger */ 01029 if( !lgQuoteFound ) 01030 { 01031 fprintf( ioQQQ, " Name of file must appear on TABLE READ.\n"); 01032 cdEXIT(EXIT_FAILURE); 01033 } 01034 01035 /* store file name for later reading */ 01036 rfield.ioTableRead[rfield.nspec] = string( chFile ); 01037 01038 /* set flag saying really just read in continuum exactly as punched */ 01039 strcpy( rfield.chSpType[rfield.nspec], "READ " ); 01040 01041 /* this will be flag saying continuum not read in here */ 01042 rfield.tslop[rfield.nspec][0] = 0.; 01043 rfield.tslop[rfield.nspec][1] = 0.; 01044 /* nothing left to do here, so return 01045 * the file will actually read in by routine ReadTable (in ffun.c) 01046 * called by ffun when continuum is being defined */ 01047 01048 /* number of spectra shapes that have been specified */ 01049 ++rfield.nspec; 01050 return; 01051 } 01052 01053 else if( nMatch("TLUSTY",chCard) && !nMatch("STAR",chCard) ) 01054 { 01055 /* >>chng 04 nov 30, retired TABLE TLUSTY command, PvH */ 01056 fprintf( ioQQQ, " The TABLE TLUSTY command is no longer supported.\n" ); 01057 fprintf( ioQQQ, " Please use TABLE STAR TLUSTY instead. See Hazy for details.\n" ); 01058 cdEXIT(EXIT_FAILURE); 01059 } 01060 01061 else if( nMatch("RUBI",chCard) ) 01062 { 01063 /* do Rubin modified theta 1 Ori c */ 01064 for( i=0; i < NRUBIN; i++ ) 01065 { 01066 rfield.tNuRyd[rfield.nspec][i] = (realnum)tnurbn[i]; 01067 rfield.tslop[rfield.nspec][i] = (realnum)log10(fnurbn[i] /tnurbn[i] ); 01068 } 01069 ncont = NRUBIN - 1; 01070 } 01071 01072 /* >>chng 06 jul 10, retired TABLE STARBURST command, PvH */ 01073 01074 else if( nMatch("STAR",chCard) ) 01075 { 01076 char *ptr, chMetalicity[5] = "", chODFNew[8], chVaryFlag[6] = ""; 01077 bool lgPG1159 = false, lgHCa = false, lgHNi = false, lgSolar, lgHalo, lgList; 01078 long nval, ndim=0; 01079 double Tlow = -1., Thigh = -1.; 01080 double val[MDIM]; 01081 01082 /* >>chng 06 jun 22, add support for 3 and 4-dimensional grids, PvH */ 01083 if( (ptr = strstr(chCard,"1-DI")) != NULL ) 01084 ndim = 1; 01085 else if( (ptr = strstr(chCard,"2-DI")) != NULL ) 01086 ndim = 2; 01087 else if( (ptr = strstr(chCard,"3-DI")) != NULL ) 01088 ndim = 3; 01089 else if( (ptr = strstr(chCard,"4-DI")) != NULL ) 01090 ndim = 4; 01091 01092 if( ptr != NULL ) 01093 { 01094 /* remember keyword for possible use in optimization command */ 01095 strncpy(chVaryFlag,ptr,5); 01096 chVaryFlag[5] = '\0'; 01097 /* erase this keyword, it upsets FFmtRead */ 01098 strncpy(ptr," ",5); 01099 } 01100 01101 if( nMatch( "TIME" , chCard ) ) 01102 { 01103 /* time option to vary only some continua with time */ 01104 rfield.lgTimeVary[*nqh] = true; 01105 } 01106 01107 /* >>chng 05 nov 19, add support for non-solar metalicities, PvH */ 01108 if( (ptr = strstr(chCard,"Z+1.0")) != NULL ) 01109 strncpy( chMetalicity, "p10", sizeof(chMetalicity) ); 01110 else if( (ptr = strstr(chCard,"Z+0.75")) != NULL ) 01111 strncpy( chMetalicity, "p075", sizeof(chMetalicity) ); 01112 else if( (ptr = strstr(chCard,"Z+0.5")) != NULL ) 01113 strncpy( chMetalicity, "p05", sizeof(chMetalicity) ); 01114 else if( (ptr = strstr(chCard,"Z+0.4")) != NULL ) 01115 strncpy( chMetalicity, "p04", sizeof(chMetalicity) ); 01116 else if( (ptr = strstr(chCard,"Z+0.3")) != NULL ) 01117 strncpy( chMetalicity, "p03", sizeof(chMetalicity) ); 01118 else if( (ptr = strstr(chCard,"Z+0.25")) != NULL ) 01119 strncpy( chMetalicity, "p025", sizeof(chMetalicity) ); 01120 else if( (ptr = strstr(chCard,"Z+0.2")) != NULL ) 01121 strncpy( chMetalicity, "p02", sizeof(chMetalicity) ); 01122 else if( (ptr = strstr(chCard,"Z+0.1")) != NULL ) 01123 strncpy( chMetalicity, "p01", sizeof(chMetalicity) ); 01124 else if( (ptr = strstr(chCard,"Z+0.0")) != NULL ) 01125 strncpy( chMetalicity, "p00", sizeof(chMetalicity) ); 01126 else if( (ptr = strstr(chCard,"Z-0.1")) != NULL ) 01127 strncpy( chMetalicity, "m01", sizeof(chMetalicity) ); 01128 else if( (ptr = strstr(chCard,"Z-0.2")) != NULL ) 01129 strncpy( chMetalicity, "m02", sizeof(chMetalicity) ); 01130 else if( (ptr = strstr(chCard,"Z-0.25")) != NULL ) 01131 strncpy( chMetalicity, "m025", sizeof(chMetalicity) ); 01132 else if( (ptr = strstr(chCard,"Z-0.3")) != NULL ) 01133 strncpy( chMetalicity, "m03", sizeof(chMetalicity) ); 01134 else if( (ptr = strstr(chCard,"Z-0.4")) != NULL ) 01135 strncpy( chMetalicity, "m04", sizeof(chMetalicity) ); 01136 else if( (ptr = strstr(chCard,"Z-0.5")) != NULL ) 01137 strncpy( chMetalicity, "m05", sizeof(chMetalicity) ); 01138 else if( (ptr = strstr(chCard,"Z-0.7")) != NULL ) 01139 strncpy( chMetalicity, "m07", sizeof(chMetalicity) ); 01140 else if( (ptr = strstr(chCard,"Z-0.75")) != NULL ) 01141 strncpy( chMetalicity, "m075", sizeof(chMetalicity) ); 01142 else if( (ptr = strstr(chCard,"Z-1.0")) != NULL ) 01143 strncpy( chMetalicity, "m10", sizeof(chMetalicity) ); 01144 else if( (ptr = strstr(chCard,"Z-1.3")) != NULL ) 01145 strncpy( chMetalicity, "m13", sizeof(chMetalicity) ); 01146 else if( (ptr = strstr(chCard,"Z-1.5")) != NULL ) 01147 strncpy( chMetalicity, "m15", sizeof(chMetalicity) ); 01148 else if( (ptr = strstr(chCard,"Z-1.7")) != NULL ) 01149 strncpy( chMetalicity, "m17", sizeof(chMetalicity) ); 01150 else if( (ptr = strstr(chCard,"Z-2.0")) != NULL ) 01151 strncpy( chMetalicity, "m20", sizeof(chMetalicity) ); 01152 else if( (ptr = strstr(chCard,"Z-2.5")) != NULL ) 01153 strncpy( chMetalicity, "m25", sizeof(chMetalicity) ); 01154 else if( (ptr = strstr(chCard,"Z-3.0")) != NULL ) 01155 strncpy( chMetalicity, "m30", sizeof(chMetalicity) ); 01156 else if( (ptr = strstr(chCard,"Z-3.5")) != NULL ) 01157 strncpy( chMetalicity, "m35", sizeof(chMetalicity) ); 01158 else if( (ptr = strstr(chCard,"Z-4.0")) != NULL ) 01159 strncpy( chMetalicity, "m40", sizeof(chMetalicity) ); 01160 else if( (ptr = strstr(chCard,"Z-4.5")) != NULL ) 01161 strncpy( chMetalicity, "m45", sizeof(chMetalicity) ); 01162 else if( (ptr = strstr(chCard,"Z-5.0")) != NULL ) 01163 strncpy( chMetalicity, "m50", sizeof(chMetalicity) ); 01164 else if( (ptr = strstr(chCard,"Z-INF")) != NULL ) 01165 strncpy( chMetalicity, "m99", sizeof(chMetalicity) ); 01166 else 01167 strncpy( chMetalicity, "p00", sizeof(chMetalicity) ); 01168 01169 if( ptr != NULL ) 01170 { 01171 /* remember keyword for possible use in optimization command */ 01172 strncpy(chVaryFlag,ptr,5); 01173 chVaryFlag[5] = '\0'; 01174 /* erase this keyword, it upsets FFmtRead */ 01175 strncpy(ptr," ",5); 01176 } 01177 01178 /* there are two abundance sets (solar and halo) for CoStar and Rauch H-Ca/H-Ni models. 01179 * If halo keyword appears use halo, else use solar unless other metalicity was requested */ 01180 if( nMatch("HALO",chCard) ) 01181 lgHalo = true; 01182 else 01183 lgHalo = false; 01184 lgSolar = ( !lgHalo && strcmp( chMetalicity, "p00" ) == 0 ); 01185 01186 /* >>chng 05 oct 27, added support for PG1159 Rauch models, PvH */ 01187 if( (ptr = strstr(chCard,"PG1159")) != NULL ) 01188 { 01189 lgPG1159 = true; 01190 /* erase this keyword, it upsets FFmtRead */ 01191 strncpy(ptr," ",6); 01192 } 01193 01194 if( nMatch("LIST",chCard) ) 01195 lgList = true; 01196 else 01197 lgList = false; 01198 01199 if( nMatch("AVAI",chCard) ) 01200 { 01201 AtmospheresAvail(); 01202 cdEXIT(EXIT_SUCCESS); 01203 } 01204 01205 /* now that all the keywords are out of the way, scan numbers from line image */ 01206 i = 5; 01207 for( nval=0; nval < MDIM; nval++ ) 01208 { 01209 val[nval] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 01210 if( lgEOL ) 01211 break; 01212 } 01213 01214 if( nval == 0 && !lgList ) 01215 { 01216 fprintf( ioQQQ, " The stellar temperature MUST be entered.\n" ); 01217 cdEXIT(EXIT_FAILURE); 01218 } 01219 01220 /* option for log if less than or equal to 10 */ 01221 /* Caution: For CoStar models this implicitly assumes that 01222 * the minimum ZAMS mass and youngest age is greater than 10. 01223 * As of June 1999 this is the case. However, this is not 01224 * guaranteed for possible future upgrades. */ 01225 /* match on "LOG " rather than " LOG" to avoid confusion with "LOG(G)" */ 01226 if( ( val[0] <= 10. && !nMatch("LINE",chCard) ) || nMatch("LOG ",chCard) ) 01227 val[0] = pow(10.,val[0]); 01228 01229 if( lgQuoteFound ) 01230 { 01231 nstar = GridInterpolate(val,&nval,&ndim,chFile,lgList,&Tlow,&Thigh); 01232 } 01233 01234 else if( nMatch("ATLA",chCard) ) 01235 { 01236 /* this sub-branch added by Kevin Volk, to read in large 01237 * grid of Kurucz atmospheres */ 01238 /* >>chng 05 nov 19, option TRACE is no longer supported, PvH */ 01239 01240 if( nMatch("ODFN",chCard) ) 01241 strncpy( chODFNew, "_odfnew", sizeof(chODFNew) ); 01242 else 01243 strncpy( chODFNew, "", sizeof(chODFNew) ); 01244 01245 /* >>chng 05 nov 19, add support for non-solar metalicities, PvH */ 01246 nstar = AtlasInterpolate(val,&nval,&ndim,chMetalicity,chODFNew,lgList,&Tlow,&Thigh); 01247 } 01248 01249 else if( nMatch("COST",chCard) ) 01250 { 01251 /* >>chng 99 apr 30, added CoStar stellar atmospheres */ 01252 /* this subbranch reads in CoStar atmospheres, no log(g), 01253 * but second parameter is age sequence, a number between 1 and 7, 01254 * default is 1 */ 01255 if( nMatch("INDE",chCard) ) 01256 { 01257 imode = IM_COSTAR_TEFF_MODID; 01258 if( nval == 1 ) 01259 { 01260 val[1] = 1.; 01261 nval++; 01262 } 01263 01264 /* now make sure that second parameter is within allowed range - 01265 * we do not have enough information at this time to verify temperature */ 01266 if( val[1] < 1. || val[1] > 7. ) 01267 { 01268 fprintf( ioQQQ, " The second number must be the id sequence number, 1 to 7.\n" ); 01269 fprintf( ioQQQ, " reset to 1.\n" ); 01270 val[1] = 1.; 01271 } 01272 } 01273 else if( nMatch("ZAMS",chCard) ) 01274 { 01275 imode = IM_COSTAR_MZAMS_AGE; 01276 if( nval == 1 ) 01277 { 01278 fprintf( ioQQQ, " A second number, the age of the star, must be present.\n" ); 01279 cdEXIT(EXIT_FAILURE); 01280 } 01281 } 01282 else if( nMatch(" AGE",chCard) ) 01283 { 01284 imode = IM_COSTAR_AGE_MZAMS; 01285 if( nval == 1 ) 01286 { 01287 fprintf( ioQQQ, " A second number, the ZAMS mass of the star, must be present.\n" ); 01288 cdEXIT(EXIT_FAILURE); 01289 } 01290 } 01291 else 01292 { 01293 if( nval == 1 ) 01294 { 01295 imode = IM_COSTAR_TEFF_MODID; 01296 /* default is to use ZAMS models, i.e. use index 1 */ 01297 val[1] = 1.; 01298 nval++; 01299 } 01300 else 01301 { 01302 /* Caution: the code in CoStarInterpolate implicitly 01303 * assumes that the dependence of log(g) with age is 01304 * strictly monotonic. As of June 1999 this is the case. */ 01305 imode = IM_COSTAR_TEFF_LOGG; 01306 } 01307 } 01308 01309 if( ! ( lgSolar || lgHalo ) ) 01310 { 01311 fprintf( ioQQQ, " Please choose SOLAR or HALO abundances.\n" ); 01312 cdEXIT(EXIT_FAILURE); 01313 } 01314 01315 nstar = CoStarInterpolate(val,&nval,&ndim,imode,lgHalo,lgList,&Tlow,&Thigh); 01316 } 01317 01318 else if( nMatch("KURU",chCard) ) 01319 { 01320 nstar = Kurucz79Interpolate(val,&nval,&ndim,lgList,&Tlow,&Thigh); 01321 } 01322 01323 else if( nMatch("MIHA",chCard) ) 01324 { 01325 nstar = MihalasInterpolate(val,&nval,&ndim,lgList,&Tlow,&Thigh); 01326 } 01327 01328 else if( nMatch("RAUC",chCard) ) 01329 { 01330 if( ! ( lgSolar || lgHalo ) ) 01331 { 01332 fprintf( ioQQQ, " Please choose SOLAR or HALO abundances.\n" ); 01333 cdEXIT(EXIT_FAILURE); 01334 } 01335 01336 /* >>chng 97 aug 23, K Volk added Rauch stellar atmospheres */ 01337 /* >>chng 05 oct 27, added support for PG1159 Rauch models, PvH */ 01338 /* >>chng 06 jun 26, added support for H, He, H+He Rauch models, PvH */ 01339 if( nMatch("H-CA",chCard) || nMatch(" OLD",chCard) ) 01340 { 01341 lgHCa = true; 01342 nstar = RauchInterpolateHCa(val,&nval,&ndim,lgHalo,lgList,&Tlow,&Thigh); 01343 } 01344 else if( nMatch("HYDR",chCard) ) 01345 { 01346 nstar = RauchInterpolateHydr(val,&nval,&ndim,lgList,&Tlow,&Thigh); 01347 } 01348 else if( nMatch("HELI",chCard) ) 01349 { 01350 nstar = RauchInterpolateHelium(val,&nval,&ndim,lgList,&Tlow,&Thigh); 01351 } 01352 else if( nMatch("H+HE",chCard) ) 01353 { 01354 nstar = RauchInterpolateHpHe(val,&nval,&ndim,lgList,&Tlow,&Thigh); 01355 } 01356 else if( lgPG1159 ) /* the keyword PG1159 was already matched and erased above */ 01357 { 01358 nstar = RauchInterpolatePG1159(val,&nval,&ndim,lgList,&Tlow,&Thigh); 01359 } 01360 else 01361 { 01362 lgHNi = true; 01363 nstar = RauchInterpolateHNi(val,&nval,&ndim,lgHalo,lgList,&Tlow,&Thigh); 01364 } 01365 } 01366 01367 else if( nMatch("TLUS",chCard) ) 01368 { 01369 if( nMatch("BSTA",chCard) ) 01370 { 01371 /* >>chng 06 oct 19, this sub-branch added to read in 01372 * large 2006 grid of Tlusty B-star atmospheres, PvH */ 01373 nstar = TlustyInterpolate(val,&nval,&ndim,TL_BSTAR,chMetalicity,lgList,&Tlow,&Thigh); 01374 } 01375 else if( nMatch("OSTA",chCard) ) 01376 { 01377 /* >>chng 05 nov 21, this sub-branch added to read in 01378 * large 2002 grid of Tlusty O-star atmospheres, PvH */ 01379 nstar = TlustyInterpolate(val,&nval,&ndim,TL_OSTAR,chMetalicity,lgList,&Tlow,&Thigh); 01380 } 01381 else 01382 { 01383 fprintf( ioQQQ, " There must be a third key on TABLE STAR TLUSTY;" ); 01384 fprintf( ioQQQ, " the options are BSTAR, OSTAR.\n" ); 01385 cdEXIT(EXIT_FAILURE); 01386 } 01387 } 01388 01389 else if( nMatch("WERN",chCard) ) 01390 { 01391 /* this subbranch added by Kevin Volk, to read in large 01392 * grid of hot PN atmospheres computer by Klaus Werner */ 01393 nstar = WernerInterpolate(val,&nval,&ndim,lgList,&Tlow,&Thigh); 01394 } 01395 01396 else if( nMatch("WMBA",chCard) ) 01397 { 01398 /* >>chng 06 jun 27, this subbranch added to read in 01399 * grid of hot atmospheres computer by Pauldrach */ 01400 nstar = WMBASICInterpolate(val,&nval,&ndim,lgList,&Tlow,&Thigh); 01401 } 01402 01403 else 01404 { 01405 fprintf( ioQQQ, " There must be a second key on TABLE STAR;" ); 01406 fprintf( ioQQQ, " the options are ATLAS, KURUCZ, MIHALAS, RAUCH, WERNER, and WMBASIC.\n" ); 01407 cdEXIT(EXIT_FAILURE); 01408 } 01409 01410 /* set flag saying really just read in continuum exactly as punched */ 01411 strcpy( rfield.chSpType[rfield.nspec], "VOLK " ); 01412 01413 ncont = nstar - 1; 01414 01415 /* vary option */ 01416 if( optimize.lgVarOn ) 01417 { 01418 optimize.vincr[optimize.nparm] = (realnum)0.1; 01419 01420 if( lgQuoteFound ) 01421 { 01422 /* this is number of parameters to feed onto the input line */ 01423 optimize.nvarxt[optimize.nparm] = nval; 01424 sprintf( optimize.chVarFmt[optimize.nparm], "TABLE STAR \"%s\"", chFile ); 01425 strcat( optimize.chVarFmt[optimize.nparm], " %f LOG" ); 01426 for( i=1; i < nval; i++ ) 01427 strcat( optimize.chVarFmt[optimize.nparm], " , %f" ); 01428 } 01429 01430 if( nMatch("ATLA",chCard) ) 01431 { 01432 /* this is number of parameters to feed onto the input line */ 01433 optimize.nvarxt[optimize.nparm] = ndim; 01434 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE STAR ATLAS " ); 01435 strcat( optimize.chVarFmt[optimize.nparm], chVaryFlag ); 01436 if( nMatch("ODFN",chCard) ) 01437 strcat( optimize.chVarFmt[optimize.nparm], " ODFNEW" ); 01438 01439 strcat( optimize.chVarFmt[optimize.nparm], " %f LOG , LOG(G)= %f" ); 01440 01441 if( ndim == 3 ) 01442 strcat( optimize.chVarFmt[optimize.nparm], " , LOG(Z)= %f" ); 01443 01444 } 01445 01446 else if( nMatch("COST",chCard) ) 01447 { 01448 /* this is number of parameters to feed onto the input line */ 01449 optimize.nvarxt[optimize.nparm] = 2; 01450 01451 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE STAR COSTAR " ); 01452 if( lgHalo ) 01453 strcat( optimize.chVarFmt[optimize.nparm], "HALO " ); 01454 01455 if( imode == IM_COSTAR_TEFF_MODID ) 01456 { 01457 strcat( optimize.chVarFmt[optimize.nparm], "TEFF= %f LOG , INDEX= %f" ); 01458 } 01459 else if( imode == IM_COSTAR_TEFF_LOGG ) 01460 { 01461 strcat( optimize.chVarFmt[optimize.nparm], "TEFF= %f LOG , LOG(G)= %f" ); 01462 } 01463 else if( imode == IM_COSTAR_MZAMS_AGE ) 01464 { 01465 /* don't use keyword _AGE here! */ 01466 strcat( optimize.chVarFmt[optimize.nparm], "ZAMS MASS= %f LOG , TIME= %f" ); 01467 } 01468 else if( imode == IM_COSTAR_AGE_MZAMS ) 01469 { 01470 /* don't use keyword ZAMS here! */ 01471 strcat( optimize.chVarFmt[optimize.nparm], "AGE= %f LOG , MASS= %f" ); 01472 optimize.vincr[optimize.nparm] = (realnum)0.5; 01473 } 01474 } 01475 01476 else if( nMatch("KURU",chCard) ) 01477 { 01478 /* this is number of parameters to feed onto the input line */ 01479 optimize.nvarxt[optimize.nparm] = 1; 01480 strcpy( optimize.chVarFmt[optimize.nparm], 01481 "TABLE STAR KURUCZ %f LOG" ); 01482 } 01483 01484 else if( nMatch("MIHA",chCard) ) 01485 { 01486 /* this is number of parameters to feed onto the input line */ 01487 optimize.nvarxt[optimize.nparm] = 1; 01488 strcpy( optimize.chVarFmt[optimize.nparm], 01489 "TABLE STAR MIHALAS %f LOG" ); 01490 } 01491 01492 else if( nMatch("RAUC",chCard) ) 01493 { 01494 /* this is number of parameters to feed onto the input line */ 01495 optimize.nvarxt[optimize.nparm] = ndim; 01496 01497 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE STAR RAUCH " ); 01498 01499 if( nMatch("HYDR",chCard) ) 01500 strcat( optimize.chVarFmt[optimize.nparm], "HYDR " ); 01501 else if( nMatch("HELI",chCard) ) 01502 strcat( optimize.chVarFmt[optimize.nparm], "HELIUM " ); 01503 else if( nMatch("H+HE",chCard) ) 01504 strcat( optimize.chVarFmt[optimize.nparm], "H+HE " ); 01505 else if( lgPG1159 ) 01506 strcat( optimize.chVarFmt[optimize.nparm], "PG1159 " ); 01507 else if( lgHCa ) 01508 strcat( optimize.chVarFmt[optimize.nparm], "H-CA " ); 01509 01510 if( ( lgHCa || lgHNi ) && ndim == 2 ) 01511 { 01512 if( lgHalo ) 01513 strcat( optimize.chVarFmt[optimize.nparm], "HALO " ); 01514 else 01515 strcat( optimize.chVarFmt[optimize.nparm], "SOLAR " ); 01516 } 01517 01518 strcat( optimize.chVarFmt[optimize.nparm], "%f LOG , LOG(G)= %f" ); 01519 01520 if( ndim == 3 ) 01521 { 01522 if( nMatch("H+HE",chCard) ) 01523 strcat( optimize.chVarFmt[optimize.nparm], " , F(HE)= %f" ); 01524 else 01525 strcat( optimize.chVarFmt[optimize.nparm], " , LOG(Z)= %f 3-DIM" ); 01526 } 01527 } 01528 01529 if( nMatch("TLUS",chCard) ) 01530 { 01531 /* this is number of parameters to feed onto the input line */ 01532 optimize.nvarxt[optimize.nparm] = ndim; 01533 strcpy( optimize.chVarFmt[optimize.nparm], "TABLE STAR TLUSTY " ); 01534 if( nMatch("BSTA",chCard) ) 01535 strcat( optimize.chVarFmt[optimize.nparm], "BSTAR " ); 01536 else if( nMatch("OSTA",chCard) ) 01537 strcat( optimize.chVarFmt[optimize.nparm], "OSTAR " ); 01538 else 01539 TotalInsanity(); 01540 strcat( optimize.chVarFmt[optimize.nparm], chVaryFlag ); 01541 strcat( optimize.chVarFmt[optimize.nparm], " %f LOG , LOG(G)= %f" ); 01542 01543 if( ndim == 3 ) 01544 strcat( optimize.chVarFmt[optimize.nparm], " , LOG(Z)= %f" ); 01545 } 01546 01547 else if( nMatch("WERN",chCard) ) 01548 { 01549 /* this is number of parameters to feed onto the input line */ 01550 optimize.nvarxt[optimize.nparm] = 2; 01551 strcpy( optimize.chVarFmt[optimize.nparm], 01552 "TABLE STAR WERNER %f LOG , LOG(G)= %f" ); 01553 } 01554 01555 else if( nMatch("WMBA",chCard) ) 01556 { 01557 /* this is number of parameters to feed onto the input line */ 01558 optimize.nvarxt[optimize.nparm] = 3; 01559 strcpy( optimize.chVarFmt[optimize.nparm], 01560 "TABLE STAR WMBASIC %f LOG , LOG(G)= %f , LOG(Z)= %f" ); 01561 } 01562 01563 /* pointer to where to write */ 01564 optimize.nvfpnt[optimize.nparm] = input.nRead; 01565 01566 ASSERT( nval <= LIMEXT ); 01567 01568 /* log of temp will be pointer */ 01569 optimize.vparm[0][optimize.nparm] = (realnum)log10(val[0]); 01570 for( i=1; i < nval; i++ ) 01571 optimize.vparm[i][optimize.nparm] = (realnum)val[i]; 01572 01573 /* following are upper and lower limits to temperature range */ 01574 optimize.varang[optimize.nparm][0] = (realnum)log10(Tlow); 01575 optimize.varang[optimize.nparm][1] = (realnum)log10(Thigh); 01576 01577 /* finally increment this counter */ 01578 ++optimize.nparm; 01579 } 01580 } 01581 01582 else 01583 { 01584 fprintf( ioQQQ, " There MUST be a keyword on this line. The keys are:" 01585 "AGN, AKN120, CRAB, COOL, DRAINE, HM96, HM05, _ISM, LINE, POWERlaw, " 01586 "READ, RUBIN, and STAR.\n Sorry.\n" ); 01587 cdEXIT(EXIT_FAILURE); 01588 } 01589 01590 /* table star Werner and table star atlas are special 01591 * cases put in by Kevin Volk - they are not really tables 01592 * at all, so return if chSpType is Volk */ 01593 if( strcmp(rfield.chSpType[rfield.nspec],"VOLK ") == 0 ) 01594 { 01595 ++rfield.nspec; 01596 return; 01597 } 01598 01599 /* this flag set true if we did not parse a continuum source, thus creating 01600 * the arrays that are needed - return at this point */ 01601 if( lgNoContinuum ) 01602 { 01603 return; 01604 } 01605 01606 /* ncont only positive when real continuum entered 01607 * convert from log(Hz) to Ryd if first nu>5 */ 01608 if( rfield.tNuRyd[rfield.nspec][0] >= 5. ) 01609 { 01610 for( i=0; i < (ncont + 1); i++ ) 01611 { 01612 rfield.tNuRyd[rfield.nspec][i] = 01613 (realnum)pow(10.,rfield.tNuRyd[rfield.nspec][i] - 15.51718); 01614 } 01615 } 01616 01617 else if( rfield.tNuRyd[rfield.nspec][0] < 0. ) 01618 { 01619 /* energies entered as logs */ 01620 for( i=0; i < (ncont + 1); i++ ) 01621 { 01622 rfield.tNuRyd[rfield.nspec][i] = 01623 (realnum)pow(10.,(double)rfield.tNuRyd[rfield.nspec][i]); 01624 } 01625 } 01626 01627 /* tFluxLog will be log(fnu) at that spot, read into tslop */ 01628 for( i=0; i < (ncont + 1); i++ ) 01629 { 01630 rfield.tFluxLog[rfield.nspec][i] = rfield.tslop[rfield.nspec][i]; 01631 } 01632 01633 for( i=0; i < ncont; i++ ) 01634 { 01635 rfield.tslop[rfield.nspec][i] = 01636 (realnum)((rfield.tslop[rfield.nspec][i+1] - 01637 rfield.tslop[rfield.nspec][i])/ 01638 log10(rfield.tNuRyd[rfield.nspec][i+1]/ 01639 rfield.tNuRyd[rfield.nspec][i])); 01640 } 01641 01642 if( ncont > 0 && (ncont + 2 < rfield.nupper) ) 01643 { 01644 /* zero out remainder of array */ 01645 for( i=ncont + 1; i < rfield.nupper; i++ ) 01646 { 01647 rfield.tNuRyd[rfield.nspec][i] = 0.; 01648 } 01649 } 01650 01651 if( trace.lgConBug && trace.lgTrace ) 01652 { 01653 fprintf( ioQQQ, " Table for this continuum; TNU, TFAC, TSLOP\n" ); 01654 for( i=0; i < (ncont + 1); i++ ) 01655 { 01656 fprintf( ioQQQ, "%12.4e %12.4e %12.4e\n", rfield.tNuRyd[rfield.nspec][i], 01657 rfield.tFluxLog[rfield.nspec][i], rfield.tslop[rfield.nspec][i] ); 01658 } 01659 } 01660 01661 /* renormalize continuum so that quantity passes through unity at 1 Ryd 01662 * lgHit will be set false when we get a hit in following loop, 01663 * then test made to make sure it happened*/ 01664 lgHit = false; 01665 /*following will be reset when hit occurs*/ 01666 fac = -DBL_MAX; 01667 /* >>chng 04 mar 16, chng loop so breaks when hit, previously had cycled 01668 * until last point reached, so last point always used */ 01669 for( i=0; i < (ncont + 1) && !lgHit; i++ ) 01670 { 01671 if( rfield.tNuRyd[rfield.nspec][i] > 1. ) 01672 { 01673 fac = rfield.tFluxLog[rfield.nspec][i]; 01674 lgHit = true; 01675 } 01676 } 01677 01678 if( ncont > 0 && lgHit ) 01679 { 01680 /* do the renormalization */ 01681 for( i=0; i < (ncont + 1); i++ ) 01682 { 01683 rfield.tFluxLog[rfield.nspec][i] -= (realnum)fac; 01684 } 01685 } 01686 01687 if( ncont > 0 ) 01688 ++rfield.nspec; 01689 return; 01690 } 01691 01692 01693 /* this allows the low energy point of any built in array to be reset to the 01694 * current low energy point in the code - nothing need be done if this is reset 01695 * tnu is array of energies, [0] is first, and we want it to be lower 01696 * fluxlog is flux at tnu, and may or may not be log 01697 * lgLog says whether it is */ 01698 STATIC void resetBltin( double *tnu , double *fluxlog , bool lgLog ) 01699 { 01700 /* this will multiply low-energy bounds of code and go into element[0] 01701 * ensures that energy range is fully covered */ 01702 const double RESETFACTOR = 0.98; 01703 double power; 01704 /* this makes sure we are called after emm is defined */ 01705 ASSERT( rfield.emm > 0. ); 01706 01707 if( lgLog ) 01708 { 01709 /* continuum comes in as log of flux */ 01710 /* this is current power-law slope of low-energy continuum */ 01711 power = (fluxlog[1] - fluxlog[0] ) / log10( tnu[1]/tnu[0] ); 01712 /* this will be new low energy bounds to this continuum */ 01713 tnu[0] = rfield.emm*RESETFACTOR; 01714 fluxlog[0] = fluxlog[1] + power * log10( tnu[0] /tnu[1] ); 01715 } 01716 else 01717 { 01718 /* continuum comes in as linear flux */ 01719 /* this is current power-law slope of low-energy continuum */ 01720 power = log10( fluxlog[1]/fluxlog[0]) / log10( tnu[1]/tnu[0] ); 01721 /* this will be new low energy bounds to this continuum */ 01722 tnu[0] = rfield.emm*RESETFACTOR; 01723 fluxlog[0] = log10(fluxlog[1]) + power * log10( tnu[0] /tnu[1] ); 01724 /* flux is not really log, we want linear */ 01725 fluxlog[0] = pow(10. , fluxlog[0]); 01726 } 01727 /*fprintf(ioQQQ," power is %f lgLog is %i\n", power, lgLog );*/ 01728 return; 01729 } 01730 01731 /*ZeroContin store sets of built in continua */ 01732 STATIC void ZeroContin(void) 01733 { 01734 01735 DEBUG_ENTRY( "ZeroContin()" ); 01736 01737 /* Draine 1978 ISM continuum */ 01738 /* freq in ryd */ 01739 tnudrn[0] = 0.3676; 01740 tnudrn[1] = 0.4144; 01741 tnudrn[2] = 0.4558; 01742 tnudrn[3] = 0.5064; 01743 tnudrn[4] = 0.5698; 01744 tnudrn[5] = 0.6511; 01745 tnudrn[6] = 0.7012; 01746 tnudrn[7] = 0.7597; 01747 tnudrn[8] = 0.8220; 01748 tnudrn[9] = 0.9116; 01749 tnudrn[10] = 0.9120; 01750 tnudrn[11] = 0.9306; 01751 tnudrn[12] = 0.9600; 01752 tnudrn[13] = 0.9806; 01753 /* >>chng 05 aug 03, this energy is so close to 1 ryd that it spills over 01754 * into the H-ionizing continuum - move down to 0.99 */ 01755 /* >>chng 05 aug 08, this destabilized pdr_leiden_hack_v4 (!) so put back to 01756 * original value and include extinguish command */ 01757 tnudrn[14] = 0.9999; 01758 /*tnudrn[14] = 0.99;*/ 01759 01760 /* f_nu from equation 23 of 01761 * >>refer ism field Draine, B.T. & Bertoldi, F. 1996, ApJ, 468, 269 */ 01762 int i; 01763 i= 0; 01764 tsldrn[i] = -17.8063; 01765 ++i;tsldrn[i] = -17.7575; 01766 ++i;tsldrn[i] = -17.7268; 01767 ++i;tsldrn[i] = -17.7036; 01768 ++i;tsldrn[i] = -17.6953; 01769 ++i;tsldrn[i] = -17.7182; 01770 ++i;tsldrn[i] = -17.7524; 01771 ++i;tsldrn[i] = -17.8154; 01772 ++i;tsldrn[i] = -17.9176; 01773 ++i;tsldrn[i] = -18.1675; 01774 ++i;tsldrn[i] = -18.1690; 01775 ++i;tsldrn[i] = -18.2477; 01776 ++i;tsldrn[i] = -18.4075; 01777 ++i;tsldrn[i] = -18.5624; 01778 ++i;tsldrn[i] = -18.7722; 01779 01780 /* generic AGN continuum taken from 01781 * >>refer AGN cont Mathews and Ferland ApJ Dec15 '87 01782 * except self absorption at 10 microns, nu**-2.5 below that */ 01783 tnuagn[0] = 1e-5; 01784 tnuagn[1] = 9.12e-3; 01785 tnuagn[2] = .206; 01786 tnuagn[3] = 1.743; 01787 tnuagn[4] = 4.13; 01788 tnuagn[5] = 26.84; 01789 tnuagn[6] = 7353.; 01790 tnuagn[7] = 7.4e6; 01791 01792 tslagn[0] = -3.388; 01793 tslagn[1] = 4.0115; 01794 tslagn[2] = 2.6576; 01795 tslagn[3] = 2.194; 01796 tslagn[4] = 1.819; 01797 tslagn[5] = -.6192; 01798 tslagn[6] = -2.326; 01799 tslagn[7] = -7.34; 01800 resetBltin( tnuagn , tslagn , true ); 01801 01802 /* Crab Nebula continuum from 01803 *>>refer continuum CrabNebula Davidson, K., & Fesen, 1985, ARAA, 01804 * second vector is F_nu as seen from Earth */ 01805 tnucrb[0] = 1.0e-5; 01806 tnucrb[1] = 5.2e-4; 01807 tnucrb[2] = 1.5e-3; 01808 tnucrb[3] = 0.11; 01809 tnucrb[4] = 0.73; 01810 tnucrb[5] = 7.3; 01811 tnucrb[6] = 73.; 01812 tnucrb[7] = 7300.; 01813 tnucrb[8] = 1.5e6; 01814 tnucrb[9] = 7.4e6; 01815 01816 fnucrb[0] = 3.77e-21; 01817 fnucrb[1] = 1.38e-21; 01818 fnucrb[2] = 2.10e-21; 01819 fnucrb[3] = 4.92e-23; 01820 fnucrb[4] = 1.90e-23; 01821 fnucrb[5] = 2.24e-24; 01822 fnucrb[6] = 6.42e-26; 01823 fnucrb[7] = 4.02e-28; 01824 fnucrb[8] = 2.08e-31; 01825 fnucrb[9] = 1.66e-32; 01826 resetBltin( tnucrb , fnucrb , false ); 01827 01828 /* Bob Rubin's theta 1 Ori C continuum, modified from Kurucz to 01829 * get NeIII right */ 01830 /* >>chng 02 may 02, revise continuum while working with Bob Rubin on NII revisit */ 01831 resetBltin( tnurbn , fnurbn , false ); 01832 01833 /* cooling flow continuum from Johnstone et al. MNRAS 255, 431. */ 01834 cfle[0] = 0.0000100; 01835 cflf[0] = -0.8046910; 01836 cfle[1] = 0.7354023; 01837 cflf[1] = -0.8046910; 01838 cfle[2] = 1.4708046; 01839 cflf[2] = -0.7436830; 01840 cfle[3] = 2.2062068; 01841 cflf[3] = -0.6818757; 01842 cfle[4] = 2.9416091; 01843 cflf[4] = -0.7168990; 01844 cfle[5] = 3.6770115; 01845 cflf[5] = -0.8068384; 01846 cfle[6] = 4.4124136; 01847 cflf[6] = -0.6722584; 01848 cfle[7] = 5.1478162; 01849 cflf[7] = -0.7626385; 01850 cfle[8] = 5.8832183; 01851 cflf[8] = -1.0396487; 01852 cfle[9] = 6.6186204; 01853 cflf[9] = -0.7972314; 01854 cfle[10] = 7.3540230; 01855 cflf[10] = -0.9883416; 01856 cfle[11] = 14.7080460; 01857 cflf[11] = -1.1675659; 01858 cfle[12] = 22.0620689; 01859 cflf[12] = -1.1985949; 01860 cfle[13] = 29.4160919; 01861 cflf[13] = -1.2263466; 01862 cfle[14] = 36.7701149; 01863 cflf[14] = -1.2918345; 01864 cfle[15] = 44.1241379; 01865 cflf[15] = -1.3510833; 01866 cfle[16] = 51.4781609; 01867 cflf[16] = -1.2715496; 01868 cfle[17] = 58.8321838; 01869 cflf[17] = -1.1098027; 01870 cfle[18] = 66.1862030; 01871 cflf[18] = -1.4315782; 01872 cfle[19] = 73.5402298; 01873 cflf[19] = -1.1327956; 01874 cfle[20] = 147.080459; 01875 cflf[20] = -1.6869649; 01876 cfle[21] = 220.620681; 01877 cflf[21] = -2.0239367; 01878 cfle[22] = 294.160919; 01879 cflf[22] = -2.2130392; 01880 cfle[23] = 367.701141; 01881 cflf[23] = -2.3773901; 01882 cfle[24] = 441.241363; 01883 cflf[24] = -2.5326197; 01884 cfle[25] = 514.7816162; 01885 cflf[25] = -2.5292389; 01886 cfle[26] = 588.3218384; 01887 cflf[26] = -2.8230250; 01888 cfle[27] = 661.8620605; 01889 cflf[27] = -2.9502323; 01890 cfle[28] = 735.4022827; 01891 cflf[28] = -3.0774822; 01892 cfle[29] = 1470.8045654; 01893 cflf[29] = -4.2239799; 01894 cfle[30] = 2206.2067871; 01895 cflf[30] = -5.2547927; 01896 cfle[31] = 2941.6091309; 01897 cflf[31] = -6.2353640; 01898 cfle[32] = 3677.0114746; 01899 cflf[32] = -7.1898708; 01900 cfle[33] = 4412.4135742; 01901 cflf[33] = -8.1292381; 01902 cfle[34] = 5147.8159180; 01903 cflf[34] = -9.0594845; 01904 cfle[35] = 5883.2182617; 01905 cflf[35] = -9.9830370; 01906 cfle[36] = 6618.6206055; 01907 cflf[36] = -10.9028034; 01908 cfle[37] = 7354.0229492; 01909 cflf[37] = -11.8188877; 01910 cfle[38] = 7400.0000000; 01911 cflf[38] = -30.0000000; 01912 cfle[39] = 10000000.0000000; 01913 cflf[39] = -30.0000000; 01914 resetBltin( cfle , cflf , true ); 01915 01916 /* AKN120 continuum from Brad Peterson's plot 01917 * second vector is F_nu*1E10 as seen from Earth */ 01918 tnuakn[0] = 1e-5; 01919 tnuakn[1] = 1.9e-5; 01920 tnuakn[2] = 3.0e-4; 01921 tnuakn[3] = 2.4e-2; 01922 tnuakn[4] = 0.15; 01923 tnuakn[5] = 0.30; 01924 tnuakn[6] = 0.76; 01925 tnuakn[7] = 2.0; 01926 tnuakn[8] = 76.0; 01927 tnuakn[9] = 760.; 01928 tnuakn[10] = 7.4e6; 01929 fnuakn[0] = 1.5e-16; 01930 fnuakn[1] = 1.6e-16; 01931 fnuakn[2] = 1.4e-13; 01932 fnuakn[3] = 8.0e-15; 01933 fnuakn[4] = 1.6e-15; 01934 fnuakn[5] = 1.8e-15; 01935 fnuakn[6] = 7.1e-16; 01936 fnuakn[7] = 7.9e-17; 01937 fnuakn[8] = 1.1e-18; 01938 fnuakn[9] = 7.1e-20; 01939 fnuakn[10] = 1.3e-24; 01940 resetBltin( fnuakn , fnuakn , false ); 01941 01942 /* interstellar radiation field from Black 1987, "Interstellar Processes" 01943 * table of log nu, log nu*fnu taken from his figure 2 */ 01944 /* >>chng 99 jun 14 energy range lowered to 1e-8 ryd */ 01945 tnuism[0] = 6.00; 01946 /*tnuism[0] = 9.00;*/ 01947 tnuism[1] = 10.72; 01948 tnuism[2] = 11.00; 01949 tnuism[3] = 11.23; 01950 tnuism[4] = 11.47; 01951 tnuism[5] = 11.55; 01952 tnuism[6] = 11.85; 01953 tnuism[7] = 12.26; 01954 tnuism[8] = 12.54; 01955 tnuism[9] = 12.71; 01956 tnuism[10] = 13.10; 01957 tnuism[11] = 13.64; 01958 tnuism[12] = 14.14; 01959 tnuism[13] = 14.38; 01960 tnuism[14] = 14.63; 01961 tnuism[15] = 14.93; 01962 tnuism[16] = 15.08; 01963 tnuism[17] = 15.36; 01964 tnuism[18] = 15.43; 01965 tnuism[19] = 16.25; 01966 tnuism[20] = 17.09; 01967 tnuism[21] = 18.00; 01968 tnuism[22] = 23.00; 01969 /* these are log nu*Fnu */ 01970 fnuism[0] = -16.708; 01971 /*fnuism[0] = -7.97;*/ 01972 fnuism[1] = -2.96; 01973 fnuism[2] = -2.47; 01974 fnuism[3] = -2.09; 01975 fnuism[4] = -2.11; 01976 fnuism[5] = -2.34; 01977 fnuism[6] = -3.66; 01978 fnuism[7] = -2.72; 01979 fnuism[8] = -2.45; 01980 fnuism[9] = -2.57; 01981 fnuism[10] = -3.85; 01982 fnuism[11] = -3.34; 01983 fnuism[12] = -2.30; 01984 fnuism[13] = -1.79; 01985 fnuism[14] = -1.79; 01986 fnuism[15] = -2.34; 01987 fnuism[16] = -2.72; 01988 fnuism[17] = -2.55; 01989 fnuism[18] = -2.62; 01990 fnuism[19] = -5.68; 01991 fnuism[20] = -6.45; 01992 fnuism[21] = -6.30; 01993 fnuism[22] = -11.3; 01994 01995 return; 01996 } 01997 01998 /*lines_table invoked by table lines command, check if we can find all lines in a given list */ 01999 int lines_table(void) 02000 { 02001 long int n, 02002 miss; 02003 02004 if( !nLINE_TABLE ) 02005 return 0; 02006 02007 DEBUG_ENTRY( "lines_table()" ); 02008 02009 fprintf( ioQQQ , "lines_table checking lines within data table %s\n", chLINE_LIST ); 02010 miss = 0; 02011 02012 /* \todo 2 DOS carriage return on linux screws this up. Can we overlook the CR? Skip in cdgetlinelist? */ 02013 for( n=0; n<nLINE_TABLE; ++n ) 02014 { 02015 double relative , absolute; 02016 if( (cdLine( chLabel[n], wl[n] , &relative , &absolute ))<=0 ) 02017 { 02018 ++miss; 02019 fprintf(ioQQQ,"lines_table in parse_table.cpp did not find line %4s ",chLabel[n]); 02020 prt_wl(ioQQQ,wl[n]); 02021 fprintf(ioQQQ,"\n"); 02022 } 02023 } 02024 if( miss ) 02025 { 02026 /* this is so that we pick up problem automatically */ 02027 fprintf( ioQQQ , " BOTCHED ASSERTS!!! Botched Asserts!!! lines_table could not find a total of %li lines\n\n", miss ); 02028 } 02029 else 02030 { 02031 fprintf( ioQQQ , "lines_table found all lines\n\n" ); 02032 } 02033 return miss; 02034 02035 } 02036 #if defined(__HP_aCC) 02037 #pragma OPTIMIZE OFF 02038 #pragma OPTIMIZE ON 02039 #endif