v0.3.10 changes.

nicer if the Occluder would always get removed before its model branch,
but that's not easily enforcable)

Doxyfile

@@ -22,7 +22,7 @@ PROJECT_NAME           = SimGear
 # This could be handy for archiving the generated documentation or 
 # if some version control system is used.
 
-PROJECT_NUMBER         = 0.3.9
+PROJECT_NUMBER         = 0.3.10
 
 # The OUTPUT_DIRECTORY tag is used to specify the (relative or absolute) 
 # base path where the generated documentation will be put. 
@@ -303,16 +303,19 @@ INPUT                  = \
 	simgear/compiler.h \
 	simgear/constants.h \
 	simgear/debug \
+	simgear/environment \
 	simgear/ephemeris \
 	simgear/io \
 	simgear/magvar \
 	simgear/math \
 	simgear/misc \
+	simgear/nasal \
         simgear/props \
 	simgear/route \
 	simgear/scene \
 	simgear/screen \
 	simgear/serial \
+	simgear/structure \
 	simgear/sg_inlines.h \
 	simgear/sg_traits.hxx \
         simgear/sound \

Makefile.am

@@ -4,13 +4,15 @@ EXTRA_DIST = \
 	DoxygenMain.cxx \
 	README.MSVC \
 	README.zlib \
+        projects \
 	SimGear.dsp \
 	SimGear.dsw
 
-SUBDIRS	= src-libs simgear
+SUBDIRS	= simgear
 
 dist-hook:
 	(cd $(top_srcdir); $(HOME)/Projects/FlightGear/admin/am2dsp.pl)
+	rm -rf `find $(distdir)/projects -name CVS`
 
 #
 # Rule to build RPM distribution package

NEWS

@@ -1,3 +1,31 @@
+New in 0.3.10
+* April 5, 2006
+
+* Add a small accessor function to expose local timezone offset.
+* Improved exception handling and made output more helpful in various places.
+* Better pbuffer runtime detection.
+* Add directory creation capability to file/path library.
+* Added a basic reference counting class to improve robustness of
+  memory management in places.  Use this for all scenegraph
+  references, sgmaterial references, sgmatmodel references, and
+  sgsoundsample references.
+* Add support for point sprites.
+* Updates to rain cone rendering.
+* Add a new vector library and integrate that with improved coordinate
+  system conversion code.
+* Mutex locking and cleanup improvements in the threading abstraction
+  library.
+* Add MacOS RenderTexture support.
+* Add a Nasal based io libarary that is not activated by default.
+* Added a set of MS-VC8 project files.
+
+* Various platform related bug fixes.
+* Various compiler related bug/warning fixes.
+* Clean up some things that triggered valgrind warnings.
+* Fix a Nasal cmp() bug.
+* Removed ancient version of zlib from distribution.
+
+
 New in 0.3.9
 * November 17, 2005
 

README.zlib

@@ -1,173 +1,8 @@
-For your convenience and allowed by zlib's license terms:
-http://www.gzip.org/zlib/zlib_license.html a copy of the zlib source
-is bundled with SimGear in $(top_srcdir)/src-libs/.  You must have
-zlib installed before you can build SimGear.
+zlib is prerequisite for compiling SimGear and FlightGear.  You can fetch
+the latest version of zlib from:
 
-- Most linux distributions have a zlib package.  For linux
-  developers, we recommend ysou install your distributions package
-  rather than building from source.
+    http://www.zlib.net
 
-- Cygwin installs zlib automatically.
+Most modern unix distributions (and cygwin) come with a version of zlib
+already installed or available to install as a package.
 
-- For developers on most other platforms, you will have to build
-  zlib from source and install it yourself.  For your convenience a
-  tar ball of the zlib source is included with the simgear source
-  distribution.  Untar the zlib source, and follow the included
-  build and installation instructions.
-
-Once zlib is installed you can return to configuring and building
-Simgear.
-
-We now send you to the official zlib README ... 
-
-
-=============================================================================
-
-
-zlib 1.1.4 is a general purpose data compression library.  All the code
-is thread safe.  The data format used by the zlib library
-is described by RFCs (Request for Comments) 1950 to 1952 in the files 
-http://www.ietf.org/rfc/rfc1950.txt (zlib format), rfc1951.txt (deflate
-format) and rfc1952.txt (gzip format). These documents are also available in
-other formats from ftp://ftp.uu.net/graphics/png/documents/zlib/zdoc-index.html
-
-All functions of the compression library are documented in the file zlib.h
-(volunteer to write man pages welcome, contact jloup@gzip.org). A usage
-example of the library is given in the file example.c which also tests that
-the library is working correctly. Another example is given in the file
-minigzip.c. The compression library itself is composed of all source files
-except example.c and minigzip.c.
-
-To compile all files and run the test program, follow the instructions
-given at the top of Makefile. In short "make test; make install"
-should work for most machines. For Unix: "./configure; make test; make install"
-For MSDOS, use one of the special makefiles such as Makefile.msc.
-For VMS, use Make_vms.com or descrip.mms.
-
-Questions about zlib should be sent to <zlib@gzip.org>, or to
-Gilles Vollant <info@winimage.com> for the Windows DLL version.
-The zlib home page is http://www.zlib.org or http://www.gzip.org/zlib/
-Before reporting a problem, please check this site to verify that
-you have the latest version of zlib; otherwise get the latest version and
-check whether the problem still exists or not.
-
-PLEASE read the zlib FAQ http://www.gzip.org/zlib/zlib_faq.html
-before asking for help.
-
-Mark Nelson <markn@ieee.org> wrote an article about zlib for the Jan. 1997
-issue of  Dr. Dobb's Journal; a copy of the article is available in
-http://dogma.net/markn/articles/zlibtool/zlibtool.htm
-
-The changes made in version 1.1.4 are documented in the file ChangeLog.
-The only changes made since 1.1.3 are bug corrections:
-
-- ZFREE was repeated on same allocation on some error conditions.
-  This creates a security problem described in
-  http://www.zlib.org/advisory-2002-03-11.txt
-- Returned incorrect error (Z_MEM_ERROR) on some invalid data
-- Avoid accesses before window for invalid distances with inflate window
-  less than 32K.
-- force windowBits > 8 to avoid a bug in the encoder for a window size
-  of 256 bytes. (A complete fix will be available in 1.1.5).
-
-The beta version 1.1.5beta includes many more changes. A new official
-version 1.1.5 will be released as soon as extensive testing has been
-completed on it.
-
-
-Unsupported third party contributions are provided in directory "contrib".
-
-A Java implementation of zlib is available in the Java Development Kit
-http://www.javasoft.com/products/JDK/1.1/docs/api/Package-java.util.zip.html
-See the zlib home page http://www.zlib.org for details.
-
-A Perl interface to zlib written by Paul Marquess <pmarquess@bfsec.bt.co.uk>
-is in the CPAN (Comprehensive Perl Archive Network) sites
-http://www.cpan.org/modules/by-module/Compress/
-
-A Python interface to zlib written by A.M. Kuchling <amk@magnet.com>
-is available in Python 1.5 and later versions, see
-http://www.python.org/doc/lib/module-zlib.html
-
-A zlib binding for TCL written by Andreas Kupries <a.kupries@westend.com>
-is availlable at http://www.westend.com/~kupries/doc/trf/man/man.html
-
-An experimental package to read and write files in .zip format,
-written on top of zlib by Gilles Vollant <info@winimage.com>, is
-available at http://www.winimage.com/zLibDll/unzip.html
-and also in the contrib/minizip directory of zlib.
-
-
-Notes for some targets:
-
-- To build a Windows DLL version, include in a DLL project zlib.def, zlib.rc
-  and all .c files except example.c and minigzip.c; compile with -DZLIB_DLL
-  The zlib DLL support was initially done by Alessandro Iacopetti and is
-  now maintained by Gilles Vollant <info@winimage.com>. Check the zlib DLL
-  home page at http://www.winimage.com/zLibDll
-
-  From Visual Basic, you can call the DLL functions which do not take
-  a structure as argument: compress, uncompress and all gz* functions.
-  See contrib/visual-basic.txt for more information, or get
-  http://www.tcfb.com/dowseware/cmp-z-it.zip
-
-- For 64-bit Irix, deflate.c must be compiled without any optimization.
-  With -O, one libpng test fails. The test works in 32 bit mode (with
-  the -n32 compiler flag). The compiler bug has been reported to SGI.
-
-- zlib doesn't work with gcc 2.6.3 on a DEC 3000/300LX under OSF/1 2.1   
-  it works when compiled with cc.
-
-- on Digital Unix 4.0D (formely OSF/1) on AlphaServer, the cc option -std1
-  is necessary to get gzprintf working correctly. This is done by configure.
-
-- zlib doesn't work on HP-UX 9.05 with some versions of /bin/cc. It works
-  with other compilers. Use "make test" to check your compiler.
-
-- gzdopen is not supported on RISCOS, BEOS and by some Mac compilers.
-
-- For Turbo C the small model is supported only with reduced performance to
-  avoid any far allocation; it was tested with -DMAX_WBITS=11 -DMAX_MEM_LEVEL=3
-
-- For PalmOs, see http://www.cs.uit.no/~perm/PASTA/pilot/software.html
-  Per Harald Myrvang <perm@stud.cs.uit.no>
-
-
-Acknowledgments:
-
-  The deflate format used by zlib was defined by Phil Katz. The deflate
-  and zlib specifications were written by L. Peter Deutsch. Thanks to all the
-  people who reported problems and suggested various improvements in zlib;
-  they are too numerous to cite here.
-
-Copyright notice:
-
- (C) 1995-2002 Jean-loup Gailly and Mark Adler
-
-  This software is provided 'as-is', without any express or implied
-  warranty.  In no event will the authors be held liable for any damages
-  arising from the use of this software.
-
-  Permission is granted to anyone to use this software for any purpose,
-  including commercial applications, and to alter it and redistribute it
-  freely, subject to the following restrictions:
-
-  1. The origin of this software must not be misrepresented; you must not
-     claim that you wrote the original software. If you use this software
-     in a product, an acknowledgment in the product documentation would be
-     appreciated but is not required.
-  2. Altered source versions must be plainly marked as such, and must not be
-     misrepresented as being the original software.
-  3. This notice may not be removed or altered from any source distribution.
-
-  Jean-loup Gailly        Mark Adler
-  jloup@gzip.org          madler@alumni.caltech.edu
-
-If you use the zlib library in a product, we would appreciate *not*
-receiving lengthy legal documents to sign. The sources are provided
-for free but without warranty of any kind.  The library has been
-entirely written by Jean-loup Gailly and Mark Adler; it does not
-include third-party code.
-
-If you redistribute modified sources, we would appreciate that you include
-in the file ChangeLog history information documenting your changes.

SimGear.dsp

@@ -35,7 +35,7 @@ RSC=rc.exe
 # PROP Output_Dir "Release"
 # PROP Intermediate_Dir "Release"
 # PROP Target_Dir ""
-# ADD CPP /nologo /W3 /GX /O2 /D "NDEBUG" /D "WIN32" /D "_MBCS" /FD /c  /MT /I "." /I ".." /I ".\SimGear" /I ".\SimGear\metakit-2.4.3\include" /I "..\SimGear\zlib-1.1.4" /D "HAVE_CONFIG_H"
+# ADD CPP /nologo /W3 /GX /O2 /D "NDEBUG" /D "WIN32" /D "_MBCS" /FD /c  /MT /I "." /I ".." /I ".\SimGear" /I "..\zlib-1.2.3" /I "..\OpenAL 1.0 Software Development Kit\include" /D "_USE_MATH_DEFINES" /D "_CRT_SECURE_NO_DEPRECATE" /D "HAVE_CONFIG_H"
 # ADD RSC /l 0x409 /d "NDEBUG"
 BSC32=bscmake.exe
 # ADD BASE BSC32 /nologo
@@ -51,7 +51,7 @@ LINK32=link.exe -lib
 # PROP Output_Dir "Debug"
 # PROP Intermediate_Dir "Debug"
 # PROP Target_Dir ""
-# ADD CPP /nologo /W3 /GX /ZI /Od /D "_DEBUG" /D "WIN32" /D "_MBCS" /FR /FD /GZ /c  /MT /I "." /I ".." /I ".\SimGear" /I ".\SimGear\metakit-2.4.3\include" /I "..\SimGear\zlib-1.1.4" /D "HAVE_CONFIG_H"
+# ADD CPP /nologo /W3 /GX /ZI /Od /D "_DEBUG" /D "WIN32" /D "_MBCS" /FR /FD /GZ /c  /MTd /I "." /I ".." /I ".\SimGear" /I "..\zlib-1.2.3" /I "..\OpenAL 1.0 Software Development Kit\include" /D "_USE_MATH_DEFINES" /D "_CRT_SECURE_NO_DEPRECATE" /D "HAVE_CONFIG_H"
 # ADD RSC /l 0x409 /d "_DEBUG"
 BSC32=bscmake.exe
 # ADD BASE BSC32 /nologo
@@ -568,6 +568,21 @@ SOURCE=.\simgear\math\fastmath.cxx
 
 !ENDIF 
 
+# End Source File
+# Begin Source File
+
+SOURCE=.\simgear\math\SGGeodesy.cxx
+
+!IF  "$(CFG)" == "SimGear - Win32 Release"
+
+# PROP Intermediate_Dir "Release\Lib_sgmath"
+
+!ELSEIF  "$(CFG)" == "SimGear - Win32 Debug"
+
+# PROP Intermediate_Dir "Debug\Lib_sgmath"
+
+!ENDIF 
+
 # End Source File
 # End Group
 # Begin Group "Lib_sgenvironment"
@@ -853,6 +868,51 @@ SOURCE=.\simgear\nasal\mathlib.c
 # End Source File
 # Begin Source File
 
+SOURCE=.\simgear\nasal\iolib.c
+
+!IF  "$(CFG)" == "SimGear - Win32 Release"
+
+# PROP Intermediate_Dir "Release\Lib_sgnasal"
+
+!ELSEIF  "$(CFG)" == "SimGear - Win32 Debug"
+
+# PROP Intermediate_Dir "Debug\Lib_sgnasal"
+
+!ENDIF 
+
+# End Source File
+# Begin Source File
+
+SOURCE=.\simgear\nasal\iolib.h
+
+!IF  "$(CFG)" == "SimGear - Win32 Release"
+
+# PROP Intermediate_Dir "Release\Lib_sgnasal"
+
+!ELSEIF  "$(CFG)" == "SimGear - Win32 Debug"
+
+# PROP Intermediate_Dir "Debug\Lib_sgnasal"
+
+!ENDIF 
+
+# End Source File
+# Begin Source File
+
+SOURCE=.\simgear\nasal\bitslib.c
+
+!IF  "$(CFG)" == "SimGear - Win32 Release"
+
+# PROP Intermediate_Dir "Release\Lib_sgnasal"
+
+!ELSEIF  "$(CFG)" == "SimGear - Win32 Debug"
+
+# PROP Intermediate_Dir "Debug\Lib_sgnasal"
+
+!ENDIF 
+
+# End Source File
+# Begin Source File
+
 SOURCE=.\simgear\nasal\misc.c
 
 !IF  "$(CFG)" == "SimGear - Win32 Release"

SimGear.dsw

@@ -15,30 +15,6 @@ Package=<4>
 
 ###############################################################################
 
-Project: "mklib"=".\SimGear\metakit-2.4.3\win\msvc60\mklib.dsp" - Package Owner=<4>
-
-Package=<5>
-{{{
-}}}
-
-Package=<4>
-{{{
-}}}
-
-###############################################################################
-
-Project: "zlib"=".\SimGear\zlib.dsp" - Package Owner=<4>
-
-Package=<5>
-{{{
-}}}
-
-Package=<4>
-{{{
-}}}
-
-###############################################################################
-
 Global:
 
 Package=<5>

am2dsp.cfg

@@ -4,14 +4,13 @@ exclude_dir = threads
 include_path = .
 include_path = ..
 include_path = .\SimGear
-include_path = .\SimGear\metakit-2.4.3\include
-include_path = ..\SimGear\zlib-1.1.4
+include_path = ..\zlib-1.2.3
+include_path = "..\OpenAL 1.0 Software Development Kit\include"
 
+define = _USE_MATH_DEFINES
+define = _CRT_SECURE_NO_DEPRECATE
 define = HAVE_CONFIG_H
 
-add_project = .\SimGear\metakit-2.4.3\win\msvc60\mklib.dsp,mklib
-add_project = .\SimGear\zlib.dsp,zlib
-
 # Rule to create simgear_config.h
 add_source_file = SOURCE=.\simgear\simgear_config.h.vc5\
 \

configure.ac

@@ -8,7 +8,7 @@ dnl Require at least automake 2.52
 AC_PREREQ(2.52)
 
 dnl Initialize the automake stuff
-AM_INIT_AUTOMAKE(SimGear, 0.3.9)
+AM_INIT_AUTOMAKE(SimGear, 0.3.10)
 
 dnl Specify KAI C++ compiler and flags.
 dnl Borrowed with slight modification from blitz distribution.
@@ -250,6 +250,9 @@ case "${host}" in
 
 esac
 
+AC_SEARCH_LIBS(glutGetModifiers, [ glut glut32 freeglut ], have_glut=yes, have_glut=no)
+AM_CONDITIONAL(HAVE_GLUT, test "x$have_glut" = "xyes")
+
 opengl_LIBS="$LIBS"
 LIBS="$base_LIBS"
 
@@ -387,6 +390,7 @@ fi
 dnl Checks for typedefs, structures, and compiler characteristics.
 AC_C_CONST
 AC_TYPE_SIZE_T
+AC_TYPE_MODE_T
 AC_HEADER_TIME
 AC_STRUCT_TM
 
