cloudy trunk
Loading...
Searching...
No Matches
species2.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#include "cddefines.h"
4#include "atmdat.h"
5#include "phycon.h"
6#include "taulines.h"
7#include "mole.h"
8#include "mole_priv.h"
9#include "atoms.h"
10#include "string.h"
11#include "thirdparty.h"
12#include "dense.h"
13#include "conv.h"
14#include "h2.h"
15#include "physconst.h"
16#include "secondaries.h"
17#include "thermal.h"
18#include "cooling.h"
19#include "lines_service.h"
20#include "atmdat.h"
21
22STATIC double LeidenCollRate(long, long, const TransitionProxy& ,double);
23STATIC double StoutCollRate(long ipSpecies, long ipCollider, const TransitionProxy&, double ftemp);
24STATIC double ChiantiCollRate(long ipSpecies, long ipCollider, const TransitionProxy&, double ftemp);
25
26static const bool DEBUGSTATE = false;
27
28static double *g, *ex, *pops, *depart, *source, *sink;
29static double **AulEscp, **col_str, **AulDest, **AulPump, **CollRate;
30
31/*Solving for the level populations*/
32
33void dBase_solve(void )
34{
36 DEBUG_ENTRY( "dBase_solve()" );
37 static bool lgFirstPass = true;
38 static long maxNumLevels = 1;
39 double totalHeating = 0.;
40
41 if( nSpecies==0 )
42 return;
43
44 if( lgFirstPass )
45 {
46 for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
47 maxNumLevels = MAX2( maxNumLevels, dBaseSpecies[ipSpecies].numLevels_max );
48
49 g = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
50 ex = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
51 pops = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
52 depart = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
53 source = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
54 sink = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
55
56 AulEscp = (double **)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double *));
57 col_str = (double **)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double *));
58 AulDest = (double **)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double *));
59 AulPump = (double **)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double *));
60 CollRate = (double **)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double *));
61
62 for( long j=0; j< maxNumLevels; j++ )
63 {
64 AulEscp[j] = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
65 col_str[j] = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
66 AulDest[j] = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
67 AulPump[j] = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
68 CollRate[j] = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
69 }
70
71 lgFirstPass = false;
72 }
73
74 // zero all of these values
75 memset( g, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
76 memset( ex, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
77 memset( pops, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
78 memset( depart, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
79 memset( source, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
80 memset( sink, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
81 for( long j=0; j< maxNumLevels; j++ )
82 {
83 memset( AulEscp[j], 0, (unsigned long)(maxNumLevels)*sizeof(double) );
84 memset( col_str[j], 0, (unsigned long)(maxNumLevels)*sizeof(double) );
85 memset( AulDest[j], 0, (unsigned long)(maxNumLevels)*sizeof(double) );
86 memset( AulPump[j], 0, (unsigned long)(maxNumLevels)*sizeof(double) );
87 memset( CollRate[j], 0, (unsigned long)(maxNumLevels)*sizeof(double) );
88 }
89
90
91 for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
92 {
93 const char *spName = dBaseSpecies[ipSpecies].chLabel;
94 double cooltl, coolder;
95 int nNegPop;
96 bool lgZeroPop, lgDeBug = false;
97
98 dBaseSpecies[ipSpecies].CoolTotal = 0.;
99
100#if 0
101 //limit for now to small number of levels
102 dBaseSpecies[ipSpecies].numLevels_local = MIN2( dBaseSpecies[ipSpecies].numLevels_local, 10 );
103#endif
104
105 /* first find current density (cm-3) of species */
106 if( dBaseSpecies[ipSpecies].lgMolecular )
107 {
108 molezone *SpeciesCurrent;
111 if( (SpeciesCurrent = findspecieslocal(dBaseSpecies[ipSpecies].chLabel)) == null_molezone )
112 {
113 /* did not find the species - print warning for now */
114 if( !conv.nTotalIoniz )
115 fprintf(ioQQQ," PROBLEM dBase_solve did not find molecular species %li\n",ipSpecies);
116 }
117 abund = (realnum)SpeciesCurrent->den;
118 }
119 else
120 {
121 /* an atom or ion */
122 ASSERT( dBaseStates[ipSpecies][0].nelem()<=LIMELM && dBaseStates[ipSpecies][0].IonStg()<=dBaseStates[ipSpecies][0].nelem()+1 );
123 abund = dense.xIonDense[ dBaseStates[ipSpecies][0].nelem()-1 ][ dBaseStates[ipSpecies][0].IonStg()-1 ];
124 }
125
126 abund *= dBaseSpecies[ipSpecies].fracType * dBaseSpecies[ipSpecies].fracIsotopologue;
127
128 // initialization at start of each iteration
129 if( conv.nTotalIoniz == 0)
130 dBaseSpecies[ipSpecies].lgActive = true;
131
132 bool lgMakeInActive = (abund <= 1e-20 * dense.xNucleiTotal);
133 if( lgMakeInActive && dBaseSpecies[ipSpecies].lgActive )
134 {
135 // zero out populations and intensities, if previously not set
136 dBaseStates[ipSpecies][0].Pop() = 0.;
137 for(long ipHi = 1; ipHi < dBaseSpecies[ipSpecies].numLevels_max; ipHi++ )
138 {
139 dBaseStates[ipSpecies][ipHi].Pop() = 0.;
140# ifdef USE_NLTE7
141 dBaseStates[ipSpecies][ipHi].DestCollBB() = 0.;
142 dBaseStates[ipSpecies][ipHi].DestPhotoBB() = 0.;
143# endif
144 }
145 for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin();
146 tr != dBaseTrans[ipSpecies].end(); ++tr)
147 {
148 (*tr).Emis().phots() = 0.;
149 (*tr).Emis().xIntensity() = 0.;
150 (*tr).Coll().col_str() = 0.;
151 (*tr).Coll().cool() = 0.;
152 (*tr).Coll().heat() = 0.;
153 }
154 dBaseSpecies[ipSpecies].lgActive = false;
155 }
156
157 if( !lgMakeInActive )
158 dBaseSpecies[ipSpecies].lgActive = true;
159
160 if( !dBaseSpecies[ipSpecies].lgActive )
161 continue;
162
163 // we always hit search phase first, reset number of levels
164 if( conv.lgSearch )
165 dBaseSpecies[ipSpecies].numLevels_local = dBaseSpecies[ipSpecies].numLevels_max;
166 for( long ipLo = 0; ipLo < dBaseSpecies[ipSpecies].numLevels_local; ipLo++ )
167 {
168 /* statistical weights & Excitation Energies*/
169 g[ipLo] = dBaseStates[ipSpecies][ipLo].g() ;
170 // parts of the code assert that ground is at zero energy - this is
171 // not true for the stored molecular data - so rescale to zero
172 ex[ipLo] = dBaseStates[ipSpecies][ipLo].energy().WN() -
173 dBaseStates[ipSpecies][0].energy().WN();
174 /* zero some rates */
175 source[ipLo] = 0.;
176 sink[ipLo] = 0.;
177 }
178
179 // non-zero was due to roundoff errors on 32-bit
180 if( ex[0] <= dBaseStates[ipSpecies][0].energy().WN()* 10. *DBL_EPSILON )
181 ex[0] = 0.;
182 else
184
185 for( long ipHi= 0; ipHi<dBaseSpecies[ipSpecies].numLevels_local; ipHi++)
186 {
187 for( long ipLo= 0; ipLo<dBaseSpecies[ipSpecies].numLevels_local; ipLo++)
188 {
189 if (ipHi > ipLo)
190 {
191 AulEscp[ipHi][ipLo] = SMALLFLOAT;
192 AulDest[ipHi][ipLo] = SMALLFLOAT;
193 AulPump[ipLo][ipHi] = SMALLFLOAT;
194 }
195 else
196 {
197 AulEscp[ipHi][ipLo] = 0.;
198 AulDest[ipHi][ipLo] = 0.;
199 AulPump[ipLo][ipHi] = 0.;
200 }
201 }
202 }
203
204 for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
205 tr != dBaseTrans[ipSpecies].end(); ++tr)
206 {
207 int ipHi = (*tr).ipHi();
208 if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local || (*tr).ipCont() <= 0)
209 continue;
210 int ipLo = (*tr).ipLo();
211 AulEscp[ipHi][ipLo] = (*tr).Emis().Aul()*
212 ((*tr).Emis().Pesc() + (*tr).Emis().Pelec_esc());
213 AulDest[ipHi][ipLo] = (*tr).Emis().Aul()*(*tr).Emis().Pdest();
214 AulPump[ipLo][ipHi] = (*tr).Emis().pump();
215 }
216
217 /*Setting all the collision strengths and collision rate to zero*/
218 for( long ipHi= 0; ipHi<dBaseSpecies[ipSpecies].numLevels_local; ipHi++)
219 {
220 for( long ipLo= 0; ipLo<dBaseSpecies[ipSpecies].numLevels_local; ipLo++)
221 {
222 col_str[ipHi][ipLo] = 0.;
223 CollRate[ipHi][ipLo] = 0.;
224 }
225 }
226
227 /*Setting all the collision strengths and collision rate to zero*/
228 for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
229 tr != dBaseTrans[ipSpecies].end(); ++tr)
230 {
231 int ipHi = (*tr).ipHi();
232 if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local)
233 continue;
234 CollisionZero( (*tr).Coll() );
235 for( long k=0; k<ipNCOLLIDER; ++k )
236 (*tr).Coll().rate_coef_ul_set()[k] = 0.f;
237 }
238
239 /* update the collision rates */
240 /* molecule */
241 if( dBaseSpecies[ipSpecies].lgLAMDA )
242 {
243 for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
244 tr != dBaseTrans[ipSpecies].end(); ++tr)
245 {
246 int ipHi = (*tr).ipHi();
247 if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local)
248 continue;
249 for( long intCollNo=0; intCollNo<ipNCOLLIDER; intCollNo++)
250 {
251 /*using the collision rate coefficients directly*/
252 (*tr).Coll().rate_coef_ul_set()[intCollNo] =
253 (realnum)LeidenCollRate(ipSpecies, intCollNo, *tr, phycon.te);
254 }
255 }
256 }
257 /* Chianti */
258 else if( dBaseTrans[ipSpecies].chLabel() == "Chianti" )
259 {
260 for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
261 tr != dBaseTrans[ipSpecies].end(); ++tr)
262 {
263 if ((*tr).ipHi() >= dBaseSpecies[ipSpecies].numLevels_local)
264 continue;
265 for( long intCollNo=0; intCollNo<ipNCOLLIDER; intCollNo++)
266 {
267 (*tr).Coll().rate_coef_ul_set()[intCollNo] =
268 (realnum)ChiantiCollRate(ipSpecies, intCollNo, *tr, phycon.te);
269 }
270 }
271 }
272 /* Stout */
273 else if( dBaseTrans[ipSpecies].chLabel() == "Stout" )
274 {
275 for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
276 tr != dBaseTrans[ipSpecies].end(); ++tr)
277 {
278 if ((*tr).ipHi() >= dBaseSpecies[ipSpecies].numLevels_local)
279 continue;
280 for( long intCollNo=0; intCollNo<ipNCOLLIDER; intCollNo++)
281 {
282 (*tr).Coll().rate_coef_ul_set()[intCollNo] =
283 (realnum)StoutCollRate(ipSpecies, intCollNo, *tr, phycon.te);
284 }
285 }
286 }
287 else
289
290 /* guess some missing data */
291 for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
292 tr != dBaseTrans[ipSpecies].end(); ++tr)
293 {
294 int ipHi = (*tr).ipHi();
295 if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local)
296 continue;
297 int ipLo = (*tr).ipLo();
298 const CollisionProxy &coll_temp = (*tr).Coll();
299
300 /* make educated guesses for some missing data */
301 if( dBaseSpecies[ipSpecies].lgMolecular )
302 {
303 ASSERT( dBaseSpecies[ipSpecies].lgLAMDA );
304 /*The collision rate coefficients for helium should not be present and that for molecular hydrogen should be present*/
305 if( AtmolCollRateCoeff[ipSpecies][ipATOM_HE].temps.size() == 0 &&
306 AtmolCollRateCoeff[ipSpecies][ipH2].temps.size() != 0 )
307 {
308 coll_temp.rate_coef_ul_set()[ipATOM_HE] = 0.7f * coll_temp.rate_coef_ul()[ipH2];
309 }
310
311 /* Put in something for hydrogen collisions if not in database */
312 if( AtmolCollRateCoeff[ipSpecies][ipATOM_H].temps.size() == 0 )
313 {
314 if( AtmolCollRateCoeff[ipSpecies][ipATOM_HE].temps.size() != 0 ) //He0
315 {
316 coll_temp.rate_coef_ul_set()[ipATOM_H] = 2.0f * coll_temp.rate_coef_ul()[ipATOM_HE];
317 }
318 else if( AtmolCollRateCoeff[ipSpecies][ipH2_ORTHO].temps.size() != 0 ) //ortho-H2
319 {
320 coll_temp.rate_coef_ul_set()[ipATOM_H] = 1.4f * coll_temp.rate_coef_ul()[ipH2_ORTHO];
321 }
322 else if( AtmolCollRateCoeff[ipSpecies][ipH2_PARA].temps.size() != 0 ) //para-H2
323 {
324 coll_temp.rate_coef_ul_set()[ipATOM_H] = 1.4f * coll_temp.rate_coef_ul()[ipH2_PARA];
325 }
326 else if( AtmolCollRateCoeff[ipSpecies][ipH2].temps.size() != 0 ) // total H2
327 {
328 coll_temp.rate_coef_ul_set()[ipATOM_H] = 1.4f * coll_temp.rate_coef_ul()[ipH2];
329 }
330 else
331 coll_temp.rate_coef_ul_set()[ipATOM_H] = 1e-13f * (realnum)g[ipLo];
332 }
333
334 /* Put in something for proton collisions if not in database */
335 if( AtmolCollRateCoeff[ipSpecies][ipPROTON].temps.size() == 0 )
336 {
337 if( AtmolCollRateCoeff[ipSpecies][ipHE_PLUS].temps.size() != 0 ) //He+
338 {
339 coll_temp.rate_coef_ul_set()[ipPROTON] = 2.0f * coll_temp.rate_coef_ul()[ipHE_PLUS];
340 }
341 else
342 coll_temp.rate_coef_ul_set()[ipPROTON] = 1e-13f * (realnum)g[ipLo];
343
344 }
345
346#if 0
347 /* if nothing else has been done, just put a small rate coefficient in */
348 for( long intCollNo=0; intCollNo<ipNCOLLIDER; intCollNo++)
349 {
350 if( coll_temp.rate_coef_ul()[intCollNo] == 0. )
351 coll_temp.rate_coef_ul_set()[intCollNo] = 1e-13;
352 }
353#endif
354 }
355 else
356 {
357 /* test for transitions without collision data */
358 if( (*tr).Coll().rate_coef_ul_set()[ipELECTRON] == 0. )
359 {
360 /* Specific transitions of Fe 3,4,and 5 have data that are not in the Chianti/Kurucz files*/
361 if( strcmp(dBaseSpecies[ipSpecies].chLabel,"Fe 3") == 0 && ipHi < 14 )
362 {
363 coll_temp.col_str() = Fe3_cs(ipLo,ipHi);
364 }
365 else if( strcmp(dBaseSpecies[ipSpecies].chLabel,"Fe 4") == 0 && ipHi < 12 )
366 {
367 coll_temp.col_str() = Fe4_cs(ipLo,ipHi);
368 }
369 else if( strcmp(dBaseSpecies[ipSpecies].chLabel,"Fe 5") == 0 && ipHi < 14 )
370 {
371 coll_temp.col_str() = Fe5_cs(ipLo,ipHi);
372 }
373 else if( atmdat.lgGbarOn )
374 {
375 /* All other transitions should use gbar if enabled */
376 MakeCS(*tr);
377 }
378 else
379 {
380 //If gbar is off, use the default collision strength value.
381 coll_temp.col_str() = atmdat.collstrDefault;
382 }
383
384 coll_temp.rate_coef_ul_set()[ipELECTRON] = (COLL_CONST*coll_temp.col_str())/
385 (g[ipHi]*phycon.sqrte);
386 }
387 }
388 }
389
390 /*Updating the CollRate*/
391 for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
392 tr != dBaseTrans[ipSpecies].end(); ++tr)
393 {
394 int ipHi = (*tr).ipHi();
395 if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local)
396 continue;
397 int ipLo = (*tr).ipLo();
398 CollRate[ipHi][ipLo] = (*tr).Coll().ColUL( colliders );
399
400 CollRate[ipLo][ipHi] = CollRate[ipHi][ipLo] * g[ipHi] / g[ipLo] *
401 dsexp( (*tr).EnergyK() / phycon.te );
402
403 /* now add in excitations resulting from cosmic ray secondaries */
404 // \todo 2 add branch to do forbidden transitions
405 // this g-bar only works for permitted lines
406 if( (*tr).ipCont() > 0 )
407 {
408 /* get secondaries for all permitted lines by scaling LyA
409 * excitation by ratio of cross section (oscillator strength/energy)
410 * Born approximation or plane-wave approximation */
411 (*tr).Coll().rate_lu_nontherm_set() = secondaries.x12tot *
412 ((*tr).Emis().gf()/(*tr).EnergyWN()) /
413 (iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,0).Emis().gf()/
414 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,0).EnergyWN());
415
416 CollRate[ipLo][ipHi] += (*tr).Coll().rate_lu_nontherm();
417 CollRate[ipHi][ipLo] += (*tr).Coll().rate_lu_nontherm() * (*(*tr).Lo()).g() / (*(*tr).Hi()).g();
418 }
419 }
420
421 multi_arr<double,2> Cool(dBaseSpecies[ipSpecies].numLevels_local, dBaseSpecies[ipSpecies].numLevels_local);
422
423 /* solve the n-level atom */
425 /* dBaseSpecies[ipSpecies].numLevels_local is the number of levels to compute*/
426 dBaseSpecies[ipSpecies].numLevels_local,
427 /* ABUND is total abundance of species, used for nth equation
428 * if balance equations are homogeneous */
429 abund,
430 /* g(dBaseSpecies[ipSpecies].numLevels_local) is statistical weight of levels */
431 g,
432 /* EX(dBaseSpecies[ipSpecies].numLevels_local) is excitation potential of levels, deg K or wavenumbers
433 * 0 for lowest level, all are energy rel to ground NOT d(ENER) */
434 ex,
435 /* this is 'K' for ex[] as Kelvin deg, is 'w' for wavenumbers */
436 'w',
437 /* populations [cm-3] of each level as deduced here */
438 pops,
439 /* departure coefficient, derived below */
440 depart,
441 /* net transition rate, A * esc prob, s-1 */
442 &AulEscp,
443 /* col str from up to low */
444 &col_str,
445 /* AulDest[ihi][ilo] is destruction rate, trans from ihi to ilo, A * dest prob,
446 * asserts confirm that [ihi][ilo] is zero */
447 &AulDest,
448 /* AulPump[ilo][ihi] is pumping rate from lower to upper level (s^-1), (hi,lo) must be zero */
449 &AulPump,
450 /* collision rates (s^-1), evaluated here and returned for cooling by calling function,
451 * unless following flag is true. If true then calling function has already filled
452 * in these rates. CollRate[ipSpecies][j] is rate from ipSpecies to j */
453 &CollRate,
454 /* this is an additional creation rate from continuum, normally zero, units cm-3 s-1 */
455 source,
456 /* this is an additional destruction rate to continuum, normally zero, units s-1 */
457 sink,
458 // flag saying whether CollRate already done (true), or we need to do it here (false),
459 true,
460 /* total cooling and its derivative, set here but nothing done with it*/
461 &cooltl,
462 &coolder,
463 /* string used to identify calling program in case of error */
464 spName,
465 /* nNegPop flag indicating what we have done
466 * positive if negative populations occurred
467 * zero if normal calculation done
468 * negative if too cold (for some atoms other routine will be called in this case) */
469 &nNegPop,
470 /* true if populations are zero, either due to zero abundance of very low temperature */
471 &lgZeroPop,
472 /* option to print debug information */
473 lgDeBug,
474 /* option to do the molecule in LTE */
475 dBaseSpecies[ipSpecies].lgLTE,
476 /* cooling per line */
477 &Cool);
478
479 if( nNegPop > 0 )
480 {
481 /* negative populations occurred */
482 fprintf(ioQQQ," PROBLEM in dBase_solve, atom_levelN returned negative population .\n");
484 }
485
486 // highest levels may have no population
487 while( (pops[dBaseSpecies[ipSpecies].numLevels_local-1]<=0 ) &&
488 (dBaseSpecies[ipSpecies].numLevels_local > 1) )
489 --dBaseSpecies[ipSpecies].numLevels_local;
490
491 for( long j=0;j< dBaseSpecies[ipSpecies].numLevels_local; j++ )
492 dBaseStates[ipSpecies][j].Pop() = pops[j];
493
494 for( long j=dBaseSpecies[ipSpecies].numLevels_local;
495 j< dBaseSpecies[ipSpecies].numLevels_max; j++ )
496 dBaseStates[ipSpecies][j].Pop() = 0.;
497
498# ifdef USE_NLTE7
499 // do sums of totals out (destruction) and into (creation) level
500 for( int lvl = 0; lvl < dBaseSpecies[ipSpecies].numLevels_local; lvl++)
501 {
502 for( int lvl2 = 0; lvl2 < dBaseSpecies[ipSpecies].numLevels_local; lvl2++)
503 {
504 if( lvl != lvl2 )
505 {
506 dBaseStates[ipSpecies][lvl].DestCollBB() += CollRate[lvl][lvl2];
507 dBaseStates[ipSpecies][lvl].CreatCollBB() += CollRate[lvl2][lvl]*pops[lvl2];
508 if( lvl2 < lvl)
509 dBaseStates[ipSpecies][lvl].DestPhotoBB() += AulDest[lvl][lvl2];
510 }
511 }
512 }
513# endif
514
515 /*Atmol line*/
516 for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
517 tr != dBaseTrans[ipSpecies].end(); ++tr)
518 {
519 int ipHi = (*tr).ipHi();
520 if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local)
521 continue;
522 int ipLo = (*tr).ipLo();
523 (*tr).Coll().cool() = max(Cool[ipHi][ipLo],0.);
524 (*tr).Coll().heat() = max(-Cool[ipHi][ipLo],0.);
525
526 if ( (*tr).ipCont() > 0 )
527 {
528
529 (*tr).Emis().phots() = (*tr).Emis().Aul() * (*(*tr).Hi()).Pop() *
530 ((*tr).Emis().Pesc() + (*tr).Emis().Pelec_esc());
531
532 (*tr).Emis().xIntensity() = (*tr).Emis().phots() * (*tr).EnergyErg();
533
534 /* population of lower level rel to ion, corrected for stim em */
535 (*tr).Emis().PopOpc() = (*(*tr).Lo()).Pop() - (*(*tr).Hi()).Pop()*
536 (*(*tr).Lo()).g()/(*(*tr).Hi()).g();
537
538 /* it's possible for this sum to be zero. Set ratio to zero in that case. */
539 if( CollRate[ipLo][ipHi]+AulPump[ipLo][ipHi] > 0. )
540 {
541 (*tr).Emis().ColOvTot() = CollRate[ipLo][ipHi]/
542 (CollRate[ipLo][ipHi]+AulPump[ipLo][ipHi]);
543 }
544 else
545 (*tr).Emis().ColOvTot() = 0.;
546
547 // this is only used for save line printout. Maybe colliders may be involved, but
548 // this simple approximation of a "collision strength" should be good enough for
549 // the purposes of that printout.
550 (*tr).Coll().col_str() = (realnum)( (*tr).Coll().rate_coef_ul()[ipELECTRON] *
551 (g[ipHi]*phycon.sqrte)/COLL_CONST);
552 }
553 else
554 {
555 (*tr).Coll().col_str() = 0.;
556 }
557 }
558
559 for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
560 tr != dBaseTrans[ipSpecies].end(); ++tr)
561 {
562 int ipHi = (*tr).ipHi();
563 if (ipHi < dBaseSpecies[ipSpecies].numLevels_local)
564 continue;
565 EmLineZero( (*tr).Emis() );
566 }
567
568 dBaseSpecies[ipSpecies].CoolTotal = cooltl;
569 CoolAdd( dBaseSpecies[ipSpecies].chLabel, 0., max(0.,dBaseSpecies[ipSpecies].CoolTotal) );
570 if( dBaseSpecies[ipSpecies].lgMolecular )
571 thermal.elementcool[LIMELM] += max(0.,dBaseSpecies[ipSpecies].CoolTotal);
572 else
573 thermal.elementcool[dBaseStates[ipSpecies][0].nelem()-1] += max(0.,dBaseSpecies[ipSpecies].CoolTotal);
574 totalHeating += max(0., -dBaseSpecies[ipSpecies].CoolTotal);
575 thermal.dCooldT += coolder;
576
577 /* option to print departure coefficients */
578 {
579 enum {DEBUG_LOC=false};
580
581 if( DEBUG_LOC )
582 {
583 fprintf( ioQQQ, " Departure coefficients for species %li\n", ipSpecies );
584 for( long j=0; j< dBaseSpecies[ipSpecies].numLevels_local; j++ )
585 {
586 fprintf( ioQQQ, " level %li \t Depar Coef %e\n", j, depart[j] );
587 }
588 }
589 }
590 }
591
592 // total heating for all dBase dBaseSpecies
593 thermal.heating[0][27] = totalHeating;
594
595 return;
596}
597
598/*Leiden*/
599STATIC double LeidenCollRate(long ipSpecies, long ipCollider, const TransitionProxy& tr, double ftemp)
600{
601 DEBUG_ENTRY("LeidenCollRate()");
602 double ret_collrate = InterpCollRate( AtmolCollRateCoeff[ipSpecies][ipCollider], tr.ipHi(), tr.ipLo(), ftemp);
603 return ret_collrate;
604}
605
606/*STOUT*/
607STATIC double StoutCollRate(long ipSpecies, long ipCollider, const TransitionProxy& tr, double ftemp)
608{
609 DEBUG_ENTRY("StoutCollRate()");
610
611 double rate = 0.;
612 // deexcitation rate or collision strength?
613 bool lgIsRate = StoutCollData[ipSpecies][tr.ipHi()][tr.ipLo()][ipCollider].lgIsRate;
614 int n = StoutCollData[ipSpecies][tr.ipHi()][tr.ipLo()][ipCollider].ntemps;
615 if( n < 2)
616 return 0.;
617
618 double *x = (double*)MALLOC(n*sizeof(double));
619 double *y = (double*)MALLOC(n*sizeof(double));
620
621 double fupsilon = 0.;
622 for(int j = 0; j < n; j ++)
623 {
624 x[j] = StoutCollData[ipSpecies][tr.ipHi()][tr.ipLo()][ipCollider].temps[j];
625 y[j] = StoutCollData[ipSpecies][tr.ipHi()][tr.ipLo()][ipCollider].collstrs[j];
626 ASSERT( x[j] > 0. && y[j] > 0.);
627 }
628 //If the temperature is above or below the temperature range, use the CS from the closest temperature.
629 //Otherwise, do the linear interpolation.
630 if( ftemp < x[0] )
631 {
632 fupsilon = y[0];
633 }
634 else if( ftemp > x[n-1] )
635 {
636 fupsilon = y[n-1];
637 }
638 else
639 {
640 fupsilon = linint(x,y,n,ftemp);
641 }
642
643 free(x);
644 free(y);
645 ASSERT(fupsilon > 0);
646
647 /* We can deal with derexcitation rates and collision strengths currently */
648 if( lgIsRate )
649 {
650 rate = fupsilon;
651 }
652 else
653 {
654 /* convert the collision strength to a collision rate coefficient */
655 /* This formula converting collision strength to collision rate coefficient works fine for the electrons*/
656 /* For any other collider the mass would be different*/
657 rate = (COLL_CONST*fupsilon)/((*tr.Hi()).g()*sqrt(ftemp));
658 }
659
660 return rate;
661}
662
663/*CHIANTI*/
664STATIC double ChiantiCollRate(long ipSpecies, long ipCollider, const TransitionProxy& tr, double ftemp)
665{
666 DEBUG_ENTRY("ChiantiCollRate()");
667
668 double rate = 0.;
669 double fupsilon = CHIANTI_Upsilon(ipSpecies, ipCollider, tr.ipHi(), tr.ipLo(), ftemp);
670
671 /* NB NB - if proton colliders, the upsilons returned here are actually already rate coefficients. */
672 /* these are designated by a collider index and a transition type */
673 if( ipCollider == ipPROTON && AtmolCollSplines[ipSpecies][tr.ipHi()][tr.ipLo()][ipCollider].intTranType == 6 )
674 {
675 rate = fupsilon;
676 }
677 else if( ipCollider == ipELECTRON )
678 {
679 /* convert the collision strength to a collision rate coefficient */
680 /*This formula converting collision strength to collision rate coefficient works fine for the electrons*/
681 /*For any other collider the mass would be different*/
682 rate = (COLL_CONST*fupsilon)/((*tr.Hi()).g()*sqrt(ftemp));
683 }
684 else
685 rate = 0.;
686
687 return rate;
688}
689
690double CHIANTI_Upsilon(long ipSpecies, long ipCollider, long ipHi, long ipLo, double ftemp)
691{
692 double fdeltae,fscalingparam,fkte,fxt,fsups,fups;
693 int intxs,inttype,intsplinepts;
694
695 DEBUG_ENTRY("CHIANTI_Upsilon()");
696
697 if( AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline == NULL )
698 {
699 return 0.;
700 }
701
702 intsplinepts = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].nSplinePts;
703 inttype = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].intTranType;
704 fdeltae = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].EnergyDiff;
705 fscalingparam = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].ScalingParam;
706
707 fkte = ftemp/fdeltae/1.57888e5;
708
709 /*Way the temperature is scaled*/
710 /*Burgess&Tully 1992:Paper gives only types 1 to 4*/
711 /*Found that the expressions were the same for 5 & 6 from the associated routine DESCALE_ALL*/
712 /*What about 7,8&9?*/
713 if( inttype ==1 || inttype==4 )
714 {
715 fxt = 1-(log(fscalingparam)/(log(fkte+fscalingparam)));
716 }
717 else if(inttype == 2 || inttype == 3||inttype == 5 || inttype == 6)
718 {
719 fxt = fkte/(fkte+fscalingparam);
720 }
721 else
723
724 double xs[9],*spl,*spl2;
725 /*Creating spline points array*/
726 if(intsplinepts == 5)
727 {
728 spl = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline;
729 for(intxs=0;intxs<5;intxs++)
730 {
731 xs[intxs] = 0.25*intxs;
732 if(DEBUGSTATE)
733 {
734 printf("The xs and spl values are %f and %f \n",xs[intxs],spl[intxs]);
735 getchar();
736 }
737 }
738 }
739 else if(intsplinepts == 9)
740 {
741 spl = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline;
742 for( intxs=0; intxs<9; intxs++ )
743 {
744 xs[intxs] = 0.125*intxs;
745 if(DEBUGSTATE)
746 {
747 printf("The xs and spl values are %f and %f \n",xs[intxs],spl[intxs]);
748 getchar();
749 }
750 }
751 }
752 else
753 {
755 }
756
757 /*Finding the second derivative*/
758 spl2 = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].SplineSecDer;
759
760 if(DEBUGSTATE)
761 {
762 printf("\n");
763 for(intxs=0;intxs<intsplinepts;intxs++)
764 {
765 printf("The %d value of 2nd derivative is %f \n",intxs+1,spl2[intxs]);
766 }
767 }
768
769 /*Extracting out the value*/
770 //splint(xs,spl,spl2,intsplinepts,fxt,&fsups);
771
772 fsups = linint( xs, spl, intsplinepts, fxt);
773
774 /*Finding upsilon*/
775 if(inttype == 1)
776 {
777 fups = fsups*log(fkte+exp(1.0));
778 }
779 else if(inttype == 2)
780 {
781 fups = fsups;
782 }
783 else if(inttype == 3)
784 {
785 fups = fsups/(fkte+1.0) ;
786 }
787 else if(inttype == 4)
788 {
789 fups = fsups*log(fkte+fscalingparam) ;
790 }
791 else if(inttype == 5)
792 {
793 fups = fsups/fkte ;
794 }
795 else if(inttype == 6)
796 {
797 fups = pow(10.0,fsups) ;
798 }
799 else
800 {
802 }
803
804 if( fups < 0. )
805 {
806 fprintf( ioQQQ," WARNING: Negative upsilon in species %s, collider %li, indices %4li %4li, Te = %e.\n",
807 dBaseSpecies[ipSpecies].chLabel, ipCollider, ipHi, ipLo, ftemp );
808 fups = 0.;
809 }
810 ASSERT(fups>=0);
811 return(fups);
812}
813
t_abund abund
Definition abund.cpp:5
double InterpCollRate(const CollRateCoeffArray &rate_table, const long &ipHi, const long &ipLo, const double &ftemp)
Definition atmdat.cpp:135
t_atmdat atmdat
Definition atmdat.cpp:6
#define DEBUGSTATE
void atom_levelN(long int nLevelCalled, realnum abund, const double g[], const double ex[], char chExUnits, double pops[], double depart[], double ***AulEscp, double ***col_str, double ***AulDest, double ***AulPump, double ***CollRate, const double source[], const double sink[], bool lgCollRateDone, double *cooltl, double *coolder, const char *chLabel, int *nNegPop, bool *lgZeroPop, bool lgDeBug, bool lgLTE, multi_arr< double, 2 > *Cool, multi_arr< double, 2 > *dCooldT)
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define MIN2
Definition cddefines.h:761
#define MALLOC(exp)
Definition cddefines.h:501
const int LIMELM
Definition cddefines.h:258
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
long max(int a, long b)
Definition cddefines.h:775
NORETURN void TotalInsanity(void)
Definition service.cpp:886
const int ipHYDROGEN
Definition cddefines.h:305
double dsexp(double x)
Definition service.cpp:953
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
const double * rate_coef_ul() const
Definition collision.h:176
realnum & col_str() const
Definition collision.h:167
double * rate_coef_ul_set() const
Definition collision.h:172
TransitionProxy::iterator iterator
Definition transition.h:280
int & ipLo() const
Definition transition.h:458
int & ipHi() const
Definition transition.h:466
qList::iterator Hi() const
Definition transition.h:396
double den
Definition mole.h:358
ColliderList colliders
Definition collision.cpp:57
void CollisionZero(const CollisionProxy &t)
Definition collision.cpp:81
@ ipATOM_H
Definition collision.h:13
@ ipH2_ORTHO
Definition collision.h:15
@ ipHE_PLUS
Definition collision.h:11
@ ipH2_PARA
Definition collision.h:16
@ ipH2
Definition collision.h:17
@ ipPROTON
Definition collision.h:10
@ ipNCOLLIDER
Definition collision.h:18
@ ipELECTRON
Definition collision.h:9
@ ipATOM_HE
Definition collision.h:14
t_conv conv
Definition conv.cpp:5
void CoolAdd(const char *chLabel, realnum lambda, double cool)
Definition cool_etc.cpp:13
double Fe3_cs(long ipLo, long ipHi)
double Fe5_cs(long ipLo, long ipHi)
double Fe4_cs(long ipLo, long ipHi)
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
void EmLineZero(EmissionList::reference t)
Definition emission.cpp:92
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH2p
Definition iso.h:29
const int ipH_LIKE
Definition iso.h:62
molezone * findspecieslocal(const char buf[])
molezone * null_molezone
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double COLL_CONST
Definition physconst.h:229
t_secondaries secondaries
static double ** AulEscp
Definition species2.cpp:29
double CHIANTI_Upsilon(long ipSpecies, long ipCollider, long ipHi, long ipLo, double ftemp)
Definition species2.cpp:690
static double * source
Definition species2.cpp:28
static double * ex
Definition species2.cpp:28
static double * pops
Definition species2.cpp:28
static double ** col_str
Definition species2.cpp:29
STATIC double LeidenCollRate(long, long, const TransitionProxy &, double)
Definition species2.cpp:599
static double * depart
Definition species2.cpp:28
void dBase_solve(void)
Definition species2.cpp:33
static double ** AulPump
Definition species2.cpp:29
static double * g
Definition species2.cpp:28
STATIC double StoutCollRate(long ipSpecies, long ipCollider, const TransitionProxy &, double ftemp)
Definition species2.cpp:607
static double ** AulDest
Definition species2.cpp:29
STATIC double ChiantiCollRate(long ipSpecies, long ipCollider, const TransitionProxy &, double ftemp)
Definition species2.cpp:664
static double * sink
Definition species2.cpp:28
static double ** CollRate
Definition species2.cpp:29
long int nSpecies
Definition taulines.cpp:21
vector< qList > dBaseStates
Definition taulines.cpp:15
vector< TransitionList > dBaseTrans
Definition taulines.cpp:17
multi_arr< CollRateCoeffArray, 2 > AtmolCollRateCoeff
Definition taulines.cpp:18
StoutColls **** StoutCollData
Definition taulines.cpp:20
species * dBaseSpecies
Definition taulines.cpp:14
CollSplinesArray **** AtmolCollSplines
Definition taulines.cpp:19
t_thermal thermal
Definition thermal.cpp:5
double linint(const double x[], const double y[], long n, double xval)
void MakeCS(const TransitionProxy &t)