Mercurial > octave-nkf
changeset 3452:a892190f4977
[project @ 2000-01-18 06:00:57 by jwe]
line wrap: on
line diff
--- a/doc/interpreter/Makefile.in Tue Jan 18 04:49:49 2000 +0000 +++ b/doc/interpreter/Makefile.in Tue Jan 18 06:01:04 2000 +0000 @@ -24,9 +24,9 @@ fn-idx.txi func.txi gpl.txi grammar.txi image.txi \ install.txi intro.txi io.txi linalg.txi matrix.txi \ nonlin.txi numbers.txi op-idx.txi optim.txi plot.txi \ - poly.txi preface.txi quad.txi set.txi signal.txi \ - stats.txi stmt.txi stream.txi strings.txi struct.txi \ - system.txi tips.txi var.txi vr-idx.txi + poly.txi preface.txi quad.txi quaternion.txi set.txi \ + signal.txi stats.txi stmt.txi stream.txi strings.txi \ + struct.txi system.txi tips.txi var.txi vr-idx.txi SOURCES := $(SUB_SOURCE)
--- a/doc/interpreter/audio.txi Tue Jan 18 04:49:49 2000 +0000 +++ b/doc/interpreter/audio.txi Tue Jan 18 06:01:04 2000 +0000 @@ -3,7 +3,7 @@ @c This is part of the Octave manual. @c For copying conditions, see the file gpl.texi. -@node Audio Processing, System Utilities, Image Processing, Top +@node Audio Processing, Quaternions, Image Processing, Top @chapter Audio Processing Octave provides a few functions for dealing with audio data. An audio
--- a/doc/interpreter/octave.texi Tue Jan 18 04:49:49 2000 +0000 +++ b/doc/interpreter/octave.texi Tue Jan 18 06:01:04 2000 +0000 @@ -40,7 +40,7 @@ @defindex op -@c Things like the Octave version number are defined in conf.texi. +@c Things like the Octave version number are defined in conf.txi. @c This file doesn't include a chapter, so it must not be included @c if you want to run the Emacs function texinfo-multiple-files-update. @@ -146,6 +146,7 @@ * Signal Processing:: * Image Processing:: * Audio Processing:: +* Quaternions:: * System Utilities:: * Tips:: * Trouble:: If you have trouble installing Octave. @@ -164,316 +165,6 @@ Preface -* Acknowledgements:: -* How You Can Contribute to Octave:: -* Distribution:: - -A Brief Introduction to Octave - -* Running Octave:: -* Simple Examples:: -* Conventions:: - -Conventions - -* Fonts:: -* Evaluation Notation:: -* Printing Notation:: -* Error Messages:: -* Format of Descriptions:: - -Format of Descriptions - -* A Sample Function Description:: -* A Sample Command Description:: -* A Sample Variable Description:: - -Getting Started - -* Invoking Octave:: -* Quitting Octave:: -* Getting Help:: -* Command Line Editing:: -* Errors:: -* Executable Octave Programs:: -* Comments:: - -Invoking Octave - -* Command Line Options:: -* Startup Files:: - -Command Line Editing - -* Cursor Motion:: -* Killing and Yanking:: -* Commands For Text:: -* Commands For Completion:: -* Commands For History:: -* Customizing readline:: -* Customizing the Prompt:: -* Diary and Echo Commands:: - -Data Types - -* Built-in Data Types:: -* User-defined Data Types:: -* Object Sizes:: - -Built-in Data Types - -* Numeric Objects:: -* String Objects:: -* Data Structure Objects:: - -Numeric Data Types - -* Matrices:: -* Ranges:: -* Logical Values:: -* Predicates for Numeric Objects:: - -Matrices - -* Empty Matrices:: - -Strings - -* Creating Strings:: -* Searching and Replacing:: -* String Conversions:: -* Character Class Functions:: - -Containers - -* Lists:: -* Cell Arrays:: - -Variables - -* Global Variables:: -* Status of Variables:: -* Summary of Built-in Variables:: -* Defaults from the Environment:: - -Expressions - -* Index Expressions:: -* Calling Functions:: -* Arithmetic Ops:: -* Comparison Ops:: -* Boolean Expressions:: -* Assignment Ops:: -* Increment Ops:: -* Operator Precedence:: - -Calling Functions - -* Call by Value:: -* Recursion:: - -Boolean Expressions - -* Element-by-element Boolean Operators:: -* Short-circuit Boolean Operators:: - -Statements - -* The if Statement:: -* The switch Statement:: -* The while Statement:: -* The for Statement:: -* The break Statement:: -* The continue Statement:: -* The unwind_protect Statement:: -* The try Statement:: -* Continuation Lines:: - -The @code{for} Statement - -* Looping Over Structure Elements:: - -Functions and Script Files - -* Defining Functions:: -* Multiple Return Values:: -* Variable-length Argument Lists:: -* Variable-length Return Lists:: -* Returning From a Function:: -* Function Files:: -* Script Files:: -* Dynamically Linked Functions:: -* Organization of Functions:: - -Input and Output - -* Basic Input and Output:: -* C-Style I/O Functions:: - -Basic Input and Output - -* Terminal Output:: -* Terminal Input:: -* Simple File I/O:: - -C-Style I/O Functions - -* Opening and Closing Files:: -* Simple Output:: -* Line-Oriented Input:: -* Formatted Output:: -* Output Conversion for Matrices:: -* Output Conversion Syntax:: -* Table of Output Conversions:: -* Integer Conversions:: -* Floating-Point Conversions:: Other Output Conversions:: -* Other Output Conversions:: -* Formatted Input:: -* Input Conversion Syntax:: -* Table of Input Conversions:: -* Numeric Input Conversions:: -* String Input Conversions:: -* Binary I/O:: -* Temporary Files:: -* EOF and Errors:: -* File Positioning:: - -Plotting - -* Two-Dimensional Plotting:: -* Specialized Two-Dimensional Plots:: -* Three-Dimensional Plotting:: -* Plot Annotations:: -* Multiple Plots on One Page:: -* Multiple Plot Windows:: -* Interaction with gnuplot:: - -Matrix Manipulation - -* Finding Elements and Checking Conditions:: -* Rearranging Matrices:: -* Special Utility Matrices:: -* Famous Matrices:: - -Arithmetic - -* Utility Functions:: -* Complex Arithmetic:: -* Trigonometry:: -* Sums and Products:: -* Special Functions:: -* Mathematical Constants:: - -Linear Algebra - -* Basic Matrix Functions:: -* Matrix Factorizations:: -* Functions of a Matrix:: - -Quadrature - -* Functions of One Variable:: -* Orthogonal Collocation:: - -Differential Equations - -* Ordinary Differential Equations:: -* Differential-Algebraic Equations:: - -Optimization - -* Quadratic Programming:: -* Nonlinear Programming:: -* Linear Least Squares:: - -Control Theory - -* sysstruct:: -* sysinterface:: -* sysdisp:: -* blockdiag:: -* numerical:: -* sysprop:: -* systime:: -* sysfreq:: -* cacsd:: -* misc:: - -System Data Structure - -* sysstructvars:: -* sysstructtf:: -* sysstructzp:: -* sysstructss:: - -System Construction and Interface Functions - -* fir2sys:: -* ss2sys:: -* tf2sys:: -* zp2sys:: -* structaccess:: -* structintern:: - -System Utilities - -* Timing Utilities:: -* Filesystem Utilities:: -* Controlling Subprocesses:: -* Process ID Information:: -* Environment Variables:: -* Current Working Directory:: -* Password Database Functions:: -* Group Database Functions:: -* System Information:: - -Tips and Standards - -* Style Tips:: Writing clean and robust programs. -* Coding Tips:: Making code run faster. -* Documentation Tips:: Writing readable documentation strings. -* Comment Tips:: Conventions for writing comments. -* Function Headers:: Standard headers for functions. - -Known Causes of Trouble with Octave - -* Actual Bugs:: Bugs we will fix later. -* Reporting Bugs:: -* Bug Criteria:: -* Bug Lists:: -* Bug Reporting:: -* Sending Patches:: -* Service:: - -Reporting Bugs - -* Bug Criteria:: -* Where: Bug Lists. Where to send your bug report. -* Reporting: Bug Reporting. How to report a bug effectively. -* Patches: Sending Patches. How to send a patch for Octave. - -Installing Octave - -* Installation Problems:: -* Binary Distributions:: - -Binary Distributions - -* Installing Octave from a Binary Distribution:: -* Creating a Binary Distribution:: - -Emacs Octave Support - -* Installing EOS:: -* Using Octave Mode:: -* Running Octave From Within Emacs:: -* Using the Emacs Info Reader for Octave:: - -Grammar - -* Keywords:: - -@end detailmenu @end menu @include preface.texi @@ -508,6 +199,7 @@ @include signal.texi @include image.texi @include audio.texi +@include quaternion.texi @include system.texi @c maybe add again later, if anyone every writes any really interesting
--- a/doc/interpreter/system.txi Tue Jan 18 04:49:49 2000 +0000 +++ b/doc/interpreter/system.txi Tue Jan 18 06:01:04 2000 +0000 @@ -2,7 +2,7 @@ @c This is part of the Octave manual. @c For copying conditions, see the file gpl.texi. -@node System Utilities, Tips, Audio Processing, Top +@node System Utilities, Tips, Quaternions, Top @chapter System Utilities This chapter describes the functions that are available to allow you to
--- a/scripts/ChangeLog Tue Jan 18 04:49:49 2000 +0000 +++ b/scripts/ChangeLog Tue Jan 18 06:01:04 2000 +0000 @@ -61,6 +61,10 @@ * signal/yulewalker.m: Ditto. * control/util/strappend.m: Ditto. * control/base/nichols.m: Ditto. + * control/system/is_signal_list.m: Ditto. + * control/system/listidx.m: Ditto. + * control/system/sysgettsam.m: Ditto. + * control/system/sysidx.m: Ditto. 2000-01-14 John W. Eaton <jwe@bevo.che.wisc.edu> @@ -78,6 +82,17 @@ * control/system/__tf2sysl__.m: Likewise. * control/util/__zgpbal__.m: Likewise. * control/system/__zp2ssg2__.m: Likewise. + * quaternion/demoquat.m: Ditto. + * quaternion/qconj.m: Ditto. + * quaternion/qcoordinate_plot.m: Ditto. + * quaternion/qderiv.m: Ditto. + * quaternion/qderivmat.m: Ditto. + * quaternion/qinv.m: Ditto. + * quaternion/qmult.m: Ditto. + * quaternion/qtrans.m: Ditto. + * quaternion/qtransv.m: Ditto. + * quaternion/qtransvmat.m: Ditto. + * quaternion/quaternion.m: Ditto. 2000-01-13 John W. Eaton <jwe@bevo.che.wisc.edu>
--- a/scripts/control/hinf/dhinfdemo.m Tue Jan 18 04:49:49 2000 +0000 +++ b/scripts/control/hinf/dhinfdemo.m Tue Jan 18 06:01:04 2000 +0000 @@ -16,26 +16,26 @@ ## along with Octave; see the file COPYING. If not, write to the Free ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. -## ------------------------------------------------------------ -## dhinfdemo Design of a discrete H_infinity controller. -## This is not a true discrete design. The design -## is carried out in continuous time while the -## effect of sampling is described by a bilinear -## transformation of the sampled system. -## This method works quite well if the sampling -## period is "small" compared to the plant time -## constants. +## -*- texinfo -*- +## @deftypefn {Function File} {} dhinfdemo () +## Demonstrate the functions available for designining a discrete +## H_infinity controller. This is not a true discrete design. The +## design is carried out in continuous time while the effect of sampling +## is described by a bilinear transformation of the sampled system. +## This method works quite well if the sampling period is "small" +## compared to the plant time constants. ## -## This is a script file for OCTAVE. -## ------------------------------------------------------------ +## Continuous plant: ## -## continuous plant: +## @example ## 1 ## G(s) = -------------- ## (s + 2)(s + 1) +## @end example ## -## discretised plant with ZOH (Sampling period = Ts = 1 second) +## Discretised plant with ZOH (Sampling period = Ts = 1 second): ## +## @example ## 0.39958z + 0.14700 ## G(s) = -------------------------- ## (z - 0.36788)(z - 0.13533) @@ -52,9 +52,11 @@ ## | +---+ | ## -----| K |<------- ## +---+ +## @end example ## -## W1 and W2 are the robustness and performancs weighting -## functions +## @noindent +## W1 and W2 are the robustness and performancs weighting functions. +## @end deftypefn ## K. Mueller, <mueller@ifr.ing.tu-bs.de> ## Technical University of Braunschweig, IfR
--- a/scripts/control/system/is_signal_list.m Tue Jan 18 04:49:49 2000 +0000 +++ b/scripts/control/system/is_signal_list.m Tue Jan 18 06:01:04 2000 +0000 @@ -16,8 +16,10 @@ ## along with Octave; see the file COPYING. If not, write to the Free ## Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. -## function flg = is_signal_list (mylist) -## returns true if mylist is a list of individual strings. +## -*- texiofn -*- +## @deftypefn {Function File} {} is_signal_list (@var{mylist}) +## Return true if @var{mylist} is a list of individual strings. +## @end deftypefn function flg = is_signal_list(mylist)
--- a/scripts/control/system/listidx.m Tue Jan 18 04:49:49 2000 +0000 +++ b/scripts/control/system/listidx.m Tue Jan 18 06:01:04 2000 +0000 @@ -16,25 +16,22 @@ ## along with Octave; see the file COPYING. If not, write to the Free ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. -## [idxvec, errmsg] = listidx(listvar, strlist) -## return indices of string entries in listvar that match strings in strlist -## Inputs: -## listvar: list of strings to be searched -## strlist: list of strings to be located in listvar. -## Note: listvar, strlist may be passed as strings or string matrices; -## in this case, each entry is processed by deblank() prior to searching -## for the entries of strlist in listvar. -## Outputs: -## idxvec -## vector of indices in listvar; -## listvar(idxvec(k)) == strlist(kk). -## errmsg -## if strlist contains a string not in listvar, then -## an error message is returned in errmsg. If only one output -## argument is requested, e.g., idxvec = listidx(listvar, strlist), -## then listidx prints errmsg to the screen and exits with -## an error. +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{idxvec}, @var{errmsg}] =} listidx (@var{listvar}, @var{strlist}) +## Return indices of string entries in @var{listvar} that match strings +## in @var{strlist}. ## +## Both @var{listvar} and @var{strlist} may be passed as strings or +## string matrices. If they are passed as string matrices, each entry +## is processed by @code{deblank} prior to searching for the entries. +## +## The first output is the vector of indices in @var{listvar}. +## +## If @var{strlist} contains a string not in @var{listvar}, then +## an error message is returned in @var{errmsg}. If only one output +## argument is requested, then @var{listidx} prints @var{errmsg} to the +## screen and exits with an error. +## @end deftypefn function [idxvec,errmsg] = listidx(listvar,strlist)
--- a/scripts/control/system/sysgettsam.m Tue Jan 18 04:49:49 2000 +0000 +++ b/scripts/control/system/sysgettsam.m Tue Jan 18 06:01:04 2000 +0000 @@ -16,8 +16,10 @@ ## along with Octave; see the file COPYING. If not, write to the Free ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. -## T = sysgettsam(sys) -## return the sampling time of the system +## -*- texinfo -*- +## @deftypefn {Function File} {} sysgettsam (@var{sys}) +## Return the sampling time of the system @var{sys}. +## @end deftypefn function T = sysgettsam (sys)
--- a/scripts/control/system/sysidx.m Tue Jan 18 04:49:49 2000 +0000 +++ b/scripts/control/system/sysidx.m Tue Jan 18 06:01:04 2000 +0000 @@ -16,14 +16,13 @@ ## along with Octave; see the file COPYING. If not, write to the Free ## Software Foundation, 59 Temple Place, Suite 330, Boston, MA 02111 USA. -## idxvec = sysidx (sys, sigtype, signamelist) -## return indices of signals with specified signal names -## inputs: -## sys: OCST system data structure -## sigtype: signal type to be selected: "in", "out", "st" -## signamelist: list of desired signal names -## outputs: -## idxvec: vector of signal indices (appropriate for use with sysprune) +## -*- texinfo -*- +## @deftypefn {Function File} {} sysidx (@var{sys}, @var{sigtype}, @var{signamelist}) +## Return indices of signals with specified signal names +## inputs given a system data structure @var{sys}, a signal type to be +## selected @var{sigtype} (@code{"in"}, @code{"out"}, @code{"st"}), and +## a list of desired signal names @var{signamelist}. +## @end deftypefn function idxvec = sysidx (sys, sigtype, signamelist)
--- a/scripts/quaternion/demoquat.m Tue Jan 18 04:49:49 2000 +0000 +++ b/scripts/quaternion/demoquat.m Tue Jan 18 06:01:04 2000 +0000 @@ -1,233 +1,267 @@ - -# Thanks to Mr. Charles Hall, Dr. Don Krupp and Dr. Larry Mullins at -# NASA's Marshall Space Flight Center for notes and instruction in -# use and conventions with quaternions. - A. S. Hodel - -function opt = demoquat() +## Copyright (C) 1998 Auburn University. All rights reserved. +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, or (at your option) +## any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, 59 Temple Place - Suite 330, Boston, MA +## 02111-1307, USA. -opt = 0; -quitopt = 5; -while(opt != quitopt) - opt = menu("Quaternion function demo (c) 1998 A. S. Hodel, a.s.hodel@eng.auburn.edu", - "quaternion construction/data extraction", - "simple quaternion functions", - "transformation functions", - "body-inertial frame demo", - "Quit"); +## -*- texinfo -*- +## @deftypefn {Function File} {} demoquat () +## Demonstrate the functions available for manipulating quaternions. +## +## Thanks to Mr. Charles Hall, Dr. Don Krupp and Dr. Larry Mullins at +## NASA's Marshall Space Flight Center for notes and instruction in +## use and conventions with quaternions. - A. S. Hodel +## @end deftypefn + +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Adapted-By: jwe + +function opt = demoquat () + + opt = 0; + quitopt = 5; + + while (opt != quitopt) + opt = menu ("Quaternion function demo (c) 1998 A. S. Hodel, a.s.hodel@eng.auburn.edu", + "quaternion construction/data extraction", + "simple quaternion functions", + "transformation functions", + "body-inertial frame demo", + "Quit"); + + switch(opt) - switch(opt) - case(1), - printf("Quaternion construction/data extraction\n"); - help quaternion - prompt - cmd = "q = quaternion(1,2,3,4)"; - run_cmd - disp("This format stores the i,j,k parts of the quaternion first;") - disp("the real part is stored last.") - prompt - disp(" ") - disp("i, j, and k are all square roots of -1; however they do not") - disp("commute under multiplication (discussed further with the function") - disp("qmult). Therefore quaternions do not commute under multiplcation:") - disp(" q1*q2 != q2*q1 (usually)") - prompt + case(1) + printf("Quaternion construction/data extraction\n"); + help quaternion + prompt + cmd = "q = quaternion(1,2,3,4)"; + run_cmd + disp("This format stores the i,j,k parts of the quaternion first;") + disp("the real part is stored last.") + prompt + disp(" ") + disp("i, j, and k are all square roots of -1; however they do not") + disp("commute under multiplication (discussed further with the function") + disp("qmult). Therefore quaternions do not commute under multiplcation:") + disp(" q1*q2 != q2*q1 (usually)") + prompt - disp("Quaternions as rotations: unit quaternion to represent") - disp("rotation of 45 degrees about the vector [1 1 1]") - cmd = "degrees = pi/180; q1 = quaternion([1 1 1],45*degrees)"; - run_cmd - prompt - cmd = "real_q = cos(45*degrees/2)"; - run_cmd - printf("The real part of the quaternion q(4) is cos(theta/2).\n----\n\n"); - cmd = "imag_q = sin(45*degrees/2)*[1 1 1]/norm([1 1 1])" - run_cmd - disp("The imaginary part of the quaternion is sin(theta/2)*unit vector"); - disp("The constructed quaternion is a unit quaternion."); - prompt - disp("Can also extract both forms of the quaternion:") - disp("Vector/angle form of 1i + 2j + 3k + 4:") - cmd = "[vv,th] = quaternion(q)"; - run_cmd - cmd = "vv_norm = norm(vv)"; - run_cmd - disp("Returns the eigenaxis as a 3-d unit vector"); - disp("Check values: ") - cmd = "th_deg = th*180/pi"; - run_cmd - disp("") - disp("This concludes the quaternion construction/extraction demo."); - prompt + disp("Quaternions as rotations: unit quaternion to represent") + disp("rotation of 45 degrees about the vector [1 1 1]") + cmd = "degrees = pi/180; q1 = quaternion([1 1 1],45*degrees)"; + run_cmd + prompt + cmd = "real_q = cos(45*degrees/2)"; + run_cmd + printf("The real part of the quaternion q(4) is cos(theta/2).\n----\n\n"); + cmd = "imag_q = sin(45*degrees/2)*[1 1 1]/norm([1 1 1])" + run_cmd + disp("The imaginary part of the quaternion is sin(theta/2)*unit vector"); + disp("The constructed quaternion is a unit quaternion."); + prompt + disp("Can also extract both forms of the quaternion:") + disp("Vector/angle form of 1i + 2j + 3k + 4:") + cmd = "[vv,th] = quaternion(q)"; + run_cmd + cmd = "vv_norm = norm(vv)"; + run_cmd + disp("Returns the eigenaxis as a 3-d unit vector"); + disp("Check values: ") + cmd = "th_deg = th*180/pi"; + run_cmd + disp("") + disp("This concludes the quaternion construction/extraction demo."); + prompt + + case(2) + printf("Simple quaternion functions\n"); + cmd = "help qconj"; + run_cmd + cmd = "degrees = pi/180; q1 = quaternion([1 1 1],45*degrees)"; + run_cmd + cmd = "q2 = qconj(q1)"; + run_cmd + disp("The conjugate changes the sign of the complex part of the") + printf("quaternion.\n\n"); + prompt + printf("\n\n\nMultiplication of quaternions:\n"); + cmd = "help qmult"; + run_cmd + cmd = "help qinv" + run_cmd + disp("Inverse quaternion: q*qi = qi*q = 1:") + cmd = "q1i = qinv(q1)"; + run_cmd + cmd = "one = qmult(q1,q1i)"; + run_cmd + + printf("Conclusion of simple quaternion functions"); + prompt - case(2), - printf("Simple quaternion functions\n"); - cmd = "help qconj"; - run_cmd - cmd = "degrees = pi/180; q1 = quaternion([1 1 1],45*degrees)"; - run_cmd - cmd = "q2 = qconj(q1)"; - run_cmd - disp("The conjugate changes the sign of the complex part of the") - printf("quaternion.\n\n"); - prompt - printf("\n\n\nMultiplication of quaternions:\n"); - cmd = "help qmult"; - run_cmd - cmd = "help qinv" - run_cmd - disp("Inverse quaternion: q*qi = qi*q = 1:") - cmd = "q1i = qinv(q1)"; - run_cmd - cmd = "one = qmult(q1,q1i)"; - run_cmd - - printf("Conclusion of simple quaternion functions"); - prompt - case(3), - printf("Transformation functions\n"); - disp("A problem with the discussion of coordinate transformations is that"); - disp("one must be clear on what is being transformed: does a rotation of"); - disp("theta degrees mean that you're rotating the VECTOR by theta degrees,"); - disp("also called the 'active convention,' or does it mean that you rotate "); - disp("the COORDINATE FRAME by theta degrees, also called the 'passive convention,' "); - disp("which is equivalent to rotating the VECTOR by (-theta) degrees. The"); - disp("functions in this demo use the active convention. I'll point out where"); - disp("this changes the code as the demo runs."); - disp(" -- The author"); - prompt - printf("\n\n"); - disp("Sequences of rotations:") - printf("\n\nRotation of a vector by 90 degrees about the reference z axis\n"); - cmd = "qz = quaternion([0 0 1], pi/2);"; - disp(cmd) ; eval(cmd); - printf("\n\nRotation of a vector by 90 degrees about the reference y axis\n"); - cmd="qy = quaternion([0 1 0], pi/2);"; - disp(cmd) ; eval(cmd); - printf("\n\nRotation of a vector by 90 degrees about the reference x axis\n"); - cmd="qx = quaternion([1 0 0], pi/2);"; - run_cmd - printf("\n\nSequence of three rotations: 90 degrees about x, then 90 degrees\n"); - disp("about y, then 90 degrees about z (all axes specified in the reference frame):"); - qchk = qmult(qz,qmult(qy,qx)); - cmd = "[vv,th] = quaternion(qchk), th_deg = th*180/pi"; - run_cmd - disp("The sequence of the three rotations above is equivalent to a single rotation") - disp("of 90 degrees about the y axis. Check:"); - cmd = "err = norm(qchk - qy)"; - run_cmd - - disp("Transformation of a quaternion by a quaternion:") - disp("The three quaternions above were rotations specified about") - disp("a single reference frame. It is often convenient to specify the"); - disp("eigenaxis of a rotation in a different frame (e.g., when computing"); - disp("the transformation rotation in terms of the Euler angles yaw-pitch-roll)."); - cmd = "help qtrans"; - run_cmd - disp("") - disp("NOTE: If the passive convention is used, then the above"); - disp("formula changes to v = qinv(q)*v*q instead of ") - disp("v = q*v*qinv(q).") - prompt - disp("") - disp("Example: Vectors in Frame 2 are obtained by rotating them from ") - disp(" from Frame 1 by 90 degrees about the x axis (quaternion qx)") - disp(" A quaternion in Frame 2 rotates a vector by 90 degrees about") - disp(" the Frame 2 y axis (quaternion qy). The equivalent rotation") - disp(" in the reference frame is:") - cmd = "q_eq = qtrans(qy,qx); [vv,th] = quaternion(q_eq)"; - run_cmd - disp("The rotation is equivalent to rotating about the reference z axis") - disp("by 90 degrees (quaternion qz)") - prompt + case(3) + printf("Transformation functions\n"); + disp("A problem with the discussion of coordinate transformations is that"); + disp("one must be clear on what is being transformed: does a rotation of"); + disp("theta degrees mean that you're rotating the VECTOR by theta degrees,"); + disp("also called the 'active convention,' or does it mean that you rotate "); + disp("the COORDINATE FRAME by theta degrees, also called the 'passive convention,' "); + disp("which is equivalent to rotating the VECTOR by (-theta) degrees. The"); + disp("functions in this demo use the active convention. I'll point out where"); + disp("this changes the code as the demo runs."); + disp(" -- The author"); + prompt + printf("\n\n"); + disp("Sequences of rotations:") + printf("\n\nRotation of a vector by 90 degrees about the reference z axis\n"); + cmd = "qz = quaternion([0 0 1], pi/2);"; + disp(cmd) ; eval(cmd); + printf("\n\nRotation of a vector by 90 degrees about the reference y axis\n"); + cmd="qy = quaternion([0 1 0], pi/2);"; + disp(cmd) ; eval(cmd); + printf("\n\nRotation of a vector by 90 degrees about the reference x axis\n"); + cmd="qx = quaternion([1 0 0], pi/2);"; + run_cmd + printf("\n\nSequence of three rotations: 90 degrees about x, then 90 degrees\n"); + disp("about y, then 90 degrees about z (all axes specified in the reference frame):"); + qchk = qmult(qz,qmult(qy,qx)); + cmd = "[vv,th] = quaternion(qchk), th_deg = th*180/pi"; + run_cmd + disp("The sequence of the three rotations above is equivalent to a single rotation") + disp("of 90 degrees about the y axis. Check:"); + cmd = "err = norm(qchk - qy)"; + run_cmd - disp("Transformation of a vector by a quaternion"); - cmd = "help qtransv"; - run_cmd + disp("Transformation of a quaternion by a quaternion:") + disp("The three quaternions above were rotations specified about") + disp("a single reference frame. It is often convenient to specify the"); + disp("eigenaxis of a rotation in a different frame (e.g., when computing"); + disp("the transformation rotation in terms of the Euler angles yaw-pitch-roll)."); + cmd = "help qtrans"; + run_cmd + disp("") + disp("NOTE: If the passive convention is used, then the above"); + disp("formula changes to v = qinv(q)*v*q instead of ") + disp("v = q*v*qinv(q).") + prompt + disp("") + disp("Example: Vectors in Frame 2 are obtained by rotating them from ") + disp(" from Frame 1 by 90 degrees about the x axis (quaternion qx)") + disp(" A quaternion in Frame 2 rotates a vector by 90 degrees about") + disp(" the Frame 2 y axis (quaternion qy). The equivalent rotation") + disp(" in the reference frame is:") + cmd = "q_eq = qtrans(qy,qx); [vv,th] = quaternion(q_eq)"; + run_cmd + disp("The rotation is equivalent to rotating about the reference z axis") + disp("by 90 degrees (quaternion qz)") + prompt + + disp("Transformation of a vector by a quaternion"); + cmd = "help qtransv"; + run_cmd + + disp("NOTE: the above formula changes if the passive quaternion ") + disp("is used; the cross product term is subtracted instead of added."); + prompt + disp("Example: rotate the vector [1,1,1] by 90 degrees about the y axis"); + cmd = "vec_r = qtransv([1,1,1],qy)"; + run_cmd + prompt + disp("Equivalently, one may multiply by qtransvmat:") + cmd = "help qtransvmat"; + run_cmd + disp("NOTE: the passive quaternion convention would use the transpose") + disp("(inverse) of the orthogonal matrix returned by qtransvmat."); + prompt + cmd = "vec_r_2 = qtransvmat(qy)*[1;1;1]; vec_err = norm(vec_r - vec_r_2)"; + run_cmd - disp("NOTE: the above formula changes if the passive quaternion ") - disp("is used; the cross product term is subtracted instead of added."); - prompt - disp("Example: rotate the vector [1,1,1] by 90 degrees about the y axis"); - cmd = "vec_r = qtransv([1,1,1],qy)"; - run_cmd - prompt - disp("Equivalently, one may multiply by qtransvmat:") - cmd = "help qtransvmat"; - run_cmd - disp("NOTE: the passive quaternion convention would use the transpose") - disp("(inverse) of the orthogonal matrix returned by qtransvmat."); - prompt - cmd = "vec_r_2 = qtransvmat(qy)*[1;1;1]; vec_err = norm(vec_r - vec_r_2)"; - run_cmd + disp("") + disp("The last transformation function is the derivative of a quaternion") + disp("Given rotation rates about the reference x, y, and z axes."); + cmd = "help qderivmat"; + run_cmd + disp("") + disp("Example:") + disp("Frame is rotating about the z axis at 1 rad/s") + cmd = "Omega = [0,0,1]; Dmat = qderivmat(Omega)"; + run_cmd + disp("Notice that Dmat is skew symmetric, as it should be.") + disp("expm(Dmat*t) is orthogonal, so that unit quaternions remain") + disp("unit quaternions as the rotating frame precesses."); + disp(" ") + disp("This concludes the transformation demo."); + prompt; + + case(4) + printf("Body-inertial frame demo: Look at the source code for\n"); + printf("demoquat.m and qcoordinate_plot.m to see how it's done.\n"); + + # i,j,k units + iv = quaternion(1,0,0,0); jv = quaternion(0,1,0,0); + kv = quaternion(0,0,1,0); + + # construct quaternion to desired view. + degrees = pi/180; daz = 45*degrees; del = -30*degrees; + qazimuth = quaternion([0,0,1],daz); + qelevation = quaternion([cos(daz),sin(daz),0],del); + qview = qmult(qelevation,qazimuth); + + # inertial frame i, j, k axes. + iif = iv; jf = qtrans(jv,iv); kf = qtrans(kv,iv); + + # rotation steps + th = 0:5:20; ov = ones(size(th)); myth = [th,max(th)*ov ; 0*ov,th]; - disp("") - disp("The last transformation function is the derivative of a quaternion") - disp("Given rotation rates about the reference x, y, and z axes."); - cmd = "help qderivmat"; - run_cmd - disp("") - disp("Example:") - disp("Frame is rotating about the z axis at 1 rad/s") - cmd = "Omega = [0,0,1]; Dmat = qderivmat(Omega)"; - run_cmd - disp("Notice that Dmat is skew symmetric, as it should be.") - disp("expm(Dmat*t) is orthogonal, so that unit quaternions remain") - disp("unit quaternions as the rotating frame precesses."); - disp(" ") - disp("This concludes the transformation demo."); - prompt; - case(4), - printf("Body-inertial frame demo: Look at the source code for\n"); - printf("demoquat.m and qcoordinate_plot.m to see how it's done.\n"); - - # i,j,k units - iv = quaternion(1,0,0,0); jv = quaternion(0,1,0,0); - kv = quaternion(0,0,1,0); - - # construct quaternion to desired view. - degrees = pi/180; daz = 45*degrees; del = -30*degrees; - qazimuth = quaternion([0,0,1],daz); - qelevation = quaternion([cos(daz),sin(daz),0],del); - qview = qmult(qelevation,qazimuth); - - # inertial frame i, j, k axes. - iif = iv; jf = qtrans(jv,iv); kf = qtrans(kv,iv); - - # rotation steps - th = 0:5:20; ov = ones(size(th)); myth = [th,max(th)*ov ; 0*ov,th]; - - # construct yaw-pitch-roll cartoon - for kk=1:length(myth(1,:)) - figure(kk) - thy = myth(1,kk); - thp = myth(2,kk); - - qyaw = quaternion([0,0,1],thy*pi/180); - [jvy,th] = quaternion(qtrans(jf,qyaw)); - qpitch = quaternion(jvy(1:3),thp*pi/180); - qb = qmult(qpitch, qyaw); - qi = quaternion([1, 0, 0],180*degrees); - - printf("yaw=%8.4f, pitch=%8.4f, \n qbi = (%8.4f)i + (%8.4e)j + (%8.4f)k + (%8.4f)\n",thy,thp, ... - qb(1), qb(2), qb(3), qb(4)); - [vv,th] = quaternion(qb); - printf(" = (vector) = [%8.4f %8.4f %8.4f], th=%5.2f deg\n", ... - vv(1), vv(2), vv(3), th*180/pi); - - qb = qmult(qb,qi); - title(sprintf("yaw=%5.2f deg, pitch=%5.2f deg",thy,thp)) - qcoordinate_plot(qi,qb,qview); - # uncomment the next four lines to get eps output - #gset terminal postscript eps - #eval(sprintf("gset output 'fig%d.eps'",kk)); - #replot - #gset terminal x11 - #prompt - endfor - case(quitopt), - printf("Exiting quaternion demo\n"); - otherwise, - error ("invalid option %f", opt); - endswitch -endwhile + # construct yaw-pitch-roll cartoon + for kk=1:length(myth(1,:)) + figure(kk) + thy = myth(1,kk); + thp = myth(2,kk); + + qyaw = quaternion([0,0,1],thy*pi/180); + [jvy,th] = quaternion(qtrans(jf,qyaw)); + qpitch = quaternion(jvy(1:3),thp*pi/180); + qb = qmult(qpitch, qyaw); + qi = quaternion([1, 0, 0],180*degrees); + + printf("yaw=%8.4f, pitch=%8.4f, \n qbi = (%8.4f)i + (%8.4e)j + (%8.4f)k + (%8.4f)\n",thy,thp, ... + qb(1), qb(2), qb(3), qb(4)); + [vv,th] = quaternion(qb); + printf(" = (vector) = [%8.4f %8.4f %8.4f], th=%5.2f deg\n", ... + vv(1), vv(2), vv(3), th*180/pi); + + qb = qmult(qb,qi); + title(sprintf("yaw=%5.2f deg, pitch=%5.2f deg",thy,thp)) + qcoordinate_plot(qi,qb,qview); + # uncomment the next four lines to get eps output + #gset terminal postscript eps + #eval(sprintf("gset output 'fig%d.eps'",kk)); + #replot + #gset terminal x11 + #prompt + endfor + + case(quitopt) + printf ("Exiting quaternion demo\n"); + + otherwise + error ("invalid option %f", opt); + + endswitch + endwhile + endfunction
--- a/scripts/quaternion/qconj.m Tue Jan 18 04:49:49 2000 +0000 +++ b/scripts/quaternion/qconj.m Tue Jan 18 06:01:04 2000 +0000 @@ -1,11 +1,40 @@ +## Copyright (C) 1998 Auburn University. All rights reserved. +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, or (at your option) +## any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, 59 Temple Place - Suite 330, Boston, MA +## 02111-1307, USA. -function retval = qconj(q) - # function retval = qconj(q) - # conjugate of a quaternion - # q = [w,x,y,z] = w*i + x*j + y*k + z - # qconj(q) = -w*i -x*j -y*k + z +## -*- texinfo -*- +## @deftypefn {Function File} {} qconj (@var{q}) +## Conjugate of a quaternion. +## +## @example +## q = [w, x, y, z] = w*i + x*j + y*k + z +## qconj (q) = -w*i -x*j -y*k + z +## @end example +## @end deftypefn - [a,b,c,d] = quaternion(q); - retval = quaternion(-a,-b,-c,d); +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Adapted-By: jwe + +function retval = qconj (q) + + [a, b, c, d] = quaternion (q); + + retval = quaternion (-a, -b, -c, d); + endfunction
--- a/scripts/quaternion/qcoordinate_plot.m Tue Jan 18 04:49:49 2000 +0000 +++ b/scripts/quaternion/qcoordinate_plot.m Tue Jan 18 06:01:04 2000 +0000 @@ -1,64 +1,130 @@ -function qcoordinate_plot(qf,qb,qv) -# function qcoordinate_plot(qf,qb,qv) -# plot in the current figure a set of coordinate axes as viewed from -# the orientation specified by quaternion qv. Inertial axes are -# also plotted -# qf: quaternion from reference (x,y,z) to inertial -# qb: quaternion from reference to body -# qv: quaternion from reference to view angle - -degrees = pi/180; -d180 = 180*degrees; +## Copyright (C) 1998 Auburn University. All rights reserved. +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, or (at your option) +## any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, 59 Temple Place - Suite 330, Boston, MA +## 02111-1307, USA. -# construct coordinate transformation to view frame -cm = qtransvmat(qv); -p1 = [-1,-1,1]; p2 = [-1,-1,-1]; p3 = [1,-1,-1]; p4 = [ 1,-1, 1]; -p5 = [-1, 1,1]; p6 = [ 1, 1, 1]; p7 = [1, 1,-1]; p8 = [-1, 1,-1]; -# outline positive quadrant -box1 = cm*[p4; p6; p5; p6; p7]'; -# outline rest of the box -box2 =cm*[p7; p8; p5; p1; p4; p3; p7; p3; p2; p1; p2; p8]'; +## -*- texinfo -*- +## @deftypefn {Function File} {} qcoordinate_plot (@var{qf}, @var{qb}, @var{qv}) +## Plot in the current figure a set of coordinate axes as viewed from +## the orientation specified by quaternion @var{qv}. Inertial axes are +## also plotted: +## +## @table @var +## @item qf +## Quaternion from reference (x,y,z) to inertial. +## @item qb +## Quaternion from reference to body. +## @item qv +## Quaternion from reference to view angle. +## @end table +## @end deftypefn + +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Adapted-By: jwe + +function qcoordinate_plot (qf, qb, qv) -# compute inertial to body rotation eigenaxis -# qb = qbf*qf => qbf = qb/qf -# -# need to use inverse quaternion to rotate axes -qbf = qinv(qmult(qb,qinv(qf))); + degrees = pi / 180; + d180 = 180 * degrees; + + ## construct coordinate transformation to view frame + + cm = qtransvmat(qv); -[eaxv,th_eig] = quaternion(qbf); + p1 = [-1, -1, 1]; + p2 = [-1, -1, -1]; + p3 = [ 1, -1, -1]; + p4 = [ 1, -1, 1]; + p5 = [-1, 1, 1]; + p6 = [ 1, 1, 1]; + p7 = [ 1, 1, -1]; + p8 = [-1, 1, -1]; + + ## outline positive quadrant + + box1 = cm * [p4; p6; p5; p6; p7]'; + + ## outline rest of the box + + box2 = cm * [p7; p8; p5; p1; p4; p3; p7; p3; p2; p1; p2; p8]'; -# draw 1/3 circle in x-y plane around a unit z axis -th = (0:-12:-120)*degrees*sign(th_eig); lth = length(th); -cpts = [0, 0, 0.1*cos(th) ; 0, 0, 0.1*sin(th); 0, 1, 1*ones(1,lth)]; + ## compute inertial to body rotation eigenaxis + ## qb = qbf*qf => qbf = qb/qf + ## + ## need to use inverse quaternion to rotate axes + + qbf = qinv (qmult (qb, qinv (qf))); + + [eaxv, th_eig] = quaternion (qbf); + + ## draw 1/3 circle in x-y plane around a unit z axis + + th = (0:-12:-120) * degrees * sign (th_eig); + lth = length (th); -# rotate the 1/3 circle around eigenaxis of inertial to body rotation -# qez = qe/qz = rotation to get from z axis to eigenaxis. -# This rotates the 1/3 circle from x-y plane to the plane normal to -# eigenaxis -qez = qmult(qbf,qinv(quaternion(0,0,1,0))); -eig_xm = qtransvmat(qez); -cpts = cm*eig_xm*cpts; + cpts = [0, 0, 0.1*cos(th); + 0, 0, 0.1*sin(th); + 0, 1, 1*ones(1,lth)]; + + ## rotate the 1/3 circle around eigenaxis of inertial to body rotation + ## qez = qe/qz = rotation to get from z axis to eigenaxis. + ## This rotates the 1/3 circle from x-y plane to the plane normal to + ## eigenaxis -# transform inertial and body quaternions to view coordinates (rotate -# by azimuth, elevation) -qfm = qtransvmat(qf); qbm = qtransvmat(qf); -qf = qmult(qv,qf); qb = qmult(qv,qb); + qez = qmult (qbf, qinv (quaternion (0, 0, 1, 0))); + eig_xm = qtransvmat (qez); + cpts = cm*eig_xm * cpts; + + ## transform inertial and body quaternions to view coordinates (rotate + ## by azimuth, elevation) + + qfm = qtransvmat (qf); + qbm = qtransvmat (qf); -# get coordinate axes in inertial and reference frame -jnk = qtransvmat(qf); ifv = jnk(:,1); jfv = jnk(:,2); kfv = jnk(:,3); -jnk = qtransvmat(qb); ibv = jnk(:,1); jbv = jnk(:,2); kbv = jnk(:,3); + qf = qmult (qv, qf); + qb = qmult (qv, qb); + + ## get coordinate axes in inertial and reference frame + + jnk = qtransvmat (qf); + ifv = jnk(:,1); + jfv = jnk(:,2); + kfv = jnk(:,3); + + jnk = qtransvmat (qb); + ibv = jnk(:,1); + jbv = jnk(:,2); + kbv = jnk(:,3); -gset size square -axis([-2,2,-2,2]); -[vv,theta] = quaternion(qb); -xlabel(sprintf("rotate about eigenaxis %5.2f deg",th_eig/degrees)); -plot( [ibv(1),0],[ibv(3),0],"-@11;x (body);", ... - [0,jbv(1)],[0,jbv(3)],"-@21;y (body);", ... - [0,kbv(1)],[0,kbv(3)],"-@32;z (body);", ... - [ifv(1),0],[ifv(3),0],"-@13;x (inertial);", ... - [0,jfv(1)],[0,jfv(3)],"-@23;y (inertial);", ... - [0,kfv(1)],[0,kfv(3)],"-@34;z (inertial);", ... - cpts(1,:), cpts(3,:),".-6 ;eigenaxis;", ... - box2(1,:),box2(3,:),"-4;;", ... - box1(1,:),box1(3,:),"-5;;"); + gset size square + axis ([-2, 2, -2, 2]); + + [vv, theta] = quaternion (qb); + + xlabel (sprintf ("rotate about eigenaxis %5.2f deg", th_eig/degrees)); + + plot ([ibv(1), 0], [ibv(3), 0], "-@11;x (body);", + [0, jbv(1)], [0, jbv(3)], "-@21;y (body);", + [0, kbv(1)], [0, kbv(3)], "-@32;z (body);", + [ifv(1), 0], [ifv(3), 0], "-@13;x (inertial);", + [0, jfv(1)], [0, jfv(3)], "-@23;y (inertial);", + [0, kfv(1)], [0, kfv(3)], "-@34;z (inertial);", + cpts(1,:), cpts(3,:), ".-6 ;eigenaxis;", + box2(1,:), box2(3,:), "-4;;", + box1(1,:), box1(3,:), "-5;;"); + endfunction
--- a/scripts/quaternion/qderiv.m Tue Jan 18 04:49:49 2000 +0000 +++ b/scripts/quaternion/qderiv.m Tue Jan 18 06:01:04 2000 +0000 @@ -1,23 +1,56 @@ -function Dmat = qderivmat(Omega) -# function Dmat = qderivmat(Omega) -# Derivative of a quaternion. -# Let Q be a quaternion to transform a vector from a fixed frame to -# a rotating frame. If the rotating frame is rotating about the -# [x,y,z] axes at angular rates [wx, wy, wz], then the derivative -# of Q is given by -# Q' = qderivmat(Omega)*Q -# -# If the passive convention is used (rotate the frame, not the vector), -# then Q' = -qderivmat(Omega)*Q; see the accompanying document -# quaternion.ps for details. +## Copyright (C) 1998 Auburn University. All rights reserved. +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, or (at your option) +## any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, 59 Temple Place - Suite 330, Boston, MA +## 02111-1307, USA. -Omega = vec(Omega); -if(length(Omega) != 3) - error("qderivmat: Omega must be a length 3 vector."); -endif +## -*- texinfo -*- +## @deftypefn {Function File} {} qderiv (omega) +## Derivative of a quaternion. +## +## Let Q be a quaternion to transform a vector from a fixed frame to +## a rotating frame. If the rotating frame is rotating about the +## [x, y, z] axes at angular rates [wx, wy, wz], then the derivative +## of Q is given by +## +## @example +## Q' = qderivmat (omega) * Q +## @end example +## +## If the passive convention is used (rotate the frame, not the vector), +## then +## +## @example +## Q' = -qderivmat (omega) * Q +## @end example +## @end deftypefn -Dmat = 0.5*[ 0.0, Omega(3), -Omega(2), Omega(1); ... - -Omega(3), 0.0, Omega(1), Omega(2); ... - Omega(2), -Omega(1), 0.0, Omega(3); ... - -Omega(1), -Omega(2), -Omega(3), 0.0 ]; +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Adapted-By: jwe + +function Dmat = qderivmat (Omega) + + Omega = vec (Omega); + + if (length (Omega) != 3) + error ("qderivmat: Omega must be a length 3 vector."); + endif + + Dmat = 0.5 * [ 0.0, Omega(3), -Omega(2), Omega(1); + -Omega(3), 0.0, Omega(1), Omega(2); + Omega(2), -Omega(1), 0.0, Omega(3); + -Omega(1), -Omega(2), -Omega(3), 0.0 ]; endfunction
--- a/scripts/quaternion/qderivmat.m Tue Jan 18 04:49:49 2000 +0000 +++ b/scripts/quaternion/qderivmat.m Tue Jan 18 06:01:04 2000 +0000 @@ -1,23 +1,56 @@ -function Dmat = qderivmat(Omega) -# function Dmat = qderivmat(Omega) -# Derivative of a quaternion. -# Let Q be a quaternion to transform a vector from a fixed frame to -# a rotating frame. If the rotating frame is rotating about the -# [x,y,z] axes at angular rates [wx, wy, wz], then the derivative -# of Q is given by -# Q' = qderivmat(Omega)*Q -# -# If the passive convention is used (rotate the frame, not the vector), -# then Q' = -qderivmat(Omega)*Q; see the accompanying document -# quaternion.ps for details. +## Copyright (C) 1998 Auburn University. All rights reserved. +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, or (at your option) +## any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, 59 Temple Place - Suite 330, Boston, MA +## 02111-1307, USA. -Omega = vec(Omega); -if(length(Omega) != 3) - error("qderivmat: Omega must be a length 3 vector."); -endif +## -*- texinfo -*- +## @deftypefn {Function File} {} qderivmat (@var{omega}) +## Derivative of a quaternion. +## +## Let Q be a quaternion to transform a vector from a fixed frame to +## a rotating frame. If the rotating frame is rotating about the +## [x, y, z] axes at angular rates [wx, wy, wz], then the derivative +## of Q is given by +## +## @example +## Q' = qderivmat (omega) * Q +## @end example +## +## If the passive convention is used (rotate the frame, not the vector), +## then +## +## @example +## Q' = -qderivmat (omega) * Q. +## @end example +## @end deftypefn -Dmat = 0.5*[ 0.0, Omega(3), -Omega(2), Omega(1); ... - -Omega(3), 0.0, Omega(1), Omega(2); ... - Omega(2), -Omega(1), 0.0, Omega(3); ... - -Omega(1), -Omega(2), -Omega(3), 0.0 ]; +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Adapted-By: jwe + +function Dmat = qderivmat (Omega) + + Omega = vec (Omega); + + if (length (Omega) != 3) + error ("qderivmat: Omega must be a length 3 vector."); + endif + + Dmat = 0.5 * [ 0.0, Omega(3), -Omega(2), Omega(1); + -Omega(3), 0.0, Omega(1), Omega(2); + Omega(2), -Omega(1), 0.0, Omega(3); + -Omega(1), -Omega(2), -Omega(3), 0.0 ]; endfunction
--- a/scripts/quaternion/qinv.m Tue Jan 18 04:49:49 2000 +0000 +++ b/scripts/quaternion/qinv.m Tue Jan 18 06:01:04 2000 +0000 @@ -1,12 +1,41 @@ -function retval = qinv(q) -# function b = qinv(q) -# return the inverse of a quaternion -# q = [w,x,y,z] = w*i + x*j + y*k + z -# qmult(q,qinv(q)) = 1 = [0 0 0 1] +## Copyright (C) 1998 Auburn University. All rights reserved. +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, or (at your option) +## any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, 59 Temple Place - Suite 330, Boston, MA +## 02111-1307, USA. -if(norm(q) != 0) - retval = qconj(q) /sum(q .* q); -else - error("qinv: zero quaternion passed!"); -end +## -*- texinfo -*- +## @deftypefn {Function File} {} qinv (@var{q}) +## Return the inverse of a quaternion. +## +## @example +## q = [w, x, y, z] = w*i + x*j + y*k + z +## qmult (q, qinv (q)) = 1 = [0 0 0 1] +## @end example +## @end deftypefn + +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Adapted-By: jwe + +function retval = qinv (q) + + if (norm (q) != 0) + retval = qconj (q) / sum (q .* q); + else + error ("qinv: zero quaternion passed!"); + endif + endfunction
--- a/scripts/quaternion/qmult.m Tue Jan 18 04:49:49 2000 +0000 +++ b/scripts/quaternion/qmult.m Tue Jan 18 06:01:04 2000 +0000 @@ -1,21 +1,55 @@ -function retval = qmult(a,b) - # function retval = qmult(a,b) - # multiply two quaternions - # [w,x,y,z] = w*i + x*j + y*k + z - # identities: - # i^2 = j^2 = k^2 = -1 - # ij = k jk = i - # ki = j kj = -i - # ji = -k ik = -j +## Copyright (C) 1998 Auburn University. All rights reserved. +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, or (at your option) +## any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, 59 Temple Place - Suite 330, Boston, MA +## 02111-1307, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {} qmult (@var{a}, @var{b}) +## Multiply two quaternions. +## +## @example +## [w, x, y, z] = w*i + x*j + y*k + z +## @end example +## +## @noindent +## identities: +## +## @example +## i^2 = j^2 = k^2 = -1 +## ij = k jk = i +## ki = j kj = -i +## ji = -k ik = -j +## @end example +## @end deftypefn + +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Adapted-By: jwe + +function retval = qmult (a, b) - [a1,b1,c1,d1] = quaternion(a); - [a2,b2,c2,d2] = quaternion(b); + [a1, b1, c1, d1] = quaternion (a); + [a2, b2, c2, d2] = quaternion (b); ri = b1*c2 - c1*b2 + d1*a2 + a1*d2; rj = c1*a2 - a1*c2 + d1*b2 + b1*d2; rk = a1*b2 - b1*a2 + d1*c2 + c1*d2; rr = -(a1*a2 + b1*b2 + c1*c2) + d1*d2; - retval = quaternion(ri,rj,rk,rr); + retval = quaternion (ri, rj, rk, rr); + endfunction
--- a/scripts/quaternion/qtrans.m Tue Jan 18 04:49:49 2000 +0000 +++ b/scripts/quaternion/qtrans.m Tue Jan 18 06:01:04 2000 +0000 @@ -1,13 +1,39 @@ -function v = qtrans(v,q) -# function v = qtrans(v,q) -# transform the unit quaternion v by the unit quaternion q; -# v = [w x y z], q = transformation quaternion -# returns v = q*v/q +## Copyright (C) 1998 Auburn University. All rights reserved. +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, or (at your option) +## any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, 59 Temple Place - Suite 330, Boston, MA +## 02111-1307, USA. -if(!is_vector(v) | length(v) != 4) - error(sprintf("qtrans: v(%d,%d) must be a quaternion",rows(v),columns(v))) -elseif(!is_vector(q) | length(q) != 4) - error(sprintf("qtrans: q(%d,%d) must be a quaternion",rows(q),columns(q))) -endif +## -*- texinfo -*- +## @deftypefn {Function File} {} qtrans (@var{v}, @var{q}) +## Transform the unit quaternion @var{v} by the unit quaternion @var{q}. +## Returns @code{@var{v} = @var{q}*@var{v}/@var{q}}. +## @end deftypefn + +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Adapted-By: jwe -v = qmult(q,qmult(v,qinv(q))); +function v = qtrans (v, q) + + if (! is_vector (v) || length (v) != 4) + error ("qtrans: v(%d,%d) must be a quaternion", rows (v), columns (v)); + elseif (! is_vector (q) || length (q) != 4) + error ("qtrans: q(%d,%d) must be a quaternion", rows (q), columns (q)); + endif + + v = qmult (q, qmult (v, qinv (q))); + +endfunction
--- a/scripts/quaternion/qtransv.m Tue Jan 18 04:49:49 2000 +0000 +++ b/scripts/quaternion/qtransv.m Tue Jan 18 06:01:04 2000 +0000 @@ -1,23 +1,52 @@ -function vi = qtransv(vb,qib) -# function vi = qtransv(vb,q) -# transform the 3-D vector v by the unit quaternion q; -# v = [w x y z], q = transformation quaternion -# returns vi = column vector -# vi = (2*real(q)^2 - 1)*vb + 2*imag(q)*(imag(q)'*vb) -# + 2*real(q)*cross(imag(q),vb) -# where imag(q) is a column vector of length 3. +## Copyright (C) 1998 Auburn University. All rights reserved. +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, or (at your option) +## any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, 59 Temple Place - Suite 330, Boston, MA +## 02111-1307, USA. -if(!is_vector(vb) | length(vb) != 3) - error(sprintf("qtransv: v(%d,%d) must be a 3-D vector",rows(vb),columns(vb))) -elseif(!is_vector(qib) | length(qib) != 4) - error(sprintf("qtransv: q(%d,%d) must be a quaternion",rows(qib),columns(qib))) -elseif(max(abs(imag(vb))) + max(abs(imag(qib))) != 0) - vb - qib - error("qtransv: input values must be real."); -endif +## -*- texinfo -*- +## @deftypefn {Function File} {} qtransv (@var{v}, @var{q}) +## Transform the 3-D vector @var{v} by the unit quaternion @var{q}. +## Return a column vector. +## +## @example +## vi = (2*real(q)^2 - 1)*vb + 2*imag(q)*(imag(q)'*vb) +## + 2*real(q)*cross(imag(q),vb) +## @end example +## +## @noindent +## Where imag(q) is a column vector of length 3. +## @end deftypefn -qr = qib(4); qimag = vec(qib(1:3)); vb = vec(vb); -vi = (2*qr^2 - 1)*vb + 2*qimag*(qimag'*vb) + 2*qr*cross(qimag,vb); +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Adapted-By: jwe + +function vi = qtransv (vb, qib) + + if (! is_vector (vb) || length (vb) != 3) + error ("qtransv: v(%d,%d) must be a 3-D vector", rows (vb), columns (vb)); + elseif (! is_vector (qib) || length (qib) != 4) + error ("qtransv: q(%d,%d) must be a quaternion", rows (qib), columns (qib)); + elseif (max (abs (imag (vb))) + max (abs (imag (qib))) != 0) + error ("qtransv: input values must be real."); + endif + + qr = qib(4); + qimag = vec (qib(1:3)); + vb = vec (vb); + vi = (2*qr^2 - 1)*vb + 2*qimag*(qimag'*vb) + 2*qr*cross (qimag, vb); endfunction
--- a/scripts/quaternion/qtransvmat.m Tue Jan 18 04:49:49 2000 +0000 +++ b/scripts/quaternion/qtransvmat.m Tue Jan 18 06:01:04 2000 +0000 @@ -1,24 +1,49 @@ -function Aib = qtransvmat(qib) -# function Aib = qtransvmat(qib) -# construct 3x3 transformation matrix from quaternion qib -# Aib is equivalent to rotation of th radians about axis vv, where -# [vv,th] = quaternion(qib) +## Copyright (C) 1998 Auburn University. All rights reserved. +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, or (at your option) +## any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, 59 Temple Place - Suite 330, Boston, MA +## 02111-1307, USA. -if(!is_vector(qib) | length(qib) != 4) - error(sprintf("qtransvmat: q(%d,%d) must be a quaternion",rows(qib),columns(qib))) -elseif(max(abs(imag(qib))) != 0) - qib - error("qtransvmat: input values must be real."); -endif +## -*- texinfo -*- +## @deftypefn {Function File} {} qtransvmat (@var{qib}) +## Construct a 3x3 transformation matrix from quaternion @var{qib} that +## is equivalent to rotation of th radians about axis @var{vv}, where +## @code{[@var{vv}, @var{th}] = quaternion (@var{qib})}. +## @end deftypefn + +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Adapted-By: jwe + +function Aib = qtransvmat (qib) -Aib = [ (2.*(qib(1)^2 + qib(4)^2) -1.), ... - (2.*(qib(1)*qib(2)-qib(3)*qib(4))), ... - (2.*(qib(1)*qib(3)+qib(2)*qib(4))); ... - (2.*(qib(1)*qib(2)+qib(3)*qib(4))), ... - (2.*(qib(2)*qib(2)+qib(4)*qib(4))-1.), ... - (2.*(qib(2)*qib(3)-qib(1)*qib(4))); ... - (2.*(qib(1)*qib(3)-qib(2)*qib(4))), ... - (2.*(qib(2)*qib(3)+qib(1)*qib(4))), ... - (2.*(qib(3)*qib(3)+qib(4)*qib(4))-1.)]; + if (! is_vector(qib) || length (qib) != 4) + error ("qtransvmat: q(%d,%d) must be a quaternion", rows (qib), \ + columns (qib)); + elseif (max (abs (imag (qib))) != 0) + error("qtransvmat: input values must be real."); + endif + + Aib = [(2.*(qib(1)^2 + qib(4)^2) -1.), + (2.*(qib(1)*qib(2)-qib(3)*qib(4))), + (2.*(qib(1)*qib(3)+qib(2)*qib(4))); + (2.*(qib(1)*qib(2)+qib(3)*qib(4))), + (2.*(qib(2)*qib(2)+qib(4)*qib(4))-1.), + (2.*(qib(2)*qib(3)-qib(1)*qib(4))); + (2.*(qib(1)*qib(3)-qib(2)*qib(4))), + (2.*(qib(2)*qib(3)+qib(1)*qib(4))), + (2.*(qib(3)*qib(3)+qib(4)*qib(4))-1.)]; + endfunction -
--- a/scripts/quaternion/quaternion.m Tue Jan 18 04:49:49 2000 +0000 +++ b/scripts/quaternion/quaternion.m Tue Jan 18 06:01:04 2000 +0000 @@ -1,86 +1,116 @@ -function [a,b,c,d] = quaternion(w,x,y,z) -# quaternion: construct or extract a quaternion -# w = a*i + b*j + c*k + d from given data. -# -# calling formats: -# [a,b,c,d] = quaternion(w) -or- -# [vv,theta] = quaternion(w) -# w = quaternion(a,b,c,d) -# w = quaternion(vv,theta) +## Copyright (C) 1998 Auburn University. All rights reserved. +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2, or (at your option) +## any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +## General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, write to the Free +## Software Foundation, 59 Temple Place - Suite 330, Boston, MA +## 02111-1307, USA. + +## -*- texinfo -*- +## @deftypefn {Function File} {[@var{a}, @var{b}, @var{c}, @var{d}] =} quaternion (w) +## @deftypefnx {Function File} {[@var{vv}, @var{theta}] =} quaternion (w) +## @deftypefnx {Function File} {@var{w} =} quaternion (@var{a}, @var{b}, @var{c}, @var{d}) +## @deftypefnx {Function File} {@var{w} =} quaternion (@var{vv}, @var{theta}) +## Construct or extract a quaternion +## +## @example +## w = a*i + b*j + c*k + d +## @end example +## +## @noindent +## from given data. +## @end deftypefn + +## Author: A. S. Hodel <a.s.hodel@eng.auburn.edu> +## Adapted-By: jwe + +function [a, b, c, d] = quaternion (w, x, y, z) + + switch (nargin) + case(1) + if (! (is_vector (w) && length (w) == 4)) + error ("input vector must be of length 4)"); + endif + ## extract data + switch (nargout) + case(4) + a = w(1); + b = w(2); + c = w(3); + d = w(4); -switch(nargin) -case(1), # extract information - if(!(is_vector(w) & length(w) == 4) ) - error("input vector must be of length 4)") - endif - # extract data - switch(nargout) - case(4), - a = w(1); - b = w(2); - c = w(3); - d = w(4); - case(2), - if(abs(norm(w)-1) > 1e-12) - warning(sprintf("quaternion: ||w||=%e, setting=1 for vv, theta",norm(w))) - w = w/norm(w); + case(2) + if (abs (norm (w) - 1) > 1e-12) + warning ("quaternion: ||w||=%e, setting=1 for vv, theta", norm(w)); + w = w/norm(w); + endif + [a, b, c, d] = quaternion (w); + theta = acos (d) * 2; + if (abs (theta) > pi) + theta = theta - sign (theta) * pi; + endif + sin_th_2 = norm ([a, b, c]); + + if (sin_th_2 != 0) + vv = [a, b, c] / sin_th_2; + else + vv = [a, b, c]; + endif + a = vv; + b = theta; + otherwise + usage ("[a, b, c, d] = quaternion (w) or [vv, theta] = quaternion (w)"); + endswitch + + case(2) + if (nargout != 1) + usage ("w = quaterion (vv, theta)"); endif - [a,b,c,d] = quaternion(w); - theta = acos(d)*2; - if(abs(theta) > pi) - theta = theta - sign(theta)*pi; + vv = w; + theta = x; + + if (! is_vector (vv) || length (vv) != 3) + error ("vv must be a length three vector"); + elseif (! is_scalar (theta)) + error ("theta must be a scalar"); + elseif (norm (vv) == 0) + error ("quaternion: vv is zero."); + elseif (abs (norm (vv) - 1) > 1e-12) + warning ("quaternion: ||vv|| != 1, normalizing") + vv = vv / norm (vv); endif - sin_th_2 = norm([a, b, c]); - if(sin_th_2 != 0) - vv = [a,b,c]/sin_th_2; - else - vv = [a, b, c]; + if (abs (theta) > 2*pi) + warning ("quaternion: |theta| > 2 pi, normalizing") + theta = rem (theta, 2*pi); endif - a = vv; - b = theta; - otherwise, - usage("[a,b,c,d] = quaternion(w) or [vv,theta] = quaternion(w)"); + vv = vv * sin (theta / 2); + d = cos (theta / 2); + a = quaternion (vv(1), vv(2), vv(3), d); + + case(4) + if (nargout != 1) + usage ("w = quaterion (a, b, c, d)"); + endif + if (! (is_scalar (w) && is_scalar (x) && is_scalar (y) && is_scalar (z))) + error ("input values must be scalars."); + endif + a = [w, x, y, z]; + + otherwise + error ("usage: [a, b, c, d] = quaternion (w), a = quaternion (w, x, y, z)"); + endswitch -case(2), - if(nargout != 1) - usage("w = quaterion(vv,theta)"); - endif - vv = w; - theta = x; - - if(!is_vector(vv) | length(vv) != 3) - error("vv must be a length three vector") - elseif(!is_scalar(theta)) - error("theta must be a scalar"); - elseif(norm(vv) == 0) - error("quaternion: vv is zero.") - elseif(abs(norm(vv)-1) > 1e-12) - warning("quaternion: ||vv|| != 1, normalizing") - vv = vv/norm(vv); - endif - - if(abs(theta) > 2*pi) - warning("quaternion: |theta| > 2 pi, normalizing") - theta = rem(theta,2*pi); - endif - vv = vv*sin(theta/2); - d = cos(theta/2); - a = quaternion(vv(1), vv(2), vv(3), d); - -case(4), - if(nargout != 1) - usage("w = quaterion(a,b,c,d)"); - endif - if ( !(is_scalar(w) & is_scalar(x) & is_scalar(y) & is_scalar(z)) ) - error("input values must be scalars.") - endif - a = [w, x, y, z]; - -otherwise, - error("usage: [a,b,c,d] = quaternion(w), a = quaternion(w,x,y,z)") - -endswitch - endfunction