COSMOS core  1.0.2 (beta)
Comprehensive Open-architecture Solution for Mission Operations Systems
Collaboration diagram for Math library functions:

Classes

class  LsFit
 

Functions

double gaussian_random (double mean, double stdev)
 Normal Distribution random number. More...
 
quaternion q_drotate_between_rv (rvector from, rvector to)
 Create rotation quaternion from 2 row vectors. More...
 
quaternion q_change_around_rv (rvector around, double angle)
 Create rotation quaternion from row vector axis and angle. More...
 
quaternion q_irotate_for (rvector sourcea, rvector sourceb, rvector targeta, rvector targetb)
 Create irotate quaternion from two orthogonal vectors. More...
 
rvector rv_quaternion2axis (quaternion q)
 Quaternion to row vector axis and angle. More...
 
rmatrix rm_quaternion2dcm (quaternion q)
 Quaternion to row matrix Direction Cosine Matrix. More...
 
quaternion q_axis2quaternion_rv (rvector v)
 Row vector axis and angle to Quaternion. More...
 
double distance_rv (rvector p0, rvector p1, rvector p2)
 Distance from a line. More...
 
double area_rv (rvector p0, rvector p1, rvector p2)
 Area of a triangle. More...
 
rvector rv_normalto (rvector p0, rvector p1, rvector p2)
 Normal to a polygon. More...
 
ByteOrder local_byte_order ()
 Determine local byte order. More...
 
uint16_t uint16from (uint8_t *pointer, ByteOrder order)
 Memory to 16 bit unsigned integer. More...
 
int16_t int16from (uint8_t *pointer, ByteOrder order)
 Memory to 16 bit signed integer. More...
 
uint32_t uint32from (uint8_t *pointer, ByteOrder order)
 Memory to 32 bit unsigned integer. More...
 
int32_t int32from (uint8_t *pointer, ByteOrder order)
 Memory to 32 bit signed integer. More...
 
float floatfrom (uint8_t *pointer, ByteOrder order)
 Memory to 32 bit float. More...
 
double doublefrom (uint8_t *pointer, ByteOrder order)
 Memory to 64 bit float. More...
 
void uint32to (uint32_t value, uint8_t *pointer, ByteOrder order)
 32 bit unsigned integer to memory More...
 
void int32to (int32_t value, uint8_t *pointer, ByteOrder order)
 32 bit signed integer to memory More...
 
void uint16to (uint16_t value, uint8_t *pointer, ByteOrder order)
 16 bit unsigned integer to memory More...
 
void int16to (int16_t value, uint8_t *pointer, ByteOrder order)
 16 bit signed integer to memory More...
 
void floatto (float value, uint8_t *pointer, ByteOrder order)
 32 bit floating point to memory More...
 
void doubleto (double value, uint8_t *pointer, ByteOrder order)
 64 bit floating point to memory More...
 
void open_estimate (estimatorhandle *estimate, uint32_t size, uint32_t degree)
 Initialize estimator. More...
 
int16_t set_estimate (estimatorhandle *estimate, double independent, double dependent)
 Set estimator. More...
 
estimatorstruc get_estimate (estimatorhandle *estimate, double x)
 Get estimate. More...
 
void multisolve (vector< vector< double > > x, vector< double > y, vector< double > &a)
 Perform N equation solution. More...
 
double evaluate_poly (double x, vector< double > parms)
 Evaluate polynomial. More...
 
double evaluate_poly_slope (double x, vector< double > parms)
 Evaluate polynomial slope. More...
 
double evaluate_poly_accel (double x, vector< double > parms)
 Evaluate polynomial acceleration. More...
 
double evaluate_poly_jerk (double x, vector< double > parms)
 Evaluate polynomial jerk. More...
 
rvector rv_evaluate_poly (double x, vector< vector< double > > parms)
 Evaluate vector polynomial. More...
 
rvector rv_evaluate_poly_slope (double x, vector< vector< double > > parms)
 Evaluate vector polynomial slope. More...
 
rvector rv_evaluate_poly_accel (double x, vector< vector< double > > parms)
 Evaluate vector polynomial acceleration. More...
 
rvector rv_evaluate_poly_jerk (double x, vector< vector< double > > parms)
 Evaluate vector polynomial jerk. More...
 
gvector gv_evaluate_poly (double x, vector< vector< double > > parms)
 Evaluate vector polynomial. More...
 
gvector gv_evaluate_poly_slope (double x, vector< vector< double > > parms)
 Evaluate vector polynomial slope. More...
 
gvector gv_evaluate_poly_accel (double x, vector< vector< double > > parms)
 Evaluate vector polynomial acceleration. More...
 
gvector gv_evaluate_poly_jerk (double x, vector< vector< double > > parms)
 Evaluate vector polynomial jerk. More...
 
quaternion q_evaluate_poly (double x, vector< vector< double > > parms)
 Evaluate quaternion polynomial. More...
 
quaternion q_evaluate_poly_slope (double x, vector< vector< double > > parms)
 Evaluate quaternion polynomial slope. More...
 
quaternion q_evaluate_poly_accel (double x, vector< vector< double > > parms)
 Evaluate quaternion polynomial acceleration. More...
 
quaternion q_evaluate_poly_jerk (double x, vector< vector< double > > parms)
 Evaluate quaternion polynomial jerk. More...
 
vector< double > polyfit (vector< double > &x, vector< double > &y)
 Perform general order polynomial fit. More...
 
uvector rv_fitpoly (uvector x, uvector y, uint32_t order)
 Perform nth order polynomial fit. More...
 
gj_kernelgauss_jackson_kernel (int32_t order, double dvi)
 Create Gauss-Jackson Integration Kernel. More...
 
void gauss_jackson_dekernel (gj_kernel *gjk)
 Free Gauss-Jackson Integration Kernel. More...
 
gj_instancegauss_jackson_instance (gj_kernel *kern, int32_t axes, void(*calc_vd2)(double vi, double *vd0, double *vd2, int32_t axes))
 Initialize an Instance of a Gauss-Jackson Integrator. More...
 
gj_stepgauss_jackson_step (gj_kernel *kern)
 Initialize a Step of a Gauss-Jackson integrator. More...
 
void gauss_jackson_destep (gj_kernel *kern, gj_step *step)
 Destroy a Step of a Gauss-Jackson integrator. More...
 
int gauss_jackson_setstep (gj_instance *gji, double vi, double *vd0, double *vd1, double *vd2, int32_t istep)
 Set Independent and Dependent variables for Gauss-Jackson step. More...
 
int gauss_jackson_getstep (gj_instance *gji, double *vi, double *vd0, double *vd1, double *vd2, int32_t istep)
 Get Independent and Dependent variables for Gauss-Jackson step. More...
 
void gauss_jackson_preset (gj_instance *gji)
 Converge all axes of a Gauss-Jackson integrator prior to propagation. More...
 
void gauss_jackson_extrapolate (gj_instance *gji, double target)
 Propagate Gauss-Jackson integration. More...
 
double fixangle (double angle)
 Limit angle to range 0-2PI. More...
 
double actan (double y, double x)
 ArcTan, limited to range 0-2PI. More...
 
double fixprecision (double number, double prec)
 Limit precision. More...
 
uint16_t calc_crc16ccitt (uint8_t *buf, int size, bool lsb)
 Calculate CRC-16-CCITT. More...
 
rvector drotate (quaternion q, rvector v)
 Rotate a row vector using a quaternion. More...
 
rvector rotate_q (quaternion q, rvector v)
 
cvector drotate (quaternion q, cvector v)
 Rotate a cartesian vector using a quaternion. More...
 
cvector rotate_q (quaternion q, cvector v)
 
rvector irotate (quaternion q, rvector v)
 Indirectly rotate a row vector using a quaternion. More...
 
rvector transform_q (quaternion q, rvector v)
 
cvector irotate (quaternion q, cvector v)
 Indirectly rotate a cartesian vector using a quaternion. More...
 
quaternion q_change_between_cv (cvector from, cvector to)
 Create rotation quaternion from 2 vectors. More...
 
quaternion q_change_between_rv (rvector from, rvector to)
 
cmatrix cm_change_between_cv (cvector from, cvector to)
 Create rotation matrix from 2 vectors. More...
 
double evaluate_poly (double x, rvector parms)
 
double evaluate_poly_slope (double x, rvector parms)
 
double evaluate_poly_accel (double x, rvector parms)
 
double evaluate_poly_jerk (double x, rvector parms)
 
uint16_t calc_crc16ccitt_lsb (vector< uint8_t >buf)
 
uint16_t calc_crc16ccitt_msb (vector< uint8_t >buf)
 
rvector rv_mmult (rmatrix m, rvector v)
 Multiply rmatrix by rvector. More...
 
rvector operator* (rmatrix m, rvector v)
 
rvector rv_diag (rmatrix a)
 Matrix diagonal. More...
 
cvector cv_mmult (cmatrix m, cvector v)
 Multiply cartesian vector by cartesian matrix. More...
 
rmatrix rm_diag (rvector a)
 Diagonal rmatrix Creates an rmatrix whose diagonal is filled with the supplied rvector. More...
 
rmatrix rm_eye ()
 Identity rmatrix. More...
 
rmatrix rm_zero ()
 Zero filled rmatrix. More...
 
double norm_rm (rmatrix mat)
 rmatrix norm. Calculates the Norm of the supplied rmatrix More...
 
double trace_rm (rmatrix mat)
 rmatrix Trace Calculates the trace of the supplied rmatrix. More...
 
rmatrix rm_transpose (rmatrix a)
 rmatrix Transpose. Calculate the transpose of the supplied rmatrix. More...
 
rmatrix rm_mmult (rmatrix a, rmatrix b)
 rmatrix Matrix Product More...
 
rmatrix rm_mult (rmatrix a, rmatrix b)
 Element-wise rmatrix multiplication. More...
 
rmatrix rm_smult (double a, rmatrix b)
 Scalar rmatrix multiplication. More...
 
rmatrix rm_add (rmatrix a, rmatrix b)
 rmatrix addition. Sum of two rmatrix values. More...
 
rmatrix rm_sub (rmatrix a, rmatrix b)
 rmatrix subtraction. Subtract two rmatrix values. More...
 
rmatrix rm_square (rmatrix a)
 Square rmatrix. More...
 
rmatrix rm_change_around_x (double angle)
 Rotation matrix for X axis. More...
 
rmatrix rm_change_around_y (double angle)
 Rotation matrix for Y axis. More...
 
rmatrix rm_change_around_z (double angle)
 Rotation matrix for Z axis. More...
 
rmatrix rm_change_around (int axis, double angle)
 Rotation matrix for indicated axis. More...
 
cmatrix cm_from_rm (rmatrix matrix)
 cmatrix from rmatrix More...
 
cmatrix cm_diag (cvector a)
 Diagonal cmatrix Creates an cmatrix whose diagonal is filled with the supplied cvector. More...
 
cmatrix cm_eye ()
 Identity cmatrix. More...
 
cmatrix cm_zero ()
 Zero filled cmatrix. More...
 
double norm_cm (cmatrix mat)
 cmatrix norm. Calculates the Norm of the supplied cmatrix More...
 
double trace_cm (cmatrix mat)
 cmatrix Trace Calculates the trace of the supplied cmatrix. More...
 
cmatrix cm_transpose (cmatrix a)
 cmatrix Transpose. Calculate the transpose of the supplied cmatrix. More...
 
cvector cv_diag (cmatrix a)
 Matrix diagonal. More...
 
cmatrix cm_mmult (cmatrix a, cmatrix b)
 cmatrix Matrix Product More...
 
cmatrix cm_mult (cmatrix a, cmatrix b)
 Element-wise cmatrix multiplication. More...
 
cmatrix cm_smult (double a, cmatrix b)
 Scalar cmatrix multiplication. More...
 
cmatrix cm_add (cmatrix a, cmatrix b)
 cmatrix addition. Sum of two cmatrix values. More...
 
cmatrix cm_sub (cmatrix a, cmatrix b)
 cmatrix subtraction. Subtract two cmatrix values. More...
 
cmatrix cm_square (cmatrix a)
 Square cmatrix. More...
 
cmatrix cm_change_around_x (double angle)
 Rotation matrix for X axis. More...
 
cmatrix cm_change_around_y (double angle)
 Rotation matrix for Y axis. More...
 
cmatrix cm_change_around_z (double angle)
 Rotation matrix for Z axis. More...
 
cmatrix cm_change_around (int axis, double angle)
 Rotation matrix for indicated axis. More...
 
rmatrix rm_from_cm (cmatrix matrix)
 rmatrix from cmatrix More...
 
rmatrix rm_from_rv (rvector vector, int direction)
 rvector to rmatrix. More...
 
rmatrix rm_skew (rvector row)
 Create skew symmetric rmatrix from rvector. More...
 
rvector rv_unskew (rmatrix matrix)
 Unskew 3x3 row matrix. More...
 
rmatrix rm_inverse (rmatrix m)
 Inverse of rmatrix. More...
 
rmatrix rm_from_m2 (matrix2d matrix)
 rmatrix from matrix2d More...
 
double determinant_rm (rmatrix m)
 Determinant of row column matrix. More...
 
matrix1d m1_zero (uint16_t cols)
 Fill 1D matrix with zeros. More...
 
matrix1d m1_smult (double number, matrix1d row)
 Multiply 1D matrix by a scalar. More...
 
matrix1d m1_add (matrix1d row1, matrix1d row2)
 Add one 1D matrix to another. More...
 
matrix1d m1_sub (matrix1d row1, matrix1d row2)
 Subtract one 1D matrix from another. More...
 
matrix1d m1_mmult (matrix2d Matrix, matrix1d Vector)
 Multiply matrix1d by matrix2d. More...
 
matrix1d m1_cross (matrix1d Vector1, matrix1d Vector2)
 matrix1d cross product More...
 
double m1_dot (matrix1d a, matrix1d b)
 matrix1d dot product More...
 
matrix2d m2_skew (matrix1d row)
 Create skew symmetric matrix2d from matrix1d. More...
 
matrix2d m2_diag (matrix1d row)
 Create diagonal matrix2d from matrix1d. More...
 
matrix2d m2_inverse (matrix2d m)
 Inverse of matrix2d. More...
 
double m2_determinant (matrix2d m)
 Determinant of a 2D matrix. More...
 
double m1_norm (matrix1d row)
 Compute the Euclidean norm of a 1D matrix. More...
 
matrix2d m2_zero (uint16_t rows, uint16_t cols)
 Create 2D zero matrix. More...
 
matrix2d m2_eye (uint16_t rows)
 Create 2D identity matrix. More...
 
matrix2d m2_smult (double number, matrix2d matrix)
 Multiply 2D matrix by a scalar. More...
 
matrix2d m2_add (matrix2d matrix1, matrix2d matrix2)
 Add one matrix2d to another. More...
 
matrix2d m2_sub (matrix2d matrix1, matrix2d matrix2)
 Subtract one matrix2d from another. More...
 
matrix2d m2_transpose (matrix2d matrix)
 Return transpose of a 2D matrix. More...
 
matrix1d m2_unskew (matrix2d matrix)
 Unskew 3x3 2D matrix. More...
 
double m2_trace (matrix2d matrix)
 Calculate the trace of a 2D matrix. More...
 
matrix2d m2_mmult (matrix2d matrix1, matrix2d matrix2)
 Matrix product. More...
 
matrix2d m2_from_rm (rmatrix matrix)
 rmatrix from rmatrix More...
 
matrix2d cm3x3_to_m2 (cmatrix matrix)
 
matrix1d cv_to_m1 (cvector vector)
 
matrix2d cv_to_m2 (cvector vector, int direction)
 
matrix2d m1_to_m2 (matrix1d vector, int direction)
 Matrix1d to matrix2d. More...
 
matrix1d m2_eig2x2 (matrix2d matrix)
 Eigen values of a 2x2 square matrix. More...
 
double m2_snorm2x2 (matrix2d matrix)
 Spectral norm of a 2x2 matrix. More...
 
cmatrix cm_quaternion2dcm (quaternion q)
 
quaternion q_dcm2quaternion_cm (cmatrix dcm)
 Direction Cosine Matrix to Quaternion. More...
 
rmatrix rm_change_between_rv (rvector from, rvector to)
 Create rotation matrix from 2 row vectors. More...
 
quaternion q_dcm2quaternion_rm (rmatrix m)
 Row matrix DCM to Quaternion. More...
 
double sep_rv (rvector v1, rvector v2)
 Angular separation between row vectors. More...
 
svector s_convert (rvector from)
 Convert rvector to svector. More...
 
rvector rv_convert (svector from)
 Convert svector to rvector. More...
 
rvector rv_zero ()
 Zero row order vector. More...
 
rvector rv_unitx (double scale)
 Scaled x row vector. More...
 
rvector rv_unity (double scale)
 Scaled y row vector. More...
 
rvector rv_unitz (double scale)
 Scaled z row vector. More...
 
rvector rv_one ()
 Row vector of ones. More...
 
rvector rv_one (double x, double y, double z)
 Row vector of values. More...
 
rvector rv_shortest (rvector v)
 Shortest vector. More...
 
rvector rv_shortest2 (rvector v)
 
rvector rv_normal (rvector v)
 Normalize row order vector. More...
 
void normalize_rv (rvector &v)
 Normalize row order vector in place. More...
 
double normVector3 (double x, double y, double z)
 basic function to compute the L2 norm of a 3d generic vector with separate entries More...
 
void normalizeVector3 (double &x, double &y, double &z)
 basic function to normalize any 3d vector with separate entries More...
 
rvector rv_smult (double a, rvector b)
 Multiply row vector by scalar. More...
 
rvector rv_sadd (double a, rvector b)
 Add scalar to each element of vector. More...
 
rvector rv_add (rvector a, rvector b)
 Add two row vectors. More...
 
rvector rv_sub (rvector a, rvector b)
 Subtract two vectors. More...
 
rvector rv_div (rvector a, rvector b)
 Divide two row vectors. More...
 
rvector rv_mult (rvector a, rvector b)
 Multiply two row vectors. More...
 
rvector rv_cross (rvector a, rvector b)
 Take cross product of two row vectors. More...
 
double dot_rv (rvector a, rvector b)
 Dot product of two row vectors. More...
 
double sep_cv (cvector v1, cvector v2)
 Angular separation between vectors. More...
 
cvector cv_zero ()
 Zero cartesian vector. More...
 
cvector cv_unitx ()
 Unit x vector. More...
 
cvector cv_unity ()
 Unit y vector. More...
 
cvector cv_unitz ()
 Unit z vector. More...
 
cvector cv_one ()
 Vector of ones. More...
 
cvector cv_normal (cvector v)
 Normalize cartesian vector. More...
 
void normalize_cv (cvector &v)
 Normalize cartesian vector in place, i.e. divides it by its own norm. More...
 
cvector cv_smult (double a, cvector b)
 Multiply vector by scalar. More...
 
cvector cv_sadd (double a, cvector b)
 Add scalar to each element of vector. More...
 
cvector cv_add (cvector a, cvector b)
 Add two vectors. More...
 
cvector cv_sub (cvector a, cvector b)
 Subtract two vectors. More...
 
cvector cv_div (cvector a, cvector b)
 Divide two vectors. More...
 
cvector cv_mult (cvector a, cvector b)
 Multiply two vectors. More...
 
cvector cv_cross (cvector a, cvector b)
 Take cross product of two vectors. More...
 
double dot_cv (cvector a, cvector b)
 
double length_cv (cvector v)
 
double cv_norm (cvector v)
 
double norm_cv (cvector v)
 
double sum_cv (cvector vec)
 
cvector cv_sqrt (cvector a)
 
bool equal_rv (rvector v1, rvector v2)
 Boolean equate of row vetor. More...
 
double length_rv (rvector v)
 Length of row vector. More...
 
double norm_rv (rvector vec)
 Infinite norm of row vector. More...
 
double sum_rv (rvector vec)
 Sum elements of a row vector. More...
 
rvector rv_sqrt (rvector vec)
 Row vector square root. More...
 
 LsFit::LsFit (uint16_t cnt=10, uint16_t ord=2)
 Multi element, variable order least squares fit. More...
 
void LsFit::initialize (uint16_t cnt=10, uint16_t ord=2)
 Initialize Least Squares Fit. More...
 
void LsFit::update (double x, double y)
 Update scalar Least Squares Fit. More...
 
void LsFit::update (double x, rvector y)
 Update rvector Least Squares Fit. More...
 
void LsFit::update (double x, gvector y)
 Update gvector Least Squares Fit. More...
 
void LsFit::update (double x, quaternion y)
 Update quaternion Least Squares Fit. More...
 
void LsFit::update (fitelement cfit, uint16_t dep)
 Update generic Least Squares Fit. More...
 
void LsFit::fit ()
 Calculate least squares fit. More...
 
double LsFit::lastx ()
 Least squares last independent value. More...
 
double LsFit::firstx ()
 Least squares first independent value. More...
 
size_t LsFit::size ()
 Least squares number of values. More...
 
double LsFit::eval (double x)
 Least squares dependent scalar value. More...
 
rvector LsFit::evalrvector (double x)
 Least squares dependent rvector value. More...
 
gvector LsFit::evalgvector (double x)
 Least squares dependent gvector value. More...
 
quaternion LsFit::evalquaternion (double x)
 Least squares dependent quaternion value. More...
 
double LsFit::slope (double x)
 Least squares dependent scalar 1st derivative. More...
 
rvector LsFit::slopervector (double x)
 Least squares dependent rvector 1st derivative. More...
 
gvector LsFit::slopegvector (double x)
 Least squares dependent gvector 1st derivative. More...
 
quaternion LsFit::slopequaternion (double x)
 Least squares dependent quaternion 1st derivative. More...
 
double LsFit::accel (double x)
 Least squares dependent scalar 2nd derivative. More...
 
rvector LsFit::accelrvector (double x)
 Least squares dependent rvector 2nd derivative. More...
 
gvector LsFit::accelgvector (double x)
 Least squares dependent gvector 2nd derivative. More...
 
quaternion LsFit::accelquaternion (double x)
 Least squares dependent quaternion 2nd derivative. More...
 
double LsFit::jerk (double x)
 Least squares dependent scalar 3rd derivative. More...
 
rvector LsFit::jerkrvector (double x)
 Least squares dependent rvector 3rd derivative. More...
 
gvector LsFit::jerkgvector (double x)
 Least squares dependent gvector 3rd derivative. More...
 
quaternion LsFit::jerkquaternion (double x)
 Least squares dependent quaternion 3rd derivative. More...
 
vector< vector< double > > LsFit::getparms ()
 Least Squares parameters. More...
 
double LsFit::getbasex ()
 Least Squares base. More...
 
double DCM::dotProduct (cvector a, cvector b)
 
cmatrix DCM::transposeMatrix (cmatrix a)
 
cmatrix DCM::base2_from_base1 (basisOrthonormal base2, basisOrthonormal base1)
 
cmatrix DCM::base1_from_base2 (basisOrthonormal base1, basisOrthonormal base2)
 
void basisOrthonormal::normalize ()
 
void cvector::normalize (double scale=1.)
 Normalize cartesian vector in place, i.e. divides it by its own norm. More...
 
cvector cvector::normalized (double scale=1.)
 Normalize cartesian vector. More...
 
double cvector::norm2 ()
 
double cvector::norm ()
 
double cvector::length ()
 
double & cvector::operator[] (const int index)
 Index into cvector. More...
 

Detailed Description

Function Documentation

double gaussian_random ( double  mean,
double  stdev 
)

Normal Distribution random number.

Random number generated using the Central Value Theorem to approximate a Gaussian distribution. Twelve random numbers between -1 and 1 are averaged to approximate the final random number. This will come from a distribution whose mean is 0 and variance is .5. The desired input parameters are then used to scale the output.

Parameters
meanDesired central value of the gaussian.
stdevDesired Standard Deviation of the gaussian.
Returns
Random number from desired distribution.
53 {
54  double trand;
55 
56  stdev *= stdev;
57 
58  trand = 0.;
59  for (uint32_t i=0; i<12; ++i)
60  {
61  trand += rand();
62  }
63 
64  trand = trand / (3. * RAND_MAX) - 2.;
65  trand = trand * stdev + mean;
66 
67  return (trand);
68 }
int i
Definition: rw_test.cpp:37
quaternion q_drotate_between_rv ( rvector  from,
rvector  to 
)

Create rotation quaternion from 2 row vectors.

Generate the quaternion that represents a direct rotation from one row order vector to a second row order vector.

