00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "dense.h"
00007 #include "prt.h"
00008 #include "helike.h"
00009 #include "iso.h"
00010 #include "atmdat.h"
00011 #include "lines.h"
00012 #include "lines_service.h"
00013 #include "phycon.h"
00014 #include "physconst.h"
00015 #include "taulines.h"
00016 #include "thirdparty.h"
00017 #include "trace.h"
00018
00019 #define NUMTEMPS 21
00020 #define NUMDENS 14
00021 typedef struct
00022 {
00023
00024 long int ipHi;
00025 long int ipLo;
00026
00027 char label[5];
00028
00029 } stdLines;
00030
00031 STATIC void GetStandardHeLines(void);
00032 STATIC double TempInterp2( double* TempArray , double* ValueArray, long NumElements, double Te );
00033 STATIC void DoSatelliteLines( long nelem );
00034
00035 static bool lgFirstRun = true;
00036 static double CaBDensities[NUMDENS];
00037 static double CaBTemps[NUMTEMPS];
00038 static long NumLines;
00039 static double ****CaBIntensity;
00040 static stdLines **CaBLines;
00041
00042 void lines_helium(void)
00043 {
00044 long ipISO = ipHE_LIKE;
00045 long int i, nelem, ipHi, ipLo;
00046 char chLabel[5]=" ";
00047
00048 double
00049 sum,
00050 photons_3889_plus_7065 = 0.;
00051
00052 DEBUG_ENTRY( "lines_helium()" );
00053
00054 if( trace.lgTrace )
00055 fprintf( ioQQQ, " prt_lines_helium called\n" );
00056
00057
00058
00059 ASSERT( iso_sp[ipHE_LIKE][ipHELIUM].n_HighestResolved_max >= 3 );
00060
00061 i = StuffComment( "He-like iso-sequence" );
00062 linadd( 0., (realnum)i , "####", 'i',
00063 " start He-like iso sequence");
00064
00065 linadd(MAX2(0.,iso_sp[ipHE_LIKE][ipHELIUM].xLineTotCool),506,"Clin",'c',
00066 " total collisional cooling due to all HeI lines ");
00067
00068 linadd(MAX2(0.,-iso_sp[ipHE_LIKE][ipHELIUM].xLineTotCool),506,"Hlin",'h' ,
00069 " total collisional heating due to all HeI lines ");
00070
00071
00072 if( lgFirstRun )
00073 {
00074 GetStandardHeLines();
00075 lgFirstRun = false;
00076 }
00077
00078
00079 for( nelem=ipISO; nelem < LIMELM; nelem++ )
00080 {
00081 if( dense.lgElmtOn[nelem] )
00082 {
00083 t_iso_sp* sp = &iso_sp[ipHE_LIKE][nelem];
00084
00085 ASSERT( sp->n_HighestResolved_max >= 3 );
00086
00087 if( nelem == ipHELIUM )
00088 {
00089 for( i=0; i< NumLines; i++ )
00090 {
00091 ipHi = CaBLines[nelem][i].ipHi;
00092 ipLo = CaBLines[nelem][i].ipLo;
00093 double intens_at_Te[NUMDENS];
00094 for( long ipDens = 0; ipDens < NUMDENS; ++ipDens )
00095 intens_at_Te[ipDens] = TempInterp2( CaBTemps, CaBIntensity[nelem][i][ipDens], NUMTEMPS, phycon.te );
00096 double intens = linint( CaBDensities, intens_at_Te, NUMDENS, log10(dense.eden) );
00097 intens = pow( 10., intens ) * dense.xIonDense[nelem][nelem+1-ipISO]*dense.eden;
00098 ASSERT( intens >= 0. );
00099 linadd( intens, atmdat.CaseBWlHeI[i], CaBLines[nelem][i].label,'i', "Case B intensity ");
00100 }
00101 }
00102
00103
00104 if( LineSave.ipass == 0 )
00105 {
00106
00107
00108 chIonLbl(chLabel, nelem+1, nelem+1-ipISO);
00109 }
00110 for( vector<two_photon>::iterator tnu = sp->TwoNu.begin(); tnu != sp->TwoNu.end(); ++tnu )
00111 {
00112 fixit();
00113 fixit();
00114 linadd( tnu->AulTotal * tnu->E2nu * EN1RYD * (*tnu->Pop),
00115 0, chLabel, 'r',
00116 " two photon continuum ");
00117
00118 linadd( tnu->induc_dn * tnu->E2nu * EN1RYD * (*tnu->Pop),
00119 22, chLabel ,'i',
00120 " induced two photon emission ");
00121 }
00122
00123
00124
00125
00126 sum = 0.;
00127 for( i=ipHe2p3P0; i <= ipHe2p3P2; i++ )
00128 {
00129 if( sp->trans(i,ipHe1s1S).ipCont() <= 0 )
00130 continue;
00131
00132 sum +=
00133 sp->trans(i,ipHe1s1S).Emis().Aul()*
00134 sp->st[i].Pop()*
00135 (sp->trans(i,ipHe1s1S).Emis().Pesc() +
00136 sp->trans(i,ipHe1s1S).Emis().Pelec_esc() ) *
00137 sp->trans(i,ipHe1s1S).EnergyErg();
00138 }
00139
00140 linadd(sum,sp->trans(ipHe2p3P1,ipHe1s1S).WLAng(),"TOTL",'i' ,
00141 " total emission in He-like intercombination lines from 2p3P to ground ");
00142
00143
00144
00145
00146 long int nLoop = sp->numLevels_max - sp->nCollapsed_max;
00147 if( prt.lgPrnIsoCollapsed )
00148 nLoop = sp->numLevels_max;
00149
00150
00151
00152
00153
00154 for( ipLo=0; ipLo < ipHe2p3P0; ipLo++ )
00155 {
00156 vector<long> EnterTheseLast;
00157 for( ipHi=ipLo+1; ipHi < nLoop; ipHi++ )
00158 {
00159
00160
00161 if( sp->trans(ipHi,ipLo).ipCont() < 1 )
00162 continue;
00163
00164
00165
00166 if( iso_ctrl.lgFSM[ipISO] && ( abs(sp->st[ipHi].l() -
00167 sp->st[ipLo].l())==1 )
00168 && (sp->st[ipLo].l()>1)
00169 && (sp->st[ipHi].l()>1)
00170 && ( sp->st[ipHi].n() ==
00171 sp->st[ipLo].n() ) )
00172 {
00173
00174 if( (sp->st[ipHi].S()==1)
00175 && (sp->st[ipLo].S()==1) )
00176 {
00177 continue;
00178 }
00179 else if( (sp->st[ipHi].S()==3)
00180 && (sp->st[ipLo].S()==3) )
00181 {
00182
00183
00184 sp->trans(ipHi+1,ipLo+1).Emis().phots() =
00185 sp->trans(ipHi,ipLo+1).Emis().Aul()*
00186 sp->st[ipHi].Pop()*
00187 (sp->trans(ipHi,ipLo+1).Emis().Pesc() +
00188 sp->trans(ipHi,ipLo+1).Emis().Pelec_esc() ) +
00189 sp->trans(ipHi+1,ipLo+1).Emis().Aul()*
00190 sp->st[ipHi+1].Pop()*
00191 (sp->trans(ipHi+1,ipLo+1).Emis().Pesc() +
00192 sp->trans(ipHi+1,ipLo+1).Emis().Pelec_esc() );
00193
00194 sp->trans(ipHi+1,ipLo+1).Emis().xIntensity() =
00195 sp->trans(ipHi+1,ipLo+1).Emis().phots() *
00196 sp->trans(ipHi+1,ipLo+1).EnergyErg();
00197
00198 PutLine(sp->trans(ipHi+1,ipLo+1), " ");
00199
00200
00201 sp->trans(ipHi,ipLo).Emis().phots() =
00202 sp->trans(ipHi,ipLo).Emis().Aul()*
00203 sp->st[ipHi].Pop()*
00204 (sp->trans(ipHi,ipLo).Emis().Pesc() +
00205 sp->trans(ipHi,ipLo).Emis().Pelec_esc() ) +
00206 sp->trans(ipHi+1,ipLo).Emis().Aul()*
00207 sp->st[ipHi+1].Pop()*
00208 (sp->trans(ipHi+1,ipLo).Emis().Pesc() +
00209 sp->trans(ipHi+1,ipLo).Emis().Pelec_esc() );
00210
00211 sp->trans(ipHi,ipLo).Emis().xIntensity() =
00212 sp->trans(ipHi,ipLo).Emis().phots() *
00213 sp->trans(ipHi,ipLo).EnergyErg();
00214
00215 PutLine(sp->trans(ipHi,ipLo), " ");
00216 }
00217 }
00218
00219 else if( ipLo==ipHe2s3S && ipHi == ipHe2p3P0 )
00220 {
00221
00222
00223
00224
00225
00226 realnum av_wl = 0.;
00227 sum = 0.;
00228 for( i=ipHe2p3P0; i <= ipHe2p3P2; i++ )
00229 {
00230 sum +=
00231 sp->trans(i,ipLo).Emis().Aul()*
00232 sp->st[i].Pop()*
00233 (sp->trans(i,ipLo).Emis().Pesc() +
00234 sp->trans(i,ipLo).Emis().Pelec_esc() )*
00235 sp->trans(i,ipLo).EnergyErg();
00236 av_wl += sp->trans(i,ipLo).WLAng();
00237 }
00238 av_wl /= 3.;
00239 # if 0
00240 {
00241 # include "elementnames.h"
00242 # include "prt.h"
00243 fprintf(ioQQQ,"DEBUG 2P - 2S multiplet wl %s ",
00244 elementnames.chElementSym[nelem] );
00245 prt_wl( ioQQQ , av_wl );
00246 fprintf(ioQQQ,"\n" );
00247 }
00248 # endif
00249
00250 linadd(sum,av_wl,"TOTL",'i',
00251 "total emission in He-like lines, use average of three line wavelengths " );
00252
00253
00254 linadd(sum,av_wl,chLabel,'i',
00255 "total emission in He-like lines, use average of three line wavelengths " );
00256
00257
00258
00259
00260
00261 sp->trans(ipHi,ipLo).Emis().phots() =
00262 sp->trans(ipHi,ipLo).Emis().Aul()*
00263 sp->st[ipHi].Pop()*
00264 (sp->trans(ipHi,ipLo).Emis().Pesc() +
00265 sp->trans(ipHi,ipLo).Emis().Pelec_esc() );
00266
00267
00268 sp->trans(ipHi,ipLo).Emis().xIntensity() =
00269 sp->trans(ipHi,ipLo).Emis().phots()*
00270 sp->trans(ipHi,ipLo).EnergyErg();
00271
00272 if( iso_ctrl.lgRandErrGen[ipISO] )
00273 {
00274 sp->trans(ipHi,ipLo).Emis().phots() *=
00275 sp->ex[ipHi][ipLo].ErrorFactor[IPRAD];
00276 sp->trans(ipHi,ipLo).Emis().xIntensity() *=
00277 sp->ex[ipHi][ipLo].ErrorFactor[IPRAD];
00278 }
00279 PutLine(sp->trans(ipHi,ipLo), " ");
00280 }
00281
00282 else
00283 {
00284
00285 sp->trans(ipHi,ipLo).Emis().phots() =
00286 sp->trans(ipHi,ipLo).Emis().Aul()*
00287 sp->st[ipHi].Pop()*
00288 (sp->trans(ipHi,ipLo).Emis().Pesc() +
00289 sp->trans(ipHi,ipLo).Emis().Pelec_esc() );
00290
00291
00292
00293 sp->trans(ipHi,ipLo).Emis().xIntensity() =
00294 sp->trans(ipHi,ipLo).Emis().phots()*
00295 sp->trans(ipHi,ipLo).EnergyErg();
00296
00297 if( iso_ctrl.lgRandErrGen[ipISO] )
00298 {
00299 sp->trans(ipHi,ipLo).Emis().phots() *=
00300 sp->ex[ipHi][ipLo].ErrorFactor[IPRAD];
00301 sp->trans(ipHi,ipLo).Emis().xIntensity() *=
00302 sp->ex[ipHi][ipLo].ErrorFactor[IPRAD];
00303 }
00304
00305 if( abs( L_(ipHi) - L_(ipLo) ) != 1 )
00306 {
00307 EnterTheseLast.push_back( ipHi );
00308 continue;
00309 }
00310
00311
00312
00313
00314 PutLine(sp->trans(ipHi,ipLo),
00315 "total intensity of He-like line");
00316 {
00317
00318
00319 enum {DEBUG_LOC=false};
00320 if( DEBUG_LOC )
00321 {
00322 if( nelem==1 && ipLo==0 && ipHi==1 )
00323 {
00324 fprintf(ioQQQ,"he1 626 %.2e %.2e \n",
00325 sp->trans(ipHi,ipLo).Emis().TauIn(),
00326 sp->trans(ipHi,ipLo).Emis().TauTot()
00327 );
00328 }
00329 }
00330 }
00331 }
00332 }
00333
00334
00335 for( vector<long>::iterator it = EnterTheseLast.begin(); it != EnterTheseLast.end(); it++ )
00336 PutLine( sp->trans(*it,ipLo),
00337 "predicted line, all processes included");
00338 }
00339
00340
00341 for( ipHi=ipHe2p3P2+1; ipHi < nLoop; ipHi++ )
00342 {
00343 double sumcool , sumheat ,
00344 save , savecool , saveheat;
00345
00346 sum = 0;
00347 sumcool = 0.;
00348 sumheat = 0.;
00349 for( ipLo=ipHe2p3P0; ipLo <= ipHe2p3P2; ++ipLo )
00350 {
00351 if( sp->trans(ipHi,ipLo).ipCont() <= 0 )
00352 continue;
00353
00354
00355 sp->trans(ipHi,ipLo).Emis().phots() =
00356 sp->trans(ipHi,ipLo).Emis().Aul()*
00357 sp->st[ipHi].Pop()*
00358 (sp->trans(ipHi,ipLo).Emis().Pesc() +
00359 sp->trans(ipHi,ipLo).Emis().Pelec_esc() );
00360
00361
00362
00363 sp->trans(ipHi,ipLo).Emis().xIntensity() =
00364 sp->trans(ipHi,ipLo).Emis().phots()*
00365 sp->trans(ipHi,ipLo).EnergyErg();
00366
00367 if( iso_ctrl.lgRandErrGen[ipISO] )
00368 {
00369 sp->trans(ipHi,ipLo).Emis().phots() *=
00370 sp->ex[ipHi][ipLo].ErrorFactor[IPRAD];
00371 sp->trans(ipHi,ipLo).Emis().xIntensity() *=
00372 sp->ex[ipHi][ipLo].ErrorFactor[IPRAD];
00373 }
00374
00375 sumcool += sp->trans(ipHi,ipLo).Coll().cool();
00376 sumheat += sp->trans(ipHi,ipLo).Coll().heat();
00377 sum += sp->trans(ipHi,ipLo).Emis().xIntensity();
00378 }
00379
00380
00381 if( sp->trans(ipHi,ipHe2p3P2).ipCont() < 1 )
00382 continue;
00383
00384
00385 save = sp->trans(ipHi,ipHe2p3P2).Emis().xIntensity();
00386 savecool = sp->trans(ipHi,ipHe2p3P2).Coll().cool();
00387 saveheat = sp->trans(ipHi,ipHe2p3P2).Coll().heat();
00388
00389 sp->trans(ipHi,ipHe2p3P2).Emis().xIntensity() = sum;
00390 sp->trans(ipHi,ipHe2p3P2).Coll().cool() = sumcool;
00391 sp->trans(ipHi,ipHe2p3P2).Coll().heat() = sumheat;
00392
00393
00394
00395 PutLine(sp->trans(ipHi,ipHe2p3P2),
00396 "predicted line, all processes included");
00397
00398 sp->trans(ipHi,ipHe2p3P2).Emis().xIntensity() = save;
00399 sp->trans(ipHi,ipHe2p3P2).Coll().cool() = savecool;
00400 sp->trans(ipHi,ipHe2p3P2).Coll().heat() = saveheat;
00401 }
00402 for( ipLo=ipHe2p3P2+1; ipLo < nLoop-1; ipLo++ )
00403 {
00404 vector<long> EnterTheseLast;
00405 for( ipHi=ipLo+1; ipHi < nLoop; ipHi++ )
00406 {
00407
00408 if( sp->trans(ipHi,ipLo).ipCont() < 1 )
00409 continue;
00410
00411
00412 sp->trans(ipHi,ipLo).Emis().phots() =
00413 sp->trans(ipHi,ipLo).Emis().Aul()*
00414 sp->st[ipHi].Pop()*
00415 (sp->trans(ipHi,ipLo).Emis().Pesc() +
00416 sp->trans(ipHi,ipLo).Emis().Pelec_esc() );
00417
00418
00419
00420 sp->trans(ipHi,ipLo).Emis().xIntensity() =
00421 sp->trans(ipHi,ipLo).Emis().phots()*
00422 sp->trans(ipHi,ipLo).EnergyErg();
00423
00424 if( iso_ctrl.lgRandErrGen[ipISO] )
00425 {
00426 sp->trans(ipHi,ipLo).Emis().phots() *=
00427 sp->ex[ipHi][ipLo].ErrorFactor[IPRAD];
00428 sp->trans(ipHi,ipLo).Emis().xIntensity() *=
00429 sp->ex[ipHi][ipLo].ErrorFactor[IPRAD];
00430 }
00431
00432 if( abs( L_(ipHi) - L_(ipLo) ) != 1 )
00433 {
00434 EnterTheseLast.push_back( ipHi );
00435 continue;
00436 }
00437
00438
00439 PutLine(sp->trans(ipHi,ipLo),
00440 "predicted line, all processes included");
00441 }
00442
00443
00444 for( vector<long>::iterator it = EnterTheseLast.begin(); it != EnterTheseLast.end(); it++ )
00445 PutLine(sp->trans(*it,ipLo),
00446 "predicted line, all processes included");
00447 }
00448
00449
00450 if( iso_ctrl.lgDielRecom[ipISO] )
00451 DoSatelliteLines(nelem);
00452 }
00453 }
00454
00455 if( iso_sp[ipHE_LIKE][ipHELIUM].n_HighestResolved_max >= 4 &&
00456 ( iso_sp[ipH_LIKE][ipHYDROGEN].n_HighestResolved_max + iso_sp[ipH_LIKE][ipHYDROGEN].nCollapsed_max ) >=8 )
00457 {
00458 t_iso_sp* sp = &iso_sp[ipHE_LIKE][ipHELIUM];
00459 const long ipHe4s3S = iso_sp[ipHE_LIKE][ipHELIUM].QuantumNumbers2Index[4][0][3];
00460 const long ipHe4p3P = iso_sp[ipHE_LIKE][ipHELIUM].QuantumNumbers2Index[4][1][3];
00461
00462
00463
00464 photons_3889_plus_7065 =
00465
00466 sp->trans(ipHe3s3S,ipHe2p3P2).Emis().xIntensity()/
00467 sp->trans(ipHe3s3S,ipHe2p3P2).EnergyErg() +
00468 sp->trans(ipHe3d3D,ipHe2p3P2).Emis().xIntensity()/
00469 sp->trans(ipHe3d3D,ipHe2p3P2).EnergyErg() +
00470 sp->trans(ipHe4s3S,ipHe2p3P2).Emis().xIntensity()/
00471 sp->trans(ipHe4s3S,ipHe2p3P2).EnergyErg() +
00472
00473 sp->trans(ipHe3s3S,ipHe2p3P1).Emis().xIntensity()/
00474 sp->trans(ipHe3s3S,ipHe2p3P1).EnergyErg() +
00475 sp->trans(ipHe3d3D,ipHe2p3P1).Emis().xIntensity()/
00476 sp->trans(ipHe3d3D,ipHe2p3P1).EnergyErg() +
00477 sp->trans(ipHe4s3S,ipHe2p3P1).Emis().xIntensity()/
00478 sp->trans(ipHe4s3S,ipHe2p3P1).EnergyErg() +
00479
00480 sp->trans(ipHe3s3S,ipHe2p3P0).Emis().xIntensity()/
00481 sp->trans(ipHe3s3S,ipHe2p3P0).EnergyErg() +
00482 sp->trans(ipHe3d3D,ipHe2p3P0).Emis().xIntensity()/
00483 sp->trans(ipHe3d3D,ipHe2p3P0).EnergyErg() +
00484 sp->trans(ipHe4s3S,ipHe2p3P0).Emis().xIntensity()/
00485 sp->trans(ipHe4s3S,ipHe2p3P0).EnergyErg() +
00486
00487 sp->trans(ipHe3p3P,ipHe2s3S).Emis().xIntensity()/
00488 sp->trans(ipHe3p3P,ipHe2s3S).EnergyErg() +
00489 sp->trans(ipHe4p3P,ipHe2s3S).Emis().xIntensity()/
00490 sp->trans(ipHe4p3P,ipHe2s3S).EnergyErg();
00491
00492 long upperIndexofH8 = iso_sp[ipH_LIKE][ipHYDROGEN].QuantumNumbers2Index[8][1][2];
00493
00494
00495 photons_3889_plus_7065 +=
00496 iso_sp[ipH_LIKE][ipHYDROGEN].trans(upperIndexofH8,1).Emis().xIntensity()/
00497 iso_sp[ipH_LIKE][ipHYDROGEN].trans(upperIndexofH8,1).EnergyErg();
00498
00499
00500
00501 if( LineSave.WavLNorm > 0 )
00502 photons_3889_plus_7065 *= (ERG1CM*1.e8)/LineSave.WavLNorm;
00503 linadd( photons_3889_plus_7065, 3889., "Pho+", 'i',
00504 "photon sum given in Porter et al. 2007 (astro-ph/0611579)");
00505 }
00506
00507
00508
00509
00510
00511 if( trace.lgTrace )
00512 {
00513 fprintf( ioQQQ, " lines_helium returns\n" );
00514 }
00515 return;
00516 }
00517
00518
00519 STATIC void GetStandardHeLines(void)
00520 {
00521 FILE *ioDATA;
00522 bool lgEOL, lgHIT;
00523 long i, i1, i2;
00524
00525 # define chLine_LENGTH 1000
00526 char chLine[chLine_LENGTH];
00527
00528 DEBUG_ENTRY( "GetStandardHeLines()" );
00529
00530 if( trace.lgTrace )
00531 fprintf( ioQQQ," lines_helium opening he1_case_b.dat\n");
00532
00533 ioDATA = open_data( "he1_case_b.dat", "r" );
00534
00535
00536 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00537 {
00538 fprintf( ioQQQ, " lines_helium could not read first line of he1_case_b.dat.\n");
00539 cdEXIT(EXIT_FAILURE);
00540 }
00541 i = 1;
00542
00543 i1 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00544 i2 = (long)FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00545 NumLines = i2;
00546
00547
00548 if( i1 !=CASEBMAGIC )
00549 {
00550 fprintf( ioQQQ,
00551 " lines_helium: the version of he1_case_ab.dat is not the current version.\n" );
00552 fprintf( ioQQQ,
00553 " lines_helium: I expected to find the number %i and got %li instead.\n" ,
00554 CASEBMAGIC, i1 );
00555 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00556 cdEXIT(EXIT_FAILURE);
00557 }
00558
00559
00560 lgHIT = false;
00561 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00562 {
00563
00564 if( chLine[0] != '#')
00565 {
00566 lgHIT = true;
00567 long j = 0;
00568 lgEOL = false;
00569 i = 1;
00570 while( !lgEOL && j < NUMDENS)
00571 {
00572 CaBDensities[j] = FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00573 ++j;
00574 }
00575 break;
00576 }
00577 }
00578
00579
00580 lgHIT = false;
00581 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00582 {
00583
00584 if( chLine[0] != '#')
00585 {
00586 lgHIT = true;
00587 long j = 0;
00588 lgEOL = false;
00589 i = 1;
00590 while( !lgEOL && j < NUMTEMPS)
00591 {
00592 CaBTemps[j] = FFmtRead(chLine,&i,sizeof(chLine),&lgEOL);
00593 ++j;
00594 }
00595 break;
00596 }
00597 }
00598
00599 if( !lgHIT )
00600 {
00601 fprintf( ioQQQ, " lines_helium could not find line of temperatures in he1_case_ab.dat.\n");
00602 cdEXIT(EXIT_FAILURE);
00603 }
00604
00605
00606 CaBIntensity = (double ****)MALLOC(sizeof(double ***)*(unsigned)LIMELM );
00607
00608 CaBLines = (stdLines **)MALLOC(sizeof(stdLines *)*(unsigned)LIMELM );
00609
00610 for( long nelem=ipHELIUM; nelem<LIMELM; ++nelem )
00611 {
00615 if( nelem != ipHELIUM )
00616 continue;
00617
00618
00619 if( nelem == ipHELIUM || dense.lgElmtOn[nelem])
00620 {
00621
00622 CaBIntensity[nelem] = (double ***)MALLOC(sizeof(double **)*(unsigned)(i2) );
00623 CaBLines[nelem] = (stdLines *)MALLOC(sizeof(stdLines )*(unsigned)(i2) );
00624
00625
00626
00627 for( long j = 0; j < i2; ++j )
00628 {
00629 CaBIntensity[nelem][j] = (double **)MALLOC(sizeof(double*)*(unsigned)NUMDENS );
00630
00631 CaBLines[nelem][j].ipHi = -1;
00632 CaBLines[nelem][j].ipLo = -1;
00633 strcpy( CaBLines[nelem][j].label , " " );
00634
00635 for( long k = 0; k < NUMDENS; ++k )
00636 {
00637 CaBIntensity[nelem][j][k] = (double *)MALLOC(sizeof(double)*(unsigned)NUMTEMPS );
00638 for( long l = 0; l < NUMTEMPS; ++l )
00639 {
00640 CaBIntensity[nelem][j][k][l] = 0.;
00641 }
00642 }
00643 }
00644 }
00645 }
00646
00647
00648 lgHIT = false;
00649 long nelem = ipHELIUM;
00650 int line = 0;
00651 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00652 {
00653 char *chTemp = NULL;
00654
00655
00656 if( (chLine[0] == ' ') || (chLine[0]=='\n') )
00657 break;
00658 if( chLine[0] != '#')
00659 {
00660 lgHIT = true;
00661
00662
00663 long j = 1;
00664
00665
00666 realnum wl = (realnum)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
00667 long ipLo = (long)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
00668 long ipHi = (long)FFmtRead(chLine,&j,sizeof(chLine),&lgEOL);
00669 CaBLines[nelem][line].ipLo = ipLo;
00670 CaBLines[nelem][line].ipHi = ipHi;
00671
00672 ASSERT( CaBLines[nelem][line].ipLo < CaBLines[nelem][line].ipHi );
00673
00674 strncpy( CaBLines[nelem][line].label, "Ca B" , 4 );
00675 CaBLines[nelem][line].label[4] = 0;
00676
00677 t_iso_sp* sp = &iso_sp[ipHE_LIKE][nelem];
00678 if( ipHi < sp->numLevels_max - sp->nCollapsed_max )
00679 atmdat.CaseBWlHeI.push_back( sp->trans(ipHi,ipLo).WLAng() );
00680 else
00681 atmdat.CaseBWlHeI.push_back( wl );
00682
00683 for( long ipDens = 0; ipDens < NUMDENS; ++ipDens )
00684 {
00685 if( read_whole_line( chLine, (int)sizeof(chLine) , ioDATA ) == NULL )
00686 {
00687 fprintf( ioQQQ, " lines_helium could not scan case B lines, current indices: %li %li\n",
00688 CaBLines[nelem][line].ipHi,
00689 CaBLines[nelem][line].ipLo );
00690 cdEXIT(EXIT_FAILURE);
00691 }
00692
00693 chTemp = chLine;
00694 j = 1;
00695 long den = (long)FFmtRead(chTemp,&j,sizeof(chTemp),&lgEOL);
00696 ASSERT( den == ipDens + 1 );
00697
00698 for( long ipTe = 0; ipTe < NUMTEMPS; ++ipTe )
00699 {
00700 double b;
00701 if( (chTemp = strchr_s( chTemp, '\t' )) == NULL )
00702 {
00703 fprintf( ioQQQ, " lines_helium could not scan case B lines, current indices: %li %li\n",
00704 CaBLines[nelem][line].ipHi,
00705 CaBLines[nelem][line].ipLo );
00706 cdEXIT(EXIT_FAILURE);
00707 }
00708 ++chTemp;
00709 sscanf( chTemp, "%le" , &b );
00710 CaBIntensity[nelem][line][ipDens][ipTe] = b;
00711 }
00712 }
00713 line++;
00714 }
00715 }
00716
00717 ASSERT( line == NumLines );
00718 ASSERT( atmdat.CaseBWlHeI.size() == (unsigned)line );
00719
00720
00721 fclose( ioDATA );
00722 return;
00723 }
00724
00726 STATIC double TempInterp2( double* TempArray , double* ValueArray, long NumElements, double Te )
00727 {
00728 long int ipTe=-1;
00729 double rate = 0.;
00730 long i0;
00731
00732 DEBUG_ENTRY( "TempInterp2()" );
00733
00734
00735 if( Te <= TempArray[0] )
00736 {
00737 return ValueArray[0] + log10( TempArray[0] ) - log10( Te );
00738 }
00739 else if( Te >= TempArray[NumElements-1] )
00740 {
00741 return ValueArray[NumElements-1];
00742 }
00743
00744
00745 ipTe = hunt_bisect( TempArray, NumElements, Te );
00746
00747 ASSERT( (ipTe >=0) && (ipTe < NumElements-1) );
00748
00749
00750 const int ORDER = 3;
00751 i0 = max(min(ipTe-ORDER/2,NumElements-ORDER-1),0);
00752 rate = lagrange( &TempArray[i0], &ValueArray[i0], ORDER+1, Te );
00753
00754 return rate;
00755 }
00756
00758
00759
00760 STATIC void DoSatelliteLines( long nelem )
00761 {
00762 long ipISO = ipHE_LIKE;
00763
00764 DEBUG_ENTRY( "DoSatelliteLines()" );
00765
00766 ASSERT( iso_ctrl.lgDielRecom[ipISO] && dense.lgElmtOn[nelem] );
00767
00768 for( long i=0; i < iso_sp[ipISO][nelem].numLevels_max; i++ )
00769 {
00770 double dr_rate = iso_sp[ipISO][nelem].fb[i].DielecRecomb;
00771 TransitionProxy tr = SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][i]];
00772
00773 tr.Emis().phots() =
00774 dr_rate * dense.eden * dense.xIonDense[nelem][nelem+1-ipISO];
00775
00776 tr.Emis().xIntensity() =
00777 tr.Emis().phots() * ERG1CM * tr.EnergyWN();
00778 tr.Emis().pump() = 0.;
00779
00780 PutLine( tr, "iso satellite line" );
00781 }
00782
00783 return;
00784 }