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