00001
00002
00003
00004
00005 #include "cddefines.h"
00006 #include "thirdparty.h"
00007 #include "trace.h"
00008 #include "secondaries.h"
00009 #include "atmdat.h"
00010 #include "phycon.h"
00011 #include "ionbal.h"
00012 #include "dense.h"
00013 #include "conv.h"
00014 #include "rfield.h"
00015 #include "he.h"
00016 #include "iso.h"
00017 #include "helike.h"
00018 #include "mole.h"
00019 #include "taulines.h"
00020 #include "physconst.h"
00021 #include "dynamics.h"
00022 #include "continuum.h"
00023
00024 #if defined (__ICC) && defined(__ia64) && __INTEL_COMPILER < 910
00025 #pragma optimize("", off)
00026 #endif
00027
00028
00029
00030
00031 void HeLikeLevel( long nelem)
00032 {
00033 int32 nerror;
00034 long int n, i, ipHi, ipLo, j,
00035 level;
00036 double HighestPColOut[2] = {0.,0.},
00037 sum_popn_ov_ion,
00038 TooSmall;
00039 bool lgNegPop;
00040 static bool lgAMAT_alloc=false;
00041 static double *amat;
00042 static double SaveHe23S_photorate=0.;
00043 static int32 *ipiv;
00044
00045 static double
00046 *creation ,
00047 *error,
00048 *work,
00049 **SaveZ,
00050 **z;
00051 static double
00052 *CollisionsGoingUp, *CollisionsGoingDown,
00053 *OutOfNGoingUp, *OutOfNGoingDown;
00054 static bool lgSpaceMalloc=false;
00055 static FILE *ioMATRIX;
00056
00057
00058 # define PARALLEL false
00059
00060 DEBUG_ENTRY( "HeLikeLevel()" );
00061
00062 ASSERT( nelem < LIMELM );
00063
00064
00065 if( PARALLEL )
00066 {
00067 ipiv = (int32 *)MALLOC(sizeof(int32)*(unsigned)iso.numLevels_local[ipHE_LIKE][nelem] );
00068
00069 creation = (double *)MALLOC(sizeof(double)*(unsigned)iso.numLevels_local[ipHE_LIKE][nelem] );
00070
00071 CollisionsGoingUp = (double *)MALLOC(sizeof(double)*(unsigned)max_num_levels );
00072
00073 CollisionsGoingDown = (double *)MALLOC(sizeof(double)*(unsigned)max_num_levels );
00074
00075 OutOfNGoingUp = (double *)MALLOC(sizeof(double)*(unsigned)(max_n + 1) );
00076
00077 OutOfNGoingDown = (double *)MALLOC(sizeof(double)*(unsigned)(max_n + 1) );
00078
00079 error = (double *)MALLOC(sizeof(double)*(unsigned)iso.numLevels_local[ipHE_LIKE][nelem] );
00080
00081 work = (double *)MALLOC(sizeof(double)*(unsigned)iso.numLevels_local[ipHE_LIKE][nelem] );
00082
00083
00084 SaveZ = (double **)MALLOC(sizeof(double *)*(unsigned)iso.numLevels_local[ipHE_LIKE][nelem] );
00085
00086 z = (double **)MALLOC(sizeof(double *)*(unsigned)iso.numLevels_local[ipHE_LIKE][nelem] );
00087
00088
00089 for( i=0; i<(iso.numLevels_local[ipHE_LIKE][nelem]); ++i )
00090 {
00091 SaveZ[i] = (double *)MALLOC(sizeof(double)*(unsigned)iso.numLevels_local[ipHE_LIKE][nelem] );
00092
00093 z[i] = (double *)MALLOC(sizeof(double)*(unsigned)iso.numLevels_local[ipHE_LIKE][nelem] );
00094 }
00095 }
00096 else if( !lgSpaceMalloc )
00097 {
00098
00099 lgSpaceMalloc = true;
00100
00101 ipiv = (int32 *)MALLOC(sizeof(int32)*(unsigned)(max_num_levels) );
00102
00103 creation = (double *)MALLOC(sizeof(double)*(unsigned)(max_num_levels) );
00104
00105 CollisionsGoingUp = (double *)MALLOC(sizeof(double)*(unsigned) (max_num_levels) );
00106
00107 CollisionsGoingDown = (double *)MALLOC(sizeof(double)*(unsigned) (max_num_levels) );
00108
00109 OutOfNGoingUp = (double *)MALLOC(sizeof(double)*(unsigned)(max_n + 1) );
00110
00111 OutOfNGoingDown = (double *)MALLOC(sizeof(double)*(unsigned)(max_n + 1) );
00112
00113 error = (double *)MALLOC(sizeof(double)*(unsigned)(max_num_levels) );
00114
00115 work = (double *)MALLOC(sizeof(double)*(unsigned)(max_num_levels) );
00116
00117
00118 SaveZ = (double **)MALLOC(sizeof(double *)*(unsigned)max_num_levels );
00119
00120 z = (double **)MALLOC(sizeof(double *)*(unsigned)max_num_levels );
00121
00122
00123 for( i=0; i<max_num_levels; ++i )
00124 {
00125 SaveZ[i] = (double *)MALLOC(sizeof(double)*(unsigned)(max_num_levels+2) );
00126
00127 z[i] = (double *)MALLOC(sizeof(double)*(unsigned)(max_num_levels+2) );
00128 }
00129 }
00130
00131 for( i=0; i < max_num_levels; ++i )
00132 {
00133 CollisionsGoingUp[i] = 0.;
00134 CollisionsGoingDown[i] = 0.;
00135 }
00136 sum_popn_ov_ion = 0.;
00137
00138 for( i=0; i <= max_n; ++i )
00139 {
00140 OutOfNGoingUp[i] = 0.;
00141 OutOfNGoingDown[i] = 0.;
00142 }
00143
00144
00145
00146
00147 if( nelem == ipHELIUM )
00148 {
00149
00150
00151
00152 if( nzone>1 )
00153 {
00154 double frac=0.5;
00155 iso.gamnc[ipHE_LIKE][nelem][ipHe2s3S] = frac*iso.gamnc[ipHE_LIKE][nelem][ipHe2s3S] +
00156 (1.-frac)*SaveHe23S_photorate;
00157 }
00158 SaveHe23S_photorate = iso.gamnc[ipHE_LIKE][nelem][ipHe2s3S];
00159 }
00160
00161 for( n=0; n < iso.numLevels_local[ipHE_LIKE][nelem]; n++ )
00162 {
00163
00164 creation[n] = iso.RateCont2Level[ipHE_LIKE][nelem][n];
00165 }
00166
00167
00168
00169 ASSERT( ionbal.CollIonRate_Ground[nelem][nelem-1][0]< SMALLFLOAT ||
00170 fabs( iso.ColIoniz[ipHE_LIKE][nelem][0]* dense.EdenHCorr /
00171 SDIV(ionbal.CollIonRate_Ground[nelem][nelem-1][0] ) - 1.) < 0.001 );
00172
00173 {
00174
00175 enum {DEBUG_LOC=false};
00176
00177 if( DEBUG_LOC && nelem==1 )
00178 {
00179 fprintf(ioQQQ,"hebugggg\t%li\t%.2e\t%.2e\t%.2e\n",
00180 nzone,
00181 iso.xIonSimple[ipHE_LIKE][nelem],
00182 ionbal.RateIonizTot[nelem][nelem-1] ,
00183 ionbal.RateRecomTot[nelem][nelem-1] );
00184 }
00185 }
00186
00187
00188
00189 if( nelem==ipHELIUM )
00190 TooSmall = 1e-20;
00191 else
00192 TooSmall = 1e-10;
00193
00194
00195
00196 if( (strcmp( iso.chTypeAtomSet[ipHE_LIKE] , "LOWT" )==0) ||
00197 iso.xIonSimple[ipHE_LIKE][nelem] < TooSmall )
00198 {
00199 strcpy( iso.chTypeAtomUsed[ipHE_LIKE][nelem], "Low T" );
00200
00201 SaveZ[0][0] = -1.;
00202
00203 ionbal.RateIonizTot[nelem][nelem-1] = iso.RateLevel2Cont[ipHE_LIKE][nelem][0];
00204
00205 lgNegPop = false;
00206
00207
00208
00209
00210 iso.pop_ion_ov_neut[ipHE_LIKE][nelem] = iso.xIonSimple[ipHE_LIKE][nelem];
00211 iso.Pop2Ion[ipHE_LIKE][nelem][0] = 1./SDIV(iso.pop_ion_ov_neut[ipHE_LIKE][nelem]);
00212 for( n=1; n < iso.numLevels_local[ipHE_LIKE][nelem]; n++ )
00213 {
00214 iso.Pop2Ion[ipHE_LIKE][nelem][n] = 0.;
00215 }
00216 helike.qTot2TripS[nelem] = 0.;
00217 }
00218 else
00219 {
00220 strcpy( iso.chTypeAtomUsed[ipHE_LIKE][nelem], "Popul" );
00221
00222 for( level=ipHe1s1S; level < iso.numLevels_local[ipHE_LIKE][nelem]; level++ )
00223 {
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233 z[level][level] = iso.RateLevel2Cont[ipHE_LIKE][nelem][level];
00234
00235
00236
00237
00238
00239
00240 if( level == ipHe2s3S )
00241
00242 helike.qTot2TripS[nelem] = iso.ColIoniz[ipHE_LIKE][nelem][level];
00243
00244
00245 for( ipLo=ipHe1s1S; ipLo < level; ipLo++ )
00246 {
00247 double coll_down = EmisLines[ipHE_LIKE][nelem][level][ipLo].ColUL * dense.EdenHCorr;
00248 double pump = (double)EmisLines[ipHE_LIKE][nelem][level][ipLo].pump *
00249 KillIfBelowPlasma(EmisLines[ipHE_LIKE][nelem][level][ipLo].EnergyWN/WAVNRYD);
00250 double RadDecay = EmisLines[ipHE_LIKE][nelem][level][ipLo].Aul*
00251 (EmisLines[ipHE_LIKE][nelem][level][ipLo].Pesc +
00252 EmisLines[ipHE_LIKE][nelem][level][ipLo].Pelec_esc +
00253 EmisLines[ipHE_LIKE][nelem][level][ipLo].Pdest)*
00254 KillIfBelowPlasma(EmisLines[ipHE_LIKE][nelem][level][ipLo].EnergyWN/WAVNRYD);
00255
00256 if( helike.lgRandErrGen )
00257 {
00258 coll_down *= helike.ErrorFactor[nelem][level][ipLo][IPCOLLIS];
00259 RadDecay *= helike.ErrorFactor[nelem][level][ipLo][IPRAD];
00260 pump *= helike.ErrorFactor[nelem][level][ipLo][IPRAD];
00261 }
00262
00263 z[ipLo][level] =
00264 -coll_down*
00265 (double)EmisLines[ipHE_LIKE][nelem][level][ipLo].gHi/
00266 (double)EmisLines[ipHE_LIKE][nelem][level][ipLo].gLo*
00267 iso.Boltzmann[ipHE_LIKE][nelem][level][ipLo]
00268 -pump;
00269
00270
00271
00272 z[level][level] += pump *
00273 (double)EmisLines[ipHE_LIKE][nelem][level][ipLo].gLo/
00274 (double)EmisLines[ipHE_LIKE][nelem][level][ipLo].gHi;
00275
00276
00277 z[level][level] += coll_down;
00278
00279 if( iso.quant_desig[ipHE_LIKE][nelem][level].n >
00280 iso.quant_desig[ipHE_LIKE][nelem][ipLo].n )
00281 CollisionsGoingDown[level] += coll_down;
00282
00283
00284 z[level][level] += RadDecay;
00285
00286 if( (nelem == ipHELIUM) && iso.quant_desig[ipHE_LIKE][nelem][level].n == iso.n_HighestResolved_local[ipHE_LIKE][nelem]
00287 && iso.quant_desig[ipHE_LIKE][nelem][level].l == 1 )
00288 {
00289 HighestPColOut[iso.quant_desig[ipHE_LIKE][nelem][level].s] += coll_down;
00290 }
00291
00292 if( (level == ipHe2s3S) && (iso.quant_desig[ipHE_LIKE][nelem][ipLo].s==0))
00293 {
00294 helike.qTot2TripS[nelem] += coll_down/dense.EdenHCorr;
00295 }
00296
00297 }
00298
00299
00300
00301 for( ipHi=level + 1; ipHi < iso.numLevels_local[ipHE_LIKE][nelem]; ipHi++ )
00302 {
00303 double coll_down = EmisLines[ipHE_LIKE][nelem][ipHi][level].ColUL * dense.EdenHCorr;
00304 double coll_up;
00305 double pump = (double)EmisLines[ipHE_LIKE][nelem][ipHi][level].pump *
00306 KillIfBelowPlasma(EmisLines[ipHE_LIKE][nelem][ipHi][level].EnergyWN/WAVNRYD);
00307 double RadDecay = EmisLines[ipHE_LIKE][nelem][ipHi][level].Aul*
00308 (EmisLines[ipHE_LIKE][nelem][ipHi][level].Pesc +
00309 EmisLines[ipHE_LIKE][nelem][ipHi][level].Pelec_esc +
00310 EmisLines[ipHE_LIKE][nelem][ipHi][level].Pdest)*
00311 KillIfBelowPlasma(EmisLines[ipHE_LIKE][nelem][ipHi][level].EnergyWN/WAVNRYD);
00312
00313 if( helike.lgRandErrGen )
00314 {
00315 coll_down *= helike.ErrorFactor[nelem][ipHi][level][IPCOLLIS];
00316 RadDecay *= helike.ErrorFactor[nelem][ipHi][level][IPRAD];
00317 pump *= helike.ErrorFactor[nelem][ipHi][level][IPRAD];
00318 }
00319
00320 coll_up = coll_down * iso.Boltzmann[ipHE_LIKE][nelem][ipHi][level] *
00321 (double)EmisLines[ipHE_LIKE][nelem][ipHi][level].gHi /
00322 (double)EmisLines[ipHE_LIKE][nelem][ipHi][level].gLo;
00323
00324 z[ipHi][level] =
00325 - RadDecay
00326 - pump*(double)EmisLines[ipHE_LIKE][nelem][ipHi][level].gLo/
00327 (double)EmisLines[ipHE_LIKE][nelem][ipHi][level].gHi
00328 - coll_down;
00329
00330
00331 z[level][level] += pump;
00332
00333
00334 z[level][level] += coll_up;
00335
00336 if( iso.quant_desig[ipHE_LIKE][nelem][level].n <
00337 iso.quant_desig[ipHE_LIKE][nelem][ipHi].n )
00338 {
00339 CollisionsGoingUp[level] += coll_up;
00340 }
00341
00342 if( (nelem == ipHELIUM) && iso.quant_desig[ipHE_LIKE][nelem][level].n == iso.n_HighestResolved_local[ipHE_LIKE][nelem]
00343 && iso.quant_desig[ipHE_LIKE][nelem][level].l == 1 )
00344 {
00345 HighestPColOut[iso.quant_desig[ipHE_LIKE][nelem][level].s] += coll_up;
00346 }
00347
00348 if( (level == ipHe2s3S) && (iso.quant_desig[ipHE_LIKE][nelem][ipHi].s==0))
00349 {
00350 helike.qTot2TripS[nelem] += coll_up/dense.EdenHCorr;
00351 }
00352
00353 }
00354 }
00355
00356
00357
00358
00359
00360 z[ipHe2s1S][ipHe1s1S] -= iso.TwoNu_induc_dn[ipHE_LIKE][nelem]*iso.lgInd2nu_On;
00361 z[ipHe1s1S][ipHe2s1S] -= iso.TwoNu_induc_up[ipHE_LIKE][nelem]*iso.lgInd2nu_On;
00362
00363
00364 z[ipHe1s1S][ipHe1s1S] += iso.TwoNu_induc_up[ipHE_LIKE][nelem]*iso.lgInd2nu_On;
00365 z[ipHe2s1S][ipHe2s1S] += iso.TwoNu_induc_dn[ipHE_LIKE][nelem]*iso.lgInd2nu_On;
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382 if( dynamics.lgAdvection )
00383 {
00384
00385
00386 creation[ipHe1s1S] += dynamics.Source[nelem][nelem-1]/SDIV(dense.xIonDense[nelem][nelem])*
00387 dynamics.lgISO[ipHE_LIKE];
00388
00389
00390 for( ipLo=ipHe1s1S; ipLo<iso.numLevels_local[ipHE_LIKE][nelem]; ++ipLo )
00391 {
00392 z[ipLo][ipLo] += dynamics.Rate*dynamics.lgISO[ipHE_LIKE];
00393 }
00394
00395
00396 }
00397
00398
00399
00405 # if 0
00406 if( Secondaries.Hx12[0][1] > 0. )
00407 {
00408
00409
00410 for( level=1; level < iso.numLevels_local[ipHE_LIKE][nelem]; level++ )
00411 {
00412 double RateUp , RateDown;
00413
00414 RateUp = Secondaries.Hx12[MIN2(nelem,1)][level];
00415 RateDown = RateUp * (double)EmisLines[ipHE_LIKE][nelem][level][ipH1s].gLo /
00416 (double)EmisLines[ipHE_LIKE][nelem][level][ipH1s].gHi;
00417
00418
00419
00421 z[ipH1s][ipH1s] += RateUp;
00422
00423
00424 z[level][ipH1s] -= RateDown;
00425
00426
00427 z[ipH1s][level] -= RateUp;
00428
00429 z[level][level] += RateDown;
00430 }
00431 }
00432 # endif
00433
00434 if( trace.lgTrace &&
00435 trace.lgIsoTraceFull[ipHE_LIKE] && (nelem == trace.ipIsoTrace[ipHE_LIKE]) )
00436 {
00437 fprintf( ioQQQ, " pop level others => (HeLikeLevel)\n" );
00438
00439 for( n=0; n < iso.numLevels_local[ipHE_LIKE][nelem]; n++ )
00440 {
00441 fprintf( ioQQQ, " HeI%2ld", n );
00442
00443
00444 for( j=0; j < iso.numLevels_local[ipHE_LIKE][nelem]; j++ )
00445 {
00446 fprintf( ioQQQ,"\t%.9e", z[j][n] );
00447 }
00448 fprintf( ioQQQ, "\n" );
00449 }
00450 fprintf(ioQQQ," recomb ");
00451 for( n=0; n < iso.numLevels_local[ipHE_LIKE][nelem]; n++ )
00452 {
00453 fprintf( ioQQQ,"\t%.9e", creation[n] );
00454 }
00455 fprintf( ioQQQ, "\n" );
00456 fprintf(ioQQQ," recomb ct %.2e co %.2e hectr %.2e hctr %.2e\n",
00457 atmdat.HeCharExcRecTotal,
00458 findspecies("CO")->hevmol ,
00459 atmdat.HeCharExcRecTo[nelem][nelem-1]*dense.xIonDense[ipHELIUM][0] ,
00460 atmdat.HCharExcRecTo[nelem][nelem-1]*dense.xIonDense[ipHYDROGEN][0] );
00461 }
00462
00463
00464 for( ipLo=ipHe1s1S; ipLo < iso.numLevels_local[ipHE_LIKE][nelem]; ipLo++ )
00465 {
00466 for( n=ipHe1s1S; n < iso.numLevels_local[ipHE_LIKE][nelem]; n++ )
00467 {
00468 SaveZ[n][ipLo] = z[n][ipLo];
00469 }
00470 }
00471
00472
00473 # ifdef AMAT
00474 # undef AMAT
00475 # endif
00476
00477 # define AMAT(I_,J_) (*(amat+(I_)*(iso.numLevels_local[ipHE_LIKE][nelem])+(J_)))
00478
00479
00480 if( PARALLEL )
00481 {
00482 amat=(double*)MALLOC( (sizeof(double)*(unsigned)(iso.numLevels_local[ipHE_LIKE][nelem]*iso.numLevels_local[ipHE_LIKE][nelem]) ));
00483 }
00484 else if( !lgAMAT_alloc )
00485 {
00486
00487 lgAMAT_alloc = true;
00488 amat=(double*)MALLOC( (sizeof(double)*(unsigned)((max_num_levels)*(max_num_levels)) ));
00489 }
00490
00491
00492 for( j=0; j < iso.numLevels_local[ipHE_LIKE][nelem]; j++ )
00493 {
00494 for( n=0; n < iso.numLevels_local[ipHE_LIKE][nelem]; n++ )
00495 {
00496 AMAT(n,j) = z[n][j];
00497 }
00498 }
00499
00500 nerror = 0;
00501
00502 getrf_wrapper(iso.numLevels_local[ipHE_LIKE][nelem],iso.numLevels_local[ipHE_LIKE][nelem],
00503 amat,iso.numLevels_local[ipHE_LIKE][nelem],ipiv,&nerror);
00504
00505
00506
00507
00508 getrs_wrapper('N',iso.numLevels_local[ipHE_LIKE][nelem],1,amat,iso.numLevels_local[ipHE_LIKE][nelem],ipiv,
00509 creation,iso.numLevels_local[ipHE_LIKE][nelem],&nerror);
00510
00511 if( nerror != 0 )
00512 {
00513 fprintf( ioQQQ, " HeLikeLevel: dgetrs finds singular or ill-conditioned matrix\n" );
00514 puts( "[Stop in helike]" );
00515 cdEXIT(EXIT_FAILURE);
00516 }
00517
00518 if( PARALLEL )
00519 free( amat);
00520
00521 lgNegPop = false;
00522 for( level=0; level < iso.numLevels_local[ipHE_LIKE][nelem]; level++ )
00523 {
00524 iso.Pop2Ion[ipHE_LIKE][nelem][level] = creation[level];
00525
00526
00527 if( iso.Pop2Ion[ipHE_LIKE][nelem][level] < 0. )
00528 lgNegPop = true;
00529 }
00530
00531
00532 for( level=iso.numLevels_local[ipHE_LIKE][nelem]; level < iso.numLevels_max[ipHE_LIKE][nelem]; level++ )
00533 {
00534 iso.Pop2Ion[ipHE_LIKE][nelem][level] = 0.;
00535
00536
00537
00538
00539 }
00540
00541
00542
00543
00544 ionbal.RateIonizTot[nelem][nelem-1] = 0.;
00545 sum_popn_ov_ion = 0.;
00546
00547
00548
00549 iso.pop_ion_ov_neut[ipHE_LIKE][nelem] = 0.;
00550
00551
00552
00553 ionbal.RateRecomTot[nelem][nelem-1] = 0.;
00554
00555 for( level=0; level < iso.numLevels_local[ipHE_LIKE][nelem]; level++ )
00556 {
00557
00558
00559 ionbal.RateRecomTot[nelem][nelem-1] += iso.RateCont2Level[ipHE_LIKE][nelem][level];
00560
00561
00562 ionbal.RateIonizTot[nelem][nelem-1] +=
00563 iso.Pop2Ion[ipHE_LIKE][nelem][level] * iso.RateLevel2Cont[ipHE_LIKE][nelem][level];
00564
00565
00566
00567 sum_popn_ov_ion += iso.Pop2Ion[ipHE_LIKE][nelem][level];
00568
00569 if( iso.PopLTE[ipHE_LIKE][nelem][level] > 0. )
00570 {
00571
00572 iso.DepartCoef[ipHE_LIKE][nelem][level] =
00573 (iso.Pop2Ion[ipHE_LIKE][nelem][level]/
00574 (iso.PopLTE[ipHE_LIKE][nelem][level]* dense.eden) );
00575 }
00576 else
00577 {
00578 iso.DepartCoef[ipHE_LIKE][nelem][level] = 0.;
00579 }
00580 }
00581
00582
00583
00584 iso.xIonSimple[ipHE_LIKE][nelem] = ionbal.RateIonizTot[nelem][nelem-1] / ionbal.RateRecomTot[nelem][nelem-1];
00585
00586 {
00587
00588 enum {DEBUG_LOC=false};
00589
00590 if( DEBUG_LOC && nzone > 700 && nelem==ipHELIUM )
00591 {
00592 fprintf(ioQQQ,"DEBUG hedest\t%.2f",fnzone);
00593 for( level=0; level < iso.numLevels_local[ipHE_LIKE][nelem]; level++ )
00594 {
00595 if( iso.Pop2Ion[ipHE_LIKE][nelem][level] * iso.RateLevel2Cont[ipHE_LIKE][nelem][level] /
00596 SDIV(ionbal.RateIonizTot[nelem][nelem-1] ) > 0.1 )
00597 fprintf(ioQQQ,"\t%li\t%.3f",level,
00598 iso.Pop2Ion[ipHE_LIKE][nelem][level] * iso.RateLevel2Cont[ipHE_LIKE][nelem][level] /
00599 SDIV(ionbal.RateIonizTot[nelem][nelem-1] ));
00600 }
00601 fprintf(ioQQQ,"\n");
00602 }
00603 }
00604
00605 if( ( ionbal.RateIonizTot[nelem][nelem-1]/MAX2(SMALLDOUBLE , sum_popn_ov_ion) ) > BIGDOUBLE )
00606 {
00607 fprintf( ioQQQ, "RateIonizTot for Z%li, ion%li is larger than BIGDOUBLE. This is a problem.",
00608 nelem+1, nelem-1);
00609 cdEXIT(EXIT_FAILURE);
00610 }
00611 else
00612 ionbal.RateIonizTot[nelem][nelem-1] /= MAX2(SMALLDOUBLE , sum_popn_ov_ion);
00613
00614
00615
00616
00617 if( sum_popn_ov_ion > SMALLDOUBLE )
00618 {
00619 iso.pop_ion_ov_neut[ipHE_LIKE][nelem] = 1./sum_popn_ov_ion;
00620 }
00621 else
00622 {
00623 iso.pop_ion_ov_neut[ipHE_LIKE][nelem] = 0.;
00624 }
00625 }
00626
00627 ASSERT( ionbal.RateIonizTot[nelem][nelem-1] >= 0. );
00628 {
00629
00630 enum {DEBUG_LOC=false};
00631
00632 if( DEBUG_LOC && nelem==1 )
00633 {
00634 fprintf(ioQQQ,"hebugggg2\t%li\t%.2e\t%.2e\t%.2e\n",
00635 nzone,
00636 iso.pop_ion_ov_neut[ipHE_LIKE][nelem],
00637 ionbal.RateIonizTot[nelem][nelem-1] ,
00638 ionbal.RateRecomTot[nelem][nelem-1] );
00639 }
00640 }
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674 ASSERT( conv.lgSearch || helike.lgFSM || dynamics.lgAdvection ||
00675
00676
00677 iso.pop_ion_ov_neut[ipHE_LIKE][nelem]<1e-15 ||
00678 (strcmp(iso.chTypeAtomUsed[ipHE_LIKE][nelem] ,"Low T")==0 && iso.pop_ion_ov_neut[ipHE_LIKE][nelem]<1e-10 ) ||
00679 fabs(iso.pop_ion_ov_neut[ipHE_LIKE][nelem] -
00680 ionbal.RateIonizTot[nelem][nelem-1] /
00681 SDIV(ionbal.RateRecomTot[nelem][nelem-1]) )/
00682 SDIV(iso.pop_ion_ov_neut[ipHE_LIKE][nelem]) < 0.001);
00683
00684
00685 if( lgNegPop || iso.pop_ion_ov_neut[ipHE_LIKE][nelem] < 0.)
00686 {
00687 fprintf( ioQQQ,
00688 " DISASTER HeLikeLevel finds negative He-like ion fraction for nelem=%ld "\
00689 "IonFrac= %.3e solver=%s simple=%.3e TE=%.3e nzone=%4ld\n",
00690 nelem,
00691 iso.pop_ion_ov_neut[ipHE_LIKE][nelem],
00692 iso.chTypeAtomUsed[ipHE_LIKE][nelem] ,
00693 iso.xIonSimple[ipHE_LIKE][nelem],
00694 phycon.te,
00695 nzone );
00696
00697 fprintf( ioQQQ, " level pop are:\n" );
00698 for( i=0; i < iso.numLevels_local[ipHE_LIKE][nelem]; i++ )
00699 {
00700 fprintf( ioQQQ," %8.1e;", iso.Pop2Ion[ipHE_LIKE][nelem][i] );
00701 if( (i!=0) && !(i%10) ) fprintf( ioQQQ,"\n" );
00702 }
00703 fprintf( ioQQQ, "\n" );
00704 ContNegative();
00705 ShowMe();
00706 puts( "[Stop in helike]" );
00707 cdEXIT(EXIT_FAILURE);
00708 }
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719 # if 0
00720 if( false && iso.pop_ion_ov_neut[ipHE_LIKE][nelem] >= 0. &&
00721 iso.pop_ion_ov_neut[ipHE_LIKE][nelem] < TooSmall &&
00722 dense.IonHigh[nelem] == nelem )
00723 {
00724
00725
00726
00727
00728
00729 if( iso.pop_ion_ov_neut[ipHE_LIKE][nelem] == 0. && !dense.lgSetIoniz[nelem] )
00730 {
00731 iso.pop_ion_ov_neut[ipHE_LIKE][nelem] = 0.;
00732
00733
00734 dense.IonHigh[nelem] = nelem-1;
00735
00736
00737
00738 dense.IonLow[nelem] = MIN2( dense.IonLow[nelem] , dense.IonHigh[nelem]-1 );
00739 }
00740
00741
00742 for( ipHi=ipHe2s3S; ipHi < iso.numLevels_local[ipHE_LIKE][nelem]; ipHi++ )
00743 {
00744 iso.Pop2Ion[ipHE_LIKE][nelem][ipHi] = 0.;
00745 for( ipLo=ipHe1s1S; ipLo < ipHi; ipLo++ )
00746 {
00747
00748 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopLo = 0.;
00749
00750
00751 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi = 0.;
00752
00753
00754 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopOpc = 0.;
00755
00756
00757 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].ots = 0.;
00758 }
00759 }
00760
00761 ionbal.RateIonizTot[nelem][nelem-1] = 0.;
00762
00763 if( trace.lgTrace && trace.lgHeBug )
00764 {
00765 fprintf(ioQQQ," HeLike for Z=%li finds small population set to zerozero.\n",nelem);
00766 }
00767 }
00768
00769 else
00770 # endif
00771
00772
00773
00774
00775
00776
00777 for( ipHi=0; ipHi<iso.numLevels_local[ipHE_LIKE][nelem]; ++ipHi )
00778 {
00779 for( ipLo=0; ipLo<ipHi; ++ipLo )
00780 {
00781 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopLo = iso.Pop2Ion[ipHE_LIKE][nelem][ipLo];
00782 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopHi = iso.Pop2Ion[ipHE_LIKE][nelem][ipHi];
00783 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopOpc = (iso.Pop2Ion[ipHE_LIKE][nelem][ipLo] - iso.Pop2Ion[ipHE_LIKE][nelem][ipHi]*
00784 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gLo/EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].gHi);
00785
00786 if( ipHi < iso.numLevels_local[ipHE_LIKE][nelem] - 1 )
00787 EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopOpc = MAX2(SMALLDOUBLE,EmisLines[ipHE_LIKE][nelem][ipHi][ipLo].PopOpc );
00788 }
00789 }
00790
00791 if( !helike.lgSetBenjamin )
00792 {
00793
00794
00795
00796
00797 for( ipHi=ipHe2p1P; ipHi < iso.numLevels_local[ipHE_LIKE][nelem]; ipHi++ )
00798 {
00799 double sum;
00800 long OneOfDegenerate, OtherDegenerate;
00801
00802 OneOfDegenerate = ipHe2p3P1;
00803
00804 if( nelem == ipHELIUM )
00805 OtherDegenerate = ipHe2p3P2;
00806 else
00807 OtherDegenerate = ipHe2p3P0;
00808
00809 sum = EmisLines[ipHE_LIKE][nelem][ipHi][OneOfDegenerate].PopOpc +
00810 EmisLines[ipHE_LIKE][nelem][ipHi][OtherDegenerate].PopOpc;
00811
00812 EmisLines[ipHE_LIKE][nelem][ipHi][OneOfDegenerate].PopOpc = sum;
00813 EmisLines[ipHE_LIKE][nelem][ipHi][OtherDegenerate].PopOpc = sum;
00814 }
00815 }
00816
00817
00818
00819
00820 if( ( (HighestPColOut[0] < 1./helike.Lifetime[ipHELIUM][QuantumNumbers2Index[ipHELIUM][iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM]][1][0]] )
00821 || (HighestPColOut[1] < 1./helike.Lifetime[ipHELIUM][QuantumNumbers2Index[ipHELIUM][iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM]][1][1]] ) )
00822 && nelem==ipHELIUM )
00823 {
00824 helike.lgCritDensLMix = false;
00825 }
00826
00827 for( n=0; n < iso.numLevels_local[ipHE_LIKE][nelem]; n++ )
00828 {
00829 OutOfNGoingUp[iso.quant_desig[ipHE_LIKE][nelem][n].n] +=
00830 CollisionsGoingUp[n]*iso.stat[ipHE_LIKE][nelem][n]/(4.*N_(n)*N_(n));
00831 OutOfNGoingDown[iso.quant_desig[ipHE_LIKE][nelem][n].n] +=
00832 CollisionsGoingDown[n]*iso.stat[ipHE_LIKE][nelem][n]/(4.*N_(n)*N_(n));
00833 }
00834
00835 #if 0
00836
00837 if( iso.lgColl_excite[ipHE_LIKE] && !helike.lgSetBenjamin && helike.lgCS_Vriens
00838 && !helike.lgRandErrGen && dense.eden > 100. )
00839 {
00840 for( n = 4; n < iso.n_HighestResolved_local[ipHE_LIKE][nelem]; n++ )
00841 {
00842 double DecayNPlusOne = 1./helike.Lifetime[nelem][QuantumNumbers2Index[nelem][n+1][1][0]];
00843
00844
00845
00846
00847
00848
00849
00850
00851 if( (OutOfNGoingUp[n]>(1E-8*DecayNPlusOne)) && n!=iso.n_HighestResolved_local[ipHE_LIKE][nelem] )
00852 {
00853
00854 ASSERT( (OutOfNGoingUp[n-1]<OutOfNGoingUp[n]) );
00855 ASSERT( (OutOfNGoingDown[n]<OutOfNGoingDown[n+1]) );
00856 }
00857 }
00858 }
00859 #endif
00860
00861 {
00862
00863 enum {DEBUG_LOC=false};
00864
00865 if( DEBUG_LOC && nelem == ipHELIUM )
00866 {
00867
00868 static long RUNTHISONCE = false;
00869
00870 if( !RUNTHISONCE )
00871 {
00872 #if 0
00873 fprintf( ioQQQ,"n\tl\ts\tCollisionsGoingUp\tCollisionsGoingDown\n");
00874
00875 for( n=0; n < iso.numLevels_local[ipHE_LIKE][nelem]; n++ )
00876 {
00877 fprintf( ioQQQ,"%2ld\t%2ld\t%2ld\t%.9e\t%.9e\n",
00878 iso.quant_desig[ipHE_LIKE][nelem][n].n,
00879 iso.quant_desig[ipHE_LIKE][nelem][n].l,
00880 iso.quant_desig[ipHE_LIKE][nelem][n].s,
00881 CollisionsGoingUp[n]*iso.Pop2Ion[ipHE_LIKE][nelem][n],
00882 CollisionsGoingDown[n]*iso.Pop2Ion[ipHE_LIKE][nelem][n] );
00883 }
00884 #endif
00885
00886 fprintf( ioQQQ,"n\tOutOfNGoingUp\tOutOfNGoingDown\n");
00887
00888 for( n=1; n <= iso.n_HighestResolved_local[ipHE_LIKE][nelem] + iso.nCollapsed_local[ipHE_LIKE][nelem]; n++ )
00889 {
00890 fprintf( ioQQQ,"%2ld\t%.9e\t%.9e\n",
00891 n, OutOfNGoingUp[n], OutOfNGoingDown[n] );
00892 }
00893
00894 fprintf( ioQQQ, "\n" );
00895 }
00896
00897 RUNTHISONCE = true;
00898 }
00899 }
00900
00901 {
00902
00903 enum {DEBUG_LOC=false};
00904
00905 if( DEBUG_LOC && nelem == ipHELIUM )
00906 {
00907
00908
00909 static long RUNONCE = false, nlo, nhi;
00910
00911
00912 if( !RUNONCE )
00913 {
00914 if( (ioMATRIX = fopen("RateMatrix.txt","w")) == NULL )
00915 return;
00916
00917 #if 0
00918 zPerN = (double **)MALLOC(sizeof(double*)*(unsigned) (iso.n_HighestResolved_local[ipHE_LIKE][nelem] + iso.nCollapsed_local[ipHE_LIKE][nelem] + 1) );
00919
00920 for( i=0; i <= (iso.n_HighestResolved_local[ipHE_LIKE][nelem] + iso.nCollapsed_local[ipHE_LIKE][nelem]); ++i )
00921 {
00922 zPerN[i] = (double *)MALLOC(sizeof(double)*(unsigned) (iso.n_HighestResolved_local[ipHE_LIKE][nelem] + iso.nCollapsed_local[ipHE_LIKE][nelem] + 1) );
00923 }
00924
00925
00926 for( nlo = 1; nlo <= iso.n_HighestResolved_local[ipHE_LIKE][nelem] + iso.nCollapsed_local[ipHE_LIKE][nelem]; nlo++ )
00927 {
00928 for( nhi = 1; nhi <= iso.n_HighestResolved_max[ipHE_LIKE][ipHELIUM] + iso.nCollapsed_max[ipHE_LIKE][ipHELIUM]; nhi++ )
00929 {
00930 zPerN[nlo][nhi] = 0.;
00931 }
00932 }
00933
00934
00935 for( level=ipHe1s1S; level < iso.numLevels_local[ipHE_LIKE][nelem]; level++ )
00936 {
00937 for( ipLo=ipHe1s1S; ipLo < level; ipLo++ )
00938 {
00939 zPerN[iso.quant_desig[ipHE_LIKE][nelem][ipLo].n][iso.quant_desig[ipHE_LIKE][nelem][level].n]
00940 += z[ipLo][level]*iso.Pop2Ion[ipHE_LIKE][nelem][ipLo];
00941 }
00942
00943 for( ipHi=level; ipHi < iso.numLevels_local[ipHE_LIKE][nelem]; ipHi++ )
00944 {
00945 zPerN[iso.quant_desig[ipHE_LIKE][nelem][ipHi].n][iso.quant_desig[ipHE_LIKE][nelem][level].n]
00946 += z[ipHi][level]*iso.Pop2Ion[ipHE_LIKE][nelem][ipHi];
00947 }
00948 }
00949
00950
00951
00952 for( nlo = 1; nlo <= iso.n_HighestResolved_local[ipHE_LIKE][nelem] + iso.nCollapsed_local[ipHE_LIKE][nelem]; nlo++ )
00953 {
00954 for( nhi = 1; nhi <= iso.n_HighestResolved_local[ipHE_LIKE][nelem] + iso.nCollapsed_local[ipHE_LIKE][nelem]; nhi++ )
00955 {
00956 double f=fabs(zPerN[nlo][nhi]);
00957 fprintf( ioMATRIX, "%.2e\t", log10(SDIV(f)) );
00958 }
00959 fprintf( ioMATRIX, "\n" );
00960 }
00961 #endif
00962
00963 for( nlo = 0; nlo < iso.numLevels_local[ipHE_LIKE][nelem]; nlo++ )
00964 {
00965 for( nhi = 0; nhi < iso.numLevels_local[ipHE_LIKE][nelem]; nhi++ )
00966 {
00967 fprintf( ioMATRIX, "%.2e\t", z[nlo][nhi] );
00968 #if 0
00969
00970 if( z[nlo][nhi] == SaveZ[nlo][nhi] )
00971 fprintf( ioMATRIX, "%.2e\t", 1.0 );
00972 else if( SaveZ[nlo][nhi] == 0. )
00973 fprintf( ioMATRIX, "%.2e\t", 0.0 );
00974 else
00975 fprintf( ioMATRIX, "%.2e\t", z[nlo][nhi]/SaveZ[nlo][nhi] );
00976
00977 #endif
00978 }
00979 fprintf( ioMATRIX, "\n" );
00980 }
00981
00982 fprintf( ioMATRIX, "\n" );
00983 fclose( ioMATRIX );
00984 }
00985
00986 RUNONCE = true;
00987 }
00988 }
00989
00990 {
00991
00992
00993 enum {DEBUG_LOC=false};
00994
00995 if( DEBUG_LOC )
00996 {
00997 if( nelem==ipHELIUM )
00998 fprintf(ioQQQ,"\n");
00999 fprintf(ioQQQ,"%li\t%.4e", nelem,
01000 iso.pop_ion_ov_neut[ipHE_LIKE][nelem] );
01001 if( dense.xIonDense[nelem][nelem]>SMALLFLOAT )
01002 {
01003 float ratio = dense.xIonDense[nelem][nelem]/ dense.xIonDense[nelem][nelem-1];
01004
01005 fprintf(ioQQQ,"\t%.4e\t%.4e", ratio ,
01006 (ratio-iso.pop_ion_ov_neut[ipHE_LIKE][nelem])/ratio );
01007 }
01008 fprintf(ioQQQ,"\n");
01009 }
01010
01011 }
01012 {
01013
01014
01015 enum {DEBUG_LOC=false};
01016
01017 if( DEBUG_LOC && (nelem==ipCARBON) )
01018 {
01019 {
01020 fprintf( ioQQQ, " He-like Z=%2ld H2OVH1=",nelem);
01021 fprintf( ioQQQ,"%10.3e", iso.pop_ion_ov_neut[ipHE_LIKE][nelem]);
01022 fprintf( ioQQQ, " simple=");
01023 fprintf( ioQQQ,"%10.3e", iso.xIonSimple[ipHE_LIKE][nelem]);
01024 fprintf( ioQQQ," dest, creat rates %.2e %.2e", ionbal.RateIonizTot[nelem][nelem-1] , ionbal.RateRecomTot[nelem][nelem-1]);
01025 fprintf( ioQQQ, "\n");
01026 }
01027 }
01028 }
01029
01030
01031
01032
01033
01034
01035
01036
01037 if( nzone>0 && nelem==ipHELIUM )
01038 {
01039
01040 double ratio = iso.Pop2Ion[ipHE_LIKE][nelem][ipHe2s3S] * iso.RateLevel2Cont[ipHE_LIKE][nelem][ipHe2s3S] /
01041 SDIV(iso.Pop2Ion[ipHE_LIKE][nelem][ipHe1s1S]*ionbal.RateIonizTot[nelem][nelem-1]);
01042 if( ratio > he.frac_he0dest_23S )
01043 {
01044
01045 he.nzone = nzone;
01046 he.frac_he0dest_23S = (float)ratio;
01047 he.frac_he0dest_23S_photo = (float)(iso.Pop2Ion[ipHE_LIKE][nelem][ipHe2s3S] *iso.gamnc[ipHE_LIKE][nelem][ipHe2s3S]/
01048 SDIV(iso.Pop2Ion[ipHE_LIKE][nelem][ipHe1s1S]*ionbal.RateIonizTot[nelem][nelem-1]));
01049 }
01050 }
01051
01052
01053
01054
01055
01056 if( trace.lgTrace )
01057 {
01058 fprintf( ioQQQ, " He-like %.2f Z:%2ld H2OVH1:",fnzone,nelem);
01059 fprintf( ioQQQ,"%10.3e", iso.pop_ion_ov_neut[ipHE_LIKE][nelem]);
01060 fprintf( ioQQQ, " simple: %.3e", iso.xIonSimple[ipHE_LIKE][nelem]);
01061 fprintf( ioQQQ," rate dest:%.2e creat:%.2e", SaveZ[0][0] , ionbal.RateRecomTot[nelem][nelem-1]);
01062 fprintf(ioQQQ," photo:%.3e coll:%.3e sec:%.3e ",
01063 iso.gamnc[ipHE_LIKE][nelem][0],
01064 iso.ColIoniz[ipHE_LIKE][nelem][0]* dense.EdenHCorr,
01065 secondaries.csupra[nelem][nelem-1]);
01066 fprintf( ioQQQ, "\n");
01067 }
01068
01069
01070 if( PARALLEL )
01071 {
01072 free( ipiv );
01073 free( creation );
01074 free( OutOfNGoingUp );
01075 free( OutOfNGoingDown );
01076 free( CollisionsGoingUp );
01077 free( CollisionsGoingDown );
01078 free( error );
01079 free( work );
01080
01081
01082 for( i=0; i<iso.numLevels_local[ipHE_LIKE][nelem]; ++i )
01083 {
01084 free( SaveZ[i] );
01085
01086 free( z[i] );
01087 }
01088
01089 free( SaveZ );
01090 free( z );
01091 }
01092
01093 DEBUG_EXIT( "HeLikeLevel()" );
01094
01095 return;
01096 }
01097
01098 void HeLikeError( long nelem )
01099 {
01100
01101 long ipHi, ipLo, typeOfRate;
01102 float ErrorToPut;
01103
01104
01105 for( ipHi=ipHe1s1S; ipHi<iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
01106 {
01107 if( N_(ipHi)<=5 )
01108 {
01109 ErrorToPut = 0.02f;
01110 }
01111 else if( N_(ipHi)>iso.n_HighestResolved_max[ipHE_LIKE][nelem] )
01112 {
01113 ErrorToPut = 0.02f;
01114 }
01115 else if( L_(ipHi)>=3 )
01116 {
01117 ErrorToPut = 0.005f;
01118 }
01119 else if( L_(ipHi)==2 )
01120 {
01121 ErrorToPut = 0.01f;
01122 }
01123 else if( L_(ipHi)==1 )
01124 {
01125 ErrorToPut = 0.025f;
01126 }
01127 else
01128 {
01129 ErrorToPut = 0.05f;
01130 }
01131 putError(nelem,iso.numLevels_max[ipHE_LIKE][nelem],ipHi,IPRAD,ErrorToPut);
01132 }
01133
01134 putError(nelem,iso.numLevels_max[ipHE_LIKE][nelem],iso.numLevels_max[ipHE_LIKE][nelem],IPRAD,0.02f);
01135
01136 for( ipHi=ipHe2s3S; ipHi<= iso.numLevels_max[ipHE_LIKE][nelem]; ipHi++ )
01137 {
01138
01139 for( ipLo=ipHe1s1S; ipLo<= ipHi; ipLo++ )
01140 {
01141 for( typeOfRate=0; typeOfRate<=1; typeOfRate++ )
01142 {
01143 if( helike.lgHugeCaseB == true )
01144 {
01145
01146
01147 helike.ErrorFactor[nelem][ipHi][ipLo][typeOfRate] = 1.0f;
01148 }
01149 #if 0
01150 else if( typeOfRate == 0 )
01151 {
01153
01154
01155 helike.ErrorFactor[nelem][ipHi][ipLo][typeOfRate] = 1.0f;
01156 }
01157 #endif
01158 else if( helike.Error[nelem][ipHi][ipLo][typeOfRate] >= 0. )
01159 {
01160 helike.ErrorFactor[nelem][ipHi][ipLo][typeOfRate] =
01161 (float)MyGaussRand( helike.Error[nelem][ipHi][ipLo][typeOfRate] );
01162 ASSERT( helike.ErrorFactor[nelem][ipHi][ipLo][typeOfRate] > 0. );
01163 }
01164 else
01165 {
01166 helike.ErrorFactor[nelem][ipHi][ipLo][typeOfRate] = 1.0f;
01167 }
01168 }
01169 }
01170 }
01171
01172 return;
01173 }
01174
01175
01176