cloudy  trunk
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 
22 STATIC double LeidenCollRate(long, long, const TransitionProxy& ,double);
23 STATIC double StoutCollRate(long ipSpecies, long ipCollider, const TransitionProxy&, double ftemp);
24 STATIC double ChiantiCollRate(long ipSpecies, long ipCollider, const TransitionProxy&, double ftemp);
25 
26 static const bool DEBUGSTATE = false;
27 
28 static double *g, *ex, *pops, *depart, *source, *sink;
29 static double **AulEscp, **col_str, **AulDest, **AulPump, **CollRate;
30 
31 /*Solving for the level populations*/
32 
33 void dBase_solve(void )
34 {
35  realnum abund;
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
183  TotalInsanity();
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
288  TotalInsanity();
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 */
424  atom_levelN(
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");
483  cdEXIT( EXIT_FAILURE );
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*/
599 STATIC 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*/
607 STATIC 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*/
664 STATIC 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 
690 double 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
722  TotalInsanity();
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  {
754  TotalInsanity();
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  {
801  TotalInsanity();
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 
thermal.h
findspecieslocal
molezone * findspecieslocal(const char buf[])
Definition: mole_species.cpp:833
AtmolCollSplines
CollSplinesArray **** AtmolCollSplines
Definition: taulines.cpp:19
StoutCollRate
STATIC double StoutCollRate(long ipSpecies, long ipCollider, const TransitionProxy &, double ftemp)
Definition: species2.cpp:607
colliders
ColliderList colliders
Definition: collision.cpp:57
dense
t_dense dense
Definition: dense.cpp:24
t_CollSplinesArray::intTranType
long intTranType
Definition: cddefines.h:1286
secondaries.h
atoms.h
ipHE_PLUS
@ ipHE_PLUS
Definition: collision.h:11
MakeCS
void MakeCS(const TransitionProxy &t)
Definition: transition.cpp:613
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
dBaseTrans
vector< TransitionList > dBaseTrans
Definition: taulines.cpp:17
CollisionProxy::rate_coef_ul
const double * rate_coef_ul() const
Definition: collision.h:176
ipH2_PARA
@ ipH2_PARA
Definition: collision.h:16
realnum
float realnum
Definition: cddefines.h:103
conv.h
STATIC
#define STATIC
Definition: cddefines.h:97
mole.h
TransitionProxy::ipLo
int & ipLo() const
Definition: transition.h:458
multi_arr< double, 2 >
linint
double linint(const double x[], const double y[], long n, double xval)
Definition: thirdparty_interpolate.cpp:456
thirdparty.h
phycon
t_phycon phycon
Definition: phycon.cpp:6
t_StoutColls::lgIsRate
bool lgIsRate
Definition: cddefines.h:1303
Fe3_cs
double Fe3_cs(long ipLo, long ipHi)
Definition: cool_iron.cpp:1036
ProxyIterator
Definition: proxy_iterator.h:58
ex
static double * ex
Definition: species2.cpp:28
ipH2_ORTHO
@ ipH2_ORTHO
Definition: collision.h:15
t_atmdat::collstrDefault
double collstrDefault
Definition: atmdat.h:254
lines_service.h
t_secondaries::x12tot
realnum x12tot
Definition: secondaries.h:53
CoolAdd
void CoolAdd(const char *chLabel, realnum lambda, double cool)
Definition: cool_etc.cpp:13
DEBUGSTATE
static const bool DEBUGSTATE
Definition: species2.cpp:26
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
Fe4_cs
double Fe4_cs(long ipLo, long ipHi)
Definition: cool_iron.cpp:1159
ipATOM_H
@ ipATOM_H
Definition: collision.h:13
atom_levelN
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)
Definition: atom_leveln.cpp:15
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
atmdat.h
null_molezone
molezone * null_molezone
Definition: mole_species.cpp:66
t_StoutColls::temps
double * temps
Definition: cddefines.h:1299
MIN2
#define MIN2
Definition: cddefines.h:761
TransitionProxy
Definition: transition.h:23
ipPROTON
@ ipPROTON
Definition: collision.h:10
t_species::chLabel
char * chLabel
Definition: cddefines.h:1235
t_StoutColls::collstrs
double * collstrs
Definition: cddefines.h:1301
source
static double * source
Definition: species2.cpp:28
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
t_species::numLevels_local
long numLevels_local
Definition: cddefines.h:1241
CollisionProxy::rate_coef_ul_set
double * rate_coef_ul_set() const
Definition: collision.h:172
AulPump
static double ** AulPump
Definition: species2.cpp:29
col_str
static double ** col_str
Definition: species2.cpp:29
dense.h
t_CollSplinesArray::collspline
double * collspline
Definition: cddefines.h:1282
CollisionZero
void CollisionZero(const CollisionProxy &t)
Definition: collision.cpp:81
cooling.h
cddefines.h
CollRate
static double ** CollRate
Definition: species2.cpp:29
ipH2
@ ipH2
Definition: collision.h:17
thermal
t_thermal thermal
Definition: thermal.cpp:5
TransitionProxy::Hi
qList::iterator Hi() const
Definition: transition.h:396
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
ipATOM_HE
@ ipATOM_HE
Definition: collision.h:14
t_thermal::heating
double heating[LIMELM][LIMELM]
Definition: thermal.h:158
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
t_CollSplinesArray::ScalingParam
double ScalingParam
Definition: cddefines.h:1288
t_species::lgActive
bool lgActive
Definition: cddefines.h:1255
ipNCOLLIDER
@ ipNCOLLIDER
Definition: collision.h:18
CollisionProxy
Definition: collision.h:79
MAX2
#define MAX2
Definition: cddefines.h:782
CollisionProxy::col_str
realnum & col_str() const
Definition: collision.h:167
LIMELM
const int LIMELM
Definition: cddefines.h:258
ipELECTRON
@ ipELECTRON
Definition: collision.h:9
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
StoutCollData
StoutColls **** StoutCollData
Definition: taulines.cpp:20
t_CollSplinesArray::SplineSecDer
double * SplineSecDer
Definition: cddefines.h:1283
t_atmdat::lgGbarOn
bool lgGbarOn
Definition: atmdat.h:249
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
ChiantiCollRate
STATIC double ChiantiCollRate(long ipSpecies, long ipCollider, const TransitionProxy &, double ftemp)
Definition: species2.cpp:664
abund
t_abund abund
Definition: abund.cpp:5
LeidenCollRate
STATIC double LeidenCollRate(long, long, const TransitionProxy &, double)
Definition: species2.cpp:599
ipH2p
const int ipH2p
Definition: iso.h:29
t_species::fracType
realnum fracType
Definition: cddefines.h:1249
TransitionProxy::ipHi
int & ipHi() const
Definition: transition.h:466
t_conv::nTotalIoniz
long int nTotalIoniz
Definition: conv.h:166
TransitionProxy::EnergyWN
realnum & EnergyWN() const
Definition: transition.h:438
AulEscp
static double ** AulEscp
Definition: species2.cpp:29
COLL_CONST
const UNUSED double COLL_CONST
Definition: physconst.h:229
t_species::CoolTotal
double CoolTotal
Definition: cddefines.h:1253
dBase_solve
void dBase_solve(void)
Definition: species2.cpp:33
physconst.h
secondaries
t_secondaries secondaries
Definition: secondaries.cpp:5
molezone
Definition: mole.h:326
t_dense::xNucleiTotal
realnum xNucleiTotal
Definition: dense.h:104
conv
t_conv conv
Definition: conv.cpp:5
t_CollSplinesArray::nSplinePts
long nSplinePts
Definition: cddefines.h:1285
t_thermal::elementcool
double elementcool[LIMELM+1]
Definition: thermal.h:95
dBaseSpecies
species * dBaseSpecies
Definition: taulines.cpp:14
InterpCollRate
double InterpCollRate(const CollRateCoeffArray &rate_table, const long &ipHi, const long &ipLo, const double &ftemp)
Definition: atmdat.cpp:135
t_thermal::dCooldT
double dCooldT
Definition: thermal.h:119
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
taulines.h
pops
static double * pops
Definition: species2.cpp:28
phycon.h
dBaseStates
vector< qList > dBaseStates
Definition: taulines.cpp:15
t_CollSplinesArray::EnergyDiff
double EnergyDiff
Definition: cddefines.h:1287
atmdat
t_atmdat atmdat
Definition: atmdat.cpp:6
t_species::fracIsotopologue
realnum fracIsotopologue
Definition: cddefines.h:1251
depart
static double * depart
Definition: species2.cpp:28
Fe5_cs
double Fe5_cs(long ipLo, long ipHi)
Definition: cool_iron.cpp:1187
t_phycon::sqrte
double sqrte
Definition: phycon.h:48
AtmolCollRateCoeff
multi_arr< CollRateCoeffArray, 2 > AtmolCollRateCoeff
Definition: taulines.cpp:18
t_species::numLevels_max
long numLevels_max
Definition: cddefines.h:1239
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
mole_priv.h
EmLineZero
void EmLineZero(EmissionList::reference t)
Definition: emission.cpp:92
dsexp
double dsexp(double x)
Definition: service.cpp:953
h2.h
CHIANTI_Upsilon
double CHIANTI_Upsilon(long ipSpecies, long ipCollider, long ipHi, long ipLo, double ftemp)
Definition: species2.cpp:690
molezone::den
double den
Definition: mole.h:358
nSpecies
long int nSpecies
Definition: taulines.cpp:21
t_conv::lgSearch
bool lgSearch
Definition: conv.h:175
t_phycon::te
double te
Definition: phycon.h:11
AulDest
static double ** AulDest
Definition: species2.cpp:29
sink
static double * sink
Definition: species2.cpp:28
max
long max(int a, long b)
Definition: cddefines.h:775
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:191
t_StoutColls::ntemps
long ntemps
Definition: cddefines.h:1297
g
static double * g
Definition: species2.cpp:28
Generated by   doxygen 1.8.17