@@ -405,7 +409,6 @@ AM_CONFIG_HEADER(simgear/simgear_config.h)
 AC_CONFIG_FILES([ \
 	Makefile \
 	SimGear.spec \
-	src-libs/Makefile \
 	simgear/Makefile \
 	simgear/version.h \
 	simgear/compatibility/Makefile \

projects/VC8/SimGear.sln

@@ -0,0 +1,20 @@
+
+Microsoft Visual Studio Solution File, Format Version 9.00
+# Visual C++ Express 2005
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "SimGear", "SimGear.vcproj", "{952B5B5B-7FC8-4AE9-A664-333322BEEEB0}"
+EndProject
+Global
+	GlobalSection(SolutionConfigurationPlatforms) = preSolution
+		Debug|Win32 = Debug|Win32
+		Release|Win32 = Release|Win32
+	EndGlobalSection
+	GlobalSection(ProjectConfigurationPlatforms) = postSolution
+		{952B5B5B-7FC8-4AE9-A664-333322BEEEB0}.Debug|Win32.ActiveCfg = Debug|Win32
+		{952B5B5B-7FC8-4AE9-A664-333322BEEEB0}.Debug|Win32.Build.0 = Debug|Win32
+		{952B5B5B-7FC8-4AE9-A664-333322BEEEB0}.Release|Win32.ActiveCfg = Release|Win32
+		{952B5B5B-7FC8-4AE9-A664-333322BEEEB0}.Release|Win32.Build.0 = Release|Win32
+	EndGlobalSection
+	GlobalSection(SolutionProperties) = preSolution
+		HideSolutionNode = FALSE
+	EndGlobalSection
+EndGlobal

simgear/bucket/newbucket.cxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/bucket/newbucket.hxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/
@@ -115,10 +114,10 @@ class SGBucket {
 
 private:
     double cx, cy;  // centerpoint (lon, lat) in degrees of bucket
-    int lon;        // longitude index (-180 to 179)
-    int lat;        // latitude index (-90 to 89)
-    int x;          // x subdivision (0 to 7)
-    int y;          // y subdivision (0 to 7)
+    short lon;        // longitude index (-180 to 179)
+    short lat;        // latitude index (-90 to 89)
+    char x;          // x subdivision (0 to 7)
+    char y;          // y subdivision (0 to 7)
 
 public:
 

simgear/compiler.h

@@ -11,10 +11,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  *

simgear/constants.h

@@ -14,10 +14,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 
@@ -32,9 +31,6 @@
 
 #include <simgear/compiler.h>
 
-#ifdef _MSC_VER
-#  define  _USE_MATH_DEFINES
-#endif
 #ifdef SG_HAVE_STD_INCLUDES
 #  include <cmath>
 #else

simgear/debug/logstream.cxx

@@ -14,10 +14,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/debug/logstream.hxx

@@ -16,10 +16,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 
@@ -194,7 +193,7 @@ logbuf::overflow( int c )
             freopen("conout$", "w", stderr);
             has_console = true;
         }
-        sbuf->sputc(c);
+        return sbuf->sputc(c);
     }
     else
         return EOF == 0 ? 1: 0;

simgear/environment/metar.cxx

@@ -16,7 +16,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 
@@ -107,7 +107,7 @@ SGMetar::SGMetar(const string& m, const string& proxy, const string& port,
 	scanType();
 	if (!scanId() || !scanDate()) {
 		delete[] _data;
-		throw sg_io_exception("metar data bogus (" + _url + ')');
+		throw sg_io_exception("metar data bogus ", sg_location(_url));
 	}
 	scanModifier();
 
@@ -133,7 +133,7 @@ SGMetar::SGMetar(const string& m, const string& proxy, const string& port,
 
 	if (_grpcount < 4) {
 		delete[] _data;
-		throw sg_io_exception("metar data incomplete (" + _url + ')');
+		throw sg_io_exception("metar data incomplete ", sg_location(_url));
 	}
 
 	_url = "";
@@ -196,7 +196,7 @@ char *SGMetar::loadData(const char *id, const string& proxy, const string& port,
 	sock->set_timeout(10000);
 	if (!sock->open(SG_IO_OUT)) {
 		delete sock;
-		throw sg_io_exception("cannot connect to " + host);
+		throw sg_io_exception("cannot connect to ", sg_location(host));
 	}
 
 	string get = "GET ";
@@ -218,7 +218,7 @@ char *SGMetar::loadData(const char *id, const string& proxy, const string& port,
 	while ((i = sock->readline(buf, buflen))) {
 		if (i <= 2 && isspace(buf[0]) && (!buf[1] || isspace(buf[1])))
 			break;
-		if (!strncmp(buf, "X-MetarProxy: ", 9))
+		if (!strncmp(buf, "X-MetarProxy: ", 13))
 			_x_proxy = true;
 	}
 	if (i) {
@@ -233,7 +233,8 @@ char *SGMetar::loadData(const char *id, const string& proxy, const string& port,
 	char *b = buf;
 	scanBoundary(&b);
 	if (*b == '<')
-		throw sg_io_exception("no metar data available from " + _url);
+		throw sg_io_exception("no metar data available from ", 
+				sg_location(_url));
 
 	char *metar = new char[strlen(b) + 2];	// make room for " \0"
 	strcpy(metar, b);
@@ -608,56 +609,56 @@ bool SGMetar::scanRwyVisRange()
 
 
 static const struct Token special[] = {
-	"NSW",	"no significant weather",
-	"VCSH",	"showers in the vicinity",
-	"VCTS",	"thunderstorm in the vicinity",
-	0, 0
+	{ "NSW",  "no significant weather" },
+	{ "VCSH", "showers in the vicinity" },
+	{ "VCTS", "thunderstorm in the vicinity" },
+	{ 0, 0 }
 };
 
 
 static const struct Token description[] = {
-	"SH",	"showers of",
-	"TS",	"thunderstorm with",
-	"BC",	"patches of",
-	"BL",	"blowing",
-	"DR",	"low drifting",
-	"FZ",	"freezing",
-	"MI",	"shallow",
-	"PR",	"partial",
-	0, 0
+	{ "SH", "showers of" },
+	{ "TS", "thunderstorm with" },
+	{ "BC", "patches of" },
+	{ "BL", "blowing" },
+	{ "DR", "low drifting" },
+	{ "FZ", "freezing" },
+	{ "MI", "shallow" },
+	{ "PR", "partial" },
+	{ 0, 0 }
 };
 
 
 static const struct Token phenomenon[] = {
-	"DZ",	"drizzle",
-	"GR",	"hail",
-	"GS",	"small hail and/or snow pellets",
-	"IC",	"ice crystals",
-	"PE",	"ice pellets",
-	"RA",	"rain",
-	"SG",	"snow grains",
-	"SN",	"snow",
-	"UP",	"unknown precipitation",
-	"BR",	"mist",
-	"DU",	"widespread dust",
-	"SG",	"fog",
-	"SGBR",	"fog bank",
-	"FU",	"smoke",
-	"HZ",	"haze",
-	"PY",	"spray",
-	"SA",	"sand",
-	"VA",	"volcanic ash",
-	"DS",	"duststorm",
-	"FC",	"funnel cloud/tornado waterspout",
-	"PO",	"well-developed dust/sand whirls",
-	"SQ",	"squalls",
-	"SS",	"sandstorm",
-	"UP",	"unknown",	// ... due to failed automatic acquisition
-	0, 0
+	{ "DZ",   "drizzle" },
+	{ "GR",   "hail" },
+	{ "GS",   "small hail and/or snow pellets" },
+	{ "IC",   "ice crystals" },
+	{ "PE",   "ice pellets" },
+	{ "RA",   "rain" },
+	{ "SG",   "snow grains" },
+	{ "SN",   "snow" },
+	{ "UP",   "unknown precipitation" },
+	{ "BR",   "mist" },
+	{ "DU",   "widespread dust" },
+	{ "FG",   "fog" },
+	{ "FGBR", "fog bank" },
+	{ "FU",   "smoke" },
+	{ "HZ",   "haze" },
+	{ "PY",   "spray" },
+	{ "SA",   "sand" },
+	{ "VA",   "volcanic ash" },
+	{ "DS",   "duststorm" },
+	{ "FC",   "funnel cloud/tornado waterspout" },
+	{ "PO",   "well-developed dust/sand whirls" },
+	{ "SQ",   "squalls" },
+	{ "SS",   "sandstorm" },
+	{ "UP",   "unknown" },	// ... due to failed automatic acquisition
+	{ 0, 0 }
 };
 
 
-// (+|-|VC)?(NSW|MI|PR|BC|DR|BL|SH|TS|FZ)?((DZ|RA|SN|SG|IC|PE|GR|GS|UP){0,3})(BR|SG|FU|VA|DU|SA|HZ|PY|PO|SQ|FC|SS|DS){0,3}
+// (+|-|VC)?(NSW|MI|PR|BC|DR|BL|SH|TS|FZ)?((DZ|RA|SN|SG|IC|PE|GR|GS|UP){0,3})(BR|FG|FU|VA|DU|SA|HZ|PY|PO|SQ|FC|SS|DS){0,3}
 bool SGMetar::scanWeather()
 {
 	char *m = _m;
@@ -713,26 +714,26 @@ bool SGMetar::scanWeather()
 
 
 static const struct Token cloud_types[] = {
-	"AC",	 "altocumulus",
-	"ACC",	 "altocumulus castellanus",
-	"ACSL",	 "altocumulus standing lenticular",
-	"AS",	 "altostratus",
-	"CB",	 "cumulonimbus",
-	"CBMAM", "cumulonimbus mammatus",
-	"CC",	 "cirrocumulus",
-	"CCSL",	 "cirrocumulus standing lenticular",
-	"CI",	 "cirrus",
-	"CS",	 "cirrostratus",
-	"CU",	 "cumulus",
-	"CUFRA", "cumulus fractus",
-	"NS",	 "nimbostratus",
-	"SAC",	 "stratoaltocumulus",		// guessed
-	"SC",	 "stratocumulus",
-	"SCSL",	 "stratocumulus standing lenticular",
-	"ST",	 "stratus",
-	"STFRA", "stratus fractus",
-	"TCU",	 "towering cumulus",
-	0, 0
+	{ "AC",    "altocumulus" },
+	{ "ACC",   "altocumulus castellanus" },
+	{ "ACSL",  "altocumulus standing lenticular" },
+	{ "AS",    "altostratus" },
+	{ "CB",    "cumulonimbus" },
+	{ "CBMAM", "cumulonimbus mammatus" },
+	{ "CC",    "cirrocumulus" },
+	{ "CCSL",  "cirrocumulus standing lenticular" },
+	{ "CI",    "cirrus" },
+	{ "CS",    "cirrostratus" },
+	{ "CU",    "cumulus" },
+	{ "CUFRA", "cumulus fractus" },
+	{ "NS",    "nimbostratus" },
+	{ "SAC",   "stratoaltocumulus" },		// guessed
+	{ "SC",    "stratocumulus" },
+	{ "SCSL",  "stratocumulus standing lenticular" },
+	{ "ST",    "stratus" },
+	{ "STFRA", "stratus fractus" },
+	{ "TCU",   "towering cumulus" },
+	{ 0, 0 }
 };
 
 
@@ -957,11 +958,13 @@ bool SGMetar::scanRunwayReport()
 			m++;
 		else
 			return false;
+
 		if (*m == '1' || *m == '2' || *m == '5' || *m == '9') {	// extent of deposit
 			r._extent = *m - '0';
 			r._extent_string = runway_deposit_extent[*m - '0'];
 		} else if (*m != '/')
 			return false;
+
 		m++;
 		i = -1;
 		if (!strncmp(m, "//", 2))
@@ -970,7 +973,7 @@ bool SGMetar::scanRunwayReport()
 			return false;
 
 		if (i == 0)
-			r._depth = 0.5;					// < 1 mm deep (let's say 0.5 :-)
+			r._depth = 0.0005;				// < 1 mm deep (let's say 0.5 :-)
 		else if (i > 0 && i <= 90)
 			r._depth = i / 1000.0;				// i mm deep
 		else if (i >= 92 && i <= 98)
@@ -996,6 +999,7 @@ bool SGMetar::scanRunwayReport()
 		return false;
 
 	_runways[id]._deposit = r._deposit;
+	_runways[id]._deposit_string = r._deposit_string;
 	_runways[id]._extent = r._extent;
 	_runways[id]._extent_string = r._extent_string;
 	_runways[id]._depth = r._depth;
@@ -1075,13 +1079,13 @@ bool SGMetar::scanTrendForecast()
 
 // (BLU|WHT|GRN|YLO|AMB|RED)
 static const struct Token colors[] = {
-	"BLU", "Blue",		// 2500 ft,  8.0 km
-	"WHT", "White",		// 1500 ft,  5.0 km
-	"GRN", "Green",		//  700 ft,  3.7 km
-	"YLO", "Yellow",	//  300 ft,  1.6 km
-	"AMB", "Amber",		//  200 ft,  0.8 km
-	"RED", "Red",		// <200 ft, <0.8 km
-	0, 0
+	{ "BLU", "Blue" },      // 2500 ft,  8.0 km
+	{ "WHT", "White" },     // 1500 ft,  5.0 km
+	{ "GRN", "Green" },     //  700 ft,  3.7 km
+	{ "YLO", "Yellow" },	//  300 ft,  1.6 km
+	{ "AMB", "Amber" },     //  200 ft,  0.8 km
+	{ "RED", "Red" },       // <200 ft, <0.8 km
+	{ 0, 0 }
 };
 
 

simgear/environment/metar.hxx

@@ -16,7 +16,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/environment/visual_enviro.cxx

@@ -16,16 +16,17 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 //
-
 #ifdef HAVE_CONFIG_H
 #  include <simgear_config.h>
 #endif
 
 #include <plib/sg.h>
 #include <simgear/constants.h>
+#include <simgear/structure/SGReferenced.hxx>
+#include <simgear/structure/SGSharedPtr.hxx>
 #include <simgear/math/sg_random.h>
 #include <simgear/math/sg_geodesy.hxx>
 #include <simgear/math/point3d.hxx>
@@ -33,6 +34,7 @@
 #include <simgear/sound/soundmgr_openal.hxx>
 #include <simgear/scene/sky/cloudfield.hxx>
 #include <simgear/scene/sky/newcloud.hxx>
+#include <simgear/props/props.hxx>
 #include "visual_enviro.hxx"
 
 #include <vector>
@@ -50,6 +52,85 @@ typedef struct {
 static float rainpos[MAX_RAIN_SLICE];
 #define MAX_LT_TREE_SEG	400
 
+#define DFL_MIN_LIGHT 0.35
+sgVec3 SGEnviro::min_light = {DFL_MIN_LIGHT, DFL_MIN_LIGHT, DFL_MIN_LIGHT};
+#define DFL_STREAK_BRIGHT_NEARMOST_LAYER 0.9
+SGfloat SGEnviro::streak_bright_nearmost_layer = DFL_STREAK_BRIGHT_NEARMOST_LAYER;
+#define DFL_STREAK_BRIGHT_FARMOST_LAYER 0.5
+SGfloat SGEnviro::streak_bright_farmost_layer = DFL_STREAK_BRIGHT_FARMOST_LAYER;
+#define DFL_STREAK_PERIOD_MAX 2.5
+SGfloat SGEnviro::streak_period_max = DFL_STREAK_PERIOD_MAX;
+#define DFL_STREAK_PERIOD_CHANGE_PER_KT 0.005
+SGfloat SGEnviro::streak_period_change_per_kt = DFL_STREAK_PERIOD_CHANGE_PER_KT;
+#define DFL_STREAK_PERIOD_MIN 1.0
+SGfloat SGEnviro::streak_period_min = DFL_STREAK_PERIOD_MIN;
+#define DFL_STREAK_LENGTH_MIN 0.03
+SGfloat SGEnviro::streak_length_min = DFL_STREAK_LENGTH_MIN;
+#define DFL_STREAK_LENGTH_CHANGE_PER_KT 0.0005
+SGfloat SGEnviro::streak_length_change_per_kt = DFL_STREAK_LENGTH_CHANGE_PER_KT;
+#define DFL_STREAK_LENGTH_MAX 0.1
+SGfloat SGEnviro::streak_length_max = DFL_STREAK_LENGTH_MAX;
+#define DFL_STREAK_COUNT_MIN 40
+int SGEnviro::streak_count_min = DFL_STREAK_COUNT_MIN;
+#define DFL_STREAK_COUNT_MAX 190
+#if (DFL_STREAK_COUNT_MAX > MAX_RAIN_SLICE)
+#error "Bad default!"
+#endif
+int SGEnviro::streak_count_max = DFL_STREAK_COUNT_MAX;
+#define DFL_CONE_BASE_RADIUS 15.0
+SGfloat SGEnviro::cone_base_radius = DFL_CONE_BASE_RADIUS;
+#define DFL_CONE_HEIGHT 30.0
+SGfloat SGEnviro::cone_height = DFL_CONE_HEIGHT;
+
+
+void SGEnviro::config(const SGPropertyNode* n)
+{
+	if (!n)
+		return;
+
+	const float ml = n->getFloatValue("min-light", DFL_MIN_LIGHT);
+	sgSetVec3(min_light, ml, ml, ml);
+
+	streak_bright_nearmost_layer = n->getFloatValue(
+			"streak-brightness-nearmost-layer",
+			DFL_STREAK_BRIGHT_NEARMOST_LAYER);
+	streak_bright_farmost_layer = n->getFloatValue(
+			"streak-brightness-farmost-layer",
+			DFL_STREAK_BRIGHT_FARMOST_LAYER);
+
+	streak_period_max = n->getFloatValue(
+			"streak-period-max",
+			DFL_STREAK_PERIOD_MAX);
+	streak_period_min = n->getFloatValue(
+			"streak-period-min",
+			DFL_STREAK_PERIOD_MIN);
+	streak_period_change_per_kt = n->getFloatValue(
+			"streak-period-change-per-kt",
+			DFL_STREAK_PERIOD_CHANGE_PER_KT);
+
+	streak_length_max = n->getFloatValue(
+			"streak-length-max",
+			DFL_STREAK_LENGTH_MAX);
+	streak_length_min = n->getFloatValue(
+			"streak-length-min",
+			DFL_STREAK_LENGTH_MIN);
+	streak_length_change_per_kt = n->getFloatValue(
+			"streak-length-change-per-kt",
+			DFL_STREAK_LENGTH_CHANGE_PER_KT);
+
+	streak_count_min = n->getIntValue(
+			"streak-count-min", DFL_STREAK_COUNT_MIN);
+	streak_count_max = n->getIntValue(
+			"streak-count-max", DFL_STREAK_COUNT_MAX);
+	if (streak_count_max > MAX_RAIN_SLICE)
+		streak_count_max = MAX_RAIN_SLICE;
+
+	cone_base_radius = n->getFloatValue(
+			"cone-base-radius", DFL_CONE_BASE_RADIUS);
+	cone_height = n->getFloatValue("cone_height", DFL_CONE_HEIGHT);
+}
+
+
 /**
  * A class to render lightnings.
  */
@@ -84,23 +165,23 @@ static list_of_lightning lightnings;
 
 SGEnviro sgEnviro;
 
-SGEnviro::SGEnviro(void) :
+SGEnviro::SGEnviro() :
 	view_in_cloud(false),
-	turbulence_enable_state(false),
 	precipitation_enable_state(true),
+	precipitation_density(100.0),
+	precipitation_max_alt(0.0),
+	turbulence_enable_state(false),
+	last_cloud_turbulence(0.0),
+	cloud_turbulence(0.0),
 	lightning_enable_state(false),
+	elapsed_time(0.0),
+	dt(0.0),
 	soundMgr(NULL),
 	snd_active(false),
 	snd_dist(0.0),
-	last_cloud_turbulence(0.0),
-	cloud_turbulence(0.0),
-	elapsed_time(0.0),
-	dt(0.0),
 	min_time_before_lt(0.0),
 	fov_width(55.0),
-	fov_height(55.0),
-	precipitation_max_alt(0.0),
-	precipitation_density(100.0)
+	fov_height(55.0)
 
 {
 	for(int i = 0; i < MAX_RAIN_SLICE ; i++)
@@ -245,7 +326,7 @@ void SGEnviro::setLight(sgVec4 adj_fog_color) {
 	}
 }
 
-void SGEnviro::callback_cloud(float heading, float alt, float radius, int familly, float dist, int cloudId) {
+void SGEnviro::callback_cloud(float heading, float alt, float radius, int family, float dist, int cloudId) {
 	// send data to wx radar
 	// compute turbulence
 	// draw precipitation
@@ -255,7 +336,7 @@ void SGEnviro::callback_cloud(float heading, float alt, float radius, int famill
 	// http://www.pilotfriend.com/flight_training/weather/THUNDERSTORM%20HAZARDS1.htm
 	double turbulence = 0.0;
 	if( dist < radius * radius * 2.25f ) {
-		switch(familly) {
+		switch(family) {
 			case SGNewCloud::CLFamilly_st:
 				turbulence = 0.2;
 				break;
@@ -290,7 +371,7 @@ void SGEnviro::callback_cloud(float heading, float alt, float radius, int famill
 	// convert to LWC for radar (experimental)
 	// http://www-das.uwyo.edu/~geerts/cwx/notes/chap08/moist_cloud.html
 	double LWC = 0.0;
-	switch(familly) {
+	switch(family) {
 		case SGNewCloud::CLFamilly_st:
 			LWC = 0.29;
 			break;
@@ -324,7 +405,7 @@ void SGEnviro::callback_cloud(float heading, float alt, float radius, int famill
 	// NB:data valid only from cockpit view
 
 	// spawn a new lightning
-	if(lightning_enable_state && min_time_before_lt <= 0.0 && (familly == SGNewCloud::CLFamilly_cb) &&
+	if(lightning_enable_state && min_time_before_lt <= 0.0 && (family == SGNewCloud::CLFamilly_cb) &&
 		dist < 15000.0 * 15000.0 && sg_random() > 0.9f) {
 		double lat, lon;
 		Point3D orig, dest;
@@ -343,7 +424,7 @@ void SGEnviro::callback_cloud(float heading, float alt, float radius, int famill
 //		min_time_before_lt = 5.0;
 	}
 	if( (alt - radius * 0.1) > precipitation_max_alt )
-		switch(familly) {
+		switch(family) {
 			case SGNewCloud::CLFamilly_st:
 			case SGNewCloud::CLFamilly_cu:
 			case SGNewCloud::CLFamilly_cb:
@@ -360,26 +441,25 @@ list_of_SGWxRadarEcho *SGEnviro::get_radar_echo(void) {
 
 // precipitation rendering code
 void SGEnviro::DrawCone2(float baseRadius, float height, int slices, bool down, double rain_norm, double speed) {
-
 	sgVec3 light;
-	sgVec3 min_light = {0.35, 0.35, 0.35};
 	sgAddVec3( light, fog_color, min_light );
 	float da = SG_PI * 2.0f / (float) slices;
 	// low number = faster
-	float speedf = 2.5f - speed / 200.0;
-	if( speedf < 1.0f )
-		speedf = 1.0f;
-	float lenf = 0.03f + speed / 2000.0;
-	if( lenf > 0.10f )
-		lenf = 0.10f;
+	float speedf = streak_period_max - speed * streak_period_change_per_kt;
+	if( speedf < streak_period_min )
+		speedf = streak_period_min;
+	float lenf = streak_length_min + speed * streak_length_change_per_kt;
+	if( lenf > streak_length_max )
+		lenf = streak_length_max;
     float t = fmod((float) elapsed_time, speedf) / speedf;
 //	t = 0.1f;
 	if( !down )
 		t = 1.0f - t;
 	float angle = 0.0f;
-	glColor4f(1.0f, 0.7f, 0.7f, 0.9f);
+	//glColor4f(1.0f, 0.7f, 0.7f, 0.9f); // XXX unneeded? overriden below
 	glBegin(GL_LINES);
-	int rainpos_indice = 0;
+	if (slices >  MAX_RAIN_SLICE)
+		slices = MAX_RAIN_SLICE; // should never happen
 	for( int i = 0 ; i < slices ; i++ ) {
 		float x = cos(angle) * baseRadius;
 		float y = sin(angle) * baseRadius;
@@ -387,12 +467,14 @@ void SGEnviro::DrawCone2(float baseRadius, float height, int slices, bool down,
 		sgVec3 dir = {x, -height, y};
 
 		// rain drops at 2 different speed to simulate depth
-		float t1 = (i & 1 ? t : t + t) + rainpos[rainpos_indice];
+		float t1 = (i & 1 ? t : t + t) + rainpos[i];
 		if(t1 > 1.0f)	t1 -= 1.0f;
 		if(t1 > 1.0f)	t1 -= 1.0f;
 
 		// distant raindrops are more transparent
-		float c = (i & 1 ? t1 * 0.5f : t1 * 0.9f);
+		float c = t1 * (i & 1 ?
+				streak_bright_farmost_layer
+				: streak_bright_nearmost_layer);
 		glColor4f(c * light[0], c * light[1], c * light[2], c);
 		sgVec3 p1, p2;
 		sgScaleVec3(p1, dir, t1);
@@ -402,13 +484,28 @@ void SGEnviro::DrawCone2(float baseRadius, float height, int slices, bool down,
 
 		glVertex3f(p1[0], p1[1] + height, p1[2]);
 		glVertex3f(p2[0], p2[1] + height, p2[2]);
-		if( ++rainpos_indice >= MAX_RAIN_SLICE )
-			rainpos_indice = 0;
 	}
 	glEnd();
 }
 
-void SGEnviro::drawRain(double pitch, double roll, double heading, double speed, double rain_norm) {
+void SGEnviro::drawRain(double pitch, double roll, double heading, double hspeed, double rain_norm) {
+
+#if 0
+	static int debug_period = 0;
+	if (debug_period++ == 50) {
+		debug_period = 0;
+		cout << "drawRain("
+			<< pitch << ", "
+			<< roll  << ", "
+			<< heading << ", "
+			<< hspeed << ", "
+			<< rain_norm << ");"
+			//" angle = " << angle
+			//<< " raindrop(KTS) = " << raindrop_speed_kts
+			<< endl;
+	}
+#endif
+
 
 	glBindTexture(GL_TEXTURE_2D, 0);
 
@@ -419,25 +516,33 @@ void SGEnviro::drawRain(double pitch, double roll, double heading, double speed,
 	glDisable( GL_FOG );
 	glDisable(GL_LIGHTING);
 
-	int slice_count = (40.0 + rain_norm*150.0)* precipitation_density / 100.0;
+	int slice_count = static_cast<int>(
+			(streak_count_min + rain_norm*(streak_count_max-streak_count_min))
+				* precipitation_density / 100.0);
 
-	float angle = speed;
-	if( angle > 90.0 )
-		angle = 90.0;
+	// www.wonderquest.com/falling-raindrops.htm says that
+	// Raindrop terminal velocity is 5 to 20mph
+	// Rather than model it accurately (temp, pressure, diameter), and make it
+	// smaller than terminal when closer to the precipitation cloud base,
+	// we interpolate in the 5-20mph range according to rain_norm.
+	double raindrop_speed_kts
+		= (5.0 + rain_norm*15.0) * SG_MPH_TO_MPS * SG_MPS_TO_KT;
 
+	float angle = atanf(hspeed / raindrop_speed_kts) * SG_RADIANS_TO_DEGREES;
 	glPushMatrix();
-		// TODO:find the real view orientation, not the AC one
-		// the cone rotate with speed
+		// the cone rotate with hspeed
 		angle = -pitch - angle;
+		glRotatef(roll, 0.0, 0.0, 1.0);
+		glRotatef(heading, 0.0, 1.0, 0.0);
 		glRotatef(angle, 1.0, 0.0, 0.0);
-		glRotatef(roll, 0.0, 1.0, 0.0);
-		glRotatef(heading, 0.0, 0.0, 1.0);
 
 		// up cone
-		DrawCone2(15.0, 30.0, slice_count, true, rain_norm, speed);
+		DrawCone2(cone_base_radius, cone_height, 
+				slice_count, true, rain_norm, hspeed);
 		// down cone (usually not visible)
 		if(angle > 0.0 || heading != 0.0)
-			DrawCone2(15.0, -30.0, slice_count, false, rain_norm, speed);
+			DrawCone2(cone_base_radius, -cone_height, 
+					slice_count, false, rain_norm, hspeed);
 
 	glPopMatrix();
 
@@ -452,19 +557,19 @@ void SGEnviro::set_soundMgr(SGSoundMgr *mgr) {
 	soundMgr = mgr;
 }
 
-void SGEnviro::drawPrecipitation(double rain_norm, double snow_norm, double hail_norm, double pitch, double roll, double heading, double speed) {
+void SGEnviro::drawPrecipitation(double rain_norm, double snow_norm, double hail_norm, double pitch, double roll, double heading, double hspeed) {
 	if( precipitation_enable_state && rain_norm > 0.0)
 	  if( precipitation_max_alt >= last_alt )
-		drawRain(pitch, roll, heading, speed, rain_norm);
+		drawRain(pitch, roll, heading, hspeed, rain_norm);
 }
 
 
 SGLightning::SGLightning(double _lon, double _lat, double _alt) :
+	nb_tree(0),
 	lon(_lon),
 	lat(_lat),
 	alt(_alt),
-	age(1.0 + sg_random() * 4.0),
-	nb_tree(0)
+	age(1.0 + sg_random() * 4.0)
 {
 //	sequence_count = 1 + sg_random() * 5.0;
 	lt_build();
@@ -500,7 +605,7 @@ void SGLightning::lt_build_tree_branch(int tree_nr, Point3D &start, float energy
         nseg++;
 		// add a branch
         if( energy * sg_random() > 0.8f )
-			lt_build_tree_branch(tree_nr + 1, pt, energy * 0.9f, nbseg == 50 ? 10 : nbseg * 0.4f, segsize * 0.7f);
+			lt_build_tree_branch(tree_nr + 1, pt, energy * 0.9f, nbseg == 50 ? 10 : static_cast<int>(nbseg * 0.4f), segsize * 0.7f);
 
 		if( nb_tree >= MAX_LT_TREE_SEG )
 			return;
@@ -555,11 +660,12 @@ void SGLightning::lt_Render(void) {
 	sgVec4 c;
 
 #define DRAW_SEG() \
-			{glBegin(GL_LINES); \
+			{glColorMaterial(GL_FRONT, GL_EMISSION);  \
+			glDisable(GL_LINE_SMOOTH); glBegin(GL_LINES); \
 				glColor4fv(c); \
                 glVertex3f(lt_tree[n].pt[PX], lt_tree[n].pt[PZ], lt_tree[n].pt[PY]); \
                 glVertex3f(lt_tree[lt_tree[n].prev].pt[PX], lt_tree[lt_tree[n].prev].pt[PZ], lt_tree[lt_tree[n].prev].pt[PY]); \
-			glEnd();}
+			glEnd(); glEnable(GL_LINE_SMOOTH);}
 
 	glDepthMask( GL_FALSE );
 	glEnable(GL_BLEND);
@@ -654,7 +760,7 @@ void SGEnviro::drawLightning(void) {
 				double ax = 0.0, ay = 0.0;
 				ax = cos(course) * dist;
 				ay = sin(course) * dist;
-				SGSoundSample *snd = soundMgr->find("thunder");
+				SGSharedPtr<SGSoundSample> snd = soundMgr->find("thunder");
 				if( snd ) {
 					ALfloat pos[3]={ax, ay, -sgEnviro.last_alt };
 					snd->set_source_pos(pos);

simgear/environment/visual_enviro.hxx

@@ -16,7 +16,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 //
 #ifndef _VISUAL_ENVIRO_HXX
@@ -90,11 +90,31 @@ private:
 
 	/** a list of all the radar echo. */
 	list_of_SGWxRadarEcho radarEcho;
+	static sgVec3 min_light;
+	static SGfloat streak_bright_nearmost_layer,
+				   streak_bright_farmost_layer,
+				   streak_period_max,
+				   streak_period_change_per_kt,
+				   streak_period_min,
+				   streak_length_min,
+				   streak_length_change_per_kt,
+				   streak_length_max;
+	static int streak_count_min, streak_count_max;
+	static SGfloat cone_base_radius,
+				   cone_height;
 
 public:
 	SGEnviro();
 	~SGEnviro();
 
+	/** Read the config from the precipitation rendering config properties.
+	 * @param precip_rendering_cfgNode "/sim/rendering/precipitation" in fg
+	 * Set from whatever info present in the
+	 * subnodes passed, substituting hardwired defaults for missing fields.
+	 * If NULL is given, do nothing.
+	 */
+	void config(const class SGPropertyNode* precip_rendering_cfgNode);
+
     /**
      * Forward a few states used for renderings.
      */
@@ -107,12 +127,12 @@ public:
      * @param heading direction of cloud in radians
      * @param alt asl of cloud in meters
      * @param radius radius of cloud in meters
-     * @param familly cloud familly
+     * @param family cloud family
      * @param dist  squared dist to cloud in meters
      */
-	void callback_cloud(float heading, float alt, float radius, int familly, float dist, int cloudId);
+	void callback_cloud(float heading, float alt, float radius, int family, float dist, int cloudId);
 
-	void drawRain(double pitch, double roll, double heading, double speed, double rain_norm);
+	void drawRain(double pitch, double roll, double heading, double hspeed, double rain_norm);
     /**
      * Draw rain or snow precipitation around the viewer.
      * @param rain_norm rain normalized intensity given by metar class
@@ -120,10 +140,10 @@ public:
      * @param hail_norm hail normalized intensity given by metar class
      * @param pitch pitch rotation of viewer
      * @param roll roll rotation of viewer
-     * @param speed moving speed of viewer in kt
+     * @param hspeed moving horizontal speed of viewer in kt
      */
 	void drawPrecipitation(double rain_norm, double snow_norm, double hail_norm,
-							double pitch, double roll, double heading, double speed);
+							double pitch, double roll, double heading, double hspeed);
 
     /**
      * Draw the lightnings spawned by cumulo nimbus.
@@ -186,6 +206,12 @@ public:
 	float get_precipitation_density(void) const;
 	bool get_precipitation_enable_state(void) const;
 
+	/** 
+	 * Decrease the precipitation density to the given percentage.
+	 * (Only show the given percentage of rain streaks etc.)
+	 * Default precipitation density upon construction is 100.0.
+	 * @param density 0.0 to 100.0
+	 */
 	void set_precipitation_density(float density);
     /**
      * Enable or disable the rendering of precipitation around the viewer.

simgear/ephemeris/celestialBody.cxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/celestialBody.hxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/
@@ -27,9 +26,9 @@
 #ifndef _CELESTIALBODY_H_
 #define _CELESTIALBODY_H_
 
-#ifndef __cplusplus                                                          
+#ifndef __cplusplus
 # error This library requires C++
-#endif                                   
+#endif
 
 
 #include <simgear/constants.h>

simgear/ephemeris/ephemeris.cxx

@@ -15,10 +15,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/ephemeris/ephemeris.hxx

@@ -19,10 +19,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/ephemeris/jupiter.cxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/jupiter.hxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/mars.cxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/mars.hxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/mercury.cxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/mercury.hxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/moonpos.cxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/moonpos.hxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/neptune.cxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/neptune.hxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/pluto.hxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/saturn.cxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/saturn.hxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/star.cxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/star.hxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/stardata.cxx

@@ -14,10 +14,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/ephemeris/stardata.hxx

@@ -14,10 +14,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/ephemeris/uranus.cxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/uranus.hxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/venus.cxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/ephemeris/venus.hxx

@@ -15,10 +15,9 @@
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  * Library General Public License for more details.
  *
- * You should have received a copy of the GNU Library General Public
- * License along with this library; if not, write to the
- * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
- * Boston, MA  02111-1307, USA.
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  *
  * $Id$
  **************************************************************************/

simgear/io/iochannel.cxx

@@ -16,7 +16,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/io/iochannel.hxx

@@ -19,7 +19,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/io/lowlevel.cxx

@@ -18,7 +18,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/io/lowlevel.hxx

@@ -18,7 +18,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 //

simgear/io/sg_binobj.cxx

@@ -16,7 +16,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 //
@@ -36,6 +36,7 @@
 #include STL_STRING
 
 #include <simgear/bucket/newbucket.hxx>
+#include <simgear/misc/sg_path.hxx>
 
 #include "lowlevel.hxx"
 #include "sg_binobj.hxx"
@@ -45,7 +46,7 @@ SG_USING_STD( string );
 SG_USING_STD( vector );
 
 
-enum {
+enum sgObjectTypes {
     SG_BOUNDING_SPHERE = 0,
 
     SG_VERTEX_LIST = 1,
@@ -58,20 +59,19 @@ enum {
     SG_TRIANGLE_FACES = 10,
     SG_TRIANGLE_STRIPS = 11,
     SG_TRIANGLE_FANS = 12
-} sgObjectTypes;
+};
 
-enum {
+enum sgIndexTypes {
     SG_IDX_VERTICES =  0x01,
     SG_IDX_NORMALS =   0x02,
     SG_IDX_COLORS =    0x04,
     SG_IDX_TEXCOORDS = 0x08
-} sgIndexTypes;
+};
 
-enum {
+enum sgPropertyTypes {
     SG_MATERIAL = 0,
     SG_INDEX_TYPES = 1
-} sgPropertyTypes;
-
+};
 
 
 class sgSimpleBuffer {
@@ -620,20 +620,13 @@ bool SGBinObject::write_bin( const string& base, const string& name,
     unsigned char idx_mask;
     int idx_size;
 
-    string dir = base + "/" + b.gen_base_path();
-    string command = "mkdir -p " + dir;
-#if defined(_MSC_VER) || defined(__MINGW32__)
-    system( (string("mkdir ") + dir).c_str() );
-#else
-    system(command.c_str());
-#endif
-
-    string file = dir + "/" + name + ".gz";
-    cout << "Output file = " << file << endl;
+    SGPath file = base + "/" + b.gen_base_path() + "/" + name + ".gz";
+    file.create_dir( 0755 );
+    cout << "Output file = " << file.str() << endl;
 
     gzFile fp;
     if ( (fp = gzopen( file.c_str(), "wb9" )) == NULL ) {
-	cout << "ERROR: opening " << file << " for writing!" << endl;
+	cout << "ERROR: opening " << file.str() << " for writing!" << endl;
 	return false;
     }
 
@@ -858,7 +851,8 @@ bool SGBinObject::write_bin( const string& base, const string& name,
 	    // find next group
 	    material = tri_materials[start];
 	    while ( (end < (int)tri_materials.size()) && 
-		    (material == tri_materials[end]) )
+		    (material == tri_materials[end]) &&
+		    3*(end-start) < 32760 )
 	    {
 		// cout << "end = " << end << endl;
 		end++;
@@ -1051,21 +1045,13 @@ bool SGBinObject::write_ascii( const string& base, const string& name,
     Point3D p;
     int i, j;
 
-    string dir = base + "/" + b.gen_base_path();
-    string command = "mkdir -p " + dir;
-#if defined(_MSC_VER) || defined(__MINGW32__)
-    system( (string("mkdir ") + dir).c_str() );
-#else
-    system(command.c_str());
-#endif
-
-    // string file = dir + "/" + b.gen_index_str();
-    string file = dir + "/" + name;
-    cout << "Output file = " << file << endl;
+    SGPath file = base + "/" + b.gen_base_path() + "/" + name;
+    file.create_dir( 0755 );
+    cout << "Output file = " << file.str() << endl;
 
     FILE *fp;
     if ( (fp = fopen( file.c_str(), "w" )) == NULL ) {
-	cout << "ERROR: opening " << file << " for writing!" << endl;
+	cout << "ERROR: opening " << file.str() << " for writing!" << endl;
 	return false;
     }
 
@@ -1224,7 +1210,7 @@ bool SGBinObject::write_ascii( const string& base, const string& name,
     // close the file
     fclose(fp);
 
-    command = "gzip --force --best " + file;
+    string command = "gzip --force --best " + file.str();
     system(command.c_str());
 
     return true;

simgear/io/sg_binobj.hxx

@@ -19,7 +19,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/io/sg_file.cxx

@@ -16,7 +16,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/io/sg_file.hxx

@@ -18,7 +18,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/io/sg_serial.cxx

@@ -16,7 +16,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/io/sg_serial.hxx

@@ -19,7 +19,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/io/sg_socket.cxx

@@ -17,10 +17,13 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 
+#ifdef HAVE_CONFIG_H
+#  include <simgear_config.h>
+#endif
 
 #include <simgear/compiler.h>
 

simgear/io/sg_socket.hxx

@@ -19,7 +19,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/io/sg_socket_udp.cxx

@@ -16,7 +16,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/io/sg_socket_udp.hxx

@@ -19,7 +19,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/magvar/coremag.cxx

@@ -57,10 +57,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/magvar/coremag.hxx

@@ -26,10 +26,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/magvar/magvar.cxx

@@ -16,7 +16,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/magvar/magvar.hxx

@@ -19,7 +19,7 @@
 //
 // You should have received a copy of the GNU General Public License
 // along with this program; if not, write to the Free Software
-// Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 
@@ -28,9 +28,9 @@
 #define _MAGVAR_HXX
 
 
-#ifndef __cplusplus                                                          
+#ifndef __cplusplus
 # error This library requires C++
-#endif                                   
+#endif
 
 
 /**

simgear/math/.cvsignore

@@ -1,3 +1,4 @@
 .deps
 Makefile
 Makefile.in
+SGMathTest

simgear/math/Makefile.am

@@ -1,11 +1,18 @@
 includedir = @includedir@/math
 
+noinst_PROGRAMS = SGMathTest
+check_PROGRAMS  = SGMathTest
+TESTS           = $(check_PROGRAMS)
+
+SGMathTest_SOURCES = SGMathTest.cxx
+SGMathTest_LDADD = libsgmath.a $(base_LIBS)
+
+
 lib_LIBRARIES = libsgmath.a
 
 include_HEADERS = \
 	interpolater.hxx \
 	leastsqs.hxx \
-	localconsts.hxx \
 	point3d.hxx \
 	polar3d.hxx \
 	sg_geodesy.hxx \
@@ -13,7 +20,18 @@ include_HEADERS = \
 	sg_random.h \
 	sg_types.hxx \
 	vector.hxx \
-	fastmath.hxx
+	fastmath.hxx \
+	SGGeoc.hxx \
+	SGGeod.hxx \
+	SGGeodesy.hxx \
+	SGLimits.hxx \
+	SGMatrix.hxx \
+	SGMath.hxx \
+	SGMisc.hxx \
+	SGQuat.hxx \
+	SGVec4.hxx \
+	SGVec3.hxx
+
 
 EXTRA_DIST = linintp2.h linintp2.inl sphrintp.h sphrintp.inl
 
@@ -24,6 +42,7 @@ libsgmath_a_SOURCES = \
 	sg_geodesy.cxx \
 	sg_random.c \
 	vector.cxx \
-	fastmath.cxx
+	fastmath.cxx \
+	SGGeodesy.cxx
 
 INCLUDES = -I$(top_srcdir)

simgear/math/SGGeoc.hxx

@@ -0,0 +1,295 @@
+// Copyright (C) 2006  Mathias Froehlich - Mathias.Froehlich@web.de
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Library General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// This library 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
+// Library General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+
+#ifndef SGGeoc_H
+#define SGGeoc_H
+
+#include <simgear/constants.h>
+ 
+// #define SG_GEOC_NATIVE_DEGREE
+
+/// Class representing a geocentric location
+class SGGeoc {
+public:
+  /// Default constructor, initializes the instance to lat = lon = lat = 0
+  SGGeoc(void);
+  /// Initialize from a cartesian vector assumed to be in meters
+  /// Note that this conversion is relatively expensive to compute
+  SGGeoc(const SGVec3<double>& cart);
+  /// Initialize from a geodetic position
+  /// Note that this conversion is relatively expensive to compute
+  SGGeoc(const SGGeod& geod);
+
+  /// Factory from angular values in radians and radius in ft
+  static SGGeoc fromRadFt(double lon, double lat, double radius);
+  /// Factory from angular values in degrees and radius in ft
+  static SGGeoc fromDegFt(double lon, double lat, double radius);
+  /// Factory from angular values in radians and radius in m
+  static SGGeoc fromRadM(double lon, double lat, double radius);
+  /// Factory from angular values in degrees and radius in m
+  static SGGeoc fromDegM(double lon, double lat, double radius);
+
+  /// Return the geocentric longitude in radians
+  double getLongitudeRad(void) const;
+  /// Set the geocentric longitude from the argument given in radians
+  void setLongitudeRad(double lon);
+
+  /// Return the geocentric longitude in degrees
+  double getLongitudeDeg(void) const;
+  /// Set the geocentric longitude from the argument given in degrees
+  void setLongitudeDeg(double lon);
+
+  /// Return the geocentric latitude in radians
+  double getLatitudeRad(void) const;
+  /// Set the geocentric latitude from the argument given in radians
+  void setLatitudeRad(double lat);
+
+  /// Return the geocentric latitude in degrees
+  double getLatitudeDeg(void) const;
+  /// Set the geocentric latitude from the argument given in degrees
+  void setLatitudeDeg(double lat);
+
+  /// Return the geocentric radius in meters
+  double getRadiusM(void) const;
+  /// Set the geocentric radius from the argument given in meters
+  void setRadiusM(double radius);
+
+  /// Return the geocentric radius in feet
+  double getRadiusFt(void) const;
+  /// Set the geocentric radius from the argument given in feet
+  void setRadiusFt(double radius);
+
+private:
+  /// This one is private since construction is not unique if you do
+  /// not know the units of the arguments, use the factory methods for
+  /// that purpose
+  SGGeoc(double lon, double lat, double radius);
+
+  /// The actual data, angles in degree, radius in meters
+  /// The rationale for storing the values in degrees is that most code places
+  /// in flightgear/terragear use degrees as a nativ input and output value.
+  /// The places where it makes sense to use radians is when we convert
+  /// to other representations or compute rotation matrices. But both tasks
+  /// are computionally intensive anyway and that additional 'toRadian'
+  /// conversion does not hurt too much
+  double _lon;
+  double _lat;
+  double _radius;
+};
+
+inline
+SGGeoc::SGGeoc(void) :
+  _lon(0), _lat(0), _radius(0)
+{
+}
+
+inline
+SGGeoc::SGGeoc(double lon, double lat, double radius) :
+  _lon(lon), _lat(lat), _radius(radius)
+{
+}
+
+inline
+SGGeoc::SGGeoc(const SGVec3<double>& cart)
+{
+  SGGeodesy::SGCartToGeoc(cart, *this);
+}
+
+inline
+SGGeoc::SGGeoc(const SGGeod& geod)
+{
+  SGVec3<double> cart;
+  SGGeodesy::SGGeodToCart(geod, cart);
+  SGGeodesy::SGCartToGeoc(cart, *this);
+}
+
+inline
+SGGeoc
+SGGeoc::fromRadFt(double lon, double lat, double radius)
+{
+#ifdef SG_GEOC_NATIVE_DEGREE
+  return SGGeoc(lon*SGD_RADIANS_TO_DEGREES, lat*SGD_RADIANS_TO_DEGREES,
+                radius*SG_FEET_TO_METER);
+#else
+  return SGGeoc(lon, lat, radius*SG_FEET_TO_METER);
+#endif
+}
+
+inline
+SGGeoc
+SGGeoc::fromDegFt(double lon, double lat, double radius)
+{
+#ifdef SG_GEOC_NATIVE_DEGREE
+  return SGGeoc(lon, lat, radius*SG_FEET_TO_METER);
+#else
+  return SGGeoc(lon*SGD_DEGREES_TO_RADIANS, lat*SGD_DEGREES_TO_RADIANS,
+                radius*SG_FEET_TO_METER);
+#endif
+}
+
+inline
+SGGeoc
+SGGeoc::fromRadM(double lon, double lat, double radius)
+{
+#ifdef SG_GEOC_NATIVE_DEGREE
+  return SGGeoc(lon*SGD_RADIANS_TO_DEGREES, lat*SGD_RADIANS_TO_DEGREES,
+                radius);
+#else
+  return SGGeoc(lon, lat, radius);
+#endif
+}
+
+inline
+SGGeoc
+SGGeoc::fromDegM(double lon, double lat, double radius)
+{
+#ifdef SG_GEOC_NATIVE_DEGREE
+  return SGGeoc(lon, lat, radius);
+#else
+  return SGGeoc(lon*SGD_DEGREES_TO_RADIANS, lat*SGD_DEGREES_TO_RADIANS,
+                radius);
+#endif
+}
+
+inline
+double
+SGGeoc::getLongitudeRad(void) const
+{
+#ifdef SG_GEOC_NATIVE_DEGREE
+  return _lon*SGD_DEGREES_TO_RADIANS;
+#else
+  return _lon;
+#endif
+}
+
+inline
+void
+SGGeoc::setLongitudeRad(double lon)
+{
+#ifdef SG_GEOC_NATIVE_DEGREE
+  _lon = lon*SGD_RADIANS_TO_DEGREES;
+#else
+  _lon = lon;
+#endif
+}
+
+inline
+double
+SGGeoc::getLongitudeDeg(void) const
+{
+#ifdef SG_GEOC_NATIVE_DEGREE
+  return _lon;
+#else
+  return _lon*SGD_DEGREES_TO_RADIANS;
+#endif
+}
+
+inline
+void
+SGGeoc::setLongitudeDeg(double lon)
+{
+#ifdef SG_GEOC_NATIVE_DEGREE
+  _lon = lon;
+#else
+  _lon = lon*SGD_RADIANS_TO_DEGREES;
+#endif
+}
+
+inline
+double
+SGGeoc::getLatitudeRad(void) const
+{
+#ifdef SG_GEOC_NATIVE_DEGREE
+  return _lat*SGD_DEGREES_TO_RADIANS;
+#else
+  return _lat;
+#endif
+}
+
+inline
+void
+SGGeoc::setLatitudeRad(double lat)
+{
+#ifdef SG_GEOC_NATIVE_DEGREE
+  _lat = lat*SGD_RADIANS_TO_DEGREES;
+#else
+  _lat = lat;
+#endif
+}
+
+inline
+double
+SGGeoc::getLatitudeDeg(void) const
+{
+#ifdef SG_GEOC_NATIVE_DEGREE
+  return _lat;
+#else
+  return _lat*SGD_DEGREES_TO_RADIANS;
+#endif
+}
+
+inline
+void
+SGGeoc::setLatitudeDeg(double lat)
+{
+#ifdef SG_GEOC_NATIVE_DEGREE
+  _lat = lat;
+#else
+  _lat = lat*SGD_RADIANS_TO_DEGREES;
+#endif
+}
+
+inline
+double
+SGGeoc::getRadiusM(void) const
+{
+  return _radius;
+}
+
+inline
+void
+SGGeoc::setRadiusM(double radius)
+{
+  _radius = radius;
+}
+
+inline
+double
+SGGeoc::getRadiusFt(void) const
+{
+  return _radius*SG_METER_TO_FEET;
+}
+
+inline
+void
+SGGeoc::setRadiusFt(double radius)
+{
+  _radius = radius*SG_FEET_TO_METER;
+}
+
+/// Output to an ostream
+template<typename char_type, typename traits_type>
+inline
+std::basic_ostream<char_type, traits_type>&
+operator<<(std::basic_ostream<char_type, traits_type>& s, const SGGeoc& g)
+{
+  return s << "lon = " << g.getLongitudeDeg()
+           << ", lat = " << g.getLatitudeDeg()
+           << ", radius = " << g.getRadiusM();
+}
+
+#endif

simgear/math/SGGeod.hxx

@@ -0,0 +1,322 @@
+// Copyright (C) 2006  Mathias Froehlich - Mathias.Froehlich@web.de
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Library General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// This library 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
+// Library General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+
+#ifndef SGGeod_H
+#define SGGeod_H
+
+#include <simgear/constants.h>
+
+// #define SG_GEOD_NATIVE_DEGREE
+
+/// Class representing a geodetic location
+class SGGeod {
+public:
+  /// Default constructor, initializes the instance to lat = lon = elev = 0
+  SGGeod(void);
+  /// Initialize from a cartesian vector assumed to be in meters
+  /// Note that this conversion is relatively expensive to compute
+  SGGeod(const SGVec3<double>& cart);
+  /// Initialize from a geocentric position
+  /// Note that this conversion is relatively expensive to compute
+  SGGeod(const SGGeoc& geoc);
+
+  /// Factory from angular values in radians and elevation is 0
+  static SGGeod fromRad(double lon, double lat);
+  /// Factory from angular values in degrees and elevation is 0
+  static SGGeod fromDeg(double lon, double lat);
+  /// Factory from angular values in radians and elevation in ft
+  static SGGeod fromRadFt(double lon, double lat, double elevation);
+  /// Factory from angular values in degrees and elevation in ft
+  static SGGeod fromDegFt(double lon, double lat, double elevation);
+  /// Factory from angular values in radians and elevation in m
+  static SGGeod fromRadM(double lon, double lat, double elevation);
+  /// Factory from angular values in degrees and elevation in m
+  static SGGeod fromDegM(double lon, double lat, double elevation);
+
+  /// Return the geodetic longitude in radians
+  double getLongitudeRad(void) const;
+  /// Set the geodetic longitude from the argument given in radians
+  void setLongitudeRad(double lon);
+
+  /// Return the geodetic longitude in degrees
+  double getLongitudeDeg(void) const;
+  /// Set the geodetic longitude from the argument given in degrees
+  void setLongitudeDeg(double lon);
+
+  /// Return the geodetic latitude in radians
+  double getLatitudeRad(void) const;
+  /// Set the geodetic latitude from the argument given in radians
+  void setLatitudeRad(double lat);
+
+  /// Return the geodetic latitude in degrees
+  double getLatitudeDeg(void) const;
+  /// Set the geodetic latitude from the argument given in degrees
+  void setLatitudeDeg(double lat);
+
+  /// Return the geodetic elevation in meters
+  double getElevationM(void) const;
+  /// Set the geodetic elevation from the argument given in meters
+  void setElevationM(double elevation);
+
+  /// Return the geodetic elevation in feet
+  double getElevationFt(void) const;
+  /// Set the geodetic elevation from the argument given in feet
+  void setElevationFt(double elevation);
+
+private:
+  /// This one is private since construction is not unique if you do
+  /// not know the units of the arguments. Use the factory methods for
+  /// that purpose
+  SGGeod(double lon, double lat, double elevation);
+
+  /// The actual data, angles in degree, elevation in meters
+  /// The rationale for storing the values in degrees is that most code places
+  /// in flightgear/terragear use degrees as a nativ input and output value.
+  /// The places where it makes sense to use radians is when we convert
+  /// to other representations or compute rotation matrices. But both tasks
+  /// are computionally intensive anyway and that additional 'toRadian'
+  /// conversion does not hurt too much
+  double _lon;
+  double _lat;
+  double _elevation;
+};
+
+inline
+SGGeod::SGGeod(void) :
+  _lon(0), _lat(0), _elevation(0)
+{
+}
+
+inline
+SGGeod::SGGeod(double lon, double lat, double elevation) :
+  _lon(lon), _lat(lat), _elevation(elevation)
+{
+}
+
+inline
+SGGeod::SGGeod(const SGVec3<double>& cart)
+{
+  SGGeodesy::SGCartToGeod(cart, *this);
+}
+
+inline
+SGGeod::SGGeod(const SGGeoc& geoc)
+{
+  SGVec3<double> cart;
+  SGGeodesy::SGGeocToCart(geoc, cart);
+  SGGeodesy::SGCartToGeod(cart, *this);
+}
+
+inline
+SGGeod
+SGGeod::fromRad(double lon, double lat)
+{
+#ifdef SG_GEOD_NATIVE_DEGREE
+  return SGGeod(lon*SGD_RADIANS_TO_DEGREES, lat*SGD_RADIANS_TO_DEGREES, 0);
+#else
+  return SGGeod(lon, lat, 0);
+#endif
+}
+
+inline
+SGGeod
+SGGeod::fromDeg(double lon, double lat)
+{
+#ifdef SG_GEOD_NATIVE_DEGREE
+  return SGGeod(lon, lat, 0);
+#else
+  return SGGeod(lon*SGD_DEGREES_TO_RADIANS, lat*SGD_DEGREES_TO_RADIANS, 0);
+#endif
+}
+
+inline
+SGGeod
+SGGeod::fromRadFt(double lon, double lat, double elevation)
+{
+#ifdef SG_GEOD_NATIVE_DEGREE
+  return SGGeod(lon*SGD_RADIANS_TO_DEGREES, lat*SGD_RADIANS_TO_DEGREES,
+                elevation*SG_FEET_TO_METER);
+#else
+  return SGGeod(lon, lat, elevation*SG_FEET_TO_METER);
+#endif
+}
+
+inline
+SGGeod
+SGGeod::fromDegFt(double lon, double lat, double elevation)
+{
+#ifdef SG_GEOD_NATIVE_DEGREE
+  return SGGeod(lon, lat, elevation*SG_FEET_TO_METER);
+#else
+  return SGGeod(lon*SGD_DEGREES_TO_RADIANS, lat*SGD_DEGREES_TO_RADIANS,
+                elevation*SG_FEET_TO_METER);
+#endif
+}
+
+inline
+SGGeod
+SGGeod::fromRadM(double lon, double lat, double elevation)
+{
+#ifdef SG_GEOD_NATIVE_DEGREE
+  return SGGeod(lon*SGD_RADIANS_TO_DEGREES, lat*SGD_RADIANS_TO_DEGREES,
+                elevation);
+#else
+  return SGGeod(lon, lat, elevation);
+#endif
+}
+
+inline
+SGGeod
+SGGeod::fromDegM(double lon, double lat, double elevation)
+{
+#ifdef SG_GEOD_NATIVE_DEGREE
+  return SGGeod(lon, lat, elevation);
+#else
+  return SGGeod(lon*SGD_DEGREES_TO_RADIANS, lat*SGD_DEGREES_TO_RADIANS,
+                elevation);
+#endif
+}
+
+inline
+double
+SGGeod::getLongitudeRad(void) const
+{
+#ifdef SG_GEOD_NATIVE_DEGREE
+  return _lon*SGD_DEGREES_TO_RADIANS;
+#else
+  return _lon;
+#endif
+}
+
+inline
+void
+SGGeod::setLongitudeRad(double lon)
+{
+#ifdef SG_GEOD_NATIVE_DEGREE
+  _lon = lon*SGD_RADIANS_TO_DEGREES;
+#else
+  _lon = lon;
+#endif
+}
+
+inline
+double
+SGGeod::getLongitudeDeg(void) const
+{
+#ifdef SG_GEOD_NATIVE_DEGREE
+  return _lon;
+#else
+  return _lon*SGD_RADIANS_TO_DEGREES;
+#endif
+}
+
+inline
+void
+SGGeod::setLongitudeDeg(double lon)
+{
+#ifdef SG_GEOD_NATIVE_DEGREE
+  _lon = lon;
+#else
+  _lon = lon*SGD_DEGREES_TO_RADIANS;
+#endif
+}
+
+inline
+double
+SGGeod::getLatitudeRad(void) const
+{
+#ifdef SG_GEOD_NATIVE_DEGREE
+  return _lat*SGD_DEGREES_TO_RADIANS;
+#else
+  return _lat;
+#endif
+}
+
+inline
+void
+SGGeod::setLatitudeRad(double lat)
+{
+#ifdef SG_GEOD_NATIVE_DEGREE
+  _lat = lat*SGD_RADIANS_TO_DEGREES;
+#else
+  _lat = lat;
+#endif
+}
+
+inline
+double
+SGGeod::getLatitudeDeg(void) const
+{
+#ifdef SG_GEOD_NATIVE_DEGREE
+  return _lat;
+#else
+  return _lat*SGD_RADIANS_TO_DEGREES;
+#endif
+}
+
+inline
+void
+SGGeod::setLatitudeDeg(double lat)
+{
+#ifdef SG_GEOD_NATIVE_DEGREE
+  _lat = lat;
+#else
+  _lat = lat*SGD_DEGREES_TO_RADIANS;
+#endif
+}
+
+inline
+double
+SGGeod::getElevationM(void) const
+{
+  return _elevation;
+}
+
+inline
+void
+SGGeod::setElevationM(double elevation)
+{
+  _elevation = elevation;
+}
+
+inline
+double
+SGGeod::getElevationFt(void) const
+{
+  return _elevation*SG_METER_TO_FEET;
+}
+
+inline
+void
+SGGeod::setElevationFt(double elevation)
+{
+  _elevation = elevation*SG_FEET_TO_METER;
+}
+
+/// Output to an ostream
+template<typename char_type, typename traits_type>
+inline
+std::basic_ostream<char_type, traits_type>&
+operator<<(std::basic_ostream<char_type, traits_type>& s, const SGGeod& g)
+{
+  return s << "lon = " << g.getLongitudeDeg()
+           << "deg, lat = " << g.getLatitudeDeg()
+           << "deg, elev = " << g.getElevationM()
+           << "m";
+}
+
+#endif

simgear/math/SGGeodesy.cxx

@@ -0,0 +1,157 @@
+// Copyright (C) 2006  Mathias Froehlich - Mathias.Froehlich@web.de
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Library General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// This library 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
+// Library General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+
+#ifdef HAVE_CONFIG_H
+#  include <simgear_config.h>
+#endif
+
+#include <cmath>
+
+#include "SGMath.hxx"
+
+// These are hard numbers from the WGS84 standard.  DON'T MODIFY
+// unless you want to change the datum.
+#define _EQURAD 6378137.0
+#define _FLATTENING 298.257223563
+
+// These are derived quantities more useful to the code:
+#if 0
+#define _SQUASH (1 - 1/_FLATTENING)
+#define _STRETCH (1/_SQUASH)
+#define _POLRAD (EQURAD * _SQUASH)
+#else
+// High-precision versions of the above produced with an arbitrary
+// precision calculator (the compiler might lose a few bits in the FPU
+// operations).  These are specified to 81 bits of mantissa, which is
+// higher than any FPU known to me:
+#define _SQUASH    0.9966471893352525192801545
+#define _STRETCH   1.0033640898209764189003079
+#define _POLRAD    6356752.3142451794975639668
+#endif
+
+// The constants from the WGS84 standard
+const double SGGeodesy::EQURAD = _EQURAD;
+const double SGGeodesy::iFLATTENING = _FLATTENING;
+const double SGGeodesy::SQUASH = _SQUASH;
+const double SGGeodesy::STRETCH = _STRETCH;
+const double SGGeodesy::POLRAD = _POLRAD;
+
+// additional derived and precomputable ones
+// for the geodetic conversion algorithm
+
+#define E2 fabs(1 - _SQUASH*_SQUASH)
+static double a = _EQURAD;
+static double ra2 = 1/(_EQURAD*_EQURAD);
+static double e = sqrt(E2);
+static double e2 = E2;
+static double e4 = E2*E2;
+
+#undef _EQURAD
+#undef _FLATTENING
+#undef _SQUASH
+#undef _STRETCH
+#undef _POLRAD
+#undef E2
+
+void
+SGGeodesy::SGCartToGeod(const SGVec3<double>& cart, SGGeod& geod)
+{
+  // according to
+  // H. Vermeille,
+  // Direct transformation from geocentric to geodetic ccordinates,
+  // Journal of Geodesy (2002) 76:451-454
+  double X = cart(0);
+  double Y = cart(1);
+  double Z = cart(2);
+  double XXpYY = X*X+Y*Y;
+  double sqrtXXpYY = sqrt(XXpYY);
+  double p = XXpYY*ra2;
+  double q = Z*Z*(1-e2)*ra2;
+  double r = 1/6.0*(p+q-e4);
+  double s = e4*p*q/(4*r*r*r);
+  double t = pow(1+s+sqrt(s*(2+s)), 1/3.0);
+  double u = r*(1+t+1/t);
+  double v = sqrt(u*u+e4*q);
+  double w = e2*(u+v-q)/(2*v);
+  double k = sqrt(u+v+w*w)-w;
+  double D = k*sqrtXXpYY/(k+e2);
+  geod.setLongitudeRad(2*atan2(Y, X+sqrtXXpYY));
+  double sqrtDDpZZ = sqrt(D*D+Z*Z);
+  geod.setLatitudeRad(2*atan2(Z, D+sqrtDDpZZ));
+  geod.setElevationM((k+e2-1)*sqrtDDpZZ/k);
+}
+
+void
+SGGeodesy::SGGeodToCart(const SGGeod& geod, SGVec3<double>& cart)
+{
+  // according to
+  // H. Vermeille,
+  // Direct transformation from geocentric to geodetic ccordinates,
+  // Journal of Geodesy (2002) 76:451-454
+  double lambda = geod.getLongitudeRad();
+  double phi = geod.getLatitudeRad();
+  double h = geod.getElevationM();
+  double sphi = sin(phi);
+  double n = a/sqrt(1-e2*sphi*sphi);
+  double cphi = cos(phi);
+  double slambda = sin(lambda);
+  double clambda = cos(lambda);
+  cart(0) = (h+n)*cphi*clambda;
+  cart(1) = (h+n)*cphi*slambda;
+  cart(2) = (h+n-e2*n)*sphi;
+}
+
+double
+SGGeodesy::SGGeodToSeaLevelRadius(const SGGeod& geod)
+{
+  // this is just a simplified version of the SGGeodToCart function above,
+  // substitute h = 0, take the 2-norm of the cartesian vector and simplify
+  double phi = geod.getLatitudeRad();
+  double sphi = sin(phi);
+  double sphi2 = sphi*sphi;
+  return a*sqrt((1 + (e4 - 2*e2)*sphi2)/(1 - e2*sphi2));
+}
+
+void
+SGGeodesy::SGCartToGeoc(const SGVec3<double>& cart, SGGeoc& geoc)
+{
+  double minVal = SGLimits<double>::min();
+  if (fabs(cart(0)) < minVal && fabs(cart(1)) < minVal)
+    geoc.setLongitudeRad(0);
+  else
+    geoc.setLongitudeRad(atan2(cart(1), cart(0)));
+
+  double nxy = sqrt(cart(0)*cart(0) + cart(1)*cart(1));
+  if (fabs(nxy) < minVal && fabs(cart(2)) < minVal)
+    geoc.setLatitudeRad(0);
+  else
+    geoc.setLatitudeRad(atan2(cart(2), nxy));
+
+  geoc.setRadiusM(norm(cart));
+}
+
+void
+SGGeodesy::SGGeocToCart(const SGGeoc& geoc, SGVec3<double>& cart)
+{
+  double lat = geoc.getLatitudeRad();
+  double lon = geoc.getLongitudeRad();
+  double slat = sin(lat);
+  double clat = cos(lat);
+  double slon = sin(lon);
+  double clon = cos(lon);
+  cart = geoc.getRadiusM()*SGVec3<double>(clat*clon, clat*slon, slat);
+}

simgear/math/SGGeodesy.hxx

@@ -0,0 +1,56 @@
+// Copyright (C) 2006  Mathias Froehlich - Mathias.Froehlich@web.de
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Library General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// This library 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
+// Library General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+
+#ifndef SGGeodesy_H
+#define SGGeodesy_H
+
+class SGGeoc;
+class SGGeod;
+
+template<typename T>
+class SGVec3;
+
+class SGGeodesy {
+public:
+  // Hard numbers from the WGS84 standard.
+  static const double EQURAD;
+  static const double iFLATTENING;
+  static const double SQUASH;
+  static const double STRETCH;
+  static const double POLRAD;
+
+  /// Takes a cartesian coordinate data and returns the geodetic
+  /// coordinates.
+  static void SGCartToGeod(const SGVec3<double>& cart, SGGeod& geod);
+  
+  /// Takes a geodetic coordinate data and returns the cartesian
+  /// coordinates.
+  static void SGGeodToCart(const SGGeod& geod, SGVec3<double>& cart);
+  
+  /// Takes a geodetic coordinate data and returns the sea level radius.
+  static double SGGeodToSeaLevelRadius(const SGGeod& geod);
+
+  /// Takes a cartesian coordinate data and returns the geocentric
+  /// coordinates.
+  static void SGCartToGeoc(const SGVec3<double>& cart, SGGeoc& geoc);
+  
+  /// Takes a geocentric coordinate data and returns the cartesian
+  /// coordinates.
+  static void SGGeocToCart(const SGGeoc& geoc, SGVec3<double>& cart);
+};
+
+#endif

simgear/math/SGLimits.hxx

@@ -0,0 +1,32 @@
+// Copyright (C) 2006  Mathias Froehlich - Mathias.Froehlich@web.de
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Library General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// This library 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
+// Library General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+
+#ifndef SGLimits_H
+#define SGLimits_H
+
+#include <limits>
+
+/// Helper class for epsilon and so on
+/// This is the possible place to hook in for machines not
+/// providing numeric_limits ...
+template<typename T>
+class SGLimits : public std::numeric_limits<T> {};
+
+typedef SGLimits<float> SGLimitsf;
+typedef SGLimits<double> SGLimitsd;
+
+#endif

simgear/math/SGMath.hxx

@@ -0,0 +1,35 @@
+// Copyright (C) 2006  Mathias Froehlich - Mathias.Froehlich@web.de
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Library General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// This library 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
+// Library General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+
+#ifndef SGMath_H
+#define SGMath_H
+
+/// Just include them all
+
+#include <iosfwd>
+
+#include "SGLimits.hxx"
+#include "SGMisc.hxx"
+#include "SGGeodesy.hxx"
+#include "SGVec3.hxx"
+#include "SGVec4.hxx"
+#include "SGQuat.hxx"
+#include "SGMatrix.hxx"
+#include "SGGeoc.hxx"
+#include "SGGeod.hxx"
+
+#endif

simgear/math/SGMathTest.cxx

@@ -0,0 +1,360 @@
+// Copyright (C) 2006  Mathias Froehlich - Mathias.Froehlich@web.de
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Library General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// This library 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
+// Library General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+
+#ifdef HAVE_CONFIG_H
+#  include <simgear_config.h>
+#endif
+
+#include <cstdlib>
+#include <iostream>
+
+#include <plib/sg.h>
+
+#include "SGMath.hxx"
+
+template<typename T>
+bool
+Vec3Test(void)
+{
+  SGVec3<T> v1, v2, v3;
+
+  // Check if the equivalent function works
+  v1 = SGVec3<T>(1, 2, 3);
+  v2 = SGVec3<T>(3, 2, 1);
+  if (equivalent(v1, v2))
+    return false;
+
+  // Check the unary minus operator
+  v3 = SGVec3<T>(-1, -2, -3);
+  if (!equivalent(-v1, v3))
+    return false;
+
+  // Check the unary plus operator
+  v3 = SGVec3<T>(1, 2, 3);
+  if (!equivalent(+v1, v3))
+    return false;
+
+  // Check the addition operator
+  v3 = SGVec3<T>(4, 4, 4);
+  if (!equivalent(v1 + v2, v3))
+    return false;
+
+  // Check the subtraction operator
+  v3 = SGVec3<T>(-2, 0, 2);
+  if (!equivalent(v1 - v2, v3))
+    return false;
+
+  // Check the scaler multiplication operator
+  v3 = SGVec3<T>(2, 4, 6);
+  if (!equivalent(2*v1, v3))
+    return false;
+
+  // Check the dot product
+  if (fabs(dot(v1, v2) - 10) > 10*SGLimits<T>::epsilon())
+    return false;
+
+  // Check the cross product
+  v3 = SGVec3<T>(-4, 8, -4);
+  if (!equivalent(cross(v1, v2), v3))
+    return false;
+
+  // Check the euclidean length
+  if (fabs(14 - length(v1)*length(v1)) > 14*SGLimits<T>::epsilon())
+    return false;
+
+  return true;
+}
+
+template<typename T>
+bool
+QuatTest(void)
+{
+  const SGVec3<T> e1(1, 0, 0);
+  const SGVec3<T> e2(0, 1, 0);
+  const SGVec3<T> e3(0, 0, 1);
+  SGVec3<T> v1, v2;
+  SGQuat<T> q1, q2, q3, q4;
+  // Check a rotation around the x axis
+  q1 = SGQuat<T>::fromAngleAxis(SGMisc<T>::pi(), e1);
+  v1 = SGVec3<T>(1, 2, 3);
+  v2 = SGVec3<T>(1, -2, -3);
+  if (!equivalent(q1.transform(v1), v2))
+    return false;
+  
+  // Check a rotation around the x axis
+  q1 = SGQuat<T>::fromAngleAxis(0.5*SGMisc<T>::pi(), e1);
+  v2 = SGVec3<T>(1, 3, -2);
+  if (!equivalent(q1.transform(v1), v2))
+    return false;
+
+  // Check a rotation around the y axis
+  q1 = SGQuat<T>::fromAngleAxis(SGMisc<T>::pi(), e2);
+  v2 = SGVec3<T>(-1, 2, -3);
+  if (!equivalent(q1.transform(v1), v2))
+    return false;
+  
+  // Check a rotation around the y axis
+  q1 = SGQuat<T>::fromAngleAxis(0.5*SGMisc<T>::pi(), e2);
+  v2 = SGVec3<T>(-3, 2, 1);
+  if (!equivalent(q1.transform(v1), v2))
+    return false;
+
+  // Check a rotation around the z axis
+  q1 = SGQuat<T>::fromAngleAxis(SGMisc<T>::pi(), e3);
+  v2 = SGVec3<T>(-1, -2, 3);
+  if (!equivalent(q1.transform(v1), v2))
+    return false;
+
+  // Check a rotation around the z axis
+  q1 = SGQuat<T>::fromAngleAxis(0.5*SGMisc<T>::pi(), e3);
+  v2 = SGVec3<T>(2, -1, 3);
+  if (!equivalent(q1.transform(v1), v2))
+    return false;
+
+  // Now check some successive transforms
+  // We can reuse the prevously tested stuff
+  q1 = SGQuat<T>::fromAngleAxis(0.5*SGMisc<T>::pi(), e1);
+  q2 = SGQuat<T>::fromAngleAxis(0.5*SGMisc<T>::pi(), e2);
+  q3 = q1*q2;
+  v2 = q2.transform(q1.transform(v1));
+  if (!equivalent(q3.transform(v1), v2))
+    return false;
+
+  /// Test from Euler angles
+  float x = 0.2*SGMisc<T>::pi();
+  float y = 0.3*SGMisc<T>::pi();
+  float z = 0.4*SGMisc<T>::pi();
+  q1 = SGQuat<T>::fromAngleAxis(z, e3);
+  q2 = SGQuat<T>::fromAngleAxis(y, e2);
+  q3 = SGQuat<T>::fromAngleAxis(x, e1);
+  v2 = q3.transform(q2.transform(q1.transform(v1)));
+  q4 = SGQuat<T>::fromEulerRad(z, y, x);
+  if (!equivalent(q4.transform(v1), v2))
+    return false;
+
+  /// Test angle axis forward and back transform
+  q1 = SGQuat<T>::fromAngleAxis(0.2*SGMisc<T>::pi(), e1);
+  q2 = SGQuat<T>::fromAngleAxis(0.7*SGMisc<T>::pi(), e2);
+  q3 = q1*q2;
+  SGVec3<T> angleAxis;
+  q1.getAngleAxis(angleAxis);
+  q4 = SGQuat<T>::fromAngleAxis(angleAxis);
+  if (!equivalent(q1, q4))
+    return false;
+  q2.getAngleAxis(angleAxis);
+  q4 = SGQuat<T>::fromAngleAxis(angleAxis);
+  if (!equivalent(q2, q4))
+    return false;
+  q3.getAngleAxis(angleAxis);
+  q4 = SGQuat<T>::fromAngleAxis(angleAxis);
+  if (!equivalent(q3, q4))
+    return false;
+
+  return true;
+}
+
+template<typename T>
+bool
+MatrixTest(void)
+{
+  // Create some test matrix
+  SGVec3<T> v0(2, 7, 17);
+  SGQuat<T> q0 = SGQuat<T>::fromAngleAxis(SGMisc<T>::pi(), normalize(v0));
+  SGMatrix<T> m0(q0, v0);
+
+  // Check the tqo forms of the inverse for that kind of special matrix
+  SGMatrix<T> m1, m2;
+  invert(m1, m0);
+  m2 = transNeg(m0);
+  if (!equivalent(m1, m2))
+    return false;
+
+  // Check matrix multiplication and inversion
+  if (!equivalent(m0*m1, SGMatrix<T>::unit()))
+    return false;
+  if (!equivalent(m1*m0, SGMatrix<T>::unit()))
+    return false;
+  if (!equivalent(m0*m2, SGMatrix<T>::unit()))
+    return false;
+  if (!equivalent(m2*m0, SGMatrix<T>::unit()))
+    return false;
+  
+  return true;
+}
+
+bool
+GeodesyTest(void)
+{
+  // We know that the values are on the order of 1
+  double epsDeg = 10*SGLimits<double>::epsilon();
+  // For the altitude values we need to tolerate relative errors in the order
+  // of the radius
+  double epsM = 1e6*SGLimits<double>::epsilon();
+
+  SGVec3<double> cart0, cart1;
+  SGGeod geod0, geod1;
+  SGGeoc geoc0;
+
+  // create some geodetic position
+  geod0 = SGGeod::fromDegM(30, 20, 17);
+
+  // Test the conversion routines to cartesian coordinates
+  cart0 = geod0;
+  geod1 = cart0;
+  if (epsDeg < fabs(geod0.getLongitudeDeg() - geod1.getLongitudeDeg()) ||
+      epsDeg < fabs(geod0.getLatitudeDeg() - geod1.getLatitudeDeg()) ||
+      epsM < fabs(geod0.getElevationM() - geod1.getElevationM()))
+    return false;
+
+  // Test the conversion routines to radial coordinates
+  geoc0 = cart0;
+  cart1 = geoc0;
+  if (!equivalent(cart0, cart1))
+    return false;
+
+  return true;
+}
+
+
+bool
+sgInterfaceTest(void)
+{
+  SGVec3f v3f = SGVec3f::e2();
+  SGVec4f v4f = SGVec4f::e2();
+  SGQuatf qf = SGQuatf::fromEulerRad(1.2, 1.3, -0.4);
+  SGMatrixf mf(qf, v3f);
+
+  // Copy to and from plibs types check if result is equal,
+  // test for exact equality
+  SGVec3f tv3f;
+  sgVec3 sv3f;
+  sgCopyVec3(sv3f, v3f.sg());
+  sgCopyVec3(tv3f.sg(), sv3f);
+  if (tv3f != v3f)
+    return false;
+  
+  // Copy to and from plibs types check if result is equal,
+  // test for exact equality
+  SGVec4f tv4f;
+  sgVec4 sv4f;
+  sgCopyVec4(sv4f, v4f.sg());
+  sgCopyVec4(tv4f.sg(), sv4f);
+  if (tv4f != v4f)
+    return false;
+
+  // Copy to and from plibs types check if result is equal,
+  // test for exact equality
+  SGQuatf tqf;
+  sgQuat sqf;
+  sgCopyQuat(sqf, qf.sg());
+  sgCopyQuat(tqf.sg(), sqf);
+  if (tqf != qf)
+    return false;
+
+  // Copy to and from plibs types check if result is equal,
+  // test for exact equality
+  SGMatrixf tmf;
+  sgMat4 smf;
+  sgCopyMat4(smf, mf.sg());
+  sgCopyMat4(tmf.sg(), smf);
+  if (tmf != mf)
+    return false;
+
+  return true;
+}
+
+bool
+sgdInterfaceTest(void)
+{
+  SGVec3d v3d = SGVec3d::e2();
+  SGVec4d v4d = SGVec4d::e2();
+  SGQuatd qd = SGQuatd::fromEulerRad(1.2, 1.3, -0.4);
+  SGMatrixd md(qd, v3d);
+
+  // Copy to and from plibs types check if result is equal,
+  // test for exact equality
+  SGVec3d tv3d;
+  sgdVec3 sv3d;
+  sgdCopyVec3(sv3d, v3d.sg());
+  sgdCopyVec3(tv3d.sg(), sv3d);
+  if (tv3d != v3d)
+    return false;
+  
+  // Copy to and from plibs types check if result is equal,
+  // test for exact equality
+  SGVec4d tv4d;
+  sgdVec4 sv4d;
+  sgdCopyVec4(sv4d, v4d.sg());
+  sgdCopyVec4(tv4d.sg(), sv4d);
+  if (tv4d != v4d)
+    return false;
+
+  // Copy to and from plibs types check if result is equal,
+  // test for exact equality
+  SGQuatd tqd;
+  sgdQuat sqd;
+  sgdCopyQuat(sqd, qd.sg());
+  sgdCopyQuat(tqd.sg(), sqd);
+  if (tqd != qd)
+    return false;
+
+  // Copy to and from plibs types check if result is equal,
+  // test for exact equality
+  SGMatrixd tmd;
+  sgdMat4 smd;
+  sgdCopyMat4(smd, md.sg());
+  sgdCopyMat4(tmd.sg(), smd);
+  if (tmd != md)
+    return false;
+
+  return true;
+}
+
+int
+main(void)
+{
+  // Do vector tests
+  if (!Vec3Test<float>())
+    return EXIT_FAILURE;
+  if (!Vec3Test<double>())
+    return EXIT_FAILURE;
+
+  // Do quaternion tests
+  if (!QuatTest<float>())
+    return EXIT_FAILURE;
+  if (!QuatTest<double>())
+    return EXIT_FAILURE;
+
+  // Do matrix tests
+  if (!MatrixTest<float>())
+    return EXIT_FAILURE;
+  if (!MatrixTest<double>())
+    return EXIT_FAILURE;
+
+  // Check geodetic/geocentric/cartesian conversions
+//   if (!GeodesyTest())
+//     return EXIT_FAILURE;
+
+  // Check interaction with sg*/sgd*
+  if (!sgInterfaceTest())
+    return EXIT_FAILURE;
+  if (!sgdInterfaceTest())
+    return EXIT_FAILURE;
+  
+  std::cout << "Successfully passed all tests!" << std::endl;
+  return EXIT_SUCCESS;
+}

simgear/math/SGMatrix.hxx

@@ -0,0 +1,598 @@
+// Copyright (C) 2006  Mathias Froehlich - Mathias.Froehlich@web.de
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Library General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// This library 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
+// Library General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+
+#ifndef SGMatrix_H
+#define SGMatrix_H
+
+/// Expression templates for poor programmers ... :)
+template<typename T>
+struct TransNegRef;
+
+template<typename T>
+class SGMatrix;
+
+/// 3D Matrix Class
+template<typename T>
+class SGMatrix {
+public:
+  enum { nCols = 4, nRows = 4, nEnts = 16 };
+  typedef T value_type;
+
+  /// Default constructor. Does not initialize at all.
+  /// If you need them zero initialized, use SGMatrix::zeros()
+  SGMatrix(void)
+  {
+    /// Initialize with nans in the debug build, that will guarantee to have
+    /// a fast uninitialized default constructor in the release but shows up
+    /// uninitialized values in the debug build very fast ...
+#ifndef NDEBUG
+    for (unsigned i = 0; i < nEnts; ++i)
+      _data.flat[i] = SGLimits<T>::quiet_NaN();
+#endif
+  }
+  /// Constructor. Initialize by the content of a plain array,
+  /// make sure it has at least 16 elements
+  explicit SGMatrix(const T* data)
+  { for (unsigned i = 0; i < nEnts; ++i) _data.flat[i] = data[i]; }
+
+  /// Constructor, build up a SGMatrix from given elements
+  SGMatrix(T m00, T m01, T m02, T m03,
+           T m10, T m11, T m12, T m13,
+           T m20, T m21, T m22, T m23,
+           T m30, T m31, T m32, T m33)
+  {
+    _data.flat[0] = m00; _data.flat[1] = m10;
+    _data.flat[2] = m20; _data.flat[3] = m30;
+    _data.flat[4] = m01; _data.flat[5] = m11;
+    _data.flat[6] = m21; _data.flat[7] = m31;
+    _data.flat[8] = m02; _data.flat[9] = m12;
+    _data.flat[10] = m22; _data.flat[11] = m32;
+    _data.flat[12] = m03; _data.flat[13] = m13;
+    _data.flat[14] = m23; _data.flat[15] = m33;
+  }
+
+  /// Constructor, build up a SGMatrix from a translation
+  SGMatrix(const SGVec3<T>& trans)
+  { set(trans); }
+
+  /// Constructor, build up a SGMatrix from a rotation and a translation
+  SGMatrix(const SGQuat<T>& quat, const SGVec3<T>& trans)
+  { set(quat, trans); }
+  /// Constructor, build up a SGMatrix from a rotation and a translation
+  SGMatrix(const SGQuat<T>& quat)
+  { set(quat); }
+
+  /// Copy constructor for a transposed negated matrix
+  SGMatrix(const TransNegRef<T>& tm)
+  { set(tm); }
+
+  /// Set from a tranlation
+  void set(const SGVec3<T>& trans)
+  {
+    _data.flat[0] = 1; _data.flat[4] = 0;
+    _data.flat[8] = 0; _data.flat[12] = -trans(0);
+    _data.flat[1] = 0; _data.flat[5] = 1;
+    _data.flat[9] = 0; _data.flat[13] = -trans(1);
+    _data.flat[2] = 0; _data.flat[6] = 0;
+    _data.flat[10] = 1; _data.flat[14] = -trans(2);
+    _data.flat[3] = 0; _data.flat[7] = 0;
+    _data.flat[11] = 0; _data.flat[15] = 1;
+  }
+
+  /// Set from a scale/rotation and tranlation
+  void set(const SGQuat<T>& quat, const SGVec3<T>& trans)
+  {
+    T w = quat.w(); T x = quat.x(); T y = quat.y(); T z = quat.z();
+    T xx = x*x; T yy = y*y; T zz = z*z;
+    T wx = w*x; T wy = w*y; T wz = w*z;
+    T xy = x*y; T xz = x*z; T yz = y*z;
+    _data.flat[0] = 1-2*(yy+zz); _data.flat[1] = 2*(xy-wz);
+    _data.flat[2] = 2*(xz+wy); _data.flat[3] = 0;
+    _data.flat[4] = 2*(xy+wz); _data.flat[5] = 1-2*(xx+zz);
+    _data.flat[6] = 2*(yz-wx); _data.flat[7] = 0;
+    _data.flat[8] = 2*(xz-wy); _data.flat[9] = 2*(yz+wx);
+    _data.flat[10] = 1-2*(xx+yy); _data.flat[11] = 0;
+    // Well, this one is ugly here, as that xform method on the current
+    // object needs the above data to be already set ...
+    SGVec3<T> t = xformVec(trans);
+    _data.flat[12] = -t(0); _data.flat[13] = -t(1);
+    _data.flat[14] = -t(2); _data.flat[15] = 1;
+  }
+  /// Set from a scale/rotation and tranlation
+  void set(const SGQuat<T>& quat)
+  {
+    T w = quat.w(); T x = quat.x(); T y = quat.y(); T z = quat.z();
+    T xx = x*x; T yy = y*y; T zz = z*z;
+    T wx = w*x; T wy = w*y; T wz = w*z;
+    T xy = x*y; T xz = x*z; T yz = y*z;
+    _data.flat[0] = 1-2*(yy+zz); _data.flat[1] = 2*(xy-wz);
+    _data.flat[2] = 2*(xz+wy); _data.flat[3] = 0;
+    _data.flat[4] = 2*(xy+wz); _data.flat[5] = 1-2*(xx+zz);
+    _data.flat[6] = 2*(yz-wx); _data.flat[7] = 0;
+    _data.flat[8] = 2*(xz-wy); _data.flat[9] = 2*(yz+wx);
+    _data.flat[10] = 1-2*(xx+yy); _data.flat[11] = 0;
+    _data.flat[12] = 0; _data.flat[13] = 0;
+    _data.flat[14] = 0; _data.flat[15] = 1;
+  }
+
+  /// set from a transposed negated matrix
+  void set(const TransNegRef<T>& tm)
+  {
+    const SGMatrix& m = tm.m;
+    _data.flat[0] = m(0,0);
+    _data.flat[1] = m(0,1);
+    _data.flat[2] = m(0,2);
+    _data.flat[3] = m(3,0);
+
+    _data.flat[4] = m(1,0);
+    _data.flat[5] = m(1,1);
+    _data.flat[6] = m(1,2);
+    _data.flat[7] = m(3,1);
+
+    _data.flat[8] = m(2,0);
+    _data.flat[9] = m(2,1);
+    _data.flat[10] = m(2,2);
+    _data.flat[11] = m(3,2);
+
+    // Well, this one is ugly here, as that xform method on the current
+    // object needs the above data to be already set ...
+    SGVec3<T> t = xformVec(SGVec3<T>(m(0,3), m(1,3), m(2,3)));
+    _data.flat[12] = -t(0);
+    _data.flat[13] = -t(1);
+    _data.flat[14] = -t(2);
+    _data.flat[15] = m(3,3);
+  }
+
+  /// Access by index, the index is unchecked
+  const T& operator()(unsigned i, unsigned j) const
+  { return _data.flat[i + 4*j]; }
+  /// Access by index, the index is unchecked
+  T& operator()(unsigned i, unsigned j)
+  { return _data.flat[i + 4*j]; }
+
+  /// Access raw data by index, the index is unchecked
+  const T& operator[](unsigned i) const
+  { return _data.flat[i]; }
+  /// Access by index, the index is unchecked
+  T& operator[](unsigned i)
+  { return _data.flat[i]; }
+
+  /// Get the data pointer
+  const T* data(void) const
+  { return _data.flat; }
+  /// Get the data pointer
+  T* data(void)
+  { return _data.flat; }
+
+  /// Readonly interface function to ssg's sgMat4/sgdMat4
+  const T (&sg(void) const)[4][4]
+  { return _data.carray; }
+  /// Interface function to ssg's sgMat4/sgdMat4
+  T (&sg(void))[4][4]
+  { return _data.carray; }
+
+  /// Inplace addition
+  SGMatrix& operator+=(const SGMatrix& m)
+  { for (unsigned i = 0; i < nEnts; ++i) _data.flat[i] += m._data.flat[i]; return *this; }
+  /// Inplace subtraction
+  SGMatrix& operator-=(const SGMatrix& m)
+  { for (unsigned i = 0; i < nEnts; ++i) _data.flat[i] -= m._data.flat[i]; return *this; }
+  /// Inplace scalar multiplication
+  template<typename S>
+  SGMatrix& operator*=(S s)
+  { for (unsigned i = 0; i < nEnts; ++i) _data.flat[i] *= s; return *this; }
+  /// Inplace scalar multiplication by 1/s
+  template<typename S>
+  SGMatrix& operator/=(S s)
+  { return operator*=(1/T(s)); }
+  /// Inplace matrix multiplication, post multiply
+  SGMatrix& operator*=(const SGMatrix<T>& m2);
+
+  SGVec3<T> xformPt(const SGVec3<T>& pt) const
+  {
+    SGVec3<T> tpt;
+    tpt(0) = (*this)(0,3);
+    tpt(1) = (*this)(1,3);
+    tpt(2) = (*this)(2,3);
+    for (unsigned i = 0; i < SGMatrix<T>::nCols-1; ++i) {
+      T tmp = pt(i);
+      tpt(0) += tmp*(*this)(0,i);
+      tpt(1) += tmp*(*this)(1,i);
+      tpt(2) += tmp*(*this)(2,i);
+    }
+    return tpt;
+  }
+  SGVec3<T> xformVec(const SGVec3<T>& v) const
+  {
+    SGVec3<T> tv;
+    T tmp = v(0);
+    tv(0) = tmp*(*this)(0,0);
+    tv(1) = tmp*(*this)(1,0);
+    tv(2) = tmp*(*this)(2,0);
+    for (unsigned i = 1; i < SGMatrix<T>::nCols-1; ++i) {
+      T tmp = v(i);
+      tv(0) += tmp*(*this)(0,i);
+      tv(1) += tmp*(*this)(1,i);
+      tv(2) += tmp*(*this)(2,i);
+    }
+    return tv;
+  }
+
+  /// Return an all zero matrix
+  static SGMatrix zeros(void)
+  { return SGMatrix(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0); }
+
+  /// Return a unit matrix
+  static SGMatrix unit(void)
+  { return SGMatrix(1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1); }
+
+private:
+  /// Required to make that alias safe.
+  union Data {
+    T flat[16];
+    T carray[4][4];
+  };
+
+  /// The actual data, the matrix is stored in column major order,
+  /// that matches the storage format of OpenGL
+  Data _data;
+};
+
+/// Class to distinguish between a matrix and the matrix with a transposed
+/// rotational part and a negated translational part
+template<typename T>
+struct TransNegRef {
+  TransNegRef(const SGMatrix<T>& _m) : m(_m) {}
+  const SGMatrix<T>& m;
+};
+
+/// Unary +, do nothing ...
+template<typename T>
+inline
+const SGMatrix<T>&
+operator+(const SGMatrix<T>& m)
+{ return m; }
+
+/// Unary -, do nearly nothing
+template<typename T>
+inline
+SGMatrix<T>
+operator-(const SGMatrix<T>& m)
+{
+  SGMatrix<T> ret;
+  for (unsigned i = 0; i < SGMatrix<T>::nEnts; ++i)
+    ret[i] = -m[i];
+  return ret;
+}
+
+/// Binary +
+template<typename T>
+inline
+SGMatrix<T>
+operator+(const SGMatrix<T>& m1, const SGMatrix<T>& m2)
+{
+  SGMatrix<T> ret;
+  for (unsigned i = 0; i < SGMatrix<T>::nEnts; ++i)
+    ret[i] = m1[i] + m2[i];
+  return ret;
+}
+
+/// Binary -
+template<typename T>
+inline
+SGMatrix<T>
+operator-(const SGMatrix<T>& m1, const SGMatrix<T>& m2)
+{
+  SGMatrix<T> ret;
+  for (unsigned i = 0; i < SGMatrix<T>::nEnts; ++i)
+    ret[i] = m1[i] - m2[i];
+  return ret;
+}
+
+/// Scalar multiplication
+template<typename S, typename T>
+inline
+SGMatrix<T>
+operator*(S s, const SGMatrix<T>& m)
+{
+  SGMatrix<T> ret;
+  for (unsigned i = 0; i < SGMatrix<T>::nEnts; ++i)
+    ret[i] = s*m[i];
+  return ret;
+}
+
+/// Scalar multiplication
+template<typename S, typename T>
+inline
+SGMatrix<T>
+operator*(const SGMatrix<T>& m, S s)
+{
+  SGMatrix<T> ret;
+  for (unsigned i = 0; i < SGMatrix<T>::nEnts; ++i)
+    ret[i] = s*m[i];
+  return ret;
+}
+
+/// Vector multiplication
+template<typename T>
+inline
+SGVec4<T>
+operator*(const SGMatrix<T>& m, const SGVec4<T>& v)
+{
+  SGVec4<T> mv;
+  T tmp = v(0);
+  mv(0) = tmp*m(0,0);
+  mv(1) = tmp*m(1,0);
+  mv(2) = tmp*m(2,0);
+  mv(3) = tmp*m(3,0);
+  for (unsigned i = 1; i < SGMatrix<T>::nCols; ++i) {
+    T tmp = v(i);
+    mv(0) += tmp*m(0,i);
+    mv(1) += tmp*m(1,i);
+    mv(2) += tmp*m(2,i);
+    mv(3) += tmp*m(3,i);
+  }
+  return mv;
+}
+
+/// Vector multiplication
+template<typename T>
+inline
+SGVec4<T>
+operator*(const TransNegRef<T>& tm, const SGVec4<T>& v)
+{
+  const SGMatrix<T>& m = tm.m;
+  SGVec4<T> mv;
+  SGVec3<T> v2;
+  T tmp = v(3);
+  mv(0) = v2(0) = -tmp*m(0,3);
+  mv(1) = v2(1) = -tmp*m(1,3);
+  mv(2) = v2(2) = -tmp*m(2,3);
+  mv(3) = tmp*m(3,3);
+  for (unsigned i = 0; i < SGMatrix<T>::nCols - 1; ++i) {
+    T tmp = v(i) + v2(i);
+    mv(0) += tmp*m(i,0);
+    mv(1) += tmp*m(i,1);
+    mv(2) += tmp*m(i,2);
+    mv(3) += tmp*m(3,i);
+  }
+  return mv;
+}
+
+/// Matrix multiplication
+template<typename T>
+inline
+SGMatrix<T>
+operator*(const SGMatrix<T>& m1, const SGMatrix<T>& m2)
+{
+  SGMatrix<T> m;
+  for (unsigned j = 0; j < SGMatrix<T>::nCols; ++j) {
+    T tmp = m2(0,j);
+    m(0,j) = tmp*m1(0,0);
+    m(1,j) = tmp*m1(1,0);
+    m(2,j) = tmp*m1(2,0);
+    m(3,j) = tmp*m1(3,0);
+    for (unsigned i = 1; i < SGMatrix<T>::nCols; ++i) {
+      T tmp = m2(i,j);
+      m(0,j) += tmp*m1(0,i);
+      m(1,j) += tmp*m1(1,i);
+      m(2,j) += tmp*m1(2,i);
+      m(3,j) += tmp*m1(3,i);
+    }
+  }
+  return m;
+}
+
+/// Inplace matrix multiplication, post multiply
+template<typename T>
+inline
+SGMatrix<T>&
+SGMatrix<T>::operator*=(const SGMatrix<T>& m2)
+{ (*this) = operator*(*this, m2); return *this; }
+
+/// Return a reference to the matrix typed to make sure we use the transposed
+/// negated matrix
+template<typename T>
+inline
+TransNegRef<T>
+transNeg(const SGMatrix<T>& m)
+{ return TransNegRef<T>(m); }
+
+/// Compute the inverse if the matrix src. Store the result in dst.
+/// Return if the matrix is nonsingular. If it is singular dst contains
+/// undefined values
+template<typename T>
+inline
+bool
+invert(SGMatrix<T>& dst, const SGMatrix<T>& src)
+{
+  // Do a LU decomposition with row pivoting and solve into dst
+  SGMatrix<T> tmp = src;
+  dst = SGMatrix<T>::unit();
+
+  for (unsigned i = 0; i < 4; ++i) {
+    T val = tmp(i,i);
+    unsigned ind = i;
+
+    // Find the row with the maximum value in the i-th colum
+    for (unsigned j = i + 1; j < 4; ++j) {
+      if (fabs(tmp(j, i)) > fabs(val)) {
+        ind = j;
+        val = tmp(j, i);
+      }
+    }
+
+    // Do row pivoting
+    if (ind != i) {
+      for (unsigned j = 0; j < 4; ++j) {
+        T t;
+        t = dst(i,j); dst(i,j) = dst(ind,j); dst(ind,j) = t;
+        t = tmp(i,j); tmp(i,j) = tmp(ind,j); tmp(ind,j) = t;
+      }
+    }
+
+    // Check for singularity
+    if (fabs(val) <= SGLimits<T>::min())
+      return false;
+
+    T ival = 1/val;
+    for (unsigned j = 0; j < 4; ++j) {
+      tmp(i,j) *= ival;
+      dst(i,j) *= ival;
+    }
+
+    for (unsigned j = 0; j < 4; ++j) {
+      if (j == i)
+        continue;
+
+      val = tmp(j,i);
+      for (unsigned k = 0; k < 4; ++k) {
+        tmp(j,k) -= tmp(i,k) * val;
+        dst(j,k) -= dst(i,k) * val;
+      }
+    }
+  }
+  return true;
+}
+
+/// The 1-norm of the matrix, this is the largest column sum of
+/// the absolute values of A.
+template<typename T>
+inline
+T
+norm1(const SGMatrix<T>& m)
+{
+  T nrm = 0;
+  for (unsigned i = 0; i < SGMatrix<T>::nRows; ++i) {
+    T sum = fabs(m(i, 0)) + fabs(m(i, 1)) + fabs(m(i, 2)) + fabs(m(i, 3));
+    if (nrm < sum)
+      nrm = sum;
+  }
+  return nrm;
+}
+
+/// The inf-norm of the matrix, this is the largest row sum of
+/// the absolute values of A.
+template<typename T>
+inline
+T
+normInf(const SGMatrix<T>& m)
+{
+  T nrm = 0;
+  for (unsigned i = 0; i < SGMatrix<T>::nCols; ++i) {
+    T sum = fabs(m(0, i)) + fabs(m(1, i)) + fabs(m(2, i)) + fabs(m(3, i));
+    if (nrm < sum)
+      nrm = sum;
+  }
+  return nrm;
+}
+
+/// Return true if exactly the same
+template<typename T>
+inline
+bool
+operator==(const SGMatrix<T>& m1, const SGMatrix<T>& m2)
+{
+  for (unsigned i = 0; i < SGMatrix<T>::nEnts; ++i)
+    if (m1[i] != m2[i])
+      return false;
+  return true;
+}
+
+/// Return true if not exactly the same
+template<typename T>
+inline
+bool
+operator!=(const SGMatrix<T>& m1, const SGMatrix<T>& m2)
+{ return ! (m1 == m2); }
+
+/// Return true if equal to the relative tolerance tol
+template<typename T>
+inline
+bool
+equivalent(const SGMatrix<T>& m1, const SGMatrix<T>& m2, T rtol, T atol)
+{ return norm1(m1 - m2) < rtol*(norm1(m1) + norm1(m2)) + atol; }
+
+/// Return true if equal to the relative tolerance tol
+template<typename T>
+inline
+bool
+equivalent(const SGMatrix<T>& m1, const SGMatrix<T>& m2, T rtol)
+{ return norm1(m1 - m2) < rtol*(norm1(m1) + norm1(m2)); }
+
+/// Return true if about equal to roundoff of the underlying type
+template<typename T>
+inline
+bool
+equivalent(const SGMatrix<T>& m1, const SGMatrix<T>& m2)
+{
+  T tol = 100*SGLimits<T>::epsilon();
+  return equivalent(m1, m2, tol, tol);
+}
+
+#ifndef NDEBUG
+template<typename T>
+inline
+bool
+isNaN(const SGMatrix<T>& m)
+{
+  for (unsigned i = 0; i < SGMatrix<T>::nEnts; ++i) {
+    if (SGMisc<T>::isNaN(m[i]))
+      return true;
+  }
+  return false;
+}
+#endif
+
+/// Output to an ostream
+template<typename char_type, typename traits_type, typename T> 
+inline
+std::basic_ostream<char_type, traits_type>&
+operator<<(std::basic_ostream<char_type, traits_type>& s, const SGMatrix<T>& m)
+{
+  s << "[ " << m(0,0) << ", " << m(0,1) << ", " << m(0,2) << ", " << m(0,3) << "\n";
+  s << "  " << m(1,0) << ", " << m(1,1) << ", " << m(1,2) << ", " << m(1,3) << "\n";
+  s << "  " << m(2,0) << ", " << m(2,1) << ", " << m(2,2) << ", " << m(2,3) << "\n";
+  s << "  " << m(3,0) << ", " << m(3,1) << ", " << m(3,2) << ", " << m(3,3) << " ]";
+  return s;
+}
+
+/// Two classes doing actually the same on different types
+typedef SGMatrix<float> SGMatrixf;
+typedef SGMatrix<double> SGMatrixd;
+
+inline
+SGMatrixf
+toMatrixf(const SGMatrixd& m)
+{
+  return SGMatrixf((float)m(0,0), (float)m(0,1), (float)m(0,2), (float)m(0,3),
+                   (float)m(1,0), (float)m(1,1), (float)m(1,2), (float)m(1,3),
+                   (float)m(3,0), (float)m(2,1), (float)m(2,2), (float)m(2,3),
+                   (float)m(4,0), (float)m(4,1), (float)m(4,2), (float)m(4,3));
+}
+
+inline
+SGMatrixd
+toMatrixd(const SGMatrixf& m)
+{
+  return SGMatrixd(m(0,0), m(0,1), m(0,2), m(0,3),
+                   m(1,0), m(1,1), m(1,2), m(1,3),
+                   m(3,0), m(2,1), m(2,2), m(2,3),
+                   m(4,0), m(4,1), m(4,2), m(4,3));
+}
+
+#endif

simgear/math/SGMisc.hxx

@@ -0,0 +1,77 @@
+// Copyright (C) 2006  Mathias Froehlich - Mathias.Froehlich@web.de
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Library General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// This library 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
+// Library General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+
+#ifndef SGMisc_H
+#define SGMisc_H
+
+#include <cmath>
+
+template<typename T>
+class SGMisc {
+public:
+  static T pi() { return T(3.1415926535897932384626433832795029L); }
+  static T min(const T& a, const T& b)
+  { return a < b ? a : b; }
+  static T min(const T& a, const T& b, const T& c)
+  { return min(min(a, b), c); }
+  static T min(const T& a, const T& b, const T& c, const T& d)
+  { return min(min(min(a, b), c), d); }
+  static T max(const T& a, const T& b)
+  { return a > b ? a : b; }
+  static T max(const T& a, const T& b, const T& c)
+  { return max(max(a, b), c); }
+  static T max(const T& a, const T& b, const T& c, const T& d)
+  { return max(max(max(a, b), c), d); }
+  static int sign(const T& a)
+  {
+    if (a < -SGLimits<T>::min())
+      return -1;
+    else if (SGLimits<T>::min() < a)
+      return 1;
+    else
+      return 0;
+  }
+
+  static T rad2deg(const T& val)
+  { return val*180/pi(); }
+  static T deg2rad(const T& val)
+  { return val*pi()/180; }
+
+#ifndef NDEBUG
+  /// Returns true if v is a NaN value
+  /// Use with care: allways code that you do not need to use that!
+  static bool isNaN(const T& v)
+  {
+#ifdef HAVE_ISNAN
+    return isnan(v);
+#elif defined HAVE_STD_ISNAN
+    return std::isnan(v);
+#else
+    // Use that every compare involving a NaN returns false
+    // But be careful, some usual compiler switches like for example
+    // -fast-math from gcc might optimize that expression to v != v which
+    // behaves exactly like the opposite ...
+    return !(v == v);
+#endif
+  }
+#endif
+};
+
+typedef SGMisc<float> SGMiscf;
+typedef SGMisc<double> SGMiscd;
+
+#endif

simgear/math/SGQuat.hxx

@@ -0,0 +1,569 @@
+// Copyright (C) 2006  Mathias Froehlich - Mathias.Froehlich@web.de
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Library General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// This library 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
+// Library General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+
+#ifndef SGQuat_H
+#define SGQuat_H
+
+/// 3D Vector Class
+template<typename T>
+class SGQuat {
+public:
+  typedef T value_type;
+
+  /// Default constructor. Does not initialize at all.
+  /// If you need them zero initialized, SGQuat::zeros()
+  SGQuat(void)
+  {
+    /// Initialize with nans in the debug build, that will guarantee to have
+    /// a fast uninitialized default constructor in the release but shows up
+    /// uninitialized values in the debug build very fast ...
+#ifndef NDEBUG
+    for (unsigned i = 0; i < 4; ++i)
+      _data[i] = SGLimits<T>::quiet_NaN();
+#endif
+  }
+  /// Constructor. Initialize by the given values
+  SGQuat(T _x, T _y, T _z, T _w)
+  { x() = _x; y() = _y; z() = _z; w() = _w; }
+  /// Constructor. Initialize by the content of a plain array,
+  /// make sure it has at least 4 elements
+  explicit SGQuat(const T* d)
+  { _data[0] = d[0]; _data[1] = d[1]; _data[2] = d[2]; _data[3] = d[3]; }
+
+  /// Return a unit quaternion
+  static SGQuat unit(void)
+  { return fromRealImag(1, SGVec3<T>(0)); }
+
+  /// Return a quaternion from euler angles
+  static SGQuat fromEulerRad(T z, T y, T x)
+  {
+    SGQuat q;
+    T zd2 = T(0.5)*z; T yd2 = T(0.5)*y; T xd2 = T(0.5)*x;
+    T Szd2 = sin(zd2); T Syd2 = sin(yd2); T Sxd2 = sin(xd2);
+    T Czd2 = cos(zd2); T Cyd2 = cos(yd2); T Cxd2 = cos(xd2);
+    T Cxd2Czd2 = Cxd2*Czd2; T Cxd2Szd2 = Cxd2*Szd2;
+    T Sxd2Szd2 = Sxd2*Szd2; T Sxd2Czd2 = Sxd2*Czd2;
+    q.w() = Cxd2Czd2*Cyd2 + Sxd2Szd2*Syd2;
+    q.x() = Sxd2Czd2*Cyd2 - Cxd2Szd2*Syd2;
+    q.y() = Cxd2Czd2*Syd2 + Sxd2Szd2*Cyd2;
+    q.z() = Cxd2Szd2*Cyd2 - Sxd2Czd2*Syd2;
+    return q;
+  }
+
+  /// Return a quaternion from euler angles in degrees
+  static SGQuat fromEulerDeg(T z, T y, T x)
+  {
+    return fromEulerRad(SGMisc<T>::deg2rad(z), SGMisc<T>::deg2rad(y),
+                        SGMisc<T>::deg2rad(x));
+  }
+
+  /// Return a quaternion from euler angles
+  static SGQuat fromYawPitchRoll(T y, T p, T r)
+  { return fromEulerRad(y, p, r); }
+
+  /// Return a quaternion from euler angles
+  static SGQuat fromYawPitchRollDeg(T y, T p, T r)
+  { return fromEulerDeg(y, p, r); }
+
+  /// Return a quaternion from euler angles
+  static SGQuat fromHeadAttBank(T h, T a, T b)
+  { return fromEulerRad(h, a, b); }
+
+  /// Return a quaternion from euler angles
+  static SGQuat fromHeadAttBankDeg(T h, T a, T b)
+  { return fromEulerDeg(h, a, b); }
+
+  /// Return a quaternion rotation the the horizontal local frame from given
+  /// longitude and latitude
+  static SGQuat fromLonLatRad(T lon, T lat)
+  {
+    SGQuat q;
+    T zd2 = T(0.5)*lon;
+    T yd2 = T(-0.25)*SGMisc<value_type>::pi() - T(0.5)*lat;
+    T Szd2 = sin(zd2);
+    T Syd2 = sin(yd2);
+    T Czd2 = cos(zd2);
+    T Cyd2 = cos(yd2);
+    q.w() = Czd2*Cyd2;
+    q.x() = -Szd2*Syd2;
+    q.y() = Czd2*Syd2;
+    q.z() = Szd2*Cyd2;
+    return q;
+  }
+
+  /// Return a quaternion rotation the the horizontal local frame from given
+  /// longitude and latitude
+  static SGQuat fromLonLatDeg(T lon, T lat)
+  { return fromLonLatRad(SGMisc<T>::deg2rad(lon), SGMisc<T>::deg2rad(lat)); }
+
+  /// Create a quaternion from the angle axis representation
+  static SGQuat fromAngleAxis(T angle, const SGVec3<T>& axis)
+  {
+    T angle2 = 0.5*angle;
+    return fromRealImag(cos(angle2), T(sin(angle2))*axis);
+  }
+
+  /// Create a quaternion from the angle axis representation
+  static SGQuat fromAngleAxisDeg(T angle, const SGVec3<T>& axis)
+  { return fromAngleAxis(SGMisc<T>::deg2rad(angle), axis); }
+
+  /// Create a quaternion from the angle axis representation where the angle
+  /// is stored in the axis' length
+  static SGQuat fromAngleAxis(const SGVec3<T>& axis)
+  {
+    T nAxis = norm(axis);
+    if (nAxis <= SGLimits<T>::min())
+      return SGQuat(1, 0, 0, 0);
+    T angle2 = 0.5*nAxis;
+    return fromRealImag(cos(angle2), T(sin(angle2)/nAxis)*axis);
+  }
+
+  /// Return a quaternion from real and imaginary part
+  static SGQuat fromRealImag(T r, const SGVec3<T>& i)
+  {
+    SGQuat q;
+    q.w() = r;
+    q.x() = i(0);
+    q.y() = i(1);
+    q.z() = i(2);
+    return q;
+  }
+
+  /// Return an all zero vector
+  static SGQuat zeros(void)
+  { return SGQuat(0, 0, 0, 0); }
+
+  /// write the euler angles into the references
+  void getEulerRad(T& zRad, T& yRad, T& xRad) const
+  {
+    value_type sqrQW = w()*w();
+    value_type sqrQX = x()*x();
+    value_type sqrQY = y()*y();
+    value_type sqrQZ = z()*z();
+
+    value_type num = 2*(y()*z() + w()*x());
+    value_type den = sqrQW - sqrQX - sqrQY + sqrQZ;
+    if (fabs(den) < SGLimits<value_type>::min() &&
+        fabs(num) < SGLimits<value_type>::min())
+      xRad = 0;
+    else
+      xRad = atan2(num, den);
+    
+    value_type tmp = 2*(x()*z() - w()*y());
+    if (tmp < -1)
+      yRad = 0.5*SGMisc<value_type>::pi();
+    else if (1 < tmp)
+      yRad = -0.5*SGMisc<value_type>::pi();
+    else
+      yRad = -asin(tmp);
+   
+    num = 2*(x()*y() + w()*z()); 
+    den = sqrQW + sqrQX - sqrQY - sqrQZ;
+    if (fabs(den) < SGLimits<value_type>::min() &&
+        fabs(num) < SGLimits<value_type>::min())
+      zRad = 0;
+    else {
+      value_type psi = atan2(num, den);
+      if (psi < 0)
+        psi += 2*SGMisc<value_type>::pi();
+      zRad = psi;
+    }
+  }
+
+  /// write the euler angles in degrees into the references
+  void getEulerDeg(T& zDeg, T& yDeg, T& xDeg) const
+  {
+    getEulerRad(zDeg, yDeg, xDeg);
+    zDeg = SGMisc<T>::rad2deg(zDeg);
+    yDeg = SGMisc<T>::rad2deg(yDeg);
+    xDeg = SGMisc<T>::rad2deg(xDeg);
+  }
+
+  /// write the angle axis representation into the references
+  void getAngleAxis(T& angle, SGVec3<T>& axis) const
+  {
+    T nrm = norm(*this);
+    if (nrm < SGLimits<T>::min()) {
+      angle = 0;
+      axis = SGVec3<T>(0, 0, 0);
+    } else {
+      T rNrm = 1/nrm;
+      angle = acos(SGMisc<T>::max(-1, SGMisc<T>::min(1, rNrm*w())));
+      T sAng = sin(angle);
+      if (fabs(sAng) < SGLimits<T>::min())
+        axis = SGVec3<T>(1, 0, 0);
+      else 
+        axis = (rNrm/sAng)*imag(*this);
+      angle *= 2;
+    }
+  }
+
+  /// write the angle axis representation into the references
+  void getAngleAxis(SGVec3<T>& axis) const
+  {
+    T angle;
+    getAngleAxis(angle, axis);
+    axis *= angle;
+  }
+
+  /// Access by index, the index is unchecked
+  const T& operator()(unsigned i) const
+  { return _data[i]; }
+  /// Access by index, the index is unchecked
+  T& operator()(unsigned i)
+  { return _data[i]; }
+
+  /// Access raw data by index, the index is unchecked
+  const T& operator[](unsigned i) const
+  { return _data[i]; }
+  /// Access raw data by index, the index is unchecked
+  T& operator[](unsigned i)
+  { return _data[i]; }
+
+  /// Access the x component
+  const T& x(void) const
+  { return _data[0]; }
+  /// Access the x component
+  T& x(void)
+  { return _data[0]; }
+  /// Access the y component
+  const T& y(void) const
+  { return _data[1]; }
+  /// Access the y component
+  T& y(void)
+  { return _data[1]; }
+  /// Access the z component
+  const T& z(void) const
+  { return _data[2]; }
+  /// Access the z component
+  T& z(void)
+  { return _data[2]; }
+  /// Access the w component
+  const T& w(void) const
+  { return _data[3]; }
+  /// Access the w component
+  T& w(void)
+  { return _data[3]; }
+
+  /// Get the data pointer, usefull for interfacing with plib's sg*Vec
+  const T* data(void) const
+  { return _data; }
+  /// Get the data pointer, usefull for interfacing with plib's sg*Vec
+  T* data(void)
+  { return _data; }
+
+  /// Readonly interface function to ssg's sgQuat/sgdQuat
+  const T (&sg(void) const)[4]
+  { return _data; }
+  /// Interface function to ssg's sgQuat/sgdQuat
+  T (&sg(void))[4]
+  { return _data; }
+
+  /// Inplace addition
+  SGQuat& operator+=(const SGQuat& v)
+  { _data[0]+=v(0);_data[1]+=v(1);_data[2]+=v(2);_data[3]+=v(3);return *this; }
+  /// Inplace subtraction
+  SGQuat& operator-=(const SGQuat& v)
+  { _data[0]-=v(0);_data[1]-=v(1);_data[2]-=v(2);_data[3]-=v(3);return *this; }
+  /// Inplace scalar multiplication
+  template<typename S>
+  SGQuat& operator*=(S s)
+  { _data[0] *= s; _data[1] *= s; _data[2] *= s; _data[3] *= s; return *this; }
+  /// Inplace scalar multiplication by 1/s
+  template<typename S>
+  SGQuat& operator/=(S s)
+  { return operator*=(1/T(s)); }
+  /// Inplace quaternion multiplication
+  SGQuat& operator*=(const SGQuat& v);
+
+  /// Transform a vector from the current coordinate frame to a coordinate
+  /// frame rotated with the quaternion
+  SGVec3<T> transform(const SGVec3<T>& v) const
+  {
+    value_type r = 2/dot(*this, *this);
+    SGVec3<T> qimag = imag(*this);
+    value_type qr = real(*this);
+    return (r*qr*qr - 1)*v + (r*dot(qimag, v))*qimag - (r*qr)*cross(qimag, v);
+  }
+  /// Transform a vector from the coordinate frame rotated with the quaternion
+  /// to the current coordinate frame
+  SGVec3<T> backTransform(const SGVec3<T>& v) const
+  {
+    value_type r = 2/dot(*this, *this);
+    SGVec3<T> qimag = imag(*this);
+    value_type qr = real(*this);
+    return (r*qr*qr - 1)*v + (r*dot(qimag, v))*qimag + (r*qr)*cross(qimag, v);
+  }
+
+  /// Rotate a given vector with the quaternion
+  SGVec3<T> rotate(const SGVec3<T>& v) const
+  { return backTransform(v); }
+  /// Rotate a given vector with the inverse quaternion
+  SGVec3<T> rotateBack(const SGVec3<T>& v) const
+  { return transform(v); }
+
+  /// Return the time derivative of the quaternion given the angular velocity
+  SGQuat
+  derivative(const SGVec3<T>& angVel)
+  {
+    SGQuat deriv;
+
+    deriv.w() = 0.5*(-x()*angVel(0) - y()*angVel(1) - z()*angVel(2));
+    deriv.x() = 0.5*( w()*angVel(0) - z()*angVel(1) + y()*angVel(2));
+    deriv.y() = 0.5*( z()*angVel(0) + w()*angVel(1) - x()*angVel(2));
+    deriv.z() = 0.5*(-y()*angVel(0) + x()*angVel(1) + w()*angVel(2));
+    
+    return deriv;
+  }
+
+private:
+  /// The actual data
+  T _data[4];
+};
+
+/// Unary +, do nothing ...
+template<typename T>
+inline
+const SGQuat<T>&
+operator+(const SGQuat<T>& v)
+{ return v; }
+
+/// Unary -, do nearly nothing
+template<typename T>
+inline
+SGQuat<T>
+operator-(const SGQuat<T>& v)
+{ return SGQuat<T>(-v(0), -v(1), -v(2), -v(3)); }
+
+/// Binary +
+template<typename T>
+inline
+SGQuat<T>
+operator+(const SGQuat<T>& v1, const SGQuat<T>& v2)
+{ return SGQuat<T>(v1(0)+v2(0), v1(1)+v2(1), v1(2)+v2(2), v1(3)+v2(3)); }
+
+/// Binary -
+template<typename T>
+inline
+SGQuat<T>
+operator-(const SGQuat<T>& v1, const SGQuat<T>& v2)
+{ return SGQuat<T>(v1(0)-v2(0), v1(1)-v2(1), v1(2)-v2(2), v1(3)-v2(3)); }
+
+/// Scalar multiplication
+template<typename S, typename T>
+inline
+SGQuat<T>
+operator*(S s, const SGQuat<T>& v)
+{ return SGQuat<T>(s*v(0), s*v(1), s*v(2), s*v(3)); }
+
+/// Scalar multiplication
+template<typename S, typename T>
+inline
+SGQuat<T>
+operator*(const SGQuat<T>& v, S s)
+{ return SGQuat<T>(s*v(0), s*v(1), s*v(2), s*v(3)); }
+
+/// Quaterion multiplication
+template<typename T>
+inline
+SGQuat<T>
+operator*(const SGQuat<T>& v1, const SGQuat<T>& v2)
+{
+  SGQuat<T> v;
+  v.x() = v1.w()*v2.x() + v1.x()*v2.w() + v1.y()*v2.z() - v1.z()*v2.y();
+  v.y() = v1.w()*v2.y() - v1.x()*v2.z() + v1.y()*v2.w() + v1.z()*v2.x();
+  v.z() = v1.w()*v2.z() + v1.x()*v2.y() - v1.y()*v2.x() + v1.z()*v2.w();
+  v.w() = v1.w()*v2.w() - v1.x()*v2.x() - v1.y()*v2.y() - v1.z()*v2.z();
+  return v;
+}
+
+/// Now define the inplace multiplication
+template<typename T>
+inline
+SGQuat<T>&
+SGQuat<T>::operator*=(const SGQuat& v)
+{ (*this) = (*this)*v; return *this; }
+
+/// The conjugate of the quaternion, this is also the
+/// inverse for normalized quaternions
+template<typename T>
+inline
+SGQuat<T>
+conj(const SGQuat<T>& v)
+{ return SGQuat<T>(-v(0), -v(1), -v(2), v(3)); }
+
+/// The conjugate of the quaternion, this is also the
+/// inverse for normalized quaternions
+template<typename T>
+inline
+SGQuat<T>
+inverse(const SGQuat<T>& v)
+{ return (1/dot(v, v))*SGQuat<T>(-v(0), -v(1), -v(2), v(3)); }
+
+/// The imagniary part of the quaternion
+template<typename T>
+inline
+T
+real(const SGQuat<T>& v)
+{ return v.w(); }
+
+/// The imagniary part of the quaternion
+template<typename T>
+inline
+SGVec3<T>
+imag(const SGQuat<T>& v)
+{ return SGVec3<T>(v.x(), v.y(), v.z()); }
+
+/// Scalar dot product
+template<typename T>
+inline
+T
+dot(const SGQuat<T>& v1, const SGQuat<T>& v2)
+{ return v1(0)*v2(0) + v1(1)*v2(1) + v1(2)*v2(2) + v1(3)*v2(3); }
+
+/// The euclidean norm of the vector, that is what most people call length
+template<typename T>
+inline
+T
+norm(const SGQuat<T>& v)
+{ return sqrt(dot(v, v)); }
+
+/// The euclidean norm of the vector, that is what most people call length
+template<typename T>
+inline
+T
+length(const SGQuat<T>& v)
+{ return sqrt(dot(v, v)); }
+
+/// The 1-norm of the vector, this one is the fastest length function we
+/// can implement on modern cpu's
+template<typename T>
+inline
+T
+norm1(const SGQuat<T>& v)
+{ return fabs(v(0)) + fabs(v(1)) + fabs(v(2)) + fabs(v(3)); }
+
+/// The euclidean norm of the vector, that is what most people call length
+template<typename T>
+inline
+SGQuat<T>
+normalize(const SGQuat<T>& q)
+{ return (1/norm(q))*q; }
+
+/// Return true if exactly the same
+template<typename T>
+inline
+bool
+operator==(const SGQuat<T>& v1, const SGQuat<T>& v2)
+{ return v1(0)==v2(0) && v1(1)==v2(1) && v1(2)==v2(2) && v1(3)==v2(3); }
+
+/// Return true if not exactly the same
+template<typename T>
+inline
+bool
+operator!=(const SGQuat<T>& v1, const SGQuat<T>& v2)
+{ return ! (v1 == v2); }
+
+/// Return true if equal to the relative tolerance tol
+/// Note that this is not the same than comparing quaternions to represent
+/// the same rotation
+template<typename T>
+inline
+bool
+equivalent(const SGQuat<T>& v1, const SGQuat<T>& v2, T tol)
+{ return norm1(v1 - v2) < tol*(norm1(v1) + norm1(v2)); }
+
+/// Return true if about equal to roundoff of the underlying type
+/// Note that this is not the same than comparing quaternions to represent
+/// the same rotation
+template<typename T>
+inline
+bool
+equivalent(const SGQuat<T>& v1, const SGQuat<T>& v2)
+{ return equivalent(v1, v2, 100*SGLimits<T>::epsilon()); }
+
+#ifndef NDEBUG
+template<typename T>
+inline
+bool
+isNaN(const SGQuat<T>& v)
+{
+  return SGMisc<T>::isNaN(v(0)) || SGMisc<T>::isNaN(v(1))
+    || SGMisc<T>::isNaN(v(2)) || SGMisc<T>::isNaN(v(3));
+}
+#endif
+
+/// quaternion interpolation for t in [0,1] interpolate between src (=0)
+/// and dst (=1)
+template<typename T>
+inline
+SGQuat<T>
+interpolate(T t, const SGQuat<T>& src, const SGQuat<T>& dst)
+{
+  T cosPhi = dot(src, dst);
+  // need to take the shorter way ...
+  int signCosPhi = SGMisc<T>::sign(cosPhi);
+  // cosPhi must be corrected for that sign
+  cosPhi = fabs(cosPhi);
+
+  // first opportunity to fail - make sure acos will succeed later -
+  // result is correct
+  if (1 <= cosPhi)
+    return dst;
+
+  // now the half angle between the orientations
+  T o = acos(cosPhi);
+
+  // need the scales now, if the angle is very small, do linear interpolation
+  // to avoid instabilities
+  T scale0, scale1;
+  if (fabs(o) < SGLimits<T>::epsilon()) {
+    scale0 = 1 - t;
+    scale1 = t;
+  } else {
+    // note that we can give a positive lower bound for sin(o) here
+    T sino = sin(o);
+    T so = 1/sino;
+    scale0 = sin((1 - t)*o)*so;
+    scale1 = sin(t*o)*so;
+  }
+
+  return scale0*src + signCosPhi*scale1*dst;
+}
+
+/// Output to an ostream
+template<typename char_type, typename traits_type, typename T>
+inline
+std::basic_ostream<char_type, traits_type>&
+operator<<(std::basic_ostream<char_type, traits_type>& s, const SGQuat<T>& v)
+{ return s << "[ " << v(0) << ", " << v(1) << ", " << v(2) << ", " << v(3) << " ]"; }
+
+/// Two classes doing actually the same on different types
+typedef SGQuat<float> SGQuatf;
+typedef SGQuat<double> SGQuatd;
+
+inline
+SGQuatf
+toQuatf(const SGQuatd& v)
+{ return SGQuatf((float)v(0), (float)v(1), (float)v(2), (float)v(3)); }
+
+inline
+SGQuatd
+toQuatd(const SGQuatf& v)
+{ return SGQuatd(v(0), v(1), v(2), v(3)); }
+
+#endif

simgear/math/SGVec3.hxx

@@ -0,0 +1,292 @@
+// Copyright (C) 2006  Mathias Froehlich - Mathias.Froehlich@web.de
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Library General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// This library 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
+// Library General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+
+#ifndef SGVec3_H
+#define SGVec3_H
+
+/// 3D Vector Class
+template<typename T>
+class SGVec3 {
+public:
+  typedef T value_type;
+
+  /// Default constructor. Does not initialize at all.
+  /// If you need them zero initialized, use SGVec3::zeros()
+  SGVec3(void)
+  {
+    /// Initialize with nans in the debug build, that will guarantee to have
+    /// a fast uninitialized default constructor in the release but shows up
+    /// uninitialized values in the debug build very fast ...
+#ifndef NDEBUG
+    for (unsigned i = 0; i < 3; ++i)
+      _data[i] = SGLimits<T>::quiet_NaN();
+#endif
+  }
+  /// Constructor. Initialize by the given values
+  SGVec3(T x, T y, T z)
+  { _data[0] = x; _data[1] = y; _data[2] = z; }
+  /// Constructor. Initialize by the content of a plain array,
+  /// make sure it has at least 3 elements
+  explicit SGVec3(const T* data)
+  { _data[0] = data[0]; _data[1] = data[1]; _data[2] = data[2]; }
+  /// Constructor. Initialize by a geodetic coordinate
+  /// Note that this conversion is relatively expensive to compute
+  SGVec3(const SGGeod& geod)
+  { SGGeodesy::SGGeodToCart(geod, *this); }
+  /// Constructor. Initialize by a geocentric coordinate
+  /// Note that this conversion is relatively expensive to compute
+  SGVec3(const SGGeoc& geoc)
+  { SGGeodesy::SGGeocToCart(geoc, *this); }
+
+  /// Access by index, the index is unchecked
+  const T& operator()(unsigned i) const
+  { return _data[i]; }
+  /// Access by index, the index is unchecked
+  T& operator()(unsigned i)
+  { return _data[i]; }
+
+  /// Access raw data by index, the index is unchecked
+  const T& operator[](unsigned i) const
+  { return _data[i]; }
+  /// Access raw data by index, the index is unchecked
+  T& operator[](unsigned i)
+  { return _data[i]; }
+
+  /// Access the x component
+  const T& x(void) const
+  { return _data[0]; }
+  /// Access the x component
+  T& x(void)
+  { return _data[0]; }
+  /// Access the y component
+  const T& y(void) const
+  { return _data[1]; }
+  /// Access the y component
+  T& y(void)
+  { return _data[1]; }
+  /// Access the z component
+  const T& z(void) const
+  { return _data[2]; }
+  /// Access the z component
+  T& z(void)
+  { return _data[2]; }
+
+  /// Get the data pointer
+  const T* data(void) const
+  { return _data; }
+  /// Get the data pointer
+  T* data(void)
+  { return _data; }
+
+  /// Readonly interface function to ssg's sgVec3/sgdVec3
+  const T (&sg(void) const)[3]
+  { return _data; }
+  /// Interface function to ssg's sgVec3/sgdVec3
+  T (&sg(void))[3]
+  { return _data; }
+
+  /// Inplace addition
+  SGVec3& operator+=(const SGVec3& v)
+  { _data[0] += v(0); _data[1] += v(1); _data[2] += v(2); return *this; }
+  /// Inplace subtraction
+  SGVec3& operator-=(const SGVec3& v)
+  { _data[0] -= v(0); _data[1] -= v(1); _data[2] -= v(2); return *this; }
+  /// Inplace scalar multiplication
+  template<typename S>
+  SGVec3& operator*=(S s)
+  { _data[0] *= s; _data[1] *= s; _data[2] *= s; return *this; }
+  /// Inplace scalar multiplication by 1/s
+  template<typename S>
+  SGVec3& operator/=(S s)
+  { return operator*=(1/T(s)); }
+
+  /// Return an all zero vector
+  static SGVec3 zeros(void)
+  { return SGVec3(0, 0, 0); }
+  /// Return unit vectors
+  static SGVec3 e1(void)
+  { return SGVec3(1, 0, 0); }
+  static SGVec3 e2(void)
+  { return SGVec3(0, 1, 0); }
+  static SGVec3 e3(void)
+  { return SGVec3(0, 0, 1); }
+
+private:
+  /// The actual data
+  T _data[3];
+};
+
+/// Unary +, do nothing ...
+template<typename T>
+inline
+const SGVec3<T>&
+operator+(const SGVec3<T>& v)
+{ return v; }
+
+/// Unary -, do nearly nothing
+template<typename T>
+inline
+SGVec3<T>
+operator-(const SGVec3<T>& v)
+{ return SGVec3<T>(-v(0), -v(1), -v(2)); }
+
+/// Binary +
+template<typename T>
+inline
+SGVec3<T>
+operator+(const SGVec3<T>& v1, const SGVec3<T>& v2)
+{ return SGVec3<T>(v1(0)+v2(0), v1(1)+v2(1), v1(2)+v2(2)); }
+
+/// Binary -
+template<typename T>
+inline
+SGVec3<T>
+operator-(const SGVec3<T>& v1, const SGVec3<T>& v2)
+{ return SGVec3<T>(v1(0)-v2(0), v1(1)-v2(1), v1(2)-v2(2)); }
+
+/// Scalar multiplication
+template<typename S, typename T>
+inline
+SGVec3<T>
+operator*(S s, const SGVec3<T>& v)
+{ return SGVec3<T>(s*v(0), s*v(1), s*v(2)); }
+
+/// Scalar multiplication
+template<typename S, typename T>
+inline
+SGVec3<T>
+operator*(const SGVec3<T>& v, S s)
+{ return SGVec3<T>(s*v(0), s*v(1), s*v(2)); }
+
+/// Scalar dot product
+template<typename T>
+inline
+T
+dot(const SGVec3<T>& v1, const SGVec3<T>& v2)
+{ return v1(0)*v2(0) + v1(1)*v2(1) + v1(2)*v2(2); }
+
+/// The euclidean norm of the vector, that is what most people call length
+template<typename T>
+inline
+T
+norm(const SGVec3<T>& v)
+{ return sqrt(dot(v, v)); }
+
+/// The euclidean norm of the vector, that is what most people call length
+template<typename T>
+inline
+T
+length(const SGVec3<T>& v)
+{ return sqrt(dot(v, v)); }
+
+/// The 1-norm of the vector, this one is the fastest length function we
+/// can implement on modern cpu's
+template<typename T>
+inline
+T
+norm1(const SGVec3<T>& v)
+{ return fabs(v(0)) + fabs(v(1)) + fabs(v(2)); }
+
+/// Vector cross product
+template<typename T>
+inline
+SGVec3<T>
+cross(const SGVec3<T>& v1, const SGVec3<T>& v2)
+{
+  return SGVec3<T>(v1(1)*v2(2) - v1(2)*v2(1),
+                   v1(2)*v2(0) - v1(0)*v2(2),
+                   v1(0)*v2(1) - v1(1)*v2(0));
+}
+
+/// The euclidean norm of the vector, that is what most people call length
+template<typename T>
+inline
+SGVec3<T>
+normalize(const SGVec3<T>& v)
+{ return (1/norm(v))*v; }
+
+/// Return true if exactly the same
+template<typename T>
+inline
+bool
+operator==(const SGVec3<T>& v1, const SGVec3<T>& v2)
+{ return v1(0) == v2(0) && v1(1) == v2(1) && v1(2) == v2(2); }
+
+/// Return true if not exactly the same
+template<typename T>
+inline
+bool
+operator!=(const SGVec3<T>& v1, const SGVec3<T>& v2)
+{ return ! (v1 == v2); }
+
+/// Return true if equal to the relative tolerance tol
+template<typename T>
+inline
+bool
+equivalent(const SGVec3<T>& v1, const SGVec3<T>& v2, T rtol, T atol)
+{ return norm1(v1 - v2) < rtol*(norm1(v1) + norm1(v2)) + atol; }
+
+/// Return true if equal to the relative tolerance tol
+template<typename T>
+inline
+bool
+equivalent(const SGVec3<T>& v1, const SGVec3<T>& v2, T rtol)
+{ return norm1(v1 - v2) < rtol*(norm1(v1) + norm1(v2)); }
+
+/// Return true if about equal to roundoff of the underlying type
+template<typename T>
+inline
+bool
+equivalent(const SGVec3<T>& v1, const SGVec3<T>& v2)
+{
+  T tol = 100*SGLimits<T>::epsilon();
+  return equivalent(v1, v2, tol, tol);
+}
+
+#ifndef NDEBUG
+template<typename T>
+inline
+bool
+isNaN(const SGVec3<T>& v)
+{
+  return SGMisc<T>::isNaN(v(0)) ||
+    SGMisc<T>::isNaN(v(1)) || SGMisc<T>::isNaN(v(2));
+}
+#endif
+
+/// Output to an ostream
+template<typename char_type, typename traits_type, typename T>
+inline
+std::basic_ostream<char_type, traits_type>&
+operator<<(std::basic_ostream<char_type, traits_type>& s, const SGVec3<T>& v)
+{ return s << "[ " << v(0) << ", " << v(1) << ", " << v(2) << " ]"; }
+
+/// Two classes doing actually the same on different types
+typedef SGVec3<float> SGVec3f;
+typedef SGVec3<double> SGVec3d;
+
+inline
+SGVec3f
+toVec3f(const SGVec3d& v)
+{ return SGVec3f((float)v(0), (float)v(1), (float)v(2)); }
+
+inline
+SGVec3d
+toVec3d(const SGVec3f& v)
+{ return SGVec3d(v(0), v(1), v(2)); }
+
+#endif

simgear/math/SGVec4.hxx

@@ -0,0 +1,283 @@
+// Copyright (C) 2006  Mathias Froehlich - Mathias.Froehlich@web.de
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Library General Public
+// License as published by the Free Software Foundation; either
+// version 2 of the License, or (at your option) any later version.
+//
+// This library 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
+// Library General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+//
+
+#ifndef SGVec4_H
+#define SGVec4_H
+
+/// 4D Vector Class
+template<typename T>
+class SGVec4 {
+public:
+  typedef T value_type;
+
+  /// Default constructor. Does not initialize at all.
+  /// If you need them zero initialized, use SGVec4::zeros()
+  SGVec4(void)
+  {
+    /// Initialize with nans in the debug build, that will guarantee to have
+    /// a fast uninitialized default constructor in the release but shows up
+    /// uninitialized values in the debug build very fast ...
+#ifndef NDEBUG
+    for (unsigned i = 0; i < 4; ++i)
+      _data[i] = SGLimits<T>::quiet_NaN();
+#endif
+  }
+  /// Constructor. Initialize by the given values
+  SGVec4(T x, T y, T z, T w)
+  { _data[0] = x; _data[1] = y; _data[2] = z; _data[3] = w; }
+  /// Constructor. Initialize by the content of a plain array,
+  /// make sure it has at least 3 elements
+  explicit SGVec4(const T* d)
+  { _data[0] = d[0]; _data[1] = d[1]; _data[2] = d[2]; _data[3] = d[3]; }
+
+
+  /// Access by index, the index is unchecked
+  const T& operator()(unsigned i) const
+  { return _data[i]; }
+  /// Access by index, the index is unchecked
+  T& operator()(unsigned i)
+  { return _data[i]; }
+
+  /// Access raw data by index, the index is unchecked
+  const T& operator[](unsigned i) const
+  { return _data[i]; }
+  /// Access raw data by index, the index is unchecked
+  T& operator[](unsigned i)
+  { return _data[i]; }
+
+  /// Access the x component
+  const T& x(void) const
+  { return _data[0]; }
+  /// Access the x component
+  T& x(void)
+  { return _data[0]; }
+  /// Access the y component
+  const T& y(void) const
+  { return _data[1]; }
+  /// Access the y component
+  T& y(void)
+  { return _data[1]; }
+  /// Access the z component
+  const T& z(void) const
+  { return _data[2]; }
+  /// Access the z component
+  T& z(void)
+  { return _data[2]; }
+  /// Access the x component
+  const T& w(void) const
+  { return _data[3]; }
+  /// Access the x component
+  T& w(void)
+  { return _data[3]; }
+
+
+  /// Get the data pointer, usefull for interfacing with plib's sg*Vec
+  const T* data(void) const
+  { return _data; }
+  /// Get the data pointer, usefull for interfacing with plib's sg*Vec
+  T* data(void)
+  { return _data; }
+
+  /// Readonly interface function to ssg's sgVec3/sgdVec3
+  const T (&sg(void) const)[4]
+  { return _data; }
+  /// Interface function to ssg's sgVec3/sgdVec3
+  T (&sg(void))[4]
+  { return _data; }
+
+  /// Inplace addition
+  SGVec4& operator+=(const SGVec4& v)
+  { _data[0]+=v(0);_data[1]+=v(1);_data[2]+=v(2);_data[3]+=v(3);return *this; }
+  /// Inplace subtraction
+  SGVec4& operator-=(const SGVec4& v)
+  { _data[0]-=v(0);_data[1]-=v(1);_data[2]-=v(2);_data[3]-=v(3);return *this; }
+  /// Inplace scalar multiplication
+  template<typename S>
+  SGVec4& operator*=(S s)
+  { _data[0] *= s; _data[1] *= s; _data[2] *= s; _data[3] *= s; return *this; }
+  /// Inplace scalar multiplication by 1/s
+  template<typename S>
+  SGVec4& operator/=(S s)
+  { return operator*=(1/T(s)); }
+
+  /// Return an all zero vector
+  static SGVec4 zeros(void)
+  { return SGVec4(0, 0, 0, 0); }
+  /// Return unit vectors
+  static SGVec4 e1(void)
+  { return SGVec4(1, 0, 0, 0); }
+  static SGVec4 e2(void)
+  { return SGVec4(0, 1, 0, 0); }
+  static SGVec4 e3(void)
+  { return SGVec4(0, 0, 1, 0); }
+  static SGVec4 e4(void)
+  { return SGVec4(0, 0, 0, 1); }
+
+private:
+  /// The actual data
+  T _data[4];
+};
+
+/// Unary +, do nothing ...
+template<typename T>
+inline
+const SGVec4<T>&
+operator+(const SGVec4<T>& v)
+{ return v; }
+
+/// Unary -, do nearly nothing
+template<typename T>
+inline
+SGVec4<T>
+operator-(const SGVec4<T>& v)
+{ return SGVec4<T>(-v(0), -v(1), -v(2), -v(3)); }
+
+/// Binary +
+template<typename T>
+inline
+SGVec4<T>
+operator+(const SGVec4<T>& v1, const SGVec4<T>& v2)
+{ return SGVec4<T>(v1(0)+v2(0), v1(1)+v2(1), v1(2)+v2(2), v1(3)+v2(3)); }
+
+/// Binary -
+template<typename T>
+inline
+SGVec4<T>
+operator-(const SGVec4<T>& v1, const SGVec4<T>& v2)
+{ return SGVec4<T>(v1(0)-v2(0), v1(1)-v2(1), v1(2)-v2(2), v1(3)-v2(3)); }
+
+/// Scalar multiplication
+template<typename S, typename T>
+inline
+SGVec4<T>
+operator*(S s, const SGVec4<T>& v)
+{ return SGVec4<T>(s*v(0), s*v(1), s*v(2), s*v(3)); }
+
+/// Scalar multiplication
+template<typename S, typename T>
+inline
+SGVec4<T>
+operator*(const SGVec4<T>& v, S s)
+{ return SGVec4<T>(s*v(0), s*v(1), s*v(2), s*v(3)); }
+
+/// Scalar dot product
+template<typename T>
+inline
+T
+dot(const SGVec4<T>& v1, const SGVec4<T>& v2)
+{ return v1(0)*v2(0) + v1(1)*v2(1) + v1(2)*v2(2) + v1(3)*v2(3); }
+
+/// The euclidean norm of the vector, that is what most people call length
+template<typename T>
+inline
+T
+norm(const SGVec4<T>& v)
+{ return sqrt(dot(v, v)); }
+
+/// The euclidean norm of the vector, that is what most people call length
+template<typename T>
+inline
+T
+length(const SGVec4<T>& v)
+{ return sqrt(dot(v, v)); }
+
+/// The 1-norm of the vector, this one is the fastest length function we
+/// can implement on modern cpu's
+template<typename T>
+inline
+T
+norm1(const SGVec4<T>& v)
+{ return fabs(v(0)) + fabs(v(1)) + fabs(v(2)) + fabs(v(3)); }
+
+/// The euclidean norm of the vector, that is what most people call length
+template<typename T>
+inline
+SGVec4<T>
+normalize(const SGVec4<T>& v)
+{ return (1/norm(v))*v; }
+
+/// Return true if exactly the same
+template<typename T>
+inline
+bool
+operator==(const SGVec4<T>& v1, const SGVec4<T>& v2)
+{ return v1(0)==v2(0) && v1(1)==v2(1) && v1(2)==v2(2) && v1(3)==v2(3); }
+
+/// Return true if not exactly the same
+template<typename T>
+inline
+bool
+operator!=(const SGVec4<T>& v1, const SGVec4<T>& v2)
+{ return ! (v1 == v2); }
+
+/// Return true if equal to the relative tolerance tol
+template<typename T>
+inline
+bool
+equivalent(const SGVec4<T>& v1, const SGVec4<T>& v2, T rtol, T atol)
+{ return norm1(v1 - v2) < rtol*(norm1(v1) + norm1(v2)) + atol; }
+
+/// Return true if equal to the relative tolerance tol
+template<typename T>
+inline
+bool
+equivalent(const SGVec4<T>& v1, const SGVec4<T>& v2, T rtol)
+{ return norm1(v1 - v2) < rtol*(norm1(v1) + norm1(v2)); }
+
+/// Return true if about equal to roundoff of the underlying type
+template<typename T>
+inline
+bool
+equivalent(const SGVec4<T>& v1, const SGVec4<T>& v2)
+{
+  T tol = 100*SGLimits<T>::epsilon();
+  return equivalent(v1, v2, tol, tol);
+}
+
+#ifndef NDEBUG
+template<typename T>
+inline
+bool
+isNaN(const SGVec4<T>& v)
+{
+  return SGMisc<T>::isNaN(v(0)) || SGMisc<T>::isNaN(v(1))
+    || SGMisc<T>::isNaN(v(2)) || SGMisc<T>::isNaN(v(3));
+}
+#endif
+
+/// Output to an ostream
+template<typename char_type, typename traits_type, typename T>
+inline
+std::basic_ostream<char_type, traits_type>&
+operator<<(std::basic_ostream<char_type, traits_type>& s, const SGVec4<T>& v)
+{ return s << "[ " << v(0) << ", " << v(1) << ", " << v(2) << ", " << v(3) << " ]"; }
+
+/// Two classes doing actually the same on different types
+typedef SGVec4<float> SGVec4f;
+typedef SGVec4<double> SGVec4d;
+
+inline
+SGVec4f
+toVec4f(const SGVec4d& v)
+{ return SGVec4f((float)v(0), (float)v(1), (float)v(2), (float)v(3)); }
+
+inline
+SGVec4d
+toVec4d(const SGVec4f& v)
+{ return SGVec4d(v(0), v(1), v(2), v(3)); }
+
+#endif

simgear/math/fastmath.hxx

@@ -2,7 +2,7 @@
  * \file fastmath.hxx
  * fast mathematics routines.
  *
- * Refferences:
+ * References:
  *
  * A Fast, Compact Approximation of the Exponential Function
  * Nicol N. Schraudolph
@@ -25,9 +25,6 @@
 # error This library requires C++
 #endif
 
-#ifdef _MSC_VER
-#define  _USE_MATH_DEFINES
-#endif
 #include <math.h>
 
 
@@ -54,16 +51,18 @@ void fast_BSR(float &x, register unsigned long shiftAmount);
 
 inline float fast_log2 (float val)
 {
-    int * const  exp_ptr = reinterpret_cast <int *> (&val);
-    int          x = *exp_ptr;
-    const int    log_2 = ((x >> 23) & 255) - 128;
-    x &= ~(255 << 23);
-    x += 127 << 23;
-    *exp_ptr = x;
+    union {
+        float f;
+        int i;
+    } v;
+    v.f = val;
+    const int log_2 = ((v.i >> 23) & 255) - 128;
+    v.i &= ~(255 << 23);
+    v.i += 127 << 23;
 
-    val = ((-1.0f/3) * val + 2) * val - 2.0f/3;   // (1)
+    v.f = ((-1.0f/3) * v.f + 2) * v.f - 2.0f/3;   // (1)
 
-    return (val + log_2);
+    return (v.f + log_2);
 }
 
 

simgear/math/interpolater.cxx

@@ -16,10 +16,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/math/interpolater.hxx

@@ -18,10 +18,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 
@@ -30,9 +29,9 @@
 #define _INTERPOLATER_H
 
 
-#ifndef __cplusplus                                                          
+#ifndef __cplusplus
 # error This library requires C++
-#endif                                   
+#endif
 
 #include <simgear/compiler.h>
 

simgear/math/leastsqs.cxx

@@ -14,10 +14,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 //

simgear/math/leastsqs.hxx

@@ -17,10 +17,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 
@@ -29,9 +28,9 @@
 #define _LEASTSQS_H
 
 
-#ifndef __cplusplus                                                          
+#ifndef __cplusplus
 # error This library requires C++
-#endif                                   
+#endif
 
 
 /**

simgear/math/localconsts.hxx

@@ -1,42 +0,0 @@
-// localconsts.hxx -- various constant that are shared 
-//
-// Written by Curtis Olson, started September 2000.
-//
-// Copyright (C) 2000  Curtis L. Olson  - http://www.flightgear.org/~curt
-//
-// This library is free software; you can redistribute it and/or
-// modify it under the terms of the GNU Library General Public
-// License as published by the Free Software Foundation; either
-// version 2 of the License, or (at your option) any later version.
-//
-// This library 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
-// Library General Public License for more details.
-//
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
-//
-// $Id$
-
-
-#ifndef _SG_LOCAL_CONSTS_HXX
-#define _SG_LOCAL_CONSTS_HXX
-
-
-// Value of earth flattening parameter from ref [8] 
-//
-//      Note: FP = f
-//            E  = 1-f
-//            EPS = sqrt(1-(1-f)^2)
-//
-
-static const double FP =   0.003352813178;
-static const double E  =   0.996647186;
-static const double EPS =  0.081819221;
-static const double INVG = 0.031080997;
-
-
-#endif // _SG_LOCAL_CONSTS_HXX

simgear/math/point3d.hxx

@@ -20,10 +20,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 
@@ -32,9 +31,9 @@
 #define _POINT3D_HXX
 
 
-#ifndef __cplusplus                                                          
+#ifndef __cplusplus
 # error This library requires C++
-#endif                                   
+#endif
 
 #include <simgear/compiler.h>
 
@@ -52,7 +51,7 @@
 # include <math.h>
 #endif
 
-#include <simgear/math/localconsts.hxx>
+#include "SGMath.hxx"
 
 // I don't understand ... <math.h> or <cmath> should be included
 // already depending on how you defined SG_HAVE_STD_INCLUDES, but I
@@ -95,6 +94,10 @@ public:
     explicit Point3D(const double d);
     Point3D(const Point3D &p);
 
+    static Point3D fromSGGeod(const SGGeod& geod);
+    static Point3D fromSGGeoc(const SGGeoc& geoc);
+    static Point3D fromSGVec3(const SGVec3<double>& cart);
+
     // Assignment operators
 
     Point3D& operator = ( const Point3D& p );	 // assignment of a Point3D
@@ -126,6 +129,9 @@ public:
     double radius() const; // polar radius
     double elev() const;   // geodetic elevation (if specifying a surface point)
 
+    SGGeod toSGGeod(void) const;
+    SGGeoc toSGGeoc(void) const;
+
     // friends
     friend Point3D operator - (const Point3D& p);	            // -p1
     friend bool operator == (const Point3D& a, const Point3D& b);  // p1 == p2?
@@ -206,6 +212,33 @@ inline Point3D::Point3D(const Point3D& p)
     n[PX] = p.n[PX]; n[PY] = p.n[PY]; n[PZ] = p.n[PZ];
 }
 
+inline Point3D Point3D::fromSGGeod(const SGGeod& geod)
+{
+  Point3D pt;
+  pt.setlon(geod.getLongitudeRad());
+  pt.setlat(geod.getLatitudeRad());
+  pt.setelev(geod.getElevationM());
+  return pt;
+}
+
+inline Point3D Point3D::fromSGGeoc(const SGGeoc& geoc)
+{
+  Point3D pt;
+  pt.setlon(geoc.getLongitudeRad());
+  pt.setlat(geoc.getLatitudeRad());
+  pt.setradius(geoc.getRadiusM());
+  return pt;
+}
+
+inline Point3D Point3D::fromSGVec3(const SGVec3<double>& cart)
+{
+  Point3D pt;
+  pt.setx(cart.x());
+  pt.sety(cart.y());
+  pt.setz(cart.z());
+  return pt;
+}
+
 // ASSIGNMENT OPERATORS
 
 inline Point3D& Point3D::operator = (const Point3D& p)
@@ -290,6 +323,23 @@ inline double Point3D::radius() const { return n[PZ]; }
 
 inline double Point3D::elev() const { return n[PZ]; }
 
+inline SGGeod Point3D::toSGGeod(void) const
+{
+  SGGeod geod;
+  geod.setLongitudeRad(lon());
+  geod.setLatitudeRad(lat());
+  geod.setElevationM(elev());
+  return geod;
+}
+
+inline SGGeoc Point3D::toSGGeoc(void) const
+{
+  SGGeoc geoc;
+  geoc.setLongitudeRad(lon());
+  geoc.setLatitudeRad(lat());
+  geoc.setRadiusM(radius());
+  return geoc;
+}
 
 // FRIENDS
 

simgear/math/polar3d.cxx

@@ -14,83 +14,22 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 
+#ifdef HAVE_CONFIG_H
+#  include <simgear_config.h>
+#endif
 
 #include <math.h>
-#include <stdio.h>
 
 #include <simgear/constants.h>
 
 #include "polar3d.hxx"
 
-
-/**
- * Find the Altitude above the Ellipsoid (WGS84) given the Earth
- * Centered Cartesian coordinate vector Distances are specified in
- * meters.
- * @param cp point specified in cartesian coordinates
- * @return altitude above the (wgs84) earth in meters
- */
-double sgGeodAltFromCart(const Point3D& cp)
-{
-    double t_lat, x_alpha, mu_alpha;
-    double lat_geoc, radius;
-    double result;
-
-    lat_geoc = SGD_PI_2 - atan2( sqrt(cp.x()*cp.x() + cp.y()*cp.y()), cp.z() );
-    radius = sqrt( cp.x()*cp.x() + cp.y()*cp.y() + cp.z()*cp.z() );
-	
-    if( ( (SGD_PI_2 - lat_geoc) < SG_ONE_SECOND )        // near North pole
-	|| ( (SGD_PI_2 + lat_geoc) < SG_ONE_SECOND ) )   // near South pole
-    {
-	result = radius - SG_EQUATORIAL_RADIUS_M*E;
-    } else {
-	t_lat = tan(lat_geoc);
-	x_alpha = E*SG_EQUATORIAL_RADIUS_M/sqrt(t_lat*t_lat + E*E);
-	mu_alpha = atan2(sqrt(SG_EQ_RAD_SQUARE_M - x_alpha*x_alpha),E*x_alpha);
-	if (lat_geoc < 0) {
-	    mu_alpha = - mu_alpha;
-	}
-	result = (radius - x_alpha/cos(lat_geoc))*cos(mu_alpha - lat_geoc);
-    }
-
-    return(result);
-}
-
-/**
- * Convert a polar coordinate to a cartesian coordinate.  Lon and Lat
- * must be specified in radians.  The SG convention is for distances
- * to be specified in meters
- * @param p point specified in polar coordinates
- * @return the same point in cartesian coordinates
- */
-Point3D sgPolarToCart3d(const Point3D& p) {
-    double tmp = cos( p.lat() ) * p.radius();
-
-    return Point3D( cos( p.lon() ) * tmp,
-                    sin( p.lon() ) * tmp,
-                    sin( p.lat() ) * p.radius() );
-}
-
-/**
- * Convert a cartesian coordinate to polar coordinates (lon/lat
- * specified in radians.  Distances are specified in meters.
- * @param cp point specified in cartesian coordinates
- * @return the same point in polar coordinates
- */
-Point3D sgCartToPolar3d(const Point3D& cp) {
-    return Point3D( atan2( cp.y(), cp.x() ),
-                    SGD_PI_2 -
-                    atan2( sqrt(cp.x()*cp.x() + cp.y()*cp.y()), cp.z() ),
-                    sqrt(cp.x()*cp.x() + cp.y()*cp.y() + cp.z()*cp.z()) );
-}
-
 /**
  * Calculate new lon/lat given starting lon/lat, and offset radial, and
  * distance.  NOTE: starting point is specifed in radians, distance is

simgear/math/polar3d.hxx

@@ -17,10 +17,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 
@@ -29,16 +28,15 @@
 #define _POLAR3D_HXX
 
 
-#ifndef __cplusplus                                                          
+#ifndef __cplusplus
 # error This library requires C++
-#endif                                   
+#endif
 
 
-#include <math.h>
-
 #include <simgear/constants.h>
 #include <simgear/math/point3d.hxx>
 
+#include "SGMath.hxx"
 
 /** 
  * Find the Altitude above the Ellipsoid (WGS84) given the Earth
@@ -47,7 +45,12 @@
  * @param cp point specified in cartesian coordinates
  * @return altitude above the (wgs84) earth in meters
  */
-double sgGeodAltFromCart(const Point3D& cp);
+inline double sgGeodAltFromCart(const Point3D& p)
+{
+  SGGeod geod;
+  SGGeodesy::SGCartToGeod(SGVec3<double>(p.x(), p.y(), p.z()), geod);
+  return geod.getElevationM();
+}
 
 
 /**
@@ -57,7 +60,12 @@ double sgGeodAltFromCart(const Point3D& cp);
  * @param p point specified in polar coordinates
  * @return the same point in cartesian coordinates
  */
-Point3D sgPolarToCart3d(const Point3D& p);
+inline Point3D sgPolarToCart3d(const Point3D& p)
+{
+  SGVec3<double> cart;
+  SGGeodesy::SGGeocToCart(SGGeoc::fromRadM(p.lon(), p.lat(), p.radius()), cart);
+  return Point3D::fromSGVec3(cart);
+}
 
 
 /**
@@ -66,7 +74,12 @@ Point3D sgPolarToCart3d(const Point3D& p);
  * @param cp point specified in cartesian coordinates
  * @return the same point in polar coordinates
  */
-Point3D sgCartToPolar3d(const Point3D& cp);
+inline Point3D sgCartToPolar3d(const Point3D& p)
+{
+  SGGeoc geoc;
+  SGGeodesy::SGCartToGeoc(SGVec3<double>(p.x(), p.y(), p.z()), geoc);
+  return Point3D::fromSGGeoc(geoc);
+}
 
 
 /**

simgear/math/sg_geodesy.cxx

@@ -1,4 +1,10 @@
+
+#ifdef HAVE_CONFIG_H
+#  include <simgear_config.h>
+#endif
+
 #include <simgear/constants.h>
+#include "SGMath.hxx"
 #include "sg_geodesy.hxx"
 
 // Notes:
@@ -40,171 +46,6 @@ static const double STRETCH = 1.0033640898209764189003079;
 static const double POLRAD  = 6356752.3142451794975639668;
 #endif
 
-// Returns a "local" geodetic latitude: an approximation that will be
-// correct only at zero altitude.
-static double localLat(double r, double z)
-{
-    // Squash to a spherical earth, compute a tangent vector to the
-    // surface circle (in squashed space, the surface is a perfect
-    // sphere) by swapping the components and negating one, stretch to
-    // real coordinates, and take an inverse-tangent/perpedicular
-    // vector to get a local geodetic "up" vector.  (Those steps all
-    // cook down to just a few multiplies).  Then just turn it into an
-    // angle.
-    double upr = r * SQUASH;
-    double upz = z * STRETCH;
-    return atan2(upz, upr);
-}
-
-// This is the inverse of the algorithm in localLat().  It returns the
-// (cylindrical) coordinates of a surface latitude expressed as an
-// "up" unit vector.
-static void surfRZ(double upr, double upz, double* r, double* z)
-{
-    // We are
-    // converting a (2D, cylindrical) "up" vector defined by the
-    // geodetic latitude into unitless R and Z coordinates in
-    // cartesian space.
-    double R = upr * STRETCH;
-    double Z = upz * SQUASH;
-
-    // Now we need to turn R and Z into a surface point.  That is,
-    // pick a coefficient C for them such that the point is on the
-    // surface when converted to "squashed" space:
-    //  (C*R*SQUASH)^2 + (C*Z)^2 = POLRAD^2
-    //   C^2 = POLRAD^2 / ((R*SQUASH)^2 + Z^2)
-    double sr = R * SQUASH;
-    double c = POLRAD / sqrt(sr*sr + Z*Z);
-    R *= c;
-    Z *= c;
-
-    *r = R; *z = Z;
-}
-
-// Returns the insersection of the line joining the center of the
-// earth and the specified cylindrical point with the surface of the
-// WGS84 ellipsoid.  Works by finding a normalization constant (in
-// squashed space) that places the squashed point on the surface of
-// the sphere.
-static double seaLevelRadius(double r, double z)
-{
-    double sr = r * SQUASH;
-    double norm = POLRAD/sqrt(sr*sr + z*z);
-    r *= norm;
-    z *= norm;
-    return sqrt(r*r + z*z);
-}
-
-// Convert a cartexian XYZ coordinate to a geodetic lat/lon/alt.  This
-// is a "recursion relation".  In essence, it iterates over the 2D
-// part of sgGeodToCart refining its approximation at each step.  The
-// MAX_LAT_ERROR threshold is picked carefully to allow us to reach
-// the full precision of an IEEE double.  While this algorithm might
-// look slow, it's not.  It actually converges very fast indeed --
-// I've never seen it take more than six iterations under normal
-// conditions.  Three or four is more typical. (It gets slower as the
-// altitude/error gets larger; at 50000m altitude, it starts to need
-// seven loops.)  One caveat is that at *very* large altitudes, it
-// starts making very poor guesses as to latitude.  As altitude
-// approaches infinity, it should be guessing with geocentric
-// coordinates, not "local geodetic up" ones.
-void sgCartToGeod(const double* xyz, double* lat, double* lon, double* alt)
-{
-    // The error is expressed as a radian angle, and we want accuracy
-    // to 1 part in 2^50 (an IEEE double has between 51 and 52
-    // significant bits of magnitude due to the "hidden" digit; leave
-    // at least one bit free for potential slop).  In real units, this
-    // works out to about 6 nanometers.
-    static const double MAX_LAT_ERROR = 8.881784197001252e-16;
-    double x = xyz[0], y = xyz[1], z = xyz[2];
-
-    // Longitude is trivial.  Convert to cylindrical "(r, z)"
-    // coordinates while we're at it.
-    *lon = atan2(y, x);
-    double r = sqrt(x*x + y*y);
-
-    double lat1, lat2 = localLat(r, z);
-    double r2, z2, dot;
-    do {
-        lat1 = lat2;
-
-        // Compute an "up" vector
-        double upr = cos(lat1);
-        double upz = sin(lat1);
-
-        // Find the surface point with that latitude
-        surfRZ(upr, upz, &r2, &z2);
-
-        // Convert r2z2 to the vector pointing from the surface to rz
-        r2 = r - r2;
-        z2 = z - z2;
-
-        // Dot it with "up" to get an approximate altitude
-        dot = r2*upr + z2*upz;
-
-        // And compute an approximate geodetic surface coordinate
-        // using that altitude, so now: R2Z2 = RZ - ((RZ - SURF) dot
-        // UP)
-        r2 = r - dot * upr;
-        z2 = z - dot * upz;
-
-        // Find the latitude of *that* point, and iterate
-        lat2 = localLat(r2, z2);
-    } while(fabs(lat2 - lat1) > MAX_LAT_ERROR);
-
-    // All done!  We have an accurate geodetic lattitude, now
-    // calculate the altitude as a cartesian distance between the
-    // final geodetic surface point and the initial r/z coordinate.
-    *lat = lat1;
-    double dr = r - r2;
-    double dz = z - z2;
-    double altsign = (dot > 0) ? 1 : -1;
-    *alt = altsign * sqrt(dr*dr + dz*dz);
-}
-
-void sgGeodToCart(double lat, double lon, double alt, double* xyz)
-{
-    // This is the inverse of the algorithm in localLat().  We are
-    // converting a (2D, cylindrical) "up" vector defined by the
-    // geodetic latitude into unitless R and Z coordinates in
-    // cartesian space.
-    double upr = cos(lat);
-    double upz = sin(lat);
-    double r, z;
-    surfRZ(upr, upz, &r, &z);
-
-    // Add the altitude using the "up" unit vector we calculated
-    // initially.
-    r += upr * alt;
-    z += upz * alt;
-
-    // Finally, convert from cylindrical to cartesian
-    xyz[0] = r * cos(lon);
-    xyz[1] = r * sin(lon);
-    xyz[2] = z;
-}
-
-void sgGeocToGeod(double lat_geoc, double radius,
-                  double *lat_geod, double *alt, double *sea_level_r)
-{
-    // Build a fake cartesian point, and run it through CartToGeod
-    double lon_dummy, xyz[3];
-    xyz[0] = cos(lat_geoc) * radius;
-    xyz[1] = 0;
-    xyz[2] = sin(lat_geoc) * radius;
-    sgCartToGeod(xyz, lat_geod, &lon_dummy, alt);
-    *sea_level_r = seaLevelRadius(xyz[0], xyz[2]);
-}
-
-void sgGeodToGeoc(double lat_geod, double alt,
-                  double *sl_radius, double *lat_geoc)
-{
-    double xyz[3];
-    sgGeodToCart(lat_geod, 0, alt, xyz);
-    *lat_geoc = atan2(xyz[2], xyz[0]);
-    *sl_radius = seaLevelRadius(xyz[0], xyz[2]);
-}
-
 ////////////////////////////////////////////////////////////////////////
 //
 // Direct and inverse distance functions 

simgear/math/sg_geodesy.hxx

@@ -2,6 +2,19 @@
 #define _SG_GEODESY_HXX
 
 #include <simgear/math/point3d.hxx>
+#include "SGMath.hxx"
+
+// Returns the insersection of the line joining the center of the
+// earth and the specified cylindrical point with the surface of the
+// WGS84 ellipsoid.  Works by finding a normalization constant (in
+// squashed space) that places the squashed point on the surface of
+// the sphere.
+inline double seaLevelRadius(double r, double z)
+{
+    double sr = r * SGGeodesy::SQUASH;
+    double zz = z*z;
+    return SGGeodesy::POLRAD*sqrt((r*r + zz)/(sr*sr + zz));
+}
 
 /** 
  * Convert from geocentric coordinates to geodetic coordinates
@@ -12,9 +25,17 @@
  * @param sea_level_r (out) radius from earth center to sea level at
  *        local vertical (surface normal) of C.G. (meters)
  */
-void sgGeocToGeod(double lat_geoc, double radius,
-                  double *lat_geod, double *alt, double *sea_level_r);
-
+inline void sgGeocToGeod(double lat_geoc, double radius,
+                         double *lat_geod, double *alt, double *sea_level_r)
+{
+  SGVec3<double> cart;
+  SGGeodesy::SGGeocToCart(SGGeoc::fromRadM(0, lat_geoc, radius), cart);
+  SGGeod geod;
+  SGGeodesy::SGCartToGeod(cart, geod);
+  *lat_geod = geod.getLatitudeRad();
+  *alt = geod.getElevationM();
+  *sea_level_r = SGGeodesy::SGGeodToSeaLevelRadius(geod);
+}
 
 /**
  * Convert from geodetic coordinates to geocentric coordinates.
@@ -31,8 +52,17 @@ void sgGeocToGeod(double lat_geoc, double radius,
  * @param sl_radius (out) SEA LEVEL radius to earth center (meters)
  * @param lat_geoc (out) Geocentric latitude, radians, + = North
  */
-void sgGeodToGeoc(double lat_geod, double alt,
-                  double *sl_radius, double *lat_geoc );
+inline void sgGeodToGeoc(double lat_geod, double alt,
+                         double *sl_radius, double *lat_geoc)
+{
+  SGVec3<double> cart;
+  SGGeodesy::SGGeodToCart(SGGeod::fromRadM(0, lat_geod, alt), cart);
+  SGGeoc geoc;
+  SGGeodesy::SGCartToGeoc(cart, geoc);
+  *lat_geoc = geoc.getLatitudeRad();
+  *sl_radius = seaLevelRadius(cart(0), cart(2));
+}
+
 
 /**
  * Convert a cartesian point to a geodetic lat/lon/altitude.
@@ -42,7 +72,14 @@ void sgGeodToGeoc(double lat_geod, double alt,
  * @param lon (out) Longitude, in radians
  * @param alt (out) Altitude, in meters above the WGS84 ellipsoid
  */
-void sgCartToGeod(const double* xyz, double* lat, double* lon, double* alt);
+inline void sgCartToGeod(const double* xyz, double* lat, double* lon, double* alt)
+{
+  SGGeod geod;
+  SGGeodesy::SGCartToGeod(SGVec3<double>(xyz), geod);
+  *lat = geod.getLatitudeRad();
+  *lon = geod.getLongitudeRad();
+  *alt = geod.getElevationM();
+}
 
 /**
  * Convert a cartesian point to a geodetic lat/lon/altitude.
@@ -53,10 +90,9 @@ void sgCartToGeod(const double* xyz, double* lat, double* lon, double* alt);
  */
 inline Point3D sgCartToGeod(const Point3D& p)
 {
-    double lat, lon, alt, xyz[3];
-    xyz[0] = p.x(); xyz[1] = p.y(); xyz[2] = p.z(); 
-    sgCartToGeod(xyz, &lat, &lon, &alt);
-    return Point3D(lon, lat, alt);
+  SGGeod geod;
+  SGGeodesy::SGCartToGeod(SGVec3<double>(p.x(), p.y(), p.z()), geod);
+  return Point3D::fromSGGeod(geod);
 }
 
 
@@ -68,7 +104,14 @@ inline Point3D sgCartToGeod(const Point3D& p)
  * @param alt (in) Altitude, in meters above the WGS84 ellipsoid
  * @param xyz (out) Pointer to cartesian point.
  */
-void sgGeodToCart(double lat, double lon, double alt, double* xyz);
+inline void sgGeodToCart(double lat, double lon, double alt, double* xyz)
+{
+  SGVec3<double> cart;
+  SGGeodesy::SGGeodToCart(SGGeod::fromRadM(lon, lat, alt), cart);
+  xyz[0] = cart(0);
+  xyz[1] = cart(1);
+  xyz[2] = cart(2);
+}
 
 /**
  * Convert a geodetic lat/lon/altitude to a cartesian point.
@@ -79,9 +122,9 @@ void sgGeodToCart(double lat, double lon, double alt, double* xyz);
  */
 inline Point3D sgGeodToCart(const Point3D& geod)
 {
-    double xyz[3];
-    sgGeodToCart(geod.lat(), geod.lon(), geod.elev(), xyz);
-    return Point3D(xyz[0], xyz[1], xyz[2]);
+  SGVec3<double> cart;
+  SGGeodesy::SGGeodToCart(SGGeod::fromRadM(geod.lon(), geod.lat(), geod.elev()), cart);
+  return Point3D::fromSGVec3(cart);
 }
 
 /**

simgear/math/sg_memory.h

@@ -14,10 +14,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 
@@ -26,7 +25,7 @@
 #define _SG_MEMORY_H
 
 #ifdef HAVE_CONFIG_H
-#  include <config.h>
+#  include <simgear_config.h>
 #endif
 
 #ifdef HAVE_MEMCPY

simgear/math/sg_random.c

@@ -14,10 +14,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/math/sg_random.h

@@ -17,10 +17,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/math/sg_types.hxx

@@ -17,10 +17,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 
@@ -29,9 +28,9 @@
 #define _SG_TYPES_HXX
 
 
-#ifndef __cplusplus                                                          
+#ifndef __cplusplus
 # error This library requires C++
-#endif                                   
+#endif
 
 
 #include <simgear/compiler.h>

simgear/math/vector.cxx

@@ -14,10 +14,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 

simgear/math/vector.hxx

@@ -17,10 +17,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 
@@ -29,9 +28,9 @@
 #define _VECTOR_HXX
 
 
-#ifndef __cplusplus                                                          
+#ifndef __cplusplus
 # error This library requires C++
-#endif                                   
+#endif
 
 #include <simgear/compiler.h>
 

simgear/misc/.cvsignore

@@ -3,3 +3,4 @@ Makefile
 Makefile.in
 props_test
 tabbed_value_test
+swap_test

simgear/misc/interpolator.cxx

@@ -12,10 +12,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 
 #include "interpolator.hxx"
 

simgear/misc/interpolator.hxx

@@ -19,10 +19,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 
 #include <simgear/props/props.hxx>
 #include <simgear/structure/subsystem_mgr.hxx>

simgear/misc/sg_path.cxx

@@ -15,10 +15,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 
@@ -26,8 +25,13 @@
 #include <simgear/compiler.h>
 
 #include <simgear_config.h>
+#include <simgear/debug/logstream.hxx>
 #include <stdio.h>
-
+#include <sys/stat.h>
+#include <sys/stat.h>
+#if defined( _MSC_VER) || defined(__MINGW32__)
+#  include <direct.h>
+#endif
 #include "sg_path.hxx"
 
 
@@ -182,6 +186,67 @@ bool SGPath::exists() const {
     return true;
 }
 
+#if defined( _MSC_VER) || defined(__MINGW32__)
+#  define sgMkDir(d,m)       _mkdir(d)
+#else
+#  define sgMkDir(d,m)       mkdir(d,m)
+#endif
+
+
+void SGPath::create_dir( mode_t mode ) {
+    string_list dirlist = sgPathSplit(dir());
+    string path = dirlist[0];
+    string_list path_elements = sgPathBranchSplit(path);
+    bool absolute = !path.empty() && path[0] == sgDirPathSep;
+
+    unsigned int i = 1;
+    SGPath dir = absolute ? string( 1, sgDirPathSep ) : "";
+    dir.concat( path_elements[0] );
+#if defined( _MSC_VER) || defined(__MINGW32__)
+    if ( dir.str().find(':') != string::npos && path_elements.size() >= 2 ) {
+        dir.append( path_elements[1] );
+        i = 2;
+    }
+#endif
+    struct stat info;
+    int r;
+    for(; ( r = stat( dir.c_str(), &info ) ) == 0 && i < path_elements.size(); i++) {
+        dir.append(path_elements[i]);
+    }
+    if ( r == 0 ) {
+        return; // Directory already exists
+    }
+    if ( sgMkDir( dir.c_str(), mode) ) {
+        SG_LOG( SG_IO, SG_ALERT, "Error creating directory: " + dir.str() );
+        return;
+    }
+    for(;i < path_elements.size(); i++) {
+        dir.append(path_elements[i]);
+        if ( sgMkDir( dir.c_str(), mode) ) {
+            SG_LOG( SG_IO, SG_ALERT, "Error creating directory: " + dir.str() );
+            break;
+        }
+    }
+}
+
+string_list sgPathBranchSplit( const string &dirpath ) {
+    string_list path_elements;
+    string element, path = dirpath;
+    while ( path.size() ) {
+        size_t p = path.find( sgDirPathSep );
+        if ( p != string::npos ) {
+            element = path.substr( 0, p );
+            path.erase( 0, p + 1 );
+        } else {
+            element = path;
+            path = "";
+        }
+        if ( element.size() )
+            path_elements.push_back( element );
+    }
+    return path_elements;
+}
+
 
 string_list sgPathSplit( const string &search_path ) {
     string tmp = search_path;

simgear/misc/sg_path.hxx

@@ -18,10 +18,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$
 
@@ -29,6 +28,7 @@
 #ifndef _SG_PATH_HXX
 #define _SG_PATH_HXX
 
+#include <sys/types.h>
 
 #include <simgear/compiler.h>
 #include STL_STRING
@@ -37,6 +37,9 @@
 
 SG_USING_STD(string);
 
+#ifdef _MSC_VER
+  typedef int mode_t;
+#endif
 
 /**
  * A class to hide path separator difference across platforms and assist
@@ -116,12 +119,14 @@ public:
      */
     string extension() const;
 
-    /** Get the path string
+    /**
+     * Get the path string
      * @return path string
      */
     string str() const { return path; }
 
-    /** Get the path string
+    /**
+     * Get the path string
      * @return path in "C" string (ptr to char array) form.
      */
     const char* c_str() { return path.c_str(); }
@@ -132,6 +137,11 @@ public:
      */
     bool exists() const;
 
+    /**
+     * Create the designated directory.
+     */
+    void create_dir(mode_t mode);
+
 private:
 
     void fix();
@@ -139,6 +149,11 @@ private:
 };
 
 
+/**
+ * Split a directory string into a list of it's parent directories.
+ */
+string_list sgPathBranchSplit( const string &path );
+
 /**
  * Split a directory search path into a vector of individual paths
  */

simgear/misc/sgstream.cxx

@@ -14,10 +14,9 @@
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 // Library General Public License for more details.
 //
-// You should have received a copy of the GNU Library General Public
-// License along with this library; if not, write to the
-// Free Software Foundation, Inc., 59 Temple Place - Suite 330,
-// Boston, MA  02111-1307, USA.
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 //
 // $Id$