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>

main()

{

typedef mac::plane P;

const double a = 2;

const double b = 3;

const double c = 4;

const double d = 5;

P p(a, b, c, d); // general form

cout << p << endl;

assert (!p.poorDefined());

assert (P(TEENY).poorDefined());

// P(TEENY) is a poorly defined plane.

// assert (!P(0)); // P(0) is not a plane.

assert (~p != p);

assert (~p && p);

assert ((~p).normal() == -p.normal()); // flipping side

assert ((~p).orientationCosine() == -p.pcos());

assert (mac::approx((~p).distance(), -p.distance()));

assert (P(mac::vector3D<double>(-1, -1, 0), p.normal()) == p); // equal

// point-normal form

assert (P(p.pcos(), p.distance()) &= p); // approximately equal

// normal form

assert (P(mac::vector3D<double>( a,  b,  c)) || // parallel

P(mac::vector3D<double>(-a, -b, -c)));

// intercept form.

assert (P(mac::vector3D<double>(1, 0, 0), 0) *

P(mac::vector3D<double>(0, 1, 0), 0));

// approximately perpendicular

assert ((p ^ ~p) == M_PI);

const mac::line l = p.intersection(P());

assert ((l || p)   && mac::approx(p.distance(l.getPoint())));

assert ((l || P()) && mac::approx(P().distance(l.getPoint())));

const int n = 100;

const double amp = 2;

std::vector<mac::vector3D<double> > data;

for (int i = 0; i < n; i++)

for (int j = 0; j < n; j++)

data.push_back(mac::vector3D<double>(i, j, a * i + b * j + c) +

mac::vector3D<double>(mac::rand_fractional() - 0.5,

mac::rand_fractional() - 0.5,

mac::rand_fractional() - 0.5) *

amp);

p = P(data);

cout << p << endl;

return 0;

}

Listing 4.3.1.0.1 Testing class mac::plane.

4.3.2 Class cpl.mac.geo.plane

In this section, p and pl are planes.  v1, v2, v3, v and n are vector3DDBLs.  a, b, c, and d are float values.

plane(a = 1, b = 0, c = 0, d = 0)

plane(v, n)

plane(v)

plane(n, d)

plane(v1, v2, v3)

plane(p)

These are constructors of general form, point-normal form, intercept form, normal form, three-point form, and copy constructor, respectively.  See 4.3.1 for detail.

plane(data)

Constructs a plane least-squares fitted to the given data.  data is an object of cpl.mac.vm.vectorV3DDBL.

p.getA()

p.getB()

p.getC()

p.getD()

Return coefficients of the plane equation.

p.x(y, z)

p.y(z, x)

p.z(x, y)

Return coordinates of a point on the plane, respectively, given other two coordinates of the point.

p.is_zero()

p.poorDefined()

Return true if p is not a plane or is a poorly defined plane.

p.normal()

Returns plane normal of p.

p.pcos()

p.orientationCosine()

Return orientation cosine of plane normal.

p.orientationAngle()

Returns orientation angles in radian of plane normal.

p.distance(v)

Returns the distance from v to p.  The returned distance could be negative.  See 4.3.1 for detail.

~p

Returns flipped plane.  p remains.

p == pl

p != pl

Return true if p is equal or not equal to pl, respectively.  See 4.3.1 for detail.

p.ap_eq(pl)

Returns true if p is approximately equal to pl.

p.parallel(pl)

p | pl

Return true if p is parallel or approximately parallel to pl.

p.antiparallel(pl)

p & pl

Return true if p is anti-parallel or approximately anti-parallel to pl.

p.perpendicular(pl)

p + pl

Return true if p is perpendicular to pl.

p * pl

Returns true if p is approximately perpendicular to pl.

p.cosAngle(pl)

Returns cosine of inter-plane angle.

p.angle(pl)

p ^ pl

Return inter-plane angle in radian.

p.intersection(pl)

Returns a straight line cpl.mac.geo.line at the intersection of p and pl.

p.__str__()

Prints equation of p.

p.__nonzero__()

Returns true if p is a valid plane.

p.__hash__()

Returns a long integer for hash key.

The following testing program prints an equation of the plane if everything is OK.

#!/usr/local/bin/python

import mac

P = mac.geo.plane

V = mac.vm.vector3DDBL

a = 2.

b = 3.

c = 4.

d = 5.

p = P(a, b, c, d) # general form

print p

assert not p.poorDefined()

assert P(mac.utility.TEENY).poorDefined()

# This is a poorly defined plane.

# assert P(0).is_zero() # P(0) is not a plane.

assert ~p != p

assert (~p).normal() == -p.normal() # flipping side

assert (~p).orientationCosine() == -p.pcos()

assert mac.utility.approx((~p).distance(), -p.distance())

assert P(V(-1, -1, 0), p.normal()) == p # equal

# point-normal form

assert P(p.pcos(), p.distance()).ap_eq(p) # approximately equal

# normal form

assert P(V( a,  b,  c)) | P(V(-a, -b, -c)) # parallel

# intercept form.

assert P(V(1, 0, 0), 0) * P(V(0, 1, 0), 0)

# approximately perpendicular

assert (p ^ ~p) == mac.utility.M_PI

l = p.intersection(P());

assert (l | p)   and mac.utility.approx(p.distance(l.getPoint()))

assert (l | P()) and mac.utility.approx(P().distance(l.getPoint()))

n = 100

amp = 2

data = mac.vm.vectorV3DDBL()

for i in range(n):

