45 ipLineTypePradMax=-1 ,
47 ipLinePradMax=-LONG_MAX,
55 static double Piso_seq[
NISO]={0.,0.},
78 ipLinePradMax=-LONG_MAX;
92 bool lgMustEvaluate =
false;
96 bool lgMustEvaluateTrivial =
false;
106 lgMustEvaluate =
true;
107 lgMustEvaluateTrivial =
true;
111 static long int nzoneEvaluated=0, iterationEvaluated=0;
114 lgMustEvaluate =
true;
115 lgMustEvaluateTrivial =
true;
117 ipLineTypePradMax = -1;
124 lgMustEvaluate =
true;
130 nzoneEvaluated =
nzone;
154 for( ion=0; ion<=nelem+1; ++ion )
171 fprintf(
ioQQQ,
"\n\n PROBLEM DISASTER PressureTotal: the chemical species "
172 "are not conserved.\n");
173 fprintf(
ioQQQ,
" The sum of the densities of the atoms and ions for "
174 "%s is %.3e which is greater than the gas-phase abundance "
175 "of the element, %.3e.\n",
181 fprintf(
ioQQQ,
" The chemistry network is known to fail this way"
182 " in cold molecular environments when cosmic rays are not"
183 " included. They were not - try including them with the"
184 " COSMIC RAY BACKBROUND command.\n");
185 fprintf(
ioQQQ,
" The calculation is aborting. conv.nTotalIoniz=%li\n "
193 DenseAtomsIons += DenseOne;
200 ASSERT( DenseAtomsIons > 0. );
217 if( ipISO >= 0 && ipISO<
NISO )
219 for( i=1; i<=ion; ++i )
221 long int nelec = nelem+1 - i;
248 if(
COmole[i]->n_nuclei > 1)
265 fprintf(
ioQQQ,
"PROBLEM DISASTER pressure_total has found "
266 "dense.xNucleiTotal with an insane density.\n");
267 fprintf(
ioQQQ,
"The density was %.2e\n",
284 for( i=0; i <
LIMELM; i++ )
348 fprintf(
ioQQQ,
"DEBUG pressure_total updates AccelTot to %.2e grav %.2e\n",
358 fprintf(
ioQQQ,
"DEBUG widLya\t%li\t%.2f\t%.3e",
368 fprintf(
ioQQQ,
"\t%.3e\n",
381 " PresTotCurrent nzone %li iteration %li lgMustEvaluate:%c "
382 "lgMustEvaluateTrivial:%c "
383 "pressure.lgLineRadPresOn:%c "
384 "rfield.lgDoLineTrans:%c \n",
397 if( lgMustEvaluateTrivial || Piso_seq[ipISO] > TrivialLineRadiationPressure )
399 Piso_seq[ipISO] = 0.;
400 for( nelem=ipISO; nelem <
LIMELM; nelem++ )
409 for( ipLo=0; ipLo < ipHi; ipLo++ )
411 if(
Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 )
419 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc > FLT_EPSILON*100. )
425 ipLineTypePradMax = 2;
430 ipLinePradMax = ipLo;
432 Piso_seq[ipISO] += RadPres1;
435 enum {DEBUG_LOC=
false};
436 if( DEBUG_LOC && ipISO==
ipH_LIKE && ipLo==3 && ipHi==5 &&
nzone > 260 )
439 "DEBUG prad1 \t%.2f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t\n",
442 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc,
453 ASSERT( Piso_seq[ipISO] >= 0. );
459 enum {DEBUG_LOC=
false};
461 if( DEBUG_LOC &&
nzone > 260 )
464 "DEBUG prad2 \t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
467 Transitions[ipISO][ip3][ipLinePradMax][ip2].Emis->PopOpc,
479 "DEBUG prad3\t%.2e\t%li\t%li\t%li\t%li\t%.2e\t%.2e\t%.2e\n",
485 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->PopOpc,
487 1.-
Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->Pesc );
491 if( lgMustEvaluateTrivial || PLevel1 > TrivialLineRadiationPressure )
503 ipLineTypePradMax = 3;
514 if( lgMustEvaluateTrivial || PLevel2 > TrivialLineRadiationPressure )
529 ipLineTypePradMax = 4;
541 if( lgMustEvaluateTrivial || PHfs > TrivialLineRadiationPressure )
546 if(
HFLines[i].Hi->Pop > 1e-30 )
553 ipLineTypePradMax = 8;
564 if( lgMustEvaluateTrivial || P_H2 > TrivialLineRadiationPressure )
571 ipLineTypePradMax = 9;
578 if( lgMustEvaluateTrivial || P_FeII > TrivialLineRadiationPressure )
584 ipLineTypePradMax = 7;
591 if( lgMustEvaluateTrivial || PCO > TrivialLineRadiationPressure )
603 ipLineTypePradMax = 5;
615 ipLineTypePradMax = 6;
629 ipLineTypePradMax = 0;
639 "PresTotCurrent: negative pressure, constituents follow %e %e %e %e %e \n",
654 " PresTotCurrent, pressure.pbeta = %e, ipLineTypePradMax%li ipLinePradMax=%li \n",
673 for( nelem=ipISO; nelem <
LIMELM; nelem++ )
715 double Energy12 = 0.;
716 double Energy13 = 0.;
743 " PresTotCurrent zn %.2f Ptot:%.5e Pgas:%.3e Prad:%.3e Pmag:%.3e Pram:%.3e "
744 "gas parts; H:%.2e He:%.2e L1:%.2e L2:%.2e CO:%.2e fs%.2e h2%.2e turb%.2e\n",
815 if( ipLineTypePradMax == 2 )
819 ASSERT( ip4 < NISO && ip3<LIMELM );
820 ASSERT( ipLinePradMax>=0 && ip2>=0 && ip3>=0 && ip4>=0 );
824 enum {DEBUG_LOC=
false};
829 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",
831 ip4,ip3,ip2,ipLinePradMax,
832 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->TauIn,
833 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->TauTot,
834 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->Pesc,
835 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->Pelec_esc,
836 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->Pdest);
837 if( ip2==5 && ipLinePradMax==4 )
840 fprintf(
ioQQQ,
"hit it\n");
842 fprintf(
ioQQQ,
"DEBUG width %e\n", width);
847 else if( ipLineTypePradMax == 3 )
850 ASSERT( ipLinePradMax>=0 );
854 else if( ipLineTypePradMax == 4 )
857 ASSERT( ipLinePradMax>=0 );
861 else if( ipLineTypePradMax == 5 )
864 ASSERT( ipLinePradMax>=0 );
868 else if( ipLineTypePradMax == 6 )
871 ASSERT( ipLinePradMax>=0 );
875 else if( ipLineTypePradMax == 7 )
880 else if( ipLineTypePradMax == 8 )
884 ASSERT( ipLinePradMax>=0 );
887 else if( ipLineTypePradMax == 9 )
894 fprintf(
ioQQQ,
" PresTotCurrent ipLineTypePradMax set to %li, this is impossible.\n", ipLineTypePradMax);
903 " PresTotCurrent Pbeta:%.2f due to %s\n",
933 enum{DEBUG_LOC=
false};
936 fprintf(
ioQQQ,
"pressureee%li\t%.4e\t%.4e\t%.4e\t%.3f\t%.3f\t%.3f\n",
951 fprintf(
ioQQQ,
" The negative pressure due to ordered magnetic field overwhelms the total pressure - please reconsider the geometry & field.\n");
955 ASSERT( TotalPressure_v > 0. );
977 double PresTotlInitSave;
978 double PresRamInitSave;
983 PresRamInitSave =
pressure.PresRamInit;
987 PresTotlInitSave = 0.;
988 PresRamInitSave = 0.;
995 " PresTotCurrent 1st zn reset PresTotlInit to PresTotlCurr=%.3e "
996 "PresRamInit to PresRamCurr=%.3e old tot=%.3e old ram %.3e hden=%.3e\n",
999 PresTotlInitSave,PresRamInitSave,