cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
opacity_addtotal.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 /*OpacityAddTotal derive total opacity for this position,
4  * called by ConvBase */
5 #include "cddefines.h"
6 #include "physconst.h"
7 #include "iso.h"
8 #include "grainvar.h"
9 #include "ca.h"
10 #include "rfield.h"
11 #include "oxy.h"
12 #include "mole.h"
13 #include "dense.h"
14 #include "atoms.h"
15 #include "conv.h"
16 #include "ionbal.h"
17 #include "trace.h"
18 #include "hmi.h"
19 #include "phycon.h"
20 #include "opacity.h"
21 
22 void OpacityAddTotal(void)
23 {
24  long int i,
25  ion,
26  ipHi,
27  ipISO,
28  ipop,
29  limit,
30  low,
31  nelem,
32  n;
33  double DepartCoefInv ,
34  fac,
35  sum;
36  realnum SaveOxygen1 ,
37  SaveCarbon1;
38 
39  DEBUG_ENTRY( "OpacityAddTotal()" );
40 
41  /* OpacityZero will zero out scattering and absorption opacities,
42  * and set OldOpacSave to opac to save it */
43  OpacityZero();
44 
45  /* free electron scattering opacity, Compton recoil energy loss */
46  /* >>chng 02 mar 26, make this loop only one over electrons alone */
47  /*for( i=0; i < (ionbal.ipCompRecoil[ipHYDROGEN][0] - 1); i++ )*/
48  for( i=0; i < rfield.nflux; i++ )
49  {
50  /* scattering part of total opacity */
52  dense.eden;
53  }
54 
55  /* opacity due to compton bound recoil ionization */
56  /* >>chng 01 dec 19, rewrite as loop over all elements */
57  /* >>chng 02 mar 26, add in ions */
58  for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
59  {
60  if( dense.lgElmtOn[nelem] )
61  {
62  for( ion=0; ion<nelem+1; ++ion )
63  {
64  realnum factor = dense.xIonDense[nelem][ion];
65  /*>>chng 05 nov 26, add molecular hydrogen assuming same
66  * as two free hydrogen atoms - hence mult density by two
67  *>>KEYWORD H2 bound Compton ionization */
68  if( nelem==ipHYDROGEN )
69  factor += hmi.H2_total*2.f;
70  if( factor > 0. )
71  {
72  // loop_min and loop_max are needed to work around a bug in icc 10.0
73  long loop_min = ionbal.ipCompRecoil[nelem][ion]-1;
74  long loop_max = rfield.nflux;
75  /* ionbal.nCompRecoilElec number of electrons in valence shell
76  * that can compton recoil ionize */
77  factor *= ionbal.nCompRecoilElec[nelem-ion];
78  for( i=loop_min; i < loop_max; i++ )
79  {
80  /* add in bound hydrogen electron scattering, treated as absorption */
81  opac.opacity_abs[i] += opac.OpacStack[i-1+opac.iopcom]*factor;
82  }
83  }
84  }
85  }
86  }
87 
88  /* opacity due to pair production - does not matter what form these
89  * elements are in */
93 
94  /* free-free free free brems emission for all ions */
95 
96  /* >>chng 02 jul 21, use full set of ions and gaunt factor */
97  /* ipEnergyBremsThin is index to energy where gas goes optically thin to brems,
98  * so this loop is over optically thick frequencies */
99  static double *TotBremsAllIons;
100  static bool lgFirstTime=true;
101  double BremsThisEner,bfac, bhfac,sumion[LIMELM+1];
102  long int ion_lo , ion_hi;
103 
104  if( lgFirstTime )
105  {
106  /* rfield.nupper will not change in one coreload, so just malloc this once */
107  TotBremsAllIons = (double *)MALLOC((unsigned long)rfield.nupper*sizeof(double));
108  lgFirstTime = false;
109  }
110  bfac = (dense.eden/1e20)/phycon.sqrte/1e10;
111  /* gaunt factors depend only on photon energy and ion charge, so do
112  * sum of ions here before entering into loop over photon energy */
113  sumion[0] = 0.;
114  for(ion=1; ion<=LIMELM; ++ion )
115  {
116  sumion[ion] = 0.;
117  for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
118  {
119  if( dense.lgElmtOn[nelem] && ion<=nelem+1 )
120  {
121  sumion[ion] += dense.xIonDense[nelem][ion];
122  }
123  }
124  /* now include the charge, density, and temperature */
125  sumion[ion] *= POW2((double)ion)*bfac;
126  }
127  /* now find lowest and highest ion we need to consider - following loop is over
128  * full continuum and eats time
129  * >>chng 05 oct 20, following had tests on ion being within bounds, bug,
130  * changed to ion_lo and ion_hi */
131  ion_lo = 1;
132  while( sumion[ion_lo]==0 && ion_lo<LIMELM-1 )
133  ++ion_lo;
134  ion_hi = LIMELM;
135  while( sumion[ion_hi]==0 && ion_hi>0 )
136  --ion_hi;
137 
138  bhfac = bfac*(dense.xIonDense[ipHYDROGEN][1]+hmi.Hmolec[ipMH2p]+hmi.Hmolec[ipMH3p]);
139  for( i=0; i < rfield.nflux; i++ )
140  {
141  /* do hydrogen first, before main loop since want to add on H- brems */
142  nelem = ipHYDROGEN;
143  ion = 1;/* >>chng 02 nov 07, add charged ions as H+ */
144 
145  /* for case of hydrogen, add H- brems - OpacStack contains the ratio
146  * of the H- to H brems cross section - multiply by this and the fraction in ground */
147  BremsThisEner = bhfac * rfield.gff[ion][i]*
149  /*TotBremsAllIons[i] = BremsThisIon;*/
150  /* there is at least one standard test (grainlte) which has ZERO ionization -
151  * this assert will pass that test (the ==0) and also the usual case */
152  /* ASSERT( (dense.xIonDense[nelem][ion]==0.) || (TotBremsAllIons[i] > 0.) );*/
153 
154  /* chng 02 may 16, by Ryan...do all brems for all ions in one fell swoop,
155  * using gaunt factors from rfield.gff. */
156  /* >>chng 05 jul 11 reorganize loop for speed */
157  for(ion=ion_lo; ion<=ion_hi; ++ion )
158  {
159  BremsThisEner += sumion[ion]*rfield.gff[ion][i];
160  }
161  TotBremsAllIons[i] = BremsThisEner;
162 
163  /* >>>chng 05 jul 12, bring these two loops together */
164  if( rfield.ContBoltz[i] < 0.995 )
165  {
166  TotBremsAllIons[i] *= opac.OpacStack[i-1+opac.ipBrems]*
167  (1. - rfield.ContBoltz[i]);
168  }
169  else
170  {
171  TotBremsAllIons[i] *= opac.OpacStack[i-1+opac.ipBrems]*
173  }
174  opac.FreeFreeOpacity[i] = TotBremsAllIons[i];
175  /* >>chng 02 jul 23, from >0 to >=0, some models have no ionization,
176  * like grainlte.in */
177  /*ASSERT( (opac.FreeFreeOpacity[i] > 0.) || (dense.xIonDense[ipHYDROGEN][1] == 0.) );*/
179  }
180 
181 
182  /* H minus absorption, with correction for stimulated emission */
183  if( hmi.hmidep > SMALLFLOAT )
184  {
185  DepartCoefInv = 1./hmi.hmidep;
186  }
187  else
188  {
189  /* the hmidep departure coef can become vastly small in totally
190  * neutral gas (no electrons) */
191  DepartCoefInv = 1.;
192  }
193 
194  limit = iso.ipIsoLevNIonCon[ipHE_LIKE][ipHELIUM][0];
195  if(limit > rfield.nflux)
196  limit = rfield.nflux;
197 
198  for( i=hmi.iphmin-1; i < limit; i++ )
199  {
200  double factor;
201  factor = 1. - rfield.ContBoltz[i]*DepartCoefInv;
202  if(factor > 0)
204  hmi.Hmolec[ipMHm]*factor;
205  }
206 
207  /* H2P h2plus bound free opacity */
208  limit = opac.ih2pnt[1];
209  if(limit > rfield.nflux)
210  limit = rfield.nflux;
211  for( i=opac.ih2pnt[0]-1; i < limit; i++ )
212  {
214  }
215 
216  /* get total population of hydrogen ground to do Rayleigh scattering */
217  if( dense.xIonDense[ipHYDROGEN][1] <= 0. )
218  {
219  fac = dense.xIonDense[ipHYDROGEN][0];
220  }
221  else
222  {
224  }
225 
226  /* Ly a damp wing opac (Rayleigh scattering) */
228  if(limit > rfield.nflux)
229  limit = rfield.nflux;
230  for( i=0; i < limit; i++ )
231  {
232  opac.opacity_sct[i] += (fac*opac.OpacStack[i-1+opac.ipRayScat]);
233  }
234 
235  /* remember largest correction for stimulated emission */
236  if( iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH1s] > 1e-30 && !conv.lgSearch )
237  {
238  realnum factor;
240  if(opac.stimax[0] < factor)
241  opac.stimax[0] = factor;
242  }
243 
244  if( iso.DepartCoef[ipH_LIKE][ipHYDROGEN][ipH2p] > 1e-30 && !conv.lgSearch )
245  {
246  realnum factor;
248  if(opac.stimax[1] < factor)
249  opac.stimax[1] = factor;
250  }
251 
252 # if 0
253  /* check whether hydrogen or Helium singlets mased, if not in search mode */
254  if( !conv.lgSearch )
255  {
256  if( Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].PopOpc < 0. )
257  {
258  hydro.lgHLyaMased = true;
259  }
260  }
261 # endif
262 
263  /* >>chng 05 nov 25, use Yan et al. H2 photo cs
264  * following reference gives cross section for all energies
265  * >>refer H2 photo cs Yan, M., Sadeghpour, H.R., & Dalgarno, A., 1998, ApJ, 496, 1044
266  * Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914 */
267  /* >>chng 02 jan 16, approximate inclusion of H_2 photoelectric opacity
268  * include H_2 in total photoelectric opacity as twice H0 cs */
269  /* set lower and upper limits to this range */
270  /*>>KEYWORD H2 photoionization opacity */
271  low = opac.ipH2_photo_thresh;
272  ipHi = rfield.nupper;
274  /* OpacityAdd1Subshell just returns for static opacities if opac.lgRedoStatic not set*/
275  /* >>chng 05 nov 27, change on nov 25 had left 2*density from H0, so
276  * twice the H2 density was used - * also changed to static opacity
277  * this assumes that all v,J levels contribute the same opacity */
278  OpacityAdd1Subshell( ipop , low , ipHi , hmi.H2_total , 's' );
279 
280  /*>>KEYWORD CO photoionization opacity */
281  /* include photoionization of CO - assume C and O in CO each have
282  * same photo cs as atom - this should only be significant in highly
283  * shielded regions where only very hard photons penetrate
284  * also H2O condensed onto grain surfaces - very important deep in cloud */
285  SaveOxygen1 = dense.xIonDense[ipOXYGEN][0];
286  SaveCarbon1 = dense.xIonDense[ipCARBON][0];
287  dense.xIonDense[ipOXYGEN][0] += findspecies("CO")->hevmol + findspecies("H2Ogrn")->hevmol;
288  dense.xIonDense[ipCARBON][0] += findspecies("CO")->hevmol;
289 
290  /* following loop adds standard opacities for first 30 elements
291  * most heavy element opacity added here */
292  for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
293  {
294  /* this element may be turned off */
295  if( dense.lgElmtOn[nelem] )
296  {
297  OpacityAdd1Element(nelem);
298  }
299  }
300 
301  /* now reset the abundances */
302  dense.xIonDense[ipOXYGEN][0] = SaveOxygen1;
303  dense.xIonDense[ipCARBON][0] = SaveCarbon1;
304 
305  /* following are opacities due to specific excited levels */
306 
307  /* nitrogen opacity
308  * excited level of N+ */
310  dense.xIonDense[ipNITROGEN][0]*atoms.p2nit , 'v' );
311 
312  /* oxygen opacity
313  * excited level of Oo */
315  dense.xIonDense[ipOXYGEN][0]*oxy.poiexc,'v');
316 
317  /* O2+ excited states */
319  dense.xIonDense[ipOXYGEN][2]*oxy.poiii2,'v');
320 
322  dense.xIonDense[ipOXYGEN][2]*oxy.poiii3,'v');
323 
324  /* magnesium opacity
325  * excited level of Mg+ */
328 
329  /* calcium opacity
330  * photoionization of excited levels of Ca+ */
332  ca.popca2ex,'v');
333 
334  /*******************************************************************
335  *
336  * complete evaluation of total opacity by adding in the static part and grains
337  *
338  *******************************************************************/
339 
340  /* this loop defines the variable iso.ConOpacRatio[ipH_LIKE][nelem][n],
341  * the ratio of not H to Hydrogen opacity. for grain free environments
342  * at low densities this is nearly zero. The correction includes
343  * stimulated emission correction */
344  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
345  {
346  for( nelem=ipISO; nelem < LIMELM; nelem++ )
347  {
348  /* this element may be turned off */
349  if( dense.lgElmtOn[nelem] )
350  {
351  /* this branch is for startup only */
352  if( nzone < 1 )
353  {
354  /* >>chng 06 aug 17, should go to numLevels_local instead of _max */
355  for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ )
356  {
357  if(iso.ipIsoLevNIonCon[ipISO][nelem][n] < rfield.nflux )
358  {
359  /*>>chng 04 dec 12, had tested against 1e-30, set to zero if not
360  * greater than this - caused oscillations as opacity fell below
361  * and around this value - change to greater than 0 */
362  /*if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 1e-30 )*/
363  if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 0. )
364  {
365  /* >>chng 02 may 06, use general form of threshold cs */
366  /*double t1 = atmdat_sth(n)/(POW2(nelem+1.-ipISO));*/
367  long int ip = iso.ipIsoLevNIonCon[ipISO][nelem][n];
368  double t2 = csphot(
369  ip ,
370  ip ,
371  iso.ipOpac[ipISO][nelem][n] );
372 
373  iso.ConOpacRatio[ipISO][nelem][n] = 1.f-(realnum)(
374  (StatesElem[ipISO][nelem][n].Pop*dense.xIonDense[nelem][nelem+1-ipISO]*t2)/
375  opac.opacity_abs[ip-1]);
376  }
377  else
378  {
379  iso.ConOpacRatio[ipISO][nelem][n] = 0.;
380  }
381  }
382  }
383  }
384  /* end branch for startup only, start branch for all zones including startup */
385  /* >>chng 06 aug 17, should go to numLevels_local instead of _max */
386  for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ )
387  {
388  /* ratios of other to total opacity for continua of all atoms done with iso model */
389  if(iso.ipIsoLevNIonCon[ipISO][nelem][n] < rfield.nflux )
390  {
391  /*>>chng 04 dec 12, had tested against 1e-30, set to zero if not
392  * greater than this - caused oscillations as opacity fell below
393  * and around this value - change to greater than 0 */
394  /*if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 1e-30 )*/
395  if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 0. )
396  {
397  /* first get departure coef */
398  if( iso.DepartCoef[ipISO][nelem][n] > 1e-30 && (!conv.lgSearch ) )
399  {
400  /* this is the usual case, use inverse of departure coef */
401  fac = 1./iso.DepartCoef[ipISO][nelem][n];
402  }
403  else if( conv.lgSearch )
404  {
405  /* do not make correction for stim emission during search
406  * for initial temperature solution, since trys are very non-equil */
407  fac = 0.;
408  }
409  else
410  {
411  fac = 1.;
412  }
413 
416  /* now get opaicty ratio with correction for stimulated emission */
417  /*>>chng 04 dec 12, had tested against 1e-30, set to zero if not
418  * greater than this - caused oscillations as opacity fell below
419  * and around this value - change to greater than 0 */
420  /*if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 1e-30 )*/
421  if( opac.opacity_abs[iso.ipIsoLevNIonCon[ipISO][nelem][n]-1] > 0. )
422  {
423  /* >>chng 02 may 06, use general form of threshold cs */
424  long int ip = iso.ipIsoLevNIonCon[ipISO][nelem][n];
425 
426  double t2 = csphot(
427  ip ,
428  ip ,
429  iso.ipOpac[ipISO][nelem][n] );
430 
431  double opacity_this_species =
432  StatesElem[ipISO][nelem][n].Pop*dense.xIonDense[nelem][nelem+1-ipISO]*t2*
433  (1. - fac*rfield.ContBoltz[ip-1]);
434 
435  double opacity_fraction = 1. - opacity_this_species / opac.opacity_abs[ip-1];
436  if(opacity_fraction < 0)
437  opacity_fraction = 0.;
438 
439  /* use mean of old and new ratios */
440  iso.ConOpacRatio[ipISO][nelem][n] = (realnum)(
441  iso.ConOpacRatio[ipISO][nelem][n]* 0.75 + 0.25*opacity_fraction );
442 
443  if(iso.ConOpacRatio[ipISO][nelem][n] < 0.)
444  iso.ConOpacRatio[ipISO][nelem][n] = 0.;
445  }
446  else
447  {
448  iso.ConOpacRatio[ipISO][nelem][n] = 0.;
449  }
450  }
451  else
452  {
453  iso.ConOpacRatio[ipISO][nelem][n] = 0.;
454  }
455  }
456  else
457  {
458  iso.ConOpacRatio[ipISO][nelem][n] = 0.;
459  }
460  }
461  }
462  }
463  }
464 
465  /* add dust grain opacity if dust present */
466  if( gv.lgDustOn )
467  {
468  /* generate current grain opacities since may be function of depth */
469  /* >>chng 01 may 11, removed code to update grain opacities, already done by GrainChargeTemp */
470  for( i=0; i < rfield.nflux; i++ )
471  {
472  /* units cm-1 */
475  }
476  }
477 
478  /* check that opacity is sane */
479  for( i=0; i < rfield.nflux; i++ )
480  {
481  /* OpacStatic was zeroed in OpacityZero, incremented in opacityadd1subshell */
482  opac.opacity_abs[i] += opac.OpacStatic[i];
483  /* make sure that opacity is positive */
484  /*ASSERT( opac.opacity_abs[i] > 0. );*/
485  }
486 
487  /* compute gas albedo here */
488  for( i=0; i < rfield.nflux; i++ )
489  {
490  opac.albedo[i] = opac.opacity_sct[i]/
491  (opac.opacity_sct[i] + opac.opacity_abs[i]);
492  }
493 
494  /* during search phase set old opacity array to current value */
495  if( conv.lgSearch )
496  OpacityZeroOld();
497 
498  if( trace.lgTrace )
499  {
500  fprintf( ioQQQ, " OpacityAddTotal returns; grd rec eff (opac) for Hn=1,4%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e%10.2e HeI,II:%10.2e%10.2e\n",
512  iso.ConOpacRatio[ipH_LIKE][ipHELIUM][ipH1s] );
513  }
514  {
515  /* following should be set true to print strongest ots contributors */
516  /*@-redef@*/
517  enum {DEBUG_LOC=false};
518  /*@+redef@*/
519  if( DEBUG_LOC && (nzone>=378) )
520  {
521  if( nzone > 380 )
522  cdEXIT( EXIT_FAILURE );
523  for( i=0; i<rfield.nflux; ++i )
524  {
525  fprintf(ioQQQ,"rtotsbugggg\t%li\t%.3e\t%.3e\t%.3e\t%.3e\n",
527  rfield.anu[i],
528  opac.opacity_abs[i],
529  rfield.otscon[i],
530  rfield.otslin[i]);
531  }
532  }
533  }
534  return;
535 }

Generated for cloudy by doxygen 1.8.1.1