for j in range(n):

x = i + amp * (mac.utility.rand_fractional() - 0.5)

y = j + amp * (mac.utility.rand_fractional() - 0.5)

z = a * i + b * j + c + amp * (mac.utility.rand_fractional() - 0.5)

data.append(mac.vm.vector3DDBL(x, y, z))

print P(data)

Listing 4.3.2.0.1 Testing class cpl.mac.geo.plane.

4.4 Straight Line

Straight line, either in 3-dimensional space or on a 2-dimensional plane, is no doubt very useful.  MAC includes two classes that model these types of straight line, which contain some common operations.

4.4.1 Straight line in 3-dimensional space

Straight line in 3-dimensional space can be defined by two elements, an orientation vector and a point on the line.  The former must be a vector with finite length, otherwise the orientation of the straight line is not defined.  If the length is too small to be practically represented by a machine, we say the straight line is poorly defined.  The latter can also be thought as a vector in this 3-dimensional space.

4.4.1.1 Class mac::line

This model of straight line in fact is more than a straight line, since it has a direction associate with a simple straight line.  It is more like a straight line with an arrow head.

#### Constructors, destructor, and assignment

line(const mac::vector3D<double>& o =

mac::vector3D<double>(1, 0, 0),

const mac::vector3D<double>& p =

mac::vector3D<double>());

Constructs a straight line which passes point p and whose orientation is o, where o must be a vector with non-zero length.  The default constructor returns a straight line of the x-axis.

line(const mac::vector3D<double>&,

const mac::vector3D<double>&, const int);

Constructs a straight line which passes two distinct points and whose orientation is from the first point to the second.  This is the two-point form of constructor.  Because of the overloading rule, the last int is necessary for distinguishing this form with the previous, general form.  The int can be any value.

line(const std::vector<mac::vector3D<double> >&);

Constructs a straight line least-squares fitted to the data given.

line(const line&);

Copy constructor.

~line();

Destructor.

const line &operator=(const line&);

Assignment operator.

#### Properties

bool operator!()   const;

bool poorDefined() const;

Return true if the object constructed is not a straight line or is a poorly defined straight line, respectively.

mac::vector3D<double> getOrientation() const;

mac::vector3D<double> getPoint()       const;

Return the orientation vector of and a point on the straight line, respectively.

mac::vector3D<double> lcos() const;

mac::vector3D<double> orientationAngle() const;

Return the orientation cosine or angles in radian of the straight line, respectively.

line operator~() const;

Returns a straight line that coincides with the original one, but with its orientation flipped towards the opposite direction.

#### Equality and Inequality

bool operator==(const line&) const;

bool operator&=(const line&) const;

bool operator!=(const line&) const;

Return true if two straight lines are equal, approximately equal, or not equal to each other, respectively.  If two straight lines coincide with each other and they have the same orientation, they are equal.  If two straight lines have nearly the same orientation, and the distance between them are very small, they are approximately equal.  See 2.1.1 and 3.1.1.19 for detail.

#### Parallelism and perpendicularity

bool operator||(const line&) const;

bool parallel(const line&) const;

Return true if two straight lines are parallel to each other.  The orientations of these lines must be towards the same direction, if they are parallel.  If they are opposite, these lines are said to be anti-parallel.  See below.

bool operator|(const line&) const;

Returns true if two straight lines are approximately parallel to each other.

bool operator&&(const line&) const;

bool antiparallel(const line&) const;

Return true if two straight lines are anti-parallel.

bool operator&(const line&) const;

Returns true if two straight lines are approximately anti-parallel.

bool operator||(const mac::plane&) const;

bool parallel(const mac::plane&) const;

Return true if the straight line is parallel to the plane.

bool operator|(const mac::plane&) const;

Returns true if the straight line is approximately parallel to the plane.

bool operator+(const line&) const;

bool perpendicular(const line&) const;

Return true if two straight lines are perpendicular to each other.

bool operator*(const line&) const;

Returns true if two straight lines are approximately perpendicular to each other.

bool operator+(const mac::plane&) const;

bool perpendicular(const mac::plane&) const;

Return true if the straight line is perpendicular to the plane.

bool operator*(const mac::plane&) const;

Returns true if the straight line is approximately perpendicular to the plane.

#### Distance

double distance(

const mac::vector3D<double>& =

mac::vector3D<double>()) const;

double distance(const line&) const;

Return distance from the straight line to a given point or another straight line, respectively.

#### Angle

double sinAngle(const line&) const;

double cosAngle(const line&) const;

double tanAngle(const line&) const;

Return sine, cosine, and tangent of the angle between two lines, respectively.

double angle(const line&) const;

double operator^(const line&) const;

Return angle between two straight lines in radian.

double sinAngle(const mac::plane&) const;

double cosAngle(const mac::plane&) const;

double tanAngle(const mac::plane&) const;

Return sine, cosine, and tangent of the angle between the straight line and the given plane, respectively.

double angle(const mac::plane&) const;

double operator^(const mac::plane&) const;

Return angle between the straight line and the given plane in radian.

#### Plane and straight line

mac::vector3D<double> intersection(

const mac::plane&) const;

Returns intersection of the straight line and the given plane.

bool coplanar(const line&,

mac::plane* = NULL) const;

Returns true if two straight lines are coplanar.  If so and the last pointer is valid, the plane is replaced by the coplane.

#### Output

friend ostream &operator<<(ostream&,

const line&);

Returns ostream to allow printing the equations of the straight line.

#include <assert.h>

