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 /*PresTotCurrent determine the gas and line radiation pressures for current conditions, 00004 * this sets values of pressure.PresTotlCurr, also calls tfidle */ 00005 #include "cddefines.h" 00006 #include "physconst.h" 00007 #include "taulines.h" 00008 #include "opacity.h" 00009 #include "hextra.h" 00010 #include "elementnames.h" 00011 #include "hydrogenic.h" 00012 #include "conv.h" 00013 #include "iso.h" 00014 #include "wind.h" 00015 #include "magnetic.h" 00016 #include "doppvel.h" 00017 #include "rfield.h" 00018 #include "phycon.h" 00019 #include "thermal.h" 00020 #include "hmi.h" 00021 #include "h2.h" 00022 #include "dense.h" 00023 #include "atomfeii.h" 00024 #include "mole.h" 00025 #include "dynamics.h" 00026 #include "trace.h" 00027 #include "rt.h" 00028 #include "atmdat.h" 00029 #include "lines_service.h" 00030 #include "pressure.h" 00031 #include "radius.h" 00032 00033 /* this sets values of pressure.PresTotlCurr, also calls tfidle */ 00034 void PresTotCurrent(void) 00035 { 00036 long int i, 00037 ion, 00038 ipISO , 00039 ipHi, 00040 ipLo, 00041 nelem; 00042 00043 static long int 00044 /* used in debug print statement to see which line adds most pressure */ 00045 ipLineTypePradMax=-1 , 00046 /* points to line causing radiation pressure */ 00047 ipLinePradMax=-LONG_MAX, 00048 ip2=-LONG_MAX, 00049 ip3=-LONG_MAX, 00050 ip4=-LONG_MAX; 00051 00052 /* the line radiation pressure variables that must be preserved since 00053 * a particular line radiation pressure may not be evaluated if it is 00054 * not important */ 00055 static double Piso_seq[NISO]={0.,0.}, 00056 PLevel1=0., 00057 PLevel2=0., 00058 PHfs=0., 00059 PCO=0., 00060 P_H2=0., 00061 P_FeII=0.; 00062 00063 double 00064 RadPres1, 00065 TotalPressure_v, 00066 pmx; 00067 realnum den_hmole , 00068 den_comole , 00069 DenseAtomsIons, 00070 DensElements; 00071 00072 DEBUG_ENTRY( "PresTotCurrent()" ); 00073 00074 if( !conv.nTotalIoniz ) 00075 { 00076 /* resetting ipLinePradMax, necessary for 00077 * multiple cloudy calls in single coreload. */ 00078 ipLinePradMax=-LONG_MAX; 00079 //pressure.PresTotlCurr = 0.; 00080 } 00081 00082 /* PresTotCurrent - evaluate total pressure, dyne cm^-2 00083 * and radiative acceleration */ 00084 00085 /* several loops over the emission lines for radiation pressure and 00086 * radiative acceleration are expensive due to the number of lines and 00087 * the memory they occupy. Many equations of state do not include 00088 * radiation pressure or radiative acceleration. Only update these 00089 * when included in EOS - when not included only evaluate once per zone. 00090 * do it once per zone since we will still report these quantities. 00091 * this flag says whether we must update all terms */ 00092 bool lgMustEvaluate = false; 00093 00094 /* this says we already have an evaluation in this zone so do not 00095 * evaluate terms known to be trivial, even when reevaluating major terms */ 00096 bool lgMustEvaluateTrivial = false; 00097 /* if an individual agent is larger than this fraction of the total current 00098 * radiation pressure then it is not trivial 00099 * conv.PressureErrorAllowed is relative error allowed in pressure */ 00100 double TrivialLineRadiationPressure = conv.PressureErrorAllowed/10. * 00101 pressure.pres_radiation_lines_curr; 00102 00103 /* reevaluate during search phase when conditions are changing dramatically */ 00104 if( conv.lgSearch ) 00105 { 00106 lgMustEvaluate = true; 00107 lgMustEvaluateTrivial = true; 00108 } 00109 00110 /* reevaluate if zone or iteration has changed */ 00111 static long int nzoneEvaluated=0, iterationEvaluated=0; 00112 if( nzone!=nzoneEvaluated || iteration!=iterationEvaluated ) 00113 { 00114 lgMustEvaluate = true; 00115 lgMustEvaluateTrivial = true; 00116 /* this flag says which source of radiation pressure was strongest */ 00117 ipLineTypePradMax = -1; 00118 } 00119 00120 /* constant pressure or dynamical sim - we must reevaluate terms 00121 * but do not need to reevaluate trivial contributors */ 00122 if( (strcmp(dense.chDenseLaw,"WIND") == 0 ) || 00123 (strcmp(dense.chDenseLaw,"CPRE") == 0 ) ) 00124 lgMustEvaluate = true; 00125 00126 if( lgMustEvaluate ) 00127 { 00128 /* we are reevaluating quantities in this zone and iteration, 00129 * remember that we did it */ 00130 nzoneEvaluated = nzone; 00131 iterationEvaluated = iteration; 00132 } 00133 00134 /* update density - temperature variables which may have changed */ 00135 /* must call TempChange since ionization has changed, there are some 00136 * terms that affect collision rates (H0 term in electron collision) */ 00137 TempChange(phycon.te , false); 00138 00139 /* find total number of atoms and ions */ 00140 DenseAtomsIons = 0.; 00141 DensElements = 0.; 00142 realnum factor = 1.1f; 00143 if( conv.lgSearch ) 00144 factor = 2.f; 00145 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 00146 { 00147 /* only check element solution if it is on, and ionization 00148 * has not been set with TABLE ELEMENT command */ 00149 if( dense.lgElmtOn[nelem] ) 00150 { 00151 /* gas_phase is density of all atoms and ions, but not molecules */ 00152 DensElements += dense.gas_phase[nelem]; 00153 realnum DenseOne = 0; 00154 for( ion=0; ion<=nelem+1; ++ion ) 00155 DenseOne += dense.xIonDense[nelem][ion]; 00156 00157 /* during search phase sums may not add up correctly, and during 00158 * early parts of search phase may be badly off. This test 00159 * was introduced as result off fpe due to another 00160 * reason - Te had changed too much during initial search 00161 * for sim in which chemistry was important - fix was to not 00162 * let Te change. During resulting insanity caused by large 00163 * change, linearization did not work, CO and ionization solvers 00164 * could diverge leading to molecule or ionization population 00165 * larger than 1e38 limit to a realnum, fpe followed. this is to 00166 * catch that in its early stages */ 00167 if( conv.nTotalIoniz>5 && !dense.lgSetIoniz[nelem] && 00168 DenseOne > dense.gas_phase[nelem]*factor ) 00169 { 00170 /* species is not conserved */ 00171 fprintf(ioQQQ,"\n\n PROBLEM DISASTER PressureTotal: the chemical species " 00172 "are not conserved.\n"); 00173 fprintf(ioQQQ," The sum of the densities of the atoms and ions for " 00174 "%s is %.3e which is greater than the gas-phase abundance " 00175 "of the element, %.3e.\n", 00176 elementnames.chElementName[nelem],DenseOne, 00177 dense.gas_phase[nelem]); 00178 /* not including cosmic rays is the most likely cause of this */ 00179 if( hextra.cryden == 0. ) 00180 { 00181 fprintf(ioQQQ," The chemistry network is known to fail this way" 00182 " in cold molecular environments when cosmic rays are not" 00183 " included. They were not - try including them with the" 00184 " COSMIC RAY BACKBROUND command.\n"); 00185 fprintf(ioQQQ," The calculation is aborting. conv.nTotalIoniz=%li\n " 00186 "Sorry.\n\n", conv.nTotalIoniz ); 00187 lgAbort = true; 00188 } 00189 /* if they were included then something seriously wrong has happened*/ 00190 else 00191 TotalInsanity(); 00192 } 00193 DenseAtomsIons += DenseOne; 00194 } 00195 } 00196 /* DensElements is sum of all gas-phase densities of all elements, 00197 * DenseAtomsIons is sum of density of atoms and ions, does not 00198 * include molecules */ 00199 ASSERT( conv.nTotalIoniz<5 || DenseAtomsIons <= DensElements*factor ); 00200 ASSERT( DenseAtomsIons > 0. ); 00201 00202 /* evaluate the total ionization energy on a scale where neutral atoms 00203 * correspond to an energy of zero, so the ions have a positive energy */ 00204 phycon.EnergyIonization = 0.; 00205 #if 0 00206 for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ ) 00207 { 00208 for( ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion ) 00209 { 00210 /* lgMETALS is true, set false with "no advection metals" command */ 00211 int kadvec = dynamics.lgMETALS; 00212 /* this gives the iso sequence for this ion - should it be included in the 00213 * advected energy equation? lgISO is true, set false with 00214 * "no advection H-like" and "he-like" - for testing*/ 00215 ipISO = nelem-ion; 00216 fixit(); /* should this be kadvec = kadvec && dynamics.lgISO[ipISO]; ? */ 00217 if( ipISO >= 0 && ipISO<NISO ) 00218 kadvec = dynamics.lgISO[ipISO]; 00219 for( i=1; i<=ion; ++i ) 00220 { 00221 long int nelec = nelem+1 - i; 00222 /* this is the sum of all the energy needed to bring the atom up 00223 * to the ion+1 level of ionization - at this point a positive number */ 00224 phycon.EnergyIonization += dense.xIonDense[nelem][ion] * 00225 t_ADfA::Inst().ph1(Heavy.nsShells[nelem][i-1]-1,nelec,nelem,0)/EVRYD* 0.9998787*kadvec; 00226 } 00227 } 00228 } 00229 /* convert to ergs/cm^3 */ 00230 phycon.EnergyIonization *= EN1RYD; 00231 #endif 00232 00236 phycon.EnergyBinding = 0.; 00237 00238 /* find number of molecules of the heavy elements in the gas phase. 00239 * Atoms and ions were counted above. Do not include ices, solids 00240 * on grains */ 00241 den_comole = 0.; 00242 for( i=0; i < mole.num_comole_calc; i++ ) 00243 { 00244 /* term COmole[i]->lgGas_Phase is 0 or 1 depending on whether molecule 00245 * is on grain or in gas phase 00246 * n_nuclei is number of nuclei in molecule, this tests whether 00247 * actually an atom or molecule - atoms were already counted above */ 00248 if(COmole[i]->n_nuclei > 1) 00249 den_comole += COmole[i]->hevmol * COmole[i]->lgGas_Phase; 00250 } 00251 00252 /* number of hydrogen molecules, loop over all H molecular species, 00253 * do not include H0, H+ */ 00254 den_hmole = 0.; 00255 for( i=0; i<N_H_MOLEC; ++i ) 00256 { 00257 if( i!=ipMH && i!=ipMHp ) 00258 den_hmole += hmi.Hmolec[i]; 00259 } 00260 00261 /* total number of atoms, ions, and molecules in gas phase per unit vol */ 00262 dense.xNucleiTotal = den_hmole + den_comole + DenseAtomsIons; 00263 if( dense.xNucleiTotal > BIGFLOAT ) 00264 { 00265 fprintf(ioQQQ, "PROBLEM DISASTER pressure_total has found " 00266 "dense.xNucleiTotal with an insane density.\n"); 00267 fprintf(ioQQQ, "The density was %.2e\n", 00268 dense.xNucleiTotal); 00269 TotalInsanity(); 00270 } 00271 ASSERT( dense.xNucleiTotal > 0. ); 00272 00273 /* particle density that enters into the pressure includes electrons 00274 * cm-3 */ 00275 dense.pden = (realnum)(dense.eden + dense.xNucleiTotal); 00276 00277 /* the current gas pressure */ 00278 pressure.PresGasCurr = dense.pden*phycon.te*BOLTZMANN; 00279 /*fprintf(ioQQQ,"DEBUG gassss %.2f %.4e %.4e %.4e \n", 00280 fnzone, pressure.PresGasCurr , dense.pden , pressure.PresInteg );*/ 00281 00282 /* dense.wmole is molecular weight, AtomicMassUnit per particle */ 00283 dense.wmole = 0.; 00284 for( i=0; i < LIMELM; i++ ) 00285 { 00286 dense.wmole += dense.gas_phase[i]*(realnum)dense.AtomicWeight[i]; 00287 } 00288 00289 /* dense.wmole is now molecular weight, AtomicMassUnit per particle */ 00290 dense.wmole /= dense.pden; 00291 ASSERT( dense.wmole > 0. && dense.pden > 0. ); 00292 00293 /* xMassDensity is density in grams per cc */ 00296 dense.xMassDensity = (realnum)(dense.wmole*ATOMIC_MASS_UNIT*dense.pden); 00297 00298 /*>>chng 04 may 25, introduce following test on xMassDensity0 00299 * WJH 21 may 04, this is the mass density that corresponds to the hden 00300 * specified in the init file. It is used to calculate the mass flux in the 00301 * dynamics models. It may not necessarily be the same as struc.DenMass[0], 00302 * since the pressure solver may have jumped to a different density at the 00303 * illuminated face from that specified.*/ 00304 if( dense.xMassDensity0 < 0.0 ) 00305 dense.xMassDensity0 = dense.xMassDensity; 00306 00307 /* >>chng 01 nov 02, move to here from dynamics routines */ 00308 /* >>chng 02 jun 18 (rjrw), fix to use local values */ 00309 /* WJH 21 may 04, now recalculate wind v for the first zone too. 00310 * This is necessary when we are forcing the sub or supersonic branch */ 00311 if(wind.windv < 0) 00312 wind.windv = DynaFlux(radius.depth)/(dense.xMassDensity); 00313 00314 /* this is the current ram pressure, at this location */ 00315 pressure.PresRamCurr = dense.xMassDensity*POW2( wind.windv ); 00316 00317 /* this is the current turbulent pressure, not now added to equation of state 00318 * >>chng 06 mar 29, add Heiles_Troland_F / 6. term*/ 00319 pressure.PresTurbCurr = dense.xMassDensity*POW2( DoppVel.TurbVel ) * 00320 DoppVel.Heiles_Troland_F / 6.; 00321 00324 /* radiative acceleration, lines and continuum */ 00325 if( lgMustEvaluate ) 00326 { 00327 /* continuous radiative acceleration */ 00328 double rforce = 0.; 00329 for( i=0; i < rfield.nflux; i++ ) 00330 { 00331 rforce += (rfield.flux[i] + rfield.outlin[i] + rfield.outlin_noplot[i]+ rfield.ConInterOut[i])* 00332 rfield.anu[i]*(opac.opacity_abs[i] + opac.opacity_sct[i]); 00333 } 00334 00335 /* line acceleration; xMassDensity is gm per cc */ 00336 wind.AccelLine = (realnum)(RT_line_driving()/SPEEDLIGHT/dense.xMassDensity); 00337 wind.AccelCont = (realnum)(rforce*EN1RYD/SPEEDLIGHT/dense.xMassDensity); 00338 /* this is numerically unstable */ 00339 wind.AccelPres = 0.; 00340 /* total acceleration */ 00341 wind.AccelTot = wind.AccelCont + wind.AccelLine + wind.AccelPres; 00342 /* G, COMASS = mass of star in solar masses */ 00343 double reff = radius.Radius-radius.drad/2.; 00344 wind.AccelGravity = (realnum)(GRAV_CONST*wind.comass*SOLAR_MASS/POW2(reff)* 00345 (1.-wind.DiskRadius/reff) ); 00346 # if 0 00347 if( fudge(-1) ) 00348 fprintf(ioQQQ,"DEBUG pressure_total updates AccelTot to %.2e grav %.2e\n", 00349 wind.AccelTot , wind.AccelGravity ); 00350 # endif 00351 } 00352 00353 /* must always evaluate H La width since used elsewhere */ 00354 if( Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->PopOpc > 0. ) 00355 { 00356 hydro.HLineWidth = (realnum)RT_LineWidth(&Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s]); 00357 # if 0 00358 fprintf(ioQQQ,"DEBUG widLya\t%li\t%.2f\t%.3e", 00359 iteration, 00360 fnzone, 00361 hydro.HLineWidth); 00362 hydro.HLineWidth = (realnum)( 00363 RT_LyaWidth( 00364 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauIn, 00365 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauTot, 00366 4.72e-2/phycon.sqrte, 00367 DoppVel.doppler[ipHYDROGEN] ) ); 00368 fprintf(ioQQQ,"\t%.3e\n", 00369 hydro.HLineWidth); 00370 # endif 00371 } 00372 else 00373 hydro.HLineWidth = 0.; 00374 00375 00376 /* find max rad pressure even if capped 00377 * lgLineRadPresOn is turned off by NO RADIATION PRESSURE command */ 00378 if( trace.lgTrace ) 00379 { 00380 fprintf(ioQQQ, 00381 " PresTotCurrent nzone %li iteration %li lgMustEvaluate:%c " 00382 "lgMustEvaluateTrivial:%c " 00383 "pressure.lgLineRadPresOn:%c " 00384 "rfield.lgDoLineTrans:%c \n", 00385 nzone , iteration , TorF(lgMustEvaluate) , TorF(lgMustEvaluateTrivial), 00386 TorF(pressure.lgLineRadPresOn), TorF(rfield.lgDoLineTrans) ); 00387 } 00388 00389 if( lgMustEvaluate && pressure.lgLineRadPresOn && rfield.lgDoLineTrans ) 00390 { 00391 /* RadPres is pressure due to lines, lgPres_radiation_ON turns off or on */ 00392 pressure.pres_radiation_lines_curr = 0.; 00393 /* used to remember largest radiation pressure source */ 00394 pmx = 0.; 00395 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00396 { 00397 if( lgMustEvaluateTrivial || Piso_seq[ipISO] > TrivialLineRadiationPressure ) 00398 { 00399 Piso_seq[ipISO] = 0.; 00400 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00401 { 00402 /* does this ion stage exist? */ 00403 if( dense.IonHigh[nelem] >= nelem + 1 - ipISO ) 00404 { 00405 /* do not include highest levels since maser can occur due to topoff, 00406 * and pops are set to small number in this case */ 00407 for( ipHi=1; ipHi <iso.numLevels_local[ipISO][nelem] - iso.nCollapsed_local[ipISO][nelem]; ipHi++ ) 00408 { 00409 for( ipLo=0; ipLo < ipHi; ipLo++ ) 00410 { 00411 if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 ) 00412 continue; 00413 00414 ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul > iso.SmallA ); 00415 00416 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc > SMALLFLOAT && 00417 /* test that have not overrun optical depth scale */ 00418 ( (Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauTot - Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn) > SMALLFLOAT ) && 00419 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc > FLT_EPSILON*100. ) 00420 { 00421 RadPres1 = PressureRadiationLine( &Transitions[ipISO][nelem][ipHi][ipLo], dense.xIonDense[nelem][nelem+1-ipISO] ); 00422 00423 if( RadPres1 > pmx ) 00424 { 00425 ipLineTypePradMax = 2; 00426 pmx = RadPres1; 00427 ip4 = ipISO; 00428 ip3 = nelem; 00429 ip2 = ipHi; 00430 ipLinePradMax = ipLo; 00431 } 00432 Piso_seq[ipISO] += RadPres1; 00433 { 00434 /* option to print particulars of some line when called */ 00435 enum {DEBUG_LOC=false}; 00436 if( DEBUG_LOC && ipISO==ipH_LIKE && ipLo==3 && ipHi==5 && nzone > 260 ) 00437 { 00438 fprintf(ioQQQ, 00439 "DEBUG prad1 \t%.2f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t\n", 00440 fnzone, 00441 RadPres1, 00442 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc, 00443 StatesElem[ipISO][nelem][ipLo].Pop, 00444 StatesElem[ipISO][nelem][ipHi].Pop, 00445 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc); 00446 } 00447 } 00448 } 00449 } 00450 } 00451 } 00452 } 00453 ASSERT( Piso_seq[ipISO] >= 0. ); 00454 } 00455 pressure.pres_radiation_lines_curr += Piso_seq[ipISO]; 00456 } 00457 { 00458 /* option to print particulars of some line when called */ 00459 enum {DEBUG_LOC=false}; 00460 # if 0 00461 if( DEBUG_LOC /*&& iteration > 1*/ && nzone > 260 ) 00462 { 00463 fprintf(ioQQQ, 00464 "DEBUG prad2 \t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n", 00465 nzone, 00466 pmx, 00467 Transitions[ipISO][ip3][ipLinePradMax][ip2].Emis->PopOpc, 00468 StatesElem[ipISO][ip3][0].Pop, 00469 StatesElem[ipISO][ip3][2].Pop, 00470 StatesElem[ipISO][ip3][3].Pop, 00471 StatesElem[ipISO][ip3][4].Pop, 00472 StatesElem[ipISO][ip3][5].Pop, 00473 StatesElem[ipISO][ip3][6].Pop); 00474 } 00475 # endif 00476 if( DEBUG_LOC /*&& iteration > 1 && nzone > 150 */) 00477 { 00478 fprintf(ioQQQ, 00479 "DEBUG prad3\t%.2e\t%li\t%li\t%li\t%li\t%.2e\t%.2e\t%.2e\n", 00480 pmx, 00481 ip4, 00482 ip3, 00483 ip2, 00484 ipLinePradMax, 00485 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->PopOpc, 00486 StatesElem[ip4][ip3][ip2].Pop, 00487 1.-Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->Pesc ); 00488 } 00489 } 00490 00491 if( lgMustEvaluateTrivial || PLevel1 > TrivialLineRadiationPressure ) 00492 { 00493 PLevel1 = 0.; 00494 /* line radiation pressure from large set of level 1 lines */ 00495 for( i=1; i <= nLevel1; i++ ) 00496 { 00497 if( TauLines[i].Hi->Pop > 1e-30 ) 00498 { 00499 RadPres1 = PressureRadiationLine( &TauLines[i], 1. ); 00500 00501 if( RadPres1 > pmx ) 00502 { 00503 ipLineTypePradMax = 3; 00504 pmx = RadPres1; 00505 ipLinePradMax = i; 00506 } 00507 PLevel1 += RadPres1; 00508 } 00509 } 00510 ASSERT( PLevel1 >= 0. ); 00511 } 00512 pressure.pres_radiation_lines_curr += PLevel1; 00513 00514 if( lgMustEvaluateTrivial || PLevel2 > TrivialLineRadiationPressure ) 00515 { 00516 /* level 2 lines */ 00517 PLevel2 = 0.; 00518 for( i=0; i < nWindLine; i++ ) 00519 { 00520 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 00521 { 00522 if( TauLine2[i].Hi->Pop > 1e-30 ) 00523 { 00524 RadPres1 = PressureRadiationLine( &TauLine2[i], 1. ); 00525 00526 PLevel2 += RadPres1; 00527 if( RadPres1 > pmx ) 00528 { 00529 ipLineTypePradMax = 4; 00530 pmx = RadPres1; 00531 ipLinePradMax = i; 00532 } 00533 } 00534 } 00535 } 00536 ASSERT( PLevel2 >= 0. ); 00537 } 00538 pressure.pres_radiation_lines_curr += PLevel2; 00539 00540 /* fine structure lines */ 00541 if( lgMustEvaluateTrivial || PHfs > TrivialLineRadiationPressure ) 00542 { 00543 PHfs = 0.; 00544 for( i=0; i < nHFLines; i++ ) 00545 { 00546 if( HFLines[i].Hi->Pop > 1e-30 ) 00547 { 00548 RadPres1 = PressureRadiationLine( &HFLines[i], 1. ); 00549 00550 PHfs += RadPres1; 00551 if( RadPres1 > pmx ) 00552 { 00553 ipLineTypePradMax = 8; 00554 pmx = RadPres1; 00555 ipLinePradMax = i; 00556 } 00557 } 00558 } 00559 ASSERT( PHfs >= 0. ); 00560 } 00561 pressure.pres_radiation_lines_curr += PHfs; 00562 00563 /* radiation pressure due to H2 lines */ 00564 if( lgMustEvaluateTrivial || P_H2 > TrivialLineRadiationPressure ) 00565 { 00566 P_H2 = H2_RadPress(); 00567 /* flag to remember H2 radiation pressure */ 00568 if( P_H2 > pmx ) 00569 { 00570 pmx = P_H2; 00571 ipLineTypePradMax = 9; 00572 } 00573 ASSERT( P_H2 >= 0. ); 00574 } 00575 pressure.pres_radiation_lines_curr += P_H2; 00576 00577 /* radiation pressure due to FeII lines and large atom */ 00578 if( lgMustEvaluateTrivial || P_FeII > TrivialLineRadiationPressure ) 00579 { 00580 P_FeII = FeIIRadPress(); 00581 if( P_FeII > pmx ) 00582 { 00583 pmx = P_FeII; 00584 ipLineTypePradMax = 7; 00585 } 00586 ASSERT( P_FeII >= 0. ); 00587 } 00588 pressure.pres_radiation_lines_curr += P_FeII; 00589 00590 /* co carbon monoxide lines */ 00591 if( lgMustEvaluateTrivial || PCO > TrivialLineRadiationPressure ) 00592 { 00593 PCO = 0.; 00594 for( i=0; i < nCORotate; i++ ) 00595 { 00596 if( C12O16Rotate[i].Hi->Pop > 1e-30 ) 00597 { 00598 RadPres1 = PressureRadiationLine( &C12O16Rotate[i], 1. ); 00599 00600 PCO += RadPres1; 00601 if( RadPres1 > pmx ) 00602 { 00603 ipLineTypePradMax = 5; 00604 pmx = RadPres1; 00605 ipLinePradMax = i; 00606 } 00607 } 00608 if( C13O16Rotate[i].Hi->Pop > 1e-30 ) 00609 { 00610 RadPres1 = PressureRadiationLine( &C13O16Rotate[i], 1. ); 00611 00612 PCO += RadPres1; 00613 if( RadPres1 > pmx ) 00614 { 00615 ipLineTypePradMax = 6; 00616 pmx = RadPres1; 00617 ipLinePradMax = i; 00618 } 00619 } 00620 } 00621 ASSERT( PCO >= 0. ); 00622 } 00623 pressure.pres_radiation_lines_curr += PCO; 00624 } 00625 else if( !pressure.lgLineRadPresOn || !rfield.lgDoLineTrans ) 00626 { 00627 /* case where radiation pressure is not turned on */ 00628 ipLinePradMax = -1; 00629 ipLineTypePradMax = 0; 00630 } 00631 00632 /* the ratio of radiation to total (gas + continuum + etc) pressure */ 00633 pressure.pbeta = (realnum)(pressure.pres_radiation_lines_curr/SDIV(pressure.PresTotlCurr)); 00634 00635 /* this would be a major logical error */ 00636 if( pressure.pres_radiation_lines_curr < 0. ) 00637 { 00638 fprintf(ioQQQ, 00639 "PresTotCurrent: negative pressure, constituents follow %e %e %e %e %e \n", 00640 Piso_seq[ipH_LIKE], 00641 Piso_seq[ipHE_LIKE], 00642 PLevel1, 00643 PLevel2, 00644 PCO); 00645 ShowMe(); 00646 cdEXIT(EXIT_FAILURE); 00647 } 00648 00649 /* following test will never succeed, here to trick lint, ipLineTypePradMax is only used 00650 * when needed for debug */ 00651 if( trace.lgTrace && ipLineTypePradMax <0 ) 00652 { 00653 fprintf(ioQQQ, 00654 " PresTotCurrent, pressure.pbeta = %e, ipLineTypePradMax%li ipLinePradMax=%li \n", 00655 pressure.pbeta,ipLineTypePradMax,ipLinePradMax ); 00656 } 00657 00658 /* this is the total internal energy of the gas, the amount of energy needed 00659 * to produce the current level populations, relative to ground, 00660 * only used for energy terms in advection */ 00661 phycon.EnergyExcitation = 0.; 00662 #if 0 00663 fixit(); /* the quantities phycon.EnergyExcitation, phycon.EnergyIonization, phycon.EnergyBinding 00664 * are not used anywhere, except in print statements... */ 00665 broken(); /* the code below contains serious bugs! It is supposed to implement loops 00666 * over quantum states in order to evaluate the potential energy stored in 00667 * excited states of all atoms, ions, and molecules. The code below however 00668 * often implements loops over all combinations of upper and lower levels! 00669 * This code needs to be rewritten once quantumstates are fully implemented. */ 00670 ipLo = 0; 00671 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00672 { 00673 for( nelem=ipISO; nelem < LIMELM; nelem++ ) 00674 { 00675 if( dense.IonHigh[nelem] == nelem + 1 - ipISO ) 00676 { 00677 /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */ 00678 for( ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ ) 00679 { 00680 phycon.EnergyExcitation += 00681 Transitions[ipISO][nelem][ipHi][ipLo].Hi->Pop * 00682 Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg* 00683 /* last term is option to turn off energy term, no advection hlike, he-like */ 00684 dense.xIonDense[nelem][nelem+1-ipISO]*dynamics.lgISO[ipISO]; 00685 } 00686 } 00687 } 00688 } 00689 00690 if( dynamics.lgISO[ipH_LIKE] ) 00691 /* internal energy of H2 */ 00692 phycon.EnergyExcitation += H2_InterEnergy(); 00693 00694 /* this is option to turn off energy effects of advection, no advection metals */ 00695 if( dynamics.lgMETALS ) 00696 { 00697 for( i=1; i <= nLevel1; i++ ) 00698 { 00699 phycon.EnergyExcitation += 00700 TauLines[i].Hi->Pop* TauLines[i].EnergyErg; 00701 } 00702 for( i=0; i < nWindLine; i++ ) 00703 { 00704 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO ) 00705 { 00706 phycon.EnergyExcitation += 00707 TauLine2[i].Hi->Pop* TauLine2[i].EnergyErg; 00708 } 00709 } 00710 for( i=0; i < nHFLines; i++ ) 00711 { 00712 phycon.EnergyExcitation += 00713 HFLines[i].Hi->Pop* HFLines[i].EnergyErg; 00714 } 00715 double Energy12 = 0.; 00716 double Energy13 = 0.; 00717 for( i=0; i < nCORotate; i++ ) 00718 { 00719 Energy12 += C12O16Rotate[i].EnergyErg; 00720 phycon.EnergyExcitation += 00721 C12O16Rotate[i].Hi->Pop* Energy12; 00722 00723 Energy13 += C13O16Rotate[i].EnergyErg; 00724 phycon.EnergyExcitation += 00725 C13O16Rotate[i].Hi->Pop* Energy13; 00726 } 00727 00728 /* internal energy of large FeII atom */ 00729 phycon.EnergyExcitation += FeII_InterEnergy(); 00730 } 00731 #endif 00732 00733 /* ================================================================== 00734 * end internal energy of atoms and molecules */ 00735 00736 /* evaluate some parameters to do with magnetic field */ 00737 Magnetic_evaluate(); 00738 00739 /*lint -e644 Piso_seq not init */ 00740 if( trace.lgTrace && (pressure.PresTotlCurr > SMALLFLOAT) ) 00741 { 00742 fprintf(ioQQQ, 00743 " PresTotCurrent zn %.2f Ptot:%.5e Pgas:%.3e Prad:%.3e Pmag:%.3e Pram:%.3e " 00744 "gas parts; H:%.2e He:%.2e L1:%.2e L2:%.2e CO:%.2e fs%.2e h2%.2e turb%.2e\n", 00745 fnzone, 00746 pressure.PresTotlCurr, 00747 pressure.PresGasCurr/pressure.PresTotlCurr, 00748 pressure.pres_radiation_lines_curr*pressure.lgPres_radiation_ON/pressure.PresTotlCurr, 00749 magnetic.pressure/pressure.PresTotlCurr, 00750 pressure.PresRamCurr*pressure.lgPres_ram_ON/pressure.PresTotlCurr, 00751 Piso_seq[ipH_LIKE]/pressure.PresTotlCurr, 00752 Piso_seq[ipHE_LIKE]/pressure.PresTotlCurr, 00753 PLevel1/pressure.PresTotlCurr, 00754 PLevel2/pressure.PresTotlCurr, 00755 PCO/pressure.PresTotlCurr, 00756 PHfs/pressure.PresTotlCurr, 00757 P_H2/pressure.PresTotlCurr, 00758 pressure.PresTurbCurr*DoppVel.lgTurb_pressure/pressure.PresTotlCurr); 00759 /*lint +e644 Piso_seq not initialized */ 00760 } 00761 00762 /* Conserved quantity in steady-state energy equation for dynamics: 00763 * thermal energy, since recombination is treated as cooling 00764 * (i.e. is loss of electron k.e. to emitted photon), so don't 00765 * include 00766 * ...phycon.EnergyExcitation + phycon.EnergyIonization + phycon.EnergyBinding 00767 * */ 00768 00769 /* constant density case is also hypersonic case */ 00770 if( dynamics.lgStatic && dense.lgDenseInitConstant ) 00771 { 00772 /* this branch is time dependent AND constant density */ 00773 /*fprintf(ioQQQ,"DEBUG enthalpy HIT1\n");*/ 00774 /* this is the time-varying case where density is constant */ 00775 /* \todo 1 this has 3/2 on the PresGasCurr while true dynamics case below 00776 * has 5/2 - so this is not really the enthalpy density - better 00777 * would be to always use this term and add the extra PresGasCurr 00778 * when the enthalpy is actually needed */ 00779 phycon.EnthalpyDensity = 00780 0.5*POW2(wind.windv)*dense.xMassDensity + /* KE */ 00781 3./2.*pressure.PresGasCurr + /* thermal plus PdV work */ 00782 magnetic.EnthalpyDensity; /* magnetic terms */ 00783 } 00784 else 00785 { 00786 /* this branch either advective or constant pressure */ 00787 /*fprintf(ioQQQ,"DEBUG enthalpy HIT2\n");*/ 00788 /* this is usual dynamics case, or time-varying case where pressure 00789 * is kept constant and PdV work is added to the cell */ 00790 phycon.EnthalpyDensity = 00791 0.5*POW2(wind.windv)*dense.xMassDensity + /* KE */ 00792 5./2.*pressure.PresGasCurr + /* thermal plus PdV work */ 00793 magnetic.EnthalpyDensity; /* magnetic terms */ 00794 } 00795 00796 /* stop model from running away on first iteration, when line optical 00797 * depths are not defined correctly anyway. 00798 * if( iter.le.itermx .and. RadPres.ge.GasPres ) then 00799 * >>chng 97 jul 23, only cap radiation pressure on first iteration */ 00800 if( iteration <= 1 && pressure.pres_radiation_lines_curr >= pressure.PresGasCurr ) 00801 { 00802 /* stop RadPres from exceeding the gas pressure on first iteration */ 00803 pressure.pres_radiation_lines_curr = 00804 MIN2(pressure.pres_radiation_lines_curr,pressure.PresGasCurr); 00805 pressure.lgPradCap = true; 00806 } 00807 00808 /* remember the globally most important line, in the entire model 00809 * test on nzone so we only do test after solution is stable */ 00810 if( pressure.pbeta > pressure.RadBetaMax && nzone ) 00811 { 00812 pressure.RadBetaMax = pressure.pbeta; 00813 pressure.ipPradMax_line = ipLinePradMax; 00814 pressure.ipPradMax_nzone = nzone; 00815 if( ipLineTypePradMax == 2 ) 00816 { 00817 /* hydrogenic */ 00818 strcpy( pressure.chLineRadPres, "ISO "); 00819 ASSERT( ip4 < NISO && ip3<LIMELM ); 00820 ASSERT( ipLinePradMax>=0 && ip2>=0 && ip3>=0 && ip4>=0 ); 00821 strcat( pressure.chLineRadPres, chLineLbl(&Transitions[ip4][ip3][ip2][ipLinePradMax]) ); 00822 { 00823 /* option to print particulars of some line when called */ 00824 enum {DEBUG_LOC=false}; 00825 /*lint -e644 Piso_seq not initialized */ 00826 /* trace serious constituents in radiation pressure */ 00827 if( DEBUG_LOC ) 00828 { 00829 fprintf(ioQQQ,"DEBUG iso prad\t%li\t%li\tISO,nelem=\t%li\t%li\tnu, no=\t%li\t%li\t%.4e\t%.4e\t%e\t%e\t%e\n", 00830 iteration, nzone, 00831 ip4,ip3,ip2,ipLinePradMax, 00832 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->TauIn, 00833 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->TauTot, 00834 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->Pesc, 00835 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->Pelec_esc, 00836 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->Pdest); 00837 if( ip2==5 && ipLinePradMax==4 ) 00838 { 00839 double width; 00840 fprintf(ioQQQ,"hit it\n"); 00841 width = RT_LineWidth(&Transitions[ip4][ip3][ip2][ipLinePradMax]); 00842 fprintf(ioQQQ,"DEBUG width %e\n", width); 00843 } 00844 } 00845 } 00846 } 00847 else if( ipLineTypePradMax == 3 ) 00848 { 00849 /* level 1 lines */ 00850 ASSERT( ipLinePradMax>=0 ); 00851 strcpy( pressure.chLineRadPres, "Level1 "); 00852 strcat( pressure.chLineRadPres, chLineLbl(&TauLines[ipLinePradMax]) ); 00853 } 00854 else if( ipLineTypePradMax == 4 ) 00855 { 00856 /* level 2 lines */ 00857 ASSERT( ipLinePradMax>=0 ); 00858 strcpy( pressure.chLineRadPres, "Level2 "); 00859 strcat( pressure.chLineRadPres, chLineLbl(&TauLine2[ipLinePradMax]) ); 00860 } 00861 else if( ipLineTypePradMax == 5 ) 00862 { 00863 /* c12o16 carbon monoxide lines */ 00864 ASSERT( ipLinePradMax>=0 ); 00865 strcpy( pressure.chLineRadPres, "12CO "); 00866 strcat( pressure.chLineRadPres, chLineLbl(&C12O16Rotate[ipLinePradMax]) ); 00867 } 00868 else if( ipLineTypePradMax == 6 ) 00869 { 00870 /* c13o16 carbon monoxide lines */ 00871 ASSERT( ipLinePradMax>=0 ); 00872 strcpy( pressure.chLineRadPres, "13CO "); 00873 strcat( pressure.chLineRadPres, chLineLbl(&C13O16Rotate[ipLinePradMax]) ); 00874 } 00875 else if( ipLineTypePradMax == 7 ) 00876 { 00877 /* FeII lines */ 00878 strcpy( pressure.chLineRadPres, "Fe II "); 00879 } 00880 else if( ipLineTypePradMax == 8 ) 00881 { 00882 /* hyperfine struct lines */ 00883 strcpy( pressure.chLineRadPres, "hyperf "); 00884 ASSERT( ipLinePradMax>=0 ); 00885 strcat( pressure.chLineRadPres, chLineLbl(&HFLines[ipLinePradMax]) ); 00886 } 00887 else if( ipLineTypePradMax == 9 ) 00888 { 00889 /* large H2 lines */ 00890 strcpy( pressure.chLineRadPres, "large H2 "); 00891 } 00892 else 00893 { 00894 fprintf(ioQQQ," PresTotCurrent ipLineTypePradMax set to %li, this is impossible.\n", ipLineTypePradMax); 00895 ShowMe(); 00896 cdEXIT(EXIT_FAILURE); 00897 } 00898 } 00899 00900 if( trace.lgTrace && pressure.pbeta > 0.5 ) 00901 { 00902 fprintf(ioQQQ, 00903 " PresTotCurrent Pbeta:%.2f due to %s\n", 00904 pressure.pbeta , 00905 pressure.chLineRadPres 00906 ); 00907 } 00908 00909 /* >>chng 02 jun 27 - add in magnetic pressure */ 00910 /* start to bring total pressure together */ 00911 TotalPressure_v = pressure.PresGasCurr; 00912 00913 /* these flags are set false by default since constant density is default, 00914 * set true for constant pressure or dynamics */ 00915 TotalPressure_v += pressure.PresRamCurr * pressure.lgPres_ram_ON; 00916 00917 /* magnetic pressure, evaluated in magnetic.c - this can be negative for an ordered field! 00918 * option is on by default, turned off with constant density, or constant gas pressure, cases */ 00921 TotalPressure_v += magnetic.pressure * pressure.lgPres_magnetic_ON; 00922 00923 /* turbulent pressure 00924 * >>chng 06 mar 24, include this by default */ 00925 TotalPressure_v += pressure.PresTurbCurr * DoppVel.lgTurb_pressure; 00926 00927 /* radiation pressure only included in total eqn of state when this flag is 00928 * true, set with constant pressure command */ 00929 /* option to add in internal line radiation pressure */ 00930 TotalPressure_v += pressure.pres_radiation_lines_curr * pressure.lgPres_radiation_ON; 00931 00932 { 00933 enum{DEBUG_LOC=false}; 00934 if( DEBUG_LOC && pressure.PresTotlCurr > SMALLFLOAT /*&& iteration > 1*/ ) 00935 { 00936 fprintf(ioQQQ,"pressureee%li\t%.4e\t%.4e\t%.4e\t%.3f\t%.3f\t%.3f\n", 00937 nzone, 00938 pressure.PresTotlCorrect, 00939 pressure.PresTotlCurr, 00940 TotalPressure_v , 00941 pressure.PresGasCurr/pressure.PresTotlCurr, 00942 pressure.pres_radiation_lines_curr/pressure.PresTotlCurr , 00943 pressure.PresRamCurr/pressure.PresTotlCurr 00944 ); 00945 } 00946 } 00947 00948 if( TotalPressure_v< 0. && magnetic.pressure < 0. ) 00949 { 00950 /* negative pressure due to ordered field overwhelms total pressure - punt */ 00951 fprintf(ioQQQ," The negative pressure due to ordered magnetic field overwhelms the total pressure - please reconsider the geometry & field.\n"); 00952 cdEXIT(EXIT_FAILURE); 00953 } 00954 00955 ASSERT( TotalPressure_v > 0. ); 00956 00957 /* remember highest pressure anywhere */ 00958 pressure.PresMax = MAX2(pressure.PresMax,(realnum)TotalPressure_v); 00959 00960 /* this is what we came for - set the current pressure */ 00961 pressure.PresTotlCurr = TotalPressure_v; 00962 00963 # if 0 00964 /* this is special case where we are working on first zone, at 00965 * illuminated face of the cloud. Here we want to remember the 00966 * initial pressure in case constant pressure is needed */ 00967 /* >>chng 05 jan 10, chng from nzone==1 to nzone<=1 so pressure not changed 00968 * during search phase */ 00969 /*>>chng 06 jun 20, add test on first iteration, or we are holding 00970 * density constant - flag dense.lgDenseInitConstant false if 00971 * constant pressure reset is used - default is true, after first iteration 00972 * initial density is kept constant, when set false with reset option on 00973 * constant pressure density on first iteration is allowed to change to keep 00974 * pressure constant */ 00975 if( nzone <= 1 && (iteration==1 || dense.lgDenseInitConstant) ) 00976 { 00977 double PresTotlInitSave; 00978 double PresRamInitSave; 00979 /* this is first zone, lock onto pressure */ 00980 if( conv.nTotalIoniz ) 00981 { 00982 PresTotlInitSave = pressure.PresTotlInit; 00983 PresRamInitSave = pressure.PresRamInit; 00984 } 00985 else 00986 { 00987 PresTotlInitSave = 0.; 00988 PresRamInitSave = 0.; 00989 } 00990 pressure.PresTotlInit = pressure.PresTotlCurr; 00991 pressure.PresRamInit = pressure.PresRamCurr; 00992 if( trace.lgTrace ) 00993 { 00994 fprintf( ioQQQ, 00995 " PresTotCurrent 1st zn reset PresTotlInit to PresTotlCurr=%.3e " 00996 "PresRamInit to PresRamCurr=%.3e old tot=%.3e old ram %.3e hden=%.3e\n", 00997 pressure.PresTotlInit, 00998 pressure.PresRamInit, 00999 PresTotlInitSave,PresRamInitSave, 01000 dense.gas_phase[ipHYDROGEN] ); 01001 } 01002 } 01003 # endif 01004 return; 01005 }