@@ -65,7 +65,7 @@ class WaasKirfFormula
6565 }
6666
6767 // methods
68- double atstol (double stol) const
68+ double xrayatstol (double stol) const
6969 {
7070 double rv = c;
7171 for (int i = 0 ; i < WKTerms; ++i)
@@ -75,6 +75,23 @@ class WaasKirfFormula
7575 return rv;
7676 }
7777
78+ double Z () const
79+ {
80+ double rv = c;
81+ for (int i = 0 ; i < WKTerms; ++i) rv += a[i];
82+ return round (rv);
83+ }
84+
85+ double electronatstol (double stol) const
86+ {
87+ using diffpy::mathutils::eps_eq;
88+ using diffpy::mathutils::DOUBLE_MAX;
89+ if (eps_eq (stol, 0.0 )) return DOUBLE_MAX;
90+ double rv = 0.023934 *
91+ (this ->Z () - this ->xrayatstol (stol)) / (stol * stol);
92+ return rv;
93+ }
94+
7895 // data
7996 string symbol;
8097 double a[WKTerms];
@@ -164,6 +181,36 @@ const SetOfWKFormulas& getWKFormulasSet()
164181 return getWKFormulasSet ();
165182}
166183
184+
185+ SetOfWKFormulas::const_iterator findWKFormula (const string& smbl)
186+ {
187+ using diffpy::srreal::atomBareSymbol;
188+ using diffpy::srreal::atomValence;
189+ WaasKirfFormula wk;
190+ wk.symbol = smbl;
191+ const SetOfWKFormulas& swk = getWKFormulasSet ();
192+ SetOfWKFormulas::const_iterator wkit = swk.find (wk);
193+ if (wkit != swk.end ()) return wkit;
194+ // try again with an adjusted charge suffix
195+ wk.symbol = atomBareSymbol (smbl);
196+ int v = atomValence (smbl);
197+ if (v == 0 ) wkit = swk.find (wk);
198+ else if (abs (v) == 1 )
199+ {
200+ wk.symbol += ((v > 0 ) ? " 1+" : " 1-" );
201+ wkit = swk.find (wk);
202+ }
203+ // raise exception if no match
204+ if (wkit == swk.end ())
205+ {
206+ string emsg (" Unknown atom or ion symbol '" );
207+ emsg += smbl + " '." ;
208+ throw invalid_argument (emsg);
209+ }
210+ return wkit;
211+ }
212+
213+
167214// Electron numbers for elements and ions
168215
169216typedef boost::unordered_map<string,int > ElectronNumberStorage;
@@ -305,51 +352,25 @@ namespace srreal {
305352// / X-ray scattering factor of an element or ion a given Q
306353double fxrayatq (const string& smbl, double q)
307354{
308- double stol = q / (4 * M_PI);
355+ const double stol = q / (4 * M_PI);
309356 return fxrayatstol (smbl, stol);
310357}
311358
312359
313360// / X-ray scattering factor of an element or ion a given sin(theta)/lambda
314361double fxrayatstol (const string& smbl, double stol)
315362{
316- WaasKirfFormula wk;
317- wk.symbol = smbl;
318- const SetOfWKFormulas& swk = getWKFormulasSet ();
319- SetOfWKFormulas::const_iterator wkit = swk.find (wk);
320- if (wkit == swk.end ())
321- {
322- wk.symbol = atomBareSymbol (smbl);
323- int v = atomValence (smbl);
324- if (v == 0 ) wkit = swk.find (wk);
325- if (abs (v) == 1 )
326- {
327- wk.symbol += ((v > 0 ) ? " 1+" : " 1-" );
328- wkit = swk.find (wk);
329- }
330- }
331- // raise exception if no match
332- if (wkit == swk.end ())
333- {
334- string emsg (" Unknown atom or ion symbol '" );
335- emsg += smbl + " '." ;
336- throw invalid_argument (emsg);
337- }
338- double rv = wkit->atstol (stol);
339- return rv;
363+ SetOfWKFormulas::const_iterator wkit = findWKFormula (smbl);
364+ return wkit->xrayatstol (stol);
340365}
341366
342367
343368// / Electron scattering factor of an element or ion a given Q
344369double felectronatq (const string& smbl, double q)
345370{
346- using namespace diffpy ::mathutils;
347- // resolve Z first so that invalid symbol can throw an exception
348- double Z = round (fxrayatstol (smbl, 0.0 ));
349- if (eps_eq (q, 0.0 )) return DOUBLE_MAX;
350- double stol = q / (4 * M_PI);
351- double rv = 0.023934 * (Z - fxrayatstol (smbl, stol)) / (stol * stol);
352- return rv;
371+ const double stol = q / (4 * M_PI);
372+ SetOfWKFormulas::const_iterator wkit = findWKFormula (smbl);
373+ return wkit->electronatstol (stol);
353374}
354375
355376
0 commit comments