#include <mac/utility.h>

#include <mac/vector3D.h>

#include <mac/plane.h>

#include <mac/line.h>

main()

{

typedef mac::line             L;

typedef mac::plane            P;

typedef mac::vector3D<double> V;

const V u(1, 2, 3);

const V v(4, 5, 6);

const V w(3, 6, 4);

const L l1(u, v);    // general form

const L l2(u, v, 1); // two-point form

cout << "line 1:" << endl;

cout << l1 << endl << endl;

cout << "line 2:" << endl;

cout << l2 << endl;

assert (l1 != l2);

assert (l1 && ~l1);

assert (l1 || P(v ->* u, 0));

assert (l1 +  P(u, 0));

assert (l1.distance(v)  == 0);

assert (l2.distance(v)  == 0);

assert (l1.distance(l2) == 0);

assert (l1.distance(l2) == l2.distance(l1));

assert (l1.intersection(P(w, v, V())) == v);

P p;

assert (L(u, v).coplanar(L(u, w), &p));

cout << p << endl << endl;

const int n = 100;

const double a =  5;

const double b =  7;

const double c =  3;

const double d = -6;

const double amp = 1;

std::vector<mac::vector3D<double> > data(n);

for (int i = 0; i < n; i++)

{

const double x = i + amp * (mac::rand_fractional() - 0.5);

const double y = a * x + b + amp * (mac::rand_fractional() - 0.5);

const double z = c * x + d + amp * (mac::rand_fractional() - 0.5);

data[i] = mac::vector3D<double>(x, y, z);

}

cout << L(data) << endl;

return 0;

}

Listing 4.4.1.1.1 Testing class mac::line.

4.4.1.2 Class cpl.mac.geo.line

In this section, u and v are two vectors cpl.mac.vm.vector3DDBL.  l and k are objects of this class.  p is a plane.

#### Constructors

line(u, v)

Constructs a straight line passing v with orientation u.

line(u, v, i)

Constructs a straight line passing two distinct points u and v.  The orientation of the line is from the first point to the second.  i is necessary to distinguish this constructor from the previous one, but it can be any integer.

line(data)

Constructs a straight line least-squares fitted to the given data.  data is an object of cpl.mac.vm.vectorV3DDBL.

line(l)

Constructs a straight line same as the given one l.

#### Properties

l.is_zero()

l.poorDefined()

Return true if l is not a valid straight line or is a poorly defined straight line.

l.__nonzero__()

Returns true if l is a valid straight line.

l.getOrientation()

l.getPoint()

Return the orientation vector of and a point on l, respectively.

l.lcos()

l.orientationAngle()

Return the orientation cosine and angles of l, respectively.

~l

Returns a straight line that coincides with l, but with its orientation flipped towards the opposite direction.

#### Equality and inequality

l == k

l != k

l.ap_eq(k)

Return true of l and k are equal, not equal, or approximately equal to each other.

#### Parallelism and perpendicularity

l.parallel(k)

l | k

Return true if l is parallel or approximately parallel to straight line k, respectively.

l.parallel(p)

l | p

Return true if l is parallel or approximately parallel to plane p, respectively.

l.antiparallel(k)

l & k

Return true if l is anti-parallel or approximately anti-parallel to k, respectively.

l.perpendicular(k)

l + k

Return true if l and k are perpendicular to each other.

l * k

Returns true if l and k are approximately perpendicular to each other.

l.perpendicular(p)

l + p

Return true if l is perpendicular to plane p.

l * p

Returns true if l is approximately perpendicular to plane p.

#### Distance

l.distance(v)

l.distance(k)

Return distance to a given point v or another straight line k.

#### Angle

l.sinAngle(k)

l.cosAngle(k)

l.tanAngle(k)

Return sine, cosine, and tangent of the angle between l and k.

l.angle(k)

l ^ k

Return the angle in radian between l and k.

l.sinAngle(p)

l.cosAngle(p)

l.tanAngle(p)

Return sine, cosine, and tangent of the angle between the straight line l and plane p.

l.angle(p)

l ^ p

Return the angle in radian between the straight line l and plane p.

#### Plane and straight line

l.intersection(p)

Returns the intersection of line l with given plane p as a vector.

l.coplanar(k, p = None)

Returns true if straight lines l and k are coplanar.  If so and the given p will be replaced by the coplane.

#### Output

l.__str__()

print l prints the equations of straight line l.

#!/usr/local/bin/python

import mac

L = mac.geo.line

P = mac.geo.plane

V = mac.vm.vector3DDBL

u = V(1, 2, 3)

v = V(4, 5, 6)

w = V(3, 6, 4)

l1 = L(u, v)    # general form

l2 = L(u, v, 1) # two-point form

print "line 1:"

print l1

print

print "line 2:"

print l2

assert l1 != l2

assert l1 & ~l1

assert l1 | P(v ** u, 0)

assert l1 + P(u, 0)

assert l1.distance(v)  == 0

assert l2.distance(v)  == 0

assert l1.distance(l2) == 0

assert l1.distance(l2) == l2.distance(l1)

assert l1.intersection(P(w, v, V())) == v

p = P()

assert l1.coplanar(L(u, w), p)

print p

n = 100

a =  5

b =  7

c =  3

d = -6

amp = 1

data = mac.vm.vectorV3DDBL(n)

for i in range(n):

data[i] = mac.vm.vector3DDBL(i, a * i + b, c * i + d) +\

