cloudy  trunk
ion_solver.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*ion_solver solve the bi-diagonal matrix for ionization balance */
4 #include "cddefines.h"
5 #include "yield.h"
6 #include "prt.h"
7 #include "continuum.h"
8 #include "iso.h"
9 #include "dynamics.h"
10 #include "grainvar.h"
11 #include "hmi.h"
12 #include "mole.h"
13 #include "thermal.h"
14 #include "newton_step.h"
15 #include "thirdparty.h"
16 #include "conv.h"
17 #include "secondaries.h"
18 #include "phycon.h"
19 #include "atmdat.h"
20 #include "heavy.h"
21 #include "elementnames.h"
22 #include "dense.h"
23 #include "radius.h"
24 #include "ionbal.h"
25 #include "taulines.h"
26 #include "trace.h"
27 
28 #define SOLVE_TWO 0
29 //bool lgOH_ChargeTransferDominant( void );
30 
31 STATIC bool lgTrivialSolution( long nelem, double abund_total );
32 
33 STATIC void find_solution( long nelem, long ion_range, valarray<double> &xmat, valarray<double> &source);
34 
35 STATIC void fill_array( long int nelem, long ion_range, valarray<double> &xmat, valarray<double> &source, valarray<double> &auger, double *abund_total );
36 
37 #if SOLVE_TWO
38 STATIC void combine_arrays( valarray<double> &xmat, const valarray<double> &xmat1, const valarray<double> &xmat2, long ion_range1, long ion_range2 );
39 #endif
40 
41 STATIC double get_total_abundance_ions( long int nelem );
42 
43 STATIC void HomogeneousSource( long nelem, long ion_low, long ion_range, valarray<double> &xmat, valarray<double> &source, double abund_total );
44 
45 STATIC void store_new_densities( long nelem, long ion_range, long ion_low, double *source, double abund_total, bool *lgNegPop ) ;
46 
47 STATIC void PrintRates( long nelem, bool lgNegPop, double abund_total, valarray<double> &auger, bool lgPrintIt );
48 
49 void solveions(double *ion, double *rec, double *snk, double *src,
50  long int nlev, long int nmax);
51 
52  /* this will be used to address the 2d arrays */
53 # undef MAT
54 # define MAT(M_,I_,J_) ((M_)[(I_)*(ion_range)+(J_)])
55 
56 # undef MAT1
57 # define MAT1(M_,I_,J_) ((M_)[(I_)*(ion_range1)+(J_)])
58 
59 # undef MAT2
60 # define MAT2(M_,I_,J_) ((M_)[(I_)*(ion_range2)+(J_)])
61 
62 void ion_solver( long int nelem, bool lgPrintIt)
63 {
64  double abund_total = 0.0;
65  long ion_low = dense.IonLow[nelem];
66  bool lgNegPop = false;
67  double error = 0.0;
68 
69  DEBUG_ENTRY( "ion_solver()" );
70 
71  iso_charge_transfer_update(nelem);
72 
73  long ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
74  valarray<double> xmat(ion_range*ion_range);
75  valarray<double> source(ion_range);
76  valarray<double> auger(LIMELM+1);
77 
78  // V/W cycle would do ion solution *first*, then iso, then
79  // ion again if iso made any changes
80  // Seem to need to break out after iso at the moment. Why?
81 
82  conv.incrementCounter(ION_SOLVES);
83  for (long it=0; it<4; ++it)
84  {
85  conv.incrementCounter(ISO_LOOPS);
86  for( long ipISO=ipH_LIKE; ipISO<MIN2(NISO,nelem+1); ipISO++ )
87  {
88  if( (dense.IonHigh[nelem] >= nelem - ipISO) &&
89  (dense.IonLow[nelem] <= nelem - ipISO))
90  {
91  iso_set_ion_rates(ipISO, nelem);
92  }
93  }
94 
95  //if ( it > 1 && error > 0.0)
96  // fprintf(ioQQQ,"ion nelem %ld loop %ld error %g nzone %ld\n",nelem,it,error,nzone);
97 
98  abund_total = get_total_abundance_ions( nelem );
99 
100  bool lgTrivial = lgTrivialSolution( nelem, abund_total );
101  if( !lgTrivial )
102  {
103  // fill xmat and source with appropriate terms
104  fill_array( nelem, ion_range, xmat, source, auger, &abund_total );
105 
106  // decide if matrix is homogeneous
107  HomogeneousSource( nelem, ion_low, ion_range, xmat, source, abund_total );
108 
109  // Now find the solution
110  find_solution( nelem, ion_range, xmat, source);
111 
112  // save the results in the global density variables
113  store_new_densities( nelem, ion_range, ion_low, &source[0], abund_total, &lgNegPop );
114 
115  ASSERT(ion_range >= dense.IonHigh[nelem]-dense.IonLow[nelem]+1);
116  if (ion_range != dense.IonHigh[nelem]-dense.IonLow[nelem]+1)
117  {
118  ion_range = dense.IonHigh[nelem]-dense.IonLow[nelem]+1;
119  xmat.resize(ion_range*ion_range);
120  source.resize(ion_range);
121  }
122 
123  }
124 
125  error = 0.0;
126  /* update and solve iso levels */
127  for( long ipISO=ipH_LIKE; ipISO<MIN2(NISO,nelem+1); ipISO++ )
128  {
129  double thiserror;
130  iso_solve( ipISO, nelem, thiserror );
131  if (thiserror > error)
132  error = thiserror;
133  }
134  if ( error < 1e-4 )
135  break;
136  }
137 
138  iso_satellite_update(nelem);
139 
140  for( long ipISO=ipH_LIKE; ipISO<MIN2(NISO,nelem+1); ipISO++ )
141  {
142  /* now evaluate departure coefficients */
143  iso_departure_coefficients( ipISO, nelem );
144  }
145 
146  if( prt.lgPrtArry[nelem] || lgPrintIt )
147  PrintRates( nelem, lgNegPop, abund_total, auger, lgPrintIt );
148 
149  return;
150 }
151 
152 #if SOLVE_TWO
153 // ion_solver is overloaded - this version looks for a simultaneous solution to the ionization balance of two elements
154 void ion_solver( long int nelem1, long int nelem2, bool lgPrintIt)
155 {
156  bool lgHomogeneous = true;
157  bool lgNegPop1 = false;
158  bool lgNegPop2 = false;
159 
160  DEBUG_ENTRY( "ion_solver()" );
161 
162  long ion_low1 = dense.IonLow[nelem1];
163  long ion_low2 = dense.IonLow[nelem2];
164  long ion_range1 = dense.IonHigh[nelem1]-dense.IonLow[nelem1]+1;
165  long ion_range2 = dense.IonHigh[nelem2]-dense.IonLow[nelem2]+1;
166  long ion_range = ion_range1 + ion_range2;
167  double abund_total1 = get_total_abundance_ions( nelem1 );
168  double abund_total2 = get_total_abundance_ions( nelem2 );
169  valarray<double> xmat(ion_range*ion_range);
170  valarray<double> xmat1(ion_range1*ion_range1);
171  valarray<double> xmat2(ion_range2*ion_range2);
172  valarray<double> source(ion_range);
173  valarray<double> auger1(LIMELM+1);
174  valarray<double> auger2(LIMELM+1);
175 
176 #define ENABLE_SIMULTANEOUS_SOLUTION 0
177 
178  if( lgTrivialSolution( nelem1, abund_total1 ) )
179  {
180  // nelem2 has only an ancillary role in solving nelem1. We do not need to solve nelem2 balance if nelem1 is trivial
181  return;
182  }
183  else if( lgTrivialSolution( nelem2, abund_total2 ) || !ENABLE_SIMULTANEOUS_SOLUTION )
184  {
185  // if nelem2 has a trivial solution, we still need to solve nelem1, just do it by itself
186  ion_solver( nelem1, lgPrintIt );
187  return;
188  }
189 
190  /* don't do simultaneous if one or both does not include the neutral stage. */
191  if( dense.IonLow[nelem1]>0 || dense.IonLow[nelem2]>0 || nelem1!=ipOXYGEN || nelem2!=ipHYDROGEN )
192  {
193  fprintf( ioQQQ, "This routine is currently intended to do only O and H when both have significant neutral fractions.\n" );
194  fprintf( ioQQQ, "It should be generalized further for other cases. Exiting. Sorry.\n" );
195  cdEXIT( EXIT_FAILURE );
196  }
197 
198  ASSERT( nelem1 != nelem2 );
199 
200  // fill active element's array
201  fill_array( nelem1, ion_range1, xmat1, source, auger1, &abund_total1 );
202  // fill passive element's array
203  fill_array( nelem2, ion_range2, xmat2, source, auger2, &abund_total2 );
204 
205  combine_arrays( xmat, xmat1, xmat2, ion_range1, ion_range2 );
206 
207  // force homogeneity in simultaneous solutions
208  for( long i=0; i< ion_range; ++i )
209  source[i] = 0.;
210 
211  // replace neutral stages with abundance total
212  for( long i=0; i< ion_range1; ++i )
213  MAT( xmat, i, 0 ) = 1.;
214 
215  for( long i=ion_range1; i< ion_range; ++i )
216  MAT( xmat, i, ion_range1 ) = 1.;
217 
218  source[0] = dense.xIonDense[nelem1][0] + dense.xIonDense[nelem1][1];
219  source[ion_range1] = dense.xIonDense[nelem2][0] + dense.xIonDense[nelem2][1];
220 
221  // declare homogeneous.
222  lgHomogeneous = true;
223 
224  // Now find the solution
225  find_solution( nelem1, ion_range, xmat, source);
226 
227  // save the results in the global density variables
228  store_new_densities( nelem1, ion_range1, ion_low1, &source[0], abund_total1, &lgNegPop1 );
229  store_new_densities( nelem2, ion_range2, ion_low2, &source[ion_range1], abund_total2, &lgNegPop2 );
230  ASSERT( lgNegPop2 == false );
231 
232  if( prt.lgPrtArry[nelem1] || lgPrintIt )
233  PrintRates( nelem1, lgNegPop1, abund_total1, auger1, lgPrintIt );
234 
235  return;
236 }
237 #endif
238 
239 STATIC bool lgTrivialSolution( long nelem, double abund_total )
240 {
241  double renorm;
242  /* return if IonHigh is zero, since no ionization at all */
243  if( dense.IonHigh[nelem] == dense.IonLow[nelem] )
244  {
245  /* set the atom to the total gas phase abundance */
246  dense.xIonDense[nelem][dense.IonHigh[nelem]] = abund_total;
247  for ( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
248  iso_renorm(nelem, ipISO, renorm);
249  return true;
250  }
251  else if( dense.lgSetIoniz[nelem] )
252  {
253  /* option to force ionization distribution with element name ioniz */
254  for( long ion=0; ion<nelem+2; ++ion )
255  dense.xIonDense[nelem][ion] = dense.SetIoniz[nelem][ion]*abund_total;
256  for ( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
257  iso_renorm(nelem, ipISO, renorm);
258  return true;
259  }
260  else
261  return false;
262 }
263 
264 STATIC void find_solution( long nelem, long ion_range, valarray<double> &xmat, valarray<double> &source)
265 {
266  int32 nerror;
267  valarray<double> xmatsave(ion_range*ion_range);
268  valarray<double> sourcesave(ion_range);
269  valarray<int32> ipiv(ion_range);
270 
271  DEBUG_ENTRY("find_solution()");
272 
273  // save copy of xmat before continuing.
274  for( long i=0; i< ion_range; ++i )
275  {
276  sourcesave[i] = source[i];
277  for( long j=0; j< ion_range; ++j )
278  {
279  MAT( xmatsave, i, j ) = MAT( xmat, i, j );
280  }
281  }
282 
283  if (1)
284  {
285  nerror = solve_system(xmat,source,ion_range,NULL);
286 
287  if (nerror)
288  {
289  // solve_system doesn't return a valid solution, so just
290  // provide initial values to keep things running
291  fprintf(ioQQQ,"Error for nelem %ld, active ion range %ld--%ld\n",
292  nelem,dense.IonLow[nelem],dense.IonHigh[nelem]);
293  fprintf(ioQQQ,"Initial ion abundances:");
294  for( long j=0; j<nelem+2; ++j )
295  fprintf(ioQQQ," %g",dense.xIonDense[nelem][j]);
296  fprintf(ioQQQ,"\n");
297  for( long j=0; j<ion_range; ++j )
298  source[j] = dense.xIonDense[nelem][dense.IonLow[nelem]+j];
299  }
300  }
301  else
302  {
303  /* this is the default solver - now get new solution */
304  nerror = 0;
305  /* Use general matrix solver */
306  getrf_wrapper(ion_range, ion_range, &xmat[0], ion_range, &ipiv[0], &nerror);
307  if( nerror != 0 )
308  {
309  fprintf( ioQQQ,
310  " DISASTER ion_solver: dgetrf finds singular or ill-conditioned matrix nelem=%li %s ion_range=%li, nConv %li IonLow %li IonHi %li\n",
311  nelem ,
312  elementnames.chElementSym[nelem],
313  ion_range,
314  conv.nTotalIoniz ,
315  dense.IonLow[nelem], dense.IonHigh[nelem]);
316  fprintf( ioQQQ, " xmat follows\n");
317  for( long i=0; i<ion_range; ++i )
318  {
319  for( long j=0;j<ion_range;j++ )
320  {
321  fprintf(ioQQQ,"%e\t",MAT(xmatsave,j,i));
322  }
323  fprintf(ioQQQ,"\n");
324  }
325  fprintf(ioQQQ,"source follows\n");
326  for( long i=0; i<ion_range;i++ )
327  {
328  fprintf(ioQQQ,"%e\t",sourcesave[i]);
329  }
330  fprintf(ioQQQ,"\n");
331  cdEXIT(EXIT_FAILURE);
332  }
333  getrs_wrapper('N', ion_range, 1, &xmat[0], ion_range, &ipiv[0], &source[0], ion_range, &nerror);
334  if( nerror != 0 )
335  {
336  fprintf( ioQQQ, " DISASTER ion_solver: dgetrs finds singular or ill-conditioned matrix nelem=%li ionrange=%li\n",
337  nelem , ion_range );
338  cdEXIT(EXIT_FAILURE);
339  }
340  }
341 
342  {
343  /* this is to debug following failed assert */
344  enum {DEBUG_LOC=false};
345  if( DEBUG_LOC && nelem == ipHYDROGEN )
346  {
347  fprintf(ioQQQ,"debuggg ion_solver1 %ld\t%.2f\t%.4e\t%.4e\tIon\t%.3e\tRec\t%.3e\n",
348  nelem,
349  fnzone,
350  phycon.te,
351  dense.eden,
352  ionbal.RateIonizTot(nelem,0) ,
353  ionbal.RateRecomTot[nelem][0]);
354  fprintf(ioQQQ," Msrc %.3e %.3e\n", mole.source[nelem][0], mole.source[nelem][1]);
355  fprintf(ioQQQ," Msnk %.3e %.3e\n", mole.sink[nelem][0], mole.sink[nelem][1]);
356  fprintf(ioQQQ," Poprat %.3e nomol %.3e\n",source[1]/source[0],
357  ionbal.RateIonizTot(nelem,0)/ionbal.RateRecomTot[nelem][0]);
358  }
359  }
360 
361  for( long i=0; i<ion_range; i++ )
362  {
363  ASSERT( !isnan( source[i] ) );
364  ASSERT( source[i] < MAX_DENSITY );
365  }
366 
367  return;
368 }
369 
370 bool lgOH_ChargeTransferDominant( void )
371 {
372 #define THRESHOLD 0.25
373 
374  bool lgDominant = false;
375  double denominator;
376 
377  if( (denominator = ionbal.RateIonizTot(ipHYDROGEN,0) + atmdat.CharExcIonTotal[ipHYDROGEN] + mole.sink[ipHYDROGEN][0] ) > 0. )
378  {
379  lgDominant = lgDominant ||
380  ( atmdat.CharExcRecTo[ipHYDROGEN][ipOXYGEN][0]*dense.xIonDense[ipOXYGEN][1]/denominator > THRESHOLD );
381  }
382 
383  if( (denominator = ionbal.RateRecomTot[ipHYDROGEN][0] + atmdat.CharExcRecTotal[ipHYDROGEN] ) > 0. )
384  {
385  lgDominant = lgDominant ||
386  ( atmdat.CharExcIonOf[ipHYDROGEN][ipOXYGEN][0]*dense.xIonDense[ipOXYGEN][0]/denominator > THRESHOLD );
387  }
388 
389  return lgDominant;
390 }
391 
392 #if SOLVE_TWO
393 STATIC void combine_arrays( valarray<double> &xmat, const valarray<double> &xmat1, const valarray<double> &xmat2, long ion_range1, long ion_range2 )
394 {
395  DEBUG_ENTRY( "combine_arrays()" );
396 
397  long ion_range = ion_range1 + ion_range2;
398 
399  for( long i=0; i<ion_range1; i++ )
400  for( long j=0; j<ion_range1; j++ )
401  MAT( xmat, i, j ) = MAT1( xmat1, i, j );
402 
403  for( long i=0; i<ion_range2; i++ )
404  for( long j=0; j<ion_range2; j++ )
405  MAT( xmat, i+ion_range1, j+ion_range1 ) = MAT2( xmat2, i, j );
406 
407 #if 0
408  bool lgNoDice = false;
409  for( long i=0; i<ion_range1; i++ )
410  {
411  if( !fp_equal( MAT( xmat, i, 0), 1.0, 1 ) )
412  {
413  lgNoDice = true;
414  break;
415  }
416  }
417 
418  if( !lgNoDice )
419  {
420  for( long i=ion_range1; i<ion_range; i++ )
421  MAT( xmat, i, 0 ) = 1.0;
422  }
423 #endif
424 
425  return;
426 }
427 #endif
428 
429 STATIC void store_new_densities( long nelem, long ion_range, long ion_low, double *source, double abund_total, bool *lgNegPop )
430 {
431  DEBUG_ENTRY( "store_new_densities()" );
432 
433  ASSERT( nelem >= 0 );
434  ASSERT( nelem < LIMELM );
435  ASSERT( ion_range <= nelem + 2 );
436  ASSERT( ion_low >= 0 );
437  ASSERT( ion_low <= nelem + 1 );
438 
439 #if 0
440 # define RJRW 0
441  if( RJRW && 0 )
442  {
443  /* verify that the rates are sensible */
444  double test;
445  for(long i=0; i<ion_range; i++) {
446  test = 0.;
447  for(long j=0; j<ion_range; j++) {
448  test = test+source[j]*MAT(xmatsave,j,i);
449  }
450  fprintf(ioQQQ,"%e\t",test);
451  }
452  fprintf(ioQQQ,"\n");
453 
454  test = 0.;
455  fprintf(ioQQQ," ion %li abundance %.3e\n",nelem,abund_total);
456  for( long ion=dense.IonLow[nelem]; ion < dense.IonHigh[nelem]; ion++ )
457  {
458  if( ionbal.RateRecomTot[nelem][ion] != 0 && source[ion-ion_low] != 0 )
459  fprintf(ioQQQ," %li %.3e %.3e : %.3e\n",ion,source[ion-ion_low],
460  source[ion-ion_low+1]/source[ion-ion_low],
461  ionbal.RateIonizTot(nelem,ion)/ionbal.RateRecomTot[nelem][ion]);
462  else
463  fprintf(ioQQQ," %li %.3e [One ratio infinity]\n",ion,source[ion-ion_low]);
464  test += source[ion-ion_low];
465  }
466  }
467 #endif
468 
469  /*
470  * >> chng 03 jan 15 rjrw:- terms are now included for
471  * molecular sources and sinks.
472  *
473  * When the network is not in equilibrium, this will lead to a
474  * change in the derived abundance of coupled ions after the matrix
475  * solution.
476  *
477  * We therefore renormalize to keep the total abundance of the
478  * states treated by the ionization ladder constant -- only the
479  * molecular network is allowed to change this.
480  *
481  * The difference between `renorm' and 1. is a measure of the
482  * quality of the solution (it will be 1. if the rate of transfer
483  * into the ionization ladder species balances the rate of transfer
484  * out, for the consistent relative abundances)
485  *
486  */
487 
488  /* source[i] contains new solution for ionization populations
489  * save resulting abundances into main ionization density array,
490  * while checking whether any negative level populations occurred */
491  *lgNegPop = false;
492  for( long i=0; i < ion_range; i++ )
493  {
494  long ion = i+ion_low;
495 
496  if( source[i] < 0. )
497  {
498  /* >>chng 04 dec 04, put test on neg abund here, don't print unless value is very -ve */
499  /* >>chng 06 feb 28, from -1e-10 to -1e-9, sim func_t10 had several negative
500  * during initial search, due to extremely high ionization */
501  /* >>chng 06 mar 11, from 1e-9 to 2e-9 make many struc elements floats from double */
502  if( source[i]<-2e-9 )
503  {
504  fprintf(ioQQQ,
505  " PROBLEM negative ion population in ion_solver, nelem=%li, %s ion=%li val=%.3e Search?%c zone=%li iteration=%li\n",
506  nelem ,
507  elementnames.chElementSym[nelem],
508  ion ,
509  source[i] ,
510  TorF(conv.lgSearch) ,
511  nzone ,
512  iteration );
513  *lgNegPop = true;
514  fixit(); //break PrintRates into one NegPop case and one trace? No auger defined here.
515  //PrintRates( nelem, *lgNegPop, abund_total, auger );
516  }
517 
518  fixit(); // Kill this bit and force exit on negative populations.
519 #if 1
520  source[i] = 0.;
521  /* if this is one of the iso seq model atoms then must also zero out pops */
522  //if( ion == nelem+1-NISO ) //newmole had this should have been next line?
523  if( ion > nelem-NISO && ion < nelem + 1 )
524  {
525  long int ipISO = nelem - ion;
526  ASSERT( ipISO>=0 && ipISO<NISO );
527  for( long level = 0; level < iso_sp[ipISO][nelem].numLevels_max; level++ )
528  iso_sp[ipISO][nelem].st[level].Pop() = 0.;
529  }
530 #endif
531  }
532  }
533 
534  double renormnew = 1.;
535  {
536  double dennew = 0.;
537 
538  /* find total population to renorm - also here check that negative pops did not occur */
539  for( long i=0;i < ion_range; i++ )
540  {
541  dennew += source[i];
542  }
543 
544  if(dennew > 0.)
545  {
546  renormnew = abund_total / dennew;
551  }
552  else
553  {
554  renormnew = 1.;
555  }
556 
557  dennew = 0.;
558  for( long i=0;i < ion_range; i++ )
559  {
560  source[i] *= renormnew;
561  dennew += source[i];
562  }
563  }
564  /* check not negative, should be +ve, can be zero if species has become totally molecular.
565  * this happens for hydrogen if no cosmic rays, or cr ion set very low */
566  if( renormnew < 0)
567  {
568  fprintf(ioQQQ,"PROBLEM impossible value of renorm \n");
569  }
570  ASSERT( renormnew>=0 );
571 
572  for( long i=0; i < ion_range; i++ )
573  {
574  long ion = i+ion_low;
575  dense.xIonDense[nelem][ion] = source[i];
576  if( dense.xIonDense[nelem][ion] >= MAX_DENSITY )
577  {
578  fprintf( ioQQQ, "PROBLEM DISASTER Huge density in ion_solver, nelem %ld ion %ld source %e renormnew %e\n",
579  nelem, ion, source[i], renormnew );
580  }
581  ASSERT( dense.xIonDense[nelem][ion] < MAX_DENSITY );
582  }
583 
584  fixit(); // this should only be done if trimming is not disabled?
585 
586  /* Zero levels with abundances < 1e-25 which which will suffer numerical noise */
587  while( dense.IonHigh[nelem] > dense.IonLow[nelem] &&
588  dense.xIonDense[nelem][dense.IonHigh[nelem]] < 1e-25*abund_total &&
589  dense.IonHigh[nelem] > 1)
590  {
591  ASSERT( dense.xIonDense[nelem][dense.IonHigh[nelem]] >= 0. );
592  /* zero out abundance and heating due to stage of ionization we are about to zero out */
593  dense.xIonDense[nelem][dense.IonHigh[nelem]] = 0.;
594  thermal.heating[nelem][dense.IonHigh[nelem]-1] = 0.;
595  /* decrement counter */
596  --dense.IonHigh[nelem];
597  }
598 
599  for ( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
600  {
601  double renorm;
602  iso_renorm(nelem, ipISO, renorm);
603  }
604 
605  // sanity check, either offset stages of low and high ionization
606  ASSERT( (dense.IonLow[nelem] <= dense.IonHigh[nelem]) ||
607  // both totally neutral
608  (dense.IonLow[nelem]==0 && dense.IonHigh[nelem]==0 ) ||
609  // both fully stripped
610  ( dense.IonLow[nelem]==nelem+1 && dense.IonHigh[nelem]==nelem+1 ) );
611 
612  fixit(); //this routine does not ever set *lgNegPop true. Think it does on newmole though.
613  //return *lgNegPop;
614  return;
615 }
616 
617 STATIC double get_total_abundance_ions( long int nelem )
618 {
619  DEBUG_ENTRY( "get_total_abundance_ions()" );
620 
621  ASSERT( nelem >= 0 );
622  ASSERT( nelem < LIMELM );
623 
624  ionbal.elecsnk[nelem] = 0.;
625  ionbal.elecsrc[nelem] = 0.;
626 
627  double abund_total = 0.;
628  for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
629  {
630  abund_total += dense.xIonDense[nelem][ion];
631  }
632 
633  realnum tot1 = dense.gas_phase[nelem];
634  realnum tot2 = (realnum)(dense.xMolecules(nelem)+abund_total);
635  if (0)
636  {
637  if (fabs(tot2-tot1) > conv.GasPhaseAbundErrorAllowed*tot1)
638  fprintf(ioQQQ,"%ld %13.8g %13.8g %13.8g %13.8g\n",nelem,tot1,tot2,abund_total,tot2/tot1 - 1.0);
639  }
640 
641  ASSERT( fp_equal_tol(tot1, tot2, realnum(conv.GasPhaseAbundErrorAllowed*tot1 + 100.f*FLT_MIN)) );
642 
643  ASSERT( abund_total < MAX_DENSITY );
644 
645  return abund_total;
646 }
647 
648 STATIC void fill_array( long int nelem, long ion_range, valarray<double> &xmat, valarray<double> &source, valarray<double> &auger, double *abund_total )
649 {
650  long int limit,
651  IonProduced;
652  double rateone;
653  long ion_low;
654 
655  valarray<double> sink(ion_range);
656  valarray<int32> ipiv(ion_range);
657 
658  DEBUG_ENTRY( "fill_array()" );
659 
660  /* this is on the c scale, so H is 0 */
661  ASSERT( nelem >= 0);
662  ASSERT( dense.IonLow[nelem] >= 0 );
663  ASSERT( dense.IonHigh[nelem] >= 0 );
664 
665  /* impossible for HIonFrac[nelem] to be zero if IonHigh(nelem)=nelem+1
666  * HIonFrac(nelem) is stripped to hydrogen */
667  /* >>chng 01 oct 30, to assert */
668  //ASSERT( (dense.IonHigh[nelem] < nelem + 1) || dense.xIonDense[nelem][nelem+1-ipH_LIKE] > 0. );
669 
670  /* zero out the ionization and recombination rates that we will modify here,
671  * but not the iso-electronic stages which are done elsewhere,
672  * the nelem stage of ionization is he-like,
673  * the nelem+1 stage of ionization is h-like */
674 
675  /* loop over stages of ionization that we solve for here,
676  * up through and including one less than nelem-NISO,
677  * never actually do highest NISO stages of ionization since they
678  * come from the ionization ratio from the next lower stage */
679  limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);
680 
681  /* the full range of ionization - this is number of ionization stages */
682  ASSERT( ion_range <= nelem+2 );
683 
684  ion_low = dense.IonLow[nelem];
685  for( long i=0; i<ion_range;i++ )
686  {
687  source[i] = 0.;
688  }
689 
690  /* zero-out loop comes before main loop since there are off-diagonal
691  * elements in the main ionization loop, due to multi-electron processes,
692  * TotIonizRate and TotRecom were already set in h-like and he-like solvers
693  * other recombination rates were already set by routines responsible for them */
694  for( long ion_from=0; ion_from <= limit; ion_from++ )
695  {
696  for( long ion_to=0; ion_to < nelem+2; ion_to++ )
697  {
698  ionbal.RateIoniz[nelem][ion_from][ion_to] = 0.;
699  }
700  }
701 
702  /* auger is used only for debug printout - it is special because with multi-electron
703  * Auger ejection, very high stages of ionization can be produced, well beyond
704  * the highest stage considered here. so we malloc to the limit */
705  for( long i=0; i< LIMELM+1; ++i )
706  {
707  auger[i] = 0.;
708  }
709 
710  /* zero out xmat */
711  for( long i=0; i< ion_range; ++i )
712  {
713  for( long j=0; j< ion_range; ++j )
714  {
715  MAT( xmat, i, j ) = 0.;
716  }
717  }
718 
719  {
720  /* this sets up a fake ionization balance problem, with a trivial solution,
721  * for debugging the ionization solver */
722  enum {DEBUG_LOC=false};
723  if( DEBUG_LOC && nelem==ipCARBON )
724  {
725  broken();/* within optional debug print statement */
726  dense.IonLow[nelem] = 0;
727  dense.IonHigh[nelem] = 3;
728  *abund_total = 1.;
729  /* make up ionization and recombination rates */
730  for( long ion=dense.IonLow[nelem]; ion <= limit; ion++ )
731  {
732  double fac=1;
733  if(ion)
734  fac = 1e-10;
735  ionbal.RateRecomTot[nelem][ion] = 100.;
736  for( long ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
737  {
738  /* direct photoionization of this shell */
739  ionbal.PhotoRate_Shell[nelem][ion][ns][0] = fac;
740  }
741  }
742  }
743  }
744 
745  bool lgGrainsOn = gv.lgDustOn() && ionbal.lgGrainIonRecom && gv.lgGrainPhysicsOn;
746 
747  /* Now put in all recombination and ionization terms from CO_mole() that
748  * come from molecular reactions. this traces molecular process that
749  * change ionization stages with this ladder - but do not remove from
750  * the ladder */
751  for( long ion_to=dense.IonLow[nelem]; ion_to <= dense.IonHigh[nelem]; ion_to++ )
752  {
753  for( long ion_from=dense.IonLow[nelem]; ion_from <= dense.IonHigh[nelem]; ++ion_from )
754  {
755  /* do not do ion onto itself */
756  if( ion_to != ion_from )
757  {
758  /* CT with molecules */
759  rateone = mole.xMoleChTrRate[nelem][ion_from][ion_to] * atmdat.lgCTOn;
760  MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
761  MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
762  /* CT with grains */
763  rateone = gv.GrainChTrRate[nelem][ion_from][ion_to]*lgGrainsOn;
764  MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= rateone;
765  MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += rateone;
766  }
767  }
768  }
769 
770  for( long ion=dense.IonLow[nelem]; ion <= limit; ion++ )
771  {
772  /* thermal & secondary collisional ionization */
773  ionbal.RateIoniz[nelem][ion][ion+1] +=
774  ionbal.CollIonRate_Ground[nelem][ion][0] +
775  secondaries.csupra[nelem][ion] +
776  /* inner shell ionization by UTA lines */
777  ionbal.UTA_ionize_rate[nelem][ion];
778 
779  /* loop over all atomic sub-shells to include photoionization */
780  for( long ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
781  {
782  /* this is the primary ionization rate - add to diagonal element,
783  * test on ion stage is so that we don't include ionization from the very highest
784  * ionization stage to even higher - since those even higher stages are not considered
785  * this would appear as a sink - but populations of this highest level is ensured to
786  * be nearly trivial and neglecting it production of even higher ionization OK */
787  /* >>chng 04 nov 29 RJRW, include following in this branch so only
788  * evaluated when below ions done with iso-sequence */
789  if( ion+1-ion_low < ion_range )
790  {
791  /* this will be redistributed into charge states in following loop */
792 
793  /* t_yield::Inst().nelec_eject(nelem,ion,ns) is total number of electrons that can
794  * possibly be freed
795  * loop over nej, the number of electrons ejected including the primary,
796  * nej = 1 is primary, nej > 1 includes primary plus Auger
797  * t_yield::Inst().elec_eject_frac is prob of nej electrons */
798  for( long nej=1; nej <= t_yield::Inst().nelec_eject(nelem,ion,ns); nej++ )
799  {
800  /* this is the ion that is produced by this ejection,
801  * limited by highest possible stage of ionization -
802  * do not want to ignore ionization that go beyond this */
803  IonProduced = MIN2(ion+nej,dense.IonHigh[nelem]);
804  rateone = ionbal.PhotoRate_Shell[nelem][ion][ns][0]*
805  t_yield::Inst().elec_eject_frac(nelem,ion,ns,nej-1);
806 
807  /* direct photoionization of this shell */
808  ionbal.RateIoniz[nelem][ion][IonProduced] += rateone;
809 
810  /* only used for possible printout - multiple electron Auger rate -
811  * do not count one-electron as Auger */
812  if( nej>1 )
813  auger[IonProduced-1] += rateone;
814  }
815  }
816  }
817  }
818 
819  for( long ion=dense.IonLow[nelem]; ion < dense.IonHigh[nelem]; ion++ )
820  {
821  /* this is charge transfer ionization of this species by hydrogen and helium */
822  double ction;
823  long ipISO=nelem-ion;
824 
825  if ( nelem < t_atmdat::NCX && nelem == ipISO )
826  {
827  ction = atmdat.CharExcIonTotal[nelem] * iso_sp[ipISO][nelem].st[0].Pop() / SDIV(dense.xIonDense[nelem][nelem-ipISO]);
828  }
829  else
830  {
831  ction=0;
832  for (long nelem1=0; nelem1 < t_atmdat::NCX; ++nelem1)
833  ction += atmdat.CharExcIonOf[nelem1][nelem][ion]*dense.xIonDense[nelem1][1];
834  }
835  //ionbal.elecsrc[nelem] += ction*dense.xIonDense[nelem][ion];
836  /* depopulation processes enter with negative sign */
837  MAT( xmat, ion-ion_low, ion-ion_low ) -= ction;
838  MAT( xmat, ion-ion_low, ion+1-ion_low ) += ction;
839  }
840 
841  for( long ion=dense.IonLow[nelem]; ion < dense.IonHigh[nelem]; ion++ )
842  {
843  /* this is charge transfer ionization of this species by hydrogen and helium */
844  double ctrec;
845  long ipISO=nelem-ion;
846 
847  if ( nelem==ipHYDROGEN )
848  ctrec = atmdat.CharExcRecTotal[nelem];
849  else if( nelem==ipHELIUM && ipISO==ipHE_LIKE )
850  ctrec = atmdat.CharExcRecTotal[nelem];
851  else
852  {
853  ctrec = 0.;
854  for (long nelem1=0; nelem1<t_atmdat::NCX; ++nelem1)
855  {
856  long ipISO = nelem1;
857  ctrec +=
858  atmdat.CharExcRecTo[nelem1][nelem][ion]*iso_sp[ipISO][nelem1].st[0].Pop();
859  }
860  }
861 
862  MAT( xmat, ion+1-ion_low, ion+1-ion_low ) -= ctrec;
863  MAT( xmat, ion+1-ion_low, ion-ion_low ) += ctrec;
864  //ionbal.elecsnk[nelem] += ctrec;
865  }
866 
867 
868  for( long ion_from=dense.IonLow[nelem]; ion_from < dense.IonHigh[nelem]; ion_from++ )
869  {
870  for( long ion_to=ion_from+1; ion_to <= dense.IonHigh[nelem]; ion_to++ )
871  {
872  ionbal.elecsrc[nelem] += ionbal.RateIoniz[nelem][ion_from][ion_to]*dense.xIonDense[nelem][ion_from]*
873  (ion_to-ion_from);
874  /* depopulation processes enter with negative sign */
875  MAT( xmat, ion_from-ion_low, ion_from-ion_low ) -= ionbal.RateIoniz[nelem][ion_from][ion_to];
876  MAT( xmat, ion_from-ion_low, ion_to-ion_low ) += ionbal.RateIoniz[nelem][ion_from][ion_to];
877  }
878  }
879 
880  for( long ion=dense.IonLow[nelem]; ion<dense.IonHigh[nelem]; ion++ )
881  {
882  /* loss of next higher ion due to recombination to this ion stage */
883  MAT( xmat, ion+1-ion_low, ion+1-ion_low ) -= ionbal.RateRecomTot[nelem][ion];
884  MAT( xmat, ion+1-ion_low, ion-ion_low ) += ionbal.RateRecomTot[nelem][ion];
885  ionbal.elecsnk[nelem] += ionbal.RateRecomTot[nelem][ion]*dense.xIonDense[nelem][ion+1];
886  }
887 
888  for( long i=0; i<ion_range;i++ )
889  {
890  long ion = i+ion_low;
891 
892  /* these are the external source and sink terms */
893  /* need negative sign to get positive pops */
894  source[i] -= mole.source[nelem][ion];
895  MAT( xmat, i, i ) -= mole.sink[nelem][ion];
896  }
897 
898  return;
899 }
900 
901 STATIC void HomogeneousSource( long nelem, long ion_low, long ion_range, valarray<double> &xmat, valarray<double> &source, double abund_total )
902 {
903  bool lgHomogeneous = true;
904  double totsrc,
905  totscl,
906  maxdiag;
907 
908  DEBUG_ENTRY( "lgHomogeneousSource()" );
909 
910  /* this will be sum of source and sink terms, will be used to decide if
911  * matrix is singular */
912  totsrc = totscl = maxdiag = 0.;
913  for( long i=0; i<ion_range;i++ )
914  {
915  long ion = i + ion_low;
916 
917  totsrc -= source[i];
918  fixit();
919  // In old newmole, totscl was calculated before
920  // mole.sink terms are added into xmat, but those terms are now already done.
921  // the above hack takes care of that, but a cleaner solution is needed.
922  double diag = -(MAT( xmat, i, i )+mole.sink[nelem][ion]);
923  if (diag > maxdiag)
924  maxdiag = diag;
925  totscl += diag*dense.xIonDense[nelem][ion];
926  }
927 
928  // Largest condition number before we decide that conservation will get
929  // a better answer than the raw linear solver
930  const double CONDITION_SCALE = 1e8;
931  double conserve_scale = maxdiag/CONDITION_SCALE;
932  ASSERT( conserve_scale > 0.0 );
933 
934  /* matrix is not homogeneous if source is non-zero */
935  if( totsrc > 1e-10*totscl )
936  lgHomogeneous = false;
937 
938  fixit(); // dynamics rates need to be moved into fill_array?
939  /* chng 03 jan 13 rjrw, add in dynamics if required here,
940  * - only do advection if we have not overrun the radius scale */
941  if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection && !dynamics.lgEquilibrium
942  && dynamics.Rate != 0. )
943  {
944  for( long i=0; i<ion_range;i++ )
945  {
946  long ion = i+ion_low;
947  MAT( xmat, i, i ) -= dynamics.Rate;
948  source[i] -= dynamics.Source[nelem][ion];
949  /* fprintf(ioQQQ," %li %li %.3e (%.3e %.3e)\n",i,i,MAT(*xmat,i,i),
950  dynamics.Rate, dynamics.Source[nelem][ion]);*/
951  }
952  lgHomogeneous = false;
953  }
954 
955  /* >>chng 06 nov 21, for very high ionization parameter sims the H molecular
956  * fraction can become so small that atom = atom + molecule. In this case we
957  * will not count system as an inhomogeneous case since linear algebra package
958  * will fail. If molecules are very small, count as homogeneous.
959  * Note that we have already added sink terms to the main matrix and they
960  * will not be removed, a possible source of error, but they must not
961  * have been significant, given that the molecular fraction is so small */
962  if( !lgHomogeneous && ion_range==2 )
963  {
964  /* solve algebraically */
965  double a = MAT( xmat, 0, 0 ),
966  b = MAT( xmat, 1, 0 ) ,
967  c = MAT( xmat, 0, 1 ) ,
968  d = MAT( xmat, 1, 1 );
969  b = SDIV(b);
970  d = SDIV(d);
971  double ratio1 = a/b , ratio2 = c/d , fratio1=fabs(a/b),fratio2=fabs(c/d);
972  if( fabs(ratio1-ratio2) <= DBL_EPSILON*1e4*MAX2(fratio1,fratio2) )
973  {
974  lgHomogeneous = true;
975  }
976  }
977 
978 #if 1
979  if(
980  fabs(dense.xMolecules(nelem)) <DBL_EPSILON*100.*dense.gas_phase[nelem] )
981  {
982  lgHomogeneous = true;
983  }
984 #endif
985 
986  /* this is true if no source terms
987  * we will use total population and species conservation to replace one
988  * set of balance equations since overdetermined */
989  if( 1 || lgHomogeneous )
990  {
991  double rate_ioniz=1., rate_recomb /*, scale = 0.*/;
992  /* Simple estimate of most abundant ion */
993  long jmax = 0;
994  for( long i=0; i<ion_range-1;i++)
995  {
996  long ion = i+ion_low;
997  rate_ioniz *= ionbal.RateIonizTot(nelem,ion);
998  rate_recomb = ionbal.RateRecomTot[nelem][ion];
999  /* trips if ion rate zero, so all the gas will be more neutral than this */
1000  if( rate_ioniz == 0)
1001  break;
1002  /* rec rate is zero */
1003  if( rate_recomb <= 0.)
1004  break;
1005 
1006  rate_ioniz /= rate_recomb;
1007  if( rate_ioniz > 1.)
1008  {
1009  /* this is peak ionization stage */
1010  jmax = i;
1011  rate_ioniz = 1.;
1012  }
1013  }
1014 
1015  /* replace its matrix elements with population sum */
1016  for( long i=0; i<ion_range;i++ )
1017  {
1018  MAT(xmat,i,jmax) -= conserve_scale;
1019  }
1020  source[jmax] -= abund_total*conserve_scale;
1021  }
1022 
1023 #if 0
1024  if( false && nelem == ipHYDROGEN && dynamics.lgAdvection&& iteration>1 )
1025  {
1026  fprintf(ioQQQ,"DEBUGG Rate %.2f %.3e \n",fnzone,dynamics.Rate);
1027  fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateIonizTot(nelem,0), ionbal.RateIonizTot(nelem,1) );
1028  fprintf(ioQQQ," %.3e %.3e\n", ionbal.RateRecomTot[nelem][0], ionbal.RateRecomTot[nelem][1]);
1029  fprintf(ioQQQ," %.3e %.3e %.3e\n\n", dynamics.Source[nelem][0], dynamics.Source[nelem][1], dynamics.Source[nelem][2]);
1030  }
1031 
1032  {
1033  /* option to print matrix */
1034  enum {DEBUG_LOC=false};
1035  if( DEBUG_LOC && nzone==1 && lgPrintIt )
1036  {
1037  fprintf( ioQQQ,
1038  " DEBUG ion_solver: nelem=%li ion_range=%li, limit=%li, nConv %li xmat follows\n",
1039  nelem , ion_range,limit , conv.nTotalIoniz );
1040  if( lgHomogeneous )
1041  fprintf(ioQQQ , "Homogeneous \n");
1042  for( long i=0; i<ion_range; ++i )
1043  {
1044  for( long j=0;j<ion_range;j++ )
1045  {
1046  fprintf(ioQQQ,"%e\t",MAT(xmat,i,j));
1047  }
1048  fprintf(ioQQQ,"\n");
1049  }
1050  fprintf(ioQQQ,"source follows\n");
1051  for( long i=0; i<ion_range;i++ )
1052  {
1053  fprintf(ioQQQ,"%e\t",source[i]);
1054  }
1055  fprintf(ioQQQ,"\n");
1056  fprintf(ioQQQ,"totsrc/totscl %g %g\n",totsrc,totscl);
1057  }
1058  }
1059 #endif
1060 
1061 }
1062 
1063 STATIC void PrintRates( long nelem, bool lgNegPop, double abund_total, valarray<double> &auger, bool lgPrintIt )
1064 {
1065  DEBUG_ENTRY( "PrintRates()" );
1066 
1067  long ion;
1068 
1069  /* this should not happen */
1070  if( lgNegPop )
1071  {
1072  fprintf( ioQQQ, " PROBLEM Negative population found for abundance of ionization stage of element %4.4s, ZONE=%4ld\n",
1073  elementnames.chElementNameShort[nelem], nzone );
1074 
1075  fprintf( ioQQQ, " Populations were" );
1076  for( ion=1; ion <= dense.IonHigh[nelem]+1; ion++ )
1077  {
1078  fprintf( ioQQQ, "%9.1e", dense.xIonDense[nelem][ion-1] );
1079  }
1080  fprintf( ioQQQ, "\n" );
1081 
1082  fprintf( ioQQQ, " destroy vector =" );
1083  for( ion=1; ion <= dense.IonHigh[nelem]; ion++ )
1084  {
1085  fprintf( ioQQQ, "%9.1e", ionbal.RateIonizTot(nelem,ion-1) );
1086  }
1087  fprintf( ioQQQ, "\n" );
1088 
1089  fprintf( ioQQQ, " CTHeavy vector =" );
1090  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1091  {
1092  fprintf( ioQQQ, "%9.1e", atmdat.CharExcIonOf[ipHELIUM][nelem][ion] );
1093  }
1094  fprintf( ioQQQ, "\n" );
1095 
1096  fprintf( ioQQQ, " CharExcIonOf[ipHYDROGEN] vtr=" );
1097  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1098  {
1099  fprintf( ioQQQ, "%9.1e", atmdat.CharExcIonOf[ipHYDROGEN][nelem][ion] );
1100  }
1101  fprintf( ioQQQ, "\n" );
1102 
1103  fprintf( ioQQQ, " CollidRate vtr=" );
1104  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1105  {
1106  fprintf( ioQQQ, "%9.1e", ionbal.CollIonRate_Ground[nelem][ion][0] );
1107  }
1108  fprintf( ioQQQ, "\n" );
1109 
1110  /* photo rates per subshell */
1111  fprintf( ioQQQ, " photo rates per subshell, ion\n" );
1112  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1113  {
1114  fprintf( ioQQQ, "%3ld", ion );
1115  for( long ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
1116  {
1117  fprintf( ioQQQ, "%9.1e", ionbal.PhotoRate_Shell[nelem][ion][ns][0] );
1118  }
1119  fprintf( ioQQQ, "\n" );
1120  }
1121 
1122  /* now check out creation vector */
1123  fprintf( ioQQQ, " create vector =" );
1124  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1125  {
1126  fprintf( ioQQQ, "%9.1e", ionbal.RateRecomTot[nelem][ion] );
1127  }
1128  fprintf( ioQQQ, "\n" );
1129  }
1130 
1131  /* option to print ionization and recombination arrays
1132  * prt flag set with print array print arrays command */
1133  if( lgPrintIt || prt.lgPrtArry[nelem] || lgNegPop )
1134  {
1135  const char* format = " %9.2e";
1136  /* say who we are, what we are doing .... */
1137  fprintf( ioQQQ,
1138  "\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",
1139  elementnames.chElementSym[nelem],
1140  elementnames.chElementName[nelem],
1141  fnzone,
1142  phycon.te ,
1143  dense.eden,
1144  dense.gas_phase[nelem],
1145  abund_total ,
1146  dense.xMolecules(nelem) );
1147  /* total ionization rate, all processes */
1148  fprintf( ioQQQ, " %s Ioniz total " ,elementnames.chElementSym[nelem]);
1149  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1150  {
1151  fprintf( ioQQQ, format, ionbal.RateIonizTot(nelem,ion) );
1152  }
1153  fprintf( ioQQQ, "\n" );
1154 
1155  /* sum of all creation processes */
1156  fprintf( ioQQQ, " %s Recom total " ,elementnames.chElementSym[nelem]);
1157  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1158  {
1159  fprintf( ioQQQ, format, ionbal.RateRecomTot[nelem][ion] );
1160  }
1161  fprintf( ioQQQ, "\n" );
1162 
1163  /* collisional ionization */
1164  fprintf( ioQQQ, " %s Coll ioniz " ,elementnames.chElementSym[nelem] );
1165  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1166  {
1167  double ColIoniz = ionbal.CollIonRate_Ground[nelem][ion][0];
1168  if( ion > nelem - NISO )
1169  {
1170  long ipISO = nelem-ion;
1171  ASSERT( ipISO >=0 && ipISO < NISO );
1172  if( dense.xIonDense[nelem][nelem-ipISO] > 0. )
1173  {
1174  ColIoniz *= iso_sp[ipISO][nelem].st[0].Pop();
1175  for( long ipLevel=1; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
1176  {
1177  ColIoniz += iso_sp[ipISO][nelem].fb[ipLevel].ColIoniz * dense.EdenHCorr * iso_sp[ipISO][nelem].st[ipLevel].Pop();
1178  }
1179  ColIoniz /= dense.xIonDense[nelem][nelem-ipISO];
1180  }
1181  }
1182  fprintf( ioQQQ, format, ColIoniz );
1183  }
1184  fprintf( ioQQQ, "\n" );
1185 
1186  /* UTA ionization */
1187  fprintf( ioQQQ, " %s UTA ioniz " ,elementnames.chElementSym[nelem] );
1188  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1189  {
1190  fprintf( ioQQQ, format, ionbal.UTA_ionize_rate[nelem][ion] );
1191  }
1192  fprintf( ioQQQ, "\n" );
1193 
1194  /* photo ionization */
1195  fprintf( ioQQQ, " %s Photoion snk" ,elementnames.chElementSym[nelem] );
1196  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1197  {
1198  double PhotIoniz = 0.;
1199  for( long ipShell = 0; ipShell < Heavy.nsShells[nelem][ion]; ipShell++ )
1200  PhotIoniz += ionbal.PhotoRate_Shell[nelem][ion][ipShell][0];
1201 
1202  // still don't have the total if stage is in one of the iso-sequences
1203  if( ion > nelem - NISO )
1204  {
1205  long ipISO = nelem-ion;
1206  ASSERT( ipISO >=0 && ipISO < NISO );
1207  if( dense.xIonDense[nelem][nelem-ipISO]>0 )
1208  {
1209  PhotIoniz *= iso_sp[ipISO][nelem].st[0].Pop();
1210  for( long ipLevel=1; ipLevel < iso_sp[ipISO][nelem].numLevels_local; ipLevel++ )
1211  {
1212  PhotIoniz += iso_sp[ipISO][nelem].fb[ipLevel].gamnc * iso_sp[ipISO][nelem].st[ipLevel].Pop();
1213  }
1214  PhotIoniz /= dense.xIonDense[nelem][nelem-ipISO];
1215  }
1216  }
1217  fprintf( ioQQQ, format, PhotIoniz );
1218  }
1219  fprintf( ioQQQ, "\n" );
1220 
1221  /* photoionization source (source of this stage due to ionization of next stage) */
1222  fprintf( ioQQQ, " %s Photoion src" ,elementnames.chElementSym[nelem]);
1223  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1224  {
1225  double source = 0.;
1226  if( ion>0 )
1227  source = ionbal.RateIoniz[nelem][ion-1][ion];
1228 
1229  fprintf( ioQQQ, format, source );
1230  }
1231  fprintf( ioQQQ, "\n" );
1232 
1233  /* auger production (of this stage due to auger ionization of another) */
1234  fprintf( ioQQQ, " %s Auger src " ,elementnames.chElementSym[nelem]);
1235  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1236  {
1237  double source = 0.;
1238  if( ion>0 )
1239  source = auger[ion-1];
1240 
1241  fprintf( ioQQQ, format, source );
1242  }
1243  fprintf( ioQQQ, "\n" );
1244 
1245  /* secondary ionization */
1246  fprintf( ioQQQ, " %s Secon ioniz " ,elementnames.chElementSym[nelem]);
1247  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1248  {
1249  fprintf( ioQQQ, format,
1250  secondaries.csupra[nelem][ion] );
1251  }
1252  fprintf( ioQQQ, "\n" );
1253 
1254  /* grain ionization - not total rate but should be dominant process */
1255  fprintf( ioQQQ, " %s ion on grn " ,elementnames.chElementSym[nelem]);
1256  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1257  {
1258  fprintf( ioQQQ, format, gv.GrainChTrRate[nelem][ion][ion+1] );
1259  }
1260  fprintf( ioQQQ, "\n" );
1261 
1262  /* charge transfer ionization from chemistry */
1263  fprintf( ioQQQ, " %s ion xfr mol " ,elementnames.chElementSym[nelem]);
1264  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1265  {
1266  fprintf( ioQQQ, format, mole.xMoleChTrRate[nelem][ion][ion+1] );
1267  }
1268  fprintf( ioQQQ, "\n" );
1269 
1270  /* charge exchange ionization */
1271  fprintf( ioQQQ, " %s chr trn ion " ,elementnames.chElementSym[nelem] );
1272  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1273  {
1274  /* sum has units s-1 */
1275  double sum;
1276  long ipISO=nelem-ion;
1277 
1278  if( nelem < t_atmdat::NCX && nelem == ipISO )
1279  {
1280  sum = atmdat.CharExcIonTotal[nelem] * iso_sp[ipISO][nelem].st[0].Pop() / SDIV(dense.xIonDense[nelem][nelem-ipISO]);
1281  }
1282  else
1283  {
1284  sum = 0.0;
1285  for (long nelem1=0; nelem1 < t_atmdat::NCX; ++nelem1)
1286  sum += atmdat.CharExcIonOf[nelem1][nelem][ion] * dense.xIonDense[nelem1][1];
1287  }
1288 
1289  fprintf( ioQQQ, format, sum );
1290  }
1291  fprintf( ioQQQ, "\n" );
1292 
1293  /* radiative recombination */
1294  fprintf( ioQQQ, " %s radiati rec " ,elementnames.chElementSym[nelem]);
1295  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1296  {
1297  fprintf( ioQQQ, format, dense.eden*ionbal.RR_rate_coef_used[nelem][ion] );
1298  }
1299  fprintf( ioQQQ, "\n" );
1300 
1301  /* Badnell DR recombination */
1302  fprintf( ioQQQ, " %s drBadnel rec" ,elementnames.chElementSym[nelem]);
1303  for( ion=0; ion < min(nelem-1,dense.IonHigh[nelem]); ion++ )
1304  {
1305  fprintf( ioQQQ, format, dense.eden*ionbal.DR_Badnell_rate_coef[nelem][ion] );
1306  }
1307  fprintf( ioQQQ, "\n" );
1308 
1309  /* Cota rate */
1310  fprintf( ioQQQ, " %s CotaRate rec" ,elementnames.chElementSym[nelem]);
1311  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1312  {
1313  fprintf( ioQQQ, format, dense.eden*ionbal.CotaRate[ion] );
1314  }
1315  fprintf( ioQQQ, "\n" );
1316 
1317  /* grain recombination - not total but from next higher ion, should
1318  * be dominant */
1319  fprintf( ioQQQ, " %s rec on grn " ,elementnames.chElementSym[nelem]);
1320  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1321  {
1322  fprintf( ioQQQ, format, gv.GrainChTrRate[nelem][ion+1][ion] );
1323  }
1324  fprintf( ioQQQ, "\n" );
1325 
1326  /* charge transfer recombination from chemistry */
1327  fprintf( ioQQQ, " %s rec xfr mol " ,elementnames.chElementSym[nelem]);
1328  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1329  {
1330  fprintf( ioQQQ, format, mole.xMoleChTrRate[nelem][ion+1][ion] );
1331  }
1332  fprintf( ioQQQ, "\n" );
1333 
1334  /* charge exchange recombination */
1335  fprintf( ioQQQ, " %s chr trn rec " ,elementnames.chElementSym[nelem]);
1336  for( ion=0; ion < dense.IonHigh[nelem]; ion++ )
1337  {
1338  double sum;
1339 
1340  if( nelem==ipHELIUM && ion==0 )
1341  {
1342  sum = atmdat.CharExcRecTotal[nelem];
1343  }
1344  else if( nelem==ipHYDROGEN && ion==0 )
1345  {
1346  sum = atmdat.CharExcRecTotal[nelem];
1347  }
1348  else
1349  {
1350  sum = 0.0;
1351  for (long nelem1=0; nelem1<t_atmdat::NCX; ++nelem1)
1352  {
1353  long ipISO=nelem1;
1354  sum += atmdat.CharExcRecTo[nelem1][nelem][ion] * iso_sp[ipISO][nelem1].st[0].Pop();
1355  }
1356  }
1357  fprintf( ioQQQ, format, sum );
1358  }
1359  fprintf( ioQQQ, "\n" );
1360 
1361  {
1362  /* sources from the chemistry network */
1363  fprintf(ioQQQ," %s molecule src",
1364  elementnames.chElementSym[nelem]);
1365  for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
1366  {
1367  fprintf(ioQQQ,format, mole.source[nelem][ion]/SDIV(dense.xIonDense[nelem][ion]) );
1368  }
1369  fprintf( ioQQQ, "\n" );
1370 
1371  /* sinks from the chemistry network */
1372  fprintf(ioQQQ," %s molecule snk",
1373  elementnames.chElementSym[nelem]);
1374  for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
1375  {
1376  fprintf(ioQQQ,format, mole.sink[nelem][ion] );
1377  }
1378  fprintf( ioQQQ, "\n" );
1379 
1380  if( dynamics.lgAdvection )
1381  {
1382  /* source from the dynamcs */
1383  fprintf(ioQQQ," %s dynamics src",
1384  elementnames.chElementSym[nelem]);
1385  for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
1386  {
1387  fprintf(ioQQQ,format, dynamics.Source[nelem][ion]/SDIV(dense.xIonDense[nelem][ion]) );
1388  }
1389  fprintf( ioQQQ, "\n" );
1390 
1391  /* sinks from the dynamics */
1392  fprintf(ioQQQ," %s dynamics snk",
1393  elementnames.chElementSym[nelem]);
1394  for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
1395  {
1396  fprintf(ioQQQ,format, dynamics.Rate );
1397  }
1398  fprintf( ioQQQ, "\n" );
1399  }
1400  }
1401 
1402  /* the "new" abundances the resulted from the previous ratio */
1403  fprintf( ioQQQ, " %s Abun [cm-3] " ,elementnames.chElementSym[nelem] );
1404  for( ion=0; ion <= dense.IonHigh[nelem]; ion++ )
1405  {
1406  fprintf( ioQQQ, format, dense.xIonDense[nelem][ion] );
1407  }
1408  fprintf( ioQQQ, "\n" );
1409  }
1410 
1411  if( lgNegPop )
1412  {
1413  ContNegative();
1414  ShowMe();
1415  cdEXIT(EXIT_FAILURE);
1416  }
1417 
1418  return;
1419 }
1420 
1421 /*
1422 
1423  Solve an ionization level system with specified ionization and
1424  recombination rates between neighboring ions, and additional sink
1425  and source terms. The sink array is overwritten, and the results
1426  appear in the source array.
1427 
1428  Written in matrix form, the algorithm is equivalent to the
1429  tridiagonal algorithm in Numerical Recipes applied to:
1430 
1431  / i_0+a_0 -r_0 . . . \ / x_0 \ / s_0 \
1432  | -i_0 i_1+a_1+r_0 -r_1 . . | | x_1 | | s_1 |
1433  | . -i_1 i_2+a_2+r_1 -r_2 . | | x_2 | | s_2 |
1434  | . . (etc....) | | ... | = | ... |
1435  \ . . . / \ / \ /
1436 
1437  where i, r are the ionization and recombination rates, s is the
1438  source rate and a is the sink rate.
1439 
1440  This matrix is diagonally dominant only when the sink terms are
1441  large -- the alternative method coded here prevents rounding error
1442  in the diagonal terms disturbing the solution when this is not the
1443  case.
1444 
1445 */
1446 
1447 /* solveions tridiagonal solver but optimized for structure of balance matrix */
1448 void solveions(double *ion, double *rec, double *snk, double *src,
1449  long int nlev, long int nmax)
1450 {
1451  double kap, bet;
1452  long int i;
1453 
1454  DEBUG_ENTRY( "solveions()" );
1455 
1456  if(nmax != -1)
1457  {
1458  /* Singular case */
1459  src[nmax] = 1.;
1460  for(i=nmax;i<nlev-1;i++)
1461  src[i+1] = src[i]*ion[i]/rec[i];
1462  for(i=nmax-1;i>=0;i--)
1463  src[i] = src[i+1]*rec[i]/ion[i];
1464  }
1465  else
1466  {
1467  i = 0;
1468  kap = snk[0];
1469  for(i=0;i<nlev-1;i++)
1470  {
1471  bet = ion[i]+kap;
1472  if(bet == 0.)
1473  {
1474  fprintf(ioQQQ,"Ionization solver error\n");
1475  cdEXIT(EXIT_FAILURE);
1476  }
1477  bet = 1./bet;
1478  src[i] *= bet;
1479  src[i+1] += ion[i]*src[i];
1480  snk[i] = bet*rec[i];
1481  kap = kap*snk[i]+snk[i+1];
1482  }
1483  bet = kap;
1484  if(bet == 0.)
1485  {
1486  fprintf(ioQQQ,"Ionization solver error\n");
1487  cdEXIT(EXIT_FAILURE);
1488  }
1489  src[i] /= bet;
1490 
1491  for(i=nlev-2;i>=0;i--)
1492  {
1493  src[i] += snk[i]*src[i+1];
1494  }
1495  }
1496 }
1497 
1498 void ion_wrapper( long nelem )
1499 {
1500 
1501  DEBUG_ENTRY( "ion_wrapper()" );
1502 
1503  ASSERT( nelem >= ipHYDROGEN );
1504  ASSERT( nelem < LIMELM );
1505 
1506  switch( nelem )
1507  {
1508  case ipHYDROGEN:
1509  IonHydro();
1510  break;
1511  case ipHELIUM:
1512  IonHelium();
1513  break;
1514  default:
1515  IonNelem(false,nelem);
1516  break;
1517  }
1518 
1519  if( trace.lgTrace && dense.lgElmtOn[nelem] && trace.lgHeavyBug )
1520  {
1521  fprintf( ioQQQ, " ion_wrapper returns; %s fracs = ", elementnames.chElementSym[nelem] );
1522  for( long ion = 0; ion<=nelem+1; ion++ )
1523  fprintf( ioQQQ,"%10.3e ", dense.xIonDense[nelem][ion]/dense.gas_phase[nelem] );
1524  fprintf( ioQQQ, "\n" );
1525  }
1526 
1527  ASSERT(lgElemsConserved());
1528 
1529  return;
1530 }
t_elementnames::chElementName
char chElementName[LIMELM][CHARS_ELEMENT_NAME]
Definition: elementnames.h:17
thermal.h
t_atmdat::CharExcIonOf
double CharExcIonOf[NCX][LIMELM][LIMELM+1]
Definition: atmdat.h:152
t_prt::lgPrtArry
bool lgPrtArry[LIMELM]
Definition: prt.h:195
fill_array
STATIC void fill_array(long int nelem, long ion_range, valarray< double > &xmat, valarray< double > &source, valarray< double > &auger, double *abund_total)
Definition: ion_solver.cpp:648
TorF
char TorF(bool l)
Definition: cddefines.h:710
ipOXYGEN
const int ipOXYGEN
Definition: cddefines.h:312
solveions
void solveions(double *ion, double *rec, double *snk, double *src, long int nlev, long int nmax)
Definition: ion_solver.cpp:1448
iso_solve
void iso_solve(long ipISO, long nelem, double &maxerr)
Definition: iso_solve.cpp:102
yield.h
t_ionbal::PhotoRate_Shell
double **** PhotoRate_Shell
Definition: ionbal.h:111
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:63
IonHelium
void IonHelium(void)
Definition: ion_helium.cpp:12
t_dense::eden
double eden
Definition: dense.h:190
t_dense::EdenHCorr
double EdenHCorr
Definition: dense.h:216
dense
t_dense dense
Definition: dense.cpp:24
elementnames.h
secondaries.h
t_dynamics::lgAdvection
bool lgAdvection
Definition: dynamics.h:60
ion_wrapper
void ion_wrapper(long nelem)
Definition: ion_solver.cpp:1498
Singleton< t_yield >::Inst
static t_yield & Inst()
Definition: cddefines.h:175
t_ionbal::CotaRate
realnum CotaRate[LIMELM]
Definition: ionbal.h:242
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
newton_step.h
t_atmdat::CharExcIonTotal
double CharExcIonTotal[NCX]
Definition: atmdat.h:162
realnum
float realnum
Definition: cddefines.h:103
conv.h
STATIC
#define STATIC
Definition: cddefines.h:97
t_ionbal::RateRecomTot
double ** RateRecomTot
Definition: ionbal.h:198
ipCARBON
const int ipCARBON
Definition: cddefines.h:310
mole.h
get_total_abundance_ions
STATIC double get_total_abundance_ions(long int nelem)
Definition: ion_solver.cpp:617
PrintRates
STATIC void PrintRates(long nelem, bool lgNegPop, double abund_total, valarray< double > &auger, bool lgPrintIt)
Definition: ion_solver.cpp:1063
t_trace::lgHeavyBug
bool lgHeavyBug
Definition: trace.h:21
t_dynamics::Rate
double Rate
Definition: dynamics.h:71
t_Heavy::nsShells
long int nsShells[LIMELM][LIMELM]
Definition: heavy.h:28
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
find_solution
STATIC void find_solution(long nelem, long ion_range, valarray< double > &xmat, valarray< double > &source)
Definition: ion_solver.cpp:264
t_ionbal::RR_rate_coef_used
double ** RR_rate_coef_used
Definition: ionbal.h:212
thirdparty.h
t_iso_sp::numLevels_local
long int numLevels_local
Definition: iso.h:498
phycon
t_phycon phycon
Definition: phycon.cpp:6
trace.h
t_conv::GasPhaseAbundErrorAllowed
realnum GasPhaseAbundErrorAllowed
Definition: conv.h:285
GrainVar::lgGrainPhysicsOn
bool lgGrainPhysicsOn
Definition: grainvar.h:479
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:952
t_dynamics::Source
double ** Source
Definition: dynamics.h:74
t_dense::gas_phase
realnum gas_phase[LIMELM]
Definition: dense.h:71
getrs_wrapper
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
Definition: thirdparty_lapack.cpp:69
t_dynamics::n_initial_relax
long int n_initial_relax
Definition: dynamics.h:126
Heavy
t_Heavy Heavy
Definition: heavy.cpp:5
dynamics.h
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
iso_departure_coefficients
void iso_departure_coefficients(long ipISO, long nelem)
Definition: iso_solve.cpp:363
HomogeneousSource
STATIC void HomogeneousSource(long nelem, long ion_low, long ion_range, valarray< double > &xmat, valarray< double > &source, double abund_total)
Definition: ion_solver.cpp:901
t_ionbal::UTA_ionize_rate
double ** UTA_ionize_rate
Definition: ionbal.h:170
MAT1
#define MAT1(M_, I_, J_)
Definition: ion_solver.cpp:57
iso.h
ISO_LOOPS
@ ISO_LOOPS
Definition: conv.h:79
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
atmdat.h
t_ionbal::DR_Badnell_rate_coef
double ** DR_Badnell_rate_coef
Definition: ionbal.h:205
MIN2
#define MIN2
Definition: cddefines.h:761
lgOH_ChargeTransferDominant
bool lgOH_ChargeTransferDominant(void)
Definition: ion_solver.cpp:370
t_mole_local::sink
double ** sink
Definition: mole.h:394
nzone
long int nzone
Definition: cddefines.cpp:14
broken
void broken(void)
Definition: service.cpp:982
MAT2
#define MAT2(M_, I_, J_)
Definition: ion_solver.cpp:60
source
static double * source
Definition: species2.cpp:28
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
ContNegative
void ContNegative(void)
Definition: cont_negative.cpp:10
dense.h
t_ionbal::elecsrc
double elecsrc[LIMELM]
Definition: ionbal.h:251
mole
t_mole_local mole
Definition: mole.cpp:7
trace
t_trace trace
Definition: trace.cpp:5
prt
t_prt prt
Definition: prt.cpp:10
cddefines.h
GrainVar::GrainChTrRate
realnum GrainChTrRate[LIMELM][LIMELM+1][LIMELM+1]
Definition: grainvar.h:541
ion_solver
void ion_solver(long int nelem, bool lgPrintIt)
Definition: ion_solver.cpp:62
t_mole_local::xMoleChTrRate
realnum *** xMoleChTrRate
Definition: mole.h:396
t_iso_sp::numLevels_max
long int numLevels_max
Definition: iso.h:493
thermal
t_thermal thermal
Definition: thermal.cpp:5
MAX_DENSITY
const double MAX_DENSITY
Definition: cddefines.h:269
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
isnan
#define isnan
Definition: cddefines.h:620
t_atmdat::CharExcRecTotal
double CharExcRecTotal[NCX]
Definition: atmdat.h:163
ION_SOLVES
@ ION_SOLVES
Definition: conv.h:78
t_thermal::heating
double heating[LIMELM][LIMELM]
Definition: thermal.h:158
radius.h
IonHydro
void IonHydro()
Definition: iso_solve.cpp:149
t_yield::elec_eject_frac
realnum elec_eject_frac(long n, long i, long ns, long ne) const
Definition: yield.h:48
t_dense::xMolecules
realnum xMolecules(long nelem)
Definition: dense.h:83
heavy.h
hmi.h
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:452
MAX2
#define MAX2
Definition: cddefines.h:782
ionbal
t_ionbal ionbal
Definition: ionbal.cpp:5
iso_renorm
void iso_renorm(long nelem, long ipISO, double &renorm)
Definition: iso_solve.cpp:272
LIMELM
const int LIMELM
Definition: cddefines.h:258
getrf_wrapper
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
Definition: thirdparty_lapack.cpp:47
t_atmdat::NCX
@ NCX
Definition: atmdat.h:144
t_dense::IonLow
long int IonLow[LIMELM+1]
Definition: dense.h:119
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
t_iso_sp::st
qList st
Definition: iso.h:453
fp_equal_tol
bool fp_equal_tol(sys_float x, sys_float y, sys_float tol)
Definition: cddefines.h:854
t_mole_local::source
double ** source
Definition: mole.h:394
t_yield::nelec_eject
long nelec_eject(long n, long i, long ns) const
Definition: yield.h:55
t_elementnames::chElementNameShort
char chElementNameShort[LIMELM][CHARS_ELEMENT_NAME_SHORT]
Definition: elementnames.h:21
t_dense::IonHigh
long int IonHigh[LIMELM+1]
Definition: dense.h:120
iteration
long int iteration
Definition: cddefines.cpp:16
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
prt.h
grainvar.h
t_secondaries::csupra
realnum ** csupra
Definition: secondaries.h:21
t_dynamics::lgEquilibrium
bool lgEquilibrium
Definition: dynamics.h:174
t_conv::nTotalIoniz
long int nTotalIoniz
Definition: conv.h:166
ionbal.h
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:812
min
long min(int a, long b)
Definition: cddefines.h:723
secondaries
t_secondaries secondaries
Definition: secondaries.cpp:5
dynamics
t_dynamics dynamics
Definition: dynamics.cpp:44
conv
t_conv conv
Definition: conv.cpp:5
gv
GrainVar gv
Definition: grainvar.cpp:5
THRESHOLD
#define THRESHOLD
t_ionbal::RateIonizTot
double RateIonizTot(long nelem, long ion)
Definition: ionbal.h:254
IonNelem
void IonNelem(bool lgPrintIt, long int nelem)
Definition: ion_nelem.cpp:12
t_conv::incrementCounter
void incrementCounter(const counter_type type)
Definition: conv.h:308
fixit
void fixit(void)
Definition: service.cpp:991
iso_set_ion_rates
void iso_set_ion_rates(long ipISO, long nelem)
Definition: iso_level.cpp:685
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
taulines.h
t_ionbal::CollIonRate_Ground
double *** CollIonRate_Ground
Definition: ionbal.h:120
MAT
#define MAT(M_, I_, J_)
Definition: ion_solver.cpp:54
t_dense::lgSetIoniz
bool lgSetIoniz[LIMELM]
Definition: dense.h:149
lgElemsConserved
bool lgElemsConserved(void)
Definition: dense.cpp:99
phycon.h
atmdat
t_atmdat atmdat
Definition: atmdat.cpp:6
ShowMe
void ShowMe(void)
Definition: service.cpp:181
t_atmdat::lgCTOn
bool lgCTOn
Definition: atmdat.h:177
t_dense::SetIoniz
realnum SetIoniz[LIMELM][LIMELM+1]
Definition: dense.h:154
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
t_atmdat::CharExcRecTo
double CharExcRecTo[NCX][LIMELM][LIMELM+1]
Definition: atmdat.h:153
continuum.h
t_ionbal::lgGrainIonRecom
int lgGrainIonRecom
Definition: ionbal.h:228
fnzone
double fnzone
Definition: cddefines.cpp:15
store_new_densities
STATIC void store_new_densities(long nelem, long ion_range, long ion_low, double *source, double abund_total, bool *lgNegPop)
Definition: ion_solver.cpp:429
t_conv::lgSearch
bool lgSearch
Definition: conv.h:175
t_phycon::te
double te
Definition: phycon.h:11
NISO
const int NISO
Definition: cddefines.h:261
iso_satellite_update
void iso_satellite_update(long nelem)
Definition: iso_create.cpp:1381
lgTrivialSolution
STATIC bool lgTrivialSolution(long nelem, double abund_total)
Definition: ion_solver.cpp:239
GrainVar::lgDustOn
bool lgDustOn() const
Definition: grainvar.h:471
sink
static double * sink
Definition: species2.cpp:28
t_ionbal::elecsnk
double elecsnk[LIMELM]
Definition: ionbal.h:251
solve_system
int32 solve_system(const valarray< double > &a, valarray< double > &b, long int n, error_print_t error_print)
Definition: newton_step.cpp:387
t_ionbal::RateIoniz
double *** RateIoniz
Definition: ionbal.h:184
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
iso_charge_transfer_update
void iso_charge_transfer_update(long nelem)
Definition: iso_ionize_recombine.cpp:20
t_trace::lgTrace
bool lgTrace
Definition: trace.h:12
Generated by   doxygen 1.8.17