00001
00002
00003
00004 #include "cddefines.h"
00005 #include "atmdat.h"
00006
00007 double atmdat_HS_caseB(
00008
00009 long int iHi,
00010 long int iLo,
00011
00012 long int nelem,
00013
00014 double TempIn,
00015
00016 double DenIn,
00017
00018 char chCase )
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034 {
00035 long int
00036 ipTemp,
00037 ipDens,
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
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
00054 --nelem;
00055
00056
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
00074
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
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
00106 if( DenIn > atmdat.Density[iCase][nelem][atmdat.nDensity[iCase][nelem]-1] )
00107 {
00108
00109 return -1;
00110 }
00111
00112 if( DenIn <= atmdat.Density[iCase][nelem][0] )
00113 {
00114
00115
00116 ipDens = 0;
00117 }
00118 else
00119 {
00120
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
00131 if( TempIn < atmdat.ElecTemp[iCase][nelem][0] ||
00132 TempIn > atmdat.ElecTemp[iCase][nelem][atmdat.ntemp[iCase][nelem]-1] )
00133 {
00134
00135 return -1;
00136 }
00137
00138
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
00147
00148 if( ipDens+1 > atmdat.nDensity[iCase][nelem]-1 )
00149 {
00150
00151 ipDensHi = atmdat.nDensity[iCase][nelem]-1;
00152 }
00153 else if( DenIn < atmdat.Density[iCase][nelem][0])
00154 {
00155
00156 ipDensHi = 0;
00157 }
00158 else
00159 {
00160 ipDensHi = ipDens+1;
00161 }
00162
00163
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
00180 ip = (int)((((atmdat.ncut[iCase][nelem]-ipUp)*(atmdat.ncut[iCase][nelem]+ipUp-1))/2)+ipLo - 1);
00181
00182
00183 ASSERT( ip < NLINEHS );
00184 ASSERT( ip >= 0 );
00185
00186
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 }