mac.vm.vector3DDBL(mac.utility.rand_fractional() - 0.5,\

mac.utility.rand_fractional() - 0.5,\

mac.utility.rand_fractional() - 0.5) * amp

print

print L(data)

Listing 4.4.1.2.1 Testing class cpl.mac.geo.line.

4.4.2 Straight line in 2-dimensional space

Straight line in 2-dimensional space, that is, a plane, is a special case of straight line in 3-dimensional space.  A class mac::line2D is a derived class from mac::line.  A number of additional member functions special for the 2-dimensional straight line are included.  All functions available to mac::line are also available to mac::line2D, except those overridden.  For example, operators for comparison, parallelism, perpendicularity, and functions for angle between lines are inherited.  The function coplanar() always returns true for two mac::line2Ds.

4.4.2.1 Class mac::line2D

#### Constructors, destructor, and assignment

line2D(const double, const double, const double);

This is the constructor in general form that takes a, b, and c to construct a straight line ax + by + c = 0.

line2D(const mac::vector3D<double>& =

mac::vector3D<double>(1));

This is the constructor in general form too that takes a vector to construct a straight line.  This is also the default constructor that returns y-axis.

line2D(const double k, const double b);

Constructs a straight line y = kx + b with slope k and intercept b on y-axis.  This is a constructor in slope-intercept form.

line2D(const double,

const std::complex<double>&);

Constructs a straight line with its normal forming an angle as specified by the given double from the x-axis, and passing through a point given by the complex.  This is the normal-point form.

line2D(const std::complex<double>&,

const double);

Constructs a straight line with slope equal to the given double, and passing through a point given as the complex.  These are constructors in point-slope form.

line2D(const std::complex<double>& i);

Constructs a straight line with intercepts on x- and y-axes being the real and imaginary components of the given complex, respectively.  This is called intercept form.

line2D(const std::complex<double>&,

const std::complex<double>&);

Constructs a straight line passing through two points as the given complexes.  This is the two-point form constructor.

line2D(const double, const double, int)

Constructs a straight line with its normal forming an angle from x-axis as the first double in radian.  The distance from the origin to the straight line is given by the second double.  The last int is necessary to distinguish this constructor from the slope-intercept form.  This is called normal form.

line2D(const std::vector<

std::complex<double> >&);

Constructs a straight line least-squares fitted to the given points.  If the length of the vector is less than 2, an assertion error will occur.

line2D(const line2D&);

Copies the given straight line to a new one.

~line2D();

Destructor.

const line2D &operator=(const line2D&);

Assignment operator.

#### Properties

bool operator!() const;

Returns true if the object is not a valid straight line.

double getA() const;

double getB() const;

double getC() const;

Return the coefficients a, b, and c, respectively.

double slope() const;

Returns the slope.

std::complex<double> intercepts() const;

Returns the intercepts on x- and y-axes as the real and imaginary components of the complex, respectively.

double normalAngle() const;

Returns the angle in radian between the line normal and x-axis.

#### Point and lines

double distance(

const std::complex<double>& =

std::complex<double>()) const;

Returns distance from the given point to the straight line.  If no point is given, returns the distance to the origin.

std::complex<double> intersection(

const line2D&) const;

static std::complex<double> intersection(

const line2D&, const line2D&);

Return intersection of two straight lines.  If two lines are parallel to each other, an assertion error will occur.

bool concurrent(const line2D&, const line2D&);

static bool concurrent(const line2D&,

const line2D&,

const line2D&);

Return true if three straight lines are concurrent.

#### Output

friend ostream &operator<<(ostream&,

const line2D&);

Returns ostream for printing.

#include <assert.h>

#include <complex.h>

#include <mac/utility.h>

#include <mac/line2D.h>

main()

{

typedef mac::line2D L;

const L l1(1., 2., 3.);    // general form

const L l2(3., 4., 1.);

cout << "line 1:" << endl;

cout << l1 << endl << endl;

cout << "line 2:" << endl;

cout << l2 << endl;

assert (l1 != l2);

assert (l1 && ~l1);

const int n = 100;

const double a = 5;

const double b = 7;

std::vector<std::complex<double> > data(n);

const double amp = 30;

for (int i = 0; i < n; i++)

data[i] = std::complex<double>(i, a * i + b) +

std::complex<double>(mac::rand_fractional() - 0.5,

mac::rand_fractional() - 0.5) * amp;

const L l(data);

cout << l << endl;

cout << l.distance() << endl;

cout << l.slope() << endl;

for (int i = 0; i < n; i++)

{

const double x = real(data[i]);

const double y = (-l.getC() - l.getA() * x) / l.getB();

cout << x << " " << imag(data[i]) << " " << y << endl;

}

return 0;

}

Listing 4.4.2.1.1 Testing class mac::line2D.

4.4.2.2 Class cpl.mac.geo.line2D

In this section, a, b, c, and k are floats; p and q are complexes; v is a vector3DDBL; l, l1, and l2 are objects of this class.

#### Constructors

line2D(a, b, c)

Constructs a straight line ax + by + c = 0.

line2D(v = cpl.mac.vm.vector3DDBL(1))

Constructs a straight line with the coefficients collected in the vector.  This is also the default constructor that returns y-axis.

line2D(k, b)

Constructs a straight line y = kx + b with slope k and intercept b on y-axis.

line2D(k, p)

Constructs a straight line with its normal forming an angle k from the x-axis, and passing through a point p.

line2D(p, k)

Construct a straight line with slope k and passing through a point p.

line2D(p)

