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