00001
00002
00003
00004 #include "cddefines.h"
00005 #include "transition.h"
00006 #include "version.h"
00007 #include "dense.h"
00008 #include "elementnames.h"
00009 #include "lines.h"
00010 #include "lines_service.h"
00011 #include "opacity.h"
00012 #include "phycon.h"
00013 #include "radius.h"
00014 #include "rfield.h"
00015 #include "rt.h"
00016 #include "taulines.h"
00017 #include "conv.h"
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 void TransitionProxy::outline_resonance( ) const
00038 {
00039 bool lgDoChecks = true;
00040 outline(Emis().ColOvTot(), lgDoChecks);
00041 }
00042
00043
00044 void TransitionProxy::outline( double nonScatteredFraction,
00045 bool lgDoChecks ) const
00046 {
00047 long int ip = ipCont()-1;
00048
00049 DEBUG_ENTRY( "TransitionProxy::outline()" );
00050
00051 if ( 0 && lgDoChecks)
00052 {
00053
00054 ASSERT( Emis().phots() >= 0. );
00055 if( Emis().phots() < SMALLFLOAT || EnergyErg() / EN1RYD <= rfield.plsfrq )
00056 return;
00057
00058 ASSERT( Emis().FracInwd() >= 0. );
00059 ASSERT( radius.BeamInIn >= 0. );
00060 ASSERT( radius.BeamInOut >= 0. );
00061 double Ptot = Emis().Pesc() + Emis().Pelec_esc();
00062 double PhotEmit = Emis().Aul()*Ptot*(*Hi()).Pop();
00063
00064 double error = MAX2( SMALLFLOAT*1e3 , 3e-1*PhotEmit );
00065
00066 bool lgGoodSolution = conv.lgConvEden && conv.lgConvIoniz() &&
00067 conv.lgConvPops && conv.lgConvPres && conv.lgConvTemp;
00068
00069
00070
00071
00072
00073
00074 ASSERT( t_version::Inst().lgRelease || !lgGoodSolution ||
00075 Ptot < 1e-3 || fp_equal_tol( Emis().phots(), PhotEmit , error ));
00076 }
00077
00078 bool lgTransStackLine = true;
00079 outline_base(Emis().dampXvel(), Emis().damp(), lgTransStackLine, ip, Emis().phots(), Emis().FracInwd(),
00080 nonScatteredFraction);
00081 }
00082
00083
00084 double emit_frac(const TransitionProxy& t)
00085 {
00086 DEBUG_ENTRY( "emit_frac()" );
00087
00088 ASSERT( t.ipCont() > 0 );
00089
00090
00091
00092 double deexcit_loss = t.Coll().col_str() * dense.cdsqte + t.Emis().Aul()*t.Emis().Pdest();
00093
00094 double rad_deexcit = t.Emis().Aul()*(t.Emis().Pelec_esc() + t.Emis().Pesc());
00095 return rad_deexcit/(deexcit_loss + rad_deexcit);
00096 }
00097
00098
00099
00100 void DumpLine(const TransitionProxy& t)
00101 {
00102 char chLbl[110];
00103
00104 DEBUG_ENTRY( "DumpLine()" );
00105
00106 ASSERT( t.ipCont() > 0 );
00107
00108
00109 strcpy( chLbl, "DEBUG ");
00110 strcat( chLbl, chLineLbl(t) );
00111
00112 fprintf( ioQQQ,
00113 "%10.10s Te%.2e eden%.1e CS%.2e Aul%.1e Tex%.2e cool%.1e het%.1e conopc%.1e albdo%.2e\n",
00114 chLbl,
00115 phycon.te,
00116 dense.eden,
00117 t.Coll().col_str(),
00118 t.Emis().Aul(),
00119 TexcLine(t),
00120 t.Coll().cool(),
00121 t.Coll().heat() ,
00122 opac.opacity_abs[t.ipCont()-1],
00123 opac.albedo[t.ipCont()-1]);
00124
00125 fprintf( ioQQQ,
00126 "Tin%.1e Tout%.1e Esc%.1e eEsc%.1e DesP%.1e Pump%.1e OTS%.1e PopL,U %.1e %.1e PopOpc%.1e\n",
00127 t.Emis().TauIn(),
00128 t.Emis().TauTot(),
00129 t.Emis().Pesc(),
00130 t.Emis().Pelec_esc(),
00131 t.Emis().Pdest(),
00132 t.Emis().pump(),
00133 t.Emis().ots(),
00134 (*t.Lo()).Pop(),
00135 (*t.Hi()).Pop() ,
00136 t.Emis().PopOpc() );
00137 return;
00138 }
00139
00140
00141
00142 double OccupationNumberLine(const TransitionProxy& t)
00143 {
00144 double OccupationNumberLine_v;
00145
00146 DEBUG_ENTRY( "OccupationNumberLine()" );
00147
00148 ASSERT( t.ipCont() > 0 );
00149
00150
00151
00152 if( fabs(t.Emis().PopOpc()) > SMALLFLOAT )
00153 {
00154
00155 OccupationNumberLine_v = ( (*t.Hi()).Pop() / (*t.Hi()).g() ) /
00156 ( t.Emis().PopOpc() / (*t.Lo()).g() ) *
00157 (1. - t.Emis().Pesc());
00158 }
00159 else
00160 {
00161 OccupationNumberLine_v = 0.;
00162 }
00163
00164
00165 return( OccupationNumberLine_v );
00166 }
00167
00168
00169 double TexcLine(const TransitionProxy& t)
00170 {
00171 double TexcLine_v;
00172
00173 DEBUG_ENTRY( "TexcLine()" );
00174
00175
00176
00177 if( (*t.Hi()).Pop() * (*t.Lo()).Pop() > 0. )
00178 {
00179 TexcLine_v = ( (*t.Hi()).Pop() / (*t.Hi()).g() )/( (*t.Lo()).Pop() / (*t.Lo()).g() );
00180 TexcLine_v = log(TexcLine_v);
00181
00182 if( fabs(TexcLine_v) > SMALLFLOAT )
00183 {
00184 TexcLine_v = - t.EnergyK() / TexcLine_v;
00185 }
00186 }
00187 else
00188 {
00189 TexcLine_v = 0.;
00190 }
00191 return( TexcLine_v );
00192 }
00193
00194
00195 void chIonLbl(char *chIonLbl_v, const TransitionProxy& t)
00196 {
00197 DEBUG_ENTRY( "chIonLbl()" );
00198
00199
00200
00201
00202
00203 if( (*t.Hi()).nelem() < 0 )
00204 {
00205
00206 if( (*t.Hi()).chLabel()[0]=='\0' )
00207 strcpy( chIonLbl_v, "Dumy" );
00208 else
00209 strcpy( chIonLbl_v, (*t.Hi()).chLabel() );
00210 }
00211 else
00212 {
00213 chIonLbl( chIonLbl_v, (*t.Hi()).nelem(), (*t.Hi()).IonStg() );
00214 }
00215
00216 return;
00217 }
00218
00219 void chIonLbl(char *chIonLbl_v, const long& nelem, const long& IonStg)
00220 {
00221 DEBUG_ENTRY( "chIonLbl()" );
00222
00223
00224 ASSERT( nelem >= 1 );
00225 ASSERT( nelem <= LIMELM );
00226
00227
00228 strcpy( chIonLbl_v , elementnames.chElementSym[nelem-1] );
00229
00230 strcat( chIonLbl_v, elementnames.chIonStage[IonStg-1] );
00231 return;
00232 }
00233
00234
00235
00236
00237 char* chLineLbl(const TransitionProxy& t)
00238 {
00239 static char chLineLbl_v[11];
00240 static char chSpecies[5];
00241
00242 DEBUG_ENTRY( "chLineLbl()" );
00243
00244 if( (*t.Hi()).nelem() < 1 && (*t.Hi()).IonStg() < 1 )
00245 {
00246 sprintf( chSpecies, "%4.4s", (*t.Hi()).chLabel() );
00247 }
00248 else
00249 {
00250 ASSERT( (*t.Hi()).nelem() >= 1 );
00251 ASSERT( (*t.Hi()).IonStg() >= 1 && (*t.Hi()).IonStg() <= (*t.Hi()).nelem() + 1 );
00252 sprintf( chSpecies, "%2.2s%2.2s",
00253 elementnames.chElementSym[(*t.Hi()).nelem() -1],
00254 elementnames.chIonStage[(*t.Hi()).IonStg()-1] );
00255 }
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265 if( t.WLAng() > (realnum)INT_MAX )
00266 {
00267 sprintf( chLineLbl_v, "%4.4s%5i%c", chSpecies,
00268 (int)(t.WLAng()/1e8), 'c' );
00269 }
00270 else if( t.WLAng() > 9999999. )
00271 {
00272
00273 sprintf( chLineLbl_v, "%4.4s%5.2f%c", chSpecies,
00274 t.WLAng()/1e8, 'c' );
00275 }
00276 else if( t.WLAng() > 999999. )
00277 {
00278
00279 sprintf( chLineLbl_v, "%4.4s%5i%c", chSpecies,
00280 (int)(t.WLAng()/1e4), 'm' );
00281 }
00282 else if( t.WLAng() > 99999. )
00283 {
00284
00285 sprintf( chLineLbl_v, "%4.4s%5.1f%c", chSpecies,
00286 t.WLAng()/1e4, 'm' );
00287 }
00288 else if( t.WLAng() > 9999. )
00289 {
00290 sprintf( chLineLbl_v, "%4.4s%5.2f%c", chSpecies,
00291 t.WLAng()/1e4, 'm' );
00292 }
00293 else if( t.WLAng() >= 100. )
00294 {
00295 sprintf( chLineLbl_v, "%4.4s%5i%c", chSpecies,
00296 (int)t.WLAng(), 'A' );
00297 }
00298
00299
00300 else if( t.WLAng() >= 10. )
00301 {
00302 sprintf( chLineLbl_v, "%4.4s%5.1f%c", chSpecies,
00303 t.WLAng(), 'A' );
00304 }
00305 else
00306 {
00307 sprintf( chLineLbl_v, "%4.4s%5.2f%c", chSpecies,
00308 t.WLAng(), 'A' );
00309 }
00310
00311
00312 ASSERT( chLineLbl_v[10]=='\0' );
00313 return( chLineLbl_v );
00314 }
00315
00316
00317 void PutCS(double cs,
00318 const TransitionProxy& t)
00319 {
00320 DEBUG_ENTRY( "PutCS()" );
00321
00322
00323 ASSERT( cs > 0. );
00324
00325 t.Coll().col_str() = (realnum)cs;
00326
00327 return;
00328 }
00329
00330 void GenerateTransitionConfiguration( const TransitionProxy &t, char *chComment )
00331 {
00332 strcpy( chComment, (*t.Lo()).chConfig() );
00333 strcat( chComment, " - " );
00334 strcat( chComment, (*t.Hi()).chConfig() );
00335 return;
00336 }
00337
00338 static realnum ExtraInten;
00339
00340
00341 STATIC void PutLine_base(const TransitionProxy& t, const char *chComment, const char *chLabelTemp, bool lgLabel)
00342 {
00343 DEBUG_ENTRY( "PutLine_base()" );
00344
00345 char chLabel[5];
00346 double xIntensity,
00347 other,
00348 xIntensity_in;
00349
00350
00351
00352 ASSERT( t.ipCont() > 0. );
00353
00354 if (lgLabel == true)
00355 {
00356 strncpy( chLabel, chLabelTemp, 4 );
00357 chLabel[4] = '\0';
00358 }
00359
00360
00361
00362 if( LineSave.ipass == 0 )
00363 {
00364 if (lgLabel == false)
00365 {
00366
00367
00368 chIonLbl(chLabel, t);
00369 }
00370 xIntensity = 0.;
00371 }
00372 else
00373 {
00374
00375
00376 chLabel[0] = '\0';
00377
00378
00379 xIntensity = t.Emis().xIntensity() + ExtraInten;
00380 }
00381
00382
00383
00384
00385
00386
00387
00388 ExtraInten = 0.;
00389
00390
00391 rt.fracin = t.Emis().FracInwd();
00392 lindst(xIntensity,
00393 t.WLAng(),
00394 chLabel,
00395 t.ipCont(),
00396
00397 'i',
00398
00399 false,
00400 chComment);
00401 rt.fracin = 0.5;
00402
00403
00404
00405 xIntensity_in = xIntensity*t.Emis().FracInwd();
00406 ASSERT( xIntensity_in>=0. );
00407 linadd(xIntensity_in,t.WLAng(),"Inwd",'i',chComment);
00408
00409
00410 other = t.Coll().cool();
00411 linadd(other,t.WLAng(),"Coll",'i',chComment);
00412
00413
00414 double radiative_branching;
00415 enum { lgNEW = true };
00416 if (lgNEW)
00417 {
00418
00419 const double AulEscp = t.Emis().Aul()*(t.Emis().Pesc()+t.Emis().Pelec_esc());
00420
00421
00422 const double sinkrate = AulEscp + t.Emis().Aul()*t.Emis().Pdest() + t.Coll().ColUL( colliders );
00423 if (sinkrate > 0.0)
00424 {
00425 radiative_branching = AulEscp/sinkrate;
00426 }
00427 else
00428 {
00429 radiative_branching = 0.;
00430 }
00431 }
00432 else
00433 {
00434
00435
00436 radiative_branching = (1.-t.Emis().ColOvTot());
00437 }
00438
00439 other = (*t.Lo()).Pop() * t.Emis().pump() * radiative_branching * t.EnergyErg();
00440 linadd(other,t.WLAng(),"Pump",'i',chComment);
00441
00442
00443
00444 other = t.Coll().heat();
00445 linadd(other,t.WLAng(),"Heat",'i',chComment);
00446 }
00447
00448
00449 void PutLine(const TransitionProxy& t, const char *chComment, const char *chLabelTemp)
00450 {
00451 const bool lgLabel = true;
00452 DEBUG_ENTRY( "PutLine()" );
00453 PutLine_base(t, chComment, chLabelTemp, lgLabel);
00454 }
00455
00456
00457 void PutLine(const TransitionProxy& t,
00458 const char *chComment)
00459 {
00460 const bool lgLabel = false;
00461 const char *chLabelTemp = NULL;
00462 DEBUG_ENTRY( "PutLine()" );
00463 PutLine_base(t, chComment, chLabelTemp, lgLabel);
00464 }
00465
00466
00467
00468 void PutExtra(double Extra)
00469 {
00470
00471 DEBUG_ENTRY( "PutExtra()" );
00472
00473 ExtraInten = (realnum)Extra;
00474 return;
00475 }
00476
00477 void TransitionProxy::Junk( void ) const
00478 {
00479
00480 DEBUG_ENTRY( "TransitionProxy::Junk()" );
00481
00482
00483 WLAng() = -FLT_MAX;
00484
00485
00486 EnergyWN() = -FLT_MAX;
00487
00488
00489
00490 ipCont() = -10000;
00491
00492 CollisionJunk( Coll() );
00493
00494
00495
00496 ipEmis() = -1;
00497
00498 setLo(-1);
00499 setHi(-1);
00500 return;
00501 }
00502
00503
00504
00505 void TransitionProxy::Zero( void ) const
00506 {
00507
00508 DEBUG_ENTRY( "TransitionProxy::Zero()" );
00509
00510 CollisionZero( Coll() );
00511
00512 ::Zero( *Lo() );
00513 ::Zero( *Hi() );
00514 EmLineZero( Emis() );
00515 TauZero( Emis() );
00516
00517 return;
00518 }
00519
00520
00521 void LineConvRate2CS( const TransitionProxy& t , realnum rate )
00522 {
00523
00524 DEBUG_ENTRY( "LineConvRate2CS()" );
00525
00526
00527
00528
00529 t.Coll().col_str() = rate * (*t.Hi()).g() / (realnum)dense.cdsqte;
00530
00531
00532
00533 ASSERT( t.Coll().col_str() >= 0. );
00534 return;
00535 }
00536
00537
00538 STATIC void gbar0(double ex,
00539 realnum *g)
00540 {
00541 double a,
00542 b,
00543 c,
00544 d,
00545 y;
00546
00547 DEBUG_ENTRY( "gbar0()" );
00548
00549
00550
00551
00552
00553
00554
00555
00556
00557
00558
00559
00560
00561
00562 y = ex/phycon.te;
00563 if( y < 0.01 )
00564 {
00565 *g = (realnum)(0.29*(log(1.0+1.0/y) - 0.4/POW2(y + 1.0))/exp(y));
00566 }
00567 else
00568 {
00569 if( y > 10.0 )
00570 {
00571 *g = (realnum)(0.066/sqrt(y));
00572 }
00573 else
00574 {
00575 a = 1.5819068e-02;
00576 b = 1.3018207e00;
00577 c = 2.6896230e-03;
00578 d = 2.5486007e00;
00579 d = log(y/c)/d;
00580 *g = (realnum)(a + b*exp(-0.5*POW2(d)));
00581 }
00582 }
00583 return;
00584 }
00585
00586
00587 STATIC void gbar1(double ex,
00588 realnum *g)
00589 {
00590 double y;
00591
00592 DEBUG_ENTRY( "gbar1()" );
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607 y = ex/phycon.te;
00608 *g = (realnum)(0.6 + 0.28*(log(1.0+1.0/y) - 0.4/POW2(y + 1.0)));
00609 return;
00610 }
00611
00612
00613 void MakeCS(const TransitionProxy& t)
00614 {
00615 long int ion;
00616 double Abun,
00617 cs;
00618 realnum
00619 gbar;
00620
00621 DEBUG_ENTRY( "MakeCS()" );
00622
00623
00624
00625
00626 ion = (*t.Hi()).IonStg();
00627
00628 Abun = dense.xIonDense[ (*t.Hi()).nelem() -1 ][ ion-1 ];
00629 if( Abun <= 0. )
00630 {
00631 gbar = 1.;
00632 }
00633 else
00634 {
00635
00636 if( ion == 1 )
00637 {
00638
00639 gbar0(t.EnergyK(),&gbar);
00640 }
00641 else
00642 {
00643
00644 gbar1(t.EnergyK(),&gbar);
00645 }
00646 }
00647
00648
00649 cs = gbar*(14.5104/WAVNRYD)*t.Emis().gf()/t.EnergyWN();
00650
00651
00652 t.Coll().col_str() = (realnum)cs;
00653 return;
00654 }
00655
00656 void TransitionProxy::AddLine2Stack() const
00657 {
00658 DEBUG_ENTRY( "AddLine2Stack()" );
00659
00660 ASSERT( lgLinesAdded == false );
00661
00662 size_t newsize = m_list->Emis.size()+1;
00663 m_list->Emis.resize(newsize);
00664 ipEmis() = newsize-1;
00665 this->resetEmis();
00666 }
00667
00668 void TransitionProxy::AddLoState() const
00669 {
00670 DEBUG_ENTRY( "AddLoState()" );
00671
00672 ASSERT( !lgStatesAdded );
00673
00674 m_list->states->resize(m_list->states->size()+1);
00675
00676 setLo(m_list->states->size()-1);
00677 }
00678
00679 void TransitionProxy::AddHiState() const
00680 {
00681 DEBUG_ENTRY( "AddHiState()" );
00682
00683 ASSERT( !lgStatesAdded );
00684
00685 m_list->states->resize(m_list->states->size()+1);
00686
00687 setHi(m_list->states->size()-1);
00688 }