00001
00002
00003
00004 #include "cddefines.h"
00005 #include "yield.h"
00006 #include "prt.h"
00007 #include "continuum.h"
00008 #include "iso.h"
00009 #include "dynamics.h"
00010 #include "grainvar.h"
00011 #include "hmi.h"
00012 #include "mole.h"
00013 #include "thermal.h"
00014 #include "newton_step.h"
00015 #include "thirdparty.h"
00016 #include "conv.h"
00017 #include "secondaries.h"
00018 #include "phycon.h"
00019 #include "atmdat.h"
00020 #include "heavy.h"
00021 #include "elementnames.h"
00022 #include "dense.h"
00023 #include "radius.h"
00024 #include "ionbal.h"
00025 #include "taulines.h"
00026 #include "trace.h"
00027
00028 #define SOLVE_TWO 0
00029
00030
00031 STATIC bool lgTrivialSolution( long nelem, double abund_total );
00032
00033 STATIC void find_solution( long nelem, long ion_range, valarray<double> &xmat, valarray<double> &source);
00034
00035 STATIC void fill_array( long int nelem, long ion_range, valarray<double> &xmat, valarray<double> &source, valarray<double> &auger, double *abund_total );
00036
00037 #if SOLVE_TWO
00038 STATIC void combine_arrays( valarray<double> &xmat, const valarray<double> &xmat1, const valarray<double> &xmat2, long ion_range1, long ion_range2 );
00039 #endif
00040
00041 STATIC double get_total_abundance_ions( long int nelem );
00042
00043 STATIC void HomogeneousSource( long nelem, long ion_low, long ion_range, valarray<double> &xmat, valarray<double> &source, double abund_total );
00044
00045 STATIC void store_new_densities( long nelem, long ion_range, long ion_low, double *source, double abund_total, bool *lgNegPop ) ;
00046
00047 STATIC void PrintRates( long nelem, bool lgNegPop, double abund_total, valarray<double> &auger, bool lgPrintIt );
00048
00049 void solveions(double *ion, double *rec, double *snk, double *src,
00050 long int nlev, long int nmax);
00051
00052
00053 # undef MAT
00054 # define MAT(M_,I_,J_) ((M_)[(I_)*(ion_range)+(J_)])
00055
00056 # undef MAT1
00057 # define MAT1(M_,I_,J_) ((M_)[(I_)*(ion_range1)+(J_)])
00058
00059 # undef MAT2
00060 # define MAT2(M_,I_,J_) ((M_)[(I_)*(ion_range2)+(J_)])
00061
00062 void ion_solver( long int nelem, bool lgPrintIt)
00063 {
00064 double abund_total = 0.0;
00065 long ion_low = dense.IonLow[nelem];
00066 bool lgNegPop = false;
00067 double error = 0.0;
00068
00069 DEBUG_ENTRY( "ion_solver()" );
00070
00071 iso_charge_transfer_update(nelem);
00072
00073 long ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
00074 valarray<double> xmat(ion_range*ion_range);
00075 valarray<double> source(ion_range);
00076 valarray<double> auger(LIMELM+1);
00077
00078
00079
00080
00081
00082 conv.incrementCounter(ION_SOLVES);
00083 for (long it=0; it<4; ++it)
00084 {
00085 conv.incrementCounter(ISO_LOOPS);
00086 for( long ipISO=ipH_LIKE; ipISO<MIN2(NISO,nelem+1); ipISO++ )
00087 {
00088 if( (dense.IonHigh[nelem] >= nelem - ipISO) &&
00089 (dense.IonLow[nelem] <= nelem - ipISO))
00090 {
00091 iso_set_ion_rates(ipISO, nelem);
00092 }
00093 }
00094
00095
00096
00097
00098 abund_total = get_total_abundance_ions( nelem );
00099
00100 bool lgTrivial = lgTrivialSolution( nelem, abund_total );
00101 if( !lgTrivial )
00102 {
00103
00104 fill_array( nelem, ion_range, xmat, source, auger, &abund_total );
00105
00106
00107 HomogeneousSource( nelem, ion_low, ion_range, xmat, source, abund_total );
00108
00109
00110 find_solution( nelem, ion_range, xmat, source);
00111
00112
00113 store_new_densities( nelem, ion_range, ion_low, &source[0], abund_total, &lgNegPop );
00114
00115 ASSERT(ion_range >= dense.IonHigh[nelem]-dense.IonLow[nelem]+1);
00116 if (ion_range != dense.IonHigh[nelem]-dense.IonLow[nelem]+1)
00117 {
00118 ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
00119 xmat.resize(ion_range*ion_range);
00120 source.resize(ion_range);
00121 }
00122
00123 }
00124
00125 error = 0.0;
00126
00127 for( long ipISO=ipH_LIKE; ipISO<MIN2(NISO,nelem+1); ipISO++ )
00128 {
00129 double thiserror;
00130 iso_solve( ipISO, nelem, thiserror );
00131 if (thiserror > error)
00132 error = thiserror;
00133 }
00134 if ( error < 1e-4 )
00135 break;
00136 }
00137
00138 iso_satellite_update(nelem);
00139
00140 for( long ipISO=ipH_LIKE; ipISO<MIN2(NISO,nelem+1); ipISO++ )
00141 {
00142
00143 iso_departure_coefficients( ipISO, nelem );
00144 }
00145
00146 if( prt.lgPrtArry[nelem] || lgPrintIt )
00147 PrintRates( nelem, lgNegPop, abund_total, auger, lgPrintIt );
00148
00149 return;
00150 }
00151
00152 #if SOLVE_TWO
00153
00154 void ion_solver( long int nelem1, long int nelem2, bool lgPrintIt)
00155 {
00156 bool lgHomogeneous = true;
00157 bool lgNegPop1 = false;
00158 bool lgNegPop2 = false;
00159
00160 DEBUG_ENTRY( "ion_solver()" );
00161
00162 long ion_low1 = dense.IonLow[nelem1];
00163 long ion_low2 = dense.IonLow[nelem2];
00164 long ion_range1 = dense.IonHigh[nelem1]-dense.IonLow[nelem1]+1;
00165 long ion_range2 = dense.IonHigh[nelem2]-dense.IonLow[nelem2]+1;
00166 long ion_range = ion_range1 + ion_range2;
00167 double abund_total1 = get_total_abundance_ions( nelem1 );
00168 double abund_total2 = get_total_abundance_ions( nelem2 );
00169 valarray<double> xmat(ion_range*ion_range);
00170 valarray<double> xmat1(ion_range1*ion_range1);
00171 valarray<double> xmat2(ion_range2*ion_range2);
00172 valarray<double> source(ion_range);
00173 valarray<double> auger1(LIMELM+1);
00174 valarray<double> auger2(LIMELM+1);
00175
00176 #define ENABLE_SIMULTANEOUS_SOLUTION 0
00177
00178 if( lgTrivialSolution( nelem1, abund_total1 ) )
00179 {
00180
00181 return;
00182 }
00183 else if( lgTrivialSolution( nelem2, abund_total2 ) || !ENABLE_SIMULTANEOUS_SOLUTION )
00184 {
00185
00186 ion_solver( nelem1, lgPrintIt );
00187 return;
00188 }
00189
00190
00191 if( dense.IonLow[nelem1]>0 || dense.IonLow[nelem2]>0 || nelem1!=ipOXYGEN || nelem2!=ipHYDROGEN )
00192 {
00193 fprintf( ioQQQ, "This routine is currently intended to do only O and H when both have significant neutral fractions.\n" );
00194 fprintf( ioQQQ, "It should be generalized further for other cases. Exiting. Sorry.\n" );
00195 cdEXIT( EXIT_FAILURE );
00196 }
00197
00198 ASSERT( nelem1 != nelem2 );
00199
00200
00201 fill_array( nelem1, ion_range1, xmat1, source, auger1, &abund_total1 );
00202
00203 fill_array( nelem2, ion_range2, xmat2, source, auger2, &abund_total2 );
00204
00205 combine_arrays( xmat, xmat1, xmat2, ion_range1, ion_range2 );
00206
00207
00208 for( long i=0; i< ion_range; ++i )
00209 source[i] = 0.;
00210
00211
00212 for( long i=0; i< ion_range1; ++i )
00213 MAT( xmat, i, 0 ) = 1.;
00214
00215 for( long i=ion_range1; i< ion_range; ++i )
00216 MAT( xmat, i, ion_range1 ) = 1.;
00217
00218 source[0] = dense.xIonDense[nelem1][0] + dense.xIonDense[nelem1][1];
00219 source[ion_range1] = dense.xIonDense[nelem2][0] + dense.xIonDense[nelem2][1];
00220
00221
00222 lgHomogeneous = true;
00223
00224
00225 find_solution( nelem1, ion_range, xmat, source);
00226
00227
00228 store_new_densities( nelem1, ion_range1, ion_low1, &source[0], abund_total1, &lgNegPop1 );
00229 store_new_densities( nelem2, ion_range2, ion_low2, &source[ion_range1], abund_total2, &lgNegPop2 );
00230 ASSERT( lgNegPop2 == false );
00231
00232 if( prt.lgPrtArry[nelem1] || lgPrintIt )
00233 PrintRates( nelem1, lgNegPop1, abund_total1, auger1, lgPrintIt );
00234
00235 return;
00236 }
00237 #endif
00238
00239 STATIC bool lgTrivialSolution( long nelem, double abund_total )
00240 {
00241 double renorm;
00242
00243 if( dense.IonHigh[nelem] == dense.IonLow[nelem] )
00244 {
00245
00246 dense.xIonDense[nelem][dense.IonHigh[nelem]] = abund_total;
00247 for ( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00248 iso_renorm(nelem, ipISO, renorm);
00249 return true;
00250 }
00251 else if( dense.lgSetIoniz[nelem] )
00252 {
00253
00254 for( long ion=0; ion<nelem+2; ++ion )
00255 dense.xIonDense[nelem][ion] = dense.SetIoniz[nelem][ion]*abund_total;
00256 for ( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00257 iso_renorm(nelem, ipISO, renorm);
00258 return true;
00259 }
00260 else
00261 return false;
00262 }
00263
00264 STATIC void find_solution( long nelem, long ion_range, valarray<double> &xmat, valarray<double> &source)
00265 {
00266 int32 nerror;
00267 valarray<double> xmatsave(ion_range*ion_range);
00268 valarray<double> sourcesave(ion_range);
00269 valarray<int32> ipiv(ion_range);
00270
00271 DEBUG_ENTRY("find_solution()");
00272
00273
00274 for( long i=0; i< ion_range; ++i )
00275 {
00276 sourcesave[i] = source[i];
00277 for( long j=0; j< ion_range; ++j )
00278 {
00279 MAT( xmatsave, i, j ) = MAT( xmat, i, j );
00280 }
00281 }
00282
00283 if (1)
00284 {
00285 nerror = solve_system(xmat,source,ion_range,NULL);
00286
00287 if (nerror)
00288 {
00289
00290
00291 fprintf(ioQQQ,"Error for nelem %ld, active ion range %ld--%ld\n",
00292 nelem,dense.IonLow[nelem],dense.IonHigh[nelem]);
00293 fprintf(ioQQQ,"Initial ion abundances:");
00294 for( long j=0; j<nelem+2; ++j )
00295 fprintf(ioQQQ," %g",dense.xIonDense[nelem][j]);
00296 fprintf(ioQQQ,"\n");
00297 for( long j=0; j<ion_range; ++j )
00298 source[j] = dense.xIonDense[nelem][dense.IonLow[nelem]+j];
00299 }
00300 }
00301 else
00302 {
00303
00304 nerror = 0;
00305
00306 getrf_wrapper(ion_range, ion_range, &xmat[0], ion_range, &ipiv[0], &nerror);
00307 if( nerror != 0 )
00308 {
00309 fprintf( ioQQQ,
00310 " DISASTER ion_solver: dgetrf finds singular or ill-conditioned matrix nelem=%li %s ion_range=%li, nConv %li IonLow %li IonHi %li\n",
00311 nelem ,
00312 elementnames.chElementSym[nelem],
00313 ion_range,
00314 conv.nTotalIoniz ,
00315 dense.IonLow[nelem], dense.IonHigh[nelem]);
00316 fprintf( ioQQQ, " xmat follows\n");
00317 for( long i=0; i<ion_range; ++i )
00318 {
00319 for( long j=0;j<ion_range;j++ )
00320 {
00321 fprintf(ioQQQ,"%e\t",MAT(xmatsave,j,i));
00322 }
00323 fprintf(ioQQQ,"\n");
00324 }
00325 fprintf(ioQQQ,"source follows\n");
00326 for( long i=0; i<ion_range;i++ )
00327 {
00328 fprintf(ioQQQ,"%e\t",sourcesave[i]);
00329 }
00330 fprintf(ioQQQ,"\n");
00331 cdEXIT(EXIT_FAILURE);
00332 }
00333 getrs_wrapper('N', ion_range, 1, &xmat[0], ion_range, &ipiv[0], &source[0], ion_range, &nerror);
00334 if( nerror != 0 )
00335 {
00336 fprintf( ioQQQ, " DISASTER ion_solver: dgetrs finds singular or ill-conditioned matrix nelem=%li ionrange=%li\n",
00337 nelem , ion_range );
00338 cdEXIT(EXIT_FAILURE);
00339 }
00340 }
00341
00342 {
00343
00344 enum {DEBUG_LOC=false};
00345 if( DEBUG_LOC && nelem == ipHYDROGEN )
00346 {
00347 fprintf(ioQQQ,"debuggg ion_solver1 %ld\t%.2f\t%.4e\t%.4e\tIon\t%.3e\tRec\t%.3e\n",
00348 nelem,
00349 fnzone,
00350 phycon.te,
00351 dense.eden,
00352 ionbal.RateIonizTot(nelem,0) ,
00353 ionbal.RateRecomTot[nelem][0]);
00354 fprintf(ioQQQ," Msrc %.3e %.3e\n", mole.source[nelem][0], mole.source[nelem][1]);
00355 fprintf(ioQQQ," Msnk %.3e %.3e\n", mole.sink[nelem][0], mole.sink[nelem][1]);
00356 fprintf(ioQQQ," Poprat %.3e nomol %.3e\n",source[1]/source[0],
00357 ionbal.RateIonizTot(nelem,0)/ionbal.RateRecomTot[nelem][0]);
00358 }
00359 }
00360
00361 for( long i=0; i<ion_range; i++ )
00362 {
00363 ASSERT( !isnan( source[i] ) );
00364 ASSERT( source[i] < MAX_DENSITY );
00365 }
00366
00367 return;
00368 }
00369
00370 bool lgOH_ChargeTransferDominant( void )
00371 {
00372 #define THRESHOLD 0.25
00373
00374 bool lgDominant = false;
00375 double denominator;
00376
00377 if( (denominator = ionbal.RateIonizTot(ipHYDROGEN,0) + atmdat.HCharExcIonTotal + mole.sink[ipHYDROGEN][0] ) > 0. )
00378 {
00379 lgDominant = lgDominant ||
00380 ( atmdat.HCharExcRecTo[ipOXYGEN][0]*dense.xIonDense[ipOXYGEN][1]/denominator > THRESHOLD );
00381 }
00382
00383 if( (denominator = ionbal.RateRecomTot[ipHYDROGEN][0] + atmdat.HCharExcRecTotal ) > 0. )
00384 {
00385 lgDominant = lgDominant ||
00386 ( atmdat.HCharExcIonOf[ipOXYGEN][0]*dense.xIonDense[ipOXYGEN][0]/denominator > THRESHOLD );
00387 }
00388
00389 return lgDominant;
00390 }
00391
00392 #if SOLVE_TWO
00393 STATIC void combine_arrays( valarray<double> &xmat, const valarray<double> &xmat1, const valarray<double> &xmat2, long ion_range1, long ion_range2 )
00394 {
00395 DEBUG_ENTRY( "combine_arrays()" );
00396
00397 long ion_range = ion_range1 + ion_range2;
00398
00399 for( long i=0; i<ion_range1; i++ )
00400 for( long j=0; j<ion_range1; j++ )
00401 MAT( xmat, i, j ) = MAT1( xmat1, i, j );
00402
00403 for( long i=0; i<ion_range2; i++ )
00404 for( long j=0; j<ion_range2; j++ )
00405 MAT( xmat, i+ion_range1, j+ion_range1 ) = MAT2( xmat2, i, j );
00406
00407 #if 0
00408 bool lgNoDice = false;
00409 for( long i=0; i<ion_range1; i++ )
00410 {
00411 if( !fp_equal( MAT( xmat, i, 0), 1.0, 1 ) )
00412 {
00413 lgNoDice = true;
00414 break;
00415 }
00416 }
00417
00418 if( !lgNoDice )
00419 {
00420 for( long i=ion_range1; i<ion_range; i++ )
00421 MAT( xmat, i, 0 ) = 1.0;
00422 }
00423 #endif
00424
00425 return;
00426 }
00427 #endif
00428
00429 STATIC void store_new_densities( long nelem, long ion_range, long ion_low, double *source, double abund_total, bool *lgNegPop )
00430 {
00431 DEBUG_ENTRY( "store_new_densities()" );
00432
00433 ASSERT( nelem >= 0 );
00434 ASSERT( nelem < LIMELM );
00435 ASSERT( ion_range <= nelem + 2 );
00436 ASSERT( ion_low >= 0 );
00437 ASSERT( ion_low <= nelem + 1 );
00438
00439 #if 0
00440 # define RJRW 0
00441 if( RJRW && 0 )
00442 {
00443
00444 double test;
00445 for(long i=0; i<ion_range; i++) {
00446 test = 0.;
00447 for(long j=0; j<ion_range; j++) {
00448 test = test+source[j]*MAT(xmatsave,j,i);
00449 }
00450 fprintf(ioQQQ,"%e\t",test);
00451 }
00452 fprintf(ioQQQ,"\n");
00453
00454 test = 0.;
00455 fprintf(ioQQQ," ion %li abundance %.3e\n",nelem,abund_total);
00456 for( long ion=dense.IonLow[nelem]; ion < dense.IonHigh[nelem]; ion++ )
00457 {
00458 if( ionbal.RateRecomTot[nelem][ion] != 0 && source[ion-ion_low] != 0 )
00459 fprintf(ioQQQ," %li %.3e %.3e : %.3e\n",ion,source[ion-ion_low],
00460 source[ion-ion_low+1]/source[ion-ion_low],
00461 ionbal.RateIonizTot(nelem,ion)/ionbal.RateRecomTot[nelem][ion]);
00462 else
00463 fprintf(ioQQQ," %li %.3e [One ratio infinity]\n",ion,source[ion-ion_low]);
00464 test += source[ion-ion_low];
00465 }
00466 }
00467 #endif
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491 *lgNegPop = false;
00492 for( long i=0; i < ion_range; i++ )
00493 {
00494 long ion = i+ion_low;
00495
00496 if( source[i] < 0. )
00497 {
00498
00499
00500
00501
00502 if( source[i]<-2e-9 )
00503 {
00504 fprintf(ioQQQ,
00505 " PROBLEM negative ion population in ion_solver, nelem=%li, %s ion=%li val=%.3e Search?%c zone=%li iteration=%li\n",
00506 nelem ,
00507 elementnames.chElementSym[nelem],
00508 ion ,
00509 source[i] ,
00510 TorF(conv.lgSearch) ,
00511 nzone ,
00512 iteration );
00513 *lgNegPop = true;
00514 fixit();
00515
00516 }
00517
00518 fixit();
00519 #if 1
00520 source[i] = 0.;
00521
00522
00523 if( ion > nelem-NISO && ion < nelem + 1 )
00524 {
00525 long int ipISO = nelem - ion;
00526 ASSERT( ipISO>=0 && ipISO<NISO );
00527 for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
00528 iso_sp[ipISO][nelem].st[level].Pop() = 0.;
00529 }
00530 #endif
00531 }
00532 }
00533
00534 double renormnew = 1.;
00535 {
00536 double dennew = 0.;
00537
00538
00539 for( long i=0;i < ion_range; i++ )
00540 {
00541 dennew += source[i];
00542 }
00543
00544 if(dennew > 0.)
00545 {
00546 renormnew = abund_total / dennew;
00551 }
00552 else
00553 {
00554 renormnew = 1.;
00555 }
00556
00557 dennew = 0.;
00558 for( long i=0;i < ion_range; i++ )
00559 {
00560 source[i] *= renormnew;
00561 dennew += source[i];
00562 }
00563 }
00564
00565
00566 if( renormnew < 0)
00567 {
00568 fprintf(ioQQQ,"PROBLEM impossible value of renorm \n");
00569 }
00570 ASSERT( renormnew>=0 );
00571
00572 for( long i=0; i < ion_range; i++ )
00573 {
00574 long ion = i+ion_low;
00575 dense.xIonDense[nelem][ion] = source[i];
00576 if( dense.xIonDense[nelem][ion] >= MAX_DENSITY )
00577 {
00578 fprintf( ioQQQ, "PROBLEM DISASTER Huge density in ion_solver, nelem %ld ion %ld source %e renormnew %e\n",
00579 nelem, ion, source[i], renormnew );
00580 }
00581 ASSERT( dense.xIonDense[nelem][ion] < MAX_DENSITY );
00582 }
00583
00584 fixit();
00585
00586
00587 while( dense.IonHigh[nelem] > dense.IonLow[nelem] &&
00588 dense.xIonDense[nelem][dense.IonHigh[nelem]] < 1e-25*abund_total )
00589 {
00590 ASSERT( dense.xIonDense[nelem][dense.IonHigh[nelem]] >= 0. );
00591
00592 dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.;
00593 thermal.heating[nelem][dense.IonHigh[nelem]-1] = 0.;
00594
00595 --dense.IonHigh[nelem];
00596 }
00597
00598 for ( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
00599 {
00600 double renorm;
00601 iso_renorm(nelem, ipISO, renorm);
00602 }
00603
00604
00605 ASSERT( (dense.IonLow[nelem] <= dense.IonHigh[nelem]) ||
00606
00607 (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) ||
00608
00609 ( dense.IonLow[nelem]==nelem+1 && dense.IonHigh[nelem]==nelem+1 ) );
00610
00611 fixit();
00612
00613 return;
00614 }
00615
00616 STATIC double get_total_abundance_ions( long int nelem )
00617 {
00618 DEBUG_ENTRY( "get_total_abundance_ions()" );
00619
00620 ASSERT( nelem >= 0 );
00621 ASSERT( nelem < LIMELM );
00622
00623 ionbal.elecsnk[nelem] = 0.;
00624 ionbal.elecsrc[nelem] = 0.;
00625
00626 double abund_total = 0.;
00627 for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
00628 {
00629 abund_total += dense.xIonDense[nelem][ion];
00630 }
00631
00632 realnum tot1 = dense.gas_phase[nelem];
00633 realnum tot2 = (realnum)(dense.xMolecules(nelem)+abund_total);
00634 if (0)
00635 {
00636 if (fabs(tot2-tot1) > conv.GasPhaseAbundErrorAllowed*tot1)
00637 fprintf(ioQQQ,"%ld %13.8g %13.8g %13.8g %13.8g\n",nelem,tot1,tot2,abund_total,tot2/tot1 - 1.0);
00638 }
00639
00640 ASSERT( fp_equal_tol(tot1, tot2, realnum(conv.GasPhaseAbundErrorAllowed*tot1 + 100.f*FLT_MIN)) );
00641
00642 ASSERT( abund_total < MAX_DENSITY );
00643
00644 return abund_total;
00645 }
00646
00647 STATIC void fill_array( long int nelem, long ion_range, valarray<double> &xmat, valarray<double> &source, valarray<double> &auger, double *abund_total )
00648 {
00649 long int limit,
00650 IonProduced;
00651 double rateone;
00652 long ion_low;
00653
00654 valarray<double> sink(ion_range);
00655 valarray<int32> ipiv(ion_range);
00656
00657 DEBUG_ENTRY( "fill_array()" );
00658
00659
00660 ASSERT( nelem >= 0);
00661 ASSERT( dense.IonLow[nelem] >= 0 );
00662 ASSERT( dense.IonHigh[nelem] >= 0 );
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678 limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);
00679
00680
00681 ASSERT( ion_range <= nelem+2 );
00682
00683 ion_low = dense.IonLow[nelem];
00684 for( long i=0; i<ion_range;i++ )
00685 {
00686 source[i] = 0.;
00687 }
00688
00689
00690
00691
00692
00693 for( long ion_from=0; ion_from <= limit; ion_from++ )
00694 {
00695 for( long ion_to=0; ion_to < nelem+2; ion_to++ )
00696 {
00697 ionbal.RateIoniz[nelem][ion_from][ion_to] = 0.;
00698 }
00699 }
00700
00701
00702
00703
00704 for( long i=0; i< LIMELM+1; ++i )
00705 {
00706 auger[i] = 0.;
00707 }
00708
00709
00710 for( long i=0; i< ion_range; ++i )
00711 {
00712 for( long j=0; j< ion_range; ++j )
00713 {
00714 MAT( xmat, i, j ) = 0.;
00715 }
00716 }
00717
00718 {
00719
00720
00721 enum {DEBUG_LOC=false};
00722 if( DEBUG_LOC && nelem==ipCARBON )
00723 {
00724 broken();
00725 dense.IonLow[nelem] = 0;
00726 dense.IonHigh[nelem] = 3;
00727 *abund_total = 1.;
00728
00729 for( long ion=dense.IonLow[nelem]; ion <= limit; ion++ )
00730 {
00731 double fac=1;
00732 if(ion)
00733 fac = 1e-10;
00734 ionbal.RateRecomTot[nelem][ion] = 100.;
00735 for( long ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
00736 {
00737
00738 ionbal.PhotoRate_Shell[nelem][ion][ns][0] = fac;
00739 }
00740 }
00741 }
00742 }
00743
00744 bool lgGrainsOn = gv.lgDustOn() && ionbal.lgGrainIonRecom && gv.lgGrainPhysicsOn;
00745
00746
00747
00748
00749
00750 for( long ion_to=dense.IonLow[nelem]; ion_to <= dense.IonHigh[nelem]; ion_to++ )
00751 {
00752 for( long ion_from=dense.IonLow[nelem]; ion_from <= dense.IonHigh[nelem]; ++ion_from )
00753 {
00754
00755 if( ion_to != ion_from )
00756 {
00757
00758 rateone = mole.xMoleChTrRate[nelem][ion_from][ion_to] * atmdat.lgCTOn;
00759 MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
00760 MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
00761
00762 rateone = gv.GrainChTrRate[nelem][ion_from][ion_to]*lgGrainsOn;
00763 MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
00764 MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
00765 }
00766 }
00767 }
00768
00769 for( long ion=dense.IonLow[nelem]; ion <= limit; ion++ )
00770 {
00771
00772 ionbal.RateIoniz[nelem][ion][ion+1] +=
00773 ionbal.CollIonRate_Ground[nelem][ion][0] +
00774 secondaries.csupra[nelem][ion] +
00775
00776 ionbal.UTA_ionize_rate[nelem][ion];
00777
00778
00779 for( long ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
00780 {
00781
00782
00783
00784
00785
00786
00787
00788 if( ion+1-ion_low < ion_range )
00789 {
00790
00791
00792
00793
00794
00795
00796
00797 for( long nej=1; nej <= t_yield::Inst().nelec_eject(nelem,ion,ns); nej++ )
00798 {
00799
00800
00801
00802 IonProduced = MIN2(ion+nej,dense.IonHigh[nelem]);
00803 rateone = ionbal.PhotoRate_Shell[nelem][ion][ns][0]*
00804 t_yield::Inst().elec_eject_frac(nelem,ion,ns,nej-1);
00805
00806
00807 ionbal.RateIoniz[nelem][ion][IonProduced] += rateone;
00808
00809
00810
00811 if( nej>1 )
00812 auger[IonProduced-1] += rateone;
00813 }
00814 }
00815 }
00816 }
00817
00818 for( long ion=dense.IonLow[nelem]; ion < dense.IonHigh[nelem]; ion++ )
00819 {
00820
00821 double ction;
00822 long ipISO=nelem-ion;
00823
00824 if ( nelem==ipHYDROGEN )
00825 ction = atmdat.HCharExcIonTotal;
00826 else if( nelem==ipHELIUM && ipISO==ipHE_LIKE )
00827 ction = atmdat.HeCharExcIonTotal;
00828 else
00829 ction = atmdat.HeCharExcIonOf[nelem][ion]*dense.xIonDense[ipHELIUM][1]+
00830 atmdat.HCharExcIonOf[nelem][ion]*dense.xIonDense[ipHYDROGEN][1];
00831
00832
00833 MAT( xmat, ion-ion_low, ion-ion_low ) -= ction;
00834 MAT( xmat, ion-ion_low, ion+1-ion_low ) += ction;
00835 }
00836
00837 for( long ion=dense.IonLow[nelem]; ion < dense.IonHigh[nelem]; ion++ )
00838 {
00839
00840 double ctrec;
00841 long ipISO=nelem-ion;
00842
00843 if ( nelem==ipHYDROGEN )
00844 ctrec = atmdat.HCharExcRecTotal;
00845 else if( nelem==ipHELIUM && ipISO==ipHE_LIKE )
00846 ctrec = atmdat.HeCharExcRecTotal;
00847 else
00848 ctrec =
00849 atmdat.HeCharExcRecTo[nelem][ion]*dense.xIonDense[ipHELIUM][0] +
00850 atmdat.HCharExcRecTo[nelem][ion]*dense.xIonDense[ipHYDROGEN][0];
00851 MAT( xmat, ion+1-ion_low, ion+1-ion_low ) -= ctrec;
00852 MAT( xmat, ion+1-ion_low, ion-ion_low ) += ctrec;
00853
00854 }
00855
00856
00857 for( long ion_from=dense.IonLow[nelem]; ion_from < dense.IonHigh[nelem]; ion_from++ )
00858 {
00859 for( long ion_to=ion_from+1; ion_to <= dense.IonHigh[nelem]; ion_to++ )
00860 {
00861 ionbal.elecsrc[nelem] += ionbal.RateIoniz[nelem][ion_from][ion_to]*dense.xIonDense[nelem][ion_from]*
00862 (ion_to-ion_from);
00863
00864 MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= ionbal.RateIoniz[nelem][ion_from][ion_to];
00865 MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += ionbal.RateIoniz[nelem][ion_from][ion_to];
00866 }
00867 }
00868
00869 for( long ion=dense.IonLow[nelem]; ion<dense.IonHigh[nelem]; ion++ )
00870 {
00871
00872 MAT( xmat, ion+1-ion_low, ion+1-ion_low ) -= ionbal.RateRecomTot[nelem][ion];
00873 MAT( xmat, ion+1-ion_low, ion-ion_low ) += ionbal.RateRecomTot[nelem][ion];
00874 ionbal.elecsnk[nelem] += ionbal.RateRecomTot[nelem][ion]*dense.xIonDense[nelem][ion+1];
00875 }
00876
00877 for( long i=0; i<ion_range;i++ )
00878 {
00879 long ion = i+ion_low;
00880
00881
00882
00883 source[i] -= mole.source[nelem][ion];
00884 MAT( xmat, i, i ) -= mole.sink[nelem][ion];
00885 }
00886
00887 return;
00888 }
00889
00890 STATIC void HomogeneousSource( long nelem, long ion_low, long ion_range, valarray<double> &xmat, valarray<double> &source, double abund_total )
00891 {
00892 bool lgHomogeneous = true;
00893 double totsrc,
00894 totscl,
00895 maxdiag;
00896
00897 DEBUG_ENTRY( "lgHomogeneousSource()" );
00898
00899
00900
00901 totsrc = totscl = maxdiag = 0.;
00902 for( long i=0; i<ion_range;i++ )
00903 {
00904 long ion = i + ion_low;
00905
00906 totsrc -= source[i];
00907 fixit();
00908
00909
00910
00911 double diag = -(MAT( xmat, i, i )+mole.sink[nelem][ion]);
00912 if (diag > maxdiag)
00913 maxdiag = diag;
00914 totscl += diag*dense.xIonDense[nelem][ion];
00915 }
00916
00917
00918
00919 const double CONDITION_SCALE = 1e8;
00920 double conserve_scale = maxdiag/CONDITION_SCALE;
00921 ASSERT( conserve_scale > 0.0 );
00922
00923
00924 if( totsrc > 1e-10*totscl )
00925 lgHomogeneous = false;
00926
00927 fixit();
00928
00929
00930 if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection && !dynamics.lgEquilibrium
00931 && dynamics.Rate != 0. )
00932 {
00933 for( long i=0; i<ion_range;i++ )
00934 {
00935 long ion = i+ion_low;
00936 MAT( xmat, i, i ) -= dynamics.Rate;
00937 source[i] -= dynamics.Source[nelem][ion];
00938
00939
00940 }
00941 lgHomogeneous = false;
00942 }
00943
00944
00945
00946
00947
00948
00949
00950
00951 if( !lgHomogeneous && ion_range==2 )
00952 {
00953
00954 double a = MAT( xmat, 0, 0 ),
00955 b = MAT( xmat, 1, 0 ) ,
00956 c = MAT( xmat, 0, 1 ) ,
00957 d = MAT( xmat, 1, 1 );
00958 b = SDIV(b);
00959 d = SDIV(d);
00960 double ratio1 = a/b , ratio2 = c/d , fratio1=fabs(a/b),fratio2=fabs(c/d);
00961 if( fabs(ratio1-ratio2) <= DBL_EPSILON*1e4*MAX2(fratio1,fratio2) )
00962 {
00963 lgHomogeneous = true;
00964 }
00965 }
00966
00967 #if 1
00968 if(
00969 fabs(dense.xMolecules(nelem)) <DBL_EPSILON*100.*dense.gas_phase[nelem] )
00970 {
00971 lgHomogeneous = true;
00972 }
00973 #endif
00974
00975
00976
00977
00978 if( 1 || lgHomogeneous )
00979 {
00980 double rate_ioniz=1., rate_recomb ;
00981
00982 long jmax = 0;
00983 for( long i=0; i<ion_range-1;i++)
00984 {
00985 long ion = i+ion_low;
00986 rate_ioniz *= ionbal.RateIonizTot(nelem,ion);
00987 rate_recomb = ionbal.RateRecomTot[nelem][ion];
00988
00989 if( rate_ioniz == 0)
00990 break;
00991
00992 if( rate_recomb <= 0.)
00993 break;
00994
00995 rate_ioniz /= rate_recomb;
00996 if( rate_ioniz > 1.)
00997 {
00998
00999 jmax = i;
01000 rate_ioniz = 1.;
01001 }
01002 }
01003
01004
01005 for( long i=0; i<ion_range;i++ )
01006 {
01007 MAT(xmat,i,jmax) -= conserve_scale;
01008 }
01009 source[jmax] -= abund_total*conserve_scale;
01010 }
01011
01012 #if 0
01013 if( false && nelem == ipHYDROGEN && dynamics.lgAdvection&& iteration>1 )
01014 {
01015 fprintf(ioQQQ,"DEBUGG Rate %.2f %.3e \n",fnzone,dynamics.Rate);
01016 fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateIonizTot(nelem,0), ionbal.RateIonizTot(nelem,1) );
01017 fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateRecomTot[nelem][0], ionbal.RateRecomTot[nelem][1]);
01018 fprintf(ioQQQ," %.3e %.3e %.3e\n\n", dynamics.Source[nelem][0], dynamics.Source[nelem][1], dynamics.Source[nelem][2]);
01019 }
01020
01021 {
01022
01023 enum {DEBUG_LOC=false};
01024 if( DEBUG_LOC && nzone==1 && lgPrintIt )
01025 {
01026 fprintf( ioQQQ,
01027 " DEBUG ion_solver: nelem=%li ion_range=%li, limit=%li, nConv %li xmat follows\n",
01028 nelem , ion_range,limit , conv.nTotalIoniz );
01029 if( lgHomogeneous )
01030 fprintf(ioQQQ , "Homogeneous \n");
01031 for( long i=0; i<ion_range; ++i )
01032 {
01033 for( long j=0;j<ion_range;j++ )
01034 {
01035 fprintf(ioQQQ,"%e\t",MAT(xmat,i,j));
01036 }
01037 fprintf(ioQQQ,"\n");
01038 }
01039 fprintf(ioQQQ,"source follows\n");
01040 for( long i=0; i<ion_range;i++ )
01041 {
01042 fprintf(ioQQQ,"%e\t",source[i]);
01043 }
01044 fprintf(ioQQQ,"\n");
01045 fprintf(ioQQQ,"totsrc/totscl %g %g\n",totsrc,totscl);
01046 }
01047 }
01048 #endif
01049
01050 }
01051
01052 STATIC void PrintRates( long nelem, bool lgNegPop, double abund_total, valarray<double> &auger, bool lgPrintIt )
01053 {
01054 DEBUG_ENTRY( "PrintRates()" );
01055
01056 long ion;
01057
01058
01059 if( lgNegPop )
01060 {
01061 fprintf( ioQQQ, " PROBLEM Negative population found for abundance of ionization stage of element %4.4s, ZONE=%4ld\n",
01062 elementnames.chElementNameShort[nelem], nzone );
01063
01064 fprintf( ioQQQ, " Populations were" );
01065 for( ion=1; ion <= dense.IonHigh[nelem]+1; ion++ )
01066 {
01067 fprintf( ioQQQ, "%9.1e", dense.xIonDense[nelem][ion-1] );
01068 }
01069 fprintf( ioQQQ, "\n" );
01070
01071 fprintf( ioQQQ, " destroy vector =" );
01072 for( ion=1; ion <= dense.IonHigh[nelem]; ion++ )
01073 {
01074 fprintf( ioQQQ, "%9.1e", ionbal.RateIonizTot(nelem,ion-1) );
01075 }
01076 fprintf( ioQQQ, "\n" );
01077
01078 fprintf( ioQQQ, " CTHeavy vector =" );
01079 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01080 {
01081 fprintf( ioQQQ, "%9.1e", atmdat.HeCharExcIonOf[nelem][ion] );
01082 }
01083 fprintf( ioQQQ, "\n" );
01084
01085 fprintf( ioQQQ, " HCharExcIonOf vtr=" );
01086 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01087 {
01088 fprintf( ioQQQ, "%9.1e", atmdat.HCharExcIonOf[nelem][ion] );
01089 }
01090 fprintf( ioQQQ, "\n" );
01091
01092 fprintf( ioQQQ, " CollidRate vtr=" );
01093 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01094 {
01095 fprintf( ioQQQ, "%9.1e", ionbal.CollIonRate_Ground[nelem][ion][0] );
01096 }
01097 fprintf( ioQQQ, "\n" );
01098
01099
01100 fprintf( ioQQQ, " photo rates per subshell, ion\n" );
01101 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01102 {
01103 fprintf( ioQQQ, "%3ld", ion );
01104 for( long ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
01105 {
01106 fprintf( ioQQQ, "%9.1e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] );
01107 }
01108 fprintf( ioQQQ, "\n" );
01109 }
01110
01111
01112 fprintf( ioQQQ, " create vector =" );
01113 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01114 {
01115 fprintf( ioQQQ, "%9.1e", ionbal.RateRecomTot[nelem][ion] );
01116 }
01117 fprintf( ioQQQ, "\n" );
01118 }
01119
01120
01121
01122 if( lgPrintIt || prt.lgPrtArry[nelem] || lgNegPop )
01123 {
01124 const char* format = " %9.2e";
01125
01126 fprintf( ioQQQ,
01127 "\n %s ion_solver ion/rec rt [s-1] %s nz%.2f Te%.4e ne%.4e Tot abun:%.3e ion abun%.2e mole%.2e\n",
01128 elementnames.chElementSym[nelem],
01129 elementnames.chElementName[nelem],
01130 fnzone,
01131 phycon.te ,
01132 dense.eden,
01133 dense.gas_phase[nelem],
01134 abund_total ,
01135 dense.xMolecules(nelem) );
01136
01137 fprintf( ioQQQ, " %s Ioniz total " ,elementnames.chElementSym[nelem]);
01138 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01139 {
01140 fprintf( ioQQQ, format, ionbal.RateIonizTot(nelem,ion) );
01141 }
01142 fprintf( ioQQQ, "\n" );
01143
01144
01145 fprintf( ioQQQ, " %s Recom total " ,elementnames.chElementSym[nelem]);
01146 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01147 {
01148 fprintf( ioQQQ, format, ionbal.RateRecomTot[nelem][ion] );
01149 }
01150 fprintf( ioQQQ, "\n" );
01151
01152
01153 fprintf( ioQQQ, " %s Coll ioniz " ,elementnames.chElementSym[nelem] );
01154 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01155 {
01156 double ColIoniz = ionbal.CollIonRate_Ground[nelem][ion][0];
01157 if( ion > nelem - NISO )
01158 {
01159 long ipISO = nelem-ion;
01160 ASSERT( ipISO >=0 && ipISO < NISO );
01161 if( dense.xIonDense[nelem][nelem-ipISO] > 0. )
01162 {
01163 ColIoniz *= iso_sp[ipISO][nelem].st[0].Pop();
01164 for( long ipLevel=1; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
01165 {
01166 ColIoniz += iso_sp[ipISO][nelem].fb[ipLevel].ColIoniz * dense.EdenHCorr * iso_sp[ipISO][nelem].st[ipLevel].Pop();
01167 }
01168 ColIoniz /= dense.xIonDense[nelem][nelem-ipISO];
01169 }
01170 }
01171 fprintf( ioQQQ, format, ColIoniz );
01172 }
01173 fprintf( ioQQQ, "\n" );
01174
01175
01176 fprintf( ioQQQ, " %s UTA ioniz " ,elementnames.chElementSym[nelem] );
01177 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01178 {
01179 fprintf( ioQQQ, format, ionbal.UTA_ionize_rate[nelem][ion] );
01180 }
01181 fprintf( ioQQQ, "\n" );
01182
01183
01184 fprintf( ioQQQ, " %s Photoion snk" ,elementnames.chElementSym[nelem] );
01185 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01186 {
01187 double PhotIoniz = 0.;
01188 for( long ipShell = 0; ipShell < Heavy.nsShells[nelem][ion]; ipShell++ )
01189 PhotIoniz += ionbal.PhotoRate_Shell[nelem][ion][ipShell][0];
01190
01191
01192 if( ion > nelem - NISO )
01193 {
01194 long ipISO = nelem-ion;
01195 ASSERT( ipISO >=0 && ipISO < NISO );
01196 if( dense.xIonDense[nelem][nelem-ipISO]>0 )
01197 {
01198 PhotIoniz *= iso_sp[ipISO][nelem].st[0].Pop();
01199 for( long ipLevel=1; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
01200 {
01201 PhotIoniz += iso_sp[ipISO][nelem].fb[ipLevel].gamnc * iso_sp[ipISO][nelem].st[ipLevel].Pop();
01202 }
01203 PhotIoniz /= dense.xIonDense[nelem][nelem-ipISO];
01204 }
01205 }
01206 fprintf( ioQQQ, format, PhotIoniz );
01207 }
01208 fprintf( ioQQQ, "\n" );
01209
01210
01211 fprintf( ioQQQ, " %s Photoion src" ,elementnames.chElementSym[nelem]);
01212 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01213 {
01214 double source = 0.;
01215 if( ion>0 )
01216 source = ionbal.RateIoniz[nelem][ion-1][ion];
01217
01218 fprintf( ioQQQ, format, source );
01219 }
01220 fprintf( ioQQQ, "\n" );
01221
01222
01223 fprintf( ioQQQ, " %s Auger src " ,elementnames.chElementSym[nelem]);
01224 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01225 {
01226 double source = 0.;
01227 if( ion>0 )
01228 source = auger[ion-1];
01229
01230 fprintf( ioQQQ, format, source );
01231 }
01232 fprintf( ioQQQ, "\n" );
01233
01234
01235 fprintf( ioQQQ, " %s Secon ioniz " ,elementnames.chElementSym[nelem]);
01236 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01237 {
01238 fprintf( ioQQQ, format,
01239 secondaries.csupra[nelem][ion] );
01240 }
01241 fprintf( ioQQQ, "\n" );
01242
01243
01244 fprintf( ioQQQ, " %s ion on grn " ,elementnames.chElementSym[nelem]);
01245 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01246 {
01247 fprintf( ioQQQ, format, gv.GrainChTrRate[nelem][ion][ion+1] );
01248 }
01249 fprintf( ioQQQ, "\n" );
01250
01251
01252 fprintf( ioQQQ, " %s ion xfr mol " ,elementnames.chElementSym[nelem]);
01253 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01254 {
01255 fprintf( ioQQQ, format, mole.xMoleChTrRate[nelem][ion][ion+1] );
01256 }
01257 fprintf( ioQQQ, "\n" );
01258
01259
01260 fprintf( ioQQQ, " %s chr trn ion " ,elementnames.chElementSym[nelem] );
01261 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01262 {
01263
01264 double sum;
01265
01266 if( nelem==ipHELIUM && ion==0 )
01267 {
01268 sum = atmdat.HeCharExcIonTotal;
01269 }
01270 else if( nelem==ipHYDROGEN && ion==0 )
01271 {
01272 sum = atmdat.HCharExcIonTotal;
01273 }
01274 else
01275 sum = atmdat.HeCharExcIonOf[nelem][ion] * dense.xIonDense[ipHELIUM][1]+
01276 atmdat.HCharExcIonOf[nelem][ion] * dense.xIonDense[ipHYDROGEN][1];
01277
01278 fprintf( ioQQQ, format, sum );
01279 }
01280 fprintf( ioQQQ, "\n" );
01281
01282
01283 fprintf( ioQQQ, " %s radiati rec " ,elementnames.chElementSym[nelem]);
01284 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01285 {
01286 fprintf( ioQQQ, format, dense.eden*ionbal.RR_rate_coef_used[nelem][ion] );
01287 }
01288 fprintf( ioQQQ, "\n" );
01289
01290
01291 fprintf( ioQQQ, " %s drBadnel rec" ,elementnames.chElementSym[nelem]);
01292 for( ion=0; ion < min(nelem-1,dense.IonHigh[nelem]); ion++ )
01293 {
01294 fprintf( ioQQQ, format, dense.eden*ionbal.DR_Badnell_rate_coef[nelem][ion] );
01295 }
01296 fprintf( ioQQQ, "\n" );
01297
01298
01299 fprintf( ioQQQ, " %s CotaRate rec" ,elementnames.chElementSym[nelem]);
01300 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01301 {
01302 fprintf( ioQQQ, format, dense.eden*ionbal.CotaRate[ion] );
01303 }
01304 fprintf( ioQQQ, "\n" );
01305
01306
01307
01308 fprintf( ioQQQ, " %s rec on grn " ,elementnames.chElementSym[nelem]);
01309 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01310 {
01311 fprintf( ioQQQ, format, gv.GrainChTrRate[nelem][ion+1][ion] );
01312 }
01313 fprintf( ioQQQ, "\n" );
01314
01315
01316 fprintf( ioQQQ, " %s rec xfr mol " ,elementnames.chElementSym[nelem]);
01317 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01318 {
01319 fprintf( ioQQQ, format, mole.xMoleChTrRate[nelem][ion+1][ion] );
01320 }
01321 fprintf( ioQQQ, "\n" );
01322
01323
01324 fprintf( ioQQQ, " %s chr trn rec " ,elementnames.chElementSym[nelem]);
01325 for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
01326 {
01327 double sum;
01328
01329 if( nelem==ipHELIUM && ion==0 )
01330 {
01331 sum = atmdat.HeCharExcRecTotal;
01332 }
01333 else if( nelem==ipHYDROGEN && ion==0 )
01334 {
01335 sum = atmdat.HCharExcRecTotal;
01336 }
01337 else
01338 sum = atmdat.HCharExcRecTo[nelem][ion] * dense.xIonDense[ipHYDROGEN][0] +
01339 atmdat.HeCharExcRecTo[nelem][ion] * dense.xIonDense[ipHELIUM][0];
01340
01341 fprintf( ioQQQ, format, sum );
01342 }
01343 fprintf( ioQQQ, "\n" );
01344
01345 {
01346
01347 fprintf(ioQQQ," %s molecule src",
01348 elementnames.chElementSym[nelem]);
01349 for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
01350 {
01351 fprintf(ioQQQ,format, mole.source[nelem][ion]/SDIV(dense.xIonDense[nelem][ion]) );
01352 }
01353 fprintf( ioQQQ, "\n" );
01354
01355
01356 fprintf(ioQQQ," %s molecule snk",
01357 elementnames.chElementSym[nelem]);
01358 for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
01359 {
01360 fprintf(ioQQQ,format, mole.sink[nelem][ion] );
01361 }
01362 fprintf( ioQQQ, "\n" );
01363
01364 if( dynamics.lgAdvection )
01365 {
01366
01367 fprintf(ioQQQ," %s dynamics src",
01368 elementnames.chElementSym[nelem]);
01369 for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
01370 {
01371 fprintf(ioQQQ,format, dynamics.Source[nelem][ion]/SDIV(dense.xIonDense[nelem][ion]) );
01372 }
01373 fprintf( ioQQQ, "\n" );
01374
01375
01376 fprintf(ioQQQ," %s dynamics snk",
01377 elementnames.chElementSym[nelem]);
01378 for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
01379 {
01380 fprintf(ioQQQ,format, dynamics.Rate );
01381 }
01382 fprintf( ioQQQ, "\n" );
01383 }
01384 }
01385
01386
01387 fprintf( ioQQQ, " %s Abun [cm-3] " ,elementnames.chElementSym[nelem] );
01388 for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
01389 {
01390 fprintf( ioQQQ, format, dense.xIonDense[nelem][ion] );
01391 }
01392 fprintf( ioQQQ, "\n" );
01393 }
01394
01395 if( lgNegPop )
01396 {
01397 ContNegative();
01398 ShowMe();
01399 cdEXIT(EXIT_FAILURE);
01400 }
01401
01402 return;
01403 }
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432 void solveions(double *ion, double *rec, double *snk, double *src,
01433 long int nlev, long int nmax)
01434 {
01435 double kap, bet;
01436 long int i;
01437
01438 DEBUG_ENTRY( "solveions()" );
01439
01440 if(nmax != -1)
01441 {
01442
01443 src[nmax] = 1.;
01444 for(i=nmax;i<nlev-1;i++)
01445 src[i+1] = src[i]*ion[i]/rec[i];
01446 for(i=nmax-1;i>=0;i--)
01447 src[i] = src[i+1]*rec[i]/ion[i];
01448 }
01449 else
01450 {
01451 i = 0;
01452 kap = snk[0];
01453 for(i=0;i<nlev-1;i++)
01454 {
01455 bet = ion[i]+kap;
01456 if(bet == 0.)
01457 {
01458 fprintf(ioQQQ,"Ionization solver error\n");
01459 cdEXIT(EXIT_FAILURE);
01460 }
01461 bet = 1./bet;
01462 src[i] *= bet;
01463 src[i+1] += ion[i]*src[i];
01464 snk[i] = bet*rec[i];
01465 kap = kap*snk[i]+snk[i+1];
01466 }
01467 bet = kap;
01468 if(bet == 0.)
01469 {
01470 fprintf(ioQQQ,"Ionization solver error\n");
01471 cdEXIT(EXIT_FAILURE);
01472 }
01473 src[i] /= bet;
01474
01475 for(i=nlev-2;i>=0;i--)
01476 {
01477 src[i] += snk[i]*src[i+1];
01478 }
01479 }
01480 }
01481
01482 void ion_wrapper( long nelem )
01483 {
01484
01485 DEBUG_ENTRY( "ion_wrapper()" );
01486
01487 ASSERT( nelem >= ipHYDROGEN );
01488 ASSERT( nelem < LIMELM );
01489
01490 switch( nelem )
01491 {
01492 case ipHYDROGEN:
01493 IonHydro();
01494 break;
01495 case ipHELIUM:
01496 IonHelium();
01497 break;
01498 default:
01499 IonNelem(false,nelem);
01500 break;
01501 }
01502
01503 if( trace.lgTrace && dense.lgElmtOn[nelem] && trace.lgHeavyBug )
01504 {
01505 fprintf( ioQQQ, " ion_wrapper returns; %s fracs = ", elementnames.chElementSym[nelem] );
01506 for( long ion = 0; ion<=nelem+1; ion++ )
01507 fprintf( ioQQQ,"%10.3e ", dense.xIonDense[nelem][ion]/dense.gas_phase[nelem] );
01508 fprintf( ioQQQ, "\n" );
01509 }
01510
01511 ASSERT(lgElemsConserved());
01512
01513 return;
01514 }