Constructs a straight line with intercepts on x- and y-axes being the real and imaginary components of the given complex, respectively.

line2D(p, q)

Constructs a straight line passing through two points.

line2D(a, b, i)

Constructs a straight line with its normal forming an angle in radian a from x-axis.  The distance from the origin to the straight line is given by b.  The last argument i is an integer necessary to distinguish this constructor from the slope-intercept form.

line2D(data)

Constructs a straight line least-squares fitted to the given points in data.  If the length of the vector is less than 2, an assertion error will occur.  data is an object of cpl.std.vectorDCX.

line2D(l)

Copies the given straight line to a new one.

#### Properties

l.is_zero()

l.__nonzero__()

Return true and false if l is not a valid straight line, respectively.

l.getA()

l.getB()

l.getC()

Return the coefficients a, b, and c, respectively.

l.slope()

Returns the slope.

l.intercepts()

Returns the intercepts on x- and y-axes as the real and imaginary components of the complex, respectively.

l.normalAngle()

Returns the angle in radian between the line normal and x-axis.

#### Point and lines

l.distance(p = complex())

Returns distance from the given point p to the straight line.  If no point is given, returns the distance to the origin.

l1.intersection(l2)

line2D.intersection(l1, l2)

line2D.inter(l1, l2)

line2D_inter(l1, l2)

Return intersection of two straight lines l1 and l2.  If two lines are parallel to each other, an assertion error will occur.

l.concurrent(l1, l2)

line2D.concurrent(l, l1, l2)

line2D.conc(l, l1, l2)

line2D_conc(l, l1, l2)

Return true if three straight lines are concurrent.

#### Output

l.__str__()

print l prints the equation of straight line.

#!/usr/local/bin/python

import cpl

L = cpl.mac.geo.line2D

l1 = L(1., 2., 3.)

l2 = L(3., 4., 1.)

print "line 1:"

print l1

print

print "line 2:"

print l2

assert l1 != l2

assert l1 & ~l1

n    = 100

a    = 5

b    = 7

amp  = 1

data = cpl.std.vectorDCX(n)

for i in range(n):

data[i] = complex(i, a * i + b) +\

complex(cpl.mac.utility.rand_fractional() - 0.5,\

cpl.mac.utility.rand_fractional() - 0.5) * amp

print L(data)

Listing 4.4.2.2.1 Testing cpl.mac.geo.line2D.

4.5 Quadric Curves

Straight line has a linear equation.  Quadric curves have equations in general form of Ax2 + Bxy + Cy2 + Dx + Ey + F = 0.  This equation can represent circle, ellipse, parabola, hyperbola, or some imaginary curves.  These curves are also known as conic curves, since they are intersections of a cone with a plane.  The most commonly used quadric curve is certainly circle.  In Laue crystallography, centric reciprocal lattice zones are transformed into some of the conic curves on a plane detector, e.g., ellipse.

4.5.1 The base class

The base classes mac::quadricCurve and cpl.mac.geo.quadricCurve are not abstract.  One may instantiate them in order to model a conic curve if the specific type is unknown.  These base classes contain the common functionalities.

4.5.1.1 Class mac::quadricCurve

#### Constructors, destructor, and assignment

quadricCurve(const double = 1,

const double = 0,

const double = 1,

const double = 0,

const double = 0,

const double = 0);

Constructs a quadric curve from the coefficients A to F.

quadricCurve(const quadricCurve&);

Copy constructor.

virtual ~quadricCurve();

Virtual destructor.

const quadricCurve& operator=(

const quadricCurve&);

Assignment operator.

#### Coefficients and comparison

double getA() const;

double getB() const;

double getC() const;

double getD() const;

double getE() const;

double getF() const;

Return coefficients A to F, respectively.

bool operator==(const quadricCurve&) const;

bool operator&=(const quadricCurve&) const;

bool operator!=(const quadricCurve&) const;

Return true if two quadric curves are equal, approximately equal, or not equal to each other, respectively.

#### Locus type

std::string type() const;

Returns a string that indicates the type of the quadric curve.  All available types are listed below.

 Type String ellipse ellipse circle circle parabola parabola hyperbola hyperbola imaginary ellipse imaginaryEllipse conjugate complex intersecting lines conjugateComplexIntersectingLines intersecting lines intersectingLines conjugate complex parallel lines conjugateComplexParallelLines distinct parallel lines distinctParallelLines coincident lines coincidentLines single line singleLine nothing nothing

Table 4.5.1.1.1 Types of quadric curve.

bool operator!() const;

Returns true if the constructed object represents nothing.

#### Rotation and translation

The general equation of a quadric curve may represent the curve at any orientation.  At a specific orientation, the cross term of the equation Bxy can be eliminated.  This orientation can be thought as a home orientation.  An arbitrary orientation is related to home orientation by a rotation.  Similarly, a specific home position of a quadric curve allows elimination of the linear terms Dx and Ey.  An arbitrary position is related to home position by a translation.  Several functions are related to the rotation and translation.

void setRotation(const double);

void setTranslation(const std::complex<double>&);

These two functions set the rotation and translation of a quadric curve to the specified values, respectively.  Thus the original rotation and translation of the curve will be lost after calling these functions.

double getRotation() const;

std::complex<double> getTranslation() const;

Return rotation and translation of a quadric curve, respectively.

quadricCurve eliminateCross();

quadricCurve eliminateLinear();

quadricCurve eliminate();

Eliminate the cross term, linear terms, or both cross and linear terms by appropriate rotation and translation of a quadric curve, respectively.  The original curve is not only modified, but also returned.

