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 /*ParseDrive parse the drive command - drive calls to various subs */ 00004 /*DrvCaseBHS allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */ 00005 /*DrvHyas allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */ 00006 /*dgaunt drive gaunt factor routines by letting user query values */ 00007 #include "cddefines.h" 00008 #include "trace.h" 00009 #include "hydro_bauman.h" 00010 #include "hydroeinsta.h" 00011 #include "phycon.h" 00012 #include "atmdat.h" 00013 #include "taulines.h" 00014 #include "thermal.h" 00015 #include "abund.h" 00016 #include "rt.h" 00017 #include "continuum.h" 00018 #include "parse.h" 00019 #include "physconst.h" 00020 00021 /*dgaunt drive gaunt factor routines by letting user query values */ 00022 STATIC void dgaunt(void); 00023 00024 /*DrvHyas allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */ 00025 STATIC void DrvHyas(void); 00026 00027 /* drive escape probability routines */ 00028 STATIC void DrvEscP( void ); 00029 00030 /*DrvCaseBHS allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */ 00031 STATIC void DrvCaseBHS(void); 00032 00033 void ParseDrive(char *chCard ) 00034 { 00035 bool lgEOL; 00036 long int n, 00037 i; 00038 double fac, 00039 zed; 00040 00041 DEBUG_ENTRY( "ParseDrive()" ); 00042 00043 /* NB evolve all following names to style DrvSomething */ 00044 00045 /* option to drive cloudy, which one? */ 00046 if( nMatch("FFMT",chCard) ) 00047 { 00048 /* free format parser */ 00049 char chInput[INPUT_LINE_LENGTH]; 00050 fprintf( ioQQQ, " FFmtRead ParseDrive entered. Enter number.\n" ); 00051 lgEOL = false; 00052 while( !lgEOL ) 00053 { 00054 if( read_whole_line( chInput , (int)sizeof(chInput) , ioStdin ) == NULL ) 00055 { 00056 fprintf( ioQQQ, " ParseDrive.dat error getting magic number\n" ); 00057 cdEXIT(EXIT_FAILURE); 00058 } 00059 i = 1; 00060 fac = FFmtRead(chInput,&i,INPUT_LINE_LENGTH,&lgEOL); 00061 if( lgEOL ) 00062 { 00063 fprintf( ioQQQ, " FFmtRead hit the EOL with no value, return=%10.2e\n", 00064 fac ); 00065 break; 00066 } 00067 else if( fac == 0. ) 00068 { 00069 break; 00070 } 00071 else 00072 { 00073 fprintf( ioQQQ, " FFmtRead returned with value%11.4e\n", 00074 fac ); 00075 } 00076 fprintf( ioQQQ, " Enter 0 to stop, or another value.\n" ); 00077 } 00078 fprintf( ioQQQ, " FFmtRead ParseDrive exits.\n" ); 00079 } 00080 00081 else if( nMatch("CASE",chCard) ) 00082 { 00083 /* option to interpolate on Hummer and Storey case b hydrogenic emission routines */ 00084 DrvCaseBHS( ); 00085 } 00086 00087 else if( nMatch("CDLI",chCard) ) 00088 { 00089 /* drive cdLine to check that it finds all the right lines, routine is in lines.c */ 00090 trace.lgDrv_cdLine = true; 00091 } 00092 00093 else if( nMatch(" E1 ",chCard) ) 00094 { 00095 /* option to drive exponential integral routines */ 00096 for( i=0; i<30; ++i ) 00097 { 00098 double tau = pow( 10. , ((double)i/4. - 5.) ); 00099 fprintf(ioQQQ,"tau\t%.3e\t exp-tau\t%.5e\t e1 tau\t%.5e \t ee2 tau\t%.5e \n", 00100 tau , sexp(tau) , ee1(tau) , e2(tau ) ); 00101 } 00102 cdEXIT(EXIT_SUCCESS); 00103 } 00104 00105 else if( nMatch("ESCA",chCard) ) 00106 { 00107 /* option to drive escape probability routines */ 00108 DrvEscP( ); 00109 } 00110 00111 else if( nMatch("HYAS",chCard) ) 00112 { 00113 /* option to drive Jason's hydrogen transition probabilities */ 00114 DrvHyas(); 00115 } 00116 00117 else if( nMatch("GAUN",chCard) ) 00118 { 00119 /* drive gaunt factor routine */ 00120 dgaunt(); 00121 } 00122 00123 else if( nMatch("POIN",chCard) ) 00124 { 00125 /* later on, check cell pointers, centers, widths */ 00126 fprintf( ioQQQ, " Define entire model first, later I will ask for energies.\n" ); 00127 trace.lgPtrace = true; 00128 } 00129 00130 else if( nMatch("PUMP",chCard) ) 00131 { 00132 char chInput[INPUT_LINE_LENGTH]; 00133 lgEOL = false; 00134 fprintf( ioQQQ, " Continuum pump ParseDrive entered - Enter log tau\n" ); 00135 while( !lgEOL ) 00136 { 00137 if( read_whole_line( chInput , (int)sizeof(chInput) , ioStdin ) == NULL ) 00138 { 00139 fprintf( ioQQQ, " ParseDrive.dat error getting magic number\n" ); 00140 cdEXIT(EXIT_FAILURE); 00141 } 00142 /* get tau */ 00143 i = 1; 00144 fac = FFmtRead(chInput,&i,INPUT_LINE_LENGTH,&lgEOL); 00145 if( lgEOL ) 00146 break; 00147 fac = pow(10.,fac); 00148 fprintf( ioQQQ, " Tau =%11.4e\n", fac ); 00149 TauDummy.Emis->TauIn = (realnum)fac; 00150 fac = DrvContPump(&TauDummy); 00151 fprintf( ioQQQ, " ContPump =%11.4e\n", fac ); 00152 fprintf( ioQQQ, " Enter null to stop, or another value.\n" ); 00153 } 00154 fprintf( ioQQQ, " ContPump ParseDrive exits.\n" ); 00155 } 00156 00157 else if( nMatch("STAR",chCard) ) 00158 { 00159 char chInput[INPUT_LINE_LENGTH]; 00160 /* get starburst abundances */ 00161 for( i=0; i < 40; i++ ) 00162 { 00163 zed = ((double)i+1.)/4. + 0.01; 00164 sprintf( chInput, "starburst, zed=%10.4f", zed ); 00165 abund_starburst(chInput); 00166 fprintf( ioQQQ, "%8.1e", zed ); 00167 for(n=0; n < LIMELM; n++) 00168 fprintf( ioQQQ, "%8.1e", abund.solar[n] ); 00169 fprintf( ioQQQ, "\n" ); 00170 } 00171 } 00172 00173 else 00174 { 00175 fprintf( ioQQQ, 00176 " Unrecognized key; keys are FFmtRead, CASE, GAUNT, HYAS, POINT, MOLE, STAR, " 00177 "PUMP, and ESCApe. Sorry.\n" ); 00178 cdEXIT(EXIT_FAILURE); 00179 } 00180 return; 00181 } 00182 00183 /*DrvEscP user queries escape probability routines, which return values */ 00184 STATIC void DrvEscP( void ) 00185 { 00186 char chCard[INPUT_LINE_LENGTH]; 00187 bool lgEOL; 00188 long i; 00189 double tau; 00190 00191 DEBUG_ENTRY( "DrvEscP()" ); 00192 00193 /* this routine is enterd with the command escape probability, and 00194 * drives the escape probability routine to compare answers */ 00195 fprintf( ioQQQ, " Enter the log of the one-sided optical depth; line with no number to stop.\n" ); 00196 00197 lgEOL = false; 00198 while( !lgEOL ) 00199 { 00200 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL ) 00201 { 00202 break; 00203 } 00204 00205 i = 1; 00206 tau = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00207 if( lgEOL ) 00208 { 00209 break; 00210 } 00211 00212 tau = pow(10.,tau); 00213 fprintf( ioQQQ, "tau was %e\n", tau ); 00214 fprintf( ioQQQ, " ESCINC=%11.3e\n", esc_PRD_1side(tau,1e-4) ); 00215 fprintf( ioQQQ, " ESCCOM=%11.3e\n", esc_CRDwing_1side(tau,1e-4 ) ); 00216 fprintf( ioQQQ, " ESCA0K2=%11.3e\n", esca0k2(tau) ); 00217 00218 } 00219 return; 00220 } 00221 00222 /*DrvCaseBHS allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */ 00223 STATIC void DrvCaseBHS(void) 00224 { 00225 char chCard[INPUT_LINE_LENGTH]; 00226 bool lgEOL, 00227 lgHit; 00228 long int i, 00229 n1, 00230 nelem , 00231 n2; 00232 double Temperature, 00233 Density; 00234 00235 DEBUG_ENTRY( "DrvCaseBHS()" ); 00236 00237 /* this routine is entered with the command DRIVE CASEB, and 00238 * interpolates on the Hummer & Storey case b data set */ 00239 00240 /* read in some external data files, but only if this is first call */ 00241 fprintf(ioQQQ," I will get needed H data files. This will take a second.\n"); 00242 atmdat_readin(); 00243 00244 { 00245 /* following should be set true to print input lines */ 00246 /*@-redef@*/ 00247 enum {DEBUG_LOC=false}; 00248 /*@+redef@*/ 00249 if( DEBUG_LOC ) 00250 { 00251 double xLyman , alpha; 00252 long int ipHi; 00253 nelem = 2; 00254 Temperature = 2e4; 00255 Density = 1e2; 00256 for( ipHi=3; ipHi<25; ++ipHi ) 00257 { 00258 double photons = (1./POW2(ipHi-1.)-1./POW2((double)ipHi) ) /(1.-1./ipHi/ipHi ); 00259 xLyman = atmdat_HS_caseB(1,ipHi, nelem,Temperature , Density , 'A' ); 00260 alpha = atmdat_HS_caseB(ipHi-1,ipHi, nelem,Temperature , Density , 'A' ); 00261 fprintf(ioQQQ,"%li\t%.3e\t%.3e\n", ipHi, xLyman/alpha*photons, photons ); 00262 } 00263 cdEXIT(EXIT_SUCCESS); 00264 } 00265 } 00266 00267 /* first get the charge, only H and He at present */ 00268 lgHit = false; 00269 nelem = 0; 00270 while( !lgHit ) 00271 { 00272 fprintf( ioQQQ, " Enter atomic number of species, either 1(H) or 2(He).\n" ); 00273 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL ) 00274 { 00275 fprintf( ioQQQ, " error getting species \n" ); 00276 cdEXIT(EXIT_FAILURE); 00277 } 00278 00279 i = 1; 00280 nelem = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00281 if( lgEOL || nelem< 1 || nelem > 2 ) 00282 { 00283 fprintf( ioQQQ, " must be either 1 or 2!\n" ); 00284 } 00285 else 00286 { 00287 lgHit = true; 00288 } 00289 } 00290 00291 fprintf(ioQQQ," In the following temperatures <10 are log, >=10 linear.\n"); 00292 fprintf(ioQQQ," The density is always a log.\n"); 00293 fprintf(ioQQQ," The order of the quantum numbers do not matter.\n"); 00294 fprintf(ioQQQ," The smallest must not be smaller than 2,\n"); 00295 fprintf(ioQQQ," and the largest must not be larger than 25.\n"); 00296 fprintf(ioQQQ," Units of emissivity are erg cm^3 s^-1\n\n"); 00297 fprintf(ioQQQ," The limits of the HS tables are 2 <= n <= 25.\n"); 00298 00299 lgHit = true; 00300 /* this is always true */ 00301 while( lgHit ) 00302 { 00303 fprintf( ioQQQ, " Enter 4 numbers, temperature, density, 2 quantum numbers, null line stop.\n" ); 00304 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL ) 00305 { 00306 fprintf( ioQQQ, " Thanks for interpolating on the Hummer & Storey data set!\n" ); 00307 cdEXIT(EXIT_FAILURE); 00308 } 00309 00310 i = 1; 00311 Temperature = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00312 if( lgEOL ) 00313 { 00314 fprintf( ioQQQ, " error getting temperature!\n" ); 00315 break; 00316 } 00317 00318 /* log if less than 10 */ 00319 if( Temperature < 10. ) 00320 { 00321 Temperature = pow(10., Temperature ); 00322 } 00323 fprintf(ioQQQ," Temperature is %g\n", Temperature ); 00324 00325 Density = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00326 if( lgEOL ) 00327 { 00328 fprintf( ioQQQ, " error getting density!\n" ); 00329 break; 00330 } 00331 Density = pow(10., Density ); 00332 fprintf(ioQQQ," Density is %g\n", Density ); 00333 00334 /* these quantum numbers can be in any order */ 00335 n1 = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00336 if( lgEOL ) 00337 { 00338 fprintf( ioQQQ, " error getting quantum number!\n" ); 00339 break; 00340 } 00341 00342 n2 = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00343 if( lgEOL ) 00344 { 00345 fprintf( ioQQQ, " error getting quantum number!\n" ); 00346 break; 00347 } 00348 00349 if( MAX2( n1 , n2 ) > 25 ) 00350 { 00351 fprintf( ioQQQ," The limits of the HS tables are 2 <= n <= 25. Sorry.\n"); 00352 break; 00353 } 00354 00355 fprintf( ioQQQ, 00356 " 4pJ(%ld,%ld)/n_e n_p=%11.3e\n", 00357 n1, n2, 00358 atmdat_HS_caseB(n1,n2, nelem,Temperature , Density , 'B' ) ); 00359 00360 /* this is check that we were in bounds */ 00361 #if 0 00362 { 00363 long j; 00364 double tempTable[33] = { 00365 11870.,12490.,12820., 00366 11060.,17740.,12560., 00367 16390.,16700.,11360., 00368 10240.,20740.,12030., 00369 14450.,19510.,12550., 00370 16470.,16560.,12220., 00371 15820.,12960.,10190., 00372 12960.,14060.,12560., 00373 11030.,10770.,13360., 00374 10780.,11410.,11690., 00375 12500.,13190.,21120. }; 00376 double edenTable[33] = { 00377 10.,270.,80.,10.,70., 00378 110.,200.,10.,40.,90., 00379 340.,80.,60.,340.,30., 00380 120.,10.,50.,450.,30., 00381 180.,20.,170.,60.,20., 00382 40.,30.,20.,100.,130., 00383 10.,10.,110. }; 00384 00385 00386 for( j=0; j<33; j++ ) 00387 { 00388 double halpha, hbeta, hgamma; 00389 00390 halpha = atmdat_HS_caseB(2,3, 1,tempTable[j] , edenTable[j] , 'B' ); 00391 hbeta = atmdat_HS_caseB(2,4, 1,tempTable[j] , edenTable[j] , 'B' ); 00392 hgamma = atmdat_HS_caseB(2,5, 1,tempTable[j] , edenTable[j] , 'B' ); 00393 00394 fprintf( ioQQQ, "%e\t%e\t%e\t%e\n", 00395 tempTable[j], 00396 edenTable[j], 00397 halpha/hbeta, 00398 hgamma/hbeta ); 00399 } 00400 } 00401 #endif 00402 } 00403 00404 fprintf( ioQQQ, " Thanks for interpolating on the Hummer & Storey data set!\n" ); 00405 cdEXIT(EXIT_FAILURE); 00406 00407 } 00408 00409 /*DrvHyas allow user to query hydrogen A's, asks for up, low level, gives A, drive hyas */ 00410 STATIC void DrvHyas(void) 00411 { 00412 char chCard[INPUT_LINE_LENGTH]; 00413 bool lgEOL; 00414 long int i, nHi, lHi, nLo, lLo; 00415 00416 DEBUG_ENTRY( "DrvHyas()" ); 00417 00418 /* this routine is entered with the command DRIVE HYAS, and 00419 * drives Jason's hydrogen einstein A routines */ 00420 00421 nHi = 1; 00422 /* nHi never lt 1 */ 00423 while( nHi != 0 ) 00424 { 00425 fprintf( ioQQQ, " Enter four quantum numbers (n, l, n', l'), null line to stop.\n" ); 00426 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL ) 00427 { 00428 fprintf( ioQQQ, " error getting drvhyas \n" ); 00429 cdEXIT(EXIT_FAILURE); 00430 } 00431 00432 i = 1; 00433 nHi = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00434 if( lgEOL ) 00435 break; 00436 00437 lHi = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00438 if( lgEOL ) 00439 { 00440 fprintf( ioQQQ, " must be four numbers!\n" ); 00441 break; 00442 } 00443 00444 if( lHi >= nHi ) 00445 { 00446 fprintf( ioQQQ, " l must be less than n!\n" ); 00447 break; 00448 } 00449 00450 nLo = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00451 if( lgEOL ) 00452 { 00453 fprintf( ioQQQ, " must be four numbers!\n" ); 00454 break; 00455 } 00456 00457 lLo = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00458 if( lgEOL ) 00459 { 00460 fprintf( ioQQQ, " must be four numbers!\n" ); 00461 break; 00462 } 00463 00464 if( lLo >= nLo ) 00465 { 00466 fprintf( ioQQQ, " l must be less than n!\n" ); 00467 break; 00468 } 00469 00470 if( nLo > nHi ) 00471 { 00472 long nTemp, lTemp; 00473 00474 /* swap hi and lo */ 00475 nTemp = nLo; 00476 lTemp = lLo; 00477 nLo = nHi; 00478 lLo = lHi; 00479 nHi = nTemp; 00480 lHi = lTemp; 00481 } 00482 00483 fprintf( ioQQQ, " A(%3ld,%3ld->%3ld,%3ld)=%11.3e\n", 00484 nHi, lHi, nLo, lLo, 00485 H_Einstein_A( nHi, lHi, nLo, lLo, 1 ) ); 00486 00487 } 00488 fprintf( ioQQQ, " Driver exits, enter next line.\n" ); 00489 00490 return; 00491 } 00492 00493 /*dgaunt drive gaunt factor routines by letting user query values */ 00494 STATIC void dgaunt(void) 00495 { 00496 char chCard[INPUT_LINE_LENGTH]; 00497 bool lgEOL; 00498 int inputflag; 00499 long int i, 00500 ierror; 00501 realnum enerlin[1]; 00502 double SaveTemp; 00503 double z,mygaunt=0.; 00504 double loggamma2, logu; 00505 00506 DEBUG_ENTRY( "dgaunt()" ); 00507 00508 SaveTemp = phycon.te; 00509 00510 /* this routine is entered with the command DRIVE GAUNT, and 00511 * drives the gaunt factor routine to check range 00512 * */ 00513 fprintf( ioQQQ, " Enter 0 to input temp, energy, and net charge, or 1 for u and gamma**2.\n" ); 00514 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL ) 00515 { 00516 fprintf( ioQQQ, " dgaunt error getting magic number\n" ); 00517 cdEXIT(EXIT_FAILURE); 00518 } 00519 i = 1; 00520 inputflag = (int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00521 00522 if( inputflag == 0 ) 00523 { 00524 fprintf( ioQQQ, " Enter the temperature (log if <=10), energy (Ryd), and net charge. Null line to stop.\n" ); 00525 /* >>chng 96 july 07, got rid of statement labels replacing with do while 00526 * */ 00527 ierror = 0; 00528 while( ierror == 0 ) 00529 { 00530 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL ) 00531 { 00532 fprintf( ioQQQ, " dgaunt error getting magic number\n" ); 00533 cdEXIT(EXIT_FAILURE); 00534 } 00535 i = 1; 00536 phycon.alogte = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00537 /* the line may be trash but ierror will pick it up */ 00538 if( lgEOL ) 00539 { 00540 fprintf( ioQQQ, " Gaunt driver exits, enter next line.\n" ); 00541 break; 00542 } 00543 /* numbers less than or equal to 10 are the log of the temperature */ 00544 double TeNew; 00545 if( phycon.alogte > 10. ) 00546 { 00547 TeNew = phycon.alogte; 00548 } 00549 else 00550 { 00551 TeNew = pow(10.,phycon.alogte); 00552 } 00553 TempChange(TeNew , false); 00554 00555 enerlin[0] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00556 if( lgEOL || enerlin[0] == 0. ) 00557 { 00558 fprintf( ioQQQ, " Sorry, but there should be two more numbers, energy and charge.\n" ); 00559 } 00560 00561 z = (double)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00562 if( lgEOL || z == 0. ) 00563 { 00564 fprintf( ioQQQ, " Sorry, but there should be a third number, charge.\n" ); 00565 } 00566 00567 /* This is non-thermally averaged gaunt factors. */ 00568 mygaunt = cont_gaunt_calc( (double)phycon.te, z, enerlin[0] ); 00569 00570 fprintf( ioQQQ, " Using my routine, Gff= \t" ); 00571 fprintf( ioQQQ, "%11.3e\n", mygaunt ); 00572 00573 } 00574 } 00575 else 00576 { 00577 /* this routine is entered with the command DRIVE GAUNT, and 00578 * drives the gaunt factor routine to check range 00579 * */ 00580 fprintf( ioQQQ, " Enter log u and log gamma2. Null line to stop.\n" ); 00581 ierror = 0; 00582 while( ierror == 0 ) 00583 { 00584 if( read_whole_line( chCard , (int)sizeof(chCard) , ioStdin ) == NULL ) 00585 { 00586 fprintf( ioQQQ, " dgaunt error getting magic number\n" ); 00587 cdEXIT(EXIT_FAILURE); 00588 } 00589 i = 1; 00590 logu = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00591 /* the line may be trash but ierror will pick it up */ 00592 if( lgEOL ) 00593 { 00594 fprintf( ioQQQ, " Gaunt driver exits, enter next line.\n" ); 00595 break; 00596 } 00597 00598 loggamma2 = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00599 if( lgEOL ) 00600 { 00601 fprintf( ioQQQ, " Sorry, but there should be another numbers, log gamma2.\n" ); 00602 } 00603 00604 /* This is my attempt to calculate non-thermally averaged gaunt factors. */ 00605 mygaunt = cont_gaunt_calc( TE1RYD/pow(10.,loggamma2), 1., pow(10.,logu-loggamma2) ); 00606 00607 TempChange(TE1RYD/pow(10.,loggamma2) , false); 00608 00609 fprintf( ioQQQ, " Using my routine, Gaunt factor is" ); 00610 fprintf( ioQQQ, "%11.3e\n", mygaunt ); 00611 } 00612 } 00613 00614 TempChange(SaveTemp , false); 00615 return; 00616 }