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 /* atmdat_HS_caseB - interpolate on line emissivities from Storey & Hummer tables for hydrogen */ 00004 #include "cddefines.h" 00005 #include "atmdat.h" 00006 00007 double atmdat_HS_caseB( 00008 /* upper and lower quantum numbers*/ 00009 long int iHi, 00010 long int iLo, 00011 /* element number, 1 or 2 at this point, but decremented to C scale later */ 00012 long int nelem, 00013 /* temperature to interpolate, for H between 500-30,000K*/ 00014 double TempIn, 00015 /* density to interpolate */ 00016 double DenIn, 00017 /* case - 'a' or 'b' */ 00018 char chCase ) 00019 00020 /* general utility to interpolate line emissivities from the Storey & Hummer tables 00021 of case B emissivities. 00022 iHi<25, iLo, the principal quantum 00023 numbers, and are upper and lower levels in any order. 00024 nelem element number on physicial scale, 1 or 2 have data 00025 TempIn = temperature, and must lie within the range of the table, which depends on 00026 the ion charge, and is 500 - 30,000K for hydrogen. 00027 DenIn is the density and must not exceed the high density limit to the table. 00028 00029 routine returns -1 if conditions outside temperature range, or 00030 above density range of tabulated case b results 00031 If desired density is low limit, lower limit is used instead 00032 */ 00033 00034 { 00035 long int 00036 ipTemp, /*pointer to temperature*/ 00037 ipDens, /*pointer to density*/ 00038 ipDensHi , 00039 ipTempHi; 00040 int ipUp , ipLo , ip; 00041 double x1 , x2 , yy1 , yy2 , z1 , z2 , za , zb ,z; 00042 int iCase; 00043 00044 DEBUG_ENTRY( "atmdat_HS_caseB()" ); 00045 00046 /*make sure nelem is 1 or 2*/ 00047 if( nelem<ipHYDROGEN || nelem> HS_NZ ) 00048 { 00049 printf("atmdat_HS_caseB called with improper nelem, was%li and must be 1 or 2",nelem); 00050 cdEXIT(EXIT_FAILURE); 00051 00052 } 00053 /* now back nelem back one, to be on c scale */ 00054 --nelem; 00055 00056 /* case A or B? */ 00057 if( chCase == 'a' || chCase=='A' ) 00058 { 00059 iCase = 0; 00060 } 00061 else if( chCase == 'b' || chCase=='B' ) 00062 { 00063 iCase = 1; 00064 } 00065 else 00066 { 00067 iCase = false; 00068 printf("atmdat_HS_caseB called with improper case, was %c and must be A or B",chCase); 00069 cdEXIT(EXIT_FAILURE); 00070 } 00071 00072 /*===========================================================*/ 00073 /* following is option to have principal quantum number given in either order, 00074 * final result is that ipUp and ipLo will be the upper and lower levels */ 00075 if( iHi > iLo ) 00076 { 00077 ipUp = (int)iHi; ipLo = (int)iLo; 00078 } 00079 else if( iHi < iLo ) 00080 { 00081 ipUp = (int)iLo; ipLo = (int)iHi; 00082 } 00083 else 00084 { 00085 printf("atmdat_HS_caseB called with indices equal, %ld %ld \n",iHi,iLo); 00086 cdEXIT(EXIT_FAILURE); 00087 } 00088 00089 /* now check that they are in range of the predicted levels of their model atom*/ 00090 if( ipLo <1 ) 00091 { 00092 printf(" atmdat_HS_caseB called with lower quantum number less than 1, = %i\n", 00093 ipLo); 00094 cdEXIT(EXIT_FAILURE); 00095 } 00096 00097 if( ipUp >25 ) 00098 { 00099 printf(" atmdat_HS_caseB called with upper quantum number greater than 25, = %i\n", 00100 ipUp); 00101 cdEXIT(EXIT_FAILURE); 00102 } 00103 00104 /*===========================================================*/ 00105 /*bail if above high density limit */ 00106 if( DenIn > atmdat.Density[iCase][nelem][atmdat.nDensity[iCase][nelem]-1] ) 00107 { 00108 /* this is flag saying bogus results */ 00109 return -1; 00110 } 00111 00112 if( DenIn <= atmdat.Density[iCase][nelem][0] ) 00113 { 00114 /* this case, desired density is below lower limit to table, 00115 * just use the lower limit */ 00116 ipDens = 0; 00117 } 00118 else 00119 { 00120 /* this case find where within table density lies */ 00121 for( ipDens=0; ipDens < atmdat.nDensity[iCase][nelem]-1; ipDens++ ) 00122 { 00123 if( DenIn >= atmdat.Density[iCase][nelem][ipDens] && 00124 DenIn < atmdat.Density[iCase][nelem][ipDens+1] ) break; 00125 } 00126 } 00127 00128 00129 /*===========================================================*/ 00130 /* confirm within temperature range*/ 00131 if( TempIn < atmdat.ElecTemp[iCase][nelem][0] || 00132 TempIn > atmdat.ElecTemp[iCase][nelem][atmdat.ntemp[iCase][nelem]-1] ) 00133 { 00134 /* this is flag saying bogus results */ 00135 return -1; 00136 } 00137 00138 /* find where within grid this temperature lies */ 00139 for( ipTemp=0; ipTemp < atmdat.ntemp[iCase][nelem]-1; ipTemp++ ) 00140 { 00141 if( TempIn >= atmdat.ElecTemp[iCase][nelem][ipTemp] && 00142 TempIn < atmdat.ElecTemp[iCase][nelem][ipTemp+1] ) break; 00143 } 00144 00145 /*===========================================================*/ 00146 /*we now have the array indices within the temperature array*/ 00147 00148 if( ipDens+1 > atmdat.nDensity[iCase][nelem]-1 ) 00149 { 00150 /* special case, when cell is highest density point */ 00151 ipDensHi = atmdat.nDensity[iCase][nelem]-1; 00152 } 00153 else if( DenIn < atmdat.Density[iCase][nelem][0]) 00154 { 00155 /* or density below lower limit to table, set both bounds to 0 */ 00156 ipDensHi = 0; 00157 } 00158 else 00159 { 00160 ipDensHi = ipDens+1; 00161 } 00162 00163 /*special case, if cell is highest temperature point*/ 00164 if( ipTemp+1 > atmdat.ntemp[iCase][nelem]-1 ) 00165 { 00166 ipTempHi = atmdat.ntemp[iCase][nelem]-1; 00167 } 00168 else 00169 { 00170 ipTempHi = ipTemp+1; 00171 } 00172 00173 x1 = log10( atmdat.Density[iCase][nelem][ipDens] ); 00174 x2 = log10( atmdat.Density[iCase][nelem][ipDensHi] ); 00175 00176 yy1 = log10( atmdat.ElecTemp[iCase][nelem][ipTemp] ); 00177 yy2 = log10( atmdat.ElecTemp[iCase][nelem][ipTempHi] ); 00178 00179 /*now generate the index to the array, expression from Storey code -1 for c*/ 00180 ip = (int)((((atmdat.ncut[iCase][nelem]-ipUp)*(atmdat.ncut[iCase][nelem]+ipUp-1))/2)+ipLo - 1); 00181 00182 /*pointer must lie within line array*/ 00183 ASSERT( ip < NLINEHS ); 00184 ASSERT( ip >= 0 ); 00185 00186 /* interpolate on emission rate*/ 00187 z1 = log10( atmdat.Emiss[iCase][nelem][ipTemp][ipDens][ip]); 00188 z2 = log10( atmdat.Emiss[iCase][nelem][ipTemp][ipDensHi][ip]); 00189 00190 if( fp_equal( x2, x1 ) ) 00191 { 00192 za = z2; 00193 } 00194 else 00195 { 00196 za = z1 + log10( DenIn / atmdat.Density[iCase][nelem][ipDens] ) * (z2-z1)/(x2-x1); 00197 } 00198 00199 z1 = log10( atmdat.Emiss[iCase][nelem][ipTempHi][ipDens][ip]); 00200 z2 = log10( atmdat.Emiss[iCase][nelem][ipTempHi][ipDensHi][ip]); 00201 00202 if( fp_equal( x2, x1 ) ) 00203 { 00204 zb = z2; 00205 } 00206 else 00207 { 00208 zb = z1 + log10( DenIn / atmdat.Density[iCase][nelem][ipDens] ) * (z2-z1)/(x2-x1); 00209 } 00210 00211 if( fp_equal( yy2, yy1 ) ) 00212 { 00213 z = zb; 00214 } 00215 else 00216 { 00217 z = za + log10( TempIn / atmdat.ElecTemp[iCase][nelem][ipTemp] ) * (zb-za)/(yy2-yy1); 00218 } 00219 00220 return ( pow(10.,z) ); 00221 }