cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
iso_continuum_lower.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 /*iso_continuum_lower - limit max prin. quan. no. due to continuum lowering processes */
4 #include "cddefines.h"
5 #include "phycon.h"
6 #include "dense.h"
7 #include "iso.h"
8 #include "hydrogenic.h"
9 #include "trace.h"
10 
11 void iso_continuum_lower( long ipISO, long nelem )
12 {
13  double a;
14  long int np, nd, ns, nc;
15  long eff_charge;
16 
17  /* size of rate matrices will be defined according to the n calculated here */
18 
19  ASSERT( nelem < LIMELM );
20  /* this may change at a future date. */
21  ASSERT( ipISO <= 1 );
22 
23  eff_charge = nelem + 1 - ipISO;
24 
25  /* Particle packing - the density here should be density of all nuclei in the plasma */
26  /* This one is just nuclear charge, which is independent of iso, and always nelem+1. */
27  a = sqrt( 1.8887E8 * (nelem+1.) / pow((double)dense.xNucleiTotal, 0.333) );
28  ASSERT( a > 0. );
29  ASSERT( a < 1.e15 );
30  if( a > (double)iso.n_HighestResolved_max[ipISO][nelem]+(double)iso.nCollapsed_max[ipISO][nelem] )
31  {
32  np = iso.n_HighestResolved_max[ipISO][nelem]+iso.nCollapsed_max[ipISO][nelem] + 1;
33  }
34  else
35  np = (long)a;
36 
37  /* Debye shielding - the density here is electron density */
38  /* This one depends on effective charge. */
39  a = 2.6E7 * eff_charge * eff_charge * pow( phycon.te/dense.eden, 0.25);
40  ASSERT( a > 0. );
41  ASSERT( a < 1.e15 );
42  if( a > (double)iso.n_HighestResolved_max[ipISO][nelem]+(double)iso.nCollapsed_max[ipISO][nelem] )
43  {
44  nd = iso.n_HighestResolved_max[ipISO][nelem]+iso.nCollapsed_max[ipISO][nelem] + 1;
45  }
46  else
47  nd = (long)a;
48 
49  /* Stark broadening - this should be the density of singly charged ions,
50  * both positive and negative. The sum of protons, electrons, and HeII should be
51  * good enough. */
52  /* This one depends on effective charge. */
53  a = 3171. * pow( (double)eff_charge, 0.8 ) * pow( dense.eden + (double)dense.xIonDense[ipHYDROGEN][1]
54  + (double)dense.xIonDense[ipHELIUM][1], -0.1333);
55  ASSERT( a > 0. );
56  ASSERT( a < 1.e15 );
57  if( a > (double)iso.n_HighestResolved_max[ipISO][nelem]+(double)iso.nCollapsed_max[ipISO][nelem] )
58  {
59  ns = iso.n_HighestResolved_max[ipISO][nelem]+iso.nCollapsed_max[ipISO][nelem] + 1;
60  }
61  else
62  ns = (long)a;
63 
64  ASSERT( np > 3 );
65  ASSERT( nd > 3 );
66  ASSERT( ns > 3 );
67 
68  nc = MIN3(np, nd, ns);
69 
70  /* I assert greater than three because the code depends upon having at least up to n=3, and
71  * because it would take EXTREMELY dense conditions to lower the continuum that much, and
72  * the code would probably get a wrong answer then anyway. */
73  ASSERT( nc > 3 );
74 
75  if( nc < iso.n_HighestResolved_max[ipISO][nelem])
76  {
77  iso.lgLevelsLowered[ipISO][nelem] = true;
78  iso.lgLevelsEverLowered[ipISO][nelem] = true;
79  hydro.lgReevalRecom = true;
80  iso.n_HighestResolved_local[ipISO][nelem] = nc;
81  iso.nCollapsed_local[ipISO][nelem] = 0;
82  iso.numLevels_local[ipISO][nelem] = iso_get_total_num_levels( ipISO, nc, 0 );
83  }
84  /* Here is the case where the critical n lies among the one or more collapsed levels */
85  /* we just get rid of any that are too high. */
86  else if( nc <= iso.n_HighestResolved_max[ipISO][nelem] + iso.nCollapsed_max[ipISO][nelem] )
87  {
88  iso.lgLevelsLowered[ipISO][nelem] = true;
89  iso.lgLevelsEverLowered[ipISO][nelem] = true;
90  hydro.lgReevalRecom = true;
91  iso.n_HighestResolved_local[ipISO][nelem] = iso.n_HighestResolved_max[ipISO][nelem];
92  iso.nCollapsed_local[ipISO][nelem] = nc - iso.n_HighestResolved_local[ipISO][nelem];
93  iso.numLevels_local[ipISO][nelem] =
94  iso_get_total_num_levels( ipISO, iso.n_HighestResolved_max[ipISO][nelem], iso.nCollapsed_local[ipISO][nelem] );
95  }
96  /* This is usually where control will flow, because in most conditions the continuum will not be lowered.
97  * Nothing changes in this case. */
98  else
99  {
100  iso.numLevels_local[ipISO][nelem] = iso.numLevels_max[ipISO][nelem];
101  iso.nCollapsed_local[ipISO][nelem] = iso.nCollapsed_max[ipISO][nelem];
102  iso.n_HighestResolved_local[ipISO][nelem] = iso.n_HighestResolved_max[ipISO][nelem];
103 
104  /* if levels were lowered on last pass but are not now, must reeval */
105  if( iso.lgLevelsLowered[ipISO][nelem] )
106  {
107  hydro.lgReevalRecom = true;
108  }
109  else
110  {
111  hydro.lgReevalRecom = false;
112  }
113 
114  iso.lgLevelsLowered[ipISO][nelem] = false;
115  }
116 
117  /* None of these can be greater than that which was originally malloc'd. */
118  ASSERT( iso.numLevels_local[ipISO][nelem] <= iso.numLevels_max[ipISO][nelem] );
119  ASSERT( iso.nCollapsed_local[ipISO][nelem] <= iso.nCollapsed_max[ipISO][nelem] );
120  ASSERT( iso.n_HighestResolved_local[ipISO][nelem] <= iso.n_HighestResolved_max[ipISO][nelem] );
121 
122  /* don't print levels in continuum. */
123  iso.numPrintLevels[ipISO][nelem] = MIN2( iso.numPrintLevels[ipISO][nelem], iso.numLevels_local[ipISO][nelem] );
124 
125  /* Lyman lines can not be greater than original malloc or critical pqn. */
126  iso.nLyman[ipISO] = MIN2( nc, iso.nLyman[ipISO]);
127 
128  if( trace.lgTrace )
129  {
130  fprintf( ioQQQ," iso_continuum_lower: ipISO %li nelem %li nc %li numLevels %li nCollapsed %li n_HighestResolved %li \n",
131  ipISO,
132  nelem,
133  nc,
134  iso.numLevels_local[ipISO][nelem],
135  iso.nCollapsed_local[ipISO][nelem],
136  iso.n_HighestResolved_local[ipISO][nelem]
137  );
138  }
139 
140  return;
141 }

Generated for cloudy by doxygen 1.8.1.1