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