quadricCurve operator+()  const;

quadricCurve operator-()  const;

quadricCurve operator++() const;

Eliminate the cross term, linear terms, or both cross and linear terms by appropriate rotation and translation of a quadric curve, respectively.  The modified curve is returned, however, the original curve remains.

static quadricCurve eliminateCross(

const quadricCurve&);

static quadricCurve eliminateLinear(

const quadricCurve&);

static quadricCurve eliminate(

const quadricCurve&);

Eliminate the cross term, linear terms, or both cross and linear terms of the given quadric curve, and return the modified curve.

quadricCurve rotateCoordinates(const double);

quadricCurve translateCoordinates(

const std::complex<double>&);

Rotates and translates coordinate system by the specified values, respectively, and return quadric curve after the rotation or translation.

quadricCurve rotateQuadricCurve(const double);

quadricCurve translateQuadricCurve(

const std::complex<double>&);

Rotates and translates quadric curve by the specified values, respectively, and return them after the rotation or translation.

int intersection(

const mac::line2D&,

std::complex<double>& = root1,

std::complex<double>& = root2) const;

Finds intersections of the given straight line with the quadric curve, and replaces the given complexes by the solutions found, if any.  The returned int is the number of intersections found.  The defaults root1 and root2 are static data members of this class that allow one or no complex given.

void intersect(const mac::line2D&,

std::complex<double>*,

std::complex<double>*) const;

Same as above except no number of intersections returned.  This function is deprecated by the above.

#include <assert.h>

#include <mac/quadricCurve.h>

main()

{

typedef mac::quadricCurve QC;

cout << QC() << " is " << QC().type() << endl; // Default constructor

QC qc(1, 2, 3, 4, 5, 6);                       // General constructor

cout << qc << " is " << qc.type() << endl;

cout << +qc << endl;

cout << -qc << endl;

cout << ++qc << endl;

assert (qc.intersection(mac::line2D()) == 0);

// No intersection with imaginary curve.

return 0;

}

Listing 4.5.1.1.1 Testing the base class of quadric curve.

4.5.1.2 Class cpl.mac.geo.quadricCurve

In this section, r, A, B, C, D, E, and F are floats; t, p1, and p2 are complexes; l is a straight line; qc and curve are objects of quadric curve.

#### Constructors

quadricCurve(A = 1, B = 0, C = 1,

D = 0, E = 0, F = 0)

Constructs a quadric curve from the coefficients A to F.

quadricCurve(qc)

Copy constructor.

#### Coefficients and comparison

qc.getA()

qc.getB()

qc.getC()

qc.getD()

qc.getE()

qc.getF()

Return coefficients A to F, respectively.

qc == curve

qc.ap_eq(curve)

qc != curve

Return true if two quadric curves are equal, approximately equal, or not equal to each other, respectively.

#### Locus type

qc.type()

Returns a string that indicates the type of the quadric curve.  All available types are listed in Table 4.5.1.1.1.

qc.is_zero()

Returns true if qc represents nothing.

#### Rotation and translation

qc.setRotation(r);

qc.setTranslation(t);

These two functions set the rotation and translation of a quadric curve to the specified values, respectively.

qc.getRotation()

qc.getTranslation()

Return rotation and translation of a quadric curve, respectively.

qc.eliminateCross();

qc.eliminateLinear();

qc.eliminate();

Eliminate the cross term, linear terms, or both cross and linear terms, respectively.  qc is modified after calling these functions, and the modified curve is also returned.

+qc

-qc

~qc

Eliminate the cross term, linear terms, or both cross and linear terms, respectively.  The modified curve is returned, however, qc remains unchanged.

quadricCurve.elim_Cross(qc)

quadricCurve_elim_Cross(qc)

quadricCurve.elim_Linear(qc)

quadricCurve_elim_Linear(qc)

quadricCurve.elim(qc)

quadricCurve_elim(qc)

Eliminate the cross term, linear terms, or both cross and linear terms of qc, and return the modified curve.

qc.rotateCoordinates(r)

qc.translateCoordinates(t)

Rotates and translates coordinate system by the specified values, respectively, and return quadric curve after the rotation or translation.

qc.rotateQuadricCurve(r)

qc.translateQuadricCurve(t)

Rotates and translates quadric curve by the specified values, respectively, and return them after the rotation or translation.

qc.intersection(l)

Finds intersections of the given straight line l with the quadric curve qc, and returns a tuple.  The first element of the tuple is the number of intersections found.  The second and third are the intersections.

#!/usr/local/bin/python

import mac

QC = mac.geo.quadricCurve

print QC(), "is", QC().type() # Default constructor

qc = QC(1, 2, 3, 4, 5, 6)     # General constructor

print qc, "is", qc.type()

assert (+qc).getB() == 0

assert (-qc).getD() == (-qc).getE() == 0

assert qc != ~qc

assert qc.eliminate() == qc

assert qc.intersection(mac.geo.line2D())[0] == 0

# No intersection with imaginary curve.

Listing 4.5.1.2.1 Testing cpl.mac.geo.quadricCurve.

4.5.2 Circle

4.5.2.1 Class mac::circle

circle(const double = 1,

const std::complex<double>& =

std::complex<double>());

circle(const std::complex<double>&,

const double);

Construct a circle centered at the given complex with radius as the given double.  The first constructor is also the default constructor.

circle(const std::complex<double>&,

const std::complex<double>&);

Constructs a circle with one of its diameter located from and to the given complexes.

