text:
\subsection{Mathematical Model}

Since the depth information ($z$-coordinate) in this experiment is to be obtained from the equation of the surface using the $x$ and $y$ coordinates extracted from Tracker, the knowledge of the exact mathematical form of the printed surface is imperative.
However, determining the same is a slightly non-trivial issue, as one may reason from simple dimensional considerations. This $z = 1/r$ surface designed as above, must actually be of the form $z = \kappa/r$, where $\kappa$ has dimensions of area, and is just numerically equal to 1 in our surface specification at the beginning. In order to avoid distorting the shape of the sample while scaling the prototype of the 3D-print to larger dimensions, an important precaution we observed was to lock the aspect ratio so that all three coordinates scale equally through some $L_{\rm scale}$ factor. But even when lengths scale equally, area would scale as $L_{\rm scale}^2$, so that $z = \kappa/r$ correctly scales by $L_{\rm scale}$ itself, hence preserving the aspect ratio as it should. However, on careful reflection, this also entails that on scaling the prototype, while the surface still remains a $(1/r)$ surface, the mathematical form of the surface does not remain the same $z = \kappa/r$ (with the original value of $\kappa = 1$) due to unequal scaling of $\kappa$ and $r$ in this expression.
Therefore,
we model this printed surface to be
mathematically represented
as
$z = C_1 +
(C_2/r) $,
where $C_2 = - \vert C_2 \vert$ is a different constant than the original $\kappa$, and $C_1$ is another constant introduced owing to spatial translation due to the adoption of any particular choice of origin, or coordinate frame. These constants can be evaluated in two ways: mathematically, using sample geometry as an input, and experimentally, through a sampling of the surface heights against a calibrated horizontal axis. We describe the latter (experimental) approach in the next section. As for the former approach, we adopt a coordinate system shown in Fig. \ref{3dprinttheoryannotations} and choose $z = 0$ at the highest elevated points of the surface, which are the four diagonal vertices of the sample. With reference to this choice of coordinates, the constants can be evaluated to be:
\begin{equation}
C_2 = - \vert z_{\rm depth}\vert \frac{R_{\rm out} \ R_{\rm in}}{\left( R_{\rm out} -R_{\rm in} \right)}
\end{equation}
\begin{equation}
C_1 = - \frac{C_2}{R_{\rm out}} = \vert z_{\rm depth}\vert \frac{R_{\rm in}}{\left( R_{\rm out} -R_{\rm in} \right)}
\end{equation}
Substituting the values of $R_{\rm in} = 2$ cm, and $R_{\rm out} = 15.5/\sqrt{2} = 10.96$ cm, and $z_{\rm depth} = 4.756$ cm, we obtain the following theoretical estimates for the surface constants.
\begin{equation}
C_1 = 1.0616 \ {\rm cm}, \ {\rm and} \ C_2 = -11.6357 \ {\rm cm}^2
\end{equation}
Complete details of the calculation can be found in the Supplementary Information \cite{Supplement} to this article.

\begin{figure}
\begin{center}
\scalebox{0.57}{\includegraphics{Fig03_Type02_Inkscape.pdf}}\caption{ \label{3dprinttheoryannotations} (Color Online) The geometry and coordinate choice for surface calibration. As before, for clarity, only the inner surface is shown; the complete sample region with the underlying fill is represented as dotted lines. $z=0$ is chosen at the highest elevated point of the sample, corresponding to the diagonal edges. At the bottom surface $(z = z_{\rm depth} = - \vert z_{\rm depth} \vert)$, the surface carves out a circle of radius $R_{\rm in}$.}
\end{center}
\end{figure}

\subsection{Experimental Calibration}

