CHAPTER
4
Function, Equation, and Geometry
This chapter discusses several functions useful in crystallography. MAC library supplies these functions for other parts of CCL and CPL. When two functions are said equal to each other, an equation is formed. The condition under which the equation is satisfied is called solution of the equation. Solving equations is a common task. MAC library includes a few types of equations. Both function and equation have their geometric representation. It is often easier to think with some geometric objects in mind or to illustrate on paper. Some of these objects are coded in form of classes in MAC library, and wrapped into Python as module cpl.mac.geo. This chapter offers some more tools for crystallographic programming.
4.1 Generalized Gaussian Function
This is a generalized Gaussian function defined as below:
,
where a, b, c, d, e, and f are parameters. a is the center of the function, or the average; g(a) = f. b is the base of the exponential expression; and b ¡Ý 1. The natural Gaussian function has the same base as the natural logarithm, 2.71828¡ c affects the kurtosis (or curtosis) of the peak; and c ¡Ý 0. A normal kurtosis is reached when c = 1. d is deviation that affects the peak width; and d > 0. e makes the peak skew to one side or the other. A normal, non-skew function has e = 0. Positive e makes the peak skew to the positive side. f is a coefficient that controls the amplitude of the peak.
To make sure the parameters b, c, and d stay in their defined range, they may be replaced by 1 + exp(b), exp(c), and exp(d), respectively. Therefore, b, c, and d become the precursors of base, kurtosis, and deviation.
.
4.1.1 Generalized Gaussian function and partial derivatives to its parameters
double
gGaussian(
const double
x,
const
std::vector<double>& para,
std::vector<double> *deri = NULL,
const bool prec
= 0);
Given a vector para of parameters a to f, this function returns the function value of the generalized Gaussian at x. If the Boolean prec is set, all of the given parameters are treated as the precursors. If a valid pointer deri is given, the vector pointed by the pointer will be modified to contain partial derivatives of the function to the parameters or their precursors.
double gGaussian(
const double& x,
const
std::vector<double>&,
const
std::vector<double>&,
const
std::vector<double>& para,
const std::vector<bool>&
= std::vector<bool>(),
std::vector<double> *deri = NULL);
The previous function does not match the prototype of an evaluation function required by mac::nonlinear<Y, X> (see 12.4.1.2). This overloaded function is simply a wrapper of the previous one, so that this one can be directly used by the non-linear least-squares class. All parameters are precursors. This function has another name gGauss.
std::vector<double>
gGaussianAll(
const
std::vector<double>&,
const
std::vector<double>&,
const std::vector<double>&,
const std::vector<bool>&
= std::vector<bool>(),
std::vector<std::vector<double> >* = NULL);
This is an evaluation function required by mac::nonlinear<Y, X> (see 12.4.1.2) that calculates at all data points in a single call. All parameters are precursors.
See Listing 6.4.5.1.1 for testing of these functions.
4.1.2 Generalized Gaussian and call-back functions in Python
gGaussian(x,
p, d = None, precursor = 0)
Returns the function value of generalized Gaussian at the given x. If d is a 6-member cpl.std.vectorDBL, it contains partial derivatives of the function to all parameters at the given x after the function returns. If the Boolean precursor is set, p is a cpl.std.vectorDBL that contains the precursors of the 6 parameters of the function. If precursor is unset, p must contain the 6 parameters.
nonlinearDBLDBL_cb1
nonlinearDBLDBL_cb2
Pointers to call-back functions that can be used by class cpl.mac.op.nonlinearDBLDBL.
See 6.4.6.8 for testing of these functions. In addition, the following example also provides a test. This program calculates and prints out the function values and partial derivatives at various parameter settings.
#!/usr/bin/python
import cpl
import math
# Setting up parameters for a
normal Gaussian
p = cpl.std.vectorDBL(6)
p[0] = 0 # a, center
p[1] = math.e # b, base
p[2] = 1 # c, kurtosis
p[3] = 1 # d, deviation
p[4] = 0 # e, skew
p[5] = 1 # f, amplitude
k = 2 # Varying
kurtosis
s = 0.2 # Step
m = 100
n = 4
q = cpl.std.vectorDBL(p) #
Copy parameter set
d = cpl.std.vectorDBL(6) #
Derivatives
for j in range(-n, n + 1):
for i in range(-m, m + 1):
x = 0.1 * i
if k == 4: q[k]
= p[k] + (j + n) * s
else: q[k] = p[k] + j
* s
g = cpl.mac.geo.gGaussian(x,
q, d)
print
'%8.3f%8.3f' % (x, g),
for l in
range(6): print '%8.3f' % d[l],
print
Listing 4.1.2.0.1 Testing gGaussian function.
Using the function values printed by this program, generalized Gaussian functions are plotted in Figures 4.1.2.0.1 through 3. All curves have the following set of common parameters except those stated in the Figure captions: a = 0; b = 2.71828¡; c = 1; d = 1; e = 0; f = 1.
Figure 4.1.2.0.1 Generalized Gaussian of various bases b from 1.11828¡ to 4.31828¡ The wider curves have smaller b values.
Figure 4.1.2.0.2 Generalized Gaussian of various kurtosis c from 0.2 to 1.8. The sharper curves have smaller c values.
Listing 4.1.2.0.3 Generalized Gaussian of various skew e from 0 to 4. Larger e value makes the curve more asymmetric. Negative e value will cause a curve to skew in the opposite direction.
4.2 Quadratic, Cubic, and Quartic Equations
Quadratic, cubic, and quartic equations with a single unknown are commonly occurring problems to solve. MAC includes these trivial classes mac::quadraticEquation, mac::cubicEquation, and mac::quarticEquation to make the solution of these problems easy.
4.2.1 Quadratic equation
quadraticEquation(const double
a = 1,
const double b = 0,
const double c = 0);
Constructs a quadratic equation ax2 + bx + c = 0.
quadraticEquation(const
quadraticEquation&);
Copy constructor.
~quadraticEquation();
Nothing to do.
const quadraticEquation&
operator=(
const quadraticEquation&);
Assigns one equation to the other.
bool operator==(const
quadraticEquation&) const;
bool operator!=(const
quadraticEquation&) const;
Test whether two equations are equivalent or not. When the ratios between all coefficients of two equations are the same, that is, a1 : b1 : c1 = a2 : b2 : c2, these equations are equivalent.
void setAto1();
Simplifies the equation by setting a to 1, and keeps the equation equivalent.
double getA() const;
double getB() const;
double getC() const;
Return the coefficients.
int
roots(std::complex<double>& = root1,
std::complex<double>&
= root2) const;
Finds the roots of the equation, and returns them in the given complexes. The defaults root1 and root2 are static members of the class, which are only space holders, and allow calling this function with only one or no argument. The returned int is the number of real roots.
std::vector<std::complex<double>
> roots();
Finds the roots of the equation, and returns them in a
2-menbered vector. This is a
uniform prototype for all equation classes.
friend ostream
&operator<<(ostream&,
const quadraticEquation&);
Returns ostream for printing this equation.
All functions above are also available in the other two
classes mac::cubicEquation and mac::quarticEquation. The only differences are the number of
arguments required and obviously the equation type.
4.2.2 Cubic equation
All functions of mac::quadraticEquation are also available for cubic equation with similar prototypes.
cubicEquation(const
double a = 1,
const double b = 0,
const double c = 0,
const double d = 0);
Constructs a cubic equation.
int roots(double& = root0,
std::complex<double>& = root1,
std::complex<double>&
= root2) const;
Finds one real root and two other roots that could be complexes. The defaults root0 and root1, root2 are static members of this class, which are only space holders, and allow calling this function with only one or two or no arguments. For example, if only the real root is needed, one may use eq.roots(realroot). The returned int is the number of real roots.
std::vector<std::complex<double>
> roots();
Finds the roots of the equation, and returns them in a
3-menbered vector. The first member
must be a real root. This is a
uniform prototype for all equation classes.
4.2.3 Quartic equation
All functions of mac::quadraticEquation are also available for quartic equation with similar prototypes.
quarticEquation(const double a
= 1,
const double b = 0,
const double c = 0,
const double d = 0,
const double e = 0);
Constructs a quartic equation.
int roots(std::complex<double>&
= root1,
std::complex<double>&
= root2,
std::complex<double>&
= root3,
std::complex<double>&
= root4) const;
Finds four roots that may be complexes. The defaults are static members of the class, which are only space holders, and allow calling this function with less than four arguments. The returned int is the number of real roots.
#include <assert.h>
#include
<mac/utility.h>
#include
<mac/equation.h>
main()
{
typedef mac::quadraticEquation E2;
typedef mac::cubicEquation E3;
typedef mac::quarticEquation E4;
const double a = 2;
const double b = 3;
const double c = 4;
const double d = 5;
const double e = 6;
// Quadratic equation
E2 e2(a, b, c);
cout << e2 << endl;
std::complex<double> s1, s2;
e2.roots(s1, s2);
assert (mac::approxComplex(a * s1
* s1 + b * s1 + c));
assert (mac::approxComplex(a * s2
* s2 + b * s2 + c));
E2 eq(e2);
eq.setAto1();
assert (eq.getA() == 1);
assert (eq.getB() == b / a);
assert (eq.getC() == c / a);
assert (eq == e2);
std::complex<double> s3, s4;
eq.roots(s3); // Only need one
solution.
assert (mac::approxComplex(s3,
s1));
std::vector<std::complex<double> > s = eq.roots();
assert (mac::approxComplex(s1, s[0]));
assert (mac::approxComplex(s2,
s[1]));
// Cubic equation
E3 e3(a, b, c, d);
cout << e3 << endl;
double s0;
e3.roots(s0, s1, s2);
assert (mac::approx (a * s0 *
s0 * s0 + b * s0 * s0 + c * s0 + d));
assert (mac::approxComplex(a * s1
* s1 * s1 + b * s1 * s1 + c * s1 + d));
assert (mac::approxComplex(a * s2
* s2 * s2 + b * s2 * s2 + c * s2 + d));
E3 equ(e3);
equ.setAto1();
assert (equ.getA() == 1);
assert (equ.getB() == b / a);
assert (equ.getC() == c / a);
assert (equ.getD() == d / a);
assert (equ == e3);
double t0;
equ.roots(t0); // Only need the
real root.
assert (mac::approx(t0, s0));
s = equ.roots();
assert (mac::approx(real(s[0]),
s0)); // The first root must be real.
assert (mac::approxComplex(s[1],
s1));
assert (mac::approxComplex(s[2],
s2));
// Quartic equation
E4 e4(a, b, c, d, e);
cout << e4 << endl;
e4.roots(s1, s2, s3, s4);
assert (mac::approxComplex(a * s1
* s1 * s1 * s1
+ b * s1 * s1 * s1
+ c * s1 * s1
+ d * s1
+ e));
assert (mac::approxComplex(a * s2
* s2 * s2 * s2
+ b * s2 * s2 * s2
+ c
* s2 * s2
+ d * s2
+ e));
assert (mac::approxComplex(a * s3
* s3 * s3 * s3
+ b * s3 * s3 * s3
+ c * s3 * s3
+ d * s3
+ e));
assert (mac::approxComplex(a * s4
* s4 * s4 * s4
+ b * s4 * s4 * s4
+ c * s4 * s4
+ d * s4
+ e));
E4 equa(e4);
equa.setAto1();
assert (equa.getA() == 1);
assert (equa.getB() == b / a);
assert (equa.getC() == c / a);
assert (equa.getD() == d / a);
assert (equa.getE() == e / a);
assert (equa == e4);
std::complex<double> t1, t2,
t3, t4;
equa.roots(t1, t2, t3, t4);
assert (mac::approxComplex(t1,
s1));
assert (mac::approxComplex(t2,
s2));
assert (mac::approxComplex(t3,
s3));
assert (mac::approxComplex(t4,
s4));
s = equa.roots();
assert (mac::approxComplex(s[0],
s1));
assert (mac::approxComplex(s[1],
s2));
assert (mac::approxComplex(s[2],
s3));
assert (mac::approxComplex(s[3],
s4));
return 0;
}
Listing 4.2.0.0.1 Testing equation classes in C++.
4.2.4 Quadratic, cubic, and quartic equations in module cpl.mac.geo
In this section, a, b, c, d, e are real values; eq and equa are equations of same type.
quadraticEquation(a
= 1, b = 0, c = 0)
quadraticEquation(eq)
cubicEquation(a
= 1, b = 0, c = 0, d = 0)
cubicEquation(eq)
quarticEquation(a
= 1, b = 0, c = 0, d = 0, e = 0)
quarticEquation(eq)
Construct equations.
eq
== equa
eq
!= equa
Test whether two equations are equivalent or not. See 4.2.1.
eq.getA()
eq.getB()
eq.getC()
¡
Return the coefficients of the equation.
eq.setAto1()
Sets the coefficient a to 1 and keep the equation equivalent.
eq.roots()
Returns all roots in a tuple. The first element of the tuple is the number of real roots. If eq is a cubic equation, the second element of the tuple must be a real root. All other roots will be complexes.
eq.__str__()
print eq prints the equation.
The following listing tests these equation classes in Python. It will print three equations without other messages, if everything is OK.
#!/usr/local/bin/python
import cpl
E2 = cpl.mac.geo.quadraticEquation
E3 =
cpl.mac.geo.cubicEquation
E4 =
cpl.mac.geo.quarticEquation
a = 2.
b = 3.
c = 4.
d = 5.
e = 6.
# Quadratic equation
e2 = E2(a, b, c)
print e2
s = e2.roots()
for i in range(1, len(s)):
assert
cpl.mac.utility.approxComplex(a * s[i] * s[i] + b * s[i] + c)
eq = E2(e2)
eq.setAto1()
assert eq.getA() == 1
assert eq.getB() == b / a
assert eq.getC() == c / a
assert eq == e2
r, s1, s2 = eq.roots()
assert
cpl.mac.utility.approxComplex(s1, s[1])
assert
cpl.mac.utility.approxComplex(s2, s[2])
# Cubic equation
e3 = E3(a, b, c, d)
print e3
r, s0, s1, s2 = e3.roots();
assert
cpl.mac.utility.approx (a * s0 **
3 + b * s0 ** 2 + c * s0 + d)
assert
cpl.mac.utility.approxComplex(a * s1 ** 3 + b * s1 ** 2 + c * s1 + d)
assert
cpl.mac.utility.approxComplex(a * s2 ** 3 + b * s2 ** 2 + c * s2 + d)
equ = E3(e3)
equ.setAto1();
assert equ.getA() == 1
assert equ.getB() == b / a
assert equ.getC() == c / a
assert equ.getD() == d / a
assert equ == e3
r, t0, s3, s4 = equ.roots()
assert
cpl.mac.utility.approx(t0, s0)
assert
cpl.mac.utility.approxComplex(s3, s1)
assert
cpl.mac.utility.approxComplex(s4, s2)
# Quartic equation
e4 = E4(a, b, c, d, e)
print e4
s = e4.roots()
for i in range(1, len(s)):
assert
cpl.mac.utility.approxComplex(\
a * s[i] ** 4 +
b * s[i] ** 3 + c * s[i] ** 2 + d * s[i] + e)
equa = E4(e4)
equa.setAto1()
assert equa.getA() == 1
assert equa.getB() == b / a
assert equa.getC() == c / a
assert equa.getD() == d / a
assert equa.getE() == e / a
assert equa == e4
r, s1, s2, s3, s4 = equa.roots()
assert
cpl.mac.utility.approxComplex(s1, s[1])
assert
cpl.mac.utility.approxComplex(s2, s[2])
assert
cpl.mac.utility.approxComplex(s3, s[3])
assert
cpl.mac.utility.approxComplex(s4, s[4])
Listing 4.2.0.0.2 Testing the equation classes in Python.
4.3 Plane
A plane in 3-dimensional space is also very useful in crystallography, for example, a reciprocal lattice zone is a plane. A class mac::plane includes some basic operations of a plane, and this class is wrapped into module cpl.mac.geo.
4.3.1 Class mac::plane
plane(const
double a = 1,
const double b =
0,
const double c =
0,
const double d =
0);
This is a constructor of general form that returns a plane ax + by + cz + d = 0. This is also a default constructor.
plane(const mac::vector3D<double>&
p,
const
mac::vector3D<double>& n);
This is point-normal form of constructor. This constructor returns a plane that passes point p and whose normal is n.
plane(const
mac::vector3D<double>&);
Intercept form of constructor returns a plane that intercepts with 3 axes at x, y, z, respectively. These three values forms the required vector.
plane(const mac::vector3D<double>&
n, const double d);
Normal form constructs a plane whose normal is orientation cosine of n, and the distance from the origin to the plane is d. The distance can positive or negative. Positive distance means that the vector from the origin to the plane is along the plane normal.
plane(const
mac::vector3D<double>&,
const
mac::vector3D<double>&,
const mac::vector3D<double>&);
Three-point form of constructor returns a plane that passes three points.
plane(const
std::vector<mac::vector3D<double> >&);
Constructs a plane least-squares fitted to the given data
points.
plane(const plane&);
This is copy constructor.
~plane();
Nothing to do in destructor.
const
plane &operator=(const plane&);
Assignment operator.
double getA() const;
double getB() const;
double getC() const;
double getD() const;
Return coefficients a, b, c, and d, respectively.
double x(const double y, const double z) const;
double y(const double z, const double x) const;
double z(const double x, const double y) const;
Return coordinates x, y, z of a point on the plane, respectively, given the other two coordinates of the point.
bool
operator!() const;
bool
poorDefined() const;
Returns true if the object constructed is not a plane or a poorly defined plane, respectively.
mac::vector3D<double>
normal() const;
Returns the plane normal as a vector.
mac::vector3D<double>
pcos() const;
mac::vector3D<double>
orientationCosine() const;
Return orientation cosine of the plane normal.
mac::vector3D<double>
orientationAngle() const;
Returns orientation angles in radian of the plane normal.
double distance(const
mac::vector3D<double> =
mac::vector3D<double>())
const;
Returns distance from the given point to the plane. If the direction from the given point to the plane is along the plane normal, returns a positive distance, otherwise, a negative distance.
plane operator~() const;
Returns a plane that is superimposed on the original plane, except the plane normal is flipped to the opposite direction.
bool operator==(const
plane&) const;
bool operator&=(const
plane&) const;
bool operator!=(const
plane&) const;
Return true if two planes are equal, approximately equal, and not equal to each other, respectively.
If two planes have a same orientation cosine of their plane normal, and their distances to the origin are equal, these two planes are equal to each other. If two objects constructed by the constructors are not plane (see operator !), they are equal to each other. Two planes are not equal to each other, if none of the above is true.
If two planes have approximately the same orientation cosine of their plane normal (see 3.1.1.19 equality and inequality of vectors), and distances to the origin are approximately equal (see 2.1.1 for function approx()), these two planes are approximately equal to each other. If two planes are equal, they are approximately equal.
bool
operator||(const plane&) const;
bool
operator| (const plane&) const;
Return true if two planes are parallel or approximately parallel to each other, respectively.
bool operator&&(const
plane&) const;
bool operator& (const
plane&) const;
Return true if two planes are anti-parallel or approximately anti-parallel to each other, respectively.
bool operator+(const plane&)
const;
bool operator*(const plane&)
const;
Return true if two planes are perpendicular or approximately perpendicular to each other, respectively.
bool parallel (const plane&)
const;
bool antiparallel (const
plane&) const;
bool perpendicular(const
plane&) const;
Return true if two planes are parallel, anti-parallel, or perpendicular to each other, respectively.
double cosAngle(const plane&) const;
Returns cosine of inter-plane angle.
double angle(const plane&)
const;
double operator^(const plane&)
const;
Return inter-plane angle in radian.
mac::line intersection(const
plane&) const;
Returns a straight line (see 4.4.1) at the intersection of two
planes. If two planes are parallel,
an assertion error will happen.
friend ostream
&operator<<(ostream&, const plane&);
Prints equation of the plane.
The following testing program shall print an equation of a
plane only, if everything is OK.
#include <assert.h>
#include
<mac/utility.h>
#include
<mac/vector3D.h>
#include <mac/plane.h>
#include <mac/line.h>