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