-+-+-+-+-+-+-+-+ START OF PART 264 -+-+-+-+-+-+-+-+
X  ** A X B :=  sin(theta) K
X  ** The sign of sinetheta is not important
X  ** because of the way it appears in the subsequent formulae
X  */
X
X  sinetheta = ptk_modv3(&nvec);
X  if (ptk_equal(fabs(sinetheta), 0.0))`20
X  `7B
X    *error = -1;
X    sinetheta = ptkcpceps;
X  `7D
X  nvec = ptk_unitv3(&nvec);
X  ptk_unitmatrix3(temp);   /* get a unit matrix */
X  t1 = 1.0 - costheta;
X  temp`5B0`5D`5B0`5D = (nvec.x * nvec.x) + (1.0 - (nvec.x * nvec.x)) * costh
Veta;
X  temp`5B1`5D`5B1`5D = (nvec.y * nvec.y) + (1.0 - (nvec.y * nvec.y)) * costh
Veta;
X  temp`5B2`5D`5B2`5D = (nvec.z * nvec.z) + (1.0 - (nvec.z * nvec.z)) * costh
Veta;
X  temp`5B1`5D`5B0`5D = (nvec.x * nvec.y * t1) + (nvec.z * sinetheta);
X  temp`5B2`5D`5B0`5D = (nvec.x * nvec.z * t1) - (nvec.y * sinetheta);`20
X  temp`5B0`5D`5B1`5D = (nvec.x * nvec.y * t1) - (nvec.z * sinetheta);
X  temp`5B2`5D`5B1`5D = (nvec.y * nvec.z * t1) + (nvec.x * sinetheta);
X  temp`5B0`5D`5B2`5D = (nvec.x * nvec.z * t1) + (nvec.y * sinetheta);
X  temp`5B1`5D`5B2`5D = (nvec.y * nvec.z * t1) - (nvec.x * sinetheta);
X  ptk_concatenatematrix3(operation, temp, matrix, matrix);
X`7D  /* ptk_rotatevv3 */
X
X/*--------------------------------------------------------------------------
V*/
X
X/*function:external*/
Xextern void ptk_rotatevv(C(Ppoint *) v1, C(Ppoint *) v2,`20
X                          C(Pcomptype) operation,`20
X                          C(Pmatrix) matrix, C(Pint *)error)
XPreANSI(Ppoint *v1)
XPreANSI(Ppoint *v2)
XPreANSI(Pcomptype operation)`20
XPreANSI(Pmatrix matrix)
XPreANSI(Pint *error)
X/*
X** \parambegin
X** \param`7BPpoint *`7D`7Bv1`7D`7B2D vector`7D`7BIN`7D
X** \param`7BPpoint *`7D`7Bv2`7D`7B2D vector`7D`7BIN`7D
X** \param`7BPcomptype`7D`7Boperation`7D`7Bconcatenation operation`7D`7BIN`7D
X** \param`7BPmatrix`7D`7Bmatrix`7D`7B3x3 matrix`7D`7BOUT`7D
X** \param`7BPint *`7D`7Berror`7D`7Berror code`7D`7BOUT`7D
X** \paramend
X** \blurb`7BThis function computes
X**  a matrix to perform the rotation (about the origin) of the 3D vector
X** \pardesc`7Bv1`7D to the 3D vector`20
X** \pardesc`7Bv2`7D, and concatenates this matrix with`20
X** \pardesc`7Bmatrix`7D on the basis of \pardesc`7Boperation`7D (\cite`7Brog
V:mecg`7D,`20
X**  pages 55--59). If the parameters are invalid, \pardesc`7Berror`7D is set
V to`20
X** -1. Otherwise, its value is \pardesc`7Bptkcpcok`7D.`7D
X*/
X`7B
X  Pmatrix temp;
X  Ppoint dv1, dv2, nvec;
X  Pfloat t1, costheta, sinetheta;
X
X  /* Create a rotation matrix (about origin) to rotate from
X  ** vector dv1 to vector dv2
X  ** Cf ROGERS and ADAMS page 55-59
X  */
X
X  *error = ptkcpcok;   /* Assume OK */
X  dv1 = ptk_unitv(v1);
X  dv2 = ptk_unitv(v2);
X
X  /* Get sine theta (cross product) and cos theta (dot product) */
X
X  costheta = dv1.x * dv2.x + dv1.y * dv2.y;
X  nvec = ptk_point(0.0, 0.0);
X
X  /* How do I choose sign of sinetheta?
X  ** let A, B and K be unit vectors, then
X  ** A X B :=  sin(theta) K
X  ** The sign of sinetheta is not important
X  ** because of the way it appears in the subsequent formulae
X  */
X
X  sinetheta = ptk_modv(&nvec);
X  if (ptk_equal(fabs(sinetheta), 0.0))`20
X  `7B
X    *error = -1;
X    sinetheta = ptkcpceps;
X  `7D
X  nvec = ptk_unitv(&nvec);
X  ptk_unitmatrix(temp);   /* get a unit matrix */
X  t1 = 1.0 - costheta;
X  temp`5B0`5D`5B0`5D = (nvec.x * nvec.x) + (1.0 - (nvec.x * nvec.x)) * costh
Veta;
X  temp`5B1`5D`5B1`5D = (nvec.y * nvec.y) + (1.0 - (nvec.y * nvec.y)) * costh
Veta;
X  temp`5B2`5D`5B2`5D = costheta;
X  temp`5B1`5D`5B0`5D = (nvec.x * nvec.y * t1);
X  temp`5B2`5D`5B0`5D = (nvec.y * sinetheta);`20
X  temp`5B0`5D`5B1`5D = (nvec.x * nvec.y * t1);
X  temp`5B2`5D`5B1`5D = (nvec.x * sinetheta);
X  temp`5B0`5D`5B2`5D = (nvec.y * sinetheta);
X  temp`5B1`5D`5B2`5D = (nvec.x * sinetheta);
X  ptk_concatenatematrix(operation, temp, matrix, matrix);
X`7D  /* ptk_rotatevv */
X
X/*--------------------------------------------------------------------------
V*/
X
X/*function:external*/
Xextern void ptk_rotateline3(C(Ppoint3 *) p1, C(Ppoint3 *) p2,`20
X                            C(Pfloat) theta, C(Pcomptype) operation,`20
X                            C(Pmatrix3) matrix, C(Pint *)error)
XPreANSI(Ppoint3 *p1)
XPreANSI(Ppoint3 *p2)
XPreANSI(Pfloat theta)`20
XPreANSI(Pcomptype operation)
XPreANSI(Pmatrix3 matrix)
XPreANSI(Pint *error)
X/*
X** \parambegin
X** \param`7BPpoint3 *`7D`7Bp1`7D`7B3D point`7D`7BIN`7D
X** \param`7BPpoint3 *`7D`7Bp2`7D`7B3D point`7D`7BIN`7D
X** \param`7BPfloat`7D`7Btheta`7D`7Bangle in degrees`7D`7BIN`7D
X** \param`7BPcomptype`7D`7Boperation`7D`7Bconcatenation operation`7D`7BIN`7D
X** \param`7BPmatrix3`7D`7Bmatrix`7D`7B4x4 matrix`7D`7BOUT`7D
X** \param`7BPint *`7D`7Berror`7D`7Berror code`7D`7BOUT`7D
X** \paramend
X** \blurb`7BThis function computes
X** a matrix to perform a 3D rotation of
X** \pardesc`7Btheta`7D degrees
X** about the line connecting \pardesc`7Bp1`7D to \pardesc`7Bp2`7D,
X** and concatenates this matrix with`20
X** \pardesc`7Bmatrix`7D on the basis of \pardesc`7Boperation`7D (\cite`7Brog
V:mecg`7D,`20
X**  pages 55--59).
X** If the parameters are invalid, \pardesc`7Berror`7D is set to`20
X** -1. Otherwise, its value is \pardesc`7Bptkcpcok`7D.`7D
X*/
X`7B
X  Pmatrix3 temp;
X  Ppoint3 nvec, mnvec;
X  Pfloat t1, costheta, sinetheta;
X
X  *error = ptkcpcok;
X  /* make a unit vector, i.e. direction cosines
X  ** nvec = p2 - p1.
X  */
X  nvec = ptk_subv3(p2, p1);
X  if (ptk_nullv3(&nvec)) `7B*error = -1;`7D
X  nvec = ptk_unitv3(&nvec);
X  costheta = cos(theta * ptkcdtor);
X  sinetheta = sin(theta * ptkcdtor);
X  t1 = 1.0 - costheta;
X  ptk_unitmatrix3(temp); /* get a unit matrix */
X  temp`5B0`5D`5B0`5D = (nvec.x * nvec.x) + (1.0 - (nvec.x * nvec.x)) * costh
Veta;
X  temp`5B1`5D`5B1`5D = (nvec.y * nvec.y) + (1.0 - (nvec.y * nvec.y)) * costh
Veta;
X  temp`5B2`5D`5B2`5D = (nvec.z * nvec.z) + (1.0 - (nvec.z * nvec.z)) * costh
Veta;
X  temp`5B1`5D`5B0`5D = (nvec.x * nvec.y * t1) + (nvec.z * sinetheta);
X  temp`5B2`5D`5B0`5D = (nvec.x * nvec.z * t1) - (nvec.y * sinetheta);`20
X  temp`5B0`5D`5B1`5D = (nvec.x * nvec.y * t1) - (nvec.z * sinetheta);
X  temp`5B2`5D`5B1`5D = (nvec.y * nvec.z * t1) + (nvec.x * sinetheta);
X  temp`5B0`5D`5B2`5D = (nvec.x * nvec.z * t1) + (nvec.y * sinetheta);
X  temp`5B1`5D`5B2`5D = (nvec.y * nvec.z * t1) - (nvec.x * sinetheta);
X  mnvec = ptk_scalev3(p1, -1.0);
X  ptk_shift3(&mnvec, PPRECONCATENATE, temp);
X  ptk_shift3(p1, PPOSTCONCATENATE, temp);
X  ptk_concatenatematrix3(operation, temp, matrix, matrix);
X`7D  /* ptk_rotateline3 */
X
X/*--------------------------------------------------------------------------
V*/
X
X/*function:external*/
Xextern void ptk_rotateline(C(Ppoint *) p1, C(Ppoint *) p2,`20
X                            C(Pfloat) theta, C(Pcomptype) operation,`20
X                            C(Pmatrix) matrix, C(Pint *)error)
XPreANSI(Ppoint *p1)
XPreANSI(Ppoint *p2)
XPreANSI(Pfloat theta)`20
XPreANSI(Pcomptype operation)
XPreANSI(Pmatrix matrix)
XPreANSI(Pint *error)
X/*
X** \parambegin
X** \param`7BPpoint *`7D`7Bp1`7D`7B2D point`7D`7BIN`7D
X** \param`7BPpoint *`7D`7Bp2`7D`7B2D point`7D`7BIN`7D
X** \param`7BPfloat`7D`7Btheta`7D`7Bangle in degrees`7D`7BIN`7D
X** \param`7BPcomptype`7D`7Boperation`7D`7Bconcatenation operation`7D`7BIN`7D
X** \param`7BPmatrix`7D`7Bmatrix`7D`7B3x3 matrix`7D`7BOUT`7D
X** \param`7BPint *`7D`7Berror`7D`7Berror code`7D`7BOUT`7D
X** \paramend
X** \blurb`7BThis function computes
X** a matrix to perform a 2D rotation of
X** \pardesc`7Btheta`7D degrees
X** about the line connecting \pardesc`7Bp1`7D to \pardesc`7Bp2`7D,
X** and concatenates this matrix with`20
X** \pardesc`7Bmatrix`7D on the basis of \pardesc`7Boperation`7D.
X** If the parameters are invalid, \pardesc`7Berror`7D is set to`20
X** -1. Otherwise, its value is \pardesc`7Bptkcpcok`7D.`7D
X*/
X`7B
X  Pmatrix temp;
X  Ppoint nvec, mnvec;
X  Pfloat t1, costheta, sinetheta;
X
X  *error = ptkcpcok;
X  /* make a unit vector, i.e. direction cosines
X  ** nvec = p2 - p1.
X  */
X  nvec = ptk_subv(p2, p1);
X  if (ptk_nullv(&nvec)) `7B*error = -1;`7D
X  nvec = ptk_unitv(&nvec);
X  costheta = cos(theta * ptkcdtor);
X  sinetheta = sin(theta * ptkcdtor);
X  t1 = 1.0 - costheta;
X  ptk_unitmatrix(temp); /* get a unit matrix */
X  temp`5B0`5D`5B0`5D = (nvec.x * nvec.x) + (1.0 - (nvec.x * nvec.x)) * costh
Veta;
X  temp`5B1`5D`5B1`5D = (nvec.y * nvec.y) + (1.0 - (nvec.y * nvec.y)) * costh
Veta;
X  temp`5B2`5D`5B2`5D = costheta;
X  temp`5B1`5D`5B0`5D = (nvec.x * nvec.y * t1);
X  temp`5B2`5D`5B0`5D = (nvec.y * sinetheta);`20
X  temp`5B0`5D`5B1`5D = (nvec.x * nvec.y * t1);
X  temp`5B2`5D`5B1`5D = (nvec.x * sinetheta);
X  temp`5B0`5D`5B2`5D = (nvec.y * sinetheta);
X  temp`5B1`5D`5B2`5D = (nvec.x * sinetheta);
X  mnvec = ptk_scalev(p1, -1.0);
X  ptk_shift(&mnvec, PPRECONCATENATE, temp);
X  ptk_shift(p1, PPOSTCONCATENATE, temp);
X  ptk_concatenatematrix(operation, temp, matrix, matrix);
X`7D  /* ptk_rotateline */
X
X/*--------------------------------------------------------------------------
V*/
X
X/*function:external*/
Xextern Ppoint4 ptk_pt3topt4(C(Ppoint3 *) point)
XPreANSI(Ppoint3 *point)
X/*
X** \parambegin
X** \param`7BPpoint3 *`7D`7Bpoint`7D`7B3D point`7D`7BIN`7D
X** \paramend
X** \blurb`7BThis function converts the 3D point \pardesc`7Bpoint`7D to a 4D
V point,
X** which is returned as the result of the function. The $w$ coordinate of th
Ve`20
X** 4D point is set to $1.0$.`7D
X*/
X`7B
X  Ppoint4 temp;
X
X  temp.x = point->x;
X  temp.y = point->y;
X  temp.z = point->z;
X  temp.w = 1.0;
X  return temp;
X`7D  /* ptk_pt3topt4 */
X
X/*--------------------------------------------------------------------------
V*/
X
X/*function:external*/
Xextern Ppoint3 ptk_pt4topt3(C(Ppoint4 *) point)
XPreANSI(Ppoint4 *point)
X/*
X** \parambegin
X** \param`7BPpoint4 *`7D`7Bpoint`7D`7B4D point`7D`7BIN`7D
X** \paramend
X** \blurb`7BThis function converts the 4D point \pardesc`7Bpoint`7D to a 3D
V point, by
X** dividing by $w$. The 3D point is returned as the result of the function.`
V7D
X*/
X`7B
X  Ppoint3 temp;
X  Pfloat w;
X
X  if (ptk_equal(fabs(point->w), 0.0))
X    w = 1.0 / ptkcpceps;
X  else
X    w = 1.0 / point->w;
X  temp.x = point->x * w;
X  temp.y = point->y * w;
X  temp.z = point->z * w;
X  return temp;
X`7D  /* ptk_pt4topt3 */
X
X/*--------------------------------------------------------------------------
V*/
X
X/*function:external*/
Xextern Ppoint4 ptk_transform4(C(Pmatrix3) matrix, C(Ppoint4 *) point)
XPreANSI(Pmatrix3 matrix)
XPreANSI(Ppoint4 *point)
X/*
X** \parambegin
X** \param`7BPmatrix3`7D`7Bmatrix`7D`7B4x4 matrix`7D`7BIN`7D
X** \param`7BPpoint4 *`7D`7Bpoint`7D`7B4D point`7D`7BIN`7D
X** \paramend
X** \blurb`7BThis function performs the 4D transformation`20
X** (that is, with no homogeneous division) of the
X** point \pardesc`7Bpt`7D by the $4 \times 4$ matrix
X** \pardesc`7Bmatrix`7D, and returns the result as
X** the value of the function.`7D
X*/
X`7B
X  Ppoint4 temp;
X
X  temp.x = matrix`5B0`5D`5B0`5D * point->x + matrix`5B0`5D`5B1`5D * point->y
V +
X`09     matrix`5B0`5D`5B2`5D * point->z + matrix`5B0`5D`5B3`5D * point->w;
X  temp.y = matrix`5B1`5D`5B0`5D * point->x + matrix`5B1`5D`5B1`5D * point->y
V +
X`09     matrix`5B1`5D`5B2`5D * point->z + matrix`5B1`5D`5B3`5D * point->w;
X  temp.z = matrix`5B2`5D`5B0`5D * point->x + matrix`5B2`5D`5B1`5D * point->y
V +
X`09     matrix`5B2`5D`5B2`5D * point->z + matrix`5B2`5D`5B3`5D * point->w;
X  temp.w = matrix`5B3`5D`5B0`5D * point->x + matrix`5B3`5D`5B1`5D * point->y
V +
X`09     matrix`5B3`5D`5B2`5D * point->z + matrix`5B3`5D`5B3`5D * point->w;
X  return temp;
X`7D  /* ptk_transform4 */
X
X/*--------------------------------------------------------------------------
V*/
X
X/*function:external*/
Xextern Ppoint3 ptk_transform3(C(Pmatrix3) matrix, C(Ppoint3 *) point)
XPreANSI(Pmatrix3 matrix)
XPreANSI(Ppoint3 *point)
X/*
X** \parambegin
X** \param`7BPmatrix3`7D`7Bmatrix`7D`7B4x4 matrix`7D`7BIN`7D
X** \param`7BPpoint3 *`7D`7Bpoint`7D`7B3D point`7D`7BIN`7D
X** \paramend
X** \blurb`7BThis function transforms the 3D point \pardesc`7Bpoint`7D by
X** the $4 \times 4$ matrix \pardesc`7Bmatrix`7D,
X** including homogeneous division.`7D
X*/
X`7B
X  Ppoint3 temp;
X  Pfloat w;
X
X  temp.x = matrix`5B0`5D`5B0`5D * point->x + matrix`5B0`5D`5B1`5D * point->y
V +
X`09     matrix`5B0`5D`5B2`5D * point->z + matrix`5B0`5D`5B3`5D;
X  temp.y = matrix`5B1`5D`5B0`5D * point->x + matrix`5B1`5D`5B1`5D * point->y
V +
X`09     matrix`5B1`5D`5B2`5D * point->z + matrix`5B1`5D`5B3`5D;
X  temp.z = matrix`5B2`5D`5B0`5D * point->x + matrix`5B2`5D`5B1`5D * point->y
V +
X`09     matrix`5B2`5D`5B2`5D * point->z + matrix`5B2`5D`5B3`5D;
X  w = matrix`5B3`5D`5B0`5D * point->x + matrix`5B3`5D`5B1`5D * point->y +`20
X        matrix`5B3`5D`5B2`5D * point->z + matrix`5B3`5D`5B3`5D;
X  if (ptk_equal(fabs(w), 0.0))
X    w = 1.0 / ptkcpceps;
X  else
X    w = 1.0 / w;
X  temp.x *= w;
X  temp.y *= w;
X  temp.z *= w;
X  return temp;
X`7D  /* ptk_transform3 */
X
X/*-------------------------------------------------------------------*/
X
X/*function:external*/
Xextern Ppoint ptk_transform(C(Pmatrix) matrix, C(Ppoint *) point)
XPreANSI(Pmatrix matrix)
XPreANSI(Ppoint *point)
X/*
X** \parambegin
X** \param`7BPmatrix`7D`7Bmatrix`7D`7B3x3 matrix`7D`7BIN`7D
X** \param`7BPpoint *`7D`7Bpoint`7D`7B2D point`7D`7BIN`7D
X** \paramend
X** \blurb`7BThis function transforms the 2D point \pardesc`7Bpoint`7D by the
X** $3 \times 3$ matrix \pardesc`7Bmatrix`7D.`7D
X*/
X`7B
X  Ppoint temp;
X  Pfloat w;
X
X  temp.x = matrix`5B0`5D`5B0`5D * point->x + matrix`5B0`5D`5B1`5D * point->y
V + matrix`5B0`5D`5B2`5D;
X  temp.y = matrix`5B1`5D`5B0`5D * point->x + matrix`5B1`5D`5B1`5D * point->y
V + matrix`5B1`5D`5B2`5D;
X  return temp;
X`7D  /* ptk_transform */
X
X/*-------------------------------------------------------------------*/
X
X/*function:external*/
Xextern void ptk_matrixtomatrix3(C(Pmatrix) mat, C(Pmatrix3) mat3)
XPreANSI(Pmatrix mat)
XPreANSI(Pmatrix mat3)
X/*
X** \parambegin
X** \param`7BPmatrix`7D`7Bmat`7D`7B3x3 matrix`7D`7BIN`7D
X** \param`7BPmatrix3`7D`7Bmat3`7D`7B4x4 matrix`7D`7BOUT`7D
X** \paramend
X** \blurb`7BThis function converts the $3 \times 3$ matrix \pardesc`7Bmat`7D
X** to the $4 \times 4$ matrix \pardesc`7Bmat3`7D, as follows:
X** $$  \left( \begin`7Barray`7D`7Bccc`7D
X           a & b & c\\
X           d & e & f\\
X           g & h & j
X         \end`7Barray`7D \right)
X
+-+-+-+-+-+-+-+-  END  OF PART 264 +-+-+-+-+-+-+-+-