MbDMath correction

.github/workflows/build.yml

@@ -0,0 +1,68 @@
+name: cmake-multiplatform-build
+on:
+  push:
+  pull_request:
+    branches: [ "main" ]
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.ref }}
+  cancel-in-progress: true
+
+
+jobs:
+  build:
+    runs-on: ${{ matrix.os }}
+
+    strategy:
+      fail-fast: false
+
+      matrix:
+        os: [ubuntu-latest, windows-latest, macos-latest]
+        build_type: [Release]
+        c_compiler: [gcc, clang, cl]
+        include:
+          - os: windows-latest
+            c_compiler: cl
+            cpp_compiler: cl
+          - os: ubuntu-latest
+            c_compiler: gcc
+            cpp_compiler: g++
+          - os: ubuntu-latest
+            c_compiler: clang
+            cpp_compiler: clang++
+          - os: macos-latest
+            c_compiler: clang
+            cpp_compiler: clang++
+        exclude:
+          - os: windows-latest
+            c_compiler: gcc
+          - os: windows-latest
+            c_compiler: clang
+          - os: ubuntu-latest
+            c_compiler: cl
+          - os: macos-latest
+            c_compiler: gcc
+          - os: macos-latest
+            c_compiler: cl
+
+    steps:
+    - uses: actions/checkout@v3
+
+    - name: Set reusable strings
+      id: strings
+      shell: bash
+      run: |
+        echo "build-output-dir=${{ github.workspace }}/build" >> "$GITHUB_OUTPUT"
+
+    - name: Configure CMake
+      run: >
+        cmake -B ${{ steps.strings.outputs.build-output-dir }}
+        -DCMAKE_CXX_COMPILER=${{ matrix.cpp_compiler }}
+        -DCMAKE_C_COMPILER=${{ matrix.c_compiler }}
+        -DCMAKE_BUILD_TYPE=${{ matrix.build_type }}
+        -S ${{ github.workspace }}
+
+    - name: Build
+      # Note that --config is needed because the default Windows generator is a multi-config generator (Visual Studio generator).
+      run: cmake --build ${{ steps.strings.outputs.build-output-dir }} --config ${{ matrix.build_type }} -j 2
+

.github/workflows/push-freecad.yml

@@ -0,0 +1,25 @@
+name: push-freecad-submodule
+on:
+  push:
+
+jobs:
+  push-submodule:
+    runs-on: ubuntu-latest
+    steps:
+    - uses: actions/checkout@v3
+      with:
+        repository: Ondsel-Development/FreeCAD
+        ref: OndselSolver-testing
+        submodules: true
+        token: ${{ secrets.TOKEN }}
+
+    - run: |
+        cd src/3rdParty/OndselSolver
+        git fetch origin $GITHUB_REF
+        git checkout $GITHUB_SHA
+    
+    - run: |
+        git config --global user.name 'github-actions[bot]'
+        git config --global user.email '41898282+github-actions[bot]@users.noreply.github.com'
+        git commit -am "Update OndselSolver submodule to $GITHUB_SHA"
+        git push

.gitignore

@@ -35,5 +35,6 @@
 x64/
 *.bak
 
+build
 cmake-build-debug
 .idea

CMakeLists.txt

@@ -7,7 +7,30 @@ set(CMAKE_CXX_STANDARD_REQUIRED True)
 
 include(GNUInstallDirs)
 
