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