cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ion_carbo.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 /*IonCarbo compute ionization balance for carbon */
4 #include "cddefines.h"
5 #include "thermal.h"
6 #include "carb.h"
7 #include "trace.h"
8 #include "dense.h"
9 #include "phycon.h"
10 #include "hmi.h"
11 #include "mole.h"
12 #include "rfield.h"
13 #include "punch.h"
14 #include "ionbal.h"
15 
16 void IonCarbo(void)
17 {
18  const int NDIM = ipCARBON+1;
19 
20  static const double dicoef[2][NDIM] = {
21  {2.54e-3,6.15e-3,1.62e-3,4.78e-2,3.22e-2,0.}, {4.42e-2,5.88e-2,0.343,0.362,0.315,0.}
22  };
23  static const double dite[2][NDIM] = {
24  {1.57e5,1.41e5,8.19e4,3.44e6,4.06e6,0.}, {3.74e5,1.41e5,1.59e5,5.87e5,8.31e5,0.}
25  };
26  static const double ditcrt[NDIM] = {1.2e4,1.2e4,1.1e4,4.4e5,7.0e5,1e20};
27  /* >>refer C DR Nussbaumer H. & Storey, P 1983, A&A, 126, 75 */
28  static const double aa[NDIM] = {.0108,1.8267,2.3196,0.,0.,0.};
29  static const double bb[NDIM] = {-0.1075,4.1012,10.7328,0.,0.,0.};
30  static const double cc[NDIM] = {.2810,4.8443,6.8830,0.,0.,0.};
31  static const double dd[NDIM] = {-0.0193,.2261,-0.1824,0.,0.,0.};
32  /* NB C+ - C0 recom numerically unstable below 2000K, use est instead */
33  static const double ff[NDIM] = {0.000,0.5960,0.4101,0.1,0.1,0.};
34 
35  double save_rec;
36 
37  DEBUG_ENTRY( "IonCarbo()" );
38 
39  if( trace.lgTrace && trace.lgCarBug )
40  {
41  fprintf( ioQQQ, " IonCarbo called.\n" );
42  }
43 
44  if( !dense.lgElmtOn[ipCARBON] )
45  {
46  carb.p1909 = 0.;
47  carb.p2326 = 0.;
48  thermal.heating[ipCARBON][9] = 0.;
49  return;
50  }
51 
52  /* zero out ionization balance arrays */
54 
55  ion_photo(ipCARBON,false);
56 
57  /* >>chng 05 aug 10, add Leiden hack here to get same C0 photo rate
58  * as UMIST - negates difference in grain opacities */
59  if(!co.lgUMISTrates)
60  {
61  int nelem=ipCARBON , ion=0 , ns=2;
62  ionbal.PhotoRate_Shell[nelem][ion][ns][0] =
63  (HMRATE((1e-10)*3.0,0,0)*(hmi.UV_Cont_rel2_Habing_TH85_face*
64  exp(-(3.0*rfield.extin_mag_V_point))/1.66));
65  /* heating rates */
66  ionbal.PhotoRate_Shell[nelem][ion][ns][1] = 0.;
67  ionbal.PhotoRate_Shell[nelem][ion][ns][2] = 0.;
68  /*fprintf(ioQQQ,"DEBUG hack %li\tav\t%.2e\tgo\t%.2e\trate\t%.2e\tabun\t%.2e\n",
69  nzone ,
70  rfield.extin_mag_V_point,
71  hmi.UV_Cont_rel2_Habing_TH85_face,
72  ionbal.PhotoRate_Shell[nelem][ion][ns][0] ,
73  dense.xIonDense[ipCARBON][1]/SDIV(dense.xIonDense[ipCARBON][0]) );*/
74  }
75 
76  /* find collisional ionization rates */
78 
79  /* get recombination coefficients */
80  /*lint -e64 double = pointer */
81  ion_recomb(false,(const double*)dicoef,(const double*)dite,ditcrt,aa,bb,cc,dd,ff,ipCARBON);
82  /*lint +e64 double = pointer */
83 
84  /* Photoproduction of 3P */
86 
87  /* photoproduction of C II] 2326 */
89 
90  /* >>chng 03 sep 29, synch up ion and co solvers.
91  * the co solver will have a different C/C+ balance that
92  * is strongly affected by the chemistry if we are in neutral gas
93  * in this case use co solver's C/C+ balance */
94  save_rec = ionbal.RateRecomTot[ipCARBON][0];
95 
96  bool lgDEBUG=false;
97  /*if( iteration>1 )*/
98  /*if( nzone > 258 )
99  lgDEBUG = true;
100  else
101  lgDEBUG = false;*/
102  ion_solver(ipCARBON,lgDEBUG);
103 
104  /* reset the var we just hosed */
105  if( save_rec > 0. )
106  ionbal.RateRecomTot[ipCARBON][0] = save_rec;
107 
108  /* rate will be cm-3 s-1, into triplets only
109  * >>chng 96 nov 22, rid of carb() */
110  carb.p1909 *= dense.xIonDense[ipCARBON][1]*0.62;
111  carb.p2326 *= dense.xIonDense[ipCARBON][0]*0.1;
112 
113  if( trace.lgTrace )
114  {
115  fprintf( ioQQQ, " IonCarbo returns; fracs=" );
116  for( int i=0; i < 7; i++ )
117  {
118  fprintf( ioQQQ, " %10.3e", dense.xIonDense[ipCARBON][i]/
120  }
121  fprintf( ioQQQ, "\n" );
122  }
123 
124  /* We have already checked the validity of cont_gaunt_calc in sanitycheck.c.
125  * Now we check to see if the InterpolateGff routine also works correctly. */
126  enum {AGN=false};
127  /* if set true there will be two problems at very low temps */
128  if( AGN )
129  {
130  phycon.te=10.;
131  /* this tells ion_recomb to punch data */
132  punch.lgioRecom = true;
133  /* this is where it will be punched */
134  punch.ioRecom = ioQQQ;
135  while( phycon.te<1e7 )
136  {
137  /* get recombination coefficients */
138  /*lint -e64 double = pointer */
139  ion_recomb(false,(double*)dicoef,(double*)dite,ditcrt,aa,bb,cc,
140  dd,ff,ipCARBON);
141  /*lint +e64 double = pointer */
142  TempChange(phycon.te *2.f , true);
143  }
144  /* bail out */
145  cdEXIT(EXIT_SUCCESS);
146  }
147  return;
148 }

Generated for cloudy by doxygen 1.8.3.1