00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #include "cddefines.h"
00016 #include "lines_service.h"
00017 #include "dense.h"
00018 #include "geometry.h"
00019 #include "hydrogenic.h"
00020 #include "ipoint.h"
00021 #include "iso.h"
00022 #include "lines.h"
00023 #include "trace.h"
00024 #include "opacity.h"
00025 #include "physconst.h"
00026 #include "radius.h"
00027 #include "rfield.h"
00028 #include "rt.h"
00029 #include "taulines.h"
00030 #include "thirdparty.h"
00031
00032 void LineStackCreate()
00033 {
00034 DEBUG_ENTRY( "LineStackCreate()" );
00035
00036
00037
00038
00039
00040
00041 LineSave.ipass = -1;
00042 lines();
00043 ASSERT( LineSave.nsum > 0 );
00044
00045
00046
00047
00048 if( LineSv != NULL )
00049 free( LineSv );
00050 if( LineSvSortWL != NULL )
00051 free( LineSvSortWL );
00052
00053
00054 LineSv = (LinSv*)MALLOC((unsigned)LineSave.nsum*sizeof( LinSv ) );
00055 LineSvSortWL = (LinSv*)MALLOC((unsigned)LineSave.nsum*sizeof( LinSv ));
00056
00057
00058 LineSave.nsumAllocated = LineSave.nsum;
00059
00060
00061 for( long i=0; i < LineSave.nsum; i++ )
00062 {
00063 for( long j =0; j<4; ++j )
00064 LineSv[i].SumLine[j] = 0;
00065 }
00066
00067
00068
00069
00070
00071 LineSave.ipass = 0;
00072 lines();
00073
00074 ASSERT( LineSave.nsum > 0);
00075
00076 LineSave.ipass = 1;
00077
00078 if( trace.lgTrace )
00079 fprintf( ioQQQ, "%7ld lines printed in main line array\n",
00080 LineSave.nsum );
00081 }
00082
00083
00084 double eina(double gf,
00085 double enercm,
00086 double gup)
00087 {
00088 double eina_v;
00089
00090 DEBUG_ENTRY( "eina()" );
00091
00092
00093
00094
00095
00096 eina_v = (gf/gup)*TRANS_PROB_CONST*POW2(enercm);
00097 return( eina_v );
00098 }
00099
00100
00101 double GetGF(double trans_prob,
00102 double enercm,
00103 double gup)
00104 {
00105 double GetGF_v;
00106
00107 DEBUG_ENTRY( "GetGF()" );
00108
00109 ASSERT( enercm > 0. );
00110 ASSERT( trans_prob > 0. );
00111 ASSERT( gup > 0.);
00112
00113
00114
00115
00116
00117 GetGF_v = trans_prob*gup/TRANS_PROB_CONST/POW2(enercm);
00118 return( GetGF_v );
00119 }
00120
00121
00122 double abscf(double gf,
00123 double enercm,
00124 double gl)
00125 {
00126 double abscf_v;
00127
00128 DEBUG_ENTRY( "abscf()" );
00129
00130 ASSERT(gl > 0. && enercm > 0. && gf > 0. );
00131
00132
00133
00134
00135 abscf_v = 1.4974e-6*(gf/gl)*(1e4/enercm);
00136 return( abscf_v );
00137 }
00138
00139
00140
00141 double RefIndex(double EnergyWN )
00142 {
00143 double RefIndex_v,
00144 WaveMic,
00145 xl,
00146 xn;
00147
00148 DEBUG_ENTRY( "RefIndex()" );
00149
00150 ASSERT( EnergyWN > 0. );
00151
00152
00153 WaveMic = 1.e4/EnergyWN;
00154
00155
00156 if( WaveMic > 0.2 )
00157 {
00158
00159
00160 xl = 1.0/WaveMic/WaveMic;
00161
00162
00163
00164 xn = 255.4/(41. - xl);
00165 xn += 29498.1/(146. - xl);
00166 xn += 64.328;
00167 RefIndex_v = xn/1.e6 + 1.;
00168 }
00169 else
00170 {
00171 RefIndex_v = 1.;
00172 }
00173 ASSERT( RefIndex_v >= 1. );
00174 return( RefIndex_v );
00175 }
00176
00177
00178
00179
00180
00181
00182 realnum WavlenErrorGet( realnum wavelength )
00183 {
00184 double a;
00185 realnum errorwave;
00186
00187 DEBUG_ENTRY( "WavlenErrorGet()" );
00188
00189 ASSERT( LineSave.sig_figs <= 6 );
00190
00191 if( wavelength > 0. )
00192 {
00193
00194 a = log10( wavelength+FLT_EPSILON);
00195 a = floor(a);
00196 }
00197 else
00198 {
00199
00200
00201 a = 0.;
00202 }
00203
00204 errorwave = 5.f * (realnum)pow( 10., a - (double)LineSave.sig_figs );
00205 return errorwave;
00206 }
00207
00208
00209 STATIC void lincom(
00210 double xInten,
00211 realnum wavelength,
00212 const char *chLab,
00213
00214 long int ipnt,
00215 char chInfo,
00216
00217
00218 const char *chComment,
00219
00220 bool lgAdd)
00221 {
00222 DEBUG_ENTRY( "lincom()" );
00223
00224
00225
00226
00227
00228
00229
00230 if( LineSave.ipass > 0 )
00231 {
00232
00233
00234 if (lgAdd || xInten > 0.)
00235 {
00236
00237
00238
00239
00240 LineSv[LineSave.nsum].SumLine[0] += xInten*radius.dVeffAper;
00241
00242
00243 LineSv[LineSave.nsum].emslin[0] = xInten;
00244 }
00245
00246 if (lgAdd)
00247 {
00248 if (wavelength > 0 )
00249 {
00250
00251
00252
00253
00254
00255 LineSv[LineSave.nsum].emslin[1] = LineSv[LineSave.nsum].emslin[0];
00256 LineSv[LineSave.nsum].SumLine[1] = LineSv[LineSave.nsum].SumLine[0];
00257 }
00258 }
00259 else
00260 {
00261 if ( xInten > 0. && ipnt <= rfield.nflux )
00262 {
00263
00264
00265 const double saveemis = emergent_line(
00266 xInten*rt.fracin , xInten*(1.-rt.fracin) , ipnt );
00267 LineSv[LineSave.nsum].emslin[1] = saveemis;
00268 LineSv[LineSave.nsum].SumLine[1] += saveemis*radius.dVeffAper;
00269 }
00270 }
00271 }
00272
00273 else if( LineSave.ipass == 0 )
00274 {
00275
00276
00277 ASSERT( (chInfo == 'c') || (chInfo == 'h') || (chInfo == 'i') || (chInfo == 'r' ) );
00278
00279 LineSv[LineSave.nsum].chSumTyp = (char)chInfo;
00280 LineSv[LineSave.nsum].emslin[0] = 0.;
00281 LineSv[LineSave.nsum].emslin[1] = 0.;
00282 LineSv[LineSave.nsum].chComment = chComment;
00283
00284
00285 ASSERT( strlen( chLab )<5 );
00286 strcpy( LineSv[LineSave.nsum].chALab, chLab );
00287
00288 if (lgAdd)
00289 {
00290 LineSv[LineSave.nsum].wavelength = wavelength;
00291 }
00292 else
00293 {
00294
00295
00296 LineSv[LineSave.nsum].wavelength = fabs(wavelength);
00297 LineSv[LineSave.nsum].SumLine[0] = 0.;
00298 LineSv[LineSave.nsum].SumLine[1] = 0.;
00299
00300
00301
00302
00303
00304 ASSERT( ipnt > 0 );
00305 # ifndef NDEBUG
00306 double error = MAX2(0.1*rfield.AnuOrg[ipnt-1] , rfield.widflx[ipnt-1] );
00307 ASSERT( wavelength<=0 ||
00308 fabs( rfield.AnuOrg[ipnt-1] - RYDLAM / wavelength) < error );
00309 # endif
00310 }
00311 }
00312
00313
00314 ++LineSave.nsum;
00315
00316
00317
00318 }
00319
00320
00321 void linadd(
00322 double xInten,
00323 realnum wavelength,
00324 const char *chLab,
00325 char chInfo,
00326
00327 const char *chComment )
00328 {
00329 DEBUG_ENTRY( "linadd()" );
00330
00331
00332 const long int ipnt = LONG_MAX;
00333
00334 lincom( xInten, wavelength, chLab, ipnt, chInfo, chComment, true );
00335 }
00336
00337
00338
00339
00340 double emergent_line(
00341
00342 double emissivity_in ,
00343
00344 double emissivity_out ,
00345
00346 long int ipCont )
00347 {
00348
00349 double emergent_in , emergent_out;
00350 long int i = ipCont-1;
00351
00352 DEBUG_ENTRY( "emergent_line()" );
00353
00354 ASSERT( i >= 0 && i < rfield.nupper-1 );
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365 if( iteration == 1 )
00366 {
00367
00368
00369 emergent_in = emissivity_in*opac.E2TauAbsFace[i];
00370 emergent_out = emissivity_out;
00371 }
00372 else
00373 {
00374 if( geometry.lgSphere )
00375 {
00376
00377
00378
00379 emergent_in = emissivity_in * opac.E2TauAbsFace[i] *opac.E2TauAbsTotal[i];
00380
00381
00382 emergent_out = emissivity_out * opac.E2TauAbsOut[i];
00383 }
00384 else
00385 {
00386
00387
00388
00389 double reflected = emissivity_out * opac.albedo[i] * (1.-opac.E2TauAbsOut[i]);
00390
00391 emergent_in = (emissivity_in + reflected) * opac.E2TauAbsFace[i];
00392
00393 emergent_out = (emissivity_out - reflected) * opac.E2TauAbsOut[i];
00394 }
00395 }
00396
00397 return( emergent_in + emergent_out );
00398 }
00399
00400
00401 void outline_base(double dampXvel, double damp, bool lgTransStackLine, long int ip, double phots, realnum inwd,
00402 double nonScatteredFraction)
00403 {
00404 DEBUG_ENTRY( "outline_base()" );
00405
00406 static const bool DO_PROFILE = false;
00407
00408 if( !DO_PROFILE || !rfield.lgDoLineTrans )
00409 outline_base_bin(lgTransStackLine, ip, phots, inwd, nonScatteredFraction);
00410 else
00411 {
00412 ASSERT( damp > 0. );
00413 double LineWidth = dampXvel/damp;
00414 LineWidth = MIN2( 0.1 * SPEEDLIGHT, LineWidth );
00415 double sigma = (LineWidth/SPEEDLIGHT);
00416 long ip3SigmaRed = ipoint( MAX2( rfield.emm, rfield.anu[ip] - 3.*sigma*rfield.anu[ip] ) );
00417 long ip3SigmaBlue = ipoint( MIN2( rfield.egamry, rfield.anu[ip] + 3.*sigma*rfield.anu[ip] ) );
00418 ASSERT( ip3SigmaBlue >= ip3SigmaRed );
00419 long numBins = ip3SigmaBlue - ip3SigmaRed + 1;
00420
00421 if( numBins < 3 )
00422 outline_base_bin(lgTransStackLine, ip, phots, inwd, nonScatteredFraction);
00423 else
00424 {
00425 valarray<realnum> x(numBins);
00426 valarray<realnum> profile(numBins);
00427
00428 for( long ipBin=ip3SigmaRed; ipBin<=ip3SigmaBlue; ipBin++ )
00429 x[ipBin] = (rfield.anu[ip] - rfield.anu[ipBin])/rfield.anu[ip]/sigma;
00430
00431 VoigtU(damp,get_ptr(x),get_ptr(profile),numBins);
00432
00433 for( long ipBin=ip3SigmaRed; ipBin<=ip3SigmaBlue; ipBin++ )
00434 outline_base_bin(lgTransStackLine, ipBin, phots*profile[ipBin-ip3SigmaRed], inwd, nonScatteredFraction);
00435 }
00436 }
00437 }
00438
00439
00440 void outline_base_bin(bool lgTransStackLine, long int ip, double phots, realnum inwd,
00441 double nonScatteredFraction)
00442 {
00443 DEBUG_ENTRY( "outline_base_bin()" );
00444
00445 if (lgTransStackLine)
00446 {
00447 rfield.DiffuseLineEmission[ip] +=
00448 (realnum)phots;
00449
00450
00451 rfield.reflin[0][ip] +=
00452 (realnum)(inwd*phots*radius.BeamInIn);
00453
00454
00455 rfield.outlin[0][ip] +=
00456 (realnum)(inwd*phots*radius.BeamInOut*opac.tmn[ip]*nonScatteredFraction);
00457
00458
00459 rfield.outlin[0][ip] +=
00460 (realnum)((1.-inwd)*phots*radius.BeamOutOut*opac.tmn[ip]*nonScatteredFraction);
00461 }
00462 else
00463 {
00464 rfield.reflin[0][ip] +=
00465 (realnum)(phots*radius.dVolReflec);
00466
00467 rfield.outlin[0][ip] +=
00468 (realnum)(phots*radius.dVolOutwrd*opac.ExpZone[ip]);
00469 }
00470 }
00471
00472
00473 void lindst(
00474
00475 double xInten,
00476
00477 realnum wavelength,
00478
00479 const char *chLab,
00480
00481 long int ipnt,
00482
00483 char chInfo,
00484
00485 bool lgOutToo,
00486
00487 const char *chComment )
00488 {
00489 DEBUG_ENTRY( "lindst()" );
00490
00491
00492 ASSERT( !lgOutToo || chInfo!='i' );
00493
00494 lincom(xInten, wavelength, chLab, ipnt, chInfo, chComment, false );
00495
00496 if( LineSave.ipass > 0 )
00497 {
00498
00499
00500 if (lgOutToo && xInten > 0.)
00501 {
00502
00503
00504
00505
00506
00507
00508 const bool lgTransStackLine = false;
00509 const long int ip = ipnt - 1;
00510 const double phots = xInten/(rfield.anu[ipnt-1]*EN1RYD);
00511 const realnum inwd = (realnum)(1.0-(1.+geometry.covrt)/2.);
00512 const double nonScatteredFraction = 1.;
00513
00514 outline_base_bin(lgTransStackLine, ip, phots, inwd, nonScatteredFraction);
00515 }
00516 }
00517 }
00518
00519
00520 void lindst(
00521 double dampXvel,
00522 double damp,
00523
00524 double xInten,
00525
00526 realnum wavelength,
00527
00528 const char *chLab,
00529
00530 long int ipnt,
00531
00532 char chInfo,
00533
00534 bool lgOutToo,
00535
00536 const char *chComment )
00537 {
00538 DEBUG_ENTRY( "lindst()" );
00539
00540
00541 ASSERT( !lgOutToo || chInfo!='i' );
00542
00543 lincom(xInten, wavelength, chLab, ipnt, chInfo, chComment, false );
00544
00545 if( LineSave.ipass > 0 )
00546 {
00547
00548
00549 if (lgOutToo && xInten > 0.)
00550 {
00551
00552
00553
00554
00555
00556
00557 const bool lgTransStackLine = false;
00558 const long int ip = ipnt - 1;
00559 const double phots = xInten/(rfield.anu[ipnt-1]*EN1RYD);
00560 const realnum inwd = (realnum)(1.0-(1.+geometry.covrt)/2.);
00561 const double nonScatteredFraction = 1.;
00562
00563 outline_base(dampXvel, damp, lgTransStackLine, ip, phots, inwd, nonScatteredFraction);
00564 }
00565 }
00566 }
00567
00568
00569 void lindst(
00570 const TransitionProxy& t,
00571
00572 const char *chLab,
00573
00574 char chInfo,
00575
00576 bool lgOutToo,
00577
00578 const char *chComment )
00579 {
00580 DEBUG_ENTRY( "lindst()" );
00581
00582 lindst( t.Emis().dampXvel(), t.Emis().damp(), t.Emis().xIntensity(), t.WLAng(), chLab, t.ipCont(), chInfo,
00583 lgOutToo, chComment );
00584
00585 }
00586
00587
00588 void PntForLine(
00589
00590 double wavelength,
00591
00592 const char *chLabel,
00593
00594
00595 long int *ipnt)
00596 {
00597
00598
00599
00600
00601
00602 const int MAXFORLIN = 1000;
00603 static long int ipForLin[MAXFORLIN]={0};
00604
00605
00606 static long int nForLin;
00607
00608 DEBUG_ENTRY( "PntForLine()" );
00609
00610
00611 ASSERT( wavelength >= 0. );
00612
00613 if( wavelength == 0. )
00614 {
00615
00616 nForLin = 0;
00617 }
00618 else
00619 {
00620
00621 if( LineSave.ipass > 0 )
00622 {
00623
00624 *ipnt = ipForLin[nForLin];
00625 }
00626 else if( LineSave.ipass == 0 )
00627 {
00628
00629 if( nForLin >= MAXFORLIN )
00630 {
00631 fprintf( ioQQQ, "PROBLEM %5ld lines is too many for PntForLine.\n",
00632 nForLin );
00633 fprintf( ioQQQ, " Increase the value of maxForLine everywhere in the code.\n" );
00634 cdEXIT(EXIT_FAILURE);
00635 }
00636
00637
00638 const double EnergyRyd = RYDLAM/wavelength;
00639 ipForLin[nForLin] = ipLineEnergy(EnergyRyd,chLabel , 0);
00640 *ipnt = ipForLin[nForLin];
00641 }
00642 else
00643 {
00644
00645 *ipnt = 0;
00646 }
00647 ++nForLin;
00648 }
00649 return;
00650 }
00651
00652
00653 double ConvRate2CS( realnum gHi , realnum rate )
00654 {
00655
00656 double cs;
00657
00658 DEBUG_ENTRY( "ConvRate2CS()" );
00659
00660
00661
00662
00663 cs = rate * gHi / dense.cdsqte;
00664
00665
00666
00667 ASSERT( cs >= 0. );
00668 return cs;
00669 }
00670
00671
00672 double ConvCrossSect2CollStr( double CrsSectCM2, double gLo, double E_ProjectileRyd, double reduced_mass_grams )
00673 {
00674 double CollisionStrength;
00675
00676 DEBUG_ENTRY( "ConvCrossSect2CollStr()" );
00677
00678 ASSERT( CrsSectCM2 >= 0. );
00679 ASSERT( gLo >= 0. );
00680 ASSERT( E_ProjectileRyd >= 0. );
00681 ASSERT( reduced_mass_grams >= 0. );
00682
00683 CollisionStrength = CrsSectCM2 * gLo * E_ProjectileRyd / (PI*BOHR_RADIUS_CM*BOHR_RADIUS_CM);
00684
00685
00686 #if 0
00687 CollisionStrength *= reduced_mass_grams / ELECTRON_MASS;
00688 #endif
00689
00690 ASSERT( CollisionStrength >= 0. );
00691 return CollisionStrength;
00692 }
00693
00694
00695 double totlin(
00696
00697
00698
00699
00700 int chInfo)
00701 {
00702 long int i;
00703 double totlin_v;
00704
00705 DEBUG_ENTRY( "totlin()" );
00706
00707
00708
00709
00710
00711
00712 if( (chInfo != 'i' && chInfo != 'r') && chInfo != 'c' )
00713 {
00714 fprintf( ioQQQ, " TOTLIN does not understand chInfo=%c\n",
00715 chInfo );
00716 cdEXIT(EXIT_FAILURE);
00717 }
00718
00719
00720
00721 totlin_v = 0.;
00722 for( i=0; i < LineSave.nsum; i++ )
00723 {
00724 if( LineSv[i].chSumTyp == chInfo )
00725 {
00726 totlin_v += LineSv[i].SumLine[0];
00727 }
00728 }
00729 return( totlin_v );
00730 }
00731
00732
00733
00734 const TransitionProxy FndLineHt(long int *level)
00735 {
00736 long int i;
00737 TransitionProxy t;
00738 DEBUG_ENTRY( "FndLineHt()" );
00739
00740 double Strong = -1.;
00741 *level = 0;
00742
00743
00744 for( i=1; i <= nLevel1; i++ )
00745 {
00746
00747 if( TauLines[i].Coll().heat() > Strong )
00748 {
00749 *level = 1;
00750 t = TauLines[i];
00751 Strong = TauLines[i].Coll().heat();
00752 }
00753 }
00754
00755
00756 for( i=0; i < nWindLine; i++ )
00757 {
00758 if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
00759 {
00760
00761 if( TauLine2[i].Coll().heat() > Strong )
00762 {
00763 *level = 2;
00764 t = TauLine2[i];
00765 Strong = TauLine2[i].Coll().heat();
00766 }
00767 }
00768 }
00769
00770
00771 for( i=0; i < nHFLines; i++ )
00772 {
00773
00774 if( HFLines[i].Coll().heat() > Strong )
00775 {
00776 *level = 3;
00777 t = HFLines[i];
00778 Strong = HFLines[i].Coll().heat();
00779 }
00780 }
00781
00782
00783 for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
00784 {
00785 for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
00786 em != dBaseTrans[ipSpecies].Emis().end(); ++em)
00787 {
00788
00789 if( (*em).Tran().Coll().heat() > Strong )
00790 {
00791 *level = 4;
00792 t = (*em).Tran();
00793 Strong = t.Coll().heat();
00794 }
00795 }
00796 }
00797
00798 fixit();
00799
00800
00801 ASSERT( t.associated() );
00802 return t;
00803 }
00804