9 #include <libnova/julian_day.h>
11 #include <gsl/gsl_blas.h>
12 #include <gsl/gsl_permutation.h>
13 #include <gsl/gsl_linalg.h>
23 namespace AlignmentSubsystem
55 switch (SyncPoints.size())
91 DummyActualDirectionCosine2.
x = 0.0;
92 DummyActualDirectionCosine2.
y = 0.0;
93 DummyActualDirectionCosine2.
z = 1.0;
94 DummyApparentDirectionCosine2 = DummyActualDirectionCosine2;
105 DummyApparentDirectionCosine2 = DummyActualDirectionCosine2;
116 DummyApparentDirectionCosine2 = DummyActualDirectionCosine2;
120 DummyActualDirectionCosine3 = ActualDirectionCosine1 * DummyActualDirectionCosine2;
122 DummyApparentDirectionCosine3 = Entry1.
TelescopeDirection * DummyApparentDirectionCosine2;
123 DummyApparentDirectionCosine3.
Normalise();
164 DummyActualDirectionCosine3 = ActualDirectionCosine1 * ActualDirectionCosine2;
167 DummyApparentDirectionCosine3.
Normalise();
236 int VertexNumber = 1;
238 for (InMemoryDatabase::AlignmentDatabaseType::const_iterator Itr = SyncPoints.begin();
239 Itr != SyncPoints.end(); Itr++)
258 ActualDirectionCosine.
z, VertexNumber);
260 (*Itr).TelescopeDirection.z, VertexNumber);
277 #ifdef CONVEX_HULL_DEBUGGING
280 if (
nullptr != CurrentFace)
284 #ifdef CONVEX_HULL_DEBUGGING
290 #ifdef CONVEX_HULL_DEBUGGING
291 ASSDEBUGF(
"Initialise - Ignoring actual face %d", ActualFaces);
296 #ifdef CONVEX_HULL_DEBUGGING
297 ASSDEBUGF(
"Initialise - Processing actual face %d v1 %d v2 %d v3 %d", ActualFaces,
304 SyncPoints[CurrentFace->
vertex[0]->
vnum - 1].TelescopeDirection,
305 SyncPoints[CurrentFace->
vertex[1]->
vnum - 1].TelescopeDirection,
306 SyncPoints[CurrentFace->
vertex[2]->
vnum - 1].TelescopeDirection,
307 CurrentFace->
pMatrix,
nullptr);
309 CurrentFace = CurrentFace->
next;
316 #ifdef CONVEX_HULL_DEBUGGING
317 int ApparentFaces = 0;
319 if (
nullptr != CurrentFace)
323 #ifdef CONVEX_HULL_DEBUGGING
329 #ifdef CONVEX_HULL_DEBUGGING
330 ASSDEBUGF(
"Initialise - Ignoring apparent face %d", ApparentFaces);
335 #ifdef CONVEX_HULL_DEBUGGING
336 ASSDEBUGF(
"Initialise - Processing apparent face %d v1 %d v2 %d v3 %d", ApparentFaces,
341 SyncPoints[CurrentFace->
vertex[1]->
vnum - 1].TelescopeDirection,
342 SyncPoints[CurrentFace->
vertex[2]->
vnum - 1].TelescopeDirection,
346 CurrentFace->
pMatrix,
nullptr);
348 CurrentFace = CurrentFace->
next;
353 #ifdef CONVEX_HULL_DEBUGGING
354 ASSDEBUGF(
"Initialise - ActualFaces %d ApparentFaces %d", ActualFaces, ApparentFaces);
379 switch (SyncPoints.size())
421 gsl_vector *pGSLActualVector = gsl_vector_alloc(3);
422 gsl_vector_set(pGSLActualVector, 0, ActualVector.
x);
423 gsl_vector_set(pGSLActualVector, 1, ActualVector.
y);
424 gsl_vector_set(pGSLActualVector, 2, ActualVector.
z);
425 gsl_vector *pGSLApparentVector = gsl_vector_alloc(3);
427 ApparentTelescopeDirectionVector.
x = gsl_vector_get(pGSLApparentVector, 0);
428 ApparentTelescopeDirectionVector.
y = gsl_vector_get(pGSLApparentVector, 1);
429 ApparentTelescopeDirectionVector.
z = gsl_vector_get(pGSLApparentVector, 2);
430 ApparentTelescopeDirectionVector.
Normalise();
431 gsl_vector_free(pGSLActualVector);
432 gsl_vector_free(pGSLApparentVector);
450 gsl_matrix *pTransform;
451 gsl_matrix *pComputedTransform =
nullptr;
457 #ifdef CONVEX_HULL_DEBUGGING
460 if (
nullptr != CurrentFace)
464 #ifdef CONVEX_HULL_DEBUGGING
471 #ifdef CONVEX_HULL_DEBUGGING
472 ASSDEBUGF(
"Celestial to telescope - Ignoring actual face %d", ActualFaces);
477 #ifdef CONVEX_HULL_DEBUGGING
478 ASSDEBUGF(
"Celestial to telescope - Processing actual face %d v1 %d v2 %d v3 %d", ActualFaces,
488 CurrentFace = CurrentFace->
next;
494 std::map<double, const AlignmentDatabaseEntry *> NearestMap;
495 for (InMemoryDatabase::AlignmentDatabaseType::const_iterator Itr = SyncPoints.begin();
496 Itr != SyncPoints.end(); Itr++)
512 NearestMap[(ActualDirectionCosine - ActualVector).Length()] = &(*Itr);
515 std::map<double, const AlignmentDatabaseEntry *>::const_iterator Nearest = NearestMap.begin();
555 pComputedTransform = gsl_matrix_alloc(3, 3);
559 pTransform = pComputedTransform;
562 pTransform = CurrentFace->
pMatrix;
568 gsl_vector *pGSLActualVector = gsl_vector_alloc(3);
569 gsl_vector_set(pGSLActualVector, 0, ActualVector.
x);
570 gsl_vector_set(pGSLActualVector, 1, ActualVector.
y);
571 gsl_vector_set(pGSLActualVector, 2, ActualVector.
z);
572 gsl_vector *pGSLApparentVector = gsl_vector_alloc(3);
574 ApparentTelescopeDirectionVector.
x = gsl_vector_get(pGSLApparentVector, 0);
575 ApparentTelescopeDirectionVector.
y = gsl_vector_get(pGSLApparentVector, 1);
576 ApparentTelescopeDirectionVector.
z = gsl_vector_get(pGSLApparentVector, 2);
577 ApparentTelescopeDirectionVector.
Normalise();
578 gsl_vector_free(pGSLActualVector);
579 gsl_vector_free(pGSLApparentVector);
580 if (
nullptr != pComputedTransform)
581 gsl_matrix_free(pComputedTransform);
594 double &RightAscension,
double &Declination)
608 ASSDEBUG(
"No database or no position in database");
612 switch (SyncPoints.size())
624 ASSDEBUGF(
"ApparentVector x %lf y %lf z %lf", ApparentTelescopeDirectionVector.
x,
625 ApparentTelescopeDirectionVector.
y, ApparentTelescopeDirectionVector.
z);
648 gsl_vector *pGSLApparentVector = gsl_vector_alloc(3);
649 gsl_vector_set(pGSLApparentVector, 0, ApparentTelescopeDirectionVector.
x);
650 gsl_vector_set(pGSLApparentVector, 1, ApparentTelescopeDirectionVector.
y);
651 gsl_vector_set(pGSLApparentVector, 2, ApparentTelescopeDirectionVector.
z);
652 gsl_vector *pGSLActualVector = gsl_vector_alloc(3);
655 Dump3(
"ApparentVector", pGSLApparentVector);
656 Dump3(
"ActualVector", pGSLActualVector);
659 ActualTelescopeDirectionVector.
x = gsl_vector_get(pGSLActualVector, 0);
660 ActualTelescopeDirectionVector.
y = gsl_vector_get(pGSLActualVector, 1);
661 ActualTelescopeDirectionVector.
z = gsl_vector_get(pGSLActualVector, 2);
662 ActualTelescopeDirectionVector.
Normalise();
674 gsl_vector_free(pGSLActualVector);
675 gsl_vector_free(pGSLApparentVector);
681 gsl_matrix *pTransform;
682 gsl_matrix *pComputedTransform =
nullptr;
688 #ifdef CONVEX_HULL_DEBUGGING
689 int ApparentFaces = 0;
691 if (
nullptr != CurrentFace)
695 #ifdef CONVEX_HULL_DEBUGGING
702 #ifdef CONVEX_HULL_DEBUGGING
703 ASSDEBUGF(
"Celestial to telescope - Ignoring apparent face %d", ApparentFaces);
708 #ifdef CONVEX_HULL_DEBUGGING
709 ASSDEBUGF(
"TelescopeToCelestial - Processing apparent face %d v1 %d v2 %d v3 %d", ApparentFaces,
714 SyncPoints[CurrentFace->
vertex[0]->
vnum - 1].TelescopeDirection,
715 SyncPoints[CurrentFace->
vertex[1]->
vnum - 1].TelescopeDirection,
716 SyncPoints[CurrentFace->
vertex[2]->
vnum - 1].TelescopeDirection))
719 CurrentFace = CurrentFace->
next;
725 std::map<double, const AlignmentDatabaseEntry *> NearestMap;
726 for (InMemoryDatabase::AlignmentDatabaseType::const_iterator Itr = SyncPoints.begin();
727 Itr != SyncPoints.end(); Itr++)
729 NearestMap[((*Itr).TelescopeDirection - ApparentTelescopeDirectionVector).Length()] = &(*Itr);
732 std::map<double, const AlignmentDatabaseEntry *>::const_iterator Nearest = NearestMap.begin();
771 pComputedTransform = gsl_matrix_alloc(3, 3);
774 ActualDirectionCosine2, ActualDirectionCosine3, pComputedTransform,
776 pTransform = pComputedTransform;
779 pTransform = CurrentFace->
pMatrix;
785 gsl_vector *pGSLApparentVector = gsl_vector_alloc(3);
786 gsl_vector_set(pGSLApparentVector, 0, ApparentTelescopeDirectionVector.
x);
787 gsl_vector_set(pGSLApparentVector, 1, ApparentTelescopeDirectionVector.
y);
788 gsl_vector_set(pGSLApparentVector, 2, ApparentTelescopeDirectionVector.
z);
789 gsl_vector *pGSLActualVector = gsl_vector_alloc(3);
792 ActualTelescopeDirectionVector.
x = gsl_vector_get(pGSLActualVector, 0);
793 ActualTelescopeDirectionVector.
y = gsl_vector_get(pGSLActualVector, 1);
794 ActualTelescopeDirectionVector.
z = gsl_vector_get(pGSLActualVector, 2);
795 ActualTelescopeDirectionVector.
Normalise();
808 gsl_vector_free(pGSLActualVector);
809 gsl_vector_free(pGSLApparentVector);
810 if (
nullptr != pComputedTransform)
811 gsl_matrix_free(pComputedTransform);
824 ASSDEBUGF(
"%lf %lf %lf", gsl_vector_get(pVector, 0), gsl_vector_get(pVector, 1), gsl_vector_get(pVector, 2));
830 ASSDEBUGF(
"Row 0 %lf %lf %lf", gsl_matrix_get(pMatrix, 0, 0), gsl_matrix_get(pMatrix, 0, 1),
831 gsl_matrix_get(pMatrix, 0, 2));
832 ASSDEBUGF(
"Row 1 %lf %lf %lf", gsl_matrix_get(pMatrix, 1, 0), gsl_matrix_get(pMatrix, 1, 1),
833 gsl_matrix_get(pMatrix, 1, 2));
834 ASSDEBUGF(
"Row 2 %lf %lf %lf", gsl_matrix_get(pMatrix, 2, 0), gsl_matrix_get(pMatrix, 2, 1),
835 gsl_matrix_get(pMatrix, 2, 2));
841 gsl_permutation *pPermutation = gsl_permutation_alloc(3);
842 gsl_matrix *pDecomp = gsl_matrix_alloc(3, 3);
846 gsl_matrix_memcpy(pDecomp, pMatrix);
848 gsl_linalg_LU_decomp(pDecomp, pPermutation, &Signum);
850 Determinant = gsl_linalg_LU_det(pDecomp, Signum);
852 gsl_matrix_free(pDecomp);
853 gsl_permutation_free(pPermutation);
862 gsl_permutation *pPermutation = gsl_permutation_alloc(3);
863 gsl_matrix *pDecomp = gsl_matrix_alloc(3, 3);
866 gsl_matrix_memcpy(pDecomp, pInput);
868 gsl_linalg_LU_decomp(pDecomp, pPermutation, &Signum);
871 if (0 == gsl_linalg_LU_det(pDecomp, Signum))
876 gsl_linalg_LU_invert(pDecomp, pPermutation, pInversion);
878 gsl_matrix_free(pDecomp);
879 gsl_permutation_free(pPermutation);
889 gsl_matrix_set_zero(pC);
891 gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, pA, pB, 0.0, pC);
899 gsl_vector_set_zero(pC);
901 gsl_blas_dgemv(CblasNoTrans, 1.0, pA, pB, 0.0, pC);
915 double Determinant = Edge1 ^ P;
916 double InverseDeterminant = 1.0 / Determinant;
920 if ((Determinant > -std::numeric_limits<double>::epsilon()) &&
921 (Determinant < std::numeric_limits<double>::epsilon()))
928 double u = (T ^ P) * InverseDeterminant;
930 if (u < 0.0 || u > 1.0)
938 double v = (Ray ^ Q) * InverseDeterminant;
940 if (v < 0.0 || u + v > 1.0)
944 double t = (Edge2 ^ Q) * InverseDeterminant;
946 if (t > std::numeric_limits<double>::epsilon())
#define ASSDEBUGF(msg,...)
bool MatrixInvert3x3(gsl_matrix *pInput, gsl_matrix *pInversion)
Calculate the inverse of the supplied matrix.
gsl_matrix * pApparentToActualTransform
gsl_matrix * pActualToApparentTransform
double Matrix3x3Determinant(gsl_matrix *pMatrix)
Caluclate the determinant of the supplied matrix.
void MatrixMatrixMultiply(gsl_matrix *pA, gsl_matrix *pB, gsl_matrix *pC)
Multiply matrix A by matrix B and put the result in C.
void Dump3(const char *Label, gsl_vector *pVector)
Print out a 3 vector to debug.
virtual void CalculateTransformMatrices(const TelescopeDirectionVector &Alpha1, const TelescopeDirectionVector &Alpha2, const TelescopeDirectionVector &Alpha3, const TelescopeDirectionVector &Beta1, const TelescopeDirectionVector &Beta2, const TelescopeDirectionVector &Beta3, gsl_matrix *pAlphaToBeta, gsl_matrix *pBetaToAlpha)=0
Calculate tranformation matrices from the supplied vectors.
std::vector< TelescopeDirectionVector > ActualDirectionCosines
ConvexHull ActualConvexHull
virtual bool TransformTelescopeToCelestial(const TelescopeDirectionVector &ApparentTelescopeDirectionVector, double &RightAscension, double &Declination)
Override for the base class virtual function.
void MatrixVectorMultiply(gsl_matrix *pA, gsl_vector *pB, gsl_vector *pC)
Multiply matrix A by vector B and put the result in vector C.
BasicMathPlugin()
Default constructor.
virtual ~BasicMathPlugin()
Virtual destructor.
bool RayTriangleIntersection(TelescopeDirectionVector &Ray, TelescopeDirectionVector &TriangleVertex1, TelescopeDirectionVector &TriangleVertex2, TelescopeDirectionVector &TriangleVertex3)
Test if a ray intersects a triangle in 3d space.
virtual bool TransformCelestialToTelescope(const double RightAscension, const double Declination, double JulianOffset, TelescopeDirectionVector &ApparentTelescopeDirectionVector)
Override for the base class virtual function.
ConvexHull ApparentConvexHull
virtual bool Initialise(InMemoryDatabase *pInMemoryDatabase)
Override for the base class virtual function.
void Dump3x3(const char *Label, gsl_matrix *pMatrix)
Print out a 3x3 matrix to debug.
void EdgeOrderOnFaces()
EdgeOrderOnFaces: puts e0 between v0 and v1, e1 between v1 and v2, e2 between v2 and v0 on each face....
void Reset()
Frees the vertices edges and faces lists and resets the debug and check flags.
void MakeNewVertex(double x, double y, double z, int VertexId)
Makes a vertex from the supplied information and adds it to the vertices list.
void PrintOut(const char *FileName, tVertex v)
Prints vertices, edges and faces to the standard error output.
bool DoubleTriangle()
DoubleTriangle builds the initial double triangle. It first finds 3 noncollinear points and makes two...
void ConstructHull()
ConstructHull adds the vertices to the hull one at a time. The hull vertices are those in the list ma...
void PrintObj(const char *FileName="chull.obj")
Outputs the faces in Lightwave obj format for 3d viewing. The files chull.obj and chull....
This class provides the driver side API to the in memory alignment database.
AlignmentDatabaseType & GetAlignmentDatabase()
Get a reference to the in memory database.
std::vector< AlignmentDatabaseEntry > AlignmentDatabaseType
bool GetDatabaseReferencePosition(IGeographicCoordinates &Position)
Get the database reference position.
InMemoryDatabase * pInMemoryDatabase
virtual bool Initialise(InMemoryDatabase *pInMemoryDatabase)
Initialise or re-initialise the math plugin. Re-reading the in memory database as necessary.
MountAlignment_t ApproximateMountAlignment
Describe the approximate alignment of the mount. This information is normally used in a one star alig...
const TelescopeDirectionVector TelescopeDirectionVectorFromAltitudeAzimuth(INDI::IHorizontalCoordinates HorizontalCoordinates)
Calculates a normalised direction vector from the supplied altitude and azimuth.
void EquatorialCoordinatesFromTelescopeDirectionVector(const TelescopeDirectionVector TelescopeDirectionVector, INDI::IEquatorialCoordinates &EquatorialCoordinates)
Calculates equatorial coordinates from the supplied telescope direction vector and declination.
const TelescopeDirectionVector TelescopeDirectionVectorFromEquatorialCoordinates(INDI::IEquatorialCoordinates EquatorialCoordinates)
Calculates a telescope direction vector from the supplied equatorial coordinates.
void AltitudeAzimuthFromTelescopeDirectionVector(const TelescopeDirectionVector TelescopeDirectionVector, INDI::IHorizontalCoordinates &HorizontalCoordinates)
Calculates an altitude and azimuth from the supplied normalised direction vector and declination.
Implementations for common driver routines.
Namespace to encapsulate INDI client, drivers, and mediator classes.
void EquatorialToHorizontal(IEquatorialCoordinates *object, IGeographicCoordinates *observer, double JD, IHorizontalCoordinates *position)
EquatorialToHorizontal Calculate horizontal coordinates from equatorial coordinates.
void HorizontalToEquatorial(IHorizontalCoordinates *object, IGeographicCoordinates *observer, double JD, IEquatorialCoordinates *position)
HorizontalToEquatorial Calculate Equatorial EOD Coordinates from horizontal coordinates.
Entry in the in memory alignment database.
double RightAscension
Right ascension in decimal hours. N.B. libnova works in decimal degrees so conversion is always neede...
double ObservationJulianDate
TelescopeDirectionVector TelescopeDirection
Normalised vector giving telescope pointing direction. This is referred to elsewhere as the "apparent...
double Declination
Declination in decimal degrees.
Holds a nomalised direction vector (direction cosines)
void Normalise()
Normalise the vector.