circle(const std::complex<double>&,

const std::complex<double>&,

const std::complex<double>&);

Constructs a circle through three points as the given complexes.  If all points are collinear, the returned object is a straight line.

circle(const circle&);

Copy constructor.

virtual ~circle();

Destructor.

const circle& operator=(const circle&);

Assignment operator.

void setRadius(const double);

double getRadius() const;

Sets and gets radius, respectively.

void center(const std::complex<double>&);

std::complex<double> center() const;

Sets and gets center, respectively.

double area() const;

Returns area.

int tangent(const std::complex<double>&,

mac::line2D&,

mac::line2D&) const;

Replaces the straight lines by tangent lines of the circle through the given point.  The returned int is the number of tangent lines found.  Due to floating operation, it may be rare to find one and only one tangent when the given point is on the circle.

int intersection(

const mac::circle&,

std::complex<double>& = root1,

std::complex<double>& = root2) const;

Replaces the complexes by intersecting points with another given circle.  The returned int is the number of intersection found.  If two circles are identical, the number of intersections is -1.

int intersections(std::complex<double>*,

std::complex<double>*,

const double = 1) const;

Same as above except that the intersecting circle is centered at the origin.  This function is deprecated by the above.

#include <assert.h>

#include <mac/utility.h>

#include <mac/quadricCurve.h>

main()

{

typedef mac::circle C;

const std::complex<double> p1(1, 2), p2(3, 4), p3(6, 5);

cout << C() << " is " << C().type() << endl; // Default constructor

const C c1(p1, p2, p3);

assert (c1.type() == mac::quadricCurve::Circle);

const C c2(p1, p2);

std::complex<double> q1, q2;

assert (c1.intersection(c2, q1, q2) == 2);

assert (mac::approxComplex(p1, q1) || mac::approxComplex(p1, q2));

assert (mac::approxComplex(p2, q1) || mac::approxComplex(p2, q2));

mac::line2D t1, t2;

assert (c2.tangent(p3, t1, t2) == 2);

cout << t1 << endl;

cout << t2 << endl;

return 0;

}

Listing 4.5.2.1.1 Testing mac::circle.

4.5.2.2 Class cpl.mac.geo.circle

In this section, r is a float; p1, p2, and p3 are complexes; c and circle are two objects of circle.

circle(r = 1, p1 = complex())

circle(p1, r)

Construct a circle centered at the given complex p1 with radius r.  The first constructor is also the default constructor.

circle(p1, p2)

Constructs a circle with one of its diameter located from p1 to p2.

circle(p1, p1, p3)

Constructs a circle through three points p1, p2, and p3.  If all points are collinear, the returned object is a straight line.

circle(c)

Copy constructor.

c.setRadius(r)

c.getRadius()

Sets and gets radius, respectively.

c.center(p1)

c.center()

Sets and gets center, respectively.

c.area()

Returns area.

c.tangent(p1, l1, l2)

Replaces the straight lines l1 and l2 by tangent lines of the circle c through point p1, and returns the number of tangent lines found.

c.intersection(circle)

Returns a tuple containing the number of intersections found between two circles and the intersecting points.  If two circles are identical, the number of intersections is -1.

#!/usr/local/bin/python

import mac

C  = mac.geo.circle

p1 = complex(1, 2)

p2 = complex(3, 4)

p3 = complex(6, 5)

print C(), "is", C().type() # Default constructor

c1 = C(p1, p2, p3)

assert c1.type() == "circle"

c2 = C(p1, p2)

q1 = complex()

q2 = complex()

n, q1, q2 = c1.intersection(c2, q1, q2)

assert n == 2

assert mac.utility.approxComplex(p1, q1) or mac.utility.approxComplex(p1, q2)

assert mac.utility.approxComplex(p2, q1) or mac.utility.approxComplex(p2, q2)

t1 = mac.geo.line2D()

t2 = mac.geo.line2D()

assert c2.tangent(p3, t1, t2) == 2

print t1

print t2

Listing 4.5.2.2.1 Testing cpl.mac.geo.circle.

4.5.3 Ellipse

4.5.3.1 Class mac::ellipse

ellipse(const std::complex<double>& =

std::complex<double>(1, 1),

const double = 0,

const std::complex<double>& =

std::complex<double>(0, 0));

Constructs an ellipse with semi-major and -minor axes specified by the first complex.  The rotation and translation of the ellipse can be specified by the given double and the second complex.

ellipse(const double,

const double,

const double = 0,

const std::complex<double>& =

std::complex<double>());

Constructs an ellipse with its eccentricity and focus-to-directrix being the first and second doubles.  By default, the ellipse has one of its focus at the origin and one of its directrix being x = -p, where p is the second double.  Additional rotation and translation of the ellipse can be set by the third double and the complex.

ellipse(const ellipse&);

Copy constructor.

virtual ~ellipse();

Destructor;

const ellipse& operator=(const ellipse&);

Assignment operator.

void setAxes(const std::complex<double>&);

std::complex<double> getAxes() const;

Sets and gets semi-major and -minor axes as the real and imaginary components of complex.

double getSemimajor() const;

double getSemiminor() const;

Return the lengths of semi-major and minor axes, respectively.

mac::line2D majorAxis() const;

mac::line2D minorAxis() const;

Return the equation of major and minor axes, respectively.

void center(const std::complex<double>&);

std::complex<double> center() const;

Sets and gets the center, respectively.

double area() const;

double eccentricity() const;

double focalLength() const;

