00001 /* This file is part of Cloudy and is copyright (C)1978-2010 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*OpacityAdd1Element enter total photo cross section for all subshells into opacity array */ 00004 #include "cddefines.h" 00005 #include "iso.h" 00006 #include "rfield.h" 00007 #include "dense.h" 00008 #include "heavy.h" 00009 #include "opacity.h" 00010 #include "taulines.h" 00011 00012 void OpacityAdd1Element( 00013 /* nelem is 0 for H, 1 for He, etc */ 00014 long int nelem) 00015 { 00016 long int ipHi, 00017 ipop, 00018 limit, 00019 low, 00020 n, 00021 ion, 00022 nshell; 00023 char chStat; 00024 double abundance; 00025 00026 DEBUG_ENTRY( "OpacityAdd1Element()" ); 00027 00028 /* this routine drives OpacityAdd1Subshell to put in total opacities for all shells*/ 00029 00030 /*begin sanity check */ 00031 ASSERT( (nelem >=0 ) && (nelem < LIMELM) ); 00032 00033 /* first do simple two-level systems - 00034 * this is number of species that are not treated on common iso-electronic series */ 00035 limit = nelem + 1 - NISO; 00036 /* this can be called with hydrogen itself, in which case nelem is 0, and limit is 00037 * -1 - do not do any of the simple ions */ 00038 limit = MAX2( 0 , limit ); 00039 00040 /* do not include the ion stages that have complete atoms, 00041 * currently H and He like iso sequences */ 00042 for( ion=0; ion < limit; ion++ ) 00043 { 00044 if( dense.xIonDense[nelem][ion] > 0. ) 00045 { 00046 /*start with static opacities, then do volatile*/ 00047 00048 chStat = 's'; 00049 /* number of bound electrons */ 00050 for( nshell=0; nshell < Heavy.nsShells[nelem][ion]; nshell++ ) 00051 { 00052 /* highest shell will be volatile*/ 00053 if( nshell== Heavy.nsShells[nelem][ion]-1 ) 00054 chStat = 'v'; 00055 /* set lower and upper limits to this range */ 00056 low = opac.ipElement[nelem][ion][nshell][0]; 00057 ipHi = opac.ipElement[nelem][ion][nshell][1]; 00058 ipop = opac.ipElement[nelem][ion][nshell][2]; 00059 /* OpacityAdd1Subshell will not do anything if static opacities do not need to be reset*/ 00060 OpacityAdd1Subshell(ipop,low,ipHi,dense.xIonDense[nelem][ion] , chStat ); 00061 } 00062 } 00063 } 00064 00065 /* now loop over all species done as large multi-level systems */ 00066 /* >>chng 02 jan 17, add loop over H and He like */ 00067 /* ion is on the c scale, =0 for HI, =1 for HeII */ 00068 for( ion=limit; ion<nelem+1; ++ion ) 00069 { 00070 /* ipISO is 0 for H-like, 1 for He-like */ 00071 long int ipISO = nelem-ion; 00072 00073 /* do multi level systems, but only if present 00074 * test for nelem+1 in case atom present but not ion, test is whether the 00075 * abundance of the recombined species is present */ 00076 /* >>chng 02 jan 17, sec dim had been nelem+1, change to ion+1 */ 00077 /*if( dense.xIonDense[nelem][nelem] > 0. )*/ 00078 if( dense.xIonDense[nelem][ion] > 0. ) 00079 { 00080 ASSERT( ipISO < NISO ); 00081 00082 /* do ground first, then all excited states */ 00083 n = 0; 00084 /* abundance of recombined species, which can be zero if no ion present */ 00085 abundance = StatesElemNEW[nelem][nelem-ipISO][n].Pop; 00086 00087 /* >>chng 02 may 06, add second test, had been just the chck on helium, 00088 * with no option to use new soln */ 00089 if( abundance == 0. ) 00090 { 00091 /* no ionized species, assume everything in ground */ 00092 abundance = dense.xIonDense[nelem][ion]; 00093 } 00094 00095 /* >>chng 02 jan 17, to arbitrary iso sequence */ 00096 /* use computed opacities and departure coef for level */ 00097 OpacityAdd1SubshellInduc( 00098 iso.ipOpac[ipISO][nelem][n], 00099 iso.ipIsoLevNIonCon[ipISO][nelem][n], 00100 /* the upper limit to the integration, 00101 * ground opacity goes up to the high-energy limit of code*/ 00102 rfield.nflux, 00103 /* the abundance of the ion */ 00104 abundance, 00105 /* departure coef, volatile opac, always reevaluate */ 00106 iso.DepartCoef[ipISO][nelem][n] , 'v' ); 00107 00108 /* do excited levvels, 00109 * this loop only if upper levels have finite population*/ 00110 if( StatesElemNEW[nelem][nelem-ipISO][3].Pop > 0. ) 00111 { 00112 char chType = 'v'; 00113 /* always want to evaluate all opacities for n=3, 4, use static opacities for higher levels */ 00114 /* >>chng 06 aug 17, should go to numLevels_local instead of _max */ 00115 for( long level =1; level < iso.numLevels_local[ipISO][nelem]; level++ ) 00116 { 00117 /* above 4 is static */ 00118 if( StatesElemNEW[nelem][nelem-ipISO][level].n >= 5 ) 00119 chType = 's'; 00120 00121 /* include correction for stimulated emission */ 00122 OpacityAdd1SubshellInduc( 00123 iso.ipOpac[ipISO][nelem][level], 00124 iso.ipIsoLevNIonCon[ipISO][nelem][level], 00125 /* the high energy bound of excited states is the 00126 * edge of the Lyman continuum */ 00127 iso.ipIsoLevNIonCon[ipISO][nelem][0], 00128 StatesElemNEW[nelem][nelem-ipISO][level].Pop, 00129 /* departure coef, volitile opacities */ 00130 iso.DepartCoef[ipISO][nelem][level] , chType ); 00131 } 00132 } 00133 } 00134 } 00135 return; 00136 }