00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "taulines.h"
00007 #include "dense.h"
00008 #include "phycon.h"
00009 #include "lines_service.h"
00010 #include "iso.h"
00011 #include "trace.h"
00012 #include "lines.h"
00013 #include "path.h"
00014 #include "helike.h"
00015 #include "thirdparty.h"
00016 #include "physconst.h"
00017
00018
00019
00020
00021
00022
00023 #define NUMTEMPS 22
00024
00025 typedef struct
00026 {
00027
00028 long int ipHi;
00029 long int ipLo;
00030
00031 char label[5];
00032
00033 } stdLines;
00034
00035 static void GetStandardHeLines(void);
00036
00037 static double TempInterp( double* TempArray , double* ValueArray, long NumElements );
00038
00039 static bool lgFirstRun = true;
00040 static double CaABTemps[NUMTEMPS];
00041 static long NumLines;
00042 static double ***CaABIntensity;
00043 static stdLines **CaABLines;
00044
00045 void lines_helium(void)
00046 {
00047 long int i, j,
00048 ipHi,
00049 ipLo,
00050 nelem;
00051 char chLabel[5]=" ";
00052 double
00053
00054
00055
00056
00057
00058 sum,
00059 Pop2_3S,
00060 photons_3889_plus_7065 = 0.;
00061
00062 DEBUG_ENTRY( "lines_helium()" );
00063
00064 if( trace.lgTrace )
00065 {
00066 fprintf( ioQQQ, " lines_helium called\n" );
00067 }
00068
00069
00070 i = StuffComment( "He-like iso-sequence" );
00071 linadd( 0., (float)i , "####", 'i');
00072
00073 #if 0
00074
00075 if( phycon.te <= 1e4 )
00076 {
00077 em = 2.42e-22/(phycon.te/phycon.te10);
00078 }
00079 else
00080 {
00081 em = 8.79e-22/phycon.te/phycon.te03/phycon.te01;
00082 }
00083
00084
00085
00086 rec = dense.eden*dense.xIonDense[ipHELIUM][1]*em;
00087 linadd( rec , 4471 , "Ca B",'i');
00088
00089
00090 cb5876 = dense.eden*dense.xIonDense[ipHELIUM][1]*4.22e-21/(phycon.te*phycon.te10);
00091 linadd( cb5876 , 5876 , "Ca B" , 'i' );
00092
00093
00094
00095 if( dense.eden < 1e8 )
00096 {
00097 double D = 1. + 3131.*100./phycon.sqrte / dense.eden;
00098 em = (6.78*0.52581*phycon.te07*sexp(3.776e4/phycon.te) +
00099 1.67*3.9811/(phycon.te10*phycon.te05)*sexp(4.545e4/phycon.te ) +
00100 0.60*22.9087/(phycon.te30*phycon.te04)*sexp(4.901e4/phycon.te) ) / D;
00101 }
00102 else
00103 {
00104 em = 0.;
00105 }
00106
00107 linadd( cb5876*(1.+em) , 5876 , "wCol" , 'i' );
00108
00109
00110 rec = dense.eden*dense.xIonDense[ipHELIUM][1]*1.32e-21/(phycon.te*phycon.te10*phycon.te01);
00111 linadd( rec , 6678 , "Ca B" ,'i' );
00112
00113 em = dense.eden*dense.xIonDense[ipHELIUM][1];
00114 cb206 = 1.064e-22/phycon.te70/phycon.te10/phycon.te03/phycon.te03*
00115 em;
00116
00117
00118 linadd(cb206,20600.,"Ca B",'i');
00119
00120
00121 cb5016 = 5.43e-23/phycon.te70/phycon.te10*em;
00122 linadd( cb5016 , 5016 , "Ca B", 'i' );
00123 #endif
00124
00125 #if 0
00126 {
00127 float DoublyIonizedWavelengths[18] = {
00128 18.96f, 19.2418f, 19.306f, 19.2831f,19.311f, 19.5747f, 19.2134f, 19.2186f, 19.439f,
00129 19.3217f,19.3126f,19.3217f,19.307f, 19.3179f,19.3269f, 19.5405f, 19.5498f, 19.5348f
00130 };
00131 double energyRyd[18] = {
00132 34.77502169,34.77502169,33.91347631,33.96960771,33.90125237,33.96960771,
00133 34.74903112,34.74903112,34.74903112,34.48322096,34.50539991,34.48322096,
00134 34.51915012,34.50539991,34.48322096,34.50539991,34.48322096,34.51915012};
00135 double stat_weight_hi[18] = {
00136 3.,3.,5.,3.,1.,3.,5.,5.,5.,
00137 3.,5.,3.,1.,5.,3.,5.,3.,1.
00138 };
00139 double stat_weight_lo[18] = {
00140 3.,1.,3.,3.,3.,1.,3.,5.,3.,
00141 1.,3.,3.,3.,5.,5.,3.,3.,3.
00142 };
00143 double branch_ratio[18] = {
00144 5.61E-05,1.00E+00,1.00E+00,1.00E+00,1.00E+00,5.60E-05,1.14E-06,4.36E-04,1.00E+00,
00145 1.11E-01,8.36E-02,8.36E-02,3.33E-01,2.50E-01,1.39E-01,1.02E-04,1.37E-05,2.77E-05
00146 };
00147 double sum_of_As[18] = {
00148 2.44014E+12,2.44014E+12,2.32013E+12,2.32013E+12,2.32013E+12,2.32013E+12,
00149 4.70206E+12,4.70206E+12,4.70206E+12,1.4002E+13,1.4002E+13,1.4002E+13,
00150 1.4002E+13,1.4002E+13,1.4002E+13,1.4002E+13,1.4002E+13,1.4002E+13
00151 };
00152
00153 long nelem = ipOXYGEN;
00154 double total_dr_rate = 0.;
00155
00156 for( i=0; i<18; i++ )
00157 {
00158 double ConBoltz, LTE_pop, factor1, ConvLTEPOP;
00159 double xInWrd;
00160 long ip;
00161 double temp;
00162 EmLine line1;
00163
00164 EmLineZero( &line1 );
00165
00166 factor1 = HION_LTE_POP*dense.AtomicWeight[nelem]/
00167 (dense.AtomicWeight[nelem]+ELECTRON_MASS/ATOMIC_MASS_UNIT);
00168
00169
00170 ConvLTEPOP = pow(factor1,1.5)/(2.*iso.stat_ion[ipHE_LIKE])/phycon.te32;
00171
00172
00173 ConBoltz = dsexp(energyRyd[i]/phycon.te_ryd);
00174
00175 if( ConBoltz >= SMALLDOUBLE )
00176 {
00177
00178
00179
00180
00181 LTE_pop = stat_weight_hi[i] * ConBoltz * ConvLTEPOP;
00182 }
00183
00184 temp = sum_of_As[i] * branch_ratio[i] * LTE_pop;
00185 total_dr_rate += temp;
00186 fprintf( ioQQQ, "%li\t%e\t%e\n", i, temp, total_dr_rate );
00187
00188 line1.WLAng = DoublyIonizedWavelengths[i];
00189 line1.EnergyWN = (float)(RYD_INF*RYDLAM)/DoublyIonizedWavelengths[i];
00190
00191 line1.phots = LTE_pop * sum_of_As[i] * branch_ratio[i] *
00192 dense.eden * dense.eden * dense.xIonDense[nelem][nelem];
00193
00194
00195
00196 line1.xIntensity = line1.phots * EN1RYD * (RYDLAM/line1.WLAng);
00197
00198 line1.ipCont = (int)ipLineEnergy(line1.EnergyWN * WAVNRYD , "O 7" ,
00199 iso.ipIsoLevNIonCon[ipHE_LIKE][ipOXYGEN][ipHe1s1S]);
00200
00201 line1.ipFine = ipFineCont(line1.EnergyWN * WAVNRYD );
00202
00203 line1.iRedisFun = ipCRDW;
00204 line1.nelem = ipOXYGEN + 1;
00205 line1.IonStg = 7;
00206 line1.Aul = (float)(sum_of_As[i] * branch_ratio[i]);
00207 line1.EnergyErg = (float)ERG1CM * line1.EnergyWN;
00208 line1.EnergyK = (float)T1CM * line1.EnergyWN;
00209 line1.gHi = (float)stat_weight_hi[i];
00210 line1.gLo = (float)stat_weight_lo[i];
00211 line1.gf = (float)GetGF(line1.Aul, line1.EnergyWN, line1.gHi);
00212 ASSERT(line1.gf > 0.);
00213 line1.FracInwd = 0.;
00214
00215 PutLine(&line1);
00216
00217 ip = line1.ipCont-1;
00218
00219
00220 xInWrd = line1.phots*line1.FracInwd;
00221
00222
00223 rfield.reflin[ip] += (float)(xInWrd*radius.BeamInIn);
00224
00225
00226
00227
00228
00229 rfield.outlin[ip] += (float)(xInWrd*radius.BeamInOut*opac.tmn[ip]);
00230
00231
00232 rfield.outlin[ip] +=
00233 (float)(line1.phots*(1. - line1.FracInwd)*
00234 radius.BeamOutOut*opac.tmn[ip]);
00235 }
00236 fprintf( ioQQQ, "dr rate\t%e\t%e\t%e\n\n",
00237 total_dr_rate, dense.eden, dense.xIonDense[nelem][nelem] );
00238 }
00239 #endif
00240
00241
00242 if( lgFirstRun )
00243 {
00244 GetStandardHeLines();
00245 lgFirstRun = false;
00246 }
00247
00248
00249 for( nelem=1; nelem < LIMELM; nelem++ )
00250 {
00251 if( dense.lgElmtOn[nelem] )
00252 {
00253 if( nelem == ipHELIUM )
00254 {
00255 double *qTotEff;
00256
00257
00258
00259 qTotEff = (double*)MALLOC(sizeof(double)*(unsigned)(iso.numLevels_max[ipHE_LIKE][nelem]) );
00260
00261 qTotEff[0] = 0.;
00262 qTotEff[1] = 0.;
00263
00264 for( i=ipHe2s3S+1; i<iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem]; i++ )
00265 {
00266 qTotEff[i] = 0.;
00267 for( j = i; j<iso.numLevels_max[ipHE_LIKE][nelem]-iso.nCollapsed_max[ipHE_LIKE][nelem]; j++ )
00268 {
00269
00270
00271 qTotEff[i] +=
00272 EmisLines[ipHE_LIKE][nelem][j][ipHe2s3S].ColUL*dense.EdenHCorr*
00273 iso.Boltzmann[ipHE_LIKE][nelem][j][ipHe2s3S] *
00274 (double)EmisLines[ipHE_LIKE][nelem][j][ipHe2s3S].gHi /
00275 (double)EmisLines[ipHE_LIKE][nelem][j][ipHe2s3S].gLo*
00276 helike.CascadeProb[nelem][j][i];
00277
00278 }
00279 }
00280
00281
00282
00283
00284 Pop2_3S = dense.eden*iso.RadRec_caseB[ipHE_LIKE][nelem]/
00285 ( EmisLines[ipHE_LIKE][nelem][ipHe2s3S][ipHe1s1S].Aul+ dense.eden*helike.qTot2TripS[nelem]);
00286
00287 for( i=0; i< NumLines; i++ )
00288 {
00289 ipHi = CaABLines[nelem][i].ipHi;
00290 ipLo = CaABLines[nelem][i].ipLo;
00291
00292
00293
00294 if( ipHi <= iso.n_HighestResolved_max[ipHE_LIKE][nelem]*(iso.n_HighestResolved_max[ipHE_LIKE][nelem]+1))
00295 {
00296 double intens = TempInterp( CaABTemps , CaABIntensity[nelem][i], NUMTEMPS )*
00297 dense.xIonDense[nelem][nelem]*dense.eden;
00298
00299 linadd( intens,
00300 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng,
00301 CaABLines[nelem][i].label,
00302 'i');
00303
00304 if( nMatch("Ca B",CaABLines[nelem][i].label) )
00305 {
00306
00307
00308
00309 linadd( intens +
00310 Pop2_3S*qTotEff[ipHi]*dense.xIonDense[nelem][nelem]*
00311 helike.BranchRatio[nelem][ipHi][ipLo]*
00312 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg,
00313 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].WLAng,
00314 "+Col",
00315 'i');
00316 }
00317 }
00318
00319 else if( ipLo==ipHe2p3P1 && ipHi==ipHe4d3D && nMatch("Ca B",CaABLines[nelem][i].label) )
00320 {
00321 double intens = TempInterp( CaABTemps , CaABIntensity[nelem][i], NUMTEMPS )*
00322 dense.xIonDense[nelem][nelem]*dense.eden;
00323
00324 linadd( intens, 4471, CaABLines[nelem][i].label, 'i');
00325 }
00326
00327 }
00328 free( qTotEff );
00329 }
00330
00331
00332
00333
00334
00335
00336 EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].xIntensity =
00337 (double)EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].Aul*
00338 (double)EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].PopHi*
00339 (double)EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].Pesc*
00340 (double)dense.xIonDense[nelem][nelem]*
00341 (double)EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].EnergyErg;
00342 if( LineSave.ipass == 0 )
00343 {
00344
00345
00346
00347 chIonLbl(chLabel, &EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S]);
00348 }
00349
00350 linadd( EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].xIntensity , 0,chLabel,'r');
00351
00352
00353 linadd(
00354 iso.Pop2Ion[ipHE_LIKE][nelem][ipHe2s1S]*
00355 dense.xIonDense[nelem][nelem]*
00356 iso.TwoNu_induc_dn[ipHE_LIKE][nelem]*
00357 EmisLines[ipHE_LIKE][nelem][ipHe2s1S][ipHe1s1S].EnergyErg,
00358 22, chLabel ,'i');
00359
00360
00361
00362
00363 sum = 0.;
00364 for( i=ipHe2p3P0; i <= ipHe2p3P2; i++ )
00365 {
00366 sum +=
00367 (double)EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].Aul*
00368 (double)EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].PopHi*
00369 (double)EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].Pesc*
00370 (double)dense.xIonDense[nelem][nelem]*
00371 (double)EmisLines[ipHE_LIKE][nelem][i][ipHe1s1S].EnergyErg;
00372 }
00373
00374
00375 linadd(sum,EmisLines[ipHE_LIKE][nelem][ipHe2p3P1][ipHe1s1S].WLAng,"TOTL",'i' );
00376
00377
00378 for( ipLo=0; ipLo < ipHe2p3P0; ipLo++ )
00379 {
00380 for( ipHi=ipLo+1; ipHi < iso.numPrintLevels[ipHE_LIKE][nelem]; ipHi++ )
00381 {
00382
00383
00384 if( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ipCont < 1 )
00385 continue;
00386
00387
00388
00389 if( helike.lgFSM && ( abs(iso.quant_desig[ipHE_LIKE][nelem][ipHi].l -
00390 iso.quant_desig[ipHE_LIKE][nelem][ipLo].l)==1 )
00391 && (iso.quant_desig[ipHE_LIKE][nelem][ipLo].l>1)
00392 && (iso.quant_desig[ipHE_LIKE][nelem][ipHi].l>1)
00393 && ( iso.quant_desig[ipHE_LIKE][nelem][ipHi].n ==
00394 iso.quant_desig[ipHE_LIKE][nelem][ipLo].n ) )
00395 {
00396 if( (iso.quant_desig[ipHE_LIKE][nelem][ipHi].s==0)
00397 && (iso.quant_desig[ipHE_LIKE][nelem][ipLo].s==0) )
00398 {
00399 continue;
00400 }
00401 else if( (iso.quant_desig[ipHE_LIKE][nelem][ipHi].s==1)
00402 && (iso.quant_desig[ipHE_LIKE][nelem][ipLo].s==1) )
00403 {
00404
00405
00406 EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo+1].phots =
00407 ( (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo+1].Aul*
00408 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo+1].PopHi*
00409 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo+1].Pesc +
00410 (double)EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo+1].Aul*
00411 (double)EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo+1].PopHi*
00412 (double)EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo+1].Pesc )*
00413 (double)dense.xIonDense[nelem][nelem];
00414
00415 EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo+1].xIntensity =
00416 EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo+1].phots *
00417 (double)EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo+1].EnergyErg;
00418
00419 PutLine(&EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo+1]);
00420
00421
00422 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots =
00423 ( (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
00424 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
00425 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc +
00426 (double)EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo].Aul*
00427 (double)EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo].PopHi*
00428 (double)EmisLines[ipHE_LIKE][nelem][ipHi+1][ipLo].Pesc )*
00429 (double)dense.xIonDense[nelem][nelem];
00430
00431 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity =
00432 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots *
00433 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;
00434
00435 PutLine(&EmisLines[ipHE_LIKE][nelem][ipHi][ipLo]);
00436 }
00437 }
00438
00439 else if( ipLo==ipHe2s3S && ipHi == ipHe2p3P0 )
00440 {
00441
00442
00443
00444
00445
00446 float av_wl = 0.;
00447 sum = 0.;
00448 for( i=ipHe2p3P0; i <= ipHe2p3P2; i++ )
00449 {
00450 sum +=
00451 (double)EmisLines[ipHE_LIKE][nelem][i][ipLo].Aul*
00452 (double)EmisLines[ipHE_LIKE][nelem][i][ipLo].PopHi*
00453 (double)EmisLines[ipHE_LIKE][nelem][i][ipLo].Pesc*
00454 (double)dense.xIonDense[nelem][nelem]*
00455 (double)EmisLines[ipHE_LIKE][nelem][i][ipLo].EnergyErg;
00456 av_wl += EmisLines[ipHE_LIKE][nelem][i][ipLo].WLAng;
00457 }
00458 av_wl /= 3.;
00459 # if 0
00460 {
00461 # include "elementnames.h"
00462 # include "prt.h"
00463 fprintf(ioQQQ,"DEBUG 2P - 2S multiplet wl %s ",
00464 elementnames.chElementSym[nelem] );
00465 prt_wl( ioQQQ , av_wl );
00466 fprintf(ioQQQ,"\n" );
00467 }
00468 # endif
00469
00470
00471
00472 linadd(sum,av_wl,"TOTL",'i' );
00473 # ifdef PRT_HELIKE
00474 # endif
00475
00476
00477
00478
00479
00480 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots =
00481 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
00482 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
00483 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc*
00484 dense.xIonDense[nelem][nelem];
00485
00486
00487 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity =
00488 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
00489 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
00490 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc*
00491 (double)dense.xIonDense[nelem][nelem]*
00492 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;
00493
00494 if( helike.lgRandErrGen )
00495 {
00496 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots *=
00497 helike.ErrorFactor[nelem][ipHi][ipLo][IPRAD];
00498 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity *=
00499 helike.ErrorFactor[nelem][ipHi][ipLo][IPRAD];
00500 }
00501 PutLine(&EmisLines[ipHE_LIKE][nelem][ipHi][ipLo]);
00502 }
00503
00504 else
00505 {
00506
00507
00508 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots =
00509 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
00510 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
00511 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc*
00512 dense.xIonDense[nelem][nelem];
00513
00514
00515
00516 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity =
00517 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
00518 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
00519 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc*
00520 (double)dense.xIonDense[nelem][nelem]*
00521 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;
00522
00523 if( helike.lgRandErrGen )
00524 {
00525 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots *=
00526 helike.ErrorFactor[nelem][ipHi][ipLo][IPRAD];
00527 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity *=
00528 helike.ErrorFactor[nelem][ipHi][ipLo][IPRAD];
00529 }
00530
00531
00532
00533
00534 PutLine(&EmisLines[ipHE_LIKE][nelem][ipHi][ipLo]);
00535 {
00536
00537
00538
00539 enum {DEBUG_LOC=false};
00540
00541 if( DEBUG_LOC )
00542 {
00543 if( nelem==1 && ipLo==0 && ipHi==1 )
00544 {
00545 fprintf(ioQQQ,"he1 626 %.2e %.2e \n",
00546 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].TauIn,
00547 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].TauTot
00548 );
00549 }
00550 }
00551 }
00552 }
00553 }
00554 }
00555
00556
00557 for( ipHi=ipHe2p3P2+1; ipHi < iso.numPrintLevels[ipHE_LIKE][nelem]; ipHi++ )
00558 {
00559 double sumcool , sumheat ,
00560 save , savecool , saveheat;
00561
00562 sum = 0;
00563 sumcool = 0.;
00564 sumheat = 0.;
00565 for( ipLo=ipHe2p3P0; ipLo <= ipHe2p3P2; ++ipLo )
00566 {
00567
00568 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots =
00569 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
00570 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
00571 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc*
00572 dense.xIonDense[nelem][nelem];
00573
00574
00575
00576 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity =
00577 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
00578 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
00579 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc*
00580 (double)dense.xIonDense[nelem][nelem]*
00581 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;
00582
00583 if( helike.lgRandErrGen )
00584 {
00585 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots *=
00586 helike.ErrorFactor[nelem][ipHi][ipLo][IPRAD];
00587 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity *=
00588 helike.ErrorFactor[nelem][ipHi][ipLo][IPRAD];
00589 }
00590
00591 sumcool += EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].cool;
00592 sumheat += EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].heat;
00593 sum += EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity;
00594 }
00595
00596
00597
00598 if( EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].ipCont < 1 )
00599 continue;
00600
00601 save = EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].xIntensity;
00602 savecool = EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].cool;
00603 saveheat = EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].heat;
00604
00605 EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].xIntensity = sum;
00606 EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].cool = sumcool;
00607 EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].heat = sumheat;
00608
00609
00610
00611 PutLine(&EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2]);
00612
00613 EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].xIntensity = save;
00614 EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].cool = savecool;
00615 EmisLines[ipHE_LIKE][nelem][ipHi][ipHe2p3P2].heat = saveheat;
00616 }
00617 for( ipLo=ipHe2p3P2+1; ipLo < iso.numPrintLevels[ipHE_LIKE][nelem]-1; ipLo++ )
00618 {
00619 for( ipHi=ipLo+1; ipHi < iso.numPrintLevels[ipHE_LIKE][nelem]; ipHi++ )
00620 {
00621
00622
00623
00624
00625 if( EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ipCont < 1 )
00626 continue;
00627
00628
00629 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots =
00630 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
00631 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
00632 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc*
00633 dense.xIonDense[nelem][nelem];
00634
00635
00636
00637 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity =
00638 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Aul*
00639 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi*
00640 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].Pesc*
00641 (double)dense.xIonDense[nelem][nelem]*
00642 (double)EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].EnergyErg;
00643
00644 if( helike.lgRandErrGen )
00645 {
00646 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].phots *=
00647 helike.ErrorFactor[nelem][ipHi][ipLo][IPRAD];
00648 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].xIntensity *=
00649 helike.ErrorFactor[nelem][ipHi][ipLo][IPRAD];
00650 }
00651
00652
00653
00654
00655 PutLine(&EmisLines[ipHE_LIKE][nelem][ipHi][ipLo]);
00656 }
00657 }
00658 }
00659 }
00660
00661 if( iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM] >= 4 &&
00662 iso.n_HighestResolved_max[ipH_LIKE][ipHYDROGEN] >=8 )
00663 {
00664
00665
00666 photons_3889_plus_7065 =
00667
00668 EmisLines[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P2].xIntensity/
00669 EmisLines[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P2].EnergyErg +
00670 EmisLines[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P2].xIntensity/
00671 EmisLines[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P2].EnergyErg +
00672 EmisLines[ipHE_LIKE][ipHELIUM][ipHe4s3S][ipHe2p3P2].xIntensity/
00673 EmisLines[ipHE_LIKE][ipHELIUM][ipHe4s3S][ipHe2p3P2].EnergyErg +
00674
00675 EmisLines[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P1].xIntensity/
00676 EmisLines[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P1].EnergyErg +
00677 EmisLines[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P1].xIntensity/
00678 EmisLines[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P1].EnergyErg +
00679 EmisLines[ipHE_LIKE][ipHELIUM][ipHe4s3S][ipHe2p3P1].xIntensity/
00680 EmisLines[ipHE_LIKE][ipHELIUM][ipHe4s3S][ipHe2p3P1].EnergyErg +
00681
00682 EmisLines[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P0].xIntensity/
00683 EmisLines[ipHE_LIKE][ipHELIUM][ipHe3s3S][ipHe2p3P0].EnergyErg +
00684 EmisLines[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P0].xIntensity/
00685 EmisLines[ipHE_LIKE][ipHELIUM][ipHe3d3D][ipHe2p3P0].EnergyErg +
00686 EmisLines[ipHE_LIKE][ipHELIUM][ipHe4s3S][ipHe2p3P0].xIntensity/
00687 EmisLines[ipHE_LIKE][ipHELIUM][ipHe4s3S][ipHe2p3P0].EnergyErg +
00688
00689 EmisLines[ipHE_LIKE][ipHELIUM][ipHe3p3P][ipHe2s3S].xIntensity/
00690 EmisLines[ipHE_LIKE][ipHELIUM][ipHe3p3P][ipHe2s3S].EnergyErg +
00691 EmisLines[ipHE_LIKE][ipHELIUM][ipHe4p3P][ipHe2s3S].xIntensity/
00692 EmisLines[ipHE_LIKE][ipHELIUM][ipHe4p3P][ipHe2s3S].EnergyErg;
00693
00694
00695 photons_3889_plus_7065 +=
00696 EmisLines[ipH_LIKE][ipHYDROGEN][8][1].xIntensity/
00697 EmisLines[ipH_LIKE][ipHYDROGEN][8][1].EnergyErg;
00698
00699
00700 photons_3889_plus_7065 *= (ERG1CM*1.e8)/LineSave.WavLNorm;
00701 linadd( photons_3889_plus_7065, 3889., "Pho+", 'i');
00702 }
00703
00704
00705
00706
00707
00708 #if 0
00709
00710
00711 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
00712 {
00716 if( nelem != ipHELIUM )
00717 continue;
00718
00719
00720 if( nelem == ipHELIUM || dense.lgElmtOn[nelem])
00721 {
00722 for( j = 0; j < NumLines; ++j )
00723 {
00724 free( CaABIntensity[nelem][j] );
00725 }
00726 free( CaABIntensity[nelem] );
00727 free( CaABLines[nelem] );
00728 }
00729 }
00730 free( CaABIntensity );
00731 free( CaABLines );
00732 #endif
00733
00734 if( trace.lgTrace )
00735 {
00736 fprintf( ioQQQ, " lines_helium returns\n" );
00737 }
00738
00739 DEBUG_EXIT( "lines_helium()" );
00740 return;
00741 }
00742
00743
00744 static void GetStandardHeLines(void)
00745 {
00746 FILE *ioDATA;
00747 bool lgEOL, lgHIT;
00748 long i, i1, i2, j, nelem;
00749
00750 # define chLine_LENGTH 1000
00751 char chLine[chLine_LENGTH] ,
00752
00753 chFilename[FILENAME_PATH_LENGTH_2];
00754
00755
00756
00757 if( lgDataPathSet == true )
00758 {
00759
00760 strcpy( chFilename , chDataPath );
00761 strcat( chFilename , "he1_case_ab.dat" );
00762 }
00763 else
00764 {
00765
00766 strcpy( chFilename , "he1_case_ab.dat" );
00767 }
00768
00769 if( trace.lgTrace )
00770 fprintf( ioQQQ," lines_helium opening he1_case_ab.dat:");
00771
00772 if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
00773 {
00774 fprintf( ioQQQ, " lines_helium could not open he1_case_ab.dat\n" );
00775 if( lgDataPathSet == true )
00776 fprintf( ioQQQ, " even tried path\n" );
00777
00778 if( lgDataPathSet == true )
00779 {
00780 fprintf( ioQQQ, " lines_helium could not open he1_case_ab.dat\n");
00781 fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
00782 fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
00783 }
00784
00785 puts( "[Stop in lines_helium]" );
00786 cdEXIT(EXIT_FAILURE);
00787 }
00788
00789
00790 if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
00791 {
00792 fprintf( ioQQQ, " lines_helium could not read first line of he1_case_ab.dat.\n");
00793 puts( "[Stop in lines_helium]" );
00794 cdEXIT(EXIT_FAILURE);
00795 }
00796 i = 1;
00797
00798 i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00799 i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00800 NumLines = i2;
00801
00802
00803 if( i1 !=CASEABMAGIC )
00804 {
00805 fprintf( ioQQQ,
00806 " lines_helium: the version of he1_case_ab.dat is not the current version.\n" );
00807 fprintf( ioQQQ,
00808 " lines_helium: I expected to find the number %i and got %li instead.\n" ,
00809 CASEABMAGIC, i1 );
00810 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
00811 puts( "[Stop in lines_helium]" );
00812 cdEXIT(EXIT_FAILURE);
00813 }
00814
00815
00816 lgHIT = false;
00817 while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00818 {
00819
00820 if( chLine[0] != '#')
00821 {
00822 lgHIT = true;
00823 j = 0;
00824 lgEOL = false;
00825 i = 1;
00826 while( !lgEOL && j < NUMTEMPS)
00827 {
00828 CaABTemps[j] = FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
00829 ++j;
00830 }
00831 break;
00832 }
00833 }
00834
00835 if( !lgHIT )
00836 {
00837 fprintf( ioQQQ, " lines_helium could not find line of temperatures in he1_case_ab.dat.\n");
00838 puts( "[Stop in lines_helium]" );
00839 cdEXIT(EXIT_FAILURE);
00840 }
00841
00842
00843 CaABIntensity = (double ***)MALLOC(sizeof(double **)*(unsigned)LIMELM );
00844
00845 CaABLines = (stdLines **)MALLOC(sizeof(stdLines *)*(unsigned)LIMELM );
00846
00847 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
00848 {
00852 if( nelem != ipHELIUM )
00853 continue;
00854
00855
00856 if( nelem == ipHELIUM || dense.lgElmtOn[nelem])
00857 {
00858
00859 CaABIntensity[nelem] = (double **)MALLOC(sizeof(double *)*(unsigned)(i2) );
00860 CaABLines[nelem] = (stdLines *)MALLOC(sizeof(stdLines )*(unsigned)(i2) );
00861
00862
00863
00864 for( j = 0; j < i2; ++j )
00865 {
00866 CaABIntensity[nelem][j] = (double *)MALLOC(sizeof(double)*(unsigned)NUMTEMPS );
00867
00868 CaABLines[nelem][j].ipHi = -1;
00869 CaABLines[nelem][j].ipLo = -1;
00870 strcpy( CaABLines[nelem][j].label , " " );
00871
00872 for( i=0; i<NUMTEMPS; ++i )
00873 {
00874 CaABIntensity[nelem][j][i] = 0.;
00875 }
00876 }
00877 }
00878 }
00879
00880
00881 lgHIT = false;
00882 nelem = ipHELIUM;
00883 while( fgets( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00884 {
00885 static long line = 0;
00886 char *chTemp;
00887
00888
00889 if( (chLine[0] == ' ') || (chLine[0]=='\n') )
00890 break;
00891 if( chLine[0] != '#')
00892 {
00893 lgHIT = true;
00894
00895
00896 j = 1;
00897
00898
00899 i1 = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00900 CaABLines[nelem][line].ipLo = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00901 CaABLines[nelem][line].ipHi = (long)FFmtRead(chLine,&j,INPUT_LINE_LENGTH,&lgEOL);
00902
00903 ASSERT( CaABLines[nelem][line].ipLo < CaABLines[nelem][line].ipHi );
00904
00905
00906
00907 chTemp = chLine;
00908
00909 for( i=0; i<3; ++i )
00910 {
00911 if( (chTemp = strchr( chTemp, '\t' )) == NULL )
00912 {
00913 fprintf( ioQQQ, " lines_helium no init case A and B\n" );
00914 puts( "[Stop in lines_helium]" );
00915 cdEXIT(EXIT_FAILURE);
00916 }
00917 ++chTemp;
00918 }
00919
00920 strncpy( CaABLines[nelem][line].label, chTemp , 4 );
00921 CaABLines[nelem][line].label[4] = 0;
00922
00923 for( i=0; i<NUMTEMPS; ++i )
00924 {
00925 float b;
00926 if( (chTemp = strchr( chTemp, '\t' )) == NULL )
00927 {
00928 fprintf( ioQQQ, " lines_helium could not scan case A and B lines, current indices: %li %li\n",
00929 CaABLines[nelem][line].ipHi,
00930 CaABLines[nelem][line].ipLo );
00931 puts( "[Stop in lines_helium]" );
00932 cdEXIT(EXIT_FAILURE);
00933 }
00934 ++chTemp;
00935 sscanf( chTemp , "%e" , &b );
00936 CaABIntensity[nelem][line][i] = pow(10.,(double)b);
00937 }
00938 line++;
00939 }
00940 }
00941
00942
00943 fclose( ioDATA );
00944
00945 return;
00946 }
00947
00949 static double TempInterp( double* TempArray , double* ValueArray, long NumElements )
00950 {
00951 long int ipTe=-1;
00952 double rate = 0.;
00953 long i0;
00954
00955 DEBUG_ENTRY( "TempInterp()" );
00956
00957 ASSERT( fabs( 1. - (double)phycon.alogte/log10(phycon.te) ) < 0.0001 );
00958
00959
00960 if( phycon.alogte <= TempArray[0] )
00961 {
00962 return ValueArray[0];
00963 }
00964 else if( phycon.alogte >= TempArray[NumElements-1] )
00965 {
00966 return ValueArray[NumElements-1];
00967 }
00968
00969
00970 ipTe = hunt_bisect( TempArray, NumElements, phycon.alogte );
00971
00972 ASSERT( (ipTe >=0) && (ipTe < NumElements-1) );
00973
00974 #if 0
00975 rate = ((double)phycon.alogte - TempArray[ipTe]) / (TempArray[ipTe+1] - TempArray[ipTe]) *
00976 (ValueArray[ipTe+1] - ValueArray[ipTe]) + ValueArray[ipTe];
00977 #endif
00978
00979
00980 const int ORDER = 3;
00981 i0 = max(min(ipTe-ORDER/2,NumElements-ORDER-1),0);
00982 rate = lagrange( &TempArray[i0], &ValueArray[i0], ORDER+1, phycon.alogte );
00983
00984 DEBUG_EXIT( "TempInterp()" );
00985
00986 return rate;
00987 }