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
00062 DEBUG_ENTRY( "PresTotCurrent()" );
00063
00064 if( !conv.nTotalIoniz )
00065 {
00066
00067
00068 ipLinePradMax=-LONG_MAX;
00069
00070 }
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082 bool lgMustEvaluate = false;
00083
00084
00085
00086 bool lgMustEvaluateTrivial = false;
00087
00088
00089
00090 double TrivialLineRadiationPressure = conv.PressureErrorAllowed/10. *
00091 pressure.pres_radiation_lines_curr;
00092
00093
00094 if( conv.lgSearch )
00095 {
00096 lgMustEvaluate = true;
00097 lgMustEvaluateTrivial = true;
00098 }
00099
00100
00101 static long int nzoneEvaluated=0, iterationEvaluated=0;
00102 if( nzone!=nzoneEvaluated || iteration!=iterationEvaluated )
00103 {
00104 lgMustEvaluate = true;
00105 lgMustEvaluateTrivial = true;
00106
00107 ipLineTypePradMax = -1;
00108 }
00109
00110
00111
00112 if( (strcmp(dense.chDenseLaw,"WIND") == 0 ) ||
00113 (strcmp(dense.chDenseLaw,"DYNA") == 0 ) ||
00114 (strcmp(dense.chDenseLaw,"CPRE") == 0 ) )
00115 lgMustEvaluate = true;
00116
00117 if( lgMustEvaluate )
00118 {
00119
00120
00121 nzoneEvaluated = nzone;
00122 iterationEvaluated = iteration;
00123 }
00124
00125
00126
00127
00128 TempChange(phycon.te , false);
00129
00130 ASSERT(lgElemsConserved());
00131
00132
00133
00134 phycon.EnergyIonization = 0.;
00135 #if 0
00136 for( long nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
00137 {
00138 for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
00139 {
00140
00141 int kadvec = dynamics.lgMETALS;
00142
00143
00144
00145 ipISO = nelem-ion;
00146 fixit();
00147 if( ipISO >= 0 && ipISO<NISO )
00148 kadvec = dynamics.lgISO[ipISO];
00149 for( long i=1; i<=ion; ++i )
00150 {
00151 long int nelec = nelem+1 - i;
00152
00153
00154 phycon.EnergyIonization += dense.xIonDense[nelem][ion] *
00155 t_ADfA::Inst().ph1(Heavy.nsShells[nelem][i-1]-1,nelec,nelem,0)/EVRYD* 0.9998787*kadvec;
00156 }
00157 }
00158 }
00159
00160 phycon.EnergyIonization *= EN1RYD;
00161 #endif
00162
00166 phycon.EnergyBinding = 0.;
00167
00168
00169 SumDensities();
00170
00171
00172 pressure.PresGasCurr = dense.pden*phycon.te*BOLTZMANN;
00173
00174
00175
00176
00177
00178
00179
00180 if(!(wind.lgBallistic() || wind.lgStatic()))
00181 {
00182 wind.windv = DynaFlux(radius.depth)/(dense.xMassDensity);
00183 }
00184
00185
00186 pressure.PresRamCurr = dense.xMassDensity*POW2( wind.windv );
00187
00188
00189
00190 pressure.PresTurbCurr = dense.xMassDensity*POW2( DoppVel.TurbVel ) *
00191 DoppVel.Heiles_Troland_F / 6.;
00192
00195
00196 if( lgMustEvaluate )
00197 {
00198
00199 double rforce = 0.;
00200 double relec = 0.;
00201 for( long i=0; i < rfield.nflux; i++ )
00202 {
00203 rforce += (rfield.flux[0][i] + rfield.outlin[0][i] + rfield.outlin_noplot[i]+ rfield.ConInterOut[i])*
00204 rfield.anu[i]*(opac.opacity_abs[i] + opac.opacity_sct[i]);
00205
00206
00207 relec +=
00208 (rfield.flux[0][i] + rfield.outlin[0][i] + rfield.outlin_noplot[i]+
00209 rfield.ConInterOut[i]) *
00210 opac.OpacStack[i-1+opac.iopcom]*dense.eden*rfield.anu[i];
00211 }
00212
00213 wind.AccelCont = (realnum)(rforce*EN1RYD/SPEEDLIGHT/dense.xMassDensity);
00214
00215 wind.AccelElectron = (realnum)(relec*EN1RYD/SPEEDLIGHT/dense.xMassDensity);
00216
00217
00218 wind.AccelLine = (realnum)(RT_line_driving()/SPEEDLIGHT/dense.xMassDensity);
00219
00220
00221 wind.AccelTotalOutward = wind.AccelCont + wind.AccelLine;
00222
00223
00224
00225
00226
00227 if( wind.AccelElectron > SMALLFLOAT )
00228 wind.fmul = (realnum)( (wind.AccelLine + wind.AccelCont) / wind.AccelElectron);
00229 else
00230 wind.fmul = 0.;
00231
00232 double reff = radius.Radius-radius.drad/2.;
00233
00234 wind.AccelGravity = (realnum)(
00235
00236 GRAV_CONST*wind.comass*SOLAR_MASS/POW2(reff)*
00237
00238 (1.-wind.DiskRadius/reff) );
00239 # if 0
00240 if( fudge(-1) )
00241 fprintf(ioQQQ,"DEBUG pressure_total updates AccelTotalOutward to %.2e grav %.2e\n",
00242 wind.AccelTotalOutward , wind.AccelGravity );
00243 # endif
00244 }
00245
00246
00247 if( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().PopOpc() > 0. )
00248 {
00249 hydro.HLineWidth = (realnum)RT_LineWidth(iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s),GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN]));
00250 }
00251 else
00252 hydro.HLineWidth = 0.;
00253
00254
00255
00256
00257 if( trace.lgTrace )
00258 {
00259 fprintf(ioQQQ,
00260 " PresTotCurrent nzone %li iteration %li lgMustEvaluate:%c "
00261 "lgMustEvaluateTrivial:%c "
00262 "pressure.lgLineRadPresOn:%c "
00263 "rfield.lgDoLineTrans:%c \n",
00264 nzone , iteration , TorF(lgMustEvaluate) , TorF(lgMustEvaluateTrivial),
00265 TorF(pressure.lgLineRadPresOn), TorF(rfield.lgDoLineTrans) );
00266 }
00267
00268 if( lgMustEvaluate && pressure.lgLineRadPresOn && rfield.lgDoLineTrans )
00269 {
00270
00271 pressure.pres_radiation_lines_curr = 0.;
00272
00273 pmx = 0.;
00274 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00275 {
00276 if( lgMustEvaluateTrivial || Piso_seq[ipISO] > TrivialLineRadiationPressure )
00277 {
00278 Piso_seq[ipISO] = 0.;
00279 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
00280 {
00281
00282 if( dense.IonHigh[nelem] >= nelem + 1 - ipISO )
00283 {
00284
00285
00286 for( long ipHi=1; ipHi <iso_sp[ipISO][nelem].numLevels_local - iso_sp[ipISO][nelem].nCollapsed_local; ipHi++ )
00287 {
00288 for( long ipLo=0; ipLo < ipHi; ipLo++ )
00289 {
00290 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() <= 0 )
00291 continue;
00292
00293 ASSERT( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() > iso_ctrl.SmallA );
00294
00295 if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc() > SMALLFLOAT &&
00296
00297 ( (iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauTot() -
00298 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn()) > SMALLFLOAT ) &&
00299 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pesc() > FLT_EPSILON*100. )
00300 {
00301 RadPres1 = PressureRadiationLine( iso_sp[ipISO][nelem].trans(ipHi,ipLo), GetDopplerWidth(dense.AtomicWeight[nelem]) );
00302
00303 if( RadPres1 > pmx )
00304 {
00305 ipLineTypePradMax = 2;
00306 pmx = RadPres1;
00307 ip4 = ipISO;
00308 ip3 = nelem;
00309 ip2 = ipHi;
00310 ipLinePradMax = ipLo;
00311 }
00312 Piso_seq[ipISO] += RadPres1;
00313 {
00314
00315 enum {DEBUG_LOC=false};
00316 if( DEBUG_LOC && ipISO==ipH_LIKE && ipLo==3 && ipHi==5 && nzone > 260 )
00317 {
00318 fprintf(ioQQQ,
00319 "DEBUG prad1 \t%.2f\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t\n",
00320 fnzone,
00321 RadPres1,
00322 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc(),
00323 iso_sp[ipISO][nelem].st[ipLo].Pop(),
00324 iso_sp[ipISO][nelem].st[ipHi].Pop(),
00325 iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Pesc());
00326 }
00327 }
00328 }
00329 }
00330 }
00331 }
00332 }
00333 ASSERT( Piso_seq[ipISO] >= 0. );
00334 }
00335 pressure.pres_radiation_lines_curr += Piso_seq[ipISO];
00336 }
00337 {
00338
00339 enum {DEBUG_LOC=false};
00340 # if 0
00341 if( DEBUG_LOC && nzone > 260 )
00342 {
00343 fprintf(ioQQQ,
00344 "DEBUG prad2 \t%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
00345 nzone,
00346 pmx,
00347 iso_sp[ipISO][ip3].trans(ipLinePradMax,ip2).Emis().PopOpc(),
00348 iso_sp[ipISO][ip3].st[0].Pop(),
00349 iso_sp[ipISO][ip3].st[2].Pop(),
00350 iso_sp[ipISO][ip3].st[3].Pop(),
00351 iso_sp[ipISO][ip3].st[4].Pop(),
00352 iso_sp[ipISO][ip3].st[5].Pop(),
00353 iso_sp[ipISO][ip3].st[6].Pop());
00354 }
00355 # endif
00356 if( DEBUG_LOC )
00357 {
00358 fprintf(ioQQQ,
00359 "DEBUG prad3\t%.2e\t%li\t%li\t%li\t%li\t%.2e\t%.2e\t%.2e\n",
00360 pmx,
00361 ip4,
00362 ip3,
00363 ip2,
00364 ipLinePradMax,
00365 iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().PopOpc(),
00366 iso_sp[ip4][ip3].st[ip2].Pop(),
00367 1.-iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().Pesc() );
00368 }
00369 }
00370
00371 if( lgMustEvaluateTrivial || PLevel1 > TrivialLineRadiationPressure )
00372 {
00373 PLevel1 = 0.;
00374
00375 for( long i=1; i <= nLevel1; i++ )
00376 {
00377 if( (*TauLines[i].Hi()).Pop() > 1e-30 )
00378 {
00379 RadPres1 = PressureRadiationLine( TauLines[i], GetDopplerWidth(dense.AtomicWeight[(*TauLines[i].Hi()).nelem()-1]) );
00380
00381 if( RadPres1 > pmx )
00382 {
00383 ipLineTypePradMax = 3;
00384 pmx = RadPres1;
00385 ipLinePradMax = i;
00386 }
00387 PLevel1 += RadPres1;
00388 }
00389 }
00390 ASSERT( PLevel1 >= 0. );
00391 }
00392 pressure.pres_radiation_lines_curr += PLevel1;
00393
00394 if( lgMustEvaluateTrivial || PLevel2 > TrivialLineRadiationPressure )
00395 {
00396
00397 PLevel2 = 0.;
00398 for( long i=0; i < nWindLine; i++ )
00399 {
00400 if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
00401 {
00402 if( (*TauLine2[i].Hi()).Pop() > 1e-30 )
00403 {
00404 RadPres1 = PressureRadiationLine( TauLine2[i], GetDopplerWidth(dense.AtomicWeight[(*TauLine2[i].Hi()).nelem()-1]) );
00405
00406 PLevel2 += RadPres1;
00407 if( RadPres1 > pmx )
00408 {
00409 ipLineTypePradMax = 4;
00410 pmx = RadPres1;
00411 ipLinePradMax = i;
00412 }
00413 }
00414 }
00415 }
00416 ASSERT( PLevel2 >= 0. );
00417 }
00418 pressure.pres_radiation_lines_curr += PLevel2;
00419
00420
00421 if( lgMustEvaluateTrivial || PHfs > TrivialLineRadiationPressure )
00422 {
00423 PHfs = 0.;
00424 for( long i=0; i < nHFLines; i++ )
00425 {
00426 if( (*HFLines[i].Hi()).Pop() > 1e-30 )
00427 {
00428 RadPres1 = PressureRadiationLine( HFLines[i], GetDopplerWidth(dense.AtomicWeight[(*HFLines[i].Hi()).nelem()-1]) );
00429
00430 PHfs += RadPres1;
00431 if( RadPres1 > pmx )
00432 {
00433 ipLineTypePradMax = 8;
00434 pmx = RadPres1;
00435 ipLinePradMax = i;
00436 }
00437 }
00438 }
00439 ASSERT( PHfs >= 0. );
00440 }
00441 pressure.pres_radiation_lines_curr += PHfs;
00442
00443
00444 if( lgMustEvaluateTrivial || P_H2 > TrivialLineRadiationPressure )
00445 {
00446 P_H2 = 0.;
00447 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
00448 {
00449 P_H2 += (*diatom)->H2_RadPress();
00450
00451 if( P_H2 > pmx )
00452 {
00453 pmx = P_H2;
00454 ipLineTypePradMax = 9;
00455 }
00456 ASSERT( P_H2 >= 0. );
00457 }
00458 }
00459 pressure.pres_radiation_lines_curr += P_H2;
00460
00461
00462 if( lgMustEvaluateTrivial || P_FeII > TrivialLineRadiationPressure )
00463 {
00464 P_FeII = FeIIRadPress();
00465 if( P_FeII > pmx )
00466 {
00467 pmx = P_FeII;
00468 ipLineTypePradMax = 7;
00469 }
00470 ASSERT( P_FeII >= 0. );
00471 }
00472 pressure.pres_radiation_lines_curr += P_FeII;
00473
00474
00475 if( lgMustEvaluateTrivial || P_dBase > TrivialLineRadiationPressure )
00476 {
00477 P_dBase = 0.;
00478 for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
00479 {
00480 if( dBaseSpecies[ipSpecies].lgActive )
00481 {
00482 realnum DopplerWidth = GetDopplerWidth( dBaseSpecies[ipSpecies].fmolweight );
00483 for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin();
00484 tr != dBaseTrans[ipSpecies].end(); ++tr)
00485 {
00486 int ipHi = (*tr).ipHi();
00487 if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local)
00488 continue;
00489 int ipLo = (*tr).ipLo();
00490 if( (*tr).ipCont() > 0 && (*(*tr).Hi()).Pop() > 1e-30 )
00491 {
00492 RadPres1 = PressureRadiationLine( *tr, DopplerWidth );
00493
00494 if( RadPres1 > pmx )
00495 {
00496 ipLineTypePradMax = 10;
00497 pmx = RadPres1;
00498 ip3 = ipSpecies;
00499 ip2 = ipHi;
00500 ipLinePradMax = ipLo;
00501 }
00502 P_dBase += RadPres1;
00503 }
00504 }
00505 }
00506 }
00507 ASSERT( P_dBase >= 0. );
00508 }
00509 pressure.pres_radiation_lines_curr += P_dBase;
00510
00511 }
00512 else if( !pressure.lgLineRadPresOn || !rfield.lgDoLineTrans )
00513 {
00514
00515 ipLinePradMax = -1;
00516 ipLineTypePradMax = 0;
00517 }
00518
00519 fixit();
00520
00521
00522
00523 pressure.pbeta = (realnum)(pressure.pres_radiation_lines_curr/SDIV(pressure.PresTotlCurr));
00524
00525
00526 if( pressure.pres_radiation_lines_curr < 0. )
00527 {
00528 fprintf(ioQQQ,
00529 "PresTotCurrent: negative pressure, constituents follow %e %e %e %e %e \n",
00530 Piso_seq[ipH_LIKE],
00531 Piso_seq[ipHE_LIKE],
00532 PLevel1,
00533 PLevel2,
00534 PCO);
00535 ShowMe();
00536 cdEXIT(EXIT_FAILURE);
00537 }
00538
00539
00540
00541 if( trace.lgTrace && ipLineTypePradMax <0 )
00542 {
00543 fprintf(ioQQQ,
00544 " PresTotCurrent, pressure.pbeta = %e, ipLineTypePradMax%li ipLinePradMax=%li \n",
00545 pressure.pbeta,ipLineTypePradMax,ipLinePradMax );
00546 }
00547
00548
00549
00550
00551 phycon.EnergyExcitation = 0.;
00552 #if 0
00553 fixit();
00554
00555 broken();
00556
00557
00558
00559
00560 ipLo = 0;
00561 for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00562 {
00563 for( long nelem=ipISO; nelem < LIMELM; nelem++ )
00564 {
00565 if( dense.IonHigh[nelem] == nelem + 1 - ipISO )
00566 {
00567
00568 for( long ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
00569 {
00570 phycon.EnergyExcitation +=
00571 iso_sp[ipISO][nelem].st[ipHi].Pop() *
00572 iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyErg() *
00573
00574 dynamics.lgISO[ipISO];
00575 }
00576 }
00577 }
00578 }
00579
00580 if( dynamics.lgISO[ipH_LIKE] )
00581
00582 phycon.EnergyExcitation += H2_InterEnergy();
00583
00584
00585 if( dynamics.lgMETALS )
00586 {
00587 for( long i=1; i <= nLevel1; i++ )
00588 {
00589 phycon.EnergyExcitation +=
00590 (*TauLines[i].Hi()).Pop()* TauLines[i].EnergyErg;
00591 }
00592 for( long i=0; i < nWindLine; i++ )
00593 {
00594 if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
00595 {
00596 phycon.EnergyExcitation +=
00597 (*TauLine2[i].Hi()).Pop()* TauLine2[i].EnergyErg;
00598 }
00599 }
00600 for( long i=0; i < nHFLines; i++ )
00601 {
00602 phycon.EnergyExcitation +=
00603 (*HFLines[i].Hi()).Pop()* HFLines[i].EnergyErg;
00604 }
00605
00606
00607 phycon.EnergyExcitation += FeII_InterEnergy();
00608 }
00609 #endif
00610
00611
00612
00613
00614
00615 Magnetic_evaluate();
00616
00617
00618 if( trace.lgTrace && (pressure.PresTotlCurr > SMALLFLOAT) )
00619 {
00620 fprintf(ioQQQ,
00621 " PresTotCurrent zn %.2f Ptot:%.5e Pgas:%.3e Prad:%.3e Pmag:%.3e Pram:%.3e "
00622 "gas parts; H:%.2e He:%.2e L1:%.2e L2:%.2e CO:%.2e fs%.2e h2%.2e turb%.2e\n",
00623 fnzone,
00624 pressure.PresTotlCurr,
00625 pressure.PresGasCurr/pressure.PresTotlCurr,
00626 pressure.pres_radiation_lines_curr*pressure.lgPres_radiation_ON/pressure.PresTotlCurr,
00627 magnetic.pressure/pressure.PresTotlCurr,
00628 pressure.PresRamCurr*pressure.lgPres_ram_ON/pressure.PresTotlCurr,
00629 Piso_seq[ipH_LIKE]/pressure.PresTotlCurr,
00630 Piso_seq[ipHE_LIKE]/pressure.PresTotlCurr,
00631 PLevel1/pressure.PresTotlCurr,
00632 PLevel2/pressure.PresTotlCurr,
00633 PCO/pressure.PresTotlCurr,
00634 PHfs/pressure.PresTotlCurr,
00635 P_H2/pressure.PresTotlCurr,
00636 pressure.PresTurbCurr*DoppVel.lgTurb_pressure/pressure.PresTotlCurr);
00637
00638 }
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648 if( dynamics.lgTimeDependentStatic )
00649 {
00650
00651
00652
00653
00654
00655 phycon.EnthalpyDensity =
00656 0.5*POW2(wind.windv)*dense.xMassDensity +
00657 3./2.*pressure.PresGasCurr +
00658 magnetic.EnthalpyDensity;
00659
00660 }
00661 else
00662 {
00663
00664
00665
00666
00667 phycon.EnthalpyDensity =
00668 0.5*POW2(wind.windv)*dense.xMassDensity +
00669 5./2.*pressure.PresGasCurr +
00670 magnetic.EnthalpyDensity;
00671
00672 }
00673
00674
00675
00676
00677
00678 if( iteration <= 1 && pressure.pres_radiation_lines_curr >= pressure.PresGasCurr )
00679 {
00680
00681 pressure.pres_radiation_lines_curr =
00682 MIN2(pressure.pres_radiation_lines_curr,pressure.PresGasCurr);
00683 pressure.lgPradCap = true;
00684 }
00685
00686
00687
00688 if( pressure.pbeta > pressure.RadBetaMax && nzone )
00689 {
00690 pressure.RadBetaMax = pressure.pbeta;
00691 pressure.ipPradMax_line = ipLinePradMax;
00692 pressure.ipPradMax_nzone = nzone;
00693 if( ipLineTypePradMax == 2 )
00694 {
00695
00696 strcpy( pressure.chLineRadPres, "ISO ");
00697 ASSERT( ip4 < NISO && ip3<LIMELM );
00698 ASSERT( ipLinePradMax>=0 && ip2>=0 && ip3>=0 && ip4>=0 );
00699 strcat( pressure.chLineRadPres, chLineLbl(iso_sp[ip4][ip3].trans(ip2,ipLinePradMax)) );
00700 {
00701
00702 enum {DEBUG_LOC=false};
00703
00704
00705 if( DEBUG_LOC )
00706 {
00707 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",
00708 iteration, nzone,
00709 ip4,ip3,ip2,ipLinePradMax,
00710 iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().TauIn(),
00711 iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().TauTot(),
00712 iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().Pesc(),
00713 iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().Pelec_esc(),
00714 iso_sp[ip4][ip3].trans(ip2,ipLinePradMax).Emis().Pdest());
00715 if( ip2==5 && ipLinePradMax==4 )
00716 {
00717 double width;
00718 fprintf(ioQQQ,"hit it\n");
00719 width = RT_LineWidth(iso_sp[ip4][ip3].trans(ip2,ipLinePradMax),GetDopplerWidth(dense.AtomicWeight[ip3]));
00720 fprintf(ioQQQ,"DEBUG width %e\n", width);
00721 }
00722 }
00723 }
00724 }
00725 else if( ipLineTypePradMax == 3 )
00726 {
00727
00728 ASSERT( ipLinePradMax>=0 );
00729 strcpy( pressure.chLineRadPres, "Level1 ");
00730 strcat( pressure.chLineRadPres, chLineLbl(TauLines[ipLinePradMax]) );
00731 }
00732 else if( ipLineTypePradMax == 4 )
00733 {
00734
00735 ASSERT( ipLinePradMax>=0 );
00736 strcpy( pressure.chLineRadPres, "Level2 ");
00737 strcat( pressure.chLineRadPres, chLineLbl(TauLine2[ipLinePradMax]) );
00738 }
00739 else if( ipLineTypePradMax == 5 )
00740 {
00741 cdEXIT( EXIT_FAILURE );
00742 }
00743 else if( ipLineTypePradMax == 6 )
00744 {
00745 cdEXIT( EXIT_FAILURE );
00746 }
00747 else if( ipLineTypePradMax == 7 )
00748 {
00749
00750 strcpy( pressure.chLineRadPres, "Fe II ");
00751 }
00752 else if( ipLineTypePradMax == 8 )
00753 {
00754
00755 strcpy( pressure.chLineRadPres, "hyperf ");
00756 ASSERT( ipLinePradMax>=0 );
00757 strcat( pressure.chLineRadPres, chLineLbl(HFLines[ipLinePradMax]) );
00758 }
00759 else if( ipLineTypePradMax == 9 )
00760 {
00761
00762 strcpy( pressure.chLineRadPres, "large H2 ");
00763 }
00764 else if( ipLineTypePradMax == 10 )
00765 {
00766
00767 strcpy( pressure.chLineRadPres, "dBaseLin " );
00768 strcat( pressure.chLineRadPres, dBaseSpecies[ip3].chLabel );
00769 }
00770 else
00771 {
00772 fprintf(ioQQQ," PresTotCurrent ipLineTypePradMax set to %li, this is impossible.\n", ipLineTypePradMax);
00773 ShowMe();
00774 cdEXIT(EXIT_FAILURE);
00775 }
00776 }
00777
00778 if( trace.lgTrace && pressure.pbeta > 0.5 )
00779 {
00780 fprintf(ioQQQ,
00781 " PresTotCurrent Pbeta:%.2f due to %s\n",
00782 pressure.pbeta ,
00783 pressure.chLineRadPres
00784 );
00785 }
00786
00787
00788
00789 TotalPressure_v = pressure.PresGasCurr;
00790
00791
00792
00793 TotalPressure_v += pressure.PresRamCurr * pressure.lgPres_ram_ON;
00794
00795
00796
00799 TotalPressure_v += magnetic.pressure * pressure.lgPres_magnetic_ON;
00800
00801
00802
00803 TotalPressure_v += pressure.PresTurbCurr * DoppVel.lgTurb_pressure;
00804
00805
00806
00807
00808 TotalPressure_v += pressure.pres_radiation_lines_curr * pressure.lgPres_radiation_ON;
00809
00810 {
00811 enum{DEBUG_LOC=false};
00812 if( DEBUG_LOC && pressure.PresTotlCurr > SMALLFLOAT )
00813 {
00814 fprintf(ioQQQ,"pressureee%li\t%.4e\t%.4e\t%.4e\t%.3f\t%.3f\t%.3f\n",
00815 nzone,
00816 pressure.PresTotlError,
00817 pressure.PresTotlCurr,
00818 TotalPressure_v ,
00819 pressure.PresGasCurr/pressure.PresTotlCurr,
00820 pressure.pres_radiation_lines_curr/pressure.PresTotlCurr ,
00821 pressure.PresRamCurr/pressure.PresTotlCurr
00822 );
00823 }
00824 }
00825
00826 if( TotalPressure_v< 0. )
00827 {
00828 ASSERT( magnetic.pressure < 0. );
00829
00830
00831 fprintf(ioQQQ," The negative pressure due to ordered magnetic field overwhelms the total outward pressure - please reconsider the geometry & field.\n");
00832 cdEXIT(EXIT_FAILURE);
00833 }
00834
00835 ASSERT( TotalPressure_v > 0. );
00836
00837
00838 pressure.PresMax = MAX2(pressure.PresMax,(realnum)TotalPressure_v);
00839
00840
00841 pressure.PresTotlCurr = TotalPressure_v;
00842
00843 return;
00844 }