diff --git a/docs/Smoldyn/SmoldynCodeDoc.pdf b/docs/Smoldyn/SmoldynCodeDoc.pdf index c419d871..dfd8d44e 100644 Binary files a/docs/Smoldyn/SmoldynCodeDoc.pdf and b/docs/Smoldyn/SmoldynCodeDoc.pdf differ diff --git a/docs/Smoldyn/SmoldynCodeDoc.tex b/docs/Smoldyn/SmoldynCodeDoc.tex index 509b531a..fbf654d6 100644 --- a/docs/Smoldyn/SmoldynCodeDoc.tex +++ b/docs/Smoldyn/SmoldynCodeDoc.tex @@ -3970,7 +3970,7 @@ \subsection{Data structures declared in smoldyn.h} \subsubsection{Enumerations} \begin{lstlisting} -enum FilamentDynamics {FDnone, FDeuler}; +enum FilamentDynamics {FDnone, FDeuler,FDRK2,FDRK4}; \end{lstlisting} The filament dynamics enumeration describes what algorithms to use for its dynamics. They are summarized below; this list will be expanded. @@ -3979,7 +3979,9 @@ \subsubsection{Enumerations} Dynamics & Description\\ \hline none & filament is rigid and does not move\\ -euler & use Euler ODE integration\\ +euler & Euler ODE integration\\ +RK2 & 2nd order Runge-Kutta integration \\ +RK4 & 4th order Runge-Kutta integration \end{longtable} \subsubsection{Segments} @@ -3993,8 +3995,10 @@ \subsubsection{Segments} double len; // segment length double thk; // thickness of segment double ypr[3]; // relative ypr angles - double dcm[9]; // relative dcm - double adcm[9]; // absolute segment orientation + double qrel[4]; // relative rotation quaternion + double qabs[4]; // absolute rotation quaternion +// double dcm[9]; // relative dcm +// double adcm[9]; // absolute segment orientation } * segmentptr; \end{lstlisting} @@ -4002,13 +4006,11 @@ \subsubsection{Segments} \ttt{xyzfront} is the $x, y, z$ coordinate of the segment starting point and \ttt{xyzback} is the $x, y, z$ coordinate of the segment end point. These vectors are not owned by the segment but point to vectors in the \ttt{nodes} element of the owning filament. The back of one segment is the identical vector (same memory address) as the front of the next segment. \ttt{len} is the segment length and \ttt{thk} is the segment thickness. -\ttt{ypr} is the relative yaw, pitch, and roll angle for the segment relative to the orientation of the prior segment. \ttt{dcm} is the relative direction cosine matrix for the segment relative to the orientation of the one before it. Identical information is contained in \ttt{ypr} and \ttt{dcm}, but the latter is here for faster computation. \ttt{adcm} is the direction cosine matrix for the absolute orientation; again, this information is redundant but included for faster computation. Restated, \ttt{ypr} and \ttt{dcm} both represent the bend at the front of the segment, and \ttt{adcm} represents the rotation matrix from the absolute coordinate system to this segment's orientation. +\ttt{ypr} is the relative yaw, pitch, and roll angle for the segment relative to the orientation of the prior segment. It is largely ignored for the first segment, whose orientiation is instead given with \ttt{seg0dcm} in the filament data structure. \ttt{qrel} and \ttt{qabs} are rotation quaternions that express the segment's orientation, relative to the prior segment and in the system coordinate system respectively. Although commented out now, \ttt{dcm} and \ttt{adcm} were the direction cosine matrix for the segment relative to the orientation of the one before it, and in the system reference frame, respectively. Identical information is contained in \ttt{ypr}, \ttt{qrel}, and \ttt{qabs} (and, previously, in \ttt{dcm} and \ttt{adcm}), but the redundancy is here for faster computation. -For 2D simulations, only the first 2 values are used in the coordinates vectors (\ttt{xyzfront} and \ttt{xyzback}) and the third coordinate should always equal 0; only the first value (yaw) is used in the ypr angles and the other two should always equal 0; and only the first 2 rows and columns are used in the direction cosine matrices, which still have 9 total elements, while the others equal 0 for unequal indices and 1 for equal indices. +For 2D simulations, only the first 2 values are used in the coordinates vectors (\ttt{xyzfront} and \ttt{xyzback}) and the third coordinate should always equal 0; only the first value (yaw) is used in the ypr angles and the other two should always equal 0; and only the first and last elements are used in the quaternions. -Segment lengths can vary but standard lengths are constant over an entire filament. This means that segments all have roughly equal lengths. I initially thought that it could be useful to have large standard length variations to focus computation on more important regions, but this is sufficiently hard to program that it's not worth it. - -A worthwhile change would be to replace the \ttt{dcm} and \ttt{adcm} matrices with quaternions because that would avoid gimble lock and would enable interpolation and other mathematics. However, the analytical mathematics, presented below, is probably still easier with matrices.\\ +Segment lengths can vary but standard lengths are constant over an entire filament. This means that segments all have roughly equal lengths. I initially thought that it could be useful to have large standard length variations to focus computation on more important regions, but this is sufficiently hard to program that it's not worth it.\\ \subsubsection{Filaments} @@ -4019,9 +4021,15 @@ \subsubsection{Filaments} int maxseg; // number of segments allocated int nseg; // number of segments segmentptr* segments; // array of segments if any - double **nodes; // list of nodes - double **forces; // list of forces on nodes (nseg+1) + double **nodes; // list of nodes (nseg+1) + double **nodes1; // nodes for RK dynamics (nseg+1) + double **nodes2; // nodes for RK4 dynamics (nseg+1) + double **nodesx; // old node locations (nseg+1) + double *roll1; // roll for RK dynamics (nseg) + double *roll2; // roll for RK dynamics (nseg) + double **forces; // forces on nodes (nseg+1) double *torques; // list of segment torques (nseg) + double seg0dcm[9]; // segment 0 dcm struct filamentstruct* frontend; // what front attaches to struct filamentstruct* backend; // what back attaches to int maxbranch; // max branches off this filament @@ -4037,7 +4045,7 @@ \subsubsection{Filaments} Each filament data structure represents a single unbranched chain of segments. Multiple filaments of the same type (e.g. two microtubules) are stored in multiple data structures. In the filament data structure, \ttt{filtype} points to the filament type that owns this filament and \ttt{filname} points to the name of this filament; names are assigned by default if not assigned by the user. -Next, \ttt{maxseg} and \ttt{nseg} give the allocated and actual numbers of segments. These are stored in an array, where the index of the front is equal to 0. Segments are stored in \ttt{segments}. The segment endpoints are stored in \ttt{nodes}, which is allocated to be \ttt{maxseg+1} long by 3 elements wide for the $x$, $y$, and $z$ values. +Next, \ttt{maxseg} and \ttt{nseg} give the allocated and actual numbers of segments. These are stored in an array, where the index of the front is equal to 0. Segments are stored in \ttt{segments}. The segment endpoints are stored in \ttt{nodes}, which is allocated to be \ttt{maxseg+1} long by 3 elements wide for the $x$, $y$, and $z$ values. The \ttt{nodes1} and \ttt{nodes2} elements are similar to nodes, but are used solely for multi-step (Runge-Kutta) integration purposes. The \ttt{nodesx} vector contains node positions at the end of the previous time step, enabling determination of surface crossing or other events during the time step. The \ttt{roll1} and \ttt{roll2} vectors contain the roll values for individual segments and are only used for multi-step integration. The \ttt{forces} element is used to store the cumulative forces on each node. Similarly, the \ttt{torques} element is used to store the torque acting on each segment (along its axis). The \ttt{seg0dcm} element stores the direction cosine matrix for segment 0. Filaments are allowed to branch. If the current filament is a branch off another filament, then \ttt{frontend} points to the filament that the front end attaches to and/or \ttt{backend} points to the filament that the back end attaches to; these are \ttt{NULL} if they are free ends. Additionally, this filament knows about the filaments that branch off this one, which are listed in \ttt{maxbranch}, \ttt{nbranch}, and \ttt{branches}. The branched filaments branch off at the locations listed in \ttt{branchspots}. @@ -4103,11 +4111,14 @@ \subsubsection{Filament superstructure} The superstructure contains a list of filament types. It starts with its condition element in \ttt{condition} and a pointer up the heirarchy to the simulation structure in \ttt{sim}. The list of filament types is allocated with \ttt{maxtype} entries, has \ttt{ntype} entries, and has list \ttt{ftnames} for the type names and \ttt{filtypes} is the actual list of types.\\ +% Filament math - 3D conformations \subsection{Filament math - 3D conformations} -This subsection relates the filament relative angles ($\bm{A}$, \ttt{dcm}), absolute angles ($\bm{B}$, \ttt{adcm}), and positions ($\bm{x}$, \ttt{xyz}). See also Goldstein Chapter 4 (page 128) and Appendix B (page 608), Andrews, 2004 (rotational averaging paper), and Andrews, 2014 (filament paper). Equations are given here for 2D and 3D systems. +There are several ways to describe 3D rotations, including with various conventions of Euler angles, with direction cosine matrices, and with quaternions. Smoldyn uses the yaw-pitch-roll Euler angle convention for relative angles between segments, which it also uses for force and energy computations. I developed most of the filament theory using direction cosine matrices, which are relatively straightforward and widely used in the physics literature, and are described in this subsection. However, I then replaced the direction cosine matrices with quaternions, described below, because they are likely to be faster computationally. At this point, I am continuing to do theory with dcms and code with quaternions. + +This subsection uses direction cosine matrices to relate filament relative angles ($\bm{A}$, was \ttt{dcm} in the code), absolute angles ($\bm{B}$, was \ttt{adcm} in the code), and positions ($\bm{x}$, \ttt{xyz}). See also Goldstein Chapter 4 (page 128) and Appendix B (page 608), Andrews, 2004 (rotational averaging paper), and Andrews, 2014 (filament paper). Equations are given here for 2D and 3D systems. -Filament bending is described by the Tait-Bryan convention: +Filament bending is described by the Tait-Bryan convention for Euler angles: \begin{longtable}[c]{lcccc} name & symbol & perpendicular & positive & positive \\ & & axis & direction & unit vector \\ @@ -4118,7 +4129,7 @@ \subsection{Filament math - 3D conformations} \hspace{3.5ex} pitch & $\theta$ & $y$ & turn down & $\hat{y}$ \\ \hspace{3.5ex} roll & $\psi$ & $x$ & tilt right & $\hat{x}$ \end{longtable} -Note that pitch is defined so that positive pitch is rotation downward. This is the opposite of many conventions, but is used here so that the rotation vector is along the positive $\hat{y}$ vector while using a right-handed coordinate system. In contrast, other conventions use a left-handed coordinate system or other definitions. Figure \ref{fig:CoordinateRotation} shows these rotations using a standard right-handed coordinate system. The pitch rotation angle is negative, with the result that the rotation vector is along the rotated $-\hat{y}$. +Note that pitch is defined so that positive pitch is rotation downward. This is the opposite of many conventions, but is used here so that the rotation vector is along the positive $\hat{y}$ vector while using a right-handed coordinate system. Figure \ref{fig:CoordinateRotation} shows these rotations using a standard right-handed coordinate system. The pitch rotation angle is negative, with the result that the rotation vector is along the rotated $-\hat{y}$. \begin{figure}[h] \centering \includegraphics[height=5cm]{figures/FilamentBend3D.png} @@ -4152,18 +4163,18 @@ \subsection{Filament math - 3D conformations} \label{eq:CoordinateRotationYPR} \end{align} -Right-multiplying the direction cosine matrix with a vector that's expressed in the lab frame returns the coordinates of the same vector but now using the molecule frame. For example, suppose $\bm{r}_L$ is a vector expressed in the lab frame. Then, +Right-multiplying the direction cosine matrix with a vector that's expressed in the lab frame returns the coordinates of the same vector but now using the molecule frame (i.e. this is passive rotation for a molecule frame that is rotated from the lab frame by the given angles). For example, suppose $\bm{r}_L$ is a vector expressed in the lab frame. Then, \begin{equation} \bm{r} = \bm{\Phi} \bm{r}_L \end{equation} gives the same vector but expressed in the molecule frame. As a particularly simple vector (shown with the red arrow in Figure \ref{fig:CoordinateRotation}), consider 2D and suppose $\bm{r}_L=[1,0]^T$. In the molecule frame, which is rotated counter-clockwise by $\phi$ from the lab frame, then $\bm{r} = [c \phi,- s\phi]^T$. -In the code and most of my math, the absolute direction cosine matrix, meaning for conversion between the simulation reference frame and the local segment frame, is called $\bm{B}$. Right multiplying the dcm with a lab-frame vector, here $\bm{p}$, gives the molecule frame, leading to $\bm{p} = \bm{B} \bm{p}_L$. This expression can be transposed and/or multiplied on both sides by $\bm{B}^T$ to yield +In the code and most of my math, the absolute direction cosine matrix, meaning for conversion between the simulation reference frame and the local segment frame, is called $\bm{B}$. Right multiplying the dcm with a lab-frame vector, here $\bm{p}_L$, gives the molecule frame, leading to $\bm{p} = \bm{B} \bm{p}_L$. This expression can be transposed and/or multiplied on both sides by $\bm{B}^T$ to yield (non-transposed vectors are column vectors and transposed ones are row vectors) \begin{align} \bm{p} &= \bm{B} \bm{p}_L \hspace{2cm} \bm{p}_L = \bm{B}^T \bm{p} \\ \bm{p}^T &= \bm{p}_L^T \bm{B}^T \hspace{1.5cm} \bm{p}_L^T = \bm{p}^T \bm{B} \end{align} -The Sphere.c library has functions for these multiplications. The transpose of the dcm listed above performs an active rotation of a vector in a constant coordinate system, in which case it's often called a rotation matrix. +The Sphere.c library has functions for these multiplications. The transpose of the dcm listed above performs an active rotation of a vector in a constant coordinate system. The front end of a filament is the front end of segment 0, which is at $\bm{x}_0$. In its own coordinate system, this segment lies along direction $\hat{x} = [1,0,0]^T$. In the lab coordinate system, this same vector is $\bm{B}_0^T \hat{x}$. From this, the back of segment 0, which is also the front of segment 1, is at $$\bm{x}_1=\bm{x}_0 + l_0 \bm{B}^T_0 \bm{\hat{x}}$$ @@ -4176,42 +4187,125 @@ \subsection{Filament math - 3D conformations} \bm{B}^T_i &= \bm{B}^T_{i-1} \bm{A}^T_i = \bm{A}^T_0 \bm{A}^T_1 \cdots \bm{A}^T_i \label{eq:AbsoluteDCM} \end{align} -The locations of the segment front ends for subsequent segments are +To make sense of this, consider the orientation of the $i$th segment. It is $\hat{\bm{x}}$ in its own reference frame and $\bm{A}_i^T\hat{\bm{x}}$ in the reference frame of the prior segment; continuing this process, it is $\bm{A}^T_0 \bm{A}^T_1 \cdots \bm{A}^T_i\hat{\bm{x}}$, which is $\bm{B}^T_i\hat{\bm{x}}$, in the lab frame. The locations of the segment front ends for subsequent segments are \begin{align} \bm{x}_{i+1} &= \bm{x}_i + l_i \bm{B}^T_i \bm{\hat{x}} \nonumber \\ \Delta \bm{x}_{i} &= \bm{x}_{i+1} - \bm{x}_i = l_i \bm{B}^T_i \bm{\hat{x}} \label{eq:NodePositions} \end{align} -Thus, if we know the ypr angles, we can compute the $\bm{A}$ matrices from eq. \ref{eq:CoordinateRotationYPR}. From those, we can find the $\bm{B}$ matrices from eq. \ref{eq:AbsoluteDCM}. Then, we can find the node positions from eq. \ref{eq:NodePositions}. These computations are the same for 2D and 3D. + +To compute a filament's node positions from angle data, we start with segment 0. This segment's orientation can be set initially with ypr angles, but those are not used when computing dynamics because the yaw and roll angles become undefined when the segment is vertical, and are extremely sensitive to the segment's orientation when it is nearly vertical. Instead, the segment 0 orientation is stored in the \ttt{seg0dcm} direction cosine matrix; this is both $\bm{A}_0$ and $\bm{B}_0$ in the analytical math. For subsequent segments, we know the ypr angles, so we compute the $\bm{A}_i$ matrices from eq. \ref{eq:CoordinateRotationYPR}. From those, we find the $\bm{B}_i$ matrices from eq. \ref{eq:AbsoluteDCM}. Then, we find the node positions from eq. \ref{eq:NodePositions}. These computations are the same for 2D and 3D. Going backward from node positions to ypr angles is more difficult. Starting with 2D systems, the only things that are needed are absolute and relative yaw angles, since pitch and roll are always zero. For segment, $i$, which goes from $\bm{x}_i$ and $\bm{x}_{i+1}$ and has spatial difference vector $\Delta \bm{x}_i=\bm{x}_{i+1} - \bm{x}_{i}$, the segment's absolute angle is $$\phi = \atan2(\Delta x_{i,y},\Delta x_{i,x})$$ This uses the $\atan2$ function, which is the same as the arctangent, but keeps track of the quadrant; this uses $\Delta x_{i,x}$ and $\Delta x_{i,y}$ for the $x$ and $y$ axes. The result is between $-\pi$ and $\pi$, (inclusive at both ends). Note that this fails if both arguments of the function are equal to zero. The code doesn't store the absolute angle, but this is converted to an absolute direction cosine matrix that is stored. The difference between this angle and the angle for the prior segment is the relative angle, but that is no longer necessarily within $[-\pi,\pi]$, so that angle needs to be clamped to this range. Testing shows that this works well, so long as the mobility value isn't so high as to create unstable integration. -For 3D, the best approach does not use absolute angles. Consider segment number $i$, which has displacement $\Delta \bm{x}_i$ and roll $\psi_i$. This displacement is in the lab coordinate system and the roll is in the coordinate system of the prior segment. From \ref{eq:NodePositions}, +The 3D computation does not use absolute angles. Also, the node positions are not sufficient input information but one also needs to know the roll value for each segment, except for segment 0. These roll values give the roll amount relative to the prior segment's orientation. For segment 0, I tried using a roll value relative to the absolute $\hat{\bm{x}}$ vector, but that creates discontinuities when the segment is nearly vertical. Instead, the solution I'm choosing is to record an upward vector, which is in the system reference frame, for segment 0. Ideally, this vector would be perpendicular to the segment and would have unit length, in which case it is the $\hat{\bm{z}}$ vector in the segment 0 reference frame. Deviations from this ideal are not a problem because they can be easily fixed. + +For other segments, consider segment $i$, which has displacement $\Delta \bm{x}_i$ and roll $\psi_i$. This displacement is in the lab coordinate system and the roll is in the coordinate system of the prior segment. From \ref{eq:NodePositions}, $$\Delta \bm{x}_i = l_i \bm{B}^T_i \hat{\bm{x}} = l_i \bm{B}^T_{i-1} \bm{A}^T_i \hat{\bm{x}}$$ Rearrange to \begin{equation} \bm{A}^T_i \hat{\bm{x}} = \frac{1}{l_i} \bm{B}_{i-1} \Delta \bm{x}_i \label{eq:AFromB} \end{equation} -We know the right hand side and want to compute the left hand side. The matrix product on the left expands to local ypr angles, which are in the frame of the $i-1$ segment. Also, define $\Delta \bm{x}'_i$ as the rotated version of $\Delta \bm{x}_i$ into the $i-1$ segment frame, with $\Delta \bm{x}' = \bm{B}_{i-1} \Delta \bm{x}$. These yield +We know the right hand side, and also the roll $\psi_i$, and want to compute the left hand side. The matrix product on the left expands to local ypr angles, which are in the frame of the $i-1$ segment. Define $\Delta \bm{x}'_i$ as the normalized and rotated version of $\Delta \bm{x}_i$ into the $i-1$ segment frame, with $\Delta \bm{x}'_i = \bm{B}_{i-1} \Delta \bm{x}_i/l_i$. These yield $$\bm{A}^T(\phi_i,\theta_i,\psi_i) \hat{\bm{x}} = \left[ \begin{array}{c} \textrm{c} \theta_i \textrm{c} \phi_i \\ \textrm{c} \theta_i \textrm{s} \phi_i \\ -\textrm{s} \theta_i \end{array} \right] = \frac{1}{l_i} \bm{B}_{i-1} \Delta \bm{x}_i -= \frac{1}{l_i} \Delta \bm{x}'_i -= \frac{1}{l_i} \left[ \begin{array}{c} \Delta x'_{i,x} \\ \Delta x'_{i,y} \\ \Delta x'_{i,z} \end{array} \right]$$ += \Delta \bm{x}'_i += \left[ \begin{array}{c} \Delta x'_{i,x} \\ \Delta x'_{i,y} \\ \Delta x'_{i,z} \end{array} \right]$$ Solving the three scalar equations for the two unknowns, $\theta_i$ and $\phi_i$ yield \begin{align} \phi_i &= \atan2(\Delta x'_{i,y}, \Delta x'_{i,x}) \nonumber \\ \theta_i &= -\asin \frac{\Delta x'_{i,z}}{l_i} \label{eq:ThetaPhiFromX} \end{align} -We still have $\psi_i$ from beforehand, still in the coordinate system of the $i-1$ segment. Combine these three values, $\phi_i$, $\theta_i$, and $\psi_i$, with eq. \ref{eq:CoordinateRotationYPR} to get $\bm{A}_i$. Then, +We still have $\psi_i$ from beforehand, still in the coordinate system of the $i-1$ segment. Combine these three values, $\phi_i$, $\theta_i$, and $\psi_i$, with eq. \ref{eq:CoordinateRotationYPR} to get $\bm{A}_i$. This can be done by first computing the ypr angles and substituting them into eq. \ref{eq:CoordinateRotationYPR}. Alternatively, one can substitute in the analytical values from eq. \ref{eq:ThetaPhiFromX} and then simplify. Using the latter approach (and abbreviating $\Delta x'_{i,x}$, $\Delta x'_{i,y}$, and $\Delta x'_{i,z}$ as $x'$, $y'$, and $z'$), +$$\bm{A} = \left[ \begin{array}{ccc} +x' & y' & z' \\ +\frac{-y' \textrm{c} \psi - x'z' \textrm{s} \psi}{\sqrt{x'^2+y'^2}} & +\frac{x' \textrm{c} \psi - y'z' \textrm{s} \psi}{\sqrt{x'^2+y'^2}} & +\sqrt{x'^2+y'^2}\textrm{s} \psi \\ +\frac{y' \textrm{s} \psi - x'z' \textrm{c} \psi}{\sqrt{x'^2+y'^2}} & +\frac{- x' \textrm{s} \psi - y'z' \textrm{c} \psi}{\sqrt{x'^2+y'^2}} & +\sqrt{x'^2+y'^2}\textrm{c} \psi + \end{array} \right]$$ +While we used ypr angles to find this solution, they do not appear here (except for $\psi$) because this is the unique solution to the stated problem of finding $\bm{A}_i$ from $\Delta \bm{x}'_i$ and $\psi_i$, independent of approach. The top row is equal to $x'$, $y'$, and $z'$ because this comes directly from the equation $\bm{A}^T\hat{\bm{x}} = \Delta \bm{x}'_i$. This top row not hypersensitive to particular values of $x'$, $y'$, or $z'$. However, the first two columns of the other two rows are highly sensitive to small changes in $x'$ or $y'$ if the segment is nearly vertical (in the reference frame of the prior segment), so that $x'^2 + y'^2 \approx 0$. This is an unavoidable consequence of how the roll value is defined. + +Next, $$\bm{B}_i = \bm{A}_i \bm{B}_{i-1}$$ This finishes the calculations for segment $i$, and then walk forward to the next segment. + +\subsection{Filament math - quaternions} + +Quaternions work in relatively similar ways (see the online article ``Rotation Quaternions, and How to Use Them'' by Rose, 2015). A quaternion, $\bm{q}$ has 4 values, $(q_0,q_1,q_2,q_3)$. If it represents a 3D point in space, then $\bm{q}=(0,x,y,z)$ and if it represents a rotation, then all four values can be non-zero and they have total magnitude of 1. The quaternion $\bm{q}=(1,0,0,0)$ is the rotation identity quaternion, which produces no rotation. The inverse or complex conjugate quaternion, $\bm{q}^* = (q_0,-q_1,-q_2,-q_3)$, performs the opposite rotation as the original one. + +Quaternions correspond to direction cosine matrices. For quaternion $\bm{q}$, the corresponding direction cosine matrix is as follows (from Rose). Note that the unitary nature of a rotation quaternion ($q_0^2+q_1^2+q_2^2+q_3^2=1$) results in different ways of expressing the diagonal elements. +\begin{equation} +\bm{Q} = \left[ \begin{array}{ccc} +1-2q_2^2-2q_3^2 & 2q_1q_2-2q_0q_3 & 2q_1q_3+2q_0q_2\\ +2q_1q_2+2q_0q_3 & 1-2q_1^2-2q_3^2 & 2q_2q_3-2q_0q_1\\ +2q_1q_3-2q_0q_2 & 2q_2q_3+2q_0q_1 & 1-2q_1^2-2q_2^2 +\end{array} \right] +\label{eq:qtn2dcm} +\end{equation} + +Consider rotation by ypr angles $(\phi,\theta,\psi)$. This converts to quaternion $\bm{q}$ with the following equation (this equation is from Rose but with reversed signs because theirs is for active rotation, whereas this work assumes passive rotation). The equation on the right performs the opposite conversion (again, equation from Rose, but with sign changes) +\begin{equation} +\bm{q} = \left[ \begin{array}{c} +\textrm{c}\frac{\psi}{2} \textrm{c}\frac{\theta}{2} \textrm{c}\frac{\phi}{2} + +\textrm{s}\frac{\psi}{2} \textrm{s}\frac{\theta}{2} \textrm{s}\frac{\phi}{2} \\ +-\textrm{s}\frac{\psi}{2} \textrm{c}\frac{\theta}{2} \textrm{c}\frac{\phi}{2} + +\textrm{c}\frac{\psi}{2} \textrm{s}\frac{\theta}{2} \textrm{s}\frac{\phi}{2} \\ +-\textrm{c}\frac{\psi}{2} \textrm{s}\frac{\theta}{2} \textrm{c}\frac{\phi}{2} - +\textrm{s}\frac{\psi}{2} \textrm{c}\frac{\theta}{2} \textrm{s}\frac{\phi}{2} \\ +-\textrm{c}\frac{\psi}{2} \textrm{c}\frac{\theta}{2} \textrm{s}\frac{\phi}{2} + +\textrm{s}\frac{\psi}{2} \textrm{s}\frac{\theta}{2} \textrm{c}\frac{\phi}{2} +\end{array} \right] +\hspace{1cm} +\bm{ypr}=\left[\begin{array}{c} +\atan2(-2q_0q_3+2q_1q_2,q_0^2+q_1^2-q_2^2-q_3^2) \\ +\asin(-2q_0q_2-2q_1q_3) \\ +\atan2(-2q_0q_1+2q_2q_3,q_0^2-q_1^2-q_2^2+q_3^2) +\end{array} \right] +\label{eq:ypr2qtn} +\end{equation} +To check these equations, I computed from ypr angles $(\phi,\theta,\psi)$ to quaternion and then to dcm, and arrived at eq. \ref{eq:CoordinateRotationYPR}. I also converted from ypr angles to quaternion and back again. + +For passive rotation of position quaternion $\bm{p}$ between the lab frame and molecule frame using rotation quaternion $\bm{q}$, +$$\bm{p} = \bm{q}^* \bm{p}_L \bm{q} \hspace{1cm} +\bm{p}_L = \bm{q} \bm{p} \bm{q}^*$$ +The former equation is for ``forward'' rotation and is equivalent to the matrix equation $\bm{p}=\bm{B}\bm{p}_L$, whereas the latter equation is for ``inverse'' rotation and is equivalent to the matrix equation $\bm{p}_L=\bm{B}^T\bm{p}$. I checked this by converting ypr angles $(\phi,\theta,\psi)$ to a dcm with eq. \ref{eq:CoordinateRotationYPR} and using that to rotate vector $[x,y,z]^T$, and also converting ypr angles to a quaternion with eq. \ref{eq:ypr2qtn} and rotating the same vector; the results were the same. + +For a filament, suppose $\hat{\bm{x}}$ is the quaternion version of the unit vector on the $x$-axis, and $\bm{b}_i$ is the absolute rotation quaternion for segment $i$. Then, the segment displacement is +$$\Delta \bm{x}_i = \bm{x}_{i+1}-\bm{x}_i = l_i \bm{b}_i \hat{\bm{x}} \bm{b}^*_i$$ +From this, it can be seen that quaternions for multiple rotations combine as +$$\bm{b}_i = \bm{b}_{i-1} \bm{a}_i$$ +$$\bm{b}_i^* = \bm{a}_i^* \bm{b}_{i-1}^*$$ +These are opposite the equations for matrices, where $\bm{B}_i=\bm{A}_i\bm{B}_{i-1}$. + +Converting from node positions and roll values to quaternions is similar to the case with direction cosine matrices. The direct approach, which will be shown to fail, starts with the equation +$$\Delta \bm{x}_i = l_i \bm{b}_i \hat{\bm{x}} \bm{b}^*_i = l_i \bm{b}_{i-1} \bm{a}_i \hat{\bm{x}} \bm{a}^*_i \bm{b}^*_{i-1}$$ +Defining $\Delta \bm{x}'_i$ as the vector for the $i$th segment in the reference frame of the $i-1$ segment, and also normalized, this equation rearranges and expands to +$$\Delta \bm{x}'_i +=\left[ \begin{array}{c} \Delta x'_{i,x} \\ \Delta x'_{i,y} \\ \Delta x'_{i,z} \end{array} \right] += \frac{1}{l_i} \bm{b}^*_{i-1} \Delta \bm{x}_i \bm{b}_{i-1} = \bm{a}_i \hat{\bm{x}} \bm{a}^*_i += \left[ \begin{array}{c} 1-2a_2^2-2a_3^2 \\ 2a_1a_2 + 2a_0a_3 \\ 2a_1a_3-2a_0a_2 \end{array} \right]$$ +Note that the vector on the right is the same as the first column of eq. \ref{eq:qtn2dcm}, much as was seen before with dcms. In principle, one could solve this equation for the four $a$ values, and also include the roll value somehow. However, the solution seems difficult (and is beyond Mathematica), and I'm also not sure how to include the roll value. + +Another approach, which does work, is to use the prior solution for the ypr values, from eq. \ref{eq:ThetaPhiFromX} and the previously known roll value. Then substitute these solutions into eq. \ref{eq:ypr2qtn} to compute the quaternion. Doing so is straightforward numerically. It's also easy analytically, although the equations are very long and don't clean up. Thus, one needs to compute the ypr angles $\phi$, $\theta$, and $\psi$ first, and then use those to compute the quaternion. This works well, for the most part, but is problematic when $x$ and $y$ are small, meaning that the filament segment is nearly vertical. + +\subsection{Filament math - segment 0 orientation} + +The segment 0 orientation cannot be described with ypr angles because its orientation values become highly sensitive on the precise orientation when the segment is nearly vertical. Furthermore, the segment orientation also cannot be described as a $\Delta \bm{x}_0$ displacement in combination with a roll value, because the roll angle is still highly sensitive when the segment is nearly vertical. + +The solution is as follows. First, compute all forces. Move nodes in response. Apply rolls in response. However, for segment 0, the roll is not just an angle value. Instead, this segment has a perpendicular vector that defines the local up direction. This vector is rotated for segment 0 rolling. + + + \subsection{Filament math - filament stretching forces} For filament stretching forces, consider a segment, $i$. Its segment vector is $\Delta \bm{x}_i = \bm{x}_{i+1}-\bm{x}_i$, as usual. Also, its length is $l_i = | \bm{x}_i |$. For standard length $l\degree$, the extension force is @@ -4414,6 +4508,14 @@ \subsection{Functions.} \item[\underline{Low level utilities}] +\item[\ttt{int}] +\ttt{filReadFilName(simptr sim,const char *str,filamenttypeptr *filtypeptr,filamentptr *filptr,char *filname)} +\hfill \\ +Reads a filament name in \ttt{str} and returns the filament type in \ttt{filtypeptr}, the filament in \ttt{filptr}, and the filament name in \ttt{filname}. This accepts a name in several possible formats. (1) It allows the format ``filtype:filname'', where a colon separates the type name from the filament name; if the filament type is recognized but not the filament name, then this sends back the filament type and returns -5 (see below for return codes). (2) If the filament type is known, then it can be sent in using \ttt{filtypeptr}, and only that type will be searched for the filament name. (3) If the filament type is not known, then this searches all filament types and, if found, returns the type in \ttt{filtypeptr} and the filament in \ttt{filptr}. + +Note that \ttt{filtypeptr} and \ttt{filptr} cannot be \ttt{NULL} and, the value of \ttt{filtypeptr} must equal \ttt{NULL} if the filament type is unknown or to a filament type if the type is known. Returns 0 for success, -1 for no \ttt{str}, -2 for no filaments have been defined yet, -3 for \ttt{str} is unreadable, -4 for unrecognized filament type name, or -5 for unrecognized filament name. + + \item[\ttt{double}] \ttt{filRandomLength(const filamenttypeptr filtype, double thickness, double sigmamult)} \hfill \\ @@ -4424,12 +4526,14 @@ \subsection{Functions.} \item[\ttt{double}] \ttt{*filRandomAngle(const filamenttypeptr filtype,int dim,int n,double thickness,double sigmamult,double *angle))} \hfill \\ -Local. Returns a random bending angle, which is a relative ypr angle, in \ttt{angle} (which needs to be allocated with 3 elements) and directly. This is for filament type \ttt{filtype}, system dimensionality \ttt{dim} and segment number \ttt{n}; this \ttt{n} number is only checked to see if it's zero, in which case angles are chosen uniformly over all possibilities or non-zero, in which case angles are chosen from angle bending forces. Angles are for a segment with thickness \ttt{thickness} and the computed standard deviation is multiplied by \ttt{sigmamult}. For \ttt{n>0}, if $k_{ypr}$ is less than 0 for a particular coordinate, this returns the standard bending angle, if it is equal to 0, this returns a uniformly distributed bending angle, and if it is greater than 0, this returns a mean of $ypr_{std}$ and standard deviation of +Local. Returns a random bending angle, which is a relative ypr angle, in \ttt{angle} (which needs to be allocated with 3 elements) and directly. This is for filament type \ttt{filtype}, system dimensionality \ttt{dim} and segment number \ttt{n}; this \ttt{n} number is only checked to see if it's zero or non-zero; in the latter case, angles are chosen from angle bending forces. Angles are for a segment with thickness \ttt{thickness} and the computed standard deviation is multiplied by \ttt{sigmamult}. For \ttt{n>0}, if $k_{ypr}$ is less than 0 for a particular coordinate, this returns the standard bending angle, if it is equal to 0, this returns a uniformly distributed bending angle, and if it is greater than 0, this returns a mean of $ypr_{std}$ and standard deviation of $$\sigma_{mult} \sqrt{\frac{kT}{R*k_{ypr}}}$$ For 2D space, this sets the last two elements of \ttt{angle} to 0. +For segment 0, angles are chosen uniformly over all possibilities. For 2D, the returned angle, which is in the first element only (yaw) is relative to the system $\hat{\bm{x}}$ vector. + \item[\underline{Computations on filaments}] \item[\ttt{double filStretchEnergy(const filamentptr fil, int seg1, int seg2)}] @@ -4448,7 +4552,7 @@ \subsection{Functions.} \item[\ttt{void filBendTorque(const filamentptr fil,int node,double *torque)}] \hfill \\ -Computes the bending torque vector for node number \ttt{node}, which is the front of segment \ttt{node}, in filament \ttt{fil}. The \ttt{torque} vector needs to be pre-allocated. This uses eq. \ref{eq:BendTorque} and then rotates the torque with the absolute direction cosine matrix of the prior segment to express the torque in absolute coordinates, which is how it's returned. +Computes the bending torque vector for node number \ttt{node}, which is the front of segment \ttt{node}, in filament \ttt{fil}. The \ttt{torque} vector needs to be pre-allocated. This uses eq. \ref{eq:BendTorque} and then rotates the torque with the absolute quaternion of the prior segment to express the torque in absolute coordinates, which is how it's returned ($\bm{\tau}_L = \bm{B}^T\bm{\tau}$). \item[\underline{Memory management}] @@ -4529,15 +4633,15 @@ \subsection{Functions.} $\bm{x}_1=\bm{x}_0 + l_0 \bm{B}^T_0 \cdot \bm{\hat{x}}$ For segments added to the back, with index $i$: -$\bm{a}_i = \ttt{angle}$, -$\bm{A}_i = DCM(\bm{a}_i)$, +$\bm{ypr}_i = \ttt{angle}$, +$\bm{A}_i = DCM(\bm{ypr}_i)$, $\bm{B}_i = \bm{A}_i \cdot \bm{B}_{i-1}$, $\bm{x}_{i+1} = \bm{x}_i + l_i \bm{B}^T_i \cdot \bm{\hat{x}}$ For segments added to the front, with index $i$ (typically equal to 0): $\bm{B}_i = \bm{A}_{i+1}^T \cdot \bm{B}_{i+1}$, -$\bm{A}_i = \bm{B}_i$ -$\bm{a}_i = XYZ(\bm{B}_i)$, +$\bm{A}_i = \bm{B}_i$, +$\bm{ypr}_i = YPR(\bm{B}_i)$, $\bm{x}_i = \bm{x}_{i+1} - l_i \bm{B}^T_i \cdot \bm{\hat{x}}$ \item[\ttt{int}] @@ -4551,26 +4655,35 @@ \subsection{Functions.} \hfill \\ Removes one segment or one bead from either the front or back end of a filament. Specify the end in \ttt{endchar} with `f' for front and `b' for back. Returns 0 for success or -1 for a filament that already had no segments. +Removing a segment from the back requires only shortening the filament by one segment. Removing one from the front requires copying the absolute dcm or quaternion into the relative dcm or quaternion for the new front segment, and then translating this into the ypr angles. + \item[\ttt{void filTranslate(filamentptr fil, const double *vect, char func)}] \hfill \\ Translates an entire filament in space without rescaling or rotation. Enter \ttt{vect} with a 3D vector and \ttt{func} with `=' for translate to the given posiition, with `-' for to subtract the value of \ttt{vect} from the current position, and with `+' to add the value of \ttt{vect} to the current position. For 2D simulations, \ttt{vect[2]} needs to equal 0. -\item[\ttt{void}] +\item[\ttt{int}] \ttt{filLengthenSegment(filamentptr fil, int seg, double length, char endchar, char func)} \hfill \\ -This function modifies the length of a single segment of index \ttt{seg} in filament \ttt{fil}. The \ttt{endchar} input can be \ttt{b} or \ttt{f} and gives which end of the segment moves, along with the entire rest of that end of the filament. The \ttt{func} input can be \ttt{=}, \ttt{+}, or \ttt{-} and gives whether the segment length should be made equal to \ttt{length}, increased by \ttt{length}, or decreased by \ttt{length}. This function is designed to run fairly fast, so doesn't loop over dimensions and it doesn't call functions from elsewhere. +This function modifies the length of a single segment of index \ttt{seg} in filament \ttt{fil}. The \ttt{endchar} input can be \ttt{b} or \ttt{f} and gives which end of the segment moves, along with the entire rest of that end of the filament. The \ttt{func} input can be \ttt{=}, \ttt{+}, or \ttt{-} and gives whether the segment length should be made equal to \ttt{length}, increased by \ttt{length}, or decreased by \ttt{length}. Returns 0 for success or an error code of 1 if the new segment length would be 0 or negative. + +This function computes the change of length, and then multiplies that by the local segment orientation to create the shift vector \ttt{xdelta}. This is added to the following or preceding nodes, depending on \ttt{endchar}. + +\item[\ttt{int}] +\ttt{filChangeThickness(filamentptr fil, int seg, double thick, char func)} +\hfill \\ +This function modifies the thickness of a single segment of index \ttt{seg} in filament \ttt{fil}. The \ttt{func} input can be \ttt{=}, \ttt{+}, or \ttt{-} and gives whether the segment thickness should be made equal to \ttt{thick}, increased by \ttt{thick}, or decreased by \ttt{thick}. Returns 0 for success or an error code of 1 if the new segment thickness would be negative. \item[\ttt{void}] \ttt{filRotateVertex(filamentptr fil, int seg, const double *angle, char endchar, char func)} \hfill \\ This function rotates part of the filament about the vertex that is at the front of segment \ttt{seg} by ypr angle \ttt{angle}. Enter \ttt{endchar} as \ttt{b} to have the back of the filament move and as \ttt{f} to have the front move. Enter \ttt{func} as \ttt{=} to set the current relative rotation angles (\ttt{segment->dcm}) to the input ones, as \ttt{+} for increased rotation by the given angles and as \ttt{-} for decreased rotation by these angles. -This function makes the requested modification to the \ttt{dcm} angles and then propagates the changes either backward or forward as needed. During propagation, new segment positions are computed from neighbor positions, rather than being rotated independently, in order to avoid round-off errors. +This function makes the requested modification to the segment's relative quaternion. Then, it propagates the changes either backward or forward as needed. During propagation, new segment positions are computed from neighbor positions, rather than being rotated independently, in order to avoid round-off errors. \item[\ttt{int}] -\ttt{filCopyFilament(filamentptr filto, const filamentptr filfrom)} +\ttt{filCopyFilament(filamentptr filto, const filamentptr filfrom,const filamenttypeptr filtype)} \hfill \\ -This function copies all of the values from \ttt{filfrom} to \ttt{filto} except for the name, completely overwriting any prior data in \ttt{filto} in the process. Any elements in bead, segment, and monomer lists are shifted back so that they start at index number 0. Returns 0 for success, 1 for out of memory, or 2 for illegal function inputs. +This function copies nearly all of the values from \ttt{filfrom} to \ttt{filto}, overwriting any prior data in \ttt{filto} in the process. It does not copy the name. Also, it assigns the filament type to \ttt{filtype} if this value is not \ttt{NULL}, and to the \ttt{filfrom} filament type otherwise. Returns 0 for success, 1 for out of memory, or 2 for illegal function inputs. \item[\ttt{filamentptr}] \ttt{filAddFilament(filamenttypeptr filtype, const char *filname)} diff --git a/docs/Smoldyn/SmoldynManual.pdf b/docs/Smoldyn/SmoldynManual.pdf index 4be853f4..33fa8f82 100644 Binary files a/docs/Smoldyn/SmoldynManual.pdf and b/docs/Smoldyn/SmoldynManual.pdf differ diff --git a/docs/Smoldyn/SmoldynManual.tex b/docs/Smoldyn/SmoldynManual.tex index 1d197908..a5ca4d2f 100644 --- a/docs/Smoldyn/SmoldynManual.tex +++ b/docs/Smoldyn/SmoldynManual.tex @@ -3585,69 +3585,65 @@ \section{Statements about filaments} Start of filament type definition block. The filament type name may be given with $name$, or it may be given afterward with the \ttt{name} statement. If the name has not been used yet for a filament type, then a new filament type is started. Between this instruction and \ttt{end\_filament\_type}, all lines need to pertain to filament types. Parameters of one filament type can be listed in multiple blocks, or parameters for many filament types can be listed in one block. -\item{\ttt{name} $name$} +\item{\ttt{* name} $name$} From within a filament type definition block, this switches to a new filament type and creates it if needed. -\item{\ttt{dynamics} $dynamics$} +\item{\ttt{* dynamics} $dynamics$} Sets the dynamics for the filament type. Current options are: ``none'', ``Euler''. These are case-insensitive. -\item{\ttt{biology} $biology$} - -Sets the biology value for the filament type. Options are: ``actin'', ``microtubule'', ``intermediate'', ``dsDNA'', ``ssDNA'', and ``other''. These are case-insensitive. - -\item{\ttt{color} $color$} +\item{\ttt{* color} $color$} Sets the drawing color for the filament type. Enter the color as either a color description as a word or as RGB$\alpha$ values. -\item{\ttt{thickness} $thickness$} +\item{\ttt{* thickness} $thickness$} Sets the drawing thickness for the filament type. Enter this in pixels if the graphics type is ``opengl'' or as an actual length value for the better drawing options. -\item{\ttt{stipple} $factor$ $pattern$} +\item{\ttt{* stipple} $factor$ $pattern$} Sets the stipple pattern and factor for filament drawing. This is only relevant for simulations with ``opengl\_good'' or better display method. In $factor$, which is an integer, enter the repeat distance for the entire stippling pattern (1 is a good choice). In $pattern$, which is a hexadecimal integer, enter the stippling pattern between 0x0000 and 0xFFFF. 0x00FF has long dashes, 0x0F0F has medium dashes, 0x5555 has dots, etc. Turn stippling off with 0xFFFF. -\item{\ttt{polygon} $drawmode$} +\item{\ttt{* polygon} $drawmode$} Drawing mode for filament surfaces. Options are: ``vertex'', ``edge'', ``face'', or combinations of these by using single letters for vertex, edge, and face, such as ``ve'', ``vef'', etc. -\item{\ttt{shininess} $value$} +\item{\ttt{* shininess} $value$} Shininess of the filament for drawing purposes. This value can range from 0 for visually flat surfaces to 128 for very shiny surfaces. This is only relevant for some simulations. -\item{\ttt{kT} $value$} +\item{\ttt{* kT} $value$} Set the thermal energy value, equal to Boltzmann's constant times temperature. Enter this in energy units that are consistent with other units in the input file. This is only used for some of the dynamics options. -\item{\ttt{treadmill\_rate} $value$} +\item{\ttt{* treadmill\_rate} $value$} Some filaments treadmill, in which they drop off monomers at one end and add them to the other end. Set the treadmilling rate here. -\item{\ttt{mobility} $value$} +\item{\ttt{* mobility} $value$} Mobility of the nodes within the surrounding medium. -\item{\ttt{standard\_length} $length$} +\item{\ttt{* standard\_length} $length$} Relaxed length of a filament segment. It can change through stretching or compression. -\item{\ttt{standard\_angle} $yaw$\\ -\ttt{standard\_angle} $yaw$ $pitch$ $roll$} +\item{\ttt{* standard\_angle} $yaw$\\ +\ttt{* standard\_angle} $yaw$ $pitch$ $roll$} Relaxed angles between adjacent filament segments. When facing toward the filament's front end, yaw represents left-right bending, pitch represents up-down bending, and roll represents rotation about the filament axis. For 2D simulations, only enter a single bending angle. -\item{\ttt{force\_length} $value$} +\item{\ttt{* force\_length} $value$} Stretching force constant for filament segments. -\item{\ttt{force\_angle} $yaw$\\ -\ttt{force\_angle} $yaw$ $pitch$ $roll$} +\item{\ttt{* force\_angle} $yaw$\\ +\ttt{* force\_angle} $yaw$ $pitch$ $roll$} Force constants for bending and torsion. -\item{\ttt{end\_filament\_type}} +\item{\ttt{* end\_filament\_type}} End of a block of filament type definitions. Filament type statements are no longer recognized but other simulation statements are. @@ -3657,36 +3653,42 @@ \section{Statements about filaments} \begin{description} -\item{\ttt{random\_filament} $name$ $type$ $segments$ $[x$ $y$ $z$ $\theta$ $\phi$ $\chi]$ $[thickness]$} +\item{\ttt{random\_filament} $type$:$name$ $segments$ $[x$ $y$ $\phi]$ $[thickness]$}\\ +\ttt{random\_filament} $type$:$name$ $segments$ $[x$ $y$ $z$ $\phi$ $\theta$ $\psi]$ $[thickness]$ -Create a new filament with random segments. It is named $name$, is of type $type$, has $segments$ number of segments, and, optionally, has its starting location at $(x,y,z)$ and initial yaw-pitch-roll angle $(\theta, \phi, \chi)$. It also has optional thickness $thickness$. +Create a new filament with random segments. It is named $name$, is of type $type$, has $segments$ number of segments, and, optionally, has its starting location at $(x,y,z)$ and initial yaw-pitch-roll angle $(\phi, \theta, \psi)$. For 2D, only enter $x,y$ and $\phi$. It also has optional thickness $thickness$. -\item{\ttt{start\_filament} $type$ $name$} +\item{\ttt{start\_filament} $type$:$name$} Start of filament definition block, for a filament of type $type$ and named $name$. If the name has not been used yet for a filament of this type, then a new filament is started. Between this instruction and \ttt{end\_filament}, all lines need to pertain to filaments. Parameters of one filament can be listed in multiple blocks, or parameters for many filaments can be listed in one block. -\item{\ttt{new\_filament} $type$ $name$} +\item{\ttt{new\_filament} $type$:$name$} Defines a new filament of type $type$ that is called $name$, but does not start a filament block. This statement is largely redundant with \ttt{start\_filament}. -\item{\ttt{* name} $type$ $name$} +\item{\ttt{* name} $type$:$name$} \\ +\ttt{* name} $name$ Name of the filament for editing. This statement is not required because the filament name can also be given with \ttt{start\_filament}. This statement gives the name of the current filament for editing, and creates a new filament if needed. -\item{\ttt{* first\_segment} $x$ $y$ $length$ $angle_0$ $[thickness]$\\ -\ttt{* first\_segment} $x$ $y$ $z$ $length$ $angle_0$ $angle_1$ $angle_2$ $[thickness]$} +\item{\ttt{* first\_segment} $x$ $y$ $length$ $angle$ $[thickness]$\\ +\ttt{* first\_segment} $x$ $y$ $z$ $length$ $yaw$ $pitch$ $roll$ $[thickness]$} -Defines the first segment of a filament. Enter the starting location in $x$, $y$, and $z$, the segment length in $length$, and its absolute orientation with the angle values using spherical coordinates. The thickness value is optional. +Defines the first segment of a filament. Enter the starting location in $x$, $y$, and $z$, the segment length in $length$, and its absolute orientation with the angle values. For 2D, the angle is the angle from the $x$-axis and in 3D these values are the yaw-pitch-roll values away from the $x$-axis. The thickness value is optional. -\item{\ttt{* add\_segment} $length$ $angle_0$ $[thickness\ [end]]$\\ -\ttt{* add\_segment} $length$ $angle_0$ $angle_1$ $angle_2$ $[thickness\ [end]]$} +\item{\ttt{* add\_segment} $length$ $angle$ $[thickness\ [end]]$\\ +\ttt{* add\_segment} $length$ $yaw$ $pitch$ $roll$ $[thickness\ [end]]$} -Add a segment to the filament, with the given length. The angles are relative angles using yaw-pitch-roll values, and the thickness is optional. The segment can be added to either end, given in $end$, where the options are $front$ or $back$. +Add a segment to the filament, with the given length. The angles are relative angles from the prior segment using yaw-pitch-roll values, and the thickness is optional. The segment can be added to either end, given in $end$, where the options are $front$ or $back$. \item{\ttt{* remove\_segment} $end$} Remove one segment from end $end$ of the filament (``front'' or ``back''). +\item{\ttt{* modify\_segment} $segment$ length/angle/thickness +/-/= $value$ $end$} + +This modifies the length, relative angle, or thickness parameters for the single segment number $segment$. After the segment number, enter one of ``length'', ``angle'', or ``thickness''. Then, enter ``+'' to increase the current value by the given amount, ``-'' to decrease the current value by the given amount, or ``='' to set the value to the given amount, where the given amount is entered in $value$. For and angle option with 3D, enter 3 values for yaw, pitch, and roll. Finally, for the length and angle options, enter $end$ as either ``front'' or ``back'' for which filament end is moved. + \item{\ttt{* random\_segments} $number$ $[x$ $y$ $z]$ $[thickness]$} Add $number$ of random segments to the filament. @@ -3695,9 +3697,9 @@ \section{Statements about filaments} Translate the filament. Set $symbol$ to ``='' for absolute coordinates, ``+'' for relative translation in which the values are added to the current filament location, and ``-'' for relative translation in the other direction. -\item{\ttt{* copy} $filament$} +\item{\ttt{* copy\_to} $filament$} -Copy monomers to the current filament from $filament$. +Copy the current filament to $filament$, which is created here if it doesn't already exist; any prior data in $filament$ is overwritten. \item{\ttt{* end\_filament}} diff --git a/docs/libSteve/Rn_doc.doc b/docs/libSteve/Rn_doc.doc index 3ede92f0..205e16e0 100644 Binary files a/docs/libSteve/Rn_doc.doc and b/docs/libSteve/Rn_doc.doc differ diff --git a/docs/libSteve/Sphere.pdf b/docs/libSteve/Sphere.pdf index c13a9a7b..46661dc5 100644 Binary files a/docs/libSteve/Sphere.pdf and b/docs/libSteve/Sphere.pdf differ diff --git a/docs/libSteve/Sphere.tex b/docs/libSteve/Sphere.tex index f7b56e2a..242a9852 100644 --- a/docs/libSteve/Sphere.tex +++ b/docs/libSteve/Sphere.tex @@ -43,6 +43,7 @@ \lstset{style=mystyle} \newcommand {\ttt} {\texttt} +\DeclareMathOperator{\atan2} {atan2} \begin{document} \maketitle @@ -52,7 +53,7 @@ \section{Header file: Sphere.h} \begin{lstlisting}[language=C] /* Steven Andrews 2/05. See documentation called Sphere_doc.doc. -Copyright 2005-2007 by Steven Andrews. This work is distributed under the terms +Copyright 2005-2007 by Steven Andrews. This work is distributed under the terms of the Gnu Lesser General Public License (LGPL). */ @@ -65,38 +66,56 @@ \section{Header file: Sphere.h} Dcm = Direction cosine matrix (A00, A01, A02, A10, A11, A12, A20, A21, A22) Eax = Euler angles x-convention (theta, phi, psi) Eay = Euler angles y-convention (theta, phi, chi) -Ep = Euler parameters (e0, e1, e2, e3) -Xyz = xyz- or ypr- or Tait-Bryan convention (yaw, pitch, roll): rotate on z, then y, then x +Qtn = Quaternions (q0, q1, q2, q3) +Ypr = xyz- or ypr- or Tait-Bryan convention (yaw, pitch, roll): rotate on z, then y, then x */ -void Sph\_Cart2Sc(const double *Cart,double *Sc); -void Sph\_Sc2Cart(const double *Sc,double *Cart); -void Sph\_Eay2Ep(const double *Eay,double *Ep); -void Sph\_Xyz2Xyz(const double *Xyz1,double *Xyz2); -void Sph\_Eax2Xyz(const double *Eax,double *Xyz); - -void Sph\_Eax2Dcm(const double *Eax,double *Dcm); -void Sph\_Eay2Dcm(const double *Eay,double *Dcm); -void Sph\_Xyz2Dcm(const double *Xyz,double *Dcm); -void Sph\_Xyz2Dcmt(const double *Xyz,double *Dcmt); -void Sph\_Dcm2Xyz(const double *Dcm,double *Xyz); -void Sph\_Dcm2Dcm(const double *Dcm1,double *Dcm2); -void Sph\_Dcm2Dcmt(const double *Dcm1,double *Dcm2); - -void Sph\_DcmxDcm(const double *Dcm1,const double *Dcm2,double *Dcm3); -void Sph\_DcmxDcmt(const double *Dcm1,const double *Dcmt,double *Dcm3); -void Sph\_DcmtxDcm(const double *Dcmt,const double *Dcm2,double *Dcm3); -void Sph\_DcmxCart(const double *Dcm,const double *Cart,double *Cart2); - -void Sph\_One2Dcm(double *Dcm); -void Sph\_Xyz2Xyzr(const double *Xyz,double *Xyzr); -void Sph\_Dcm2Dcmr(const double *Dcm,double *Dcmr); -void Sph\_Rot2Dcm(char axis,double angle,double *Dcm); -void Sph\_Newz2Dcm(const double *Newz,double psi,double *Dcm); - -void Sph\_DcmtxUnit(const double *Dcmt,char unit,double *vect,const double *add,double mult); - -double Sph\_RotateVectWithNormals3D(const double *pt1,const double *pt2,double *newpt2,double *oldnorm,double *newnorm,int sign); +void Sph_Cart2Sc(const double *Cart,double *Sc); +void Sph_Sc2Cart(const double *Sc,double *Cart); +void Sph_Cart2Cart(const double *Cart1,double *Cart2); +void Sph_Ypr2Ypr(const double *Ypr1,double *Ypr2); +void Sph_Eax2Ypr(const double *Eax,double *Ypr); + +void Sph_Eax2Dcm(const double *Eax,double *Dcm); +void Sph_Eay2Dcm(const double *Eay,double *Dcm); +void Sph_Ypr2Dcm(const double *Ypr,double *Dcm); +void Sph_Ypr2Dcmt(const double *Ypr,double *Dcmt); +void Sph_Dcm2Ypr(const double *Dcm,double *Ypr); +void Sph_Dcm2Dcm(const double *Dcm1,double *Dcm2); +void Sph_Dcm2Dcmt(const double *Dcm1,double *Dcm2); + +void Sph_DcmxDcm(const double *Dcm1,const double *Dcm2,double *Dcm3); +void Sph_DcmxDcmt(const double *Dcm1,const double *Dcmt,double *Dcm3); +void Sph_DcmtxDcm(const double *Dcmt,const double *Dcm2,double *Dcm3); +void Sph_DcmxCart(const double *Dcm,const double *Cart,double *Cart2); +void Sph_DcmtxCart(const double *Dcm,const double *Cart,double *Cart2); + +void Sph_One2Dcm(double *Dcm); +void Sph_Ypr2Yprr(const double *Ypr,double *Yprr); +void Sph_Dcm2Dcmr(const double *Dcm,double *Dcmr); +void Sph_Rot2Dcm(char axis,double angle,double *Dcm); +void Sph_Newz2Dcm(const double *Newz,double psi,double *Dcm); + +void Sph_DcmtxUnit(const double *Dcmt,char unit,double *vect,const double *add,double mult); + +void Sph_One2Qtn(double *Qtn); +void Sph_Qtn2Qtn(const double *Qtn1,double *Qtn2); +void Sph_Ypr2Qtn(const double *Ypr,double *Qtn); +void Sph_Ypr2Qtni(const double *Ypr,double *Qtni); +void Sph_Qtn2Ypr(const double *Qtn,double *Ypr); +void Sph_Dcm2Qtn(const double *Dcm,double *Qtn); +void Sph_Qtn2Dcm(const double *Qtn,double *Dcm); +void Sph_XZ2Qtni(const double *x,const double *z,double *Qtn); +void Sph_QtnxQtn(const double *Qtn1,const double *Qtn2,double *Qtn3); +void Sph_QtnixQtn(const double *Qtn1,const double *Qtn2,double *Qtn3); +void Sph_QtnxQtni(const double *Qtn1,const double *Qtn2,double *Qtn3); +void Sph_QtnRotate(const double *Qtn,const double *Cart,double *Cart2); +void Sph_QtniRotate(const double *Qtn,const double *Cart,double *Cart2); +void Sph_QtniRotateUnitx(const double *Qtni,double *vect,const double *add,double mult); +void Sph_QtniRotateUnitz(const double *Qtni,double *vect,const double *add,double mult); + +double Sph_RotateVectWithNormals3D(const double *pt1,const double *pt2,double *newpt2,double *oldnorm,double *newnorm,int sign); +void Sph_RotateVectorAxisAngle(const double *vect,const double *axis,double angle,double *rotated); #endif \end{lstlisting} @@ -111,36 +130,36 @@ \section{Description} \begin{description} \item[Cartesian coordinates (Cart)] -\hfill \\ -Vector is $[x,y,z]$, all of which are on $(-\infty,\infty)$. - -\item[Spherical coordinates (Sc)] -\hfill \\ -Vector is $[r,\theta,\phi]$. $r$ is on $[0,\infty)$, $\theta$ is on $[0,\pi]$, and $\phi$ is on $(-\pi,\pi]$. - -\item[Direction cosine matrix (Dcm)] -\hfill \\ -Matrix is given as a 9 element array, which lists the matrix row by row. This is useful for all coordinate transformations and is not associated with any particular convention. - -\item[Direction cosine matrix transpose (Dcmt)] -\hfill \\ -This is entered as a normal, non-transposed, direction cosine matrix. However, it is interpreted as a transposed direction cosine matrix in the code. - -\item[Euler angle $x$-convention (Eax)] -\hfill \\ -Vector is $[\theta,\phi,\psi]$. $\phi$ is on $[0,2\pi)$, $\theta$ is on $[0,\pi]$, $\psi$ is on $[0,2\pi)$. $\bm{A} = \bm{Z}(\psi)\bm{X}(\theta)\bm{Z}(\phi)$. - -\item[Euler angle $y$-convention (Eay)] -\hfill \\ -Vector is $[\theta,\phi,\chi]$. $\phi$ is on $[0,2\pi)$, $\theta$ is on $[0,\pi]$, $\chi$ is on $[0,2\pi)$. $\bm{A} = \bm{Z}(\chi)\bm{Y}(\theta)\bm{Z}(\phi)$. - -\item[Euler parameters (Eap)] -\hfill \\ -Vector is $[e_0,e_1,e_2,e_3]$. - -\item[Yaw-pitch-roll (Xyz or Ypr)] -\hfill \\ -Vector is $[\phi,\theta,\psi]$. All are on $(-\pi,\pi]$. These angles are sometimes not clamped because multiple rotations are possible. $\bm{A} = \bm{X}(\psi)\bm{Y}(\theta)\bm{Z}(\phi)$. +\hfill \\ +Vector is $[x,y,z]$, all of which are on $(-\infty,\infty)$. + +\item[Spherical coordinates (Sc)] +\hfill \\ +Vector is $[r,\theta,\phi]$. $r$ is on $[0,\infty)$, $\theta$ is on $[0,\pi]$, and $\phi$ is on $(-\pi,\pi]$. + +\item[Direction cosine matrix (Dcm)] +\hfill \\ +Matrix is given as a 9 element array, which lists the matrix row by row. This is useful for all coordinate transformations and is not associated with any particular convention. + +\item[Direction cosine matrix transpose (Dcmt)] +\hfill \\ +This is entered as a normal, non-transposed, direction cosine matrix. However, it is interpreted as a transposed direction cosine matrix in the code. + +\item[Euler angle $x$-convention (Eax)] +\hfill \\ +Vector is $[\theta,\phi,\psi]$. $\phi$ is on $[0,2\pi)$, $\theta$ is on $[0,\pi]$, $\psi$ is on $[0,2\pi)$. $\bm{A} = \bm{Z}(\psi)\bm{X}(\theta)\bm{Z}(\phi)$. + +\item[Euler angle $y$-convention (Eay)] +\hfill \\ +Vector is $[\theta,\phi,\chi]$. $\phi$ is on $[0,2\pi)$, $\theta$ is on $[0,\pi]$, $\chi$ is on $[0,2\pi)$. $\bm{A} = \bm{Z}(\chi)\bm{Y}(\theta)\bm{Z}(\phi)$. + +\item[Quaternions (Qtn)] +\hfill \\ +Vector is $[q_0,q_1,q_2,q_3]$. + +\item[Yaw-pitch-roll (Xyz or Ypr)] +\hfill \\ +Vector is $[\phi,\theta,\psi]$. All are on $(-\pi,\pi]$. These angles are sometimes not clamped because multiple rotations are possible. $\bm{A} = \bm{X}(\psi)\bm{Y}(\theta)\bm{Z}(\phi)$. \end{description} @@ -163,6 +182,7 @@ \section{History} \item[3/13/24] Added \ttt{Sph\_Eax2Xyz}. \item[3/15/24] Added \ttt{Sph\_DcmxCart}. \item[3/16/24] Rewrote docs in LaTeX. +\item[10/5/24] Renamed Xyz to Ypr. Added math text. Removed Euler parameters. Added quaternions. \end{description} @@ -170,12 +190,48 @@ \section{History} % Math \section{Math} -This library performs passive rotations, meaning that the coordinate system rotates while the vector stays fixed. This is the standard textbook method, such as in Goldstein. Elementary rotation matrices are: -\begin{align*} -\bm{X}(\psi) &= \left[ \begin{array}{ccc} 1 & 0 & 0 \\ 0 & \textrm{c} \psi & \textrm{s} \psi \\ 0 & -\textrm{s} \psi & \textrm{c} \psi \end{array} \right]\\ -\bm{Y}(\theta) &= \left[ \begin{array}{ccc} \textrm{c} \theta & 0 & -\textrm{s} \theta \\ 0 & 1 & 0 \\ \textrm{s} \theta & 0 & \textrm{c} \theta \end{array} \right] \\ -\bm{Z}(\phi) &= \left[ \begin{array}{ccc} \textrm{c} \phi & \textrm{s} \phi & 0 \\ -\textrm{s} \phi & \textrm{c} \phi & 0 \\ 0 & 0 & 1 \end{array} \right] -\end{align*} +The following text was copied from SmoldynCodeDoc, and then edited. See also Goldstein Chapter 4 (page 128) and Appendix B (page 608), Andrews, 2004 (rotational averaging paper), and Andrews, 2014 (filament paper). + +Yaw-pitch-roll angles are described by the Tait-Bryan convention: +\begin{longtable}[c]{lcccc} +name & symbol & perpendicular & positive & positive \\ + & & axis & direction & unit vector \\ +\hline +yaw & $\phi$ & $z$ & turn left & $\hat{z}$ \\ +pitch & $\theta$ & $y$ & turn down & $\hat{y}$ \\ +roll & $\psi$ & $x$ & tilt right & $\hat{x}$ +\end{longtable} +Note that pitch is defined so that positive pitch is rotation downward. This is the opposite of many conventions, but is used here so that the rotation vector is along the positive $\hat{y}$ vector while using a right-handed coordinate system. + +Rotation can be given with a direction cosine matrix (dcm), $\bm{\Phi}$, which can be expanded as +\begin{align} +\bm{\Phi} &= \left[ \begin{array}{ccc} \Phi_{Xx} & \Phi_{Xy} & \Phi_{Xz} \\ \Phi_{Yx} & \Phi_{Yy} & \Phi_{Yz} \\ \Phi_{Zx} & \Phi_{Zy} & \Phi_{Zz} \end{array} \right] +\label{eq:DirectionCosineMatrix} +\end{align} +This dcm can convert between two frames of reference for a static object, called passive rotation, or can rotate an object in a fixed reference frame, called active rotation. Here, $X$, $Y$, and $Z$ are the unit vectors of the ``lab frame'', meaning the absolute coordinates of the simulation, and $x$, $y$, and $z$ are the unit vectors of the ``molecule frame''. The dcm expresses the dot products of these unit vectors, so each column gives the lab frame coordinates of each unit vector of the molecule frame, and each row gives the molecule frame coordinates of each unit vector of the lab frame. It is unitary, meaning that all its eigenvalues equal 1, and its transpose is its inverse, so $\bm{\Phi}^T \bm{\Phi} = \bm{1}$, where $\bm{1}$ is the identity matrix. + +The dcm for a product of yaw, pitch, and roll transformations is (Goldstein pp. 135 and 609) +\begin{align} +\bm{\Phi} &= \bm{X}(\psi) \bm{Y}(\theta) \bm{Z}(\phi) \nonumber \\ +&= \left[ \begin{array}{ccc} 1 & 0 & 0 \\ 0 & \textrm{c} \psi & \textrm{s} \psi \\ 0 & -\textrm{s} \psi & \textrm{c} \psi \end{array} \right] +\left[ \begin{array}{ccc} \textrm{c} \theta & 0 & -\textrm{s} \theta \\ 0 & 1 & 0 \\ \textrm{s} \theta & 0 & \textrm{c} \theta \end{array} \right] +\left[ \begin{array}{ccc} \textrm{c} \phi & \textrm{s} \phi & 0 \\ -\textrm{s} \phi & \textrm{c} \phi & 0 \\ 0 & 0 & 1 \end{array} \right] \nonumber \\ +&= \left[ \begin{array}{ccc} \textrm{c} \theta \textrm{c} \phi & \textrm{c} \theta \textrm{s} \phi & -\textrm{s} \theta \\ \textrm{s} \psi \textrm{s} \theta \textrm{c} \phi - \textrm{c} \psi \textrm{s} \phi & \textrm{s} \psi \textrm{s} \theta \textrm{s} \phi + \textrm{c} \psi \textrm{c} \phi & \textrm{c} \theta \textrm{s} \psi \\ \textrm{c} \psi \textrm{s} \theta \textrm{c} \phi + \textrm{s} \psi \textrm{s} \phi & \textrm{c} \psi \textrm{s} \theta \textrm{s} \phi - \textrm{s} \psi \textrm{c} \phi & \textrm{c} \theta \textrm{c} \psi \end{array} \right] +\label{eq:CoordinateRotationYPR} +\end{align} + +Right-multiplying the direction cosine matrix with a vector that's expressed in the lab frame returns the coordinates of the same vector but now using the molecule frame. For example, suppose $\bm{r}_L$ is a vector expressed in the lab frame. Then, +\begin{equation} +\bm{r} = \bm{\Phi} \bm{r}_L +\end{equation} +gives the same vector but expressed in the molecule frame. (As a particularly simple vector, suppose $\bm{r}_L=[1,0,0]^T$; in the molecule frame, which is rotated counter-clockwise by $\phi$ from the lab frame, $\bm{r} = [c \phi,- s\phi, 0]^T$.) This expression can be transposed and/or multiplied on both sides by $\bm{\Phi}^T$ to yield +\begin{align} +\bm{r} &= \bm{\Phi} \bm{r}_L \hspace{2cm} \bm{r}_L = \bm{\Phi}^T \bm{r} \\ +\bm{r}^T &= \bm{r}_L^T \bm{\Phi}^T \hspace{1.5cm} \bm{r}_L^T = \bm{r}^T \bm{\Phi} +\end{align} +This library has functions for these multiplications. The transpose of the dcm listed above performs an active rotation of a vector in a constant coordinate system. + +Quaternions work in relatively similar ways (see the online article ``Rotation Quaternions, and How to Use Them'' by Rose, 2015). A quaternion, $\bm{q}$ has 4 values, $(q_0,q_1,q_2,q_3)$. If it represents a 3D point in space, then $\bm{q}=(0,x,y,z)$ and if it represents a rotation, then all four values can be non-zero and they have total magnitude of 1. The quaternion $\bm{q}=(1,0,0,0)$ is the rotation identity quaternion, which produces no rotation. \section{Code documentation} @@ -197,7 +253,7 @@ \subsection{Internal macros and global variables} \begin{description} \item[\ttt{\#define PI 3.14159265358979323846}] -hfill \\ +\hfill \\ $\pi$ \item[\ttt{double Work[9],Work2[9];}] @@ -212,117 +268,309 @@ \subsection{Externally accessible functions} \item[\underline{Vector conversion functions}] + \item[\ttt{void Sph\_Cart2Sc(double *Cart,double *Sc);}] \hfill \\ Converts Cartesian coordinates to spherical coordinates. + \item[\ttt{void Sph\_Sc2Cart(double *Sc,double *Cart);}] \hfill \\ -Converts spherical coordinates to Cartesian coordinates. - -\item[\ttt{void Sph\_Eay2Ep(double *Eay,double *Ep);}] -\hfill \\ -Converts Euler angle $y$-convention transformation to Euler parameters. Equations from Goldstein p. 608. - -\item[\ttt{void Sph\_Xyz2Xyz(double *Xyz1,double *Xyz2);}] -\hfill \\ -Copies yaw-pitch-roll vector (or any other 3D vector) \ttt{Xyz1} to \ttt{Xyz2}. This does not clamp angles. - -\item[\ttt{void Sph\_Eax2Xyz(double *Eax,double *Xyz);}] -\hfill \\ -Converts Euler angle $x$-convention to yaw-pitch-roll vector. Equations are from Sphere.nb. - - +Converts spherical coordinates to Cartesian coordinates. + + +\item[\ttt{void Sph\_Cart2Cart(const double *Cart1,double *Cart2);}] +\hfill \\ +Copies the 3D vector \ttt{Cart1} to the 3D vector \ttt{Cart2}. + + +\item[\ttt{void Sph\_Ypr2Ypr(double *Ypr1,double *Ypr2);}] +\hfill \\ +Copies yaw-pitch-roll vector (or any other 3D vector) \ttt{Ypr1} to \ttt{Ypr2}. This does not clamp angles. + + +\item[\ttt{void Sph\_Eax2Ypr(double *Eax,double *Ypr);}] +\hfill \\ +Converts Euler angle $x$-convention to yaw-pitch-roll vector. Equations are from Sphere.nb. + + \item[\ttt{void Sph\_Eax2Dcm(double *Eax,double *Dcm);}] -\hfill \\ -Calculates direction cosine matrix from Euler angle $x$-convention vector. Equations from Wolfram Mathworld and Sphere.nb. - +\hfill \\ +Calculates direction cosine matrix from Euler angle $x$-convention vector. Equations from Wolfram Mathworld and Sphere.nb. + + \item[\ttt{void Sph\_Eay2Dcm(double *Eay,double *Dcm);}] -\hfill \\ -Calculates direction cosine matrix from Euler angle $y$-convention vector. Equations from Wolfram Mathworld and Sphere.nb. - -\item[\ttt{void Sph\_Xyz2Dcm(double *Xyz,double *Dcm);}] -\hfill \\ -Calculates direction cosine matrix from yaw-pitch-roll vector. Equations from Goldstein p. 609 and Sphere.nb. $\bm{A} = \bm{X}(\psi)\bm{Y}(\theta)\bm{Z}(\phi)$. - -\item[\ttt{void Sph\_Xyz2Dcmt(double *Xyz,double *Dcmt);}] -\hfill \\ -Calculates transposed direction cosine matrix from yaw-pitch-roll vector. This is just \ttt{Sph\_Xyz2Dcm}, but for a transposed result. Rather than using this function, it's usually best to use the untransposed version in a \ttt{...Dcmt...} function. - -\item[\ttt{void Sph\_Dcm2Xyz(double *Dcm,double *Xyz);}] -\hfill \\ -Calculates yaw-pitch-roll vector from a direction cosine matrix. Equations derived from Goldstein p. 609 and from Sphere.nb. - +\hfill \\ +Calculates direction cosine matrix from Euler angle $y$-convention vector. Equations from Wolfram Mathworld and Sphere.nb. + + +\item[\ttt{void Sph\_Ypr2Dcm(double *Ypr,double *Dcm);}] +\hfill \\ +Calculates direction cosine matrix from yaw-pitch-roll vector. Equations from Goldstein p. 609 and Sphere.nb. $\bm{A} = \bm{X}(\psi)\bm{Y}(\theta)\bm{Z}(\phi)$. + + +\item[\ttt{void Sph\_Ypr2Dcmt(double *Ypr,double *Dcmt);}] +\hfill \\ +Calculates transposed direction cosine matrix from yaw-pitch-roll vector. This is just \ttt{Sph\_Ypr2Dcm}, but for a transposed result. Rather than using this function, it's usually best to use the untransposed version in a \ttt{...Dcmt...} function. + + +\item[\ttt{void Sph\_Dcm2Ypr(double *Dcm,double *Ypr);}] +\hfill \\ +Calculates yaw-pitch-roll vector from a direction cosine matrix. Equations derived from Goldstein p. 609 and from Sphere.nb. + + \item[\ttt{void Sph\_Dcm2Dcm(double *Dcm1,double *Dcm2);}] -\hfill \\ -Copies direction cosine matrix (or any other 9 element vector) \ttt{Dcm1} to a new one in \ttt{Dcm2}. - +\hfill \\ +Copies direction cosine matrix (or any other 9 element vector) \ttt{Dcm1} to a new one in \ttt{Dcm2}. + + \item[\ttt{void Sph\_Dcm2Dcmt(double *Dcm1,double *Dcm2);}] -\hfill \\ -Transposes direction cosine matrix \ttt{Dcm1} to yield matrix inverse in \ttt{Dcm2}. $\bm{A}_2 = \bm{A}_1^{-1}$. - +\hfill \\ +Transposes direction cosine matrix \ttt{Dcm1} to yield matrix inverse in \ttt{Dcm2}. $\bm{A}_2 = \bm{A}_1^{-1}$. + + \item[\ttt{void Sph\_DcmxDcm(double *Dcm1,double *Dcm2,double *Dcm3);}] -\hfill \\ -Matrix multiplies \ttt{Dcm1} by \ttt{Dcm2} and returns result in \ttt{Dcm3}. Note that the transformation is \ttt{Dcm2} first, then \ttt{Dcm1}. $\bm{A}_3 = \bm{A}_1\bm{A}_2$. - - +\hfill \\ +Matrix multiplies \ttt{Dcm1} by \ttt{Dcm2} and returns result in \ttt{Dcm3}. Note that the transformation is \ttt{Dcm2} first, then \ttt{Dcm1}. $\bm{A}_3 = \bm{A}_1\bm{A}_2$. + + \item[\ttt{void Sph\_DcmxDcmt(double *Dcm1,double *Dcmt,double *Dcm3);}] -\hfill \\ -Matrix multiplies \ttt{Dcm1} by the transpose of \ttt{Dcmt} and returns result in \ttt{Dcm3} (\ttt{Dcmt} is entered as an untransposed matrix). Essentially, this is a negative rotation of \ttt{Dcmt} followed by a positive rotation of \ttt{Dcm1}. $\bm{A}_3 = \bm{A}_1\bm{A}_2^{-1}$. - +\hfill \\ +Matrix multiplies \ttt{Dcm1} by the transpose of \ttt{Dcmt} and returns result in \ttt{Dcm3} (\ttt{Dcmt} is entered as an untransposed matrix). Essentially, this is a negative rotation of \ttt{Dcmt} followed by a positive rotation of \ttt{Dcm1}. $\bm{A}_3 = \bm{A}_1\bm{A}_2^{-1}$. + + \item[\ttt{void Sph\_DcmtxDcm(double *Dcmt,double *Dcm2,double *Dcm3);}] -\hfill \\ -Matrix multiplies the transpose of \ttt{Dcmt} by \ttt{Dcm2} and returns the result in \ttt{Dcm3} (\ttt{Dcmt} is entered as an untransposed matrix). Essentially, this is a positive rotation of \ttt{Dcm2} followed by a negative rotation of \ttt{Dcmt}. $\bm{A}_3 = \bm{A}_1^{-1}\bm{A}_2$. - +\hfill \\ +Matrix multiplies the transpose of \ttt{Dcmt} by \ttt{Dcm2} and returns the result in \ttt{Dcm3} (\ttt{Dcmt} is entered as an untransposed matrix). Essentially, this is a positive rotation of \ttt{Dcm2} followed by a negative rotation of \ttt{Dcmt}. $\bm{A}_3 = \bm{A}_1^{-1}\bm{A}_2$. + + \item[\ttt{void Sph\_DcmxCart(const double *Dcm,const double *Cart,double *Cart2);}] -\hfill \\ -Multiplies matrix \ttt{Dcm} by Cartesian vector \ttt{Cart} and returns the result in \ttt{Cart2}. This is a passive rotation of the coordinate system for the vector, $\bm{x}_2 = \bm{A}\cdot \bm{x}$. This is just a matrix times a vector, but is specifically for 3D. - +\hfill \\ +Multiplies matrix \ttt{Dcm} by Cartesian vector \ttt{Cart} and returns the result in \ttt{Cart2}. This is a passive rotation of the coordinate system for the vector, $\bm{x}_2 = \bm{A}\cdot \bm{x}$. This is just a matrix times a vector, but is specifically for 3D. + + \item[\ttt{void Sph\_One2Dcm(double *Dcm);}] -\hfill \\ -Returns the identity direction cosine matrix. $A_{ij} = \delta_{ij}$. - -\item[\ttt{void Sph\_Xyz2Xyzr(double *Xyz,double *Xyzr);}] -\hfill \\ -Converts the forward-direction yaw-pitch-roll vector \ttt{Xyz} to a relative direction change, but for travel in the reverse direction. For example, suppose an airplane performs the direction change that corresponds to \ttt{Xyz}. If it then turns around, with the local $z$-vector as it was initially, but with both $x$- and $y$-vectors reversed (180\degree yaw), then it needs to execute rotation \ttt{Xyzr} to retrace its original track. $\bm{A} = \bm{Z}^{-1}(\phi)\bm{Y}(\theta)\bm{X}(\psi)$. Note that this reverses a relative direction change between two vectors and does not reverse an absolute vector (the airplane traveling west being converted to it traveling east). Equations from Sphere.nb. - +\hfill \\ +Returns the identity direction cosine matrix. $A_{ij} = \delta_{ij}$. + + +\item[\ttt{void Sph\_Ypr2Yprr(double *Ypr,double *Yprr);}] +\hfill \\ +Converts the forward-direction yaw-pitch-roll vector \ttt{Ypr} to a relative direction change, but for travel in the reverse direction. For example, suppose an airplane performs the direction change that corresponds to \ttt{Ypr}. If it then turns around, with the local $z$-vector as it was initially, but with both $x$- and $y$-vectors reversed (180\degree yaw), then it needs to execute rotation \ttt{Yprr} to retrace its original track. $\bm{A} = \bm{Z}^{-1}(\phi)\bm{Y}(\theta)\bm{X}(\psi)$. Note that this reverses a relative direction change between two vectors and does not reverse an absolute vector (the airplane traveling west being converted to it traveling east). Equations from Sphere.nb. + + \item[\ttt{void Sph\_Dcm2Dcmr(double *Dcm,double *Dcmr);}] -\hfill \\ -Converts an absolute dcm to a dcm in the reverse direction. This reverses the local $x$ and $y$ directions, while preserving the local $z$ direction. This is unlike \ttt{Sph\_Xyz2Xyzr} in that this is for absolute directions while that one was for relative directions. $\bm{A}_r = \bm{Z}(\pi)\bm{A}$. - +\hfill \\ +Converts an absolute dcm to a dcm in the reverse direction. This reverses the local $x$ and $y$ directions, while preserving the local $z$ direction. This is unlike \ttt{Sph\_Ypr2Yprr} in that this is for absolute directions while that one was for relative directions. $\bm{A}_r = \bm{Z}(\pi)\bm{A}$. + + \item[\ttt{void Sph\_Rot2Dcm(char axis,double angle,double *Dcm);}] -\hfill \\ -Returns the direction cosine matrix that corresponds to rotation by angle \ttt{angle} about axis \ttt{axis}, where this latter parameter is the character \ttt{x}, \ttt{y}, or \ttt{z} (or upper-case). $\bm{A} = \bm{X}(a)$ or $\bm{A} = \bm{Y}(a)$ or $\bm{A} = \bm{Z}(a)$. - +\hfill \\ +Returns the direction cosine matrix that corresponds to rotation by angle \ttt{angle} about axis \ttt{axis}, where this latter parameter is the character \ttt{x}, \ttt{y}, or \ttt{z} (or upper-case). $\bm{A} = \bm{X}(a)$ or $\bm{A} = \bm{Y}(a)$ or $\bm{A} = \bm{Z}(a)$. + + \item[\ttt{void Sph\_Newz2Dcm(double *Newz,double psi,double *Dcm);}] -\hfill \\ -Returns the direction cosine matrix that can be used to rotate the coordinate system such that the original $z$-axis will line up with the vector \ttt{Newz}. The length of \ttt{Newz} is irrelevent; it does not need to be normalized. Additional rotation about the new $z$-axis is entered with \ttt{psi}. This works as follows: \ttt{Newz} is converted to spherical coordinates $\theta$ and $\phi$, then the dcm is $\bm{A} = \bm{Z}(\psi-\phi) \bm{X}(\theta) \bm{Z}(\phi)$, which is transposed to yield the active matrix. - +\hfill \\ +Returns the direction cosine matrix that can be used to rotate the coordinate system such that the original $z$-axis will line up with the vector \ttt{Newz}. The length of \ttt{Newz} is irrelevent; it does not need to be normalized. Additional rotation about the new $z$-axis is entered with \ttt{psi}. This works as follows: \ttt{Newz} is converted to spherical coordinates $\theta$ and $\phi$, then the dcm is $\bm{A} = \bm{Z}(\psi-\phi) \bm{X}(\theta) \bm{Z}(\phi)$, which is transposed to yield the active matrix. + + \item[\ttt{void Sph\_DcmtxUnit(double *Dcmt,char axis,double *vect,double *add,double mult);}] -\hfill \\ -Multiplies the transpose of \ttt{Dcmt} (entered as a non-transposed direction cosine matrix) with the unit vector for axis axis (entered as \ttt{x}, \ttt{y}, or \ttt{z}, or upper case) and returns the result in the 3-dimensional vector \ttt{vect}. This multiplies the result by the scalar \ttt{mult}. If \ttt{add} is non-\ttt{NULL}, this adds \ttt{add} to \ttt{vect} before returning the result. - +\hfill \\ +Multiplies the transpose of \ttt{Dcmt} (entered as a non-transposed direction cosine matrix) with the unit vector for axis \ttt{axis} (entered as \ttt{x}, \ttt{y}, or \ttt{z}, or upper case) and returns the result in the 3-dimensional vector \ttt{vect}. This multiplies the result by the scalar \ttt{mult}. If \ttt{add} is non-\ttt{NULL}, this adds \ttt{add} to \ttt{vect} before returning the result. + +% Quaternions +\item{\underline{Quaternions}} + +\item[\ttt{void Sph\_One2Qtn(double *Qtn)}] +\hfill \\ +Returns the rotation identity quaternion, which is $\bm{q} = (1,0,0,0)$. + + +\item[\ttt{void Sph\_Qtn2Qtn(const double *Qtn1,double *Qtn2)}] +\hfill \\ +Copies quaternion \ttt{Qtn1} to \ttt{Qtn2}. + + +\item[\ttt{void Sph\_Ypr2Qtn(const double Ypr,double *Qtn)}] +\hfill \\ +Converts Ypr angles to quaternion \ttt{Qtn}. This uses the following equation, which is from Rose 2015 but with sign changes because theirs is for active rotation and this is for passive rotation. +$$\bm{q} = \left[ \begin{array}{c} +\textrm{c}\frac{\psi}{2} \textrm{c}\frac{\theta}{2} \textrm{c}\frac{\phi}{2} + +\textrm{s}\frac{\psi}{2} \textrm{s}\frac{\theta}{2} \textrm{s}\frac{\phi}{2} \\ +-\textrm{s}\frac{\psi}{2} \textrm{c}\frac{\theta}{2} \textrm{c}\frac{\phi}{2} + +\textrm{c}\frac{\psi}{2} \textrm{s}\frac{\theta}{2} \textrm{s}\frac{\phi}{2} \\ +-\textrm{c}\frac{\psi}{2} \textrm{s}\frac{\theta}{2} \textrm{c}\frac{\phi}{2} - +\textrm{s}\frac{\psi}{2} \textrm{c}\frac{\theta}{2} \textrm{s}\frac{\phi}{2} \\ +-\textrm{c}\frac{\psi}{2} \textrm{c}\frac{\theta}{2} \textrm{s}\frac{\phi}{2} + +\textrm{s}\frac{\psi}{2} \textrm{s}\frac{\theta}{2} \textrm{c}\frac{\phi}{2} +\end{array} \right]$$ + + +\item[\ttt{void Sph\_Ypr2Qtni(const double Ypr,double *Qtni)}] +\hfill \\ +Converts Ypr angles to the inverse quaternion \ttt{Qtn}, but returned as a normal quaternion. This uses the following equation, from Rose, 2015. The signs are the same as in Rose because this is for an inverse passive rotation and theirs is for a forward active rotation. +$$\bm{q} = \left[ \begin{array}{c} +\textrm{c}\frac{\psi}{2} \textrm{c}\frac{\theta}{2} \textrm{c}\frac{\phi}{2} + +\textrm{s}\frac{\psi}{2} \textrm{s}\frac{\theta}{2} \textrm{s}\frac{\phi}{2} \\ +\textrm{s}\frac{\psi}{2} \textrm{c}\frac{\theta}{2} \textrm{c}\frac{\phi}{2} - +\textrm{c}\frac{\psi}{2} \textrm{s}\frac{\theta}{2} \textrm{s}\frac{\phi}{2} \\ +\textrm{c}\frac{\psi}{2} \textrm{s}\frac{\theta}{2} \textrm{c}\frac{\phi}{2} + +\textrm{s}\frac{\psi}{2} \textrm{c}\frac{\theta}{2} \textrm{s}\frac{\phi}{2} \\ +\textrm{c}\frac{\psi}{2} \textrm{c}\frac{\theta}{2} \textrm{s}\frac{\phi}{2} - +\textrm{s}\frac{\psi}{2} \textrm{s}\frac{\theta}{2} \textrm{c}\frac{\phi}{2} +\end{array} \right]$$ + + +\item[\ttt{void Sph\_Qtn2Ypr(const double *Qtn,double *Ypr)}] +\hfill \\ +Converts quaternion \ttt{Qtn} to ypr angles in \ttt{Ypr}. This uses the following equation, which is from Rose, 2015. +$$\bm{ypr}=\left[\begin{array}{c} +\atan2(2q_0q_3+2q_1q_2,q_0^2+q_1^2-q_2^2-q_3^2) \\ +\asin(2q_0q_2-2q_1q_3) \\ +\atan2(2q_0q_1+2q_2q_3,q_0^2-q_1^2-q_2^2+q_3^2) +\end{array} \right]$$ + + +\item[\ttt{void Sph\_Dcm2Qtn(const double *Dcm,double *Qtn)}] +\hfill \\ +Converts direction cosine matrix \ttt{Dcm} to quaternion \ttt{Qtn}. This is based on equations from Rose, 2015, but modified. Those equations are slightly slower to compute and also unstable to roundoff error (one can end up taking the square root of a negative number). My equations are as follows. First, compute +\begin{align*} +q_0 &= r_{11}+r_{22}+r_{33} \\ +q_1 &= r_{11}-r_{22}-r_{33} \\ +q_2 &= -r_{11}+r_{22}-r_{33} \\ +q_3 &= -r_{11}-r_{22}+r_{33} +\end{align*} +Next, determine which of these $q$ values is largest and recompute each one according the following equations: + +\begin{tabular}{llll} +$q_0$ largest & $q_1$ largest & $q_2$ largest & $q_3$ largest \\ +\hline +$q_0 = 0.5 \sqrt{1+q_0}$ & $q_1 = 0.5 \sqrt{1+q_1}$ & $q_2 = 0.5 \sqrt{1+q_2}$ & $q_3 = 0.5 \sqrt{1+q_3}$ \\ +$f = 0.25/q_0$ & $f = 0.25/q_1$ & $f = 0.25/q_2$ & $f = 0.25/q_3$ \\ +$q_1=f(r_{32}-r_{23})$ & $q_0=f(r_{32}-r_{23})$ & $q_0=f(r_{13}-r_{31})$ & $q_0=f(r_{21}-r_{12})$ \\ +$q_2=f(r_{13}-r_{31})$ & $q_2=f(r_{12}+r_{21})$ & $q_1=f(r_{12}+r_{21})$ & $q_1=f(r_{13}+r_{31})$ \\ +$q_3=f(r_{21}-r_{12})$ & $q_3=f(r_{13}+r_{31})$ & $q_3=f(r_{23}+r_{32})$ & $q_2=f(r_{23}+r_{32})$ +\end{tabular} + + +\item[\ttt{void Sph\_Qtn2Dcm(const double *Qtn,double *Dcm)}] +\hfill \\ +Converts quaternion \ttt{Qtn} to direction cosine matrix \ttt{Dcm}. This uses the equation, from Rose, 2015, +$$\bm{Q} = \left[ \begin{array}{ccc} +1-2q_2^2-2q_3^2 & 2q_1q_2-2q_0q_3 & 2q_1q_3+2q_0q_2\\ +2q_1q_2+2q_0q_3 & 1-2q_1^2-2q_3^2 & 2q_2q_3-2q_0q_1\\ +2q_1q_3-2q_0q_2 & 2q_2q_3+2q_0q_1 & 1-2q_1^2-2q_2^2 +\end{array} \right]$$ + + +\item[\ttt{void Sph\_QtnxQtn(const double Qtn1,const double Qtn2,double *Qtn3)}] +\hfill \\ +Multiplies quaternion \ttt{Qtn1} by \ttt{Qtn2} to yield \ttt{Qtn3}. Equation is from Rose, 2015. For the product $\bm{t}=\bm{r}\bm{s}$, +$$\bm{t}=\left[\begin{array}{c} +r_0s_0 - r_1s_1 - r_2s_2 - r_3s_3 \\ +r_0s_1 + r_1s_0 - r_2s_3 + r_3s_2 \\ +r_0s_2 + r_1s_3 + r_2s_0 - r_3s_1 \\ +r_0s_3 - r_1s_2 + r_2s_1 + r_3s_0 +\end{array} \right]$$ + + +\item[\ttt{void Sph\_QtnixQtn(const double Qtn1,const double Qtn2,double *Qtn3)}] +\hfill \\ +Multiplies the inverse of quaternion \ttt{Qtn1} (entered in its non-inverse version) by \ttt{Qtn2} to yield \ttt{Qtn3}. For the product $\bm{t}=\bm{r}^{-1}\bm{s}$, +$$\bm{t}=\left[\begin{array}{c} +r_0s_0 + r_1s_1 + r_2s_2 + r_3s_3 \\ +r_0s_1 - r_1s_0 + r_2s_3 - r_3s_2 \\ +r_0s_2 - r_1s_3 - r_2s_0 + r_3s_1 \\ +r_0s_3 + r_1s_2 - r_2s_1 - r_3s_0 +\end{array} \right]$$ + + +\item[\ttt{void Sph\_QtnxQtni(const double Qtn1,const double Qtn2,double *Qtn3)}] +\hfill \\ +Multiplies quaternion \ttt{Qtn1} by the inverse of \ttt{Qtn2} (entered in its non-inverse version) to yield \ttt{Qtn3}. For the product $\bm{t}=\bm{r}^{-1}\bm{s}$, +$$\bm{t}=\left[\begin{array}{c} +r_0s_0 + r_1s_1 + r_2s_2 + r_3s_3 \\ +-r_0s_1 + r_1s_0 + r_2s_3 - r_3s_2 \\ +-r_0s_2 - r_1s_3 + r_2s_0 + r_3s_1 \\ +-r_0s_3 + r_1s_2 - r_2s_1 + r_3s_0 +\end{array} \right]$$ + + +\item[\ttt{void}] +\ttt{Sph\_QtnRotate(const double *Qtn,const double *Cart,double *Cart2);} +\hfill \\ +Rotates 3D Cartesian vector \ttt{Cart} using quaternion \ttt{Qtn}, returning the answer in 3D vector \ttt{Cart2}. This computes the product +$$\bm{c}_2 = \bm{q}^* \bm{c} \bm{q}$$ +where $\bm{c}_2$ is \ttt{Cart2}, $\bm{q}$ is \ttt{Qtn}, and $\bm{c}$ is \ttt{Cart}. The math is derived in the Mathematica document SphereQuat.nb. This equation performs passive rotation from the $\bm{c}$ frame to the $\bm{c}_2$ frame, where the rotation angles represent the rotation of $\bm{c}_2$ from $\bm{c}$ (e.g. $\bm{c}$ is the lab frame and $\bm{c}_2$ is the molecule frame). + + +\item[\ttt{void}] +\ttt{Sph\_QtniRotate(const double *Qtn,const double *Cart,double *Cart2);} +\hfill \\ +Performs the inverse rotation of 3D Cartesian vector \ttt{Cart} using quaternion \ttt{Qtn}, returning the answer in 3D vector \ttt{Cart2}. This computes the product +$$\bm{c}_2 = \bm{q} \bm{c} \bm{q}^*$$ +where $\bm{c}_2$ is \ttt{Cart2}, $\bm{q}$ is \ttt{Qtn}, and $\bm{c}$ is \ttt{Cart}. The math is derived in the Mathematica document SphereQuat.nb. See the discussion for \ttt{Sph\_QtnRotate}. + + +\item[\ttt{void}] +\ttt{Sph\_QtniRotateUnitx(const double *Qtni, double *vect, const double *add, double mult)} +\hfill \\ +Rotates the $\hat{\bm{x}}$ unit vector with the inverse quaternion \ttt{Qtni}, which is entered as a non-inverted quaternion, and then multiplies the result by \ttt{mult} and adds it to \ttt{add}, returning the answer in \ttt{vect}. The \ttt{add} vector is required and cannot be the same as \ttt{vect}. Using $\bm{a}$ for \ttt{add}, $m$ for \ttt{mult}, $\bm{v}$ for \ttt{vect}, and $\bm{q}$ for \ttt{Qtni}, the equation for this is +$$\bm{v} = \bm{a} + m \left[ \begin{array}{c} +q_0^2 + q_1^2 - q_2^2 - q_3^2 \\ +2 q_1 q_2 - 2 q_0 q_3 \\ +2 q_0 q_2 + 2 q_1 q_3 +\end{array} \right]$$ +This is derived in the Mathematica document SphereQuat.nb. + + +\item[\ttt{void}] +\ttt{Sph\_QtniRotateUnitz(const double *Qtni, double *vect, const double *add, double mult)} +\hfill \\ +Rotates the $\hat{\bm{z}}$ unit vector with the inverse quaternion \ttt{Qtni}, which is entered as a non-inverted quaternion, and then multiplies the result by \ttt{mult} and adds it to \ttt{add}, returning the answer in \ttt{vect}. The \ttt{add} vector is required and cannot be the same as \ttt{vect}. Using $\bm{a}$ for \ttt{add}, $m$ for \ttt{mult}, $\bm{v}$ for \ttt{vect}, and $\bm{q}$ for \ttt{Qtni}, the equation for this is +$$\bm{v} = \bm{a} + m \left[ \begin{array}{c} +-2 q_0 q_2 + 2 q_1 q_3 \\ +2 q_0 q_1 + 2 q_2 q_3 \\ +q_0^2 - q_1^2 - q_2^2 + q_3^2 +\end{array} \right]$$ +This is derived in the Mathematica document SphereQuat.nb. + +% More rotations +\item{\underline{More rotations}} + +\item[\ttt{void}] +\ttt{Sph\_RotateAxisAngle(const double *pt1,const double *pt2,double theta,const double *oldvect,double *newvect)} +\hfill \\ +Text. + + \item[\ttt{double}] \ttt{Sph\_RotateVectWithNormals3D(double *pt1, double *pt2, double *newpt2, double *oldnorm, double *newnorm, int sign);} -\hfill \\ -This is for the case where the line from \ttt{pt1} to \ttt{pt2} is in the plane that has normal \ttt{oldnorm}, and then the plane is rotated about point \ttt{pt1} to so that its normal becomes \ttt{newnorm}. This function calculates the new value for \ttt{pt2}, returned in \ttt{newpt2}. \ttt{newpt2} and \ttt{pt2} are allowed to point to the same memory. Both \ttt{oldnorm} and \ttt{newnorm} need to have unit length. This returns the cosine of the angle between the two normals, which is also the dot product of the two normal vectors. If this cosine is 1, then the two normals are parallel to each other and \ttt{newpt2} is set equal to \ttt{pt2} because no rotation takes place. If this cosine is -1, then the two normals are anti-parallel to each other, in which case the problem is ill-determined because the rotation axis cannot be determined; if that's the case, then this function assumes that the rotation axis is perpendicular to the vector from \ttt{pt1} to \ttt{pt2}, with the result that the new vector is in the opposite direction as the original vector. New function Sept. 2015. - -The \ttt{sign} input is here to allow the normals to internally inconsistent. That is, it is good practice for all normals to points toward the same face of a surface, such as the outside or inside. If this is the case, then enter \ttt{sign} as 0. However, if this is not done, then enter \ttt{sign} as 1 if the total rotation should be less than 90\degree and as -1 if the total rotation should be more than 90\degree. - -It is permitted to enter \ttt{oldnorm} as \ttt{NULL}. In this case, the vector is rotated around a random rotation axis that is perpendicular to \ttt{newnorm}. In other words, \ttt{newpt2} is still placed in the new plane and it is still the correct distance from \ttt{pt1}, but the rotation direction to this new position is random. - +\hfill \\ +This is for the case where the line from \ttt{pt1} to \ttt{pt2} is in the plane that has normal \ttt{oldnorm}, and then the plane is rotated about point \ttt{pt1} to so that its normal becomes \ttt{newnorm}. This function calculates the new value for \ttt{pt2}, returned in \ttt{newpt2}. \ttt{newpt2} and \ttt{pt2} are allowed to point to the same memory. Both \ttt{oldnorm} and \ttt{newnorm} need to have unit length. This returns the cosine of the angle between the two normals, which is also the dot product of the two normal vectors. If this cosine is 1, then the two normals are parallel to each other and \ttt{newpt2} is set equal to \ttt{pt2} because no rotation takes place. If this cosine is -1, then the two normals are anti-parallel to each other, in which case the problem is ill-determined because the rotation axis cannot be determined; if that's the case, then this function assumes that the rotation axis is perpendicular to the vector from \ttt{pt1} to \ttt{pt2}, with the result that the new vector is in the opposite direction as the original vector. New function Sept. 2015. + +The \ttt{sign} input is here to allow the normals to internally inconsistent. That is, it is good practice for all normals to points toward the same face of a surface, such as the outside or inside. If this is the case, then enter \ttt{sign} as 0. However, if this is not done, then enter \ttt{sign} as 1 if the total rotation should be less than 90\degree and as -1 if the total rotation should be more than 90\degree. + +It is permitted to enter \ttt{oldnorm} as \ttt{NULL}. In this case, the vector is rotated around a random rotation axis that is perpendicular to \ttt{newnorm}. In other words, \ttt{newpt2} is still placed in the new plane and it is still the correct distance from \ttt{pt1}, but the rotation direction to this new position is random. + The math is as follows. Define $\bm{p}_1$ as \ttt{pt1}, $\bm{p}_2$ as \ttt{pt2}, $\bm{o}$ as \ttt{oldnorm}, and $\bm{n}$ as \ttt{newnorm}. Also, define $\bm{p}$ as the vector from $\bm{p}_1$ to $\bm{p}_2$, meaning that $\bm{p} = \bm{p}_2 - \bm{p}_1$. Also define $\bm{a}$ as the unit vector for the axis about which the rotation takes place; it is the line that is shared by the old plane and the new plane. Define $\theta$ as the rotation angle about this axis. These values are $$\bm{a} = \frac{\bm{o} \times \bm{n}}{\sqrt{(\bm{o} \times \bm{n})\cdot (\bm{o} \times \bm{n})}}$$ $$\cos \theta = \bm{o} \cdot \bm{n}$$ - -The $\theta$ equation relies on the requirement that $\bm{o}$ and $\bm{n}$ have unit length. The direction cosine matrix for rotation by angle $\theta$ about axis $\bm{a}$ is (from Wikipedia ``Rotation matrix'') - + +The $\theta$ equation relies on the requirement that $\bm{o}$ and $\bm{n}$ have unit length. The direction cosine matrix for rotation by angle $\theta$ about axis $\bm{a}$ is (from Wikipedia ``Rotation matrix'') + $$\left[ \begin{array}{ccc} c\theta+ a_x^2(1-c\theta) & a_x a_y(1-c\theta)-a_z s\theta & a_x a_z (1-c\theta)+a_y s\theta \\ a_y a_x(1-c\theta)+a_zs\theta & c\theta + a_y^2(1-c\theta) & a_y z_z(1-c\theta)-a_x s\theta \\ a_z a_x(1-c\theta)-a_y s\theta & a_z a_y(1-c\theta) + a_x s\theta & c\theta + a_z^2(1-c\theta) \end{array} \right ]$$ + +\item[\ttt{void}] +\ttt{Sph\_RotateVectorAxisAngle(const double *vect, const double *axis, double angle, double *rotated)} +\hfill \\ +Rotates vector \ttt{vect} about axis \ttt{axis} by angle \ttt{angle}, returning the answer in \ttt{rotated}. This uses Rodrigues's rotation formula. I didn't write this function, but told ChatGPT to write it for me, and then I edited the result. + \end{description} \bibliographystyle{plain} diff --git a/examples/S13_filaments/FilDefineTest.txt b/examples/S13_filaments/FilDefineTest.txt new file mode 100644 index 00000000..392368ea --- /dev/null +++ b/examples/S13_filaments/FilDefineTest.txt @@ -0,0 +1,55 @@ +# Filaments + +graphics opengl +random_seed 2 + +dim 3 +boundaries x -5 5 r +boundaries y -5 5 r +boundaries z -5 5 r +frame_color lightblue + +species red + +difc red 3 + +color red red + +time_start 0 +time_stop 1000 +time_step 0.01 + + +start_filament_type black +dynamics none +color black +thickness 2 +polygon ve +kT 1 +standard_length 1 +standard_angle 0 0 0 +force_length 1 +force_angle 1 1 1 +end_filament_type + +random_filament black:fil1 5 0 0 0 0 0 0 + +start_filament black:fil2 +first_segment -3 0 0 1 0 0 0 +add_segment 1 1 0 0 2 back +add_segment 1 0 1 0 2 back +add_segment 1 1 0 0 2 front +copy_to fil3 +end_filament + +filament fil3 translate + 0 2 0 +filament fil3 remove_segment front +filament fil3 remove_segment back + +filament fil2 modify_segment 2 angle + 1 0 0 back + + +end_file + + + diff --git a/examples/S13_filaments/filament.txt b/examples/S13_filaments/filament.txt index 5f66e92f..04d35f9f 100644 --- a/examples/S13_filaments/filament.txt +++ b/examples/S13_filaments/filament.txt @@ -1,7 +1,7 @@ # Filaments graphics opengl -#random_seed 2 +random_seed 2 #graphics opengl_better dim 3 @@ -37,218 +37,218 @@ force_angle 100 0 0 end_filament_type -start_filament black template +start_filament black:template random_segments 15 10 50 1 u u u -#copy fil2 +#copy_to fil2 #name black fil2 #translate + 15 15 0 end_filament -filament template copy fil2 +filament template copy_to fil2 filament fil2 translate + 15 15 0 -filament template copy fil3 -filament template copy fil4 -filament template copy fil5 -filament template copy fil6 -filament template copy fil7 -filament template copy fil8 -filament template copy fil9 +filament template copy_to fil3 +filament template copy_to fil4 +filament template copy_to fil5 +filament template copy_to fil6 +filament template copy_to fil7 +filament template copy_to fil8 +filament template copy_to fil9 /* -filament template copy fil10 -filament template copy fil11 -filament template copy fil12 -filament template copy fil13 -filament template copy fil14 -filament template copy fil15 -filament template copy fil16 -filament template copy fil17 -filament template copy fil18 -filament template copy fil19 -filament template copy fil20 -filament template copy fil21 -filament template copy fil22 -filament template copy fil23 -filament template copy fil24 -filament template copy fil25 -filament template copy fil26 -filament template copy fil27 -filament template copy fil28 -filament template copy fil29 -filament template copy fil30 -filament template copy fil31 -filament template copy fil32 -filament template copy fil33 -filament template copy fil34 -filament template copy fil35 -filament template copy fil36 -filament template copy fil37 -filament template copy fil38 -filament template copy fil39 -filament template copy fil40 -filament template copy fil41 -filament template copy fil42 -filament template copy fil43 -filament template copy fil44 -filament template copy fil45 -filament template copy fil46 -filament template copy fil47 -filament template copy fil48 -filament template copy fil49 -filament template copy fil50 -filament template copy fil51 -filament template copy fil52 -filament template copy fil53 -filament template copy fil54 -filament template copy fil55 -filament template copy fil56 -filament template copy fil57 -filament template copy fil58 -filament template copy fil59 -filament template copy fil60 -filament template copy fil61 -filament template copy fil62 -filament template copy fil63 -filament template copy fil64 -filament template copy fil65 -filament template copy fil66 -filament template copy fil67 -filament template copy fil68 -filament template copy fil69 -filament template copy fil70 -filament template copy fil71 -filament template copy fil72 -filament template copy fil73 -filament template copy fil74 -filament template copy fil75 -filament template copy fil76 -filament template copy fil77 -filament template copy fil78 -filament template copy fil79 -filament template copy fil80 -filament template copy fil81 -filament template copy fil82 -filament template copy fil83 -filament template copy fil84 -filament template copy fil85 -filament template copy fil86 -filament template copy fil87 -filament template copy fil88 -filament template copy fil89 -filament template copy fil90 -filament template copy fil91 -filament template copy fil92 -filament template copy fil93 -filament template copy fil94 -filament template copy fil95 -filament template copy fil96 -filament template copy fil97 -filament template copy fil98 -filament template copy fil99 +filament template copy_to fil10 +filament template copy_to fil11 +filament template copy_to fil12 +filament template copy_to fil13 +filament template copy_to fil14 +filament template copy_to fil15 +filament template copy_to fil16 +filament template copy_to fil17 +filament template copy_to fil18 +filament template copy_to fil19 +filament template copy_to fil20 +filament template copy_to fil21 +filament template copy_to fil22 +filament template copy_to fil23 +filament template copy_to fil24 +filament template copy_to fil25 +filament template copy_to fil26 +filament template copy_to fil27 +filament template copy_to fil28 +filament template copy_to fil29 +filament template copy_to fil30 +filament template copy_to fil31 +filament template copy_to fil32 +filament template copy_to fil33 +filament template copy_to fil34 +filament template copy_to fil35 +filament template copy_to fil36 +filament template copy_to fil37 +filament template copy_to fil38 +filament template copy_to fil39 +filament template copy_to fil40 +filament template copy_to fil41 +filament template copy_to fil42 +filament template copy_to fil43 +filament template copy_to fil44 +filament template copy_to fil45 +filament template copy_to fil46 +filament template copy_to fil47 +filament template copy_to fil48 +filament template copy_to fil49 +filament template copy_to fil50 +filament template copy_to fil51 +filament template copy_to fil52 +filament template copy_to fil53 +filament template copy_to fil54 +filament template copy_to fil55 +filament template copy_to fil56 +filament template copy_to fil57 +filament template copy_to fil58 +filament template copy_to fil59 +filament template copy_to fil60 +filament template copy_to fil61 +filament template copy_to fil62 +filament template copy_to fil63 +filament template copy_to fil64 +filament template copy_to fil65 +filament template copy_to fil66 +filament template copy_to fil67 +filament template copy_to fil68 +filament template copy_to fil69 +filament template copy_to fil70 +filament template copy_to fil71 +filament template copy_to fil72 +filament template copy_to fil73 +filament template copy_to fil74 +filament template copy_to fil75 +filament template copy_to fil76 +filament template copy_to fil77 +filament template copy_to fil78 +filament template copy_to fil79 +filament template copy_to fil80 +filament template copy_to fil81 +filament template copy_to fil82 +filament template copy_to fil83 +filament template copy_to fil84 +filament template copy_to fil85 +filament template copy_to fil86 +filament template copy_to fil87 +filament template copy_to fil88 +filament template copy_to fil89 +filament template copy_to fil90 +filament template copy_to fil91 +filament template copy_to fil92 +filament template copy_to fil93 +filament template copy_to fil94 +filament template copy_to fil95 +filament template copy_to fil96 +filament template copy_to fil97 +filament template copy_to fil98 +filament template copy_to fil99 */ /* -filament template copy fil102 -filament template copy fil103 -filament template copy fil104 -filament template copy fil105 -filament template copy fil106 -filament template copy fil107 -filament template copy fil108 -filament template copy fil109 -filament template copy fil110 -filament template copy fil111 -filament template copy fil112 -filament template copy fil113 -filament template copy fil114 -filament template copy fil115 -filament template copy fil116 -filament template copy fil117 -filament template copy fil118 -filament template copy fil119 -filament template copy fil120 -filament template copy fil121 -filament template copy fil122 -filament template copy fil123 -filament template copy fil124 -filament template copy fil125 -filament template copy fil126 -filament template copy fil127 -filament template copy fil128 -filament template copy fil129 -filament template copy fil130 -filament template copy fil131 -filament template copy fil132 -filament template copy fil133 -filament template copy fil134 -filament template copy fil135 -filament template copy fil136 -filament template copy fil137 -filament template copy fil138 -filament template copy fil139 -filament template copy fil140 -filament template copy fil141 -filament template copy fil142 -filament template copy fil143 -filament template copy fil144 -filament template copy fil145 -filament template copy fil146 -filament template copy fil147 -filament template copy fil148 -filament template copy fil149 -filament template copy fil150 -filament template copy fil151 -filament template copy fil152 -filament template copy fil153 -filament template copy fil154 -filament template copy fil155 -filament template copy fil156 -filament template copy fil157 -filament template copy fil158 -filament template copy fil159 -filament template copy fil160 -filament template copy fil161 -filament template copy fil162 -filament template copy fil163 -filament template copy fil164 -filament template copy fil165 -filament template copy fil166 -filament template copy fil167 -filament template copy fil168 -filament template copy fil169 -filament template copy fil170 -filament template copy fil171 -filament template copy fil172 -filament template copy fil173 -filament template copy fil174 -filament template copy fil175 -filament template copy fil176 -filament template copy fil177 -filament template copy fil178 -filament template copy fil179 -filament template copy fil180 -filament template copy fil181 -filament template copy fil182 -filament template copy fil183 -filament template copy fil184 -filament template copy fil185 -filament template copy fil186 -filament template copy fil187 -filament template copy fil188 -filament template copy fil189 -filament template copy fil190 -filament template copy fil191 -filament template copy fil192 -filament template copy fil193 -filament template copy fil194 -filament template copy fil195 -filament template copy fil196 -filament template copy fil197 -filament template copy fil198 -filament template copy fil199 +filament template copy_to fil102 +filament template copy_to fil103 +filament template copy_to fil104 +filament template copy_to fil105 +filament template copy_to fil106 +filament template copy_to fil107 +filament template copy_to fil108 +filament template copy_to fil109 +filament template copy_to fil110 +filament template copy_to fil111 +filament template copy_to fil112 +filament template copy_to fil113 +filament template copy_to fil114 +filament template copy_to fil115 +filament template copy_to fil116 +filament template copy_to fil117 +filament template copy_to fil118 +filament template copy_to fil119 +filament template copy_to fil120 +filament template copy_to fil121 +filament template copy_to fil122 +filament template copy_to fil123 +filament template copy_to fil124 +filament template copy_to fil125 +filament template copy_to fil126 +filament template copy_to fil127 +filament template copy_to fil128 +filament template copy_to fil129 +filament template copy_to fil130 +filament template copy_to fil131 +filament template copy_to fil132 +filament template copy_to fil133 +filament template copy_to fil134 +filament template copy_to fil135 +filament template copy_to fil136 +filament template copy_to fil137 +filament template copy_to fil138 +filament template copy_to fil139 +filament template copy_to fil140 +filament template copy_to fil141 +filament template copy_to fil142 +filament template copy_to fil143 +filament template copy_to fil144 +filament template copy_to fil145 +filament template copy_to fil146 +filament template copy_to fil147 +filament template copy_to fil148 +filament template copy_to fil149 +filament template copy_to fil150 +filament template copy_to fil151 +filament template copy_to fil152 +filament template copy_to fil153 +filament template copy_to fil154 +filament template copy_to fil155 +filament template copy_to fil156 +filament template copy_to fil157 +filament template copy_to fil158 +filament template copy_to fil159 +filament template copy_to fil160 +filament template copy_to fil161 +filament template copy_to fil162 +filament template copy_to fil163 +filament template copy_to fil164 +filament template copy_to fil165 +filament template copy_to fil166 +filament template copy_to fil167 +filament template copy_to fil168 +filament template copy_to fil169 +filament template copy_to fil170 +filament template copy_to fil171 +filament template copy_to fil172 +filament template copy_to fil173 +filament template copy_to fil174 +filament template copy_to fil175 +filament template copy_to fil176 +filament template copy_to fil177 +filament template copy_to fil178 +filament template copy_to fil179 +filament template copy_to fil180 +filament template copy_to fil181 +filament template copy_to fil182 +filament template copy_to fil183 +filament template copy_to fil184 +filament template copy_to fil185 +filament template copy_to fil186 +filament template copy_to fil187 +filament template copy_to fil188 +filament template copy_to fil189 +filament template copy_to fil190 +filament template copy_to fil191 +filament template copy_to fil192 +filament template copy_to fil193 +filament template copy_to fil194 +filament template copy_to fil195 +filament template copy_to fil196 +filament template copy_to fil197 +filament template copy_to fil198 +filament template copy_to fil199 */ #filament fil2 color red diff --git a/examples/S13_filaments/filament1.txt b/examples/S13_filaments/filament1.txt index a564f3f9..804b4a68 100644 --- a/examples/S13_filaments/filament1.txt +++ b/examples/S13_filaments/filament1.txt @@ -34,11 +34,11 @@ force_length 10 force_angle 5 -1 -1 end_filament_type -random_filament fil1 tread 10 10 50 1 0 0 0 -random_filament fil2 tread 10 20 50 1 1 0 0 -random_filament fil3 tread 10 10 40 1 -1 0 0 -random_filament fil4 tread 10 30 50 1 2 0 0 -random_filament fil5 tread 10 30 30 1 -1 0 0 +random_filament tread:fil1 10 10 50 1 0 0 0 +random_filament tread:fil2 10 20 50 1 1 0 0 +random_filament tread:fil3 10 10 40 1 -1 0 0 +random_filament tread:fil4 10 30 50 1 2 0 0 +random_filament tread:fil5 10 30 30 1 -1 0 0 start_surface bounds diff --git a/examples/S13_filaments/filament2.txt b/examples/S13_filaments/filament2.txt index cd140836..4bd9f56f 100644 --- a/examples/S13_filaments/filament2.txt +++ b/examples/S13_filaments/filament2.txt @@ -33,24 +33,24 @@ force_length 1 force_angle 100 80 -1 end_filament_type -new_filament redtype red -new_filament redtype red2 -new_filament redtype red3 -new_filament redtype red4 -new_filament redtype red5 -new_filament redtype red6 -new_filament redtype red7 -new_filament redtype red8 -new_filament redtype red9 - -filament red copy red2 -filament red copy red3 -filament red copy red4 -filament red copy red5 -filament red copy red6 -filament red copy red7 -filament red copy red8 -filament red copy red9 +new_filament redtype:red +new_filament redtype:red2 +new_filament redtype:red3 +new_filament redtype:red4 +new_filament redtype:red5 +new_filament redtype:red6 +new_filament redtype:red7 +new_filament redtype:red8 +new_filament redtype:red9 + +filament red copy_to red2 +filament red copy_to red3 +filament red copy_to red4 +filament red copy_to red5 +filament red copy_to red6 +filament red copy_to red7 +filament red copy_to red8 +filament red copy_to red9 filament red first_segment 46.19222236 50.60101155 1 5 4.656180038 0 0 filament red2 first_segment 76.90132296 26.70912275 1 5 4.351468544 0 0 @@ -83,16 +83,16 @@ force_length -1 force_angle 10 0 0 end_filament_type -new_filament greentype green +new_filament greentype:green -filament green copy green2 -filament green copy green3 -filament green copy green4 -filament green copy green5 -filament green copy green6 -filament green copy green7 -filament green copy green8 -filament green copy green9 +filament green copy_to green2 +filament green copy_to green3 +filament green copy_to green4 +filament green copy_to green5 +filament green copy_to green6 +filament green copy_to green7 +filament green copy_to green8 +filament green copy_to green9 filament green first_segment 68.21505086 46.64608437 1 2 1.691868734 0 0 filament green2 first_segment 51.732133 55.70341596 1 2 0.764343705 0 0 filament green3 first_segment 99.77730256 10.79488604 1 2 1.517531223 0 0 diff --git a/examples/S13_filaments/filament2d.txt b/examples/S13_filaments/filament2d.txt index ef2ae6d7..d3d481b1 100644 --- a/examples/S13_filaments/filament2d.txt +++ b/examples/S13_filaments/filament2d.txt @@ -37,20 +37,22 @@ force_angle 10 mobility 10 end_filament_type -#start_filament template fil1 +#start_filament template:fil1 #first_segment 0 0 5 0 #add_segment 5 0.3 #add_segment 5 0.3 #end_filament -random_filament fil1 template 8 10 50 u -random_filament fil2 template 49 20 50 u -random_filament fil3 template 20 30 50 u -#random_filament fil4 template 49 40 50 u -#random_filament fil5 template 49 50 50 u +random_filament template:fil1 8 10 50 u +random_filament template:fil2 49 20 50 u +random_filament template:fil3 20 30 50 u +#random_filament template:fil4 49 40 50 u +#random_filament template:fil5 49 50 50 u +filament_type template standard_angle 0.3 + cmd B pause #cmd E printFilament template fil1 xabf stdout diff --git a/examples/S13_filaments/filament3.txt b/examples/S13_filaments/filament3.txt index 41e94160..392c5697 100644 --- a/examples/S13_filaments/filament3.txt +++ b/examples/S13_filaments/filament3.txt @@ -1,5 +1,12 @@ # Filaments +# Current issue. Run this for a few steps and stop. +# The computed YPR values are not correct. The segment 0 ones have the wrong sign +# and the segment 1 ones never change. Don't change the qtn2ypr signs. See SmoldynCodeDoc +# eq. 5.12. +# I suspect that XY2Qtn should return the inverse qtn instead. + + graphics opengl #random_seed 2 #graphics opengl_better @@ -16,8 +23,8 @@ difc red 3 color red white time_start 0 -time_stop 50000 -time_step 0.3 +time_stop 10000 +time_step 0.5 frame_thickness 0 @@ -29,7 +36,7 @@ force_arrows 10 red thickness 2 polygon ve kT 0.1 -dynamics Euler +dynamics RK4 standard_length 5 standard_angle 0 0 0 force_length 1 @@ -37,26 +44,28 @@ force_angle 2 2 2 mobility 0.5 end_filament_type -#start_filament template fil1 +#start_filament template:fil1 #first_segment 0 0 0 5 0 0 0 #add_segment 5 1 0 0 #add_segment 5 0 1 0 #add_segment 5 0 0 1 #end_filament -random_filament fil1 template 10 0 0 0 u u u -random_filament fil2 template 20 20 50 1 u u u -random_filament fil3 template 15 30 50 1 u u u -random_filament fil4 template 30 40 50 1 u u u -random_filament fil5 template 25 50 50 1 u u u +random_filament template:fil1 10 0 0 0 u u u +#random_filament template:fil1 10 0 0 0 0 0 0 + + +random_filament template:fil2 20 20 50 1 u u u +random_filament template:fil3 15 30 50 1 u u u +random_filament template:fil4 30 40 50 1 u u u +random_filament template:fil5 25 50 50 1 u u u + +filament_type template standard_angle 0.3 0 0 -start_filament_type template -#standard_angle 0.3 0.3 0.1 -standard_angle 0.3 0 0 -end_filament_type #cmd B printFilament template fil1 xabyft stdout cmd B pause +#cmd i 0 10000 1000 printFilament template fil1 xabyft stdout end_file diff --git a/examples/S13_filaments/filament4.txt b/examples/S13_filaments/filament4.txt index 72b5b22d..5cc3362a 100644 --- a/examples/S13_filaments/filament4.txt +++ b/examples/S13_filaments/filament4.txt @@ -32,16 +32,16 @@ force_length -1 force_angle 100 80 0 end_filament_type -random_filament r0 red 10 u u u u u u -random_filament r1 red 10 u u u u u u -random_filament r2 red 10 u u u u u u -random_filament r3 red 10 u u u u u u -random_filament r4 red 10 u u u u u u -random_filament r5 red 10 u u u u u u -random_filament r6 red 10 u u u u u u -random_filament r7 red 10 u u u u u u -random_filament r8 red 10 u u u u u u -random_filament r9 red 10 u u u u u u +random_filament red:r0 10 u u u u u u +random_filament red:r1 10 u u u u u u +random_filament red:r2 10 u u u u u u +random_filament red:r3 10 u u u u u u +random_filament red:r4 10 u u u u u u +random_filament red:r5 10 u u u u u u +random_filament red:r6 10 u u u u u u +random_filament red:r7 10 u u u u u u +random_filament red:r8 10 u u u u u u +random_filament red:r9 10 u u u u u u start_filament_type green @@ -56,16 +56,16 @@ force_length -1 force_angle 10 10 10 end_filament_type -random_filament g0 green 80 u u u u u u -random_filament g1 green 80 u u u u u u -random_filament g2 green 80 u u u u u u -random_filament g3 green 80 u u u u u u -random_filament g4 green 80 u u u u u u -random_filament g5 green 80 u u u u u u -random_filament g6 green 80 u u u u u u -random_filament g7 green 80 u u u u u u -random_filament g8 green 80 u u u u u u -random_filament g9 green 80 u u u u u u +random_filament green:g0 80 u u u u u u +random_filament green:g1 80 u u u u u u +random_filament green:g2 80 u u u u u u +random_filament green:g3 80 u u u u u u +random_filament green:g4 80 u u u u u u +random_filament green:g5 80 u u u u u u +random_filament green:g6 80 u u u u u u +random_filament green:g7 80 u u u u u u +random_filament green:g8 80 u u u u u u +random_filament green:g9 80 u u u u u u end_file diff --git a/examples/S13_filaments/filament5.txt b/examples/S13_filaments/filament5.txt index 588fab79..4908a5a8 100644 --- a/examples/S13_filaments/filament5.txt +++ b/examples/S13_filaments/filament5.txt @@ -42,16 +42,16 @@ force_length -1 force_angle 100 80 0 end_filament_type -random_filament r0 red 10 u u u u u u -random_filament r1 red 10 u u u u u u -random_filament r2 red 10 u u u u u u -random_filament r3 red 10 u u u u u u -random_filament r4 red 10 u u u u u u -random_filament r5 red 10 u u u u u u -random_filament r6 red 10 u u u u u u -random_filament r7 red 10 u u u u u u -random_filament r8 red 10 u u u u u u -random_filament r9 red 10 u u u u u u +random_filament red:r0 10 u u u u u u +random_filament red:r1 10 u u u u u u +random_filament red:r2 10 u u u u u u +random_filament red:r3 10 u u u u u u +random_filament red:r4 10 u u u u u u +random_filament red:r5 10 u u u u u u +random_filament red:r6 10 u u u u u u +random_filament red:r7 10 u u u u u u +random_filament red:r8 10 u u u u u u +random_filament red:r9 10 u u u u u u start_filament_type green @@ -66,16 +66,16 @@ force_length -1 force_angle 10 10 10 end_filament_type -random_filament g0 green 80 u u u u u u -random_filament g1 green 80 u u u u u u -random_filament g2 green 80 u u u u u u -random_filament g3 green 80 u u u u u u -random_filament g4 green 80 u u u u u u -random_filament g5 green 80 u u u u u u -random_filament g6 green 80 u u u u u u -random_filament g7 green 80 u u u u u u -random_filament g8 green 80 u u u u u u -random_filament g9 green 80 u u u u u u +random_filament green:0 80 u u u u u u +random_filament green:1 80 u u u u u u +random_filament green:2 80 u u u u u u +random_filament green:3 80 u u u u u u +random_filament green:4 80 u u u u u u +random_filament green:5 80 u u u u u u +random_filament green:6 80 u u u u u u +random_filament green:7 80 u u u u u u +random_filament green:8 80 u u u u u u +random_filament green:9 80 u u u u u u end_file diff --git a/examples/S13_filaments/integration2d.txt b/examples/S13_filaments/integration2d.txt new file mode 100644 index 00000000..a800ae19 --- /dev/null +++ b/examples/S13_filaments/integration2d.txt @@ -0,0 +1,98 @@ +# Test of 2D filament integration methods + +#units nm us zJ + +#random_seed 2 +graphics opengl_good + +dim 2 +boundaries x -100 100 r +boundaries y -100 100 r + +#species red + +#difc red 3 + +#color red red + +time_start 0 +time_stop 50000 +time_step 0.03 + +#frame_thickness 1 + +#mol 1 red u u + +start_filament_type Euler +color blue +thickness 2 +polygon ve +force_arrows 1 red +kT 1 +dynamics Euler +standard_length 2 +standard_angle 0 +force_length 1 +force_angle 2 +mobility 10 +end_filament_type + +start_filament_type RK2 +color blue +thickness 2 +polygon ve +force_arrows 1 red +kT 1 +dynamics RK2 +standard_length 2 +standard_angle 0 +force_length 1 +force_angle 2 +mobility 10 +end_filament_type + +start_filament_type RK4 +color blue +thickness 2 +polygon ve +force_arrows 1 red +kT 1 +dynamics RK4 +standard_length 2 +standard_angle 0 +force_length 1 +force_angle 2 +mobility 10 +end_filament_type + + +random_filament Euler:e1 10 -50 -50 u +random_filament Euler:e2 20 0 -50 u +random_filament Euler:e3 30 50 -50 u + + +filament e1 copy_to RK2:t1 +filament e2 copy_to RK2:t2 +filament e3 copy_to RK2:t3 +filament t1 translate + 0 50 +filament t2 translate + 0 50 +filament t3 translate + 0 50 + + +filament e1 copy_to RK4:f1 +filament e2 copy_to RK4:f2 +filament e3 copy_to RK4:f3 +filament f1 translate + 0 100 +filament f2 translate + 0 100 +filament f3 translate + 0 100 + + +#filament_type template standard_angle 0.3 + +cmd B pause +#cmd E printFilament template fil1 xabf stdout + +end_file + + + diff --git a/source/Smoldyn/smolcmd.c b/source/Smoldyn/smolcmd.c index 42107e68..6bbd5089 100644 --- a/source/Smoldyn/smolcmd.c +++ b/source/Smoldyn/smolcmd.c @@ -3153,15 +3153,9 @@ enum CMDcode cmdprintFilament(simptr sim,cmdptr cmd,char *line2) { else scmdfprintf(cmd->cmds,fptr," %i len=%1.3g thick=%1.3g pos.=(%1.3g,%1.3g,%1.3g)->(%1.3g,%1.3g,%1.3g) angle=(%1.3g,%1.3g,%1.3g)",segment->index,segment->len,segment->thk,segment->xyzfront[0],segment->xyzfront[1],segment->xyzfront[2],segment->xyzback[0],segment->xyzback[1],segment->xyzback[2],segment->ypr[0],segment->ypr[1],segment->ypr[2]); } if(strchr(code,'a')) { - if(sim->dim==2) - scmdfprintf(cmd->cmds,fptr," A=(%1.3g,%1.3g;%1.3g,%1.3g)",segment->dcm[0],segment->dcm[1],segment->dcm[3],segment->dcm[4]); - else - scmdfprintf(cmd->cmds,fptr," A=(%1.3g,%1.3g,%1.3g;%1.3g,%1.3g,%1.3g;%1.3g,%1.3g,%1.3g)",segment->dcm[0],segment->dcm[1],segment->dcm[2],segment->dcm[3],segment->dcm[4],segment->dcm[5],segment->dcm[6],segment->dcm[7],segment->dcm[8]); } + scmdfprintf(cmd->cmds,fptr," a=(%1.3g,%1.3g,%1.3g,%1.3g)",segment->qrel[0],segment->qrel[1],segment->qrel[2],segment->qrel[3]); } if(strchr(code,'b')) { - if(sim->dim==2) - scmdfprintf(cmd->cmds,fptr," B=(%1.3g,%1.3g;%1.3g,%1.3g)",segment->adcm[0],segment->adcm[1],segment->adcm[3],segment->adcm[4]); - else - scmdfprintf(cmd->cmds,fptr," B=(%1.3g,%1.3g,%1.3g;%1.3g,%1.3g,%1.3g;%1.3g,%1.3g,%1.3g)",segment->adcm[0],segment->adcm[1],segment->adcm[2],segment->adcm[3],segment->adcm[4],segment->adcm[5],segment->adcm[6],segment->adcm[7],segment->adcm[8]); } + scmdfprintf(cmd->cmds,fptr," b=(%1.3g,%1.3g,%1.3g,%1.3g)",segment->qabs[0],segment->qabs[1],segment->qabs[2],segment->qabs[3]); } if(strchr(code,'f')) { if(sim->dim==2) scmdfprintf(cmd->cmds,fptr," F=(%1.3g,%1.3g)",fil->forces[seg][0],fil->forces[seg][1]); diff --git a/source/Smoldyn/smoldyn.h b/source/Smoldyn/smoldyn.h index 127a5a69..1705543b 100644 --- a/source/Smoldyn/smoldyn.h +++ b/source/Smoldyn/smoldyn.h @@ -119,7 +119,7 @@ typedef struct qstruct int n; int f; int b; -} * queue; +} * q_queue; #endif @@ -185,8 +185,8 @@ typedef struct cmdstruct typedef struct cmdsuperstruct { - queue cmd; // queue of normal run-time commands - queue cmdi; // queue of integer time commands + q_queue cmd; // queue of normal run-time commands + q_queue cmdi; // queue of integer time commands enum CMDcode (*cmdfn)(void*, cmdptr, char*); // function that runs commands void* cmdfnarg; // function argument (e.g. sim) int iter; // number of times integer commands have run @@ -773,8 +773,10 @@ typedef struct segmentstruct double len; // segment length double thk; // thickness of segment double ypr[3]; // relative ypr angles - double dcm[9]; // relative dcm - double adcm[9]; // absolute segment orientation + double qrel[4]; // relative rotation quaternion + double qabs[4]; // absolute rotation quaternion +// double dcm[9]; // relative dcm +// double adcm[9]; // absolute segment orientation } * segmentptr; typedef struct filamentstruct @@ -788,10 +790,13 @@ typedef struct filamentstruct double **nodes1; // nodes for RK dynamics (nseg+1) double **nodes2; // nodes for RK4 dynamics (nseg+1) double **nodesx; // old node locations (nseg+1) - double *roll1; // roll for RK dynamics (nseg) - double *roll2; // roll for RK dynamics (nseg) - double **forces; // list of forces on nodes (nseg+1) - double *torques; // list of segment torques (nseg) + double *roll1; // roll for RK dynamics (nseg) + double *roll2; // roll for RK dynamics (nseg) + double **forces; // forces on nodes (nseg+1) + double *torques; // list of segment torques (nseg) + double seg0up[3]; // segment 0 up vector + double seg0up1[3]; // segment 0 up vector for RK + double seg0up2[3]; // segment 0 up vector for RK4 struct filamentstruct* frontend; // what front attaches to struct filamentstruct* backend; // what back attaches to int maxbranch; // max branches off this filament @@ -1031,21 +1036,14 @@ typedef struct simstruct boxssptr boxs; // box superstructure compartssptr cmptss; // compartment superstructure - portssptr portss; // port superstructure - latticessptr latticess; // lattice superstructure - bngssptr bngss; // bionetget superstructure - filamentssptr filss; // filament superstructure - cmdssptr cmds; // command superstructure - graphicsssptr graphss; // graphics superstructure diffusefnptr diffusefn; // function for molecule diffusion - surfaceboundfnptr surfaceboundfn; // function for surface-bound molecules surfacecollisionsfnptr surfacecollisionsfn; // function for surface collisons assignmols2boxesfnptr assignmols2boxesfn; // function that assigns molecs to boxes diff --git a/source/Smoldyn/smoldynfuncs.h b/source/Smoldyn/smoldynfuncs.h index ceae468e..1bc8464f 100644 --- a/source/Smoldyn/smoldynfuncs.h +++ b/source/Smoldyn/smoldynfuncs.h @@ -401,6 +401,7 @@ int latticeruntimestep(simptr sim); // enumerated types // low level utilities +int filReadFilName(simptr sim,const char *str,filamenttypeptr *filtypeptr,filamentptr *filptr,char *filname); // memory management void filssfree(filamentssptr filss); diff --git a/source/Smoldyn/smolfilament.c b/source/Smoldyn/smolfilament.c index 2105b9f6..67913a7f 100644 --- a/source/Smoldyn/smolfilament.c +++ b/source/Smoldyn/smolfilament.c @@ -67,9 +67,10 @@ int filAddSegment(filamentptr fil,const double *x,double length,const double *an int filRemoveSegment(filamentptr fil,char endchar); void filTranslate(filamentptr fil,const double *vect,char func); void filRotateVertex(filamentptr fil,int seg,const double *angle,char endchar,char func); -void filLengthenSegment(filamentptr fil,int seg,double length,char endchar,char func); +int filLengthenSegment(filamentptr fil,int seg,double length,char endchar,char func); +int filChangeThickness(filamentptr fil,int seg,double thick,char func); void filReverseFilament(filamentptr fil); -int filCopyFilament(filamentptr filto,const filamentptr filfrom); +int filCopyFilament(filamentptr filto,const filamentptr filfrom,const filamenttypeptr filtype); // filament type int filtypeSetParam(filamenttypeptr filtype,const char *param,int index,double value); @@ -127,6 +128,55 @@ enum FilamentDynamics filstring2FD(const char *string) { /******************************************************************************/ +/* filReadFilName */ +int filReadFilName(simptr sim,const char *str,filamenttypeptr *filtypeptr,filamentptr *filptr,char *filname) { + int itct,ft,f; + char nm[STRCHAR],*fnm,*colon; + filamenttypeptr filtype; + + if(!str) return -1; // str missing + if(!sim->filss || !sim->filss->ntype) return -2; // no filaments + + itct=sscanf(str,"%s",nm); + if(itct!=1) return -3; // can't read str + colon=strchr(nm,':'); + + if(colon) { // format is filtype:fil + fnm=colon+1; + *colon='\0'; + if(filname) strcpy(filname,fnm); + ft=stringfind(sim->filss->ftnames,sim->filss->ntype,nm); + if(ft<0) return -4; // filtype not recognized + filtype=sim->filss->filtypes[ft]; + *filtypeptr=filtype; + f=stringfind(filtype->filnames,filtype->nfil,fnm); + if(f<0) return -5; // filament not recognized + *filptr=filtype->fillist[f]; } + + else if(*filtypeptr) { // known filtype, looking for fil + filtype=*filtypeptr; + fnm=nm; + if(filname) strcpy(filname,fnm); + f=stringfind(filtype->filnames,filtype->nfil,fnm); + if(f<0) return -5; // filament not recognized + *filptr=filtype->fillist[f]; } + + else { // unknown filtype + fnm=nm; + if(filname) strcpy(filname,fnm); + for(ft=0;ftfilss->ntype;ft++) { + filtype=sim->filss->filtypes[ft]; + f=stringfind(filtype->filnames,filtype->nfil,fnm); + if(f>=0) break; } + if(f>=0) { + *filtypeptr=filtype; + *filptr=filtype->fillist[f]; } + else + return -5; } // filament not recognized + + return 0; } + + /* filRandomLength */ double filRandomLength(const filamenttypeptr filtype,double thickness,double sigmamult) { double lsigma; @@ -164,7 +214,7 @@ double *filRandomAngle(const filamenttypeptr filtype,int dim,int n,double thickn angle[0]=thetarandCCD(); angle[1]=unirandCOD(0,2*PI); angle[2]=unirandCOD(0,2*PI); - Sph_Eax2Xyz(angle,angle); } + Sph_Eax2Ypr(angle,angle); } else { // 2D first segment angle[0]=unirandOCD(-PI,PI); angle[1]=angle[2]=0; } @@ -225,7 +275,7 @@ double filBendEnergy(const filamentptr fil,int seg1,int seg2) { // filBendTorque void filBendTorque(const filamentptr fil,int node,double *torque) { filamenttypeptr filtype; - double *kypr,*stdypr,*ypr,cphi,sphi,ctht,stht,deltaypr[3]; + double *kypr,*stdypr,*ypr,cphi,sphi,ctht,stht,deltaypr[3],reltorque[3]; int dim; if(node<=0 || node>=fil->nseg) { @@ -252,11 +302,14 @@ void filBendTorque(const filamentptr fil,int node,double *torque) { deltaypr[1]=kypr[1]*(ypr[1]-stdypr[1]); deltaypr[2]=kypr[2]*(ypr[2]-stdypr[2]); - torque[0]=-deltaypr[2]*ctht*cphi+deltaypr[1]*sphi; // compute torque in segment reference frame - torque[1]=-deltaypr[1]*cphi+deltaypr[2]*sphi; - torque[2]=-deltaypr[0]-deltaypr[2]*cphi*stht; + reltorque[0]=-deltaypr[2]*ctht*cphi+deltaypr[1]*sphi; // compute torque in segment reference frame + reltorque[1]=-deltaypr[1]*cphi-deltaypr[2]*sphi; + reltorque[2]=-deltaypr[0]+deltaypr[2]*cphi*stht; - Sph_DcmtxCart(fil->segments[node-1]->adcm,torque,torque); } // convert torque to system reference frame + Sph_QtniRotate(fil->segments[node-1]->qabs,reltorque,torque); // convert torque to system reference frame +// Sph_DcmtxCart(fil->segments[node-1]->adcm,reltorque,torque); +// printf("*** filBendTorque. node=%i reltorque=(%1.3g %1.3g %1.3g) abstorque=(%1.3g %1.3g %1.3g)\n",node,reltorque[0],reltorque[1],reltorque[2],torque[0],torque[1],torque[2]);//?? + } return; } @@ -277,8 +330,10 @@ segmentptr segmentalloc() { segment->len=0; segment->thk=1; segment->ypr[0]=segment->ypr[1]=segment->ypr[2]=0; - Sph_One2Dcm(segment->dcm); - Sph_One2Dcm(segment->adcm); + Sph_One2Qtn(segment->qrel); + Sph_One2Qtn(segment->qabs); +// Sph_One2Dcm(segment->dcm); +// Sph_One2Dcm(segment->adcm); return segment; failure: @@ -316,6 +371,8 @@ filamentptr filalloc(filamentptr fil,int maxseg,int maxbranch,int maxmonomer) { fil->roll2=NULL; fil->forces=NULL; fil->torques=NULL; + fil->seg0up[0]=fil->seg0up[1]=0; + fil->seg0up[2]=1; fil->frontend=NULL; fil->backend=NULL; fil->maxbranch=0; @@ -651,41 +708,53 @@ void filoutput(const filamentptr fil) { sim=NULL; dim=3; } - simLog(sim,2," Filament: %s\n",fil->filname); - simLog(sim,1," type: %s\n",fil->filtype?fil->filtype->ftname:"None (assuming dim=3)"); - simLog(sim,1," allocated segments: %i\n",fil->maxseg); - simLog(sim,2," number of segments: %i\n",fil->nseg); + simLog(sim,2," Filament: %s\n",fil->filname); + simLog(sim,1," type: %s\n",fil->filtype?fil->filtype->ftname:"None (assuming dim=3)"); + simLog(sim,1," allocated segments: %i\n",fil->maxseg); + simLog(sim,2," number of segments: %i\n",fil->nseg); + if(dim==3) + simLog(sim,2," segment 0 up vector: (%1.3g %1.3g %1.3g)\n",fil->seg0up[0],fil->seg0up[1],fil->seg0up[2]); - simLog(sim,2," segment, length, thickness, front position, relative angle:\n"); - for(seg=0;segnseg;seg++) { - segment=fil->segments[seg]; - if(dim==2) - simLog(sim,seg>5?1:2," %i length=%1.3g, thick=%1.3g, front pos.=(%1.3g %1.3g), rel. angle=%1.3g\n",segment->index,segment->len,segment->thk,segment->xyzfront[0],segment->xyzfront[1],segment->ypr[0]); - else - simLog(sim,seg>5?1:2," %i length=%1.3g, thick=%1.3g, front pos.=(%1.3g %1.3g %1.3g), rel. angle=(%1.3g %1.3g %1.3g)\n",segment->index,segment->len,segment->thk,segment->xyzfront[0],segment->xyzfront[1],segment->xyzfront[2],segment->ypr[0],segment->ypr[1],segment->ypr[2]); } + simLog(sim,2," segments:\n"); + if(dim==2) { + for(seg=0;segnseg;seg++) { + segment=fil->segments[seg]; + simLog(sim,seg>5?1:2," %i length=%1.3g, thick=%1.3g, front pos.=(%1.3g %1.3g), rel. angle=%1.3g",segment->index,segment->len,segment->thk,segment->xyzfront[0],segment->xyzfront[1],segment->ypr[0]); + simLog(sim,1,", back pos.=(%1.3g %1.3g), qrel=(%1.3g %1.3g %1.3g %1.3g), qabs=(%1.3g %1.3g %1.3g %1.3g)",segment->xyzback[0],segment->xyzback[1],segment->qrel[0],segment->qrel[1],segment->qrel[2],segment->qrel[3],segment->qabs[0],segment->qabs[1],segment->qabs[2],segment->qabs[3]); + simLog(sim,seg>5?1:2,"\n"); } + segment=fil->segments[fil->nseg-1]; + simLog(sim,2," back pos.=(%1.3g %1.3g)\n",segment->xyzback[0],segment->xyzback[1]); } + else { + for(seg=0;segnseg;seg++) { + segment=fil->segments[seg]; + simLog(sim,seg>5?1:2," %i length=%1.3g, thick=%1.3g, front pos.=(%1.3g %1.3g %1.3g), rel. angle=(%1.3g %1.3g %1.3g)",segment->index,segment->len,segment->thk,segment->xyzfront[0],segment->xyzfront[1],segment->xyzfront[2],segment->ypr[0],segment->ypr[1],segment->ypr[2]); + simLog(sim,1,", back pos.=(%1.3g %1.3g %1.3g), qrel=(%1.3g %1.3g %1.3g %1.3g), qabs=(%1.3g %1.3g %1.3g %1.3g)",segment->xyzback[0],segment->xyzback[1],segment->xyzback[2],segment->qrel[0],segment->qrel[1],segment->qrel[2],segment->qrel[3],segment->qabs[0],segment->qabs[1],segment->qabs[2],segment->qabs[3]); + simLog(sim,seg>5?1:2,"\n"); } + segment=fil->segments[fil->nseg-1]; + simLog(sim,2," back pos.=(%1.3g %1.3g %1.3g)\n",segment->xyzback[0],segment->xyzback[1],segment->xyzback[2]); } if(fil->frontend) - simLog(sim,2," front branched from: %s\n",fil->frontend->filname); + simLog(sim,2," front branched from: %s\n",fil->frontend->filname); if(fil->backend) - simLog(sim,2," back branched from: %s\n",fil->backend->filname); + simLog(sim,2," back branched from: %s\n",fil->backend->filname); - simLog(sim,1," allocated branches: %i\n",fil->maxbranch); - simLog(sim,fil->nbranch>0?2:1," number of branches: %i\n",fil->nbranch); + simLog(sim,1," allocated branches: %i\n",fil->maxbranch); + simLog(sim,fil->nbranch>0?2:1," number of branches: %i\n",fil->nbranch); for(br=0;brnbranch;br++) - simLog(sim,2," %s at %i\n",fil->branches[br]->filname,fil->branchspots[br]); + simLog(sim,2," %s at %i\n",fil->branches[br]->filname,fil->branchspots[br]); - simLog(sim,1," monomer codes: %i of %i allocated,",fil->nmonomer,fil->maxmonomer); - simLog(sim,1," starting index: %i\n",fil->frontmonomer); + simLog(sim,1," monomer codes: %i of %i allocated,",fil->nmonomer,fil->maxmonomer); + simLog(sim,1," starting index: %i\n",fil->frontmonomer); if(fil->nmonomer) { - simLog(sim,2," monomer code: "); + simLog(sim,2," monomer code: "); for(mn=0;mnnmonomer;mn++) simLog(sim,2,"%c",fil->monomers[mn]); simLog(sim,2,"\n"); } if(fil->filtype->klen>0) - simLog(sim,2," stretching energy: %g|E\n",filStretchEnergy(fil,-1,-1)); + simLog(sim,2," stretching energy: %g|E\n",filStretchEnergy(fil,-1,-1)); if(fil->filtype->kypr[0]>0 || fil->filtype->kypr[1]>0 || fil->filtype->kypr[2]>0) - simLog(sim,2," bending energy: %g|E\n",filBendEnergy(fil,-1,-1)); + simLog(sim,2," bending energy: %g|E\n",filBendEnergy(fil,-1,-1)); return; } @@ -831,8 +900,8 @@ void filArrayShift(filamentptr fil,int shift) { double *newnode,*newforce,newtorque; if(shift>0) { - if(fil->nseg+shift>=fil->maxseg) - simLog(fil->filtype->filss->sim,10,"BUG in filArrayShift. Memory overwrite."); + if(fil->nseg+shift>fil->maxseg) + simLog(fil->filtype->filss->sim,10,"BUG in filArrayShift. Memory overwrite.\n"); for(i=fil->nseg+shift-1;i>=shift;i--) { // i is new index newsegment = fil->segments[i]; newtorque = fil->torques[i]; @@ -846,7 +915,10 @@ void filArrayShift(filamentptr fil,int shift) { fil->nodes[i]=fil->nodes[i-shift]; fil->forces[i]=fil->forces[i-shift]; fil->nodes[i-shift]=newnode; - fil->forces[i-shift]=newforce; }} + fil->forces[i-shift]=newforce; } + for(i=0;isegments[i]->xyzfront=fil->nodes[i]; + fil->segments[i]->xyzback=fil->nodes[i+1]; }} else if(shift<0) { shift=-shift; // now shift is positive @@ -863,7 +935,10 @@ void filArrayShift(filamentptr fil,int shift) { fil->nodes[i]=fil->nodes[i+shift]; fil->forces[i]=fil->forces[i+shift]; fil->nodes[i+shift]=newnode; - fil->forces[i+shift]=newforce; }} + fil->forces[i+shift]=newforce; } + for(i=fil->nseg-shift;inseg;i++) { + fil->segments[i]->xyzfront=fil->nodes[i]; + fil->segments[i]->xyzback=fil->nodes[i+1]; }} for(i=0;inseg;i++) fil->segments[i]->index=i; @@ -875,12 +950,7 @@ void filArrayShift(filamentptr fil,int shift) { int filAddSegment(filamentptr fil,const double *x,double length,const double *angle,double thickness,char endchar) { int seg; segmentptr segment,segmentm1,segmentp1; - double dcm[9]; - -// printf("filAddSegment. end=%c",endchar);//?? -// if(x) printf(" x=(%g,%g,%g)",x[0],x[1],x[2]);//?? -// else printf(" x=NULL");//?? -// printf(" length=%g, angle=(%g,%g,%g), thickness=%g\n",length,angle[0],angle[1],angle[2],thickness);//?? + double zero[3]={0,0,0}; if(fil->nseg==fil->maxseg) { fil=filalloc(fil,fil->maxseg*2+1,0,0); @@ -891,19 +961,22 @@ int filAddSegment(filamentptr fil,const double *x,double length,const double *an segment=fil->segments[seg]; segment->len=length; segment->thk=thickness; - segment->xyzfront=fil->nodes[seg]; - segment->xyzback=fil->nodes[seg+1]; - Sph_Xyz2Xyz(angle,segment->ypr); // ypr = angle - Sph_Xyz2Dcm(angle,segment->dcm); // A = Dcm(angle) + Sph_Ypr2Ypr(angle,segment->ypr); // ypr = angle +// Sph_Ypr2Dcm(angle,segment->dcm); // A = Dcm(angle) + Sph_Ypr2Qtn(angle,segment->qrel); // a = Qtn(angle) if(seg==0) { segment->xyzfront[0]=x[0]; // x_i = input value segment->xyzfront[1]=x[1]; segment->xyzfront[2]=x[2]; - Sph_Dcm2Dcm(segment->dcm,segment->adcm); } // B = A +// Sph_Dcm2Dcm(segment->dcm,segment->adcm); // B = A + Sph_Qtn2Qtn(segment->qrel,segment->qabs); // b = a + Sph_QtniRotateUnitz(segment->qabs,fil->seg0up,zero,1); } // rotate segment z to lab frame else { segmentm1=fil->segments[seg-1]; - Sph_DcmxDcm(segment->dcm,segmentm1->adcm,segment->adcm); } // B_i = A_i . B_{i-1} - Sph_DcmtxUnit(segment->adcm,'x',segment->xyzback,segment->xyzfront,segment->len); // x_{i+1} = x_i + l_i * BT_i . xhat +// Sph_DcmxDcm(segment->dcm,segmentm1->adcm,segment->adcm); // B_i = A_i . B_{i-1} + Sph_QtnxQtn(segmentm1->qabs,segment->qrel,segment->qabs); } // b_i = b_{i-1} . a_i +// Sph_DcmtxUnit(segment->adcm,'x',segment->xyzback,segment->xyzfront,segment->len); // x_{i+1} = x_i + l_i * BT_i . xhat + Sph_QtniRotateUnitx(segment->qabs,segment->xyzback,segment->xyzfront,segment->len); // x_{i+1} = x_i + l_i * b.xhat.bT fil->nseg++; } else { @@ -912,21 +985,26 @@ int filAddSegment(filamentptr fil,const double *x,double length,const double *an segment=fil->segments[seg]; segment->len=length; segment->thk=thickness; - segment->xyzfront=fil->nodes[0]; - segment->xyzback=fil->nodes[1]; if(fil->nseg==0) { - Sph_Xyz2Dcmt(angle,segment->adcm); // B_0 = Dcmt(angle) +// Sph_Ypr2Dcmt(angle,segment->adcm); // B_0 = Dcmt(angle) + Sph_Ypr2Qtni(angle,segment->qabs); // b_0 segment->xyzback[0]=x[0]; // back of segment = input value segment->xyzback[1]=x[1]; segment->xyzback[2]=x[2]; } else { segmentp1=fil->segments[1]; - Sph_Xyz2Dcm(angle,dcm); - Sph_DcmtxDcm(dcm,segmentp1->adcm,segment->adcm); - Sph_Dcm2Dcm(dcm,segmentp1->dcm); } - Sph_Dcm2Dcm(segment->adcm,segment->dcm); // A_i = B_i - Sph_Dcm2Xyz(segment->dcm,segment->ypr); // a_0 = Xyz(B_0) - Sph_DcmtxUnit(segment->adcm,'x',segment->xyzfront,segment->xyzback,-segment->len); // x_i = x_{i+1} - l_i * BT_i . xhat + Sph_Ypr2Ypr(angle,segmentp1->ypr); +// Sph_Ypr2Dcm(angle,segmentp1->dcm); + Sph_Ypr2Qtn(angle,segmentp1->qrel); +// Sph_DcmtxDcm(segmentp1->dcm,segmentp1->adcm,segment->adcm); // B_i = A_{i+1}T . B_{i+1} + Sph_QtnxQtni(segmentp1->qabs,segmentp1->qrel,segment->qabs); } // b_i = b_{i+1} . a_{i+1}* +// Sph_Dcm2Dcm(segment->adcm,segment->dcm); // A_i = B_i + Sph_Qtn2Qtn(segment->qabs,segment->qrel); // a_i = b_i +// Sph_Dcm2Ypr(segment->dcm,segment->ypr); // a_0 = Ypr(B_0) + Sph_Qtn2Ypr(segment->qrel,segment->ypr); // a_0 = Ypr(B_0) +// Sph_DcmtxUnit(segment->adcm,'x',segment->xyzfront,segment->xyzback,-segment->len); // x_i = x_{i+1} - l_i * BT_i . xhat + Sph_QtniRotateUnitx(segment->qabs,segment->xyzfront,segment->xyzback,-segment->len); + Sph_QtniRotateUnitz(segment->qabs,fil->seg0up,zero,1); fil->nseg++; } return 0; } @@ -962,13 +1040,13 @@ int filAddRandomSegments(filamentptr fil,int number,const char *xstr,const char else return 2; if(dim==2) { - angle[0]=unirandCCD(-PI,PI); + angle[0]=unirandOCD(-PI,PI); angle[1]=angle[2]=0; } else { angle[0]=thetarandCCD(); angle[1]=unirandCOD(0,2*PI); angle[2]=unirandCOD(0,2*PI); - Sph_Eax2Xyz(angle,angle); } + Sph_Eax2Ypr(angle,angle); } if(!strcmp(phistr,"u")); else if(strmathsscanf(phistr,"%mlg|",varnames,varvalues,nvar,&f1)==1) angle[0]=f1; else return 3; @@ -998,6 +1076,7 @@ int filAddRandomSegments(filamentptr fil,int number,const char *xstr,const char int filRemoveSegment(filamentptr fil,char endchar) { int seg; segmentptr segment; + double zero[3]={0,0,0}; if(fil->nseg==0) return -1; @@ -1006,8 +1085,11 @@ int filRemoveSegment(filamentptr fil,char endchar) { else { seg=1; // new front segment segment=fil->segments[seg]; - Sph_Dcm2Dcm(segment->adcm,segment->dcm); - Sph_Dcm2Xyz(segment->dcm,segment->ypr); +// Sph_Dcm2Dcm(segment->adcm,segment->dcm); + Sph_Qtn2Qtn(segment->qabs,segment->qrel); +// Sph_Dcm2Ypr(segment->dcm,segment->ypr); + Sph_Qtn2Ypr(segment->qrel,segment->ypr); + Sph_QtniRotateUnitz(segment->qrel,fil->seg0up,zero,1); filArrayShift(fil,-1); fil->nseg--; } @@ -1020,12 +1102,13 @@ void filTranslate(filamentptr fil,const double *vect,char func) { double shift[3],*node; segmentptr segment; - seg=0; + if(fil->nseg==0) return; + if(func=='=') { - segment=fil->segments[seg]; - shift[0]=segment->xyzfront[0]-vect[0]; - shift[1]=segment->xyzfront[1]-vect[1]; - shift[2]=segment->xyzfront[2]-vect[2];} + segment=fil->segments[0]; + shift[0]=vect[0]-segment->xyzfront[0]; + shift[1]=vect[1]-segment->xyzfront[1]; + shift[2]=vect[2]-segment->xyzfront[2]; } else if(func=='-') { shift[0]=-vect[0]; shift[1]=-vect[1]; @@ -1045,20 +1128,21 @@ void filTranslate(filamentptr fil,const double *vect,char func) { /* filLengthenSegment */ -void filLengthenSegment(filamentptr fil,int seg,double length,char endchar,char func) { +int filLengthenSegment(filamentptr fil,int seg,double length,char endchar,char func) { int i; - double lenold,lendelta,xdelta[3],**nodes; + double lenold,lendelta,xdelta[3],**nodes,zero[3]={0,0,0}; segmentptr segment; segment=fil->segments[seg]; lenold=segment->len; if(func=='=') lendelta=length-lenold; else if(func=='+') lendelta=length; - else lendelta=lenold-length; + else lendelta=-length; + + if(lenold+lendelta<=0) return 1; - xdelta[0]=lendelta*segment->adcm[0]; // same as SphDcmtxUnit(segment->adcm,'x',xdelta,NULL,lendelta); - xdelta[1]=lendelta*segment->adcm[1]; - xdelta[2]=lendelta*segment->adcm[2]; +// Sph_DcmtxUnit(segment->adcm,'x',xdelta,NULL,lendelta); // Delta Delta x = Delta l * B_iT xhat + Sph_QtniRotateUnitx(segment->qabs,xdelta,zero,lendelta); // rotate to lab frame. Delta Delta x = Delta l * b_i xhat b_iT nodes=fil->nodes; if(endchar=='b') { for(i=seg+1;i<=fil->nseg;i++) { @@ -1071,55 +1155,86 @@ void filLengthenSegment(filamentptr fil,int seg,double length,char endchar,char nodes[i][1]-=xdelta[1]; nodes[i][2]-=xdelta[2]; }} - return; } + return 0; } + + +int filChangeThickness(filamentptr fil,int seg,double thick,char func) { + segmentptr segment; + double thkold,thkdelta; + + segment=fil->segments[seg]; + thkold=segment->thk; + if(func=='=') thkdelta=thick-thkold; + else if(func=='+') thkdelta=thick; + else thkdelta=-thick; + + if(thkold+thkdelta<0) return 1; + segment->thk+=thkdelta; + + return 0; } /* filRotateVertex */ void filRotateVertex(filamentptr fil,int seg,const double *angle,char endchar,char func) { int i; - double dcmdelta[9]; +//double dcmdelta[9]; + double qtndelta[4],zero[3]={0,0,0}; segmentptr segment,segmentm1,segmentp1; segment=fil->segments[seg]; - Sph_Xyz2Dcm(angle,dcmdelta); - if(func=='=') Sph_Dcm2Dcm(dcmdelta,segment->dcm); - else if(func=='+') Sph_DcmxDcm(dcmdelta,segment->dcm,segment->dcm); - else Sph_DcmtxDcm(dcmdelta,segment->dcm,segment->dcm); +//Sph_Ypr2Dcm(angle,dcmdelta); + Sph_Ypr2Qtn(angle,qtndelta); + if(func=='=') { +// Sph_Dcm2Dcm(dcmdelta,segment->dcm); + Sph_Qtn2Qtn(qtndelta,segment->qrel); } + else if(func=='+') { +// Sph_DcmxDcm(dcmdelta,segment->dcm,segment->dcm); + Sph_QtnxQtn(segment->qrel,qtndelta,segment->qrel); } + else { +// Sph_DcmtxDcm(dcmdelta,segment->dcm,segment->dcm); + Sph_QtnxQtni(segment->qrel,qtndelta,segment->qrel); } if(endchar=='b') { for(i=seg;inseg;i++) { segment=fil->segments[i]; - if(i==0) Sph_Dcm2Dcm(segment->dcm,segment->adcm); + if(i==0) { +// Sph_Dcm2Dcm(segment->dcm,segment->adcm); // B_0 = A_0 + Sph_Qtn2Qtn(segment->qrel,segment->qabs); // b_0 = a_0 + Sph_QtniRotateUnitz(segment->qabs,fil->seg0up,zero,1); } else { segmentm1=fil->segments[i-1]; - Sph_DcmxDcm(segment->dcm,segmentm1->adcm,segment->adcm); } - Sph_Dcm2Xyz(segment->dcm,segment->ypr); - segment->xyzback[0]=segment->xyzfront[0]+segment->len*segment->adcm[0]; - segment->xyzback[1]=segment->xyzfront[1]+segment->len*segment->adcm[1]; - segment->xyzback[2]=segment->xyzfront[2]+segment->len*segment->adcm[2]; }} +// Sph_DcmxDcm(segment->dcm,segmentm1->adcm,segment->adcm); // B_i = A_i B_{i-1} + Sph_QtnxQtn(segmentm1->qabs,segment->qrel,segment->qabs); } // b_i = b_{i-1} a_i +// Sph_Dcm2Ypr(segment->dcm,segment->ypr); + Sph_Qtn2Ypr(segment->qrel,segment->ypr); +// Sph_DcmtxUnit(segment->adcm,'x',segment->xyzback,segment->xyzfront,segment->len); // x_{i+1} = x_i + l_i B_iT xhat + Sph_QtniRotateUnitx(segment->qabs,segment->xyzback,segment->xyzfront,segment->len); + }} else { for(i=seg;i>=0;i--) { segment=fil->segments[i]; if(i==fil->nseg-1); else { segmentp1=fil->segments[i+1]; - Sph_DcmtxDcm(segmentp1->dcm,segmentp1->adcm,segment->adcm); } - Sph_Dcm2Xyz(segment->dcm,segment->ypr); - segment->xyzfront[0]=segment->xyzback[0]-segment->len*segment->adcm[0]; - segment->xyzfront[1]=segment->xyzback[1]-segment->len*segment->adcm[1]; - segment->xyzfront[2]=segment->xyzback[2]-segment->len*segment->adcm[2]; }} +// Sph_DcmtxDcm(segmentp1->dcm,segmentp1->adcm,segment->adcm); // B_i = A_{i+1}T B_i + Sph_QtnxQtni(segmentp1->qabs,segmentp1->qrel,segment->qabs); } // b_i = b_i a_{i+1)T +// Sph_Dcm2Ypr(segment->dcm,segment->ypr); + Sph_Qtn2Ypr(segment->qrel,segment->ypr); +// Sph_DcmtxUnit(segment->adcm,'x',segment->xyzfront,segment->xyzback,-segment->len); // x_i = x_{i+1} - l_i B_{i+1}T xhat + Sph_QtniRotateUnitx(segment->qabs,segment->xyzfront,segment->xyzback,-segment->len); } + Sph_QtniRotateUnitz(segment->qabs,fil->seg0up,zero,1); } return; } /* filCopyFilament */ -int filCopyFilament(filamentptr filto,const filamentptr filfrom) { +int filCopyFilament(filamentptr filto,const filamentptr filfrom,const filamenttypeptr filtype) { int i; - segmentptr segmentto,segmentfrom; + segmentptr segmentfrom; if(!filto || !filfrom) return 2; - filto->filtype=filfrom->filtype; + filto->filtype=filtype?filtype:filfrom->filtype; filto->nseg=0; filto->nbranch=0; @@ -1129,24 +1244,16 @@ int filCopyFilament(filamentptr filto,const filamentptr filfrom) { if(!filto) return 1; for(i=0;inseg;i++) { - segmentto=filto->segments[i]; segmentfrom=filfrom->segments[i]; - segmentto->index=segmentfrom->index; - copyVD(segmentfrom->xyzfront,segmentto->xyzfront,3); - copyVD(filfrom->forces[i],filto->forces[i],3); - segmentto->len=segmentfrom->len; - segmentto->thk=segmentfrom->thk; - copyVD(segmentfrom->ypr,segmentto->ypr,3); - copyVD(segmentfrom->dcm,segmentto->dcm,9); - copyVD(segmentfrom->adcm,segmentto->adcm,9); } - if(filfrom->nseg>0) { - copyVD(segmentfrom->xyzback,segmentto->xyzback,3); - i=filfrom->nseg; - copyVD(filfrom->forces[i],filto->forces[i],3); } - filto->nseg=filfrom->nseg; + filAddSegment(filto,segmentfrom->xyzfront,segmentfrom->len,segmentfrom->ypr,segmentfrom->thk,'b'); } + + filto->seg0up[0]=filfrom->seg0up[0]; + filto->seg0up[1]=filfrom->seg0up[1]; + filto->seg0up[2]=filfrom->seg0up[2]; filto->frontend=filfrom->frontend; filto->backend=filfrom->backend; + for(i=0;inbranch;i++) { filto->branchspots[i]=filfrom->branchspots[i]; filto->branches[i]=filfrom->branches[i]; } @@ -1179,6 +1286,12 @@ filamentptr filAddFilament(filamenttypeptr filtype,const char *filname) { strncpy(filtype->filnames[f],filname,STRCHAR-1); filtype->filnames[f][STRCHAR-1]='\0'; fil=filtype->fillist[f]; + fil->filtype=filtype; + fil->nseg=0; + fil->frontend=fil->backend=NULL; + fil->nbranch=0; + fil->nmonomer=0; + fil->frontmonomer=0; filsetcondition(filtype->filss,SClists,0); return fil; } @@ -1413,7 +1526,7 @@ filamenttypeptr filtypereadstring(simptr sim,ParseFilePtr pfp,filamenttypeptr fi enum DrawMode dm; enum FilamentDynamics fd; - printf("%s\n",word);//?? debug + printf("%s %s\n",word,line2);//?? debug dim=sim->dim; varnames=sim->varnames; @@ -1588,13 +1701,12 @@ filamentptr filreadstring(simptr sim,ParseFilePtr pfp,filamentptr fil,filamentty int nvar,dim; filamentptr fil2; - filamentssptr filss; - int itct,er,i1; + filamenttypeptr filtype2; + int itct,er,i1,seg; char nm[STRCHAR],nm1[STRCHAR],endchar,str1[STRCHAR],str2[STRCHAR],str3[STRCHAR],str4[STRCHAR],str5[STRCHAR],str6[STRCHAR],symbol; double fltv1[9],length,angle[3],thick; - printf("%s\n",word);//?? debug - filss=sim->filss; + printf("%s %s\n",word,line2);//?? debug dim=sim->dim; varnames=sim->varnames; @@ -1602,14 +1714,22 @@ filamentptr filreadstring(simptr sim,ParseFilePtr pfp,filamentptr fil,filamentty nvar=sim->nvar; if(!strcmp(word,"name")) { - itct=sscanf(line2,"%s %s",nm1,nm); - CHECKS(itct==2,"correct format: filament_type filament_name"); - i1=stringfind(filss->ftnames,filss->ntype,nm1); - CHECKS(i1>=0,"filament type not recognized"); - filtype=filss->filtypes[i1]; - fil=filAddFilament(filtype,nm); - CHECKS(fil,"failed to add filament"); - CHECKS(!strnword(line2,3),"unexpected text following filament name"); } + itct=sscanf(line2,"%s",nm); + CHECKS(itct==1,"correct format: filament_type:filament_name"); + filtype2=filtype; + fil2=NULL; + er=filReadFilName(sim,nm,&filtype2,&fil2,nm1); + CHECKS(er!=-1 && er!=-3,"correct format: filament_type:filament_name"); + CHECKS(er!=-2,"filaments have not been defined yet"); + CHECKS(er!=-4,"filament type not recognized"); + CHECKS(filtype2,"missing filament type. Format: filament_type:filament_name"); + filtype=filtype2; + if(er==-5) { + fil=filAddFilament(filtype,nm1); + CHECKS(fil,"failed to add filament"); } + else + fil=fil2; + CHECKS(!strnword(line2,2),"unexpected text following filament name"); } else if(!strcmp(word,"first_segment")) { // first_segment CHECKS(fil,"need to enter filament name before first_segment"); @@ -1617,17 +1737,17 @@ filamentptr filreadstring(simptr sim,ParseFilePtr pfp,filamentptr fil,filamentty if(dim==2) { itct=strmathsscanf(line2,"%mlg|L %mlg|L %mlg|L %mlg|",varnames,varvalues,nvar,&fltv1[0],&fltv1[1],&length,&angle[0]); fltv1[2]=angle[1]=angle[2]=0; - CHECKS(itct==4,"first_segment format: x y length angle0 [thickness]"); } + CHECKS(itct==4,"first_segment format: x y length angle [thickness]"); } else { itct=strmathsscanf(line2,"%mlg|L %mlg|L %mlg|L %mlg|L %mlg| %mlg| %mlg|",varnames,varvalues,nvar,&fltv1[0],&fltv1[1],&fltv1[2],&length,&angle[0],&angle[1],&angle[2]); - CHECKS(itct==7,"first_segment format: x y z length angle0 angle1 angle2 [thickness]"); } + CHECKS(itct==7,"first_segment format: x y z length yaw pitch roll [thickness]"); } CHECKS(length>0,"length needs to be >0"); line2=strnword(line2,itct+1); thick=1; if(line2) { itct=strmathsscanf(line2,"%mlg|L",varnames,varvalues,nvar,&thick); - if(dim==2) {CHECKS(itct==1,"first_segment format: x y length angle0 [thickness]");} - else {CHECKS(itct==1,"first_segment format: x y z length angle0 angle1 angle2 [thickness]");} + if(dim==2) {CHECKS(itct==1,"first_segment format: x y length angle [thickness]");} + else {CHECKS(itct==1,"first_segment format: x y z length yaw pitch roll [thickness]");} CHECKS(thick>0,"thickness needs to be >0"); line2=strnword(line2,2); } er=filAddSegment(fil,fltv1,length,angle,thick,'b'); @@ -1713,23 +1833,90 @@ filamentptr filreadstring(simptr sim,ParseFilePtr pfp,filamentptr fil,filamentty else if(!strcmp(word,"translate")) { // translate CHECKS(fil,"need to enter filament name before translate"); - itct=strmathsscanf(line2,"%c %mlg|L %mlg|L %mlg|L",varnames,varvalues,nvar,&symbol,&fltv1[0],&fltv1[1],&fltv1[2]); - CHECKS(itct==4,"translate format: symbol x y z"); + if(dim==2) { + itct=strmathsscanf(line2,"%c %mlg|L %mlg|L",varnames,varvalues,nvar,&symbol,&fltv1[0],&fltv1[1]); + CHECKS(itct==3,"translate format: symbol x y"); + fltv1[2]=0; + line2=strnword(line2,4); } + else { + itct=strmathsscanf(line2,"%c %mlg|L %mlg|L %mlg|L",varnames,varvalues,nvar,&symbol,&fltv1[0],&fltv1[1],&fltv1[2]); + CHECKS(itct==4,"translate format: symbol x y z"); + line2=strnword(line2,5); } CHECKS(symbol=='=' || symbol=='+' || symbol=='-',"symbol needs to be '=', '+', or '-'"); filTranslate(fil,fltv1,symbol); - CHECKS(!strnword(line2,5),"unexpected text following translate"); } + CHECKS(!line2,"unexpected text following translate"); } + + else if(!strcmp(word,"modify_segment")) { // modify_segment + CHECKS(fil,"need to enter filament name before modify_segment"); + itct=strmathsscanf(line2,"%mi| %s %c",varnames,varvalues,nvar,&seg,nm,&symbol); + CHECKS(itct==3,"modify_segment format: segment length/angle/thickness +/-/= value [end]"); + CHECKS(seg>=0 && segnseg,"segment number is not within filament"); + line2=strnword(line2,4); + CHECKS(line2,"modify_segment format: segment length/angle/thickness +/-/= value [end]"); + if(!strcmp(nm,"length")) { + CHECKS(symbol=='=' || symbol=='+' || symbol=='-',"symbol needs to be '=', '+', or '-'"); + itct=strmathsscanf(line2,"%mlg|L",varnames,varvalues,nvar,&length); + CHECKS(itct==1,"modify_segment format: segment length/angle/thickness +/-/= value [end]"); + line2=strnword(line2,2); + endchar='b'; + if(line2) { + itct=sscanf(line2,"%s",nm1); + CHECKS(itct==1,"modify_segment format: segment length/angle/thickness +/-/= value [end]"); + line2=strnword(line2,2); + if(nm1[0]=='B' || nm1[0]=='b') endchar='b'; + else if(nm1[0]=='F' || nm1[0]=='f') endchar='f'; + else CHECKS(0,"end needs to be 'back' or 'front'"); } + er=filLengthenSegment(fil,seg,length,endchar,symbol); + CHECKS(!er,"Error while trying to change segment length - new length must be positive"); } + else if(!strcmp(nm,"angle")) { + CHECKS(symbol=='=' || symbol=='+' || symbol=='-',"symbol needs to be '=', '+', or '-'"); + angle[0]=angle[1]=angle[2]=0; + if(dim==2) { + itct=strmathsscanf(line2,"%mlg|",varnames,varvalues,nvar,&angle[0]); + CHECKS(itct==1,"modify_segment format: segment length/angle/thickness +/-/= value [end]"); + line2=strnword(line2,2); } + else { + itct=strmathsscanf(line2,"%mlg| %mlg| %mlg|",varnames,varvalues,nvar,&angle[0],&angle[1],&angle[2]); + CHECKS(itct==3,"insufficient angle values"); + line2=strnword(line2,4); } + endchar='b'; + if(line2) { + itct=sscanf(line2,"%s",nm1); + CHECKS(itct==1,"modify_segment format: segment length/angle/thickness +/-/= value [end]"); + line2=strnword(line2,2); + if(nm1[0]=='B' || nm1[0]=='b') endchar='b'; + else if(nm1[0]=='F' || nm1[0]=='f') endchar='f'; + else CHECKS(0,"end needs to be 'back' or 'front'"); } + filRotateVertex(fil,seg,angle,endchar,symbol); } + else if(!strcmp(nm,"thickness")) { + CHECKS(symbol=='=' || symbol=='+' || symbol=='-',"first symbol needs to be '=', '+', or '-'"); + itct=strmathsscanf(line2,"%mlg|L",varnames,varvalues,nvar,&thick); + CHECKS(itct==1,"modify_segment format: segment length/angle/thickness +/-/= value [end]"); + line2=strnword(line2,2); + er=filChangeThickness(fil,seg,thick,symbol); + CHECKS(!er,"Error while trying to change segment thickness - new value cannot be negative"); } + else + CHECKS(0,"modification option needs to be 'length', 'angle', or 'thickness'"); + CHECKS(!line2,"unexpected text following modify_segment"); } - else if(!strcmp(word,"copy")) { // copy - CHECKS(fil,"need to enter filament name before copy"); + else if(!strcmp(word,"copy_to")) { // copy_to + CHECKS(fil,"need to enter filament name before copy_to"); itct=sscanf(line2,"%s",nm); CHECKS(itct==1,"error reading filament name"); - fil2=filAddFilament(fil->filtype,nm); - CHECKS(fil2,"failed to add filament"); - er=filCopyFilament(fil2,fil); - CHECKS(!strnword(line2,2),"unexpected text following copy"); } + filtype2=filtype; + fil2=NULL; + er=filReadFilName(sim,nm,&filtype2,&fil2,nm1); + CHECKS(er!=-1 && er!=-3,"error reading filament name"); + CHECKS(er!=-4,"unknown filament type"); + CHECKS(fil2!=fil,"a filament cannot be copied to itself"); + if(!fil2) { + fil2=filAddFilament(filtype2,nm1); + CHECKS(fil2,"failed to add filament"); } + er=filCopyFilament(fil2,fil,filtype2); + CHECKS(!strnword(line2,2),"unexpected text following copy_to"); } else { // unknown word - CHECKS(0,"syntax error within filament_type block: statement not recognized"); } + CHECKS(0,"syntax error within filament block: statement not recognized"); } return fil; @@ -1861,30 +2048,40 @@ void filNodes2Angles(filamentptr fil) { if(segment->ypr[0]<-PI) segment->ypr[0]+=2*PI; else if(segment->ypr[0]>PI) segment->ypr[0]-=2*PI; segment->ypr[1]=segment->ypr[2]=0; - Sph_Xyz2Dcm(segment->ypr,segment->dcm); // segment->dcm - Sph_Xyz2Dcm(aypr,segment->adcm); // segment->adcm +// Sph_Ypr2Dcm(segment->ypr,segment->dcm); // segment->dcm + Sph_Ypr2Qtn(segment->ypr,segment->qrel); +// Sph_Ypr2Dcm(aypr,segment->adcm); // segment->adcm + Sph_Ypr2Qtn(aypr,segment->qabs); aphim1=aypr[0]; xm1[0]=x[0]; xm1[1]=x[1]; }} else { // 3D -// printf("\nfilNodes2Anges:");//?? - for(seg=0;segnseg;seg++) { + segment=fil->segments[0]; + x[0]=segment->xyzback[0]-segment->xyzfront[0]; // displacement of segment + x[1]=segment->xyzback[1]-segment->xyzfront[1]; + x[2]=segment->xyzback[2]-segment->xyzfront[2]; + len=sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]); + segment->len=len; + Sph_XZ2Qtni(x,fil->seg0up,segment->qrel); + Sph_Qtn2Qtn(segment->qrel,segment->qabs); + Sph_Qtn2Ypr(segment->qrel,segment->ypr); + for(seg=1;segnseg;seg++) { segment=fil->segments[seg]; x[0]=segment->xyzback[0]-segment->xyzfront[0]; // displacement of segment x[1]=segment->xyzback[1]-segment->xyzfront[1]; x[2]=segment->xyzback[2]-segment->xyzfront[2]; len=sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]); segment->len=len; - if(seg>0) Sph_DcmxCart(fil->segments[seg-1]->adcm,x,x); // rotate x into prior segment reference frame - segment->ypr[0]=atan2(x[1],x[0]); // local ypr angles +// Sph_DcmxCart(fil->segments[seg-1]->adcm,x,x); // rotate x into prior segment reference frame + Sph_QtnRotate(fil->segments[seg-1]->qabs,x,x); + segment->ypr[0]=atan2(x[1],x[0]); // local ypr angles ?? This was the problem. For segment 0, if x[1]=x[0]=0, then this blows up. segment->ypr[1]=-asin(x[2]/len); - Sph_Xyz2Dcm(segment->ypr,segment->dcm); -// if(seg==0) -// printf(" x=(%1.3g,%1.3g,%1.3g), YPR=(%1.3g,%1.3g,%1.3g), Dcm=((%1.3g,%1.3g,%1.3g),(%1.3g,%1.3g,%1.3g),(%1.3g,%1.3g,%1.3g))",x[0],x[1],x[2],segment->ypr[0],segment->ypr[1],segment->ypr[2],segment->dcm[0],segment->dcm[1],segment->dcm[2],segment->dcm[3],segment->dcm[4],segment->dcm[5],segment->dcm[6],segment->dcm[7],segment->dcm[8]);//?? - if(seg>0) Sph_DcmxDcm(segment->dcm,fil->segments[seg-1]->adcm,segment->adcm); - else Sph_Dcm2Dcm(segment->dcm,segment->adcm); }} - +// Sph_Ypr2Dcm(segment->ypr,segment->dcm); + Sph_Ypr2Qtn(segment->ypr,segment->qrel); +// Sph_DcmxDcm(segment->dcm,fil->segments[seg-1]->adcm,segment->adcm); // B_i = A_i . B_{i-1} + Sph_QtnxQtn(fil->segments[seg-1]->qabs,segment->qrel,segment->qabs); // b_i = b_{i-1} . a_i + }} return; } @@ -2096,10 +2293,8 @@ void filAddBendForces(filamentptr fil) { forces[node+1][0]+=forcep1[0]; forces[node+1][1]+=forcep1[1]; }} else { -// printf("\nTorque:");//?? for(node=1;nodenseg;node++) { // compute bending forces filBendTorque(fil,node,bendtorque); // get bending torque in system reference frame -// printf(" (%1.2g,%1.2g,%1.2g)",bendtorque[0],bendtorque[1],bendtorque[2]);//?? segmentm1=fil->segments[node-1]; segment=fil->segments[node]; @@ -2149,9 +2344,10 @@ void filComputeForces(filamentptr fil) { forces=fil->forces; torques=fil->torques; - for(node=0;node<=fil->nseg;node++) { // clear all forces + for(node=0;node<=fil->nseg;node++) // clear all forces forces[node][0]=forces[node][1]=forces[node][2]=0; - torques[node]=0; } + for(node=0;nodenseg;node++) // clear all torques + torques[node]=0; filAddStretchForces(fil); filAddBendForces(fil); // filAddThermalForces(fil); //?? commented out for now because it needs checking @@ -2162,40 +2358,44 @@ void filComputeForces(filamentptr fil) { /* filEulerDynamics */ void filEulerDynamics(simptr sim,filamenttypeptr filtype) { int node,seg,f; - double **nodes,mobility; + double **nodes,dtmu,axis[3]; segmentptr segment; filamentptr fil; for(f=0;fnfil;f++) { fil=filtype->fillist[f]; - mobility=fil->filtype->mobility*sim->dt; + dtmu=fil->filtype->mobility*sim->dt; filComputeForces(fil); nodes=fil->nodes; if(sim->dim==2) { for(node=0;node<=fil->nseg;node++) { - nodes[node][0]+=mobility*fil->forces[node][0]; - nodes[node][1]+=mobility*fil->forces[node][1]; }} + nodes[node][0]+=dtmu*fil->forces[node][0]; + nodes[node][1]+=dtmu*fil->forces[node][1]; }} else { - for(node=0;node<=fil->nseg;node++) { - nodes[node][0]+=mobility*fil->forces[node][0]; - nodes[node][1]+=mobility*fil->forces[node][1]; - nodes[node][2]+=mobility*fil->forces[node][2]; } - for(seg=0;segnseg;seg++) { + axis[0]=nodes[1][0]-nodes[0][0]; // process segment 0 torque + axis[1]=nodes[1][1]-nodes[0][1]; + axis[2]=nodes[1][2]-nodes[0][2]; + Sph_RotateVectorAxisAngle(fil->seg0up,axis,dtmu*fil->torques[0],fil->seg0up); + for(seg=1;segnseg;seg++) { // process other segment torques segment=fil->segments[seg]; - segment->ypr[2]+=mobility*fil->torques[seg]; }} + segment->ypr[2]+=dtmu*fil->torques[seg]; } + for(node=0;node<=fil->nseg;node++) { // process node forces + nodes[node][0]+=dtmu*fil->forces[node][0]; + nodes[node][1]+=dtmu*fil->forces[node][1]; + nodes[node][2]+=dtmu*fil->forces[node][2]; }} filNodes2Angles(fil); } return; } void filRK2Dynamics(simptr sim,filamenttypeptr filtype) { int node,seg,f; - double **nodes,**nodes1,*roll1,mobility; + double **nodes,**nodes1,*roll1,dtmu2,dtmu,axis[3]; segmentptr segment; filamentptr fil; for(f=0;fnfil;f++) { // store initial state in nodes1 and roll1 and then take a half step fil=filtype->fillist[f]; - mobility=fil->filtype->mobility*sim->dt*0.5; + dtmu2=fil->filtype->mobility*sim->dt*0.5; filComputeForces(fil); nodes=fil->nodes; nodes1=fil->nodes1; @@ -2204,42 +2404,50 @@ void filRK2Dynamics(simptr sim,filamenttypeptr filtype) { for(node=0;node<=fil->nseg;node++) { nodes1[node][0]=nodes[node][0]; nodes1[node][1]=nodes[node][1]; - nodes[node][0]+=mobility*fil->forces[node][0]; - nodes[node][1]+=mobility*fil->forces[node][1]; }} + nodes[node][0]+=dtmu2*fil->forces[node][0]; + nodes[node][1]+=dtmu2*fil->forces[node][1]; }} else { + Sph_Cart2Cart(fil->seg0up,fil->seg0up1); + axis[0]=nodes[1][0]-nodes[0][0]; // process segment 0 torque + axis[1]=nodes[1][1]-nodes[0][1]; + axis[2]=nodes[1][2]-nodes[0][2]; + Sph_RotateVectorAxisAngle(fil->seg0up,axis,dtmu2*fil->torques[0],fil->seg0up); + for(seg=0;segnseg;seg++) { + segment=fil->segments[seg]; + roll1[seg]=segment->ypr[2]; + segment->ypr[2]+=dtmu2*fil->torques[seg]; } for(node=0;node<=fil->nseg;node++) { nodes1[node][0]=nodes[node][0]; nodes1[node][1]=nodes[node][1]; nodes1[node][2]=nodes[node][2]; - nodes[node][0]+=mobility*fil->forces[node][0]; - nodes[node][1]+=mobility*fil->forces[node][1]; - nodes[node][2]+=mobility*fil->forces[node][2]; } - for(seg=0;segnseg;seg++) { - segment=fil->segments[seg]; - roll1[seg]=segment->ypr[2]; - segment->ypr[2]+=mobility*fil->torques[seg]; }} + nodes[node][0]+=dtmu2*fil->forces[node][0]; + nodes[node][1]+=dtmu2*fil->forces[node][1]; + nodes[node][2]+=dtmu2*fil->forces[node][2]; }} filNodes2Angles(fil); } for(f=0;fnfil;f++) { // compute force from new state and then take full step from initial state fil=filtype->fillist[f]; - mobility=fil->filtype->mobility*sim->dt; + dtmu=fil->filtype->mobility*sim->dt; filComputeForces(fil); nodes=fil->nodes; nodes1=fil->nodes1; roll1=fil->roll1; - roll1=fil->roll1; if(sim->dim==2) { for(node=0;node<=fil->nseg;node++) { - nodes[node][0]=nodes1[node][0]+mobility*fil->forces[node][0]; - nodes[node][1]=nodes1[node][1]+mobility*fil->forces[node][1]; }} + nodes[node][0]=nodes1[node][0]+dtmu*fil->forces[node][0]; + nodes[node][1]=nodes1[node][1]+dtmu*fil->forces[node][1]; }} else { - for(node=0;node<=fil->nseg;node++) { - nodes[node][0]=nodes1[node][0]+mobility*fil->forces[node][0]; - nodes[node][1]=nodes1[node][1]+mobility*fil->forces[node][1]; - nodes[node][2]=nodes1[node][2]+mobility*fil->forces[node][2]; } + axis[0]=nodes[1][0]-nodes[0][0]; // process segment 0 torque + axis[1]=nodes[1][1]-nodes[0][1]; + axis[2]=nodes[1][2]-nodes[0][2]; + Sph_RotateVectorAxisAngle(fil->seg0up1,axis,dtmu*fil->torques[0],fil->seg0up); for(seg=0;segnseg;seg++) { segment=fil->segments[seg]; - segment->ypr[2]=roll1[seg]+mobility*fil->torques[seg]; }} + segment->ypr[2]=roll1[seg]+dtmu*fil->torques[seg]; } + for(node=0;node<=fil->nseg;node++) { + nodes[node][0]=nodes1[node][0]+dtmu*fil->forces[node][0]; + nodes[node][1]=nodes1[node][1]+dtmu*fil->forces[node][1]; + nodes[node][2]=nodes1[node][2]+dtmu*fil->forces[node][2]; }} filNodes2Angles(fil); } return; } @@ -2247,11 +2455,11 @@ void filRK2Dynamics(simptr sim,filamenttypeptr filtype) { void filRK4Dynamics(simptr sim,filamenttypeptr filtype) { int node,seg,f; - double **nodes,**nodes1,**nodes2,*roll1,*roll2,mobility; + double **nodes,**nodes1,**nodes2,*roll1,*roll2,dtmu,axis[3]; segmentptr segment; filamentptr fil; - mobility=filtype->mobility*sim->dt; + dtmu=filtype->mobility*sim->dt; for(f=0;fnfil;f++) { // copy initial state to nodes1/roll1, then take a half step and add to nodes2/roll2 fil=filtype->fillist[f]; filComputeForces(fil); @@ -2264,26 +2472,32 @@ void filRK4Dynamics(simptr sim,filamenttypeptr filtype) { for(node=0;node<=fil->nseg;node++) { nodes1[node][0]=nodes[node][0]; nodes1[node][1]=nodes[node][1]; - nodes[node][0]+=0.5*mobility*fil->forces[node][0]; - nodes[node][1]+=0.5*mobility*fil->forces[node][1]; - nodes2[node][0]=nodes1[node][0]+(1.0/6)*mobility*fil->forces[node][0]; - nodes2[node][1]=nodes1[node][1]+(1.0/6)*mobility*fil->forces[node][1]; }} + nodes[node][0]+=0.5*dtmu*fil->forces[node][0]; + nodes[node][1]+=0.5*dtmu*fil->forces[node][1]; + nodes2[node][0]=nodes1[node][0]+(1.0/6)*dtmu*fil->forces[node][0]; + nodes2[node][1]=nodes1[node][1]+(1.0/6)*dtmu*fil->forces[node][1]; }} else { + Sph_Cart2Cart(fil->seg0up,fil->seg0up1); + axis[0]=nodes[1][0]-nodes[0][0]; // process segment 0 torque + axis[1]=nodes[1][1]-nodes[0][1]; + axis[2]=nodes[1][2]-nodes[0][2]; + Sph_RotateVectorAxisAngle(fil->seg0up,axis,0.5*dtmu*fil->torques[0],fil->seg0up); + Sph_RotateVectorAxisAngle(fil->seg0up1,axis,(1.0/6)*dtmu*fil->torques[0],fil->seg0up2); + for(seg=0;segnseg;seg++) { + segment=fil->segments[seg]; + roll1[seg]=segment->ypr[2]; + segment->ypr[2]+=0.5*dtmu*fil->torques[seg]; + roll2[seg]=roll1[seg]+(1.0/6)*dtmu*fil->torques[seg]; } for(node=0;node<=fil->nseg;node++) { nodes1[node][0]=nodes[node][0]; nodes1[node][1]=nodes[node][1]; nodes1[node][2]=nodes[node][2]; - nodes[node][0]+=0.5*mobility*fil->forces[node][0]; - nodes[node][1]+=0.5*mobility*fil->forces[node][1]; - nodes[node][2]+=0.5*mobility*fil->forces[node][2]; - nodes2[node][0]=nodes1[node][0]+(1.0/6)*mobility*fil->forces[node][0]; - nodes2[node][1]=nodes1[node][1]+(1.0/6)*mobility*fil->forces[node][1]; - nodes2[node][2]=nodes1[node][2]+(1.0/6)*mobility*fil->forces[node][2]; } - for(seg=0;segnseg;seg++) { - segment=fil->segments[seg]; - roll1[seg]=segment->ypr[2]; - segment->ypr[2]+=0.5*mobility*fil->torques[seg]; - roll2[seg]=roll1[seg]+(1.0/6)*mobility*fil->torques[seg]; }} + nodes[node][0]+=0.5*dtmu*fil->forces[node][0]; + nodes[node][1]+=0.5*dtmu*fil->forces[node][1]; + nodes[node][2]+=0.5*dtmu*fil->forces[node][2]; + nodes2[node][0]=nodes1[node][0]+(1.0/6)*dtmu*fil->forces[node][0]; + nodes2[node][1]=nodes1[node][1]+(1.0/6)*dtmu*fil->forces[node][1]; + nodes2[node][2]=nodes1[node][2]+(1.0/6)*dtmu*fil->forces[node][2]; }} filNodes2Angles(fil); } for(f=0;fnfil;f++) { // with new force, take half step from initial state, and add to nodes2/roll2 @@ -2296,22 +2510,30 @@ void filRK4Dynamics(simptr sim,filamenttypeptr filtype) { roll2=fil->roll2; if(sim->dim==2) { for(node=0;node<=fil->nseg;node++) { - nodes[node][0]=nodes1[node][0]+0.5*mobility*fil->forces[node][0]; - nodes[node][1]=nodes1[node][1]+0.5*mobility*fil->forces[node][1]; - nodes2[node][0]+=(1.0/3)*mobility*fil->forces[node][0]; - nodes2[node][1]+=(1.0/3)*mobility*fil->forces[node][1]; }} + nodes[node][0]=nodes1[node][0]+0.5*dtmu*fil->forces[node][0]; + nodes[node][1]=nodes1[node][1]+0.5*dtmu*fil->forces[node][1]; + nodes2[node][0]+=(1.0/3)*dtmu*fil->forces[node][0]; + nodes2[node][1]+=(1.0/3)*dtmu*fil->forces[node][1]; }} else { - for(node=0;node<=fil->nseg;node++) { - nodes[node][0]=nodes1[node][0]+0.5*mobility*fil->forces[node][0]; - nodes[node][1]=nodes1[node][1]+0.5*mobility*fil->forces[node][1]; - nodes[node][2]=nodes1[node][2]+0.5*mobility*fil->forces[node][2]; - nodes2[node][0]+=(1.0/3)*mobility*fil->forces[node][0]; - nodes2[node][1]+=(1.0/3)*mobility*fil->forces[node][1]; - nodes2[node][2]+=(1.0/3)*mobility*fil->forces[node][2]; } + axis[0]=nodes1[1][0]-nodes1[0][0]; // process segment 0 torque + axis[1]=nodes1[1][1]-nodes1[0][1]; + axis[2]=nodes1[1][2]-nodes1[0][2]; + Sph_RotateVectorAxisAngle(fil->seg0up1,axis,0.5*dtmu*fil->torques[0],fil->seg0up); + axis[0]=nodes2[1][0]-nodes2[0][0]; + axis[1]=nodes2[1][1]-nodes2[0][1]; + axis[2]=nodes2[1][2]-nodes2[0][2]; + Sph_RotateVectorAxisAngle(fil->seg0up2,axis,(1.0/3)*dtmu*fil->torques[0],fil->seg0up2); for(seg=0;segnseg;seg++) { segment=fil->segments[seg]; - segment->ypr[2]=roll1[seg]+0.5*mobility*fil->torques[seg]; - roll2[seg]+=(1.0/3)*mobility*fil->torques[seg]; }} + segment->ypr[2]=roll1[seg]+0.5*dtmu*fil->torques[seg]; + roll2[seg]+=(1.0/3)*dtmu*fil->torques[seg]; } + for(node=0;node<=fil->nseg;node++) { + nodes[node][0]=nodes1[node][0]+0.5*dtmu*fil->forces[node][0]; + nodes[node][1]=nodes1[node][1]+0.5*dtmu*fil->forces[node][1]; + nodes[node][2]=nodes1[node][2]+0.5*dtmu*fil->forces[node][2]; + nodes2[node][0]+=(1.0/3)*dtmu*fil->forces[node][0]; + nodes2[node][1]+=(1.0/3)*dtmu*fil->forces[node][1]; + nodes2[node][2]+=(1.0/3)*dtmu*fil->forces[node][2]; }} filNodes2Angles(fil); } for(f=0;fnfil;f++) { // with new force, take full step from initial state, and add to nodes2/roll2 @@ -2324,22 +2546,30 @@ void filRK4Dynamics(simptr sim,filamenttypeptr filtype) { roll2=fil->roll2; if(sim->dim==2) { for(node=0;node<=fil->nseg;node++) { - nodes[node][0]=nodes1[node][0]+mobility*fil->forces[node][0]; - nodes[node][1]=nodes1[node][1]+mobility*fil->forces[node][1]; - nodes2[node][0]+=(1.0/3)*mobility*fil->forces[node][0]; - nodes2[node][1]+=(1.0/3)*mobility*fil->forces[node][1]; }} + nodes[node][0]=nodes1[node][0]+dtmu*fil->forces[node][0]; + nodes[node][1]=nodes1[node][1]+dtmu*fil->forces[node][1]; + nodes2[node][0]+=(1.0/3)*dtmu*fil->forces[node][0]; + nodes2[node][1]+=(1.0/3)*dtmu*fil->forces[node][1]; }} else { - for(node=0;node<=fil->nseg;node++) { - nodes[node][0]=nodes1[node][0]+mobility*fil->forces[node][0]; - nodes[node][1]=nodes1[node][1]+mobility*fil->forces[node][1]; - nodes[node][2]=nodes1[node][2]+mobility*fil->forces[node][2]; - nodes2[node][0]+=(1.0/3)*mobility*fil->forces[node][0]; - nodes2[node][1]+=(1.0/3)*mobility*fil->forces[node][1]; - nodes2[node][2]+=(1.0/3)*mobility*fil->forces[node][2]; } + axis[0]=nodes1[1][0]-nodes1[0][0]; // process segment 0 torque + axis[1]=nodes1[1][1]-nodes1[0][1]; + axis[2]=nodes1[1][2]-nodes1[0][2]; + Sph_RotateVectorAxisAngle(fil->seg0up1,axis,dtmu*fil->torques[0],fil->seg0up); + axis[0]=nodes1[1][0]-nodes1[0][0]; + axis[1]=nodes1[1][1]-nodes1[0][1]; + axis[2]=nodes1[1][2]-nodes1[0][2]; + Sph_RotateVectorAxisAngle(fil->seg0up2,axis,(1.0/3)*dtmu*fil->torques[0],fil->seg0up2); for(seg=0;segnseg;seg++) { segment=fil->segments[seg]; - segment->ypr[2]=roll1[seg]+mobility*fil->torques[seg]; - roll2[seg]+=(1.0/3)*mobility*fil->torques[seg]; }} + segment->ypr[2]=roll1[seg]+dtmu*fil->torques[seg]; + roll2[seg]+=(1.0/3)*dtmu*fil->torques[seg]; } + for(node=0;node<=fil->nseg;node++) { + nodes[node][0]=nodes1[node][0]+dtmu*fil->forces[node][0]; + nodes[node][1]=nodes1[node][1]+dtmu*fil->forces[node][1]; + nodes[node][2]=nodes1[node][2]+dtmu*fil->forces[node][2]; + nodes2[node][0]+=(1.0/3)*dtmu*fil->forces[node][0]; + nodes2[node][1]+=(1.0/3)*dtmu*fil->forces[node][1]; + nodes2[node][2]+=(1.0/3)*dtmu*fil->forces[node][2]; }} filNodes2Angles(fil); } for(f=0;fnfil;f++) { // with new force, add to nodes2/roll2 @@ -2352,16 +2582,20 @@ void filRK4Dynamics(simptr sim,filamenttypeptr filtype) { roll2=fil->roll2; if(sim->dim==2) { for(node=0;node<=fil->nseg;node++) { - nodes[node][0]=nodes2[node][0]+(1.0/6)*mobility*fil->forces[node][0]; - nodes[node][1]=nodes2[node][1]+(1.0/6)*mobility*fil->forces[node][1]; }} + nodes[node][0]=nodes2[node][0]+(1.0/6)*dtmu*fil->forces[node][0]; + nodes[node][1]=nodes2[node][1]+(1.0/6)*dtmu*fil->forces[node][1]; }} else { - for(node=0;node<=fil->nseg;node++) { - nodes[node][0]=nodes2[node][0]+(1.0/6)*mobility*fil->forces[node][0]; - nodes[node][1]=nodes2[node][1]+(1.0/6)*mobility*fil->forces[node][1]; - nodes[node][2]=nodes2[node][2]+(1.0/6)*mobility*fil->forces[node][2]; } + axis[0]=nodes2[1][0]-nodes2[0][0]; // process segment 0 torque + axis[1]=nodes2[1][1]-nodes2[0][1]; + axis[2]=nodes2[1][2]-nodes2[0][2]; + Sph_RotateVectorAxisAngle(fil->seg0up2,axis,(1.0/6)*dtmu*fil->torques[0],fil->seg0up); for(seg=0;segnseg;seg++) { segment=fil->segments[seg]; - segment->ypr[2]=roll2[seg]+(1.0/6)*mobility*fil->torques[seg]; }} + segment->ypr[2]=roll2[seg]+(1.0/6)*dtmu*fil->torques[seg]; } + for(node=0;node<=fil->nseg;node++) { + nodes[node][0]=nodes2[node][0]+(1.0/6)*dtmu*fil->forces[node][0]; + nodes[node][1]=nodes2[node][1]+(1.0/6)*dtmu*fil->forces[node][1]; + nodes[node][2]=nodes2[node][2]+(1.0/6)*dtmu*fil->forces[node][2]; }} filNodes2Angles(fil); } return; } diff --git a/source/Smoldyn/smolsim.cpp b/source/Smoldyn/smolsim.cpp index dca8792e..01afc6c2 100644 --- a/source/Smoldyn/smolsim.cpp +++ b/source/Smoldyn/smolsim.cpp @@ -635,7 +635,6 @@ int simreadstring(simptr sim,ParseFilePtr pfp,const char *word,char *line2) { portptr port; filamentptr fil; filamenttypeptr filtype; - filamentssptr filss; long int li1; listptrli lilist; listptrv vlist; @@ -1542,6 +1541,19 @@ int simreadstring(simptr sim,ParseFilePtr pfp,const char *word,char *line2) { if(!fil) pfp=NULL; CHECK(fil!=NULL); } + else if(!strcmp(word,"filament_type")) { // filament_type + CHECKS(sim->filss,"individual filament types need to be defined before using filament_type"); + itct=sscanf(line2,"%s %s",nm,nm1); + CHECKS(itct==2,"filament_type format: filament_type_name statement_name statement_text"); + line2=strnword(line2,3); + CHECKS(line2,"filament_type format: filament_type_name statement_name statement_text"); + ft=stringfind(sim->filss->ftnames,sim->filss->ntype,nm); + CHECKS(ft>=0,"filament type is unrecognized"); + filtype=sim->filss->filtypes[ft]; + filtype=filtypereadstring(sim,pfp,filtype,nm1,line2); + if(!filtype) pfp=NULL; + CHECK(filtype!=NULL); } + else if(!strcmp(word,"filament")) { // filament CHECKS(sim->filss,"individual filaments need to be defined before using filament"); itct=sscanf(line2,"%s %s",nm,nm1); @@ -1559,31 +1571,31 @@ int simreadstring(simptr sim,ParseFilePtr pfp,const char *word,char *line2) { else if(!strcmp(word,"random_filament")) { // random_filament CHECKS(sim->filss,"need to enter a filament type before random_filament"); - filss=sim->filss; - itct=sscanf(line2,"%s %s",nm,nm1); - CHECKS(itct==2,"random_filament format: name type segments [x y z theta phi chi] [thickness]"); - ft=stringfind(filss->ftnames,filss->ntype,nm1); - CHECKS(ft>=0,"filament type is unknown"); - filtype=filss->filtypes[ft]; + itct=sscanf(line2,"%s",nm); + CHECKS(itct==1,"random_filament format: type:name segments [x y z phi theta psi] [thickness]"); + er=filReadFilName(sim,nm,&filtype,&fil,nm1); + CHECKS(er!=-1 && er!=-3,"cannot read filament name"); + CHECKS(er!=-4,"unknonwn filament type"); + CHECKS(filtype,"missing filament type. Format: type:name segments [x y z phi theta psi] [thickness]"); CHECKS(filtype->klen==-1 || filtype->klen>0,"cannot compute random segments because the filament type has length force constant equals 0"); - fil=filAddFilament(filtype,nm); + fil=filAddFilament(filtype,nm1); CHECKS(fil,"unable to add filament to simulation"); - line2=strnword(line2,3); + line2=strnword(line2,2); - CHECKS(line2,"random_filament format: name type segments [x y z theta phi chi] [thickness]"); + CHECKS(line2,"random_filament format: type:name segments [x y z phi theta psi] [thickness]"); itct=strmathsscanf(line2,"%mi",varnames,varvalues,nvar,&i1); // number of segments - CHECKS(itct==1,"random_filament format: name type segments [x y z theta phi chi] [thickness]"); - CHECKS(i1>0,"number needs to be >0"); + CHECKS(itct==1,"random_filament format: type:name segments [x y z phi theta psi] [thickness]"); + CHECKS(i1>0,"number of segments needs to be >0"); line2=strnword(line2,2); if(fil->nseg==0 && dim==3) { CHECKS(line2,"missing position and angle information"); itct=sscanf(line2,"%s %s %s %s %s %s",str1,str2,str3,str4,str5,str6); - CHECKS(itct==6,"random_filament format: name type number [x y z theta phi chi] [thickness]"); + CHECKS(itct==6,"random_filament format: type:name number [x y z phi theta psi] [thickness]"); line2=strnword(line2,7); } else if(fil->nseg==0 && dim==2) { CHECKS(line2,"missing position and angle information"); itct=sscanf(line2,"%s %s %s",str1,str2,str4); - CHECKS(itct==3,"random_filament format: name type number [x y theta] [thickness]"); + CHECKS(itct==3,"random_filament format: type:name number [x y phi] [thickness]"); sprintf(str3,"%i",0); sprintf(str5,"%i",0); sprintf(str6,"%i",0); @@ -1598,7 +1610,7 @@ int simreadstring(simptr sim,ParseFilePtr pfp,const char *word,char *line2) { thick=1; if(line2) { itct=strmathsscanf(line2,"%mlg|L",varnames,varvalues,nvar,&thick); - CHECKS(itct==1,"random_segments format: name type number [x y z theta phi chi] [thickness]"); + CHECKS(itct==1,"random_segments format: type:name number [x y z phi theta psi] [thickness]"); CHECKS(thick>0,"thickness needs to be >0"); line2=strnword(line2,2); } er=filAddRandomSegments(fil,i1,str1,str2,str3,str4,str5,str6,thick); diff --git a/source/libSteve/Rn.c b/source/libSteve/Rn.c index 9a4d06b3..a218d36b 100755 --- a/source/libSteve/Rn.c +++ b/source/libSteve/Rn.c @@ -443,16 +443,16 @@ double dotVVD(double *a,double *b,int n) { void crossVV(float *a,float *b,float *c) { - c[0]=a[2]*b[3]-a[3]*b[2]; - c[1]=a[3]*b[1]-a[1]*b[3]; - c[2]=a[1]*b[2]-a[2]*b[1]; + c[0]=a[1]*b[2]-a[2]*b[1]; + c[1]=a[2]*b[0]-a[0]*b[2]; + c[2]=a[0]*b[1]-a[1]*b[0]; return; } void crossVVD(const double *a,const double *b,double *c) { - c[0]=a[2]*b[3]-a[3]*b[2]; - c[1]=a[3]*b[1]-a[1]*b[3]; - c[2]=a[1]*b[2]-a[2]*b[1]; + c[0]=a[1]*b[2]-a[2]*b[1]; + c[1]=a[2]*b[0]-a[0]*b[2]; + c[2]=a[0]*b[1]-a[1]*b[0]; return; } diff --git a/source/libSteve/SimCommand.c b/source/libSteve/SimCommand.c index 0f8937ce..e76e90d3 100644 --- a/source/libSteve/SimCommand.c +++ b/source/libSteve/SimCommand.c @@ -556,7 +556,7 @@ enum CMDcode scmdcmdtype(cmdssptr cmds,cmdptr cmd) { /* scmdoutput */ void scmdoutput(cmdssptr cmds) { int fid,i,did; - queue cmdq; + q_queue cmdq; cmdptr cmd; void *voidptr,*simvd; char timing,string[STRCHAR],string2[STRCHAR]; diff --git a/source/libSteve/SimCommand.h b/source/libSteve/SimCommand.h index f16f1ff4..d5679909 100644 --- a/source/libSteve/SimCommand.h +++ b/source/libSteve/SimCommand.h @@ -46,8 +46,8 @@ typedef struct cmdsuperstruct { int maxcmdlist; // allocated size of command list int ncmdlist; // actual size of command list cmdptr *cmdlist; // list of all added commands - queue cmd; // queue of normal run-time commands - queue cmdi; // queue of integer time commands + q_queue cmd; // queue of normal run-time commands + q_queue cmdi; // queue of integer time commands enum CMDcode (*cmdfn)(void*,cmdptr,char*); // function that runs commands void *simvd; // void pointer to simulation structure int iter; // number of times integer commands have run diff --git a/source/libSteve/Sphere.c b/source/libSteve/Sphere.c index c3a10e2a..ddccc578 100644 --- a/source/libSteve/Sphere.c +++ b/source/libSteve/Sphere.c @@ -13,10 +13,21 @@ of the Gnu Lesser General Public License (LGPL). */ #define PI 3.14159265358979323846 #define EPSILON 100*DBL_EPSILON -double Work[9],Work2[9]; +/* *** Copied over from Rn.c to avoid a dependency *** */ + +void Sph_crossVVD(const double *a,const double *b,double *c) { + c[0]=a[1]*b[2]-a[2]*b[1]; + c[1]=a[2]*b[0]-a[0]*b[2]; + c[2]=a[0]*b[1]-a[1]*b[0]; + return; } + + +/* *** Start of Sphere.c library functions *** */ void Sph_Cart2Sc(const double *Cart,double *Sc) { + double Work[3]; + Work[0]=sqrt(Cart[0]*Cart[0]+Cart[1]*Cart[1]+Cart[2]*Cart[2]); if(Work[0]>0) { Work[1]=acos(Cart[2]/Work[0]); @@ -30,6 +41,8 @@ void Sph_Cart2Sc(const double *Cart,double *Sc) { void Sph_Sc2Cart(const double *Sc,double *Cart) { + double Work[3]; + Work[0]=Sc[0]*sin(Sc[1])*cos(Sc[2]); Work[1]=Sc[0]*sin(Sc[1])*sin(Sc[2]); Work[2]=Sc[0]*cos(Sc[1]); @@ -39,26 +52,21 @@ void Sph_Sc2Cart(const double *Sc,double *Cart) { return; } -void Sph_Eay2Ep(const double *Eay,double *Ep) { - Work[0]=cos(0.5*(Eay[3]+Eay[2]))*cos(0.5*Eay[1]); - Work[1]=sin(0.5*(Eay[3]-Eay[2]))*sin(0.5*Eay[1]); - Work[2]=cos(0.5*(Eay[3]-Eay[2]))*sin(0.5*Eay[1]); - Work[3]=sin(0.5*(Eay[3]+Eay[2]))*cos(0.5*Eay[1]); - Ep[0]=Work[0]; - Ep[1]=Work[1]; - Ep[2]=Work[2]; - Ep[3]=Work[3]; +void Sph_Cart2Cart(const double *Cart1,double *Cart2) { + Cart2[0]=Cart1[0]; + Cart2[1]=Cart1[1]; + Cart2[2]=Cart1[2]; return; } -void Sph_Xyz2Xyz(const double *Xyz1,double *Xyz2) { - Xyz2[0]=Xyz1[0]; - Xyz2[1]=Xyz1[1]; - Xyz2[2]=Xyz1[2]; +void Sph_Ypr2Ypr(const double *Ypr1,double *Ypr2) { + Ypr2[0]=Ypr1[0]; + Ypr2[1]=Ypr1[1]; + Ypr2[2]=Ypr1[2]; return; } -void Sph_Eax2Xyz(const double *Eax,double *Xyz) { +void Sph_Eax2Ypr(const double *Eax,double *Ypr) { double cf,cq,cy,sf,sq,sy; cq=cos(Eax[0]); @@ -67,9 +75,9 @@ void Sph_Eax2Xyz(const double *Eax,double *Xyz) { sq=sin(Eax[0]); sf=sin(Eax[1]); sy=sin(Eax[2]); - Xyz[0]=atan2(cy*sf+cq*cf*sy,cy*cf-cq*sf*sy); - Xyz[1]=asin(-sy*sq); - Xyz[2]=atan2(cy*sq,cq); + Ypr[0]=atan2(cy*sf+cq*cf*sy,cy*cf-cq*sf*sy); + Ypr[1]=asin(-sy*sq); + Ypr[2]=atan2(cy*sq,cq); return; } @@ -115,15 +123,15 @@ void Sph_Eay2Dcm(const double *Eay,double *Dcm) { return; } -void Sph_Xyz2Dcm(const double *Xyz,double *Dcm) { +void Sph_Ypr2Dcm(const double *Ypr,double *Dcm) { double cf,cq,cy,sf,sq,sy; - cf=cos(Xyz[0]); - cq=cos(Xyz[1]); - cy=cos(Xyz[2]); - sf=sin(Xyz[0]); - sq=sin(Xyz[1]); - sy=sin(Xyz[2]); + cf=cos(Ypr[0]); + cq=cos(Ypr[1]); + cy=cos(Ypr[2]); + sf=sin(Ypr[0]); + sq=sin(Ypr[1]); + sy=sin(Ypr[2]); Dcm[0]=cq*cf; Dcm[1]=cq*sf; Dcm[2]=-sq; @@ -136,15 +144,15 @@ void Sph_Xyz2Dcm(const double *Xyz,double *Dcm) { return; } -void Sph_Xyz2Dcmt(const double *Xyz,double *Dcmt) { +void Sph_Ypr2Dcmt(const double *Ypr,double *Dcmt) { double cf,cq,cy,sf,sq,sy; - cf=cos(Xyz[0]); - cq=cos(Xyz[1]); - cy=cos(Xyz[2]); - sf=sin(Xyz[0]); - sq=sin(Xyz[1]); - sy=sin(Xyz[2]); + cf=cos(Ypr[0]); + cq=cos(Ypr[1]); + cy=cos(Ypr[2]); + sf=sin(Ypr[0]); + sq=sin(Ypr[1]); + sy=sin(Ypr[2]); Dcmt[0]=cq*cf; Dcmt[3]=cq*sf; Dcmt[6]=-sq; @@ -157,13 +165,15 @@ void Sph_Xyz2Dcmt(const double *Xyz,double *Dcmt) { return; } -void Sph_Dcm2Xyz(const double *Dcm,double *Xyz) { +void Sph_Dcm2Ypr(const double *Dcm,double *Ypr) { + double Work[3]; + Work[0]=atan2(Dcm[1],Dcm[0]); Work[1]=asin(-Dcm[2]); Work[2]=atan2(Dcm[5],Dcm[8]); - Xyz[0]=Work[0]; - Xyz[1]=Work[1]; - Xyz[2]=Work[2]; + Ypr[0]=Work[0]; + Ypr[1]=Work[1]; + Ypr[2]=Work[2]; return; } @@ -181,6 +191,8 @@ void Sph_Dcm2Dcm(const double *Dcm1,double *Dcm2) { void Sph_Dcm2Dcmt(const double *Dcm1,double *Dcm2) { + double Work[9]; + Dcm2[0]=Dcm1[0]; Dcm2[4]=Dcm1[4]; Dcm2[8]=Dcm1[8]; @@ -200,6 +212,8 @@ void Sph_Dcm2Dcmt(const double *Dcm1,double *Dcm2) { void Sph_DcmxDcm(const double *Dcm1,const double *Dcm2,double *Dcm3) { + double Work[9]; + Work[0]=Dcm1[0]*Dcm2[0]+Dcm1[1]*Dcm2[3]+Dcm1[2]*Dcm2[6]; Work[1]=Dcm1[0]*Dcm2[1]+Dcm1[1]*Dcm2[4]+Dcm1[2]*Dcm2[7]; Work[2]=Dcm1[0]*Dcm2[2]+Dcm1[1]*Dcm2[5]+Dcm1[2]*Dcm2[8]; @@ -222,6 +236,8 @@ void Sph_DcmxDcm(const double *Dcm1,const double *Dcm2,double *Dcm3) { void Sph_DcmxDcmt(const double *Dcm1,const double *Dcmt,double *Dcm3) { + double Work[9]; + Work[0]=Dcm1[0]*Dcmt[0]+Dcm1[1]*Dcmt[1]+Dcm1[2]*Dcmt[2]; Work[1]=Dcm1[0]*Dcmt[3]+Dcm1[1]*Dcmt[4]+Dcm1[2]*Dcmt[5]; Work[2]=Dcm1[0]*Dcmt[6]+Dcm1[1]*Dcmt[7]+Dcm1[2]*Dcmt[8]; @@ -244,6 +260,8 @@ void Sph_DcmxDcmt(const double *Dcm1,const double *Dcmt,double *Dcm3) { void Sph_DcmtxDcm(const double *Dcmt,const double *Dcm2,double *Dcm3) { + double Work[9]; + Work[0]=Dcmt[0]*Dcm2[0]+Dcmt[3]*Dcm2[3]+Dcmt[6]*Dcm2[6]; Work[1]=Dcmt[0]*Dcm2[1]+Dcmt[3]*Dcm2[4]+Dcmt[6]*Dcm2[7]; Work[2]=Dcmt[0]*Dcm2[2]+Dcmt[3]*Dcm2[5]+Dcmt[6]*Dcm2[8]; @@ -266,6 +284,8 @@ void Sph_DcmtxDcm(const double *Dcmt,const double *Dcm2,double *Dcm3) { void Sph_DcmxCart(const double *Dcm,const double *Cart,double *Cart2) { + double Work[3]; + Work[0]=Dcm[0]*Cart[0]+Dcm[1]*Cart[1]+Dcm[2]*Cart[2]; Work[1]=Dcm[3]*Cart[0]+Dcm[4]*Cart[1]+Dcm[5]*Cart[2]; Work[2]=Dcm[6]*Cart[0]+Dcm[7]*Cart[1]+Dcm[8]*Cart[2]; @@ -276,6 +296,8 @@ void Sph_DcmxCart(const double *Dcm,const double *Cart,double *Cart2) { void Sph_DcmtxCart(const double *Dcm,const double *Cart,double *Cart2) { + double Work[3]; + Work[0]=Dcm[0]*Cart[0]+Dcm[3]*Cart[1]+Dcm[6]*Cart[2]; Work[1]=Dcm[1]*Cart[0]+Dcm[4]*Cart[1]+Dcm[7]*Cart[2]; Work[2]=Dcm[2]*Cart[0]+Dcm[5]*Cart[1]+Dcm[8]*Cart[2]; @@ -291,18 +313,18 @@ void Sph_One2Dcm(double *Dcm) { return; } -void Sph_Xyz2Xyzr(const double *Xyz,double *Xyzr) { +void Sph_Ypr2Yprr(const double *Ypr,double *Yprr) { double cf,cq,cy,sf,sq,sy; - cf=cos(Xyz[0]); - cq=cos(Xyz[1]); - cy=cos(Xyz[2]); - sf=sin(Xyz[0]); - sq=sin(Xyz[1]); - sy=sin(Xyz[2]); - Xyzr[0]=atan2(sy*sq*cf-cy*sf,cq*cf); - Xyzr[1]=-asin(-cy*sq*cf-sy*sf); - Xyzr[2]=atan2(-cy*sq*sf+sy*cf,cq*cy); + cf=cos(Ypr[0]); + cq=cos(Ypr[1]); + cy=cos(Ypr[2]); + sf=sin(Ypr[0]); + sq=sin(Ypr[1]); + sy=sin(Ypr[2]); + Yprr[0]=atan2(sy*sq*cf-cy*sf,cq*cf); + Yprr[1]=-asin(-cy*sq*cf-sy*sf); + Yprr[2]=atan2(-cy*sq*sf+sy*cf,cq*cy); return; } @@ -339,15 +361,19 @@ void Sph_Rot2Dcm(char axis,double angle,double *Dcm) { void Sph_Newz2Dcm(const double *Newz,double psi,double *Dcm) { - Sph_Cart2Sc(Newz,Work2); - Work2[2]+=PI/2.0; - Work2[3]=psi-Work2[2]; - Sph_Eax2Dcm(Work2+1,Dcm); + double Work[4]; + + Sph_Cart2Sc(Newz,Work); + Work[2]+=PI/2.0; + Work[3]=psi-Work[2]; + Sph_Eax2Dcm(Work+1,Dcm); Sph_Dcm2Dcmt(Dcm,Dcm); return; } void Sph_DcmtxUnit(const double *Dcmt,char axis,double *vect,const double *add,double mult) { + double Work[3]; + if(add) { Work[0]=add[0]; Work[1]=add[1]; @@ -371,6 +397,219 @@ void Sph_DcmtxUnit(const double *Dcmt,char axis,double *vect,const double *add,d return; } +/* *** Quaternions *** */ + +void Sph_One2Qtn(double *Qtn) { + Qtn[0]=1; + Qtn[1]=Qtn[2]=Qtn[3]=0; + return; } + + +void Sph_Qtn2Qtn(const double *Qtn1,double *Qtn2) { + Qtn2[0]=Qtn1[0]; + Qtn2[1]=Qtn1[1]; + Qtn2[2]=Qtn1[2]; + Qtn2[3]=Qtn1[3]; + return; } + + +void Sph_Ypr2Qtn(const double *Ypr,double *Qtn) { + double cf2,cq2,cy2,sf2,sq2,sy2; + + cf2=cos(Ypr[0]/2); + sf2=sin(Ypr[0]/2); + cq2=cos(Ypr[1]/2); + sq2=sin(Ypr[1]/2); + cy2=cos(Ypr[2]/2); + sy2=sin(Ypr[2]/2); + Qtn[0]=cy2*cq2*cf2+sy2*sq2*sf2; + Qtn[1]=-sy2*cq2*cf2+cy2*sq2*sf2; + Qtn[2]=-cy2*sq2*cf2-sy2*cq2*sf2; + Qtn[3]=-cy2*cq2*sf2+sy2*sq2*cf2; + return; } + + +void Sph_Ypr2Qtni(const double *Ypr,double *Qtni) { + double cf2,cq2,cy2,sf2,sq2,sy2; + + cf2=cos(Ypr[0]/2); + sf2=sin(Ypr[0]/2); + cq2=cos(Ypr[1]/2); + sq2=sin(Ypr[1]/2); + cy2=cos(Ypr[2]/2); + sy2=sin(Ypr[2]/2); + Qtni[0]=cy2*cq2*cf2+sy2*sq2*sf2; + Qtni[1]=sy2*cq2*cf2-cy2*sq2*sf2; + Qtni[2]=cy2*sq2*cf2+sy2*cq2*sf2; + Qtni[3]=cy2*cq2*sf2-sy2*sq2*cf2; + return; } + + +void Sph_Qtn2Ypr(const double *Qtn,double *Ypr) { + Ypr[0]=atan2(-2*Qtn[0]*Qtn[3]+2*Qtn[1]*Qtn[2],Qtn[0]*Qtn[0]+Qtn[1]*Qtn[1]-Qtn[2]*Qtn[2]-Qtn[3]*Qtn[3]); + Ypr[1]=asin(-2*Qtn[0]*Qtn[2]-2*Qtn[1]*Qtn[3]); + Ypr[2]=atan2(-2*Qtn[0]*Qtn[1]+2*Qtn[2]*Qtn[3],Qtn[0]*Qtn[0]-Qtn[1]*Qtn[1]-Qtn[2]*Qtn[2]+Qtn[3]*Qtn[3]); + return; } + + +void Sph_Dcm2Qtn(const double *Dcm,double *Qtn) { + double factor; + + Qtn[0]=Dcm[0]+Dcm[4]+Dcm[8]; + Qtn[1]=Dcm[0]-Dcm[4]-Dcm[8]; + Qtn[2]=-Dcm[0]+Dcm[4]-Dcm[8]; + Qtn[3]=-Dcm[0]-Dcm[4]+Dcm[8]; + if(Qtn[0]>=Qtn[1] && Qtn[0]>=Qtn[2] && Qtn[0]>=Qtn[3]) { +// printf("*** Sph_Dcm2Qtn A. Qtn=(%1.3g %1.3g %1.3g %1.3g)\n",Qtn[0],Qtn[1],Qtn[2],Qtn[3]);//?? + Qtn[0]=0.5*sqrt(1+Qtn[0]); + factor=0.25/Qtn[0]; + Qtn[1]=factor*(Dcm[7]-Dcm[5]); + Qtn[2]=factor*(Dcm[2]-Dcm[6]); + Qtn[3]=factor*(Dcm[3]-Dcm[1]); } + else if(Qtn[1]>=Qtn[2] && Qtn[1]>=Qtn[3]) { +// printf("*** Sph_Dcm2Qtn B. Qtn=(%1.3g %1.3g %1.3g %1.3g)\n",Qtn[0],Qtn[1],Qtn[2],Qtn[3]);//?? + Qtn[1]=0.5*sqrt(1+Qtn[1]); + factor=0.25/Qtn[1]; + Qtn[0]=factor*(Dcm[7]-Dcm[5]); + Qtn[2]=factor*(Dcm[1]+Dcm[3]); + Qtn[3]=factor*(Dcm[2]+Dcm[6]); } + else if(Qtn[2]>=Qtn[3]) { +// printf("*** Sph_Dcm2Qtn C. Qtn=(%1.3g %1.3g %1.3g %1.3g)\n",Qtn[0],Qtn[1],Qtn[2],Qtn[3]);//?? + Qtn[2]=0.5*sqrt(1+Qtn[2]); + factor=0.25/Qtn[2]; + Qtn[0]=factor*(Dcm[2]-Dcm[6]); + Qtn[1]=factor*(Dcm[1]+Dcm[3]); + Qtn[3]=factor*(Dcm[5]+Dcm[7]); } + else { +// printf("*** Sph_Dcm2Qtn D. Qtn=(%1.3g %1.3g %1.3g %1.3g)\n",Qtn[0],Qtn[1],Qtn[2],Qtn[3]);//?? + Qtn[3]=0.5*sqrt(1+Qtn[3]); + factor=0.25/Qtn[3]; + Qtn[0]=factor*(Dcm[3]-Dcm[1]); + Qtn[1]=factor*(Dcm[2]+Dcm[6]); + Qtn[2]=factor*(Dcm[5]+Dcm[7]); } +// printf("*** Sph_Dcm2Qtn Z. Qtn=(%1.3g %1.3g %1.3g %1.3g)\n",Qtn[0],Qtn[1],Qtn[2],Qtn[3]);//?? + return; } + + +void Sph_Qtn2Dcm(const double *Qtn,double *Dcm) { + Dcm[0]=1-2*Qtn[2]*Qtn[2]-2*Qtn[3]*Qtn[3]; + Dcm[1]=2*Qtn[1]*Qtn[2]-2*Qtn[0]*Qtn[3]; + Dcm[2]=2*Qtn[1]*Qtn[3]+2*Qtn[0]*Qtn[2]; + Dcm[3]=2*Qtn[1]*Qtn[2]+2*Qtn[0]*Qtn[3]; + Dcm[4]=1-2*Qtn[1]*Qtn[1]-2*Qtn[3]*Qtn[3]; + Dcm[5]=2*Qtn[2]*Qtn[3]-2*Qtn[0]*Qtn[1]; + Dcm[6]=2*Qtn[1]*Qtn[3]-2*Qtn[0]*Qtn[2]; + Dcm[7]=2*Qtn[2]*Qtn[3]+2*Qtn[0]*Qtn[1]; + Dcm[8]=1-2*Qtn[1]*Qtn[1]-2*Qtn[2]*Qtn[2]; + return; } + + +void Sph_XZ2Qtni(const double *x,const double *z,double *Qtn) { + double *x2,*y2,*z2,dcm[9],len; + + x2=dcm; + y2=dcm+3; + z2=dcm+6; + len=sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]); + x2[0]=x[0]/len; // normalize x and store in first row + x2[1]=x[1]/len; + x2[2]=x[2]/len; + Sph_crossVVD(z,x2,y2); // y2 = z cross x2 + len=sqrt(y2[0]*y2[0]+y2[1]*y2[1]+y2[2]*y2[2]); + y2[0]/=len; // normalize y2 + y2[1]/=len; + y2[2]/=len; + Sph_crossVVD(x2,y2,z2); // z2 = x2 cross y2 +// printf("*** Sph_XZ2Qtni A. dcm=(%1.3g %1.3g %1.3g; %1.3g %1.3g %1.3g; %1.3g %1.3g %1.3g)\n",dcm[0],dcm[1],dcm[2],dcm[3],dcm[4],dcm[5],dcm[6],dcm[7],dcm[8]);//?? + Sph_Dcm2Qtn(dcm,Qtn); // convert matrix to qtn +// printf("*** Sph_XZ2Qtni B. qtn=(%1.3g %1.3g %1.3g %1.3g)\n",Qtn[0],Qtn[1],Qtn[2],Qtn[3]);//?? + return; } + + +void Sph_QtnxQtn(const double *Qtn1,const double *Qtn2,double *Qtn3) { + double Work[4]; + + Work[0]=Qtn1[0]*Qtn2[0]-Qtn1[1]*Qtn2[1]-Qtn1[2]*Qtn2[2]-Qtn1[3]*Qtn2[3]; + Work[1]=Qtn1[0]*Qtn2[1]+Qtn1[1]*Qtn2[0]-Qtn1[2]*Qtn2[3]+Qtn1[3]*Qtn2[2]; + Work[2]=Qtn1[0]*Qtn2[2]+Qtn1[1]*Qtn2[3]+Qtn1[2]*Qtn2[0]-Qtn1[3]*Qtn2[1]; + Work[3]=Qtn1[0]*Qtn2[3]-Qtn1[1]*Qtn2[2]+Qtn1[2]*Qtn2[1]+Qtn1[3]*Qtn2[0]; + Qtn3[0]=Work[0]; + Qtn3[1]=Work[1]; + Qtn3[2]=Work[2]; + Qtn3[3]=Work[3]; + return; } + + +void Sph_QtnixQtn(const double *Qtn1,const double *Qtn2,double *Qtn3) { + double Work[4]; + + Work[0]=Qtn1[0]*Qtn2[0]+Qtn1[1]*Qtn2[1]+Qtn1[2]*Qtn2[2]+Qtn1[3]*Qtn2[3]; + Work[1]=Qtn1[0]*Qtn2[1]-Qtn1[1]*Qtn2[0]+Qtn1[2]*Qtn2[3]-Qtn1[3]*Qtn2[2]; + Work[2]=Qtn1[0]*Qtn2[2]-Qtn1[1]*Qtn2[3]-Qtn1[2]*Qtn2[0]+Qtn1[3]*Qtn2[1]; + Work[3]=Qtn1[0]*Qtn2[3]+Qtn1[1]*Qtn2[2]-Qtn1[2]*Qtn2[1]-Qtn1[3]*Qtn2[0]; + Qtn3[0]=Work[0]; + Qtn3[1]=Work[1]; + Qtn3[2]=Work[2]; + Qtn3[3]=Work[3]; + return; } + + +void Sph_QtnxQtni(const double *Qtn1,const double *Qtn2,double *Qtn3) { + double Work[4]; + + Work[0]=Qtn1[0]*Qtn2[0]+Qtn1[1]*Qtn2[1]+Qtn1[2]*Qtn2[2]+Qtn1[3]*Qtn2[3]; + Work[1]=-Qtn1[0]*Qtn2[1]+Qtn1[1]*Qtn2[0]+Qtn1[2]*Qtn2[3]-Qtn1[3]*Qtn2[2]; + Work[2]=-Qtn1[0]*Qtn2[2]-Qtn1[1]*Qtn2[3]+Qtn1[2]*Qtn2[0]+Qtn1[3]*Qtn2[1]; + Work[3]=-Qtn1[0]*Qtn2[3]+Qtn1[1]*Qtn2[2]-Qtn1[2]*Qtn2[1]+Qtn1[3]*Qtn2[0]; + Qtn3[0]=Work[0]; + Qtn3[1]=Work[1]; + Qtn3[2]=Work[2]; + Qtn3[3]=Work[3]; + return; } + + +void Sph_QtnRotate(const double *Qtn,const double *Cart,double *Cart2) { + double qxc[4]; + + qxc[0]=Cart[0]*Qtn[1]+Cart[1]*Qtn[2]+Cart[2]*Qtn[3]; + qxc[1]=Cart[0]*Qtn[0]+Cart[2]*Qtn[2]-Cart[1]*Qtn[3]; + qxc[2]=Cart[1]*Qtn[0]-Cart[2]*Qtn[1]+Cart[0]*Qtn[3]; + qxc[3]=Cart[2]*Qtn[0]+Cart[1]*Qtn[1]-Cart[0]*Qtn[2]; + Cart2[0]=Qtn[1]*qxc[0]+Qtn[0]*qxc[1]-Qtn[3]*qxc[2]+Qtn[2]*qxc[3]; + Cart2[1]=Qtn[2]*qxc[0]+Qtn[3]*qxc[1]+Qtn[0]*qxc[2]-Qtn[1]*qxc[3]; + Cart2[2]=Qtn[3]*qxc[0]-Qtn[2]*qxc[1]+Qtn[1]*qxc[2]+Qtn[0]*qxc[3]; + return; } + + +void Sph_QtniRotate(const double *Qtn,const double *Cart,double *Cart2) { + double qxc[4]; + + qxc[0]=-Cart[0]*Qtn[1]-Cart[1]*Qtn[2]-Cart[2]*Qtn[3]; + qxc[1]= Cart[0]*Qtn[0]-Cart[2]*Qtn[2]+Cart[1]*Qtn[3]; + qxc[2]= Cart[1]*Qtn[0]+Cart[2]*Qtn[1]-Cart[0]*Qtn[3]; + qxc[3]= Cart[2]*Qtn[0]-Cart[1]*Qtn[1]+Cart[0]*Qtn[2]; + Cart2[0]=-Qtn[1]*qxc[0]+Qtn[0]*qxc[1]+Qtn[3]*qxc[2]-Qtn[2]*qxc[3]; + Cart2[1]=-Qtn[2]*qxc[0]-Qtn[3]*qxc[1]+Qtn[0]*qxc[2]+Qtn[1]*qxc[3]; + Cart2[2]=-Qtn[3]*qxc[0]+Qtn[2]*qxc[1]-Qtn[1]*qxc[2]+Qtn[0]*qxc[3]; + return; } + + +void Sph_QtniRotateUnitx(const double *Qtni,double *vect,const double *add,double mult) { + vect[0]=add[0]+mult*(Qtni[0]*Qtni[0]+Qtni[1]*Qtni[1]-Qtni[2]*Qtni[2]-Qtni[3]*Qtni[3]); + vect[1]=add[1]+mult*(2*Qtni[1]*Qtni[2]-2*Qtni[0]*Qtni[3]); + vect[2]=add[2]+mult*(2*Qtni[0]*Qtni[2]+2*Qtni[1]*Qtni[3]); + return; } + + +void Sph_QtniRotateUnitz(const double *Qtni,double *vect,const double *add,double mult) { + vect[0]=add[0]+mult*(-2*Qtni[0]*Qtni[2]+2*Qtni[1]*Qtni[3]); + vect[1]=add[1]+mult*(2*Qtni[0]*Qtni[1]+2*Qtni[2]*Qtni[3]); + vect[2]=add[2]+mult*(Qtni[0]*Qtni[0]-Qtni[1]*Qtni[1]-Qtni[2]*Qtni[2]+Qtni[3]*Qtni[3]); + return; } + + +/* *** Unit normals *** */ + double Sph_RotateVectWithNormals3D(const double *pt1,const double *pt2,double *newpt2,double *oldnorm,double *newnorm,int sign) { double costheta,ax,ay,az,sintheta,oldnormint[3],*oldnormptr; double costheta1,dcmxx,dcmxy,dcmxz,dcmyx,dcmyy,dcmyz,dcmzx,dcmzy,dcmzz,x,y,z; @@ -436,3 +675,24 @@ double Sph_RotateVectWithNormals3D(const double *pt1,const double *pt2,double *n return costheta; } +void Sph_RotateVectorAxisAngle(const double *vect,const double *axis,double angle,double *rotated) { + double ax[3],len,ca,sa,ans[3]; + + len=sqrt(axis[0]*axis[0]+axis[1]*axis[1]+axis[2]*axis[2]); + ax[0]=axis[0]/len; + ax[1]=axis[1]/len; + ax[2]=axis[2]/len; + + ca=cos(angle); + sa=sin(angle); + + ans[0]=vect[0]*(ca+ax[0]*ax[0]*(1-ca)) + vect[1]*(ax[0]*ax[1]*(1-ca)-ax[2]*sa) + vect[2]*(ax[0]*ax[2]*(1-ca)+ax[1]*sa); + ans[1]=vect[0]*(ax[1]*ax[0]*(1-ca)+ax[2]*sa) + vect[1]*(ca+ax[1]*ax[1]*(1-ca)) + vect[2]*(ax[1]*ax[2]*(1-ca)-ax[0]*sa); + ans[2]=vect[0]*(ax[2]*ax[0]*(1-ca)-ax[1]*sa) + vect[1]*(ax[2]*ax[1]*(1-ca)+ax[0]*sa) + vect[2]*(ca+ax[2]*ax[2]*(1-ca)); + rotated[0]=ans[0]; + rotated[1]=ans[1]; + rotated[2]=ans[2]; + return; } + + + diff --git a/source/libSteve/Sphere.h b/source/libSteve/Sphere.h index 6668519e..f1d4b2b7 100644 --- a/source/libSteve/Sphere.h +++ b/source/libSteve/Sphere.h @@ -13,21 +13,21 @@ Sc = Spherical coordinates (r, theta, phi) Dcm = Direction cosine matrix (A00, A01, A02, A10, A11, A12, A20, A21, A22) Eax = Euler angles x-convention (theta, phi, psi) Eay = Euler angles y-convention (theta, phi, chi) -Ep = Euler parameters (e0, e1, e2, e3) -Xyz = xyz- or ypr- or Tait-Bryan convention (yaw, pitch, roll): rotate on z, then y, then x +Qtn = Quaternions (q0, q1, q2, q3) +Ypr = xyz- or ypr- or Tait-Bryan convention (yaw, pitch, roll): rotate on z, then y, then x */ void Sph_Cart2Sc(const double *Cart,double *Sc); void Sph_Sc2Cart(const double *Sc,double *Cart); -void Sph_Eay2Ep(const double *Eay,double *Ep); -void Sph_Xyz2Xyz(const double *Xyz1,double *Xyz2); -void Sph_Eax2Xyz(const double *Eax,double *Xyz); +void Sph_Cart2Cart(const double *Cart1,double *Cart2); +void Sph_Ypr2Ypr(const double *Ypr1,double *Ypr2); +void Sph_Eax2Ypr(const double *Eax,double *Ypr); void Sph_Eax2Dcm(const double *Eax,double *Dcm); void Sph_Eay2Dcm(const double *Eay,double *Dcm); -void Sph_Xyz2Dcm(const double *Xyz,double *Dcm); -void Sph_Xyz2Dcmt(const double *Xyz,double *Dcmt); -void Sph_Dcm2Xyz(const double *Dcm,double *Xyz); +void Sph_Ypr2Dcm(const double *Ypr,double *Dcm); +void Sph_Ypr2Dcmt(const double *Ypr,double *Dcmt); +void Sph_Dcm2Ypr(const double *Dcm,double *Ypr); void Sph_Dcm2Dcm(const double *Dcm1,double *Dcm2); void Sph_Dcm2Dcmt(const double *Dcm1,double *Dcm2); @@ -38,13 +38,30 @@ void Sph_DcmxCart(const double *Dcm,const double *Cart,double *Cart2); void Sph_DcmtxCart(const double *Dcm,const double *Cart,double *Cart2); void Sph_One2Dcm(double *Dcm); -void Sph_Xyz2Xyzr(const double *Xyz,double *Xyzr); +void Sph_Ypr2Yprr(const double *Ypr,double *Yprr); void Sph_Dcm2Dcmr(const double *Dcm,double *Dcmr); void Sph_Rot2Dcm(char axis,double angle,double *Dcm); void Sph_Newz2Dcm(const double *Newz,double psi,double *Dcm); void Sph_DcmtxUnit(const double *Dcmt,char unit,double *vect,const double *add,double mult); +void Sph_One2Qtn(double *Qtn); +void Sph_Qtn2Qtn(const double *Qtn1,double *Qtn2); +void Sph_Ypr2Qtn(const double *Ypr,double *Qtn); +void Sph_Ypr2Qtni(const double *Ypr,double *Qtni); +void Sph_Qtn2Ypr(const double *Qtn,double *Ypr); +void Sph_Dcm2Qtn(const double *Dcm,double *Qtn); +void Sph_Qtn2Dcm(const double *Qtn,double *Dcm); +void Sph_XZ2Qtni(const double *x,const double *z,double *Qtn); +void Sph_QtnxQtn(const double *Qtn1,const double *Qtn2,double *Qtn3); +void Sph_QtnixQtn(const double *Qtn1,const double *Qtn2,double *Qtn3); +void Sph_QtnxQtni(const double *Qtn1,const double *Qtn2,double *Qtn3); +void Sph_QtnRotate(const double *Qtn,const double *Cart,double *Cart2); +void Sph_QtniRotate(const double *Qtn,const double *Cart,double *Cart2); +void Sph_QtniRotateUnitx(const double *Qtni,double *vect,const double *add,double mult); +void Sph_QtniRotateUnitz(const double *Qtni,double *vect,const double *add,double mult); + double Sph_RotateVectWithNormals3D(const double *pt1,const double *pt2,double *newpt2,double *oldnorm,double *newnorm,int sign); +void Sph_RotateVectorAxisAngle(const double *vect,const double *axis,double angle,double *rotated); #endif diff --git a/source/libSteve/queue.c b/source/libSteve/queue.c index 6e40ab3d..2003778a 100755 --- a/source/libSteve/queue.c +++ b/source/libSteve/queue.c @@ -9,12 +9,12 @@ of the Gnu Lesser General Public License (LGPL). */ #define CHECK(A) if(!(A)) {goto failure;} else (void)0 -queue q_alloc(int n,enum Q_types type,int (*keycmp)(void *,void *)) { - queue q; +q_queue q_alloc(int n,enum Q_types type,int (*keycmp)(void *,void *)) { + q_queue q; int i; if(n<0) return NULL; - q=(queue) malloc(sizeof(struct qstruct)); + q=(q_queue) malloc(sizeof(struct qstruct)); if(!q) return NULL; q->type=type; @@ -61,7 +61,7 @@ queue q_alloc(int n,enum Q_types type,int (*keycmp)(void *,void *)) { -int q_expand(queue q,int addspace) { +int q_expand(q_queue q,int addspace) { int i,j,num,nnew; void **xnew,**kvnew; int *kinew; @@ -135,7 +135,7 @@ int q_expand(queue q,int addspace) { -void q_free(queue q,int freek,int freex) { +void q_free(q_queue q,int freek,int freex) { int i; if(!q) return; @@ -152,12 +152,12 @@ void q_free(queue q,int freek,int freex) { -void q_null(queue q) { +void q_null(q_queue q) { q->f=q->b=0; } -int q_enqueue(void *kv,int ki,double kd,Q_LONGLONG kl,void *x,queue q) { +int q_enqueue(void *kv,int ki,double kd,Q_LONGLONG kl,void *x,q_queue q) { int sp; if(q->type==Qvoid) { @@ -179,7 +179,7 @@ int q_enqueue(void *kv,int ki,double kd,Q_LONGLONG kl,void *x,queue q) { -int q_push(void *kv,int ki,double kd,Q_LONGLONG kl,void *x,queue q) { +int q_push(void *kv,int ki,double kd,Q_LONGLONG kl,void *x,q_queue q) { int sp; q->f=(q->n+q->f-1)%q->n; @@ -201,7 +201,7 @@ int q_push(void *kv,int ki,double kd,Q_LONGLONG kl,void *x,queue q) { -int q_insert(void *kv,int ki,double kd,Q_LONGLONG kl,void *x,queue q) { +int q_insert(void *kv,int ki,double kd,Q_LONGLONG kl,void *x,q_queue q) { int i,n,sp,im1; n=q->n; @@ -256,7 +256,7 @@ int q_insert(void *kv,int ki,double kd,Q_LONGLONG kl,void *x,queue q) { -void q_front(queue q,void **kvptr,int *kiptr,double *kdptr,Q_LONGLONG *klptr,void **xptr) { +void q_front(q_queue q,void **kvptr,int *kiptr,double *kdptr,Q_LONGLONG *klptr,void **xptr) { if(q->f==q->b) { if(kvptr) *kvptr=NULL; if(kiptr) *kiptr=0; @@ -274,7 +274,7 @@ void q_front(queue q,void **kvptr,int *kiptr,double *kdptr,Q_LONGLONG *klptr,voi -int q_pop(queue q,void **kvptr,int *kiptr,double *kdptr,Q_LONGLONG *klptr,void **xptr) { +int q_pop(q_queue q,void **kvptr,int *kiptr,double *kdptr,Q_LONGLONG *klptr,void **xptr) { if(q->f==q->b) { if(kvptr) *kvptr=NULL; if(kiptr) *kiptr=0; @@ -294,17 +294,17 @@ int q_pop(queue q,void **kvptr,int *kiptr,double *kdptr,Q_LONGLONG *klptr,void * -int q_length(queue q) { +int q_length(q_queue q) { return (q->n+q->b-q->f)%q->n; } -int q_maxlength(queue q) { +int q_maxlength(q_queue q) { return q->n-1; } -int q_next(int i,void **kvptr,int *kiptr,double *kdptr,Q_LONGLONG *klptr,void **xptr,queue q) { +int q_next(int i,void **kvptr,int *kiptr,double *kdptr,Q_LONGLONG *klptr,void **xptr,q_queue q) { int f,b; f=q->f; diff --git a/source/libSteve/queue.h b/source/libSteve/queue.h index 80909c64..0d37b1e8 100755 --- a/source/libSteve/queue.h +++ b/source/libSteve/queue.h @@ -30,25 +30,25 @@ typedef struct qstruct{ void **x; int n; int f; - int b; } *queue; + int b; } *q_queue; #define q_frontkeyV(q) ((q)->b==(q)->f?NULL:(q)->kv[(q)->f]) #define q_frontkeyI(q) ((q)->b==(q)->f?0:(q)->ki[(q)->f]) #define q_frontkeyD(q) ((q)->b==(q)->f?0:(q)->kd[(q)->f]) #define q_frontkeyL(q) ((q)->b==(q)->f?0:(q)->kl[(q)->f]) -queue q_alloc(int n,enum Q_types type,int (*keycmp)(void *,void *)); -int q_expand(queue q,int addspace); -void q_free(queue q,int freek,int freex); -void q_null(queue q); -int q_enqueue(void *kv,int ki,double kd,Q_LONGLONG kl,void *x,queue q); -int q_push(void *kv,int ki,double kd,Q_LONGLONG kl,void *x,queue q); -int q_insert(void *kv,int ki,double kd,Q_LONGLONG kl,void *x,queue q); -void q_front(queue q,void **kvptr,int *kiptr,double *kdptr,Q_LONGLONG *klptr,void **xptr); -int q_pop(queue q,void **kvptr,int *kiptr,double *kdptr,Q_LONGLONG *klptr,void **xptr); -int q_length(queue q); -int q_maxlength(queue q); -int q_next(int i,void **kvptr,int *kiptr,double *kdptr,Q_LONGLONG *klptr,void **xptr,queue q); +q_queue q_alloc(int n,enum Q_types type,int (*keycmp)(void *,void *)); +int q_expand(q_queue q,int addspace); +void q_free(q_queue q,int freek,int freex); +void q_null(q_queue q); +int q_enqueue(void *kv,int ki,double kd,Q_LONGLONG kl,void *x,q_queue q); +int q_push(void *kv,int ki,double kd,Q_LONGLONG kl,void *x,q_queue q); +int q_insert(void *kv,int ki,double kd,Q_LONGLONG kl,void *x,q_queue q); +void q_front(q_queue q,void **kvptr,int *kiptr,double *kdptr,Q_LONGLONG *klptr,void **xptr); +int q_pop(q_queue q,void **kvptr,int *kiptr,double *kdptr,Q_LONGLONG *klptr,void **xptr); +int q_length(q_queue q); +int q_maxlength(q_queue q); +int q_next(int i,void **kvptr,int *kiptr,double *kdptr,Q_LONGLONG *klptr,void **xptr,q_queue q); #endif