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 /*radius_first derive thickness of first zone, called after conditions in first zone 00004 * are established, sets 00005 radius.drad_x_fillfac 00006 radius.drad 00007 */ 00008 #include "cddefines.h" 00009 #define Z 1.0001 00010 #include "wind.h" 00011 #include "stopcalc.h" 00012 #include "thermal.h" 00013 #include "dynamics.h" 00014 #include "trace.h" 00015 #include "punch.h" 00016 #include "pressure.h" 00017 #include "iso.h" 00018 #include "h2.h" 00019 #include "rfield.h" 00020 #include "dense.h" 00021 #include "hmi.h" 00022 #include "geometry.h" 00023 #include "opacity.h" 00024 #include "ipoint.h" 00025 #include "radius.h" 00026 00027 void radius_first(void) 00028 { 00029 long int i , 00030 ip; 00031 00032 bool lgDoPun; 00033 00034 int indexOfSmallest = 0; 00035 00036 const int NUM_DR_TYPES = 13; 00037 00038 struct t_drValues{ 00039 double dr; 00040 char whatToSay[40]; 00041 } drValues[NUM_DR_TYPES]; 00042 00043 double accel, 00044 BigOpacity, 00045 change, 00046 dr912, 00047 drH2 , 00048 drContPres , 00049 drOpacity , 00050 drStromgren, /* used for Stromgren length */ 00051 drTabDen, 00052 dradf, 00053 drcol, 00054 dr_time_dep, 00055 drthrm, 00056 factor, 00057 winddr; 00058 static double drad_last_iteration=-1.; 00059 00060 DEBUG_ENTRY( "radius_first()" ); 00061 00062 /*********************************************************************** 00063 * 00064 * wind model, use acceleration length 00065 * 00066 ***********************************************************************/ 00067 00068 if( wind.windv > 0. ) 00069 { 00070 /* evaluate total pressure, although value not used (so stuffed into dr912) */ 00071 /* >>chng 01 nov 02, remove call, confirm vals defined with assert */ 00072 ASSERT( dense.pden > 0. && dense.wmole > 0. ); 00073 accel = 1.3e-10*dense.pden*dense.wmole; 00074 winddr = POW2(wind.windv)/25./accel; 00075 } 00076 else 00077 { 00078 winddr = 1e30; 00079 } 00080 00081 /* key off of Lyman continuum optical depth */ 00082 if( StopCalc.taunu > 0.99 && StopCalc.taunu < 3. ) 00083 { 00084 dr912 = StopCalc.tauend/6.3e-18/(dense.xIonDense[ipHYDROGEN][0]*geometry.FillFac)*Z/50.; 00085 } 00086 else 00087 { 00088 dr912 = 1e30; 00089 } 00090 00091 if( dynamics.lgStatic && iteration > 2 ) 00092 { 00093 /* when time dependent case on do not let dr change since current continuum 00094 * is not good indicator of conditions */ 00095 dr_time_dep = drad_last_iteration; 00096 } 00097 else 00098 { 00099 dr_time_dep = 1e30; 00100 } 00101 00102 /*********************************************************************** 00103 * 00104 * key off of column density; total, neutral, or ionized 00105 * 00106 ***********************************************************************/ 00107 00108 if( StopCalc.HColStop < 5e29 ) 00109 { 00110 /* this is useful for very thin columns, normally larger than 1st DR */ 00111 drcol = log10(StopCalc.HColStop) - log10(dense.gas_phase[ipHYDROGEN]*geometry.FillFac* 20.); 00112 } 00113 else if( StopCalc.colpls < 5e29 ) 00114 { 00115 /* ionized column density */ 00116 drcol = log10(StopCalc.colpls) - log10(dense.xIonDense[ipHYDROGEN][1]*geometry.FillFac* 20.); 00117 } 00118 else if( StopCalc.colnut < 5e29 ) 00119 { 00120 /* neutral column denisty */ 00121 drcol = log10(StopCalc.colnut) - log10(dense.xIonDense[ipHYDROGEN][0]*geometry.FillFac*50.); 00122 } 00123 else 00124 { 00125 /* not used */ 00126 drcol = 30.; 00127 } 00128 /* finally convert the drived column density to linear scale */ 00129 drcol = pow(10.,MIN2(35.,drcol)); 00130 00131 /*********************************************************************** 00132 * 00133 * key off of density or abundance fluctuations, must be small part of wavelength 00134 * 00135 ***********************************************************************/ 00136 00137 if( dense.flong != 0. ) 00138 { 00139 /* flong set => density fluctuations */ 00140 dradf = 6.283/dense.flong/10.; 00141 dradf = MIN4(dradf,radius.router[iteration-1]*Z,drcol,dr912); 00142 } 00143 else 00144 { 00145 dradf = FLT_MAX; 00146 } 00147 00148 /* >>>chng 99 nov 18, add check on stromgren length */ 00149 /* estimate Stromgren length, but only if there are ionizing photons 00150 * and not constant temperature model */ 00151 if( (rfield.qhtot>0.) && (rfield.qhtot> rfield.qbal*0.01) && (rfield.uh>1e-10) ) 00152 { 00153 /* >>chng 99 dec 23, double to allow lte.in to work on alphas */ 00154 /* >>chng 03 mar 15, density to double to avoid overflow, PvH */ 00155 drStromgren = (double)(rfield.qhtot)/iso.RadRec_caseB[ipH_LIKE][ipHYDROGEN]/ 00156 POW2((double)dense.gas_phase[ipHYDROGEN]); 00157 00158 /* different logic if this is a sphere */ 00159 if( drStromgren/radius.rinner > 1. ) 00160 { 00161 /* >>chng 03 mar 15, to double to avoid FP overflow, PvH */ 00162 drStromgren = (double)rfield.qhtot*3./(radius.rinner* 00163 iso.RadRec_caseB[ipH_LIKE][ipHYDROGEN]*POW2((double)dense.gas_phase[ipHYDROGEN]) ); 00164 drStromgren += 1.; 00165 /* this results in r_out / r_in */ 00166 drStromgren = pow( drStromgren , 0.33333); 00167 /* make it a physics thickness in cm */ 00168 drStromgren *= radius.rinner; 00169 } 00170 00171 /* remember the Stromgren thickness */ 00172 radius.thickness_stromgren = (realnum)drStromgren; 00173 00174 /* take one hundredth of this */ 00175 drStromgren /= 100.; 00176 } 00177 else 00178 { 00179 drStromgren = FLT_MAX; 00180 radius.thickness_stromgren = FLT_MAX; 00181 } 00182 00183 /*********************************************************************** 00184 * 00185 * find largest opacity, to keep the first zone optical depth 1 00186 * this is usually the physics that sets the first zone thickness 00187 * 00188 ***********************************************************************/ 00189 00190 /* >>>chng 99 jun 25, this is to simulate behavior of code before extension 00191 * of continuum array to 1e-8 Ryd */ 00192 ip = ipoint(1e-5); 00193 00194 /* find largest opacity */ 00195 BigOpacity = 0.; 00196 for( i=ip; i < rfield.nflux; i++ ) 00197 { 00198 /* remember largest opacity, and energy where this happened, 00199 * make sure flux at energy is gt 0, can be zero for case where 00200 * nflux increased to include emission from some ions */ 00201 if( rfield.flux[i]>0. && opac.opacity_abs[i] > BigOpacity ) 00202 { 00203 BigOpacity = opac.opacity_abs[i]; 00204 } 00205 } 00206 /* BigOpacity may be zero on very first call */ 00207 00208 /* drChange set with set didz command, is only number set with this command, 00209 * default in zerologic is 0.15 00210 * set drad to small part of*/ 00211 if( BigOpacity > SMALLFLOAT ) 00212 { 00213 drOpacity = (radius.drChange/100.)/BigOpacity/geometry.FillFac; 00214 } 00215 else 00216 { 00217 drOpacity = 1e30; 00218 } 00219 00220 /*********************************************************************** 00221 * 00222 * thermalization length of typical lines 00223 * 00224 ***********************************************************************/ 00225 00226 drthrm = 1.5e31/MAX2(1.,POW2((double)dense.gas_phase[ipHYDROGEN])); 00227 00228 /*********************************************************************** 00229 * 00230 * make sure we resolve initial structure in dense_tabden command 00231 * if interpolated table we need to make sure that we resolve the 00232 * initial changes in the structure 00233 * 00234 ***********************************************************************/ 00235 00236 if( strcmp(dense.chDenseLaw,"DLW2") == 0 ) 00237 { 00238 drTabDen = 1.; 00239 i = 1; 00240 factor = 0.; 00241 while( i < 100 && factor < 0.05 ) 00242 { 00243 /* check densities at ever larger dr's, until factor becomes more than 5% */ 00244 factor = dense.gas_phase[ipHYDROGEN]/ 00245 dense_tabden(radius.Radius+drTabDen, drTabDen ); 00246 /* density change can be positive or negative sign */ 00247 factor = fabs(factor-1.); 00248 drTabDen *= 2.; 00249 i += 1; 00250 } 00251 drTabDen /= 2.; 00252 } 00253 else 00254 { 00255 drTabDen = 1e30; 00256 } 00257 00258 /* >>chng 03 mar 20, add check on lyman band optical depth - want first zone 00259 * to be thin in H2 bands */ 00260 /* some tests are fully molecular with solomon process turned off, 00261 * do not sense this when already almost fully molecular */ 00262 if( hmi.H2_total/dense.gas_phase[ipHYDROGEN] < 0.1 ) 00263 { 00264 change = 0.1; 00265 } 00266 else 00267 { 00268 /* >>chng 04 mar 14, this branch, H is quite molecular, 00269 * still do not want large changes in solomon rate since linearization 00270 * would not work in hmole network, bu do not need such fine steps */ 00271 change = 1.; 00272 } 00273 00274 /* >>chng 04 mar 14 go back to original logic since molecular 00275 * pdr's had big jump in conditions from 00276 * first to second zon even when most H in H2 00277 change = 0.1; */ 00278 /* >>chng 04 apr 18, change from 0.1 to 0.001, inital zones too large in 00279 * leiden test case f1 */ 00280 change = 0.001; 00281 /* >>chng 04 mar 13, not too large when big H2 is on */ 00282 if( h2.lgH2ON && hmi.lgBigH2_evaluated ) 00283 { 00284 if( fabs(hmi.HeatH2Dexc_BigH2)/thermal.ctot > 0.05 ) 00285 { 00286 /* changes in H2 heating caused by changes in solomon rate 00287 * would drive temperature failures */ 00288 /* >>chng 04 apr 18, change from 0.001 to 0.0001, inital zones too large in 00289 * leiden test case f1 */ 00290 change = 0.0001; 00291 } 00292 else 00293 { 00294 /* >>chng 04 apr 18, change from 0.01 to 0.001, inital zones too large in 00295 * leiden test case f1 */ 00296 change = 0.001; 00297 } 00298 } 00299 drH2 = change / SDIV( 00300 hmi.H2_total * geometry.FillFac * hmi.H2Opacity ); 00301 00302 /* >>chng 06 feb 01, very high U ulirg models had dramatic increase in 00303 * cont pre in first few zones, 00304 * in constant total pressure case, don't want acceleration across first zone to 00305 * be large compared with current gas pressure */ 00306 if( (strcmp( dense.chDenseLaw, "CPRE" )==0) && pressure.lgContRadPresOn ) 00307 { 00308 /* radiative acceleration was evaluated in PressureTotal */ 00309 drContPres = 0.05 * pressure.PresTotlCurr / 00310 (wind.AccelTot*dense.xMassDensity*geometry.FillFac*geometry.DirectionalCosin); 00311 } 00312 else if( wind.windv != 0. ) 00313 { 00314 /* acceleration and change in v in wind */ 00315 double g = fabs(wind.AccelTot-wind.AccelGravity); 00316 /* wind - do not let velocity change by too much */ 00317 drContPres = 0.05*POW2(wind.windv)/(2.*SDIV(g)); 00318 } 00319 else 00320 drContPres = 1e30; 00321 00322 drValues[0].dr = drOpacity; 00323 drValues[1].dr = radius.Radius/20.; 00324 drValues[2].dr = drStromgren; 00325 drValues[3].dr = radius.router[iteration-1]/10.; 00326 drValues[4].dr = drcol; 00327 drValues[5].dr = dr912; 00328 drValues[6].dr = drthrm; 00329 drValues[7].dr = winddr; 00330 drValues[8].dr = dradf; 00331 drValues[9].dr = drTabDen; 00332 drValues[10].dr = drH2; 00333 drValues[11].dr = drContPres; 00334 drValues[12].dr = dr_time_dep; 00335 00336 strcpy( drValues[0].whatToSay, "drOpacity" ); 00337 strcpy( drValues[1].whatToSay, "radius.Radius/20."); 00338 strcpy( drValues[2].whatToSay, "drStromgren"); 00339 strcpy( drValues[3].whatToSay, "radius.router[iteration-1]/10."); 00340 strcpy( drValues[4].whatToSay, "drcol"); 00341 strcpy( drValues[5].whatToSay, "dr912"); 00342 strcpy( drValues[6].whatToSay, "drthrm"); 00343 strcpy( drValues[7].whatToSay, "winddr"); 00344 strcpy( drValues[8].whatToSay, "dradf"); 00345 strcpy( drValues[9].whatToSay, "drTabDen"); 00346 strcpy( drValues[10].whatToSay, "drH2"); 00347 strcpy( drValues[11].whatToSay, "drContPres"); 00348 strcpy( drValues[12].whatToSay, "dr_time_dep"); 00349 00350 for( i=0; i<NUM_DR_TYPES; i++ ) 00351 { 00352 if( drValues[i].dr < drValues[indexOfSmallest].dr ) 00353 { 00354 indexOfSmallest = i; 00355 } 00356 } 00357 00358 radius.drad = drValues[indexOfSmallest].dr; 00359 00360 /* reset if radius.drad is less than radius.sdrmin */ 00361 if( radius.sdrmin >= radius.drad ) 00362 { 00363 radius.drad = radius.sdrmin; 00364 /* set flag for comment if the previous line forced a larger dr than 00365 * would otherwise have been chosen. will cause comment to be generated 00366 * in PrtComment if set true*/ 00367 radius.lgDR2Big = true; 00368 } 00369 else 00370 { 00371 radius.lgDR2Big = false; 00372 } 00373 00374 /* this min had been in the big min set above, but caused a false alarm 00375 * on the lgDR2Big test above since the set dr command sets both in and max */ 00376 radius.drad = MIN2( radius.sdrmax , radius.drad ); 00377 00378 #if 0 00379 /*********************************************************************** 00380 * 00381 * we have now generated range of estimates of first thickness, 00382 * now choose smallest of the group 00383 * 00384 ***********************************************************************/ 00385 00386 /* radius div by 20, to prevent big change in iron ionization for high ioniz gas, 00387 * this is also the ONLY place that sphericity comes in */ 00388 radius.drad = MIN4( MIN3( drOpacity, radius.Radius/20., drStromgren ), 00389 MIN3( radius.router[iteration-1]/10., drcol, dr912 ), 00390 MIN4( drthrm, winddr, dradf, drTabDen ), 00391 MIN3( drH2, drContPres, dr_time_dep ) ); 00392 00393 /* option to set lower limit to zone thickness, with set drmin command*/ 00394 radius.drad = MAX2( radius.drad , radius.sdrmin ); 00395 00396 /* set flag for comment if the previous line forced a larger dr than 00397 * would otherwise have been chosen. will cause comment to be generated 00398 * in PrtComment if set true*/ 00399 if( fp_equal( radius.drad, radius.sdrmin ) ) 00400 { 00401 radius.lgDR2Big = true; 00402 } 00403 else 00404 { 00405 radius.lgDR2Big = false; 00406 } 00407 00408 /* this min had been in the big min set above, but caused a false alarm 00409 * on the lgDR2Big test above since the set dr command sets both in and max */ 00410 radius.drad = MIN2( radius.sdrmax , radius.drad ); 00411 #endif 00412 00413 radius.drad_x_fillfac = radius.drad * geometry.FillFac; 00414 00415 /* save dr for this iteration */ 00416 drad_last_iteration = radius.drad; 00417 00418 /* drMinimum is smallest acceptable DRAD, and is 1/100 OF DRAD(1) */ 00419 /* this can be turned off by GLOB command */ 00420 if( radius.lgDrMnOn ) 00421 { 00422 /* >>chng 05 mar 05, drMinimum is now drad * hden, to make propro to optical depth 00423 * avoid false trigger across thermal fronts 00424 * add * dense.gas_phase */ 00425 /* NB - drMinimum not used in code - delete? */ 00426 radius.drMinimum = (realnum)(radius.drad * dense.gas_phase[ipHYDROGEN]/1e7); 00427 } 00428 else 00429 { 00430 radius.drMinimum = 0.; 00431 } 00432 00433 /* if set drmin is used, make sure drMinimum (which will cause an abort) is 00434 * smaller than drmin */ 00435 if( radius.lgSMinON ) 00436 { 00437 /* >>chng 05 mar 05, drMinimum is now drad * hden, to make propro to optical depth 00438 * avoid false trigger across thermal fronts 00439 * add * dense.gas_phase */ 00440 /* NB - drMinimum not used in code - delete? */ 00441 radius.drMinimum = MIN2(radius.drMinimum * dense.gas_phase[ipHYDROGEN],(realnum)(radius.sdrmin/10.f) ); 00442 } 00443 00444 if( trace.lgTrace ) 00445 { 00446 fprintf( ioQQQ, 00447 " radius_first called, finds dr=%13.5e drMinimum=%12.3e sdrmin=%10.2e sdrmax=%10.2e\n", 00448 radius.drad, radius.drMinimum/ dense.gas_phase[ipHYDROGEN], radius.sdrmin, radius.sdrmax ); 00449 } 00450 00451 if( radius.drad < SMALLFLOAT*1.1 ) 00452 { 00453 fprintf( ioQQQ, 00454 " PROBLEM radius_first detected likely insanity, found dr=%13.5e \n", radius.drad); 00455 fprintf( ioQQQ, 00456 " radius_first: calculation continuing but crash is likely. \n"); 00457 /* this sets flag that insanity has occurred */ 00458 TotalInsanity(); 00459 } 00460 00461 /* all this is to only punch on last iteration 00462 * the punch dr command is not really a punch command, making this necessary 00463 * lgDRon is set true if "punch dr" entered */ 00464 if( punch.lgDROn ) 00465 { 00466 lgDoPun = true; 00467 } 00468 else 00469 { 00470 lgDoPun = false; 00471 } 00472 00473 /* punch what we decided up? */ 00474 if( lgDoPun ) 00475 { 00476 /* create hash marks on second and later iterations */ 00477 if( iteration > 1 && punch.lgDRHash ) 00478 { 00479 static int iter_punch=-1; 00480 if( iteration !=iter_punch ) 00481 fprintf( punch.ipDRout, "%s\n",punch.chHashString ); 00482 iter_punch = iteration; 00483 } 00484 /* this is common part of each line, the zone count, depth, chosen dr, and depth2go */ 00485 /* >>chng 05 aug 15, had printed drNext, always zero, rather the drad, which is set here */ 00486 fprintf( punch.ipDRout , "%ld\t%.5e\t%.3e\t%.3e\t", nzone, radius.depth, radius.drad, radius.Depth2Go ); 00487 00488 if( radius.lgDR2Big ) 00489 { 00490 fprintf( punch.ipDRout, 00491 "radius_first keys from radius.sdrmin\n"); 00492 00493 } 00494 else if( fp_equal( radius.drad, radius.sdrmax ) ) 00495 { 00496 00497 fprintf( punch.ipDRout, 00498 "radius_first keys from radius.sdrmax\n"); 00499 } 00500 else 00501 { 00502 ASSERT( indexOfSmallest < NUM_DR_TYPES - 1 ); 00503 fprintf( punch.ipDRout, "radius_first keys from %s\n", 00504 drValues[indexOfSmallest].whatToSay); 00505 } 00506 00507 /* \todo 1 improve this printout and the drValues treatment above. */ 00508 00509 #if 0 00510 if( fp_equal( radius.drad, drOpacity ) ) 00511 { 00512 fprintf( punch.ipDRout, 00513 "radius_first keys from drOpacity, opac was %.2e at %.2e Ryd\n", 00514 BigOpacity , BigOpacityAnu ); 00515 } 00516 else if( fp_equal( radius.drad, radius.Radius/20. ) ) 00517 { 00518 fprintf( punch.ipDRout, 00519 "radius_first keys from radius.Radius\n" ); 00520 } 00521 else if( fp_equal( radius.drad, drStromgren ) ) 00522 { 00523 fprintf( punch.ipDRout, 00524 "radius_first keys from drStromgren\n"); 00525 } 00526 else if( fp_equal( radius.drad, dr_time_dep ) ) 00527 { 00528 fprintf( punch.ipDRout, 00529 "radius_first keys from time dependent\n"); 00530 } 00531 else if( fp_equal( radius.drad, radius.router[iteration-1]/10. ) ) 00532 { 00533 fprintf( punch.ipDRout, 00534 "radius_first keys from radius.router[iteration-1]\n"); 00535 } 00536 else if( fp_equal( radius.drad, drcol ) ) 00537 { 00538 fprintf( punch.ipDRout, 00539 "radius_first keys from drcol\n"); 00540 } 00541 else if( fp_equal( radius.drad, radius.sdrmin ) ) 00542 { 00543 fprintf( punch.ipDRout, 00544 "radius_first keys from radius.sdrmin\n"); 00545 } 00546 else if( fp_equal( radius.drad, dr912 ) ) 00547 { 00548 fprintf( punch.ipDRout, 00549 "radius_first keys from dr912\n"); 00550 } 00551 else if( fp_equal( radius.drad, radius.sdrmax ) ) 00552 { 00553 fprintf( punch.ipDRout, 00554 "radius_first keys from radius.sdrmax\n"); 00555 } 00556 else if( fp_equal( radius.drad, drthrm ) ) 00557 { 00558 fprintf( punch.ipDRout, 00559 "radius_first keys from drthrm\n"); 00560 } 00561 else if( fp_equal( radius.drad, winddr ) ) 00562 { 00563 fprintf( punch.ipDRout, 00564 "radius_first keys from winddr\n"); 00565 } 00566 else if( fp_equal( radius.drad, drH2 ) ) 00567 { 00568 fprintf( punch.ipDRout, 00569 "radius_first keys from H2 lyman lines\n"); 00570 } 00571 else if( fp_equal( radius.drad, dradf ) ) 00572 { 00573 fprintf( punch.ipDRout, 00574 "radius_first keys from dradf\n"); 00575 } 00576 else if( fp_equal( radius.drad, drTabDen ) ) 00577 { 00578 fprintf( punch.ipDRout, 00579 "radius_first keys from drTabDen\n"); 00580 } 00581 else if( fp_equal( radius.drad, drContPres ) ) 00582 { 00583 fprintf( punch.ipDRout, 00584 "radius_first keys from radiative acceleration across zone\n"); 00585 } 00586 else 00587 { 00588 fprintf( punch.ipDRout, "radius_first insanity\n" ); 00589 fprintf( ioQQQ, "radius_first insanity, radius is %e\n" , 00590 radius.drad); 00591 ShowMe(); 00592 } 00593 #endif 00594 00595 } 00596 return; 00597 }