-set(ONDSEL_EXPORT
+if( UNIX )
+    set( ONDSELSOLVER_BUILD_SHARED ON )
+ELSEIF ( APPLE )
+    set( ONDSELSOLVER_BUILD_SHARED ON )
+ELSE()
+    set( ONDSELSOLVER_BUILD_SHARED OFF )
+ENDIF ()
+
+if ( ONDSELSOLVER_BUILD_SHARED )
+    message( STATUS "[OndselSolver] Building shared library" )
+    add_library(OndselSolver SHARED)
+else()
+    message( STATUS "[OndselSolver] Building static library" )
+    add_library(OndselSolver STATIC)
+endif()
+
+
+set(ONDSELSOLVER_SRC
+        OndselSolver/Array.cpp
+        OndselSolver/FullVector.cpp
+        OndselSolver/RowTypeMatrix.cpp
+        OndselSolver/FullColumn.cpp
+        OndselSolver/FullMatrix.cpp
+        OndselSolver/FullRow.cpp
         OndselSolver/Abs.cpp
         OndselSolver/AbsConstraint.cpp
         OndselSolver/AccICKineNewtonRaphson.cpp
@@ -20,9 +43,9 @@ set(ONDSEL_EXPORT
         OndselSolver/AngleZIeqcJeqc.cpp
         OndselSolver/AnyGeneralSpline.cpp
         OndselSolver/AnyPosICNewtonRaphson.cpp
+        OndselSolver/ArcSine.cpp
         OndselSolver/ArcTan.cpp
         OndselSolver/ArcTan2.cpp
-        OndselSolver/Array.cpp
         OndselSolver/ASMTAnimationParameters.cpp
         OndselSolver/ASMTAssembly.cpp
         OndselSolver/ASMTConstantGravity.cpp
@@ -38,7 +61,9 @@ set(ONDSEL_EXPORT
         OndselSolver/ASMTKinematicIJ.cpp
         OndselSolver/ASMTMarker.cpp
         OndselSolver/ASMTMotion.cpp
+        OndselSolver/ASMTNoRotationJoint.cpp
         OndselSolver/ASMTPart.cpp
+        OndselSolver/ASMTPointInLineJoint.cpp
         OndselSolver/ASMTPointInPlaneJoint.cpp
         OndselSolver/ASMTPrincipalMassMarker.cpp
         OndselSolver/ASMTRefCurve.cpp
@@ -52,8 +77,8 @@ set(ONDSEL_EXPORT
         OndselSolver/ASMTSpatialItem.cpp
         OndselSolver/ASMTSphericalJoint.cpp
         OndselSolver/ASMTTime.cpp
-        OndselSolver/ASMTTranslationalMotion.cpp
         OndselSolver/ASMTTranslationalJoint.cpp
+        OndselSolver/ASMTTranslationalMotion.cpp
         OndselSolver/ASMTUniversalJoint.cpp
         OndselSolver/AtPointConstraintIJ.cpp
         OndselSolver/AtPointConstraintIqcJc.cpp
@@ -78,24 +103,9 @@ set(ONDSEL_EXPORT
         OndselSolver/CREATE.cpp
         OndselSolver/CylindricalJoint.cpp
         OndselSolver/CylSphJoint.cpp
-        OndselSolver/DifferentiatedGeneralSpline.cpp
-        OndselSolver/EndFrameqct2.cpp
-        OndselSolver/EulerAngles.cpp
-        OndselSolver/EulerAnglesDDot.cpp
-        OndselSolver/EulerAnglesDot.cpp
-        OndselSolver/Exponential.cpp
-        OndselSolver/ExternalSystem.cpp
-        OndselSolver/FunctionFromData.cpp
-        OndselSolver/FunctionXcParameter.cpp
-        OndselSolver/FunctionXY.cpp
-        OndselSolver/GeneralSpline.cpp
-        OndselSolver/Ln.cpp
-        OndselSolver/Log10.cpp
-        OndselSolver/LogN.cpp
-        OndselSolver/Negative.cpp
-        OndselSolver/PosVelAccData.cpp
         OndselSolver/DiagonalMatrix.cpp
         OndselSolver/DifferenceOperator.cpp
+        OndselSolver/DifferentiatedGeneralSpline.cpp
         OndselSolver/DirectionCosineConstraintIJ.cpp
         OndselSolver/DirectionCosineConstraintIqcJc.cpp
         OndselSolver/DirectionCosineConstraintIqcJqc.cpp
@@ -134,9 +144,14 @@ set(ONDSEL_EXPORT
         OndselSolver/DistxyIeqcJec.cpp
         OndselSolver/DistxyIeqcJeqc.cpp
         OndselSolver/DistxyIeqctJeqc.cpp
+        OndselSolver/EigenDecomposition.cpp
         OndselSolver/EndFramec.cpp
         OndselSolver/EndFrameqc.cpp
         OndselSolver/EndFrameqct.cpp
+        OndselSolver/EndFrameqct2.cpp
+        OndselSolver/EulerAngles.cpp
+        OndselSolver/EulerAnglesDDot.cpp
+        OndselSolver/EulerAnglesDot.cpp
         OndselSolver/EulerAngleszxz.cpp
         OndselSolver/EulerAngleszxzDDot.cpp
         OndselSolver/EulerAngleszxzDot.cpp
@@ -145,17 +160,19 @@ set(ONDSEL_EXPORT
         OndselSolver/EulerParameters.cpp
         OndselSolver/EulerParametersDDot.cpp
         OndselSolver/EulerParametersDot.cpp
+        OndselSolver/Exponential.cpp
         OndselSolver/ExpressionX.cpp
+        OndselSolver/ExternalSystem.cpp
         OndselSolver/FixedJoint.cpp
         OndselSolver/ForceTorqueData.cpp
         OndselSolver/ForceTorqueItem.cpp
-        OndselSolver/FullColumn.cpp
-        OndselSolver/FullMatrix.cpp
         OndselSolver/FullMotion.cpp
-        OndselSolver/FullRow.cpp
         OndselSolver/Function.cpp
+        OndselSolver/FunctionFromData.cpp
         OndselSolver/FunctionWithManyArgs.cpp
         OndselSolver/FunctionX.cpp
+        OndselSolver/FunctionXcParameter.cpp
+        OndselSolver/FunctionXY.cpp
         OndselSolver/GearConstraintIJ.cpp
         OndselSolver/GearConstraintIqcJc.cpp
         OndselSolver/GearConstraintIqcJqc.cpp
@@ -163,6 +180,7 @@ set(ONDSEL_EXPORT
         OndselSolver/GEFullMat.cpp
         OndselSolver/GEFullMatFullPv.cpp
         OndselSolver/GEFullMatParPv.cpp
+        OndselSolver/GeneralSpline.cpp
         OndselSolver/GESpMat.cpp
         OndselSolver/GESpMatFullPv.cpp
         OndselSolver/GESpMatFullPvPosIC.cpp
@@ -180,6 +198,7 @@ set(ONDSEL_EXPORT
         OndselSolver/Joint.cpp
         OndselSolver/KineIntegrator.cpp
         OndselSolver/KinematicIeJe.cpp
+        OndselSolver/LDUFullMat.cpp
         OndselSolver/LDUFullMatParPv.cpp
         OndselSolver/LDUSpMat.cpp
         OndselSolver/LDUSpMatParPv.cpp
@@ -187,6 +206,9 @@ set(ONDSEL_EXPORT
         OndselSolver/LDUSpMatParPvPrecise.cpp
         OndselSolver/LinearMultiStepMethod.cpp
         OndselSolver/LineInPlaneJoint.cpp
+        OndselSolver/Ln.cpp
+        OndselSolver/Log10.cpp
+        OndselSolver/LogN.cpp
         OndselSolver/MarkerFrame.cpp
         OndselSolver/MatrixDecomposition.cpp
         OndselSolver/MatrixGaussElimination.cpp
@@ -194,6 +216,21 @@ set(ONDSEL_EXPORT
         OndselSolver/MatrixSolver.cpp
         OndselSolver/MaximumIterationError.cpp
         OndselSolver/MbdMath.cpp
+        OndselSolver/MBDynBlock.cpp
+        OndselSolver/MBDynBody.cpp
+        OndselSolver/MBDynControlData.cpp
+        OndselSolver/MBDynData.cpp
+        OndselSolver/MBDynElement.cpp
+        OndselSolver/MBDynInitialValue.cpp
+        OndselSolver/MBDynItem.cpp
+        OndselSolver/MBDynJoint.cpp
+        OndselSolver/MBDynMarker.cpp
+        OndselSolver/MBDynNode.cpp
+        OndselSolver/MBDynReference.cpp
+        OndselSolver/MBDynStructural.cpp
+        OndselSolver/MBDynSystem.cpp
+        OndselSolver/MomentOfInertiaSolver.cpp
+        OndselSolver/Negative.cpp
         OndselSolver/NewtonRaphson.cpp
         OndselSolver/NewtonRaphsonError.cpp
         OndselSolver/NoRotationJoint.cpp
@@ -214,6 +251,7 @@ set(ONDSEL_EXPORT
         OndselSolver/PosICNewtonRaphson.cpp
         OndselSolver/PosKineNewtonRaphson.cpp
         OndselSolver/PosNewtonRaphson.cpp
+        OndselSolver/PosVelAccData.cpp
         OndselSolver/Power.cpp
         OndselSolver/PrescribedMotion.cpp
         OndselSolver/Product.cpp
@@ -226,7 +264,6 @@ set(ONDSEL_EXPORT
         OndselSolver/RedundantConstraint.cpp
         OndselSolver/RevCylJoint.cpp
         OndselSolver/RevoluteJoint.cpp
-        OndselSolver/RowTypeMatrix.cpp
         OndselSolver/ScalarNewtonRaphson.cpp
         OndselSolver/ScrewConstraintIJ.cpp
         OndselSolver/ScrewConstraintIqcJc.cpp
@@ -243,7 +280,6 @@ set(ONDSEL_EXPORT
         OndselSolver/SphericalJoint.cpp
         OndselSolver/SphSphJoint.cpp
         OndselSolver/StableBackwardDifference.cpp
-        OndselSolver/LDUFullMat.cpp
         OndselSolver/StateData.cpp
         OndselSolver/Sum.cpp
         OndselSolver/Symbolic.cpp
@@ -265,7 +301,6 @@ set(ONDSEL_EXPORT
         OndselSolver/UniversalJoint.cpp
         OndselSolver/UserFunction.cpp
         OndselSolver/Variable.cpp
-        OndselSolver/FullVector.cpp
         OndselSolver/VectorNewtonRaphson.cpp
         OndselSolver/VelICKineSolver.cpp
         OndselSolver/VelICSolver.cpp
@@ -273,25 +308,312 @@ set(ONDSEL_EXPORT
         OndselSolver/VelSolver.cpp
         OndselSolver/ZRotation.cpp
         OndselSolver/ZTranslation.cpp
-        OndselSolver/ASMTAssembly.cpp
 )
 
-add_executable(ondsel-solver OndselSolver/OndselSolver.cpp)
-
-add_library(OndselSolver
-        STATIC ${ONDSEL_EXPORT}
+set(ONDSELSOLVER_HEADERS
+        OndselSolver/Array.h
+        OndselSolver/FullVector.h
+        OndselSolver/RowTypeMatrix.h
+        OndselSolver/FullRow.h
+        OndselSolver/FullColumn.h
+        OndselSolver/FullMatrix.h
+        OndselSolver/Abs.h
+        OndselSolver/AbsConstraint.h
+        OndselSolver/AccICKineNewtonRaphson.h
+        OndselSolver/AccICNewtonRaphson.h
+        OndselSolver/AccKineNewtonRaphson.h
+        OndselSolver/AccNewtonRaphson.h
+        OndselSolver/AngleJoint.h
+        OndselSolver/AngleZIecJec.h
+        OndselSolver/AngleZIeqcJec.h
+        OndselSolver/AngleZIeqcJeqc.h
+        OndselSolver/AnyGeneralSpline.h
+        OndselSolver/AnyPosICNewtonRaphson.h
+        OndselSolver/APIExport.h
+        OndselSolver/ArcSine.h
+        OndselSolver/ArcTan.h
+        OndselSolver/ArcTan2.h
+        OndselSolver/ASMTAnimationParameters.h
+        OndselSolver/ASMTAssembly.h
+        OndselSolver/ASMTConstantGravity.h
+        OndselSolver/ASMTConstraintSet.h
+        OndselSolver/ASMTCylindricalJoint.h
+        OndselSolver/ASMTExtrusion.h
+        OndselSolver/ASMTFixedJoint.h
+        OndselSolver/ASMTForceTorque.h
+        OndselSolver/ASMTGeneralMotion.h
+        OndselSolver/ASMTItem.h
+        OndselSolver/ASMTItemIJ.h
+        OndselSolver/ASMTJoint.h
+        OndselSolver/ASMTKinematicIJ.h
+        OndselSolver/ASMTMarker.h
+        OndselSolver/ASMTMotion.h
+        OndselSolver/ASMTNoRotationJoint.h
+        OndselSolver/ASMTPart.h
+        OndselSolver/ASMTPointInLineJoint.h
+        OndselSolver/ASMTPointInPlaneJoint.h
+        OndselSolver/ASMTPrincipalMassMarker.h
+        OndselSolver/ASMTRefCurve.h
+        OndselSolver/ASMTRefItem.h
+        OndselSolver/ASMTRefPoint.h
+        OndselSolver/ASMTRefSurface.h
+        OndselSolver/ASMTRevoluteJoint.h
+        OndselSolver/ASMTRotationalMotion.h
+        OndselSolver/ASMTSimulationParameters.h
+        OndselSolver/ASMTSpatialContainer.h
+        OndselSolver/ASMTSpatialItem.h
+        OndselSolver/ASMTSphericalJoint.h
+        OndselSolver/ASMTTime.h
+        OndselSolver/ASMTTranslationalJoint.h
+        OndselSolver/ASMTTranslationalMotion.h
+        OndselSolver/ASMTUniversalJoint.h
+        OndselSolver/AtPointConstraintIJ.h
+        OndselSolver/AtPointConstraintIqcJc.h
+        OndselSolver/AtPointConstraintIqcJqc.h
+        OndselSolver/AtPointConstraintIqctJqc.h
+        OndselSolver/AtPointJoint.h
+        OndselSolver/BasicIntegrator.h
+        OndselSolver/BasicQuasiIntegrator.h
+        OndselSolver/BasicUserFunction.h
+        OndselSolver/CADSystem.h
+        OndselSolver/CartesianFrame.h
+        OndselSolver/CompoundJoint.h
+        OndselSolver/Constant.h
+        OndselSolver/ConstantGravity.h
+        OndselSolver/ConstantVelocityJoint.h
+        OndselSolver/Constraint.h
+        OndselSolver/ConstraintIJ.h
+        OndselSolver/ConstVelConstraintIJ.h
+        OndselSolver/ConstVelConstraintIqcJc.h
+        OndselSolver/ConstVelConstraintIqcJqc.h
+        OndselSolver/corecrt_math_defines.h
+        OndselSolver/Cosine.h
+        OndselSolver/CREATE.h
+        OndselSolver/CylindricalJoint.h
+        OndselSolver/CylSphJoint.h
+        OndselSolver/DiagonalMatrix.h
+        OndselSolver/DifferenceOperator.h
+        OndselSolver/DifferentiatedGeneralSpline.h
+        OndselSolver/DirectionCosineConstraintIJ.h
+        OndselSolver/DirectionCosineConstraintIqcJc.h
+        OndselSolver/DirectionCosineConstraintIqcJqc.h
+        OndselSolver/DirectionCosineConstraintIqctJqc.h
+        OndselSolver/DirectionCosineIecJec.h
+        OndselSolver/DirectionCosineIeqcJec.h
+        OndselSolver/DirectionCosineIeqcJeqc.h
+        OndselSolver/DirectionCosineIeqctJeqc.h
+        OndselSolver/DiscontinuityError.h
+        OndselSolver/DispCompIecJecIe.h
+        OndselSolver/DispCompIecJecKec.h
+        OndselSolver/DispCompIecJecKeqc.h
+        OndselSolver/DispCompIecJecO.h
+        OndselSolver/DispCompIeqcJecIe.h
+        OndselSolver/DispCompIeqcJecKeqc.h
+        OndselSolver/DispCompIeqcJecO.h
+        OndselSolver/DispCompIeqcJeqcIe.h
+        OndselSolver/DispCompIeqcJeqcKeqc.h
+        OndselSolver/DispCompIeqcJeqcKeqct.h
+        OndselSolver/DispCompIeqcJeqcO.h
+        OndselSolver/DispCompIeqctJeqcIe.h
+        OndselSolver/DispCompIeqctJeqcKeqct.h
+        OndselSolver/DispCompIeqctJeqcO.h
+        OndselSolver/DistanceConstraintIJ.h
+        OndselSolver/DistanceConstraintIqcJc.h
+        OndselSolver/DistanceConstraintIqcJqc.h
+        OndselSolver/DistanceConstraintIqctJqc.h
+        OndselSolver/DistancexyConstraintIJ.h
+        OndselSolver/DistancexyConstraintIqcJc.h
+        OndselSolver/DistancexyConstraintIqcJqc.h
+        OndselSolver/DistIecJec.h
+        OndselSolver/DistIeqcJec.h
+        OndselSolver/DistIeqcJeqc.h
+        OndselSolver/DistIeqctJeqc.h
+        OndselSolver/DistxyIecJec.h
+        OndselSolver/DistxyIeqcJec.h
+        OndselSolver/DistxyIeqcJeqc.h
+        OndselSolver/DistxyIeqctJeqc.h
+        OndselSolver/EigenDecomposition.h
+        OndselSolver/EndFramec.h
+        OndselSolver/EndFrameqc.h
+        OndselSolver/EndFrameqct.h
+        OndselSolver/EndFrameqct2.h
+        OndselSolver/enum.h
+        OndselSolver/EulerAngles.h
+        OndselSolver/EulerAnglesDDot.h
+        OndselSolver/EulerAnglesDot.h
+        OndselSolver/EulerAngleszxz.h
+        OndselSolver/EulerAngleszxzDDot.h
+        OndselSolver/EulerAngleszxzDot.h
+        OndselSolver/EulerArray.h
+        OndselSolver/EulerConstraint.h
+        OndselSolver/EulerParameters.h
+        OndselSolver/EulerParametersDDot.h
+        OndselSolver/EulerParametersDot.h
+        OndselSolver/Exponential.h
+        OndselSolver/ExpressionX.h
+        OndselSolver/ExternalSystem.h
+        OndselSolver/FixedJoint.h
+        OndselSolver/ForceTorqueData.h
+        OndselSolver/ForceTorqueItem.h
+        OndselSolver/FullMotion.h
+        OndselSolver/FullRow.h
+        OndselSolver/Function.h
+        OndselSolver/FunctionFromData.h
+        OndselSolver/FunctionWithManyArgs.h
+        OndselSolver/FunctionX.h
+        OndselSolver/FunctionXcParameter.h
+        OndselSolver/FunctionXY.h
+        OndselSolver/GearConstraintIJ.h
+        OndselSolver/GearConstraintIqcJc.h
+        OndselSolver/GearConstraintIqcJqc.h
+        OndselSolver/GearJoint.h
+        OndselSolver/GEFullMat.h
+        OndselSolver/GEFullMatFullPv.h
+        OndselSolver/GEFullMatParPv.h
+        OndselSolver/GeneralSpline.h
+        OndselSolver/GESpMat.h
+        OndselSolver/GESpMatFullPv.h
+        OndselSolver/GESpMatFullPvPosIC.h
+        OndselSolver/GESpMatParPv.h
+        OndselSolver/GESpMatParPvMarko.h
+        OndselSolver/GESpMatParPvMarkoFast.h
+        OndselSolver/GESpMatParPvPrecise.h
+        OndselSolver/ICKineIntegrator.h
+        OndselSolver/IndependentVariable.h
+        OndselSolver/InLineJoint.h
+        OndselSolver/InPlaneJoint.h
+        OndselSolver/Integrator.h
+        OndselSolver/IntegratorInterface.h
+        OndselSolver/Item.h
+        OndselSolver/Joint.h
+        OndselSolver/KineIntegrator.h
+        OndselSolver/KinematicIeJe.h
+        OndselSolver/LDUFullMat.h
+        OndselSolver/LDUFullMatParPv.h
+        OndselSolver/LDUSpMat.h
+        OndselSolver/LDUSpMatParPv.h
+        OndselSolver/LDUSpMatParPvMarko.h
+        OndselSolver/LDUSpMatParPvPrecise.h
+        OndselSolver/LinearMultiStepMethod.h
+        OndselSolver/LineInPlaneJoint.h
+        OndselSolver/Ln.h
+        OndselSolver/Log10.h
+        OndselSolver/LogN.h
+        OndselSolver/MarkerFrame.h
+        OndselSolver/MatrixDecomposition.h
+        OndselSolver/MatrixGaussElimination.h
+        OndselSolver/MatrixLDU.h
+        OndselSolver/MatrixSolver.h
+        OndselSolver/MaximumIterationError.h
+        OndselSolver/MbdMath.h
+        OndselSolver/MBDynBlock.h
+        OndselSolver/MBDynBody.h
+        OndselSolver/MBDynControlData.h
+        OndselSolver/MBDynData.h
+        OndselSolver/MBDynElement.h
+        OndselSolver/MBDynInitialValue.h
+        OndselSolver/MBDynItem.h
+        OndselSolver/MBDynJoint.h
+        OndselSolver/MBDynMarker.h
+        OndselSolver/MBDynNode.h
+        OndselSolver/MBDynReference.h
+        OndselSolver/MBDynStructural.h
+        OndselSolver/MBDynSystem.h
+        OndselSolver/MomentOfInertiaSolver.h
+        OndselSolver/Negative.h
+        OndselSolver/NewtonRaphson.h
+        OndselSolver/NewtonRaphsonError.h
+        OndselSolver/NoRotationJoint.h
+        OndselSolver/NotKinematicError.h
+        OndselSolver/Numeric.h
+        OndselSolver/OrbitAnglezIecJec.h
+        OndselSolver/OrbitAnglezIeqcJec.h
+        OndselSolver/OrbitAnglezIeqcJeqc.h
+        OndselSolver/Orientation.h
+        OndselSolver/ParallelAxesJoint.h
+        OndselSolver/Part.h
+        OndselSolver/PartFrame.h
+        OndselSolver/PerpendicularJoint.h
+        OndselSolver/PlanarJoint.h
+        OndselSolver/PointInLineJoint.h
+        OndselSolver/PointInPlaneJoint.h
+        OndselSolver/PosICKineNewtonRaphson.h
+        OndselSolver/PosICNewtonRaphson.h
+        OndselSolver/PosKineNewtonRaphson.h
+        OndselSolver/PosNewtonRaphson.h
+        OndselSolver/PosVelAccData.h
+        OndselSolver/Power.h
+        OndselSolver/PrescribedMotion.h
+        OndselSolver/Product.h
+        OndselSolver/QuasiIntegrator.h
+        OndselSolver/RackPinConstraintIJ.h
+        OndselSolver/RackPinConstraintIqcJc.h
+        OndselSolver/RackPinConstraintIqcJqc.h
+        OndselSolver/RackPinJoint.h
+        OndselSolver/Reciprocal.h
+        OndselSolver/RedundantConstraint.h
+        OndselSolver/resource.h
+        OndselSolver/RevCylJoint.h
+        OndselSolver/RevoluteJoint.h
+        OndselSolver/ScalarNewtonRaphson.h
+        OndselSolver/ScrewConstraintIJ.h
+        OndselSolver/ScrewConstraintIqcJc.h
+        OndselSolver/ScrewConstraintIqcJqc.h
+        OndselSolver/ScrewJoint.h
+        OndselSolver/SimulationStoppingError.h
+        OndselSolver/Sine.h
+        OndselSolver/SingularMatrixError.h
+        OndselSolver/Solver.h
+        OndselSolver/SparseColumn.h
+        OndselSolver/SparseMatrix.h
+        OndselSolver/SparseRow.h
+        OndselSolver/SparseVector.h
+        OndselSolver/SphericalJoint.h
+        OndselSolver/SphSphJoint.h
+        OndselSolver/StableBackwardDifference.h
+        OndselSolver/StateData.h
+        OndselSolver/Sum.h
+        OndselSolver/Symbolic.h
+        OndselSolver/SymbolicParser.h
+        OndselSolver/SyntaxError.h
+        OndselSolver/System.h
+        OndselSolver/SystemNewtonRaphson.h
+        OndselSolver/SystemSolver.h
+        OndselSolver/Time.h
+        OndselSolver/TooManyTriesError.h
+        OndselSolver/TooSmallStepSizeError.h
+        OndselSolver/Translation.h
+        OndselSolver/TranslationalJoint.h
+        OndselSolver/TranslationConstraintIJ.h
+        OndselSolver/TranslationConstraintIqcJc.h
+        OndselSolver/TranslationConstraintIqcJqc.h
+        OndselSolver/TranslationConstraintIqctJqc.h
+        OndselSolver/Units.h
+        OndselSolver/UniversalJoint.h
+        OndselSolver/UserFunction.h
+        OndselSolver/Variable.h
+        OndselSolver/VectorNewtonRaphson.h
+        OndselSolver/VelICKineSolver.h
+        OndselSolver/VelICSolver.h
+        OndselSolver/VelKineSolver.h
+        OndselSolver/VelSolver.h
+        OndselSolver/ZRotation.h
+        OndselSolver/ZTranslation.h
 )
 
-add_subdirectory(OndselSolver)
+target_sources(OndselSolver PRIVATE
+        "${ONDSELSOLVER_SRC}"
+        "${ONDSELSOLVER_HEADERS}")
 
 set_target_properties(OndselSolver
         PROPERTIES VERSION ${PROJECT_VERSION}
         SOVERSION 1
+        PUBLIC_HEADER "${ONDSELSOLVER_HEADERS}"
 )
 
 configure_file(OndselSolver.pc.in ${CMAKE_BINARY_DIR}/OndselSolver.pc @ONLY)
 install(TARGETS OndselSolver
         LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
-        PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR})
+        PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/OndselSolver)
 install(FILES ${CMAKE_BINARY_DIR}/OndselSolver.pc
         DESTINATION ${CMAKE_INSTALL_DATAROOTDIR}/pkgconfig)

OndselSolver/ASMTAnimationParameters.h

@@ -18,7 +18,7 @@ namespace MbD {
         void parseASMT(std::vector<std::string>& lines) override;
         void storeOnLevel(std::ofstream& os, int level) override;
 
-		int nframe = 1000000, icurrent = 0, istart = 0, iend = 1000000, framesPerSecond = 30;
+		int nframe = 1000000, icurrent = 1, istart = 1, iend = 1000000, framesPerSecond = 30;
         bool isForward = true;
 
 

OndselSolver/ASMTAssembly.cpp

@@ -248,7 +248,7 @@ void MbD::ASMTAssembly::runSinglePendulum()
 	assembly->setName(str);
 	auto pos3D = std::make_shared<FullColumn<double>>(ListD{ 0, 0, 0 });
 	assembly->setPosition3D(pos3D);
-	auto rotMat = std::make_shared<FullMatrix<double>>(ListListD{
+	auto rotMat = std::make_shared<FullMatrixDouble>(ListListD{
 		{ 1, 0, 0 },
 		{ 0, 1, 0 },
 		{ 0, 0, 1 }
@@ -262,11 +262,11 @@ void MbD::ASMTAssembly::runSinglePendulum()
 	auto massMarker = std::make_shared<ASMTPrincipalMassMarker>();
 	massMarker->setMass(0.0);
 	massMarker->setDensity(0.0);
-	auto aJ = std::make_shared<DiagonalMatrix<double>>(ListD{ 0, 0, 0 });
+	auto aJ = std::make_shared<DiagonalMatrix>(ListD{ 0, 0, 0 });
 	massMarker->setMomentOfInertias(aJ);
 	pos3D = std::make_shared<FullColumn<double>>(ListD{ 0, 0, 0 });
 	massMarker->setPosition3D(pos3D);
-	rotMat = std::make_shared<FullMatrix<double>>(ListListD{
+	rotMat = std::make_shared<FullMatrixDouble>(ListListD{
 		{ 1, 0, 0 },
 		{ 0, 1, 0 },
 		{ 0, 0, 1 }
@@ -279,7 +279,7 @@ void MbD::ASMTAssembly::runSinglePendulum()
 	mkr->setName(str);
 	pos3D = std::make_shared<FullColumn<double>>(ListD{ 0, 0, 0 });
 	mkr->setPosition3D(pos3D);
-	rotMat = std::make_shared<FullMatrix<double>>(ListListD{
+	rotMat = std::make_shared<FullMatrixDouble>(ListListD{
 		{ 1, 0, 0 },
 		{ 0, 1, 0 },
 		{ 0, 0, 1 }
@@ -292,7 +292,7 @@ void MbD::ASMTAssembly::runSinglePendulum()
 	part->setName(str);
 	pos3D = std::make_shared<FullColumn<double>>(ListD{ -0.1, -0.1, -0.1 });
 	part->setPosition3D(pos3D);
-	rotMat = std::make_shared<FullMatrix<double>>(ListListD{
+	rotMat = std::make_shared<FullMatrixDouble>(ListListD{
 		{ 1, 0, 0 },
 		{ 0, 1, 0 },
 		{ 0, 0, 1 }
@@ -307,11 +307,11 @@ void MbD::ASMTAssembly::runSinglePendulum()
 	massMarker = std::make_shared<ASMTPrincipalMassMarker>();
 	massMarker->setMass(0.2);
 	massMarker->setDensity(10.0);
-	aJ = std::make_shared<DiagonalMatrix<double>>(ListD{ 8.3333333333333e-4, 0.016833333333333, 0.017333333333333 });
+	aJ = std::make_shared<DiagonalMatrix>(ListD{ 8.3333333333333e-4, 0.016833333333333, 0.017333333333333 });
 	massMarker->setMomentOfInertias(aJ);
 	pos3D = std::make_shared<FullColumn<double>>(ListD{ 0.5, 0.1, 0.05 });
 	massMarker->setPosition3D(pos3D);
-	rotMat = std::make_shared<FullMatrix<double>>(ListListD{
+	rotMat = std::make_shared<FullMatrixDouble>(ListListD{
 		{ 1, 0, 0 },
 		{ 0, 1, 0 },
 		{ 0, 0, 1 }
@@ -324,7 +324,7 @@ void MbD::ASMTAssembly::runSinglePendulum()
 	mkr->setName(str);
 	pos3D = std::make_shared<FullColumn<double>>(ListD{ 0.1, 0.1, 0.1 });
 	mkr->setPosition3D(pos3D);
-	rotMat = std::make_shared<FullMatrix<double>>(ListListD{
+	rotMat = std::make_shared<FullMatrixDouble>(ListListD{
 		{ 1, 0, 0 },
 		{ 0, 1, 0 },
 		{ 0, 0, 1 }
@@ -940,7 +940,7 @@ void MbD::ASMTAssembly::createMbD(std::shared_ptr<System> mbdSys, std::shared_pt
 void MbD::ASMTAssembly::outputFile(std::string filename)
 {
 	std::ofstream os(filename);
-	os << std::setprecision(std::numeric_limits<double>::digits10 + 1);
+	os << std::setprecision(static_cast<std::streamsize>(std::numeric_limits<double>::digits10) + 1);
 	//	try {
 	os << "OndselSolver" << std::endl;
 	storeOnLevel(os, 0);
@@ -995,9 +995,9 @@ void MbD::ASMTAssembly::initprincipalMassMarker()
 	principalMassMarker = std::make_shared<ASMTPrincipalMassMarker>();
 	principalMassMarker->mass = 0.0;
 	principalMassMarker->density = 0.0;
-	principalMassMarker->momentOfInertias = std::make_shared<DiagonalMatrix<double>>(3, 0);
+	principalMassMarker->momentOfInertias = std::make_shared<DiagonalMatrix>(3, 0);
 	//principalMassMarker->position3D = std::make_shared<FullColumn<double>>(3, 0);
-	//principalMassMarker->rotationMatrix = FullMatrix<double>>::identitysptr(3);
+	//principalMassMarker->rotationMatrix = FullMatrixDouble>::identitysptr(3);
 }
 
 std::shared_ptr<ASMTSpatialContainer> MbD::ASMTAssembly::spatialContainerAt(std::shared_ptr<ASMTAssembly> self, std::string& longname)
@@ -1144,7 +1144,7 @@ std::shared_ptr<ASMTPart> MbD::ASMTAssembly::partPartialNamed(std::string partia
 		auto fullName = prt->fullName("");
 		return fullName.find(partialName) != std::string::npos;
 		});
-	auto part = *it;
+	auto& part = *it;
 	return part;
 }
 
@@ -1212,6 +1212,7 @@ void MbD::ASMTAssembly::storeOnLevelGeneralConstraintSets(std::ofstream& os, int
 
 void MbD::ASMTAssembly::storeOnTimeSeries(std::ofstream& os)
 {
+	if (times->empty()) return;
 	os << "TimeSeries" << std::endl;
 	os << "Number\tInput\t";
 	for (int i = 1; i < times->size(); i++)

OndselSolver/ASMTAssembly.h

@@ -121,7 +121,7 @@ namespace MbD {
 		void storeOnLevelGeneralConstraintSets(std::ofstream& os, int level);
 		void storeOnTimeSeries(std::ofstream& os) override;
 
-		std::string notes;
+		std::string notes = "(Text string: '' runs: (Core.RunArray runs: #() values: #()))";
 		std::shared_ptr<std::vector<std::shared_ptr<ASMTPart>>> parts = std::make_shared<std::vector<std::shared_ptr<ASMTPart>>>();
 		std::shared_ptr<std::vector<std::shared_ptr<ASMTKinematicIJ>>> kinematicIJs = std::make_shared<std::vector<std::shared_ptr<ASMTKinematicIJ>>>();
 		std::shared_ptr<std::vector<std::shared_ptr<ASMTConstraintSet>>> constraintSets = std::make_shared<std::vector<std::shared_ptr<ASMTConstraintSet>>>();
@@ -129,11 +129,11 @@ namespace MbD {
 		std::shared_ptr<std::vector<std::shared_ptr<ASMTMotion>>> motions = std::make_shared<std::vector<std::shared_ptr<ASMTMotion>>>();
 		std::shared_ptr<std::vector<std::shared_ptr<ASMTForceTorque>>> forcesTorques = std::make_shared<std::vector<std::shared_ptr<ASMTForceTorque>>>();
 		std::shared_ptr<ASMTConstantGravity> constantGravity = std::make_shared<ASMTConstantGravity>();
-		std::shared_ptr<ASMTSimulationParameters> simulationParameters;
-		std::shared_ptr<ASMTAnimationParameters> animationParameters;
-		std::shared_ptr<std::vector<double>> times;
+		std::shared_ptr<ASMTSimulationParameters> simulationParameters = std::make_shared<ASMTSimulationParameters>();
+		std::shared_ptr<ASMTAnimationParameters> animationParameters = std::make_shared<ASMTAnimationParameters>();
+		std::shared_ptr<std::vector<double>> times = std::make_shared<std::vector<double>>();
 		std::shared_ptr<ASMTTime> asmtTime = std::make_shared<ASMTTime>();
-		std::shared_ptr<Units> mbdUnits;
+		std::shared_ptr<Units> mbdUnits = std::make_shared<Units>();
 		MBDynSystem* mbdynItem = nullptr;
 	};
 }

OndselSolver/ASMTConstraintSet.cpp

@@ -10,6 +10,7 @@
 #include "ASMTAssembly.h"
 #include "ASMTMarker.h"
 #include "Joint.h"
+#include "FullMatrix.h"
 
 using namespace MbD;
 

OndselSolver/ASMTGeneralMotion.cpp

@@ -92,21 +92,21 @@ void MbD::ASMTGeneralMotion::createMbD(std::shared_ptr<System> mbdSys, std::shar
 	userFunc = std::make_shared<BasicUserFunction>(rIJI->at(0), 1.0);
 	parser->parseUserFunction(userFunc);
 	auto geoX = parser->stack->top();
-	geoX = Symbolic::times(geoX, std::make_shared<Constant>(1.0 / mbdUnits->length));
+	geoX = Symbolic::times(geoX, sptrConstant(1.0 / mbdUnits->length));
 	geoX->createMbD(mbdSys, mbdUnits);
 	auto xBlk = geoX->simplified(geoX);
 
 	userFunc = std::make_shared<BasicUserFunction>(rIJI->at(1), 1.0);
 	parser->parseUserFunction(userFunc);
 	auto geoY = parser->stack->top();
-	geoY = Symbolic::times(geoY, std::make_shared<Constant>(1.0 / mbdUnits->length));
+	geoY = Symbolic::times(geoY, sptrConstant(1.0 / mbdUnits->length));
 	geoY->createMbD(mbdSys, mbdUnits);
 	auto yBlk = geoY->simplified(geoY);
 
 	userFunc = std::make_shared<BasicUserFunction>(rIJI->at(2), 1.0);
 	parser->parseUserFunction(userFunc);
 	auto geoZ = parser->stack->top();
-	geoZ = Symbolic::times(geoZ, std::make_shared<Constant>(1.0 / mbdUnits->length));
+	geoZ = Symbolic::times(geoZ, sptrConstant(1.0 / mbdUnits->length));
 	geoZ->createMbD(mbdSys, mbdUnits);
 	auto zBlk = geoZ->simplified(geoZ);
 
@@ -117,21 +117,21 @@ void MbD::ASMTGeneralMotion::createMbD(std::shared_ptr<System> mbdSys, std::shar
 	userFunc = std::make_shared<BasicUserFunction>(angIJJ->at(0), 1.0);
 	parser->parseUserFunction(userFunc);
 	auto geoPhi = parser->stack->top();
-	geoPhi = Symbolic::times(geoPhi, std::make_shared<Constant>(1.0 / mbdUnits->angle));
+	geoPhi = Symbolic::times(geoPhi, sptrConstant(1.0 / mbdUnits->angle));
 	geoPhi->createMbD(mbdSys, mbdUnits);
 	auto phiBlk = geoPhi->simplified(geoPhi);
 
 	userFunc = std::make_shared<BasicUserFunction>(angIJJ->at(1), 1.0);
 	parser->parseUserFunction(userFunc);
 	auto geoThe = parser->stack->top();
-	geoThe = Symbolic::times(geoThe, std::make_shared<Constant>(1.0 / mbdUnits->angle));
+	geoThe = Symbolic::times(geoThe, sptrConstant(1.0 / mbdUnits->angle));
 	geoThe->createMbD(mbdSys, mbdUnits);
 	auto theBlk = geoThe->simplified(geoThe);
 
 	userFunc = std::make_shared<BasicUserFunction>(angIJJ->at(2), 1.0);
 	parser->parseUserFunction(userFunc);
 	auto geoPsi = parser->stack->top();
-	geoPsi = Symbolic::times(geoPsi, std::make_shared<Constant>(1.0 / mbdUnits->angle));
+	geoPsi = Symbolic::times(geoPsi, sptrConstant(1.0 / mbdUnits->angle));
 	geoPsi->createMbD(mbdSys, mbdUnits);
 	auto psiBlk = geoPsi->simplified(geoPsi);
 

OndselSolver/ASMTItem.cpp

@@ -11,6 +11,7 @@
 #include "ASMTSpatialContainer.h"
 #include "ASMTAssembly.h"
 #include "Constant.h"
+#include <algorithm>
 
 using namespace MbD;
 

OndselSolver/ASMTMarker.cpp

@@ -7,6 +7,7 @@
  ***************************************************************************/
 
 #include "ASMTMarker.h"
+#include "FullMatrix.h"
 #include "ASMTRefItem.h"
 #include "ASMTPart.h"
 #include "Part.h"
@@ -14,57 +15,57 @@
 #include "MarkerFrame.h"
 #include "ASMTPrincipalMassMarker.h"
 
-using namespace MbD;
+namespace MbD {
+    void ASMTMarker::parseASMT(std::vector<std::string>& lines)
+    {
+        readName(lines);
+        readPosition3D(lines);
+        readRotationMatrix(lines);
+    }
 
-void MbD::ASMTMarker::parseASMT(std::vector<std::string>& lines)
-{
-	readName(lines);
-	readPosition3D(lines);
-	readRotationMatrix(lines);
-}
-
-FColDsptr MbD::ASMTMarker::rpmp()
-{
-	//p is cm
-	auto refItem = static_cast<ASMTRefItem*>(owner);
-	auto& rPrefP = refItem->position3D;
-	auto& aAPref = refItem->rotationMatrix;
-	auto& rrefmref = position3D;
-	auto rPmP = rPrefP->plusFullColumn(aAPref->timesFullColumn(rrefmref));
-	auto& principalMassMarker = static_cast<ASMTPart*>(refItem->owner)->principalMassMarker;
-	auto& rPcmP = principalMassMarker->position3D;
-	auto& aAPcm = principalMassMarker->rotationMatrix;
-	auto rpmp = aAPcm->transposeTimesFullColumn(rPmP->minusFullColumn(rPcmP));
-	return rpmp;
-}
-
-FMatDsptr MbD::ASMTMarker::aApm()
-{
-	//p is cm
-	auto refItem = static_cast<ASMTRefItem*>(owner);
-	auto& aAPref = refItem->rotationMatrix;
-	auto& aArefm = rotationMatrix;
-	auto& principalMassMarker = static_cast<ASMTPart*>(refItem->owner)->principalMassMarker;
-	auto& aAPcm = principalMassMarker->rotationMatrix;
-	auto aApm = aAPcm->transposeTimesFullMatrix(aAPref->timesFullMatrix(aArefm));
-	return aApm;
-}
-
-void MbD::ASMTMarker::createMbD(std::shared_ptr<System> mbdSys, std::shared_ptr<Units> mbdUnits)
-{
-	auto mkr = CREATE<MarkerFrame>::With(name.c_str());
-	auto prt = std::static_pointer_cast<Part>(partOrAssembly()->mbdObject);
-	prt->partFrame->addMarkerFrame(mkr);
-
-	mkr->rpmp = rpmp()->times(1.0 / mbdUnits->length);
-	mkr->aApm = aApm();
-	mbdObject = mkr->endFrames->at(0);
-}
-
-void MbD::ASMTMarker::storeOnLevel(std::ofstream& os, int level)
-{
-	storeOnLevelString(os, level, "Marker");
-	storeOnLevelString(os, level + 1, "Name");
-	storeOnLevelString(os, level + 2, name);
-	ASMTSpatialItem::storeOnLevel(os, level);
+    FColDsptr ASMTMarker::rpmp()
+    {
+        //p is cm
+        auto refItem = static_cast<ASMTRefItem*>(owner);
+        auto& rPrefP = refItem->position3D;
+        auto& aAPref = refItem->rotationMatrix;
+        auto& rrefmref = position3D;
+        auto rPmP = rPrefP->plusFullColumn(aAPref->timesFullColumn(rrefmref));
+        auto& principalMassMarker = static_cast<ASMTPart*>(refItem->owner)->principalMassMarker;
+        auto& rPcmP = principalMassMarker->position3D;
+        auto& aAPcm = principalMassMarker->rotationMatrix;
+        auto rpmp = aAPcm->transposeTimesFullColumn(rPmP->minusFullColumn(rPcmP));
+        return rpmp;
+    }
+
+    FMatDsptr ASMTMarker::aApm()
+    {
+        //p is cm
+        auto refItem = static_cast<ASMTRefItem*>(owner);
+        auto& aAPref = refItem->rotationMatrix;
+        auto& aArefm = rotationMatrix;
+        auto& principalMassMarker = static_cast<ASMTPart*>(refItem->owner)->principalMassMarker;
+        auto& aAPcm = principalMassMarker->rotationMatrix;
+        auto aApm = aAPcm->transposeTimesFullMatrix(aAPref->timesFullMatrix(aArefm));
+        return aApm;
+    }
+
+    void ASMTMarker::createMbD(std::shared_ptr<System> mbdSys, std::shared_ptr<Units> mbdUnits)
+    {
+        auto mkr = CREATE<MarkerFrame>::With(name.c_str());
+        auto prt = std::static_pointer_cast<Part>(partOrAssembly()->mbdObject);
+        prt->partFrame->addMarkerFrame(mkr);
+
+        mkr->rpmp = rpmp()->times(1.0 / mbdUnits->length);
+        mkr->aApm = aApm();
+        mbdObject = mkr->endFrames->at(0);
+    }
+
+    void ASMTMarker::storeOnLevel(std::ofstream& os, int level)
+    {
+        storeOnLevelString(os, level, "Marker");
+        storeOnLevelString(os, level + 1, "Name");
+        storeOnLevelString(os, level + 2, name);
+        ASMTSpatialItem::storeOnLevel(os, level);
+    }
 }

OndselSolver/ASMTMarker.h

@@ -23,7 +23,6 @@ namespace MbD {
         FMatDsptr aApm();
         void createMbD(std::shared_ptr<System> mbdSys, std::shared_ptr<Units> mbdUnits) override;
         void storeOnLevel(std::ofstream& os, int level) override;
-
     };
 }
 

OndselSolver/ASMTMotion.cpp

@@ -9,36 +9,36 @@
 
 #include "ASMTMotion.h"
 
-using namespace MbD;
+namespace MbD {
+    void ASMTMotion::readMotionSeries(std::vector<std::string>& lines)
+    {
+        std::string str = lines[0];
+        std::string substr = "MotionSeries";
+        auto pos = str.find(substr);
+        assert(pos != std::string::npos);
+        str.erase(0, pos + substr.length());
+        auto seriesName = readString(str);
+        assert(fullName("") == seriesName);
+        lines.erase(lines.begin());
+        readFXonIs(lines);
+        readFYonIs(lines);
+        readFZonIs(lines);
+        readTXonIs(lines);
+        readTYonIs(lines);
+        readTZonIs(lines);
+    }
 
-void MbD::ASMTMotion::readMotionSeries(std::vector<std::string>& lines)
-{
-	std::string str = lines[0];
-	std::string substr = "MotionSeries";
-	auto pos = str.find(substr);
-	assert(pos != std::string::npos);
-	str.erase(0, pos + substr.length());
-	auto seriesName = readString(str);
-	assert(fullName("") == seriesName);
-	lines.erase(lines.begin());
-	readFXonIs(lines);
-	readFYonIs(lines);
-	readFZonIs(lines);
-	readTXonIs(lines);
-	readTYonIs(lines);
-	readTZonIs(lines);
-}
+    void ASMTMotion::initMarkers()
+    {
+    }
 
-void MbD::ASMTMotion::initMarkers()
-{
-}
+    void ASMTMotion::storeOnLevel(std::ofstream& os, int level)
+    {
+        assert(false);
+    }
 
-void MbD::ASMTMotion::storeOnLevel(std::ofstream& os, int level)
-{
-	assert(false);
-}
-
-void MbD::ASMTMotion::storeOnTimeSeries(std::ofstream& os)
-{
-	assert(false);
+    void ASMTMotion::storeOnTimeSeries(std::ofstream& os)
+    {
+        assert(false);
+    }
 }

OndselSolver/ASMTPart.cpp

@@ -11,7 +11,7 @@
 #include "CREATE.h"
 #include "ASMTPrincipalMassMarker.h"
 #include "Part.h"
-
+#include <algorithm>
 
 using namespace MbD;
 

OndselSolver/ASMTPointInPlaneJoint.cpp

@@ -25,11 +25,14 @@ void MbD::ASMTPointInPlaneJoint::parseASMT(std::vector<std::string>& lines)
 
 void MbD::ASMTPointInPlaneJoint::readOffset(std::vector<std::string>& lines)
 {
-	assert(lines[0].find("offset") != std::string::npos);
-	lines.erase(lines.begin());
-	offset = readDouble(lines[0]);
-	lines.erase(lines.begin());
-
+	if (lines[0].find("offset") == std::string::npos) {
+		offset = 0.0;
+	}
+	else {
+		lines.erase(lines.begin());
+		offset = readDouble(lines[0]);
+		lines.erase(lines.begin());
+	}
 }
 
 void MbD::ASMTPointInPlaneJoint::createMbD(std::shared_ptr<System> mbdSys, std::shared_ptr<Units> mbdUnits)

OndselSolver/ASMTPrincipalMassMarker.cpp

@@ -12,6 +12,11 @@
 
 using namespace MbD;
 
+MbD::ASMTPrincipalMassMarker::ASMTPrincipalMassMarker()
+{
+	name = "MassMarker";
+}
+
 void MbD::ASMTPrincipalMassMarker::parseASMT(std::vector<std::string>& lines)
 {
 	size_t pos = lines[0].find_first_not_of("\t");
@@ -26,7 +31,7 @@ void MbD::ASMTPrincipalMassMarker::parseASMT(std::vector<std::string>& lines)
 	lines.erase(lines.begin());
 	assert(lines[0] == (leadingTabs + "RotationMatrix"));
 	lines.erase(lines.begin());
-	rotationMatrix = std::make_shared<FullMatrix<double>>(3);
+	rotationMatrix = std::make_shared<FullMatrixDouble>(3);
 	for (int i = 0; i < 3; i++)
 	{
 		auto row = readRowOfDoubles(lines[0]);
@@ -39,7 +44,7 @@ void MbD::ASMTPrincipalMassMarker::parseASMT(std::vector<std::string>& lines)
 	lines.erase(lines.begin());
 	assert(lines[0] == (leadingTabs + "MomentOfInertias"));
 	lines.erase(lines.begin());
-	momentOfInertias = std::make_shared<DiagonalMatrix<double>>(3);
+	momentOfInertias = std::make_shared<DiagonalMatrix>(3);
 	auto row = readRowOfDoubles(lines[0]);
 	lines.erase(lines.begin());
 	for (int i = 0; i < 3; i++)
@@ -70,7 +75,7 @@ void MbD::ASMTPrincipalMassMarker::setMomentOfInertias(DiagMatDsptr mat)
 // Overloads to simplify syntax.
 void MbD::ASMTPrincipalMassMarker::setMomentOfInertias(double a, double b, double c)
 {
-	momentOfInertias = std::make_shared<DiagonalMatrix<double>>(ListD{ a, b, c });
+	momentOfInertias = std::make_shared<DiagonalMatrix>(ListD{ a, b, c });
 }
 
 void MbD::ASMTPrincipalMassMarker::storeOnLevel(std::ofstream& os, int level)

OndselSolver/ASMTPrincipalMassMarker.h

@@ -5,29 +5,30 @@
  *                                                                         *
  *   See LICENSE file for details about copyright.                         *
  ***************************************************************************/
- 
+
 #pragma once
 
 #include "ASMTSpatialItem.h"
 
 namespace MbD {
-    class EXPORT ASMTPrincipalMassMarker : public ASMTSpatialItem
-    {
-        //
-    public:
-        void parseASMT(std::vector<std::string>& lines) override;
-        void setMass(double mass);
-        void setDensity(double density);
-        void setMomentOfInertias(DiagMatDsptr momentOfInertias);
+	class EXPORT ASMTPrincipalMassMarker : public ASMTSpatialItem
+	{
+		//
+	public:
+		ASMTPrincipalMassMarker();
+		void parseASMT(std::vector<std::string>& lines) override;
+		void setMass(double mass);
+		void setDensity(double density);
+		void setMomentOfInertias(DiagMatDsptr momentOfInertias);
 
-        // Overloads to simplify syntax.
-        void setMomentOfInertias(double a, double b, double c);
-        void storeOnLevel(std::ofstream& os, int level) override;
+		// Overloads to simplify syntax.
+		void setMomentOfInertias(double a, double b, double c);
+		void storeOnLevel(std::ofstream& os, int level) override;
 
-        double mass = 0.0;
-        double density = 0.0;
-        DiagMatDsptr momentOfInertias = std::make_shared<DiagonalMatrix<double>>(ListD{ 0.,0.,0. });
+		double mass = 1.0;
+		double density = 10.0;
+		DiagMatDsptr momentOfInertias = std::make_shared<DiagonalMatrix>(ListD{ 1.0, 2.0, 3.0 });
 
-    };
+	};
 }
 

OndselSolver/ASMTRefItem.cpp

@@ -8,6 +8,7 @@
  
 #include "ASMTRefItem.h"
 #include "CREATE.h"
+#include <algorithm>
 
 using namespace MbD;
 

OndselSolver/ASMTRotationalMotion.cpp

@@ -68,9 +68,13 @@ void MbD::ASMTRotationalMotion::createMbD(std::shared_ptr<System> mbdSys, std::s
 	auto userFunc = std::make_shared<BasicUserFunction>(rotationZ, 1.0);
 	parser->parseUserFunction(userFunc);
 	auto geoPhi = parser->stack->top();
-	geoPhi = Symbolic::times(geoPhi, std::make_shared<Constant>(1.0 / mbdUnits->angle));
+	std::cout << *geoPhi << std::endl;
+	geoPhi = Symbolic::times(geoPhi, sptrConstant(1.0 / mbdUnits->angle));
 	geoPhi->createMbD(mbdSys, mbdUnits);
-	std::static_pointer_cast<ZRotation>(mbdObject)->phiBlk = geoPhi->simplified(geoPhi);
+	std::cout << *geoPhi << std::endl;
+	auto simple = geoPhi->simplified(geoPhi);
+	std::cout << *simple << std::endl;
+	std::static_pointer_cast<ZRotation>(mbdObject)->phiBlk = simple;
 }
 
 std::shared_ptr<Joint> MbD::ASMTRotationalMotion::mbdClassNew()

OndselSolver/ASMTSpatialContainer.cpp

@@ -550,13 +550,25 @@ void MbD::ASMTSpatialContainer::setOmega3D(double a, double b, double c)
 void MbD::ASMTSpatialContainer::storeOnLevelVelocity(std::ofstream& os, int level)
 {
 	storeOnLevelString(os, level, "Velocity3D");
-	storeOnLevelArray(os, level + 1, *velocity3D);
+	if (vxs == nullptr || vxs->empty()) {
+		storeOnLevelArray(os, level + 1, *velocity3D);
+	}
+	else {
+		auto array = getVelocity3D(0);
+		storeOnLevelArray(os, level + 1, *array);
+	}
 }
 
 void MbD::ASMTSpatialContainer::storeOnLevelOmega(std::ofstream& os, int level)
 {
 	storeOnLevelString(os, level, "Omega3D");
-	storeOnLevelArray(os, level + 1, *omega3D);
+	if (omexs == nullptr || omexs->empty()) {
+		storeOnLevelArray(os, level + 1, *omega3D);
+	}
+	else {
+		auto array = getOmega3D(0);
+		storeOnLevelArray(os, level + 1, *array);
+	}
 }
 
 void MbD::ASMTSpatialContainer::storeOnLevelRefPoints(std::ofstream& os, int level)
@@ -697,3 +709,21 @@ void MbD::ASMTSpatialContainer::storeOnTimeSeries(std::ofstream& os)
 	}
 	os << std::endl;
 }
+
+FColDsptr MbD::ASMTSpatialContainer::getVelocity3D(size_t i)
+{
+	auto vec3 = std::make_shared<FullColumn<double>>(3);
+	vec3->atiput(0, vxs->at(i));
+	vec3->atiput(1, vys->at(i));
+	vec3->atiput(2, vzs->at(i));
+	return vec3;
+}
+
+FColDsptr MbD::ASMTSpatialContainer::getOmega3D(size_t i)
+{
+	auto vec3 = std::make_shared<FullColumn<double>>(3);
+	vec3->atiput(0, omexs->at(i));
+	vec3->atiput(1, omeys->at(i));
+	vec3->atiput(2, omezs->at(i));
+	return vec3;
+}

OndselSolver/ASMTSpatialContainer.h

@@ -86,16 +86,16 @@ namespace MbD {
         void storeOnLevelRefCurves(std::ofstream& os, int level);
         void storeOnLevelRefSurfaces(std::ofstream& os, int level);
         void storeOnTimeSeries(std::ofstream& os) override;
+        FColDsptr getVelocity3D(size_t i);
+        FColDsptr getOmega3D(size_t i);
 
         FColDsptr velocity3D = std::make_shared<FullColumn<double>>(3);
         FColDsptr omega3D = std::make_shared<FullColumn<double>>(3);
         std::shared_ptr<std::vector<std::shared_ptr<ASMTRefPoint>>> refPoints;
         std::shared_ptr<std::vector<std::shared_ptr<ASMTRefCurve>>> refCurves;
         std::shared_ptr<std::vector<std::shared_ptr<ASMTRefSurface>>> refSurfaces;
-        FRowDsptr xs, ys, zs, bryxs, bryys, bryzs;
         FRowDsptr vxs, vys, vzs, omexs, omeys, omezs;
         FRowDsptr axs, ays, azs, alpxs, alpys, alpzs;
-        FRowDsptr inxs, inys, inzs, inbryxs, inbryys, inbryzs;
         FRowDsptr invxs, invys, invzs, inomexs, inomeys, inomezs;
         FRowDsptr inaxs, inays, inazs, inalpxs, inalpys, inalpzs;
         std::shared_ptr<ASMTPrincipalMassMarker> principalMassMarker = std::make_shared<ASMTPrincipalMassMarker>();

OndselSolver/ASMTSpatialItem.cpp

@@ -10,6 +10,7 @@
 #include "Units.h"
 #include "Part.h"
 #include "ASMTSpatialContainer.h"
+#include "EulerAngles.h"
 
 using namespace MbD;
 
@@ -47,7 +48,7 @@ void MbD::ASMTSpatialItem::readRotationMatrix(std::vector<std::string>& lines)
 {
 	assert(lines[0].find("RotationMatrix") != std::string::npos);
 	lines.erase(lines.begin());
-	rotationMatrix = std::make_shared<FullMatrix<double>>(3, 0);
+	rotationMatrix = std::make_shared<FullMatrixDouble>(3, 0);
 	for (int i = 0; i < 3; i++)
 	{
 		auto& row = rotationMatrix->at(i);
@@ -86,7 +87,7 @@ void MbD::ASMTSpatialItem::setRotationMatrix(double v11, double v12, double v13,
 	double v21, double v22, double v23,
 	double v31, double v32, double v33)
 {
-	rotationMatrix = std::make_shared<FullMatrix<double>>(ListListD{
+	rotationMatrix = std::make_shared<FullMatrixDouble>(ListListD{
 		{ v11, v12, v13 },
 		{ v21, v22, v23 },
 		{ v31, v32, v33 }
@@ -102,14 +103,50 @@ void MbD::ASMTSpatialItem::storeOnLevel(std::ofstream& os, int level)
 void MbD::ASMTSpatialItem::storeOnLevelPosition(std::ofstream& os, int level)
 {
 	storeOnLevelString(os, level, "Position3D");
-	storeOnLevelArray(os, level + 1, *position3D);
+	if (xs == nullptr || xs->empty()) {
+		storeOnLevelArray(os, level + 1, *position3D);
+	}
+	else {
+		auto array = getPosition3D(0);
+		storeOnLevelArray(os, level + 1, *array);
+	}
 }
 
 void MbD::ASMTSpatialItem::storeOnLevelRotationMatrix(std::ofstream& os, int level)
 {
 	storeOnLevelString(os, level, "RotationMatrix");
-	for (int i = 0; i < 3; i++)
-	{
-		storeOnLevelArray(os, level + 1, *rotationMatrix->at(i));
+	if (xs == nullptr || xs->empty()) {
+		for (int i = 0; i < 3; i++)
+		{
+			storeOnLevelArray(os, level + 1, *rotationMatrix->at(i));
+		}
 	}
+	else {
+		auto rotMat = getRotationMatrix(0);
+		for (int i = 0; i < 3; i++)
+		{
+			storeOnLevelArray(os, level + 1, *rotMat->at(i));
+		}
+	}
+
+}
+
+FColDsptr MbD::ASMTSpatialItem::getPosition3D(size_t i)
+{
+	auto vec3 = std::make_shared<FullColumn<double>>(3);
+	vec3->atiput(0, xs->at(i));
+	vec3->atiput(1, ys->at(i));
+	vec3->atiput(2, zs->at(i));
+	return vec3;
+}
+
+FMatDsptr MbD::ASMTSpatialItem::getRotationMatrix(size_t i)
+{
+	auto bryantAngles = std::make_shared<EulerAngles<double>>();
+	bryantAngles->setRotOrder(1, 2, 3);
+	bryantAngles->at(0) = bryxs->at(i);
+	bryantAngles->at(1) = bryys->at(i);
+	bryantAngles->at(2) = bryzs->at(i);
+	bryantAngles->calc();
+	return bryantAngles->aA;
 }

OndselSolver/ASMTSpatialItem.h

@@ -31,13 +31,18 @@ namespace MbD {
 		void storeOnLevel(std::ofstream& os, int level) override;
 		void storeOnLevelPosition(std::ofstream& os, int level);
 		void storeOnLevelRotationMatrix(std::ofstream& os, int level);
+		FColDsptr getPosition3D(size_t i);
+		FMatDsptr getRotationMatrix(size_t i);
 
 		FColDsptr position3D = std::make_shared<FullColumn<double>>(3);
-		FMatDsptr rotationMatrix = std::make_shared<FullMatrix<double>>(ListListD{
+		FMatDsptr rotationMatrix = std::make_shared<FullMatrixDouble>(ListListD{
 				{ 1, 0, 0 },
 				{ 0, 1, 0 },
 				{ 0, 0, 1 }
 			});
+		FRowDsptr xs, ys, zs, bryxs, bryys, bryzs;
+		FRowDsptr inxs, inys, inzs, inbryxs, inbryys, inbryzs;
+
 	};
 }
 

OndselSolver/ASMTTime.cpp

@@ -25,7 +25,27 @@ void MbD::ASMTTime::createMbD(std::shared_ptr<System> mbdSys, std::shared_ptr<Un
 {
 	auto mbdTime = mbdSys->time;
 	if (xx == mbdTime) return;
-	auto timeScale = std::make_shared<Constant>(mbdUnits->time);
+	auto timeScale = sptrConstant(mbdUnits->time);
 	auto geoTime = std::make_shared<Product>(timeScale, mbdTime);
 	this->xexpression(mbdTime, geoTime->simplified(geoTime));
 }
+
+Symsptr MbD::ASMTTime::expandUntil(Symsptr sptr, std::shared_ptr<std::unordered_set<Symsptr>> set)
+{
+	return sptr;
+}
+
+Symsptr MbD::ASMTTime::simplifyUntil(Symsptr sptr, std::shared_ptr<std::unordered_set<Symsptr>> set)
+{
+	return sptr;
+}
+
+bool MbD::ASMTTime::isVariable()
+{
+	return true;
+}
+
+void MbD::ASMTTime::setValue(double val)
+{
+	xx->setValue(val);
+}

OndselSolver/ASMTTime.h

@@ -19,6 +19,10 @@ namespace MbD {
     public:
         void deleteMbD();
         void createMbD(std::shared_ptr<System> mbdSys, std::shared_ptr<Units> mbdUnits);
+        Symsptr expandUntil(Symsptr sptr, std::shared_ptr<std::unordered_set<Symsptr>> set) override;
+        Symsptr simplifyUntil(Symsptr sptr, std::shared_ptr<std::unordered_set<Symsptr>> set) override;
+        bool isVariable() override;
+        void setValue(double val) override;
 
     };
 }

OndselSolver/ASMTTranslationalMotion.cpp

@@ -42,7 +42,7 @@ void MbD::ASMTTranslationalMotion::createMbD(std::shared_ptr<System> mbdSys, std
 	auto userFunc = std::make_shared<BasicUserFunction>(translationZ, 1.0);
 	parser->parseUserFunction(userFunc);
 	auto zIJ = parser->stack->top();
-	zIJ = Symbolic::times(zIJ, std::make_shared<Constant>(1.0 / mbdUnits->length));
+	zIJ = Symbolic::times(zIJ, sptrConstant(1.0 / mbdUnits->length));
 	zIJ->createMbD(mbdSys, mbdUnits);
 	std::static_pointer_cast<ZTranslation>(mbdObject)->zBlk = zIJ->simplified(zIJ);
 }

OndselSolver/Abs.cpp

@@ -19,8 +19,13 @@ double MbD::Abs::getValue()
     return std::abs(xx->getValue());
 }
 
+Symsptr MbD::Abs::copyWith(Symsptr arg)
+{
+	return std::make_shared<Abs>(arg);
+}
+
 std::ostream& MbD::Abs::printOn(std::ostream& s) const
 {
-	s << "abs(" << xx << ")";
+	s << "abs(" << *xx << ")";
 	return s;
 }

OndselSolver/Abs.h

@@ -18,6 +18,7 @@ namespace MbD {
         Abs() = default;
         Abs(Symsptr arg);
         double getValue() override;
+        Symsptr copyWith(Symsptr arg) override;
 
         std::ostream& printOn(std::ostream& s) const override;
 

OndselSolver/AngleZIeqcJec.cpp

@@ -65,7 +65,7 @@ void MbD::AngleZIeqcJec::initialize()
 {
 	AngleZIecJec::initialize();
 	pthezpEI = std::make_shared<FullRow<double>>(4);
-	ppthezpEIpEI = std::make_shared<FullMatrix<double>>(4, 4);
+	ppthezpEIpEI = std::make_shared<FullMatrixDouble>(4, 4);
 }
 
 FMatDsptr MbD::AngleZIeqcJec::ppvaluepEIpEI()

OndselSolver/AngleZIeqcJeqc.cpp

@@ -85,8 +85,8 @@ void MbD::AngleZIeqcJeqc::initialize()
 {
 	AngleZIeqcJec::initialize();
 	pthezpEJ = std::make_shared<FullRow<double>>(4);
-	ppthezpEIpEJ = std::make_shared<FullMatrix<double>>(4, 4);
-	ppthezpEJpEJ = std::make_shared<FullMatrix<double>>(4, 4);
+	ppthezpEIpEJ = std::make_shared<FullMatrixDouble>(4, 4);
+	ppthezpEJpEJ = std::make_shared<FullMatrixDouble>(4, 4);
 }
 
 FMatDsptr MbD::AngleZIeqcJeqc::ppvaluepEIpEJ()

OndselSolver/AnyPosICNewtonRaphson.cpp

@@ -32,7 +32,7 @@ void AnyPosICNewtonRaphson::initializeGlobally()
 void AnyPosICNewtonRaphson::createVectorsAndMatrices()
 {
 	qsuOld = std::make_shared<FullColumn<double>>(nqsu);
-	qsuWeights = std::make_shared<DiagonalMatrix<double>>(nqsu);
+	qsuWeights = std::make_shared<DiagonalMatrix>(nqsu);
 	SystemNewtonRaphson::createVectorsAndMatrices();
 }
 

OndselSolver/ArcSine.cpp

@@ -19,8 +19,13 @@ double MbD::ArcSine::getValue()
 	return std::asin(xx->getValue());
 }
 
+Symsptr MbD::ArcSine::copyWith(Symsptr arg)
+{
+	return std::make_shared<ArcSine>(arg);
+}
+
 std::ostream& MbD::ArcSine::printOn(std::ostream& s) const
 {
-	s << "arcsin(" << xx << ")";
+	s << "arcsin(" << *xx << ")";
 	return s;
 }

OndselSolver/ArcSine.h

@@ -18,6 +18,7 @@ namespace MbD {
         ArcSine() = default;
         ArcSine(Symsptr arg);
         double getValue() override;
+        Symsptr copyWith(Symsptr arg) override;
 
         std::ostream& printOn(std::ostream& s) const override;
 

OndselSolver/ArcTan.cpp

@@ -19,8 +19,13 @@ double MbD::ArcTan::getValue()
     return std::atan(xx->getValue());
 }
 
+Symsptr MbD::ArcTan::copyWith(Symsptr arg)
+{
+	return std::make_shared<ArcTan>(arg);
+}
+
 std::ostream& MbD::ArcTan::printOn(std::ostream& s) const
 {
-	s << "arctan(" << xx << ")";
+	s << "arctan(" << *xx << ")";
 	return s;
 }

OndselSolver/ArcTan.h

@@ -18,6 +18,7 @@ namespace MbD {
         ArcTan() = default;
         ArcTan(Symsptr arg);
         double getValue() override;
+        Symsptr copyWith(Symsptr arg) override;
 
         std::ostream& printOn(std::ostream& s) const override;
 

OndselSolver/ArcTan2.cpp

@@ -16,5 +16,5 @@ MbD::ArcTan2::ArcTan2(Symsptr arg, Symsptr arg1) : FunctionXY(arg, arg1)
 
 double MbD::ArcTan2::getValue()
 {
-    return std::atan2(y->getValue(), x->getValue());;
+    return std::atan2(y->getValue(), x->getValue());
 }

OndselSolver/Array.h

@@ -121,7 +121,7 @@ namespace MbD {
 	{
 		for (int ii = 0; ii < this->size(); ii++)
 		{
-			this->at(ii) = array->at(i + ii);
+			this->at(ii) = array->at((size_t)i + ii);
 		}
 	}
 	//template<>
@@ -182,4 +182,4 @@ namespace MbD {
 	{
 		this->at(i) *= factor;
 	}
-}
+}

OndselSolver/AtPointConstraintIqcJc.cpp

@@ -45,7 +45,7 @@ void AtPointConstraintIqcJc::useEquationNumbers()
 void AtPointConstraintIqcJc::fillPosICError(FColDsptr col)
 {
 	Constraint::fillPosICError(col);
-	col->at(iqXIminusOnePlusAxis) -= lam;
+	col->atiminusNumber(iqXIminusOnePlusAxis, lam);
 	col->atiplusFullVectortimes(iqEI, pGpEI, lam);
 }
 

OndselSolver/AtPointConstraintIqcJqc.cpp

@@ -46,7 +46,7 @@ void AtPointConstraintIqcJqc::useEquationNumbers()
 void AtPointConstraintIqcJqc::fillPosICError(FColDsptr col)
 {
 	AtPointConstraintIqcJc::fillPosICError(col);
-	col->at(iqXJminusOnePlusAxis) += lam;
+	col->atiplusNumber(iqXJminusOnePlusAxis, lam);
 	col->atiplusFullVectortimes(iqEJ, pGpEJ, lam);
 }
 

OndselSolver/BasicIntegrator.cpp

@@ -5,7 +5,7 @@
  *                                                                         *
  *   See LICENSE file for details about copyright.                         *
  ***************************************************************************/
- 
+
 #include "BasicIntegrator.h"
 #include "CREATE.h"
 #include "StableBackwardDifference.h"

OndselSolver/BasicQuasiIntegrator.cpp

@@ -14,15 +14,15 @@ using namespace MbD;
 void BasicQuasiIntegrator::firstStep()
 {
 	istep = 0;
-	this->preFirstStep();
+	preFirstStep();
 	iTry = 1;
 	orderNew = 1;
-	this->selectFirstStepSize();
-	this->incrementTime();
-	this->runInitialConditionTypeSolution();
-	//this->reportTrialStepStats();
-	this->postFirstStep();
-	//this->reportStepStats();
+	selectFirstStepSize();
+	incrementTime();
+	runInitialConditionTypeSolution();
+	//reportTrialStepStats();
+	postFirstStep();
+	//reportStepStats();
 }
 
 bool BasicQuasiIntegrator::isRedoingFirstStep()
@@ -32,15 +32,15 @@ bool BasicQuasiIntegrator::isRedoingFirstStep()
 
 void BasicQuasiIntegrator::nextStep()
 {
-	this->preStep();
+	preStep();
 	iTry = 1;
-	this->selectOrder();
-	this->selectStepSize();
-	this->incrementTime();
-	this->runInitialConditionTypeSolution();
-	//this->reportTrialStepStats();
-	this->postStep();
-	//this->reportStepStats();
+	selectOrder();
+	selectStepSize();
+	incrementTime();
+	runInitialConditionTypeSolution();
+	//reportTrialStepStats();
+	postStep();
+	//reportStepStats();
 }
 
 void BasicQuasiIntegrator::runInitialConditionTypeSolution()

OndselSolver/CADSystem.cpp

@@ -87,9 +87,9 @@ void CADSystem::runOndselSinglePendulum()
 	auto assembly1 = CREATE<Part>::With("/Assembly1");
 	std::cout << "assembly1->name " << assembly1->name << std::endl;
 	assembly1->m = 0.0;
-	assembly1->aJ = std::make_shared<DiagonalMatrix<double>>(ListD{ 0, 0, 0 });
+	assembly1->aJ = std::make_shared<DiagonalMatrix>(ListD{ 0, 0, 0 });
 	qX = std::make_shared<FullColumn<double>>(ListD{ 0, 0, 0 });
-	aAap = std::make_shared<FullMatrix<double>>(ListListD{
+	aAap = std::make_shared<FullMatrixDouble>(ListListD{
 		{ 1, 0, 0 },
 		{ 0, 1, 0 },
 		{ 0, 0, 1 }
@@ -106,7 +106,7 @@ void CADSystem::runOndselSinglePendulum()
 		auto marker2 = CREATE<MarkerFrame>::With("/Assembly1/Marker2");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ 0.0, 0.0, 0.0 });
 		marker2->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{ 1, 0, 0 },
 			{ 0, 1, 0 },
 			{ 0, 0, 1 }
@@ -117,7 +117,7 @@ void CADSystem::runOndselSinglePendulum()
 		auto marker1 = CREATE<MarkerFrame>::With("/Assembly1/Marker1");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ 0.0, 3.0, 0.0 });
 		marker1->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{ 1, 0, 0 },
 			{ 0, 0, 1 },
 			{ 0, -1, 0 }
@@ -130,9 +130,9 @@ void CADSystem::runOndselSinglePendulum()
 	auto crankPart1 = CREATE<Part>::With("/Assembly1/Part1");
 	std::cout << "crankPart1->name " << crankPart1->name << std::endl;
 	crankPart1->m = 1.0;
-	crankPart1->aJ = std::make_shared<DiagonalMatrix<double>>(ListD{ 1, 1, 1 });
+	crankPart1->aJ = std::make_shared<DiagonalMatrix>(ListD{ 1, 1, 1 });
 	qX = std::make_shared<FullColumn<double>>(ListD{ 0.4, 0.0, -0.05 });
-	aAap = std::make_shared<FullMatrix<double>>(ListListD{
+	aAap = std::make_shared<FullMatrixDouble>(ListListD{
 		{ 1, 0, 0 },
 		{ 0, 1, 0 },
 		{ 0, 0, 1 }
@@ -149,7 +149,7 @@ void CADSystem::runOndselSinglePendulum()
 		auto marker1 = CREATE<MarkerFrame>::With("/Assembly1/Part1/Marker1");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ -0.4, 0.0, 0.05 });
 		marker1->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{ 1, 0, 0 },
 			{ 0, 1, 0 },
 			{ 0, 0, 1 }
@@ -160,7 +160,7 @@ void CADSystem::runOndselSinglePendulum()
 		auto marker2 = CREATE<MarkerFrame>::With("/Assembly1/Part1/Marker2");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ 0.4, 0.0, 0.05 });
 		marker2->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{ 1, 0, 0 },
 			{ 0, 1, 0 },
 			{ 0, 0, 1 }
@@ -222,9 +222,9 @@ void CADSystem::runOndselDoublePendulum()
 	auto assembly1 = CREATE<Part>::With("/Assembly1");
 	std::cout << "assembly1->name " << assembly1->name << std::endl;
 	assembly1->m = 0.0;
-	assembly1->aJ = std::make_shared<DiagonalMatrix<double>>(ListD{ 0, 0, 0 });
+	assembly1->aJ = std::make_shared<DiagonalMatrix>(ListD{ 0, 0, 0 });
 	qX = std::make_shared<FullColumn<double>>(ListD{ 0, 0, 0 });
-	aAap = std::make_shared<FullMatrix<double>>(ListListD{
+	aAap = std::make_shared<FullMatrixDouble>(ListListD{
 		{ 1, 0, 0 },
 		{ 0, 1, 0 },
 		{ 0, 0, 1 }
@@ -241,7 +241,7 @@ void CADSystem::runOndselDoublePendulum()
 		auto marker2 = CREATE<MarkerFrame>::With("/Assembly1/Marker2");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ 0.0, 0.0, 0.0 });
 		marker2->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{ 1, 0, 0 },
 			{ 0, 1, 0 },
 			{ 0, 0, 1 }
@@ -252,7 +252,7 @@ void CADSystem::runOndselDoublePendulum()
 		auto marker1 = CREATE<MarkerFrame>::With("/Assembly1/Marker1");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ 0.0, 3.0, 0.0 });
 		marker1->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{ 1, 0, 0 },
 			{ 0, 0, 1 },
 			{ 0, -1, 0 }
@@ -265,9 +265,9 @@ void CADSystem::runOndselDoublePendulum()
 	auto crankPart1 = CREATE<Part>::With("/Assembly1/Part1");
 	std::cout << "crankPart1->name " << crankPart1->name << std::endl;
 	crankPart1->m = 1.0;
-	crankPart1->aJ = std::make_shared<DiagonalMatrix<double>>(ListD{ 1, 1, 1 });
+	crankPart1->aJ = std::make_shared<DiagonalMatrix>(ListD{ 1, 1, 1 });
 	qX = std::make_shared<FullColumn<double>>(ListD{ 0.4, 0.0, -0.05 });
-	aAap = std::make_shared<FullMatrix<double>>(ListListD{
+	aAap = std::make_shared<FullMatrixDouble>(ListListD{
 		{ 1, 0, 0 },
 		{ 0, 1, 0 },
 		{ 0, 0, 1 }
@@ -284,7 +284,7 @@ void CADSystem::runOndselDoublePendulum()
 		auto marker1 = CREATE<MarkerFrame>::With("/Assembly1/Part1/Marker1");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ -0.4, 0.0, 0.05 });
 		marker1->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{ 1, 0, 0 },
 			{ 0, 1, 0 },
 			{ 0, 0, 1 }
@@ -295,7 +295,7 @@ void CADSystem::runOndselDoublePendulum()
 		auto marker2 = CREATE<MarkerFrame>::With("/Assembly1/Part1/Marker2");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ 0.4, 0.0, 0.05 });
 		marker2->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{ 1, 0, 0 },
 			{ 0, 1, 0 },
 			{ 0, 0, 1 }
@@ -307,7 +307,7 @@ void CADSystem::runOndselDoublePendulum()
 	auto conrodPart2 = CREATE<Part>::With("/Assembly1/Part2");
 	std::cout << "conrodPart2->name " << conrodPart2->name << std::endl;
 	conrodPart2->m = 1.0;
-	conrodPart2->aJ = std::make_shared<DiagonalMatrix<double>>(ListD{ 1, 1, 1 });
+	conrodPart2->aJ = std::make_shared<DiagonalMatrix>(ListD{ 1, 1, 1 });
 	qX = std::make_shared<FullColumn<double>>(ListD{ 0.15, 0.1, 0.05 });
 	qE = std::make_shared<FullColumn<double>>(ListD{ 0.0, 0.0, 1.0, 0.0 });
 	auto eulerParameters = CREATE<EulerParameters<double>>::With(ListD{ 0.0, 0.0, 1.0, 0.0 });
@@ -325,7 +325,7 @@ void CADSystem::runOndselDoublePendulum()
 		auto marker1 = CREATE<MarkerFrame>::With("/Assembly1/Part2/Marker1");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ -0.65, 0.0, -0.05 });
 		marker1->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{1.0, 0.0, 0.0},
 			{0.0, 1.0, 0.0},
 			{0.0, 0.0, 1.0}
@@ -336,7 +336,7 @@ void CADSystem::runOndselDoublePendulum()
 		auto marker2 = CREATE<MarkerFrame>::With("/Assembly1/Part2/Marker2");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ 0.65, 0.0, -0.05 });
 		marker2->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{1.0, 0.0, 0.0},
 			{0.0, 1.0, 0.0},
 			{0.0, 0.0, 1.0}
@@ -396,7 +396,7 @@ void CADSystem::runOndselPiston()
 	auto assembly1 = CREATE<Part>::With("/Assembly1");
 	std::cout << "assembly1->name " << assembly1->name << std::endl;
 	assembly1->m = 0.0;
-	assembly1->aJ = std::make_shared<DiagonalMatrix<double>>(ListD{ 0, 0, 0 });
+	assembly1->aJ = std::make_shared<DiagonalMatrix>(ListD{ 0, 0, 0 });
 	qX = std::make_shared<FullColumn<double>>(ListD{ 0, 0, 0 });
 	qE = std::make_shared<FullColumn<double>>(ListD{ 0, 0, 0, 1 });
 	assembly1->setqX(qX);
@@ -417,7 +417,7 @@ void CADSystem::runOndselPiston()
 		auto marker2 = CREATE<MarkerFrame>::With("/Assembly1/Marker2");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ 0.0, 0.0, 0.0 });
 		marker2->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{ 1, 0, 0 },
 			{ 0, 1, 0 },
 			{ 0, 0, 1 }
@@ -428,7 +428,7 @@ void CADSystem::runOndselPiston()
 		auto marker1 = CREATE<MarkerFrame>::With("/Assembly1/Marker1");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ 0.0, 3.0, 0.0 });
 		marker1->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{ 1, 0, 0 },
 			{ 0, 0, 1 },
 			{ 0, -1, 0 }
@@ -441,7 +441,7 @@ void CADSystem::runOndselPiston()
 	auto crankPart1 = CREATE<Part>::With("/Assembly1/Part1");
 	std::cout << "crankPart1->name " << crankPart1->name << std::endl;
 	crankPart1->m = 1.0;
-	crankPart1->aJ = std::make_shared<DiagonalMatrix<double>>(ListD{ 1, 1, 1 });
+	crankPart1->aJ = std::make_shared<DiagonalMatrix>(ListD{ 1, 1, 1 });
 	qX = std::make_shared<FullColumn<double>>(ListD{ 0.4, 0.0, -0.05 });
 	qE = std::make_shared<FullColumn<double>>(ListD{ 0.0, 0.0, 0.0, 1.0 });
 	crankPart1->setqX(qX);
@@ -460,7 +460,7 @@ void CADSystem::runOndselPiston()
 		auto marker1 = CREATE<MarkerFrame>::With("/Assembly1/Part1/Marker1");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ -0.4, 0.0, 0.05 });
 		marker1->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{ 1, 0, 0 },
 			{ 0, 1, 0 },
 			{ 0, 0, 1 }
@@ -471,7 +471,7 @@ void CADSystem::runOndselPiston()
 		auto marker2 = CREATE<MarkerFrame>::With("/Assembly1/Part1/Marker2");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ 0.4, 0.0, 0.05 });
 		marker2->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{ 1, 0, 0 },
 			{ 0, 1, 0 },
 			{ 0, 0, 1 }
@@ -483,7 +483,7 @@ void CADSystem::runOndselPiston()
 	auto conrodPart2 = CREATE<Part>::With("/Assembly1/Part2");
 	std::cout << "conrodPart2->name " << conrodPart2->name << std::endl;
 	conrodPart2->m = 1.0;
-	conrodPart2->aJ = std::make_shared<DiagonalMatrix<double>>(ListD{ 1, 1, 1 });
+	conrodPart2->aJ = std::make_shared<DiagonalMatrix>(ListD{ 1, 1, 1 });
 	qX = std::make_shared<FullColumn<double>>(ListD{ 0.15, 0.1, 0.05 });
 	qE = std::make_shared<FullColumn<double>>(ListD{ 0.0, 0.0, 1.0, 0.0 });
 	conrodPart2->setqX(qX);
@@ -502,7 +502,7 @@ void CADSystem::runOndselPiston()
 		auto marker1 = CREATE<MarkerFrame>::With("/Assembly1/Part2/Marker1");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ -0.65, 0.0, -0.05 });
 		marker1->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{1.0, 0.0, 0.0},
 			{0.0, 1.0, 0.0},
 			{0.0, 0.0, 1.0}
@@ -513,7 +513,7 @@ void CADSystem::runOndselPiston()
 		auto marker2 = CREATE<MarkerFrame>::With("/Assembly1/Part2/Marker2");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ 0.65, 0.0, -0.05 });
 		marker2->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{1.0, 0.0, 0.0},
 			{0.0, 1.0, 0.0},
 			{0.0, 0.0, 1.0}
@@ -525,7 +525,7 @@ void CADSystem::runOndselPiston()
 	auto pistonPart3 = CREATE<Part>::With("/Assembly1/Part3");
 	std::cout << "pistonPart3->name " << pistonPart3->name << std::endl;
 	pistonPart3->m = 1.0;
-	pistonPart3->aJ = std::make_shared<DiagonalMatrix<double>>(ListD{ 1, 1, 1 });
+	pistonPart3->aJ = std::make_shared<DiagonalMatrix>(ListD{ 1, 1, 1 });
 	qX = std::make_shared<FullColumn<double>>(ListD{ -0.0, 1.5, 0.0 });
 	qE = std::make_shared<FullColumn<double>>(ListD{ 0.70710678118655, 0.70710678118655, 0.0, 0.0 });
 	pistonPart3->setqX(qX);
@@ -544,7 +544,7 @@ void CADSystem::runOndselPiston()
 		auto marker1 = CREATE<MarkerFrame>::With("/Assembly1/Part3/Marker1");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ -0.5, 0.0, 0.0 });
 		marker1->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{0.0, 1.0, 0.0},
 			{1.0, 0.0, 0.0},
 			{0.0, 0.0, -1.0}
@@ -555,7 +555,7 @@ void CADSystem::runOndselPiston()
 		auto marker2 = CREATE<MarkerFrame>::With("/Assembly1/Part3/Marker2");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ 0.5, 0.0, 0.0 });
 		marker2->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{0.0, 0.0, 1.0},
 			{1.0, 0.0, 0.0},
 			{0.0, 1.0, 0.0}
@@ -633,7 +633,7 @@ void CADSystem::runPiston()
 	auto assembly1 = CREATE<Part>::With("/Assembly1");
 	std::cout << "assembly1->name " << assembly1->name << std::endl;
 	assembly1->m = 0.0;
-	assembly1->aJ = std::make_shared<DiagonalMatrix<double>>(ListD{ 0, 0, 0 });
+	assembly1->aJ = std::make_shared<DiagonalMatrix>(ListD{ 0, 0, 0 });
 	qX = std::make_shared<FullColumn<double>>(ListD{ 0, 0, 0 });
 	qE = std::make_shared<FullColumn<double>>(ListD{ 0, 0, 0, 1 });
 	assembly1->setqX(qX);
@@ -654,7 +654,7 @@ void CADSystem::runPiston()
 		auto marker2 = CREATE<MarkerFrame>::With("/Assembly1/Marker2");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ 0.0, 0.0, 0.0 });
 		marker2->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{ 1, 0, 0 },
 			{ 0, 1, 0 },
 			{ 0, 0, 1 }
@@ -665,7 +665,7 @@ void CADSystem::runPiston()
 		auto marker1 = CREATE<MarkerFrame>::With("/Assembly1/Marker1");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ 0.0, 2.8817526385684, 0.0 });
 		marker1->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{ 1, 0, 0 },
 			{ 0, 0, 1 },
 			{ 0, -1, 0 }
@@ -678,7 +678,7 @@ void CADSystem::runPiston()
 	auto crankPart1 = CREATE<Part>::With("/Assembly1/Part1");
 	std::cout << "crankPart1->name " << crankPart1->name << std::endl;
 	crankPart1->m = 0.045210530089461;
-	crankPart1->aJ = std::make_shared<DiagonalMatrix<double>>(ListD{ 1.7381980042084e-4, 0.003511159968501, 0.0036154518487535 });
+	crankPart1->aJ = std::make_shared<DiagonalMatrix>(ListD{ 1.7381980042084e-4, 0.003511159968501, 0.0036154518487535 });
 	qX = std::make_shared<FullColumn<double>>(ListD{ 0.38423368514246, 2.6661567755108e-17, -0.048029210642807 });
 	qE = std::make_shared<FullColumn<double>>(ListD{ 0.0, 0.0, 0.0, 1.0 });
 	crankPart1->setqX(qX);
@@ -697,7 +697,7 @@ void CADSystem::runPiston()
 		auto marker1 = CREATE<MarkerFrame>::With("/Assembly1/Part1/Marker1");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ -0.38423368514246, -2.6661567755108e-17, 0.048029210642807 });
 		marker1->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{ 1, 0, 0 },
 			{ 0, 1, 0 },
 			{ 0, 0, 1 }
@@ -708,7 +708,7 @@ void CADSystem::runPiston()
 		auto marker2 = CREATE<MarkerFrame>::With("/Assembly1/Part1/Marker2");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ 0.38423368514246, -2.6661567755108e-17, 0.048029210642807 });
 		marker2->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{ 1, 0, 0 },
 			{ 0, 1, 0 },
 			{ 0, 0, 1 }
@@ -720,7 +720,7 @@ void CADSystem::runPiston()
 	auto conrodPart2 = CREATE<Part>::With("/Assembly1/Part2");
 	std::cout << "conrodPart2->name " << conrodPart2->name << std::endl;
 	conrodPart2->m = 0.067815795134192;
-	conrodPart2->aJ = std::make_shared<DiagonalMatrix<double>>(ListD{ 2.6072970063126e-4, 0.011784982468533, 0.011941420288912 });
+	conrodPart2->aJ = std::make_shared<DiagonalMatrix>(ListD{ 2.6072970063126e-4, 0.011784982468533, 0.011941420288912 });
 	qX = std::make_shared<FullColumn<double>>(ListD{ 0.38423368514246, 0.49215295678475, 0.048029210642807 });
 	qE = std::make_shared<FullColumn<double>>(ListD{ 0.0, 0.0, 0.89871703427292, 0.43852900965351 });
 	conrodPart2->setqX(qX);
@@ -739,7 +739,7 @@ void CADSystem::runPiston()
 		auto marker1 = CREATE<MarkerFrame>::With("/Assembly1/Part2/Marker1");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ -0.6243797383565, 1.1997705489799e-16, -0.048029210642807 });
 		marker1->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{1.0, 2.7755575615629e-16, 0.0},
 			{-2.7755575615629e-16, 1.0, 0.0},
 			{0.0, 0.0, 1.0}
@@ -750,7 +750,7 @@ void CADSystem::runPiston()
 		auto marker2 = CREATE<MarkerFrame>::With("/Assembly1/Part2/Marker2");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ 0.6243797383565, -2.1329254204087e-16, -0.048029210642807 });
 		marker2->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{1.0, 2.4980018054066e-16, 2.2204460492503e-16},
 			{-2.4980018054066e-16, 1.0, 4.1633363423443e-17},
 			{-2.2204460492503e-16, -4.1633363423443e-17, 1.0}
@@ -762,7 +762,7 @@ void CADSystem::runPiston()
 	auto pistonPart3 = CREATE<Part>::With("/Assembly1/Part3");
 	std::cout << "pistonPart3->name " << pistonPart3->name << std::endl;
 	pistonPart3->m = 1.730132083368;
-	pistonPart3->aJ = std::make_shared<DiagonalMatrix<double>>(ListD{ 0.19449049546716, 0.23028116340971, 0.23028116340971 });
+	pistonPart3->aJ = std::make_shared<DiagonalMatrix>(ListD{ 0.19449049546716, 0.23028116340971, 0.23028116340971 });
 	qX = std::make_shared<FullColumn<double>>(ListD{ -1.283972762056e-18, 1.4645980199976, -4.7652385308244e-17 });
 	qE = std::make_shared<FullColumn<double>>(ListD{ 0.70710678118655, 0.70710678118655, 0.0, 0.0 });
 	pistonPart3->setqX(qX);
@@ -781,7 +781,7 @@ void CADSystem::runPiston()
 		auto marker1 = CREATE<MarkerFrame>::With("/Assembly1/Part3/Marker1");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ -0.48029210642807, 7.6201599718927e-18, -2.816737703896e-17 });
 		marker1->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{9.2444637330587e-33, 1.0, 2.2204460492503e-16},
 			{1.0, -9.2444637330587e-33, -1.0785207688569e-32},
 			{-1.0785207688569e-32, 2.2204460492503e-16, -1.0}
@@ -792,7 +792,7 @@ void CADSystem::runPiston()
 		auto marker2 = CREATE<MarkerFrame>::With("/Assembly1/Part3/Marker2");
 		rpmp = std::make_shared<FullColumn<double>>(ListD{ 0.48029210642807, 1.7618247880058e-17, 2.5155758471256e-17 });
 		marker2->setrpmp(rpmp);
-		aApm = std::make_shared<FullMatrix<double>>(ListListD{
+		aApm = std::make_shared<FullMatrixDouble>(ListListD{
 			{6.9388939039072e-18, -6.4146353042213e-50, 1.0},
 			{1.0, -6.9388939039072e-18, 6.9388939039072e-18},
 			{-6.9388939039072e-18, 1.0, -7.4837411882581e-50}

OndselSolver/CMakeLists.txt

@@ -1,570 +0,0 @@
-target_sources(ondsel-solver PUBLIC
-        Abs.cpp
-        AbsConstraint.cpp
-        AccICKineNewtonRaphson.cpp
-        AccICNewtonRaphson.cpp
-        AccKineNewtonRaphson.cpp
-        AccNewtonRaphson.cpp
-        AngleJoint.cpp
-        AngleZIecJec.cpp
-        AngleZIeqcJec.cpp
-        AngleZIeqcJeqc.cpp
-        AnyGeneralSpline.cpp
-        AnyPosICNewtonRaphson.cpp
-        ArcSine.cpp
-        ArcTan.cpp
-        ArcTan2.cpp
-        Array.cpp
-        ASMTAnimationParameters.cpp
-        ASMTAssembly.cpp
-        ASMTConstantGravity.cpp
-        ASMTConstraintSet.cpp
-        ASMTCylindricalJoint.cpp
-        ASMTExtrusion.cpp
-        ASMTFixedJoint.cpp
-        ASMTForceTorque.cpp
-        ASMTGeneralMotion.cpp
-        ASMTItem.cpp
-        ASMTItemIJ.cpp
-        ASMTJoint.cpp
-        ASMTKinematicIJ.cpp
-        ASMTMarker.cpp
-        ASMTMotion.cpp
-        ASMTNoRotationJoint.cpp
-        ASMTPart.cpp
-        ASMTPointInLineJoint.cpp
-        ASMTPointInPlaneJoint.cpp
-        ASMTPrincipalMassMarker.cpp
-        ASMTRefCurve.cpp
-        ASMTRefItem.cpp
-        ASMTRefPoint.cpp
-        ASMTRefSurface.cpp
-        ASMTRevoluteJoint.cpp
-        ASMTRotationalMotion.cpp
-        ASMTSimulationParameters.cpp
-        ASMTSpatialContainer.cpp
-        ASMTSpatialItem.cpp
-        ASMTSphericalJoint.cpp
-        ASMTTime.cpp
-        ASMTTranslationalMotion.cpp
-        ASMTTranslationalJoint.cpp
-        ASMTUniversalJoint.cpp
-        AtPointConstraintIJ.cpp
-        AtPointConstraintIqcJc.cpp
-        AtPointConstraintIqcJqc.cpp
-        AtPointConstraintIqctJqc.cpp
-        AtPointJoint.cpp
-        BasicIntegrator.cpp
-        BasicQuasiIntegrator.cpp
-        BasicUserFunction.cpp
-        CADSystem.cpp
-        CartesianFrame.cpp
-        CompoundJoint.cpp
-        Constant.cpp
-        ConstantGravity.cpp
-        ConstantVelocityJoint.cpp
-        Constraint.cpp
-        ConstraintIJ.cpp
-        ConstVelConstraintIJ.cpp
-        ConstVelConstraintIqcJc.cpp
-        ConstVelConstraintIqcJqc.cpp
-        Cosine.cpp
-        CREATE.cpp
-        CylindricalJoint.cpp
-        CylSphJoint.cpp
-        DifferentiatedGeneralSpline.cpp
-        EndFrameqct2.cpp
-        EulerAngles.cpp
-        EulerAnglesDDot.cpp
-        EulerAnglesDot.cpp
-        Exponential.cpp
-        ExternalSystem.cpp
-        FunctionFromData.cpp
-        FunctionXcParameter.cpp
-        FunctionXY.cpp
-        GeneralSpline.cpp
-        Ln.cpp
-        Log10.cpp
-        LogN.cpp
-        Negative.cpp
-        PosVelAccData.cpp
-        DiagonalMatrix.cpp
-        DifferenceOperator.cpp
-        DirectionCosineConstraintIJ.cpp
-        DirectionCosineConstraintIqcJc.cpp
-        DirectionCosineConstraintIqcJqc.cpp
-        DirectionCosineConstraintIqctJqc.cpp
-        DirectionCosineIecJec.cpp
-        DirectionCosineIeqcJec.cpp
-        DirectionCosineIeqcJeqc.cpp
-        DirectionCosineIeqctJeqc.cpp
-        DiscontinuityError.cpp
-        DispCompIecJecIe.cpp
-        DispCompIecJecKec.cpp
-        DispCompIecJecKeqc.cpp
-        DispCompIecJecO.cpp
-        DispCompIeqcJecIe.cpp
-        DispCompIeqcJecKeqc.cpp
-        DispCompIeqcJecO.cpp
-        DispCompIeqcJeqcIe.cpp
-        DispCompIeqcJeqcKeqc.cpp
-        DispCompIeqcJeqcKeqct.cpp
-        DispCompIeqcJeqcO.cpp
-        DispCompIeqctJeqcIe.cpp
-        DispCompIeqctJeqcKeqct.cpp
-        DispCompIeqctJeqcO.cpp
-        DistanceConstraintIJ.cpp
-        DistanceConstraintIqcJc.cpp
-        DistanceConstraintIqcJqc.cpp
-        DistanceConstraintIqctJqc.cpp
-        DistancexyConstraintIJ.cpp
-        DistancexyConstraintIqcJc.cpp
-        DistancexyConstraintIqcJqc.cpp
-        DistIecJec.cpp
-        DistIeqcJec.cpp
-        DistIeqcJeqc.cpp
-        DistIeqctJeqc.cpp
-        DistxyIecJec.cpp
-        DistxyIeqcJec.cpp
-        DistxyIeqcJeqc.cpp
-        DistxyIeqctJeqc.cpp
-        EndFramec.cpp
-        EndFrameqc.cpp
-        EndFrameqct.cpp
-        EulerAngleszxz.cpp
-        EulerAngleszxzDDot.cpp
-        EulerAngleszxzDot.cpp
-        EulerArray.cpp
-        EulerConstraint.cpp
-        EulerParameters.cpp
-        EulerParametersDDot.cpp
-        EulerParametersDot.cpp
-        ExpressionX.cpp
-        FixedJoint.cpp
-        ForceTorqueData.cpp
-        ForceTorqueItem.cpp
-        FullColumn.cpp
-        FullMatrix.cpp
-        FullMotion.cpp
-        FullRow.cpp
-        Function.cpp
-        FunctionWithManyArgs.cpp
-        FunctionX.cpp
-        GearConstraintIJ.cpp
-        GearConstraintIqcJc.cpp
-        GearConstraintIqcJqc.cpp
-        GearJoint.cpp
-        GEFullMat.cpp
-        GEFullMatFullPv.cpp
-        GEFullMatParPv.cpp
-        GESpMat.cpp
-        GESpMatFullPv.cpp
-        GESpMatFullPvPosIC.cpp
-        GESpMatParPv.cpp
-        GESpMatParPvMarko.cpp
-        GESpMatParPvMarkoFast.cpp
-        GESpMatParPvPrecise.cpp
-        ICKineIntegrator.cpp
-        IndependentVariable.cpp
-        InLineJoint.cpp
-        InPlaneJoint.cpp
-        Integrator.cpp
-        IntegratorInterface.cpp
-        Item.cpp
-        Joint.cpp
-        KineIntegrator.cpp
-        KinematicIeJe.cpp
-        LDUFullMatParPv.cpp
-        LDUSpMat.cpp
-        LDUSpMatParPv.cpp
-        LDUSpMatParPvMarko.cpp
-        LDUSpMatParPvPrecise.cpp
-        LinearMultiStepMethod.cpp
-        LineInPlaneJoint.cpp
-        MarkerFrame.cpp
-        MatrixDecomposition.cpp
-        MatrixGaussElimination.cpp
-        MatrixLDU.cpp
-        MatrixSolver.cpp
-        MaximumIterationError.cpp
-        MBDynBlock.cpp
-        MBDynBody.cpp
-        MBDynControlData.cpp
-        MBDynData.cpp
-        MBDynElement.cpp
-        MBDynInitialValue.cpp
-        MBDynItem.cpp
-        MBDynJoint.cpp
-        MBDynMarker.cpp
-        MBDynNode.cpp
-        MBDynReference.cpp
-        MBDynStructural.cpp
-        MBDynSystem.cpp
-        MomentOfInertiaSolver.cpp
-        NewtonRaphson.cpp
-        NewtonRaphsonError.cpp
-        NoRotationJoint.cpp
-        NotKinematicError.cpp
-        Numeric.cpp
-        OrbitAnglezIecJec.cpp
-        OrbitAnglezIeqcJec.cpp
-        OrbitAnglezIeqcJeqc.cpp
-        Orientation.cpp
-        ParallelAxesJoint.cpp
-        Part.cpp
-        PartFrame.cpp
-        PerpendicularJoint.cpp
-        PlanarJoint.cpp
-        PointInLineJoint.cpp
-        PointInPlaneJoint.cpp
-        PosICKineNewtonRaphson.cpp
-        PosICNewtonRaphson.cpp
-        PosKineNewtonRaphson.cpp
-        PosNewtonRaphson.cpp
-        Power.cpp
-        PrescribedMotion.cpp
-        Product.cpp
-        QuasiIntegrator.cpp
-        RackPinConstraintIJ.cpp
-        RackPinConstraintIqcJc.cpp
-        RackPinConstraintIqcJqc.cpp
-        RackPinJoint.cpp
-        Reciprocal.cpp
-        RedundantConstraint.cpp
-        RevCylJoint.cpp
-        RevoluteJoint.cpp
-        RowTypeMatrix.cpp
-        ScalarNewtonRaphson.cpp
-        ScrewConstraintIJ.cpp
-        ScrewConstraintIqcJc.cpp
-        ScrewConstraintIqcJqc.cpp
-        ScrewJoint.cpp
-        SimulationStoppingError.cpp
-        Sine.cpp
-        SingularMatrixError.cpp
-        Solver.cpp
-        SparseColumn.cpp
-        SparseMatrix.cpp
-        SparseRow.cpp
-        SparseVector.cpp
-        SphericalJoint.cpp
-        SphSphJoint.cpp
-        StableBackwardDifference.cpp
-        LDUFullMat.cpp
-        StateData.cpp
-        Sum.cpp
-        Symbolic.cpp
-        SymbolicParser.cpp
-        SyntaxError.cpp
-        System.cpp
-        SystemNewtonRaphson.cpp
-        SystemSolver.cpp
-        Time.cpp
-        TooManyTriesError.cpp
-        TooSmallStepSizeError.cpp
-        Translation.cpp
-        TranslationalJoint.cpp
-        TranslationConstraintIJ.cpp
-        TranslationConstraintIqcJc.cpp
-        TranslationConstraintIqcJqc.cpp
-        TranslationConstraintIqctJqc.cpp
-        Units.cpp
-        UniversalJoint.cpp
-        UserFunction.cpp
-        Variable.cpp
-        FullVector.cpp
-        VectorNewtonRaphson.cpp
-        VelICKineSolver.cpp
-        VelICSolver.cpp
-        VelKineSolver.cpp
-        VelSolver.cpp
-        ZRotation.cpp
-        ZTranslation.cpp
-        ASMTAssembly.cpp
-)
-
-target_sources(OndselSolver PUBLIC
-        Abs.h
-        AbsConstraint.h
-        AccICKineNewtonRaphson.h
-        AccICNewtonRaphson.h
-        AccKineNewtonRaphson.h
-        AccNewtonRaphson.h
-        AngleJoint.h
-        AngleZIecJec.h
-        AngleZIeqcJec.h
-        AngleZIeqcJeqc.h
-        AnyGeneralSpline.h
-        AnyPosICNewtonRaphson.h
-        ArcSine.h
-        ArcTan.h
-        ArcTan2.h
-        Array.h
-        ASMTAnimationParameters.h
-        ASMTAssembly.h
-        ASMTConstantGravity.h
-        ASMTConstraintSet.h
-        ASMTCylindricalJoint.h
-        ASMTExtrusion.h
-        ASMTFixedJoint.h
-        ASMTForceTorque.h
-        ASMTGeneralMotion.h
-        ASMTItem.h
-        ASMTItemIJ.h
-        ASMTJoint.h
-        ASMTKinematicIJ.h
-        ASMTMarker.h
-        ASMTMotion.h
-        ASMTNoRotationJoint.h
-        ASMTPart.h
-        ASMTPointInLineJoint.h
-        ASMTPointInPlaneJoint.h
-        ASMTPrincipalMassMarker.h
-        ASMTRefCurve.h
-        ASMTRefItem.h
-        ASMTRefPoint.h
-        ASMTRefSurface.h
-        ASMTRevoluteJoint.h
-        ASMTRotationalMotion.h
-        ASMTSimulationParameters.h
-        ASMTSpatialContainer.h
-        ASMTSpatialItem.h
-        ASMTSphericalJoint.h
-        ASMTTime.h
-        ASMTTranslationalMotion.h
-        ASMTTranslationalJoint.h
-        ASMTUniversalJoint.h
-        AtPointConstraintIJ.h
-        AtPointConstraintIqcJc.h
-        AtPointConstraintIqcJqc.h
-        AtPointConstraintIqctJqc.h
-        AtPointJoint.h
-        BasicIntegrator.h
-        BasicQuasiIntegrator.h
-        BasicUserFunction.h
-        CADSystem.h
-        CartesianFrame.h
-        CompoundJoint.h
-        Constant.h
-        ConstantGravity.h
-        ConstantVelocityJoint.h
-        Constraint.h
-        ConstraintIJ.h
-        ConstVelConstraintIJ.h
-        ConstVelConstraintIqcJc.h
-        ConstVelConstraintIqcJqc.h
-        Cosine.h
-        CREATE.h
-        CylindricalJoint.h
-        CylSphJoint.h
-        DifferentiatedGeneralSpline.h
-        EndFrameqct2.h
-        EulerAngles.h
-        EulerAnglesDDot.h
-        EulerAnglesDot.h
-        Exponential.h
-        ExternalSystem.h
-        FunctionFromData.h
-        FunctionXcParameter.h
-        FunctionXY.h
-        GeneralSpline.h
-        Ln.h
-        Log10.h
-        LogN.h
-        Negative.h
-        PosVelAccData.h
-        DiagonalMatrix.h
-        DifferenceOperator.h
-        DirectionCosineConstraintIJ.h
-        DirectionCosineConstraintIqcJc.h
-        DirectionCosineConstraintIqcJqc.h
-        DirectionCosineConstraintIqctJqc.h
-        DirectionCosineIecJec.h
-        DirectionCosineIeqcJec.h
-        DirectionCosineIeqcJeqc.h
-        DirectionCosineIeqctJeqc.h
-        DiscontinuityError.h
-        DispCompIecJecIe.h
-        DispCompIecJecKec.h
-        DispCompIecJecKeqc.h
-        DispCompIecJecO.h
-        DispCompIeqcJecIe.h
-        DispCompIeqcJecKeqc.h
-        DispCompIeqcJecO.h
-        DispCompIeqcJeqcIe.h
-        DispCompIeqcJeqcKeqc.h
-        DispCompIeqcJeqcKeqct.h
-        DispCompIeqcJeqcO.h
-        DispCompIeqctJeqcIe.h
-        DispCompIeqctJeqcKeqct.h
-        DispCompIeqctJeqcO.h
-        DistanceConstraintIJ.h
-        DistanceConstraintIqcJc.h
-        DistanceConstraintIqcJqc.h
-        DistanceConstraintIqctJqc.h
-        DistancexyConstraintIJ.h
-        DistancexyConstraintIqcJc.h
-        DistancexyConstraintIqcJqc.h
-        DistIecJec.h
-        DistIeqcJec.h
-        DistIeqcJeqc.h
-        DistIeqctJeqc.h
-        DistxyIecJec.h
-        DistxyIeqcJec.h
-        DistxyIeqcJeqc.h
-        DistxyIeqctJeqc.h
-        EndFramec.h
-        EndFrameqc.h
-        EndFrameqct.h
-        enum.h
-        EulerAngleszxz.h
-        EulerAngleszxzDDot.h
-        EulerAngleszxzDot.h
-        EulerArray.h
-        EulerConstraint.h
-        EulerParameters.h
-        EulerParametersDDot.h
-        EulerParametersDot.h
-        ExpressionX.h
-        FixedJoint.h
-        ForceTorqueData.h
-        ForceTorqueItem.h
-        FullColumn.h
-        FullMatrix.h
-        FullMotion.h
-        FullRow.h
-        Function.h
-        FunctionWithManyArgs.h
-        FunctionX.h
-        GearConstraintIJ.h
-        GearConstraintIqcJc.h
-        GearConstraintIqcJqc.h
-        GearJoint.h
-        GEFullMat.h
-        GEFullMatFullPv.h
-        GEFullMatParPv.h
-        GESpMat.h
-        GESpMatFullPv.h
-        GESpMatFullPvPosIC.h
-        GESpMatParPv.h
-        GESpMatParPvMarko.h
-        GESpMatParPvMarkoFast.h
-        GESpMatParPvPrecise.h
-        ICKineIntegrator.h
-        IndependentVariable.h
-        InLineJoint.h
-        InPlaneJoint.h
-        Integrator.h
-        IntegratorInterface.h
-        Item.h
-        Joint.h
-        KineIntegrator.h
-        KinematicIeJe.h
-        LDUFullMatParPv.h
-        LDUSpMat.h
-        LDUSpMatParPv.h
-        LDUSpMatParPvMarko.h
-        LDUSpMatParPvPrecise.h
-        LinearMultiStepMethod.h
-        LineInPlaneJoint.h
-        MarkerFrame.h
-        MatrixDecomposition.h
-        MatrixGaussElimination.h
-        MatrixLDU.h
-        MatrixSolver.h
-        MaximumIterationError.h
-        MBDynBlock.h
-        MBDynBody.h
-        MBDynControlData.h
-        MBDynData.h
-        MBDynElement.h
-        MBDynInitialValue.h
-        MBDynItem.h
-        MBDynJoint.h
-        MBDynMarker.h
-        MBDynNode.h
-        MBDynReference.h
-        MBDynStructural.h
-        MBDynSystem.h
-        MomentOfInertiaSolver.h
-        NewtonRaphson.h
-        NewtonRaphsonError.h
-        NoRotationJoint.h
-        NotKinematicError.h
-        Numeric.h
-        OrbitAnglezIecJec.h
-        OrbitAnglezIeqcJec.h
-        OrbitAnglezIeqcJeqc.h
-        Orientation.h
-        ParallelAxesJoint.h
-        Part.h
-        PartFrame.h
-        PerpendicularJoint.h
-        PlanarJoint.h
-        PointInLineJoint.h
-        PointInPlaneJoint.h
-        PosICKineNewtonRaphson.h
-        PosICNewtonRaphson.h
-        PosKineNewtonRaphson.h
-        PosNewtonRaphson.h
-        Power.h
-        PrescribedMotion.h
-        Product.h
-        QuasiIntegrator.h
-        RackPinConstraintIJ.h
-        RackPinConstraintIqcJc.h
-        RackPinConstraintIqcJqc.h
-        RackPinJoint.h
-        Reciprocal.h
-        RedundantConstraint.h
-        resource.h
-        RevCylJoint.h
-        RevoluteJoint.h
-        RowTypeMatrix.h
-        ScalarNewtonRaphson.h
-        ScrewConstraintIJ.h
-        ScrewConstraintIqcJc.h
-        ScrewConstraintIqcJqc.h
-        ScrewJoint.h
-        SimulationStoppingError.h
-        Sine.h
-        SingularMatrixError.h
-        Solver.h
-        SparseColumn.h
-        SparseMatrix.h
-        SparseRow.h
-        SparseVector.h
-        SphericalJoint.h
-        SphSphJoint.h
-        StableBackwardDifference.h
-        LDUFullMat.h
-        StateData.h
-        Sum.h
-        Symbolic.h
-        SymbolicParser.h
-        SyntaxError.h
-        System.h
-        SystemNewtonRaphson.h
-        SystemSolver.h
-        Time.h
-        TooManyTriesError.h
-        TooSmallStepSizeError.h
-        Translation.h
-        TranslationalJoint.h
-        TranslationConstraintIJ.h
-        TranslationConstraintIqcJc.h
-        TranslationConstraintIqcJqc.h
-        TranslationConstraintIqctJqc.h
-        Units.h
-        UniversalJoint.h
-        UserFunction.h
-        Variable.h
-        FullVector.h
-        VectorNewtonRaphson.h
-        VelICKineSolver.h
-        VelICSolver.h
-        VelKineSolver.h
-        VelSolver.h
-        ZRotation.h
-        ZTranslation.h
-)

OndselSolver/CREATE.h

@@ -6,7 +6,7 @@
  *   See LICENSE file for details about copyright.                         *
  ***************************************************************************/
  
-//This header file causes wierd problems in Visual Studio when included in subclasses of std::vector or std::map. Why?
+//This header file causes weird problems in Visual Studio when included in subclasses of std::vector or std::map. Why?
 
 #pragma once
 
@@ -78,7 +78,7 @@ namespace MbD {
 				else {
 					inst = std::make_shared<AtPointConstraintIqcJqc>(frmi, frmj, axis);
 				}
-			} // "class MbD::Tran
+			}
 			else if(str.find("TranslationConstraintIJ") != std::string::npos) {
 				if (std::dynamic_pointer_cast<EndFrameqct>(frmi)) {
 					inst = std::make_shared<TranslationConstraintIqctJqc>(frmi, frmj, axis);

OndselSolver/ClassDiagram.cd

@@ -1,2 +0,0 @@
-<?xml version="1.0" encoding="utf-8"?>
-<ClassDiagram /> 

OndselSolver/ClassDiagram1.cd

@@ -1,11 +0,0 @@
-<?xml version="1.0" encoding="utf-8"?>
-<ClassDiagram MajorVersion="1" MinorVersion="1">
-  <Class Name="MbD::Solver" Collapsed="true">
-    <Position X="0.5" Y="0.5" Width="1.5" />
-    <TypeIdentifier>
-      <HashCode>AKAAAgAAAAAAgAAAEAEAAIABAAAAAAAACIAAIAAAAAI=</HashCode>
-      <FileName>Solver.h</FileName>
-    </TypeIdentifier>
-  </Class>
-  <Font Name="Segoe UI" Size="9" />
-</ClassDiagram>

OndselSolver/ConstVelConstraintIJ.cpp

@@ -52,40 +52,40 @@ void ConstVelConstraintIJ::postInput()
 {
 	aA01IeJe->postInput();
 	aA10IeJe->postInput();
-	Constraint::postInput();
+	ConstraintIJ::postInput();
 }
 
 void ConstVelConstraintIJ::postPosICIteration()
 {
 	aA01IeJe->postPosICIteration();
 	aA10IeJe->postPosICIteration();
-	Item::postPosICIteration();
+	ConstraintIJ::postPosICIteration();
 }
 
 void ConstVelConstraintIJ::preAccIC()
 {
 	aA01IeJe->preAccIC();
 	aA10IeJe->preAccIC();
-	Constraint::preAccIC();
+	ConstraintIJ::preAccIC();
 }
 
 void ConstVelConstraintIJ::prePosIC()
 {
 	aA01IeJe->prePosIC();
 	aA10IeJe->prePosIC();
-	Constraint::prePosIC();
+	ConstraintIJ::prePosIC();
 }
 
 void ConstVelConstraintIJ::preVelIC()
 {
 	aA01IeJe->preVelIC();
 	aA10IeJe->preVelIC();
-	Item::preVelIC();
+	ConstraintIJ::preVelIC();
 }
 
 void ConstVelConstraintIJ::simUpdateAll()
 {
 	aA01IeJe->simUpdateAll();
 	aA10IeJe->simUpdateAll();
-	Item::simUpdateAll();
+	ConstraintIJ::simUpdateAll();
 }

OndselSolver/ConstVelConstraintIqcJc.cpp

@@ -5,7 +5,7 @@
  *                                                                         *
  *   See LICENSE file for details about copyright.                         *
  ***************************************************************************/
- 
+
 #include <memory>
 
 #include "ConstVelConstraintIqcJc.h"
@@ -97,7 +97,7 @@ void MbD::ConstVelConstraintIqcJc::initialize()
 {
 	ConstVelConstraintIJ::initialize();
 	pGpEI = std::make_shared<FullRow<double>>(4);
-	ppGpEIpEI = std::make_shared<FullMatrix<double>>(4, 4);
+	ppGpEIpEI = std::make_shared<FullMatrixDouble>(4, 4);
 }
 
 void MbD::ConstVelConstraintIqcJc::useEquationNumbers()

OndselSolver/ConstVelConstraintIqcJqc.cpp

@@ -118,8 +118,8 @@ void MbD::ConstVelConstraintIqcJqc::initialize()
 {
 	ConstVelConstraintIqcJc::initialize();
 	pGpEJ = std::make_shared<FullRow<double>>(4);
-	ppGpEIpEJ = std::make_shared<FullMatrix<double>>(4, 4);
-	ppGpEJpEJ = std::make_shared<FullMatrix<double>>(4, 4);
+	ppGpEIpEJ = std::make_shared<FullMatrixDouble>(4, 4);
+	ppGpEJpEJ = std::make_shared<FullMatrixDouble>(4, 4);
 }
 
 void MbD::ConstVelConstraintIqcJqc::useEquationNumbers()

OndselSolver/Constant.cpp

@@ -9,6 +9,7 @@
 #include "Constant.h"
 #include "System.h"
 #include "Units.h"
+#include "Polynomial.h"
 
 using namespace MbD;
 
@@ -22,7 +23,23 @@ Constant::Constant(double val) : Variable(val)
 
 Symsptr MbD::Constant::differentiateWRT(Symsptr var)
 {
-	return std::make_shared<Constant>(0.0);
+	return sptrConstant(0.0);
+}
+
+Symsptr MbD::Constant::integrateWRT(Symsptr var)
+{
+	if (value == 0.0) return clonesptr();
+	auto slope = sptrConstant(value);
+	auto intercept = sptrConstant(0.0);
+	auto coeffs = std::make_shared<std::vector<Symsptr>>();
+	coeffs->push_back(intercept);
+	coeffs->push_back(slope);
+	return std::make_shared<Polynomial>(var, coeffs);
+}
+
+Symsptr MbD::Constant::expandUntil(Symsptr sptr, std::shared_ptr<std::unordered_set<Symsptr>> set)
+{
+	return sptr;
 }
 
 bool Constant::isConstant()
@@ -30,11 +47,6 @@ bool Constant::isConstant()
 	return true;
 }
 
-Symsptr MbD::Constant::expandUntil(std::shared_ptr<std::unordered_set<Symsptr>> set)
-{
-	return clonesptr();
-}
-
 Symsptr MbD::Constant::clonesptr()
 {
 	return std::make_shared<Constant>(*this);
@@ -60,6 +72,11 @@ double MbD::Constant::getValue()
 	return value;
 }
 
+double MbD::Constant::getValue(double arg)
+{
+	return value;
+}
+
 std::ostream& Constant::printOn(std::ostream& s) const
 {
 	return s << this->value;

OndselSolver/Constant.h

@@ -17,13 +17,15 @@ namespace MbD {
         Constant();
         Constant(double val);
         Symsptr differentiateWRT(Symsptr var) override;
+        Symsptr integrateWRT(Symsptr var) override;
+        Symsptr expandUntil(Symsptr sptr, std::shared_ptr<std::unordered_set<Symsptr>> set) override;
         bool isConstant() override;
-        Symsptr expandUntil(std::shared_ptr<std::unordered_set<Symsptr>> set) override;
         Symsptr clonesptr() override;
         bool isZero() override;
         bool isOne() override;
         void createMbD(std::shared_ptr<System> mbdSys, std::shared_ptr<Units> mbdUnits) override;
         double getValue() override;
+        double getValue(double arg) override;
 
         std::ostream& printOn(std::ostream& s) const override;
     };

OndselSolver/ConstantGravity.cpp

@@ -15,6 +15,6 @@ using namespace MbD;
 void MbD::ConstantGravity::fillAccICIterError(FColDsptr col)
 {
 	for (auto& part : *(root()->parts)) {
-		col->atiplusFullColumn(part->iqX(), gXYZ->times(part->m));
+		col->atiplusFullColumntimes(part->iqX(), gXYZ, part->m);
 	}
 }

OndselSolver/Constraint.cpp

@@ -104,7 +104,7 @@ void Constraint::setqsudotlam(FColDsptr col)
 
 void Constraint::fillPosICError(FColDsptr col)
 {
-	col->at(iG) += aG;
+	col->atiplusNumber(iG, aG);
 }
 
 void Constraint::removeRedundantConstraints(std::shared_ptr<std::vector<int>> redundantEqnNos)

OndselSolver/Constraint.h

@@ -5,7 +5,7 @@
  *                                                                         *
  *   See LICENSE file for details about copyright.                         *
  ***************************************************************************/
- 
+
 #pragma once
 
 #include <memory>

OndselSolver/Cosine.cpp

@@ -26,8 +26,13 @@ Symsptr MbD::Cosine::differentiateWRTx()
     return std::make_shared<Negative>(std::make_shared<Sine>(xx));
 }
 
+Symsptr MbD::Cosine::copyWith(Symsptr arg)
+{
+	return std::make_shared<Cosine>(arg);
+}
+
 std::ostream& MbD::Cosine::printOn(std::ostream& s) const
 {
-	s << "cos(" << xx << ")";
+	s << "cos(" << *xx << ")";
 	return s;
 }

OndselSolver/Cosine.h

@@ -19,6 +19,7 @@ namespace MbD {
         Cosine(Symsptr arg);
         double getValue() override;
         Symsptr differentiateWRTx() override;
+        Symsptr copyWith(Symsptr arg) override;
 
         std::ostream& printOn(std::ostream& s) const override;
 

OndselSolver/CrankSlider2.mbd

@@ -0,0 +1,297 @@
+-----------------------------------------------------------------------------
+ # [Data Block]
+
+ begin: data;
+     problem: initial value;
+ end: data;
+
+ #-----------------------------------------------------------------------------
+ # [Problem Block]
+
+ begin: initial value;
+     initial time: 0.0;
+     final time: 8.0;
+     time step: 0.01;
+     max iterations: 100;
+     tolerance: 1e-06;
+     derivatives tolerance: 0.0001;
+     derivatives max iterations: 100;
+     derivatives coefficient: auto;
+ end: initial value;
+
+ #-----------------------------------------------------------------------------
+ # [Control Data Block]
+
+ begin: control data;
+     max iterations: 1000;
+     default orientation: euler321;
+     omega rotates: no;
+     print: none;
+     initial stiffness: 1.0, 1.0;
+     structural nodes: 4;
+     rigid bodies: 3;
+     joints: 6;
+ end: control data;
+
+ #-----------------------------------------------------------------------------
+ # [Design Variables]
+
+     #Generic bodies
+
+         #body: 2
+             set: integer body_2 = 2; #body label
+             set: real mass_2 = 3.805252376198168; #mass [kg]
+             set: real volume_2 = 0.0004816775159744516; #volume [m^3]
+
+         #body: 3
+             set: integer body_3 = 3; #body label
+             set: real mass_3 = 15.238784954845523; #mass [kg]
+             set: real volume_3 = 0.0019289601208665218; #volume [m^3]
+
+         #body: 4
+             set: integer body_4 = 4; #body label
+             set: real mass_4 = 2.865603331977783; #mass [kg]
+             set: real volume_4 = 0.0003627345989845295; #volume [m^3]
+
+     #Nodes
+
+         #node: 1
+             set: integer structural_node_1 = 1; #node label
+
+         #node: 2
+             set: integer structural_node_2 = 2; #node label
+
+         #node: 3
+             set: integer structural_node_3 = 3; #node label
+
+         #node: 4
+             set: integer structural_node_4 = 4; #node label
+
+     #Joints
+
+         #joint: 1
+             set: integer joint_1 = 1; #joint label
+
+         #joint: 2
+             set: integer joint_2 = 2; #joint label
+
+         #joint: 3
+             set: integer joint_3 = 3; #joint label
+
+         #joint: 4
+             set: integer joint_4 = 4; #joint label
+
+         #joint: 5
+             set: integer joint_5 = 5; #joint label
+
+         #joint: 6
+             set: integer joint_6 = 6; #joint label
+
+     #Nodes: initial conditions
+
+         #node: 1
+             set: real Px_1 = -0.06210573361337854; #X component of the absolute position [m]
+             set: real Py_1 = 0.048526435375479564; #Y component of the absolute position [m]
+             set: real Pz_1 = -4.033966837940965e-17; #Z component of the absolute position [m]
+
+             set: real Vx_1 = 0.0; #X component of the absolute velocity [m/s]
+             set: real Vy_1 = 0.0; #Y component of the absolute velocity [m/s]
+             set: real Vz_1 = 0.0; #Z component of the absolute velocity [m/s]
+
+             set: real Wx_1 = 0.0; #X component of the absolute angular velocity [rad/s]
+             set: real Wy_1 = 0.0; #Y component of the absolute angular velocity [rad/s]
+             set: real Wz_1 = 0.0; #Z component of the absolute angular velocity [rad/s]
+
+         #node: 2
+             set: real Px_2 = 0.011666006676941875; #X component of the absolute position [m]
+             set: real Py_2 = 0.15999999999999778; #Y component of the absolute position [m]
+             set: real Pz_2 = -1.2084363289349542e-19; #Z component of the absolute position [m]
+
+             set: real Vx_2 = 0.0; #X component of the absolute velocity [m/s]
+             set: real Vy_2 = 0.0; #Y component of the absolute velocity [m/s]
+             set: real Vz_2 = 0.0; #Z component of the absolute velocity [m/s]
+
+             set: real Wx_2 = 0.0; #X component of the absolute angular velocity [rad/s]
+             set: real Wy_2 = 0.0; #Y component of the absolute angular velocity [rad/s]
+             set: real Wz_2 = 0.0; #Z component of the absolute angular velocity [rad/s]
+
+         #node: 3
+             set: real Px_3 = 0.2111281366952498; #X component of the absolute position [m]
+             set: real Py_3 = 0.16; #Y component of the absolute position [m]
+             set: real Pz_3 = -2.0217697810416158e-18; #Z component of the absolute position [m]
+
+             set: real Vx_3 = 0.0; #X component of the absolute velocity [m/s]
+             set: real Vy_3 = 0.0; #Y component of the absolute velocity [m/s]
+             set: real Vz_3 = 0.0; #Z component of the absolute velocity [m/s]
+
+             set: real Wx_3 = 0.0; #X component of the absolute angular velocity [rad/s]
+             set: real Wy_3 = 0.0; #Y component of the absolute angular velocity [rad/s]
+             set: real Wz_3 = 0.0; #Z component of the absolute angular velocity [rad/s]
+
+         #node: 4
+             set: real Px_4 = -0.1812239275015207; #X component of the absolute position [m]
+             set: real Py_4 = 0.16000000169909329; #Y component of the absolute position [m]
+             set: real Pz_4 = -4.340477856936436e-12; #Z component of the absolute position [m]
+
+             set: real Vx_4 = 0.0; #X component of the absolute velocity [m/s]
+             set: real Vy_4 = 0.0; #Y component of the absolute velocity [m/s]
+             set: real Vz_4 = 0.0; #Z component of the absolute velocity [m/s]
+
+             set: real Wx_4 = 0.0; #X component of the absolute angular velocity [rad/s]
+             set: real Wy_4 = 0.0; #Y component of the absolute angular velocity [rad/s]
+             set: real Wz_4 = 0.0; #Z component of the absolute angular velocity [rad/s]
+
+ #-----------------------------------------------------------------------------
+ # [Intermediate Variables]
+
+     #Moments of inertia and relative center of mass
+
+         #body 2: 
+             set: real Ixx_2 = 0.031674420620509; #moment of inertia [kg*m^2]
+             set: real Iyy_2 = 0.029604112147595; #moment of inertia [kg*m^2]
+             set: real Izz_2 = 0.002867529429125; #moment of inertia [kg*m^2]
+
+             set: real Rx_2 = 0.0; #X component of the relative center of mass [m]
+             set: real Ry_2 = 0.0; #Y component of the relative center of mass [m]
+             set: real Rz_2 = 2.0843632893495426e-20; #Z component of the relative center of mass [m]
+
+         #body 3: 
+             set: real Ixx_3 = 0.09813066341583701; #moment of inertia [kg*m^2]
+             set: real Iyy_3 = 0.095433846761275; #moment of inertia [kg*m^2]
+             set: real Izz_3 = 0.077043262824289; #moment of inertia [kg*m^2]
+
+             set: real Rx_3 = 0.0; #X component of the relative center of mass [m]
+             set: real Ry_3 = 0.0; #Y component of the relative center of mass [m]
+             set: real Rz_3 = 2.1769781041615487e-20; #Z component of the relative center of mass [m]
+
+         #body 4: 
+             set: real Ixx_4 = 0.010133521085753; #moment of inertia [kg*m^2]
+             set: real Iyy_4 = 0.006853402672398001; #moment of inertia [kg*m^2]
+             set: real Izz_4 = 0.00669113151275; #moment of inertia [kg*m^2]
+
+             set: real Rx_4 = 0.0; #X component of the relative center of mass [m]
+             set: real Ry_4 = 0.0; #Y component of the relative center of mass [m]
+             set: real Rz_4 = -4.306356366563123e-20; #Z component of the relative center of mass [m]
+
+ #-----------------------------------------------------------------------------
+ # [Nodes Block]
+
+ begin: nodes;
+
+     structural: structural_node_1,
+                 static,
+                 Px_1, Py_1, Pz_1, #<absolute_position> [m]
+                 3, 0.0, 0.0, 1.0, 2, 0.0, 1.0, 0.0, #<absolute_orientation_matrix>
+                 Vx_1, Vy_1, Vz_1, #<absolute_velocity> [m/s]
+                 Wx_1, Wy_1, Wz_1; #<absolute_angular_velocity> [rad/s]
+
+     structural: structural_node_2,
+                 dynamic,
+                 Px_2, Py_2, Pz_2, #<absolute_position> [m]
+                 3, 0.0, 0.0, 1.0, 2, 0.0, 1.0, 0.0, #<absolute_orientation_matrix>
+                 Vx_2, Vy_2, Vz_2, #<absolute_velocity> [m/s]
+                 Wx_2, Wy_2, Wz_2; #<absolute_angular_velocity> [rad/s]
+
+     structural: structural_node_3,
+                 dynamic,
+                 Px_3, Py_3, Pz_3, #<absolute_position> [m]
+                 3, 0.0, 0.0, 1.0, 2, 0.0, 1.0, 0.0, #<absolute_orientation_matrix>
+                 Vx_3, Vy_3, Vz_3, #<absolute_velocity> [m/s]
+                 Wx_3, Wy_3, Wz_3; #<absolute_angular_velocity> [rad/s]
+
+     structural: structural_node_4,
+                 dynamic,
+                 Px_4, Py_4, Pz_4, #<absolute_position> [m]
+                 3, 0.0, 0.0, 1.0, 2, 0.0, 1.0, 0.0, #<absolute_orientation_matrix>
+                 Vx_4, Vy_4, Vz_4, #<absolute_velocity> [m/s]
+                 Wx_4, Wy_4, Wz_4; #<absolute_angular_velocity> [rad/s]
+
+ end: nodes;
+
+ #-----------------------------------------------------------------------------
+ # [Elements Block]
+
+ begin: elements;
+
+     #-----------------------------------------------------------------------------
+     # [Bodies]
+
+         body: body_2,
+               structural_node_2, #<node_label> 
+               mass_2, #<mass> [kg]
+               Rx_2, Ry_2, Rz_2, #<relative_center_of_mass> [m]
+               diag, Ixx_2, Iyy_2, Izz_2, #<inertia matrix> [kg*m^2]
+               orientation, 3, 1.0, 0.0, 0.0, 2, 0.0, 0.9999999999999999, 9.62964972193618e-35; 
+
+         body: body_3,
+               structural_node_3, #<node_label> 
+               mass_3, #<mass> [kg]
+               Rx_3, Ry_3, Rz_3, #<relative_center_of_mass> [m]
+               diag, Ixx_3, Iyy_3, Izz_3, #<inertia matrix> [kg*m^2]
+               orientation, 3, 0.0, 0.0, 1.0, 2, 1.0, 0.0, 0.0; 
+
+         body: body_4,
+               structural_node_4, #<node_label> 
+               mass_4, #<mass> [kg]
+               Rx_4, Ry_4, Rz_4, #<relative_center_of_mass> [m]
+               diag, Ixx_4, Iyy_4, Izz_4, #<inertia matrix> [kg*m^2]
+               orientation, 3, 0.0, 0.0, 1.0, 2, 0.0, 1.0, 0.0; 
+
+     #-----------------------------------------------------------------------------
+     # [Joints]
+
+     joint: joint_1,
+         clamp,
+         structural_node_1, #<node_label>
+             -0.06210573361337854, 0.048526435375479564, -4.033966837940965e-17, #<absolute_pin_position> [m]
+             3, 0.0, 0.0, 1.0, 2, 0.0, 1.0, 0.0; #<absolute_orientation_matrix> 
+
+     joint: joint_2,
+         axial rotation,
+         structural_node_1, #<node_1_label>
+             0.2521057336133785, 0.11147356462452043, 0.13500000000000004, #<relative_offset_1> [m]
+             orientation, 3, 0.0, 0.0, 1.0, 2, -1.0, 2.220446049250313e-16, 0.0, #<relative_orientation_matrix_1>
+         structural_node_3, #<node_2_label>
+             -0.021128136695249794, 0.0, 0.135, #<relative_offset_2> [m]
+             orientation, 3, 0.0, 0.0, 1.0, 2, -1.0, 2.220446049250313e-16, 0.0, #<relative_orientation_matrix_2>
+         string, "model::drive(1, Time)"; #<angular_velocity> [rad/s]
+
+     joint: joint_3,
+         revolute hinge,
+         structural_node_2, #<node_1_label>
+             0.08833399332305812, 2.2168933355715126e-15, 0.065, #<relative_position_1> [m]
+             orientation, 3, 0.0, 0.0, 1.0, 2, -1.0, 2.220446049250313e-16, 0.0, #<relative_pin_orientation_matrix_1>
+         structural_node_3, #<node_2_label>
+            -0.11112813669524979, 0.0, 0.065, #<relative_position_2> [m]
+            orientation, 3, 0.0, 0.0, 1.0, 2, -1.0, 2.220446049250313e-16, 0.0; #<relative_pin_orientation_matrix_2>
+
+     joint: joint_4,
+         prismatic,
+         structural_node_1, #<node_1_label>
+             orientation, 3, 0.0, 0.0, 1.0, 2, 0.0, 1.0, 0.0, #relative_orientation_matrix_1>
+         structural_node_4, #<node_2_label>
+             orientation, 3, 0, 0, 1, 2, 0, 1, 0; #relative_orientation_matrix_2>
+
+     joint: joint_5,
+         in line,
+         structural_node_1, #<node_1_label>
+             0.0, 0.11147355803101042, 2.4118673970779456e-17, #<relative_line_position> [m]
+             3, -1.0, 1.1102230246251565e-16, -2.220446049250313e-16, 2, -1.1102230246251565e-16, 0.0, 1.0, #<relative_orientation>
+         structural_node_4, #<node_2_label>
+             offset, 0.0, -8.292603282173139e-09, 4.340448411165876e-12; #<relative_offset> [m]
+
+     joint: joint_6,
+         in line,
+         structural_node_2, #<node_1_label>
+             -0.1616660066769419, 2.2168933355715126e-15, 0.0, #<relative_line_position> [m]
+             3, 0.0, -0.0, -1.0, 2, -1.0, -2.220446049250313e-16, 0.0, #<relative_orientation>
+         structural_node_4, #<node_2_label>
+             offset, 0.031223927501520705, -1.69909327496498e-09, 7.275957614183426e-15; #<relative_offset> [m]
+
+     #-----------------------------------------------------------------------------
+     # [Drive callers]
+
+     drive caller: 1, name,"drive:1", ramp, 10.0, 0.0, 1.0, 0.0;
+
+ end: elements;
+

OndselSolver/DiagonalMatrix.cpp

@@ -7,5 +7,90 @@
  ***************************************************************************/
  
 #include "DiagonalMatrix.h"
+#include "FullMatrix.h"
 
-using namespace MbD;
+namespace MbD {
+
+    DiagMatDsptr DiagonalMatrix::times(double factor)
+    {
+        auto nrow = (int)this->size();
+        auto answer = std::make_shared<DiagonalMatrix>(nrow);
+        for (int i = 0; i < nrow; i++)
+        {
+            answer->at(i) = this->at(i) * factor;
+        }
+        return answer;
+    }
+    void DiagonalMatrix::atiputDiagonalMatrix(int i, std::shared_ptr<DiagonalMatrix> diagMat)
+    {
+        for (int ii = 0; ii < diagMat->size(); ii++)
+        {
+            this->at(i + ii) = diagMat->at(ii);
+        }
+    }
+    FColsptr<double> DiagonalMatrix::timesFullColumn(FColsptr<double> fullCol)
+    {
+        //"a*b = a(i,j)b(j) sum j."
+
+        auto nrow = (int)this->size();
+        auto answer = std::make_shared<FullColumn<double>>(nrow);
+        for (int i = 0; i < nrow; i++)
+        {
+            answer->at(i) = this->at(i) * fullCol->at(i);
+        }
+        return answer;
+    }
+    FMatDsptr DiagonalMatrix::timesFullMatrix(FMatDsptr fullMat)
+    {
+        auto nrow = (int)this->size();
+        auto answer = std::make_shared<FullMatrixDouble>(nrow);
+        for (int i = 0; i < nrow; i++)
+        {
+            answer->at(i) = fullMat->at(i)->times(this->at(i));
+        }
+        return answer;
+    }
+    double DiagonalMatrix::sumOfSquares()
+    {
+        double sum = 0.0;
+        for (int i = 0; i < this->size(); i++)
+        {
+            double element = this->at(i);
+            sum += element * element;
+        }
+        return sum;
+    }
+    int DiagonalMatrix::numberOfElements()
+    {
+        auto n = (int)this->size();
+        return n * n;
+    }
+    void DiagonalMatrix::zeroSelf()
+    {
+        for (int i = 0; i < this->size(); i++) {
+            this->at(i) = 0.0;
+        }
+    }
+    double DiagonalMatrix::maxMagnitude()
+    {
+        double max = 0.0;
+        for (int i = 0; i < this->size(); i++)
+        {
+            double element = this->at(i);
+            if (element < 0.0) element = -element;
+            if (max < element) max = element;
+        }
+        return max;
+    }
+    std::ostream& DiagonalMatrix::printOn(std::ostream& s) const
+    {
+        s << "DiagMat[";
+        s << this->at(0);
+        for (int i = 1; i < this->size(); i++)
+        {
+            s << ", " << this->at(i);
+        }
+        s << "]";
+        return s;
+    }
+}

OndselSolver/DiagonalMatrix.h

@@ -8,30 +8,26 @@
 
 #pragma once
 
+#include "FullColumn.ref.h"
+#include "FullRow.ref.h"
+#include "FullMatrix.ref.h"
+#include "DiagonalMatrix.ref.h"
 #include "Array.h"
 #include "FullColumn.h"
-#include "FullMatrix.h"
 
 namespace MbD {
-	template<typename T>
-	class DiagonalMatrix;
-	template<typename T>
-	using DiagMatsptr = std::shared_ptr<DiagonalMatrix<T>>;
-	using DiagMatDsptr = std::shared_ptr<DiagonalMatrix<double>>;
-
-	template<typename T>
-	class DiagonalMatrix : public Array<T>
+	class DiagonalMatrix : public Array<double>
 	{
 		//
 	public:
-		DiagonalMatrix() : Array<T>() {}
-		DiagonalMatrix(int count) : Array<T>(count) {}
-		DiagonalMatrix(int count, const T& value) : Array<T>(count, value) {}
-		DiagonalMatrix(std::initializer_list<T> list) : Array<T>{ list } {}
-		void atiputDiagonalMatrix(int i, std::shared_ptr<DiagonalMatrix<T>> diagMat);
-		DiagMatsptr<T> times(T factor);
-		FColsptr<T> timesFullColumn(FColsptr<T> fullCol);
-		FMatsptr<T> timesFullMatrix(FMatsptr<T> fullMat);
+		DiagonalMatrix() : Array<double>() {}
+		explicit DiagonalMatrix(int count) : Array<double>(count) {}
+		DiagonalMatrix(int count, const double& value) : Array<double>(count, value) {}
+		DiagonalMatrix(std::initializer_list<double> list) : Array<double>{ list } {}
+		void atiputDiagonalMatrix(int i, std::shared_ptr<DiagonalMatrix> diagMat);
+		DiagMatDsptr times(double factor);
+		FColsptr<double> timesFullColumn(FColsptr<double> fullCol);
+		FMatDsptr timesFullMatrix(FMatDsptr fullMat);
 		int nrow() {
 			return (int)this->size();
 		}
@@ -44,109 +40,6 @@ namespace MbD {
 		double maxMagnitude() override;
 
 		std::ostream& printOn(std::ostream& s) const override;
-
 	};
-	template<>
-	inline DiagMatDsptr DiagonalMatrix<double>::times(double factor)
-	{
-		auto nrow = (int)this->size();
-		auto answer = std::make_shared<DiagonalMatrix<double>>(nrow);
-		for (int i = 0; i < nrow; i++)
-		{
-			answer->at(i) = this->at(i) * factor;
-		}
-		return answer;
-	}
-	template<typename T>
-	inline void DiagonalMatrix<T>::atiputDiagonalMatrix(int i, std::shared_ptr<DiagonalMatrix<T>> diagMat)
-	{
-		for (int ii = 0; ii < diagMat->size(); ii++)
-		{
-			this->at(i + ii) = diagMat->at(ii);
-		}
-	}
-	template<typename T>
-	inline DiagMatsptr<T> DiagonalMatrix<T>::times(T factor)
-	{
-		assert(false);
-	}
-	template<typename T>
-	inline FColsptr<T> DiagonalMatrix<T>::timesFullColumn(FColsptr<T> fullCol)
-	{
-		//"a*b = a(i,j)b(j) sum j."
-
-		auto nrow = (int)this->size();
-		auto answer = std::make_shared<FullColumn<T>>(nrow);
-		for (int i = 0; i < nrow; i++)
-		{
-			answer->at(i) = this->at(i) * fullCol->at(i);
-		}
-		return answer;
-	}
-	template<typename T>
-	inline FMatsptr<T> DiagonalMatrix<T>::timesFullMatrix(FMatsptr<T> fullMat)
-	{
-		auto nrow = (int)this->size();
-		auto answer = std::make_shared<FullMatrix<T>>(nrow);
-		for (int i = 0; i < nrow; i++)
-		{
-			answer->at(i) = fullMat->at(i)->times(this->at(i));
-		}
-		return answer;
-	}
-	template<>
-	inline double DiagonalMatrix<double>::sumOfSquares()
-	{
-		double sum = 0.0;
-		for (int i = 0; i < this->size(); i++)
-		{
-			double element = this->at(i);
-			sum += element * element;
-		}
-		return sum;
-	}
-	template<typename T>
-	inline int DiagonalMatrix<T>::numberOfElements()
-	{
-		auto n = (int)this->size();
-		return n * n;
-	}
-	template<>
-	inline void DiagonalMatrix<double>::zeroSelf()
-	{
-		for (int i = 0; i < this->size(); i++) {
-			this->at(i) = 0.0;
-		}
-	}
-	template<>
-	inline double DiagonalMatrix<double>::maxMagnitude()
-	{
-		double max = 0.0;
-		for (int i = 0; i < this->size(); i++)
-		{
-			double element = this->at(i);
-			if (element < 0.0) element = -element;
-			if (max < element) max = element;
-		}
-		return max;
-	}
-	template<typename T>
-	inline double DiagonalMatrix<T>::maxMagnitude()
-	{
-		assert(false);
-		return 0.0;
-	}
-	template<typename T>
-	inline std::ostream& DiagonalMatrix<T>::printOn(std::ostream& s) const
-	{
-		s << "DiagMat[";
-		s << this->at(0);
-		for (int i = 1; i < this->size(); i++)
-		{
-			s << ", " << this->at(i);
-		}
-		s << "]";
-		return s;
-	}
 }
 

OndselSolver/DiagonalMatrix.ref.h

@@ -0,0 +1,9 @@
+#pragma once
+
+#include <memory>
+
+namespace MbD {
+    class DiagonalMatrix;
+    using DiagMatsptr = std::shared_ptr<DiagonalMatrix>;
+    using DiagMatDsptr = std::shared_ptr<DiagonalMatrix>;
+}

OndselSolver/DifferenceOperator.cpp

@@ -5,7 +5,7 @@
  *                                                                         *
  *   See LICENSE file for details about copyright.                         *
  ***************************************************************************/
- 
+
 #include <cmath>
 
 #include "DifferenceOperator.h"
@@ -20,7 +20,7 @@ FRowDsptr DifferenceOperator::OneOverFactorials = []() {
 	auto oneOverFactorials = std::make_shared<FullRow<double>>(10);
 	for (int i = 0; i < oneOverFactorials->size(); i++)
 	{
-		oneOverFactorials->at(i) = 1.0 / std::tgamma(i+1);
+		oneOverFactorials->at(i) = 1.0 / std::tgamma(i + 1);
 	}
 	return oneOverFactorials;
 }();
@@ -42,6 +42,12 @@ void DifferenceOperator::calcOperatorMatrix()
 
 void DifferenceOperator::initialize()
 {
+	//Do nothing
+}
+
+void MbD::DifferenceOperator::initializeLocally()
+{
+	assert(false);
 }
 
 void DifferenceOperator::setiStep(int i)
@@ -57,7 +63,7 @@ void DifferenceOperator::setorder(int o)
 void DifferenceOperator::instantiateTaylorMatrix()
 {
 	if (taylorMatrix == nullptr || (taylorMatrix->nrow() != (order + 1))) {
-		taylorMatrix = std::make_shared<FullMatrix<double>>(order + 1, order + 1);
+		taylorMatrix = std::make_shared<FullMatrixDouble>(order + 1, order + 1);
 	}
 }
 
@@ -75,9 +81,8 @@ void DifferenceOperator::formTaylorRowwithTimeNodederivative(int i, int ii, int
 	for (int j = k + 1; j < order + 1; j++)
 	{
 		hipower = hipower * hi;
-		assert(false);
-		//aij = hipower * (OneOverFactorials at : j - k - 1);
-		//rowi at : j put : aij
+		auto aij = hipower * OneOverFactorials->at((size_t)j - k);
+		rowi->atiput(j, aij);
 	}
 }
 

OndselSolver/DifferenceOperator.h

@@ -19,6 +19,7 @@ namespace MbD {
 	public:
 		void calcOperatorMatrix();
 		virtual void initialize();
+		virtual void initializeLocally();
 		virtual void setiStep(int i);
 		virtual void setorder(int o);
 		virtual void formTaylorMatrix() = 0;
@@ -29,7 +30,7 @@ namespace MbD {
 		int iStep = 0, order = 0;
 		FMatDsptr taylorMatrix, operatorMatrix;
 		double time = 0.0;
-		std::shared_ptr<std::vector<double>> timeNodes;
+		std::shared_ptr<std::vector<double>> timeNodes;	//"Row of past times in order of increasing past."
 		static FRowDsptr OneOverFactorials;
 	};
 }

OndselSolver/DifferentiatedGeneralSpline.cpp

@@ -34,6 +34,6 @@ Symsptr MbD::DifferentiatedGeneralSpline::clonesptr()
 
 std::ostream& MbD::DifferentiatedGeneralSpline::printOn(std::ostream& s) const
 {
-	s << "deriv(" << generalSpline << ", " << derivativeOrder << ")";
+	s << "deriv(" << *generalSpline << ", " << derivativeOrder << ")";
 	return s;
 }

OndselSolver/DirectionCosineConstraintIJ.cpp

@@ -59,7 +59,7 @@ void DirectionCosineConstraintIJ::prePosIC()
 void DirectionCosineConstraintIJ::postPosICIteration()
 {
 	aAijIeJe->postPosICIteration();
-	Item::postPosICIteration();
+	ConstraintIJ::postPosICIteration();
 }
 
 ConstraintType DirectionCosineConstraintIJ::type()
@@ -70,7 +70,7 @@ ConstraintType DirectionCosineConstraintIJ::type()
 void DirectionCosineConstraintIJ::preVelIC()
 {
 	aAijIeJe->preVelIC();
-	Item::preVelIC();
+	ConstraintIJ::preVelIC();
 }
 
 void MbD::DirectionCosineConstraintIJ::simUpdateAll()

OndselSolver/DirectionCosineIeqcJec.cpp

@@ -24,7 +24,7 @@ void DirectionCosineIeqcJec::initialize()
 {
 	DirectionCosineIecJec::initialize();
 	pAijIeJepEI = std::make_shared<FullRow<double>>(4);
-	ppAijIeJepEIpEI = std::make_shared<FullMatrix<double>>(4, 4);
+	ppAijIeJepEIpEI = std::make_shared<FullMatrixDouble>(4, 4);
 }
 
 void DirectionCosineIeqcJec::initializeGlobally()

OndselSolver/DirectionCosineIeqcJeqc.cpp

@@ -24,8 +24,8 @@ void DirectionCosineIeqcJeqc::initialize()
 {
 	DirectionCosineIeqcJec::initialize();
 	pAijIeJepEJ = std::make_shared<FullRow<double>>(4);
-	ppAijIeJepEIpEJ = std::make_shared<FullMatrix<double>>(4, 4);
-	ppAijIeJepEJpEJ = std::make_shared<FullMatrix<double>>(4, 4);
+	ppAijIeJepEIpEJ = std::make_shared<FullMatrixDouble>(4, 4);
+	ppAijIeJepEJpEJ = std::make_shared<FullMatrixDouble>(4, 4);
 }
 
 void DirectionCosineIeqcJeqc::initializeGlobally()

OndselSolver/DispCompIecJecKeqc.cpp

@@ -22,7 +22,7 @@ DispCompIecJecKeqc::DispCompIecJecKeqc(EndFrmsptr frmi, EndFrmsptr frmj, EndFrms
 void DispCompIecJecKeqc::initialize()
 {
     priIeJeKepEK = std::make_shared<FullRow<double>>(4);
-    ppriIeJeKepEKpEK = std::make_shared<FullMatrix<double>>(4, 4);
+    ppriIeJeKepEKpEK = std::make_shared<FullMatrixDouble>(4, 4);
 }
 
 void DispCompIecJecKeqc::initializeGlobally()

OndselSolver/DispCompIeqcJecIe.cpp

@@ -87,8 +87,8 @@ void MbD::DispCompIeqcJecIe::initialize()
 	DispCompIecJecIe::initialize();
 	priIeJeIepXI = std::make_shared<FullRow<double>>(3);
 	priIeJeIepEI = std::make_shared<FullRow<double>>(4);
-	ppriIeJeIepXIpEI = std::make_shared<FullMatrix<double>>(3, 4);
-	ppriIeJeIepEIpEI = std::make_shared<FullMatrix<double>>(4, 4);
+	ppriIeJeIepXIpEI = std::make_shared<FullMatrixDouble>(3, 4);
+	ppriIeJeIepEIpEI = std::make_shared<FullMatrixDouble>(4, 4);
 }
 
 void MbD::DispCompIeqcJecIe::initializeGlobally()

OndselSolver/DispCompIeqcJecKeqc.cpp

@@ -24,9 +24,9 @@ void DispCompIeqcJecKeqc::initialize()
 	DispCompIecJecKeqc::initialize();
 	priIeJeKepXI = std::make_shared<FullRow<double>>(3);
 	priIeJeKepEI = std::make_shared<FullRow<double>>(4);
-	ppriIeJeKepEIpEI = std::make_shared<FullMatrix<double>>(4, 4);
-	ppriIeJeKepXIpEK = std::make_shared<FullMatrix<double>>(3, 4);
-	ppriIeJeKepEIpEK = std::make_shared<FullMatrix<double>>(4, 4);
+	ppriIeJeKepEIpEI = std::make_shared<FullMatrixDouble>(4, 4);
+	ppriIeJeKepXIpEK = std::make_shared<FullMatrixDouble>(3, 4);
+	ppriIeJeKepEIpEK = std::make_shared<FullMatrixDouble>(4, 4);
 }
 
 void DispCompIeqcJecKeqc::calcPostDynCorrectorIteration()

OndselSolver/DispCompIeqcJeqcIe.cpp

@@ -83,9 +83,9 @@ void MbD::DispCompIeqcJeqcIe::initialize()
 	DispCompIeqcJecIe::initialize();
 	priIeJeIepXJ = std::make_shared<FullRow<double>>(3);
 	priIeJeIepEJ = std::make_shared<FullRow<double>>(4);
-	ppriIeJeIepEIpXJ = std::make_shared<FullMatrix<double>>(4, 3);
-	ppriIeJeIepEIpEJ = std::make_shared<FullMatrix<double>>(4, 4);
-	ppriIeJeIepEJpEJ = std::make_shared<FullMatrix<double>>(4, 4);
+	ppriIeJeIepEIpXJ = std::make_shared<FullMatrixDouble>(4, 3);
+	ppriIeJeIepEIpEJ = std::make_shared<FullMatrixDouble>(4, 4);
+	ppriIeJeIepEJpEJ = std::make_shared<FullMatrixDouble>(4, 4);
 }
 
 FMatDsptr MbD::DispCompIeqcJeqcIe::ppvaluepEIpEJ()

OndselSolver/DispCompIeqcJeqcKeqc.cpp

@@ -24,9 +24,9 @@ void DispCompIeqcJeqcKeqc::initialize()
 	DispCompIeqcJecKeqc::initialize();
 	priIeJeKepXJ = std::make_shared<FullRow<double>>(3);
 	priIeJeKepEJ = std::make_shared<FullRow<double>>(4);
-	ppriIeJeKepEJpEJ = std::make_shared<FullMatrix<double>>(4, 4);
-	ppriIeJeKepXJpEK = std::make_shared<FullMatrix<double>>(3, 4);
-	ppriIeJeKepEJpEK = std::make_shared<FullMatrix<double>>(4, 4);
+	ppriIeJeKepEJpEJ = std::make_shared<FullMatrixDouble>(4, 4);
+	ppriIeJeKepXJpEK = std::make_shared<FullMatrixDouble>(3, 4);
+	ppriIeJeKepEJpEK = std::make_shared<FullMatrixDouble>(4, 4);
 }
 
 void DispCompIeqcJeqcKeqc::calcPostDynCorrectorIteration()

OndselSolver/DistIeqcJec.cpp

@@ -41,7 +41,7 @@ void MbD::DistIeqcJec::calcPrivate()
 			pprIeJepXIipXI->atiput(j, element / rIeJe);
 		}
 	}
-	pprIeJepXIpEI = std::make_shared<FullMatrix<double>>(3, 4);
+	pprIeJepXIpEI = std::make_shared<FullMatrixDouble>(3, 4);
 	for (int i = 0; i < 3; i++)
 	{
 		auto& pprIeJepXIipEI = pprIeJepXIpEI->at(i);
@@ -53,7 +53,7 @@ void MbD::DistIeqcJec::calcPrivate()
 			pprIeJepXIipEI->atiput(j, element / rIeJe);
 		}
 	}
-	pprIeJepEIpEI = std::make_shared<FullMatrix<double>>(4, 4);
+	pprIeJepEIpEI = std::make_shared<FullMatrixDouble>(4, 4);
 	for (int i = 0; i < 4; i++)
 	{
 		auto& pprIeJepEIipEI = pprIeJepEIpEI->at(i);
@@ -74,9 +74,9 @@ void MbD::DistIeqcJec::initialize()
 	DistIecJec::initialize();
 	prIeJepXI = std::make_shared<FullRow<double>>(3);
 	prIeJepEI = std::make_shared<FullRow<double>>(4);
-	pprIeJepXIpXI = std::make_shared<FullMatrix<double>>(3, 3);
-	pprIeJepXIpEI = std::make_shared<FullMatrix<double>>(3, 4);
-	pprIeJepEIpEI = std::make_shared<FullMatrix<double>>(4, 4);
+	pprIeJepXIpXI = std::make_shared<FullMatrixDouble>(3, 3);
+	pprIeJepXIpEI = std::make_shared<FullMatrixDouble>(3, 4);
+	pprIeJepEIpEI = std::make_shared<FullMatrixDouble>(4, 4);
 }
 
 FMatDsptr MbD::DistIeqcJec::ppvaluepEIpEI()

OndselSolver/DistIeqcJeqc.cpp

@@ -122,13 +122,13 @@ void MbD::DistIeqcJeqc::initialize()
 	DistIeqcJec::initialize();
 	prIeJepXJ = std::make_shared<FullRow<double>>(3);
 	prIeJepEJ = std::make_shared<FullRow<double>>(4);
-	pprIeJepXIpXJ = std::make_shared<FullMatrix<double>>(3, 3);
-	pprIeJepEIpXJ = std::make_shared<FullMatrix<double>>(4, 3);
-	pprIeJepXJpXJ = std::make_shared<FullMatrix<double>>(3, 3);
-	pprIeJepXIpEJ = std::make_shared<FullMatrix<double>>(3, 4);
-	pprIeJepEIpEJ = std::make_shared<FullMatrix<double>>(4, 4);
-	pprIeJepXJpEJ = std::make_shared<FullMatrix<double>>(3, 4);
-	pprIeJepEJpEJ = std::make_shared<FullMatrix<double>>(4, 4);
+	pprIeJepXIpXJ = std::make_shared<FullMatrixDouble>(3, 3);
+	pprIeJepEIpXJ = std::make_shared<FullMatrixDouble>(4, 3);
+	pprIeJepXJpXJ = std::make_shared<FullMatrixDouble>(3, 3);
+	pprIeJepXIpEJ = std::make_shared<FullMatrixDouble>(3, 4);
+	pprIeJepEIpEJ = std::make_shared<FullMatrixDouble>(4, 4);
+	pprIeJepXJpEJ = std::make_shared<FullMatrixDouble>(3, 4);
+	pprIeJepEJpEJ = std::make_shared<FullMatrixDouble>(4, 4);
 }
 
 FMatDsptr MbD::DistIeqcJeqc::ppvaluepEIpEJ()

OndselSolver/DistxyIeqcJec.cpp

@@ -161,9 +161,9 @@ void MbD::DistxyIeqcJec::initialize()
 	DistxyIecJec::initialize();
 	pdistxypXI = std::make_shared<FullRow<double>>(3);
 	pdistxypEI = std::make_shared<FullRow<double>>(4);
-	ppdistxypXIpXI = std::make_shared<FullMatrix<double>>(3, 3);
-	ppdistxypXIpEI = std::make_shared<FullMatrix<double>>(3, 4);
-	ppdistxypEIpEI = std::make_shared<FullMatrix<double>>(4, 4);
+	ppdistxypXIpXI = std::make_shared<FullMatrixDouble>(3, 3);
+	ppdistxypXIpEI = std::make_shared<FullMatrixDouble>(3, 4);
+	ppdistxypEIpEI = std::make_shared<FullMatrixDouble>(4, 4);
 }
 
 FMatDsptr MbD::DistxyIeqcJec::ppvaluepEIpEI()

OndselSolver/DistxyIeqcJeqc.cpp

@@ -293,13 +293,13 @@ void MbD::DistxyIeqcJeqc::initialize()
 	DistxyIeqcJec::initialize();
 	pdistxypXJ = std::make_shared<FullRow<double>>(3);
 	pdistxypEJ = std::make_shared<FullRow<double>>(4);
-	ppdistxypXIpXJ = std::make_shared<FullMatrix<double>>(3, 3);
-	ppdistxypXIpEJ = std::make_shared<FullMatrix<double>>(3, 4);
-	ppdistxypEIpXJ = std::make_shared<FullMatrix<double>>(4, 3);
-	ppdistxypEIpEJ = std::make_shared<FullMatrix<double>>(4, 4);
-	ppdistxypXJpXJ = std::make_shared<FullMatrix<double>>(3, 3);
-	ppdistxypXJpEJ = std::make_shared<FullMatrix<double>>(3, 4);
-	ppdistxypEJpEJ = std::make_shared<FullMatrix<double>>(4, 4);
+	ppdistxypXIpXJ = std::make_shared<FullMatrixDouble>(3, 3);
+	ppdistxypXIpEJ = std::make_shared<FullMatrixDouble>(3, 4);
+	ppdistxypEIpXJ = std::make_shared<FullMatrixDouble>(4, 3);
+	ppdistxypEIpEJ = std::make_shared<FullMatrixDouble>(4, 4);
+	ppdistxypXJpXJ = std::make_shared<FullMatrixDouble>(3, 3);
+	ppdistxypXJpEJ = std::make_shared<FullMatrixDouble>(3, 4);
+	ppdistxypEJpEJ = std::make_shared<FullMatrixDouble>(4, 4);
 }
 
 FMatDsptr MbD::DistxyIeqcJeqc::ppvaluepEIpEJ()

OndselSolver/EndFramec.h

@@ -44,7 +44,7 @@ namespace MbD {
 
 		MarkerFrame* markerFrame; //Use raw pointer when pointing backwards.
 		FColDsptr rOeO = std::make_shared<FullColumn<double>>(3);
-		FMatDsptr aAOe = FullMatrix<double>::identitysptr(3);
+		FMatDsptr aAOe = FullMatrixDouble::identitysptr(3);
 	};
 	//using EndFrmsptr = std::shared_ptr<EndFramec>;
 }

OndselSolver/EndFrameqc.cpp

@@ -25,10 +25,10 @@ EndFrameqc::EndFrameqc(const char* str) : EndFramec(str) {
 
 void EndFrameqc::initialize()
 {
-	prOeOpE = std::make_shared<FullMatrix<double>>(3, 4);
-	pprOeOpEpE = std::make_shared<FullMatrix<FColDsptr>>(4, 4);
+	prOeOpE = std::make_shared<FullMatrixDouble>(3, 4);
+	pprOeOpEpE = std::make_shared<FullMatrixFullColumnDouble>(4, 4);
 	pAOepE = std::make_shared<FullColumn<FMatDsptr>>(4);
-	ppAOepEpE = std::make_shared<FullMatrix<FMatDsptr>>(4, 4);
+	ppAOepEpE = std::make_shared<FullMatrixFullMatrixDouble>(4, 4);
 }
 
 void EndFrameqc::initializeGlobally()
@@ -59,7 +59,7 @@ void MbD::EndFrameqc::initEndFrameqct2()
 
 FMatFColDsptr EndFrameqc::ppAjOepEpE(int jj)
 {
-	auto answer = std::make_shared<FullMatrix<FColDsptr>>(4, 4);
+	auto answer = std::make_shared<FullMatrixFullColumnDouble>(4, 4);
 	for (int i = 0; i < 4; i++) {
 		auto& answeri = answer->at(i);
 		auto& ppAOepEipE = ppAOepEpE->at(i);
@@ -80,7 +80,7 @@ void EndFrameqc::calcPostDynCorrectorIteration()
 
 FMatDsptr EndFrameqc::pAjOepET(int axis)
 {
-	auto answer = std::make_shared<FullMatrix<double>>(4, 3);
+	auto answer = std::make_shared<FullMatrixDouble>(4, 3);
 	for (int i = 0; i < 4; i++) {
 		auto& answeri = answer->at(i);
 		auto& pAOepEi = pAOepE->at(i);
@@ -94,7 +94,7 @@ FMatDsptr EndFrameqc::pAjOepET(int axis)
 
 FMatDsptr EndFrameqc::ppriOeOpEpE(int ii)
 {
-	auto answer = std::make_shared<FullMatrix<double>>(4, 4);
+	auto answer = std::make_shared<FullMatrixDouble>(4, 4);
 	for (int i = 0; i < 4; i++) {
 		auto& answeri = answer->at(i);
 		auto& pprOeOpEipE = pprOeOpEpE->at(i);

OndselSolver/EndFrameqct.cpp

@@ -9,7 +9,6 @@
 #include "EndFrameqct.h"
 #include "MarkerFrame.h"
 #include "System.h"
-#include "Symbolic.h"
 #include "Time.h"
 #include "EulerParameters.h"
 #include "CREATE.h"
@@ -17,347 +16,308 @@
 #include "EulerAngleszxzDot.h"
 #include "EulerAngleszxzDDot.h"
 
-using namespace MbD;
+namespace MbD {
+    template class EulerParameters<double>;
 
-EndFrameqct::EndFrameqct() {
-}
+    EndFrameqct::EndFrameqct() {
+    }
 
-EndFrameqct::EndFrameqct(const char* str) : EndFrameqc(str) {
-}
+    EndFrameqct::EndFrameqct(const char *str) : EndFrameqc(str) {
+    }
 
-void EndFrameqct::initialize()
-{
-	EndFrameqc::initialize();
-	rmem = std::make_shared<FullColumn<double>>(3);
-	prmempt = std::make_shared<FullColumn<double>>(3);
-	pprmemptpt = std::make_shared<FullColumn<double>>(3);
-	aAme = FullMatrix<double>::identitysptr(3);
-	pAmept = std::make_shared<FullMatrix<double>>(3, 3);
-	ppAmeptpt = std::make_shared<FullMatrix<double>>(3, 3);
-	pprOeOpEpt = std::make_shared<FullMatrix<double>>(3, 4);
-	pprOeOptpt = std::make_shared<FullColumn<double>>(3);
-	ppAOepEpt = std::make_shared<FullColumn<FMatDsptr>>(4);
-	ppAOeptpt = std::make_shared<FullMatrix<double>>(3, 3);
-}
+    void EndFrameqct::initialize() {
+        EndFrameqc::initialize();
+        rmem = std::make_shared<FullColumn<double>>(3);
+        prmempt = std::make_shared<FullColumn<double>>(3);
+        pprmemptpt = std::make_shared<FullColumn<double>>(3);
+        aAme = FullMatrixDouble::identitysptr(3);
+        pAmept = std::make_shared<FullMatrixDouble>(3, 3);
+        ppAmeptpt = std::make_shared<FullMatrixDouble>(3, 3);
+        pprOeOpEpt = std::make_shared<FullMatrixDouble>(3, 4);
+        pprOeOptpt = std::make_shared<FullColumn<double>>(3);
+        ppAOepEpt = std::make_shared<FullColumn<FMatDsptr>>(4);
+        ppAOeptpt = std::make_shared<FullMatrixDouble>(3, 3);
+    }
 
-void EndFrameqct::initializeLocally()
-{
-	if (!rmemBlks) {
-		rmem->zeroSelf();
-		prmempt->zeroSelf();
-		pprmemptpt->zeroSelf();
-	}
-	if (!phiThePsiBlks) {
-		aAme->identity();
-		pAmept->zeroSelf();
-		ppAmeptpt->zeroSelf();
-	}
-}
+    void EndFrameqct::initializeLocally() {
+        if (!rmemBlks) {
+            rmem->zeroSelf();
+            prmempt->zeroSelf();
+            pprmemptpt->zeroSelf();
+        }
+        if (!phiThePsiBlks) {
+            aAme->identity();
+            pAmept->zeroSelf();
+            ppAmeptpt->zeroSelf();
+        }
+    }
 
-void EndFrameqct::initializeGlobally()
-{
-	if (rmemBlks) {
-		initprmemptBlks();
-		initpprmemptptBlks();
-	}
-	if (phiThePsiBlks) {
-		initpPhiThePsiptBlks();
-		initppPhiThePsiptptBlks();
-	}
-}
+    void EndFrameqct::initializeGlobally() {
+        if (rmemBlks) {
+            initprmemptBlks();
+            initpprmemptptBlks();
+        }
+        if (phiThePsiBlks) {
+            initpPhiThePsiptBlks();
+            initppPhiThePsiptptBlks();
+        }
+    }
 
-void EndFrameqct::initprmemptBlks()
-{
-	auto& mbdTime = this->root()->time;
-	prmemptBlks = std::make_shared< FullColumn<Symsptr>>(3);
-	for (int i = 0; i < 3; i++) {
-		auto& disp = rmemBlks->at(i);
-		auto var = disp->differentiateWRT(mbdTime);
-		auto vel = var->simplified(var);
-		prmemptBlks->at(i) = vel;
-	}
-}
+    void EndFrameqct::initprmemptBlks() {
+        auto &mbdTime = this->root()->time;
+        prmemptBlks = std::make_shared<FullColumn<Symsptr>>(3);
+        for (int i = 0; i < 3; i++) {
+            auto &disp = rmemBlks->at(i);
+            auto var = disp->differentiateWRT(mbdTime);
+            auto vel = var->simplified(var);
+            prmemptBlks->at(i) = vel;
+        }
+    }
 
-void EndFrameqct::initpprmemptptBlks()
-{
-	auto& mbdTime = this->root()->time;
-	pprmemptptBlks = std::make_shared< FullColumn<Symsptr>>(3);
-	for (int i = 0; i < 3; i++) {
-		auto& vel = prmemptBlks->at(i);
-		auto var = vel->differentiateWRT(mbdTime);
-		auto acc = var->simplified(var);
-		pprmemptptBlks->at(i) = acc;
-	}
-}
+    void EndFrameqct::initpprmemptptBlks() {
+        auto &mbdTime = this->root()->time;
+        pprmemptptBlks = std::make_shared<FullColumn<Symsptr>>(3);
+        for (int i = 0; i < 3; i++) {
+            auto &vel = prmemptBlks->at(i);
+            auto var = vel->differentiateWRT(mbdTime);
+            auto acc = var->simplified(var);
+            pprmemptptBlks->at(i) = acc;
+        }
+    }
 
-void EndFrameqct::initpPhiThePsiptBlks()
-{
-	auto& mbdTime = this->root()->time;
-	pPhiThePsiptBlks = std::make_shared< FullColumn<Symsptr>>(3);
-	for (int i = 0; i < 3; i++) {
-		auto& angle = phiThePsiBlks->at(i);
-		auto var = angle->differentiateWRT(mbdTime);
-		//std::cout << "var " << *var << std::endl;
-		auto vel = var->simplified(var);
-		//std::cout << "vel " << *vel << std::endl;
-		pPhiThePsiptBlks->at(i) = vel;
-	}
-}
+    void EndFrameqct::initpPhiThePsiptBlks() {
+        auto &mbdTime = this->root()->time;
+        pPhiThePsiptBlks = std::make_shared<FullColumn<Symsptr>>(3);
+        for (int i = 0; i < 3; i++) {
+            auto &angle = phiThePsiBlks->at(i);
+            auto var = angle->differentiateWRT(mbdTime);
+            //std::cout << "var " << *var << std::endl;
+            auto vel = var->simplified(var);
+            //std::cout << "vel " << *vel << std::endl;
+            pPhiThePsiptBlks->at(i) = vel;
+        }
+    }
 
-void EndFrameqct::initppPhiThePsiptptBlks()
-{
-	auto& mbdTime = this->root()->time;
-	ppPhiThePsiptptBlks = std::make_shared< FullColumn<Symsptr>>(3);
-	for (int i = 0; i < 3; i++) {
-		auto& angleVel = pPhiThePsiptBlks->at(i);
-		auto var = angleVel->differentiateWRT(mbdTime);
-		auto angleAcc = var->simplified(var);
-		ppPhiThePsiptptBlks->at(i) = angleAcc;
-	}
-}
+    void EndFrameqct::initppPhiThePsiptptBlks() {
+        auto &mbdTime = this->root()->time;
+        ppPhiThePsiptptBlks = std::make_shared<FullColumn<Symsptr>>(3);
+        for (int i = 0; i < 3; i++) {
+            auto &angleVel = pPhiThePsiptBlks->at(i);
+            auto var = angleVel->differentiateWRT(mbdTime);
+            auto angleAcc = var->simplified(var);
+            ppPhiThePsiptptBlks->at(i) = angleAcc;
+        }
+    }
 
-void EndFrameqct::postInput()
-{
-	this->evalrmem();
-	this->evalAme();
-	Item::postInput();
-}
+    void EndFrameqct::postInput() {
+        this->evalrmem();
+        this->evalAme();
+        Item::postInput();
+    }
 
-void EndFrameqct::calcPostDynCorrectorIteration()
-{
-	auto& rOmO = markerFrame->rOmO;
-	auto& aAOm = markerFrame->aAOm;
-	rOeO = rOmO->plusFullColumn(aAOm->timesFullColumn(rmem));
-	auto& prOmOpE = markerFrame->prOmOpE;
-	auto& pAOmpE = markerFrame->pAOmpE;
-	for (int i = 0; i < 3; i++)
-	{
-		auto& prOmOpEi = prOmOpE->at(i);
-		auto& prOeOpEi = prOeOpE->at(i);
-		for (int j = 0; j < 4; j++)
-		{
-			auto prOeOpEij = prOmOpEi->at(j) + pAOmpE->at(j)->at(i)->timesFullColumn(rmem);
-			prOeOpEi->at(j) = prOeOpEij;
-		}
-	}
-	auto rpep = markerFrame->rpmp->plusFullColumn(markerFrame->aApm->timesFullColumn(rmem));
-	pprOeOpEpE = EulerParameters<double>::ppApEpEtimesColumn(rpep);
-	aAOe = aAOm->timesFullMatrix(aAme);
-	for (int i = 0; i < 4; i++)
-	{
-		pAOepE->at(i) = pAOmpE->at(i)->timesFullMatrix(aAme);
-	}
-	auto aApe = markerFrame->aApm->timesFullMatrix(aAme);
-	ppAOepEpE = EulerParameters<double>::ppApEpEtimesMatrix(aApe);
-}
+    void EndFrameqct::calcPostDynCorrectorIteration() {
+        auto &rOmO = markerFrame->rOmO;
+        auto &aAOm = markerFrame->aAOm;
+        rOeO = rOmO->plusFullColumn(aAOm->timesFullColumn(rmem));
+        auto &prOmOpE = markerFrame->prOmOpE;
+        auto &pAOmpE = markerFrame->pAOmpE;
+        for (int i = 0; i < 3; i++) {
+            auto &prOmOpEi = prOmOpE->at(i);
+            auto &prOeOpEi = prOeOpE->at(i);
+            for (int j = 0; j < 4; j++) {
+                auto prOeOpEij = prOmOpEi->at(j) + pAOmpE->at(j)->at(i)->timesFullColumn(rmem);
+                prOeOpEi->at(j) = prOeOpEij;
+            }
+        }
+        auto rpep = markerFrame->rpmp->plusFullColumn(markerFrame->aApm->timesFullColumn(rmem));
+        pprOeOpEpE = EulerParameters<double>::ppApEpEtimesColumn(rpep);
+        aAOe = aAOm->timesFullMatrix(aAme);
+        for (int i = 0; i < 4; i++) {
+            pAOepE->at(i) = pAOmpE->at(i)->timesFullMatrix(aAme);
+        }
+        auto aApe = markerFrame->aApm->timesFullMatrix(aAme);
+        ppAOepEpE = EulerParameters<double>::ppApEpEtimesMatrix(aApe);
+    }
 
-void EndFrameqct::prePosIC()
-{
-	time = this->root()->mbdTimeValue();
-	this->evalrmem();
-	this->evalAme();
-	EndFrameqc::prePosIC();
-}
+    void EndFrameqct::prePosIC() {
+        time = this->root()->mbdTimeValue();
+        this->evalrmem();
+        this->evalAme();
+        EndFrameqc::prePosIC();
+    }
 
-void EndFrameqct::evalrmem()
-{
-	if (rmemBlks) {
-		for (int i = 0; i < 3; i++)
-		{
-			auto& expression = rmemBlks->at(i);
-			double value = expression->getValue();
-			rmem->at(i) = value;
-		}
-	}
-}
+    void EndFrameqct::evalrmem() {
+        if (rmemBlks) {
+            for (int i = 0; i < 3; i++) {
+                auto &expression = rmemBlks->at(i);
+                double value = expression->getValue();
+                rmem->at(i) = value;
+            }
+        }
+    }
 
-void EndFrameqct::evalAme()
-{
-	if (phiThePsiBlks) {
-		auto phiThePsi = CREATE<EulerAngleszxz<double>>::With();
-		for (int i = 0; i < 3; i++)
-		{
-			auto& expression = phiThePsiBlks->at(i);
-			auto value = expression->getValue();
-			phiThePsi->at(i) = value;
-		}
-		phiThePsi->calc();
-		aAme = phiThePsi->aA;
-	}
-}
+    void EndFrameqct::evalAme() {
+        if (phiThePsiBlks) {
+            auto phiThePsi = CREATE<EulerAngleszxz<double>>::With();
+            for (int i = 0; i < 3; i++) {
+                auto &expression = phiThePsiBlks->at(i);
+                auto value = expression->getValue();
+                phiThePsi->at(i) = value;
+            }
+            phiThePsi->calc();
+            aAme = phiThePsi->aA;
+        }
+    }
 
-void EndFrameqct::preVelIC()
-{
-	time = this->root()->mbdTimeValue();
-	this->evalrmem();
-	this->evalAme();
-	Item::preVelIC();
-	this->evalprmempt();
-	this->evalpAmept();
-	auto& aAOm = markerFrame->aAOm;
-	prOeOpt = aAOm->timesFullColumn(prmempt);
-	pAOept = aAOm->timesFullMatrix(pAmept);
-}
+    void EndFrameqct::preVelIC() {
+        time = this->root()->mbdTimeValue();
+        this->evalrmem();
+        this->evalAme();
+        Item::preVelIC();
+        this->evalprmempt();
+        this->evalpAmept();
+        auto &aAOm = markerFrame->aAOm;
+        prOeOpt = aAOm->timesFullColumn(prmempt);
+        pAOept = aAOm->timesFullMatrix(pAmept);
+    }
 
-void EndFrameqct::postVelIC()
-{
-	auto& pAOmpE = markerFrame->pAOmpE;
-	for (int i = 0; i < 3; i++)
-	{
-		auto& pprOeOpEpti = pprOeOpEpt->at(i);
-		for (int j = 0; j < 4; j++)
-		{
-			auto pprOeOpEptij = pAOmpE->at(j)->at(i)->dot(prmempt);
-			pprOeOpEpti->atiput(j, pprOeOpEptij);
-		}
-	}
-	for (int i = 0; i < 4; i++)
-	{
-		ppAOepEpt->atiput(i, pAOmpE->at(i)->timesFullMatrix(pAmept));
-	}
-}
+    void EndFrameqct::postVelIC() {
+        auto &pAOmpE = markerFrame->pAOmpE;
+        for (int i = 0; i < 3; i++) {
+            auto &pprOeOpEpti = pprOeOpEpt->at(i);
+            for (int j = 0; j < 4; j++) {
+                auto pprOeOpEptij = pAOmpE->at(j)->at(i)->dot(prmempt);
+                pprOeOpEpti->atiput(j, pprOeOpEptij);
+            }
+        }
+        for (int i = 0; i < 4; i++) {
+            ppAOepEpt->atiput(i, pAOmpE->at(i)->timesFullMatrix(pAmept));
+        }
+    }
 
-FColDsptr EndFrameqct::pAjOept(int j)
-{
-	return pAOept->column(j);
-}
+    FColDsptr EndFrameqct::pAjOept(int j) {
+        return pAOept->column(j);
+    }
 
-FMatDsptr EndFrameqct::ppAjOepETpt(int jj)
-{
-	auto answer = std::make_shared<FullMatrix<double>>(4, 3);
-	for (int i = 0; i < 4; i++)
-	{
-		auto& answeri = answer->at(i);
-		auto& ppAOepEipt = ppAOepEpt->at(i);
-		for (int j = 0; j < 3; j++)
-		{
-			auto& answerij = ppAOepEipt->at(j)->at(jj);
-			answeri->atiput(j, answerij);
-		}
-	}
-	return answer;
-}
+    FMatDsptr EndFrameqct::ppAjOepETpt(int jj) {
+        auto answer = std::make_shared<FullMatrixDouble>(4, 3);
+        for (int i = 0; i < 4; i++) {
+            auto &answeri = answer->at(i);
+            auto &ppAOepEipt = ppAOepEpt->at(i);
+            for (int j = 0; j < 3; j++) {
+                auto &answerij = ppAOepEipt->at(j)->at(jj);
+                answeri->atiput(j, answerij);
+            }
+        }
+        return answer;
+    }
 
-FColDsptr EndFrameqct::ppAjOeptpt(int j)
-{
-	return ppAOeptpt->column(j);
-}
+    FColDsptr EndFrameqct::ppAjOeptpt(int j) {
+        return ppAOeptpt->column(j);
+    }
 
-double EndFrameqct::priOeOpt(int i)
-{
-	return prOeOpt->at(i);
-}
+    double EndFrameqct::priOeOpt(int i) {
+        return prOeOpt->at(i);
+    }
 
-FRowDsptr EndFrameqct::ppriOeOpEpt(int i)
-{
-	return pprOeOpEpt->at(i);
-}
+    FRowDsptr EndFrameqct::ppriOeOpEpt(int i) {
+        return pprOeOpEpt->at(i);
+    }
 
-double EndFrameqct::ppriOeOptpt(int i)
-{
-	return pprOeOptpt->at(i);
-}
+    double EndFrameqct::ppriOeOptpt(int i) {
+        return pprOeOptpt->at(i);
+    }
 
-void EndFrameqct::evalprmempt()
-{
-	if (rmemBlks) {
-		for (int i = 0; i < 3; i++)
-		{
-			auto& derivative = prmemptBlks->at(i);
-			auto value = derivative->getValue();
-			prmempt->at(i) = value;
-		}
-	}
-}
+    void EndFrameqct::evalprmempt() {
+        if (rmemBlks) {
+            for (int i = 0; i < 3; i++) {
+                auto &derivative = prmemptBlks->at(i);
+                auto value = derivative->getValue();
+                prmempt->at(i) = value;
+            }
+        }
+    }
 
-void EndFrameqct::evalpAmept()
-{
-	if (phiThePsiBlks) {
-		auto phiThePsi = CREATE<EulerAngleszxz<double>>::With();
-		auto phiThePsiDot = CREATE<EulerAngleszxzDot<double>>::With();
-		phiThePsiDot->phiThePsi = phiThePsi;
-		for (int i = 0; i < 3; i++)
-		{
-			auto& expression = phiThePsiBlks->at(i);
-			auto& derivative = pPhiThePsiptBlks->at(i);
-			auto value = expression->getValue();
-			auto valueDot = derivative->getValue();
-			phiThePsi->at(i) = value;
-			phiThePsiDot->at(i) = valueDot;
-		}
-		phiThePsi->calc();
-		phiThePsiDot->calc();
-		pAmept = phiThePsiDot->aAdot;
-	}
-}
+    void EndFrameqct::evalpAmept() {
+        if (phiThePsiBlks) {
+            auto phiThePsi = CREATE<EulerAngleszxz<double>>::With();
+            auto phiThePsiDot = CREATE<EulerAngleszxzDot<double>>::With();
+            phiThePsiDot->phiThePsi = phiThePsi;
+            for (int i = 0; i < 3; i++) {
+                auto &expression = phiThePsiBlks->at(i);
+                auto &derivative = pPhiThePsiptBlks->at(i);
+                auto value = expression->getValue();
+                auto valueDot = derivative->getValue();
+                phiThePsi->at(i) = value;
+                phiThePsiDot->at(i) = valueDot;
+            }
+            phiThePsi->calc();
+            phiThePsiDot->calc();
+            pAmept = phiThePsiDot->aAdot;
+        }
+    }
 
-void EndFrameqct::evalpprmemptpt()
-{
-	if (rmemBlks) {
-		for (int i = 0; i < 3; i++)
-		{
-			auto& secondDerivative = pprmemptptBlks->at(i);
-			auto value = secondDerivative->getValue();
-			pprmemptpt->atiput(i, value);
-		}
-	}
-}
+    void EndFrameqct::evalpprmemptpt() {
+        if (rmemBlks) {
+            for (int i = 0; i < 3; i++) {
+                auto &secondDerivative = pprmemptptBlks->at(i);
+                auto value = secondDerivative->getValue();
+                pprmemptpt->atiput(i, value);
+            }
+        }
+    }
 
-void EndFrameqct::evalppAmeptpt()
-{
-	if (phiThePsiBlks) {
-		auto phiThePsi = CREATE<EulerAngleszxz<double>>::With();
-		auto phiThePsiDot = CREATE<EulerAngleszxzDot<double>>::With();
-		phiThePsiDot->phiThePsi = phiThePsi;
-		auto phiThePsiDDot = CREATE<EulerAngleszxzDDot<double>>::With();
-		phiThePsiDDot->phiThePsiDot = phiThePsiDot;
-		for (int i = 0; i < 3; i++)
-		{
-			auto& expression = phiThePsiBlks->at(i);
-			auto& derivative = pPhiThePsiptBlks->at(i);
-			auto& secondDerivative = ppPhiThePsiptptBlks->at(i);
-			auto value = expression->getValue();
-			auto valueDot = derivative->getValue();
-			auto valueDDot = secondDerivative->getValue();
-			phiThePsi->atiput(i, value);
-			phiThePsiDot->atiput(i, valueDot);
-			phiThePsiDDot->atiput(i, valueDDot);
-		}
-		phiThePsi->calc();
-		phiThePsiDot->calc();
-		phiThePsiDDot->calc();
-		ppAmeptpt = phiThePsiDDot->aAddot;
-	}
-}
+    void EndFrameqct::evalppAmeptpt() {
+        if (phiThePsiBlks) {
+            auto phiThePsi = CREATE<EulerAngleszxz<double>>::With();
+            auto phiThePsiDot = CREATE<EulerAngleszxzDot<double>>::With();
+            phiThePsiDot->phiThePsi = phiThePsi;
+            auto phiThePsiDDot = CREATE<EulerAngleszxzDDot<double>>::With();
+            phiThePsiDDot->phiThePsiDot = phiThePsiDot;
+            for (int i = 0; i < 3; i++) {
+                auto &expression = phiThePsiBlks->at(i);
+                auto &derivative = pPhiThePsiptBlks->at(i);
+                auto &secondDerivative = ppPhiThePsiptptBlks->at(i);
+                auto value = expression->getValue();
+                auto valueDot = derivative->getValue();
+                auto valueDDot = secondDerivative->getValue();
+                phiThePsi->atiput(i, value);
+                phiThePsiDot->atiput(i, valueDot);
+                phiThePsiDDot->atiput(i, valueDDot);
+            }
+            phiThePsi->calc();
+            phiThePsiDot->calc();
+            phiThePsiDDot->calc();
+            ppAmeptpt = phiThePsiDDot->aAddot;
+        }
+    }
 
-FColDsptr EndFrameqct::rmeO()
-{
-	return markerFrame->aAOm->timesFullColumn(rmem);
-}
+    FColDsptr EndFrameqct::rmeO() {
+        return markerFrame->aAOm->timesFullColumn(rmem);
+    }
 
-FColDsptr EndFrameqct::rpep()
-{
-	auto& rpmp = markerFrame->rpmp;
-	auto& aApm = markerFrame->aApm;
-	auto rpep = rpmp->plusFullColumn(aApm->timesFullColumn(rmem));
-	return rpep;
-}
+    FColDsptr EndFrameqct::rpep() {
+        auto &rpmp = markerFrame->rpmp;
+        auto &aApm = markerFrame->aApm;
+        auto rpep = rpmp->plusFullColumn(aApm->timesFullColumn(rmem));
+        return rpep;
+    }
 
-void EndFrameqct::preAccIC()
-{
-	time = this->root()->mbdTimeValue();
-	this->evalrmem();
-	this->evalAme();
-	Item::preVelIC();
-	this->evalprmempt();
-	this->evalpAmept();
-	auto& aAOm = markerFrame->aAOm;
-	prOeOpt = aAOm->timesFullColumn(prmempt);
-	pAOept = aAOm->timesFullMatrix(pAmept);
-	Item::preAccIC();
-	this->evalpprmemptpt();
-	this->evalppAmeptpt();
-	aAOm = markerFrame->aAOm;
-	pprOeOptpt = aAOm->timesFullColumn(pprmemptpt);
-	ppAOeptpt = aAOm->timesFullMatrix(ppAmeptpt);
+    void EndFrameqct::preAccIC() {
+        time = this->root()->mbdTimeValue();
+        this->evalrmem();
+        this->evalAme();
+        Item::preVelIC();
+        this->evalprmempt();
+        this->evalpAmept();
+        auto &aAOm = markerFrame->aAOm;
+        prOeOpt = aAOm->timesFullColumn(prmempt);
+        pAOept = aAOm->timesFullMatrix(pAmept);
+        Item::preAccIC();
+        this->evalpprmemptpt();
+        this->evalppAmeptpt();
+        aAOm = markerFrame->aAOm;
+        pprOeOptpt = aAOm->timesFullColumn(pprmemptpt);
+        ppAOeptpt = aAOm->timesFullMatrix(ppAmeptpt);
+    }
 }

OndselSolver/EndFrameqct.h

@@ -58,7 +58,6 @@ namespace MbD {
 		FMatDsptr aAme, pAmept, ppAmeptpt, pAOept, ppAOeptpt;
 		FMatDsptr pprOeOpEpt;
 		FColFMatDsptr ppAOepEpt;
-		
 	};
 }
 

OndselSolver/EndFrameqct2.cpp

@@ -29,7 +29,7 @@ void EndFrameqct2::initpPhiThePsiptBlks()
 {
 	auto& mbdTime = this->root()->time;
 	auto eulerAngles = std::static_pointer_cast<EulerAngles<Symsptr>>(phiThePsiBlks);
-	pPhiThePsiptBlks = eulerAngles->differentiateWRT(mbdTime);
+	pPhiThePsiptBlks = differentiateWRT(*eulerAngles, mbdTime);
 }
 
 void EndFrameqct2::initppPhiThePsiptptBlks()

OndselSolver/EulerAngles.cpp

@@ -7,3 +7,73 @@
  ***************************************************************************/
  
 #include "EulerAngles.h"
+
+namespace MbD {
+    template<typename T>
+    inline void EulerAngles<T>::initialize()
+    {
+        assert(false);
+    }
+    template<>
+    inline void EulerAngles<Symsptr>::calc()
+    {
+        cA = std::make_shared<FullColumn<FMatDsptr>>(3);
+        for (int i = 0; i < 3; i++)
+        {
+            auto axis = rotOrder->at(i);
+            auto angle = this->at(i)->getValue();
+            if (axis == 1) {
+                cA->atiput(i, FullMatrixDouble::rotatex(angle));
+            }
+            else if (axis == 2) {
+                cA->atiput(i, FullMatrixDouble::rotatey(angle));
+            }
+            else if (axis == 3) {
+                cA->atiput(i, FullMatrixDouble::rotatez(angle));
+            }
+            else {
+                throw std::runtime_error("Euler angle rotation order must be any permutation of 1,2,3 without consecutive repeats.");
+            }
+        }
+        aA = cA->at(0)->timesFullMatrix(cA->at(1)->timesFullMatrix(cA->at(2)));
+    }
+    template<>
+    inline void EulerAngles<double>::calc()
+    {
+        cA = std::make_shared<FullColumn<FMatDsptr>>(3);
+        for (int i = 0; i < 3; i++)
+        {
+            auto axis = rotOrder->at(i);
+            auto angle = this->at(i);
+            if (axis == 1) {
+                cA->atiput(i, FullMatrixDouble::rotatex(angle));
+            }
+            else if (axis == 2) {
+                cA->atiput(i, FullMatrixDouble::rotatey(angle));
+            }
+            else if (axis == 3) {
+                cA->atiput(i, FullMatrixDouble::rotatez(angle));
+            }
+            else {
+                throw std::runtime_error("Euler angle rotation order must be any permutation of 1,2,3 without consecutive repeats.");
+            }
+        }
+        aA = cA->at(0)->timesFullMatrix(cA->at(1)->timesFullMatrix(cA->at(2)));
+    }
+    template<typename T>
+    inline void EulerAngles<T>::calc()
+    {
+        assert(false);
+    }
+    // type-specific helper functions
+    std::shared_ptr<EulerAnglesDot<std::shared_ptr<MbD::Symbolic>>> differentiateWRT(EulerAngles<std::shared_ptr<MbD::Symbolic>>& ref, std::shared_ptr<MbD::Symbolic> var)
+    {
+        auto derivatives = std::make_shared<EulerAnglesDot<std::shared_ptr<MbD::Symbolic>>>();
+        std::transform(ref.begin(), ref.end(), derivatives->begin(),
+                       [var](std::shared_ptr<MbD::Symbolic> term) { return term->differentiateWRT(var); }
+        );
+        derivatives->aEulerAngles = &ref;
+        return derivatives;
+    }
+}
+

OndselSolver/EulerAngles.h

@@ -16,8 +16,6 @@
 #include "Symbolic.h"
 
 namespace MbD {
-	//template<typename T>
-	//class EulerAnglesDot;
 
 	template<typename T>
 	class EulerAngles : public EulerArray<T>
@@ -29,87 +27,23 @@ namespace MbD {
 		EulerAngles(std::initializer_list<T> list) : EulerArray<T>{ list } {}
 		void initialize() override;
 		void calc() override;
-		std::shared_ptr<EulerAnglesDot<T>> differentiateWRT(T var);
 		void setRotOrder(int i, int j, int k);
 
 		std::shared_ptr<FullColumn<int>> rotOrder;
 		FColFMatDsptr cA;
 		FMatDsptr aA;
-
 	};
-	template<typename T>
-	inline void EulerAngles<T>::initialize()
-	{
-		assert(false);
-	}
-	template<>
-	inline void EulerAngles<Symsptr>::calc()
-	{
-		cA = std::make_shared<FullColumn<FMatDsptr>>(3);
-		for (int i = 0; i < 3; i++)
-		{
-			auto axis = rotOrder->at(i);
-			auto angle = this->at(i)->getValue();
-			if (axis == 1) {
-				cA->atiput(i, FullMatrix<double>::rotatex(angle));
-			}
-			else if (axis == 2) {
-				cA->atiput(i, FullMatrix<double>::rotatey(angle));
-			}
-			else if (axis == 3) {
-				cA->atiput(i, FullMatrix<double>::rotatez(angle));
-			}
-			else {
-				throw std::runtime_error("Euler angle rotation order must be any permutation of 1,2,3 without consecutive repeats.");
-			}
-		}
-		aA = cA->at(0)->timesFullMatrix(cA->at(1)->timesFullMatrix(cA->at(2)));
-	}
-	template<>
-	inline void EulerAngles<double>::calc()
-	{
-		cA = std::make_shared<FullColumn<FMatDsptr>>(3);
-		for (int i = 0; i < 3; i++)
-		{
-			auto axis = rotOrder->at(i);
-			auto angle = this->at(i);
-			if (axis == 1) {
-				cA->atiput(i, FullMatrix<double>::rotatex(angle));
-			}
-			else if (axis == 2) {
-				cA->atiput(i, FullMatrix<double>::rotatey(angle));
-			}
-			else if (axis == 3) {
-				cA->atiput(i, FullMatrix<double>::rotatez(angle));
-			}
-			else {
-				throw std::runtime_error("Euler angle rotation order must be any permutation of 1,2,3 without consecutive repeats.");
-			}
-		}
-		aA = cA->at(0)->timesFullMatrix(cA->at(1)->timesFullMatrix(cA->at(2)));
-	}
-	template<typename T>
-	inline void EulerAngles<T>::calc()
-	{
-		assert(false);
-	}
-	template<typename T>
-	inline std::shared_ptr<EulerAnglesDot<T>> EulerAngles<T>::differentiateWRT(T var)
-	{
-		auto derivatives = std::make_shared<EulerAnglesDot<T>>();
-		std::transform(this->begin(), this->end(), derivatives->begin(),
-			[var](T term) { return term->differentiateWRT(var); }
-		);
-		derivatives->aEulerAngles = this;
-		return derivatives;
-	}
-	template<typename T>
-	inline void EulerAngles<T>::setRotOrder(int i, int j, int k)
-	{
-		rotOrder = std::make_shared<FullColumn<int>>(3);
-		rotOrder->at(0) = i;
-		rotOrder->at(1) = j;
-		rotOrder->at(2) = k;
-	}
+    template class EulerAngles<std::shared_ptr<MbD::Symbolic>>;
+    template class EulerAngles<double>;
+    std::shared_ptr<EulerAnglesDot<std::shared_ptr<MbD::Symbolic>>> differentiateWRT(EulerAngles<std::shared_ptr<MbD::Symbolic>>& ref, std::shared_ptr<MbD::Symbolic> var);
+    template<typename T>
+
+    void EulerAngles<T>::setRotOrder(int i, int j, int k)
+    {
+        rotOrder = std::make_shared<FullColumn<int>>(3);
+        rotOrder->at(0) = i;
+        rotOrder->at(1) = j;
+        rotOrder->at(2) = k;
+    }
 }
 

OndselSolver/EulerAnglesDDot.h

@@ -47,13 +47,13 @@ namespace MbD {
 			auto angleDot = aEulerAnglesDot->at(i)->getValue();
 			auto angleDDot = this->at(i)->getValue();
 			if (axis == 1) {
-				cAddot->atiput(i, FullMatrix<double>::rotatexrotDotrotDDot(angle, angleDot, angleDDot));
+				cAddot->atiput(i, FullMatrixDouble::rotatexrotDotrotDDot(angle, angleDot, angleDDot));
 			}
 			else if (axis == 2) {
-				cAddot->atiput(i, FullMatrix<double>::rotateyrotDotrotDDot(angle, angleDot, angleDDot));
+				cAddot->atiput(i, FullMatrixDouble::rotateyrotDotrotDDot(angle, angleDot, angleDDot));
 			}
 			else if (axis == 3) {
-				cAddot->atiput(i, FullMatrix<double>::rotatezrotDotrotDDot(angle, angleDot, angleDDot));
+				cAddot->atiput(i, FullMatrixDouble::rotatezrotDotrotDDot(angle, angleDot, angleDDot));
 			}
 			else {
 				throw std::runtime_error("Euler angle rotation order must be any permutation of 1,2,3 without consecutive repeats.");
@@ -79,8 +79,10 @@ namespace MbD {
 		auto term7 = phiA->timesFullMatrix(theAdot->timesFullMatrix(psiAdot));
 		auto term8 = phiA->timesFullMatrix(theA->timesFullMatrix(psiAddot));
 
-		aAddot = term->plusFullMatrix(term1)->plusFullMatrix(term2)->plusFullMatrix(term3)->plusFullMatrix(term4)
-			->plusFullMatrix(term5)->plusFullMatrix(term6)->plusFullMatrix(term7)->plusFullMatrix(term8);
+		aAddot = term->plusFullMatrix(term1)->plusFullMatrix(term2)
+                ->plusFullMatrix(term3)->plusFullMatrix(term4)
+			    ->plusFullMatrix(term5)->plusFullMatrix(term6)
+                ->plusFullMatrix(term7)->plusFullMatrix(term8);
 	}
 	template<typename T>
 	inline void EulerAnglesDDot<T>::aEulerAngles(EulerAngles<T>* eulerAngles)

OndselSolver/EulerAnglesDot.h

@@ -56,13 +56,13 @@ namespace MbD {
 			auto angle = aEulerAngles->at(i)->getValue();
 			auto angleDot = this->at(i)->getValue();
 			if (axis == 1) {
-				cAdot->atiput(i, FullMatrix<double>::rotatexrotDot(angle, angleDot));
+				cAdot->atiput(i, FullMatrixDouble::rotatexrotDot(angle, angleDot));
 			}
 			else if (axis == 2) {
-				cAdot->atiput(i, FullMatrix<double>::rotateyrotDot(angle, angleDot));
+				cAdot->atiput(i, FullMatrixDouble::rotateyrotDot(angle, angleDot));
 			}
 			else if (axis == 3) {
-				cAdot->atiput(i, FullMatrix<double>::rotatezrotDot(angle, angleDot));
+				cAdot->atiput(i, FullMatrixDouble::rotatezrotDot(angle, angleDot));
 			}
 			else {
 				throw std::runtime_error("Euler angle rotation order must be any permutation of 1,2,3 without consecutive repeats.");
@@ -79,11 +79,11 @@ namespace MbD {
 		auto psiAdot = cAdot->at(2);
 
 		aAdot = phiAdot->timesFullMatrix(theA->timesFullMatrix(psiA))
-			->plusFullMatrix(phiA->timesFullMatrix(theAdot->timesFullMatrix(psiA)))
-			->plusFullMatrix(phiA->timesFullMatrix(theA->timesFullMatrix(psiAdot)));
+                ->plusFullMatrix(phiA->timesFullMatrix(theAdot->timesFullMatrix(psiA)))
+                ->plusFullMatrix(phiA->timesFullMatrix(theA->timesFullMatrix(psiAdot)));
 		omeF = (phiA->column(0)->times(phidot)
-			->plusFullColumn(phiA->timesFullMatrix(theA)->column(1)->times(thedot))
-			->plusFullColumn(aEulerAngles->aA->column(2)->times(psidot)));
+                ->plusFullColumn(phiA->timesFullMatrix(theA)->column(1)->times(thedot))
+                ->plusFullColumn(aEulerAngles->aA->column(2)->times(psidot)));
 		omef = aEulerAngles->aA->transposeTimesFullColumn(omeF);
 	}
 }

OndselSolver/EulerAngleszxz.h

@@ -29,9 +29,9 @@ namespace MbD {
 	template<typename T>
 	inline void EulerAngleszxz<T>::initialize()
 	{
-		phiA = FullMatrix<T>::identitysptr(3);
-		theA = FullMatrix<T>::identitysptr(3);
-		psiA = FullMatrix<T>::identitysptr(3);
+		phiA = FullMatrixDouble::identitysptr(3);
+		theA = FullMatrixDouble::identitysptr(3);
+		psiA = FullMatrixDouble::identitysptr(3);
 	}
 	template<typename T>
 	inline void EulerAngleszxz<T>::calc()

OndselSolver/EulerAngleszxzDDot.h

@@ -8,6 +8,7 @@
  
 #pragma once
 
+#include "FullMatrix.h"
 #include "EulerArray.h"
 #include "EulerAngleszxzDot.h"
 
@@ -23,16 +24,17 @@ namespace MbD {
         void calc() override;
 
         std::shared_ptr<EulerAngleszxzDot<double>> phiThePsiDot;
-        FMatDsptr phiAddot, theAddot, psiAddot, aAddot;
+        FMatDsptr phiAddot, theAddot, psiAddot;
+        std::shared_ptr<FullMatrixDouble> aAddot;
     };
     template<typename T>
     inline void EulerAngleszxzDDot<T>::initialize()
     {
-        phiAddot = std::make_shared<FullMatrix<double>>(3, 3);
+        phiAddot = std::make_shared<FullMatrixDouble>(3, 3);
         phiAddot->zeroSelf();
-        theAddot = std::make_shared<FullMatrix<double>>(3, 3);
+        theAddot = std::make_shared<FullMatrixDouble>(3, 3);
         theAddot->zeroSelf();
-        psiAddot = std::make_shared<FullMatrix<double>>(3, 3);
+        psiAddot = std::make_shared<FullMatrixDouble>(3, 3);
         psiAddot->zeroSelf();
     }
     template<typename T>
@@ -90,7 +92,7 @@ namespace MbD {
 			+ *(phiAdot->timesFullMatrix(theA->timesFullMatrix(psiAdot)))
 			+ *(phiA->timesFullMatrix(theAdot->timesFullMatrix(psiAdot)))
 			+ *(phiA->timesFullMatrix(theA->timesFullMatrix(psiAddot)));
-		aAddot = std::make_shared<FullMatrix<double>>(mat);
+		aAddot = std::make_shared<FullMatrixDouble>(mat);
 	}
 }
 

OndselSolver/EulerAngleszxzDot.h

@@ -28,11 +28,11 @@ namespace MbD {
 	template<typename T>
 	inline void EulerAngleszxzDot<T>::initialize()
 	{
-		phiAdot = std::make_shared<FullMatrix<double>>(3, 3);
+		phiAdot = std::make_shared<FullMatrixDouble>(3, 3);
 		phiAdot->zeroSelf();
-		theAdot = std::make_shared<FullMatrix<double>>(3, 3);
+		theAdot = std::make_shared<FullMatrixDouble>(3, 3);
 		theAdot->zeroSelf();
-		psiAdot = std::make_shared<FullMatrix<double>>(3, 3);
+		psiAdot = std::make_shared<FullMatrixDouble>(3, 3);
 		psiAdot->zeroSelf();
 	}
 	template<typename T>

OndselSolver/EulerParameters.cpp

@@ -7,5 +7,308 @@
  ***************************************************************************/
  
 #include "EulerParameters.h"
+#include "FullColumn.h"
+#include "FullRow.h"
+#include "FullMatrix.h"
 
-using namespace MbD;
+namespace MbD {
+    template<>
+    FMatFColDsptr EulerParameters<double>::ppApEpEtimesColumn(FColDsptr col)
+    {
+        double a2c0 = 2 * col->at(0);
+        double a2c1 = 2 * col->at(1);
+        double a2c2 = 2 * col->at(2);
+        double m2c0 = 0 - a2c0;
+        double m2c1 = 0 - a2c1;
+        double m2c2 = 0 - a2c2;
+        auto col00 = std::make_shared<FullColumn<double>>(ListD{ a2c0, m2c1, m2c2 });
+        auto col01 = std::make_shared<FullColumn<double>>(ListD{ a2c1, a2c0, 0 });
+        auto col02 = std::make_shared<FullColumn<double>>(ListD{ a2c2, 0, a2c0 });
+        auto col03 = std::make_shared<FullColumn<double>>(ListD{ 0, m2c2, a2c1 });
+        auto col11 = std::make_shared<FullColumn<double>>(ListD{ m2c0, a2c1, m2c2 });
+        auto col12 = std::make_shared<FullColumn<double>>(ListD{ 0, a2c2, a2c1 });
+        auto col13 = std::make_shared<FullColumn<double>>(ListD{ a2c2, 0, m2c0 });
+        auto col22 = std::make_shared<FullColumn<double>>(ListD{ m2c0, m2c1, a2c2 });
+        auto col23 = std::make_shared<FullColumn<double>>(ListD{ m2c1, a2c0, 0 });
+        auto col33 = std::make_shared<FullColumn<double>>(ListD{ a2c0, a2c1, a2c2 });
+        auto answer = std::make_shared<FullMatrixFullColumnDouble>(4, 4);
+        auto& row0 = answer->at(0);
+        row0->at(0) = col00;
+        row0->at(1) = col01;
+        row0->at(2) = col02;
+        row0->at(3) = col03;
+        auto& row1 = answer->at(1);
+        row1->at(0) = col01;
+        row1->at(1) = col11;
+        row1->at(2) = col12;
+        row1->at(3) = col13;
+        auto& row2 = answer->at(2);
+        row2->at(0) = col02;
+        row2->at(1) = col12;
+        row2->at(2) = col22;
+        row2->at(3) = col23;
+        auto& row3 = answer->at(3);
+        row3->at(0) = col03;
+        row3->at(1) = col13;
+        row3->at(2) = col23;
+        row3->at(3) = col33;
+        return answer;
+    }
+
+    template<>
+    FMatDsptr EulerParameters<double>::pCpEtimesColumn(FColDsptr col)
+    {
+        //"col size = 4."
+        auto c0 = col->at(0);
+        auto c1 = col->at(1);
+        auto c2 = col->at(2);
+        auto mc0 = -c0;
+        auto mc1 = -c1;
+        auto mc2 = -c2;
+        auto mc3 = -col->at(3);
+        auto answer = std::make_shared<FullMatrixDouble>(3, 4);
+        auto& row0 = answer->at(0);
+        auto& row1 = answer->at(1);
+        auto& row2 = answer->at(2);
+        row0->atiput(0, mc3);
+        row0->atiput(1, mc2);
+        row0->atiput(2, c1);
+        row0->atiput(3, c0);
+        row1->atiput(0, c2);
+        row1->atiput(1, mc3);
+        row1->atiput(2, mc0);
+        row1->atiput(3, c1);
+        row2->atiput(0, mc1);
+        row2->atiput(1, c0);
+        row2->atiput(2, mc3);
+        row2->atiput(3, c2);
+        return answer;
+    }
+
+    template<typename T>
+    FMatDsptr EulerParameters<T>::pCTpEtimesColumn(FColDsptr col)
+    {
+        //"col size = 3."
+        auto c0 = col->at(0);
+        auto c1 = col->at(1);
+        auto c2 = col->at(2);
+        auto mc0 = -c0;
+        auto mc1 = -c1;
+        auto mc2 = -c2;
+        auto answer = std::make_shared<FullMatrixDouble>(4, 4);
+        auto& row0 = answer->at(0);
+        auto& row1 = answer->at(1);
+        auto& row2 = answer->at(2);
+        auto& row3 = answer->at(3);
+        row0->atiput(0, 0.0);
+        row0->atiput(1, c2);
+        row0->atiput(2, mc1);
+        row0->atiput(3, c0);
+        row1->atiput(0, mc2);
+        row1->atiput(1, 0.0);
+        row1->atiput(2, c0);
+        row1->atiput(3, c1);
+        row2->atiput(0, c1);
+        row2->atiput(1, mc0);
+        row2->atiput(2, 0.0);
+        row2->atiput(3, c2);
+        row3->atiput(0, mc0);
+        row3->atiput(1, mc1);
+        row3->atiput(2, mc2);
+        row3->atiput(3, 0.0);
+        return answer;
+    }
+
+    template<>
+    FMatFMatDsptr EulerParameters<double>::ppApEpEtimesMatrix(FMatDsptr mat)
+    {
+        FRowDsptr a2m0 = mat->at(0)->times(2.0);
+        FRowDsptr a2m1 = mat->at(1)->times(2.0);
+        FRowDsptr a2m2 = mat->at(2)->times(2.0);
+        FRowDsptr m2m0 = a2m0->negated();
+        FRowDsptr m2m1 = a2m1->negated();
+        FRowDsptr m2m2 = a2m2->negated();
+        FRowDsptr zero = std::make_shared<FullRow<double>>(3, 0.0);
+        auto mat00 = std::make_shared<FullMatrixDouble>(ListFRD{ a2m0, m2m1, m2m2 });
+        auto mat01 = std::make_shared<FullMatrixDouble>(ListFRD{ a2m1, a2m0, zero });
+        auto mat02 = std::make_shared<FullMatrixDouble>(ListFRD{ a2m2, zero, a2m0 });
+        auto mat03 = std::make_shared<FullMatrixDouble>(ListFRD{ zero, m2m2, a2m1 });
+        auto mat11 = std::make_shared<FullMatrixDouble>(ListFRD{ m2m0, a2m1, m2m2 });
+        auto mat12 = std::make_shared<FullMatrixDouble>(ListFRD{ zero, a2m2, a2m1 });
+        auto mat13 = std::make_shared<FullMatrixDouble>(ListFRD{ a2m2, zero, m2m0 });
+        auto mat22 = std::make_shared<FullMatrixDouble>(ListFRD{ m2m0, m2m1, a2m2 });
+        auto mat23 = std::make_shared<FullMatrixDouble>(ListFRD{ m2m1, a2m0, zero });
+        auto mat33 = std::make_shared<FullMatrixDouble>(ListFRD{ a2m0, a2m1, a2m2 });
+        auto answer = std::make_shared<FullMatrixFullMatrixDouble>(4, 4);
+        auto& row0 = answer->at(0);
+        row0->at(0) = mat00;
+        row0->at(1) = mat01;
+        row0->at(2) = mat02;
+        row0->at(3) = mat03;
+        auto& row1 = answer->at(1);
+        row1->at(0) = mat01;
+        row1->at(1) = mat11;
+        row1->at(2) = mat12;
+        row1->at(3) = mat13;
+        auto& row2 = answer->at(2);
+        row2->at(0) = mat02;
+        row2->at(1) = mat12;
+        row2->at(2) = mat22;
+        row2->at(3) = mat23;
+        auto& row3 = answer->at(3);
+        row3->at(0) = mat03;
+        row3->at(1) = mat13;
+        row3->at(2) = mat23;
+        row3->at(3) = mat33;
+        return answer;
+    }
+
+    template<typename T> // this is ALWAYS double; see note below.
+    void EulerParameters<T>::initialize()
+    {
+        aA = FullMatrixDouble::identitysptr(3);
+        aB = std::make_shared<FullMatrixDouble>(3, 4);
+        aC = std::make_shared<FullMatrixDouble>(3, 4);
+        pApE = std::make_shared<FullColumn<FMatDsptr>>(4);
+        for (int i = 0; i < 4; i++)
+        {
+            pApE->at(i) = std::make_shared<FullMatrixDouble>(3, 3);
+        }
+    }
+
+// the following  can't be valid as FullMatrix instatiatiates <double>, yet
+// this class needs to see the member functions of FullMatrix
+//template<>
+//void EulerParameters<double>::initialize()
+//{
+//}
+
+    template<typename T>
+    void EulerParameters<T>::calc()
+    {
+        this->calcABC();
+        this->calcpApE();
+    }
+    template<>
+    void EulerParameters<double>::calcABC()
+    {
+        double aE0 = this->at(0);
+        double aE1 = this->at(1);
+        double aE2 = this->at(2);
+        double aE3 = this->at(3);
+        double mE0 = -aE0;
+        double mE1 = -aE1;
+        double mE2 = -aE2;
+        FRowDsptr aBi;
+        aBi = aB->at(0);
+        aBi->at(0) = aE3;
+        aBi->at(1) = mE2;
+        aBi->at(2) = aE1;
+        aBi->at(3) = mE0;
+        aBi = aB->at(1);
+        aBi->at(0) = aE2;
+        aBi->at(1) = aE3;
+        aBi->at(2) = mE0;
+        aBi->at(3) = mE1;
+        aBi = aB->at(2);
+        aBi->at(0) = mE1;
+        aBi->at(1) = aE0;
+        aBi->at(2) = aE3;
+        aBi->at(3) = mE2;
+        FRowDsptr aCi;
+        aCi = aC->at(0);
+        aCi->at(0) = aE3;
+        aCi->at(1) = aE2;
+        aCi->at(2) = mE1;
+        aCi->at(3) = mE0;
+        aCi = aC->at(1);
+        aCi->at(0) = mE2;
+        aCi->at(1) = aE3;
+        aCi->at(2) = aE0;
+        aCi->at(3) = mE1;
+        aCi = aC->at(2);
+        aCi->at(0) = aE1;
+        aCi->at(1) = mE0;
+        aCi->at(2) = aE3;
+        aCi->at(3) = mE2;
+
+        aA = aB->timesTransposeFullMatrix(aC);
+    }
+    template<>
+    void EulerParameters<double>::calcpApE()
+    {
+        double a2E0 = 2.0 * (this->at(0));
+        double a2E1 = 2.0 * (this->at(1));
+        double a2E2 = 2.0 * (this->at(2));
+        double a2E3 = 2.0 * (this->at(3));
+        double m2E0 = -a2E0;
+        double m2E1 = -a2E1;
+        double m2E2 = -a2E2;
+        double m2E3 = -a2E3;
+        FMatDsptr pApEk;
+        pApEk = pApE->at(0);
+        FRowDsptr pAipEk;
+        pAipEk = pApEk->at(0);
+        pAipEk->at(0) = a2E0;
+        pAipEk->at(1) = a2E1;
+        pAipEk->at(2) = a2E2;
+        pAipEk = pApEk->at(1);
+        pAipEk->at(0) = a2E1;
+        pAipEk->at(1) = m2E0;
+        pAipEk->at(2) = m2E3;
+        pAipEk = pApEk->at(2);
+        pAipEk->at(0) = a2E2;
+        pAipEk->at(1) = a2E3;
+        pAipEk->at(2) = m2E0;
+        //
+        pApEk = pApE->at(1);
+        pAipEk = pApEk->at(0);
+        pAipEk->at(0) = m2E1;
+        pAipEk->at(1) = a2E0;
+        pAipEk->at(2) = a2E3;
+        pAipEk = pApEk->at(1);
+        pAipEk->at(0) = a2E0;
+        pAipEk->at(1) = a2E1;
+        pAipEk->at(2) = a2E2;
+        pAipEk = pApEk->at(2);
+        pAipEk->at(0) = m2E3;
+        pAipEk->at(1) = a2E2;
+        pAipEk->at(2) = m2E1;
+        //
+        pApEk = pApE->at(2);
+        pAipEk = pApEk->at(0);
+        pAipEk->at(0) = m2E2;
+        pAipEk->at(1) = m2E3;
+        pAipEk->at(2) = a2E0;
+        pAipEk = pApEk->at(1);
+        pAipEk->at(0) = a2E3;
+        pAipEk->at(1) = m2E2;
+        pAipEk->at(2) = a2E1;
+        pAipEk = pApEk->at(2);
+        pAipEk->at(0) = a2E0;
+        pAipEk->at(1) = a2E1;
+        pAipEk->at(2) = a2E2;
+        //
+        pApEk = pApE->at(3);
+        pAipEk = pApEk->at(0);
+        pAipEk->at(0) = a2E3;
+        pAipEk->at(1) = m2E2;
+        pAipEk->at(2) = a2E1;
+        pAipEk = pApEk->at(1);
+        pAipEk->at(0) = a2E2;
+        pAipEk->at(1) = a2E3;
+        pAipEk->at(2) = m2E0;
+        pAipEk = pApEk->at(2);
+        pAipEk->at(0) = m2E1;
+        pAipEk->at(1) = a2E0;
+        pAipEk->at(2) = a2E3;
+    }
+    template<typename T>
+    void EulerParameters<T>::conditionSelf()
+    {
+        EulerArray<T>::conditionSelf();
+        this->normalizeSelf();
+    }
+
+    template class EulerParameters<double>;
+}

OndselSolver/EulerParameters.h

@@ -8,9 +8,12 @@
  
 #pragma once
 
+#include "FullColumn.ref.h"
+#include "FullMatrix.ref.h"
+#include "EulerParameters.ref.h"
 #include "EulerArray.h"
-#include "FullColumn.h"
-#include "FullMatrix.h"
+// #include "FullColumn.h"
+// #include "FullMatrix.h"
 
 namespace MbD {
 
@@ -39,11 +42,10 @@ namespace MbD {
 			this->initialize();
 			this->calc();
 		}
-
-		static std::shared_ptr<FullMatrix<FColsptr<T>>> ppApEpEtimesColumn(FColDsptr col);
+		static std::shared_ptr<FullMatrixFullColumnDouble> ppApEpEtimesColumn(FColDsptr col);
 		static FMatDsptr pCpEtimesColumn(FColDsptr col);
 		static FMatDsptr pCTpEtimesColumn(FColDsptr col);
-		static std::shared_ptr<FullMatrix<FMatsptr<T>>> ppApEpEtimesMatrix(FMatDsptr mat);
+		static std::shared_ptr<FullMatrixFullMatrixDouble> ppApEpEtimesMatrix(FMatDsptr mat);
 
 
 		void initialize() override;
@@ -57,295 +59,4 @@ namespace MbD {
 		FMatDsptr aC;
 		FColFMatDsptr pApE;
 	};
-
-	template<>
-	inline FMatFColDsptr EulerParameters<double>::ppApEpEtimesColumn(FColDsptr col)
-	{
-		double a2c0 = 2 * col->at(0);
-		double a2c1 = 2 * col->at(1);
-		double a2c2 = 2 * col->at(2);
-		double m2c0 = 0 - a2c0;
-		double m2c1 = 0 - a2c1;
-		double m2c2 = 0 - a2c2;
-		auto col00 = std::make_shared<FullColumn<double>>(ListD{ a2c0, m2c1, m2c2 });
-		auto col01 = std::make_shared<FullColumn<double>>(ListD{ a2c1, a2c0, 0 });
-		auto col02 = std::make_shared<FullColumn<double>>(ListD{ a2c2, 0, a2c0 });
-		auto col03 = std::make_shared<FullColumn<double>>(ListD{ 0, m2c2, a2c1 });
-		auto col11 = std::make_shared<FullColumn<double>>(ListD{ m2c0, a2c1, m2c2 });
-		auto col12 = std::make_shared<FullColumn<double>>(ListD{ 0, a2c2, a2c1 });
-		auto col13 = std::make_shared<FullColumn<double>>(ListD{ a2c2, 0, m2c0 });
-		auto col22 = std::make_shared<FullColumn<double>>(ListD{ m2c0, m2c1, a2c2 });
-		auto col23 = std::make_shared<FullColumn<double>>(ListD{ m2c1, a2c0, 0 });
-		auto col33 = std::make_shared<FullColumn<double>>(ListD{ a2c0, a2c1, a2c2 });
-		auto answer = std::make_shared<FullMatrix<FColDsptr>>(4, 4);
-		auto& row0 = answer->at(0);
-		row0->at(0) = col00;
-		row0->at(1) = col01;
-		row0->at(2) = col02;
-		row0->at(3) = col03;
-		auto& row1 = answer->at(1);
-		row1->at(0) = col01;
-		row1->at(1) = col11;
-		row1->at(2) = col12;
-		row1->at(3) = col13;
-		auto& row2 = answer->at(2);
-		row2->at(0) = col02;
-		row2->at(1) = col12;
-		row2->at(2) = col22;
-		row2->at(3) = col23;
-		auto& row3 = answer->at(3);
-		row3->at(0) = col03;
-		row3->at(1) = col13;
-		row3->at(2) = col23;
-		row3->at(3) = col33;
-		return answer;
-	}
-
-	template<>
-	inline FMatDsptr EulerParameters<double>::pCpEtimesColumn(FColDsptr col)
-	{
-		//"col size = 4."
-		auto c0 = col->at(0);
-		auto c1 = col->at(1);
-		auto c2 = col->at(2);
-		auto mc0 = -c0;
-		auto mc1 = -c1;
-		auto mc2 = -c2;
-		auto mc3 = -col->at(3);
-		auto answer = std::make_shared<FullMatrix<double>>(3, 4);
-		auto& row0 = answer->at(0);
-		auto& row1 = answer->at(1);
-		auto& row2 = answer->at(2);
-		row0->atiput(0, mc3);
-		row0->atiput(1, mc2);
-		row0->atiput(2, c1);
-		row0->atiput(3, c0);
-		row1->atiput(0, c2);
-		row1->atiput(1, mc3);
-		row1->atiput(2, mc0);
-		row1->atiput(3, c1);
-		row2->atiput(0, mc1);
-		row2->atiput(1, c0);
-		row2->atiput(2, mc3);
-		row2->atiput(3, c2);
-		return answer;
-	}
-
-	template<typename T>
-	inline FMatDsptr EulerParameters<T>::pCTpEtimesColumn(FColDsptr col)
-	{
-		//"col size = 3."
-		auto c0 = col->at(0);
-		auto c1 = col->at(1);
-		auto c2 = col->at(2);
-		auto mc0 = -c0;
-		auto mc1 = -c1;
-		auto mc2 = -c2;
-		auto answer = std::make_shared<FullMatrix<double>>(4, 4);
-		auto& row0 = answer->at(0);
-		auto& row1 = answer->at(1);
-		auto& row2 = answer->at(2);
-		auto& row3 = answer->at(3);
-		row0->atiput(0, 0.0);
-		row0->atiput(1, c2);
-		row0->atiput(2, mc1);
-		row0->atiput(3, c0);
-		row1->atiput(0, mc2);
-		row1->atiput(1, 0.0);
-		row1->atiput(2, c0);
-		row1->atiput(3, c1);
-		row2->atiput(0, c1);
-		row2->atiput(1, mc0);
-		row2->atiput(2, 0.0);
-		row2->atiput(3, c2);
-		row3->atiput(0, mc0);
-		row3->atiput(1, mc1);
-		row3->atiput(2, mc2);
-		row3->atiput(3, 0.0);
-		return answer;
-	}
-
-	template<>
-	inline FMatFMatDsptr EulerParameters<double>::ppApEpEtimesMatrix(FMatDsptr mat)
-	{
-		FRowDsptr a2m0 = mat->at(0)->times(2.0);
-		FRowDsptr a2m1 = mat->at(1)->times(2.0);
-		FRowDsptr a2m2 = mat->at(2)->times(2.0);
-		FRowDsptr m2m0 = a2m0->negated();
-		FRowDsptr m2m1 = a2m1->negated();
-		FRowDsptr m2m2 = a2m2->negated();
-		FRowDsptr zero = std::make_shared<FullRow<double>>(3, 0.0);
-		auto mat00 = std::make_shared<FullMatrix<double>>(ListFRD{ a2m0, m2m1, m2m2 });
-		auto mat01 = std::make_shared<FullMatrix<double>>(ListFRD{ a2m1, a2m0, zero });
-		auto mat02 = std::make_shared<FullMatrix<double>>(ListFRD{ a2m2, zero, a2m0 });
-		auto mat03 = std::make_shared<FullMatrix<double>>(ListFRD{ zero, m2m2, a2m1 });
-		auto mat11 = std::make_shared<FullMatrix<double>>(ListFRD{ m2m0, a2m1, m2m2 });
-		auto mat12 = std::make_shared<FullMatrix<double>>(ListFRD{ zero, a2m2, a2m1 });
-		auto mat13 = std::make_shared<FullMatrix<double>>(ListFRD{ a2m2, zero, m2m0 });
-		auto mat22 = std::make_shared<FullMatrix<double>>(ListFRD{ m2m0, m2m1, a2m2 });
-		auto mat23 = std::make_shared<FullMatrix<double>>(ListFRD{ m2m1, a2m0, zero });
-		auto mat33 = std::make_shared<FullMatrix<double>>(ListFRD{ a2m0, a2m1, a2m2 });
-		auto answer = std::make_shared<FullMatrix<FMatDsptr>>(4, 4);
-		auto& row0 = answer->at(0);
-		row0->at(0) = mat00;
-		row0->at(1) = mat01;
-		row0->at(2) = mat02;
-		row0->at(3) = mat03;
-		auto& row1 = answer->at(1);
-		row1->at(0) = mat01;
-		row1->at(1) = mat11;
-		row1->at(2) = mat12;
-		row1->at(3) = mat13;
-		auto& row2 = answer->at(2);
-		row2->at(0) = mat02;
-		row2->at(1) = mat12;
-		row2->at(2) = mat22;
-		row2->at(3) = mat23;
-		auto& row3 = answer->at(3);
-		row3->at(0) = mat03;
-		row3->at(1) = mat13;
-		row3->at(2) = mat23;
-		row3->at(3) = mat33;
-		return answer;
-	}
-
-	template<>
-	inline void EulerParameters<double>::initialize()
-	{
-		aA = FullMatrix<double>::identitysptr(3);
-		aB = std::make_shared<FullMatrix<double>>(3, 4);
-		aC = std::make_shared<FullMatrix<double>>(3, 4);
-		pApE = std::make_shared<FullColumn<FMatDsptr>>(4);
-		for (int i = 0; i < 4; i++)
-		{
-			pApE->at(i) = std::make_shared<FullMatrix<double>>(3, 3);
-		}
-	}
-	template<typename T>
-	inline void EulerParameters<T>::calc()
-	{
-		this->calcABC();
-		this->calcpApE();
-	}
-	template<>
-	inline void EulerParameters<double>::calcABC()
-	{
-		double aE0 = this->at(0);
-		double aE1 = this->at(1);
-		double aE2 = this->at(2);
-		double aE3 = this->at(3);
-		double mE0 = -aE0;
-		double mE1 = -aE1;
-		double mE2 = -aE2;
-		FRowDsptr aBi;
-		aBi = aB->at(0);
-		aBi->at(0) = aE3;
-		aBi->at(1) = mE2;
-		aBi->at(2) = aE1;
-		aBi->at(3) = mE0;
-		aBi = aB->at(1);
-		aBi->at(0) = aE2;
-		aBi->at(1) = aE3;
-		aBi->at(2) = mE0;
-		aBi->at(3) = mE1;
-		aBi = aB->at(2);
-		aBi->at(0) = mE1;
-		aBi->at(1) = aE0;
-		aBi->at(2) = aE3;
-		aBi->at(3) = mE2;
-		FRowDsptr aCi;
-		aCi = aC->at(0);
-		aCi->at(0) = aE3;
-		aCi->at(1) = aE2;
-		aCi->at(2) = mE1;
-		aCi->at(3) = mE0;
-		aCi = aC->at(1);
-		aCi->at(0) = mE2;
-		aCi->at(1) = aE3;
-		aCi->at(2) = aE0;
-		aCi->at(3) = mE1;
-		aCi = aC->at(2);
-		aCi->at(0) = aE1;
-		aCi->at(1) = mE0;
-		aCi->at(2) = aE3;
-		aCi->at(3) = mE2;
-
-		aA = aB->timesTransposeFullMatrix(aC);
-	}
-	template<>
-	inline void EulerParameters<double>::calcpApE()
-	{
-		double a2E0 = 2.0 * (this->at(0));
-		double a2E1 = 2.0 * (this->at(1));
-		double a2E2 = 2.0 * (this->at(2));
-		double a2E3 = 2.0 * (this->at(3));
-		double m2E0 = -a2E0;
-		double m2E1 = -a2E1;
-		double m2E2 = -a2E2;
-		double m2E3 = -a2E3;
-		FMatDsptr pApEk;
-		pApEk = pApE->at(0);
-		FRowDsptr pAipEk;
-		pAipEk = pApEk->at(0);
-		pAipEk->at(0) = a2E0;
-		pAipEk->at(1) = a2E1;
-		pAipEk->at(2) = a2E2;
-		pAipEk = pApEk->at(1);
-		pAipEk->at(0) = a2E1;
-		pAipEk->at(1) = m2E0;
-		pAipEk->at(2) = m2E3;
-		pAipEk = pApEk->at(2);
-		pAipEk->at(0) = a2E2;
-		pAipEk->at(1) = a2E3;
-		pAipEk->at(2) = m2E0;
-		//
-		pApEk = pApE->at(1);
-		pAipEk = pApEk->at(0);
-		pAipEk->at(0) = m2E1;
-		pAipEk->at(1) = a2E0;
-		pAipEk->at(2) = a2E3;
-		pAipEk = pApEk->at(1);
-		pAipEk->at(0) = a2E0;
-		pAipEk->at(1) = a2E1;
-		pAipEk->at(2) = a2E2;
-		pAipEk = pApEk->at(2);
-		pAipEk->at(0) = m2E3;
-		pAipEk->at(1) = a2E2;
-		pAipEk->at(2) = m2E1;
-		//
-		pApEk = pApE->at(2);
-		pAipEk = pApEk->at(0);
-		pAipEk->at(0) = m2E2;
-		pAipEk->at(1) = m2E3;
-		pAipEk->at(2) = a2E0;
-		pAipEk = pApEk->at(1);
-		pAipEk->at(0) = a2E3;
-		pAipEk->at(1) = m2E2;
-		pAipEk->at(2) = a2E1;
-		pAipEk = pApEk->at(2);
-		pAipEk->at(0) = a2E0;
-		pAipEk->at(1) = a2E1;
-		pAipEk->at(2) = a2E2;
-		//
-		pApEk = pApE->at(3);
-		pAipEk = pApEk->at(0);
-		pAipEk->at(0) = a2E3;
-		pAipEk->at(1) = m2E2;
-		pAipEk->at(2) = a2E1;
-		pAipEk = pApEk->at(1);
-		pAipEk->at(0) = a2E2;
-		pAipEk->at(1) = a2E3;
-		pAipEk->at(2) = m2E0;
-		pAipEk = pApEk->at(2);
-		pAipEk->at(0) = m2E1;
-		pAipEk->at(1) = a2E0;
-		pAipEk->at(2) = a2E3;
-	}
-	template<typename T>
-	inline void EulerParameters<T>::conditionSelf()
-	{
-		EulerArray<T>::conditionSelf();
-		this->normalizeSelf();
-	}
 }
-

OndselSolver/EulerParameters.ref.h

@@ -0,0 +1,6 @@
+#pragma once
+
+namespace MbD {
+    template<typename T>
+    class EulerParameters;
+}

OndselSolver/EulerParametersDot.h

@@ -55,13 +55,13 @@ namespace MbD {
 	template<typename T>
 	inline void EulerParametersDot<T>::initialize()
 	{
-		aAdot = std::make_shared<FullMatrix<double>>(3, 3);
-		aBdot = std::make_shared<FullMatrix<double>>(3, 4);
-		aCdot = std::make_shared<FullMatrix<double>>(3, 4);
+		aAdot = std::make_shared<FullMatrixDouble>(3, 3);
+		aBdot = std::make_shared<FullMatrixDouble>(3, 4);
+		aCdot = std::make_shared<FullMatrixDouble>(3, 4);
 		pAdotpE = std::make_shared<FullColumn<FMatDsptr>>(4);
 		for (int i = 0; i < 4; i++)
 		{
-			pAdotpE->at(i) = std::make_shared<FullMatrix<double>>(3, 3);
+			pAdotpE->at(i) = std::make_shared<FullMatrixDouble>(3, 3);
 		}
 	}
 
@@ -185,6 +185,7 @@ namespace MbD {
 	inline FColDsptr EulerParametersDot<T>::omeOpO()
 	{
 		auto aaa = this->aB();
+//		auto bbb = aaa->timesFullColumn((MbD::FColsptr<double>)this);
 		auto bbb = aaa->timesFullColumn(this);
 		auto ccc = bbb->times(2.0);
 		return ccc;

OndselSolver/Exponential.cpp

@@ -19,8 +19,13 @@ double MbD::Exponential::getValue()
     return std::log(xx->getValue());
 }
 
+Symsptr MbD::Exponential::copyWith(Symsptr arg)
+{
+	return std::make_shared<Exponential>(arg);
+}
+
 std::ostream& MbD::Exponential::printOn(std::ostream& s) const
 {
-	s << "exp(" << xx << ")";
+	s << "exp(" << *xx << ")";
 	return s;
 }

OndselSolver/Exponential.h

@@ -18,6 +18,7 @@ namespace MbD {
         Exponential() = default;
         Exponential(Symsptr arg);
         double getValue() override;
+        Symsptr copyWith(Symsptr arg) override;
 
         std::ostream& printOn(std::ostream& s) const override;
 

OndselSolver/ExpressionX.cpp

@@ -7,6 +7,7 @@
  ***************************************************************************/
  
 #include "ExpressionX.h"
+#include "Constant.h"
 
 using namespace MbD;
 
@@ -23,8 +24,14 @@ void MbD::ExpressionX::xexpression(Symsptr arg, Symsptr func)
 	expression = func;
 }
 
+Symsptr MbD::ExpressionX::differentiateWRTx()
+{
+	return expression->differentiateWRT(xx);
+}
+
 Symsptr MbD::ExpressionX::differentiateWRT(Symsptr var)
 {
+	if (this == var.get()) return sptrConstant(1.0);
 	return expression->differentiateWRT(var);
 }
 

OndselSolver/ExpressionX.h

@@ -17,6 +17,7 @@ namespace MbD {
     public:
 
         void xexpression(Symsptr arg, Symsptr func);
+        Symsptr differentiateWRTx() override;
         Symsptr differentiateWRT(Symsptr var) override;
         double getValue() override;
 

OndselSolver/FullColumn.cpp

@@ -9,5 +9,196 @@
 #include <sstream>
 
 #include "FullColumn.h"
+#include "FullRow.h"
+#include "FullMatrix.h"
 
-using namespace MbD;
+namespace MbD {
+    template<typename T>
+    FColsptr<T> FullColumn<T>::plusFullColumn(FColsptr<T> fullCol)
+    {
+        int n = (int) this->size();
+        auto answer = std::make_shared<FullColumn<T>>(n);
+        for (int i = 0; i < n; i++) {
+            answer->at(i) = this->at(i) + fullCol->at(i);
+        }
+        return answer;
+    }
+    template<typename T>
+    FColsptr<T> FullColumn<T>::minusFullColumn(FColsptr<T> fullCol)
+    {
+        int n = (int) this->size();
+        auto answer = std::make_shared<FullColumn<T>>(n);
+        for (int i = 0; i < n; i++) {
+            answer->at(i) = this->at(i) - fullCol->at(i);
+        }
+        return answer;
+    }
+    template<>
+    FColDsptr FullColumn<double>::times(double a)
+    {
+        int n = (int)this->size();
+        auto answer = std::make_shared<FullColumn<double>>(n);
+        for (int i = 0; i < n; i++) {
+            answer->at(i) = this->at(i) * a;
+        }
+        return answer;
+    }
+    template<typename T>
+    FColsptr<T> FullColumn<T>::times(T a)
+    {
+        assert(false);
+        auto answer = std::make_shared<FullColumn<T>>();
+        return answer;
+    }
+    template<typename T>
+    FColsptr<T> FullColumn<T>::negated()
+    {
+        return this->times(-1.0);
+    }
+    template<typename T>
+    void FullColumn<T>::atiputFullColumn(int i, FColsptr<T> fullCol)
+    {
+        for (int ii = 0; ii < fullCol->size(); ii++)
+        {
+            this->at(i + ii) = fullCol->at(ii);
+        }
+    }
+    template<typename T>
+    void FullColumn<T>::atiplusFullColumn(int i, FColsptr<T> fullCol)
+    {
+        for (int ii = 0; ii < fullCol->size(); ii++)
+        {
+            this->at(i + ii) += fullCol->at(ii);
+        }
+    }
+    template<typename T>
+    void FullColumn<T>::equalSelfPlusFullColumnAt(FColsptr<T> fullCol, int ii)
+    {
+        //self is subcolumn of fullCol
+        for (int i = 0; i < this->size(); i++)
+        {
+            this->at(i) += fullCol->at(ii + i);
+        }
+    }
+    template<typename T>
+    void FullColumn<T>::atiminusFullColumn(int i1, FColsptr<T> fullCol)
+    {
+        for (int ii = 0; ii < fullCol->size(); ii++)
+        {
+            int i = i1 + ii;
+            this->at(i) -= fullCol->at(ii);
+        }
+    }
+    template<typename T>
+    void FullColumn<T>::equalFullColumnAt(FColsptr<T> fullCol, int i)
+    {
+        this->equalArrayAt(fullCol, i);
+        //for (int ii = 0; ii < this->size(); ii++)
+        //{
+        //	this->at(ii) = fullCol->at(i + ii);
+        //}
+    }
+    template<>
+    FColDsptr FullColumn<double>::copy()
+    {
+        auto n = (int) this->size();
+        auto answer = std::make_shared<FullColumn<double>>(n);
+        for (int i = 0; i < n; i++)
+        {
+            answer->at(i) = this->at(i);
+        }
+        return answer;
+    }
+    template<typename T>
+    FRowsptr<T> FullColumn<T>::transpose()
+    {
+        return std::make_shared<FullRow<T>>(*this);
+    }
+    template<typename T>
+    void FullColumn<T>::atiplusFullColumntimes(int i1, FColsptr<T> fullCol, T factor)
+    {
+        for (int ii = 0; ii < fullCol->size(); ii++)
+        {
+            int i = i1 + ii;
+            this->at(i) += fullCol->at(ii) * factor;
+        }
+    }
+    template<typename T>
+    T FullColumn<T>::transposeTimesFullColumn(const FColsptr<T> fullCol)
+    {
+        return this->dot(fullCol);
+    }
+    template<typename T>
+    void FullColumn<T>::equalSelfPlusFullColumntimes(FColsptr<T> fullCol, T factor)
+    {
+        this->equalSelfPlusFullVectortimes(fullCol, factor);
+    }
+    template<typename T>
+    FColsptr<T> FullColumn<T>::cross(FColsptr<T> fullCol)
+    {
+        auto a0 = this->at(0);
+        auto a1 = this->at(1);
+        auto a2 = this->at(2);
+        auto b0 = fullCol->at(0);
+        auto b1 = fullCol->at(1);
+        auto b2 = fullCol->at(2);
+        auto answer = std::make_shared<FullColumn<T>>(3);
+        answer->atiput(0, a1 * b2 - (a2 * b1));
+        answer->atiput(1, a2 * b0 - (a0 * b2));
+        answer->atiput(2, a0 * b1 - (a1 * b0));
+        return answer;
+    }
+    //template<>
+    //inline std::shared_ptr<FullColumn<Symsptr>> FullColumn<Symsptr>::simplified()
+    //{
+    //	auto n = this->size();
+    //	auto answer = std::make_shared<FullColumn<Symsptr>>(n);
+    //	for (int i = 0; i < n; i++)
+    //	{
+    //		auto func = this->at(i);
+    //		answer->at(i) = func->simplified(func);
+    //	}
+    //	return answer;
+    //}
+    template<typename T>
+    FColsptr<T> FullColumn<T>::simplified()
+    {
+//        assert(false);
+        return FColsptr<T>();
+    }
+    // instantiate on purpose to make visible in library api:
+    template class FullColumn<double>;
+    template class FullColumn<int>;
+    template<typename T>
+    double FullColumn<T>::dot(std::shared_ptr<FullVector<T>> vec)
+    {
+        int n = (int)this->size();
+        double answer = 0.0;
+        for (int i = 0; i < n; i++) {
+            answer += this->at(i) * vec->at(i);
+        }
+        return answer;
+    }
+    template<typename T>
+    std::shared_ptr<FullVector<T>> FullColumn<T>::dot(std::shared_ptr<std::vector<std::shared_ptr<FullColumn<T>>>> vecvec)
+    {
+        int ncol = (int)this->size();
+        auto nelem = vecvec->at(0)->size();
+        auto answer = std::make_shared<FullVector<T>>(nelem);
+        for (int k = 0; k < nelem; k++) {
+            auto sum = 0.0;
+            for (int i = 0; i < ncol; i++)
+            {
+                sum += this->at(i) * vecvec->at(i)->at(k);
+            }
+            answer->at(k) = sum;
+        }
+        return answer;
+    }
+    template<typename T>
+    std::shared_ptr<FullColumn<T>> FullColumn<T>::clonesptr()
+    {
+        return std::make_shared<FullColumn<T>>(*this);
+    }
+
+}

OndselSolver/FullColumn.h

@@ -5,24 +5,19 @@
  *                                                                         *
  *   See LICENSE file for details about copyright.                         *
  ***************************************************************************/
- 
+
 #pragma once
 
 #include <string>
 #include <sstream> 
 
 #include "FullVector.h"
+#include "FullColumn.ref.h"
+#include "FullRow.ref.h"
+#include "FullMatrix.ref.h"
+//#include "Symbolic.h"
 
 namespace MbD {
-	template<typename T>
-	class FullColumn;
-	using FColDsptr = std::shared_ptr<FullColumn<double>>;
-	template<typename T>
-	using FColsptr = std::shared_ptr<FullColumn<T>>;
-	template<typename T>
-	class FullRow;
-	template<typename T>
-	using FRowsptr = std::shared_ptr<FullRow<T>>;
 	class Symbolic;
 
 	template<typename T>
@@ -36,6 +31,7 @@ namespace MbD {
 		FullColumn(typename std::vector<T>::iterator begin, typename std::vector<T>::iterator end) : FullVector<T>(begin, end) {}
 		FullColumn(std::initializer_list<T> list) : FullVector<T>{ list } {}
 		FColsptr<T> plusFullColumn(FColsptr<T> fullCol);
+		FColsptr<T> plusFullColumntimes(FColsptr<T> fullCol, double factor);
 		FColsptr<T> minusFullColumn(FColsptr<T> fullCol);
 		FColsptr<T> times(T a);
 		FColsptr<T> negated();
@@ -47,175 +43,27 @@ namespace MbD {
 		FColsptr<T> copy();
 		FRowsptr<T> transpose();
 		void atiplusFullColumntimes(int i, FColsptr<T> fullCol, T factor);
-		T transposeTimesFullColumn(const FColsptr<T> fullCol);		
+		T transposeTimesFullColumn(const FColsptr<T> fullCol);
 		void equalSelfPlusFullColumntimes(FColsptr<T> fullCol, T factor);
 		FColsptr<T> cross(FColsptr<T> fullCol);
 		FColsptr<T> simplified();
+		std::shared_ptr<FullColumn<T>> clonesptr();
+		double dot(std::shared_ptr<FullVector<T>> vec);
+		std::shared_ptr<FullVector<T>> dot(std::shared_ptr<std::vector<std::shared_ptr<FullColumn<T>>>> vecvec);
 
 		std::ostream& printOn(std::ostream& s) const override;
 	};
-	template<typename T>
-	inline FColsptr<T> FullColumn<T>::plusFullColumn(FColsptr<T> fullCol)
-	{
-		int n = (int) this->size();
-		auto answer = std::make_shared<FullColumn<T>>(n);
-		for (int i = 0; i < n; i++) {
-			answer->at(i) = this->at(i) + fullCol->at(i);
-		}
-		return answer;
-	}
-	template<typename T>
-	inline FColsptr<T> FullColumn<T>::minusFullColumn(FColsptr<T> fullCol)
-	{
-		int n = (int) this->size();
-		auto answer = std::make_shared<FullColumn<T>>(n);
-		for (int i = 0; i < n; i++) {
-			answer->at(i) = this->at(i) - fullCol->at(i);
-		}
-		return answer;
-	}
-	template<>
-	inline FColDsptr FullColumn<double>::times(double a)
-	{
-		int n = (int)this->size();
-		auto answer = std::make_shared<FullColumn<double>>(n);
-		for (int i = 0; i < n; i++) {
-			answer->at(i) = this->at(i) * a;
-		}
-		return answer;
-	}
-	template<typename T>
-	inline FColsptr<T> FullColumn<T>::times(T a)
-	{
-		assert(false);
-	}
-	template<typename T>
-	inline FColsptr<T> FullColumn<T>::negated()
-	{
-		return this->times(-1.0);
-	}
-	template<typename T>
-	inline void FullColumn<T>::atiputFullColumn(int i, FColsptr<T> fullCol)
-	{
-		for (int ii = 0; ii < fullCol->size(); ii++)
-		{
-			this->at(i + ii) = fullCol->at(ii);
-		}
-	}
-	template<typename T>
-	inline void FullColumn<T>::atiplusFullColumn(int i, FColsptr<T> fullCol)
-	{
-		for (int ii = 0; ii < fullCol->size(); ii++)
-		{
-			this->at(i + ii) += fullCol->at(ii);
-		}
-	}
-	template<typename T>
-	inline void FullColumn<T>::equalSelfPlusFullColumnAt(FColsptr<T> fullCol, int ii)
-	{
-		//self is subcolumn of fullCol
-		for (int i = 0; i < this->size(); i++)
-		{
-			this->at(i) += fullCol->at(ii + i);
-		}
-	}
-	template<typename T>
-	inline void FullColumn<T>::atiminusFullColumn(int i1, FColsptr<T> fullCol)
-	{
-		for (int ii = 0; ii < fullCol->size(); ii++)
-		{
-			int i = i1 + ii;
-			this->at(i) -= fullCol->at(ii);
-		}
-	}
-	template<typename T>
-	inline void FullColumn<T>::equalFullColumnAt(FColsptr<T> fullCol, int i)
-	{
-		this->equalArrayAt(fullCol, i);
-		//for (int ii = 0; ii < this->size(); ii++)
-		//{
-		//	this->at(ii) = fullCol->at(i + ii);
-		//}
-	}
-	template<>
-	inline FColDsptr FullColumn<double>::copy()
-	{
-		auto n = (int) this->size();
-		auto answer = std::make_shared<FullColumn<double>>(n);
-		for (int i = 0; i < n; i++)
-		{
-			answer->at(i) = this->at(i);
-		}
-		return answer;
-	}
-	template<typename T>
-	inline FRowsptr<T> FullColumn<T>::transpose()
-	{
-		return std::make_shared<FullRow<T>>(*this);
-	}
-	template<typename T>
-	inline void FullColumn<T>::atiplusFullColumntimes(int i1, FColsptr<T> fullCol, T factor)
-	{
-		for (int ii = 0; ii < fullCol->size(); ii++)
-		{
-			int i = i1 + ii;
-			this->at(i) += fullCol->at(ii) * factor;
-		}
-	}
-	template<typename T>
-	inline T FullColumn<T>::transposeTimesFullColumn(const FColsptr<T> fullCol)
-	{
-		return this->dot(fullCol);
-	}
-	template<typename T>
-	inline void FullColumn<T>::equalSelfPlusFullColumntimes(FColsptr<T> fullCol, T factor)
-	{
-		this->equalSelfPlusFullVectortimes(fullCol, factor);
-	}
-	template<typename T>
-	inline FColsptr<T> FullColumn<T>::cross(FColsptr<T> fullCol)
-	{
-			auto a0 = this->at(0);
-			auto a1 = this->at(1);
-			auto a2 = this->at(2);
-			auto b0 = fullCol->at(0);
-			auto b1 = fullCol->at(1);
-			auto b2 = fullCol->at(2);
-			auto answer = std::make_shared<FullColumn<T>>(3);
-			answer->atiput(0, a1 * b2 - (a2 * b1));
-			answer->atiput(1, a2 * b0 - (a0 * b2));
-			answer->atiput(2, a0 * b1 - (a1 * b0));
-			return answer;
-	}
-	//template<>
-	//inline std::shared_ptr<FullColumn<Symsptr>> FullColumn<Symsptr>::simplified()
-	//{
-	//	auto n = this->size();
-	//	auto answer = std::make_shared<FullColumn<Symsptr>>(n);
-	//	for (int i = 0; i < n; i++)
-	//	{
-	//		auto func = this->at(i);
-	//		answer->at(i) = func->simplified(func);
-	//	}
-	//	return answer;
-	//}
-	template<typename T>
-	inline FColsptr<T> FullColumn<T>::simplified()
-	{
-		assert(false);
-		return FColsptr<T>();
-	}
-	template<typename T>
-	inline std::ostream& FullColumn<T>::printOn(std::ostream& s) const
-	{
-		s << "FullCol{";
-		s << this->at(0);
-		for (int i = 1; i < this->size(); i++)
-		{
-			s << ", " << this->at(i);
-		}
-		s << "}";
-		return s;
-	}
+    // the following "printOn" needs to be in the header for unknown reasons linker
+    template<typename T>
+    std::ostream& FullColumn<T>::printOn(std::ostream& s) const
+    {
+        s << "FullCol{";
+        s << this->at(0);
+        for (int i = 1; i < int(this->size()); i++)
+        {
+            s << ", " << this->at(i);
+        }
+        s << "}";
+        return s;
+    }
 }
-

OndselSolver/FullColumn.ref.h

@@ -0,0 +1,14 @@
+#pragma once
+
+#include <algorithm>
+#include <memory>
+
+namespace MbD {
+    template<typename T>
+    class FullColumn;
+
+    using FColDsptr = std::shared_ptr<FullColumn<double>>;
+
+    template<typename T>
+    using FColsptr = std::shared_ptr<FullColumn<T>>;
+}