cloudy  trunk
 All Data Structures Namespaces 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-2022 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 "conv.h"
8 #include "iso.h"
9 #include "trace.h"
10 
11 /* iso_continuum_lower - implement continuum lowering as per section 3.1 of
12  * >>refer ionization Bautista, M.A., & Kallman, T.R., 2000, ApJ, 544, 581
13  * http://adsabs.harvard.edu/abs/2000ApJ...544..581B
14  */
15 void iso_continuum_lower( long ipISO, long nelem )
16 {
17  long nc;
18  t_iso_sp* sp = &iso_sp[ipISO][nelem];
19 
20  /* size of rate matrices will be defined according to the n calculated here */
21 
23  ASSERT( nelem < LIMELM );
24  /* this may change at a future date. */
25  ASSERT( ipISO <= 1 );
26 
27  long eff_charge = nelem + 1 - ipISO;
28 
29  /* Particle packing - the density here should be density of all nuclei in the plasma */
30  /* This one is just nuclear charge, which is independent of iso, and always nelem+1. */
31  double np = sqrt( 1.8887E8 * (nelem+1.) / cbrt((double)dense.xNucleiTotal) );
32  np = floor(np);
33  ASSERT( np > 0. );
34 
35  /* Debye shielding - the density here is electron density */
36  /* This one depends on effective charge. */
37  double nd = 2.6E7 * eff_charge * eff_charge * powpq( phycon.te/dense.eden, 1, 4 );
38  nd = floor(nd);
39  ASSERT( nd > 0. );
40 
41  /* Stark broadening - this should be the density of singly charged ions,
42  * both positive and negative. The sum of protons, electrons, and HeII should be
43  * good enough. */
44  /* This one depends on effective charge. */
45  double ns = 3171. * pow( (double)eff_charge, 0.8 ) * pow( dense.eden + (double)dense.xIonDense[ipHYDROGEN][1]
46  + (double)dense.xIonDense[ipHELIUM][1], -0.1333);
47  ns = floor(ns);
48  ASSERT( ns > 0. );
49 
50  {
51  double nc_doub = MIN3(np, nd, ns);
52  if( nc_doub > INT16_MAX )
53  nc = INT16_MAX;
54  else
55  nc = (long)nc_doub;
56  }
57 
58  /* Don't allow continuum to be lowered below n=3. */
59  nc = MAX2( nc, 3 );
60 
61  if( nc <= sp->n_HighestResolved_max)
62  {
63  sp->lgLevelsLowered = true;
64  sp->lgLevelsEverLowered = true;
65  sp->lgMustReeval = true;
66  sp->n_HighestResolved_local = nc;
67  sp->nCollapsed_local = 0;
68  sp->numLevels_local = iso_get_total_num_levels( ipISO, nc, 0 );
69  }
70  /* Here is the case where the critical n lies among the one or more collapsed levels */
71  /* we just get rid of any that are too high. */
72  else if( nc <= sp->n_HighestResolved_max + sp->nCollapsed_max )
73  {
74  sp->lgLevelsLowered = true;
75  sp->lgLevelsEverLowered = true;
76  sp->lgMustReeval = true;
79  sp->numLevels_local =
81  }
82  /* This is usually where control will flow, because in most conditions the continuum will not be lowered.
83  * Nothing changes in this case. */
84  else
85  {
89 
90  /* if levels were lowered on last pass but are not now, must reeval */
91  if( sp->lgLevelsLowered )
92  {
93  sp->lgMustReeval = true;
94  }
95  else
96  {
97  sp->lgMustReeval = false;
98  }
99 
100  sp->lgLevelsLowered = false;
101  }
102 
103  if( !conv.nTotalIoniz )
104  sp->lgMustReeval = true;
105 
106  /* None of these can be greater than that which was originally malloc'd. */
107  ASSERT( sp->numLevels_local <= sp->numLevels_max );
108  ASSERT( sp->nCollapsed_local <= sp->nCollapsed_max );
110 
111  /* Lyman lines can not be greater than original malloc or critical pqn. */
112  iso_ctrl.nLyman[ipISO] = MIN2( nc, iso_ctrl.nLyman_max[ipISO]);
113 
114  // zero out cooling and heating terms involving unused levels
115  for( long ipHi=sp->numLevels_local; ipHi < sp->numLevels_max; ++ipHi )
116  {
117  for( long ipLo=0; ipLo < ipHi; ++ipLo )
118  CollisionZero( sp->trans(ipHi,ipLo).Coll() );
119  }
120 
121  if( trace.lgTrace && (trace.lgHBug||trace.lgHeBug) )
122  {
123  fprintf( ioQQQ," iso_continuum_lower: ipISO %li nelem %li nc %li (np:%.0f,nd:%.0f,ns:%.0f) numLevels %li nCollapsed %li n_HighestResolved %li \n",
124  ipISO,
125  nelem,
126  nc,
127  np, nd, ns,
128  sp->numLevels_local,
129  sp->nCollapsed_local,
131  );
132  }
133 
134  return;
135 }
void iso_continuum_lower(long ipISO, long nelem)
bool lgHeBug
Definition: trace.h:79
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
realnum xNucleiTotal
Definition: dense.h:114
long int nCollapsed_max
Definition: iso.h:518
t_conv conv
Definition: conv.cpp:5
long iso_get_total_num_levels(long ipISO, long nmaxResolved, long numCollapsed)
t_phycon phycon
Definition: phycon.cpp:6
FILE * ioQQQ
Definition: cddefines.cpp:7
#define MIN2(a, b)
Definition: cddefines.h:803
t_dense dense
Definition: global.cpp:15
const double MAX_DENSITY
Definition: cddefines.h:319
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
long int n_HighestResolved_local
Definition: iso.h:538
t_trace trace
Definition: trace.cpp:5
long int n_HighestResolved_max
Definition: iso.h:536
long int nLyman_max[NISO]
Definition: iso.h:352
bool lgLevelsLowered
Definition: iso.h:505
long int nLyman[NISO]
Definition: iso.h:352
#define INT16_MAX
Definition: cpu.h:22
bool lgTrace
Definition: trace.h:12
bool lgLevelsEverLowered
Definition: iso.h:509
long int nTotalIoniz
Definition: conv.h:159
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
Definition: iso.h:470
#define ASSERT(exp)
Definition: cddefines.h:613
bool lgHBug
Definition: trace.h:82
const int LIMELM
Definition: cddefines.h:308
CollisionProxy Coll() const
Definition: transition.h:463
double powpq(double x, int p, int q)
Definition: service.cpp:630
const int ipHELIUM
Definition: cddefines.h:350
#define MIN3(a, b, c)
Definition: cddefines.h:808
double eden
Definition: dense.h:201
#define MAX2(a, b)
Definition: cddefines.h:824
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
long int numLevels_max
Definition: iso.h:524
void CollisionZero(const CollisionProxy &t)
Definition: collision.cpp:88
double te
Definition: phycon.h:21
const int ipHYDROGEN
Definition: cddefines.h:349
long int numLevels_local
Definition: iso.h:529
long int nCollapsed_local
Definition: iso.h:519
bool lgMustReeval
Definition: iso.h:512