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 /*OpacityAdd1Subshell add opacity due to single shell to main opacity array*/ 00004 /*OpacityAdd1SubshellInduc add opacity of individual species, including stimulated emission */ 00005 #include "cddefines.h" 00006 #include "rfield.h" 00007 #include "hydrogenic.h" 00008 #include "opacity.h" 00009 00010 void OpacityAdd1Subshell( 00011 /*ipOpac is opacity index within opac opacity offset for this species */ 00012 long int ipOpac, 00013 /* lower freq limit to opacity range on energy mesh */ 00014 long int ipLowLim, 00015 /* upper limit to opacity range on energy mesh */ 00016 long int ipUpLim, 00017 /* abundance, we bail if zero */ 00018 realnum abundance, 00019 /* either static 's' or volitile 'v' */ 00020 char chStat ) 00021 { 00022 long int i , 00023 ipOffset, 00024 limit; 00025 00026 DEBUG_ENTRY( "OpacityAdd1Subshell()" ); 00027 00028 /* code spends roughly 20% of its time in this loop*/ 00029 00030 ASSERT( chStat == 's' || chStat == 'v' ); 00031 00032 ipOffset = ipOpac - ipLowLim; 00033 ASSERT( ipLowLim > 0 ); 00034 /* >>chng 02 aug 13, negative offset is ok, it is only this plus ipLowLim that matters */ 00035 /*ASSERT( ipOffset >= 0 );*/ 00036 00037 limit = MIN2(ipUpLim,rfield.nflux); 00038 00039 /* do nothing if abundance is zero, or if static opacities do not 00040 * need to be redone */ 00041 if( abundance <= 0. || (chStat=='s' && !opac.lgRedoStatic) ) 00042 { 00043 return; 00044 } 00045 00046 /* volative (outer shell, constantly reevaluated) or static opacity? */ 00047 if( chStat=='v' ) 00048 { 00049 for( i=ipLowLim-1; i < limit; i++ ) 00050 { 00051 opac.opacity_abs[i] += opac.OpacStack[i+ipOffset]*abundance; 00052 } 00053 } 00054 else 00055 { 00056 for( i=ipLowLim-1; i < limit; i++ ) 00057 { 00058 opac.OpacStatic[i] += opac.OpacStack[i+ipOffset]*abundance; 00059 } 00060 } 00061 return; 00062 } 00063 00064 /*OpacityAdd1SubshellInduc add opacity of individual species, including stimulated emission */ 00065 void OpacityAdd1SubshellInduc( 00066 /* pointer to opacity within opacity stack */ 00067 long int ipOpac, 00068 /* pointer to low energy in continuum array for this opacity band */ 00069 long int ipLowEnergy, 00070 /* pointer to high energy in continuum array for this opacity */ 00071 long int ipHiEnergy, 00072 /* this abundance of this species, may be zero */ 00073 double abundance, 00074 /* the departure coef, may be infinite or zero */ 00075 double DepartCoef , 00076 /* either 'v' for volitile or 's' for static opacities */ 00077 char chStat ) 00078 { 00079 long int i, 00080 iup, 00081 k; 00082 00083 DEBUG_ENTRY( "OpacityAdd1SubshellInduc()" ); 00084 00085 /* add opacity of individual species, including stimulated emission 00086 * abundance is the density of the lower level (cm^-3) 00087 * DepartCoef is its departure coefficient, can be zero */ 00088 00089 /* this is opacity offset, must be positive */ 00090 ASSERT( ipOpac > 0 ); 00091 00092 /* check that chStat is either 'v' or 's' */ 00093 ASSERT( chStat == 'v' || chStat == 's' ); 00094 00095 /* do nothing if abundance is zero, or if static opacities do not 00096 * need to be redone */ 00097 if( abundance <= 0. || (chStat=='s' && !opac.lgRedoStatic) ) 00098 { 00099 return; 00100 } 00101 00102 k = ipOpac - ipLowEnergy; 00103 00104 /* DepartCoef is dep coef, rfield.lgInducProcess is turned off with 'no indcued' command */ 00105 if( (DepartCoef > 1e-35 && rfield.lgInducProcess) && hydro.lgHInducImp ) 00106 { 00107 iup = MIN2(ipHiEnergy,rfield.nflux); 00108 /* >>>chng 99 apr 29, following was present, caused pdr to make opac at impossible energy*/ 00109 /*iup = MAX2(ipLowEnergy,iup);*/ 00110 double DepartCoefInv = 1./DepartCoef; 00111 if( chStat == 'v' ) 00112 { 00113 /* volitile opacities, always reevaluate */ 00114 for( i=ipLowEnergy-1; i < iup; i++ ) 00115 { 00116 opac.opacity_abs[i] += opac.OpacStack[i+k]*abundance* 00117 MAX2(0. , 1.- rfield.ContBoltz[i]*DepartCoefInv); 00118 } 00119 } 00120 else 00121 { 00122 /* static opacities, save in special array */ 00123 for( i=ipLowEnergy-1; i < iup; i++ ) 00124 { 00125 opac.OpacStatic[i] += opac.OpacStack[i+k]*abundance* 00126 MAX2(0. , 1.- rfield.ContBoltz[i]*DepartCoefInv); 00127 } 00128 } 00129 } 00130 00131 else 00132 { 00133 /* DepartCoef is the departure coef, which can go to zero, 00134 * neglect stimulated emission in this case */ 00135 iup = MIN2(ipHiEnergy,rfield.nflux); 00136 /* >>>chng 99 apr 29, following was present, caused pdr to make opac at impossible energy*/ 00137 /*iup = MAX2(ipLowEnergy,iup);*/ 00138 if( chStat == 'v' ) 00139 { 00140 for( i=ipLowEnergy-1; i < iup; i++ ) 00141 { 00142 opac.opacity_abs[i] += opac.OpacStack[i+k]*abundance; 00143 } 00144 } 00145 else 00146 { 00147 for( i=ipLowEnergy-1; i < iup; i++ ) 00148 { 00149 opac.OpacStatic[i] += opac.OpacStack[i+k]*abundance; 00150 } 00151 } 00152 } 00153 00154 return; 00155 }