41 double x1 ,
x2 , yy1 , yy2 , z1 , z2 , za , zb ,z;
47 if( nelem<ipHYDROGEN || nelem>
HS_NZ )
49 printf(
"atmdat_HS_caseB called with improper nelem, was%li and must be 1 or 2",nelem);
57 if( chCase ==
'a' || chCase==
'A' )
61 else if( chCase ==
'b' || chCase==
'B' )
68 printf(
"atmdat_HS_caseB called with improper case, was %c and must be A or B",chCase);
77 ipUp = (int)iHi; ipLo = (int)iLo;
81 ipUp = (int)iLo; ipLo = (int)iHi;
85 printf(
"atmdat_HS_caseB called with indices equal, %ld %ld \n",iHi,iLo);
92 printf(
" atmdat_HS_caseB called with lower quantum number less than 1, = %i\n",
99 printf(
" atmdat_HS_caseB called with upper quantum number greater than 25, = %i\n",
106 if( DenIn >
atmdat.Density[iCase][nelem][
atmdat.nDensity[iCase][nelem]-1] )
112 if( DenIn <=
atmdat.Density[iCase][nelem][0] )
121 for( ipDens=0; ipDens <
atmdat.nDensity[iCase][nelem]-1; ipDens++ )
123 if( DenIn >=
atmdat.Density[iCase][nelem][ipDens] &&
124 DenIn <
atmdat.Density[iCase][nelem][ipDens+1] )
break;
131 if( TempIn <
atmdat.ElecTemp[iCase][nelem][0] ||
132 TempIn >
atmdat.ElecTemp[iCase][nelem][
atmdat.ntemp[iCase][nelem]-1] )
139 for( ipTemp=0; ipTemp <
atmdat.ntemp[iCase][nelem]-1; ipTemp++ )
141 if( TempIn >=
atmdat.ElecTemp[iCase][nelem][ipTemp] &&
142 TempIn <
atmdat.ElecTemp[iCase][nelem][ipTemp+1] )
break;
148 if( ipDens+1 >
atmdat.nDensity[iCase][nelem]-1 )
151 ipDensHi =
atmdat.nDensity[iCase][nelem]-1;
153 else if( DenIn <
atmdat.Density[iCase][nelem][0])
164 if( ipTemp+1 >
atmdat.ntemp[iCase][nelem]-1 )
166 ipTempHi =
atmdat.ntemp[iCase][nelem]-1;
173 x1 = log10(
atmdat.Density[iCase][nelem][ipDens] );
174 x2 = log10(
atmdat.Density[iCase][nelem][ipDensHi] );
176 yy1 = log10(
atmdat.ElecTemp[iCase][nelem][ipTemp] );
177 yy2 = log10(
atmdat.ElecTemp[iCase][nelem][ipTempHi] );
180 ip = (int)((((
atmdat.ncut[iCase][nelem]-ipUp)*(
atmdat.ncut[iCase][nelem]+ipUp-1))/2)+ipLo - 1);
187 z1 = log10(
atmdat.Emiss[iCase][nelem][ipTemp][ipDens][ip]);
188 z2 = log10(
atmdat.Emiss[iCase][nelem][ipTemp][ipDensHi][ip]);
196 za = z1 + log10( DenIn /
atmdat.Density[iCase][nelem][ipDens] ) * (z2-z1)/(
x2-
x1);
199 z1 = log10(
atmdat.Emiss[iCase][nelem][ipTempHi][ipDens][ip]);
200 z2 = log10(
atmdat.Emiss[iCase][nelem][ipTempHi][ipDensHi][ip]);
208 zb = z1 + log10( DenIn /
atmdat.Density[iCase][nelem][ipDens] ) * (z2-z1)/(
x2-
x1);
217 z = za + log10( TempIn /
atmdat.ElecTemp[iCase][nelem][ipTemp] ) * (zb-za)/(yy2-yy1);
220 return ( pow(10.,z) );
double atmdat_HS_caseB(long int iHi, long int iLo, long int nelem, double TempIn, double DenIn, char chCase)