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
00045 if( nelem<ipHYDROGEN || nelem> HS_NZ )
00046 {
00047 printf("atmdat_HS_caseB called with improper nelem, was%li and must be 1 or 2",nelem);
00048 cdEXIT(EXIT_FAILURE);
00049
00050 }
00051
00052 --nelem;
00053
00054
00055 if( chCase == 'a' || chCase=='A' )
00056 {
00057 iCase = 0;
00058 }
00059 else if( chCase == 'b' || chCase=='B' )
00060 {
00061 iCase = 1;
00062 }
00063 else
00064 {
00065 iCase = false;
00066 printf("atmdat_HS_caseB called with improper case, was %c and must be A or B",chCase);
00067 cdEXIT(EXIT_FAILURE);
00068 }
00069
00070
00071
00072
00073 if( iHi > iLo )
00074 {
00075 ipUp = (int)iHi; ipLo = (int)iLo;
00076 }
00077 else if( iHi < iLo )
00078 {
00079 ipUp = (int)iLo; ipLo = (int)iHi;
00080 }
00081 else
00082 {
00083 printf("atmdat_HS_caseB called with indices equal, %ld %ld \n",iHi,iLo);
00084 cdEXIT(EXIT_FAILURE);
00085 }
00086
00087
00088 if( ipLo <1 )
00089 {
00090 printf(" atmdat_HS_caseB called with lower quantum number less than 1, = %i\n",
00091 ipLo);
00092 cdEXIT(EXIT_FAILURE);
00093 }
00094
00095 if( ipUp >25 )
00096 {
00097 printf(" atmdat_HS_caseB called with upper quantum number greater than 25, = %i\n",
00098 ipUp);
00099 cdEXIT(EXIT_FAILURE);
00100 }
00101
00102
00103
00104 if( DenIn > atmdat.Density[iCase][nelem][atmdat.nDensity[iCase][nelem]-1] )
00105 {
00106
00107 return -1;
00108 }
00109
00110 if( DenIn <= atmdat.Density[iCase][nelem][0] )
00111 {
00112
00113
00114 ipDens = 0;
00115 }
00116 else
00117 {
00118
00119 for( ipDens=0; ipDens < atmdat.nDensity[iCase][nelem]-1; ipDens++ )
00120 {
00121 if( DenIn >= atmdat.Density[iCase][nelem][ipDens] &&
00122 DenIn < atmdat.Density[iCase][nelem][ipDens+1] ) break;
00123 }
00124 }
00125
00126
00127
00128
00129 if( TempIn < atmdat.ElecTemp[iCase][nelem][0] ||
00130 TempIn > atmdat.ElecTemp[iCase][nelem][atmdat.ntemp[iCase][nelem]-1] )
00131 {
00132
00133 return -1;
00134 }
00135
00136
00137 for( ipTemp=0; ipTemp < atmdat.ntemp[iCase][nelem]-1; ipTemp++ )
00138 {
00139 if( TempIn >= atmdat.ElecTemp[iCase][nelem][ipTemp] &&
00140 TempIn < atmdat.ElecTemp[iCase][nelem][ipTemp+1] ) break;
00141 }
00142
00143
00144
00145
00146 if( ipDens+1 > atmdat.nDensity[iCase][nelem]-1 )
00147 {
00148
00149 ipDensHi = atmdat.nDensity[iCase][nelem]-1;
00150 }
00151 else if( DenIn < atmdat.Density[iCase][nelem][0])
00152 {
00153
00154 ipDensHi = 0;
00155 }
00156 else
00157 {
00158 ipDensHi = ipDens+1;
00159 }
00160
00161
00162 if( ipTemp+1 > atmdat.ntemp[iCase][nelem]-1 )
00163 {
00164 ipTempHi = atmdat.ntemp[iCase][nelem]-1;
00165 }
00166 else
00167 {
00168 ipTempHi = ipTemp+1;
00169 }
00170
00171 x1 = log10( atmdat.Density[iCase][nelem][ipDens] );
00172 x2 = log10( atmdat.Density[iCase][nelem][ipDensHi] );
00173
00174 yy1 = log10( atmdat.ElecTemp[iCase][nelem][ipTemp] );
00175 yy2 = log10( atmdat.ElecTemp[iCase][nelem][ipTempHi] );
00176
00177
00178 ip = (int)((((atmdat.ncut[iCase][nelem]-ipUp)*(atmdat.ncut[iCase][nelem]+ipUp-1))/2)+ipLo - 1);
00179
00180
00181 ASSERT( ip<NLINEHS );
00182 ASSERT( ip >= 0 );
00183
00184
00185 z1 = log10( atmdat.Emiss[iCase][nelem][ipTemp][ipDens][ip]);
00186 z2 = log10( atmdat.Emiss[iCase][nelem][ipTemp][ipDensHi][ip]);
00187
00188
00189 if(x2 == x1 )
00190 {
00191 za = z2;
00192 }
00193 else
00194 {
00195 za = z1 + log10( DenIn / atmdat.Density[iCase][nelem][ipDens] ) * (z2-z1)/(x2-x1);
00196 }
00197
00198
00199 z1 = log10( atmdat.Emiss[iCase][nelem][ipTempHi][ipDens][ip]);
00200 z2 = log10( atmdat.Emiss[iCase][nelem][ipTempHi][ipDensHi][ip]);
00201
00202
00203 if(x2 == x1 )
00204 {
00205 zb = z2;
00206 }
00207 else
00208 {
00209 zb = z1 + log10( DenIn / atmdat.Density[iCase][nelem][ipDens] ) * (z2-z1)/(x2-x1);
00210 }
00211
00212 if( yy2 == yy1 )
00213 {
00214 z = zb;
00215 }
00216 else
00217 {
00218 z = za + log10( TempIn / atmdat.ElecTemp[iCase][nelem][ipTemp] ) * (zb-za)/(yy2-yy1);
00219 }
00220
00221
00222 return ( pow(10.,z) );
00223 }
00224