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