00001
00002
00003
00004
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
00034 void PresTotCurrent(void)
00035 {
00036 static long int
00037
00038 ipLineTypePradMax=-1 ,
00039
00040 ipLinePradMax=-LONG_MAX,
00041 ip2=-LONG_MAX,
00042 ip3=-LONG_MAX,
00043 ip4=-LONG_MAX;
00044
00045
00046
00047
00048 static double Piso_seq[NISO]={0.,0.},
00049 PLevel1=0.,
00050 PLevel2=0.,
00051 PHfs=0.,
00052 PCO=0.,
00053 P_H2=0.,
00054 P_FeII=0.,
00055 P_dBase=0.;
00056
00057 double
00058 RadPres1,
00059 TotalPressure_v,
00060 pmx;
00061 realnum den_hmole ,
00062 den_comole ,
00063 DenseAtomsIons,
00064 DensElements;
00065
00066 DEBUG_ENTRY( "PresTotCurrent()" );
00067
00068 if( !conv.nTotalIoniz )
00069 {
00070
00071
00072 ipLinePradMax=-LONG_MAX;
00073
00074 }
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086 bool lgMustEvaluate = false;
00087
00088
00089
00090 bool lgMustEvaluateTrivial = false;
00091
00092
00093
00094 double TrivialLineRadiationPressure = conv.PressureErrorAllowed/10. *
00095 pressure.pres_radiation_lines_curr;
00096
00097
00098 if( conv.lgSearch )
00099 {
00100 lgMustEvaluate = true;
00101 lgMustEvaluateTrivial = true;
00102 }
00103
00104
00105 static long int nzoneEvaluated=0, iterationEvaluated=0;
00106 if( nzone!=nzoneEvaluated || iteration!=iterationEvaluated )
00107 {
00108 lgMustEvaluate = true;
00109 lgMustEvaluateTrivial = true;
00110
00111 ipLineTypePradMax = -1;
00112 }
00113
00114
00115
00116 if( (strcmp(dense.chDenseLaw,"WIND") == 0 ) ||
00117 (strcmp(dense.chDenseLaw,"CPRE") == 0 ) )
00118 lgMustEvaluate = true;
00119
00120 if( lgMustEvaluate )
00121 {
00122
00123
00124 nzoneEvaluated = nzone;
00125 iterationEvaluated = iteration;
00126 }
00127
00128
00129
00130
00131 TempChange(phycon.te , false);
00132
00133
00134 DenseAtomsIons = 0.;
00135 DensElements = 0.;
00136 realnum factor = 1.1f;
00137 if( conv.lgSearch )
00138 factor = 2.f;
00139 for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
00140 {
00141
00142
00143 if( dense.lgElmtOn[nelem] )
00144 {
00145
00146 DensElements += dense.gas_phase[nelem];
00147 realnum DenseOne = 0;
00148 for( long ion=0; ion<=nelem+1; ++ion )
00149 DenseOne += dense.xIonDense[nelem][ion];
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161 if( conv.nTotalIoniz>5 && !dense.lgSetIoniz[nelem] &&
00162 DenseOne > dense.gas_phase[nelem]*factor )
00163 {
00164
00165 fprintf(ioQQQ,"\n\n PROBLEM DISASTER PressureTotal: the chemical species "
00166 "are not conserved.\n");
00167 fprintf(ioQQQ," The sum of the densities of the atoms and ions for "
00168 "%s is %.3e which is greater than the gas-phase abundance "
00169 "of the element, %.3e.\n",
00170 elementnames.chElementName[nelem],DenseOne,
00171 dense.gas_phase[nelem]);
00172
00173 if( hextra.cryden == 0. )
00174 {
00175 fprintf(ioQQQ," The chemistry network is known to fail this way"
00176 " in cold molecular environments when cosmic rays are not"
00177 " included. They were not - try including them with the"
00178 " COSMIC RAY BACKBROUND command.\n");
00179 fprintf(ioQQQ," The calculation is aborting. conv.nTotalIoniz=%li\n"
00180 " Sorry.\n\n", conv.nTotalIoniz );
00181 lgAbort = true;
00182 }
00183
00184 else
00185 TotalInsanity();
00186 }
00187 DenseAtomsIons += DenseOne;
00188 }
00189 }
00190
00191
00192
00193 ASSERT( conv.nTotalIoniz<5 || DenseAtomsIons <= DensElements*factor );
00194 ASSERT( DenseAtomsIons > 0. );
00195
00196
00197
00198 phycon.EnergyIonization = 0.;
00199 #if 0
00200 for( long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
00201 {
00202 for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
00203 {
00204
00205 int kadvec = dynamics.lgMETALS;
00206
00207
00208
00209 ipISO = nelem-ion;
00210 fixit();
00211 if( ipISO >= 0 && ipISO<NISO )
00212 kadvec = dynamics.lgISO[ipISO];
00213 for( long i=1; i<=ion; ++i )
00214 {
00215 long int nelec = nelem+1 - i;
00216
00217
00218 phycon.EnergyIonization += dense.xIonDense[nelem][ion] *
00219 t_ADfA::Inst().ph1(Heavy.nsShells[nelem][i-1]-1,nelec,nelem,0)/EVRYD* 0.9998787*kadvec;
00220 }
00221 }
00222 }
00223
00224 phycon.EnergyIonization *= EN1RYD;
00225 #endif
00226
00230 phycon.EnergyBinding = 0.;
00231
00232
00233
00234
00235 den_comole = 0.;
00236 for( long i=0; i < mole.num_comole_calc; i++ )
00237 {
00238
00239
00240
00241
00242 if(COmole[i]->n_nuclei > 1)
00243 den_comole += COmole[i]->hevmol * COmole[i]->lgGas_Phase;
00244 }
00245
00246
00247
00248 den_hmole = 0.;
00249 for( long i=0; i<N_H_MOLEC; ++i )
00250 {
00251 if( i!=ipMH && i!=ipMHp )
00252 den_hmole += hmi.Hmolec[i];
00253 }
00254
00255
00256 dense.xNucleiTotal = den_hmole + den_comole + DenseAtomsIons;
00257 if( dense.xNucleiTotal > BIGFLOAT )
00258 {
00259 fprintf(ioQQQ, "PROBLEM DISASTER pressure_total has found "
00260 "dense.xNucleiTotal with an insane density.\n");
00261 fprintf(ioQQQ, "The density was %.2e\n",
00262 dense.xNucleiTotal);
00263 TotalInsanity();
00264 }
00265 ASSERT( dense.xNucleiTotal > 0. );
00266
00267
00268
00269 dense.pden = (realnum)(dense.eden + dense.xNucleiTotal);
00270
00271
00272 pressure.PresGasCurr = dense.pden*phycon.te*BOLTZMANN;
00273
00274
00275
00276
00277 dense.wmole = 0.;
00278 for( long i=0; i < LIMELM; i++ )
00279 {
00280 dense.wmole += dense.gas_phase[i]*(realnum)dense.AtomicWeight[i];
00281 }
00282
00283
00284 dense.wmole /= dense.pden;
00285 ASSERT( dense.wmole > 0. && dense.pden > 0. );
00286
00287
00290 dense.xMassDensity = (realnum)(dense.wmole*ATOMIC_MASS_UNIT*dense.pden);
00291
00292
00293
00294
00295
00296
00297
00298 if( dense.xMassDensity0 < 0.0 )
00299 dense.xMassDensity0 = dense.xMassDensity;
00300
00301
00302
00303
00304
00305 if(!(wind.lgBallistic() || wind.lgStatic()))
00306 {
00307 wind.windv = DynaFlux(radius.depth)/(dense.xMassDensity);
00308 }
00309
00310
00311 pressure.PresRamCurr = dense.xMassDensity*POW2( wind.windv );
00312
00313
00314
00315 pressure.PresTurbCurr = dense.xMassDensity*POW2( DoppVel.TurbVel ) *
00316 DoppVel.Heiles_Troland_F / 6.;
00317
00320
00321 if( lgMustEvaluate )
00322 {
00323
00324 double rforce = 0.;
00325 double relec = 0.;
00326 for( long i=0; i < rfield.nflux; i++ )
00327 {
00328 rforce += (rfield.flux[0][i] + rfield.outlin[0][i] + rfield.outlin_noplot[i]+ rfield.ConInterOut[i])*
00329 rfield.anu[i]*(opac.opacity_abs[i] + opac.opacity_sct[i]);
00330
00331
00332 relec +=
00333 (rfield.flux[0][i] + rfield.outlin[0][i] + rfield.outlin_noplot[i]+
00334 rfield.ConInterOut[i]) *
00335 opac.OpacStack[i-1+opac.iopcom]*dense.eden*rfield.anu[i];
00336 }
00337
00338 wind.AccelCont = (realnum)(rforce*EN1RYD/SPEEDLIGHT/dense.xMassDensity);
00339
00340 wind.AccelElectron = (realnum)(relec*EN1RYD/SPEEDLIGHT/dense.xMassDensity);
00341
00342
00343 wind.AccelLine = (realnum)(RT_line_driving()/SPEEDLIGHT/dense.xMassDensity);
00344
00345
00346 wind.AccelTotalOutward = wind.AccelCont + wind.AccelLine;
00347
00348
00349
00350
00351
00352 if( wind.AccelElectron > SMALLFLOAT )
00353 wind.fmul = (realnum)( (wind.AccelLine + wind.AccelCont) / wind.AccelElectron);
00354 else
00355 wind.fmul = 0.;
00356
00357 double reff = radius.Radius-radius.drad/2.;
00358
00359 wind.AccelGravity = (realnum)(
00360
00361 GRAV_CONST*wind.comass*SOLAR_MASS/POW2(reff)*
00362
00363 (1.-wind.DiskRadius/reff) );
00364 # if 0
00365 if( fudge(-1) )
00366 fprintf(ioQQQ,"DEBUG pressure_total updates AccelTotalOutward to %.2e grav %.2e\n",
00367 wind.AccelTotalOutward , wind.AccelGravity );
00368 # endif
00369 }
00370
00371
00372 if( Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->PopOpc > 0. )
00373 {
00374 hydro.HLineWidth = (realnum)RT_LineWidth(&Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s],GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN]));
00375 # if 0
00376 fprintf(ioQQQ,"DEBUG widLya\t%li\t%.2f\t%.3e",
00377 iteration,
00378 fnzone,
00379 hydro.HLineWidth);
00380 hydro.HLineWidth = (realnum)(
00381 RT_LyaWidth(
00382 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauIn,
00383 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].TauTot,
00384 4.72e-2/phycon.sqrte,
00385 GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN]) ) );
00386 fprintf(ioQQQ,"\t%.3e\n",
00387 hydro.HLineWidth);
00388 # endif
00389 }
00390 else
00391 hydro.HLineWidth = 0.;
00392
00393
00394
00395
00396 if( trace.lgTrace )
00397 {
00398 fprintf(ioQQQ,
00399 " PresTotCurrent nzone %li iteration %li lgMustEvaluate:%c "
00400 "lgMustEvaluateTrivial:%c "
00401 "pressure.lgLineRadPresOn:%c "
00402 "rfield.lgDoLineTrans:%c \n",
00403 nzone , iteration , TorF(lgMustEvaluate) , TorF(lgMustEvaluateTrivial),
00404 TorF(pressure.lgLineRadPresOn), TorF(rfield.lgDoLineTrans) );
00405 }
00406
00407 if( lgMustEvaluate && pressure.lgLineRadPresOn && rfield.lgDoLineTrans )
00408 {
00409
00410 pressure.pres_radiation_lines_curr = 0.;
00411
00412 pmx = 0.;
00413 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00414 {
00415 if( lgMustEvaluateTrivial || Piso_seq[ipISO] > TrivialLineRadiationPressure )
00416 {
00417 Piso_seq[ipISO] = 0.;
00418 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
00419 {
00420
00421 if( dense.IonHigh[nelem] >= nelem + 1 - ipISO )
00422 {
00423
00424
00425 for( long ipHi=1; ipHi <iso.numLevels_local[ipISO][nelem] - iso.nCollapsed_local[ipISO][nelem]; ipHi++ )
00426 {
00427 for( long ipLo=0; ipLo < ipHi; ipLo++ )
00428 {
00429 if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont <= 0 )
00430 continue;
00431
00432 ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul > iso.SmallA );
00433
00434 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc > SMALLFLOAT &&
00435
00436 ( (Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauTot - Transitions[ipISO][nelem][ipHi][ipLo].Emis->TauIn) > SMALLFLOAT ) &&
00437 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc > FLT_EPSILON*100. )
00438 {
00439 RadPres1 = PressureRadiationLine( &Transitions[ipISO][nelem][ipHi][ipLo], GetDopplerWidth(dense.AtomicWeight[nelem]) );
00440
00441 if( RadPres1 > pmx )
00442 {
00443 ipLineTypePradMax = 2;
00444 pmx = RadPres1;
00445 ip4 = ipISO;
00446 ip3 = nelem;
00447 ip2 = ipHi;
00448 ipLinePradMax = ipLo;
00449 }
00450 Piso_seq[ipISO] += RadPres1;
00451 {
00452
00453 enum {DEBUG_LOC=false};
00454 if( DEBUG_LOC && ipISO==ipH_LIKE && ipLo==3 && ipHi==5 && nzone > 260 )
00455 {
00456 fprintf(ioQQQ,
00457 "DEBUG prad1 \t%.2f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t\n",
00458 fnzone,
00459 RadPres1,
00460 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc,
00461 StatesElemNEW[nelem][nelem-ipISO][ipLo].Pop,
00462 StatesElemNEW[nelem][nelem-ipISO][ipHi].Pop,
00463 Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pesc);
00464 }
00465 }
00466 }
00467 }
00468 }
00469 }
00470 }
00471 ASSERT( Piso_seq[ipISO] >= 0. );
00472 }
00473 pressure.pres_radiation_lines_curr += Piso_seq[ipISO];
00474 }
00475 {
00476
00477 enum {DEBUG_LOC=false};
00478 # if 0
00479 if( DEBUG_LOC && nzone > 260 )
00480 {
00481 fprintf(ioQQQ,
00482 "DEBUG prad2 \t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00483 nzone,
00484 pmx,
00485 Transitions[ipISO][ip3][ipLinePradMax][ip2].Emis->PopOpc,
00486 StatesElemNEW[ip3][ip3-ipISO][0].Pop,
00487 StatesElemNEW[ip3][ip3-ipISO][2].Pop,
00488 StatesElemNEW[ip3][ip3-ipISO][3].Pop,
00489 StatesElemNEW[ip3][ip3-ipISO][4].Pop,
00490 StatesElemNEW[ip3][ip3-ipISO][5].Pop,
00491 StatesElemNEW[ip3][ip3-ipISO][6].Pop);
00492 }
00493 # endif
00494 if( DEBUG_LOC )
00495 {
00496 fprintf(ioQQQ,
00497 "DEBUG prad3\t%.2e\t%li\t%li\t%li\t%li\t%.2e\t%.2e\t%.2e\n",
00498 pmx,
00499 ip4,
00500 ip3,
00501 ip2,
00502 ipLinePradMax,
00503 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->PopOpc,
00504 StatesElemNEW[ip3][ip3-ip4][ip2].Pop,
00505 1.-Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->Pesc );
00506 }
00507 }
00508
00509 if( lgMustEvaluateTrivial || PLevel1 > TrivialLineRadiationPressure )
00510 {
00511 PLevel1 = 0.;
00512
00513 for( long i=1; i <= nLevel1; i++ )
00514 {
00515 if( TauLines[i].Hi->Pop > 1e-30 )
00516 {
00517 RadPres1 = PressureRadiationLine( &TauLines[i], GetDopplerWidth(dense.AtomicWeight[TauLines[i].Hi->nelem-1]) );
00518
00519 if( RadPres1 > pmx )
00520 {
00521 ipLineTypePradMax = 3;
00522 pmx = RadPres1;
00523 ipLinePradMax = i;
00524 }
00525 PLevel1 += RadPres1;
00526 }
00527 }
00528 ASSERT( PLevel1 >= 0. );
00529 }
00530 pressure.pres_radiation_lines_curr += PLevel1;
00531
00532 if( lgMustEvaluateTrivial || PLevel2 > TrivialLineRadiationPressure )
00533 {
00534
00535 PLevel2 = 0.;
00536 for( long i=0; i < nWindLine; i++ )
00537 {
00538 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
00539 {
00540 if( TauLine2[i].Hi->Pop > 1e-30 )
00541 {
00542 RadPres1 = PressureRadiationLine( &TauLine2[i], GetDopplerWidth(dense.AtomicWeight[TauLine2[i].Hi->nelem-1]) );
00543
00544 PLevel2 += RadPres1;
00545 if( RadPres1 > pmx )
00546 {
00547 ipLineTypePradMax = 4;
00548 pmx = RadPres1;
00549 ipLinePradMax = i;
00550 }
00551 }
00552 }
00553 }
00554 ASSERT( PLevel2 >= 0. );
00555 }
00556 pressure.pres_radiation_lines_curr += PLevel2;
00557
00558
00559 if( lgMustEvaluateTrivial || PHfs > TrivialLineRadiationPressure )
00560 {
00561 PHfs = 0.;
00562 for( long i=0; i < nHFLines; i++ )
00563 {
00564 if( HFLines[i].Hi->Pop > 1e-30 )
00565 {
00566 RadPres1 = PressureRadiationLine( &HFLines[i], GetDopplerWidth(dense.AtomicWeight[HFLines[i].Hi->nelem-1]) );
00567
00568 PHfs += RadPres1;
00569 if( RadPres1 > pmx )
00570 {
00571 ipLineTypePradMax = 8;
00572 pmx = RadPres1;
00573 ipLinePradMax = i;
00574 }
00575 }
00576 }
00577 ASSERT( PHfs >= 0. );
00578 }
00579 pressure.pres_radiation_lines_curr += PHfs;
00580
00581
00582 if( lgMustEvaluateTrivial || P_H2 > TrivialLineRadiationPressure )
00583 {
00584 P_H2 = H2_RadPress();
00585
00586 if( P_H2 > pmx )
00587 {
00588 pmx = P_H2;
00589 ipLineTypePradMax = 9;
00590 }
00591 ASSERT( P_H2 >= 0. );
00592 }
00593 pressure.pres_radiation_lines_curr += P_H2;
00594
00595
00596 if( lgMustEvaluateTrivial || P_FeII > TrivialLineRadiationPressure )
00597 {
00598 P_FeII = FeIIRadPress();
00599 if( P_FeII > pmx )
00600 {
00601 pmx = P_FeII;
00602 ipLineTypePradMax = 7;
00603 }
00604 ASSERT( P_FeII >= 0. );
00605 }
00606 pressure.pres_radiation_lines_curr += P_FeII;
00607
00608
00609 if( lgMustEvaluateTrivial || P_dBase > TrivialLineRadiationPressure )
00610 {
00611 P_dBase = 0.;
00612 for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
00613 {
00614 if( Species[ipSpecies].lgActive )
00615 {
00616 for( long ipHi=1; ipHi<Species[ipSpecies].numLevels_local; ipHi++ )
00617 {
00618 transition *t = &dBaseTrans[ipSpecies][ipHi][0];
00619 for(long ipLo = 0; ipLo < ipHi; ipLo++ )
00620 {
00621 if( t->ipCont > 0 && t->Hi->Pop > 1e-30 ) {
00622 RadPres1 = PressureRadiationLine( t, GetDopplerWidth( t->Hi->sp->fmolweight ) );
00623
00624 if( RadPres1 > pmx )
00625 {
00626 ipLineTypePradMax = 10;
00627 pmx = RadPres1;
00628 ip3 = ipSpecies;
00629 ip2 = ipHi;
00630 ipLinePradMax = ipLo;
00631 }
00632 P_dBase += RadPres1;
00633 }
00634 ++t;
00635 }
00636 }
00637 }
00638 }
00639 ASSERT( P_dBase >= 0. );
00640 }
00641 pressure.pres_radiation_lines_curr += P_dBase;
00642
00643 }
00644 else if( !pressure.lgLineRadPresOn || !rfield.lgDoLineTrans )
00645 {
00646
00647 ipLinePradMax = -1;
00648 ipLineTypePradMax = 0;
00649 }
00650
00651 fixit();
00652
00653
00654
00655 pressure.pbeta = (realnum)(pressure.pres_radiation_lines_curr/SDIV(pressure.PresTotlCurr));
00656
00657
00658 if( pressure.pres_radiation_lines_curr < 0. )
00659 {
00660 fprintf(ioQQQ,
00661 "PresTotCurrent: negative pressure, constituents follow %e %e %e %e %e \n",
00662 Piso_seq[ipH_LIKE],
00663 Piso_seq[ipHE_LIKE],
00664 PLevel1,
00665 PLevel2,
00666 PCO);
00667 ShowMe();
00668 cdEXIT(EXIT_FAILURE);
00669 }
00670
00671
00672
00673 if( trace.lgTrace && ipLineTypePradMax <0 )
00674 {
00675 fprintf(ioQQQ,
00676 " PresTotCurrent, pressure.pbeta = %e, ipLineTypePradMax%li ipLinePradMax=%li \n",
00677 pressure.pbeta,ipLineTypePradMax,ipLinePradMax );
00678 }
00679
00680
00681
00682
00683 phycon.EnergyExcitation = 0.;
00684 #if 0
00685 fixit();
00686
00687 broken();
00688
00689
00690
00691
00692 ipLo = 0;
00693 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00694 {
00695 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
00696 {
00697 if( dense.IonHigh[nelem] == nelem + 1 - ipISO )
00698 {
00699
00700 for( long ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ )
00701 {
00702 phycon.EnergyExcitation +=
00703 StatesElemNEW[nelem][nelem-ipISO][ipHi].Pop *
00704 Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg*
00705
00706 dynamics.lgISO[ipISO];
00707 }
00708 }
00709 }
00710 }
00711
00712 if( dynamics.lgISO[ipH_LIKE] )
00713
00714 phycon.EnergyExcitation += H2_InterEnergy();
00715
00716
00717 if( dynamics.lgMETALS )
00718 {
00719 for( long i=1; i <= nLevel1; i++ )
00720 {
00721 phycon.EnergyExcitation +=
00722 TauLines[i].Hi->Pop* TauLines[i].EnergyErg;
00723 }
00724 for( long i=0; i < nWindLine; i++ )
00725 {
00726 if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
00727 {
00728 phycon.EnergyExcitation +=
00729 TauLine2[i].Hi->Pop* TauLine2[i].EnergyErg;
00730 }
00731 }
00732 for( long i=0; i < nHFLines; i++ )
00733 {
00734 phycon.EnergyExcitation +=
00735 HFLines[i].Hi->Pop* HFLines[i].EnergyErg;
00736 }
00737
00738
00739 phycon.EnergyExcitation += FeII_InterEnergy();
00740 }
00741 #endif
00742
00743
00744
00745
00746
00747 Magnetic_evaluate();
00748
00749
00750 if( trace.lgTrace && (pressure.PresTotlCurr > SMALLFLOAT) )
00751 {
00752 fprintf(ioQQQ,
00753 " PresTotCurrent zn %.2f Ptot:%.5e Pgas:%.3e Prad:%.3e Pmag:%.3e Pram:%.3e "
00754 "gas parts; H:%.2e He:%.2e L1:%.2e L2:%.2e CO:%.2e fs%.2e h2%.2e turb%.2e\n",
00755 fnzone,
00756 pressure.PresTotlCurr,
00757 pressure.PresGasCurr/pressure.PresTotlCurr,
00758 pressure.pres_radiation_lines_curr*pressure.lgPres_radiation_ON/pressure.PresTotlCurr,
00759 magnetic.pressure/pressure.PresTotlCurr,
00760 pressure.PresRamCurr*pressure.lgPres_ram_ON/pressure.PresTotlCurr,
00761 Piso_seq[ipH_LIKE]/pressure.PresTotlCurr,
00762 Piso_seq[ipHE_LIKE]/pressure.PresTotlCurr,
00763 PLevel1/pressure.PresTotlCurr,
00764 PLevel2/pressure.PresTotlCurr,
00765 PCO/pressure.PresTotlCurr,
00766 PHfs/pressure.PresTotlCurr,
00767 P_H2/pressure.PresTotlCurr,
00768 pressure.PresTurbCurr*DoppVel.lgTurb_pressure/pressure.PresTotlCurr);
00769
00770 }
00771
00772
00773
00774
00775
00776
00777
00778
00779
00780 if( dynamics.lgTimeDependentStatic )
00781 {
00782
00783
00784
00785
00786
00787 phycon.EnthalpyDensity =
00788 0.5*POW2(wind.windv)*dense.xMassDensity +
00789 3./2.*pressure.PresGasCurr +
00790 magnetic.EnthalpyDensity;
00791
00792 }
00793 else
00794 {
00795
00796
00797
00798
00799 phycon.EnthalpyDensity =
00800 0.5*POW2(wind.windv)*dense.xMassDensity +
00801 5./2.*pressure.PresGasCurr +
00802 magnetic.EnthalpyDensity;
00803
00804 }
00805
00806
00807
00808
00809
00810 if( iteration <= 1 && pressure.pres_radiation_lines_curr >= pressure.PresGasCurr )
00811 {
00812
00813 pressure.pres_radiation_lines_curr =
00814 MIN2(pressure.pres_radiation_lines_curr,pressure.PresGasCurr);
00815 pressure.lgPradCap = true;
00816 }
00817
00818
00819
00820 if( pressure.pbeta > pressure.RadBetaMax && nzone )
00821 {
00822 pressure.RadBetaMax = pressure.pbeta;
00823 pressure.ipPradMax_line = ipLinePradMax;
00824 pressure.ipPradMax_nzone = nzone;
00825 if( ipLineTypePradMax == 2 )
00826 {
00827
00828 strcpy( pressure.chLineRadPres, "ISO ");
00829 ASSERT( ip4 < NISO && ip3<LIMELM );
00830 ASSERT( ipLinePradMax>=0 && ip2>=0 && ip3>=0 && ip4>=0 );
00831 strcat( pressure.chLineRadPres, chLineLbl(&Transitions[ip4][ip3][ip2][ipLinePradMax]) );
00832 {
00833
00834 enum {DEBUG_LOC=false};
00835
00836
00837 if( DEBUG_LOC )
00838 {
00839 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",
00840 iteration, nzone,
00841 ip4,ip3,ip2,ipLinePradMax,
00842 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->TauIn,
00843 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->TauTot,
00844 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->Pesc,
00845 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->Pelec_esc,
00846 Transitions[ip4][ip3][ip2][ipLinePradMax].Emis->Pdest);
00847 if( ip2==5 && ipLinePradMax==4 )
00848 {
00849 double width;
00850 fprintf(ioQQQ,"hit it\n");
00851 width = RT_LineWidth(&Transitions[ip4][ip3][ip2][ipLinePradMax],GetDopplerWidth(dense.AtomicWeight[ip3]));
00852 fprintf(ioQQQ,"DEBUG width %e\n", width);
00853 }
00854 }
00855 }
00856 }
00857 else if( ipLineTypePradMax == 3 )
00858 {
00859
00860 ASSERT( ipLinePradMax>=0 );
00861 strcpy( pressure.chLineRadPres, "Level1 ");
00862 strcat( pressure.chLineRadPres, chLineLbl(&TauLines[ipLinePradMax]) );
00863 }
00864 else if( ipLineTypePradMax == 4 )
00865 {
00866
00867 ASSERT( ipLinePradMax>=0 );
00868 strcpy( pressure.chLineRadPres, "Level2 ");
00869 strcat( pressure.chLineRadPres, chLineLbl(&TauLine2[ipLinePradMax]) );
00870 }
00871 else if( ipLineTypePradMax == 5 )
00872 {
00873 cdEXIT( EXIT_FAILURE );
00874 }
00875 else if( ipLineTypePradMax == 6 )
00876 {
00877 cdEXIT( EXIT_FAILURE );
00878 }
00879 else if( ipLineTypePradMax == 7 )
00880 {
00881
00882 strcpy( pressure.chLineRadPres, "Fe II ");
00883 }
00884 else if( ipLineTypePradMax == 8 )
00885 {
00886
00887 strcpy( pressure.chLineRadPres, "hyperf ");
00888 ASSERT( ipLinePradMax>=0 );
00889 strcat( pressure.chLineRadPres, chLineLbl(&HFLines[ipLinePradMax]) );
00890 }
00891 else if( ipLineTypePradMax == 9 )
00892 {
00893
00894 strcpy( pressure.chLineRadPres, "large H2 ");
00895 }
00896 else if( ipLineTypePradMax == 10 )
00897 {
00898
00899 strcpy( pressure.chLineRadPres, "dBaseLin " );
00900 strcat( pressure.chLineRadPres, dBaseTrans[ip3][ip2][ipLinePradMax].Hi->sp->chLabel );
00901 }
00902 else
00903 {
00904 fprintf(ioQQQ," PresTotCurrent ipLineTypePradMax set to %li, this is impossible.\n", ipLineTypePradMax);
00905 ShowMe();
00906 cdEXIT(EXIT_FAILURE);
00907 }
00908 }
00909
00910 if( trace.lgTrace && pressure.pbeta > 0.5 )
00911 {
00912 fprintf(ioQQQ,
00913 " PresTotCurrent Pbeta:%.2f due to %s\n",
00914 pressure.pbeta ,
00915 pressure.chLineRadPres
00916 );
00917 }
00918
00919
00920
00921 TotalPressure_v = pressure.PresGasCurr;
00922
00923
00924
00925 TotalPressure_v += pressure.PresRamCurr * pressure.lgPres_ram_ON;
00926
00927
00928
00931 TotalPressure_v += magnetic.pressure * pressure.lgPres_magnetic_ON;
00932
00933
00934
00935 TotalPressure_v += pressure.PresTurbCurr * DoppVel.lgTurb_pressure;
00936
00937
00938
00939
00940 TotalPressure_v += pressure.pres_radiation_lines_curr * pressure.lgPres_radiation_ON;
00941
00942 {
00943 enum{DEBUG_LOC=false};
00944 if( DEBUG_LOC && pressure.PresTotlCurr > SMALLFLOAT )
00945 {
00946 fprintf(ioQQQ,"pressureee%li\t%.4e\t%.4e\t%.4e\t%.3f\t%.3f\t%.3f\n",
00947 nzone,
00948 pressure.PresTotlCorrect,
00949 pressure.PresTotlCurr,
00950 TotalPressure_v ,
00951 pressure.PresGasCurr/pressure.PresTotlCurr,
00952 pressure.pres_radiation_lines_curr/pressure.PresTotlCurr ,
00953 pressure.PresRamCurr/pressure.PresTotlCurr
00954 );
00955 }
00956 }
00957
00958 if( TotalPressure_v< 0. )
00959 {
00960 ASSERT( magnetic.pressure < 0. );
00961
00962
00963 fprintf(ioQQQ," The negative pressure due to ordered magnetic field overwhelms the total outward pressure - please reconsider the geometry & field.\n");
00964 cdEXIT(EXIT_FAILURE);
00965 }
00966
00967 ASSERT( TotalPressure_v > 0. );
00968
00969
00970 pressure.PresMax = MAX2(pressure.PresMax,(realnum)TotalPressure_v);
00971
00972
00973 pressure.PresTotlCurr = TotalPressure_v;
00974
00975 return;
00976 }