38 ipLineTypePradMax=-1 ,
40 ipLinePradMax=-LONG_MAX,
48 static double Piso_seq[
NISO]={0.,0.},
64 if( !
conv.nTotalIoniz )
68 ipLinePradMax=-LONG_MAX;
82 bool lgMustEvaluate =
false;
86 bool lgMustEvaluateTrivial =
false;
90 double TrivialLineRadiationPressure =
conv.PressureErrorAllowed/10. *
96 lgMustEvaluate =
true;
97 lgMustEvaluateTrivial =
true;
101 static long int nzoneEvaluated=0, iterationEvaluated=0;
104 lgMustEvaluate =
true;
105 lgMustEvaluateTrivial =
true;
107 ipLineTypePradMax = -1;
112 if( (strcmp(
dense.chDenseLaw,
"WIND") == 0 ) ||
113 (strcmp(
dense.chDenseLaw,
"DYNA") == 0 ) ||
114 (strcmp(
dense.chDenseLaw,
"CPRE") == 0 ) )
115 lgMustEvaluate =
true;
121 nzoneEvaluated =
nzone;
134 phycon.EnergyIonization = 0.;
138 for(
long ion=
dense.IonLow[nelem]; ion<=
dense.IonHigh[nelem]; ++ion )
147 if( ipISO >= 0 && ipISO<
NISO )
149 for(
long i=1; i<=ion; ++i )
151 long int nelec = nelem+1 - i;
154 phycon.EnergyIonization +=
dense.xIonDense[nelem][ion] *
166 phycon.EnergyBinding = 0.;
180 if(!(
wind.lgBallistic() ||
wind.lgStatic()))
201 for(
long i=0; i <
rfield.nflux; i++ )
221 wind.AccelTotalOutward =
wind.AccelCont +
wind.AccelLine;
238 (1.-
wind.DiskRadius/reff) );
241 fprintf(
ioQQQ,
"DEBUG pressure_total updates AccelTotalOutward to %.2e grav %.2e\n",
242 wind.AccelTotalOutward ,
wind.AccelGravity );
252 hydro.HLineWidth = 0.;
260 " PresTotCurrent nzone %li iteration %li lgMustEvaluate:%c "
261 "lgMustEvaluateTrivial:%c "
262 "pressure.lgLineRadPresOn:%c "
263 "rfield.lgDoLineTrans:%c \n",
268 if( lgMustEvaluate &&
pressure.lgLineRadPresOn &&
rfield.lgDoLineTrans )
271 pressure.pres_radiation_lines_curr = 0.;
276 if( lgMustEvaluateTrivial || Piso_seq[ipISO] > TrivialLineRadiationPressure )
278 Piso_seq[ipISO] = 0.;
279 for(
long nelem=ipISO; nelem <
LIMELM; nelem++ )
282 if(
dense.IonHigh[nelem] >= nelem + 1 - ipISO )
286 for(
long ipHi=1; ipHi <
iso_sp[ipISO][nelem].numLevels_local -
iso_sp[ipISO][nelem].nCollapsed_local; ipHi++ )
288 for(
long ipLo=0; ipLo < ipHi; ipLo++ )
290 if(
iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() <= 0 )
297 ( (
iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauTot() -
299 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pesc() > FLT_EPSILON*100. )
305 ipLineTypePradMax = 2;
310 ipLinePradMax = ipLo;
312 Piso_seq[ipISO] += RadPres1;
315 enum {DEBUG_LOC=
false};
316 if( DEBUG_LOC && ipISO==
ipH_LIKE && ipLo==3 && ipHi==5 &&
nzone > 260 )
319 "DEBUG prad1 \t%.2f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t\n",
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());
333 ASSERT( Piso_seq[ipISO] >= 0. );
335 pressure.pres_radiation_lines_curr += Piso_seq[ipISO];
339 enum {DEBUG_LOC=
false};
341 if( DEBUG_LOC &&
nzone > 260 )
344 "DEBUG prad2 \t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
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());
359 "DEBUG prad3\t%.2e\t%li\t%li\t%li\t%li\t%.2e\t%.2e\t%.2e\n",
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() );
371 if( lgMustEvaluateTrivial || PLevel1 > TrivialLineRadiationPressure )
375 for(
long i=1; i <=
nLevel1; i++ )
377 if( (*
TauLines[i].Hi()).Pop() > 1e-30 )
383 ipLineTypePradMax = 3;
392 pressure.pres_radiation_lines_curr += PLevel1;
394 if( lgMustEvaluateTrivial || PLevel2 > TrivialLineRadiationPressure )
402 if( (*
TauLine2[i].Hi()).Pop() > 1e-30 )
409 ipLineTypePradMax = 4;
418 pressure.pres_radiation_lines_curr += PLevel2;
421 if( lgMustEvaluateTrivial || PHfs > TrivialLineRadiationPressure )
426 if( (*
HFLines[i].Hi()).Pop() > 1e-30 )
433 ipLineTypePradMax = 8;
441 pressure.pres_radiation_lines_curr += PHfs;
444 if( lgMustEvaluateTrivial || P_H2 > TrivialLineRadiationPressure )
449 P_H2 += (*diatom)->H2_RadPress();
454 ipLineTypePradMax = 9;
459 pressure.pres_radiation_lines_curr += P_H2;
462 if( lgMustEvaluateTrivial || P_FeII > TrivialLineRadiationPressure )
468 ipLineTypePradMax = 7;
472 pressure.pres_radiation_lines_curr += P_FeII;
475 if( lgMustEvaluateTrivial || P_dBase > TrivialLineRadiationPressure )
478 for(
long ipSpecies=0; ipSpecies<
nSpecies; ipSpecies++ )
486 int ipHi = (*tr).ipHi();
489 int ipLo = (*tr).ipLo();
490 if( (*tr).ipCont() > 0 && (*(*tr).Hi()).Pop() > 1e-30 )
496 ipLineTypePradMax = 10;
500 ipLinePradMax = ipLo;
509 pressure.pres_radiation_lines_curr += P_dBase;
516 ipLineTypePradMax = 0;
526 if(
pressure.pres_radiation_lines_curr < 0. )
529 "PresTotCurrent: negative pressure, constituents follow %e %e %e %e %e \n",
541 if(
trace.lgTrace && ipLineTypePradMax <0 )
544 " PresTotCurrent, pressure.pbeta = %e, ipLineTypePradMax%li ipLinePradMax=%li \n",
545 pressure.pbeta,ipLineTypePradMax,ipLinePradMax );
551 phycon.EnergyExcitation = 0.;
563 for(
long nelem=ipISO; nelem <
LIMELM; nelem++ )
565 if(
dense.IonHigh[nelem] == nelem + 1 - ipISO )
568 for(
long ipHi=1; ipHi <
iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
570 phycon.EnergyExcitation +=
571 iso_sp[ipISO][nelem].st[ipHi].Pop() *
572 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyErg() *
582 phycon.EnergyExcitation += H2_InterEnergy();
587 for(
long i=1; i <=
nLevel1; i++ )
589 phycon.EnergyExcitation +=
596 phycon.EnergyExcitation +=
602 phycon.EnergyExcitation +=
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",
648 if(
dynamics.lgTimeDependentStatic )
681 pressure.pres_radiation_lines_curr =
691 pressure.ipPradMax_line = ipLinePradMax;
693 if( ipLineTypePradMax == 2 )
696 strcpy(
pressure.chLineRadPres,
"ISO ");
698 ASSERT( ipLinePradMax>=0 && ip2>=0 && ip3>=0 && ip4>=0 );
702 enum {DEBUG_LOC=
false};
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 )
718 fprintf(
ioQQQ,
"hit it\n");
720 fprintf(
ioQQQ,
"DEBUG width %e\n", width);
725 else if( ipLineTypePradMax == 3 )
728 ASSERT( ipLinePradMax>=0 );
729 strcpy(
pressure.chLineRadPres,
"Level1 ");
732 else if( ipLineTypePradMax == 4 )
735 ASSERT( ipLinePradMax>=0 );
736 strcpy(
pressure.chLineRadPres,
"Level2 ");
739 else if( ipLineTypePradMax == 5 )
743 else if( ipLineTypePradMax == 6 )
747 else if( ipLineTypePradMax == 7 )
750 strcpy(
pressure.chLineRadPres,
"Fe II ");
752 else if( ipLineTypePradMax == 8 )
755 strcpy(
pressure.chLineRadPres,
"hyperf ");
756 ASSERT( ipLinePradMax>=0 );
759 else if( ipLineTypePradMax == 9 )
762 strcpy(
pressure.chLineRadPres,
"large H2 ");
764 else if( ipLineTypePradMax == 10 )
767 strcpy(
pressure.chLineRadPres,
"dBaseLin " );
772 fprintf(
ioQQQ,
" PresTotCurrent ipLineTypePradMax set to %li, this is impossible.\n", ipLineTypePradMax);
781 " PresTotCurrent Pbeta:%.2f due to %s\n",
789 TotalPressure_v =
pressure.PresGasCurr;
808 TotalPressure_v +=
pressure.pres_radiation_lines_curr *
pressure.lgPres_radiation_ON;
811 enum{DEBUG_LOC=
false};
814 fprintf(
ioQQQ,
"pressureee%li\t%.4e\t%.4e\t%.4e\t%.3f\t%.3f\t%.3f\n",
826 if( TotalPressure_v< 0. )
831 fprintf(
ioQQQ,
" The negative pressure due to ordered magnetic field overwhelms the total outward pressure - please reconsider the geometry & field.\n");
835 ASSERT( TotalPressure_v > 0. );
841 pressure.PresTotlCurr = TotalPressure_v;