Return area, eccentricity, and focal length, respectively.

void foci(std::complex<double>& = root1,

std::complex<double>& = root2) const;

Returns the two foci.

void directrices(mac::line2D&,

mac::line2D&) const;

Returns the two directrices.

In addition, there are several other functions no longer supported.  They may be revised and supported in future releases.

#include <assert.h>

#include <mac/utility.h>

#include <mac/quadricCurve.h>

main()

{

typedef mac::ellipse E;

const std::complex<double> p1(1, 2), p2(3, 4), p3(6, 5);

const E ell(p1, 1, p2);

assert (ell.type() == mac::quadricCurve::Ellipse);

assert (ell.getAxes() == p1);

assert (ell.center()  == p2);

assert (ell.eccentricity() < 1);

assert (ell.majorAxis() + ell.minorAxis());

std::complex<double> q1, q2;

ell.foci(q1, q2);

assert (mac::approx(abs(q1 - q2), ell.focalLength()));

mac::line2D d1, d2;

ell.directrices(d1, d2);

assert (d1 || d2);

return 0;

}

Listing 4.5.3.1.1 Testing mac::ellipse.

4.5.3.2 Class cpl.mac.geo.ellipse

In this section, r, e, p are floats ; a and t are complexes; ell is an object of ellipse.

ellipse(a = complex(1, 1), r = 0, t = complex())

Constructs an ellipse with semi-major and -minor axes specified by the first complex a.  The rotation and translation of the ellipse can be specified by r and t, respectively.

ellipse(e, p, r = 0, t = complex())

Constructs an ellipse with its eccentricity and focus-to-directrix being e and p.  By default, the ellipse has one of its focus at the origin and one of its directrix being x = -p.  Additional rotation and translation of the ellipse can be set by r and t, respectively.

ellipse(ell)

Copy constructor.

ell.setAxes(a)

ell.getAxes()

Sets and gets semi-major and -minor axes as the real and imaginary components of complex.

ell.getSemimajor()

ell.getSemiminor()

Return the lengths of semi-major and minor axes, respectively.

ell.majorAxis()

ell.minorAxis()

Return the equation of major and minor axes as object of cpl.mac.geo.line2D, respectively.

ell.center(t)

ell.center()

Sets and gets the center, respectively.

ell.area()

ell.eccentricity()

ell.focalLength()

Return area, eccentricity, and focal length, respectively.

ell.foci()

Returns the two foci.

ell.directrices(l1, l2)

Replaces l1 and l2 by the two directrices as objects of cpl.mac.geo.line2D.

#!/usr/local/bin/python

import mac

E  = mac.geo.ellipse

p1 = complex(1, 2)

p2 = complex(3, 4)

p3 = complex(6, 5)

ell = E(p1, 1, p2)

assert ell.type() == "ellipse"

assert ell.getAxes() == p1

assert ell.center()  == p2

assert ell.eccentricity() < 1

assert ell.majorAxis() + ell.minorAxis()

q1, q2 = ell.foci()

assert mac.utility.approx(abs(q1 - q2), ell.focalLength())

d1 = mac.geo.line2D()

d2 = mac.geo.line2D()

ell.directrices(d1, d2)

assert d1 | d2

Listing 4.5.4.2.1 Testing cpl.mac.geo.ellipse.

4.5.4 Parabola

This class is still underdeveloped.

4.5.4.1 Class mac::parabola

parabola(const double = 1,

const double = 0,

const std::complex<double>& =

std::complex<double>(0, 0));

Constructs a parabola with focus-to-directrix as specified by the first double.  By default, the focus is located at the origin.  Additional rotation and translation can be set by the second double and the complex.

parabola(const parabola&);

Copy constructor.

virtual ~parabola();

Destructor.

const parabola& operator=(const parabola&);

Assignment operator.

void setFocusToDirectrix(const double);

Sets focus-to-directrix to the given double.

This class is untested.

4.5.4.2 Class cpl.mac.geo.parabola

Here, p and r are floats; t is a complex; and pd is an object of parabola.

parabola(p = 1, r = 0, t = complex())

Constructs a parabola with focus-to-directrix as specified by p.  By default, the focus is located at the origin.  Additional rotation and translation can be set by r and t, respectively.

parabola(pb)

Copy constructor.

pb.setFocusToDirectrix(p)

Sets focus-to-directrix to p.

4.5.5 Hyperbola

4.5.5.1 Class mac::hyperbola

hyperbola(const std::complex<double>& =

std::complex<double>(1, 1),

const double = 0,

const std::complex<double>& =

std::complex<double>(0, 0));

Constructs a hyperbola with semi-real and -imaginary axes as given by the real and imaginary components of the first complex.  Rotation and translation can be set by the double and the second complex, respectively.

hyperbola(const hyperbola&);

Copy constructor.

virtual ~hyperbola();

Destructor.

const hyperbola& operator=(const hyperbola&);

Assignment operator.

void setAxes(const std::complex<double>&);

Sets the semi-real and -imaginary axes.

This class is untested.

4.5.5.2 Class cpl.mac.geo.hyperbola

Here, r is a float; a and t are complexes; hb is an object of hyperbola.

hyperbola(a = complex(1, 1),

r = 0,

t = complex())

Constructs a hyperbola with semi-real and -imaginary axes as given by the real and imaginary components a.  Rotation and translation can be set by r and t, respectively.

hyperbola(hb)

Copy constructor.

hb.setAxes(a)

Sets the semi-real and -imaginary axes.

4.6 Cone