cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ion_recomb.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*ion_recomb generate recombination coefficients for any species */
4 /*ion_recombAGN generate recombination coefficients for AGN table */
5 #include "cddefines.h"
6 #include "phycon.h"
7 #include "heavy.h"
8 #include "grainvar.h"
9 #include "conv.h"
10 #include "thermal.h"
11 #include "iso.h"
12 #include "abund.h"
13 #include "save.h"
14 #include "elementnames.h"
15 #include "atmdat.h"
16 #include "ionbal.h"
17 #include "dense.h"
18 
19 /*ion_recomb generate recombination coefficients for any species */
21  /* this is debug flag */
22  bool lgPrintIt,
23  /* nelem is the atomic number on the C scale, 0 for H */
24  long int nelem/*,
25  double tlow[]*/)
26 {
27  long int i,
28  ion,
29  limit;
30 
31  DEBUG_ENTRY( "ion_recomb()" );
32 
33  ASSERT( nelem < LIMELM);
34  ASSERT( nelem > 1 );
35 
36  /* check that range of ionization is correct */
37  ASSERT( dense.IonLow[nelem] >= 0 );
38  ASSERT( dense.IonLow[nelem] <= nelem+1 );
39 
41 
42  /* this routine only does simple two-level species,
43  * loop over ions will be <= limit, IonHigh is -1 since very
44  * highest stage of ionization is not recombined into.
45  * for Li, will do only atom, since ions are H and He like,
46  * so limit is zero */
47 
48  /* zero-out loop comes before main loop since there are off-diagonal
49  * elements in the main ionization loop, due to multi-electron processes */
50  /* >>chng 00 dec 07, limit changed to identical to ion_solver */
51  for( ion=0; ion <= dense.IonHigh[nelem]-1; ion++ )
52  {
53  ionbal.RateRecomTot[nelem][ion] = 0.;
54  }
55 
56  for( long ipISO=ipH_LIKE; ipISO<MIN2(NISO,nelem+1); ipISO++ )
57  {
58  if( ((dense.IonHigh[nelem] >= nelem - ipISO) &&
59  (dense.IonLow[nelem] <= nelem - ipISO)) )
60  {
61  ionbal.RateRecomTot[nelem][nelem-ipISO] = ionbal.RateRecomIso[nelem][ipISO];
62  }
63  }
64 
65  limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);
66  ASSERT( limit >= -1 );
67 
68  /* these are counted elsewhere and must not be added here */
69  Heavy.xLyaHeavy[nelem][nelem] = 0.;
70  Heavy.xLyaHeavy[nelem][nelem-1] = 0.;
71 
72  /* IonLow is 0 for the atom, limit chosen to NOT include iso sequences */
73  for( ion=dense.IonLow[nelem]; ion <= limit; ion++ )
74  {
75  /* number of bound electrons of the ion after recombination,
76  * for an atom (ion=0) this is
77  * equal to nelem+1, the element on the physical scale, since nelem is
78  * on the C scale, being 5 for carbon */
79  ASSERT(ionbal.DR_Badnell_rate_coef[nelem][ion] >= 0);
80  ASSERT(ionbal.RR_rate_coef_used[nelem][ion] >= 0);
81 
82  /* sum of recombination rates [units s-1] for radiative, three body, charge transfer */
83  ionbal.RateRecomTot[nelem][ion] =
84  dense.eden* (
85  ionbal.RR_rate_coef_used[nelem][ion] +
86  ionbal.DR_Badnell_rate_coef[nelem][ion] +
87  ionbal.CotaRate[ion] );
88 
89  /* >>chng 01 jun 30, FRAC_LINE was 0.1, not 1, did not include anything except
90  * radiative recombination, the radrec term */
91  static const double FRAC_LINE = 1.;
92  /* was 0.1 */
93  /*Heavy.xLyaHeavy[nelem][ion] = (realnum)(dense.eden*radrec*FRAC_LINE );*/
94  Heavy.xLyaHeavy[nelem][ion] = (realnum)(dense.eden*
95  (ionbal.RR_rate_coef_used[nelem][ion]+ionbal.DR_Badnell_rate_coef[nelem][ion])*FRAC_LINE );
96  }
97 
98  /* option to save rec coefficients */
99  if( save.lgioRecom || lgPrintIt )
100  {
101  /* >>chng 04 feb 22, make option to print ions for single element */
102  FILE *ioOut;
103  if( lgPrintIt )
104  ioOut = ioQQQ;
105  else
106  ioOut = save.ioRecom;
107 
108  /* print name of element */
109  fprintf( ioOut,
110  " %s recombination coefficients fnzone:%.2f \tte\t%.4e\tne\t%.4e\n",
112 
113  /*limit = MIN2(11,dense.IonHigh[nelem]);*/
114  /* >>chng 05 sep 24, just print one long line - need info */
115  limit = dense.IonHigh[nelem];
116  // give ion stage
117  for( i=0; i<limit; ++i )
118  fprintf( ioOut, "%10ld",i+1);
119  fprintf( ioOut, "\n");
120 
121  for( i=0; i < limit; i++ )
122  {
123  fprintf( ioOut, "%10.2e", ionbal.RR_rate_coef_used[nelem][i] );
124  }
125  fprintf( ioOut, " radiative used vs Z\n" );
126 
127  for( i=0; i < limit; i++ )
128  {
129  fprintf( ioOut, "%10.2e", ionbal.RR_Verner_rate_coef[nelem][i] );
130  }
131  fprintf( ioOut, " old Verner vs Z\n" );
132 
133  for( i=0; i < limit; i++ )
134  {
135  fprintf( ioOut, "%10.2e", ionbal.RR_Badnell_rate_coef[nelem][i] );
136  }
137  fprintf( ioOut, " new Badnell vs Z\n" );
138 
139  for( i=0; i < limit; i++ )
140  {
141  /* >>chng 06 jan 19, from div by eden to div by H0 - want units of cm3 s-1 but
142  * no single collider does this so not possible to get rate coefficient easily
143  * H0 is more appropriate than electron density */
144  fprintf( ioOut, "%10.2e", ionbal.CX_recomb_rate_used[nelem][i]/SDIV(dense.xIonDense[ipHYDROGEN][0]) );
145  }
146  fprintf( ioOut, " CT/n(H0)\n" );
147 
148  for( i=0; i < limit; i++ )
149  {
150  fprintf( ioOut, "%10.2e", ionbal.CotaRate[ion] );
151  }
152  fprintf( ioOut, " 3body vs Z /ne\n" );
153 
154  /* note different upper limit - this routine does grain rec for all ions */
155  for( i=0; i < dense.IonHigh[nelem]; i++ )
156  {
157  fprintf( ioOut, "%10.2e", gv.GrainChTrRate[nelem][i+1][i]/dense.eden );
158  }
159  fprintf( ioOut, " Grain vs Z /ne\n" );
160  fprintf( ioOut, " old Nussbaumer Storey DR vs Z\n" );
161 
162  for( i=0; i < limit; i++ )
163  {
164  fprintf( ioOut, "%10.2e", ionbal.DR_Badnell_rate_coef[nelem][i] );
165  }
166  fprintf( ioOut, " new Badnell DR vs Z\n" );
167 
168  /* total recombination rate, with density included - this goes into the matrix */
169  for( i=0; i < limit; i++ )
170  {
171  fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i] );
172  }
173  fprintf( ioOut,
174  " total rec rate (with density) for %s\n",
175  elementnames.chElementSym[nelem] );
176  for( i=0; i < limit; i++ )
177  {
178  fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i]/dense.eden );
179  }
180  fprintf( ioOut,
181  " total rec rate / ne for %s\n\n",
182  elementnames.chElementSym[nelem] );
183 
184  /* spill over to next line for many stages of ionization */
185  if( dense.IonHigh[nelem] > 11 )
186  {
187  limit = MIN2(29,dense.IonHigh[nelem]);
188  fprintf( ioOut, " R " );
189  for( i=11; i < limit; i++ )
190  {
191  fprintf( ioOut, "%10.2e", dense.eden*ionbal.CotaRate[ion] );
192  }
193  fprintf( ioOut, "\n" );
194 
195  fprintf( ioOut, " " );
196  for( i=11; i < limit; i++ )
197  {
198  fprintf( ioOut, "%10.2e", ionbal.RateRecomTot[nelem][i] );
199  }
200  fprintf( ioOut, "\n\n" );
201  }
202  }
203 
204  /* >>chng 02 nov 09, from -2 to -NISO */
205  /*limit = MIN2(nelem-2,dense.IonHigh[nelem]-1);*/
206  limit = MIN2(nelem-NISO,dense.IonHigh[nelem]-1);
207  for( i=dense.IonLow[nelem]; i <= limit; i++ )
208  {
209  ASSERT( Heavy.xLyaHeavy[nelem][i] > 0. );
210  ASSERT( ionbal.RateRecomTot[nelem][i] > 0. );
211  }
212  return;
213 }
214 
215 /*ion_recombAGN generate recombination coefficients for AGN table */
216 void ion_recombAGN( FILE * io )
217 {
218  static const int N1LIM = 3;
219  static const int N2LIM = 4;
220  double te1[N1LIM]={ 5000., 10000., 20000.};
221  double te2[N2LIM]={ 20000.,50000.,100000.,1e6};
222  /* this is boundary between two tables */
223  double BreakEnergy = 100./13.0;
224  long int nelem, ion , i;
225  /* this will hold element symbol + ionization */
226  char chString[100],
227  chOutput[100];
228  /* save temp here */
229  double TempSave = phycon.te;
230  /* save ne here */
231  double EdenSave = dense.eden;
232 
233  DEBUG_ENTRY( "ion_recombAGN()" );
234 
235  EdenChange( 1. );
236  /*atmdat_readin();*/
237 
238  /* first put header on file */
239  fprintf(io,"X+i\\Te");
240  for( i=0; i<N1LIM; ++i )
241  {
242  phycon.te = te1[i];
243  fprintf(io,"\t%.0f K",phycon.te);
244  }
245  fprintf(io,"\n");
246 
247  /* now do loop over temp, but add elements */
248  for( nelem=ipLITHIUM; nelem<LIMELM; ++nelem )
249  {
250  /* this list of elements included in the AGN tables is defined in zeroabun.c */
251  if( abund.lgAGN[nelem] )
252  {
253  for( ion=0; ion<=nelem; ++ion )
254  {
255  ASSERT( Heavy.Valence_IP_Ryd[nelem][ion] > 0.05 );
256 
257  if( Heavy.Valence_IP_Ryd[nelem][ion] > BreakEnergy )
258  break;
259 
260  /* print chemical symbol */
261  sprintf(chOutput,"%s",
262  elementnames.chElementSym[nelem]);
263  /* some elements have only one letter - this avoids leaving a space */
264  if( chOutput[1]==' ' )
265  chOutput[1] = chOutput[2];
266  /* now ionization stage */
267  if( ion==0 )
268  {
269  sprintf(chString,"0 ");
270  }
271  else if( ion==1 )
272  {
273  sprintf(chString,"+ ");
274  }
275  else
276  {
277  sprintf(chString,"+%li ",ion);
278  }
279  strcat( chOutput , chString );
280  fprintf(io,"%5s",chOutput );
281 
282  for( i=0; i<N1LIM; ++i )
283  {
284  TempChange(te1[i] , false);
285  dense.IonLow[nelem] = 0;
286  dense.IonHigh[nelem] = nelem+1;
287  if( ConvBase(0) )
288  fprintf(ioQQQ,"PROBLEM ConvBase returned error.\n");
289  fprintf(io,"\t%.2e",ionbal.RateRecomTot[nelem][ion]);
290  }
291  fprintf(io,"\n");
292  }
293  fprintf(io,"\n");
294  }
295  }
296 
297  /* second put header on file */
298  fprintf(io,"X+i\\Te");
299  for( i=0; i<N2LIM; ++i )
300  {
301  TempChange(te2[i] , false);
302  fprintf(io,"\t%.0f K",phycon.te);
303  }
304  fprintf(io,"\n");
305 
306  /* now do same loop over temp, but add elements */
307  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
308  {
309  /* this list of elements included in the AGN tables is defined in zeroabun.c */
310  if( abund.lgAGN[nelem] )
311  {
312  for( ion=0; ion<=nelem; ++ion )
313  {
314  ASSERT( Heavy.Valence_IP_Ryd[nelem][ion] > 0.05 );
315 
316  if( Heavy.Valence_IP_Ryd[nelem][ion] <= BreakEnergy )
317  continue;
318 
319  /* print chemical symbol */
320  fprintf(io,"%s",
321  elementnames.chElementSym[nelem]);
322  /* now ionization stage */
323  if( ion==0 )
324  {
325  fprintf(io,"0 ");
326  }
327  else if( ion==1 )
328  {
329  fprintf(io,"+ ");
330  }
331  else
332  {
333  fprintf(io,"+%li",ion);
334  }
335 
336  for( i=0; i<N2LIM; ++i )
337  {
338  TempChange(te2[i] , false);
339  dense.IonLow[nelem] = 0;
340  dense.IonHigh[nelem] = nelem+1;
341  if( ConvBase(0) )
342  fprintf(ioQQQ,"PROBLEM ConvBase returned error.\n");
343  fprintf(io,"\t%.2e",ionbal.RateRecomTot[nelem][ion]);
344  }
345  fprintf(io,"\n");
346  }
347  fprintf(io,"\n");
348  }
349  }
350 
351  TempChange(TempSave , true);
352  EdenChange( EdenSave );
353  return;
354 }
bool lgAGN[LIMELM]
Definition: abund.h:198
double ** DR_Badnell_rate_coef
Definition: ionbal.h:197
double ** RR_Badnell_rate_coef
Definition: ionbal.h:197
t_atmdat atmdat
Definition: atmdat.cpp:6
t_Heavy Heavy
Definition: heavy.cpp:5
int ConvBase(long loopi)
Definition: conv_base.cpp:188
double ** RR_Verner_rate_coef
Definition: ionbal.h:210
const int NISO
Definition: cddefines.h:310
long int IonHigh[LIMELM+1]
Definition: dense.h:130
realnum xLyaHeavy[LIMELM][LIMELM]
Definition: heavy.h:21
t_phycon phycon
Definition: phycon.cpp:6
void ion_recombAGN(FILE *io)
Definition: ion_recomb.cpp:216
FILE * ioRecom
Definition: save.h:457
FILE * ioQQQ
Definition: cddefines.cpp:7
long int nsbig
Definition: atmdat.h:345
#define MIN2(a, b)
Definition: cddefines.h:807
void TempChange(double TempNew, bool lgForceUpdate)
Definition: temp_change.cpp:31
t_dense dense
Definition: global.cpp:15
double ** RateRecomIso
Definition: ionbal.h:194
t_elementnames elementnames
Definition: elementnames.cpp:5
double ** RateRecomTot
Definition: ionbal.h:191
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
t_ionbal ionbal
Definition: ionbal.cpp:8
t_abund abund
Definition: abund.cpp:5
long int IonLow[LIMELM+1]
Definition: dense.h:129
float realnum
Definition: cddefines.h:124
void ion_recomb(bool lgPrintIt, long int nelem)
Definition: ion_recomb.cpp:20
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
#define ASSERT(exp)
Definition: cddefines.h:617
double ** CX_recomb_rate_used
Definition: ionbal.h:197
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:307
double Valence_IP_Ryd[LIMELM][LIMELM]
Definition: heavy.h:24
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
double eden
Definition: dense.h:201
#define MAX2(a, b)
Definition: cddefines.h:828
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
bool lgioRecom
Definition: save.h:458
sys_float SDIV(sys_float x)
Definition: cddefines.h:1006
char chElementName[LIMELM][CHARS_ELEMENT_NAME]
Definition: elementnames.h:17
void EdenChange(double EdenNew)
Definition: eden_change.cpp:11
GrainVar gv
Definition: grainvar.cpp:5
double fnzone
Definition: cddefines.cpp:15
t_save save
Definition: save.cpp:5
double te
Definition: phycon.h:21
const int ipHYDROGEN
Definition: cddefines.h:348
realnum GrainChTrRate[LIMELM][LIMELM+1][LIMELM+1]
Definition: grainvar.h:543
const int ipLITHIUM
Definition: cddefines.h:350
realnum CotaRate[LIMELM]
Definition: ionbal.h:236
double ** RR_rate_coef_used
Definition: ionbal.h:207