cloudy  trunk
 All Data Structures Namespaces 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-2017 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, abort if iHi>25*/
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  printf("atmdat_HS_caseB called with improper case, was %c and must be A or B",chCase);
69  }
70 
71  /*===========================================================*/
72  /* following is option to have principal quantum number given in either order,
73  * final result is that ipUp and ipLo will be the upper and lower levels */
74  if( iHi > iLo )
75  {
76  ipUp = (int)iHi; ipLo = (int)iLo;
77  }
78  else if( iHi < iLo )
79  {
80  ipUp = (int)iLo; ipLo = (int)iHi;
81  }
82  else
83  {
84  printf("atmdat_HS_caseB called with indices equal, %ld %ld \n",iHi,iLo);
86  }
87 
88  /* now check that they are in range of the predicted levels of their model atom*/
89  if( ipLo <1 )
90  {
91  printf(" atmdat_HS_caseB called with lower quantum number less than 1, = %i\n",
92  ipLo);
94  }
95 
96  if( ipUp >25 )
97  {
98  printf(" atmdat_HS_caseB called with upper quantum number greater than 25, = %i\n",
99  ipUp);
101  }
102 
103  /*===========================================================*/
104  /*bail if above high density limit */
105  if( DenIn > atmdat.Density[iCase][nelem][atmdat.nDensity[iCase][nelem]-1] )
106  {
107  /* this is flag saying bogus results */
108  return -1;
109  }
110 
111  if( DenIn <= atmdat.Density[iCase][nelem][0] )
112  {
113  /* this case, desired density is below lower limit to table,
114  * just use the lower limit */
115  ipDens = 0;
116  }
117  else
118  {
119  /* this case find where within table density lies */
120  for( ipDens=0; ipDens < atmdat.nDensity[iCase][nelem]-1; ipDens++ )
121  {
122  if( DenIn >= atmdat.Density[iCase][nelem][ipDens] &&
123  DenIn < atmdat.Density[iCase][nelem][ipDens+1] ) break;
124  }
125  }
126 
127 
128  /*===========================================================*/
129  /* confirm within temperature range*/
130  if( TempIn < atmdat.ElecTemp[iCase][nelem][0] ||
131  TempIn > atmdat.ElecTemp[iCase][nelem][atmdat.ntemp[iCase][nelem]-1] )
132  {
133  /* this is flag saying bogus results */
134  return -1;
135  }
136 
137  /* find where within grid this temperature lies */
138  for( ipTemp=0; ipTemp < atmdat.ntemp[iCase][nelem]-1; ipTemp++ )
139  {
140  if( TempIn >= atmdat.ElecTemp[iCase][nelem][ipTemp] &&
141  TempIn < atmdat.ElecTemp[iCase][nelem][ipTemp+1] ) break;
142  }
143 
144  /*===========================================================*/
145  /*we now have the array indices within the temperature array*/
146 
147  if( ipDens+1 > atmdat.nDensity[iCase][nelem]-1 )
148  {
149  /* special case, when cell is highest density point */
150  ipDensHi = atmdat.nDensity[iCase][nelem]-1;
151  }
152  else if( DenIn < atmdat.Density[iCase][nelem][0])
153  {
154  /* or density below lower limit to table, set both bounds to 0 */
155  ipDensHi = 0;
156  }
157  else
158  {
159  ipDensHi = ipDens+1;
160  }
161 
162  /*special case, if cell is highest temperature point*/
163  if( ipTemp+1 > atmdat.ntemp[iCase][nelem]-1 )
164  {
165  ipTempHi = atmdat.ntemp[iCase][nelem]-1;
166  }
167  else
168  {
169  ipTempHi = ipTemp+1;
170  }
171 
172  x1 = log10( atmdat.Density[iCase][nelem][ipDens] );
173  x2 = log10( atmdat.Density[iCase][nelem][ipDensHi] );
174 
175  yy1 = log10( atmdat.ElecTemp[iCase][nelem][ipTemp] );
176  yy2 = log10( atmdat.ElecTemp[iCase][nelem][ipTempHi] );
177 
178  /*now generate the index to the array, expression from Storey code -1 for c*/
179  ip = (int)((((atmdat.ncut[iCase][nelem]-ipUp)*(atmdat.ncut[iCase][nelem]+ipUp-1))/2)+ipLo - 1);
180 
181  /*pointer must lie within line array*/
182  ASSERT( ip < NLINEHS );
183  ASSERT( ip >= 0 );
184 
185  /* interpolate on emission rate*/
186  z1 = log10( atmdat.Emiss[iCase][nelem][ipTemp][ipDens][ip]);
187  z2 = log10( atmdat.Emiss[iCase][nelem][ipTemp][ipDensHi][ip]);
188 
189  if( fp_equal( x2, x1 ) )
190  {
191  za = z2;
192  }
193  else
194  {
195  za = z1 + log10( DenIn / atmdat.Density[iCase][nelem][ipDens] ) * (z2-z1)/(x2-x1);
196  }
197 
198  z1 = log10( atmdat.Emiss[iCase][nelem][ipTempHi][ipDens][ip]);
199  z2 = log10( atmdat.Emiss[iCase][nelem][ipTempHi][ipDensHi][ip]);
200 
201  if( fp_equal( x2, x1 ) )
202  {
203  zb = z2;
204  }
205  else
206  {
207  zb = z1 + log10( DenIn / atmdat.Density[iCase][nelem][ipDens] ) * (z2-z1)/(x2-x1);
208  }
209 
210  if( fp_equal( yy2, yy1 ) )
211  {
212  z = zb;
213  }
214  else
215  {
216  z = za + log10( TempIn / atmdat.ElecTemp[iCase][nelem][ipTemp] ) * (zb-za)/(yy2-yy1);
217  }
218 
219  return ( exp10(z) );
220 }
t_atmdat atmdat
Definition: atmdat.cpp:6
static double x2[63]
double exp10(double x)
Definition: cddefines.h:1383
static double x1[83]
long int ncut[2][HS_NZ]
Definition: atmdat.h:337
long int ntemp[2][HS_NZ]
Definition: atmdat.h:335
#define NLINEHS
Definition: atmdat.h:266
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:858
long int nDensity[2][HS_NZ]
Definition: atmdat.h:335
#define EXIT_FAILURE
Definition: cddefines.h:168
#define cdEXIT(FAIL)
Definition: cddefines.h:484
double Density[2][HS_NZ][NHSDIM]
Definition: atmdat.h:329
double Emiss[2][HS_NZ][NHSDIM][NHSDIM][NLINEHS]
Definition: atmdat.h:329
double atmdat_HS_caseB(long int iHi, long int iLo, long int iZ, double TempIn, double DenIn, char chCase)
#define HS_NZ
Definition: atmdat.h:267
#define ASSERT(exp)
Definition: cddefines.h:617
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
double ElecTemp[2][HS_NZ][NHSDIM]
Definition: atmdat.h:329