00001
00002
00003
00004 #include "cddefines.h"
00005 #include "atmdat.h"
00006 #include "continuum.h"
00007 #include "conv.h"
00008 #include "dense.h"
00009 #include "dynamics.h"
00010 #include "elementnames.h"
00011 #include "grainvar.h"
00012 #include "he.h"
00013 #include "helike.h"
00014 #include "hmi.h"
00015 #include "hydrogenic.h"
00016 #include "ionbal.h"
00017 #include "iso.h"
00018 #include "mole.h"
00019 #include "phycon.h"
00020 #include "physconst.h"
00021 #include "rfield.h"
00022 #include "secondaries.h"
00023 #include "taulines.h"
00024 #include "thirdparty.h"
00025 #include "trace.h"
00026
00027
00028 void iso_level( const long int ipISO, const long int nelem)
00029 {
00030 long int ipHi,
00031 ipLo,
00032 i,
00033 level,
00034 level_error;
00035
00036 const long int numlevels_local = iso.numLevels_local[ipISO][nelem];
00037
00038 double BigError;
00039
00040 int32 nerror;
00041 double HighestPColOut = 0.,
00042 TotalPop;
00043 bool lgNegPop=false;
00044 valarray<int32> ipiv(numlevels_local);
00045
00046 valarray<double>
00047 creation(numlevels_local),
00048 error(numlevels_local),
00049 work(numlevels_local),
00050 Save_creation(numlevels_local);
00051 double source=0.,
00052 sink=0.;
00053 valarray<double> PopPerN(iso.n_HighestResolved_local[ipISO][nelem]+1);
00054
00055 multi_arr<double,2,C_TYPE> z, SaveZ;
00056
00057 DEBUG_ENTRY( "iso_level()" );
00058
00059
00060 ASSERT( nelem >= ipISO );
00061 ASSERT( nelem < LIMELM );
00062
00063
00064 z.alloc(numlevels_local,numlevels_local);
00065
00066
00067 for( level=0; level < numlevels_local; level++ )
00068 {
00069 ASSERT( dense.xIonDense[nelem][nelem+1-ipISO] >= 0.f );
00070
00071 creation[level] = iso.RateCont2Level[ipISO][nelem][level] * dense.xIonDense[nelem][nelem+1-ipISO];
00072 }
00073
00074
00075
00076 ASSERT( ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0]< SMALLFLOAT ||
00077 fabs( (iso.ColIoniz[ipISO][nelem][0]* dense.EdenHCorr) /
00078 SDIV(ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0] ) - 1.) < 0.001 );
00079
00080
00081 if( dense.xIonDense[nelem][nelem+1-ipISO] < SMALLFLOAT || iso.xIonSimple[ipISO][nelem] < 1e-35 )
00082 {
00083
00084 strcpy( iso.chTypeAtomUsed[ipISO][nelem], "zero " );
00085 if( trace.lgTrace && (nelem == trace.ipIsoTrace[ipISO]) )
00086 {
00087 fprintf( ioQQQ, " iso_level iso=%2ld nelem=%2ld simple II/I=%10.2e so not doing equilibrium, doing %s.\n",
00088 ipISO, nelem, iso.xIonSimple[ipISO][nelem], iso.chTypeAtomUsed[ipISO][nelem] );
00089 }
00090
00091
00092 ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] = iso.RateLevel2Cont[ipISO][nelem][0];
00093 lgNegPop = false;
00094 StatesElemNEW[nelem][nelem-ipISO][0].Pop = dense.xIonDense[nelem][nelem-ipISO];
00095
00096 for( long n=1; n < numlevels_local; n++ )
00097 {
00098 StatesElemNEW[nelem][nelem-ipISO][n].Pop = 0.;
00099 }
00100 iso.qTot2S[ipISO][nelem] = 0.;
00101 }
00102 else
00103 {
00104
00105
00106 ASSERT( NISO == 2 );
00107 long SpinUsed[NISO] = {2,1};
00108 long indexNmaxP =
00109 iso.QuantumNumbers2Index[ipISO][nelem][ iso.n_HighestResolved_local[ipISO][nelem] ][1][SpinUsed[ipISO]];
00110
00111 strcpy( iso.chTypeAtomUsed[ipISO][nelem], "Popul" );
00112 if( trace.lgTrace && (nelem == trace.ipIsoTrace[ipISO]) )
00113 {
00114 fprintf( ioQQQ, " iso_level iso=%2ld nelem=%2ld doing regular matrix inversion, %s\n",
00115 ipISO, nelem, iso.chTypeAtomUsed[ipISO][nelem] );
00116 }
00117
00118 multi_arr<quantumState,3>::const_iterator StElm = StatesElemNEW.begin(nelem,nelem-ipISO);
00119
00120
00121 for( level=0; level < numlevels_local; level++ )
00122 {
00123
00124
00125 z[level][level] = iso.RateLevel2Cont[ipISO][nelem][level];
00126
00127 if( level == 1 )
00128
00129 iso.qTot2S[ipISO][nelem] = iso.ColIoniz[ipISO][nelem][level];
00130
00131 multi_arr<transition,4>::const_iterator Trans = Transitions.begin(ipISO,nelem,level);
00132 md4ci Boltz = iso.Boltzmann.begin(ipISO,nelem,level);
00133
00134
00135 for( ipLo=0; ipLo < level; ipLo++ )
00136 {
00137 double coll_down = Trans[ipLo].Coll.ColUL * dense.EdenHCorr;
00138
00139 double RadDecay = MAX2( iso.SmallA, Trans[ipLo].Emis->Aul*
00140 (Trans[ipLo].Emis->Pesc +
00141 Trans[ipLo].Emis->Pelec_esc +
00142 Trans[ipLo].Emis->Pdest)*
00143 KILL_BELOW_PLASMA(Trans[ipLo].EnergyWN*WAVNRYD) );
00144
00145 double pump = MAX2( iso.SmallA, Trans[ipLo].Emis->pump *
00146 KILL_BELOW_PLASMA(Trans[ipLo].EnergyWN*WAVNRYD) );
00147
00148 if( iso.lgRandErrGen[ipISO] )
00149 {
00150 coll_down *= iso.ErrorFactor[ipISO][nelem][level][ipLo][IPCOLLIS];
00151 RadDecay *= iso.ErrorFactor[ipISO][nelem][level][ipLo][IPRAD];
00152 pump *= iso.ErrorFactor[ipISO][nelem][level][ipLo][IPRAD];
00153 }
00154
00155 double coll_up = coll_down *
00156 (double)StElm[level].g/
00157 (double)StElm[ipLo].g*
00158 Boltz[ipLo];
00159
00160 z[ipLo][ipLo] += coll_up + pump ;
00161 z[ipLo][level] = - ( coll_up + pump );
00162
00163 double pump_down = pump *
00164 (double)StElm[ipLo].g/
00165 (double)StElm[level].g;
00166
00167 z[level][level] += RadDecay + pump_down + coll_down;
00168 z[level][ipLo] = - (RadDecay + pump_down + coll_down);
00169
00170 if( level == indexNmaxP )
00171 {
00172 HighestPColOut += coll_down;
00173 }
00174 if( ipLo == indexNmaxP )
00175 {
00176 HighestPColOut += coll_up;
00177 }
00178
00179
00180 if( (level == 1) && (ipLo==0) )
00181 {
00182 iso.qTot2S[ipISO][nelem] += coll_down/dense.EdenHCorr;
00183 }
00184 if( (ipLo == 1) && (ipISO==ipH_LIKE || StElm[level].S==1) )
00185 {
00186 iso.qTot2S[ipISO][nelem] += coll_up/dense.EdenHCorr;
00187 }
00188 }
00189 }
00190
00191 if( ipISO == nelem )
00192 {
00193
00194
00195
00196 if( HighestPColOut < 1./StatesElemNEW[nelem][nelem-ipISO][indexNmaxP].lifetime )
00197 {
00198 iso.lgCritDensLMix[ipISO] = false;
00199 }
00200 }
00201
00203 ASSERT( ipISO <= ipHE_LIKE );
00204
00205
00206
00207
00208 z[1+ipISO][0] -= iso.TwoNu_induc_dn[ipISO][nelem]*iso.lgInd2nu_On;
00209 z[0][1+ipISO] -= iso.TwoNu_induc_up[ipISO][nelem]*iso.lgInd2nu_On;
00210
00211
00212 z[0][0] += iso.TwoNu_induc_up[ipISO][nelem]*iso.lgInd2nu_On;
00213 z[1+ipISO][1+ipISO] += iso.TwoNu_induc_dn[ipISO][nelem]*iso.lgInd2nu_On;
00214
00215
00216 for( long ion=0; ion<=nelem+1; ++ion )
00217 {
00218 if( ion!=nelem-ipISO )
00219 {
00220 source += gv.GrainChTrRate[nelem][ion][nelem-ipISO] *
00221 dense.xIonDense[nelem][ion];
00222 sink += gv.GrainChTrRate[nelem][nelem-ipISO][ion];
00223 #if 0
00224 source += mole.xMoleChTrRate[nelem][ion][nelem-ipISO] *
00225 dense.xIonDense[nelem][ion] * atmdat.lgCTOn;
00226 sink += mole.xMoleChTrRate[nelem][nelem-ipISO][ion] * atmdat.lgCTOn;
00227 #endif
00228 }
00229 }
00230
00231 #if 0
00232
00233
00234 if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection &&
00235 dynamics.Rate != 0.0 &&
00236 !dynamics.lgEquilibrium && dynamics.lgISO[ipISO])
00237 {
00238
00239 source += dynamics.Source[nelem][nelem-ipISO];
00240
00241 sink += dynamics.Rate;
00242 }
00243 #else
00244
00245 if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection &&
00246 dynamics.Rate != 0.0 &&
00247 !dynamics.lgEquilibrium && dynamics.lgISO[ipISO])
00248 {
00249 for( level=0; level < numlevels_local; level++ )
00250 {
00251 creation[level] += dynamics.StatesElemNEW[nelem][nelem-ipISO][level];
00252 z[level][level] += dynamics.Rate;
00253 }
00254 }
00255 #endif
00256
00257
00258 for(long ion_from=dense.IonLow[nelem]; ion_from < MIN2( dense.IonHigh[nelem], nelem-ipISO ) ; ion_from++ )
00259 {
00260 if( ionbal.RateIoniz[nelem][ion_from][nelem-ipISO] >= 0. )
00261 {
00262 double damper = sexp( dense.gas_phase[nelem] / 1e10 );
00263
00264
00265 source += ionbal.RateIoniz[nelem][ion_from][nelem-ipISO] * dense.xIonDense[nelem][ion_from] * damper;
00266
00267 if( ion_from == nelem-1-ipISO )
00268 sink += ionbal.RateRecomTot[nelem][ion_from] * damper;
00269 }
00270 }
00271
00272 ASSERT( source >= 0.f );
00273 creation[0] += source;
00274 for( level=0; level < numlevels_local; level++ )
00275 {
00276 z[level][level] += sink;
00277 }
00278
00279
00280 if( secondaries.Hx12[ipISO][nelem][iso.nLyaLevel[ipISO]] * iso.lgColl_excite[ipISO] > 0. )
00281 {
00282
00283 for( level=1; level < numlevels_local; level++ )
00284 {
00285 double RateUp , RateDown;
00286
00287 RateUp = secondaries.Hx12[ipISO][nelem][level];
00288 RateDown = RateUp * (double)StatesElemNEW[nelem][nelem-ipISO][ipH1s].g /
00289 (double)StatesElemNEW[nelem][nelem-ipISO][level].g;
00290
00291
00292 z[ipH1s][ipH1s] += RateUp;
00293
00294
00295 z[level][ipH1s] -= RateDown;
00296
00297
00298 z[ipH1s][level] -= RateUp;
00299
00300 z[level][level] += RateDown;
00301 }
00302 }
00303
00304
00305
00306
00307
00308
00309
00310 SaveZ = z;
00311
00312 for( ipLo=0; ipLo < numlevels_local; ipLo++ )
00313 Save_creation[ipLo] = creation[ipLo];
00314
00315 if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) )
00316 {
00317 fprintf( ioQQQ, " pop level others => (iso_level)\n" );
00318 for( long n=0; n < numlevels_local; n++ )
00319 {
00320 fprintf( ioQQQ, " %s %s %2ld", iso.chISO[ipISO], elementnames.chElementNameShort[nelem], n );
00321 for( long j=0; j < numlevels_local; j++ )
00322 {
00323 fprintf( ioQQQ,"\t%.9e", z[j][n] );
00324 }
00325 fprintf( ioQQQ, "\n" );
00326 }
00327 fprintf(ioQQQ," recomb ");
00328 for( long n=0; n < numlevels_local; n++ )
00329 {
00330 fprintf( ioQQQ,"\t%.9e", creation[n] );
00331 }
00332 fprintf( ioQQQ, "\n" );
00333 fprintf(ioQQQ," recomb ct %.2e co %.2e hectr %.2e hctr %.2e\n",
00334 atmdat.HeCharExcRecTotal,
00335 findspecies("CO")->hevmol ,
00336 atmdat.HeCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHELIUM][0] ,
00337 atmdat.HCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHYDROGEN][0] );
00338 }
00339
00340 nerror = 0;
00341
00342 getrf_wrapper(numlevels_local,numlevels_local,
00343 z.data(),numlevels_local,&ipiv[0],&nerror);
00344
00345 getrs_wrapper('N',numlevels_local,1,z.data(),numlevels_local,&ipiv[0],
00346 &creation[0],numlevels_local,&nerror);
00347
00348 if( nerror != 0 )
00349 {
00350 fprintf( ioQQQ, " iso_level: dgetrs finds singular or ill-conditioned matrix\n" );
00351 cdEXIT(EXIT_FAILURE);
00352 }
00353
00354
00355
00356 for( level=ipH1s; level < numlevels_local; level++ )
00357 {
00358 double qn = 0., qx = 0.;
00359 error[level] = 0.;
00360 for( long n=ipH1s; n < numlevels_local; n++ )
00361 {
00362 double q = SaveZ[n][level]*creation[n];
00363
00364
00365 if ( q > qx )
00366 qx = q;
00367 else if (q < qn)
00368 qn = q;
00369
00370 error[level] += q;
00371 }
00372
00373 if (-qn > qx)
00374 qx = -qn;
00375
00376 if( qx > 0. )
00377 {
00378 error[level] = (error[level] - Save_creation[level])/qx;
00379 }
00380 else
00381 {
00382 error[level] = 0.;
00383 }
00384 }
00385
00386
00387 BigError = -1.;
00388 level_error = -1;
00389
00390 for( level=ipH1s; level < numlevels_local; level++ )
00391 {
00392 double abserror;
00393 abserror = fabs( error[level]);
00394
00395 if( abserror > BigError )
00396 {
00397 BigError = abserror;
00398 level_error = level;
00399 }
00400 }
00401
00402
00403
00404 if( BigError > FLT_EPSILON )
00405 {
00406 if( !conv.lgSearch )
00407 fprintf(ioQQQ," PROBLEM" );
00408
00409 fprintf(ioQQQ,
00410 " iso_level zone %.2f - largest residual in iso=%li %s nelem=%li matrix inversion is %g "
00411 "level was %li Search?%c \n",
00412 fnzone,
00413 ipISO,
00414 elementnames.chElementName[nelem],
00415 nelem ,
00416 BigError ,
00417 level_error,
00418 TorF(conv.lgSearch) );
00419 }
00420
00421
00422 for( level=0; level < numlevels_local; level++ )
00423 {
00424 StatesElemNEW[nelem][nelem-ipISO][level].Pop = creation[level];
00425
00426
00427
00428
00429 if( StatesElemNEW[nelem][nelem-ipISO][level].Pop < 0. )
00430 lgNegPop = true;
00431
00432 if( StatesElemNEW[nelem][nelem-ipISO][level].Pop <= 0 && !conv.lgSearch )
00433 {
00434 fprintf(ioQQQ,
00435 "PROBLEM non-positive level pop for iso = %li, nelem = "
00436 "%li = %s, level=%li val=%.3e nTotalIoniz %li\n",
00437 ipISO,
00438 nelem ,
00439 elementnames.chElementSym[nelem],
00440 level,
00441 StatesElemNEW[nelem][nelem-ipISO][level].Pop ,
00442 conv.nTotalIoniz);
00443 }
00444
00445 }
00446
00447
00448 for( level=numlevels_local; level < iso.numLevels_max[ipISO][nelem]; level++ )
00449 {
00450 StatesElemNEW[nelem][nelem-ipISO][level].Pop = 0.;
00451
00452
00453 }
00454
00455
00456 TotalPop = 0.;
00457
00458 for( level=0; level < numlevels_local; level++ )
00459 TotalPop += StatesElemNEW[nelem][nelem-ipISO][level].Pop;
00460
00461
00462 if( dense.lgSetIoniz[nelem] )
00463 {
00464 lgNegPop = false;
00465 for( level=0; level < numlevels_local; level++ )
00466 {
00467 StatesElemNEW[nelem][nelem-ipISO][level].Pop *=
00468 dense.xIonDense[nelem][nelem-ipISO]/TotalPop;
00469 lgNegPop = lgNegPop || (StatesElemNEW[nelem][nelem-ipISO][level].Pop < 0.);
00470 }
00471 TotalPop = dense.xIonDense[nelem][nelem-ipISO];
00472 }
00473
00474 ASSERT( TotalPop >= 0. );
00475
00476
00477
00478 ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] = 0.;
00479 for( level=0; level < numlevels_local; level++ )
00480 {
00481
00482
00483 ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] +=
00484 StatesElemNEW[nelem][nelem-ipISO][level].Pop * iso.RateLevel2Cont[ipISO][nelem][level];
00485 }
00486
00487 if( ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] > BIGDOUBLE )
00488 {
00489 fprintf( ioQQQ, "DISASTER RateIonizTot for Z=%li, ion %li is larger than BIGDOUBLE. This is a big problem.",
00490 nelem+1, nelem-ipISO);
00491 cdEXIT(EXIT_FAILURE);
00492 }
00493 else
00494 ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] /= MAX2(SMALLDOUBLE , TotalPop);
00495 }
00496
00497
00498
00499 if( ionbal.RateRecomTot[nelem][nelem-ipISO] > 0. )
00500 iso.xIonSimple[ipISO][nelem] = ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1]/ionbal.RateRecomTot[nelem][nelem-ipISO];
00501 else
00502 iso.xIonSimple[ipISO][nelem] = 0.;
00503
00504 ASSERT( ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] >= 0. );
00505
00506
00507 if( lgNegPop || dense.xIonDense[nelem][nelem-ipISO] < 0. )
00508 {
00509 fprintf( ioQQQ,
00510 " DISASTER iso_level finds negative ion fraction for iso=%2ld nelem=%2ld "\
00511 "%s using solver %s, IonFrac=%.3e simple=%.3e TE=%.3e ZONE=%4ld\n",
00512 ipISO,
00513 nelem,
00514 elementnames.chElementSym[nelem],
00515 iso.chTypeAtomUsed[ipISO][nelem],
00516 dense.xIonDense[nelem][nelem+1-ipISO]/SDIV(dense.xIonDense[nelem][nelem-ipISO]),
00517 iso.xIonSimple[ipISO][nelem],
00518 phycon.te,
00519 nzone );
00520
00521 fprintf( ioQQQ, " level pop are:\n" );
00522 for( i=0; i < numlevels_local; i++ )
00523 {
00524 fprintf( ioQQQ,PrintEfmt("%8.1e", StatesElemNEW[nelem][nelem-ipISO][i].Pop ));
00525 if( (i!=0) && !(i%10) ) fprintf( ioQQQ,"\n" );
00526 }
00527 fprintf( ioQQQ, "\n" );
00528 ContNegative();
00529 ShowMe();
00530 cdEXIT(EXIT_FAILURE);
00531 }
00532
00533 if( ipISO == ipHE_LIKE && nelem==ipHELIUM && nzone>0 )
00534 {
00535
00536 double ratio;
00537 double rateOutOf2TripS = StatesElemNEW[nelem][nelem-ipISO][ipHe2s3S].Pop * iso.RateLevel2Cont[ipISO][nelem][ipHe2s3S];
00538 if( rateOutOf2TripS > SMALLFLOAT )
00539 {
00540 ratio = rateOutOf2TripS /
00541 ( rateOutOf2TripS + StatesElemNEW[nelem][nelem-ipISO][ipHe1s1S].Pop*ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] );
00542 }
00543 else
00544 {
00545 ratio = 0.;
00546 }
00547 if( ratio > he.frac_he0dest_23S )
00548 {
00549
00550 he.nzone = nzone;
00551 he.frac_he0dest_23S = ratio;
00552 rateOutOf2TripS = StatesElemNEW[nelem][nelem-ipISO][ipHe2s3S].Pop *iso.gamnc[ipISO][nelem][ipHe2s3S];
00553 if( rateOutOf2TripS > SMALLFLOAT )
00554 {
00555 he.frac_he0dest_23S_photo = rateOutOf2TripS /
00556 ( rateOutOf2TripS + StatesElemNEW[nelem][nelem-ipISO][ipHe1s1S].Pop*ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] );
00557 }
00558 else
00559 {
00560 he.frac_he0dest_23S_photo = 0.;
00561 }
00562 }
00563 }
00564
00565 for( ipHi=1; ipHi<numlevels_local; ++ipHi )
00566 {
00567 for( ipLo=0; ipLo<ipHi; ++ipLo )
00568 {
00569 if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
00570 continue;
00571
00572
00573 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc =
00574 StatesElemNEW[nelem][nelem-ipISO][ipLo].Pop - StatesElemNEW[nelem][nelem-ipISO][ipHi].Pop*
00575 StatesElemNEW[nelem][nelem-ipISO][ipLo].g/StatesElemNEW[nelem][nelem-ipISO][ipHi].g;
00576
00577
00578 if( N_(ipHi) > iso.n_HighestResolved_local[ipISO][nelem] )
00579 Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc = StatesElemNEW[nelem][nelem-ipISO][ipLo].Pop;
00580 }
00581 }
00582
00583 return;
00584 }