cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
nemala2.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 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 "phycon.h"
5 #include "taulines.h"
6 #include "mole.h"
7 #include "atoms.h"
8 #include "string.h"
9 #include "thirdparty.h"
10 #include "dense.h"
11 #include "conv.h"
12 #include "h2.h"
13 #include "physconst.h"
14 
15 void states_popfill( void);
16 STATIC double LeidenCollRate(long, long, long, long,double);
17 STATIC double Chianti_Upsilon(long, long, long, long,double);
18 
19 #define DEBUGSTATE false
20 
21 
22 /*Solving for the level populations*/
23 
24 void atmol_popsolve(void )
25 {
26  realnum abund;
27  DEBUG_ENTRY( "atmol_popsolve()" );
28 
29  for( long i=0; i<nSpecies; i++ )
30  {
31  const char *spName = Species[i].chptrSpName;
32  double *g, *ex, *pops, *depart;
33  double **AulEscp, **col_str, **AulDest, **AulPump, **CollRate,fcolldensity[9];
34  double *source, *sink;
35  double cooltl, coolder;
36  double fupsilon;
37  long ipLo,ipHi,intCollNo;
38  int nNegPop;
39  bool lgZeroPop, lgDeBug = false;
40  long nlev = Species[i].numLevels_max;
41 
42  /* first find current density (cm-3) of species */
43  if( lgSpeciesMolecule[i] )
44  {
45  struct molecule *SpeciesCurrent;
48  if( (SpeciesCurrent = findspecies(Species[i].chptrSpName)) == &null_mole )
49  {
50  /* did not find the species - print warning for now */
51  fprintf(ioQQQ," PROBLEM atmol_popsolve did not find molecular species %li\n",i);
52  }
53  abund = SpeciesCurrent->hevmol;
54  }
55  else
56  {
57  /* an atom or ion */
58  ASSERT( Species[i].intAtNo<LIMELM && Species[i].intIonStage<=Species[i].intAtNo);
60  }
61 
62  if(abund<=SMALLFLOAT)
63  {
64  /* zero out populations and intensities */
65  /* NB - can't use states_popfill here, because it zeros ALL species */
66  atmolStates[i][0].Pop = 0.;
67  for(long ipHi = 1; ipHi < nlev; ipHi++ )
68  {
69  atmolStates[i][ipHi].Pop = 0.;
70  for(long ipLo = 0; ipLo < ipHi; ipLo++ )
71  {
72  if( atmolTrans[i][ipHi][ipLo].ipCont > 0 )
73  {
74  atmolTrans[i][ipHi][ipLo].Emis->phots = 0.;
75  atmolTrans[i][ipHi][ipLo].Emis->xIntensity = 0.;
76  }
77  }
78  }
79  continue;
80  }
81 
82  g = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double));
83  ex = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double));
84  pops = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double));
85  depart = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double));
86  source = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double));
87  sink = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double));
88 
89  AulEscp = (double **)MALLOC( (unsigned long)(nlev)*sizeof(double *));
90  col_str = (double **)MALLOC( (unsigned long)(nlev)*sizeof(double *));
91  AulDest = (double **)MALLOC( (unsigned long)(nlev)*sizeof(double *));
92  AulPump = (double **)MALLOC( (unsigned long)(nlev)*sizeof(double *));
93  CollRate = (double **)MALLOC( (unsigned long)(nlev)*sizeof(double *));
94 
95  for( long j=0; j< nlev; j++ )
96  {
97  AulEscp[j] = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double));
98  col_str[j] = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double));
99  AulDest[j] = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double));
100  AulPump[j] = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double));
101  CollRate[j] = (double *)MALLOC( (unsigned long)(nlev)*sizeof(double));
102  }
103 
104  for( ipLo = 0; ipLo < nlev; ipLo++ )
105  {
106  /*Filling in statistical weights & Excitation Energies*/
107  g[ipLo] = atmolStates[i][ipLo].g ;
108  ex[ipLo] = atmolStates[i][ipLo].energy;
109  /* zero some rates */
110  source[ipLo] = 0.;
111  sink[ipLo] = 0.;
112  AulEscp[ipLo][ipLo] = 0.;
113  AulDest[ipLo][ipLo] = 0.;
114  AulPump[ipLo][ipLo] = 0.;
115  }
116 
117  for( ipHi=1; ipHi<nlev; ipHi++ )
118  {
119  for( ipLo=0; ipLo<ipHi; ipLo++ )
120  {
121  if( atmolTrans[i][ipHi][ipLo].ipCont > 0 )
122  {
123  AulEscp[ipHi][ipLo] = atmolTrans[i][ipHi][ipLo].Emis->Aul*
124  (atmolTrans[i][ipHi][ipLo].Emis->Pesc +
125  atmolTrans[i][ipHi][ipLo].Emis->Pelec_esc);
126  AulDest[ipHi][ipLo] = atmolTrans[i][ipHi][ipLo].Emis->Aul*
127  atmolTrans[i][ipHi][ipLo].Emis->Pdest;
128  AulPump[ipLo][ipHi] = atmolTrans[i][ipHi][ipLo].Emis->pump;
129  }
130  else
131  {
132  AulEscp[ipHi][ipLo] = SMALLFLOAT;
133  AulDest[ipHi][ipLo] = SMALLFLOAT;
134  AulPump[ipLo][ipHi] = SMALLFLOAT;
135  }
136 
137  AulEscp[ipLo][ipHi] = 0.;
138  AulDest[ipLo][ipHi] = 0.;
139  AulPump[ipHi][ipLo] = 0.;
140  }
141  }
142 
143  /*Setting all the collision strengths and collision rate to zero*/
144  for( ipHi= 0; ipHi<nlev; ipHi++)
145  {
146  for( ipLo= 0; ipLo<nlev; ipLo++)
147  {
148  col_str[ipHi][ipLo] = 0.;
149  CollRate[ipHi][ipLo] = 0.;
150  }
151  }
152 
153  /*Updating the CollRatesArray*/
154  /*You need to do this for a species indexed by i*/
155  for( intCollNo=0; intCollNo<NUM_COLLIDERS; intCollNo++)
156  {
157  for( ipHi=1; ipHi<nlev; ipHi++)
158  {
159  for( ipLo=0; ipLo<ipHi; ipLo++)
160  {
161  /* molecule */
162  if( lgSpeciesMolecule[i] )
163  {
164  if(CollRatesArray[i][intCollNo]!=NULL)
165  {
166  /*using the collision rate coefficients directly*/
167  CollRatesArray[i][intCollNo][ipHi][ipLo] =
168  LeidenCollRate(i, intCollNo, ipHi, ipLo, phycon.te);
169  }
170  }
171  /* atom or ion */
172  else
173  {
174  if(CollRatesArray[i][intCollNo]!=NULL)
175  {
176  fupsilon = Chianti_Upsilon(i, intCollNo, ipHi, ipLo, phycon.te);
177  /*converting the collision strength to a collision rate coefficient*/
178  /*This formula converting collision strength to collision rate coefficent works fine for the electrons*/
179  /*For any other collider the mass would be different*/
180  CollRatesArray[i][intCollNo][ipHi][ipLo] = ((8.6291e-6)*fupsilon)/(g[ipHi]*phycon.sqrte);
181  }
182  }
183 
184  /* now get the corresponding excitation rate */
185  if(CollRatesArray[i][intCollNo]!=NULL)
186  {
187  CollRatesArray[i][intCollNo][ipLo][ipHi] =
188  CollRatesArray[i][intCollNo][ipHi][ipLo] * g[ipHi] / g[ipLo] *
189  sexp( atmolTrans[i][ipHi][ipLo].EnergyK / phycon.te );
190  }
191  }
192  }
193  }
194  /*Get the densities seperately*/
195  /*Electron*/
196  fcolldensity[0] = dense.eden;
197  /*Proton*/
198  fcolldensity[1] = dense.xIonDense[1][1];
199  /*Atomic Hydrogen*/
200  fcolldensity[2] = dense.xIonDense[1][0];
201  /*Helium*/
202  fcolldensity[3] = dense.xIonDense[2][0];
203  /*Helium+*/
204  fcolldensity[4] = dense.xIonDense[2][1];
205  /*Helium ++*/
206  fcolldensity[5] = dense.xIonDense[2][2];
207  /*Molecular Hydrogen*/
208  fcolldensity[8] = findspecies("H2")->hevmol;
209  /*Ortho(6) and Para(7) Mol Hydrogen*/
210  if(h2.lgH2ON)
211  {
212  fcolldensity[6] = h2.ortho_density;
213  fcolldensity[7] = h2.para_density;
214  }
215  else
216  {
217  fcolldensity[6] = (fcolldensity[8])*(0.75);
218  fcolldensity[7] = (fcolldensity[8])*(0.25);
219  }
220 
221  /*Updating the CollRate*/
222  for( ipHi=0; ipHi<nlev; ipHi++ )
223  {
224  for( ipLo=0; ipLo<nlev; ipLo++ )
225  {
226  if( ipHi == ipLo )
227  {
228  CollRate[ipHi][ipLo] = 0.;
229  continue;
230  }
231 
232  for( intCollNo=0; intCollNo<NUM_COLLIDERS; intCollNo++ )
233  {
234  if(CollRatesArray[i][intCollNo]!=NULL)
235  {
236  CollRate[ipHi][ipLo] += CollRatesArray[i][intCollNo][ipHi][ipLo]*fcolldensity[intCollNo];
237  }
238  }
239  /*Correction for Helium*/
240  /*To include the effects of collision with Helium,the user must multiply the density by 1.14*/
241  /*pg 5,Schoier et al 2004*/
242  if(lgSpeciesMolecule[i])
243  {
244  /*The collision rate coefficients for helium should not be present and that for molecular hydrogen should be present*/
245  if(CollRatesArray[i][3]==NULL && CollRatesArray[i][8]!=NULL )
246  {
247  CollRate[ipHi][ipLo] += 0.7 *CollRatesArray[i][8][ipHi][ipLo]*fcolldensity[3];
248  }
249  }
250  }
251  }
253  /*Level Lowering as seen in mole_co_atom.cpp*/
254  /*Excitation energies need to be converted to temperature*/
255  /*ex[nUsed]/phycon.te cannot be much larger than 25, or matrix inversion will fail*/
256  while ( ((ex[nlev-1]*T1CM) > phycon.te*25.) && (nlev > 1) )
257  {
258  --nlev;
259  }
260  Species[i].numLevels_local = nlev;
261 
262  /* build some arrays */
263  atom_levelN(
264  /* nlev is the number of levels to compute*/
265  nlev,
266  /* ABUND is total abundance of species, used for nth equation
267  * if balance equations are homogeneous */
268  abund,
269  /* G(nlev) is stat weight of levels */
270  g,
271  /* EX(nlev) is excitation potential of levels, deg K or wavenumbers
272  * 0 for lowest level, all are energy rel to ground NOT d(ENER) */
273  ex,
274  /* this is 'K' for ex[] as Kelvin deg, is 'w' for wavenumbers */
275  'w',
276  /* populations [cm-3] of each level as deduced here */
277  pops,
278  /* departure coefficient, derived below */
279  depart,
280  /* net transition rate, A * esc prob, s-1 */
281  &AulEscp,
282  /* col str from up to low */
283  &col_str,
284  /* AulDest[ihi][ilo] is destruction rate, trans from ihi to ilo, A * dest prob,
285  * asserts confirm that [ihi][ilo] is zero */
286  &AulDest,
287  /* AulPump[ilo][ihi] is pumping rate from lower to upper level (s^-1), (hi,lo) must be zero */
288  &AulPump,
289  /* collision rates (s^-1), evaluated here and returned for cooling by calling function,
290  * unless following flag is true. If true then calling function has already filled
291  * in these rates. CollRate[i][j] is rate from i to j */
292  &CollRate,
293  /* this is an additional creation rate from continuum, normally zero, units cm-3 s-1 */
294  source,
295  /* this is an additional destruction rate to continuum, normally zero, units s-1 */
296  sink,
297  /* flag saying whether CollRate already done (true), or we need to do it here (false),
298  * this is stored in data)[ihi][ilo] as either downward rate or collis strength*/
299  true,
300  /* total cooling and its derivative, set here but nothing done with it*/
301  &cooltl,
302  &coolder,
303  /* string used to identify calling program in case of error */
304  spName,
305  /* nNegPop flag indicating what we have done
306  * positive if negative populations occurred
307  * zero if normal calculation done
308  * negative if too cold (for some atoms other routine will be called in this case) */
309  &nNegPop,
310  /* true if populations are zero, either due to zero abundance of very low temperature */
311  &lgZeroPop,
312  /* option to print debug information */
313  lgDeBug );
314 
315  if( nNegPop == 0 )
316  {
317  ASSERT( lgZeroPop == false );
318 
319  for( long j=0; j< nlev; j++ )
320  {
321  atmolStates[i][j].Pop = pops[j];
322  }
323  }
324  else if( nNegPop > 0 )
325  {
326  /* something should be done here. */
327  /* try another solver? */
328  /* probably don't want to quit */
329  fprintf(ioQQQ," PROBLEM, atom_levelN returned negative population .\n");
330  cdEXIT( EXIT_FAILURE );
331  }
332  else
333  {
334  /*In search phase we dont lower the levels*/
335  if(conv.lgSearch!= true)
336  {
337  for(long j=Species[i].numLevels_max-1; j>0; j--)
338  {
339  if((atmolStates[i][j].Pop/abund)<POPTHRES)
340  {
341  nlev--;
342  }
343  else
344  break;
345  }
346 
347  ASSERT( nlev >= 1 );
348  /*Setting the number of energy levels to new limit*/
349  Species[i].numLevels_local = nlev;
350 
351  if( nlev == 1 )
352  {
353  atmolStates[i][0].Pop = abund;
354  for( long j=1; j< Species[i].numLevels_max; j++ )
355  {
356  atmolStates[i][j].Pop = 0.;
357  }
358  }
359  else
360  {
361  /*Solving for level populations*/
362  atom_levelN(nlev,abund,g,ex,'w',pops,depart,&AulEscp,&col_str,
363  &AulDest, &AulPump, &CollRate, source, sink,true,&cooltl,
364  &coolder, spName, &nNegPop, &lgZeroPop, lgDeBug );
365 
366  for( long j=0; j< nlev; j++ )
367  {
368  atmolStates[i][j].Pop = pops[j];
369  }
370  for( long j=nlev; j< Species[i].numLevels_max; j++ )
371  {
372  atmolStates[i][j].Pop = 0.;
373  }
374  }
375  }
376  }
377 
378  /*Atmol line*/
379  for(long ipHi = 1; ipHi < nlev; ipHi++ )
380  {
381  for(long ipLo = 0; ipLo < ipHi; ipLo++ )
382  {
383  if( atmolTrans[i][ipHi][ipLo].ipCont > 0 )
384  {
385  atmolTrans[i][ipHi][ipLo].Emis->phots =
386  atmolTrans[i][ipHi][ipLo].Emis->Aul *
387  atmolTrans[i][ipHi][ipLo].Hi->Pop *
388  (atmolTrans[i][ipHi][ipLo].Emis->Pesc +
389  atmolTrans[i][ipHi][ipLo].Emis->Pelec_esc);
390 
391  atmolTrans[i][ipHi][ipLo].Emis->xIntensity =
392  atmolTrans[i][ipHi][ipLo].Emis->phots *
393  atmolTrans[i][ipHi][ipLo].EnergyErg;
394 
395  /* population of lower level rel to ion, corrected for stim em */
396 
397  atmolTrans[i][ipHi][ipLo].Emis->PopOpc =
398  atmolTrans[i][ipHi][ipLo].Lo->Pop - atmolTrans[i][ipHi][ipLo].Hi->Pop*
399  atmolTrans[i][ipHi][ipLo].Lo->g/atmolTrans[i][ipHi][ipLo].Hi->g;
400  fixit(); /* need to do something with these! */
401  atmolTrans[i][ipHi][ipLo].Emis->ColOvTot = 0.;
402  atmolTrans[i][ipHi][ipLo].Coll.cool = 0.;
403  atmolTrans[i][ipHi][ipLo].Coll.heat = 0.;
404  }
405  }
406  }
407 
408  /* option to print departure coefficients */
409  {
410  enum {DEBUG_LOC=false};
411 
412  if( DEBUG_LOC )
413  {
414  fprintf( ioQQQ, " Departure coefficients for species %li\n", i );
415  for( long j=0; j< nlev; j++ )
416  {
417  fprintf( ioQQQ, " level %li \t Depar Coef %e\n", j, depart[j] );
418  }
419  }
420  }
421 
422  /* free vectors */
423  free( g );
424  free( ex );
425  free( pops );
426  free( depart );
427  free( source );
428  free( sink );
429 
430  /* free arrays */
431  for( long j=0; j< nlev; j++ )
432  {
433  free( AulEscp[j] );
434  free( col_str[j] );
435  free( AulDest[j] );
436  free( AulPump[j] );
437  free( CollRate[j] );
438  }
439  free( AulEscp );
440  free( col_str );
441  free( AulDest );
442  free( AulPump );
443  free( CollRate );
444  }
445  return;
446 }
447 
448 /*This function fills the population of the states to a dummy value of 0*/
449 void states_popfill( void)
450 {
451  DEBUG_ENTRY( "states_popfill()" );
452 
453  for( long i=0; i<nSpecies; i++)
454  {
455  for( long j=0; j<Species[i].numLevels_max; j++)
456  {
457  atmolStates[i][j].Pop = 0.;
458  }
459  }
460  return;
461 }
462 
463 /*Leiden*/
464 double LeidenCollRate(long ipSpecies, long ipCollider, long ipHi, long ipLo, double ftemp)
465 {
466  double ret_collrate;
467  int inttemps;
468  DEBUG_ENTRY("LeidenCollRate()");
469  inttemps = AtmolCollRateCoeff[ipSpecies][ipCollider].ntemps;
470 
471  if( AtmolCollRateCoeff[ipSpecies][ipCollider].temps == NULL )
472  {
473  return 0.;
474  }
475 
476  if(ftemp < AtmolCollRateCoeff[ipSpecies][ipCollider].temps[0])
477  {
478  ret_collrate = AtmolCollRateCoeff[ipSpecies][ipCollider].collrates[ipHi][ipLo][0];
479  }
480  else if(ftemp > AtmolCollRateCoeff[ipSpecies][ipCollider].temps[inttemps-1])
481  {
482  ret_collrate = AtmolCollRateCoeff[ipSpecies][ipCollider].collrates[ipHi][ipLo][inttemps-1];
483  }
484  else
485  {
486  ret_collrate = linint(AtmolCollRateCoeff[ipSpecies][ipCollider].temps,
487  AtmolCollRateCoeff[ipSpecies][ipCollider].collrates[ipHi][ipLo],
488  AtmolCollRateCoeff[ipSpecies][ipCollider].ntemps,
489  ftemp);
490  }
491 
492  if(DEBUGSTATE)
493  printf("\nThe collision rate at temperature %f is %e",ftemp,ret_collrate);
494 
495  ASSERT( !isnan( ret_collrate ) );
496  return(ret_collrate);
497 
498 }
499 
500 /*Chianti*/
501 double Chianti_Upsilon(long ipSpecies, long ipCollider, long ipHi, long ipLo, double ftemp)
502 {
503  double fdeltae,fscalingparam,fkte,fxt,fsups,fups;
504  double *xs,*spl,*spl2;
505  int i,intxs,inttype,intsplinepts;
506 
507  DEBUG_ENTRY("Chianti_Upsilon()");
508 
509  if( AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline == NULL )
510  {
511  return 0.;
512  }
513 
514  intsplinepts = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].nSplinePts;
515  inttype = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].intTranType;
516  fdeltae = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].EnergyDiff;
517  fscalingparam = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].ScalingParam;
518 
519  fkte = ftemp/fdeltae/1.57888e5;
520 
521  /*Way the temperature is scaled*/
522  /*Burgess&Tully 1992:Paper gives only types 1 to 4*/
523  /*Found that the expressions were the same for 5 & 6 from the associated routine DESCALE_ALL*/
524  /*What about 7,8&9?*/
525  if(inttype ==1 || inttype == 4)
526  {
527  fxt = 1-(log(fscalingparam)/(log(fkte+fscalingparam)));
528  }
529  else if(inttype == 2 || inttype == 3||inttype == 5 || inttype == 6)
530  {
531  fxt = fkte/(fkte+fscalingparam);
532  }
533  else
534  TotalInsanity();
535 
536  /*Creating spline points array*/
537  if(intsplinepts == 5)
538  {
539  xs = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
540  spl = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
541  spl2 = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
542  for(intxs=0;intxs<5;intxs++)
543  {
544  xs[intxs] = 0.25*intxs;
545  spl[intxs] = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline[intxs];
546  if(DEBUGSTATE)
547  {
548  printf("The xs and spl values are %f and %f \n",xs[intxs],spl[intxs]);
549  getchar();
550  }
551  }
552  }
553  else if(intsplinepts == 9)
554  {
555  xs = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
556  spl = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
557  spl2 = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
558  for( intxs=0; intxs<9; intxs++ )
559  {
560  xs[intxs] = 0.125*intxs;
561  spl[intxs] = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline[intxs];
562  if(DEBUGSTATE)
563  {
564  printf("The xs and spl values are %f and %f \n",xs[intxs],spl[intxs]);
565  getchar();
566  }
567  }
568  }
569  else
570  {
571  TotalInsanity();
572  }
573 
574  /*Finding the second derivative*/
575  for( i=0; i<intsplinepts; i++)
576  {
577  spl2[i] = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].SplineSecDer[i];
578  }
579 
580  if(DEBUGSTATE)
581  {
582  printf("\n");
583  for(intxs=0;intxs<intsplinepts;intxs++)
584  {
585  printf("The %d value of 2nd derivative is %f \n",intxs+1,spl2[intxs]);
586  }
587  }
588 
589  /*Extracting out the value*/
590  splint(xs,spl,spl2,intsplinepts,fxt,&fsups);
591 
592  /*Finding upsilon*/
593  if(inttype == 1)
594  {
595  fups = fsups*log(fkte+exp(1.0));
596  }
597  else if(inttype == 2)
598  {
599  fups = fsups;
600  }
601  else if(inttype == 3)
602  {
603  fups = fsups/(fkte+1.0) ;
604  }
605  else if(inttype == 4)
606  {
607  fups = fsups*log(fkte+fscalingparam) ;
608  }
609  else if(inttype == 5)
610  {
611  fups = fsups/fkte ;
612  }
613  else if(inttype == 6)
614  {
615  fups = pow(10.0,fsups) ;
616  }
617  else
618  {
619  TotalInsanity();
620  }
621 
622  free( xs );
623  free( spl );
624  free( spl2 );
625 
626  ASSERT(fups>=0);
627  return(fups);
628 }

Generated for cloudy by doxygen 1.8.1.1