Parameters
frominitial row order vector
tofinal row order vector
Returns
quaternion that can be used to rotate points
82 {
83  uvector rq = {{{0.,0.,0.},0.}};
84  rvector vec1, vec2;
85 
86  normalize_rv(from);
87  normalize_rv(to);
88 
89  if (length_rv(rv_add(from,to)) < 1e-14)
90  {
91  vec1.col[0] = rand();
92  vec1.col[1] = rand();
93  vec1.col[2] = rand();
94  normalize_rv(vec1);
95  vec2 = rv_cross(vec1,to);
96  normalize_rv(vec2);
97  if (length_rv(vec2)<D_SMALL)
98  {
99  vec1.col[0] = rand();
100  vec1.col[1] = rand();
101  vec1.col[2] = rand();
102  normalize_rv(vec1);
103  vec2 = rv_cross(vec1,to);
104  normalize_rv(vec2);
105  }
106  rq.r = vec2;
107  rq.q.w = 0.;
108  }
109  else
110  {
111  vec2 = rv_cross(from,to);
112  rq.r = vec2;
113  rq.q.w = 1. + dot_rv(from,to);
114  }
115 
116  normalize_q(&rq.q);
117  return (rq.q);
118 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
Definition: eci2kep_test.cpp:33
rvector rv_add(rvector a, rvector b)
Add two row vectors.
Definition: vector.cpp:299
rvector r
Definition: mathlib.h:161
3 element generic row vector
Definition: vector.h:53
double length_rv(rvector v)
Length of row vector.
Definition: vector.cpp:748
void normalize_rv(rvector &v)
Normalize row order vector in place.
Definition: vector.cpp:222
quaternion q
Definition: mathlib.h:157
double dot_rv(rvector a, rvector b)
Dot product of two row vectors.
Definition: vector.cpp:379
void normalize_q(quaternion *q)
Definition: vector.cpp:961
double w
Rotation.
Definition: vector.h:407
const double D_SMALL
Definition: math/constants.h:40
double col[3]
Definition: vector.h:55
rvector rv_cross(rvector a, rvector b)
Take cross product of two row vectors.
Definition: vector.cpp:363
quaternion q_change_around_rv ( rvector  around,
double  angle 
)

Create rotation quaternion from row vector axis and angle.

Generate the quaternion that represents a rotation of the specified angle around the specified row vector axis.

Parameters
aroundrow order vector around which the rotation will occur
angleamount of rotation in radians
Returns
quaternion that can be used to rotate points
129 {
130  double sa;
131  uvector rq{};
132 
133  angle /= 2.;
134  sa = sin(angle);
135  normalize_rv(around);
136 
137  rq.r = rv_smult(sa,around);
138  rq.q.w = cos(angle);
139  normalize_q(&rq.q);
140  return (rq.q);
141 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
void normalize_rv(rvector &v)
Normalize row order vector in place.
Definition: vector.cpp:222
rvector rv_smult(double a, rvector b)
Multiply row vector by scalar.
Definition: vector.cpp:266
void normalize_q(quaternion *q)
Definition: vector.cpp:961
quaternion q_irotate_for ( rvector  sourcea,
rvector  sourceb,
rvector  targeta,
rvector  targetb 
)

Create irotate quaternion from two orthogonal vectors.

Using two vectors, represented in both the original and target frames, calculate the quaternion that will irotate any vector from the original to the target frame.

Parameters
sourceaFirst vector in source frame
sourcebSecond vector in source frame
targetaFirst vector in target frame
targetbSecond vector in target frame
Returns
Quaternion to use with irotate to irotate from source to target.
155 {
156  quaternion qe_a;
157  quaternion qe_b;
158  quaternion fqe;
159 
160  // Determine rotation of source A into target A
161  qe_a = q_conjugate(q_drotate_between_rv(sourcea,targeta));
162 
163  // Use to irotate source B into intermediate B
164  sourceb = irotate(qe_a,sourceb);
165  normalize_rv(sourceb);
166  normalize_rv(targetb);
167  if (length_rv(rv_add(sourceb,targetb)) < 1e-14)
168  {
169  // Antiparallel - rotate 180 degrees around vector A
170  qe_b.d.x = -targeta.col[0];
171  qe_b.d.y = -targeta.col[1];
172  qe_b.d.z = -targeta.col[2];
173  qe_b.w = 0;
174  }
175  else
176  {
177  // Determine intrinsic rotation of this intermediate B into target B
178  qe_b = q_conjugate(q_drotate_between_rv(sourceb,targetb));
179  }
180 
181  // Combine to determine complete intrinsic rotation of source into target
182  fqe = q_fmult(qe_a,qe_b);
183  normalize_q(&fqe);
184 
185  return fqe;
186 }
double y
Y value.
Definition: vector.h:114
Definition: eci2kep_test.cpp:33
rvector rv_add(rvector a, rvector b)
Add two row vectors.
Definition: vector.cpp:299
cvector d
Orientation.
Definition: vector.h:405
double length_rv(rvector v)
Length of row vector.
Definition: vector.cpp:748
Quaternion, scalar last, using x, y, z.
Definition: vector.h:402
quaternion q_conjugate(quaternion q)
Definition: vector.cpp:1010
void normalize_rv(rvector &v)
Normalize row order vector in place.
Definition: vector.cpp:222
double x
X value.
Definition: vector.h:112
quaternion q_fmult(rvector r1, quaternion q2)
rvector quaternion multiply
Definition: vector.cpp:1079
quaternion q_drotate_between_rv(rvector from, rvector to)
Create rotation quaternion from 2 row vectors.
Definition: mathlib.cpp:81
void normalize_q(quaternion *q)
Definition: vector.cpp:961
double w
Rotation.
Definition: vector.h:407
double z
Z value.
Definition: vector.h:116
double col[3]
Definition: vector.h:55
rvector irotate(quaternion q, rvector v)
Indirectly rotate a row vector using a quaternion.
Definition: mathlib.cpp:2308
rvector rv_quaternion2axis ( quaternion  q)

Quaternion to row vector axis and angle.

Convert rotation quaternion to a directional vector in row order form the length of which represents the angle of rotation in radians.

Parameters
qQuaternion to be converted.
Returns
Row vector representing the quaternion.
195 {
196  double ca, sa;
197  uvector rq = {{{0.,0.,0.},0.}};
198 
199  normalize_q(&q);
200  ca = 2.*acos(q.w);
201  if (ca > 0. && ca < D2PI)
202  {
203  sa = sin(ca/2.);
204  rq.q = q;
205  rq.r = rv_smult(ca/sa,rq.r);
206  }
207  else
208  rq.r = rv_zero();
209 
210  return (rq.r);
211 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
rvector r
Definition: mathlib.h:161
quaternion q
Definition: mathlib.h:157
rvector rv_smult(double a, rvector b)
Multiply row vector by scalar.
Definition: vector.cpp:266
void normalize_q(quaternion *q)
Definition: vector.cpp:961
double w
Rotation.
Definition: vector.h:407
rvector rv_zero()
Zero row order vector.
Definition: vector.cpp:107
const double D2PI
Double precision 2*PI.
Definition: math/constants.h:16
rmatrix rm_quaternion2dcm ( quaternion  q)

Quaternion to row matrix Direction Cosine Matrix.

Convert rotation quaternion to an equivalent direction cosine matrix in row matrix form.

Parameters
qQuaternion representing rotation.
Returns
direction Cosine rotation matrix in row matrix form
220 {
221  rmatrix m;
222  double yy, xx, zz, xy, xz, xw, yz, yw, zw;
223 
224  normalize_q(&q);
225 
226  xx = 2. * q.d.x;
227  xy = xx * q.d.y;
228  xz = xx * q.d.z;
229  xw = xx * q.w;
230  xx *= q.d.x;
231  yy = 2. * q.d.y;
232  yz = yy * q.d.z;
233  yw = yy * q.w;
234  yy *= q.d.y;
235  zz = 2. * q.d.z;
236  zw = zz * q.w;
237  zz *= q.d.z;
238 
239  m.row[0].col[0] = 1. - yy - zz;
240  m.row[0].col[1] = xy - zw;
241  m.row[0].col[2] = xz + yw;
242  m.row[1].col[0] = xy + zw;
243  m.row[1].col[1] = 1. - xx - zz;
244  m.row[1].col[2] = yz - xw;
245  m.row[2].col[0] = xz - yw;
246  m.row[2].col[1] = yz + xw;
247  m.row[2].col[2] = 1. - xx - yy;
248 
249  return (m);
250 }
double y
Y value.
Definition: vector.h:114
cvector d
Orientation.
Definition: vector.h:405
double x
X value.
Definition: vector.h:112
void normalize_q(quaternion *q)
Definition: vector.cpp:961
3x3 element generic matrix
Definition: matrix.h:41
double w
Rotation.
Definition: vector.h:407
double z
Z value.
Definition: vector.h:116
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
quaternion q_axis2quaternion_rv ( rvector  v)

Row vector axis and angle to Quaternion.

Convert axis and angle orientation represented as row vector to a rotation Quaternion.

Parameters
vRow vector axis and angle.
Returns
Rotation Quaternion.
261 {
262  double length, s2;
263  quaternion q;
264 
265  length = sqrt(v.col[0]*v.col[0]+v.col[1]*v.col[1]+v.col[2]*v.col[2]);
266  s2 = sin(length/2.)/length;
267  if (length)
268  {
269  q.d.x = s2*v.col[0];
270  q.d.y = s2*v.col[1];
271  q.d.z = s2*v.col[2];
272  }
273  else
274  q.d.x = q.d.y = q.d.z = 0.;
275  q.w =cos(length/2);
276 
277  normalize_q(&q);
278  return (q);
279 }
double y
Y value.
Definition: vector.h:114
cvector d
Orientation.
Definition: vector.h:405
Quaternion, scalar last, using x, y, z.
Definition: vector.h:402
double x
X value.
Definition: vector.h:112
void normalize_q(quaternion *q)
Definition: vector.cpp:961
double w
Rotation.
Definition: vector.h:407
double z
Z value.
Definition: vector.h:116
png_uint_32 length
Definition: png.c:2173
double col[3]
Definition: vector.h:55
double distance_rv ( rvector  p0,
rvector  p1,
rvector  p2 
)

Distance from a line.

Calculates the distance of a third point from a line defined be two points.

Parameters
p0Point not on line.
p1First point defining line.
p2Second point defining line.
Returns
Distance of Third point from line. Zero if any problems.
292 {
293  rvector p21, p01, p02;
294  double d, l;
295 
296  p21 = rv_sub(p2,p1);
297  p01 = rv_sub(p0,p1);
298  p02 = rv_sub(p0,p2);
299 
300  if ((l=length_rv(p21)) == 0.)
301  return (0.);
302 
303  d = length_rv(rv_cross(p01,p02))/l;
304 
305  return(d);
306 }
3 element generic row vector
Definition: vector.h:53
double length_rv(rvector v)
Length of row vector.
Definition: vector.cpp:748
rvector rv_sub(rvector a, rvector b)
Subtract two vectors.
Definition: vector.cpp:315
rvector rv_cross(rvector a, rvector b)
Take cross product of two row vectors.
Definition: vector.cpp:363
double area_rv ( rvector  p0,
rvector  p1,
rvector  p2 
)

Area of a triangle.

Calculates the area of a triangle defined by three points.

Parameters
p0First point defining triangle.
p1Second point defining triangle.
p2Third point defining triangle.
Returns
Area of triangle. Zero if any problems.
337 {
338  double base, altitude;
339 
340  base = length_rv(rv_sub(p1,p2));
341  altitude = distance_rv(p0,p1,p2);
342 
343  return(.5*base*altitude);
344 }
double length_rv(rvector v)
Length of row vector.
Definition: vector.cpp:748
double distance_rv(rvector p0, rvector p1, rvector p2)
Distance from a line.
Definition: mathlib.cpp:291
rvector rv_sub(rvector a, rvector b)
Subtract two vectors.
Definition: vector.cpp:315
rvector rv_normalto ( rvector  p0,
rvector  p1,
rvector  p2 
)

Normal to a polygon.

Calculate the normal vector to a polygon that contains at least three sequential vertices.

Parameters
p0First point defining polygon.
p1Second point defining polygon.
p2Third point defining polygon.
Returns
rvector representing normal in a normalized form.
355 {
356  rvector normal, p10, p20;
357 
358  p10 = rv_sub(p1,p0);
359  p20 = rv_sub(p2,p0);
360 
361  normal = rv_cross(p10,p20);
362 
363  return(rv_normal(normal));
364 }
3 element generic row vector
Definition: vector.h:53
rvector rv_normal(rvector v)
Normalize row order vector.
Definition: vector.cpp:212
rvector rv_sub(rvector a, rvector b)
Subtract two vectors.
Definition: vector.cpp:315
rvector rv_cross(rvector a, rvector b)
Take cross product of two row vectors.
Definition: vector.cpp:363
ByteOrder local_byte_order ( )

Determine local byte order.

Investigate a locally stored number to determine the byte order of the local machine.

Returns
Order as provided in ByteOrder.
376 {
377  uint16_t test = 1;
378  uint8_t *check;
379 
380  check = (uint8_t *)&test;
381 
382  if (check[0] == 0)
383  return (ByteOrder::BIGENDIAN);
384  else
385  return (ByteOrder::LITTLEENDIAN);
386 }
void test()
Definition: testprofiler.c:34
Big Endian byte order.
Little Endian byte order.
uint16_t uint16from ( uint8_t *  pointer,
ByteOrder  order 
)

Memory to 16 bit unsigned integer.

Return the 16 bit unsigned integer equivalent of a location in memory, corrected for the local byte order.

Parameters
pointerlocation in memory to be cast
orderbyte order of the data in memory. Taken from ByteOrder.
Returns
16 bit unsigned integer
396 {
397  uint16_t *result;
398  uint8_t *rb;
399  double rd;
400 
401  rb = (uint8_t *)&rd;
402 
403  result = (uint16_t *)rb;
404  if (local_byte_order() == order)
405  {
406  memcpy((void *)rb,pointer,2);
407  }
408  else
409  {
410  rb[1] = pointer[0];
411  rb[0] = pointer[1];
412  }
413 
414  return (*result);
415 }
ByteOrder local_byte_order()
Determine local byte order.
Definition: mathlib.cpp:375
int16_t int16from ( uint8_t *  pointer,
ByteOrder  order 
)

Memory to 16 bit signed integer.

Return the 16 bit signed integer equivalent of a location in memory, corrected for the local byte order.

Parameters
pointerlocation in memory to be cast
orderbyte order of the data in memory. Taken from ByteOrder.
Returns
16 bit signed integer
425 {
426  int16_t *result;
427  uint16_t rb;
428 
429  result = (int16_t *)&rb;
430  rb = uint16from(pointer,order);
431 
432  return (*result);
433 }
uint16_t uint16from(uint8_t *pointer, ByteOrder order)
Memory to 16 bit unsigned integer.
Definition: mathlib.cpp:395
uint32_t uint32from ( uint8_t *  pointer,
ByteOrder  order 
)

Memory to 32 bit unsigned integer.

Return the 32 bit unsigned integer equivalent of a location in memory, corrected for the local byte order.

Parameters
pointerlocation in memory to be cast
orderbyte order of the data in memory. Taken from ByteOrder.
Returns
32 bit unsigned integer
443 {
444  uint32_t *result;
445  uint8_t *rb;
446  double rd;
447 
448  rb = (uint8_t *)&rd;
449 
450  result = (uint32_t *)rb;
451  if (local_byte_order() == order)
452  {
453  memcpy((void *)rb,pointer,4);
454  }
455  else
456  {
457  rb[3] = pointer[0];
458  rb[2] = pointer[1];
459  rb[1] = pointer[2];
460  rb[0] = pointer[3];
461  }
462 
463  return (*result);
464 }
ByteOrder local_byte_order()
Determine local byte order.
Definition: mathlib.cpp:375
int32_t int32from ( uint8_t *  pointer,
ByteOrder  order 
)

Memory to 32 bit signed integer.

Return the 32 bit signed integer equivalent of a location in memory, corrected for the local byte order.

Parameters
pointerlocation in memory to be cast
orderbyte order of the data in memory. Taken from ByteOrder.
Returns
32 bit signed integer
474 {
475  int32_t *result;
476  uint32_t rb;
477 
478  result = (int32_t *)&rb;
479  rb = uint32from(pointer,order);
480 
481  return (*result);
482 }
uint32_t uint32from(uint8_t *pointer, ByteOrder order)
Memory to 32 bit unsigned integer.
Definition: mathlib.cpp:442
float floatfrom ( uint8_t *  pointer,
ByteOrder  order 
)

Memory to 32 bit float.

Return the 32 bit float equivalent of a location in memory, corrected for the local byte order.

Parameters
pointerlocation in memory to be cast
orderbyte order of the data in memory. Taken from ByteOrder.
Returns
32 bit float
492 {
493  float result;
494  uint8_t *rb;
495 
496  rb = (uint8_t *)&result;
497  if (local_byte_order() == order)
498  {
499  memcpy((void *)rb,pointer,4);
500  }
501  else
502  {
503  rb[3] = pointer[0];
504  rb[2] = pointer[1];
505  rb[1] = pointer[2];
506  rb[0] = pointer[3];
507  }
508 
509  return (result);
510 }
ByteOrder local_byte_order()
Determine local byte order.
Definition: mathlib.cpp:375
double doublefrom ( uint8_t *  pointer,
ByteOrder  order 
)

Memory to 64 bit float.

Return the 64 bit float equivalent of a location in memory, corrected for the local byte order.

Parameters
pointerlocation in memory to be cast
orderbyte order of the data in memory. Taken from ByteOrder.
Returns
64 bit float
520 {
521  double result;
522  uint8_t *rb;
523 
524  rb = (uint8_t *)&result;
525  if (local_byte_order() == order)
526  {
527  memcpy((void *)rb,pointer,8);
528  }
529  else
530  {
531  rb[7] = pointer[0];
532  rb[6] = pointer[1];
533  rb[5] = pointer[2];
534  rb[4] = pointer[3];
535  rb[3] = pointer[4];
536  rb[2] = pointer[5];
537  rb[1] = pointer[6];
538  rb[0] = pointer[7];
539  }
540 
541  return (result);
542 }
ByteOrder local_byte_order()
Determine local byte order.
Definition: mathlib.cpp:375
void uint32to ( uint32_t  value,
uint8_t *  pointer,
ByteOrder  order 
)

32 bit unsigned integer to memory

Cast a 32 bit unsigned integer equivalent into a location in memory, corrected for the local byte order.

Parameters
valueinteger to be cast
pointerlocation in memory
orderdesired byte order of the data in memory. Taken from ByteOrder.
552 {
553  uint32_t *result;
554  uint8_t *rb;
555  double rd;
556 
557  rb = (uint8_t *)&rd;
558 
559  result = (uint32_t *)rb;
560  *result = value;
561  if (local_byte_order() == order)
562  {
563  memcpy(pointer,(void *)rb,4);
564  }
565  else
566  {
567  pointer[0] = rb[3];
568  pointer[1] = rb[2];
569  pointer[2] = rb[1];
570  pointer[3] = rb[0];
571  }
572 
573 }
ByteOrder local_byte_order()
Determine local byte order.
Definition: mathlib.cpp:375
void int32to ( int32_t  value,
uint8_t *  pointer,
ByteOrder  order 
)

32 bit signed integer to memory

Cast a 32 bit signed integer equivalent into a location in memory, corrected for the local byte order.

Parameters
valueinteger to be cast
pointerlocation in memory
orderdesired byte order of the data in memory. Taken from ByteOrder.
583 {
584  int32_t *result;
585  uint8_t *rb;
586  double rd;
587 
588  rb = (uint8_t *)&rd;
589 
590  result = (int32_t *)rb;
591  *result = value;
592  if (local_byte_order() == order)
593  {
594  memcpy(pointer,(void *)rb,4);
595  }
596  else
597  {
598  pointer[0] = rb[3];
599  pointer[1] = rb[2];
600  pointer[2] = rb[1];
601  pointer[3] = rb[0];
602  }
603 
604 }
ByteOrder local_byte_order()
Determine local byte order.
Definition: mathlib.cpp:375
void uint16to ( uint16_t  value,
uint8_t *  pointer,
ByteOrder  order 
)

16 bit unsigned integer to memory

Cast a 16 bit unsigned integer equivalent into a location in memory, corrected for the local byte order.

Parameters
valueinteger to be cast
pointerlocation in memory
orderdesired byte order of the data in memory. Taken from ByteOrder.
614 {
615  uint16_t *result;
616  uint8_t *rb;
617  double rd;
618 
619  rb = (uint8_t *)&rd;
620 
621  result = (uint16_t *)rb;
622  *result = value;
623  if (local_byte_order() == order)
624  {
625  memcpy(pointer,(void *)rb,2);
626  }
627  else
628  {
629  pointer[0] = rb[1];
630  pointer[1] = rb[0];
631  }
632 
633 }
ByteOrder local_byte_order()
Determine local byte order.
Definition: mathlib.cpp:375
void int16to ( int16_t  value,
uint8_t *  pointer,
ByteOrder  order 
)

16 bit signed integer to memory

Cast a 16 bit signed integer equivalent into a location in memory, corrected for the local byte order.

Parameters
valueinteger to be cast
pointerlocation in memory
orderdesired byte order of the data in memory. Taken from ByteOrder.
643 {
644  int16_t *result;
645  uint8_t *rb;
646  double rd;
647 
648  rb = (uint8_t *)&rd;
649 
650  result = (int16_t *)rb;
651  *result = value;
652  if (local_byte_order() == order)
653  {
654  memcpy(pointer,(void *)rb,2);
655  }
656  else
657  {
658  pointer[0] = rb[1];
659  pointer[1] = rb[0];
660  }
661 
662 }
ByteOrder local_byte_order()
Determine local byte order.
Definition: mathlib.cpp:375
void floatto ( float  value,
uint8_t *  pointer,
ByteOrder  order 
)

32 bit floating point to memory

Cast a 32 bit floating point equivalent into a location in memory, corrected for the local byte order.

Parameters
valuefloat to be cast
pointerlocation in memory
orderdesired byte order of the data in memory. Taken from ByteOrder.
672 {
673  float *result;
674  uint8_t *rb;
675  double rd;
676 
677  rb = (uint8_t *)&rd;
678 
679  result = (float *)rb;
680  *result = value;
681  if (local_byte_order() == order)
682  {
683  memcpy(pointer,(void *)rb,4);
684  }
685  else
686  {
687  pointer[0] = rb[3];
688  pointer[1] = rb[2];
689  pointer[2] = rb[1];
690  pointer[3] = rb[0];
691  }
692 
693 }
ByteOrder local_byte_order()
Determine local byte order.
Definition: mathlib.cpp:375
void doubleto ( double  value,
uint8_t *  pointer,
ByteOrder  order 
)

64 bit floating point to memory

Cast a 64 bit floating point equivalent into a location in memory, corrected for the local byte order.

Parameters
valuefloat to be cast
pointerlocation in memory
orderdesired byte order of the data in memory. Taken from ByteOrder.
703 {
704  double *result;
705  uint8_t *rb;
706  double rd;
707 
708  rb = (uint8_t *)&rd;
709 
710  result = (double *)rb;
711  *result = value;
712  if (local_byte_order() == order)
713  {
714  memcpy(pointer,(void *)rb,8);
715  }
716  else
717  {
718  pointer[0] = rb[7];
719  pointer[1] = rb[6];
720  pointer[2] = rb[5];
721  pointer[3] = rb[4];
722  pointer[4] = rb[3];
723  pointer[5] = rb[2];
724  pointer[6] = rb[1];
725  pointer[7] = rb[0];
726  }
727 }
ByteOrder local_byte_order()
Determine local byte order.
Definition: mathlib.cpp:375
void open_estimate ( estimatorhandle estimate,
uint32_t  size,
uint32_t  degree 
)

Initialize estimator.

Setup the provided estimatorhandle so that it can be fed values and provide estimates.

Parameters
estimatePointer to an estimatorhandle.
sizeThe number of estimates to be averaged for the total estimate.
degreeThe degree of the polynomial fit for the estimator.
737 {
738  if (degree > 4) degree = 4;
739  estimate->degree = degree;
740  if (size < degree) size = degree;
741  estimate->size = size;
742  estimate->r.resize(0);
743 }
vector< estimatorstruc > r
Definition: mathlib.h:272
uint32_t size
Definition: mathlib.h:274
uint32_t degree
Definition: mathlib.h:275
int16_t set_estimate ( estimatorhandle estimate,
double  independent,
double  dependent 
)

Set estimator.

Add a pair of independent and dependent values to the supplied estimatorhandle. If more than degree+1 values have been accumulated, a fit will be done for the new value. If more than size values have been added, the oldest will be dropped.

Parameters
estimatePointer to an estimatorhandle.
independentIndependent value.
dependentDependent value.
Returns
Number of stored values.
756 {
757 
758  if (estimate->r.size() == estimate->size && estimate->r[estimate->size-1].x.size() == estimate->degree+1)
759  {
760  for (uint32_t i=0; i<estimate->size-1; ++i)
761  {
762  estimate->r[i] = estimate->r[i+1];
763  }
764  for (uint32_t i=0; i<estimate->degree; ++i)
765  {
766  estimate->r[estimate->size-1].x[i] = estimate->r[estimate->size-1].x[i+1];
767  estimate->r[estimate->size-1].y[i] = estimate->r[estimate->size-1].y[i+1];
768  }
769  estimate->r[estimate->size-1].x[estimate->degree] = independent - estimate->xbase;
770  estimate->r[estimate->size-1].y[estimate->degree] = dependent - estimate->ybase;
771  estimate->r[estimate->size-1].a = polyfit(estimate->r[estimate->size-1].x,estimate->r[estimate->size-1].y);
772  }
773  else
774  {
775  if (estimate->r.size() > 0)
776  {
777  if (estimate->r.size() < estimate->size) estimate->r.resize(estimate->r.size()+1);
778  for (uint32_t i=0; i<estimate->r.size(); ++i)
779  {
780  if (estimate->r[i].x.size() < estimate->degree+1)
781  {
782  estimate->r[i].x.push_back(independent - estimate->xbase);
783  estimate->r[i].y.push_back(dependent - estimate->ybase);
784  }
785  if (estimate->r[i].x.size() == estimate->degree+1 && !estimate->r[i].a.size()) estimate->r[i].a = polyfit(estimate->r[i].x,estimate->r[i].y);
786  }
787 
788  }
789  else
790  {
791  estimate->xbase = independent;
792  estimate->ybase = dependent;
793  estimate->r.resize(1);
794  estimate->r[0].x.push_back(independent - estimate->xbase);
795  estimate->r[0].y.push_back(dependent - estimate->ybase);
796  }
797  }
798 
799  return (estimate->r.size()+estimate->r[estimate->r.size()-1].x.size());
800 }
vector< estimatorstruc > r
Definition: mathlib.h:272
double xbase
Definition: mathlib.h:276
double ybase
Definition: mathlib.h:277
int i
Definition: rw_test.cpp:37
vector< double > polyfit(vector< double > &x, vector< double > &y)
Perform general order polynomial fit.
Definition: mathlib.cpp:1434
uint32_t size
Definition: mathlib.h:274
uint32_t degree
Definition: mathlib.h:275
estimatorstruc get_estimate ( estimatorhandle estimate,
double  x 
)

Get estimate.

Return the best estimate for the supplied independent value. Estimate will be returned as an estimatorstruc, containing the value and likely error for the 0th, 1st and 2nd derivative of the dependent value.

Parameters
estimatePointer to an estimatorhandle.
xValue of independent variable for estimate.
Returns
estimatorstruc containing the estimate of the dependent value.
813 {
814  estimatorstruc result;
815  double tx;
816 
817  result.value[0] = result.value[1] = result.value[2] = 0.;
818  result.error[0] = result.error[1] = result.error[2] = 0.;
819  if (!estimate->r.size())
820  {
821  return (result);
822  }
823 
824  x -= estimate->xbase;
825 
826  if (estimate->r.size() == estimate->size && estimate->r[estimate->size-1].a.size() == estimate->degree+1)
827  {
828  for (uint32_t i=0; i<estimate->size; i++)
829  {
830  tx = 1.;
831  for (uint32_t j=0; j<estimate->degree+1; j++)
832  {
833  switch (j)
834  {
835  case 0:
836  estimate->r[i].value[0] = estimate->r[i].a[0];
837  if (estimate->degree > 0) estimate->r[i].value[1] = estimate->r[i].a[1];
838  if (estimate->degree > 1) estimate->r[i].value[2] = 2. * estimate->r[i].a[2];
839  break;
840  default:
841  tx *= x;
842  estimate->r[i].value[0] += tx * estimate->r[i].a[j];
843  if (estimate->degree > j) estimate->r[i].value[1] += tx * (j+1) * estimate->r[i].a[j+1];
844  if (estimate->degree > j+1) estimate->r[i].value[2] += tx * (j+2) * (j+1) * estimate->r[i].a[j+2];
845  break;
846  }
847  }
848 
849  for (int j=0; j<3; ++j)
850  {
851  result.value[j] += estimate->r[i].value[j];
852  result.error[j] += (estimate->r[i].value[j]*estimate->r[i].value[j]);
853  }
854 
855  }
856 
857  for (int i=0; i<3; i++)
858  {
859  result.error[i] = result.error[i] - result.value[i]*result.value[i]/(estimate->size);
860  if (result.error[i] < 0.) result.error[i] = 0.;
861  if (estimate->size > 1)
862  {
863  result.error[i] /= (estimate->size - 1);
864  }
865  result.error[i] = sqrt(result.error[i]);
866  result.value[i] /= (estimate->size);
867  }
868 
869  result.value[0] += estimate->ybase;
870  }
871  else
872  {
873  for (uint32_t i=0; i<estimate->r[0].x.size(); ++i)
874  {
875  result.value[0] += estimate->r[0].y[i];
876  result.error[0] += (estimate->r[0].y[i]*estimate->r[0].y[i]);
877  }
878  for (int i=0; i<3; ++i)
879  {
880  result.error[i] = result.error[i] - result.value[i]*result.value[i]/estimate->r[0].x.size();
881  if (estimate->r[0].x.size() > 1)
882  {
883  result.error[i] /= estimate->r[0].x.size() - 1;
884  }
885  result.error[i] = sqrt(result.error[i]);
886  result.value[i] /= estimate->r[0].x.size();
887  }
888  }
889 
890  return (result);
891 }
vector< estimatorstruc > r
Definition: mathlib.h:272
double xbase
Definition: mathlib.h:276
double ybase
Definition: mathlib.h:277
int i
Definition: rw_test.cpp:37
Estimator structure.
Definition: mathlib.h:254
uint32_t size
Definition: mathlib.h:274
uint32_t degree
Definition: mathlib.h:275
x
Definition: inputfile.py:6
double error[3]
Definition: mathlib.h:257
double value[3]
Definition: mathlib.h:256
void multisolve ( vector< vector< double > >  x,
vector< double >  y,
vector< double > &  a 
)

Perform N equation solution.

Solve a system of N equations in N unknowns. The input data is N vectors of N-1 points ( N independent and 1 dependent variable).

Parameters
xN vectors of N-1 dependent variables.
yN independent variables.
aN returned parameters
901 {
902  uint32_t order;
903  vector< vector<double> > dx(y.size()-1, vector<double>(y.size()-1));
904  vector<double> dy(y.size()-1);
905  vector<double> da(y.size()-1);
906 
907  order = y.size() - 1;
908 
909  // Order must be at least 1
910  if (order)
911  {
912  for (uint32_t r=1; r<=order; ++r)
913  {
914  dy[r-1] = y[r] * x[0][0] - y[0] * x[r][0];
915  for (uint32_t c=1; c<=order; ++c)
916  {
917  dx[r-1][c-1] = x[r][c] * x[0][0] - x[0][c] * x[r][0];
918  }
919  }
920  multisolve(dx, dy, da);
921 
922  for (uint32_t i=0; i<order; ++i)
923  {
924  a[i+1] = da[i];
925  }
926  }
927 
928  // Make this final calculation with the largest dependent in the first position
929  uint16_t bestindex = 0;
930  double bestx = 0.;
931  for (uint16_t i=0; i<y.size(); ++i)
932  {
933  if (fabs(x[i][0]) > bestx)
934  {
935  bestx = fabs(x[i][0]);
936  bestindex = i;
937  }
938  }
939  a[0] = y[bestindex];
940  for (uint32_t i=1; i<=order; ++i)
941  {
942  a[0] -= a[i] * x[bestindex][i];
943  }
944  a[0] /= x[bestindex][0];
945 
946 }
y
Definition: inputfile.py:6
void multisolve(vector< vector< double > > x, vector< double > y, vector< double > &a)
Perform N equation solution.
Definition: mathlib.cpp:900
int i
Definition: rw_test.cpp:37
x
Definition: inputfile.py:6
Definition: eci2kep_test.cpp:33
double evaluate_poly ( double  x,
vector< double >  parms 
)

Evaluate polynomial.

Return the value of the given Nth order polynomial, evaluated at the given location.

Parameters
xIndependent variable where polynomial should be evaluated.
parmsVector of parameters for Nth order polynomial to be evaluated.
Returns
Value of the variable, evaluated at the location of the independent variable.
956 {
957  double result;
958 
959  if (parms.size() < 2)
960  {
961  return 0.;
962  }
963 
964  result = parms[parms.size()-1];
965  for (uint16_t i=parms.size()-2; i<parms.size(); --i)
966  {
967  result *= x;
968  result += parms[i];
969  }
970 
971  return result;
972 }
int i
Definition: rw_test.cpp:37
x
Definition: inputfile.py:6
double evaluate_poly_slope ( double  x,
vector< double >  parms 
)

Evaluate polynomial slope.

Return the value of the 1st derivative of the given Nth order polynomial, evaluated at the given location.

Parameters
xIndependent variable where polynomial should be evaluated.
parmsVector of parameters for Nth order polynomial to be evaluated.
Returns
Value of the slope, evaluated at the location of the independent variable.
982 {
983  double result;
984 
985  if (parms.size() < 2)
986  {
987  return 0.;
988  }
989 
990  result = parms[parms.size()-1] * (parms.size()-1);
991  for (uint16_t i=parms.size()-2; i>0; --i)
992  {
993  result *= x;
994  result += i * parms[i];
995  }
996 
997  return result;
998 }
int i
Definition: rw_test.cpp:37
x
Definition: inputfile.py:6
double evaluate_poly_accel ( double  x,
vector< double >  parms 
)

Evaluate polynomial acceleration.

Return the value of the 2nd derivative of the given Nth order polynomial, evaluated at the given location.

Parameters
xIndependent variable where polynomial should be evaluated.
parmsVector of parameters for Nth order polynomial to be evaluated.
Returns
Value of the acceleration, evaluated at the location of the independent variable.
1008 {
1009  double result;
1010 
1011  if (parms.size() < 3)
1012  {
1013  return 0.;
1014  }
1015 
1016  result = parms[parms.size()-1] * (parms.size()-1) * (parms.size()-2);
1017  for (uint16_t i=parms.size()-2; i>1; --i)
1018  {
1019  result *= x;
1020  result += i * (i-1) * parms[i];
1021  }
1022 
1023  return result;
1024 }
int i
Definition: rw_test.cpp:37
x
Definition: inputfile.py:6
double evaluate_poly_jerk ( double  x,
vector< double >  parms 
)

Evaluate polynomial jerk.

Return the value of the 3rd derivative of the given Nth order polynomial, evaluated at the given location.

Parameters
xIndependent variable where polynomial should be evaluated.
parmsVector of parameters for Nth order polynomial to be evaluated.
Returns
Value of the jerk, evaluated at the location of the independent variable.
1034 {
1035  double result;
1036 
1037  if (parms.size() < 4)
1038  {
1039  return 0.;
1040  }
1041 
1042  result = parms[parms.size()-1] * (parms.size()-1) * (parms.size()-2) * (parms.size()-3);
1043  for (uint16_t i=parms.size()-2; i>2; --i)
1044  {
1045  result *= x;
1046  result += i * (i-1) * (i-2) * parms[i];
1047  }
1048 
1049  return result;
1050 }
int i
Definition: rw_test.cpp:37
x
Definition: inputfile.py:6
rvector rv_evaluate_poly ( double  x,
vector< vector< double > >  parms 
)

Evaluate vector polynomial.

Return the value of the given Nth order vector polynomial, evaluated at the given location.

Parameters
xIndependent variable where polynomial should be evaluated.
parmsVector of parameters for Nth order vector polynomial to be evaluated.
Returns
Values, evaluated at the location of the independent variable.
1060 {
1061  uvector result{};
1062 
1063  for (uint16_t ic=0; ic<parms.size(); ++ic)
1064  {
1065  if (parms[ic].size() < 2)
1066  {
1067  result.a4[ic] = 0.;
1068  }
1069  else
1070  {
1071  result.a4[ic] = parms[ic][parms[ic].size()-1];
1072  for (uint16_t i=parms[ic].size()-2; i<parms[ic].size(); --i)
1073  {
1074  result.a4[ic] *= x;
1075  result.a4[ic] += parms[ic][i];
1076  }
1077  }
1078  }
1079 
1080  return result.r;
1081 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
int i
Definition: rw_test.cpp:37
x
Definition: inputfile.py:6
double a4[4]
Definition: mathlib.h:167
rvector rv_evaluate_poly_slope ( double  x,
vector< vector< double > >  parms 
)

Evaluate vector polynomial slope.

Return the value of the 1st derivative of the given Nth order vector polynomial, evaluated at the given location.

Parameters
xIndependent variable where polynomial should be evaluated.
parmsVector of parameters for Nth order vector polynomial to be evaluated.
Returns
Values of the slope, evaluated at the location of the independent variable.
1091 {
1092  uvector result{};
1093 
1094  for (uint16_t ic=0; ic<parms.size(); ++ic)
1095  {
1096  if (parms[ic].size() < 2)
1097  {
1098  result.a4[ic] = 0.;
1099  }
1100  else
1101  {
1102  result.a4[ic] = parms[ic][parms[ic].size()-1] * (parms[ic].size()-1);
1103  for (uint16_t i=parms[ic].size()-2; i>0; --i)
1104  {
1105  result.a4[ic] *= x;
1106  result.a4[ic] += i * parms[ic][i];
1107  }
1108  }
1109  }
1110 
1111  return result.r;
1112 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
int i
Definition: rw_test.cpp:37
x
Definition: inputfile.py:6
double a4[4]
Definition: mathlib.h:167
rvector rv_evaluate_poly_accel ( double  x,
vector< vector< double > >  parms 
)

Evaluate vector polynomial acceleration.

Return the value of the 2nd derivative of the given Nth order vector polynomial, evaluated at the given location.

Parameters
xIndependent variable where polynomial should be evaluated.
parmsVector of parameters for Nth order vector polynomial to be evaluated.
Returns
Values of the acceleration, evaluated at the location of the independent variable.
1122 {
1123  uvector result{};
1124 
1125  for (uint16_t ic=0; ic<parms.size(); ++ic)
1126  {
1127  if (parms[ic].size() < 2)
1128  {
1129  result.a4[ic] = 0.;
1130  }
1131  else
1132  {
1133  result.a4[ic] = parms[ic][parms[ic].size()-1] * (parms[ic].size()-1) * (parms[ic].size()-2);
1134  for (uint16_t i=parms[ic].size()-2; i>1; --i)
1135  {
1136  result.a4[ic] *= x;
1137  result.a4[ic] += i * (i-1) * parms[ic][i];
1138  }
1139  }
1140  }
1141 
1142  return result.r;
1143 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
int i
Definition: rw_test.cpp:37
x
Definition: inputfile.py:6
double a4[4]
Definition: mathlib.h:167
rvector rv_evaluate_poly_jerk ( double  x,
vector< vector< double > >  parms 
)

Evaluate vector polynomial jerk.

Return the value of the 3rd derivative of the given Nth order vector polynomial, evaluated at the given location.

Parameters
xIndependent variable where polynomial should be evaluated.
parmsVector of parameters for Nth order vector polynomial to be evaluated.
Returns
Values of the jerk, evaluated at the location of the independent variable.
1153 {
1154  uvector result{};
1155 
1156  for (uint16_t ic=0; ic<parms.size(); ++ic)
1157  {
1158  if (parms[ic].size() < 2)
1159  {
1160  result.a4[ic] = 0.;
1161  }
1162  else
1163  {
1164  result.a4[ic] = parms[ic][parms[ic].size()-1] * (parms[ic].size()-1) * (parms[ic].size()-2) * (parms[ic].size()-3);
1165  for (uint16_t i=parms[ic].size()-2; i>2; --i)
1166  {
1167  result.a4[ic] *= x;
1168  result.a4[ic] += i * (i-1) * (i-2) * parms[ic][i];
1169  }
1170  }
1171  }
1172 
1173  return result.r;
1174 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
int i
Definition: rw_test.cpp:37
x
Definition: inputfile.py:6
double a4[4]
Definition: mathlib.h:167
gvector gv_evaluate_poly ( double  x,
vector< vector< double > >  parms 
)

Evaluate vector polynomial.

Return the value of the given Nth order vector polynomial, evaluated at the given location.

Parameters
xIndependent variable where polynomial should be evaluated.
parmsVector of parameters for Nth order vector polynomial to be evaluated.
Returns
Values, evaluated at the location of the independent variable.
1184 {
1185  uvector result{};
1186 
1187  for (uint16_t ic=0; ic<parms.size(); ++ic)
1188  {
1189  if (parms[ic].size() < 2)
1190  {
1191  result.a4[ic] = 0.;
1192  }
1193  else
1194  {
1195  result.a4[ic] = parms[ic][parms[ic].size()-1];
1196  for (uint16_t i=parms[ic].size()-2; i<parms[ic].size(); --i)
1197  {
1198  result.a4[ic] *= x;
1199  result.a4[ic] += parms[ic][i];
1200  }
1201  }
1202  }
1203 
1204  return result.g;
1205 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
int i
Definition: rw_test.cpp:37
x
Definition: inputfile.py:6
double a4[4]
Definition: mathlib.h:167
gvector gv_evaluate_poly_slope ( double  x,
vector< vector< double > >  parms 
)

Evaluate vector polynomial slope.

Return the value of the 1st derivative of the given Nth order vector polynomial, evaluated at the given location.

Parameters
xIndependent variable where polynomial should be evaluated.
parmsVector of parameters for Nth order vector polynomial to be evaluated.
Returns
Values of the slope, evaluated at the location of the independent variable.
1215 {
1216  uvector result{};
1217 
1218  for (uint16_t ic=0; ic<parms.size(); ++ic)
1219  {
1220  if (parms[ic].size() < 2)
1221  {
1222  result.a4[ic] = 0.;
1223  }
1224  else
1225  {
1226  result.a4[ic] = parms[ic][parms[ic].size()-1] * (parms[ic].size()-1);
1227  for (uint16_t i=parms[ic].size()-2; i>0; --i)
1228  {
1229  result.a4[ic] *= x;
1230  result.a4[ic] += i * parms[ic][i];
1231  }
1232  }
1233  }
1234 
1235  return result.g;
1236 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
int i
Definition: rw_test.cpp:37
x
Definition: inputfile.py:6
double a4[4]
Definition: mathlib.h:167
gvector gv_evaluate_poly_accel ( double  x,
vector< vector< double > >  parms 
)

Evaluate vector polynomial acceleration.

Return the value of the 2nd derivative of the given Nth order vector polynomial, evaluated at the given location.

Parameters
xIndependent variable where polynomial should be evaluated.
parmsVector of parameters for Nth order vector polynomial to be evaluated.
Returns
Values of the acceleration, evaluated at the location of the independent variable.
1246 {
1247  uvector result{};
1248 
1249  for (uint16_t ic=0; ic<parms.size(); ++ic)
1250  {
1251  if (parms[ic].size() < 2)
1252  {
1253  result.a4[ic] = 0.;
1254  }
1255  else
1256  {
1257  result.a4[ic] = parms[ic][parms[ic].size()-1] * (parms[ic].size()-1) * (parms[ic].size()-2);
1258  for (uint16_t i=parms[ic].size()-2; i>1; --i)
1259  {
1260  result.a4[ic] *= x;
1261  result.a4[ic] += i * (i-1) * parms[ic][i];
1262  }
1263  }
1264  }
1265 
1266  return result.g;
1267 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
int i
Definition: rw_test.cpp:37
x
Definition: inputfile.py:6
double a4[4]
Definition: mathlib.h:167
gvector gv_evaluate_poly_jerk ( double  x,
vector< vector< double > >  parms 
)

Evaluate vector polynomial jerk.

Return the value of the 3rd derivative of the given Nth order vector polynomial, evaluated at the given location.

Parameters
xIndependent variable where polynomial should be evaluated.
parmsVector of parameters for Nth order vector polynomial to be evaluated.
Returns
Values of the jerk, evaluated at the location of the independent variable.
1277 {
1278  uvector result{};
1279 
1280  for (uint16_t ic=0; ic<parms.size(); ++ic)
1281  {
1282  if (parms[ic].size() < 2)
1283  {
1284  result.a4[ic] = 0.;
1285  }
1286  else
1287  {
1288  result.a4[ic] = parms[ic][parms[ic].size()-1] * (parms[ic].size()-1) * (parms[ic].size()-2) * (parms[ic].size()-3);
1289  for (uint16_t i=parms[ic].size()-2; i>2; --i)
1290  {
1291  result.a4[ic] *= x;
1292  result.a4[ic] += i * (i-1) * (i-2) * parms[ic][i];
1293  }
1294  }
1295  }
1296 
1297  return result.g;
1298 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
int i
Definition: rw_test.cpp:37
x
Definition: inputfile.py:6
double a4[4]
Definition: mathlib.h:167
quaternion q_evaluate_poly ( double  x,
vector< vector< double > >  parms 
)

Evaluate quaternion polynomial.

Return the value of the 1st derivative of the given Nth order quaternion polynomial, evaluated at the given location.

Parameters
xIndependent variable where polynomial should be evaluated.
parmsVector of parameters for Nth order quaternion polynomial to be evaluated.
Returns
Values, evaluated at the location of the independent variable.
1308 {
1309  uvector result{};
1310 
1311  for (uint16_t ic=0; ic<parms.size(); ++ic)
1312  {
1313  if (parms[ic].size() < 2)
1314  {
1315  result.a4[ic] = 0.;
1316  }
1317  else
1318  {
1319  result.a4[ic] = parms[ic][parms[ic].size()-1];
1320  for (uint16_t i=parms[ic].size()-2; i<parms[ic].size(); --i)
1321  {
1322  result.a4[ic] *= x;
1323  result.a4[ic] += parms[ic][i];
1324  }
1325  }
1326  }
1327 
1328  return result.q;
1329 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
int i
Definition: rw_test.cpp:37
x
Definition: inputfile.py:6
double a4[4]
Definition: mathlib.h:167
quaternion q_evaluate_poly_slope ( double  x,
vector< vector< double > >  parms 
)

Evaluate quaternion polynomial slope.

Return the value of the 1st derivative of the given Nth order quaternion polynomial, evaluated at the given location.

Parameters
xIndependent variable where polynomial should be evaluated.
parmsVector of parameters for Nth order quaternion polynomial to be evaluated.
Returns
Values of the slope, evaluated at the location of the independent variable.
1339 {
1340  uvector result{};
1341 
1342  for (uint16_t ic=0; ic<parms.size(); ++ic)
1343  {
1344  if (parms[ic].size() < 2)
1345  {
1346  result.a4[ic] = 0.;
1347  }
1348  else
1349  {
1350  result.a4[ic] = parms[ic][parms[ic].size()-1] * (parms[ic].size()-1);
1351  for (uint16_t i=parms[ic].size()-2; i>0; --i)
1352  {
1353  result.a4[ic] *= x;
1354  result.a4[ic] += i * parms[ic][i];
1355  }
1356  }
1357  }
1358 
1359  return result.q;
1360 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
int i
Definition: rw_test.cpp:37
x
Definition: inputfile.py:6
double a4[4]
Definition: mathlib.h:167
quaternion q_evaluate_poly_accel ( double  x,
vector< vector< double > >  parms 
)

Evaluate quaternion polynomial acceleration.

Return the value of the 1st derivative of the given Nth order quaternion polynomial, evaluated at the given location.

Parameters
xIndependent variable where polynomial should be evaluated.
parmsVector of parameters for Nth order quaternion polynomial to be evaluated.
Returns
Values of the acceleration, evaluated at the location of the independent variable.
1370 {
1371  uvector result{};
1372 
1373  for (uint16_t ic=0; ic<parms.size(); ++ic)
1374  {
1375  if (parms[ic].size() < 2)
1376  {
1377  result.a4[ic] = 0.;
1378  }
1379  else
1380  {
1381  result.a4[ic] = parms[ic][parms[ic].size()-1] * (parms[ic].size()-1) * (parms[ic].size()-2);
1382  for (uint16_t i=parms[ic].size()-2; i>1; --i)
1383  {
1384  result.a4[ic] *= x;
1385  result.a4[ic] += i * (i-1) * parms[ic][i];
1386  }
1387  }
1388  }
1389 
1390  return result.q;
1391 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
int i
Definition: rw_test.cpp:37
x
Definition: inputfile.py:6
double a4[4]
Definition: mathlib.h:167
quaternion q_evaluate_poly_jerk ( double  x,
vector< vector< double > >  parms 
)

Evaluate quaternion polynomial jerk.

Return the value of the 1st derivative of the given Nth order quaternion polynomial, evaluated at the given location.

Parameters
xIndependent variable where polynomial should be evaluated.
parmsVector of parameters for Nth order quaternion polynomial to be evaluated.
Returns
Values of the jerk, evaluated at the location of the independent variable.
1401 {
1402  uvector result{};
1403 
1404  for (uint16_t ic=0; ic<parms.size(); ++ic)
1405  {
1406  if (parms[ic].size() < 2)
1407  {
1408  result.a4[ic] = 0.;
1409  }
1410  else
1411  {
1412  result.a4[ic] = parms[ic][parms[ic].size()-1] * (parms[ic].size()-1) * (parms[ic].size()-2) * (parms[ic].size()-3);
1413  for (uint16_t i=parms[ic].size()-2; i>2; --i)
1414  {
1415  result.a4[ic] *= x;
1416  result.a4[ic] += i * (i-1) * (i-2) * parms[ic][i];
1417  }
1418  }
1419  }
1420 
1421  return result.q;
1422 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
int i
Definition: rw_test.cpp:37
x
Definition: inputfile.py:6
double a4[4]
Definition: mathlib.h:167
vector< double > polyfit ( vector< double > &  x,
vector< double > &  y 
)

Perform general order polynomial fit.

Fit a polynomial of type: y(x) = An*x^n + An-1*x^(n-1) + ... + * A2*x^2 + A1*x + A0 where n-1 is the order of the polynomial. The input data is n points in x and n matching points in y. The fit is achieved using a Vandermonde interpolating polynomial.

Parameters
xOrder+1 number of x values
yOrder+1 number of y values
Returns
Order+1 number of polynomial coefficients
1435 {
1436  uint32_t order;
1437  vector< vector<double> > dx;
1438  vector<double> dy;
1439  vector<double> da;
1440  vector<double> a;
1441 
1442  order = x.size() - 1;
1443  a.resize(x.size());
1444 
1445  switch (order)
1446  {
1447  case 1:
1448  a[1] = (y[1]-y[0])/(x[1]-x[0]);
1449  a[0] = y[1] - a[1]*x[1];
1450  break;
1451  case 2:
1452  a[2] = ((y[0]-y[1])/(x[1]-x[0])+(y[2]-y[1])/(x[2]-x[1]))/(x[2]-x[0]);
1453  a[1] = (y[0]-y[1])/(x[0]-x[1])-(a[2]*(x[1]+x[0]));
1454  a[0] = -((x[0]*(a[1]+(a[2]*x[0])))-y[0]);
1455  break;
1456  case 3:
1457  a[3] = (((x[1] - x[0])*(x[2] - x[0])*((y[3]*(x[2] - x[1]))
1458  + (y[2]*(x[1] - x[3])) + (y[1]*(x[3] - x[2])))) - ((((y[1]
1459  - y[0])*(x[1] - x[2])) + ((x[1] - x[0])*(y[2] -
1460  y[1])))*(x[1] - x[3])*(x[2] - x[3])))/((x[1] - x[0])*(x[2]
1461  - x[0])*((x[1]*(x[0] - x[3])*((x[1]*(x[2] - x[3])) +
1462  x[3]*x[3] - x[2]*x[2])) + (x[3]*x[2]*((x[0]*x[2]) + (x[3]*(x[3] - x[0] - x[2]))))));
1463  a[2] = ((y[1] - y[0])/(x[0] - x[1])+(y[2] - y[1])/(x[2] - x[1]))/(x[2] - x[0])-(a[3]*(x[1] + x[2] + x[0]));
1464  a[1] = (y[0] - y[1] + (a[2]*(x[1]*x[1] - x[0]*x[0])) + (a[3]*(pow(x[1],3) - pow(x[0],3))))/(x[0] - x[1]);
1465  a[0] = -((x[0]*(a[1] + (x[0]*(a[2] + (a[3]*x[0]))))) - y[0]);
1466  break;
1467  default:
1468  dx.resize(order+1);
1469  for (uint32_t r=0; r<=order; ++r)
1470  {
1471  dx[r].resize(order+1);
1472  for (uint32_t c=0; c<=order; ++c)
1473  {
1474  switch (c)
1475  {
1476  case 0:
1477  dx[r][0] = 1.;
1478  break;
1479  case 1:
1480  dx[r][1] = x[r];
1481  break;
1482  case 2:
1483  dx[r][2] = x[r]*x[r];
1484  break;
1485  default:
1486  dx[r][c] = pow(x[r],c);
1487  break;
1488  }
1489  }
1490  }
1491  multisolve(dx,y, a);
1492  break;
1493  }
1494 
1495  return (a);
1496 }
y
Definition: inputfile.py:6
void multisolve(vector< vector< double > > x, vector< double > y, vector< double > &a)
Perform N equation solution.
Definition: mathlib.cpp:900
x
Definition: inputfile.py:6
Definition: eci2kep_test.cpp:33
uvector rv_fitpoly ( uvector  x,
uvector  y,
uint32_t  order 
)

Perform nth order polynomial fit.

Using the n equations: y = A0 + A1*x + A2*x^2 + ... +An*x^n, solve for the coefficients A0...An.

Parameters
xThe n values of x
yThe n corresponding values of y
orderThe order of the polynomial (< 5)
Returns
The n resulting coefficients
1508 {
1509  uvector a = {{{0.,0.,0.},0.}};
1510 
1511  if (order < 5)
1512  {
1513  switch (order)
1514  {
1515  case 1:
1516  a.a4[1] = (y.a4[1]-y.a4[0])/(x.a4[1]-x.a4[0]);
1517  a.a4[0] = y.a4[1] - a.a4[1]*x.a4[1];
1518  break;
1519  case 2:
1520  a.a4[2] = ((y.a4[0]-y.a4[1])/(x.a4[1]-x.a4[0])+(y.a4[2]-y.a4[1])/(x.a4[2]-x.a4[1]))/(x.a4[2]-x.a4[0]);
1521  a.a4[1] = (y.a4[0]-y.a4[1])/(x.a4[0]-x.a4[1])-(a.a4[2]*(x.a4[1]+x.a4[0]));
1522  a.a4[0] = -((x.a4[0]*(a.a4[1]+(a.a4[2]*x.a4[0])))-y.a4[0]);
1523  break;
1524  case 3:
1525  a.a4[3] = (((x.a4[1] - x.a4[0])*(x.a4[2] - x.a4[0])*((y.a4[3]*(x.a4[2] - x.a4[1])) + (y.a4[2]*(x.a4[1] - x.a4[3])) + (y.a4[1]*(x.a4[3] - x.a4[2])))) - ((((y.a4[1] - y.a4[0])*(x.a4[1] - x.a4[2])) + ((x.a4[1] - x.a4[0])*(y.a4[2] - y.a4[1])))*(x.a4[1] - x.a4[3])*(x.a4[2] - x.a4[3])))/((x.a4[1] - x.a4[0])*(x.a4[2] - x.a4[0])*((x.a4[1]*(x.a4[0] - x.a4[3])*((x.a4[1]*(x.a4[2] - x.a4[3])) + pow(x.a4[3],2) - pow(x.a4[2],2))) + (x.a4[3]*x.a4[2]*((x.a4[0]*x.a4[2]) + (x.a4[3]*(x.a4[3] - x.a4[0] - x.a4[2]))))));
1526  a.a4[2] = ((y.a4[1] - y.a4[0])/(x.a4[0] - x.a4[1])+(y.a4[2] - y.a4[1])/(x.a4[2] - x.a4[1]))/(x.a4[2] - x.a4[0])-(a.a4[3]*(x.a4[1] + x.a4[2] + x.a4[0]));
1527  a.a4[1] = (y.a4[0] - y.a4[1] + (a.a4[2]*(pow(x.a4[1],2) - pow(x.a4[0],2))) + (a.a4[3]*(pow(x.a4[1],3) - pow(x.a4[0],3))))/(x.a4[0] - x.a4[1]);
1528  a.a4[0] = -((x.a4[0]*(a.a4[1] + (x.a4[0]*(a.a4[2] + (a.a4[3]*x.a4[0]))))) - y.a4[0]);
1529  break;
1530  }
1531  }
1532  return (a);
1533 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
Definition: eci2kep_test.cpp:33
double a4[4]
Definition: mathlib.h:167
gj_kernel * gauss_jackson_kernel ( int32_t  order,
double  dvi 
)

Create Gauss-Jackson Integration Kernel.

Allocate space for, and initialize all the parameters common to a Gauss-Jackson integration of the requested order and step change in the independent variable.

Parameters
orderOrder at which integration will be performed.
dviStep size of independent variable.
Returns
Pointer to a structure that can be reused for multiple integrations. NULL if error.
1544 {
1545  int i, n, j,m,k;
1546  gj_kernel *gjk = nullptr;
1547 
1548  if ((gjk = (gj_kernel *)calloc(1,sizeof(gj_kernel))) == nullptr)
1549  return (gjk);
1550 
1551  gjk->dvi = dvi;
1552  gjk->order = order;
1553  gjk->dvi2 = dvi * dvi;
1554  gjk->horder = order/2;
1555 
1556  if ((gjk->binom = (int32_t **)calloc(order+2,sizeof(int32_t *))) == nullptr)
1557  {
1559  return nullptr;
1560  }
1561 
1562  for (i=0; i<order+2; i++)
1563  {
1564  if ((gjk->binom[i] = (int32_t *)calloc(order+2,sizeof(int32_t))) == nullptr)
1565  {
1567  return nullptr;
1568  }
1569  }
1570 
1571  if ((gjk->alpha = (double **)calloc(order+2,sizeof(double *))) == nullptr)
1572  {
1574  return nullptr;
1575  }
1576 
1577  for (i=0; i<order+2; i++)
1578  {
1579  if ((gjk->alpha[i] = (double *)calloc(order+1,sizeof(double))) == nullptr)
1580  {
1582  return nullptr;
1583  }
1584  }
1585 
1586  if ((gjk->beta = (double **)calloc(order+2,sizeof(double *))) == nullptr)
1587  {
1589  return nullptr;
1590  }
1591 
1592  for (i=0; i<order+2; i++)
1593  {
1594  if ((gjk->beta[i] = (double *)calloc(order+1,sizeof(double))) == nullptr)
1595  {
1597  return nullptr;
1598  }
1599  }
1600 
1601  if ((gjk->c = (double *)calloc(order+3,sizeof(double *))) == nullptr)
1602  {
1604  return nullptr;
1605  }
1606 
1607  if ((gjk->gam = (double *)calloc(order+2,sizeof(double *))) == nullptr)
1608  {
1610  return nullptr;
1611  }
1612 
1613  if ((gjk->q = (double *)calloc(order+3,sizeof(double *))) == nullptr)
1614  {
1616  return nullptr;
1617  }
1618 
1619  if ((gjk->lam = (double *)calloc(order+3,sizeof(double *))) == nullptr)
1620  {
1622  return nullptr;
1623  }
1624 
1625  for (m=0; m<order+2; m++)
1626  {
1627  for (i=0; i<order+2; i++)
1628  {
1629  if (m > i)
1630  gjk->binom[m][i] = 0;
1631  else
1632  {
1633  if (m == i)
1634  gjk->binom[m][i] = 1;
1635  else
1636  {
1637  if (m == 0)
1638  gjk->binom[m][i] = 1;
1639  else
1640  gjk->binom[m][i] = gjk->binom[m-1][i-1] + gjk->binom[m][i-1];
1641  }
1642  }
1643  }
1644  }
1645 
1646  gjk->c[0] = 1.;
1647  for (n=1; n<order+3; n++)
1648  {
1649  gjk->c[n] = 0.;
1650  for (i=0; i<=n-1; i++)
1651  {
1652  gjk->c[n] -= gjk->c[i] / (n+1-i);
1653  }
1654  }
1655 
1656  gjk->gam[0] = gjk->c[0];
1657  for (i=1; i<gjk->order+2; i++)
1658  {
1659  gjk->gam[i] = gjk->gam[i-1] + gjk->c[i];
1660  }
1661 
1662  for (i=0; i<gjk->order+1; i++)
1663  {
1664  gjk->beta[gjk->order+1][i] = gjk->gam[i+1];
1665  gjk->beta[gjk->order][i] = gjk->c[i+1];
1666  for (j=gjk->order-1; j>=0; j--)
1667  {
1668  if (!i)
1669  gjk->beta[j][i] = gjk->beta[j+1][i];
1670  else
1671  gjk->beta[j][i] = gjk->beta[j+1][i] - gjk->beta[j+1][i-1];
1672  }
1673  }
1674 
1675  gjk->q[0] = 1.;
1676  for (i=1; i<gjk->order+3; i++)
1677  {
1678  gjk->q[i] = 0.;
1679  for (k=0; k<=i; k++)
1680  {
1681  gjk->q[i] += gjk->c[k]*gjk->c[i-k];
1682  }
1683  }
1684 
1685  gjk->lam[0] = gjk->q[0];
1686  for (i=1; i<gjk->order+3; i++)
1687  {
1688  gjk->lam[i] = gjk->lam[i-1] + gjk->q[i];
1689  }
1690 
1691  for (i=0; i<gjk->order+1; i++)
1692  {
1693  gjk->alpha[gjk->order+1][i] = gjk->lam[i+2];
1694  gjk->alpha[gjk->order][i] = gjk->q[i+2];
1695  for (j=gjk->order-1; j>=0; j--)
1696  {
1697  if (!i)
1698  gjk->alpha[j][i] = gjk->alpha[j+1][i];
1699  else
1700  gjk->alpha[j][i] = gjk->alpha[j+1][i] - gjk->alpha[j+1][i-1];
1701  }
1702  }
1703 
1704  return (gjk);
1705 }
double dvi
Definition: mathlib.h:182
double dvi2
Definition: mathlib.h:183
int32_t ** binom
Definition: mathlib.h:184
int i
Definition: rw_test.cpp:37
void gauss_jackson_dekernel(gj_kernel *gjk)
Free Gauss-Jackson Integration Kernel.
Definition: mathlib.cpp:1711
pxnxm element cube
Definition: mathlib.h:178
double * q
Definition: mathlib.h:189
double * c
Definition: mathlib.h:187
int32_t order
Definition: mathlib.h:180
double * gam
Definition: mathlib.h:188
double * lam
Definition: mathlib.h:190
double ** alpha
Definition: mathlib.h:185
double ** beta
Definition: mathlib.h:186
int32_t horder
Definition: mathlib.h:181
void gauss_jackson_dekernel ( gj_kernel gjk)

Free Gauss-Jackson Integration Kernel.

Free the space fully or partially allocated represented by the supplied kernel pointer.

Parameters
gjkPointer to Gauss-Jackson Kernel
1712 {
1713  int i;
1714 
1715  if (gjk == nullptr)
1716  return;
1717 
1718  if (gjk->binom != nullptr)
1719  {
1720  for (i=0; i<gjk->order+2; i++)
1721  {
1722  if (gjk->binom[i] != nullptr)
1723  free(gjk->binom[i]);
1724  }
1725  free(gjk->binom);
1726  }
1727 
1728  if (gjk->alpha != nullptr)
1729  {
1730  for (i=0; i<gjk->order+2; i++)
1731  {
1732  if (gjk->alpha[i] != nullptr)
1733  free(gjk->alpha[i]);
1734  }
1735  free(gjk->alpha);
1736  }
1737 
1738  if (gjk->beta != nullptr)
1739  {
1740  for (i=0; i<gjk->order+2; i++)
1741  {
1742  if (gjk->beta[i] != nullptr)
1743  free(gjk->beta[i]);
1744  }
1745  free(gjk->beta);
1746  }
1747 
1748  if (gjk->c != nullptr)
1749  free(gjk->c);
1750 
1751  if (gjk->gam != nullptr)
1752  free(gjk->gam);
1753 
1754  if (gjk->q != nullptr)
1755  free(gjk->q);
1756 
1757  if (gjk->lam != nullptr)
1758  free(gjk->lam);
1759 
1760  free(gjk);
1761 }
int32_t ** binom
Definition: mathlib.h:184
int i
Definition: rw_test.cpp:37
double * q
Definition: mathlib.h:189
double * c
Definition: mathlib.h:187
int32_t order
Definition: mathlib.h:180
double * gam
Definition: mathlib.h:188
double * lam
Definition: mathlib.h:190
double ** alpha
Definition: mathlib.h:185
double ** beta
Definition: mathlib.h:186
gj_instance * gauss_jackson_instance ( gj_kernel kern,
int32_t  axes,
void(*)(double vi, double *vd0, double *vd2, int32_t axes)  calc_vd2 
)

Initialize an Instance of a Gauss-Jackson Integrator.

Set up the integration specific variables for an integration using the provided kernel. The integration will be performed simultaneously for the indicated number of axes. The 2nd derivative for a specific set of independent and dependent variables is calculated using the provided function.

Parameters
kernPointer to Gauss-Jackson Kernel previously initialized for order and step size.
axesNumber of axes to integrate simultaneously.
calc_vd2Function that will calculate 2nd derivative for all dependent variables at once.
Returns
Pointer to a Gauss-Jackson Instance, initialized for use with this kernel. NULL if error.
1775 {
1776  int i, m;
1777  gj_instance *gji;
1778 
1779  if (kern == nullptr)
1780  return nullptr;
1781 
1782  if ((gji = (gj_instance *)calloc(1,sizeof(gj_instance))) != nullptr)
1783  {
1784  gji->kern = kern;
1785  gji->axes = axes;
1786  gji->calc_vd2 = calc_vd2;
1787  if ((gji->vi = (double *)calloc(kern->order+2,sizeof(double))) != nullptr)
1788  {
1789  if ((gji->steps = (gj_step **)calloc(axes,sizeof(gj_step *))) != nullptr)
1790  {
1791  for (i=0; i<axes; i++)
1792  {
1793  if ((gji->steps[i] = gauss_jackson_step(kern)) != nullptr)
1794  {
1795  continue;
1796  }
1797  break;
1798  }
1799  if (i == axes)
1800  return (gji);
1801  for (m=0; m<i; m++)
1802  {
1803  gauss_jackson_destep(kern,gji->steps[m]);
1804  }
1805  free(gji->steps);
1806  free(gji->vi);
1807  free(gji);
1808  return nullptr;
1809  }
1810  else
1811  {
1812  free(gji->vi);
1813  free(gji);
1814  return nullptr;
1815  }
1816  }
1817  else
1818  {
1819  free(gji);
1820  return nullptr;
1821  }
1822  }
1823  else
1824  return nullptr;
1825 }
Gauss-Jackson Integration Instance.
Definition: mathlib.h:216
int i
Definition: rw_test.cpp:37
double * vi
Pointer to array of independent variables; order+2.
Definition: mathlib.h:225
int32_t order
Definition: mathlib.h:180
int32_t axes
Number of axes of integration.
Definition: mathlib.h:221
gj_step * gauss_jackson_step(gj_kernel *kern)
Initialize a Step of a Gauss-Jackson integrator.
Definition: mathlib.cpp:1833
void(* calc_vd2)(double vi, double *vd0, double *vd2, int32_t axes)
Pointer to a function that will calculate the 2nd derivative given a axes dependent and one independe...
Definition: mathlib.h:227
void gauss_jackson_destep(gj_kernel *kern, gj_step *step)
Destroy a Step of a Gauss-Jackson integrator.
Definition: mathlib.cpp:1888
gj_kernel * kern
Kernel Pointer.
Definition: mathlib.h:219
gj_step ** steps
Pointer to array of steps; one for each axis, order+2 for the orders.
Definition: mathlib.h:223
Gauss-Jackson Integration Step.
Definition: mathlib.h:196
gj_step * gauss_jackson_step ( gj_kernel kern)

Initialize a Step of a Gauss-Jackson integrator.

Allocate the space, and set up the steps for a single axis of a Gauss-Jackson integrator of specified order.

Parameters
kernKernel for this Gauss-Jackson integrator.
Returns
Pointer to the Step.
1834 {
1835  int i,j, m;
1836 
1837  gj_step *step;
1838 
1839  if ((step = (gj_step *)calloc(kern->order+2,sizeof(gj_step))) != nullptr)
1840  {
1841  for (j=0; j<kern->order+2; j++)
1842  {
1843  if ((step[j].a = (double *)calloc(kern->order+1,sizeof(double))) != nullptr)
1844  {
1845  if ((step[j].b = (double *)calloc(kern->order+1,sizeof(double))) != nullptr)
1846  continue;
1847  free(step[j].a);
1848  }
1849  break;
1850  }
1851  if (j < kern->order+2)
1852  {
1853  for (i=0; i<j; i++)
1854  {
1855  free(step[i].a);
1856  free(step[i].b);
1857  }
1858  free(step);
1859  return nullptr;
1860  }
1861  }
1862 
1863  for (j=0; j<kern->order+2; j++)
1864  {
1865  for (m=0; m<kern->order+1; m++)
1866  {
1867  step[j].a[kern->order-m] = step[j].b[kern->order-m] = 0.;
1868  for (i=m; i<=kern->order; i++)
1869  {
1870  step[j].a[kern->order-m] += kern->alpha[j][i] * kern->binom[m][i];
1871  step[j].b[kern->order-m] += kern->beta[j][i] * kern->binom[m][i];
1872  }
1873  step[j].a[kern->order-m] *= pow(-1.,m);
1874  step[j].b[kern->order-m] *= pow(-1.,m);
1875  if (kern->order-m == j)
1876  step[j].b[kern->order-m] += .5;
1877  }
1878  }
1879  return (step);
1880 }
int32_t ** binom
Definition: mathlib.h:184
int i
Definition: rw_test.cpp:37
long b
Definition: jpegint.h:371
int32_t order
Definition: mathlib.h:180
Definition: eci2kep_test.cpp:33
double ** alpha
Definition: mathlib.h:185
Gauss-Jackson Integration Step.
Definition: mathlib.h:196
double ** beta
Definition: mathlib.h:186
double * b
Definition: mathlib.h:205
double * a
Definition: mathlib.h:204
void gauss_jackson_destep ( gj_kernel kern,
gj_step step 
)

Destroy a Step of a Gauss-Jackson integrator.

Deallocate the space a single axis of a Gauss-Jackson integrator of specified order.

Parameters
kernKernel for this Gauss-Jackson integrator.
stepPointer to the step to be destroyed.
1889 {
1890  int j;
1891 
1892  if (step != nullptr)
1893  {
1894  for (j=0; j<kern->order+2; j++)
1895  {
1896  if (step[j].a != nullptr)
1897  {
1898  free(step[j].a);
1899  }
1900  if (step[j].b != nullptr)
1901  {
1902  free(step[j].b);
1903  }
1904  }
1905  free(step);
1906  }
1907 }
long b
Definition: jpegint.h:371
int32_t order
Definition: mathlib.h:180
Definition: eci2kep_test.cpp:33
int gauss_jackson_setstep ( gj_instance gji,
double  vi,
double *  vd0,
double *  vd1,
double *  vd2,
int32_t  istep 
)

Set Independent and Dependent variables for Gauss-Jackson step.

Set the Independent and Dependent variables for a particular step of a particular Gauss-Jackson integration Instance.

Parameters
gjiPointer to Gauss-Jackson Instance.
viPointer to Independent variable.
vd0Pointer to array of Dependent variable for all axes at given step.
vd1Pointer to array of First derivative of Dependent variable for all axes at given step.
vd2Pointer to array of Second derivative of Dependent variable for all axes at given step.
istepInteger number of step.
Returns
0, otherwise negative error.
1921 {
1922  int i;
1923 
1924  if (gji == nullptr)
1925  return (MATH_ERROR_GJ_UNDEFINED);
1926 
1927  if (istep > gji->kern->order+2)
1928  return (MATH_ERROR_GJ_OUTOFRANGE);
1929 
1930  gji->vi[istep] = vi;
1931  for (i=0; i<gji->axes; i++)
1932  {
1933  gji->steps[i][istep].vd0 = vd0[i];
1934  gji->steps[i][istep].vd1 = vd1[i];
1935  gji->steps[i][istep].vd2 = vd2[i];
1936  }
1937  return 0;
1938 }
#define MATH_ERROR_GJ_OUTOFRANGE
Definition: cosmos-errno.h:116
int i
Definition: rw_test.cpp:37
double vd1
Dependent variable 1st derivative.
Definition: mathlib.h:201
double * vi
Pointer to array of independent variables; order+2.
Definition: mathlib.h:225
double vd0
Dependent variable.
Definition: mathlib.h:199
#define MATH_ERROR_GJ_UNDEFINED
Definition: cosmos-errno.h:115
int32_t order
Definition: mathlib.h:180
int32_t axes
Number of axes of integration.
Definition: mathlib.h:221
double vd2
Dependent variable 2nd derivative.
Definition: mathlib.h:203
gj_kernel * kern
Kernel Pointer.
Definition: mathlib.h:219
gj_step ** steps
Pointer to array of steps; one for each axis, order+2 for the orders.
Definition: mathlib.h:223
int gauss_jackson_getstep ( gj_instance gji,
double *  vi,
double *  vd0,
double *  vd1,
double *  vd2,
int32_t  istep 
)

Get Independent and Dependent variables for Gauss-Jackson step.

Get the Independent and Dependent variables for a particular step of a particular Gauss-Jackson integration Instance.

Parameters
gjiPointer to Gauss-Jackson Instance
viPointer to Independent variable
vd0Pointer to array of Dependent variable for all axes at given step
vd1Pointer to array of First derivative of Dependent variable for all axes at given step
vd2Pointer to array of Second derivative of Dependent variable for all axes at given step
istepIndex of step to get.
Returns
0, otherwise negative error.
1952 {
1953  int i;
1954  if (gji == nullptr)
1955  return (MATH_ERROR_GJ_UNDEFINED);
1956 
1957  if (istep > gji->kern->order+2)
1958  return (MATH_ERROR_GJ_OUTOFRANGE);
1959 
1960  *vi = gji->vi[istep];
1961  for (i=0; i<gji->axes; i++)
1962  {
1963  vd0[i] = gji->steps[i][istep].vd0;
1964  vd1[i] = gji->steps[i][istep].vd1;
1965  vd2[i] = gji->steps[i][istep].vd2;
1966  }
1967  return 0;
1968 }
#define MATH_ERROR_GJ_OUTOFRANGE
Definition: cosmos-errno.h:116
int i
Definition: rw_test.cpp:37
double vd1
Dependent variable 1st derivative.
Definition: mathlib.h:201
double * vi
Pointer to array of independent variables; order+2.
Definition: mathlib.h:225
double vd0
Dependent variable.
Definition: mathlib.h:199
#define MATH_ERROR_GJ_UNDEFINED
Definition: cosmos-errno.h:115
int32_t order
Definition: mathlib.h:180
int32_t axes
Number of axes of integration.
Definition: mathlib.h:221
double vd2
Dependent variable 2nd derivative.
Definition: mathlib.h:203
gj_kernel * kern
Kernel Pointer.
Definition: mathlib.h:219
gj_step ** steps
Pointer to array of steps; one for each axis, order+2 for the orders.
Definition: mathlib.h:223
void gauss_jackson_preset ( gj_instance gji)

Converge all axes of a Gauss-Jackson integrator prior to propagation.

Converge all axes of a Gauss-Jackson integrator prior to propagation.

Parameters
gjiPointer to a Gauss-Jackson Integration Instance that has valid values
1975 {
1976  int32_t ccount, cflag=1, k, n, i;
1977  int32_t gj_order, gj_2order;
1978  double gj_dvi, gj_dvi2;
1979  static double *oldvd2 = nullptr;
1980  static double *newvd0 = nullptr;
1981  static double *newvd2 = nullptr;
1982  static int32_t axes = 0;
1983 
1984  if (newvd0 == nullptr || axes < gji->axes)
1985  {
1986  if (newvd0)
1987  {
1988  free(newvd0);
1989  free(oldvd2);
1990  free(newvd2);
1991  }
1992  axes = gji->axes;
1993  oldvd2 = static_cast <double *>(calloc(axes,sizeof(double)));
1994  newvd0 = static_cast <double *>(calloc(axes,sizeof(double)));
1995  newvd2 = static_cast <double *>(calloc(axes,sizeof(double)));
1996  }
1997 
1998  gj_2order = gji->kern->horder;
1999  gj_order = gji->kern->order;
2000  gj_dvi = gji->kern->dvi;
2001  gj_dvi2 = gji->kern->dvi2;
2002 
2003  ccount = 0;
2004  do
2005  {
2006  for (i=0; i<axes; i++)
2007  {
2008  gji->steps[i][gj_2order].s = gji->steps[i][gj_2order].vd1/gj_dvi;
2009  for (k=0; k<=gj_order; k++)
2010  gji->steps[i][gj_2order].s -= gji->steps[i][gj_2order].b[k] * gji->steps[i][k].vd2;
2011 
2012  for (n=1; n<=gj_2order; n++)
2013  {
2014  gji->steps[i][gj_2order+n].s = gji->steps[i][gj_2order+n-1].s + (gji->steps[i][gj_2order+n].vd2+gji->steps[i][gj_2order+n-1].vd2)/2;
2015  gji->steps[i][gj_2order-n].s = gji->steps[i][gj_2order-n+1].s - (gji->steps[i][gj_2order-n].vd2+gji->steps[i][gj_2order-n+1].vd2)/2;
2016  }
2017 
2018  gji->steps[i][gj_2order].ss = gji->steps[i][gj_2order].vd0/gj_dvi2;
2019 
2020  for (k=0; k<=gj_order; k++)
2021  gji->steps[i][gj_2order].ss -= gji->steps[i][gj_2order].a[k] * gji->steps[i][k].vd2;
2022 
2023  for (n=1; n<=gj_2order; n++)
2024  {
2025  gji->steps[i][gj_2order+n].ss = gji->steps[i][gj_2order+n-1].ss + gji->steps[i][gj_2order+n-1].s + (gji->steps[i][gj_2order+n-1].vd2)/2;
2026  gji->steps[i][gj_2order-n].ss = gji->steps[i][gj_2order-n+1].ss - gji->steps[i][gj_2order-n+1].s + (gji->steps[i][gj_2order-n+1].vd2)/2;
2027  }
2028 
2029  for (n=0; n<=gj_order; n++)
2030  {
2031  if (n == gj_2order)
2032  continue;
2033  gji->steps[i][n].sb = gji->steps[i][n].sa = 0.;
2034  for (k=0; k<=gj_order; k++)
2035  {
2036  gji->steps[i][n].sb += gji->steps[i][n].b[k] * gji->steps[i][k].vd2;
2037  gji->steps[i][n].sa += gji->steps[i][n].a[k] * gji->steps[i][k].vd2;
2038  }
2039  }
2040  }
2041 
2042  for (n=1; n<=gj_2order; n++)
2043  {
2044  for (i=-1; i<2; i+=2)
2045  {
2046  // Calculate new probable position and velocity
2047  for (k=0; k<axes; k++)
2048  {
2049  gji->steps[k][gj_2order+i*n].vd1 = gj_dvi * (gji->steps[k][gj_2order+i*n].s + gji->steps[k][gj_2order+i*n].sb);
2050  newvd0[k] = gji->steps[k][gj_2order+i*n].vd0 = gj_dvi2 * (gji->steps[k][gj_2order+i*n].ss + gji->steps[k][gj_2order+i*n].sa);
2051  oldvd2[k] = gji->steps[k][gj_2order+i*n].vd2;
2052  }
2053 
2054  // Calculate new acceleration
2055  gji->calc_vd2(gji->vi[gj_2order+i*n],newvd0,newvd2,axes);
2056 
2057  // Compare acceleration at new position to previous iteration
2058  cflag = 0;
2059 
2060  for (k=0; k<axes; k++)
2061  {
2062  gji->steps[k][gj_2order+i*n].vd2 = newvd2[k];
2063  if (fabs(oldvd2[k]-newvd2[k])>1e-14)
2064  cflag = 1;
2065  }
2066 
2067  }
2068  }
2069  ccount++;
2070  } while (ccount<10 && cflag);
2071 }
double dvi
Definition: mathlib.h:182
Definition: eci2kep_test.cpp:33
double dvi2
Definition: mathlib.h:183
int i
Definition: rw_test.cpp:37
double vd1
Dependent variable 1st derivative.
Definition: mathlib.h:201
double sb
Definition: mathlib.h:210
double * vi
Pointer to array of independent variables; order+2.
Definition: mathlib.h:225
double vd0
Dependent variable.
Definition: mathlib.h:199
int32_t order
Definition: mathlib.h:180
double s
Current guess of dependent variable.
Definition: mathlib.h:207
int32_t axes
Number of axes of integration.
Definition: mathlib.h:221
double sa
Definition: mathlib.h:209
void(* calc_vd2)(double vi, double *vd0, double *vd2, int32_t axes)
Pointer to a function that will calculate the 2nd derivative given a axes dependent and one independe...
Definition: mathlib.h:227
double ss
Definition: mathlib.h:208
double vd2
Dependent variable 2nd derivative.
Definition: mathlib.h:203
gj_kernel * kern
Kernel Pointer.
Definition: mathlib.h:219
gj_step ** steps
Pointer to array of steps; one for each axis, order+2 for the orders.
Definition: mathlib.h:223
double * b
Definition: mathlib.h:205
double * a
Definition: mathlib.h:204
int32_t horder
Definition: mathlib.h:181
void gauss_jackson_extrapolate ( gj_instance gji,
double  target 
)

Propagate Gauss-Jackson integration.

Propagate the Gauss-Jackson integration Instance for all axes. Will iterate a step at a time until target value of indpendent variable has been reached.

Parameters
gjiGauss-Jackson instance
targetTarget value of the independent variable
2080 {
2081  int32_t chunks, i, j, k;
2082  int32_t gj_order, gj_2order;
2083  double gj_dvi, gj_dvi2;
2084  static double *newvd0 = nullptr;
2085  static double *newvd2 = nullptr;
2086  static int32_t axes = 0;
2087 
2088  if (newvd0 == nullptr || axes < gji->axes)
2089  {
2090  if (newvd0)
2091  {
2092  free(newvd0);
2093  free(newvd2);
2094  }
2095  axes = gji->axes;
2096  newvd0 = (double *)calloc(axes,sizeof(double));
2097  newvd2 = (double *)calloc(axes,sizeof(double));
2098  }
2099 
2100 
2101  gj_order = gji->kern->order;
2102  gj_2order = gj_order/2;
2103  gj_dvi = gji->kern->dvi;
2104  gj_dvi2 = gj_dvi * gj_dvi;
2105 
2106  chunks = (int32_t)(.5 + (target-gji->vi[gj_2order])/gji->kern->dvi);
2107  for (i=0; i<chunks; i++)
2108  {
2109  gji->vi[gj_order+1] = gji->vi[gj_order] + gji->kern->dvi;
2110 
2111  // For each axis
2112  for (j=0; j<axes; j++)
2113  {
2114  // Calculate S(order/2+1)
2115  gji->steps[j][gj_order+1].ss = gji->steps[j][gj_order].ss + gji->steps[j][gj_order].s + gji->steps[j][gj_order].vd2/2.;
2116 
2117  // Calculate Sum(order/2+1) for a and b
2118  gji->steps[j][gj_order+1].sb = gji->steps[j][gj_order+1].sa = 0.;
2119  for (k=0; k<=gj_order; k++)
2120  {
2121  gji->steps[j][gj_order+1].sb += gji->steps[j][gj_order+1].b[k] * gji->steps[j][k].vd2;
2122  gji->steps[j][gj_order+1].sa += gji->steps[j][gj_order+1].a[k] * gji->steps[j][k].vd2;
2123  }
2124 
2125  // Calculate pos.v(order/2+1)
2126  gji->steps[j][gj_order+1].vd1 = gj_dvi * (gji->steps[j][gj_order].s + gji->steps[j][gj_order].vd2/2. + gji->steps[j][gj_order+1].sb);
2127 
2128  // Calculate pos.s(order/2+1)
2129  newvd0[j] = gji->steps[j][gj_order+1].vd0 = gj_dvi2 * (gji->steps[j][gj_order+1].ss + gji->steps[j][gj_order+1].sa);
2130  }
2131 
2132  // Perform positional and attitude accelerations at new position
2133  gji->calc_vd2(gji->vi[gj_order+1],newvd0,newvd2,axes);
2134 
2135  for (j=0; j<axes; j++)
2136  {
2137  // Set new 2nd derivative
2138  gji->steps[j][gj_order+1].vd2 = newvd2[j];
2139 
2140  // Calculate s(order/2+1)
2141  gji->steps[j][gj_order+1].s = gji->steps[j][gj_order].s + (gji->steps[j][gj_order].vd2+gji->steps[j][gj_order+1].vd2)/2.;
2142 
2143  // Shift everything over 1
2144  for (k=0; k<=gj_order; k++)
2145  gji->steps[j][k] = gji->steps[j][k+1];
2146  }
2147  // Shift times over 1
2148  for (k=0; k<=gj_order; k++)
2149  gji->vi[k] = gji->vi[k+1];
2150 
2151  }
2152 }
double dvi
Definition: mathlib.h:182
int i
Definition: rw_test.cpp:37
static antstruc target
Definition: agent_antenna.cpp:160
double vd1
Dependent variable 1st derivative.
Definition: mathlib.h:201
double sb
Definition: mathlib.h:210
double * vi
Pointer to array of independent variables; order+2.
Definition: mathlib.h:225
double vd0
Dependent variable.
Definition: mathlib.h:199
int32_t order
Definition: mathlib.h:180
double s
Current guess of dependent variable.
Definition: mathlib.h:207
int32_t axes
Number of axes of integration.
Definition: mathlib.h:221
double sa
Definition: mathlib.h:209
void(* calc_vd2)(double vi, double *vd0, double *vd2, int32_t axes)
Pointer to a function that will calculate the 2nd derivative given a axes dependent and one independe...
Definition: mathlib.h:227
double ss
Definition: mathlib.h:208
double vd2
Dependent variable 2nd derivative.
Definition: mathlib.h:203
gj_kernel * kern
Kernel Pointer.
Definition: mathlib.h:219
gj_step ** steps
Pointer to array of steps; one for each axis, order+2 for the orders.
Definition: mathlib.h:223
double * b
Definition: mathlib.h:205
double * a
Definition: mathlib.h:204
double fixangle ( double  angle)

Limit angle to range 0-2PI.

Ensure that angle represents equivalent value within a range of 0 to 2*PI.

Parameters
angleInput angle in radians.
Returns
Output angle adjusted to range 0 to 2*PI.
2160 {
2161  double result;
2162 
2163  result = fmod(angle,D2PI);
2164 
2165  return (result >= 0.)?result:result+D2PI;
2166 }
const double D2PI
Double precision 2*PI.
Definition: math/constants.h:16
double actan ( double  y,
double  x 
)

ArcTan, limited to range 0-2PI.

Calculate Arc Tangent of y/x, ennsuring that angle represents equivalent value within a range of 0 to 2*PI.

Parameters
yNumerator.
xDenominator.
Returns
Output angle adjusted to range 0 to 2*PI.
2175 {
2176  double actan;
2177 
2178  actan = atan2(y,x);
2179  if (actan < 0. )
2180  actan += D2PI;
2181  return (actan);
2182 }
y
Definition: inputfile.py:6
double actan(double y, double x)
ArcTan, limited to range 0-2PI.
Definition: mathlib.cpp:2174
x
Definition: inputfile.py:6
const double D2PI
Double precision 2*PI.
Definition: math/constants.h:16
double fixprecision ( double  number,
double  prec 
)

Limit precision.

Ensure that number represents equivalent value rounded to nearest precision. The calculation performed is prec * round(number / prec).

Parameters
numberOriginal value.
precPrecision to round to.
Returns
Rounded value.
2192 {
2193  return (prec*round(number/prec));
2194 }
uint16_t calc_crc16ccitt ( uint8_t *  buf,
int  size,
bool  lsb 
)

Calculate CRC-16-CCITT.

Calculate 16-bit CCITT CRC for the indicated buffer and number of bytes. For the lsb variant, the initial shift register value is 0xffff, and the calculation starts with the LSB, so the Polynomial is 0x8408. For the msb, the initial value is 0 and the Polynomial is 0x1021.

Parameters
bufbytes to calculate on
sizenumber of bytes
Returns
calculated CRC
2207 {
2208  uint16_t crc;
2209  uint8_t ch;
2210 
2211  if (lsb)
2212  {
2213  crc = 0xffff;
2214  }
2215  else
2216  {
2217  crc = 0;
2218  }
2219  for (uint16_t i=0; i<size; i++)
2220  {
2221  ch = buf[i];
2222  for (uint16_t j=0; j<8; j++)
2223  {
2224  if (lsb)
2225  {
2226  crc = (crc >> 1) ^ (((ch^crc)&0x01)?CRC16CCITTLSB:0);
2227  ch >>= 1;
2228  }
2229  else
2230  {
2231  crc = (crc << 1) ^ (((ch&0x80)^((crc&0x8000)>>8))?CRC16CCITTMSB:0);
2232  ch <<= 1;
2233  }
2234  }
2235  }
2236  return (crc);
2237 }
#define CRC16CCITTLSB
Definition: mathlib.h:131
int i
Definition: rw_test.cpp:37
png_uint_32 crc
Definition: png.c:2173
#define CRC16CCITTMSB
Definition: mathlib.h:128
char buf[128]
Definition: rw_test.cpp:40
rvector drotate ( quaternion  q,
rvector  v 
)

Rotate a row vector using a quaternion.

Rotate a row vector within one coordinate system using the provided left quaternion.

Parameters
qQuaternion representing the rotation.
vRow vector to be rotated.
Returns
Rotated row vector in the same system.
2262 {
2263  // TODO: remove uvector, use quaternion
2264  uvector t = {{{0.,0.,0.},0.}};
2265 
2266  t.q = q_fmult(q,q_fmult(v,q_conjugate(q)));
2267 
2268  // TODO: how is this supposed to work?
2269  return (t.r);
2270 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
rvector r
Definition: mathlib.h:161
quaternion q_conjugate(quaternion q)
Definition: vector.cpp:1010
quaternion q
Definition: mathlib.h:157
quaternion q_fmult(rvector r1, quaternion q2)
rvector quaternion multiply
Definition: vector.cpp:1079
rvector rotate_q ( quaternion  q,
rvector  v 
)
2273 {
2274  return drotate(q, v);
2275 }
rvector drotate(quaternion q, rvector v)
Rotate a row vector using a quaternion.
Definition: mathlib.cpp:2261
cvector drotate ( quaternion  q,
cvector  v 
)

Rotate a cartesian vector using a quaternion.

Rotate a cartesian vector from one coordinate system to another using the provided quaternion.

Parameters
qquaternion representing the rotation
vcartesian vector to be rotated
Returns
cartesian vector in the rotated system
2285 {
2286  uvector qt{};
2287 
2288  qt.c = v;
2289  qt.q.w = 0.0;
2290 
2291  qt.q = q_fmult(q,q_fmult(qt.r,q_conjugate(q)));
2292 
2293  return (qt.c);
2294 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
quaternion q_conjugate(quaternion q)
Definition: vector.cpp:1010
quaternion q_fmult(rvector r1, quaternion q2)
rvector quaternion multiply
Definition: vector.cpp:1079
cvector c
Definition: mathlib.h:162
cvector rotate_q ( quaternion  q,
cvector  v 
)
2297 {
2298  return drotate(q, v);
2299 }
rvector drotate(quaternion q, rvector v)
Rotate a row vector using a quaternion.
Definition: mathlib.cpp:2261
rvector irotate ( quaternion  q,
rvector  v 
)

Indirectly rotate a row vector using a quaternion.

Indirectly rotate a row vector from one coordinate system to another using the provided left quaternion.

Parameters
qquaternion representing the intrinsic rotation
vrow vector to be rotated
Returns
row vector in the intrinsically rotated system
2309 {
2310  uvector t = {{{0.,0.,0.},0.}};
2311 
2312  t.q = q_fmult(q_conjugate(q),q_fmult(v,q));
2313 
2314  return (t.r);
2315 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
rvector r
Definition: mathlib.h:161
quaternion q_conjugate(quaternion q)
Definition: vector.cpp:1010
quaternion q
Definition: mathlib.h:157
quaternion q_fmult(rvector r1, quaternion q2)
rvector quaternion multiply
Definition: vector.cpp:1079
rvector transform_q ( quaternion  q,
rvector  v 
)
2318 {
2319  return irotate(q, v);
2320 }
rvector irotate(quaternion q, rvector v)
Indirectly rotate a row vector using a quaternion.
Definition: mathlib.cpp:2308
cvector irotate ( quaternion  q,
cvector  v 
)

Indirectly rotate a cartesian vector using a quaternion.

Indirectly rotate a cartesian vector from one coordinate system to another using the provided left quaternion.

Parameters
qquaternion representing the intrinsic rotation
vcartesian vector to be rotated
Returns
cartesian vector in the intrinsically rotated system
2330 {
2331  uvector qt{};
2332 
2333  qt.c = v;
2334  qt.q.w = 0.0;
2335 
2336  qt.q = q_fmult(q_conjugate(q),q_fmult(qt.r,q));
2337 
2338  return (qt.c);
2339 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
quaternion q_conjugate(quaternion q)
Definition: vector.cpp:1010
quaternion q_fmult(rvector r1, quaternion q2)
rvector quaternion multiply
Definition: vector.cpp:1079
cvector c
Definition: mathlib.h:162
quaternion q_change_between_cv ( cvector  from,
cvector  to 
)

Create rotation quaternion from 2 vectors.

Generate the quaternion that represents a rotation of from one cartesian vector to a second cartesian vector.

Parameters
frominitial cartesian vector
tofinal cartesian vector
Returns
quaternion that can be used to rotate points
2351 {
2352  uvector rq{};
2353  cvector vec1, vec2;
2354 
2355  normalize_cv(from);
2356  normalize_cv(to);
2357 
2358  if (length_cv(cv_add(from,to)) < 1e-14)
2359  {
2360  vec1.x = rand();
2361  vec1.y = rand();
2362  vec1.z = rand();
2363  normalize_cv(vec1);
2364  vec2 = cv_cross(vec1,to);
2365  normalize_cv(vec2);
2366  if (length_cv(vec2)<D_SMALL)
2367  {
2368  vec1.x = rand();
2369  vec1.y = rand();
2370  vec1.z = rand();
2371  normalize_cv(vec1);
2372  vec2 = cv_cross(vec1,to);
2373  normalize_cv(vec2);
2374  }
2375  rq.c = vec2;
2376  rq.q.w = 0.;
2377  }
2378  else
2379  {
2380  rq.c = cv_cross(from,to);
2381  rq.q.w = 1. + dot_cv(from,to);
2382  }
2383 
2384  normalize_q(&rq.q);
2385  return (rq.q);
2386 }
double y
Y value.
Definition: vector.h:114
Quaternion/Rvector Union.
Definition: mathlib.h:155
Definition: eci2kep_test.cpp:33
cvector cv_cross(cvector a, cvector b)
Take cross product of two vectors.
Definition: vector.cpp:613
cvector cv_add(cvector a, cvector b)
Add two vectors.
Definition: vector.cpp:554
double x
X value.
Definition: vector.h:112
void normalize_q(quaternion *q)
Definition: vector.cpp:961
void normalize_cv(cvector &v)
Normalize cartesian vector in place, i.e. divides it by its own norm.
Definition: vector.cpp:474
double length_cv(cvector v)
Definition: vector.cpp:638
const double D_SMALL
Definition: math/constants.h:40
double z
Z value.
Definition: vector.h:116
double dot_cv(cvector a, cvector b)
Definition: vector.cpp:623
3 element cartesian vector
Definition: vector.h:107
quaternion q_change_between_rv ( rvector  from,
rvector  to 
)
2389 {
2390  uvector uvfrom{}, uvto{};
2391  quaternion result;
2392 
2393  uvfrom.r = from;
2394  uvto.r = to;
2395  result = q_change_between_cv(uvfrom.c, uvto.c);
2396 
2397  return result;
2398 }
Quaternion/Rvector Union.
Definition: mathlib.h:155
quaternion q_change_between_cv(cvector from, cvector to)
Create rotation quaternion from 2 vectors.
Definition: mathlib.cpp:2350
Quaternion, scalar last, using x, y, z.
Definition: vector.h:402
cmatrix cm_change_between_cv ( cvector  from,
cvector  to 
)

Create rotation matrix from 2 vectors.

Generate the direction cosine matrix that represents a rotation from one cartesian vector to a second cartesian vector.

Parameters
frominitial cartesian vector
tofinal cartesian vector
Returns
direction cosine matrix that can be used to rotate points
2408 {
2409  //cmatrix m;
2410  //m = cm_quaternion2dcm(q_change_between_cv(from,to));
2411  //return (m);
2412 
2413  return cm_quaternion2dcm(q_change_between_cv(from,to));
2414 }
quaternion q_change_between_cv(cvector from, cvector to)
Create rotation quaternion from 2 vectors.
Definition: mathlib.cpp:2350
cmatrix cm_quaternion2dcm(quaternion q)
Definition: rotation.cpp:51
double evaluate_poly ( double  x,
rvector  parms 
)
double evaluate_poly_slope ( double  x,
rvector  parms 
)
double evaluate_poly_accel ( double  x,
rvector  parms 
)
double evaluate_poly_jerk ( double  x,
rvector  parms 
)
uint16_t calc_crc16ccitt_lsb ( vector< uint8_t >  buf)
3020 {
3021  uint16_t crc;
3022  uint8_t ch;
3023 
3024  crc = 0xffff;
3025 
3026  for (uint16_t i=0; i<buf.size(); i++)
3027  {
3028  ch = buf[i];
3029  for (uint16_t j=0; j<8; j++)
3030  {
3031  crc = (crc >> 1) ^ (((ch^crc)&0x01)?CRC16CCITTLSB:0);
3032  ch >>= 1;
3033 
3034  }
3035  }
3036  return (crc);
3037 }
#define CRC16CCITTLSB
Definition: mathlib.h:131
int i
Definition: rw_test.cpp:37
png_uint_32 crc
Definition: png.c:2173
char buf[128]
Definition: rw_test.cpp:40
uint16_t calc_crc16ccitt_msb ( vector< uint8_t >  buf)
3040 {
3041  uint16_t crc;
3042  uint8_t ch;
3043 
3044  crc = 0;
3045 
3046  for (uint16_t i=0; i<buf.size(); i++)
3047  {
3048  ch = buf[i];
3049  for (uint16_t j=0; j<8; j++)
3050  {
3051 
3052  crc = (crc << 1) ^ (((ch&0x80)^((crc&0x8000)>>8))?CRC16CCITTMSB:0);
3053  ch <<= 1;
3054 
3055  }
3056  }
3057  return (crc);
3058 }
int i
Definition: rw_test.cpp:37
png_uint_32 crc
Definition: png.c:2173
#define CRC16CCITTMSB
Definition: mathlib.h:128
char buf[128]
Definition: rw_test.cpp:40
rvector rv_mmult ( rmatrix  m,
rvector  v 
)

Multiply rmatrix by rvector.

Multiply 3x3 rmatrix by 3 element rvector (treated as a column order vector).

Parameters
mmatrix to multiply by, in rmatrix form.
vvector to be tranformed, in rvector form.
Returns
multiplied vector, in rvector format.
42 {
43  rvector o;
44 
45  o.col[0] = m.row[0].col[0]*v.col[0] + m.row[0].col[1]*v.col[1] + m.row[0].col[2]*v.col[2];
46  o.col[1] = m.row[1].col[0]*v.col[0] + m.row[1].col[1]*v.col[1] + m.row[1].col[2]*v.col[2];
47  o.col[2] = m.row[2].col[0]*v.col[0] + m.row[2].col[1]*v.col[1] + m.row[2].col[2]*v.col[2];
48 
49  return (o);
50 }
3 element generic row vector
Definition: vector.h:53
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
rvector operator* ( rmatrix  m,
rvector  v 
)
54 {
55  return rv_mmult(m, v);
56 }
rvector rv_mmult(rmatrix m, rvector v)
Multiply rmatrix by rvector.
Definition: matrix.cpp:41
rvector rv_diag ( rmatrix  a)

Matrix diagonal.

rvector representing the diagonal of a rmatrix

Parameters
a:rmatrix to get diagonal from.
Returns
Diagonal
64 {
65  rvector b;
66 
67  b.col[0] = a.row[0].col[0];
68  b.col[1] = a.row[1].col[1];
69  b.col[2] = a.row[2].col[2];
70 
71  return (b);
72 }
3 element generic row vector
Definition: vector.h:53
long b
Definition: jpegint.h:371
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
cvector cv_mmult ( cmatrix  m,
cvector  v 
)

Multiply cartesian vector by cartesian matrix.

Multiply 3 element cartesian vector by 3x3 cartesian matrix

Parameters
vvector to be tranformed, in cvector form
mmatrix to multiply by, in cmatrix form
Returns
multiplied vector, in cvector format
81 {
82  cvector o;
83 
84  o.x = m.r1.x*v.x + m.r1.y*v.y + m.r1.z*v.z;
85  o.y = m.r2.x*v.x + m.r2.y*v.y + m.r2.z*v.z;
86  o.z = m.r3.x*v.x + m.r3.y*v.y + m.r3.z*v.z;
87 
88  return (o);
89 }
double y
Y value.
Definition: vector.h:114
cvector r1
Row 1.
Definition: matrix.h:99
double x
X value.
Definition: vector.h:112
double z
Z value.
Definition: vector.h:116
cvector r3
Row 3.
Definition: matrix.h:103
3 element cartesian vector
Definition: vector.h:107
cvector r2
Row 2.
Definition: matrix.h:101
rmatrix rm_diag ( rvector  a)

Diagonal rmatrix Creates an rmatrix whose diagonal is filled with the supplied rvector.

Parameters
a
Returns
New rmatrix.
98 {
99  rmatrix b;
100 
101  b.row[0].col[0] = a.col[0];
102  b.row[1].col[1] = a.col[1];
103  b.row[2].col[2] = a.col[2];
104 
105  b.row[0].col[1] = b.row[0].col[2] = b.row[1].col[0] = b.row[1].col[2] = b.row[2].col[0] = b.row[2].col[1] = 0.;
106 
107  return (b);
108 }
long b
Definition: jpegint.h:371
3x3 element generic matrix
Definition: matrix.h:41
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
rmatrix rm_eye ( )

Identity rmatrix.

rmatrix with diagonal elements set to one and all others set to zero

Returns
Identity matrix
115 {
116  rmatrix mat = {{1.,0.,0.}, {0.,1.,0.}, {0.,0.,1.}};
117 
118  return (mat);
119 }
3x3 element generic matrix
Definition: matrix.h:41
rmatrix rm_zero ( )

Zero filled rmatrix.

rmatrix with all elements set to zero

Returns
Zero filled matrix
126 {
127  rmatrix mat;
128  return (mat);
129 }
3x3 element generic matrix
Definition: matrix.h:41
double norm_rm ( rmatrix  mat)

rmatrix norm. Calculates the Norm of the supplied rmatrix

Parameters
mat
Returns
Calculated norm.
138 {
139  return fmax(
140  norm_rv(mat.row[0]),
141  fmax(
142  norm_rv(mat.row[1]),
143  norm_rv(mat.row[2])
144  )
145  );
146 }
double norm_rv(rvector vec)
Infinite norm of row vector.
Definition: vector.cpp:765
rvector row[3]
Definition: matrix.h:43
double trace_rm ( rmatrix  mat)

rmatrix Trace Calculates the trace of the supplied rmatrix.

Parameters
mat
Returns
Calculated trace.
155 {
156  double trace;
157  trace = mat.row[0].col[0] + mat.row[1].col[1] + mat.row[2].col[2];
158  return (trace);
159 }
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
rmatrix rm_transpose ( rmatrix  a)

rmatrix Transpose. Calculate the transpose of the supplied rmatrix.

Parameters
a
Returns
Calculated transpose.
173 {
174  rmatrix b;
175 
176  b.row[0].col[0] = a.row[0].col[0];
177  b.row[0].col[1] = a.row[1].col[0];
178  b.row[0].col[2] = a.row[2].col[0];
179 
180  b.row[1].col[0] = a.row[0].col[1];
181  b.row[1].col[1] = a.row[1].col[1];
182  b.row[1].col[2] = a.row[2].col[1];
183 
184  b.row[2].col[0] = a.row[0].col[2];
185  b.row[2].col[1] = a.row[1].col[2];
186  b.row[2].col[2] = a.row[2].col[2];
187 
188  return (b);
189 }
long b
Definition: jpegint.h:371
3x3 element generic matrix
Definition: matrix.h:41
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
rmatrix rm_mmult ( rmatrix  a,
rmatrix  b 
)

rmatrix Matrix Product

Multiply two row order matrices together.

Parameters
afirst rmatrix.
bsecond rmatrix.
Returns
product rmatrix.
199 {
200  rmatrix mat;
201 
202  mat.row[0].col[0] = a.row[0].col[0]*b.row[0].col[0] + a.row[0].col[1]*b.row[1].col[0] + a.row[0].col[2]*b.row[2].col[0];
203  mat.row[0].col[1] = a.row[0].col[0]*b.row[0].col[1] + a.row[0].col[1]*b.row[1].col[1] + a.row[0].col[2]*b.row[2].col[1];
204  mat.row[0].col[2] = a.row[0].col[0]*b.row[0].col[2] + a.row[0].col[1]*b.row[1].col[2] + a.row[0].col[2]*b.row[2].col[2];
205 
206  mat.row[1].col[0] = a.row[1].col[0]*b.row[0].col[0] + a.row[1].col[1]*b.row[1].col[0] + a.row[1].col[2]*b.row[2].col[0];
207  mat.row[1].col[1] = a.row[1].col[0]*b.row[0].col[1] + a.row[1].col[1]*b.row[1].col[1] + a.row[1].col[2]*b.row[2].col[1];
208  mat.row[1].col[2] = a.row[1].col[0]*b.row[0].col[2] + a.row[1].col[1]*b.row[1].col[2] + a.row[1].col[2]*b.row[2].col[2];
209 
210  mat.row[2].col[0] = a.row[2].col[0]*b.row[0].col[0] + a.row[2].col[1]*b.row[1].col[0] + a.row[2].col[2]*b.row[2].col[0];
211  mat.row[2].col[1] = a.row[2].col[0]*b.row[0].col[1] + a.row[2].col[1]*b.row[1].col[1] + a.row[2].col[2]*b.row[2].col[1];
212  mat.row[2].col[2] = a.row[2].col[0]*b.row[0].col[2] + a.row[2].col[1]*b.row[1].col[2] + a.row[2].col[2]*b.row[2].col[2];
213 
214  return (mat);
215 }
3x3 element generic matrix
Definition: matrix.h:41
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
rmatrix rm_mult ( rmatrix  a,
rmatrix  b 
)

Element-wise rmatrix multiplication.

Parameters
afirst rmatrix.
bsecond rmatrix.
Returns
Element-wise product rmatrix.
224 {
225  rmatrix mat;
226  rvector *va, *vb;
227 
228  va = &a.row[0];
229  vb = &b.row[0];
230  mat.row[0] = rv_mult(*va,*vb);
231  va = &a.row[1];
232  vb = &b.row[1];
233  mat.row[1] = rv_mult(*va,*vb);
234  va = &a.row[2];
235  vb = &b.row[2];
236  mat.row[2] = rv_mult(*va,*vb);
237 
238  return (mat);
239 }
3 element generic row vector
Definition: vector.h:53
3x3 element generic matrix
Definition: matrix.h:41
rvector rv_mult(rvector a, rvector b)
Multiply two row vectors.
Definition: vector.cpp:347
rvector row[3]
Definition: matrix.h:43
rmatrix rm_smult ( double  a,
rmatrix  b 
)

Scalar rmatrix multiplication.

Parameters
aScalar to multiply each element by.
brmatrix.
Returns
Scalar product rmatrix.
248 {
249  rmatrix mat;
250  rvector *vb;
251 
252  vb = &b.row[0];
253  mat.row[0] = rv_smult(a,*vb);
254  vb = &b.row[1];
255  mat.row[1] = rv_smult(a,*vb);
256  vb = &b.row[2];
257  mat.row[2] = rv_smult(a,*vb);
258 
259  return (mat);
260 }
3 element generic row vector
Definition: vector.h:53
rvector rv_smult(double a, rvector b)
Multiply row vector by scalar.
Definition: vector.cpp:266
3x3 element generic matrix
Definition: matrix.h:41
Definition: eci2kep_test.cpp:33
rvector row[3]
Definition: matrix.h:43
rmatrix rm_add ( rmatrix  a,
rmatrix  b 
)

rmatrix addition. Sum of two rmatrix values.

Parameters
aFirst term.
bSecond term.
Returns
Matrix sum.
270 {
271  rmatrix mat;
272  rvector *va, *vb;
273 
274  va = &a.row[0];
275  vb = &b.row[0];
276  mat.row[0] = rv_add(*va,*vb);
277  va = &a.row[1];
278  vb = &b.row[1];
279  mat.row[1] = rv_add(*va,*vb);
280  va = &a.row[2];
281  vb = &b.row[2];
282  mat.row[2] = rv_add(*va,*vb);
283 
284  return (mat);
285 }
rvector rv_add(rvector a, rvector b)
Add two row vectors.
Definition: vector.cpp:299
3 element generic row vector
Definition: vector.h:53
3x3 element generic matrix
Definition: matrix.h:41
rvector row[3]
Definition: matrix.h:43
rmatrix rm_sub ( rmatrix  a,
rmatrix  b 
)

rmatrix subtraction. Subtract two rmatrix values.

Parameters
aMinuend.
bSubtrahend.
Returns
Difference.
295 {
296  rmatrix mat;
297  rvector *va, *vb;
298 
299  va = &a.row[0];
300  vb = &b.row[0];
301  mat.row[0] = rv_sub(*va,*vb);
302  va = &a.row[1];
303  vb = &b.row[1];
304  mat.row[1] = rv_sub(*va,*vb);
305  va = &a.row[2];
306  vb = &b.row[2];
307  mat.row[2] = rv_sub(*va,*vb);
308 
309  return (mat);
310 }
3 element generic row vector
Definition: vector.h:53
3x3 element generic matrix
Definition: matrix.h:41
rvector rv_sub(rvector a, rvector b)
Subtract two vectors.
Definition: vector.cpp:315
rvector row[3]
Definition: matrix.h:43
rmatrix rm_square ( rmatrix  a)

Square rmatrix.

Square a rmatrix matrix by matrix multiplying it by itself.

Parameters
amatrix to be squared
Returns
squared matrix
318 {
319  rmatrix b;
320  rvector r1, r2, r3, c1, c2, c3;
321 
322  r1 = a.row[0];
323  r2 = a.row[1];
324  r3 = a.row[2];
325 
326  c1.col[0] = a.row[0].col[0];
327  c1.col[1] = a.row[1].col[0];
328  c1.col[2] = a.row[2].col[0];
329 
330  c2.col[0] = a.row[0].col[1];
331  c2.col[1] = a.row[1].col[1];
332  c2.col[2] = a.row[2].col[1];
333 
334  c3.col[0] = a.row[0].col[2];
335  c3.col[1] = a.row[1].col[2];
336  c3.col[2] = a.row[2].col[2];
337 
338  b.row[0].col[0] = sum_rv(rv_mult(r1,c1));
339  b.row[0].col[1] = sum_rv(rv_mult(r1,c2));
340  b.row[0].col[2] = sum_rv(rv_mult(r1,c3));
341 
342  b.row[1].col[0] = sum_rv(rv_mult(r2,c1));
343  b.row[1].col[1] = sum_rv(rv_mult(r2,c2));
344  b.row[1].col[2] = sum_rv(rv_mult(r2,c3));
345 
346  b.row[2].col[0] = sum_rv(rv_mult(r3,c1));
347  b.row[2].col[1] = sum_rv(rv_mult(r3,c2));
348  b.row[2].col[2] = sum_rv(rv_mult(r3,c3));
349 
350  return (b);
351 }
3 element generic row vector
Definition: vector.h:53
double sum_rv(rvector vec)
Sum elements of a row vector.
Definition: vector.cpp:775
long b
Definition: jpegint.h:371
3x3 element generic matrix
Definition: matrix.h:41
rvector rv_mult(rvector a, rvector b)
Multiply two row vectors.
Definition: vector.cpp:347
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
rmatrix rm_change_around_x ( double  angle)

Rotation matrix for X axis.

Create the DCM that represents a rotation of the given angle around the X axis.

Parameters
angleAngle of rotation in radians
Returns
Resulting DCM
359 {
360  return rm_change_around(1,angle);
361 }
rmatrix rm_change_around(int axis, double angle)
Rotation matrix for indicated axis.
Definition: matrix.cpp:397
rmatrix rm_change_around_y ( double  angle)

Rotation matrix for Y axis.

Create the DCM that represents a rotation of the given angle around the Y axis.

Parameters
angleAngle of rotation in radians
Returns
Resulting DCM
369 {
370  //rmatrix a = {{{{0.}}}};
371  //a = rm_change_around(2, angle);
372  //return (a);
373 
374  return rm_change_around(2,angle);
375 }
rmatrix rm_change_around(int axis, double angle)
Rotation matrix for indicated axis.
Definition: matrix.cpp:397
rmatrix rm_change_around_z ( double  angle)

Rotation matrix for Z axis.

Create the DCM that represents a rotation of the given angle around the Z axis.

Parameters
angleAngle of rotation in radians
Returns
Resulting DCM
383 {
384  //rmatrix a = {{{{0.}}}};
385  //a = rm_change_around(3, angle);
386  //return (a);
387 
388  return rm_change_around(3,angle);
389 }
rmatrix rm_change_around(int axis, double angle)
Rotation matrix for indicated axis.
Definition: matrix.cpp:397
rmatrix rm_change_around ( int  axis,
double  angle 
)

Rotation matrix for indicated axis.

Create the DCM that represents a rotation of the given angle around the indicated axis.

Parameters
axisAxis of rotation: 1=X, 2=Y, 3=Z
angleAngle of rotation in radians
Returns
Resulting DCM
398 {
399  rmatrix a = {{1.,0.,0.}, {0.,1.,0.}, {0.,0.,1.}};
400 
401  switch (axis)
402  {
403  case 1:
404  a.row[1].col[1] = a.row[2].col[2] = cos(angle);
405  a.row[2].col[1] = sin(angle);
406  a.row[1].col[2] = -(a.row[2].col[1]);
407  break;
408  case 2:
409  a.row[0].col[0] = a.row[2].col[2] = cos(angle);
410  a.row[0].col[2] = sin(angle);
411  a.row[2].col[0] = -(a.row[0].col[2]);
412  break;
413  case 3:
414  a.row[0].col[0] = a.row[1].col[1] = cos(angle);
415  a.row[1].col[0] = sin(angle);
416  a.row[0].col[1] = -(a.row[1].col[0]);
417  break;
418  }
419 
420  return (a);
421 }
3x3 element generic matrix
Definition: matrix.h:41
Definition: eci2kep_test.cpp:33
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
cmatrix cm_from_rm ( rmatrix  matrix)

cmatrix from rmatrix

Converts 3x3 matrix in rmatrix form to cmatrix form.

Parameters
matrixRow major matrix to convert
Returns
Converted matrix in row order form
429 {
430  cmatrix cm;
431 
432  cm.r1.x = matrix.row[0].col[0];
433  cm.r1.y = matrix.row[0].col[1];
434  cm.r1.z = matrix.row[0].col[2];
435  cm.r2.x = matrix.row[1].col[0];
436  cm.r2.y = matrix.row[1].col[1];
437  cm.r2.z = matrix.row[1].col[2];
438  cm.r3.x = matrix.row[2].col[0];
439  cm.r3.y = matrix.row[2].col[1];
440  cm.r3.z = matrix.row[2].col[2];
441 
442  return (cm);
443 }
double y
Y value.
Definition: vector.h:114
3x3 element cartesian matrix
Definition: matrix.h:96
cvector r1
Row 1.
Definition: matrix.h:99
double x
X value.
Definition: vector.h:112
double z
Z value.
Definition: vector.h:116
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
cvector r3
Row 3.
Definition: matrix.h:103
cvector r2
Row 2.
Definition: matrix.h:101
cmatrix cm_diag ( cvector  a)

Diagonal cmatrix Creates an cmatrix whose diagonal is filled with the supplied cvector.

Parameters
a
Returns
New cmatrix.
454 {
455  cmatrix b;
456 
457  b.r1.x = a.x;
458  b.r2.y = a.y;
459  b.r3.z = a.z;
460 
461  b.r1.y = b.r1.z = b.r2.x = b.r2.z = b.r3.x = b.r3.y = 0.;
462 
463  return (b);
464 }
double y
Y value.
Definition: vector.h:114
3x3 element cartesian matrix
Definition: matrix.h:96
cvector r1
Row 1.
Definition: matrix.h:99
double x
X value.
Definition: vector.h:112
long b
Definition: jpegint.h:371
double z
Z value.
Definition: vector.h:116
cvector r3
Row 3.
Definition: matrix.h:103
cvector r2
Row 2.
Definition: matrix.h:101
cmatrix cm_eye ( )

Identity cmatrix.

cmatrix with diagonal elements set to one and all others set to zero

Returns
Identity matrix
471 {
472  cmatrix mat = {{1.,0.,0.},{0.,1.,0.},{0.,0.,1.}};
473 
474  return (mat);
475 }
3x3 element cartesian matrix
Definition: matrix.h:96
cmatrix cm_zero ( )

Zero filled cmatrix.

cmatrix with all elements set to zero

Returns
Zero filled matrix
482 {
483  cmatrix mat = {{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};
484  return (mat);
485 }
3x3 element cartesian matrix
Definition: matrix.h:96
double norm_cm ( cmatrix  mat)

cmatrix norm. Calculates the Norm of the supplied cmatrix

Parameters
mat
Returns
Calculated norm.
494 {
495  double norm;
496 
497  norm = norm_cv(mat.r1);
498  norm = fmax(norm,norm_cv(mat.r2));
499  norm = fmax(norm,norm_cv(mat.r3));
500 
501  return (norm);
502 }
cvector r1
Row 1.
Definition: matrix.h:99
double norm_cv(cvector v)
Definition: vector.cpp:699
cvector r3
Row 3.
Definition: matrix.h:103
cvector r2
Row 2.
Definition: matrix.h:101
double trace_cm ( cmatrix  mat)

cmatrix Trace Calculates the trace of the supplied cmatrix.

Parameters
mat
Returns
Calculated trace.
511 {
512  double trace;
513 
514  trace = mat.r1.x + mat.r2.y + mat.r3.z;
515 
516  return (trace);
517 }
double y
Y value.
Definition: vector.h:114
cvector r1
Row 1.
Definition: matrix.h:99
double x
X value.
Definition: vector.h:112
double z
Z value.
Definition: vector.h:116
cvector r3
Row 3.
Definition: matrix.h:103
cvector r2
Row 2.
Definition: matrix.h:101
cmatrix cm_transpose ( cmatrix  a)

cmatrix Transpose. Calculate the transpose of the supplied cmatrix.

Parameters
a
Returns
Calculated transpose.
526 {
527  cmatrix b;
528 
529  b.r1.x = a.r1.x;
530  b.r1.y = a.r2.x;
531  b.r1.z = a.r3.x;
532 
533  b.r2.x = a.r1.y;
534  b.r2.y = a.r2.y;
535  b.r2.z = a.r3.y;
536 
537  b.r3.x = a.r1.z;
538  b.r3.y = a.r2.z;
539  b.r3.z = a.r3.z;
540 
541  return (b);
542 }
double y
Y value.
Definition: vector.h:114
3x3 element cartesian matrix
Definition: matrix.h:96
cvector r1
Row 1.
Definition: matrix.h:99
double x
X value.
Definition: vector.h:112
long b
Definition: jpegint.h:371
double z
Z value.
Definition: vector.h:116
cvector r3
Row 3.
Definition: matrix.h:103
cvector r2
Row 2.
Definition: matrix.h:101
cvector cv_diag ( cmatrix  a)

Matrix diagonal.

cvector representing the diagonal of a cmatrix

Parameters
a:cmatrix to get diagonal from.
Returns
Diagonal
550 {
551  cvector b;
552 
553  b.x = a.r1.x;
554  b.y = a.r2.y;
555  b.z = a.r3.z;
556 
557  return (b);
558 }
double y
Y value.
Definition: vector.h:114
cvector r1
Row 1.
Definition: matrix.h:99
double x
X value.
Definition: vector.h:112
long b
Definition: jpegint.h:371
double z
Z value.
Definition: vector.h:116
cvector r3
Row 3.
Definition: matrix.h:103
3 element cartesian vector
Definition: vector.h:107
cvector r2
Row 2.
Definition: matrix.h:101
cmatrix cm_mmult ( cmatrix  a,
cmatrix  b 
)

cmatrix Matrix Product

Multiply two cartesian matrices together.

Parameters
afirst cartesian matrix
bsecond cartesian matrix
Returns
product cartesian matrix
567 {
568  cmatrix mat;
569 
570  mat.r1.x = a.r1.x*b.r1.x + a.r1.y*b.r2.x + a.r1.z*b.r3.x;
571  mat.r1.y = a.r1.x*b.r1.y + a.r1.y*b.r2.y + a.r1.z*b.r3.y;
572  mat.r1.z = a.r1.x*b.r1.z + a.r1.y*b.r2.z + a.r1.z*b.r3.z;
573 
574  mat.r2.x = a.r2.x*b.r1.x + a.r2.y*b.r2.x + a.r2.z*b.r3.x;
575  mat.r2.y = a.r2.x*b.r1.y + a.r2.y*b.r2.y + a.r2.z*b.r3.y;
576  mat.r2.z = a.r2.x*b.r1.z + a.r2.y*b.r2.z + a.r2.z*b.r3.z;
577 
578  mat.r3.x = a.r3.x*b.r1.x + a.r3.y*b.r2.x + a.r3.z*b.r3.x;
579  mat.r3.y = a.r3.x*b.r1.y + a.r3.y*b.r2.y + a.r3.z*b.r3.y;
580  mat.r3.z = a.r3.x*b.r1.z + a.r3.y*b.r2.z + a.r3.z*b.r3.z;
581 
582  return (mat);
583 }
double y
Y value.
Definition: vector.h:114
3x3 element cartesian matrix
Definition: matrix.h:96
cvector r1
Row 1.
Definition: matrix.h:99
double x
X value.
Definition: vector.h:112
double z
Z value.
Definition: vector.h:116
cvector r3
Row 3.
Definition: matrix.h:103
cvector r2
Row 2.
Definition: matrix.h:101
cmatrix cm_mult ( cmatrix  a,
cmatrix  b 
)

Element-wise cmatrix multiplication.

Parameters
afirst cmatrix.
bsecond cmatrix.
Returns
Element-wise product cmatrix.
592 {
593  cmatrix mat;
594  cvector *va, *vb;
595 
596  va = &a.r1;
597  vb = &b.r1;
598  mat.r1 = cv_mult(*va,*vb);
599  va = &a.r2;
600  vb = &b.r2;
601  mat.r2 = cv_mult(*va,*vb);
602  va = &a.r3;
603  vb = &b.r3;
604  mat.r3 = cv_mult(*va,*vb);
605 
606  return (mat);
607 }
3x3 element cartesian matrix
Definition: matrix.h:96
cvector r1
Row 1.
Definition: matrix.h:99
cvector cv_mult(cvector a, cvector b)
Multiply two vectors.
Definition: vector.cpp:602
cvector r3
Row 3.
Definition: matrix.h:103
3 element cartesian vector
Definition: vector.h:107
cvector r2
Row 2.
Definition: matrix.h:101
cmatrix cm_smult ( double  a,
cmatrix  b 
)

Scalar cmatrix multiplication.

Parameters
aScalar to multiply each element by.
bcmatrix.
Returns
Scalar product cmatrix.
616 {
617  cmatrix mat;
618  cvector *vb;
619 
620  vb = &b.r1;
621  mat.r1 = cv_smult(a,*vb);
622  vb = &b.r2;
623  mat.r2 = cv_smult(a,*vb);
624  vb = &b.r3;
625  mat.r3 = cv_smult(a,*vb);
626 
627  return (mat);
628 }
3x3 element cartesian matrix
Definition: matrix.h:96
cvector r1
Row 1.
Definition: matrix.h:99
Definition: eci2kep_test.cpp:33
cvector cv_smult(double a, cvector b)
Multiply vector by scalar.
Definition: vector.cpp:521
cvector r3
Row 3.
Definition: matrix.h:103
3 element cartesian vector
Definition: vector.h:107
cvector r2
Row 2.
Definition: matrix.h:101
cmatrix cm_add ( cmatrix  a,
cmatrix  b 
)

cmatrix addition. Sum of two cmatrix values.

Parameters
aFirst term.
bSecond term.
Returns
Matrix sum.
638 {
639  cmatrix mat;
640  cvector *va, *vb;
641 
642  va = &a.r1;
643  vb = &b.r1;
644  mat.r1 = cv_add(*va,*vb);
645  va = &a.r2;
646  vb = &b.r2;
647  mat.r2 = cv_add(*va,*vb);
648  va = &a.r3;
649  vb = &b.r3;
650  mat.r3 = cv_add(*va,*vb);
651 
652  return (mat);
653 }
3x3 element cartesian matrix
Definition: matrix.h:96
cvector r1
Row 1.
Definition: matrix.h:99
cvector cv_add(cvector a, cvector b)
Add two vectors.
Definition: vector.cpp:554
cvector r3
Row 3.
Definition: matrix.h:103
3 element cartesian vector
Definition: vector.h:107
cvector r2
Row 2.
Definition: matrix.h:101
cmatrix cm_sub ( cmatrix  a,
cmatrix  b 
)

cmatrix subtraction. Subtract two cmatrix values.

Parameters
aMinuend.
bSubtrahend.
Returns
Difference.
663 {
664  cmatrix mat;
665  cvector *va, *vb;
666 
667  va = &a.r1;
668  vb = &b.r1;
669  mat.r1 = cv_sub(*va,*vb);
670  va = &a.r2;
671  vb = &b.r2;
672  mat.r2 = cv_sub(*va,*vb);
673  va = &a.r3;
674  vb = &b.r3;
675  mat.r3 = cv_sub(*va,*vb);
676 
677  return (mat);
678 }
3x3 element cartesian matrix
Definition: matrix.h:96
cvector cv_sub(cvector a, cvector b)
Subtract two vectors.
Definition: vector.cpp:570
cvector r1
Row 1.
Definition: matrix.h:99
cvector r3
Row 3.
Definition: matrix.h:103
3 element cartesian vector
Definition: vector.h:107
cvector r2
Row 2.
Definition: matrix.h:101
cmatrix cm_square ( cmatrix  a)

Square cmatrix.

Square a cmatrix by matrix multiplying it by itself.

Parameters
amatrix to be squared
Returns
squared matrix
686 {
687  cmatrix b;
688  cvector r1, r2, r3, c1, c2, c3;
689 
690  r1 = a.r1;
691  r2 = a.r2;
692  r3 = a.r3;
693 
694  c1.x = a.r1.x;
695  c1.y = a.r2.x;
696  c1.z = a.r3.x;
697 
698  c2.x = a.r1.y;
699  c2.y = a.r2.y;
700  c2.z = a.r3.y;
701 
702  c3.x = a.r1.z;
703  c3.y = a.r2.z;
704  c3.z = a.r3.z;
705 
706  b.r1.x = sum_cv(cv_mult(r1,c1));
707  b.r1.y = sum_cv(cv_mult(r1,c2));
708  b.r1.z = sum_cv(cv_mult(r1,c3));
709 
710  b.r2.x = sum_cv(cv_mult(r2,c1));
711  b.r2.y = sum_cv(cv_mult(r2,c2));
712  b.r2.z = sum_cv(cv_mult(r2,c3));
713 
714  b.r3.x = sum_cv(cv_mult(r3,c1));
715  b.r3.y = sum_cv(cv_mult(r3,c2));
716  b.r3.z = sum_cv(cv_mult(r3,c3));
717 
718  return (b);
719 }
double y
Y value.
Definition: vector.h:114
3x3 element cartesian matrix
Definition: matrix.h:96
cvector r1
Row 1.
Definition: matrix.h:99
double x
X value.
Definition: vector.h:112
long b
Definition: jpegint.h:371
double z
Z value.
Definition: vector.h:116
double sum_cv(cvector vec)
Definition: vector.cpp:708
cvector cv_mult(cvector a, cvector b)
Multiply two vectors.
Definition: vector.cpp:602
cvector r3
Row 3.
Definition: matrix.h:103
3 element cartesian vector
Definition: vector.h:107
cvector r2
Row 2.
Definition: matrix.h:101
cmatrix cm_change_around_x ( double  angle)

Rotation matrix for X axis.

Create the DCM that represents a rotation of the given angle around the X axis.

Parameters
angleAngle of rotation in radians
Returns
Resulting DCM
728 {
729  cmatrix a;
730 
731  a = cm_change_around(1, angle);
732 
733  return (a);
734 }
3x3 element cartesian matrix
Definition: matrix.h:96
cmatrix cm_change_around(int axis, double angle)
Rotation matrix for indicated axis.
Definition: matrix.cpp:773
Definition: eci2kep_test.cpp:33
cmatrix cm_change_around_y ( double  angle)

Rotation matrix for Y axis.

Create the DCM that represents a rotation of the given angle around the Y axis.

Parameters
angleAngle of rotation in radians
Returns
Resulting DCM
743 {
744  cmatrix a;
745 
746  a = cm_change_around(2, angle);
747 
748  return (a);
749 }
3x3 element cartesian matrix
Definition: matrix.h:96
cmatrix cm_change_around(int axis, double angle)
Rotation matrix for indicated axis.
Definition: matrix.cpp:773
Definition: eci2kep_test.cpp:33
cmatrix cm_change_around_z ( double  angle)

Rotation matrix for Z axis.

Create the DCM that represents a rotation of the given angle around the Z axis.

Parameters
angleAngle of rotation in radians
Returns
Resulting DCM
758 {
759  cmatrix a;
760 
761  a = cm_change_around(3, angle);
762 
763  return (a);
764 }
3x3 element cartesian matrix
Definition: matrix.h:96
cmatrix cm_change_around(int axis, double angle)
Rotation matrix for indicated axis.
Definition: matrix.cpp:773
Definition: eci2kep_test.cpp:33
cmatrix cm_change_around ( int  axis,
double  angle 
)

Rotation matrix for indicated axis.

Create the DCM that represents a rotation of the given angle around the indicated axis.

Parameters
axisAxis of rotation: 1=X, 2=Y, 3=Z
angleAngle of rotation in radians
Returns
Resulting DCM
774 {
775  cmatrix a = {{1.,0.,0.},{0.,1.,0.},{0.,0.,1.}};
776 
777  switch (axis)
778  {
779  case 1:
780  a.r2.y = a.r3.z = cos(angle);
781  a.r3.y = sin(angle);
782  a.r2.z = -(a.r3.y);
783  break;
784  case 2:
785  a.r1.x = a.r3.z = cos(angle);
786  a.r1.z = sin(angle);
787  a.r3.x = -(a.r1.z);
788  break;
789  case 3:
790  a.r1.x = a.r2.y = cos(angle);
791  a.r2.x = sin(angle);
792  a.r1.y = -(a.r2.x);
793  break;
794  }
795 
796  return (a);
797 }
double y
Y value.
Definition: vector.h:114
3x3 element cartesian matrix
Definition: matrix.h:96
cvector r1
Row 1.
Definition: matrix.h:99
double x
X value.
Definition: vector.h:112
Definition: eci2kep_test.cpp:33
double z
Z value.
Definition: vector.h:116
cvector r3
Row 3.
Definition: matrix.h:103
cvector r2
Row 2.
Definition: matrix.h:101
rmatrix rm_from_cm ( cmatrix  matrix)

rmatrix from cmatrix

Converts 3x3 matrix in cartesian form to row major form.

Parameters
matrixCartesian matrix to convert
Returns
Converted matrix in row major form
805 {
806  rmatrix rm;
807 
808  rm.row[0].col[0] = matrix.r1.x;
809  rm.row[0].col[1] = matrix.r1.y;
810  rm.row[0].col[2] = matrix.r1.z;
811  rm.row[1].col[0] = matrix.r2.x;
812  rm.row[1].col[1] = matrix.r2.y;
813  rm.row[1].col[2] = matrix.r2.z;
814  rm.row[2].col[0] = matrix.r3.x;
815  rm.row[2].col[1] = matrix.r3.y;
816  rm.row[2].col[2] = matrix.r3.z;
817 
818  return (rm);
819 }
double y
Y value.
Definition: vector.h:114
cvector r1
Row 1.
Definition: matrix.h:99
double x
X value.
Definition: vector.h:112
3x3 element generic matrix
Definition: matrix.h:41
double z
Z value.
Definition: vector.h:116
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
cvector r3
Row 3.
Definition: matrix.h:103
cvector r2
Row 2.
Definition: matrix.h:101
rmatrix rm_from_rv ( rvector  vector,
int  direction 
)

rvector to rmatrix.

Convert a row vector to a row ordder rmatrix

Parameters
vectorRow vector to be converted
directionAlignment, column order if 1 (DIRECTION_COLUMN), otherwise row order.
Returns
Single row or column rmatrix
829 {
830  rmatrix answer;
831  uint16_t i;
832 
833  answer = rm_zero();
834  switch (direction)
835  {
836  case DIRECTION_COLUMN:
837  for (i=0; i<3; i++)
838  {
839  answer.row[i].col[0] = vector.col[i];
840  }
841  break;
842  break;
843  default:
844  answer.row[0] = vector;
845  break;
846  }
847 
848  return (answer);
849 }
int i
Definition: rw_test.cpp:37
rmatrix rm_zero()
Zero filled rmatrix.
Definition: matrix.cpp:125
3x3 element generic matrix
Definition: matrix.h:41
#define DIRECTION_COLUMN
Definition: math/constants.h:64
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
rmatrix rm_skew ( rvector  row)

Create skew symmetric rmatrix from rvector.

Generate a new 3x3 skew symmetric row matrix from the provided 3 element vector.

Parameters
rowrow vector to be skewed
Returns
skew symmetric row matrix
858 {
859  rmatrix answer;
860 
861  answer = rm_zero();
862 
863  answer.row[0].col[1] = -row.col[2];
864  answer.row[0].col[2] = row.col[1];
865 
866  answer.row[1].col[0] = row.col[2];
867  answer.row[1].col[2] = -row.col[0];
868 
869  answer.row[2].col[0] = -row.col[1];
870  answer.row[2].col[1] = row.col[0];
871 
872  return (answer);
873 }
rmatrix rm_zero()
Zero filled rmatrix.
Definition: matrix.cpp:125
3x3 element generic matrix
Definition: matrix.h:41
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
rvector rv_unskew ( rmatrix  matrix)

Unskew 3x3 row matrix.

Create the 3 element rvector correponding to a 3x3 skew symmetric rmatrix

Parameters
matrix3x3 row skew symmetric matrix
Returns
rvector representing unskewed elements
881 {
882  rvector answer;
883 
884  answer.col[0] = -matrix.row[1].col[2];
885  answer.col[1] = matrix.row[0].col[2];
886  answer.col[2] = -matrix.row[0].col[1];
887 
888  return (answer);
889 }
3 element generic row vector
Definition: vector.h:53
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
rmatrix rm_inverse ( rmatrix  m)

Inverse of rmatrix.

Inverse of 3x3 rmatrix using algorithm at http://mathworld.wolfram.com/MatrixInverse.html

Parameters
mrmatrix to take inverse of
Returns
Inverse matrix, or
898 {
899  rmatrix wm;
900 
901  wm.row[0].col[0] = m.row[1].col[1]*m.row[2].col[2] - m.row[1].col[2]*m.row[2].col[1];
902  wm.row[0].col[1] = m.row[0].col[2]*m.row[2].col[1] - m.row[0].col[1]*m.row[2].col[2];
903  wm.row[0].col[2] = m.row[0].col[1]*m.row[1].col[2] - m.row[0].col[2]*m.row[1].col[1];
904 
905  wm.row[1].col[0] = m.row[1].col[2]*m.row[2].col[0] - m.row[1].col[0]*m.row[2].col[2];
906  wm.row[1].col[1] = m.row[0].col[0]*m.row[2].col[2] - m.row[0].col[2]*m.row[2].col[0];
907  wm.row[1].col[2] = m.row[0].col[2]*m.row[1].col[0] - m.row[0].col[0]*m.row[1].col[2];
908 
909  wm.row[2].col[0] = m.row[1].col[0]*m.row[2].col[1] - m.row[1].col[1]*m.row[2].col[0];
910  wm.row[2].col[1] = m.row[0].col[1]*m.row[2].col[0] - m.row[0].col[0]*m.row[2].col[1];
911  wm.row[2].col[2] = m.row[0].col[0]*m.row[1].col[1] - m.row[0].col[1]*m.row[1].col[0];
912 
913  return( rm_smult((1/determinant_rm(m)), wm) );
914 }
double determinant_rm(rmatrix m)
Determinant of row column matrix.
Definition: matrix.cpp:942
3x3 element generic matrix
Definition: matrix.h:41
rmatrix rm_smult(double a, rmatrix b)
Scalar rmatrix multiplication.
Definition: matrix.cpp:247
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
rmatrix rm_from_m2 ( matrix2d  matrix)

rmatrix from matrix2d

Converts 3x3 matrix in matrix2d form to row major form.

Parameters
matrixmatrix2d matrix to convert
Returns
Converted matrix in row major form
922 {
923  rmatrix rm;
924  int i, j;
925 
926  for (i=0; i<3; i++)
927  {
928  for (j=0; j<3; j++)
929  {
930  rm.row[i].col[j] = matrix.array[i][j];
931  }
932  }
933 
934  return (rm);
935 }
int i
Definition: rw_test.cpp:37
double array[4][4]
Elements.
Definition: matrix.h:160
3x3 element generic matrix
Definition: matrix.h:41
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
double determinant_rm ( rmatrix  m)

Determinant of row column matrix.

Return the determinant for a 3x3 rmatrix

Parameters
m;;rmatrix to calculate detrminant of
Returns
determinant, otherwise NaN
943 {
944  double result;
945 
946  result = m.row[0].col[0] * (m.row[1].col[1] * m.row[2].col[2] - m.row[2].col[1] * m.row[1].col[2]);
947  result -= m.row[1].col[0] * (m.row[0].col[1] * m.row[2].col[2] - m.row[0].col[2] * m.row[2].col[1]);
948  result += m.row[2].col[0] * (m.row[0].col[1] * m.row[1].col[2] - m.row[0].col[2] * m.row[1].col[1]);
949  return (result);
950 
951 }
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
matrix1d m1_zero ( uint16_t  cols)

Fill 1D matrix with zeros.

Fill the provided vector with zeros, and set its size to the provided number of columns

Parameters
colsNumber of columns
Returns
Zero filled vector
960 {
961  matrix1d answer={{0.,0.,0.,0.},0};
962  uint16_t i;
963 
964  answer.cols = cols;
965  for (i=0; i<cols; i++)
966  {
967  answer.vector[i] = 0.;
968  }
969  return (answer);
970 }
uint16_t cols
Number of elements.
Definition: matrix.h:149
int i
Definition: rw_test.cpp:37
n element row matrix
Definition: matrix.h:144
double vector[4]
Elements.
Definition: matrix.h:147
matrix1d m1_smult ( double  number,
matrix1d  row 
)

Multiply 1D matrix by a scalar.

Multiply each element of indicated 1D matrix by a scalar and return as a new 1D matrix.

Parameters
numberscalar to multiply by
row1D matrix to be multiplied
Returns
pointer to the matrix1d result, otherwise NULL
980 {
981  matrix1d answer={{0.,0.,0.,0.},0};
982  uint16_t i;
983 
984  answer.cols = row.cols;
985  for (i=0; i<row.cols; i++)
986  {
987  answer.vector[i] = row.vector[i] * number;
988  }
989 
990  return (answer);
991 }
uint16_t cols
Number of elements.
Definition: matrix.h:149
int i
Definition: rw_test.cpp:37
n element row matrix
Definition: matrix.h:144
double vector[4]
Elements.
Definition: matrix.h:147
matrix1d m1_add ( matrix1d  row1,
matrix1d  row2 
)

Add one 1D matrix to another.

Add each element in the second matrix to each element in the first matrix, returning the answer as a new matrix.

Parameters
row1matrix to be added to
row2matrix to add
Returns
sum, as newly created 1D matrix
1001 {
1002  matrix1d answer={{0.,0.,0.,0.},0};
1003  uint16_t i;
1004 
1005  answer.cols = row1.cols;
1006  for (i=0; i<row1.cols; i++)
1007  {
1008  answer.vector[i] = row1.vector[i] + row2.vector[i];
1009  }
1010 
1011  return (answer);
1012 }
uint16_t cols
Number of elements.
Definition: matrix.h:149
int i
Definition: rw_test.cpp:37
n element row matrix
Definition: matrix.h:144
double vector[4]
Elements.
Definition: matrix.h:147
matrix1d m1_sub ( matrix1d  row1,
matrix1d  row2 
)

Subtract one 1D matrix from another.

Subtract each element in the second matrix from each element in the first matrix, returning the answer as a new matrix.

Parameters
row1matrix to be subtracted from
row2matrix to subtract
Returns
difference, as newly created 1D matrix
1022 {
1023  matrix1d answer={{0.,0.,0.,0.},0};
1024  uint16_t i;
1025 
1026  answer.cols = row1.cols;
1027  for (i=0; i<row1.cols; i++)
1028  {
1029  answer.vector[i] = row1.vector[i] - row2.vector[i];
1030  }
1031 
1032  return (answer);
1033 }
uint16_t cols
Number of elements.
Definition: matrix.h:149
int i
Definition: rw_test.cpp:37
n element row matrix
Definition: matrix.h:144
double vector[4]
Elements.
Definition: matrix.h:147
matrix1d m1_mmult ( matrix2d  Matrix,
matrix1d  Vector 
)

Multiply matrix1d by matrix2d.

Perform a matrix multiplication of a matrix1d vector as though it were a column matrix.

Parameters
MatrixMatrix to multiply by
VectorVector to multiply
Returns
Vector result
1042 {
1043  matrix1d result={{0.,0.,0.,0.},0};
1044  int i, j;
1045 
1046  if (Matrix.cols == Vector.cols)
1047  {
1048  result.cols = Vector.cols;
1049  for (i=0; i<result.cols; i++)
1050  {
1051  result.vector[i] = 0.;
1052  for (j=0; j<result.cols; j++)
1053  {
1054  result.vector[i] += Vector.vector[j]*Matrix.array[i][j];
1055  }
1056  }
1057  }
1058 
1059  return(result);
1060 }
uint16_t cols
Number of elements.
Definition: matrix.h:149
int i
Definition: rw_test.cpp:37
double array[4][4]
Elements.
Definition: matrix.h:160
n element row matrix
Definition: matrix.h:144
double vector[4]
Elements.
Definition: matrix.h:147
uint16_t cols
Number of elements.
Definition: matrix.h:158
matrix1d m1_cross ( matrix1d  Vector1,
matrix1d  Vector2 
)

matrix1d cross product

Computes the vector cross product (A x B) of two 1d vectors.

Parameters
Vector1vector A
Vector2vector B
Returns
Vector cross product as 1d vector.
1068 {
1069  matrix1d result;
1070 
1071  result.cols = Vector1.cols;
1072  result.vector[0] = Vector1.vector[1]*Vector2.vector[2] - Vector1.vector[2]*Vector2.vector[1];
1073  result.vector[1] = -(Vector1.vector[0]*Vector2.vector[2] - Vector1.vector[2]*Vector2.vector[0]);
1074  result.vector[2] = Vector1.vector[0]*Vector2.vector[1] - Vector1.vector[1]*Vector2.vector[0];
1075 
1076  return(result);
1077 }
uint16_t cols
Number of elements.
Definition: matrix.h:149
n element row matrix
Definition: matrix.h:144
double vector[4]
Elements.
Definition: matrix.h:147
double m1_dot ( matrix1d  a,
matrix1d  b 
)

matrix1d dot product

Computes the vector dot product (A x B) of two 1d vectors.

Parameters
avector A
bvector B
Returns
Vector cross product as 1d vector.
1086 {
1087  double d = 0.;
1088  uint16_t i;
1089 
1090  for(i=0; i<a.cols; i++)
1091  {
1092  d += a.vector[i]*b.vector[i];
1093  }
1094  return (d);
1095 }
uint16_t cols
Number of elements.
Definition: matrix.h:149
int i
Definition: rw_test.cpp:37
double vector[4]
Elements.
Definition: matrix.h:147
matrix2d m2_skew ( matrix1d  row)

Create skew symmetric matrix2d from matrix1d.

Generate a new 3x3 skew symmetric 2D matrix from the provided 3 element vector.

Parameters
rowrow vector to be skewed
Returns
skew symmetric 2D matrix
1104 {
1105  matrix2d answer={0,0,{{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}}};
1106 
1107  if (row.cols == 3)
1108  {
1109  answer.rows = answer.cols = 3;
1110  answer.array[0][0] = answer.array[1][1] = answer.array[2][2] = 0.;
1111 
1112  answer.array[0][1] = -row.vector[2];
1113  answer.array[0][2] = row.vector[1];
1114 
1115  answer.array[1][0] = row.vector[2];
1116  answer.array[1][2] = -row.vector[0];
1117 
1118  answer.array[2][0] = -row.vector[1];
1119  answer.array[2][1] = row.vector[0];
1120  }
1121 
1122  return (answer);
1123 }
uint16_t cols
Number of elements.
Definition: matrix.h:149
double array[4][4]
Elements.
Definition: matrix.h:160
uint16_t rows
Number of rows.
Definition: matrix.h:156
nxm element 2D matrix
Definition: matrix.h:153
double vector[4]
Elements.
Definition: matrix.h:147
uint16_t cols
Number of elements.
Definition: matrix.h:158
matrix2d m2_diag ( matrix1d  row)

Create diagonal matrix2d from matrix1d.

Generate a new nxn diagonal 2D matrix from the provided n element 1D matrix.

Parameters
row1D matrix to be skewed
Returns
result as new 2D matrix, otherwise NULL
1132 {
1133  matrix2d answer={0,0,{{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}}};
1134 
1135  if (row.cols == 3)
1136  {
1137  answer.cols = answer.rows = 3;
1138  answer.array[0][0] = row.vector[0];
1139  answer.array[1][1] = row.vector[1];
1140  answer.array[2][2] = row.vector[2];
1141  }
1142 
1143  return (answer);
1144 }
uint16_t cols
Number of elements.
Definition: matrix.h:149
double array[4][4]
Elements.
Definition: matrix.h:160
uint16_t rows
Number of rows.
Definition: matrix.h:156
nxm element 2D matrix
Definition: matrix.h:153
double vector[4]
Elements.
Definition: matrix.h:147
uint16_t cols
Number of elements.
Definition: matrix.h:158
matrix2d m2_inverse ( matrix2d  m)

Inverse of matrix2d.

Inverse of 2x2 or 3x3 Matrix2D using algorithm at http://mathworld.wolfram.com/MatrixInverse.html

Parameters
mMatrix2d to take inverse of
Returns
Inverse matrix, or
1153 {
1154  matrix2d wm;
1155 
1156  wm.rows = m.rows;
1157  wm.cols = m.cols;
1158 
1159  if (m.rows == 2)
1160  {
1161  wm.array[0][0] = m.array[1][1];
1162  wm.array[1][1] = m.array[0][0];
1163  wm.array[1][0] = -m.array[1][0];
1164  wm.array[0][1] = -m.array[0][1];
1165  }
1166  else
1167  {
1168  wm.array[0][0] = m.array[1][1]*m.array[2][2] - m.array[1][2]*m.array[2][1];
1169  wm.array[0][1] = m.array[0][2]*m.array[2][1] - m.array[0][1]*m.array[2][2];
1170  wm.array[0][2] = m.array[0][1]*m.array[1][2] - m.array[0][2]*m.array[1][1];
1171 
1172  wm.array[1][0] = m.array[1][2]*m.array[2][0] - m.array[1][0]*m.array[2][2];
1173  wm.array[1][1] = m.array[0][0]*m.array[2][2] - m.array[0][2]*m.array[2][0];
1174  wm.array[1][2] = m.array[0][2]*m.array[1][0] - m.array[0][0]*m.array[1][2];
1175 
1176  wm.array[2][0] = m.array[1][0]*m.array[2][1] - m.array[1][1]*m.array[2][0];
1177  wm.array[2][1] = m.array[0][1]*m.array[2][0] - m.array[0][0]*m.array[2][1];
1178  wm.array[2][2] = m.array[0][0]*m.array[1][1] - m.array[0][1]*m.array[1][0];
1179  }
1180 
1181  return( m2_smult((1/m2_determinant(m)), wm) );
1182 }
matrix2d m2_smult(double number, matrix2d matrix)
Multiply 2D matrix by a scalar.
Definition: matrix.cpp:1281
double array[4][4]
Elements.
Definition: matrix.h:160
uint16_t rows
Number of rows.
Definition: matrix.h:156
nxm element 2D matrix
Definition: matrix.h:153
uint16_t cols
Number of elements.
Definition: matrix.h:158
double m2_determinant(matrix2d m)
Determinant of a 2D matrix.
Definition: matrix.cpp:1189
double m2_determinant ( matrix2d  m)

Determinant of a 2D matrix.

Return the determinant for a 2x2 or 3x3 2D Matrix

Parameters
mSquare 2D matrix to calculate detrminant of
Returns
determinant, otherwise NaN
1190 {
1191  double result;
1192 
1193  if (m.rows == 2 && m.cols == 2)
1194  {
1195  result = m.array[0][0]*m.array[1][1] - m.array[0][1]*m.array[1][0];
1196  return (result);
1197  }
1198 
1199  if (m.rows == 3 && m.cols == 3)
1200  {
1201  result = m.array[0][0] * (m.array[1][1] * m.array[2][2] - m.array[1][2] * m.array[2][1]);
1202  result += m.array[0][1] * (m.array[1][2] * m.array[2][0] - m.array[1][0] * m.array[2][2]);
1203  result += m.array[0][2] * (m.array[1][0] * m.array[2][1] - m.array[1][1] * m.array[2][0]);
1204  return (result);
1205  }
1206 
1207  return (NAN);
1208 }
double array[4][4]
Elements.
Definition: matrix.h:160
uint16_t rows
Number of rows.
Definition: matrix.h:156
uint16_t cols
Number of elements.
Definition: matrix.h:158
double m1_norm ( matrix1d  row)

Compute the Euclidean norm of a 1D matrix.

Calculate and return the Euclidean norm of the provided 1D matrix.

Parameters
row1D matrix to be normed
Returns
the calculated norm
1216 {
1217  double norm=0.;
1218 
1219  for (int i=0; i<row.cols; i++)
1220  {
1221  norm += row.vector[i] * row.vector[i];
1222  }
1223 
1224  return (sqrt(norm));
1225 }
uint16_t cols
Number of elements.
Definition: matrix.h:149
int i
Definition: rw_test.cpp:37
double vector[4]
Elements.
Definition: matrix.h:147
matrix2d m2_zero ( uint16_t  rows,
uint16_t  cols 
)

Create 2D zero matrix.

Create a matrix of the requested sized, set to zero.

Parameters
rowsNumber of rows in the matrix
colsNumber of columns in the matrix
Returns
New 2D matrix
1234 {
1235  int i, j;
1236  matrix2d answer={0,0,{{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}}};
1237 
1238  if (rows <= 3 && cols <=3)
1239  {
1240  answer.rows = rows;
1241  answer.cols = cols;
1242  for (i=0; i<rows; i++)
1243  {
1244  for (j=0; j<cols; j++)
1245  {
1246  answer.array[i][j] = 0.;
1247  }
1248  }
1249  }
1250 
1251  return answer;
1252 }
int i
Definition: rw_test.cpp:37
double array[4][4]
Elements.
Definition: matrix.h:160
uint16_t rows
Number of rows.
Definition: matrix.h:156
nxm element 2D matrix
Definition: matrix.h:153
uint16_t cols
Number of elements.
Definition: matrix.h:158
matrix2d m2_eye ( uint16_t  rows)

Create 2D identity matrix.

Create a square identity matrix of the requested size and return as a new matrix.

Parameters
rowsthe size, in rows and columns, of the new matrix
Returns
the new identity matrix, as a 2D matrix
1262 {
1263  matrix2d answer={0,0,{{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}}};
1264 
1265  if (rows <= 3)
1266  {
1267  answer = m2_zero(rows,rows);
1268  answer.array[0][0] = answer.array[1][1] = answer.array[2][2] = 1.;
1269  }
1270 
1271  return (answer);
1272 }
double array[4][4]
Elements.
Definition: matrix.h:160
nxm element 2D matrix
Definition: matrix.h:153
matrix2d m2_zero(uint16_t rows, uint16_t cols)
Create 2D zero matrix.
Definition: matrix.cpp:1233
matrix2d m2_smult ( double  number,
matrix2d  matrix 
)

Multiply 2D matrix by a scalar.

Multiply each element of indicated 2D matrix by a scalar and return as a new 2D matrix.

Parameters
numberscalar to multiply by
matrix2D matrix to be multiplied
Returns
2D matrix result
1282 {
1283  matrix2d answer={0,0,{{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}}};
1284  uint16_t i, j;
1285 
1286  answer = matrix;
1287  for (i=0; i<matrix.rows; i++)
1288  {
1289  for (j=0; j<matrix.cols; j++)
1290  {
1291  answer.array[i][j] *= number;
1292  }
1293  }
1294  return (answer);
1295 }
int i
Definition: rw_test.cpp:37
double array[4][4]
Elements.
Definition: matrix.h:160
uint16_t rows
Number of rows.
Definition: matrix.h:156
nxm element 2D matrix
Definition: matrix.h:153
uint16_t cols
Number of elements.
Definition: matrix.h:158
matrix2d m2_add ( matrix2d  matrix1,
matrix2d  matrix2 
)

Add one matrix2d to another.

Add each element in the second matrix to each element in the first matrix, returning the answer as a new matrix.

Parameters
matrix1matrix to be added from
matrix2matrix to add
Returns
sum, as newly created 2D matrix
1305 {
1306  matrix2d answer={0,0,{{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}}};
1307  uint16_t i, j;
1308 
1309  if (matrix1.rows == matrix2.rows && matrix1.cols == matrix2.cols)
1310  {
1311  answer.rows = matrix1.rows;
1312  answer.cols = matrix1.cols;
1313  for (i=0; i<matrix1.rows; i++)
1314  {
1315  for (j=0; j<matrix1.cols; j++)
1316  {
1317  answer.array[i][j] = matrix1.array[i][j] + matrix2.array[i][j];
1318  }
1319  }
1320  }
1321 
1322  return (answer);
1323 }
int i
Definition: rw_test.cpp:37
double array[4][4]
Elements.
Definition: matrix.h:160
uint16_t rows
Number of rows.
Definition: matrix.h:156
nxm element 2D matrix
Definition: matrix.h:153
uint16_t cols
Number of elements.
Definition: matrix.h:158
matrix2d m2_sub ( matrix2d  matrix1,
matrix2d  matrix2 
)

Subtract one matrix2d from another.

Subtract each element in the second matrix from each element in the first matrix, returning the answer as a new matrix.

Parameters
matrix1matrix to be subtracted from
matrix2matrix to subtract
Returns
difference, as newly created 2D matrix
1333 {
1334  matrix2d answer={0,0,{{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}}};
1335  uint16_t i, j;
1336 
1337  if (matrix1.rows == matrix2.rows && matrix1.cols == matrix2.cols)
1338  {
1339  answer.rows = matrix1.rows;
1340  answer.cols = matrix1.cols;
1341  for (i=0; i<matrix1.rows; i++)
1342  {
1343  for (j=0; j<matrix1.cols; j++)
1344  {
1345  answer.array[i][j] = matrix1.array[i][j] - matrix2.array[i][j];
1346  }
1347  }
1348  }
1349 
1350  return (answer);
1351 }
int i
Definition: rw_test.cpp:37
double array[4][4]
Elements.
Definition: matrix.h:160
uint16_t rows
Number of rows.
Definition: matrix.h:156
nxm element 2D matrix
Definition: matrix.h:153
uint16_t cols
Number of elements.
Definition: matrix.h:158
matrix2d m2_transpose ( matrix2d  matrix)

Return transpose of a 2D matrix.

Generate a new matrix2d that is is the transpose of the provided matrix2d

Parameters
matrixmatrix2d to be transposed
Returns
Transpose, as new matrix2d
1359 {
1360  matrix2d answer={0,0,{{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}}};
1361  uint16_t i, j;
1362 
1363  if (matrix.rows <=3 && matrix.cols <= 3)
1364  {
1365  answer.rows = matrix.cols;
1366  answer.cols = matrix.rows;
1367  for (i=0; i<matrix.rows; i++)
1368  {
1369  for (j=0; j<matrix.cols; j++)
1370  {
1371  answer.array[j][i] = matrix.array[i][j];
1372  }
1373  }
1374  }
1375 
1376  return (answer);
1377 }
int i
Definition: rw_test.cpp:37
double array[4][4]
Elements.
Definition: matrix.h:160
uint16_t rows
Number of rows.
Definition: matrix.h:156
nxm element 2D matrix
Definition: matrix.h:153
uint16_t cols
Number of elements.
Definition: matrix.h:158
matrix1d m2_unskew ( matrix2d  matrix)

Unskew 3x3 2D matrix.

Create the 3 element 1D matrix correponding to a 3x3 skew symmetric matrix

Parameters
matrix3x3 2D skew symmetric matrix
Returns
matrix1d representing unskewed elements
1385 {
1386  matrix1d answer={{0.,0.,0.,0.},0};
1387 
1388  if (matrix.rows == 3 && matrix.cols == 3)
1389  {
1390  answer.cols = 3;
1391  answer.vector[0] = -matrix.array[1][2];
1392  answer.vector[1] = matrix.array[0][2];
1393  answer.vector[2] = -matrix.array[0][1];
1394  }
1395 
1396  return (answer);
1397 }
uint16_t cols
Number of elements.
Definition: matrix.h:149
double array[4][4]
Elements.
Definition: matrix.h:160
uint16_t rows
Number of rows.
Definition: matrix.h:156
n element row matrix
Definition: matrix.h:144
double vector[4]
Elements.
Definition: matrix.h:147
uint16_t cols
Number of elements.
Definition: matrix.h:158
double m2_trace ( matrix2d  matrix)

Calculate the trace of a 2D matrix.

Add the diagonal elements of an nxn 2D matrix.

Parameters
matrixnxn matrix to take the elements from
Returns
Sum of diagonals, otherwise 0.0
1405 {
1406  double answer=0.;
1407  uint16_t i;
1408 
1409  for (i=0; i<matrix.rows; i++)
1410  {
1411  answer += matrix.array[i][i];
1412  }
1413 
1414  return (answer);
1415 }
int i
Definition: rw_test.cpp:37
double array[4][4]
Elements.
Definition: matrix.h:160
uint16_t rows
Number of rows.
Definition: matrix.h:156
matrix2d m2_mmult ( matrix2d  matrix1,
matrix2d  matrix2 
)

Matrix product.

Calculate the matrix product of two 2D matrices and return a new 2D matrix

Parameters
matrix1mxn matrix
matrix2nxm matrix
Returns
new matrix, of dimension rows x cols, otherwise NULL
1424 {
1425  matrix2d answer={0,0,{{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}}};
1426  uint16_t i, j, k;
1427 
1428  if (matrix1.cols == matrix2.rows && matrix1.rows <= 3 && matrix1.cols <=3 && matrix2.rows <= 3 && matrix2.cols <=3)
1429  {
1430  answer.rows = matrix1.rows;
1431  answer.cols = matrix2.cols;
1432  for (i=0; i<matrix1.rows; i++)
1433  {
1434  for (j=0; j<matrix2.cols; j++)
1435  {
1436  answer.array[i][j] = 0.;
1437  for (k=0; k<matrix1.cols; k++)
1438  {
1439  answer.array[i][j] += matrix1.array[i][k] * matrix2.array[k][j];
1440  }
1441  }
1442  }
1443  }
1444 
1445  return (answer);
1446 }
int i
Definition: rw_test.cpp:37
double array[4][4]
Elements.
Definition: matrix.h:160
uint16_t rows
Number of rows.
Definition: matrix.h:156
nxm element 2D matrix
Definition: matrix.h:153
uint16_t cols
Number of elements.
Definition: matrix.h:158
matrix2d m2_from_rm ( rmatrix  matrix)

rmatrix from rmatrix

Converts 3x3 matrix in row order form to cartesian form.

Parameters
matrixRow major matrix to convert
Returns
Converted matrix in row order form
1454 {
1455  matrix2d m2;
1456  int i, j;
1457 
1458  for (i=0; i<3; i++)
1459  {
1460  for (j=0; j<3; j++)
1461  {
1462  m2.array[i][j] = matrix.row[i].col[j];
1463  }
1464  }
1465 
1466  return (m2);
1467 }
int i
Definition: rw_test.cpp:37
double array[4][4]
Elements.
Definition: matrix.h:160
nxm element 2D matrix
Definition: matrix.h:153
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
matrix2d cm3x3_to_m2 ( cmatrix  matrix)

Convert a 3x3 Cartesian matrix to a matrix2d.

Parameters
matrixcmatrix to convert
Returns
matrix2d
1475 {
1476  matrix2d answer={0,0,{{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}}};
1477 
1478  answer.rows = 3;
1479  answer.cols = 3;
1480  answer.array[0][0] = matrix.r1.x;
1481  answer.array[0][1] = matrix.r1.y;
1482  answer.array[0][2] = matrix.r1.z;
1483  answer.array[1][0] = matrix.r2.x;
1484  answer.array[1][1] = matrix.r2.y;
1485  answer.array[1][2] = matrix.r2.z;
1486  answer.array[2][0] = matrix.r3.x;
1487  answer.array[2][1] = matrix.r3.y;
1488  answer.array[2][2] = matrix.r3.z;
1489 
1490  return (answer);
1491 }
double y
Y value.
Definition: vector.h:114
cvector r1
Row 1.
Definition: matrix.h:99
double x
X value.
Definition: vector.h:112
double array[4][4]
Elements.
Definition: matrix.h:160
uint16_t rows
Number of rows.
Definition: matrix.h:156
nxm element 2D matrix
Definition: matrix.h:153
double z
Z value.
Definition: vector.h:116
uint16_t cols
Number of elements.
Definition: matrix.h:158
cvector r3
Row 3.
Definition: matrix.h:103
cvector r2
Row 2.
Definition: matrix.h:101
matrix1d cv_to_m1 ( cvector  vector)

Convert a Cartesian vector to a matrix1d.

Parameters
vectorcvector to convert
Returns
Row vector
1499 {
1500  matrix1d answer={{0.,0.,0.,0.},0};
1501 
1502  answer.cols = 3;
1503  answer.vector[0] = vector.x;
1504  answer.vector[1] = vector.y;
1505  answer.vector[2] = vector.z;
1506 
1507  return (answer);
1508 }
double y
Y value.
Definition: vector.h:114
uint16_t cols
Number of elements.
Definition: matrix.h:149
double x
X value.
Definition: vector.h:112
n element row matrix
Definition: matrix.h:144
double z
Z value.
Definition: vector.h:116
double vector[4]
Elements.
Definition: matrix.h:147
matrix2d cv_to_m2 ( cvector  vector,
int  direction 
)

Convert a Cartesian vector to either a row or column matrix2d.

Parameters
vectorcvector to convert
directioncolumn if 1 (DIRECTION_COLUMN), otherwise row
Returns
Row or column matrix
1517 {
1518  matrix2d answer={0,0,{{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}}};
1519 
1520  if (direction == DIRECTION_COLUMN)
1521  {
1522  answer.rows = 3;
1523  answer.cols = 1;
1524  answer.array[0][0] = vector.x;
1525  answer.array[1][0] = vector.y;
1526  answer.array[2][0] = vector.z;
1527  }
1528  else
1529  {
1530  answer.rows = 1;
1531  answer.cols = 3;
1532  answer.array[0][0] = vector.x;
1533  answer.array[0][1] = vector.y;
1534  answer.array[0][2] = vector.z;
1535  }
1536  return (answer);
1537 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
double array[4][4]
Elements.
Definition: matrix.h:160
uint16_t rows
Number of rows.
Definition: matrix.h:156
nxm element 2D matrix
Definition: matrix.h:153
#define DIRECTION_COLUMN
Definition: math/constants.h:64
double z
Z value.
Definition: vector.h:116
uint16_t cols
Number of elements.
Definition: matrix.h:158
matrix2d m1_to_m2 ( matrix1d  vector,
int  direction 
)

Matrix1d to matrix2d.

Convert a matrix1d vector to a row element matrix2d

Parameters
vectorMatrix1d to be converted
directionAlignment, column order if 1 (DIRECTION_COLUMN), otherwise row order.
Returns
Single row or column matrix2d
1547 {
1548  matrix2d answer={0,0,{{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}}};
1549  uint16_t i;
1550 
1551  switch (direction)
1552  {
1553  case DIRECTION_COLUMN:
1554  answer.cols = 1;
1555  answer.rows = vector.cols;
1556  for (i=0; i<vector.cols; i++)
1557  {
1558  answer.array[i][0] = vector.vector[i];
1559  }
1560  break;
1561  break;
1562  default:
1563  answer.rows = 1;
1564  answer.cols = vector.cols;
1565  for (i=0; i<vector.cols; i++)
1566  {
1567  answer.array[0][i] = vector.vector[i];
1568  }
1569  break;
1570  }
1571 
1572  return (answer);
1573 }
uint16_t cols
Number of elements.
Definition: matrix.h:149
int i
Definition: rw_test.cpp:37
double array[4][4]
Elements.
Definition: matrix.h:160
uint16_t rows
Number of rows.
Definition: matrix.h:156
nxm element 2D matrix
Definition: matrix.h:153
#define DIRECTION_COLUMN
Definition: math/constants.h:64
double vector[4]
Elements.
Definition: matrix.h:147
uint16_t cols
Number of elements.
Definition: matrix.h:158
matrix1d m2_eig2x2 ( matrix2d  matrix)

Eigen values of a 2x2 square matrix.

Calculate the 2 element 1D matrix that is the set of eigenvalues of a 2x2 2D matrix

Parameters
matrix2x2 2D matrix
Returns
2 element 1D matrix of eigenvalues, otherwise NULL
1581 {
1582  matrix1d answer={{0.,0.,0.,0.},0};
1583  double trm, detm;
1584  double trm2, detm4;
1585 
1586  if (matrix.cols == 2 && matrix.rows == 2)
1587  {
1588  trm = matrix.array[0][0] + matrix.array[1][1];
1589  detm = matrix.array[0][0]*matrix.array[1][1] - matrix.array[0][1]*matrix.array[1][0];
1590 
1591  if ((trm2=trm*trm) < (detm4=4*detm))
1592  return (answer);
1593 
1594  answer.cols = 2;
1595  answer.vector[0] = (trm + sqrt(trm2-detm4))/2.;
1596  answer.vector[1] = (trm - sqrt(trm2-detm4))/2.;
1597  }
1598 
1599  return (answer);
1600 }
uint16_t cols
Number of elements.
Definition: matrix.h:149
double array[4][4]
Elements.
Definition: matrix.h:160
uint16_t rows
Number of rows.
Definition: matrix.h:156
n element row matrix
Definition: matrix.h:144
double vector[4]
Elements.
Definition: matrix.h:147
uint16_t cols
Number of elements.
Definition: matrix.h:158
double m2_snorm2x2 ( matrix2d  matrix)

Spectral norm of a 2x2 matrix.

Calculate the spectral norm of the provided 2x2 matrix.

Parameters
matrix2x2 matrix
Returns
spectral norm, otherwise NaN
1608 {
1609  double answer=NAN;
1610  matrix1d re;
1611 
1612  if (matrix.cols == 2 && matrix.rows == 2)
1613  {
1614  re = m2_eig2x2(m2_mmult(m2_transpose(matrix),matrix));
1615 
1616  if (re.cols == 0)
1617  return (NAN);
1618 
1619  if (re.vector[0] > re.vector[1])
1620  answer = sqrt(re.vector[0]);
1621  else
1622  answer = sqrt(re.vector[1]);
1623  }
1624 
1625  return (answer);
1626 }
uint16_t cols
Number of elements.
Definition: matrix.h:149
static double re
Definition: nrlmsise-00.cpp:67
matrix2d m2_transpose(matrix2d matrix)
Return transpose of a 2D matrix.
Definition: matrix.cpp:1358
uint16_t rows
Number of rows.
Definition: matrix.h:156
matrix1d m2_eig2x2(matrix2d matrix)
Eigen values of a 2x2 square matrix.
Definition: matrix.cpp:1580
n element row matrix
Definition: matrix.h:144
double vector[4]
Elements.
Definition: matrix.h:147
uint16_t cols
Number of elements.
Definition: matrix.h:158
matrix2d m2_mmult(matrix2d matrix1, matrix2d matrix2)
Matrix product.
Definition: matrix.cpp:1423
cmatrix cm_quaternion2dcm ( quaternion  q)

Quaternion to Direction Cosine Matrix This function expects a quaternion that represents the coordinate frame transformation not the rotation. If the quaternion represents the rotation from the inertial reference frame into the frame of the sensor/body then this DCM will represent the rotation from the sensor body frame (B) to the inertial frame (I) TODO: later this should be changed to be more consistent

Convert supplied quaternion to an equivalent direction cosine matrix

Parameters
qquaternion
Returns
direction cosine matrix
52 {
53  cmatrix m;
54  double yy, xx, zz, xy, xz, xw, yz, yw, zw;
55 
56  normalize_q(&q);
57 
58  xx = 2. * q.d.x;
59  xy = xx * q.d.y;
60  xz = xx * q.d.z;
61  xw = xx * q.w;
62  xx *= q.d.x;
63 
64  yy = 2. * q.d.y;
65  yz = yy * q.d.z;
66  yw = yy * q.w;
67  yy *= q.d.y;
68 
69  zz = 2. * q.d.z;
70  zw = zz * q.w;
71  zz *= q.d.z;
72 
73  // first row
74  m.r1.x = 1. - yy - zz; // 1 - 2(qy^2 + qz^2) == 1 - 2(q2^ + q3^2)
75  m.r1.y = xy - zw; // 2*(qx*qy - qw*qz)
76  m.r1.z = xz + yw; // 2*(qx*qz + qw*qy)
77 
78  // second row
79  m.r2.x = xy + zw; // 2*(qx*qy + qw*qz) == 2(q1q2 + q0q3)
80  m.r2.y = 1. - xx - zz; // 1 - 2*(qz^2 + qx^2) == 1 - 2(q3^2 + q1^2)
81  m.r2.z = yz - xw; // 2*(qy*qz - qw*qx) == 2(q2q3 - q0q1)
82 
83  // third row
84  m.r3.x = xz - yw; // 2(qx*qz - qw*qy) == 2(q1q3 - q0q2)
85  m.r3.y = yz + xw; // 2(qy*qz + qw*qx) == 2(q2q3 + q0q1)
86  m.r3.z = 1. - xx - yy; // 1 - 2(qx^2 + qy^2) == 1 - 2(q1^2 + q2^2)
87 
88  return (m);
89 }
double y
Y value.
Definition: vector.h:114
3x3 element cartesian matrix
Definition: matrix.h:96
cvector d
Orientation.
Definition: vector.h:405
cvector r1
Row 1.
Definition: matrix.h:99
double x
X value.
Definition: vector.h:112
void normalize_q(quaternion *q)
Definition: vector.cpp:961
double w
Rotation.
Definition: vector.h:407
double z
Z value.
Definition: vector.h:116
cvector r3
Row 3.
Definition: matrix.h:103
cvector r2
Row 2.
Definition: matrix.h:101
quaternion q_dcm2quaternion_cm ( cmatrix  dcm)

Direction Cosine Matrix to Quaternion.

Convert the given DCM to an equivalent quaternion

Parameters
dcmdirection cosine matrix
Returns
q quaternion
101 {
102  quaternion q;
103  double t, tr;
104 
105  // removed trace_cm(dcm) to clear the dependencies from mathlib
106  // computing the trace in place
107  tr = dcm.r1.x + dcm.r2.y + dcm.r3.z; //trace_cm(dcm)
108 
109  // TODO: explain the different cases
110  if (tr > 0.)
111  {
112  // TODO: simplify
113  // q.w = 0.5 * sqrt(1 + tr)
114  t = .5 / sqrt(1.+tr);
115  q.w = .25 / t;
116 
117  // TODO: check because Kuipers does the inverse here
118  q.d.x = t*(dcm.r3.y - dcm.r2.z);
119  q.d.y = t*(dcm.r1.z - dcm.r3.x);
120  q.d.z = t*(dcm.r2.x - dcm.r1.y);
121  }
122  else
123  {
124  if (dcm.r1.x > dcm.r2.y && dcm.r1.x > dcm.r3.z)
125  {
126  t = 2. * sqrt(1. + dcm.r1.x - dcm.r2.y - dcm.r3.z);
127  q.w = (dcm.r3.y - dcm.r2.z) / t;
128  q.d.x = .25 * t;
129  q.d.y = (dcm.r1.y + dcm.r2.x) / t;
130  q.d.z = (dcm.r1.z + dcm.r3.x) / t;
131  }
132  else
133  {
134  if (dcm.r2.x > dcm.r3.z)
135  {
136  t = 2. * sqrt(1. + dcm.r2.y - dcm.r1.x - dcm.r3.z);
137  q.w = (dcm.r1.z - dcm.r3.x) / t;
138  q.d.x = (dcm.r1.y + dcm.r2.x) / t;
139  q.d.y = .25 * t;
140  q.d.z = (dcm.r2.z + dcm.r3.y) / t;
141  }
142  else
143  {
144  t = 2. * sqrt(1. + dcm.r3.z - dcm.r1.x - dcm.r2.y);
145  q.w = (dcm.r2.x - dcm.r1.y) / t;
146  q.d.x = (dcm.r1.z + dcm.r3.x) / t;
147  q.d.y = (dcm.r2.z + dcm.r3.y) / t;
148  q.d.z = .25 * t;
149  }
150  }
151  }
152 
153  normalize_q(&q);
154  return(q);
155 }
double y
Y value.
Definition: vector.h:114
cvector d
Orientation.
Definition: vector.h:405
Quaternion, scalar last, using x, y, z.
Definition: vector.h:402
cvector r1
Row 1.
Definition: matrix.h:99
double x
X value.
Definition: vector.h:112
void normalize_q(quaternion *q)
Definition: vector.cpp:961
double w
Rotation.
Definition: vector.h:407
double z
Z value.
Definition: vector.h:116
cvector r3
Row 3.
Definition: matrix.h:103
cvector r2
Row 2.
Definition: matrix.h:101
rmatrix rm_change_between_rv ( rvector  from,
rvector  to 
)

Create rotation matrix from 2 row vectors.

Generate the direction cosine matrix that represents a rotation from one row order vector to a second row order vector.

Parameters
frominitial row order vector
tofinal row order vector
Returns
direction cosine matrix, in row matrix form, that can be used to rotate points
165 {
166  rmatrix m;
167 
169 
170  return (m);
171 }
rmatrix rm_quaternion2dcm(quaternion q)
Quaternion to row matrix Direction Cosine Matrix.
Definition: mathlib.cpp:219
quaternion q_drotate_between_rv(rvector from, rvector to)
Create rotation quaternion from 2 row vectors.
Definition: mathlib.cpp:81
3x3 element generic matrix
Definition: matrix.h:41
quaternion q_dcm2quaternion_rm ( rmatrix  m)

Row matrix DCM to Quaternion.

Convert Direction Cosine Matrix in row matrix form to a Quaternion.

Parameters
mDirection Cosine Matrix in rmatrix form
Returns
Quaternion representing DCM.
179 {
180  quaternion q={{0.,0.,0.},0.};
181  double t, tr;
182 
183  if ((tr=trace_rm(m)) > 0.)
184  {
185  t = .5 / sqrt(1.+tr);
186  q.w = .25 / t;
187  q.d.x = t*(m.row[2].col[1] - m.row[1].col[2]);
188  q.d.y = t*(m.row[0].col[2] - m.row[2].col[0]);
189  q.d.z = t*(m.row[1].col[0] - m.row[0].col[1]);
190  }
191  else
192  {
193  if (m.row[0].col[0] > m.row[1].col[1] && m.row[0].col[0] > m.row[2].col[2])
194  {
195  t = 2. * sqrt(1. + m.row[0].col[0] - m.row[1].col[1] - m.row[2].col[2]);
196  q.w = (m.row[2].col[1] - m.row[1].col[2]) / t;
197  q.d.x = .25 * t;
198  q.d.y = (m.row[0].col[1] + m.row[1].col[0]) / t;
199  q.d.z = (m.row[0].col[2] + m.row[2].col[0]) / t;
200  }
201  else
202  {
203  if (m.row[1].col[0] > m.row[2].col[2])
204  {
205  t = 2. * sqrt(1. + m.row[1].col[1] - m.row[0].col[0] - m.row[2].col[2]);
206  q.w = (m.row[0].col[2] - m.row[2].col[0]) / t;
207  q.d.x = (m.row[0].col[1] + m.row[1].col[0]) / t;
208  q.d.y = .25 * t;
209  q.d.z = (m.row[1].col[2] + m.row[2].col[1]) / t;
210  }
211  else
212  {
213  t = 2. * sqrt(1. + m.row[2].col[2] - m.row[0].col[0] - m.row[1].col[1]);
214  q.w = (m.row[1].col[0] - m.row[0].col[1]) / t;
215  q.d.x = (m.row[0].col[2] + m.row[2].col[0]) / t;
216  q.d.y = (m.row[1].col[2] + m.row[2].col[1]) / t;
217  q.d.z = .25 * t;
218  }
219  }
220  }
221 
222  normalize_q(&q);
223  return(q);
224 }
double y
Y value.
Definition: vector.h:114
double trace_rm(rmatrix mat)
rmatrix Trace Calculates the trace of the supplied rmatrix.
Definition: matrix.cpp:154
cvector d
Orientation.
Definition: vector.h:405
Quaternion, scalar last, using x, y, z.
Definition: vector.h:402
double x
X value.
Definition: vector.h:112
void normalize_q(quaternion *q)
Definition: vector.cpp:961
double w
Rotation.
Definition: vector.h:407
double z
Z value.
Definition: vector.h:116
double col[3]
Definition: vector.h:55
rvector row[3]
Definition: matrix.h:43
double sep_rv ( rvector  v1,
rvector  v2 
)

Angular separation between row vectors.

Calculates the separation angle between two row order vectors, in radians.

Parameters
v1the first vector, in rvector format
v2the second vector, in rvector format
Returns
The separation angle in radians as a double precision
43 {
44  rvector dv;
45 
46  normalize_rv(v1);
47  normalize_rv(v2);
48 
49  dv = rv_sub(v2,v1);
50  double diff = length_rv(dv);
51 
52  double sepangle = 2. * atan2(diff/2.,sqrt(fmax(0.,1.-diff*diff/4.)));
53  return (sepangle);
54 }
3 element generic row vector
Definition: vector.h:53
double length_rv(rvector v)
Length of row vector.
Definition: vector.cpp:748
void normalize_rv(rvector &v)
Normalize row order vector in place.
Definition: vector.cpp:222
rvector rv_sub(rvector a, rvector b)
Subtract two vectors.
Definition: vector.cpp:315
svector s_convert ( rvector  from)

Convert rvector to svector.

Convert vector in cartesian coordinates to vector in spherical coordinates.

Parameters
fromVector in cartesian coordinates to be converted.
Returns
Vector in spherical coordinates.
67 {
68  svector result;
69 
70  double minir2 = from.col[0] * from.col[0] + from.col[1] * from.col[1];
71  double r2 = minir2 + from.col[2] * from.col[2];
72  result.r = sqrt(r2);
73 
74  double sp = from.col[2] / result.r;
75  result.phi = asin(sp);
76  result.lambda = atan2(from.col[1], from.col[0]);
77 
78  return result;
79 }
double lambda
E/W in radians.
Definition: vector.h:172
double phi
N/S in radians.
Definition: vector.h:170
3 element spherical vector
Definition: vector.h:167
double col[3]
Definition: vector.h:55
double r
Radius in meters.
Definition: vector.h:174
rvector rv_convert ( svector  from)

Convert svector to rvector.

Convert vector in spherical coordinates to vector in cartesian coordinates.

Parameters
fromVector in spherical coordinates to be converted.
Returns
Vector in cartesian coordinates.
87 {
88  rvector result;
89 
90  double sp = sin(from.phi);
91  double cp = cos(from.phi);
92  double sl = sin(from.lambda);
93  double cl = cos(from.lambda);
94  double cpr = cp * from.r;
95 
96  result.col[0] = cpr * cl;
97  result.col[1] = cpr * sl;
98  result.col[2] = from.r * sp;
99 
100  return result;
101 }
double lambda
E/W in radians.
Definition: vector.h:172
3 element generic row vector
Definition: vector.h:53
double phi
N/S in radians.
Definition: vector.h:170
double col[3]
Definition: vector.h:55
double r
Radius in meters.
Definition: vector.h:174
rvector rv_zero ( )

Zero row order vector.

Creates a zero length row order vector.

Returns
a rvector of zero length
108 {
109  rvector v;
110 
111  return (v);
112 }
3 element generic row vector
Definition: vector.h:53
rvector rv_unitx ( double  scale)

Scaled x row vector.

Creates a row order vector with the X value set to scale.

Returns
the vector
119 {
120  rvector v={1.,0.,0.};
121  v.col[0] *= scale;
122  return (v);
123 }
3 element generic row vector
Definition: vector.h:53
double col[3]
Definition: vector.h:55
rvector rv_unity ( double  scale)

Scaled y row vector.

Creates a row order vector with the Y value set to scale.

Returns
the vector
130 {
131  rvector v={0.,1.,0.};
132  v.col[1] *= scale;
133  return (v);
134 }
3 element generic row vector
Definition: vector.h:53
double col[3]
Definition: vector.h:55
rvector rv_unitz ( double  scale)

Scaled z row vector.

Creates a row order vector with the Z value set to scale.

Returns
the vector
141 {
142  rvector v={0.,0.,1.};
143  v.col[2] *= scale;
144  return (v);
145 }
3 element generic row vector
Definition: vector.h:53
double col[3]
Definition: vector.h:55
rvector rv_one ( )

Row vector of ones.

Creates a row order vector with all values set to one.

Returns
a rvector with each dimension set to 1.
152 {
153  rvector v={1.,1.,1.};
154 
155  return (v);
156 }
3 element generic row vector
Definition: vector.h:53
rvector rv_one ( double  x,
double  y,
double  z 
)

Row vector of values.

Creates a row order vector with values set to x, y, and z.

Returns
a rvector with dimensions set to provided values.
163 {
164  rvector v={1.,1.,1.};
165  v.col[0] *= x;
166  v.col[1] *= y;
167  v.col[2] *= z;
168 
169  return (v);
170 }
y
Definition: inputfile.py:6
3 element generic row vector
Definition: vector.h:53
x
Definition: inputfile.py:6
double col[3]
Definition: vector.h:55
rvector rv_shortest ( rvector  v)

Shortest vector.

Creates a row order vector of unit length in the direction of the shortest element of the vector provided.

Parameters
vThe rvector from which to draw the shortest element.
Returns
The unit length rvector in the proper direction.
178 {
179  rvector a={1.,0.,0.};
180 
181  for (int i=1; i<3; i++)
182  {
183  if (fabs(v.col[i]) < fabs(v.col[i-1]))
184  {
185  a.col[i-1] = 0.;
186  a.col[i] = 1.;
187  }
188  }
189  return (a);
190 }
3 element generic row vector
Definition: vector.h:53
int i
Definition: rw_test.cpp:37
Definition: eci2kep_test.cpp:33
double col[3]
Definition: vector.h:55
rvector rv_shortest2 ( rvector  v)
193 {
194  rvector rx = {1,0,0};
195  rvector ry = {0,1,0};
196  rvector rz = {0,0,1};
197 
198  if (v.col[0]<=v.col[1] && v.col[0]<=v.col[2])
199  return (rx);
200  else if (v.col[1]<=v.col[0] && v.col[1]<=v.col[2])
201  return (ry);
202  else
203  return (rz);
204 }
3 element generic row vector
Definition: vector.h:53
double col[3]
Definition: vector.h:55
rvector rv_normal ( rvector  v)

Normalize row order vector.

Returns a normalized version of the requested row order vector.

Parameters
vthe rvector to be normalized
Returns
the normalized version of the vector as rvector
213 {
214  normalize_rv(v);
215  return (v);
216 }
void normalize_rv(rvector &v)
Normalize row order vector in place.
Definition: vector.cpp:222
void normalize_rv ( rvector v)

Normalize row order vector in place.

Normalizes requested row order vector.

Parameters
va pointer to the rvector to be normalized
223 {
224  double mag;
225 
226  mag = length_rv(v);
227 
228  // if the current length is not zero (or already one)
229  if (fabs(mag - (double)0.) > D_SMALL && fabs(mag - (double)1.) > D_SMALL)
230  {
231  v.col[0] /= mag;
232  v.col[1] /= mag;
233  v.col[2] /= mag;
234  }
235 }
double length_rv(rvector v)
Length of row vector.
Definition: vector.cpp:748
const double D_SMALL
Definition: math/constants.h:40
double col[3]
Definition: vector.h:55
double normVector3 ( double  x,
double  y,
double  z 
)

basic function to compute the L2 norm of a 3d generic vector with separate entries

240 {
241  return sqrt(x*x + y*y + z*z);
242 }
y
Definition: inputfile.py:6
x
Definition: inputfile.py:6
void normalizeVector3 ( double &  x,
double &  y,
double &  z 
)

basic function to normalize any 3d vector with separate entries

247 {
248  double norm = normVector3(x,y,z);
249 
250  // if the current norm is not zero (or already one)
251  if (fabs(norm - (double)0.) > D_SMALL && fabs(norm - (double)1.) > D_SMALL)
252  {
253  x /= norm;
254  y /= norm;
255  z /= norm;
256  }
257 }
y
Definition: inputfile.py:6
double normVector3(double x, double y, double z)
basic function to compute the L2 norm of a 3d generic vector with separate entries ...
Definition: vector.cpp:239
x
Definition: inputfile.py:6
const double D_SMALL
Definition: math/constants.h:40
rvector rv_smult ( double  a,
rvector  b 
)

Multiply row vector by scalar.

Multiply a row vector by a double precision scalar.

Parameters
bvector to be tranformed, in rvector form
adouble precision scalar to multiply by
Returns
the transformed vector, in rvector form
267 {
268  rvector c;
269 
270  c.col[0] = a * b.col[0];
271  c.col[1] = a * b.col[1];
272  c.col[2] = a * b.col[2];
273  return (c);
274 }
3 element generic row vector
Definition: vector.h:53
Definition: eci2kep_test.cpp:33
double col[3]
Definition: vector.h:55
rvector rv_sadd ( double  a,
rvector  b 
)

Add scalar to each element of vector.

Add a double precision scalar to each element of a 3 element vector.

Parameters
bvector to be tranformed, in rvector form
adouble precision scalar to add
Returns
the transformed vector, in rvector form
284 {
285  rvector c;
286 
287  c.col[0] = a + b.col[0];
288  c.col[1] = a + b.col[1];
289  c.col[2] = a + b.col[2];
290  return (c);
291 }
3 element generic row vector
Definition: vector.h:53
Definition: eci2kep_test.cpp:33
double col[3]
Definition: vector.h:55
rvector rv_add ( rvector  a,
rvector  b 
)

Add two row vectors.

Add two vectors in rvector form, returning an rvector.

Parameters
afirst vector to be added, in rvector form
bsecond vector to be added, in rvector form
Returns
the transformed vector, in rvector form
300 {
301  rvector c;
302 
303  c.col[0] = a.col[0] + b.col[0];
304  c.col[1] = a.col[1] + b.col[1];
305  c.col[2] = a.col[2] + b.col[2];
306  return (c);
307 }
3 element generic row vector
Definition: vector.h:53
double col[3]
Definition: vector.h:55
rvector rv_sub ( rvector  a,
rvector  b 
)

Subtract two vectors.

Subtract two vectors in rvector form, returning a rvector.

Parameters
avector to be subtracted from, in rvector form
bvector to be subtracted, in rvector form
Returns
the transformed vector, in rvector form
316 {
317  rvector c;
318 
319  c.col[0] = a.col[0] - b.col[0];
320  c.col[1] = a.col[1] - b.col[1];
321  c.col[2] = a.col[2] - b.col[2];
322  return (c);
323 }
3 element generic row vector
Definition: vector.h:53
double col[3]
Definition: vector.h:55
rvector rv_div ( rvector  a,
rvector  b 
)

Divide two row vectors.

Divide one rvector by another, returning a rvector.

Parameters
avector to be divided by, in rvector form
bvector to divide by, in rvector form
Returns
the transformed vector, in rvector form
332 {
333  rvector c;
334 
335  c.col[0] = a.col[0] / b.col[0];
336  c.col[1] = a.col[1] / b.col[1];
337  c.col[2] = a.col[2] / b.col[2];
338  return (c);
339 }
3 element generic row vector
Definition: vector.h:53
double col[3]
Definition: vector.h:55
rvector rv_mult ( rvector  a,
rvector  b 
)

Multiply two row vectors.

Multiply two vectors in rvector form, returning a rvector.

Parameters
afirst vector to be multiplied, in rvector form
bsecond vector to be multiplied, in rvector form
Returns
the transformed vector, in rvector form
348 {
349  rvector c;
350 
351  c.col[0] = a.col[0] * b.col[0];
352  c.col[1] = a.col[1] * b.col[1];
353  c.col[2] = a.col[2] * b.col[2];
354  return (c);
355 }
3 element generic row vector
Definition: vector.h:53
double col[3]
Definition: vector.h:55
rvector rv_cross ( rvector  a,
rvector  b 
)

Take cross product of two row vectors.

Take the vector cross product of two 3 element vectors in rvector form, returning an rvector result.

Parameters
aFirst vector in product
bSecond vector in product, treated as column vector
Returns
a cross b
364 {
365  rvector c;
366 
367  c.col[0] = a.col[1]*b.col[2] - a.col[2]*b.col[1];
368  c.col[1] = a.col[2]*b.col[0] - a.col[0]*b.col[2];
369  c.col[2] = a.col[0]*b.col[1] - a.col[1]*b.col[0];
370  return (c);
371 }
3 element generic row vector
Definition: vector.h:53
double col[3]
Definition: vector.h:55
double dot_rv ( rvector  a,
rvector  b 
)

Dot product of two row vectors.

Take the vector dot product of two vectors in rvector form.

Parameters
aFirst vector in product
bSecond vector in product, treated as column vector
Returns
a dot b
380 {
381  return a.col[0]*b.col[0] + a.col[1]*b.col[1] + a.col[2]*b.col[2];
382 }
double col[3]
Definition: vector.h:55
double sep_cv ( cvector  v1,
cvector  v2 
)

Angular separation between vectors.

Calculates the separation angle between two vectors, in radians.

Parameters
v1the first vector, in cvector format
v2the first vector, in cvector format
Returns
the separation angle in radians as a double precision

< Normalize first vector

391 {
392  double length, sepangle;
393 
394  length = sqrt(v1.x*v1.x+v1.y*v1.y+v1.z*v1.z);
395  v1.x /= length;
396  v1.y /= length;
397  v1.z /= length;
398  length = sqrt(v2.x*v2.x+v2.y*v2.y+v2.z*v2.z);
399  v2.x /= length;
400  v2.y /= length;
401  v2.z /= length;
402 
403  sepangle = ((v1.x-v2.x)*(v1.x-v2.x)+(v1.y-v2.y)*(v1.y-v2.y)+(v1.z-v2.z)*(v1.z-v2.z))/4.;
404  sepangle = 2. * atan2(sqrt(sepangle),sqrt(fmax(0.,1.-sepangle)));
405  return (sepangle);
406 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
double z
Z value.
Definition: vector.h:116
png_uint_32 length
Definition: png.c:2173
cvector cv_zero ( )

Zero cartesian vector.

Creates a zero length cartesian vector.

Returns
a cvector of zero length
413 {
414  cvector v={0.,0.,0.};
415  return (v);
416 }
3 element cartesian vector
Definition: vector.h:107
cvector cv_unitx ( )

Unit x vector.

Creates a cartesian vector with the X value set to one.

Returns
the vector
423 {
424  cvector v={1.,0.,0.};
425  return (v);
426 }
3 element cartesian vector
Definition: vector.h:107
cvector cv_unity ( )

Unit y vector.

Creates a cartesian vector with the Y value set to one.

Returns
the vector
433 {
434  cvector v={0.,1.,0.};
435  return (v);
436 }
3 element cartesian vector
Definition: vector.h:107
cvector cv_unitz ( )

Unit z vector.

Creates a cartesian vector with the Z value set to one.

Returns
the vector
443 {
444  cvector v={0.,0.,1.};
445  return (v);
446 }
3 element cartesian vector
Definition: vector.h:107
cvector cv_one ( )

Vector of ones.

Creates a cartesian vector with all values set to one.

Returns
a cvector of unit length, each dimension 1
453 {
454  cvector v={1.,1.,1.};
455  return (v);
456 }
3 element cartesian vector
Definition: vector.h:107
cvector cv_normal ( cvector  v)

Normalize cartesian vector.

Returns a normalized version of the requested cartesian vector.

Parameters
vthe cvector to be normalized
Returns
the normalized version of the vector as cvector
464 {
465  normalize_cv(v);
466  return (v);
467 }
void normalize_cv(cvector &v)
Normalize cartesian vector in place, i.e. divides it by its own norm.
Definition: vector.cpp:474
void normalize_cv ( cvector v)

Normalize cartesian vector in place, i.e. divides it by its own norm.

Normalizes requested cartesian vector.

Parameters
va pointer to the cvector to be normalized
475 {
476  double mag;
477 
478  mag = v.x*v.x + v.y*v.y + v.z*v.z;
479  if (fabs(mag - (double)0.) > D_SMALL && fabs(mag - (double)1.) > D_SMALL)
480  {
481  mag = sqrt(mag);
482  v.x /= mag;
483  v.y /= mag;
484  v.z /= mag;
485  }
486 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
const double D_SMALL
Definition: math/constants.h:40
double z
Z value.
Definition: vector.h:116
cvector cv_smult ( double  a,
cvector  b 
)

Multiply vector by scalar.

Multiply a 3 element vector by a double precision scalar.

Parameters
bvector to be tranformed, in cvector form
adouble precision scalar to multiply by
Returns
the transformed vector, in cvector form
522 {
523  cvector c;
524 
525  c.x = a * b.x;
526  c.y = a * b.y;
527  c.z = a * b.z;
528  return (c);
529 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
Definition: eci2kep_test.cpp:33
double z
Z value.
Definition: vector.h:116
3 element cartesian vector
Definition: vector.h:107
cvector cv_sadd ( double  a,
cvector  b 
)

Add scalar to each element of vector.

Add a double precision scalar to each element of a 3 element vector.

Parameters
bvector to be tranformed, in cvector form
adouble precision scalar to add
Returns
the transformed vector, in cvector form
539 {
540  cvector c;
541 
542  c.x = a + b.x;
543  c.y = a + b.y;
544  c.z = a + b.z;
545  return (c);
546 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
Definition: eci2kep_test.cpp:33
double z
Z value.
Definition: vector.h:116
3 element cartesian vector
Definition: vector.h:107
cvector cv_add ( cvector  a,
cvector  b 
)

Add two vectors.

Add two vectors in cvector form, returning a cvector.

Parameters
afirst vector to be added, in cvector form
bsecond vector to be added, in cvector form
Returns
the transformed vector, in cvector form
555 {
556  cvector c;
557 
558  c.x = a.x + b.x;
559  c.y = a.y + b.y;
560  c.z = a.z + b.z;
561  return (c);
562 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
double z
Z value.
Definition: vector.h:116
3 element cartesian vector
Definition: vector.h:107
cvector cv_sub ( cvector  a,
cvector  b 
)

Subtract two vectors.

Subtract two vectors in cvector form, returning a cvector.

Parameters
avector to be subtracted from, in cvector form
bvector to be subtracted, in cvector form
Returns
the transformed vector, in cvector form
571 {
572  cvector c;
573 
574  c.x = a.x - b.x;
575  c.y = a.y - b.y;
576  c.z = a.z - b.z;
577  return (c);
578 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
double z
Z value.
Definition: vector.h:116
3 element cartesian vector
Definition: vector.h:107
cvector cv_div ( cvector  a,
cvector  b 
)

Divide two vectors.

Divide one cvector by another, returning a cvector.

Parameters
avector to be divided by, in cvector form
bvector to divide by, in cvector form
Returns
the transformed vector, in cvector form
587 {
588  cvector c;
589 
590  c.x = a.x / b.x;
591  c.y = a.y / b.y;
592  c.z = a.z / b.z;
593  return (c);
594 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
double z
Z value.
Definition: vector.h:116
3 element cartesian vector
Definition: vector.h:107
cvector cv_mult ( cvector  a,
cvector  b 
)

Multiply two vectors.

Multiply two vectors in cvector form, returning a cvector.

Parameters
afirst vector to be multiplied, in cvector form
bsecond vector to be multiplied, in cvector form
Returns
the transformed vector, in cvector form
603 {
604  cvector c;
605 
606  c.x = a.x * b.x;
607  c.y = a.y * b.y;
608  c.z = a.z * b.z;
609  return (c);
610 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
double z
Z value.
Definition: vector.h:116
3 element cartesian vector
Definition: vector.h:107
cvector cv_cross ( cvector  a,
cvector  b 
)

Take cross product of two vectors.

614 {
615  cvector c;
616 
617  c.x = a.y*b.z - a.z*b.y;
618  c.y = a.z*b.x - a.x*b.z;
619  c.z = a.x*b.y - a.y*b.x;
620  return (c);
621 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
double z
Z value.
Definition: vector.h:116
3 element cartesian vector
Definition: vector.h:107
double dot_cv ( cvector  a,
cvector  b 
)
624 {
625  double d;
626 
627  d = a.x*b.x + a.y*b.y + a.z*b.z;
628  return (d);
629 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
double z
Z value.
Definition: vector.h:116
double length_cv ( cvector  v)
639 {
640  return cv_norm(v);
641 }
double cv_norm(cvector v)
Definition: vector.cpp:644
double cv_norm ( cvector  v)
645 {
646  double length;
647 
648  length = sqrt(v.x*v.x+v.y*v.y+v.z*v.z);
649  if (length < D_SMALL)
650  return (0.);
651  else
652  return (length);
653 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
const double D_SMALL
Definition: math/constants.h:40
double z
Z value.
Definition: vector.h:116
png_uint_32 length
Definition: png.c:2173
double norm_cv ( cvector  v)
700 {
701  double norm;
702 
703  norm = fmax(fabs(v.x),fmax(fabs(v.y),fabs(v.z)));
704 
705  return (norm);
706 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
double z
Z value.
Definition: vector.h:116
double sum_cv ( cvector  vec)
709 {
710  double sum;
711 
712  sum = vec.x + vec.y + vec.z;
713 
714  return (sum);
715 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
double z
Z value.
Definition: vector.h:116
cvector cv_sqrt ( cvector  a)
718 {
719  cvector s;
720 
721  s.x = sqrt(a.x);
722  s.y = sqrt(a.y);
723  s.z = sqrt(a.z);
724 
725  return (s);
726 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
double z
Z value.
Definition: vector.h:116
3 element cartesian vector
Definition: vector.h:107
bool equal_rv ( rvector  v1,
rvector  v2 
)

Boolean equate of row vetor.

Determine whether all the elements of an rvector are equal and return either true or false.

Parameters
v1First rvector.
v2Second rvector.
Returns
Boolean true or false.
736 {
737  return (v1.col[0] == v2.col[0] && v1.col[1] == v2.col[1] && v1.col[2] == v2.col[2]);
738 
739 }
double col[3]
Definition: vector.h:55
double length_rv ( rvector  v)

Length of row vector.

Calculate the length of a vector in row vector format.

Parameters
vVector to find the length of.
Returns
Length of row vector.
749 {
750  double length;
751 
752  length = sqrt(v.col[0]*v.col[0]+v.col[1]*v.col[1]+v.col[2]*v.col[2]);
753 
754  if (length < D_SMALL)
755  return (0.);
756  else
757  return (length);
758 }
const double D_SMALL
Definition: math/constants.h:40
png_uint_32 length
Definition: png.c:2173
double col[3]
Definition: vector.h:55
double norm_rv ( rvector  vec)

Infinite norm of row vector.

Find the largest value in a row vector.

Parameters
vecRow vector to take infinte norm of
Returns
Norm
766 {
767  return fmax(fabs(vec.col[0]),fmax(fabs(vec.col[1]),fabs(vec.col[2])));
768 }
double col[3]
Definition: vector.h:55
double sum_rv ( rvector  vec)

Sum elements of a row vector.

Add up the elements of a row vector and return the sum.

Parameters
vecRow vector to take sum of.
Returns
Sum of elements
776 {
777  return vec.col[0] + vec.col[1] + vec.col[2];
778 }
double col[3]
Definition: vector.h:55
rvector rv_sqrt ( rvector  vec)

Row vector square root.

rvector whose elements are the square roots of the elements of a rvector.

Parameters
vecrvector to take the square root of.
Returns
rvector of resultant square roots.
787 {
788  rvector s;
789 
790  s.col[0] = sqrt(vec.col[0]);
791  s.col[1] = sqrt(vec.col[1]);
792  s.col[2] = sqrt(vec.col[2]);
793 
794  return (s);
795 }
3 element generic row vector
Definition: vector.h:53
double col[3]
Definition: vector.h:55
LsFit::LsFit ( uint16_t  cnt = 10,
uint16_t  ord = 2 
)

Multi element, variable order least squares fit.

Constructor for cnt element, ord order least squares fit.

Parameters
cntNumber of elements to be fit.
ordOrder of fit.
Returns
Least squares fit object.
2423 {
2424  initialize(cnt, ord);
2425 // if (ord)
2426 // {
2427 // order = ord;
2428 // }
2429 // else
2430 // {
2431 // order = 1;
2432 // }
2433 // if (cnt)
2434 // {
2435 // element_cnt = cnt;
2436 // }
2437 // else
2438 // {
2439 // element_cnt = order + 1;
2440 // }
2441 // var.resize(0);
2442 // depth = 0;
2443 }
void initialize(uint16_t cnt=10, uint16_t ord=2)
Initialize Least Squares Fit.
Definition: mathlib.cpp:2450
void LsFit::initialize ( uint16_t  cnt = 10,
uint16_t  ord = 2 
)

Initialize Least Squares Fit.

Perform setting of variables in LsFit so that it can be ready for use.

Parameters
cntNumber of elements to be fit.
ordOrder of fit.
2451 {
2452  if (ord)
2453  {
2454  order = ord;
2455  }
2456  else
2457  {
2458  order = 1;
2459  }
2460  if (cnt)
2461  {
2462  element_cnt = cnt;
2463  }
2464  else
2465  {
2466  element_cnt = order + 1;
2467  }
2468  var.resize(0);
2469  depth = 0;
2470 }
deque< fitelement > var
Definition: mathlib.h:389
uint16_t depth
Definition: mathlib.h:383
uint32_t order
Definition: mathlib.h:385
uint16_t element_cnt
Least Squares Fit Structure.
Definition: mathlib.h:381
void LsFit::update ( double  x,
double  y 
)

Update scalar Least Squares Fit.

Add independent and dependent value pair to existing LsFit, updating the fit. If the number of elements in the fit has been reached, the oldest element is dropped before fitting.

Parameters
xDependent value.
yIndependent value.
2479 {
2480  fitelement cfit;
2481 
2482  // Independent variable
2483  cfit.x = x;
2484 
2485  // Dependent variable for quaternion
2486  cfit.y.a4[0] = y;
2487 
2488  update(cfit, 1);
2489 }
y
Definition: inputfile.py:6
x
Definition: inputfile.py:6
void update(double x, double y)
Update scalar Least Squares Fit.
Definition: mathlib.cpp:2478
void LsFit::update ( double  x,
rvector  y 
)

Update rvector Least Squares Fit.

Add independent and dependent value pair to existing LsFit, updating the fit. If the number of elements in the fit has been reached, the oldest element is dropped before fitting.

Parameters
xDependent value.
yIndependent values.
2498 {
2499  fitelement cfit;
2500 
2501  // Independent variable
2502  cfit.x = x;
2503 
2504  // Dependent variable for quaternion
2505  cfit.y.r = y;
2506 
2507  update(cfit, 3);
2508 }
y
Definition: inputfile.py:6
x
Definition: inputfile.py:6
void update(double x, double y)
Update scalar Least Squares Fit.
Definition: mathlib.cpp:2478
void LsFit::update ( double  x,
gvector  y 
)

Update gvector Least Squares Fit.

Add independent and dependent value pair to existing LsFit, updating the fit. If the number of elements in the fit has been reached, the oldest element is dropped before fitting.

Parameters
xDependent value.
yIndependent values.
2517 {
2518  fitelement cfit;
2519 
2520  // Independent variable
2521  cfit.x = x;
2522 
2523  // Dependent variable for quaternion
2524  cfit.y.g = y;
2525 
2526  update(cfit, 3);
2527 }
y
Definition: inputfile.py:6
x
Definition: inputfile.py:6
void update(double x, double y)
Update scalar Least Squares Fit.
Definition: mathlib.cpp:2478
void LsFit::update ( double  x,
quaternion  y 
)

Update quaternion Least Squares Fit.

Add independent and dependent value pair to existing LsFit, updating the fit. If the number of elements in the fit has been reached, the oldest element is dropped before fitting.

Parameters
xDependent value.
yIndependent values.
2536 {
2537  fitelement cfit;
2538 
2539  // Independent variable
2540  cfit.x = x;
2541 
2542  // Dependent variable for quaternion
2543  cfit.y.q = y;
2544 
2545  update(cfit, 4);
2546 
2547  // // TODO: check this problem!!!
2548  // //if (parms.size() && (std::isnan(parms[0][0]) || std::isnan(parms[0][1] || std::isnan(parms[0][2]))))
2549  // if (1)
2550  // {
2551  // for (uint16_t i=0; i<var.size(); ++i)
2552  // {
2553  // // TODO: why are we printing this?
2554  // printf("%.15g [%g %g %g %g]\n", var[i].x, var[i].y.q.w, var[i].y.q.d.x, var[i].y.q.d.y, var[i].y.q.d.z);
2555  // }
2556  // printf("\n");
2557  // }
2558 }
y
Definition: inputfile.py:6
x
Definition: inputfile.py:6
void update(double x, double y)
Update scalar Least Squares Fit.
Definition: mathlib.cpp:2478
void LsFit::update ( fitelement  cfit,
uint16_t  dep 
)

Update generic Least Squares Fit.

Add LsFit::fitelement containing independent and dependent value pair to existing LsFit, updating the fit. If the number of elements in the fit has been reached, the oldest element is dropped before fitting.

Parameters
cfitIndependent and dependent values.
depDepth of fit.
2567 {
2568  if (var.size() && cfit.x == var[var.size()-1].x)
2569  {
2570  return;
2571  }
2572 
2573  depth = dep;
2574 
2575  // Sudden switch to mirror value will wreak havoc with fit of quaternion
2576  if (depth == 4 && var.size() > 0)
2577  {
2578  // If new value is closer to mirror of last value, then switch all previous values
2579  if (length_q(q_sub(cfit.y.q, var[var.size()-1].y.q)) > length_q(q_sub(q_smult(-1., cfit.y.q), var[var.size()-1].y.q)))
2580  {
2581  for (uint16_t i=0; i<var.size(); ++i)
2582  {
2583  var[i].y.q = q_smult(-1., var[i].y.q);
2584  }
2585  }
2586  }
2587 
2588  // Save to FIFO
2589  var.push_back(cfit);
2590 
2591  // Element_cnt element collected and we can start cycling FIFO
2592  if (var.size() > element_cnt)
2593  {
2594  var.pop_front();
2595  }
2596 
2597  // More than order elements collected and we can start fitting
2598  if (var.size() > order)
2599  {
2600  LsFit::fit();
2601  }
2602 }
quaternion q_sub(quaternion a, quaternion b)
Subtract two quaternions.
Definition: vector.cpp:1186
y
Definition: inputfile.py:6
int i
Definition: rw_test.cpp:37
quaternion q_smult(double a, quaternion b)
Multiply quaternion by scalar.
Definition: vector.cpp:1151
double length_q(quaternion q)
Length of quaternion.
Definition: vector.cpp:1329
deque< fitelement > var
Definition: mathlib.h:389
uint16_t depth
Definition: mathlib.h:383
void fit()
Calculate least squares fit.
Definition: mathlib.cpp:2608
uint32_t order
Definition: mathlib.h:385
uint16_t element_cnt
Least Squares Fit Structure.
Definition: mathlib.h:381
void LsFit::fit ( )
private

Calculate least squares fit.

Calculate least squares fit for each axis of exsiting LsFit. Parameters are updated to reflect new fit.

2609 {
2610  // Minimize independent variable by zero offsetting
2611  basex = var[0].x;
2612 
2613  // For each independent, calculate sums of powers
2614  vector<double> sumx(2*order+1);
2615  sumx[0] = var.size();
2616  for (uint16_t i=0; i<sumx[0]; ++i)
2617  {
2618  double ix = 1.;
2619  double cx = var[i].x - basex;
2620  for (uint16_t j=0; j<2*order; ++j)
2621  {
2622  ix *= cx;
2623  sumx[j+1] += ix;
2624  }
2625  }
2626  meanx = sumx[1] / sumx[0];
2627  stdevx = sqrt((sumx[2] - sumx[1]*sumx[1]/sumx[0])/(sumx[0]-1));
2628 
2629  // Calculate sums of products of dependent and independent and do least squares fit
2630  parms.resize(depth);
2631  for (uint16_t i=0; i<depth; ++i)
2632  {
2633  vector<double> sumxy(order+1);
2634  stdevy.a4[i] = 0;
2635  for (uint16_t j=0; j<var.size(); ++j)
2636  {
2637  double ixy = var[j].y.a4[i];
2638  stdevy.a4[i] += ixy * ixy;
2639  double cx = var[j].x - basex;
2640  for (uint16_t k=0; k<order+1; ++k)
2641  {
2642  sumxy[k] += ixy;
2643  ixy *= cx;
2644  }
2645  }
2646  meany.a4[i] = sumxy[0] / var.size();
2647  stdevy.a4[i] = sqrt((stdevy.a4[i] - sumxy[0]*sumxy[0]/sumx[0])/(sumx[0]-1));
2648 
2649  vector< vector<double> > xs(order+1, vector<double>(order+1));
2650  vector<double> ys(order+1);
2651  vector<double> tx(order+1);
2652  for (uint16_t j=0; j<order+1; ++j)
2653  {
2654  for (uint16_t k=0; k<order+1; ++k)
2655  {
2656  tx[k] = sumx[k+j];
2657  }
2658  xs[j] = tx;
2659  ys[j] = sumxy[j];
2660  }
2661  parms[i].resize(ys.size());
2662  multisolve(xs, ys, parms[i]);
2663  // printf("%u %u %f %f %f\n", i, parms[i].size(), parms[i][0], parms[i][1], parms[i][2]);
2664  }
2665 }
vector< vector< double > > parms
Definition: mathlib.h:391
uvector meany
Definition: mathlib.h:397
void multisolve(vector< vector< double > > x, vector< double > y, vector< double > &a)
Perform N equation solution.
Definition: mathlib.cpp:900
uvector stdevy
Definition: mathlib.h:399
int i
Definition: rw_test.cpp:37
double basex
Definition: mathlib.h:387
deque< fitelement > var
Definition: mathlib.h:389
double meanx
Definition: mathlib.h:396
uint16_t depth
Definition: mathlib.h:383
uint32_t order
Definition: mathlib.h:385
double a4[4]
Definition: mathlib.h:167
double stdevx
Definition: mathlib.h:398
double LsFit::lastx ( )

Least squares last independent value.

Return the value of the independent value added at the most recent LsFit::update.

Returns
Most recently updated independent value.
2672 {
2673  if (var.size())
2674  {
2675  return var[var.size()-1].x;
2676  }
2677  else
2678  {
2679  return 0.;
2680  }
2681 }
deque< fitelement > var
Definition: mathlib.h:389
double LsFit::firstx ( )

Least squares first independent value.

Return the value of the independent value added at the least recent LsFit::update.

Returns
Least recently updated independent value.
2688 {
2689  if (var.size())
2690  {
2691  return var[0].x;
2692  }
2693  else
2694  {
2695  return 0.;
2696  }
2697 }
deque< fitelement > var
Definition: mathlib.h:389
size_t LsFit::size ( )

Least squares number of values.

Return the number of data values in the fit as of the last LsFit::update.

Returns
Last size of fit.
2704 {
2705  return var.size();
2706 }
deque< fitelement > var
Definition: mathlib.h:389
double LsFit::eval ( double  x)

Least squares dependent scalar value.

Return the value of the dependent scalar, calculated for the provided independent value, using the parameters from the latest LsFit::update.

Parameters
xIndependent value.
Returns
Calculated scalar dependent value.
2715 {
2716  if (var.size() > order)
2717  {
2718  return evaluate_poly(x - basex, parms[0]);
2719  }
2720  else
2721  {
2722  return 0.;
2723  }
2724 
2725 }
vector< vector< double > > parms
Definition: mathlib.h:391
double basex
Definition: mathlib.h:387
double evaluate_poly(double x, vector< double > parms)
Evaluate polynomial.
Definition: mathlib.cpp:955
deque< fitelement > var
Definition: mathlib.h:389
x
Definition: inputfile.py:6
uint32_t order
Definition: mathlib.h:385
rvector LsFit::evalrvector ( double  x)

Least squares dependent rvector value.

Return the value of the dependent rvector, calculated for the provided independent value, using the parameters from the latest LsFit::update.

Parameters
xIndependent value.
Returns
Calculated rvector dependent value.
2734 {
2735  if (var.size() > order)
2736  {
2737  return rv_evaluate_poly(x - basex, parms);
2738  }
2739  else
2740  {
2741  return rv_zero();
2742  }
2743 }
vector< vector< double > > parms
Definition: mathlib.h:391
double basex
Definition: mathlib.h:387
rvector rv_evaluate_poly(double x, vector< vector< double > > parms)
Evaluate vector polynomial.
Definition: mathlib.cpp:1059
deque< fitelement > var
Definition: mathlib.h:389
x
Definition: inputfile.py:6
uint32_t order
Definition: mathlib.h:385
rvector rv_zero()
Zero row order vector.
Definition: vector.cpp:107
gvector LsFit::evalgvector ( double  x)

Least squares dependent gvector value.

Return the value of the dependent gvector, calculated for the provided independent value, using the parameters from the latest LsFit::update.

Parameters
xIndependent value.
Returns
Calculated gvector dependent value.
2752 {
2753  if (var.size() > order)
2754  {
2755  return gv_evaluate_poly(x - basex, parms);
2756  }
2757  else
2758  {
2759  return gv_zero();
2760  }
2761 }
vector< vector< double > > parms
Definition: mathlib.h:391
double basex
Definition: mathlib.h:387
deque< fitelement > var
Definition: mathlib.h:389
gvector gv_zero()
Zero geodetic vector.
Definition: vector.cpp:946
x
Definition: inputfile.py:6
uint32_t order
Definition: mathlib.h:385
gvector gv_evaluate_poly(double x, vector< vector< double > > parms)
Evaluate vector polynomial.
Definition: mathlib.cpp:1183
quaternion LsFit::evalquaternion ( double  x)

Least squares dependent quaternion value.

Return the value of the dependent quaternion, calculated for the provided independent value, using the parameters from the latest LsFit::update.

Parameters
xIndependent value.
Returns
Calculated quaternion dependent value.
2770 {
2771  if (var.size() > order)
2772  {
2773  return q_evaluate_poly(x - basex, parms);
2774  }
2775  else
2776  {
2777  return q_zero();
2778  }
2779 }
vector< vector< double > > parms
Definition: mathlib.h:391
double basex
Definition: mathlib.h:387
quaternion q_zero()
Zero quaternion.
Definition: vector.cpp:1003
deque< fitelement > var
Definition: mathlib.h:389
x
Definition: inputfile.py:6
uint32_t order
Definition: mathlib.h:385
quaternion q_evaluate_poly(double x, vector< vector< double > > parms)
Evaluate quaternion polynomial.
Definition: mathlib.cpp:1307
double LsFit::slope ( double  x)

Least squares dependent scalar 1st derivative.

Return the value of the dependent scalar 1st derivative, calculated for the provided independent value, using the parameters from the latest LsFit::update.

Parameters
xIndependent value.
Returns
Calculated scalar dependent 1st derivative.
2788 {
2789  if (var.size() > order)
2790  {
2791  return evaluate_poly_slope(x - basex, parms[0]);
2792  }
2793  else
2794  {
2795  return 0.;
2796  }
2797 }
vector< vector< double > > parms
Definition: mathlib.h:391
double basex
Definition: mathlib.h:387
double evaluate_poly_slope(double x, vector< double > parms)
Evaluate polynomial slope.
Definition: mathlib.cpp:981
deque< fitelement > var
Definition: mathlib.h:389
x
Definition: inputfile.py:6
uint32_t order
Definition: mathlib.h:385
rvector LsFit::slopervector ( double  x)

Least squares dependent rvector 1st derivative.

Return the value of the dependent rvector 1st derivative, calculated for the provided independent value, using the parameters from the latest LsFit::update.

Parameters
xIndependent value.
Returns
Calculated rvector dependent 1st derivative.
2806 {
2807  if (var.size() > order)
2808  {
2809  return rv_evaluate_poly_slope(x - basex, parms);
2810  }
2811  else
2812  {
2813  return rv_zero();
2814  }
2815 }
vector< vector< double > > parms
Definition: mathlib.h:391
double basex
Definition: mathlib.h:387
deque< fitelement > var
Definition: mathlib.h:389
x
Definition: inputfile.py:6
uint32_t order
Definition: mathlib.h:385
rvector rv_zero()
Zero row order vector.
Definition: vector.cpp:107
rvector rv_evaluate_poly_slope(double x, vector< vector< double > > parms)
Evaluate vector polynomial slope.
Definition: mathlib.cpp:1090
gvector LsFit::slopegvector ( double  x)

Least squares dependent gvector 1st derivative.

Return the value of the dependent gvector 1st derivative, calculated for the provided independent value, using the parameters from the latest LsFit::update.

Parameters
xIndependent value.
Returns
Calculated gvector dependent 1st derivative.
2824 {
2825  if (var.size() > order)
2826  {
2827  return gv_evaluate_poly_slope(x - basex, parms);
2828  }
2829  else
2830  {
2831  return gv_zero();
2832  }
2833 }
vector< vector< double > > parms
Definition: mathlib.h:391
double basex
Definition: mathlib.h:387
gvector gv_evaluate_poly_slope(double x, vector< vector< double > > parms)
Evaluate vector polynomial slope.
Definition: mathlib.cpp:1214
deque< fitelement > var
Definition: mathlib.h:389
gvector gv_zero()
Zero geodetic vector.
Definition: vector.cpp:946
x
Definition: inputfile.py:6
uint32_t order
Definition: mathlib.h:385
quaternion LsFit::slopequaternion ( double  x)

Least squares dependent quaternion 1st derivative.

Return the value of the dependent quaternion 1st derivative, calculated for the provided independent value, using the parameters from the latest LsFit::update.

Parameters
xIndependent value.
Returns
Calculated quaternion dependent 1st derivative.
2842 {
2843  if (var.size() > order)
2844  {
2845  return q_evaluate_poly_slope(x - basex, parms);
2846  }
2847  else
2848  {
2849  return q_zero();
2850  }
2851 }
vector< vector< double > > parms
Definition: mathlib.h:391
double basex
Definition: mathlib.h:387
quaternion q_zero()
Zero quaternion.
Definition: vector.cpp:1003
deque< fitelement > var
Definition: mathlib.h:389
x
Definition: inputfile.py:6
uint32_t order
Definition: mathlib.h:385
quaternion q_evaluate_poly_slope(double x, vector< vector< double > > parms)
Evaluate quaternion polynomial slope.
Definition: mathlib.cpp:1338
double LsFit::accel ( double  x)

Least squares dependent scalar 2nd derivative.

Return the value of the dependent scalar 2nd derivative, calculated for the provided independent value, using the parameters from the latest LsFit::update.

Parameters
xIndependent value.
Returns
Calculated scalar dependent 2nd derivative.
2860 {
2861  if (var.size() > order)
2862  {
2863  return evaluate_poly_accel(x - basex, parms[0]);
2864  }
2865  else
2866  {
2867  return 0.;
2868  }
2869 }
vector< vector< double > > parms
Definition: mathlib.h:391
double basex
Definition: mathlib.h:387
deque< fitelement > var
Definition: mathlib.h:389
x
Definition: inputfile.py:6
uint32_t order
Definition: mathlib.h:385
double evaluate_poly_accel(double x, vector< double > parms)
Evaluate polynomial acceleration.
Definition: mathlib.cpp:1007
rvector LsFit::accelrvector ( double  x)

Least squares dependent rvector 2nd derivative.

Return the value of the dependent rvector 2nd derivative, calculated for the provided independent value, using the parameters from the latest LsFit::update.

Parameters
xIndependent value.
Returns
Calculated rvector dependent 2nd derivative.
2878 {
2879  if (var.size() > order)
2880  {
2881  return rv_evaluate_poly_accel(x - basex, parms);
2882  }
2883  else
2884  {
2885  return rv_zero();
2886  }
2887 }
vector< vector< double > > parms
Definition: mathlib.h:391
double basex
Definition: mathlib.h:387
deque< fitelement > var
Definition: mathlib.h:389
x
Definition: inputfile.py:6
uint32_t order
Definition: mathlib.h:385
rvector rv_zero()
Zero row order vector.
Definition: vector.cpp:107
rvector rv_evaluate_poly_accel(double x, vector< vector< double > > parms)
Evaluate vector polynomial acceleration.
Definition: mathlib.cpp:1121
gvector LsFit::accelgvector ( double  x)

Least squares dependent gvector 2nd derivative.

Return the value of the dependent gvector 2nd derivative, calculated for the provided independent value, using the parameters from the latest LsFit::update.

Parameters
xIndependent value.
Returns
Calculated gvector dependent 2nd derivative.
2896 {
2897  if (var.size() > order)
2898  {
2899  return gv_evaluate_poly_accel(x - basex, parms);
2900  }
2901  else
2902  {
2903  return gv_zero();
2904  }
2905 }
vector< vector< double > > parms
Definition: mathlib.h:391
double basex
Definition: mathlib.h:387
deque< fitelement > var
Definition: mathlib.h:389
gvector gv_zero()
Zero geodetic vector.
Definition: vector.cpp:946
x
Definition: inputfile.py:6
uint32_t order
Definition: mathlib.h:385
gvector gv_evaluate_poly_accel(double x, vector< vector< double > > parms)
Evaluate vector polynomial acceleration.
Definition: mathlib.cpp:1245
quaternion LsFit::accelquaternion ( double  x)

Least squares dependent quaternion 2nd derivative.

Return the value of the dependent quaternion 2nd derivative, calculated for the provided independent value, using the parameters from the latest LsFit::update.

Parameters
xIndependent value.
Returns
Calculated quaternion dependent 2nd derivative.
2914 {
2915  if (var.size() > order)
2916  {
2917  return q_evaluate_poly_accel(x - basex, parms);
2918  }
2919  else
2920  {
2921  return q_zero();
2922  }
2923 }
vector< vector< double > > parms
Definition: mathlib.h:391
double basex
Definition: mathlib.h:387
quaternion q_evaluate_poly_accel(double x, vector< vector< double > > parms)
Evaluate quaternion polynomial acceleration.
Definition: mathlib.cpp:1369
quaternion q_zero()
Zero quaternion.
Definition: vector.cpp:1003
deque< fitelement > var
Definition: mathlib.h:389
x
Definition: inputfile.py:6
uint32_t order
Definition: mathlib.h:385
double LsFit::jerk ( double  x)

Least squares dependent scalar 3rd derivative.

Return the value of the dependent scalar 3rd derivative, calculated for the provided independent value, using the parameters from the latest LsFit::update.

Parameters
xIndependent value.
Returns
Calculated scalar dependent 3rd derivative.
2932 {
2933  if (var.size() > order)
2934  {
2935  return evaluate_poly_jerk(x - basex, parms[0]);
2936  }
2937  else
2938  {
2939  return 0.;
2940  }
2941 }
vector< vector< double > > parms
Definition: mathlib.h:391
double basex
Definition: mathlib.h:387
double evaluate_poly_jerk(double x, vector< double > parms)
Evaluate polynomial jerk.
Definition: mathlib.cpp:1033
deque< fitelement > var
Definition: mathlib.h:389
x
Definition: inputfile.py:6
uint32_t order
Definition: mathlib.h:385
rvector LsFit::jerkrvector ( double  x)

Least squares dependent rvector 3rd derivative.

Return the value of the dependent rvector 3rd derivative, calculated for the provided independent value, using the parameters from the latest LsFit::update.

Parameters
xIndependent value.
Returns
Calculated rvector dependent 3rd derivative.
2950 {
2951  if (var.size() > order)
2952  {
2953  return rv_evaluate_poly_jerk(x - basex, parms);
2954  }
2955  else
2956  {
2957  return rv_zero();
2958  }
2959 }
vector< vector< double > > parms
Definition: mathlib.h:391
double basex
Definition: mathlib.h:387
deque< fitelement > var
Definition: mathlib.h:389
x
Definition: inputfile.py:6
uint32_t order
Definition: mathlib.h:385
rvector rv_evaluate_poly_jerk(double x, vector< vector< double > > parms)
Evaluate vector polynomial jerk.
Definition: mathlib.cpp:1152
rvector rv_zero()
Zero row order vector.
Definition: vector.cpp:107
gvector LsFit::jerkgvector ( double  x)

Least squares dependent gvector 3rd derivative.

Return the value of the dependent gvector 3rd derivative, calculated for the provided independent value, using the parameters from the latest LsFit::update.

Parameters
xIndependent value.
Returns
Calculated gvector dependent 3rd derivative.
2968 {
2969  if (var.size() > order)
2970  {
2971  return gv_evaluate_poly_jerk(x - basex, parms);
2972  }
2973  else
2974  {
2975  return gv_zero();
2976  }
2977 }
vector< vector< double > > parms
Definition: mathlib.h:391
double basex
Definition: mathlib.h:387
deque< fitelement > var
Definition: mathlib.h:389
gvector gv_zero()
Zero geodetic vector.
Definition: vector.cpp:946
x
Definition: inputfile.py:6
uint32_t order
Definition: mathlib.h:385
gvector gv_evaluate_poly_jerk(double x, vector< vector< double > > parms)
Evaluate vector polynomial jerk.
Definition: mathlib.cpp:1276
quaternion LsFit::jerkquaternion ( double  x)

Least squares dependent quaternion 3rd derivative.

Return the value of the dependent quaternion 3rd derivative, calculated for the provided independent value, using the parameters from the latest LsFit::update.

Parameters
xIndependent value.
Returns
Calculated quaternion dependent 3rd derivative.
2986 {
2987  if (var.size() > order)
2988  {
2989  return q_evaluate_poly_jerk(x - basex, parms);
2990  }
2991  else
2992  {
2993  return q_zero();
2994  }
2995 }
vector< vector< double > > parms
Definition: mathlib.h:391
double basex
Definition: mathlib.h:387
quaternion q_evaluate_poly_jerk(double x, vector< vector< double > > parms)
Evaluate quaternion polynomial jerk.
Definition: mathlib.cpp:1400
quaternion q_zero()
Zero quaternion.
Definition: vector.cpp:1003
deque< fitelement > var
Definition: mathlib.h:389
x
Definition: inputfile.py:6
uint32_t order
Definition: mathlib.h:385
vector< vector< double > > LsFit::getparms ( )

Least Squares parameters.

Return the values of the parameters for all axes of generated from the latest LsFit::update.

Returns
Parameters
3002 {
3003  return parms;
3004 }
vector< vector< double > > parms
Definition: mathlib.h:391
double LsFit::getbasex ( )

Least Squares base.

Return the values of the base independent variable for the latest LsFit::update.

Returns
Parameters
3011 {
3012  return basex;
3013 }
double basex
Definition: mathlib.h:387
double DCM::dotProduct ( cvector  a,
cvector  b 
)
232 {
233  return (a.x * b.x + a.y * b.y + a.z * b.z);
234 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
double z
Z value.
Definition: vector.h:116
cmatrix DCM::transposeMatrix ( cmatrix  a)
240 {
241  cmatrix b;
242 
243  b.r1.x = a.r1.x;
244  b.r1.y = a.r2.x;
245  b.r1.z = a.r3.x;
246 
247  b.r2.x = a.r1.y;
248  b.r2.y = a.r2.y;
249  b.r2.z = a.r3.y;
250 
251  b.r3.x = a.r1.z;
252  b.r3.y = a.r2.z;
253  b.r3.z = a.r3.z;
254 
255  return (b);
256 }
double y
Y value.
Definition: vector.h:114
3x3 element cartesian matrix
Definition: matrix.h:96
cvector r1
Row 1.
Definition: matrix.h:99
double x
X value.
Definition: vector.h:112
long b
Definition: jpegint.h:371
double z
Z value.
Definition: vector.h:116
cvector r3
Row 3.
Definition: matrix.h:103
cvector r2
Row 2.
Definition: matrix.h:101
cmatrix DCM::base2_from_base1 ( basisOrthonormal  base2,
basisOrthonormal  base1 
)
258  {
259 
260  // compute dcm matrix (A) to represent vector in base 2 coordinates
261 
262  // References
263  // - Quaternion and Rotation Sequences, Kuipers, pg 161 eq 7.8 (I think these
264  // formulas are wrong in the book! they are inversed, must check)
265  // - http://people.ae.illinois.edu/tbretl/ae403/handouts/06-dcm.pdf (this
266  // reference seems to be sound)
267  // - http://www.starlino.com/dcm_tutorial.html (eq 1.4)
268 
269  // example: vector_body_coordinates = dcm_base2_from_base1 * vector_inertial_coodinates
270 
271  // Notes:
272  // - matrix A represents a frame rotation that relates the initial reference
273  // frame {X,Y,Z} to a rotated frame {x,y,z}: x = AX. Example, appliying
274  // the operation AX - where X is for instance a vector in the inertial
275  // frame - results in the coordinates of that vector in the new frame {x,y,z}
276  // in general the vector x is simply the vector X expressed in a new frame coordinates
277 
278  // Example to run:
279  // define a base1 (inertial)
280  // base1.i = [1,0,0];
281  // base1.j = [0,1,0];
282  // base1.k = [0,0,1];
283  // define a frame2 rotated 90 deg aroud z axis (of the inertial)
284  // base2.i = [0,1,0];
285  // base2.j = [-1,0,0];
286  // base2.k = [0,0,1];
287  // compute the dcm
288  // dcm_base2_from_base1(base2,base1);
289  // check if it's right, base1.i vector shoud be now [0,-1,0]
290  // dcm_base2_from_base1(base2,base1)*base1.i'
291 
292  // TODOs
293  // - add validation step to verify if the bases are orthogonal
294 
295  // normalize the basis just in case
296 
297  base1.normalize();
298  base2.normalize();
299 
300 
301  return {
302  { dotProduct(base1.i, base2.i) , dotProduct(base1.j, base2.i), dotProduct(base1.k, base2.i) },
303  { dotProduct(base1.i, base2.j) , dotProduct(base1.j, base2.j), dotProduct(base1.k, base2.j) },
304  { dotProduct(base1.i, base2.k) , dotProduct(base1.j, base2.k), dotProduct(base1.k, base2.k) }
305  };
306 }
cvector k
Definition: rotation.h:54
double dotProduct(cvector a, cvector b)
Definition: rotation.cpp:231
cvector i
Definition: rotation.h:52
cvector j
Definition: rotation.h:53
void normalize()
Definition: rotation.cpp:325
cmatrix DCM::base1_from_base2 ( basisOrthonormal  base1,
basisOrthonormal  base2 
)
310  {
311  // compute dcm matrix (A) to represent vector given in base 2
312  // in base 1 coordinates
313 
314  // example:
315  // vector_inertial_coordinates = dcm_base1_from_base2 * vector_body_coodinates
316 
317  // TODOs
318  // - add validation step to verify if the bases are orthogonal
319 
320  // just compute the transpose
321  return transposeMatrix(base2_from_base1(base2,base1));
322 }
cmatrix base2_from_base1(basisOrthonormal base2, basisOrthonormal base1)
Definition: rotation.cpp:258
cmatrix transposeMatrix(cmatrix a)
Definition: rotation.cpp:239
void basisOrthonormal::normalize ( )
326 {
327  this->i.normalize();
328  this->j.normalize();
329  this->k.normalize();
330 }
cvector k
Definition: rotation.h:54
void normalize(double scale=1.)
Normalize cartesian vector in place, i.e. divides it by its own norm.
Definition: vector.cpp:489
cvector i
Definition: rotation.h:52
cvector j
Definition: rotation.h:53
void cvector::normalize ( double  scale = 1.)

Normalize cartesian vector in place, i.e. divides it by its own norm.

490 {
491  double weight = scale / norm();
492 
493  if (fabs(weight - (double)0.) > D_SMALL && fabs(weight - (double)1.) > D_SMALL)
494  {
495  x *= weight;
496  y *= weight;
497  z *= weight;
498  }
499 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
const double D_SMALL
Definition: math/constants.h:40
double z
Z value.
Definition: vector.h:116
double norm()
Definition: vector.cpp:664
cvector cvector::normalized ( double  scale = 1.)

Normalize cartesian vector.

Returns a normalized version of the requested cartesian vector.

Parameters
scalethe weight to be applied after normalizing
Returns
the normalized version of the vector as cvector
507 {
508  cvector newv = *this;
509  newv.normalize(scale);
510  return newv;
511 }
void normalize(double scale=1.)
Normalize cartesian vector in place, i.e. divides it by its own norm.
Definition: vector.cpp:489
3 element cartesian vector
Definition: vector.h:107
double cvector::norm2 ( )
656 {
657  double norm = (this->x*this->x + this->y*this->y + this->z*this->z);
658  if (norm < D_SMALL)
659  return (0.);
660  else
661  return (norm);
662 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
const double D_SMALL
Definition: math/constants.h:40
double z
Z value.
Definition: vector.h:116
double norm()
Definition: vector.cpp:664
double cvector::norm ( )
665 {
666  double norm = sqrt(this->norm2());
667  if (norm < D_SMALL)
668  return (0.);
669  else
670  return (norm);
671 }
double norm2()
Definition: vector.cpp:655
const double D_SMALL
Definition: math/constants.h:40
double norm()
Definition: vector.cpp:664
double cvector::length ( )
674 {
675  return (this->norm());
676 }
double norm()
Definition: vector.cpp:664
double & cvector::operator[] ( const int  index)

Index into cvector.

680 {
681  switch (index)
682  {
683  case 0:
684  default:
685  return x;
686  break;
687  case 1:
688  return y;
689  break;
690  case 2:
691  return z;
692  break;
693  }
694 }
double y
Y value.
Definition: vector.h:114
double x
X value.
Definition: vector.h:112
double z
Z value.
Definition: vector.h:116