cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
atmdat_HS_caseb.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 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);
50  cdEXIT(EXIT_FAILURE);
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);
69  cdEXIT(EXIT_FAILURE);
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);
86  cdEXIT(EXIT_FAILURE);
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);
94  cdEXIT(EXIT_FAILURE);
95  }
96 
97  if( ipUp >25 )
98  {
99  printf(" atmdat_HS_caseB called with upper quantum number greater than 25, = %i\n",
100  ipUp);
101  cdEXIT(EXIT_FAILURE);
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 }

Generated for cloudy by doxygen 1.8.3.1