00001 /* This file is part of Cloudy and is copyright (C)1978-2007 by Gary J. Ferland 00002 * 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 /* do ground first, then all excited states */ 00080 n = 0; 00081 /* abundance of recombined species, which can be zero if no ion present */ 00082 /* >>chng 02 jan 17, to arbitrary iso sequence */ 00083 /*abundance = iso.Pop2Ion[ipH_LIKE][nelem][n]*dense.xIonDense[nelem][nelem+1];*/ 00084 abundance = iso.Pop2Ion[ipISO][nelem][n]*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 /* >>chng 02 jan 17, logical error here, must use dense.xIonDense if abundance is 0 */ 00090 { 00091 /* no ionized species, assume everything in ground */ 00092 /* >>chng 02 jan 17, to arbitrary iso sequence */ 00093 abundance = dense.xIonDense[nelem][ion]; 00094 } 00095 00096 /* >>chng 02 jan 17, to arbitrary iso sequence */ 00097 /* use computed opacities and departure coef for level */ 00098 OpacityAdd1SubshellInduc( 00099 iso.ipOpac[ipISO][nelem][n], 00100 iso.ipIsoLevNIonCon[ipISO][nelem][n], 00101 /* the upper limit to the integration, 00102 * ground opacity goes up to the high-energy limit of code*/ 00103 rfield.nflux, 00104 /* the abundance of the ion */ 00105 abundance, 00106 /* departure coef, volatile opac, always reevaluate */ 00107 iso.DepartCoef[ipISO][nelem][n] , 'v' ); 00108 00109 /* do excited levvels, 00110 * this loop only if upper levels have finite population*/ 00111 /* >>chng 02 jan 17, to arbitrary iso sequence */ 00112 /*if( iso.Pop2Ion[ipH_LIKE][nelem][3]*dense.xIonDense[nelem][nelem+1] > 0. )*/ 00113 if( iso.Pop2Ion[ipISO][nelem][3]*dense.xIonDense[nelem][ion+1] > 0. ) 00114 { 00115 char chType = 'v'; 00116 /* always want to evaluate all opacities for n=3, 4, use static opacities for higher levels */ 00117 /* >>chng 02 jan 17, to arbitrary iso sequence, use single loop, setting type of opacity */ 00118 /* >>chng 06 aug 17, should go to numLevels_local instead of _max */ 00119 for( n=1; n < iso.numLevels_local[ipISO][nelem]; n++ ) 00120 { 00121 /* include correction for stimulated emission */ 00122 OpacityAdd1SubshellInduc( 00123 iso.ipOpac[ipISO][nelem][n], 00124 iso.ipIsoLevNIonCon[ipISO][nelem][n], 00125 /* the high energy bound of excited states is the 00126 * edge of the Lyman continuum */ 00127 iso.ipIsoLevNIonCon[ipISO][nelem][0], 00128 iso.Pop2Ion[ipISO][nelem][n]*dense.xIonDense[nelem][ion+1], 00129 /* departure coef, volitile opacities */ 00130 iso.DepartCoef[ipISO][nelem][n] , chType ); 00131 /* above 4 is static */ 00132 if( n==4) 00133 chType = 's'; 00134 } 00135 } 00136 } 00137 } 00138 00139 DEBUG_EXIT( "OpacityAdd1Element()" ); 00140 return; 00141 } 00142