The experimental surface-calibration process requires us to sample the surface at different points in the $x-z$ plane, and to fit a curve through these points. However, this requires a cross sectional view of the inner surface of the sample, which is not accessible for direct experimentation. We adapt the experimental arrangement to provide for the same, by adding an additional reference mark provided by a pin-thermocol arrangement placed underneath the sample, which allows us to circumvent these difficulties and perform the calibration experimentally. We refer the interested reader to the Supplementary Information \cite{Supplement} to this article, where complete details about the experimental assembly have been collected. After collecting various data points/sets as desired, we best fit the data to the expected curve to determine the values of the constants. This can be done through multiple procedures, all of them mutually equivalent:
\begin{enumerate}
\item Writing a simple routine for least squares fitting in any programming language/software package, upon importing this experimental $(x_i, z_i')$ data.
\item Importing this data in Scilab \cite{Scilab_site}, and using the built-in \emph{reglin} subroutine for linear regression, directly yielding the surface constants.
\item Simply entering this experimental $(x_i, z_i')$ data in Microsoft Excel, and using Excel Solver to find an optimal fit. This can be easily done for both $z'$-versus-$(1/x)$, as well as for $z'$-versus-$x$ variation. The procedure requires some initial guesses for $C_2$ and $C_3$ to be provided, for which the theoretical estimates obtained earlier  prove very useful. The Solver optimizer searches for optimal regression solutions in the vicinity of the guesses we provide, and typically returns surface constants' values only slightly different from the theoretical values obtained earlier.
\end{enumerate}

We find the experimental values of the surface constants to be as follows:
\begin{equation}
C_2 = -12.0717 \ {\rm cm}^2, \ \ \ \ C_1 = \vert = 1.0646 \ {\rm cm}
\end{equation}
Thus, we find that the theoretical and experimental values of these surface constants are in reasonably good agreement. The experimental value of $C_3$ in particular is in excellent agreement with the theoretical estimate, while the experimental value of $C_2$ is also found to be within $4 \%$ of the theoretical estimate. This agreement lends a lot of credence to this model formulated.

\end{section}

\begin{section}{Critical Assessment of Kepler's First Law for Marble Motion}
\label{Critical_Assessment_KeplerFirstLaw_section}

As discussed in Section 1, whether or not the marble motion preserves Kepler's first law is a key question, which can now be simply answered on the basis of the data recorded from this experiment. For this, one can use a simple trick from vector analysis. One can use just three consecutive data points from the experimentally recorded datasets, two compute two consecutive distance vectors whose cross product can be used to identify the normal vector to the plane constituted by them.  If the points are distributed in a plane, the normal vectors so identified will be parallel to one another, and vice-versa. Moreover, from the information available, one can also compute the angles between these normal vectors to gauge the extent of coplanarity of these points.

The typical results are shown in Figs. 4 -- 6, showing respectively the typical orbital trajectories (for an elliptical-like, closed) orbit, the distribution of the normal vectors to the triplets of consecutive data points, and the distribution of the angles between consecutive normal vectors. Fig. 4 seems to suggest that the motion of points is approximately in a plane to a very good approximation, but Figs. 5 and 6 suggest that there is indeed a small variation from complete planarity. The largest deviation angle is $\sim 30$ for a single data point, but the typical deviation angles range around $< 10^{\circ}$, which is significantly small. But howsoever small this angle of deviation may be, in a strict mathematical sense, this is still an explicit experimental verification of the theoretical argument \cite{Middleton} for this motion of the ball to be distinct from Keplerian orbits.

However, the results of Fig. 4 are still reassuring in that one can see an approximate restoration of the initial coordinates at the end of the orbit. Considering that we are using a finite sized marble, we can envisage in general two kinds of orbit closures - a strong form and a weak form. A strict ``fully-closed'' orbit would mean the complete restoration of the initial coordinates $(x, y, z)$ at the end of the orbit, a very strict condition which upon a large number of trials, we do occasionally observe to be satisfied within a single decimal place or two. A majority of the trials result in elliptical arcs, rather than these fully closed orbits, which is understandable since the dissipative effects of friction can not be completely eliminated in this experiment. However, we can also have an intermediate subset of elliptical trajectories, characterized henceforth as ``near-closed'' orbits, which we define as pertaining to the condition $\left( {\vec r_f} - {\vec r_i} \right) < R_{\rm ball}$, or even some percentage of $R_{\rm ball}$. This generous definition accommodates a much bigger subset of the observed elliptical trajectories, in the category of closed orbits. While not very rigorous, this still serves our purpose very well in the present work, since ultimately, one requires merely five points on an ellipse to uniquely determine it \cite{Anton, MSE_5pointsdetermineellipse}, so in principle, each of these elliptical trajectories is actually an overdetermined system.

\begin{figure}
\begin{center}
\scalebox{0.4}{\includegraphics{Fig5_2_3D_Points_in_Space.pdf}}\caption{
(Color Online) Scatter plot of the experimentally observed data points, in a typical dataset. One observes approximate orbit closure.}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\scalebox{0.4}{\includegraphics{Fig6_2_Normal_Vectors.pdf}}\caption{
(Color Online) The distribution of the normal vectors of consecutive data points. Each pair gives a normal vector which points in a slightly different direction. For reference, the normal vector of the best fitted plane is also shown, annotated in red.}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\scalebox{0.4}{\includegraphics{Fig6_1_Distribution_of_Angles.pdf}}\caption{
(Color Online) The distribution of the angles between consecutive normal vectors. The maximum variation between consecutive angles is $\sim 30^{\circ}$, though the typical deviation angles are small.}
\end{center}
\end{figure}

However, these results resoundingly indicate that even though the motion is strictly speaking not in a plane, the observed data points are still scattered about some best fitted plane, with marginal deviations. Thus, although this is strictly speaking, not motion in a plane, but the deviations from perfect coplanarity is small, implying therefore that these data points are not too far from being in a plane.

In the following, we suggest a remedial procedure which may be undertaken to still approximately salvage Kepler's first law for the purpose of making a computational orbital mechanics lab, in spite of all the preceding discussion in this section.

\end{section}

\begin{section}{Resurrecting Kepler's First Law}

\label{Resurrection_Section_KeplerLaw1}

\subsection{Best Fitting a Plane to Data}
\label{Plane_fitting_section}

\subsubsection{Principal Component Analysis}
\label{PCA_subsubsection}

Plane fitting can be achieved through the powerful technique of Principal Component Analysis (PCA) \cite{Pearson_PCA_orgnl, Jollife_PCA_book, Shlens_TutorialPCA, Bishop_PatternRecog_ML, Smith_OtagoPCAReport, SE_MakingSenseofPCA, Gnanadesikan}. In essence, this method
relies on the identification of
the principal components/axes for the data \cite{Gnanadesikan}, through the aid of the covariance matrix, whose each element measures the covariance/correlation between the coordinates. The PCA procedure \cite{Jollife_PCA_book, Shlens_TutorialPCA} requires this matrix to be diagonalized, and its normalized eigenvectors to be extracted and sorted according to the corresponding eigenvalues. As regards regression problems, these eigenvectors carry a special geometrical significance. For instance, for linear regression involving two variables ($x, y$), the direction (eigenvector) corresponding to the largest eigenvalue (first principal component), maximizes the variance of data projection along this component, and minimizes the residual/projection error between the original data points and their projection on to this first principal direction \cite{Bishop_PatternRecog_ML, SE_MaxVarianceAlong_MinError}
thereby effectively solving the least square fitting problem. Orthogonality of the eigenvectors \cite{Weber_Arfken} entails that the other principal direction would then give the perpendicular to this best fitting straight line \cite{Smith_OtagoPCAReport}. Plane fitting in three dimensions is
an extension
of the same idea \cite{Pearson_PCA_orgnl}, wherein the normalized eigenvectors corresponding to the first two principal components, identify as the basis vectors of the best fitting plane, and minimize the sum of square of projection error of all these 3D points to this best fitting plane so identified.
The third principal component, corresponding to the smallest eigenvalue of the covariance matrix, identifies the outward normal vector ${\hat n}$ to this best fitting plane.

Applying this technique in our context,
we first `center' the data \cite{Pearson_PCA_orgnl}, by determining the centroid and
co-ordinates of data points relative to
it.
If ${\vec r}$ and ${\vec r_{\rm C}}$ respectively denote the position vectors of any experimental data point P, and the centroid C of all $N$ data points, then the position vector of P relative to C is ${\vec R_{\rm rel}} = \left({\vec r} - {\vec r_{\rm C}} \right)$.
In terms of its components $(X, Y, Z)$,
we can define \cite{Shlens_TutorialPCA} the covariance between, e.g. $X$ and $Y$ as
\begin{equation}
\sigma_{\rm XY}^{2} = \frac{1}{N} \sum_{i = 1}^N X_i Y_i, \label{correlationmatrix}
\end{equation}
Observing that $\sigma_{ij}^{2} = \sigma_{ji}^{2}$, with $\left\{ i, j \right\} \in\left\{  X,Y,Z \right\}$, and the diagonal term $\sigma_{jj}^{2}$ reduces to the usual variance $\sigma_{j}^{2}$ \cite{Weber_Arfken}, one can thus compose the symmetric, covariance matrix $\Sigma$, as
\begin{eqnarray}
\Sigma =\left[
\begin{array}
[c]{ccc}%
\sigma_{\rm X}^{2} & \sigma_{\rm XY}^{2} & \sigma_{\rm XZ}^{2}\\
\sigma_{\rm YX}^{2} & \sigma_{\rm Y}^{2} & \sigma_{\rm YZ}^{2}\\
\sigma_{\rm ZX}^{2} & \sigma_{\rm ZY}^{2} & \sigma_{\rm Z}^{2}
\end{array}
\right]
\label{Covariance_matrix}
\end{eqnarray}
While ab-initio numerical procedures for the implementation of the rest of the PCA procedure, viz. diagonalization and eigenvector-eigenvalue extraction, are well documented \cite{Jollife_PCA_book}, these can be readily performed with the aid of Scilab  \cite{Scilab_site}, with this ease being rooted in the powerful DSYEV subroutine of the ubiquitous linear algebra package LAPACK \cite{LAPACK}, thus making this procedure
easily accessible even to novices/undergraduates.

\subsubsection{Singular Value Decomposition}

A closely related approach, which is
also frequently used
in
regression analysis
\cite{Mandel_SVD} and can also be easily adopted
here, is to construct the full Singular Value Decomposition (SVD)
of the ${\vec R_{\rm rel}}$ obtained in Section \ref{PCA_subsubsection}.
In general terms, the SVD of a rectangular $\left( m \times n \right)$ matrix $A$ \cite{Bisgard, Jollife_PCA_book, Strang, Atkinson} refers to the decomposition $A = U {\mathbb S} V^T$, where $U$ and $V$ are orthogonal matrices of sizes $(m \times m)$ and $(n \times n)$ respectively, and ${\mathbb S}$ is a $(m \times n)$ matrix populated only along the diagonal entries,
bearing the `singular values' of $A$, which are simply the square roots of the eigenvalues of
 $A A^{\rm T}$.
Applying this SVD technique in our context, for plane fitting the centered data points, our equivalent $A$ matrix here would be a $(3 \times N)$ data matrix bearing the components of ${\vec R_{\rm rel}}$. Upon performing its SVD, the normal vector ${\hat n}$ of the best-fitting plane can be readily identified as the left singular vector (i.e. the column of $U$) corresponding to the least singular value (min$({{\mathbb S}_{\rm ii}}))$. Since $U$ is orthogonal, the other two columns of $U$ represent the null space of the vector, representing geometrically the basis vectors which span the best fitting plane. Once again, this computation of SVD of our data matrix can be greatly facilitated by Scilab \cite{Scilab_site}, wherein we can use the inbuilt DGESVD routine of the Linear Algebra package LAPACK \cite{LAPACK} to readily access the same.

The
implementation-routine for this decomposition seems to suggest that this is an alternative take on this regression problem, but in reality, this constitutes a completely equivalent regression strategy \cite{Mandel_SVD, Shlens_TutorialPCA, SE_RelationPCASVD_1, SE_RelationPCASVD_2}.
The PCA approach identifies (and sorts) mutually orthogonal axes as per the variance of the data along them, thereby identifying, most notably, the first principal component, viz. the special axis along which the variance of the data is the largest, or the data is clustered the most.
Thus, the normal vector ${\hat n}$ of the best-fitting plane is, very understandably, that orthogonal axis with the lowest variance of data-projection along it. As noted before, the singular values of the matrix form of ${\vec{R}_{\rm rel}}$, are in proportion with the eigenvalues of the correlation matrix defined through Eq. (\ref{correlationmatrix}), and hence, it is only natural that identification on the basis of lowest singular value in the SVD procedure, also results in the same normal vector ${\hat n}$. Thus, in spite of appearances, both procedures are a means to the same end. Although from a computational perspective, since the creation of the covariance matrix in the PCA procedure involves square/bilinear terms, it can potentially cause precision-loss, and hence,
instead of
a direct identification of the lowest eigenvalue solution on the basis of PCA, the equivalent SVD approach is generally considered preferable. However, in all datasets analyzed in the present work, both these procedures are observed to yield perfectly identical results for the best fitted plane, owing to this interrelation between PCA and SVD.

\subsection{Projection of Points on the Plane}
\label{Projection_points_plane_subsection}
While the procedures of Section (\ref{Plane_fitting_section}) identify the best fitting plane for the data points, in general, these data points do not strictly lie in this plane, but are distributed about it. Therefore, in order to formulate orbital mechanics within this plane, we must project these data points onto the plane, and consider these images (projections onto the plane) for subsequent analysis. As long as the deviation angle identified in Section (\ref{Critical_Assessment_KeplerFirstLaw_section}) is small, these images would very nearly correspond to the actual data points. In the following, we adopt the projection convention that the aforementioned `image' point is that point on the plane that is nearest to the actual experimental data point.

In linear regression, it is a familiar fact that the best fitted straight line passes through the centroid of the data points \cite{Pearson_PCA_orgnl, Chernov_Regression_Book, Chernov_Uniqueness}. This can be readily generalized to three dimensions, and thus, one may infer that the best fitted plane, identified in Section (\ref{Plane_fitting_section}), passes through the centroid of the data points. Thus, we know both, the normal vector $\left({\hat n}\right)$ to the best fitting plane,
and one point which surely lies on the plane. This much is sufficient information to determine the image points on the plane, in accordance with the following procedure \cite{Anton, SO_Projection, GeometricToolsBook}.
Invoking the ${\vec R_{\rm rel}}$ from Section (\ref{PCA_subsubsection}), which represents the position vector of any experimental point P relative to the centroid, $({\vec R_{\rm rel}} \cdot {\hat n})$ gives its projection on ${\hat n}$, or the distance from point to plane along the normal. Subtracting the corresponding projection vector $({\vec R_{\rm rel}} \cdot {\hat n}) {\hat n}$ from ${\vec R_{\rm rel}}$ returns the component of ${\vec R_{\rm rel}}$ perpendicular to $\hat n$, thus lying in the plane and representing the position vector of the image point relative to the centroid.

Then, relative to our original $(x,y,z)$ coordinate system, the position vector of the image point ($\vec r_{\rm (Im)}$) can then be constructed by subtracting the projection vector from ${\vec r}$, as
\begin{equation}
{\vec r_{\rm (Im)}} = {\vec r} - ({\vec R_{\rm rel}}\cdot {\hat n}) {\hat n} = {\vec r} - \left( \left({\vec r} - {\vec r_{\rm C}} \right) \cdot {\hat n} \right) {\hat n}
\end{equation}

\subsection{Transforming coordinates to the plane of motion: from 3D to 2D data}
\label{from3Dto2D}
Next, we seek to transform
the original $(x, y, z)$ coordinates of these image points to new coordinates within the plane identified above. This requires a reorientation of the co-ordinate system, such that the latter coordinates may be denoted as $(x_p, y_p, z_p)$, with $z_p = 0$, and hence, we can perform the subsequent orbital mechanical analysis on the two coordinates $(x_p, y_p)$ within the plane. This transformation may be achieved through a couple of equivalent procedures.

\subsubsection{Geometric Determination of Rotation Matrix}
\label{Geom_Det_Rot_Mat}
Mathematically, this transformation from the original to the new coordinates may be realized through rotations by three Euler angles $(\alpha, \beta, \gamma)$, eventually constituting the rotation matrix $\mathbf{R} = R_1\left({\hat z}, \alpha \right) R_2\left({\hat y}, \beta \right) R_3\left({\hat z}, \gamma \right)$ \cite{Weber_Arfken, Goldstein}. While it is usual to construct this cumulative rotation matrix, by individually working out these $R_i$ matrices, in this context, we are less concerned with the individual pieces and more with the end result $\mathbf{R}$, so an alternative procedure \cite{MSE_make2Dfrom3D} is more convenient. The transformation between the original and the projected coordinates may be represented as:
\begin{eqnarray}
\left[
\begin{array}
[c]{c}%
x_p \\
y_p \\
z_p %
\end{array}
\right] = \mathbf{R} \left[
\begin{array}
[c]{c}%
x \\
y \\
z %
\end{array}
\right] =\left[
\begin{array}
[c]{ccc}%
e_{1}^{x} & e_{2}^{x} & e_{3}^{x}\\
e_{1}^{y} & e_{2}^{y} & e_{3}^{y}\\
e_{1}^{z} & e_{2}^{z} & e_{3}^{z}%
\end{array}
\right]
\left[
\begin{array}
[c]{c}%
x \\
y \\
z %
\end{array}
\right]
\label{Eq_projcoords}
\end{eqnarray}
Invoking the properties \cite{Weber_Arfken} of the rotation matrix ${\mathbf R}$,
the vectors $e_{i}^{\upsilon}$ for $\upsilon \in\left\{  x,y,z\right\}  $ make an orthonormal set, satisfying
$
\sum_{i=1}^{3}e_{i}^{\upsilon_1}e_{i}^{\upsilon_2}= \delta_{\left( \upsilon_1, \upsilon_2 \right)}$. Using this, these vectors can be determined simply through vector calculus and geometric considerations. Taking any three points $P_1$, $P_2$, and $P_3$ in the plane of motion, each written as yet in terms of its $(x, y, z)$ coordinates, one can compute and normalize the vector $\overrightarrow{P_1 P_2}$, and use it as the unit vector  $e_{i}^{x}$ along the $x_p$ axis in the plane. Since this plane is to eventually be realized as the $x_p$-$y_p$  plane, unit vector $e_{i}^{z}$ along the $z_p$ axis must be perpendicular to any vectors within this plane. Since the previously computed $e_{i}^{x}$, and the vector $\overrightarrow{P_1 P_3}$, are two such vectors at hand, $e_{i}^{z}$ can be chosen to lie along the cross product of these two vectors, duly normalized. Thereafter, $e_{i}^{y}$ must be such that it is perpendicular to both $e_{i}^{x}$ and $e_{i}^{z}$, and hence, can be computed from their cross product. To summarize this geometrical procedure:
\begin{eqnarray}
\widehat{e^{x}} = \frac{\overrightarrow{P_1 P_2}}{\vert \vert \overrightarrow{P_1 P_2} \vert \vert}, \ \widehat{e^{z}} =  \frac{\widehat{e^{x}} \times \overrightarrow{P_1 P_3}}{\vert \vert \widehat{e^{x}} \times \overrightarrow{P_1 P_3} \vert \vert},
\ \widehat{e^{y}} =  \widehat{e^{z}} \times \widehat{e^{x}}
\label{eqn_geom_3Dto2D}
\end{eqnarray}
These equations determine ${\mathbf R}$ and the projected coordinates through the use of Eq.(\ref{Eq_projcoords}). It can be readily verified computationally that this procedure returns a $z_p \approx 0$ (typical numerical values $\sim 10^{-16}$) \cite{Github_OrbM}. One must note that since the points $P_i$ in Eq.(\ref{eqn_geom_3Dto2D}) are chosen arbitrarily out of all points in the plane, the direction of the $x_p$ axis ($\widehat{e^{x}}$) is also arbitrary within the plane. Thus, even if the $(x_p, y_p)$ data obtained through Eq. (\ref{Eq_projcoords}) represents an ellipse, it is, in general, not in the principal-axis form \cite{Anton}.

\subsubsection{Using Basis Vectors from SVD/PCA}
\label{SVDPCA_Det_Rot_Mat}
As seen previously in Section (\ref{Plane_fitting_section}), the normal vector ${\hat n}$ of the best fitting plane, as well as two orthogonal basis vectors which span the best fitting plane, can all be directly extracted from both, the PCA and SVD approaches. These three vectors can be directly utilized to construct the rotation matrix $\mathbf{R}$, and  the projected coordinates can again be computed through Eq. (\ref{Eq_projcoords}). Unsurprisingly, this procedure similarly returns $z_p \approx 0$ (typical numerical values $\sim 10^{-16}$) \cite{Github_OrbM}.

A few comments regarding how this compares against the procedure of Sec.(\ref{Geom_Det_Rot_Mat}) are in order.
Since any two mutually orthogonal unit vectors in the best fitted plane can qualify as the basis vectors in that plane, an infinite number of such choices exist which correspond to the same ${\hat n}$, implying thereby that the co-ordinates  $(x_p, y_p)$ returned by the procedures of Sec.(\ref{Geom_Det_Rot_Mat}) and Sec.(\ref{SVDPCA_Det_Rot_Mat}) may disagree in general. We indeed observe this to be the case in our computations \cite{Github_OrbM}, but this totally-legitimate variation in projected coordinates is found to be inconsequential in the larger context. With either set of coordinates, the conic section identified through the procedure is found to be robust \cite{Github_OrbM} (e.g. leads to ellipses with identical eccentricities, semi-major, and semi-minor axes), in line with intuitive expectations.

\end{section}

\begin{section}{Intermission}

The above procedures turn raw experimental data into elliptical/conic section-data in  a plane, which is liable for conic section fitting, which is the logical next step. While that may be readily performed, it is natural to wonder exactly how reasonable is this resurrection of Kepler's first law? To conclusively answer this question, one measure that can be adopted is the computation of the Laplace-Runge-Lenz (LRL) vector for the motion of the marble along the elliptic orbit so defined. For Kepler orbits, this LRL vector is a constant of motion \cite{Curtis, Goldstein}, and hence, is conserved in both magnitude and direction. We analyze the same in the subsequent work \cite{OurOrbM_Art2}, and find that the same is conserved within $0.1 \%$, which may be taken as a telltale signature of the Kepler's first law being upheld near-perfectly within this formulation.

Also, since all data points ($x_p, y_p$) are separated by the same time interval $\Delta t$, dictated by the number of frames per second, of the recording camera, velocity components in the plane of the motion can be approximated through the finite difference (either of the forward difference, or the central difference) relations \cite{Weber_Arfken}. In the present work, we work with the central difference convention, and estimate
\begin{equation}
v_{x_p}(i) = \frac{\left(x_p(i+1) - x_p(i-1)\right)}{2 \Delta t}
\end{equation}
It is straightforward to verify that this is exactly how Tracker also estimates the velocity components $v_x$ and $v_y$, from the $x, y$ coordinates, within the software itself. In principle, it is also possible to import these velocity components themselves within Scilab, but then one would have to trace their transformation through all the projections and basis changes performed in Section \ref{Resurrection_Section_KeplerLaw1}, which is significantly arduous relative to the procedure outlined above.

Thus, at this point, we possess both the ${\vec r_p}$ and ${\vec v_p}$ in the plane of motion, computed from the actual motion of the marble on the 3D printed surface. As per the well-known wisdom from classical mechanics \cite{Goldstein, Curtis}, we can then calculate all other dynamical variables from this state vector information.

We explore the applications of this formalism, as well as some generalizations in a subsequent investigation \cite{OurOrbM_Art2}.

\end{section}

\begin{section}*{Acknowledgments}
This work was performed during the 2017, 2018, and 2022 editions of the annual six-week summer undergraduate research camp, \emph{Flavor of Research}, organized by the D. S. Kothari Center for Research and Innovation (DSKC), Miranda House.
The authors
gratefully acknowledge the help from Nitnui Technology,
Delhi, towards 3D printing the eventual $z = 1/r$ surface used seminally throughout this work, and to the in-house DSKC 3D printing facilities for
the rest of 3D printing used in this investigation,
which led to the final outcome.
\end{section}