cloudy trunk
Loading...
Searching...
No Matches
pressure_total.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*PresTotCurrent determine the gas and line radiation pressures for current conditions,
4 * this sets values of pressure.PresTotlCurr, also calls tfidle */
5#include "cddefines.h"
6#include "physconst.h"
7#include "taulines.h"
8#include "opacity.h"
9#include "hextra.h"
10#include "elementnames.h"
11#include "hydrogenic.h"
12#include "conv.h"
13#include "iso.h"
14#include "wind.h"
15#include "magnetic.h"
16#include "doppvel.h"
17#include "rfield.h"
18#include "phycon.h"
19#include "thermal.h"
20#include "hmi.h"
21#include "h2.h"
22#include "dense.h"
23#include "atomfeii.h"
24#include "mole.h"
25#include "dynamics.h"
26#include "trace.h"
27#include "rt.h"
28#include "atmdat.h"
29#include "lines_service.h"
30#include "pressure.h"
31#include "radius.h"
32
33/* this sets values of pressure.PresTotlCurr, also calls tfidle */
35{
36 static long int
37 /* used in debug print statement to see which line adds most pressure */
38 ipLineTypePradMax=-1 ,
39 /* points to line causing radiation pressure */
40 ipLinePradMax=-LONG_MAX,
41 ip2=-LONG_MAX,
42 ip3=-LONG_MAX,
43 ip4=-LONG_MAX;
44
45 /* the line radiation pressure variables that must be preserved since
46 * a particular line radiation pressure may not be evaluated if it is
47 * not important */
48 static double Piso_seq[NISO]={0.,0.},
49 PLevel1=0.,
50 PLevel2=0.,
51 PHfs=0.,
52 PCO=0.,
53 P_H2=0.,
54 P_FeII=0.,
55 P_dBase=0.;
56
57 double
58 RadPres1,
59 TotalPressure_v,
60 pmx;
61
62 DEBUG_ENTRY( "PresTotCurrent()" );
63
64 if( !conv.nTotalIoniz )
65 {
66 /* resetting ipLinePradMax, necessary for
67 * multiple cloudy calls in single coreload. */
68 ipLinePradMax=-LONG_MAX;
69 //pressure.PresTotlCurr = 0.;
70 }
71
72 /* PresTotCurrent - evaluate total pressure, dyne cm^-2
73 * and radiative acceleration */
74
75 /* several loops over the emission lines for radiation pressure and
76 * radiative acceleration are expensive due to the number of lines and
77 * the memory they occupy. Many equations of state do not include
78 * radiation pressure or radiative acceleration. Only update these
79 * when included in EOS - when not included only evaluate once per zone.
80 * do it once per zone since we will still report these quantities.
81 * this flag says whether we must update all terms */
82 bool lgMustEvaluate = false;
83
84 /* this says we already have an evaluation in this zone so do not
85 * evaluate terms known to be trivial, even when reevaluating major terms */
86 bool lgMustEvaluateTrivial = false;
87 /* if an individual agent is larger than this fraction of the total current
88 * radiation pressure then it is not trivial
89 * conv.PressureErrorAllowed is relative error allowed in pressure */
90 double TrivialLineRadiationPressure = conv.PressureErrorAllowed/10. *
91 pressure.pres_radiation_lines_curr;
92
93 /* reevaluate during search phase when conditions are changing dramatically */
94 if( conv.lgSearch )
95 {
96 lgMustEvaluate = true;
97 lgMustEvaluateTrivial = true;
98 }
99
100 /* reevaluate if zone or iteration has changed */
101 static long int nzoneEvaluated=0, iterationEvaluated=0;
102 if( nzone!=nzoneEvaluated || iteration!=iterationEvaluated )
103 {
104 lgMustEvaluate = true;
105 lgMustEvaluateTrivial = true;
106 /* this flag says which source of radiation pressure was strongest */
107 ipLineTypePradMax = -1;
108 }
109
110 /* constant pressure or dynamical sim - we must reevaluate terms
111 * but do not need to reevaluate trivial contributors */
112 if( (strcmp(dense.chDenseLaw,"WIND") == 0 ) ||
113 (strcmp(dense.chDenseLaw,"DYNA") == 0 ) ||
114 (strcmp(dense.chDenseLaw,"CPRE") == 0 ) )
115 lgMustEvaluate = true;
116
117 if( lgMustEvaluate )
118 {
119 /* we are reevaluating quantities in this zone and iteration,
120 * remember that we did it */
121 nzoneEvaluated = nzone;
122 iterationEvaluated = iteration;
123 }
124
125 /* update density - temperature variables which may have changed */
126 /* must call TempChange since ionization has changed, there are some
127 * terms that affect collision rates (H0 term in electron collision) */
128 TempChange(phycon.te , false);
129
131
132 /* evaluate the total ionization energy on a scale where neutral atoms
133 * correspond to an energy of zero, so the ions have a positive energy */
134 phycon.EnergyIonization = 0.;
135#if 0
136 for( long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
137 {
138 for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
139 {
140 /* lgMETALS is true, set false with "no advection metals" command */
141 int kadvec = dynamics.lgMETALS;
142 /* this gives the iso sequence for this ion - should it be included in the
143 * advected energy equation? lgISO is true, set false with
144 * "no advection H-like" and "he-like" - for testing*/
145 ipISO = nelem-ion;
146 fixit(); /* should this be kadvec = kadvec && dynamics.lgISO[ipISO]; ? */
147 if( ipISO >= 0 && ipISO<NISO )
148 kadvec = dynamics.lgISO[ipISO];
149 for( long i=1; i<=ion; ++i )
150 {
151 long int nelec = nelem+1 - i;
152 /* this is the sum of all the energy needed to bring the atom up
153 * to the ion+1 level of ionization - at this point a positive number */
154 phycon.EnergyIonization += dense.xIonDense[nelem][ion] *
155 t_ADfA::Inst().ph1(Heavy.nsShells[nelem][i-1]-1,nelec,nelem,0)/EVRYD* 0.9998787*kadvec;
156 }
157 }
158 }
159 /* convert to ergs/cm^3 */
160 phycon.EnergyIonization *= EN1RYD;
161#endif
162
166 phycon.EnergyBinding = 0.;
167
168 /* find total number of atoms and ions */
169 SumDensities();
170
171 /* the current gas pressure */
172 pressure.PresGasCurr = dense.pden*phycon.te*BOLTZMANN;
173 /*fprintf(ioQQQ,"DEBUG gassss %.2f %.4e %.4e %.4e \n",
174 fnzone, pressure.PresGasCurr , dense.pden , pressure.PresInteg );*/
175
176 /* >>chng 01 nov 02, move to here from dynamics routines */
177 /* >>chng 02 jun 18 (rjrw), fix to use local values */
178 /* WJH 21 may 04, now recalculate wind v for the first zone too.
179 * This is necessary when we are forcing the sub or supersonic branch */
180 if(!(wind.lgBallistic() || wind.lgStatic()))
181 {
182 wind.windv = DynaFlux(radius.depth)/(dense.xMassDensity);
183 }
184
185 /* this is the current ram pressure, at this location */
186 pressure.PresRamCurr = dense.xMassDensity*POW2( wind.windv );
187
188 /* this is the current turbulent pressure, not now added to equation of state
189 * >>chng 06 mar 29, add Heiles_Troland_F / 6. term*/
190 pressure.PresTurbCurr = dense.xMassDensity*POW2( DoppVel.TurbVel ) *
191 DoppVel.Heiles_Troland_F / 6.;
192
194
195 /* radiative acceleration, lines and continuum */
196 if( lgMustEvaluate )
197 {
198 /* continuous radiative acceleration */
199 double rforce = 0.;
200 double relec = 0.;
201 for( long i=0; i < rfield.nflux; i++ )
202 {
203 rforce += (rfield.flux[0][i] + rfield.outlin[0][i] + rfield.outlin_noplot[i]+ rfield.ConInterOut[i])*
204 rfield.anu[i]*(opac.opacity_abs[i] + opac.opacity_sct[i]);
205
206 /* electron scattering acceleration - used only to derive force multiplier */
207 relec +=
208 (rfield.flux[0][i] + rfield.outlin[0][i] + rfield.outlin_noplot[i]+
209 rfield.ConInterOut[i]) *
210 opac.OpacStack[i-1+opac.iopcom]*dense.eden*rfield.anu[i];
211 }
212 /* total continuum radiative acceleration */
213 wind.AccelCont = (realnum)(rforce*EN1RYD/SPEEDLIGHT/dense.xMassDensity);
214
215 wind.AccelElectron = (realnum)(relec*EN1RYD/SPEEDLIGHT/dense.xMassDensity);
216
217 /* line acceleration; xMassDensity is gm per cc */
218 wind.AccelLine = (realnum)(RT_line_driving()/SPEEDLIGHT/dense.xMassDensity);
219
220 /* total acceleration cm s^-2 */
221 wind.AccelTotalOutward = wind.AccelCont + wind.AccelLine;
222
223 /* find wind.fmul - the force multiplier; wind.AccelElectron can be zero for very low
224 * densities - fmul is of order unity - wind.AccelLine and wind.AccelCont
225 * are both floats to will underflow long before wind.AccelElectron will - fmul is only used
226 * in output, not in any physics */
227 if( wind.AccelElectron > SMALLFLOAT )
228 wind.fmul = (realnum)( (wind.AccelLine + wind.AccelCont) / wind.AccelElectron);
229 else
230 wind.fmul = 0.;
231
232 double reff = radius.Radius-radius.drad/2.;
233 /* inward acceleration of gravity cm s^-2 */
234 wind.AccelGravity = (realnum)(
235 /* wind.comass - mass of star in solar masses */
236 GRAV_CONST*wind.comass*SOLAR_MASS/POW2(reff)*
237 /* wind.DiskRadius normally zero, set with disk option on wind command */
238 (1.-wind.DiskRadius/reff) );
239# if 0
240 if( fudge(-1) )
241 fprintf(ioQQQ,"DEBUG pressure_total updates AccelTotalOutward to %.2e grav %.2e\n",
242 wind.AccelTotalOutward , wind.AccelGravity );
243# endif
244 }
245
246 /* must always evaluate H La width since used elsewhere */
247 if( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().PopOpc() > 0. )
248 {
250 }
251 else
252 hydro.HLineWidth = 0.;
253
254
255 /* find max rad pressure even if capped
256 * lgLineRadPresOn is turned off by NO RADIATION PRESSURE command */
257 if( trace.lgTrace )
258 {
259 fprintf(ioQQQ,
260 " PresTotCurrent nzone %li iteration %li lgMustEvaluate:%c "
261 "lgMustEvaluateTrivial:%c "
262 "pressure.lgLineRadPresOn:%c "
263 "rfield.lgDoLineTrans:%c \n",
264 nzone , iteration , TorF(lgMustEvaluate) , TorF(lgMustEvaluateTrivial),
265 TorF(pressure.lgLineRadPresOn), TorF(rfield.lgDoLineTrans) );
266 }
267
268 if( lgMustEvaluate && pressure.lgLineRadPresOn && rfield.lgDoLineTrans )
269 {
270 /* RadPres is pressure due to lines, lgPres_radiation_ON turns off or on */
271 pressure.pres_radiation_lines_curr = 0.;
272 /* used to remember largest radiation pressure source */
273 pmx = 0.;
274 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
275 {
276 if( lgMustEvaluateTrivial || Piso_seq[ipISO] > TrivialLineRadiationPressure )
277 {
278 Piso_seq[ipISO] = 0.;
279 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
280 {
281 /* does this ion stage exist? */
282 if( dense.IonHigh[nelem] >= nelem + 1 - ipISO )
283 {
284 /* do not include highest levels since maser can occur due to topoff,
285 * and pops are set to small number in this case */
286 for( long ipHi=1; ipHi <iso_sp[ipISO][nelem].numLevels_local - iso_sp[ipISO][nelem].nCollapsed_local; ipHi++ )
287 {
288 for( long ipLo=0; ipLo < ipHi; ipLo++ )
289 {
290 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() <= 0 )
291 continue;
292
293 ASSERT( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() > iso_ctrl.SmallA );
294
295 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc() > SMALLFLOAT &&
296 /* test that have not overrun optical depth scale */
297 ( (iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauTot() -
298 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn()) > SMALLFLOAT ) &&
299 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pesc() > FLT_EPSILON*100. )
300 {
301 RadPres1 = PressureRadiationLine( iso_sp[ipISO][nelem].trans(ipHi,ipLo), GetDopplerWidth(dense.AtomicWeight[nelem]) );
302
303 if( RadPres1 > pmx )
304 {
305 ipLineTypePradMax = 2;
306 pmx = RadPres1;
307 ip4 = ipISO;
308 ip3 = nelem;
309 ip2 = ipHi;
310 ipLinePradMax = ipLo;
311 }
312 Piso_seq[ipISO] += RadPres1;
313 {
314 /* option to print particulars of some line when called */
315 enum {DEBUG_LOC=false};
316 if( DEBUG_LOC && ipISO==ipH_LIKE && ipLo==3 && ipHi==5 && nzone > 260 )
317 {
318 fprintf(ioQQQ,
319 "DEBUG prad1 \t%.2f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t\n",
320 fnzone,
321 RadPres1,
322 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc(),
323 iso_sp[ipISO][nelem].st[ipLo].Pop(),
324 iso_sp[ipISO][nelem].st[ipHi].Pop(),
325 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pesc());
326 }
327 }
328 }
329 }
330 }
331 }
332 }
333 ASSERT( Piso_seq[ipISO] >= 0. );
334 }
335 pressure.pres_radiation_lines_curr += Piso_seq[ipISO];
336 }
337 {
338 /* option to print particulars of some line when called */
339 enum {DEBUG_LOC=false};
340# if 0
341 if( DEBUG_LOC /*&& iteration > 1*/ && nzone > 260 )
342 {
343 fprintf(ioQQQ,
344 "DEBUG prad2 \t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
345 nzone,
346 pmx,
347 iso_sp[ipISO][ip3].trans(ipLinePradMax,ip2).Emis().PopOpc(),
348 iso_sp[ipISO][ip3].st[0].Pop(),
349 iso_sp[ipISO][ip3].st[2].Pop(),
350 iso_sp[ipISO][ip3].st[3].Pop(),
351 iso_sp[ipISO][ip3].st[4].Pop(),
352 iso_sp[ipISO][ip3].st[5].Pop(),
353 iso_sp[ipISO][ip3].st[6].Pop());
354 }
355# endif
356 if( DEBUG_LOC /*&& iteration > 1 && nzone > 150 */)
357 {
358 fprintf(ioQQQ,
359 "DEBUG prad3\t%.2e\t%li\t%li\t%li\t%li\t%.2e\t%.2e\t%.2e\n",
360 pmx,
361 ip4,
362 ip3,
363 ip2,
364 ipLinePradMax,
365 iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().PopOpc(),
366 iso_sp[ip4][ip3].st[ip2].Pop(),
367 1.-iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().Pesc() );
368 }
369 }
370
371 if( lgMustEvaluateTrivial || PLevel1 > TrivialLineRadiationPressure )
372 {
373 PLevel1 = 0.;
374 /* line radiation pressure from large set of level 1 lines */
375 for( long i=1; i <= nLevel1; i++ )
376 {
377 if( (*TauLines[i].Hi()).Pop() > 1e-30 )
378 {
379 RadPres1 = PressureRadiationLine( TauLines[i], GetDopplerWidth(dense.AtomicWeight[(*TauLines[i].Hi()).nelem()-1]) );
380
381 if( RadPres1 > pmx )
382 {
383 ipLineTypePradMax = 3;
384 pmx = RadPres1;
385 ipLinePradMax = i;
386 }
387 PLevel1 += RadPres1;
388 }
389 }
390 ASSERT( PLevel1 >= 0. );
391 }
392 pressure.pres_radiation_lines_curr += PLevel1;
393
394 if( lgMustEvaluateTrivial || PLevel2 > TrivialLineRadiationPressure )
395 {
396 /* level 2 lines */
397 PLevel2 = 0.;
398 for( long i=0; i < nWindLine; i++ )
399 {
400 if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
401 {
402 if( (*TauLine2[i].Hi()).Pop() > 1e-30 )
403 {
404 RadPres1 = PressureRadiationLine( TauLine2[i], GetDopplerWidth(dense.AtomicWeight[(*TauLine2[i].Hi()).nelem()-1]) );
405
406 PLevel2 += RadPres1;
407 if( RadPres1 > pmx )
408 {
409 ipLineTypePradMax = 4;
410 pmx = RadPres1;
411 ipLinePradMax = i;
412 }
413 }
414 }
415 }
416 ASSERT( PLevel2 >= 0. );
417 }
418 pressure.pres_radiation_lines_curr += PLevel2;
419
420 /* fine structure lines */
421 if( lgMustEvaluateTrivial || PHfs > TrivialLineRadiationPressure )
422 {
423 PHfs = 0.;
424 for( long i=0; i < nHFLines; i++ )
425 {
426 if( (*HFLines[i].Hi()).Pop() > 1e-30 )
427 {
428 RadPres1 = PressureRadiationLine( HFLines[i], GetDopplerWidth(dense.AtomicWeight[(*HFLines[i].Hi()).nelem()-1]) );
429
430 PHfs += RadPres1;
431 if( RadPres1 > pmx )
432 {
433 ipLineTypePradMax = 8;
434 pmx = RadPres1;
435 ipLinePradMax = i;
436 }
437 }
438 }
439 ASSERT( PHfs >= 0. );
440 }
441 pressure.pres_radiation_lines_curr += PHfs;
442
443 /* radiation pressure due to H2 lines */
444 if( lgMustEvaluateTrivial || P_H2 > TrivialLineRadiationPressure )
445 {
446 P_H2 = 0.;
447 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
448 {
449 P_H2 += (*diatom)->H2_RadPress();
450 /* flag to remember H2 radiation pressure */
451 if( P_H2 > pmx )
452 {
453 pmx = P_H2;
454 ipLineTypePradMax = 9;
455 }
456 ASSERT( P_H2 >= 0. );
457 }
458 }
459 pressure.pres_radiation_lines_curr += P_H2;
460
461 /* radiation pressure due to FeII lines and large atom */
462 if( lgMustEvaluateTrivial || P_FeII > TrivialLineRadiationPressure )
463 {
464 P_FeII = FeIIRadPress();
465 if( P_FeII > pmx )
466 {
467 pmx = P_FeII;
468 ipLineTypePradMax = 7;
469 }
470 ASSERT( P_FeII >= 0. );
471 }
472 pressure.pres_radiation_lines_curr += P_FeII;
473
474 /* do lines from third-party databases */
475 if( lgMustEvaluateTrivial || P_dBase > TrivialLineRadiationPressure )
476 {
477 P_dBase = 0.;
478 for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
479 {
480 if( dBaseSpecies[ipSpecies].lgActive )
481 {
482 realnum DopplerWidth = GetDopplerWidth( dBaseSpecies[ipSpecies].fmolweight );
483 for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin();
484 tr != dBaseTrans[ipSpecies].end(); ++tr)
485 {
486 int ipHi = (*tr).ipHi();
487 if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local)
488 continue;
489 int ipLo = (*tr).ipLo();
490 if( (*tr).ipCont() > 0 && (*(*tr).Hi()).Pop() > 1e-30 )
491 {
492 RadPres1 = PressureRadiationLine( *tr, DopplerWidth );
493
494 if( RadPres1 > pmx )
495 {
496 ipLineTypePradMax = 10;
497 pmx = RadPres1;
498 ip3 = ipSpecies;
499 ip2 = ipHi;
500 ipLinePradMax = ipLo;
501 }
502 P_dBase += RadPres1;
503 }
504 }
505 }
506 }
507 ASSERT( P_dBase >= 0. );
508 }
509 pressure.pres_radiation_lines_curr += P_dBase;
510
511 }
512 else if( !pressure.lgLineRadPresOn || !rfield.lgDoLineTrans )
513 {
514 /* case where radiation pressure is not turned on */
515 ipLinePradMax = -1;
516 ipLineTypePradMax = 0;
517 }
518
519 fixit(); // all other line stacks need to be included here.
520 // can we just sweep over line stack? Is that ready yet?
521
522 /* the ratio of radiation to total (gas + continuum + etc) pressure */
523 pressure.pbeta = (realnum)(pressure.pres_radiation_lines_curr/SDIV(pressure.PresTotlCurr));
524
525 /* this would be a major logical error */
526 if( pressure.pres_radiation_lines_curr < 0. )
527 {
528 fprintf(ioQQQ,
529 "PresTotCurrent: negative pressure, constituents follow %e %e %e %e %e \n",
530 Piso_seq[ipH_LIKE],
531 Piso_seq[ipHE_LIKE],
532 PLevel1,
533 PLevel2,
534 PCO);
535 ShowMe();
537 }
538
539 /* following test will never succeed, here to trick lint, ipLineTypePradMax is only used
540 * when needed for debug */
541 if( trace.lgTrace && ipLineTypePradMax <0 )
542 {
543 fprintf(ioQQQ,
544 " PresTotCurrent, pressure.pbeta = %e, ipLineTypePradMax%li ipLinePradMax=%li \n",
545 pressure.pbeta,ipLineTypePradMax,ipLinePradMax );
546 }
547
548 /* this is the total internal energy of the gas, the amount of energy needed
549 * to produce the current level populations, relative to ground,
550 * only used for energy terms in advection */
551 phycon.EnergyExcitation = 0.;
552#if 0
553 fixit(); /* the quantities phycon.EnergyExcitation, phycon.EnergyIonization, phycon.EnergyBinding
554 * are not used anywhere, except in print statements... */
555 broken(); /* the code below contains serious bugs! It is supposed to implement loops
556 * over quantum states in order to evaluate the potential energy stored in
557 * excited states of all atoms, ions, and molecules. The code below however
558 * often implements loops over all combinations of upper and lower levels!
559 * This code needs to be rewritten once quantumstates are fully implemented. */
560 ipLo = 0;
561 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
562 {
563 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
564 {
565 if( dense.IonHigh[nelem] == nelem + 1 - ipISO )
566 {
567 /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
568 for( long ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
569 {
570 phycon.EnergyExcitation +=
571 iso_sp[ipISO][nelem].st[ipHi].Pop() *
572 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyErg() *
573 /* last term is option to turn off energy term, no advection hlike, he-like */
574 dynamics.lgISO[ipISO];
575 }
576 }
577 }
578 }
579
580 if( dynamics.lgISO[ipH_LIKE] )
581 /* internal energy of H2 */
582 phycon.EnergyExcitation += H2_InterEnergy();
583
584 /* this is option to turn off energy effects of advection, no advection metals */
585 if( dynamics.lgMETALS )
586 {
587 for( long i=1; i <= nLevel1; i++ )
588 {
589 phycon.EnergyExcitation +=
590 (*TauLines[i].Hi()).Pop()* TauLines[i].EnergyErg;
591 }
592 for( long i=0; i < nWindLine; i++ )
593 {
594 if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
595 {
596 phycon.EnergyExcitation +=
597 (*TauLine2[i].Hi()).Pop()* TauLine2[i].EnergyErg;
598 }
599 }
600 for( long i=0; i < nHFLines; i++ )
601 {
602 phycon.EnergyExcitation +=
603 (*HFLines[i].Hi()).Pop()* HFLines[i].EnergyErg;
604 }
605
606 /* internal energy of large FeII atom */
607 phycon.EnergyExcitation += FeII_InterEnergy();
608 }
609#endif
610
611 /* ==================================================================
612 * end internal energy of atoms and molecules */
613
614 /* evaluate some parameters to do with magnetic field */
616
617 /*lint -e644 Piso_seq not init */
618 if( trace.lgTrace && (pressure.PresTotlCurr > SMALLFLOAT) )
619 {
620 fprintf(ioQQQ,
621 " PresTotCurrent zn %.2f Ptot:%.5e Pgas:%.3e Prad:%.3e Pmag:%.3e Pram:%.3e "
622 "gas parts; H:%.2e He:%.2e L1:%.2e L2:%.2e CO:%.2e fs%.2e h2%.2e turb%.2e\n",
623 fnzone,
624 pressure.PresTotlCurr,
625 pressure.PresGasCurr/pressure.PresTotlCurr,
626 pressure.pres_radiation_lines_curr*pressure.lgPres_radiation_ON/pressure.PresTotlCurr,
627 magnetic.pressure/pressure.PresTotlCurr,
628 pressure.PresRamCurr*pressure.lgPres_ram_ON/pressure.PresTotlCurr,
629 Piso_seq[ipH_LIKE]/pressure.PresTotlCurr,
630 Piso_seq[ipHE_LIKE]/pressure.PresTotlCurr,
631 PLevel1/pressure.PresTotlCurr,
632 PLevel2/pressure.PresTotlCurr,
633 PCO/pressure.PresTotlCurr,
634 PHfs/pressure.PresTotlCurr,
635 P_H2/pressure.PresTotlCurr,
636 pressure.PresTurbCurr*DoppVel.lgTurb_pressure/pressure.PresTotlCurr);
637 /*lint +e644 Piso_seq not initialized */
638 }
639
640 /* Conserved quantity in steady-state energy equation for dynamics:
641 * thermal energy, since recombination is treated as cooling
642 * (i.e. is loss of electron k.e. to emitted photon), so don't
643 * include
644 * ...phycon.EnergyExcitation + phycon.EnergyIonization + phycon.EnergyBinding
645 * */
646
647 /* constant density case is also hypersonic case */
648 if( dynamics.lgTimeDependentStatic )
649 {
650 /* this branch is time dependent single-zone */
651 /* \todo 1 this has 3/2 on the PresGasCurr while true dynamics case below
652 * has 5/2 - so this is not really the enthalpy density - better
653 * would be to always use this term and add the extra PresGasCurr
654 * when the enthalpy is actually needed */
655 phycon.EnthalpyDensity =
656 0.5*POW2(wind.windv)*dense.xMassDensity + /* KE */
657 3./2.*pressure.PresGasCurr + /* thermal plus PdV work */
658 magnetic.EnthalpyDensity; /* magnetic terms */
659 //pressure.RhoGravity * (radius.Radius/(radius.Radius+radius.drad)); /* gravity */
660 }
661 else
662 {
663 /* this branch either advective or constant pressure */
664 /*fprintf(ioQQQ,"DEBUG enthalpy HIT2\n");*/
665 /* this is usual dynamics case, or time-varying case where pressure
666 * is kept constant and PdV work is added to the cell */
667 phycon.EnthalpyDensity =
668 0.5*POW2(wind.windv)*dense.xMassDensity + /* KE */
669 5./2.*pressure.PresGasCurr + /* thermal plus PdV work */
670 magnetic.EnthalpyDensity; /* magnetic terms */
671 //pressure.RhoGravity * (radius.Radius/(radius.Radius+radius.drad)); /* gravity */
672 }
673
674 /* stop model from running away on first iteration, when line optical
675 * depths are not defined correctly anyway.
676 * if( iter.le.itermx .and. RadPres.ge.GasPres ) then
677 * >>chng 97 jul 23, only cap radiation pressure on first iteration */
679 {
680 /* stop RadPres from exceeding the gas pressure on first iteration */
681 pressure.pres_radiation_lines_curr =
682 MIN2(pressure.pres_radiation_lines_curr,pressure.PresGasCurr);
683 pressure.lgPradCap = true;
684 }
685
686 /* remember the globally most important line, in the entire model
687 * test on nzone so we only do test after solution is stable */
688 if( pressure.pbeta > pressure.RadBetaMax && nzone )
689 {
690 pressure.RadBetaMax = pressure.pbeta;
691 pressure.ipPradMax_line = ipLinePradMax;
692 pressure.ipPradMax_nzone = nzone;
693 if( ipLineTypePradMax == 2 )
694 {
695 /* hydrogenic */
696 strcpy( pressure.chLineRadPres, "ISO ");
697 ASSERT( ip4 < NISO && ip3<LIMELM );
698 ASSERT( ipLinePradMax>=0 && ip2>=0 && ip3>=0 && ip4>=0 );
699 strcat( pressure.chLineRadPres, chLineLbl(iso_sp[ip4][ip3].trans(ip2,ipLinePradMax)) );
700 {
701 /* option to print particulars of some line when called */
702 enum {DEBUG_LOC=false};
703 /*lint -e644 Piso_seq not initialized */
704 /* trace serious constituents in radiation pressure */
705 if( DEBUG_LOC )
706 {
707 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",
709 ip4,ip3,ip2,ipLinePradMax,
710 iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().TauIn(),
711 iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().TauTot(),
712 iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().Pesc(),
713 iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().Pelec_esc(),
714 iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().Pdest());
715 if( ip2==5 && ipLinePradMax==4 )
716 {
717 double width;
718 fprintf(ioQQQ,"hit it\n");
719 width = RT_LineWidth(iso_sp[ip4][ip3].trans(ip2,ipLinePradMax),GetDopplerWidth(dense.AtomicWeight[ip3]));
720 fprintf(ioQQQ,"DEBUG width %e\n", width);
721 }
722 }
723 }
724 }
725 else if( ipLineTypePradMax == 3 )
726 {
727 /* level 1 lines */
728 ASSERT( ipLinePradMax>=0 );
729 strcpy( pressure.chLineRadPres, "Level1 ");
730 strcat( pressure.chLineRadPres, chLineLbl(TauLines[ipLinePradMax]) );
731 }
732 else if( ipLineTypePradMax == 4 )
733 {
734 /* level 2 lines */
735 ASSERT( ipLinePradMax>=0 );
736 strcpy( pressure.chLineRadPres, "Level2 ");
737 strcat( pressure.chLineRadPres, chLineLbl(TauLine2[ipLinePradMax]) );
738 }
739 else if( ipLineTypePradMax == 5 )
740 {
742 }
743 else if( ipLineTypePradMax == 6 )
744 {
746 }
747 else if( ipLineTypePradMax == 7 )
748 {
749 /* FeII lines */
750 strcpy( pressure.chLineRadPres, "Fe II ");
751 }
752 else if( ipLineTypePradMax == 8 )
753 {
754 /* hyperfine struct lines */
755 strcpy( pressure.chLineRadPres, "hyperf ");
756 ASSERT( ipLinePradMax>=0 );
757 strcat( pressure.chLineRadPres, chLineLbl(HFLines[ipLinePradMax]) );
758 }
759 else if( ipLineTypePradMax == 9 )
760 {
761 /* large H2 lines */
762 strcpy( pressure.chLineRadPres, "large H2 ");
763 }
764 else if( ipLineTypePradMax == 10 )
765 {
766 /* database lines */
767 strcpy( pressure.chLineRadPres, "dBaseLin " );
768 strcat( pressure.chLineRadPres, dBaseSpecies[ip3].chLabel );
769 }
770 else
771 {
772 fprintf(ioQQQ," PresTotCurrent ipLineTypePradMax set to %li, this is impossible.\n", ipLineTypePradMax);
773 ShowMe();
775 }
776 }
777
778 if( trace.lgTrace && pressure.pbeta > 0.5 )
779 {
780 fprintf(ioQQQ,
781 " PresTotCurrent Pbeta:%.2f due to %s\n",
782 pressure.pbeta ,
783 pressure.chLineRadPres
784 );
785 }
786
787 /* >>chng 02 jun 27 - add in magnetic pressure */
788 /* start to bring total pressure together */
789 TotalPressure_v = pressure.PresGasCurr;
790
791 /* these flags are set false by default since constant density is default,
792 * set true for constant pressure or dynamics */
793 TotalPressure_v += pressure.PresRamCurr * pressure.lgPres_ram_ON;
794
795 /* magnetic pressure, evaluated in magnetic.c - this can be negative for an ordered field!
796 * option is on by default, turned off with constant density, or constant gas pressure, cases */
799 TotalPressure_v += magnetic.pressure * pressure.lgPres_magnetic_ON;
800
801 /* turbulent pressure
802 * >>chng 06 mar 24, include this by default */
803 TotalPressure_v += pressure.PresTurbCurr * DoppVel.lgTurb_pressure;
804
805 /* radiation pressure only included in total eqn of state when this flag is
806 * true, set with constant pressure command */
807 /* option to add in internal line radiation pressure */
808 TotalPressure_v += pressure.pres_radiation_lines_curr * pressure.lgPres_radiation_ON;
809
810 {
811 enum{DEBUG_LOC=false};
812 if( DEBUG_LOC && pressure.PresTotlCurr > SMALLFLOAT /*&& iteration > 1*/ )
813 {
814 fprintf(ioQQQ,"pressureee%li\t%.4e\t%.4e\t%.4e\t%.3f\t%.3f\t%.3f\n",
815 nzone,
816 pressure.PresTotlError,
817 pressure.PresTotlCurr,
818 TotalPressure_v ,
819 pressure.PresGasCurr/pressure.PresTotlCurr,
820 pressure.pres_radiation_lines_curr/pressure.PresTotlCurr ,
821 pressure.PresRamCurr/pressure.PresTotlCurr
822 );
823 }
824 }
825
826 if( TotalPressure_v< 0. )
827 {
828 ASSERT( magnetic.pressure < 0. );
829
830 /* negative pressure due to ordered field overwhelms total pressure - punt */
831 fprintf(ioQQQ," The negative pressure due to ordered magnetic field overwhelms the total outward pressure - please reconsider the geometry & field.\n");
833 }
834
835 ASSERT( TotalPressure_v > 0. );
836
837 /* remember highest pressure anywhere */
838 pressure.PresMax = MAX2(pressure.PresMax,(realnum)TotalPressure_v);
839
840 /* this is what we came for - set the current pressure */
841 pressure.PresTotlCurr = TotalPressure_v;
842
843 return;
844}
double FeIIRadPress(void)
double FeII_InterEnergy(void)
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
long int iteration
Definition cddefines.cpp:16
double fnzone
Definition cddefines.cpp:15
#define ASSERT(exp)
Definition cddefines.h:578
#define MIN2
Definition cddefines.h:761
void broken(void)
Definition service.cpp:982
double fudge(long int ipnt)
Definition service.cpp:481
const int LIMELM
Definition cddefines.h:258
#define POW2
Definition cddefines.h:929
#define EXIT_FAILURE
Definition cddefines.h:140
char TorF(bool l)
Definition cddefines.h:710
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void ShowMe(void)
Definition service.cpp:181
void fixit(void)
Definition service.cpp:991
long nWindLine
Definition cdinit.cpp:19
static t_ADfA & Inst()
Definition cddefines.h:175
TransitionProxy::iterator iterator
Definition transition.h:280
realnum ph1(int i, int j, int k, int l) const
Definition atmdat.h:329
t_conv conv
Definition conv.cpp:5
const realnum SMALLFLOAT
Definition cpu.h:191
bool lgElemsConserved(void)
Definition dense.cpp:99
t_dense dense
Definition dense.cpp:24
void SumDensities(void)
Definition dense.cpp:200
t_DoppVel DoppVel
Definition doppvel.cpp:5
realnum GetDopplerWidth(realnum massAMU)
t_dynamics dynamics
Definition dynamics.cpp:44
realnum DynaFlux(double depth)
vector< diatomics * > diatoms
Definition h2.cpp:8
vector< diatomics * >::iterator diatom_iter
Definition h2.h:13
t_Heavy Heavy
Definition heavy.cpp:5
t_hydro hydro
Definition hydrogenic.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
const int ipH1s
Definition iso.h:27
const int ipHE_LIKE
Definition iso.h:63
const int ipH2p
Definition iso.h:29
const int ipH_LIKE
Definition iso.h:62
t_magnetic magnetic
Definition magnetic.cpp:17
void Magnetic_evaluate(void)
Definition magnetic.cpp:39
t_opac opac
Definition opacity.cpp:5
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double SOLAR_MASS
Definition physconst.h:71
UNUSED const double SPEEDLIGHT
Definition physconst.h:100
UNUSED const double BOLTZMANN
Definition physconst.h:97
UNUSED const double EN1RYD
Definition physconst.h:179
UNUSED const double EVRYD
Definition physconst.h:189
UNUSED const double GRAV_CONST
Definition physconst.h:109
t_pressure pressure
Definition pressure.cpp:5
double PressureRadiationLine(const TransitionProxy &t, realnum DopplerWidth)
Definition pressure.h:18
void PresTotCurrent(void)
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
double RT_line_driving(void)
double RT_LineWidth(const TransitionProxy &t, realnum DopplerWidth)
long int nSpecies
Definition taulines.cpp:21
TransitionList TauLine2("TauLine2", &AnonStates)
vector< TransitionList > dBaseTrans
Definition taulines.cpp:17
TransitionList HFLines("HFLines", &AnonStates)
long int nLevel1
Definition taulines.cpp:28
long int nHFLines
Definition taulines.cpp:31
TransitionList TauLines("TauLines", &AnonStates)
species * dBaseSpecies
Definition taulines.cpp:14
void TempChange(double TempNew, bool lgForceUpdate)
t_trace trace
Definition trace.cpp:5
char * chLineLbl(const TransitionProxy &t)
Wind wind
Definition wind.cpp:5