00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #include "cddefines.h"
00013
00014 #ifdef __DECC
00015 # pragma message disable (BADBOUNDCHK)
00016 #endif
00017 #include "physconst.h"
00018 #include "taulines.h"
00019 #include "iso.h"
00020 #include "ionbal.h"
00021 #include "trace.h"
00022 #include "dense.h"
00023 #include "thermal.h"
00024 #include "phycon.h"
00025 #include "rfield.h"
00026 #include "heavy.h"
00027 #include "secondaries.h"
00028 #include "opacity.h"
00029 #include "conv.h"
00030 #include "hydro_vs_rates.h"
00031 #include "atmdat.h"
00032 #include "hydrogenic.h"
00033
00034
00035 static double Fe26cs123(long int i,
00036 long int j);
00037
00038
00039 static double He2cs123(long int i,
00040 long int j);
00041
00042
00043 static double H1cs123(long int i,
00044 long int j,
00045 long int chType);
00046
00047
00048 static double Hydcs123(long int ilow,
00049 long int ihigh,
00050 long int iz,
00051 long int chType);
00052
00053
00054 static double C6cs123(long int i,
00055 long int j);
00056
00057
00058 static double Ca20cs123(long int i,
00059 long int j);
00060
00061
00062 static float HCSAR_interp( int ipLo , int ipHi );
00063
00064
00065 static double Ne10cs123(long int i,
00066 long int j);
00067
00068
00069
00070 static const int NHCSTE = 8;
00071 static const int NHCS = 6;
00072 static const float HCSTE[NHCSTE] = {5802.f,11604.f,34812.f,58020.f,116040.f,174060.f,232080.f,290100.f};
00073
00074
00075 #if 0
00076 static const float HCS[NHCSTE*NHCS*(NHCS-1)] = {
00077 0., 0., 0., 0., 0., 0., 0., 0.,
00078 2.600e-01f, 2.960e-01f, 3.250e-01f, 3.370e-01f, 3.560e-01f, 3.680e-01f, 3.750e-01f, 3.800e-01f,
00079 4.270e-01f, 5.360e-01f, 8.570e-01f, 1.150e+00f, 1.750e+00f, 2.130e+00f, 2.350e+00f, 2.460e+00f,
00080 2.372e-01f, 2.605e-01f, 3.413e-01f, 4.104e-01f, 5.063e-01f, 5.351e-01f, 5.335e-01f, 5.236e-01f,
00081 9.990e-02f, 1.144e-01f, 1.564e-01f, 1.857e-01f, 2.263e-01f, 2.422e-01f, 2.445e-01f, 2.428e-01f,
00082 7.325e-02f, 8.263e-02f, 9.641e-02f, 1.064e-01f, 1.202e-01f, 1.233e-01f, 1.219e-01f, 1.191e-01f,
00083 0., 0., 0., 0., 0., 0., 0., 0.,
00084 0., 0., 0., 0., 0., 0., 0., 0.,
00085 0., 0., 0., 0., 0., 0., 0., 0.,
00086 5.860e+00f, 7.640e+00f, 1.419e+01f, 2.015e+01f, 3.156e+01f, 3.969e+01f, 4.598e+01f, 5.114e+01f,
00087 2.276e+00f, 2.747e+00f, 4.365e+00f, 5.609e+00f, 7.803e+00f, 9.267e+00f, 1.038e+01f, 1.126e+01f,
00088 1.610e+00f, 1.897e+00f, 2.421e+00f, 2.799e+00f, 3.512e+00f, 4.023e+00f, 4.450e+00f, 4.790e+00f,
00089 0., 0., 0., 0., 0., 0., 0., 0.,
00090 0., 0., 0., 0., 0., 0., 0., 0.,
00091 0., 0., 0., 0., 0., 0., 0., 0.,
00092 2.228e+01f, 2.818e+01f, 5.026e+01f, 7.193e+01f, 1.178e+02f, 1.542e+02f, 1.828e+02f, 2.078e+02f,
00093 8.772e+00f, 1.042e+01f, 1.600e+01f, 2.058e+01f, 2.898e+01f, 3.482e+01f, 3.923e+01f, 4.280e+01f,
00094 6.161e+00f, 7.209e+00f, 9.091e+00f, 1.049e+01f, 1.308e+01f, 1.488e+01f, 1.626e+01f, 1.736e+01f,
00095 0., 0., 0., 0., 0., 0., 0., 0.,
00096 0., 0., 0., 0., 0., 0., 0., 0.,
00097 0., 0., 0., 0., 0., 0., 0., 0.,
00098 0., 0., 0., 0., 0., 0., 0., 0.,
00099 8.648e+01f, 1.366e+02f, 3.252e+02f, 4.834e+02f, 7.715e+02f, 9.692e+02f, 1.123e+03f, 1.247e+03f,
00100 6.229e+01f, 8.263e+01f, 1.255e+02f, 1.527e+02f, 1.951e+02f, 2.211e+02f, 2.400e+02f, 2.544e+02f,
00101 0., 0., 0., 0., 0., 0., 0., 0.,
00102 0., 0., 0., 0., 0., 0., 0., 0.,
00103 0., 0., 0., 0., 0., 0., 0., 0.,
00104 0., 0., 0., 0., 0., 0., 0., 0.,
00105 0., 0., 0., 0., 0., 0., 0., 0.,
00106 4.625e+02f, 8.070e+02f, 1.707e+03f, 2.269e+03f, 3.122e+03f, 3.675e+03f, 4.083e+03f, 4.414e+03f
00107 };
00108 #endif
00109
00110
00111
00112
00113 static const float HCS[NHCSTE*NHCS*(NHCS-1)] = {
00114 0., 0., 0., 0., 0., 0., 0., 0.,
00115 2.600e-01f, 2.960e-01f, 3.260e-01f, 3.390e-01f, 3.730e-01f, 4.060e-01f, 4.360e-01f, 4.610e-01f,
00116 4.290e-01f, 5.290e-01f, 8.530e-01f, 1.150e+00f, 1.810e+00f, 2.350e+00f, 2.810e+00f, 3.200e+00f,
00117 2.372e-01f, 2.605e-01f, 3.413e-01f, 4.140e-01f, 5.500e-01f, 6.510e-01f, 7.300e-01f, 7.976e-01f,
00118 9.990e-02f, 1.144e-01f, 1.564e-01f, 1.857e-01f, 2.263e-01f, 2.422e-01f, 2.445e-01f, 2.428e-01f,
00119 7.325e-02f, 8.263e-02f, 9.641e-02f, 1.064e-01f, 1.202e-01f, 1.233e-01f, 1.219e-01f, 1.191e-01f,
00120 0., 0., 0., 0., 0., 0., 0., 0.,
00121 0., 0., 0., 0., 0., 0., 0., 0.,
00122 0., 0., 0., 0., 0., 0., 0., 0.,
00123 5.860e+00f, 7.640e+00f, 1.419e+01f, 2.015e+01f, 3.156e+01f, 3.969e+01f, 4.598e+01f, 5.114e+01f,
00124 2.276e+00f, 2.747e+00f, 4.365e+00f, 5.609e+00f, 7.803e+00f, 9.267e+00f, 1.038e+01f, 1.126e+01f,
00125 1.610e+00f, 1.897e+00f, 2.421e+00f, 2.799e+00f, 3.512e+00f, 4.023e+00f, 4.450e+00f, 4.790e+00f,
00126 0., 0., 0., 0., 0., 0., 0., 0.,
00127 0., 0., 0., 0., 0., 0., 0., 0.,
00128 0., 0., 0., 0., 0., 0., 0., 0.,
00129 2.228e+01f, 2.818e+01f, 5.026e+01f, 7.193e+01f, 1.178e+02f, 1.542e+02f, 1.828e+02f, 2.078e+02f,
00130 8.772e+00f, 1.042e+01f, 1.600e+01f, 2.058e+01f, 2.898e+01f, 3.482e+01f, 3.923e+01f, 4.280e+01f,
00131 6.161e+00f, 7.209e+00f, 9.091e+00f, 1.049e+01f, 1.308e+01f, 1.488e+01f, 1.626e+01f, 1.736e+01f,
00132 0., 0., 0., 0., 0., 0., 0., 0.,
00133 0., 0., 0., 0., 0., 0., 0., 0.,
00134 0., 0., 0., 0., 0., 0., 0., 0.,
00135 0., 0., 0., 0., 0., 0., 0., 0.,
00136 8.648e+01f, 1.366e+02f, 3.252e+02f, 4.834e+02f, 7.715e+02f, 9.692e+02f, 1.123e+03f, 1.247e+03f,
00137 6.229e+01f, 8.263e+01f, 1.255e+02f, 1.527e+02f, 1.951e+02f, 2.211e+02f, 2.400e+02f, 2.544e+02f,
00138 0., 0., 0., 0., 0., 0., 0., 0.,
00139 0., 0., 0., 0., 0., 0., 0., 0.,
00140 0., 0., 0., 0., 0., 0., 0., 0.,
00141 0., 0., 0., 0., 0., 0., 0., 0.,
00142 0., 0., 0., 0., 0., 0., 0., 0.,
00143 4.625e+02f, 8.070e+02f, 1.707e+03f, 2.269e+03f, 3.122e+03f, 3.675e+03f, 4.083e+03f, 4.414e+03f
00144 };
00145
00146 void HydroCollid(
00147
00148 long int nelem)
00149 {
00150 long int i,
00151 ipHi,
00152 ipLo,
00153 n;
00154 double CStemp,
00155 factor ,
00156 ConvLTEPOP;
00157 double collider;
00158
00159 static float told[LIMELM] = {0.};
00160
00161
00162 DEBUG_ENTRY( "HydroCollid()" );
00163
00164
00165 ASSERT( nelem >= 0);
00166 ASSERT( nelem < LIMELM );
00167
00168
00169
00170
00171
00172
00173 if( nelem==ipHYDROGEN )
00174 {
00175
00176 collider = dense.EdenHontoHCorr;
00177 }
00178 else
00179 {
00180 collider = dense.EdenHCorr;
00181 }
00182
00183
00184
00185
00186 {
00187
00188 enum {AGN=false};
00189
00190 if( AGN )
00191 {
00192 static const double tear[3] = {10000.,15000.,20000.};
00193 for( i=0;i<3;++i)
00194 {
00195 phycon.te = tear[i];
00196 tfidle(false);
00197 for( ipHi=1;ipHi<4;++ipHi )
00198 {
00199 for( ipLo=0; ipLo<ipHi; ++ipLo )
00200 {
00201 if( !(ipLo==1 && ipHi==2 ) )
00202 {
00203 fprintf(ioQQQ,"%.3f\t", HCSAR_interp(ipLo,ipHi) );
00204 }
00205 }
00206 }
00207 fprintf(ioQQQ," \n");
00208 }
00209 cdEXIT(EXIT_SUCCESS);
00210 }
00211 }
00212
00213
00214
00215
00216
00217
00218
00219 if( nelem == 0 )
00220 {
00221
00222
00223 secondaries.Hx12[0][ipH2p] = (float)(secondaries.x12tot*2./3./1.261);
00224 secondaries.Hx12[0][ipH2s] = (float)(secondaries.x12tot*1./3./1.261);
00225 secondaries.Hx12[0][3] = (float)(secondaries.x12tot*0.165/1.261);
00226 secondaries.Hx12[0][4] = (float)(secondaries.x12tot*0.059/1.261);
00227
00228
00229 for( i=5; i < iso.numLevels_max[ipH_LIKE][nelem]; i++ )
00230 {
00231 secondaries.Hx12[0][i] = (float)(secondaries.Hx12[0][i-1]/2.);
00232 }
00233 }
00234 else if( nelem==1 )
00235 {
00236
00237
00238
00239
00240
00241 secondaries.Hx12[1][ipH2p] = (float)(secondaries.x12tot/3.*2./3.);
00242 secondaries.Hx12[1][ipH2s] = (float)(secondaries.x12tot/3./3.);
00243
00244 for( i=3; i < iso.numLevels_max[ipH_LIKE][nelem]; i++ )
00245 {
00246 secondaries.Hx12[1][i] = (float)(secondaries.Hx12[1][i-1]/2.);
00247 }
00248 }
00249
00250
00251
00252
00253
00254 if( fabs(phycon.te-told[nelem])/phycon.te < 0.00002 &&
00255 EmisLines[ipH_LIKE][nelem][ipH2p][ipH1s].ColUL > 0. && conv.nTotalIoniz )
00256 {
00257 if( trace.lgTrace )
00258 {
00259 fprintf( ioQQQ,
00260 " HydroCollid called nelem %li - no reeval Boltz fac, LTE dens\n",
00261 nelem );
00262 }
00263
00264
00265 }
00266 else
00267 {
00268 if( trace.lgTrace )
00269 {
00270 fprintf( ioQQQ,
00271 " HydroCollid called nelem %li - will reeval Boltz fac, LTE dens\n",
00272 nelem );
00273 }
00274 told[nelem] = (float)phycon.te;
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284 for( ipHi=ipH2s; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00285 {
00286
00287
00288
00289
00290 iso.Boltzmann[ipH_LIKE][nelem][ipHi][ipHi-1] =
00291 exp(-(double)EmisLines[ipH_LIKE][nelem][ipHi][ipHi-1].EnergyK/phycon.te);
00292 }
00293
00294
00295
00296 iso.ConBoltz[ipH_LIKE][nelem][iso.numLevels_max[ipH_LIKE][nelem]-1] =
00297 exp(-iso.xIsoLevNIonRyd[ipH_LIKE][nelem][iso.numLevels_max[ipH_LIKE][nelem]-1]/phycon.te_ryd);
00298
00299
00300 for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 2); ipLo++ )
00301 {
00302 for( ipHi=ipLo + 2; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00303 {
00304 iso.Boltzmann[ipH_LIKE][nelem][ipHi][ipLo] =
00305 iso.Boltzmann[ipH_LIKE][nelem][ipHi-1][ipLo]*
00306 iso.Boltzmann[ipH_LIKE][nelem][ipHi][ipHi-1];
00307 }
00308 }
00309
00310 for( i=ipH1s; i < (iso.numLevels_max[ipH_LIKE][nelem] - 1); i++ )
00311 {
00312 iso.ConBoltz[ipH_LIKE][nelem][i] =
00313 iso.Boltzmann[ipH_LIKE][nelem][iso.numLevels_max[ipH_LIKE][nelem]-1][i]*
00314 iso.ConBoltz[ipH_LIKE][nelem][iso.numLevels_max[ipH_LIKE][nelem]-1];
00315
00316 if( iso.ConBoltz[ipH_LIKE][nelem][i] < SMALLDOUBLE )
00317 {
00318 iso.ConBoltz[ipH_LIKE][nelem][i] = 0.;
00319 }
00320 }
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335 factor = HION_LTE_POP*dense.AtomicWeight[nelem]/
00336 (dense.AtomicWeight[nelem]+ELECTRON_MASS/ATOMIC_MASS_UNIT);
00337
00338
00339 ConvLTEPOP = pow(factor,1.5)/(2.*iso.stat_ion[ipH_LIKE])/phycon.te32;
00340
00341 for( n=ipH1s; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00342 {
00343
00344
00345 if( iso.ConBoltz[ipH_LIKE][nelem][n] > 1e-100 )
00346 {
00347 iso.PopLTE[ipH_LIKE][nelem][n] =
00348 iso.stat[ipH_LIKE][nelem][n] / iso.ConBoltz[ipH_LIKE][nelem][n]* ConvLTEPOP;
00349 }
00350 else
00351 {
00352 iso.PopLTE[ipH_LIKE][nelem][n] = 0.;
00353 }
00354 }
00355
00356
00357 iso.lgPopLTE_OK[ipH_LIKE][nelem] = true;
00358 for( n=ipH1s; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00359 {
00360 if( iso.PopLTE[ipH_LIKE][nelem][n] <= 0. )
00361 {
00362 iso.lgPopLTE_OK[ipH_LIKE][nelem] = false;
00363 break;
00364 }
00365 }
00366
00367
00368
00369
00370
00371
00372
00373
00374 if( nelem == ipHYDROGEN )
00375 {
00376
00377
00378
00379 iso.ColIoniz[ipH_LIKE][nelem][ipH2p] = iso.lgColl_ionize[ipH_LIKE]*
00380 hydro_vs_ioniz(ipH_LIKE, nelem, ipH2p);
00381
00382
00383 for( n=3; n < iso.numLevels_max[ipH_LIKE][nelem]-1; n++ )
00384 {
00385
00386 iso.ColIoniz[ipH_LIKE][nelem][n] = iso.lgColl_ionize[ipH_LIKE] *
00387 hydro_vs_ioniz(ipH_LIKE, nelem, n);
00388 }
00389
00390 iso.ColIoniz[ipH_LIKE][nelem][iso.numLevels_max[ipH_LIKE][nelem]-1] =
00391 hydro_vs_ioniz( ipH_LIKE, nelem, iso.numLevels_max[ipH_LIKE][nelem]-1 );
00392 }
00393 else
00394 {
00395
00396
00397 iso.ColIoniz[ipH_LIKE][nelem][ipH2p] = Hion_coll_ioniz_ratecoef(ipH_LIKE,nelem,ipH2p)*iso.lgColl_ionize[ipH_LIKE];
00398 for( n=3; n < iso.numLevels_max[ipH_LIKE][nelem]-1; n++ )
00399 {
00400
00401 iso.ColIoniz[ipH_LIKE][nelem][n] = Hion_coll_ioniz_ratecoef(ipH_LIKE,nelem,n)*iso.lgColl_ionize[ipH_LIKE];
00402 }
00403
00404 iso.ColIoniz[ipH_LIKE][nelem][iso.numLevels_max[ipH_LIKE][nelem]-1] = Hion_coll_ioniz_ratecoef(ipHYDROGEN,nelem,iso.numLevels_max[ipH_LIKE][nelem]-1);
00405 }
00406
00407
00408 iso.ColIoniz[ipH_LIKE][nelem][ipH2p] *= 0.75;
00409 iso.ColIoniz[ipH_LIKE][nelem][ipH2s] = iso.ColIoniz[ipH_LIKE][nelem][ipH2p]/3.;
00410
00411
00412 iso.ColIoniz[ipH_LIKE][nelem][ipH1s] = t_ADfA::Inst().coll_ion(nelem+1,1,phycon.te);
00413
00414
00415 iso.ColIoniz[ipH_LIKE][nelem][ipH1s] *= iso.lgColl_ionize[ipH_LIKE];
00416 iso.ColIoniz[ipH_LIKE][nelem][ipH2s] *= iso.lgColl_ionize[ipH_LIKE];
00417 iso.ColIoniz[ipH_LIKE][nelem][ipH2p] *= iso.lgColl_ionize[ipH_LIKE];
00418
00419
00420 if( opac.lgCaseB_HummerStorey )
00421 {
00422 iso.ColIoniz[ipH_LIKE][nelem][ipH1s] = 0.;
00423 iso.ColIoniz[ipH_LIKE][nelem][ipH2p] = 0.;
00424 iso.ColIoniz[ipH_LIKE][nelem][ipH2s] = 0.;
00425 }
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437 for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ )
00438 {
00439 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00440 {
00441 if( ipHi == 1 )
00442 {
00443
00444 if( nelem==0 )
00445 {
00446 CStemp = HCSAR_interp(ipLo,ipHi);
00447 }
00448 else
00449 {
00450 CStemp = Hydcs123(1,2,nelem,'e');
00451 }
00452 }
00453 else if( ipHi == 2 )
00454 {
00455 if( ipLo == 0 )
00456 {
00457
00458 if( nelem==0 )
00459 {
00460 CStemp = HCSAR_interp(ipLo,ipHi);
00461 }
00462 else
00463 {
00464 CStemp = Hydcs123(1,3,nelem,'e');
00465 }
00466 }
00467
00468 else if( ipLo == 1 )
00469 {
00470
00471 if( iso.lgColl_l_mixing[ipH_LIKE] )
00472 {
00473
00474
00475
00476
00477
00478
00479 CStemp = (Hydcs123(2,3,nelem,'e')*dense.eden +
00480 Hydcs123(2,3,nelem,'p')*dense.xIonDense[ipHYDROGEN][1]) / collider;
00481 # if 0
00482
00483 CStemp = Hydcs123(2,3,nelem,'e');
00484
00485 HCol2s2pp = Hydcs123(2,3,nelem,'p');
00486
00487
00488
00489
00490
00491 CStemp += HCol2s2pp*dense.xIonDense[ipHYDROGEN][1]/collider;
00492 # endif
00493 }
00494 else
00495 {
00496 CStemp = 0.;
00497 }
00498 }
00499 else
00500 {
00501
00502 CStemp = -DBL_MAX;
00503 }
00504 }
00505
00506 else if( ipHi == 3 )
00507 {
00508
00509 if( nelem==0 )
00510 {
00511 CStemp = HCSAR_interp(ipLo,ipHi);
00512 }
00513 else
00514 {
00515 if( ipLo == 0 )
00516 {
00517
00518 CStemp = Hydcs123(1,4,nelem,'e');
00519 }
00520 else if( ipLo == 1 )
00521 {
00522
00523 CStemp = Hydcs123(2,4,nelem,'e');
00524 }
00525 else if( ipLo == 2 )
00526 {
00527
00528 CStemp = Hydcs123(3,4,nelem,'e');
00529 }
00530 else
00531 {
00532
00533 CStemp = -DBL_MAX;
00534 }
00535 }
00536 }
00537 else if( nelem==ipHYDROGEN && ipHi <= 5 )
00538 {
00539 CStemp = HCSAR_interp(ipLo,ipHi);
00540 }
00541 else
00542 {
00543
00544
00545
00546 CStemp = CS_VS80(
00547 ipH_LIKE,
00548 nelem,
00549 ipHi,
00550 ipLo,
00551 phycon.te,
00552
00553 ipELECTRON );
00554
00555
00556
00557
00558
00559
00560
00561
00562
00563
00564
00565 }
00566 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].cs = (float)CStemp;
00567 }
00568 }
00569
00570
00571
00572
00573 CStemp = EmisLines[ipH_LIKE][nelem][ipH2p][ipH2s].cs;
00574
00575
00576
00577 if( iso.lgColl_excite[ipH_LIKE] == 0 )
00578 {
00579 for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipH_LIKE][nelem]-1); ipLo++ )
00580 {
00581 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00582 {
00583 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].cs = 0.;
00584 }
00585 }
00586 }
00587
00588
00589
00590 if( opac.lgCaseB_HummerStorey )
00591 {
00592 for( ipLo=ipH1s; ipLo <= ipH2p; ipLo++ )
00593 {
00594 for( ipHi=3; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00595 {
00596 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].cs = 0.;
00597 }
00598 }
00599 }
00600
00601
00602
00603
00604 if( iso.lgColl_l_mixing[ipH_LIKE] )
00605 {
00606 EmisLines[ipH_LIKE][nelem][ipH2p][ipH2s].cs = (float)CStemp;
00607 }
00608 else
00609 {
00610 EmisLines[ipH_LIKE][nelem][ipH2p][ipH2s].cs = 0.;
00611 }
00612
00613
00614
00615
00616
00617
00618
00619 for( ipLo=ipH1s; ipLo < (iso.numLevels_max[ipH_LIKE][nelem] - 1); ipLo++ )
00620 {
00621 for( ipHi=ipLo + 1; ipHi < iso.numLevels_max[ipH_LIKE][nelem]; ipHi++ )
00622 {
00623
00624 EmisLines[ipH_LIKE][nelem][ipHi][ipLo].ColUL =
00625 (float)(EmisLines[ipH_LIKE][nelem][ipHi][ipLo].cs/
00626 phycon.sqrte*COLL_CONST/EmisLines[ipH_LIKE][nelem][ipHi][ipLo].gHi );
00627
00628 }
00629 }
00630
00631 if( (trace.lgTrace && trace.lgIsoTraceFull[ipH_LIKE]) && (nelem == trace.ipIsoTrace[ipH_LIKE]) )
00632 {
00633 fprintf( ioQQQ, " HydroCollid rate n->n-1:" );
00634 for( n=ipH1s; n < (iso.numLevels_max[ipH_LIKE][nelem] - 1); n++ )
00635 {
00636 fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n+1][n].ColUL ));
00637 }
00638 fprintf( ioQQQ, "\n" );
00639
00640 fprintf( ioQQQ, " HydroCollid col ion cof:" );
00641 for( n=ipH1s; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00642 {
00643 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.ColIoniz[ipH_LIKE][nelem][n] ));
00644 }
00645 fprintf( ioQQQ, "\n" );
00646
00647 fprintf( ioQQQ, " HydroCollid rate dn2 1s:" );
00648 for( n=ipH2s; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00649 {
00650 fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n][ipH1s].ColUL ));
00651 }
00652 fprintf( ioQQQ, "\n" );
00653
00654 fprintf( ioQQQ, " HydroCollid rate dn2 2s:" );
00655 for( n=ipH2p; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00656 {
00657 fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n][ipH2s].ColUL ));
00658 }
00659 fprintf( ioQQQ, "\n" );
00660
00661 fprintf( ioQQQ, " HydroCollid rate dn2 2p:" );
00662 for( n=3; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00663 {
00664 fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n][ipH2p].ColUL ));
00665 }
00666 fprintf( ioQQQ, "\n" );
00667
00668 fprintf( ioQQQ, " HydroCollid rate dwn2 3:" );
00669 for( n=4; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00670 {
00671 fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n][3].ColUL ));
00672 }
00673 fprintf( ioQQQ, "\n" );
00674
00675 fprintf( ioQQQ, " HydroCollid rate dwn2 4:" );
00676 for( n=5; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00677 {
00678 fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n][4].ColUL ));
00679 }
00680 fprintf( ioQQQ, "\n" );
00681
00682 fprintf( ioQQQ, " HydroCollid rate dwn2 5:" );
00683 for( n=6; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00684 {
00685 fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n][5].ColUL) );
00686 }
00687 fprintf( ioQQQ, "\n" );
00688
00689 fprintf( ioQQQ, " HydroCollid rate dwn2 6:" );
00690 for( n=7; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00691 {
00692 fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n][6].ColUL ));
00693 }
00694 fprintf( ioQQQ, "\n" );
00695
00696 fprintf( ioQQQ, " HydroCollid rate dwn2 7:" );
00697 for( n=8; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00698 {
00699 fprintf( ioQQQ,PrintEfmt("%10.3e", EmisLines[ipH_LIKE][nelem][n][7].ColUL) );
00700 }
00701 fprintf( ioQQQ, "\n" );
00702
00703 fprintf( ioQQQ, " HydroCollid boltz fac c:" );
00704 for( n=ipH1s; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00705 {
00706 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.ConBoltz[ipH_LIKE][nelem][n] ));
00707 }
00708 fprintf( ioQQQ, "\n" );
00709
00710 fprintf( ioQQQ, " HydroCollid HLTE(n,Z): " );
00711 for( n=ipH1s; n < iso.numLevels_max[ipH_LIKE][nelem]; n++ )
00712 {
00713 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.PopLTE[ipH_LIKE][nelem][n] ));
00714 }
00715 fprintf( ioQQQ, "\n" );
00716 }
00717 }
00718
00719
00720
00721 ionbal.CollIonRate_Ground[nelem][nelem][0] =
00722 (iso.ColIoniz[ipH_LIKE][nelem][ipH1s]*collider);
00723
00724
00725 ionbal.CollIonRate_Ground[nelem][nelem][1] =
00726 (iso.ColIoniz[ipH_LIKE][nelem][ipH1s]*collider*
00727 rfield.anu[Heavy.ipHeavy[nelem][nelem]-1]* EN1RYD);
00728
00729 DEBUG_EXIT( "HydroCollid()" );
00730 return;
00731 }
00732
00733
00734 #define HCSAR(ilo , ihi , nte ) (*(HCS+(nte)+(ilo)*48+(ihi)*8))
00735
00736
00737 static float HCSAR_interp( int ipLo , int ipHi )
00738 {
00739
00740 static int ip=1;
00741 float cs;
00742
00743 DEBUG_ENTRY( "HCSAR_interp()" );
00744
00745 if( ipLo==1 && ipHi==2 )
00746 {
00747 fprintf(ioQQQ,"HCSAR_interp was called for the 2s-2p transition, which it cannot do\n");
00748 cdEXIT(EXIT_FAILURE);
00749 }
00750 if( phycon.te <= HCSTE[0] )
00751 {
00752 cs = HCSAR( ipLo , ipHi , 0 );
00753 }
00754 else if( phycon.te >= HCSTE[NHCSTE-1] )
00755 {
00756 cs = HCSAR( ipLo , ipHi , NHCSTE-1 );
00757 }
00758 else
00759 {
00760
00761 if( !(HCSTE[ip-1] < phycon.te ) && ( phycon.te <= HCSTE[ip]) )
00762 {
00763
00764 for( ip=1; ip<NHCSTE; ++ip )
00765 {
00766 if( (HCSTE[ip-1] < phycon.te ) && ( phycon.te <= HCSTE[ip]) )
00767 break;
00768 }
00769 }
00770
00771 cs = HCSAR( ipLo , ipHi , ip-1 ) +
00772 (HCSAR( ipLo , ipHi , ip ) - HCSAR( ipLo , ipHi , ip-1 ) ) / (HCSTE[ip]-HCSTE[ip-1] ) *
00773 ((float)phycon.te - HCSTE[ip-1] );
00774 if( cs <= 0.)
00775 {
00776 fprintf(ioQQQ," insane cs returned by HCSAR_interp, values are\n");
00777 fprintf(ioQQQ,"%.3f %.3f \n", HCSAR( ipLo , ipHi , ip-1 ),HCSAR( ipLo , ipHi , ip ) );
00778 }
00779 }
00780
00781 DEBUG_EXIT( "HCSAR_interp()" );
00782 return(cs);
00783 }
00784
00785
00786
00787
00788 static double Hydcs123(
00789
00790
00791 long int ipLow,
00792
00793 long int ipHi,
00794
00795 long int nelem,
00796
00797 long int chType)
00798 {
00799 long int i,
00800 j,
00801 k;
00802 double C,
00803 D,
00804 EE,
00805 expq ,
00806 Hydcs123_v,
00807 Ratehigh,
00808 Ratelow,
00809 TeUse,
00810 gLo,
00811 gHi,
00812 q,
00813 rate,
00814 slope,
00815 temp,
00816 temphigh,
00817 templow,
00818 tev,
00819 x,
00820 QuanNLo,
00821 QuanNUp,
00822 Charge,
00823 ChargeSquared,
00824 zhigh,
00825 zlow;
00826 static const double ap[5] = {-2113.113,729.0084,1055.397,854.632,938.9912};
00827 static const double bp[5] = {-6783.515,-377.7190,724.1936,493.1107,735.7466};
00828 static const double cp[5] = {-3049.719,226.2320,637.8630,388.5465,554.6369};
00829 static const double dp[5] = {3514.5153,88.60169,-470.4055,-329.4914,-450.8459};
00830 static const double ep[5] = {0.005251557,0.009059154,0.008725781,0.009952418,0.01098687};
00831 static const double ae[5] = {-767.5859,-643.1189,-461.6836,-429.0543,-406.5285};
00832 static const double be[5] = {-1731.9178,-1442.548,-1055.364,-980.3079,-930.9266};
00833 static const double ce[5] = {-939.1834,-789.9569,-569.1451,-530.1974,-502.0939};
00834 static const double de[5] = {927.4773,773.2008,564.3272,524.2944,497.7763};
00835 static const double ee[5] = {-0.002528027,-0.003793665,-0.002122103,-0.002234207,-0.002317720};
00836 static const double A[2] = {4.4394,0.0};
00837 static const double B[2] = {0.8949,0.8879};
00838 static const double C0[2] = {-0.6012,-0.2474};
00839 static const double C1[2] = {-3.9710,-3.7562};
00840 static const double C2[2] = {-4.2176,2.0491};
00841 static const double D0[2] = {2.930,0.0539};
00842 static const double D1[2] = {1.7990,3.4009};
00843 static const double D2[2] = {4.9347,-1.7770};
00844
00845 DEBUG_ENTRY( "Hydcs123()" );
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872 ASSERT( nelem >= 0 );
00873 ASSERT( nelem < LIMELM );
00874
00875
00876 ASSERT( ipLow > 0);
00877 ASSERT( ipLow <= 3);
00878 ASSERT( ipHi > 1 );
00879 ASSERT( ipHi <=4 );
00880
00881
00882 if( ipHi == 4 )
00883 {
00884
00885 QuanNUp = 3.;
00886 gHi = 18.;
00887
00888
00889
00890 i = -1;
00891 }
00892 else if( ipHi == 3 )
00893 {
00894
00895 QuanNUp = 2.;
00896 gHi = 6.;
00897
00898 i = 0;
00899 }
00900 else if( ipHi == 2 )
00901 {
00902
00903 QuanNUp = 2.;
00904 gHi = 2.;
00905
00906 i = 1;
00907 }
00908 else
00909 {
00910 fprintf( ioQQQ, " Insane levels in Hydcs123\n" );
00911 puts( "[Stop in hydcs123]" );
00912 cdEXIT(EXIT_FAILURE);
00913 }
00914
00915
00916 if( ipLow == 1 )
00917 {
00918
00919 QuanNLo = 1.;
00920 gLo = 2.;
00921 }
00922 else if( ipLow == 2 )
00923 {
00924
00925 QuanNLo = 2.;
00926 gLo = 2.;
00927 }
00928 else if( ipLow == 3 )
00929 {
00930 QuanNLo = 2.;
00931 gLo = 6.;
00932 }
00933 else
00934 {
00935 fprintf( ioQQQ, " Insane levels in Hydcs123\n" );
00936 puts( "[Stop in hydcs123]" );
00937 cdEXIT(EXIT_FAILURE);
00938 }
00939
00940 if( nelem == 0 )
00941 {
00942
00943 Hydcs123_v = H1cs123(ipLow,ipHi,chType);
00944
00945 DEBUG_EXIT( "Hydcs123()" );
00946 return( Hydcs123_v );
00947 }
00948
00949
00950
00951
00952 Charge = (double)(nelem + 1);
00953
00954 ChargeSquared = Charge*Charge;
00955
00956 if( ipLow == 2 && ipHi == 3 )
00957 {
00958
00959
00960 if( nelem < 1 )
00961 {
00962
00963
00964 fprintf( ioQQQ, " insane charge given to Hydcs123\n" );
00965 puts( "[Stop in hydcs123]" );
00966 cdEXIT(EXIT_FAILURE);
00967 }
00968 else if( nelem == 1 )
00969 {
00970 zlow = 2.;
00971 j = 1;
00972 zhigh = 2.;
00973 k = 1;
00974 }
00975
00976 else if( nelem <= 5 )
00977 {
00978 zlow = 2.;
00979 j = 1;
00980 zhigh = 6.;
00981 k = 2;
00982 }
00983 else if( nelem <= 11 )
00984 {
00985 zlow = 6.;
00986 j = 2;
00987 zhigh = 12.;
00988 k = 3;
00989 }
00990 else if( nelem <= 15 )
00991 {
00992 zlow = 12.;
00993 j = 3;
00994 zhigh = 16.;
00995 k = 4;
00996 }
00997 else if( nelem <= 17 )
00998 {
00999 zlow = 16.;
01000 j = 4;
01001 zhigh = 18.;
01002 k = 5;
01003 }
01004
01005
01006
01007 else
01008 {
01009 zlow = 18.;
01010 j = 5;
01011 zhigh = 18.;
01012 k = 5;
01013 }
01014
01015
01016
01017 x = EVRYD/TE1RYD*phycon.te/(27.211396*pow2(zlow));
01018 TeUse = MIN2(x,0.80);
01019 x = MAX2(0.025,TeUse);
01020
01021
01022 if( chType == 'e' )
01023 {
01024
01025 Ratelow = ae[j-1] + be[j-1]*x + ce[j-1]*pow2(x)*log(x) + de[j-1]*
01026 exp(x) + ee[j-1]*log(x)/pow2(x);
01027 }
01028 else if( chType == 'p' )
01029 {
01030 Ratelow = ap[j-1] + bp[j-1]*x + cp[j-1]*pow2(x)*log(x) + dp[j-1]*
01031 exp(x) + ep[j-1]*log(x)/pow2(x);
01032 }
01033 else
01034 {
01035
01036 fprintf( ioQQQ, " insane collision species given to Hydcs123\n" );
01037 puts( "[Stop in hydcs123]" );
01038 cdEXIT(EXIT_FAILURE);
01039 }
01040
01041
01042 x = EVRYD/TE1RYD*phycon.te/(27.211396*pow2(zhigh));
01043 TeUse = MIN2(x,0.80);
01044 x = MAX2(0.025,TeUse);
01045 if( chType == 'e' )
01046 {
01047 Ratehigh = ae[k-1] + be[k-1]*x + ce[k-1]*pow2(x)*log(x) +
01048 de[k-1]*exp(x) + ee[k-1]*log(x)/pow2(x);
01049 }
01050 else
01051 {
01052 Ratehigh = ap[k-1] + bp[k-1]*x + cp[k-1]*pow2(x)*log(x) +
01053 dp[k-1]*exp(x) + ep[k-1]*log(x)/pow2(x);
01054 }
01055
01056
01057 if( zlow == zhigh )
01058
01059 {
01060 rate = Ratelow;
01061 }
01062 else
01063 {
01064 slope = (Ratehigh - Ratelow)/(zhigh - zlow);
01065 rate = slope*(Charge - zlow) + Ratelow;
01066 }
01067 rate = rate/ChargeSquared/Charge*1.0e-7;
01068
01069 templow = 0.025*27.211396*TE1RYD/EVRYD*ChargeSquared;
01070 temphigh = 0.80*27.211396*TE1RYD/EVRYD*ChargeSquared;
01071 TeUse = MIN2((double)phycon.te,temphigh);
01072 temp = MAX2(TeUse,templow);
01073 Hydcs123_v = rate*gHi*sqrt(temp)/COLL_CONST;
01074
01075 }
01076 else if( ipHi == 4 )
01077 {
01078
01079
01080 if( nelem < 1 )
01081 {
01082 fprintf( ioQQQ, " insane charge given to Hydcs123\n" );
01083 puts( "[Stop in hydcs123]" );
01084 cdEXIT(EXIT_FAILURE);
01085 }
01086 else if( nelem == 1 )
01087 {
01088 zlow = 2.;
01089 Ratelow = He2cs123(ipLow,ipHi);
01090 zhigh = 2.;
01091 Ratehigh = Ratelow;
01092 }
01093 else if( nelem > 1 && nelem <= 5 )
01094 {
01095 zlow = 2.;
01096 Ratelow = He2cs123(ipLow,ipHi);
01097 zhigh = 6.;
01098 Ratehigh = C6cs123(ipLow,ipHi);
01099 }
01100 else if( nelem > 5 && nelem <= 9 )
01101 {
01102 zlow = 6.;
01103 Ratelow = C6cs123(ipLow,ipHi);
01104 zhigh = 10.;
01105 Ratehigh = Ne10cs123(ipLow,ipHi);
01106 }
01107 else if( nelem > 9 && nelem <= 19 )
01108 {
01109 zlow = 10.;
01110 Ratelow = Ne10cs123(ipLow,ipHi);
01111 zhigh = 20.;
01112 Ratehigh = Ca20cs123(ipLow,ipHi);
01113 }
01114 else if( nelem > 19 && nelem <= 25 )
01115 {
01116 zlow = 20.;
01117 Ratelow = Ca20cs123(ipLow,ipHi);
01118 zhigh = 26.;
01119 Ratehigh = Fe26cs123(ipLow,ipHi);
01120 }
01121
01122
01123 else
01124 {
01125 Charge = 26.;
01126 zlow = 26.;
01127 Ratelow = Fe26cs123(ipLow,ipHi);
01128 zhigh = 26.;
01129 Ratehigh = Ratelow;
01130 }
01131
01132
01133
01134 if( zlow == zhigh )
01135
01136 {
01137 rate = Ratelow;
01138 }
01139 else
01140 {
01141 slope = (Ratehigh - Ratelow)/(zhigh - zlow);
01142 rate = slope*(Charge - zlow) + Ratelow;
01143 }
01144 Hydcs123_v = rate;
01145
01146 }
01147 else
01148 {
01149
01150 if( nelem == 1 )
01151 {
01152
01153 Hydcs123_v = He2cs123(ipLow,ipHi);
01154
01155 DEBUG_EXIT( "Hydcs123()" );
01156 return( Hydcs123_v );
01157 }
01158
01159
01160 tev = phycon.te / EVDEGK;
01161
01162 EE = ChargeSquared*EVRYD*(1./QuanNLo/QuanNLo - 1./QuanNUp/QuanNUp);
01163
01164 q = EE/tev;
01165 TeUse = MIN2(q,10.);
01166
01167 q = MAX2(1.,TeUse);
01168 expq = exp(q);
01169
01170
01171 ASSERT( i==0 || i==1 );
01172 C = C0[i] + C1[i]/Charge + C2[i]/ChargeSquared;
01173 D = D0[i] + D1[i]/Charge + D2[i]/ChargeSquared;
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194
01195
01196
01197
01198
01199 rate = (B[i] + D*q)/expq + (A[i] + C*q - D*q*q)*
01200 ee1(q);
01201 rate *= 8.010e-8/2./ChargeSquared/tev*sqrt(tev);
01202
01203 rate *= expq*gLo/gHi;
01204
01205
01206 Hydcs123_v = rate*gHi*phycon.sqrte/COLL_CONST;
01207 }
01208
01209 DEBUG_EXIT( "Hydcs123()" );
01210 return( Hydcs123_v );
01211 }
01212
01213
01214 static double C6cs123(long int i,
01215 long int j)
01216 {
01217 double C6cs123_v,
01218 TeUse,
01219 t,
01220 x;
01221 static const double a[3] = {-92.23774,-1631.3878,-6326.4947};
01222 static const double b[3] = {-11.93818,-218.3341,-849.8927};
01223 static const double c[3] = {0.07762914,1.50127,5.847452};
01224 static const double d[3] = {78.401154,1404.8475,5457.9291};
01225 static const double e[3] = {332.9531,5887.4263,22815.211};
01226
01227 DEBUG_ENTRY( "C6cs123()" );
01228
01229
01230
01231
01232
01233
01234
01235
01236
01237
01238
01239
01240
01241 TeUse = MAX2(phycon.te,6310.);
01242 t = MIN2(TeUse,1.6e6);
01243 x = log10(t);
01244
01245 if( i == 1 && j == 2 )
01246 {
01247
01248 fprintf( ioQQQ, " Carbon VI 2s-1s not done in C6cs123\n" );
01249 puts( "[Stop in c6cs123]" );
01250 cdEXIT(EXIT_FAILURE);
01251 }
01252
01253 else if( i == 1 && j == 3 )
01254 {
01255
01256 fprintf( ioQQQ, " Carbon VI 2p-1s not done in C6cs123\n" );
01257 puts( "[Stop in c6cs123]" );
01258 cdEXIT(EXIT_FAILURE);
01259 }
01260
01261 else if( i == 1 && j == 4 )
01262 {
01263
01264 C6cs123_v = a[0] + b[0]*x + c[0]*pow2(x)*sqrt(x) + d[0]*log(x) +
01265 e[0]*log(x)/pow2(x);
01266
01267
01268 }
01269 else if( i == 2 && j == 4 )
01270 {
01271
01272 C6cs123_v = a[1] + b[1]*x + c[1]*pow2(x)*sqrt(x) + d[1]*log(x) +
01273 e[1]*log(x)/pow2(x);
01274
01275
01276 }
01277 else if( i == 3 && j == 4 )
01278 {
01279
01280 C6cs123_v = a[2] + b[2]*x + c[2]*pow2(x)*sqrt(x) + d[2]*log(x) +
01281 e[2]*log(x)/pow2(x);
01282
01283
01284 }
01285 else
01286 {
01287 fprintf( ioQQQ, " insane levels for C VI n=1,2,3 !!!\n" );
01288 puts( "[Stop in c6cs123]" );
01289 cdEXIT(EXIT_FAILURE);
01290 }
01291
01292 DEBUG_EXIT( "C6cs123()" );
01293 return( C6cs123_v );
01294 }
01295
01296
01297 static double Ca20cs123(long int i,
01298 long int j)
01299 {
01300 double Ca20cs123_v,
01301 TeUse,
01302 t,
01303 x;
01304 static const double a[3] = {-12.5007,-187.2303,-880.18896};
01305 static const double b[3] = {-1.438749,-22.17467,-103.1259};
01306 static const double c[3] = {0.008219688,0.1318711,0.6043752};
01307 static const double d[3] = {10.116516,153.2650,717.4036};
01308 static const double e[3] = {45.905343,685.7049,3227.2836};
01309
01310 DEBUG_ENTRY( "Ca20cs123()" );
01311
01312
01313
01314
01315
01316
01317
01318
01319
01320
01321
01322
01323
01324
01325
01326 TeUse = MAX2(phycon.te,1.0e5);
01327 t = MIN2(TeUse,1.585e7);
01328 x = log10(t);
01329
01330 if( i == 1 && j == 2 )
01331 {
01332
01333 fprintf( ioQQQ, " Ca XX 2s-1s not done in Ca20cs123\n" );
01334 puts( "[Stop in ca20cs123]" );
01335 cdEXIT(EXIT_FAILURE);
01336 }
01337
01338 else if( i == 1 && j == 3 )
01339 {
01340
01341 fprintf( ioQQQ, " Ca XX 2p-1s not done in Ca20cs123\n" );
01342 puts( "[Stop in ca20cs123]" );
01343 cdEXIT(EXIT_FAILURE);
01344 }
01345
01346 else if( i == 1 && j == 4 )
01347 {
01348
01349 Ca20cs123_v = a[0] + b[0]*x + c[0]*pow2(x)*sqrt(x) + d[0]*log(x) +
01350 e[0]*log(x)/pow2(x);
01351
01352
01353 }
01354 else if( i == 2 && j == 4 )
01355 {
01356
01357 Ca20cs123_v = a[1] + b[1]*x + c[1]*pow2(x)*sqrt(x) + d[1]*log(x) +
01358 e[1]*log(x)/pow2(x);
01359
01360
01361 }
01362 else if( i == 3 && j == 4 )
01363 {
01364
01365 Ca20cs123_v = a[2] + b[2]*x + c[2]*pow2(x)*sqrt(x) + d[2]*log(x) +
01366 e[2]*log(x)/pow2(x);
01367
01368
01369 }
01370 else
01371 {
01372 fprintf( ioQQQ, " insane levels for Ca XX n=1,2,3 !!!\n" );
01373 puts( "[Stop in ca20cs123]" );
01374 cdEXIT(EXIT_FAILURE);
01375 }
01376
01377 DEBUG_EXIT( "Ca20cs123()" );
01378 return( Ca20cs123_v );
01379 }
01380
01381
01382 static double Ne10cs123(long int i,
01383 long int j)
01384 {
01385 double Ne10cs123_v,
01386 TeUse,
01387 t,
01388 x;
01389 static const double a[3] = {3.346644,151.2435,71.7095};
01390 static const double b[3] = {0.5176036,20.05133,13.1543};
01391 static const double c[3] = {-0.00408072,-0.1311591,-0.1099238};
01392 static const double d[3] = {-3.064742,-129.8303,-71.0617};
01393 static const double e[3] = {-11.87587,-541.8599,-241.2520};
01394
01395 DEBUG_ENTRY( "Ne10cs123()" );
01396
01397
01398
01399
01400
01401
01402
01403
01404
01405
01406
01407
01408
01409 TeUse = MAX2(phycon.te,6310.);
01410 t = MIN2(TeUse,1.6e6);
01411 x = log10(t);
01412
01413 if( i == 1 && j == 2 )
01414 {
01415
01416 fprintf( ioQQQ, " Neon X 2s-1s not done in Ne10cs123\n" );
01417 puts( "[Stop in ne10cs123]" );
01418 cdEXIT(EXIT_FAILURE);
01419 }
01420
01421 else if( i == 1 && j == 3 )
01422 {
01423
01424 fprintf( ioQQQ, " Neon X 2p-1s not done in Ne10cs123\n" );
01425 puts( "[Stop in ne10cs123]" );
01426 cdEXIT(EXIT_FAILURE);
01427 }
01428
01429 else if( i == 1 && j == 4 )
01430 {
01431
01432 Ne10cs123_v = a[0] + b[0]*x + c[0]*pow2(x)*sqrt(x) + d[0]*log(x) +
01433 e[0]*log(x)/pow2(x);
01434
01435
01436 }
01437 else if( i == 2 && j == 4 )
01438 {
01439
01440 Ne10cs123_v = a[1] + b[1]*x + c[1]*pow2(x)*sqrt(x) + d[1]*log(x) +
01441 e[1]*log(x)/pow2(x);
01442
01443
01444 }
01445 else if( i == 3 && j == 4 )
01446 {
01447
01448 Ne10cs123_v = a[2] + b[2]*x + c[2]*pow2(x)*sqrt(x) + d[2]*log(x) +
01449 e[2]*log(x)/pow2(x);
01450
01451
01452 }
01453 else
01454 {
01455 fprintf( ioQQQ, " insane levels for Ne X n=1,2,3 !!!\n" );
01456 puts( "[Stop in ne10cs123]" );
01457 cdEXIT(EXIT_FAILURE);
01458 }
01459
01460 DEBUG_EXIT( "Ne10cs123()" );
01461 return( Ne10cs123_v );
01462 }
01463
01464
01465 static double H1cs123(long int i,
01466 long int j,
01467 long int chType)
01468 {
01469 long int k;
01470 double H1cs123_v,
01471 T,
01472 TeUse,
01473 Tmax1,
01474 Tmax2,
01475 Tmax3,
01476 Tmin,
01477 bk,
01478 cs,
01479 cs3d,
01480 cs3p,
01481 cs3s,
01482 ev,
01483 rate,
01484 ryd,
01485 te03,
01486 te05,
01487 te10,
01488 te20;
01489 static const double b[10][11] = {
01490 {0.25501041,0.36105477,0.06461453,0.09895325,0.05742134,1.2168,1.9985,
01491 0.9382,2.0097,6.5316,8.9833},
01492 {0.09973588,1.73103797,-0.11160826,0.23240113,-0.00609531,4.7349,13.2632,
01493 32.8948,3.7515,26.8258,126.0513},
01494 {-0.05404749,0.37297254,0.76957269,0.20112555,0.93848281,6.5942,11.912,
01495 -17.9286,-7.9734,-18.2834,-21.9017},
01496 {0.02520311,-0.43240114,-2.18918875,0.87813331,-3.40393911,-14.9419,-13.0364,
01497 -0.8772,10.9359,2.22,20.7445},
01498 {-0.00405306,0.14296095,3.14216481,-3.34239518,5.50742903,10.7255,5.672,5.8667,
01499 -7.1074,4.1958,-20.3788},
01500 {0.03508995,-0.01670572,-2.23498961,3.7075526,-4.25196238,-3.3407,-1.2155,
01501 -2.4439,2.1359,-2.0916,7.3643},
01502 {28.056,0.,0.62651998,-1.37320827,1.2687104,0.3823,0.1034,0.3147,-0.2405,0.289,
01503 -0.9082},
01504 {7.2945,0.,0.0112515,0.,0.00724136,0.,0.,0.,0.,0.,0.},
01505 {0.,0.,28.056,0.,28.056,0.,0.,0.,0.,0.,0.},
01506 {0.,0.,7.2945,0.,7.2945,0.,0.,0.,0.,0.,0.}
01507 };
01508
01509 DEBUG_ENTRY( "H1cs123()" );
01510
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529 ev = 1.60184e-12;
01530 bk = 1.380622e-16;
01531 ryd = EVRYD ;
01532
01533
01534
01535 T = bk*phycon.te/(ev*ryd);
01536
01537
01538
01539
01540
01541
01542
01543
01544
01545
01546
01547
01548
01549
01550 Tmin = 0.01;
01551
01552
01553
01554 Tmax1 = 4.0;
01555
01556
01557
01558 Tmax2 = 1.2667;
01559
01560
01561
01562 Tmax3 = 1.0;
01563
01564
01565
01566 cs = 0.;
01567
01568
01569
01570 if( i == 1 && j == 2 )
01571 {
01572
01573
01574
01575 TeUse = MIN2(T,Tmax1);
01576 T = MAX2(TeUse,Tmin);
01577
01578 cs = 0.;
01579 for( k=0; k < 5; k++ )
01580 {
01581 cs += b[k][0]*powi(T,k);
01582 }
01583
01584 H1cs123_v = cs + b[5][0]*log(b[6][0]*T)*exp(-b[7][0]*T);
01585 }
01586
01587 else if( i == 1 && j == 3 )
01588 {
01589
01590
01591
01592 TeUse = MIN2(T,Tmax1);
01593 T = MAX2(TeUse,Tmin);
01594
01595 cs = 0.;
01596 for( k=0; k < 6; k++ )
01597 {
01598 cs += b[k][1]*powi(T,k);
01599 }
01600 H1cs123_v = cs;
01601 }
01602
01603 else if( i == 1 && j == 4 )
01604 {
01605
01606
01607
01608 TeUse = MIN2(T,Tmax2);
01609 T = MAX2(TeUse,Tmin);
01610
01611 cs3s = 0.;
01612 for( k=0; k < 7; k++ )
01613 {
01614 cs3s += b[k][2]*powi(T,k);
01615 }
01616
01617 cs3s += b[7][2]*log(b[8][2]*T)*exp(-b[9][2]*T);
01618
01619
01620 cs3p = 0.;
01621 for( k=0; k < 7; k++ )
01622 {
01623 cs3p += b[k][3]*powi(T,k);
01624 }
01625
01626
01627 cs3d = 0.;
01628 for( k=0; k < 7; k++ )
01629 {
01630 cs3d += b[k][4]*powi(T,k);
01631 }
01632
01633 cs3d += b[7][4]*log(b[8][4]*T)*exp(-b[9][4]*T);
01634
01635
01636 H1cs123_v = cs3s + cs3p + cs3d;
01637 }
01638
01639 else if( i == 2 && j == 4 )
01640 {
01641
01642
01643
01644 TeUse = MIN2(T,Tmax3);
01645 T = MAX2(TeUse,Tmin);
01646
01647 cs3s = 0.;
01648 for( k=0; k < 7; k++ )
01649 {
01650 cs3s += b[k][5]*powi(T,k);
01651 }
01652
01653
01654 cs3p = 0.;
01655 for( k=0; k < 7; k++ )
01656 {
01657 cs3p += b[k][6]*powi(T,k);
01658 }
01659
01660
01661 cs3d = 0.;
01662 for( k=0; k < 7; k++ )
01663 {
01664 cs3d += b[k][7]*powi(T,k);
01665 }
01666
01667
01668 H1cs123_v = cs3s + cs3p + cs3d;
01669 }
01670
01671 else if( i == 3 && j == 4 )
01672 {
01673
01674
01675
01676 TeUse = MIN2(T,Tmax3);
01677 T = MAX2(TeUse,Tmin);
01678
01679 cs3s = 0.;
01680 for( k=0; k < 7; k++ )
01681 {
01682 cs3s += b[k][8]*powi(T,k);
01683 }
01684
01685
01686 cs3p = 0.;
01687 for( k=0; k < 7; k++ )
01688 {
01689 cs3p += b[k][9]*powi(T,k);
01690 }
01691
01692
01693 cs3d = 0.;
01694 for( k=0; k < 7; k++ )
01695 {
01696 cs3d += b[k][10]*powi(T,k);
01697 }
01698
01699
01700 H1cs123_v = cs3s + cs3p + cs3d;
01701 }
01702
01703 else if( i == 2 && j == 3 )
01704 {
01705
01706
01707 te20 = pow(phycon.te,.2);
01708 te03 = pow(phycon.te,.03);
01709 te10 = pow(phycon.te,.10);
01710 te05 = pow(phycon.te,.05);
01711 if( chType == 'e' )
01712 {
01713 rate = 5.738e-4/(te20*te20/te03);
01714 }
01715 else
01716 {
01717 rate = 6.290e-4/(te10*te05);
01718 }
01719
01720 H1cs123_v = rate*sqrt(phycon.te)*6./COLL_CONST;
01721
01722 }
01723 else
01724 {
01725 fprintf( ioQQQ, " insane levels for H I n=1,2,3 !!!\n" );
01726 fprintf( ioQQQ, "%4ld%4ld\n", i, j );
01727 puts( "[Stop in h1cs123]" );
01728 H1cs123_v = 0.;
01729 cdEXIT(EXIT_FAILURE);
01730 }
01731
01732 DEBUG_EXIT( "H1cs123()" );
01733 return( H1cs123_v );
01734 }
01735
01736
01737 static double He2cs123(long int i,
01738 long int j)
01739 {
01740 double He2cs123_v,
01741 cs3d,
01742 cs3p,
01743 cs3s,
01744 t;
01745 static const double a[11]={0.12176209,0.32916723,0.46546497,0.044501688,
01746 0.040523277,0.5234889,1.4903214,1.4215094,1.0295881,4.769306,9.7226127};
01747 static const double b[11]={0.039936166,2.9711166e-05,-0.020835863,3.0508137e-04,
01748 -2.004485e-15,4.41475e-06,1.0622666e-05,2.0538877e-06,0.80638448,2.0967075e-06,
01749 7.6089851e-05};
01750 static const double c[11]={143284.77,0.73158545,-2.159172,0.43254802,2.1338557,
01751 8.9899702e-06,-2.9001451e-12,1.762076e-05,52741.735,-2153.1219,-3.3996921e-11};
01752
01753 DEBUG_ENTRY( "He2cs123()" );
01754
01755
01756
01757
01758
01759
01760
01761
01762
01763
01764
01765
01766
01767
01768 t = phycon.te;
01769 if( t < 5000. )
01770 {
01771 t = 5000.;
01772 }
01773 else if( t > 5.0e05 )
01774 {
01775 t = 5.0e05;
01776 }
01777
01778
01779
01780 if( i == 1 && j == 2 )
01781 {
01782
01783 He2cs123_v = a[0] + b[0]*exp(-t/c[0]);
01784
01785
01786 }
01787 else if( i == 1 && j == 3 )
01788 {
01789
01790 He2cs123_v = a[1] + b[1]*pow(t,c[1]);
01791
01792
01793 }
01794 else if( i == 1 && j == 4 )
01795 {
01796
01797 cs3s = a[2] + b[2]*log(t) + c[2]/log(t);
01798
01799
01800 cs3p = a[3] + b[3]*pow(t,c[3]);
01801
01802
01803 cs3d = a[4] + b[4]*pow(t,c[4]);
01804
01805
01806 He2cs123_v = cs3s + cs3p + cs3d;
01807
01808
01809 }
01810 else if( i == 2 && j == 4 )
01811 {
01812
01813 cs3s = (a[5] + c[5]*t)/(1 + b[5]*t);
01814
01815
01816 cs3p = a[6] + b[6]*t + c[6]*t*t;
01817
01818
01819 cs3d = (a[7] + c[7]*t)/(1 + b[7]*t);
01820
01821
01822 He2cs123_v = cs3s + cs3p + cs3d;
01823
01824
01825 }
01826 else if( i == 3 && j == 4 )
01827 {
01828
01829 cs3s = a[8] + b[8]*exp(-t/c[8]);
01830
01831
01832 cs3p = a[9] + b[9]*t + c[9]/t;
01833
01834
01835 cs3d = a[10] + b[10]*t + c[10]*t*t;
01836
01837
01838 He2cs123_v = cs3s + cs3p + cs3d;
01839
01840
01841 }
01842 else
01843 {
01844
01845 fprintf( ioQQQ, " insane levels for He II n=1,2,3 !!!\n" );
01846 puts( "[Stop in he2cs123]" );
01847 cdEXIT(EXIT_FAILURE);
01848 }
01849
01850 DEBUG_EXIT( "He2cs123()" );
01851 return( He2cs123_v );
01852 }
01853
01854
01855 static double Fe26cs123(long int i,
01856 long int j)
01857 {
01858 double Fe26cs123_v,
01859 TeUse,
01860 t,
01861 x;
01862 static const double a[3] = {-4.238398,-238.2599,-1211.5237};
01863 static const double b[3] = {-0.4448177,-27.06869,-136.7659};
01864 static const double c[3] = {0.0022861,0.153273,0.7677703};
01865 static const double d[3] = {3.303775,191.7165,972.3731};
01866 static const double e[3] = {15.82689,878.1333,4468.696};
01867
01868 DEBUG_ENTRY( "Fe26cs123()" );
01869
01870
01871
01872
01873
01874
01875
01876
01877
01878
01879
01880
01881
01882 TeUse = MAX2(phycon.te,1.585e5);
01883 t = MIN2(TeUse,1.585e7);
01884 x = log10(t);
01885
01886 if( i == 1 && j == 2 )
01887 {
01888
01889 fprintf( ioQQQ, " Fe XXVI 2s-1s not done in Fe26cs123\n" );
01890 puts( "[Stop in fe26cs123]" );
01891 cdEXIT(EXIT_FAILURE);
01892 }
01893
01894 else if( i == 1 && j == 3 )
01895 {
01896
01897 fprintf( ioQQQ, " Fe XXVI 2p-1s not done in Fe26cs123\n" );
01898 puts( "[Stop in fe26cs123]" );
01899 cdEXIT(EXIT_FAILURE);
01900 }
01901
01902 else if( i == 1 && j == 4 )
01903 {
01904
01905 Fe26cs123_v = a[0] + b[0]*x + c[0]*pow2(x)*sqrt(x) + d[0]*log(x) +
01906 e[0]*log(x)/pow2(x);
01907
01908
01909 }
01910 else if( i == 2 && j == 4 )
01911 {
01912
01913 Fe26cs123_v = a[1] + b[1]*x + c[1]*pow2(x)*sqrt(x) + d[1]*log(x) +
01914 e[1]*log(x)/pow2(x);
01915
01916
01917 }
01918 else if( i == 3 && j == 4 )
01919 {
01920
01921 Fe26cs123_v = a[2] + b[2]*x + c[2]*pow2(x)*sqrt(x) + d[2]*log(x) +
01922 e[2]*log(x)/pow2(x);
01923
01924
01925 }
01926 else
01927 {
01928 fprintf( ioQQQ, " insane levels for Ca XX n=1,2,3 !!!\n" );
01929 puts( "[Stop in fe26cs123]" );
01930 cdEXIT(EXIT_FAILURE);
01931 }
01932
01933 DEBUG_EXIT( "Fe26cs123()" );
01934 return( Fe26cs123_v );
01935 }
01936
01937