00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 #include "cddefines.h"
00035 #include "physconst.h"
00036 #include "phycon.h"
00037 #include "lines.h"
00038 #include "radius.h"
00039 #include "geometry.h"
00040 #include "elementnames.h"
00041 #include "rt.h"
00042 #include "dense.h"
00043 #include "rfield.h"
00044 #include "opacity.h"
00045 #include "ipoint.h"
00046 #include "iso.h"
00047 #include "taulines.h"
00048 #include "hydrogenic.h"
00049 #include "lines_service.h"
00050
00051
00052 void outline( EmLine *t )
00053 {
00054 long int ip = t->ipCont-1;
00055 double xInWrd = t->phots*t->FracInwd;
00056
00057 DEBUG_ENTRY( "outline()" );
00058
00059
00060 rfield.reflin[ip] += (float)(xInWrd*radius.BeamInIn);
00061
00062
00063 rfield.outlin[ip] += (float)(xInWrd*radius.BeamInOut*opac.tmn[ip]*
00064 t->ColOvTot);
00065
00066
00068 rfield.outlin[ip] += (float)(t->phots*
00069 (1. - t->FracInwd)*radius.BeamOutOut* t->ColOvTot);
00070
00071 DEBUG_EXIT( "outline()" );
00072
00073 return;
00074 }
00075
00076
00077 void HLineTransOpacSet( long int nelem )
00078 {
00079 long int limit , ipHi , ipLo;
00080 double z4 , factor;
00081
00082 DEBUG_ENTRY( "HLineTransOpacSet()" );
00083
00084
00085 z4 = POW2(nelem+1.);
00086 z4 *= z4;
00087
00088
00089
00090
00091
00092
00093
00094
00095 limit = MIN2(16,iso.numLevels_local[ipH_LIKE][nelem]);
00096
00097
00098
00099 for( ipHi=4; ipHi < limit; ipHi++ )
00100 {
00101 for( ipLo=3; ipLo < ipHi; ipLo++ )
00102 {
00103 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul =
00104 (float)(hydro.HyLife[ipHi]*
00105 HydroBranch(ipHi,ipLo,nelem+1)*z4);
00106
00107 ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul > 0.);
00108
00109
00110 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].opacity =
00111 (float)(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul*
00112 2.2448e-26*iso.stat[ipH_LIKE][nelem][ipHi]/
00113 iso.stat[ipH_LIKE][nelem][ipLo]*
00114 POW3(RYDLAM/(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyWN * WAVNRYD)));
00115
00116
00117 ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].opacity > 0.);
00118 }
00119 }
00120
00121
00122
00123
00124
00125 factor = MAX2( 0.25 , 0.32 - 0.07*dense.eden/(dense.eden + 1e7) );
00126
00127
00128 for( ipHi=4; ipHi < limit; ipHi++ )
00129 {
00130
00131 EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].Aul =
00132 (float)(hydro.HyLife[ipHi]*factor*HydroBranch(ipHi,2,nelem+1)*z4);
00133
00134
00135 ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].Aul > 0.);
00136
00137
00138 EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].Aul =
00139 (float)(EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].Aul / factor *( 1. - factor ));
00140
00141
00142 ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].Aul > 0.);
00143
00144
00145 EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].opacity =
00146 (float)(EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].Aul*
00147 2.2448e-26*iso.stat[ipH_LIKE][nelem][ipHi]/
00148 iso.stat[ipH_LIKE][nelem][ipH2s]*
00149 POW3(RYDLAM/(EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].EnergyWN * WAVNRYD)));
00150
00151
00152 ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].opacity > 0.);
00153
00154
00155 EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].opacity =
00156 (float)(EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].Aul*
00157 2.2448e-26*iso.stat[ipH_LIKE][nelem][ipHi]/
00158 iso.stat[ipH_LIKE][nelem][ipH2p]*
00159 POW3(RYDLAM/ (EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].EnergyWN * WAVNRYD)));
00160
00161
00162 ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].opacity > 0.);
00163 }
00164
00165
00166
00167 for( ipHi=limit; ipHi < iso.numLevels_local[ipH_LIKE][nelem]; ipHi++ )
00168 {
00169
00170 EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].Aul =
00171 (float)(hydro.HyLife[ipHi]*factor*HydroBranch(limit-1,2,nelem+1)*z4);
00172
00173
00174 ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].Aul > 0.);
00175
00176
00177 EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].opacity =
00178 (float)(EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].Aul*
00179 2.2448e-26*iso.stat[ipH_LIKE][nelem][ipHi]/
00180 iso.stat[ipH_LIKE][nelem][ipH2s]*
00181 POW3(RYDLAM/(EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].EnergyWN * WAVNRYD)));
00182
00183
00184 EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].Aul =
00185 (float)(EmisLines[ipH_LIKE][nelem][ipHi][ipH2s].Aul / factor *( 1. - factor ));
00186
00187
00188 ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].Aul > 0.);
00189
00190
00191 EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].opacity =
00192 (float)(EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].Aul*
00193 2.2448e-26*iso.stat[ipH_LIKE][nelem][ipHi]/
00194 iso.stat[ipH_LIKE][nelem][ipH2p]*
00195 POW3(RYDLAM/ (EmisLines[ipH_LIKE][nelem][ipHi][ipH2p].EnergyWN * WAVNRYD)));
00196
00197 for( ipLo=3; ipLo < limit - 1; ipLo++ )
00198 {
00199 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul =
00200 (float)(hydro.HyLife[ipHi]*HydroBranch(limit-1,ipLo,nelem+1)*z4);
00201
00202
00203 ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul > 0.);
00204
00205
00206 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].opacity =
00207 (float)(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].Aul*
00208 2.2448e-26*iso.stat[ipH_LIKE][nelem][ipHi]/
00209 iso.stat[ipH_LIKE][nelem][ipLo]*
00210 POW3(RYDLAM/(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].EnergyWN * WAVNRYD)));
00211
00212
00213 ASSERT(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].opacity > 0.);
00214 }
00215 }
00216
00217 DEBUG_EXIT( "HLineTransOpacSet()" );
00218
00219 return;
00220
00221 }
00222
00223
00224 double emit_frac( EmLine *t )
00225 {
00226
00227
00228 double deexcit_loss = t->cs * dense.cdsqte + t->Aul*t->Pdest;
00229
00230 double rad_deexcit = t->Aul*(t->Pelec_esc + t->Pesc);
00231
00232 return rad_deexcit/(deexcit_loss + rad_deexcit);
00233 }
00234
00235
00236
00237 void DumpLine(EmLine * t)
00238 {
00239 char chLbl[11];
00240
00241 DEBUG_ENTRY( "DumpLine()" );
00242
00243
00244 strcpy( chLbl, chLineLbl(t) );
00245
00246 fprintf( ioQQQ,
00247 " %10.10s Te%.2e eden%.1e CS%.2e Aul%.1e Tex%.2e cool%.1e het%.1e conopc%.1e albdo%.2e\n",
00248 chLbl,
00249 phycon.te,
00250 dense.eden,
00251 t->cs,
00252 t->Aul,
00253 TexcLine(t),
00254 t->cool,
00255 t->heat ,
00256 opac.opacity_abs[t->ipCont-1],
00257 opac.albedo[t->ipCont-1]);
00258
00259 fprintf( ioQQQ,
00260 " Tin%.1e Tout%.1e Esc%.1e eEsc%.1e DesP%.1e Pump%.1e OTS%.1e PopL,U %.1e %.1e PopOpc%.1e\n",
00261 t->TauIn,
00262 t->TauTot,
00263 t->Pesc,
00264 t->Pelec_esc,
00265 t->Pdest,
00266 t->pump,
00267 t->ots,
00268 t->PopLo,
00269 t->PopHi ,
00270 t->PopOpc );
00271
00272 DEBUG_EXIT( "DumpLine()" );
00273
00274 return;
00275 }
00276
00277
00278
00279 double OccupationNumberLine( EmLine * t )
00280 {
00281 double OccupationNumberLine_v;
00282
00283 DEBUG_ENTRY( "OccupationNumberLine()" );
00284
00285
00286 if( t->PopLo > SMALLFLOAT )
00287 {
00288
00289 double PopLo_corr = t->PopLo / t->gLo - t->PopHi / t->gHi;
00290 OccupationNumberLine_v = ( t->PopHi / t->gHi )/SDIV( PopLo_corr ) *
00291 (1. - t->Pesc);
00292 }
00293 else
00294 {
00295 OccupationNumberLine_v = 0.;
00296 }
00297
00298 DEBUG_EXIT( "OccupationNumberLine()" );
00299
00300 return( OccupationNumberLine_v );
00301 }
00302
00303
00304 double TexcLine(EmLine * t)
00305 {
00306 double TexcLine_v;
00307
00308 DEBUG_ENTRY( "TexcLine()" );
00309
00310
00311
00312 if( t->PopHi * t->PopLo > 0. )
00313 {
00314 TexcLine_v = ( t->PopHi / t->gHi )/( t->PopLo / t->gLo );
00315 TexcLine_v = log(TexcLine_v);
00316
00317 if( fabs(TexcLine_v) > SMALLFLOAT )
00318 {
00319 TexcLine_v = - t->EnergyK / TexcLine_v;
00320 }
00321 }
00322 else
00323 {
00324 TexcLine_v = 0.;
00325 }
00326
00327 DEBUG_EXIT( "TexcLine()" );
00328 return( TexcLine_v );
00329 }
00330
00331
00332 double eina(double gf,
00333 double enercm,
00334 double gup)
00335 {
00336 double eina_v;
00337
00338 DEBUG_ENTRY( "eina()" );
00339
00340
00341
00342
00343
00344 eina_v = (gf/gup)*TRANS_PROB_CONST*POW2(enercm);
00345
00346 DEBUG_EXIT( "eina()" );
00347
00348 return( eina_v );
00349 }
00350
00351
00352 double GetGF(double trans_prob,
00353 double enercm,
00354 double gup)
00355 {
00356 double GetGF_v;
00357
00358 DEBUG_ENTRY( "GetGF()" );
00359
00360 ASSERT( enercm > 0. );
00361 ASSERT( trans_prob > 0. );
00362 ASSERT( gup > 0.);
00363
00364
00365
00366
00367
00368 GetGF_v = trans_prob*gup/TRANS_PROB_CONST/POW2(enercm);
00369
00370 DEBUG_EXIT( "GetGF()" );
00371
00372 return( GetGF_v );
00373 }
00374
00375
00376 double abscf(double gf,
00377 double enercm,
00378 double gl)
00379 {
00380 double abscf_v;
00381
00382 DEBUG_ENTRY( "abscf()" );
00383
00384 ASSERT(gl > 0. && enercm > 0. && gf > 0. );
00385
00386
00387
00388
00389 abscf_v = 1.4974e-6*(gf/gl)*(1e4/enercm);
00390
00391 DEBUG_EXIT( "abscf()" );
00392
00393 return( abscf_v );
00394 }
00395
00396
00397 void chIonLbl(char *chIonLbl_v , EmLine * t )
00398 {
00399
00400
00401 DEBUG_ENTRY( "chIonLbl()" );
00402
00403
00404
00405
00406
00407 if( t->nelem < 0 )
00408 {
00409
00410 strcpy( chIonLbl_v, "Dumy" );
00411 }
00412 else if( t->nelem-1 >= LIMELM )
00413 {
00414
00415
00416
00417 strcpy( chIonLbl_v , elementnames.chElementNameShort[t->nelem-1] );
00418
00419
00420
00421 }
00422
00423 else
00424 {
00425
00426
00427 strcpy( chIonLbl_v , elementnames.chElementSym[t->nelem -1] );
00428
00429
00430 strcat( chIonLbl_v, elementnames.chIonStage[t->IonStg-1]);
00431 }
00432
00433 DEBUG_EXIT( "chIonLbl()" );
00434
00435
00436 return;
00437 }
00438
00439
00440
00441
00442 char* chLineLbl(EmLine * t )
00443 {
00444 static char chLineLbl_v[11];
00445
00446 DEBUG_ENTRY( "chLineLbl()" );
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457 if( t->WLAng > (float)INT_MAX )
00458 {
00459 sprintf( chLineLbl_v, "%2.2s%2.2s%5i%c",
00460 elementnames.chElementSym[t->nelem -1],
00461 elementnames.chIonStage[t->IonStg-1],
00462 (int)(t->WLAng/1e8), 'c' );
00463 }
00464 else if( t->WLAng > 9999999. )
00465 {
00466
00467 sprintf( chLineLbl_v, "%2.2s%2.2s%5.2f%c",
00468 elementnames.chElementSym[t->nelem -1],
00469 elementnames.chIonStage[t->IonStg-1],
00470 t->WLAng/1e8, 'c' );
00471 }
00472 else if( t->WLAng > 999999. )
00473 {
00474
00475 sprintf( chLineLbl_v, "%2.2s%2.2s%5i%c",
00476 elementnames.chElementSym[t->nelem -1],
00477 elementnames.chIonStage[t->IonStg-1],
00478 (int)(t->WLAng/1e4), 'm' );
00479 }
00480 else if( t->WLAng > 99999. )
00481 {
00482
00483 sprintf( chLineLbl_v, "%2.2s%2.2s%5.1f%c",
00484 elementnames.chElementSym[t->nelem -1],
00485 elementnames.chIonStage[t->IonStg-1],
00486 t->WLAng/1e4, 'm' );
00487 }
00488 else if( t->WLAng > 9999. )
00489 {
00490 sprintf( chLineLbl_v, "%2.2s%2.2s%5.2f%c",
00491 elementnames.chElementSym[ t->nelem -1],
00492 elementnames.chIonStage[t->IonStg-1],
00493 t->WLAng/1e4, 'm' );
00494 }
00495 else if( t->WLAng >= 100. )
00496 {
00497 sprintf( chLineLbl_v, "%2.2s%2.2s%5i%c",
00498 elementnames.chElementSym[ t->nelem -1],
00499 elementnames.chIonStage[t->IonStg-1],
00500 (int)t->WLAng, 'A' );
00501 }
00502
00503
00504 else if( t->WLAng >= 10. )
00505 {
00506 sprintf( chLineLbl_v, "%2.2s%2.2s%5.1f%c",
00507 elementnames.chElementSym[ t->nelem -1],
00508 elementnames.chIonStage[t->IonStg-1],
00509 t->WLAng, 'A' );
00510 }
00511 else
00512 {
00513 sprintf( chLineLbl_v, "%2.2s%2.2s%5.2f%c",
00514 elementnames.chElementSym[ t->nelem -1],
00515 elementnames.chIonStage[t->IonStg-1],
00516 t->WLAng, 'A' );
00517 }
00518
00519
00520 ASSERT( chLineLbl_v[10]==0 );
00521
00522 DEBUG_EXIT( "chLineLbl()" );
00523
00524 return( chLineLbl_v );
00525 }
00526
00527
00528
00529 double RefIndex(double EnergyWN )
00530 {
00531 double RefIndex_v,
00532 WaveMic,
00533 xl,
00534 xn;
00535
00536 DEBUG_ENTRY( "RefIndex()" );
00537
00538 ASSERT( EnergyWN > 0. );
00539
00540
00541 WaveMic = 1.e4/EnergyWN;
00542
00543
00544 if( WaveMic > 0.2 )
00545 {
00546
00547
00548 xl = 1.0/WaveMic/WaveMic;
00549
00550
00551
00552 xn = 255.4/(41. - xl);
00553 xn += 29498.1/(146. - xl);
00554 xn += 64.328;
00555 RefIndex_v = xn/1.e6 + 1.;
00556 }
00557 else
00558 {
00559 RefIndex_v = 1.;
00560 }
00561
00562 DEBUG_EXIT( "RefIndex()" );
00563
00564 return( RefIndex_v );
00565 }
00566
00567
00568 void PutCS(double cs,
00569
00570 EmLine * t)
00571 {
00572
00573 DEBUG_ENTRY( "PutCS()" );
00574
00575
00576
00577 ASSERT( cs >= 0. );
00578
00579 t->cs = (float)cs;
00580
00581 DEBUG_EXIT( "PutCS()" );
00582
00583 return;
00584 }
00585
00586
00587
00588
00589
00590
00591 float WavlenErrorGet( float wavelength )
00592 {
00593 double a;
00594 float errorwave;
00595
00596 DEBUG_ENTRY( "WavlenErrorGet()" );
00597
00598 ASSERT( LineSave.sig_figs <= 6 );
00599
00600 if( wavelength > 0. )
00601 {
00602
00603 a = log10( wavelength+FLT_EPSILON);
00604 a = floor(a);
00605 }
00606 else
00607 {
00608
00609
00610 a = 0.;
00611 }
00612
00613 errorwave = 5.f * (float)pow( 10., a - (double)LineSave.sig_figs );
00614
00615 DEBUG_EXIT( "WavlenErrorGet()" );
00616 return errorwave;
00617 }
00618
00619
00620 void linadd(
00621 double xInten,
00622 float wavelength,
00623 const char *chLab,
00624 char chInfo
00625 )
00626 {
00627
00628 DEBUG_ENTRY( "linadd()" );
00629
00630
00631
00632
00633
00634
00635
00636 if( LineSave.ipass > 0 )
00637 {
00638
00639
00640 LineSv[LineSave.nsum].sumlin[0] += xInten*radius.dVeff;
00641
00642 LineSv[LineSave.nsum].emslin[0] = xInten;
00643
00644
00645 }
00646
00647 else if( LineSave.ipass == 0 )
00648 {
00649
00650
00651
00652 ASSERT( (chInfo == 'c') || (chInfo == 'h') || (chInfo == 'i') || (chInfo == 'r' ) );
00653
00654
00655 LineSv[LineSave.nsum].chSumTyp = (char)chInfo;
00656
00657
00658 LineSv[LineSave.nsum].sumlin[0] = 0.;
00659 LineSv[LineSave.nsum].emslin[0] = 0.;
00660 LineSv[LineSave.nsum].wavelength = wavelength;
00661 LineSv[LineSave.nsum].sumlin[1] = 0.;
00662 LineSv[LineSave.nsum].emslin[1] = 0.;
00663 strcpy( LineSv[LineSave.nsum].chALab, chLab );
00664 }
00665
00666
00667 ++LineSave.nsum;
00668
00669
00670
00671
00672 DEBUG_EXIT( "linadd()" );
00673
00674 return;
00675 }
00676
00677
00678 static double EnergyRyd;
00679
00680 static bool lgEnergyRydSet=false;
00681
00682
00683 double emergent_line(
00684
00685 double emissivity_in ,
00686
00687 double emissivity_out ,
00688
00689 long int ipCont )
00690 {
00691
00692 double emergent_in , emergent_out;
00693 long int i = ipCont-1;
00694
00695 DEBUG_ENTRY( "emergent_line()" );
00696
00697 ASSERT( i >= 0 && i < rfield.nupper-1 );
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708 if( iteration == 1 )
00709 {
00710
00711
00712 # if 0
00713 emergent_in = (emissivity_in + emissivity_out * opac.albedo[i]) *
00714 opac.E2TauAbsFace[i];
00715 emergent_out = 0.;
00716 # endif
00717 emergent_in = emissivity_in*opac.E2TauAbsFace[i];
00718 emergent_out = emissivity_out;
00719 }
00720 else
00721 {
00722 if( geometry.lgSphere )
00723 {
00724
00725
00726
00727 emergent_in = emissivity_in * opac.E2TauAbsFace[i] *opac.E2TauAbsTotal[i];
00728
00729
00730 emergent_out = emissivity_out * opac.E2TauAbsOut[i];
00731 }
00732 else
00733 {
00734
00735 double reflected = emissivity_out * opac.albedo[i] * (1.-opac.E2TauAbsOut[i]);
00736
00737 emergent_in = (emissivity_in + reflected) * opac.E2TauAbsFace[i];
00738
00739 emergent_out = (emissivity_out - reflected) * opac.E2TauAbsOut[i];
00740 }
00741 }
00742
00743 DEBUG_EXIT( "emergent_line()" );
00744
00745
00746 return( (emergent_in + emergent_out) );
00747 }
00748
00749
00750 void lindst(double xInten,
00751 float wavelength,
00752 const char *chLab,
00753 long int ipnt,
00754 char chInfo,
00755 bool lgOutToo)
00756 {
00757 double saveemis;
00758
00759 DEBUG_ENTRY( "lindst()" );
00760
00761
00762
00763 if( LineSave.ipass > 0 && xInten > 0. )
00764 {
00765
00766 ASSERT( xInten >= 0.);
00767
00768 LineSv[LineSave.nsum].emslin[0] = xInten;
00769
00770 LineSv[LineSave.nsum].sumlin[0] += xInten*radius.dVeff;
00771
00772
00773
00774
00775
00776
00777
00778 if( lgOutToo )
00779 {
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789 rfield.outlin[ipnt-1] += (float)(xInten/(rfield.anu[ipnt-1]*EN1RYD)*
00790 radius.dVolOutwrd*opac.ExpZone[ipnt-1]);
00791 rfield.reflin[ipnt-1] += (float)(xInten/(rfield.anu[ipnt-1]*EN1RYD)*
00792 radius.dVolReflec);
00793 }
00794
00795
00796
00797
00798
00799 if( ipnt <= rfield.nflux )
00800 {
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810 saveemis = emergent_line(
00811 xInten*rt.fracin , xInten*(1.-rt.fracin) , ipnt );
00812 LineSv[LineSave.nsum].emslin[1] = saveemis;
00813 LineSv[LineSave.nsum].sumlin[1] += saveemis*radius.dVeff;
00814 }
00815 }
00816 else if( LineSave.ipass == 0 )
00817 {
00818 ASSERT( ipnt > 0 );
00819
00820
00821
00822
00823
00824
00825
00826
00827 if( !lgEnergyRydSet )
00828 EnergyRyd = rfield.anu[ipnt-1];
00829
00830 LineSv[LineSave.nsum].xLineEnergy = (float)EnergyRyd;
00831
00832
00833
00834 ASSERT( (chInfo == 'c') || (chInfo == 'h') || (chInfo == 'i') || (chInfo == 'r' ) );
00835 LineSv[LineSave.nsum].chSumTyp = (char)chInfo;
00836
00837 LineSv[LineSave.nsum].sumlin[0] = 0.;
00838 LineSv[LineSave.nsum].emslin[0] = 0.;
00839 LineSv[LineSave.nsum].wavelength = wavelength;
00840 LineSv[LineSave.nsum].emslin[1] = 0.;
00841 LineSv[LineSave.nsum].sumlin[1] = 0.;
00842
00843 strcpy( LineSv[LineSave.nsum].chALab, chLab );
00844 }
00845
00846
00847
00848 # if 0
00849
00850 if( gv.lgDustOn )
00851 {
00852 if( LineSave.ipass > 0 )
00853 {
00854 }
00855 else if( LineSave.ipass == 0 )
00856 {
00857
00858 LineDSv[LineSave.nsum].sumlin[LineSave.lgLineEmergent] = 0.;
00859 LineDSv[LineSave.nsum].wavelength = wavelength;
00860 strcpy( LineDSv[LineSave.nsum].chALab, chLab );
00861 }
00862
00863
00864
00865 ++LineSave.nsum;
00866 }
00867 # endif
00868
00869
00870 ++LineSave.nsum;
00871
00872
00873 lgEnergyRydSet = false;
00874 EnergyRyd = 0.;
00875
00876 DEBUG_EXIT( "lindst()" );
00877
00878 return;
00879 }
00880
00881
00882 void PntForLine(
00883
00884 double wavelength,
00885
00886 const char *chLabel,
00887
00888
00889 long int *ipnt)
00890 {
00891
00892
00893
00894
00895
00896 # define MAXFORLIN 1000
00897 static long int ipForLin[MAXFORLIN]={0};
00898
00899
00900 static long int nForLin;
00901
00902 DEBUG_ENTRY( "PntForLine()" );
00903
00904
00905 ASSERT( wavelength >= 0. );
00906
00907 if( wavelength == 0. )
00908 {
00909
00910 nForLin = 0;
00911
00912 EnergyRyd = 0.;
00913
00914 lgEnergyRydSet = false;
00915 }
00916 else
00917 {
00918
00919 if( LineSave.ipass > 0 )
00920 {
00921
00922 *ipnt = ipForLin[nForLin];
00923 }
00924 else if( LineSave.ipass == 0 )
00925 {
00926
00927 if( nForLin >= MAXFORLIN )
00928 {
00929 fprintf( ioQQQ, "PROBLEM %5ld lines is too many for PntForLine.\n",
00930 nForLin );
00931 fprintf( ioQQQ, " Increase the value of maxForLine everywhere in the code.\n" );
00932 puts( "[Stop in pntforline]" );
00933 cdEXIT(EXIT_FAILURE);
00934 }
00935
00936
00937 EnergyRyd = RYDLAM/wavelength;
00938
00939 lgEnergyRydSet = true;
00940 ipForLin[nForLin] = ipLineEnergy(EnergyRyd,chLabel , 0);
00941 *ipnt = ipForLin[nForLin];
00942 }
00943 else
00944 {
00945
00946 *ipnt = 0;
00947 }
00948 ++nForLin;
00949 }
00950
00951 DEBUG_EXIT( "PntForLine()" );
00952
00953 return;
00954 }
00955
00956 static float ExtraInten;
00957
00958
00959 void PutLine(EmLine * t)
00960 {
00961 char chLabel[5];
00962 float wl;
00963 double xIntensity,
00964 other,
00965 xIntensity_in;
00966
00967 DEBUG_ENTRY( "PutLine()" );
00968
00969
00970
00971 ASSERT( t->ipCont > 0. );
00972
00973
00974
00975 if( LineSave.ipass == 0 )
00976 {
00977
00978
00979 chIonLbl(chLabel, t);
00980
00981
00982 wl = t->WLAng;
00983
00984
00985
00986 xIntensity = 0.;
00987 }
00988 else
00989 {
00990
00991
00992
00993 chLabel[0] = '\n';
00994 wl = 0.;
00995
00996
00997 xIntensity = t->xIntensity + ExtraInten;
00998 }
00999
01000
01001
01002
01003
01004
01005
01006
01007
01008 ExtraInten = 0.;
01009
01010
01011 rt.fracin = t->FracInwd;
01012 lindst(xIntensity,
01013 wl,
01014 chLabel,
01015 t->ipCont,
01016
01017 'i',
01018
01019 false);
01020 rt.fracin = 0.5;
01021
01022
01023
01024 xIntensity_in = xIntensity*t->FracInwd;
01025 linadd(xIntensity_in,wl,"Inwd",'i');
01026
01027
01028 other = (float)t->cool;
01029 linadd(other,wl,"Coll",'i');
01030
01031
01032 other = (float)(t->PopOpc * t->pump * (1.-t->ColOvTot) * t->EnergyErg);
01033 linadd(other,wl,"Pump",'i');
01034
01035
01036 other = (float)t->heat;
01037 linadd(other,wl,"Heat",'i');
01038
01039 # if 0
01040
01041 if( gv.lgDustOn )
01042 {
01043 if( LineSave.ipass > 0 )
01044 {
01045
01046
01047
01048 if( t->ipCont <= rfield.nflux )
01049 {
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061 LineDSv[LineSave.nsum].sumlin[LineSave.lgLineEmergent] += (float)emergent_line(
01062 xIntensity_in , xIntensity_out , t->ipCont );
01063 }
01064 }
01065 else if( LineSave.ipass == 0 )
01066 {
01067
01068 LineDSv[LineSave.nsum].sumlin[LineSave.lgLineEmergent] = 0.;
01069 LineDSv[LineSave.nsum].wavelength = wl;
01070 strcpy( LineDSv[LineSave.nsum].chALab, chLabel );
01071 }
01072
01073
01074
01075 ++LineSave.nsum;
01076 }
01077 # endif
01078
01079 DEBUG_EXIT( "PutLine()" );
01080
01081 return;
01082 }
01083
01084
01085
01086 void PutExtra(double Extra)
01087 {
01088
01089 DEBUG_ENTRY( "PutExtra()" );
01090
01091 ExtraInten = (float)Extra;
01092
01093 DEBUG_EXIT( "PutExtra()" );
01094
01095 return;
01096 }
01097
01098
01099
01100
01101 void EmLineJunk( EmLine * t )
01102 {
01103
01104 DEBUG_ENTRY( "EmLineJunk()" );
01105
01106
01107 t->TauCon = -FLT_MAX;
01108
01109
01110 t->TauIn = -FLT_MAX;
01111 t->TauTot = -FLT_MAX;
01112
01113
01114 t->iRedisFun = INT_MIN;
01115
01116
01117 t->ipCont = -10000;
01118
01119
01120 t->ipFine = -10000;
01121
01122
01123 t->FracInwd = -FLT_MAX;
01124
01125
01126 t->pump = -FLT_MAX;
01127
01128
01129 t->xIntensity = -FLT_MAX;
01130
01131
01132 t->phots = -FLT_MAX;
01133
01134
01135 t->gf = -FLT_MAX;
01136
01137
01138 t->Pesc = -FLT_MAX;
01139 t->Pdest = -FLT_MAX;
01140 t->Pelec_esc = -FLT_MAX;
01141
01142
01143 t->dampXvel = -FLT_MAX;
01144 t->damp = -FLT_MAX;
01145
01146
01147 t->cool = -FLT_MAX;
01148 t->heat = -FLT_MAX;
01149
01150
01151 t->ColOvTot = -FLT_MAX;
01152
01153
01154 t->cs = -FLT_MAX;
01155
01156
01157 t->IonStg = -10000;
01158
01159
01160 t->nelem = -10000;
01161
01162
01163 t->WLAng = -FLT_MAX;
01164
01165
01166 t->EnergyK = -FLT_MAX;
01167
01168 t->EnergyErg = -FLT_MAX;
01169
01170 t->EnergyWN = -FLT_MAX;
01171
01172
01173 t->opacity = -FLT_MAX;
01174
01175
01176 t->gLo = -FLT_MAX;
01177 t->gHi = -FLT_MAX;
01178
01179
01180 t->PopLo = -FLT_MAX;
01181 t->PopHi = -FLT_MAX;
01182 t->PopOpc = -FLT_MAX;
01183
01184
01185 t->Aul = -FLT_MAX;
01186
01187
01188 t->ots = -FLT_MAX;
01189
01190 DEBUG_EXIT( "EmLineJunk()" );
01191
01192 return;
01193 }
01194
01195
01196
01197
01198 void EmLineZero( EmLine * t )
01199 {
01200
01201 DEBUG_ENTRY( "EmLineZero()" );
01202
01203
01204
01205 t->TauCon = opac.taumin;
01206
01207
01208
01209 t->TauIn = opac.taumin;
01210
01211
01212 t->TauTot = 1e20f;
01213
01214
01215
01216 t->FracInwd = 1.;
01217
01218
01219 t->pump = 0.;
01220
01221
01222 t->xIntensity = 0.;
01223
01224
01225 t->phots = 0.;
01226
01227
01228
01229 t->Pesc = 1.;
01230 t->Pdest = 0.;
01231 t->Pelec_esc = 0.;
01232
01233
01234 t->cool = 0.;
01235 t->heat = 0.;
01236
01237
01238 t->ColOvTot = 0.;
01239
01240
01241 t->PopLo = 0.;
01242 t->PopHi = 0.;
01243 t->PopOpc = 0.;
01244
01245
01246 t->ots = 0.;
01247
01248 DEBUG_EXIT( "EmLineZero()" );
01249
01250 return;
01251 }
01252
01253
01254 void LineConvRate2CS( EmLine * t , float rate )
01255 {
01256
01257 DEBUG_ENTRY( "LineConvRate2CS()" );
01258
01259
01260
01261
01262 t->cs = rate * t->gHi / (float)dense.cdsqte;
01263
01264
01265
01266 ASSERT( t->cs >= 0. );
01267
01268 DEBUG_EXIT( "LineConvRate2CS()" );
01269
01270 return;
01271 }
01272
01273
01274
01275 double ConvRate2CS( float gHi , float rate )
01276 {
01277
01278 double cs;
01279
01280 DEBUG_ENTRY( "ConvRate2CS()" );
01281
01282
01283
01284
01285 cs = rate * gHi / dense.cdsqte;
01286
01287
01288
01289 ASSERT( cs >= 0. );
01290
01291 DEBUG_EXIT( "ConvRate2CS()" );
01292
01293 return cs;
01294 }
01295
01296
01297 bool lgTauGood( EmLine * t)
01298 {
01299 bool lgGoodTau;
01300
01301 DEBUG_ENTRY( "lgTauGood()" );
01302
01303
01304 if( (t->TauTot*0.9 - t->TauIn < 0. && t->TauIn > 0.) &&
01305 t->TauTot != opac.taumin )
01306
01307 {
01308
01309 lgGoodTau = false;
01310 }
01311 else
01312 {
01313 lgGoodTau = true;
01314 }
01315
01316 DEBUG_EXIT( "lgTauGood()" );
01317
01318 return lgGoodTau;
01319
01320 }
01321
01322
01323 static void gbar0(double ex,
01324 float *g)
01325 {
01326 double a,
01327 b,
01328 c,
01329 d,
01330 y;
01331
01332 DEBUG_ENTRY( "gbar0()" );
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347 y = ex/phycon.te;
01348 if( y < 0.01 )
01349 {
01350 *g = (float)(0.29*(log(1.0+1.0/y) - 0.4/POW2(y + 1.0))/exp(y));
01351 }
01352 else
01353 {
01354 if( y > 10.0 )
01355 {
01356 *g = (float)(0.066/sqrt(y));
01357 }
01358 else
01359 {
01360 a = 1.5819068e-02;
01361 b = 1.3018207e00;
01362 c = 2.6896230e-03;
01363 d = 2.5486007e00;
01364 d = log(y/c)/d;
01365 *g = (float)(a + b*exp(-0.5*POW2(d)));
01366 }
01367 }
01368
01369 DEBUG_EXIT( "gbar0()" );
01370
01371 return;
01372 }
01373
01374
01375 static void gbar1(double ex,
01376 float *g)
01377 {
01378 double y;
01379
01380 DEBUG_ENTRY( "gbar1()" );
01381
01382
01383
01384
01385
01386
01387
01388
01389
01390
01391
01392
01393
01394
01395 y = ex/phycon.te;
01396 *g = (float)(0.6 + 0.28*(log(1.0+1.0/y) - 0.4/POW2(y + 1.0)));
01397
01398 DEBUG_EXIT( "gbar1()" );
01399
01400 return;
01401 }
01402
01403
01404 void MakeCS(EmLine * t)
01405 {
01406 long int ion;
01407 double Abun,
01408 cs;
01409 float
01410 gbar;
01411
01412 DEBUG_ENTRY( "MakeCS()" );
01413
01414
01415
01416
01417 ion = t->IonStg;
01418
01419 Abun = dense.xIonDense[ t->nelem -1 ][ ion-1 ];
01420 if( Abun <= 0. )
01421 {
01422 gbar = 1.;
01423 }
01424 else
01425 {
01426
01427 if( ion == 1 )
01428 {
01429
01430 gbar0(t->EnergyK,&gbar);
01431 }
01432 else
01433 {
01434
01435 gbar1(t->EnergyK,&gbar);
01436 }
01437 }
01438
01439
01440 cs = gbar*(14.5104/WAVNRYD)*t->gf/t->EnergyWN;
01441
01442
01443 t->cs = (float)cs;
01444
01445 DEBUG_EXIT( "MakeCS()" );
01446
01447 return;
01448 }
01449
01450
01451 double totlin(
01452
01453
01454
01455
01456 int chInfo)
01457 {
01458 long int i;
01459 double totlin_v;
01460
01461 DEBUG_ENTRY( "totlin()" );
01462
01463
01464
01465
01466
01467
01468 if( (chInfo != 'i' && chInfo != 'r') && chInfo != 'c' )
01469 {
01470 fprintf( ioQQQ, " TOTLIN does not understand chInfo=%c\n",
01471 chInfo );
01472 puts( "[Stop in totlin]" );
01473 cdEXIT(EXIT_FAILURE);
01474 }
01475
01476
01477
01478 totlin_v = 0.;
01479 for( i=0; i < LineSave.nsum; i++ )
01480 {
01481 if( LineSv[i].chSumTyp == chInfo )
01482 {
01483 totlin_v += LineSv[i].sumlin[0];
01484 }
01485 }
01486
01487 DEBUG_EXIT( "totlin()" );
01488
01489 return( totlin_v );
01490 }
01491
01492
01493
01494 void FndLineHt(long int *level,
01495
01496 long int *ipStrong,
01497 double *Strong)
01498 {
01499 long int i;
01500
01501 DEBUG_ENTRY( "FndLineHt()" );
01502
01503 *Strong = 0.;
01504 *level = 0;
01505
01506
01507 for( i=1; i <= nLevel1; i++ )
01508 {
01509
01510 if( TauLines[i].heat > *Strong )
01511 {
01512 *ipStrong = i;
01513 *level = 1;
01514 *Strong = TauLines[i].heat;
01515 }
01516 }
01517
01518
01519 for( i=0; i < nWindLine; i++ )
01520 {
01521 if( TauLine2[i].IonStg < TauLine2[i].nelem+1-NISO )
01522 {
01523
01524 if( TauLine2[i].heat > *Strong )
01525 {
01526 *ipStrong = i;
01527 *level = 2;
01528 *Strong = TauLine2[i].heat;
01529 }
01530 }
01531 }
01532
01533
01534 for( i=0; i < nCORotate; i++ )
01535 {
01536
01537 if( C12O16Rotate[i].heat > *Strong )
01538 {
01539 *ipStrong = i;
01540 *level = 3;
01541 *Strong = C12O16Rotate[i].heat;
01542 }
01543 }
01544 for( i=0; i < nCORotate; i++ )
01545 {
01546
01547 if( C13O16Rotate[i].heat > *Strong )
01548 {
01549 *ipStrong = i;
01550 *level = 4;
01551 *Strong = C13O16Rotate[i].heat;
01552 }
01553 }
01554
01555
01556 for( i=0; i < nHFLines; i++ )
01557 {
01558
01559 if( HFLines[i].heat > *Strong )
01560 {
01561 *ipStrong = i;
01562 *level = 5;
01563 *Strong = HFLines[i].heat;
01564 }
01565 }
01566
01567 DEBUG_EXIT( "FndLineHt()" );
01568
01569 return;
01570 }
01571
01572