cloudy trunk
Loading...
Searching...
No Matches
atmdat_HS_caseb.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/* atmdat_HS_caseB - interpolate on line emissivities from Storey & Hummer tables for hydrogen */
4#include "cddefines.h"
5#include "atmdat.h"
6
8 /* upper and lower quantum numbers*/
9 long int iHi,
10 long int iLo,
11 /* element number, 1 or 2 at this point, but decremented to C scale later */
12 long int nelem,
13 /* temperature to interpolate, for H between 500-30,000K*/
14 double TempIn,
15 /* density to interpolate */
16 double DenIn,
17 /* case - 'a' or 'b' */
18 char chCase )
19
20/* general utility to interpolate line emissivities from the Storey & Hummer tables
21 of case B emissivities.
22 iHi<25, iLo, the principal quantum
23 numbers, and are upper and lower levels in any order.
24 nelem element number on physicial scale, 1 or 2 have data
25 TempIn = temperature, and must lie within the range of the table, which depends on
26 the ion charge, and is 500 - 30,000K for hydrogen.
27 DenIn is the density and must not exceed the high density limit to the table.
28
29 routine returns -1 if conditions outside temperature range, or
30 above density range of tabulated case b results
31 If desired density is low limit, lower limit is used instead
32*/
33
34{
35 long int
36 ipTemp, /*pointer to temperature*/
37 ipDens, /*pointer to density*/
38 ipDensHi ,
39 ipTempHi;
40 int ipUp , ipLo , ip;
41 double x1 , x2 , yy1 , yy2 , z1 , z2 , za , zb ,z;
42 int iCase;
43
44 DEBUG_ENTRY( "atmdat_HS_caseB()" );
45
46 /*make sure nelem is 1 or 2*/
47 if( nelem<ipHYDROGEN || nelem> HS_NZ )
48 {
49 printf("atmdat_HS_caseB called with improper nelem, was%li and must be 1 or 2",nelem);
51
52 }
53 /* now back nelem back one, to be on c scale */
54 --nelem;
55
56 /* case A or B? */
57 if( chCase == 'a' || chCase=='A' )
58 {
59 iCase = 0;
60 }
61 else if( chCase == 'b' || chCase=='B' )
62 {
63 iCase = 1;
64 }
65 else
66 {
67 iCase = false;
68 printf("atmdat_HS_caseB called with improper case, was %c and must be A or B",chCase);
70 }
71
72 /*===========================================================*/
73 /* following is option to have principal quantum number given in either order,
74 * final result is that ipUp and ipLo will be the upper and lower levels */
75 if( iHi > iLo )
76 {
77 ipUp = (int)iHi; ipLo = (int)iLo;
78 }
79 else if( iHi < iLo )
80 {
81 ipUp = (int)iLo; ipLo = (int)iHi;
82 }
83 else
84 {
85 printf("atmdat_HS_caseB called with indices equal, %ld %ld \n",iHi,iLo);
87 }
88
89 /* now check that they are in range of the predicted levels of their model atom*/
90 if( ipLo <1 )
91 {
92 printf(" atmdat_HS_caseB called with lower quantum number less than 1, = %i\n",
93 ipLo);
95 }
96
97 if( ipUp >25 )
98 {
99 printf(" atmdat_HS_caseB called with upper quantum number greater than 25, = %i\n",
100 ipUp);
102 }
103
104 /*===========================================================*/
105 /*bail if above high density limit */
106 if( DenIn > atmdat.Density[iCase][nelem][atmdat.nDensity[iCase][nelem]-1] )
107 {
108 /* this is flag saying bogus results */
109 return -1;
110 }
111
112 if( DenIn <= atmdat.Density[iCase][nelem][0] )
113 {
114 /* this case, desired density is below lower limit to table,
115 * just use the lower limit */
116 ipDens = 0;
117 }
118 else
119 {
120 /* this case find where within table density lies */
121 for( ipDens=0; ipDens < atmdat.nDensity[iCase][nelem]-1; ipDens++ )
122 {
123 if( DenIn >= atmdat.Density[iCase][nelem][ipDens] &&
124 DenIn < atmdat.Density[iCase][nelem][ipDens+1] ) break;
125 }
126 }
127
128
129 /*===========================================================*/
130 /* confirm within temperature range*/
131 if( TempIn < atmdat.ElecTemp[iCase][nelem][0] ||
132 TempIn > atmdat.ElecTemp[iCase][nelem][atmdat.ntemp[iCase][nelem]-1] )
133 {
134 /* this is flag saying bogus results */
135 return -1;
136 }
137
138 /* find where within grid this temperature lies */
139 for( ipTemp=0; ipTemp < atmdat.ntemp[iCase][nelem]-1; ipTemp++ )
140 {
141 if( TempIn >= atmdat.ElecTemp[iCase][nelem][ipTemp] &&
142 TempIn < atmdat.ElecTemp[iCase][nelem][ipTemp+1] ) break;
143 }
144
145 /*===========================================================*/
146 /*we now have the array indices within the temperature array*/
147
148 if( ipDens+1 > atmdat.nDensity[iCase][nelem]-1 )
149 {
150 /* special case, when cell is highest density point */
151 ipDensHi = atmdat.nDensity[iCase][nelem]-1;
152 }
153 else if( DenIn < atmdat.Density[iCase][nelem][0])
154 {
155 /* or density below lower limit to table, set both bounds to 0 */
156 ipDensHi = 0;
157 }
158 else
159 {
160 ipDensHi = ipDens+1;
161 }
162
163 /*special case, if cell is highest temperature point*/
164 if( ipTemp+1 > atmdat.ntemp[iCase][nelem]-1 )
165 {
166 ipTempHi = atmdat.ntemp[iCase][nelem]-1;
167 }
168 else
169 {
170 ipTempHi = ipTemp+1;
171 }
172
173 x1 = log10( atmdat.Density[iCase][nelem][ipDens] );
174 x2 = log10( atmdat.Density[iCase][nelem][ipDensHi] );
175
176 yy1 = log10( atmdat.ElecTemp[iCase][nelem][ipTemp] );
177 yy2 = log10( atmdat.ElecTemp[iCase][nelem][ipTempHi] );
178
179 /*now generate the index to the array, expression from Storey code -1 for c*/
180 ip = (int)((((atmdat.ncut[iCase][nelem]-ipUp)*(atmdat.ncut[iCase][nelem]+ipUp-1))/2)+ipLo - 1);
181
182 /*pointer must lie within line array*/
183 ASSERT( ip < NLINEHS );
184 ASSERT( ip >= 0 );
185
186 /* interpolate on emission rate*/
187 z1 = log10( atmdat.Emiss[iCase][nelem][ipTemp][ipDens][ip]);
188 z2 = log10( atmdat.Emiss[iCase][nelem][ipTemp][ipDensHi][ip]);
189
190 if( fp_equal( x2, x1 ) )
191 {
192 za = z2;
193 }
194 else
195 {
196 za = z1 + log10( DenIn / atmdat.Density[iCase][nelem][ipDens] ) * (z2-z1)/(x2-x1);
197 }
198
199 z1 = log10( atmdat.Emiss[iCase][nelem][ipTempHi][ipDens][ip]);
200 z2 = log10( atmdat.Emiss[iCase][nelem][ipTempHi][ipDensHi][ip]);
201
202 if( fp_equal( x2, x1 ) )
203 {
204 zb = z2;
205 }
206 else
207 {
208 zb = z1 + log10( DenIn / atmdat.Density[iCase][nelem][ipDens] ) * (z2-z1)/(x2-x1);
209 }
210
211 if( fp_equal( yy2, yy1 ) )
212 {
213 z = zb;
214 }
215 else
216 {
217 z = za + log10( TempIn / atmdat.ElecTemp[iCase][nelem][ipTemp] ) * (zb-za)/(yy2-yy1);
218 }
219
220 return ( pow(10.,z) );
221}
t_atmdat atmdat
Definition atmdat.cpp:6
#define NLINEHS
Definition atmdat.h:124
#define HS_NZ
Definition atmdat.h:125
static double x2[63]
static double x1[83]
double atmdat_HS_caseB(long int iHi, long int iLo, long int nelem, double TempIn, double DenIn, char chCase)
#define ASSERT(exp)
Definition cddefines.h:578
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684