28 files changed, 1823 insertions, 1708 deletions

diff --git a/linden/indra/llcommon/indra_constants.h b/linden/indra/llcommon/indra_constants.h
index 279d280..ae1e04a 100644
--- a/linden/indra/llcommon/indra_constants.h
+++ b/linden/indra/llcommon/indra_constants.h
@@ -151,7 +151,12 @@ const	U32		PORT_DISCOVERY_RANGE_MAX		= PORT_DISCOVERY_RANGE_MIN + 50;
 const   char    LAND_LAYER_CODE                                 = 'L';
 const   char    WATER_LAYER_CODE                                = 'W';
 const   char    WIND_LAYER_CODE                                 = '7';
-const   char    CLOUD_LAYER_CODE                                = '8';
+const   char    CLOUD_LAYER_CODE                                = '8';

+

+// Extended land layer for Aurora Sim

+const   char    AURORA_LAND_LAYER_CODE                  = 'M';

+const   char    AURORA_WIND_LAYER_CODE                  = '9';

+const   char    AURORA_CLOUD_LAYER_CODE                 = ':';
 
 // keys
 // Bit masks for various keyboard modifier keys.
diff --git a/linden/indra/llinventory/llparcel.h b/linden/indra/llinventory/llparcel.h
index 8faa673..366f589 100644
--- a/linden/indra/llinventory/llparcel.h
+++ b/linden/indra/llinventory/llparcel.h
@@ -66,7 +66,8 @@ const S32 PARCEL_PASS_PRICE_DEFAULT = 10;
 const F32 PARCEL_PASS_HOURS_DEFAULT = 1.f;
 
 // Number of "chunks" in which parcel overlay data is sent
-// Chunk 0 = southern rows, entire width
+// Chunk 0 = southern rows, entire width

+// NOTE: NOT USABLE FOR VAR SIZED REGIONS!
 const S32 PARCEL_OVERLAY_CHUNKS = 4;
 
 // Bottom three bits are a color index for the land overlay
diff --git a/linden/indra/llmessage/message_prehash.cpp b/linden/indra/llmessage/message_prehash.cpp
index ecb45bf..a8b021f 100644
--- a/linden/indra/llmessage/message_prehash.cpp
+++ b/linden/indra/llmessage/message_prehash.cpp
@@ -586,7 +586,9 @@ char* _PREHASH_NotecardItemID = LLMessageStringTable::getInstance()->getString("
 char* _PREHASH_LastName = LLMessageStringTable::getInstance()->getString("LastName");
 char* _PREHASH_From = LLMessageStringTable::getInstance()->getString("From");
 char* _PREHASH_RoleChange = LLMessageStringTable::getInstance()->getString("RoleChange");
-char* _PREHASH_Port = LLMessageStringTable::getInstance()->getString("Port");
+char* _PREHASH_Port = LLMessageStringTable::getInstance()->getString("Port");

+char* _PREHASH_RegionSizeX = LLMessageStringTable::getInstance()->getString("RegionSizeX");

+char* _PREHASH_RegionSizeY = LLMessageStringTable::getInstance()->getString("RegionSizeY");
 char* _PREHASH_MemberTitle = LLMessageStringTable::getInstance()->getString("MemberTitle");
 char* _PREHASH_LogParcelChanges = LLMessageStringTable::getInstance()->getString("LogParcelChanges");
 char* _PREHASH_AgentCachedTextureResponse = LLMessageStringTable::getInstance()->getString("AgentCachedTextureResponse");
diff --git a/linden/indra/llmessage/message_prehash.h b/linden/indra/llmessage/message_prehash.h
index 8516483..adf27ee 100644
--- a/linden/indra/llmessage/message_prehash.h
+++ b/linden/indra/llmessage/message_prehash.h
@@ -586,7 +586,9 @@ extern char * _PREHASH_NotecardItemID;
 extern char * _PREHASH_LastName;
 extern char * _PREHASH_From;
 extern char * _PREHASH_RoleChange;
-extern char * _PREHASH_Port;
+extern char * _PREHASH_Port;

+extern char * _PREHASH_RegionSizeX;

+extern char * _PREHASH_RegionSizeY;
 extern char * _PREHASH_MemberTitle;
 extern char * _PREHASH_LogParcelChanges;
 extern char * _PREHASH_AgentCachedTextureResponse;
diff --git a/linden/indra/llmessage/patch_code.cpp b/linden/indra/llmessage/patch_code.cpp
index 90fb236..4b4313d 100644
--- a/linden/indra/llmessage/patch_code.cpp
+++ b/linden/indra/llmessage/patch_code.cpp
@@ -235,7 +235,7 @@ void	decode_patch_group_header(LLBitPack &bitpack, LLGroupHeader *gopp)
         gPatchSize = gopp->patch_size; 
 }
 
-void    decode_patch_header(LLBitPack &bitpack, LLPatchHeader *ph)
+void    decode_patch_header(LLBitPack &bitpack, LLPatchHeader *ph, BOOL b_large_patch)
 {
         U8 retvalu8;
 
@@ -274,15 +274,18 @@ void	decode_patch_header(LLBitPack &bitpack, LLPatchHeader *ph)
 #endif
         ph->range = retvalu16;
 
-        retvalu16 = 0;
+        retvalu32 = 0;
 #ifdef LL_BIG_ENDIAN
         ret = (U8 *)&retvalu16;
         bitpack.bitUnpack(&(ret[1]), 8);
         bitpack.bitUnpack(&(ret[0]), 2);
 #else
-        bitpack.bitUnpack((U8 *)&retvalu16, 10);
-#endif
-        ph->patchids = retvalu16;
+        if (b_large_patch)

+                bitpack.bitUnpack((U8 *)&retvalu32, 32);

+        else

+                bitpack.bitUnpack((U8 *)&retvalu32, 10);

+#endif

+        ph->patchids = retvalu32;
 
         gWordBits = (ph->quant_wbits & 0xf) + 2;
 }
diff --git a/linden/indra/llmessage/patch_code.h b/linden/indra/llmessage/patch_code.h
index 82fa6bb..dbfdf70 100644
--- a/linden/indra/llmessage/patch_code.h
+++ b/linden/indra/llmessage/patch_code.h
@@ -46,7 +46,7 @@ void	end_patch_coding(LLBitPack &bitpack);
 
 void    init_patch_decoding(LLBitPack &bitpack);
 void    decode_patch_group_header(LLBitPack &bitpack, LLGroupHeader *gopp);
-void    decode_patch_header(LLBitPack &bitpack, LLPatchHeader *ph);
+void    decode_patch_header(LLBitPack &bitpack, LLPatchHeader *ph, BOOL b_large_patch);
 void    decode_patch(LLBitPack &bitpack, S32 *patches);
 
 #endif
diff --git a/linden/indra/llmessage/patch_dct.h b/linden/indra/llmessage/patch_dct.h
index 663e146..ba21be7 100644
--- a/linden/indra/llmessage/patch_dct.h
+++ b/linden/indra/llmessage/patch_dct.h
@@ -79,7 +79,7 @@ public:
         F32     dc_offset;              // 4 bytes
         U16     range;                  // 2 = 7 ((S16) FP range (breaks if we need > 32K meters in 1 patch)
         U8      quant_wbits;    // 1 = 8 (upper 4 bits is quant - 2, lower 4 bits is word bits - 2)
-        U16     patchids;               // 2 = 10 (actually only uses 10 bits, 5 for each)
+        U32     patchids;               // 2 = 10 (actually only uses 10 bits, 5 for each)
 };
 
 // Compression routines
diff --git a/linden/indra/llwindow/glh/glh_linear.h b/linden/indra/llwindow/glh/glh_linear.h
index 04ae1bd..bb59d7e 100755
--- a/linden/indra/llwindow/glh/glh_linear.h
+++ b/linden/indra/llwindow/glh/glh_linear.h
@@ -1,1621 +1,1621 @@
-/*
-    glh - is a platform-indepenedent C++ OpenGL helper library 
-
-
-    Copyright (c) 2000 Cass Everitt
-        Copyright (c) 2000 NVIDIA Corporation
-    All rights reserved.
-
-    Redistribution and use in source and binary forms, with or
-        without modification, are permitted provided that the following
-        conditions are met:
-
-     * Redistributions of source code must retain the above
-           copyright notice, this list of conditions and the following
-           disclaimer.
-
-     * Redistributions in binary form must reproduce the above
-           copyright notice, this list of conditions and the following
-           disclaimer in the documentation and/or other materials
-           provided with the distribution.
-
-     * The names of contributors to this software may not be used
-           to endorse or promote products derived from this software
-           without specific prior written permission. 
-
-       THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
-           ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
-           LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
-           FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
-           REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
-           INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
-           BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
-           LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
-           CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
-           LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
-           ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
-           POSSIBILITY OF SUCH DAMAGE. 
-
-
-    Cass Everitt - cass@r3.nu
-*/
-
-/*
-glh_linear.h
-*/
-
-// Author:  Cass W. Everitt
-
-#ifndef GLH_LINEAR_H
-#define GLH_LINEAR_H
-
-#include <memory.h>
-#include <math.h>
-#include <assert.h>
-
-// only supports float for now...
-#define GLH_REAL_IS_FLOAT
-
-#ifdef GLH_REAL_IS_FLOAT
-# define GLH_REAL float
-# define GLH_REAL_NAMESPACE ns_float
-#endif
-
-#define     GLH_QUATERNION_NORMALIZATION_THRESHOLD  64
-
-#define     GLH_RAD_TO_DEG      GLH_REAL(57.2957795130823208767981548141052)
-#define     GLH_DEG_TO_RAD      GLH_REAL(0.0174532925199432957692369076848861)
-#define     GLH_ZERO            GLH_REAL(0.0)
-#define     GLH_ONE             GLH_REAL(1.0)
-#define     GLH_TWO             GLH_REAL(2.0)
-#define     GLH_EPSILON         GLH_REAL(10e-6)
-#define     GLH_PI              GLH_REAL(3.1415926535897932384626433832795)    
-
-#define     equivalent(a,b)     (((a < b + GLH_EPSILON) && (a > b - GLH_EPSILON)) ? true : false)
-
-namespace glh
-{
-
-        inline GLH_REAL to_degrees(GLH_REAL radians) { return radians*GLH_RAD_TO_DEG; }
-        inline GLH_REAL to_radians(GLH_REAL degrees) { return degrees*GLH_DEG_TO_RAD; }
-
-        // forward declarations for friend template functions.
-        template <int N, class T> class vec;
-        
-        // forward declarations for friend template functions.
-        template <int N, class T>
-                bool operator == ( const vec<N,T> & v1, const vec<N,T> & v2 );
-        
-        // forward declarations for friend template functions.
-        template <int N, class T>
-                bool operator != ( const vec<N,T> & v1, const vec<N,T> & v2 );
-
-        template <int N, class T>       
-        class vec
-        {                               
-    public:
-                int size() const { return N; }
-                
-                vec(const T & t = T()) 
-                { for(int i = 0; i < N; i++) v[i] = t; }
-                vec(const T * tp)
-                { for(int i = 0; i < N; i++) v[i] = tp[i]; }
-                
-                const T * get_value() const
-                { return v; }
-                
-                
-                T dot( const vec<N,T> & rhs ) const
-                { 
-                        T r = 0;
-                        for(int i = 0; i < N; i++) r += v[i]*rhs.v[i];
-                        return r;
-                }
-                
-                T length() const
-                {
-                        T r = 0;
-                        for(int i = 0; i < N; i++) r += v[i]*v[i]; 
-                        return T(sqrt(r));
-                }       
-                
-                T square_norm() const
-                {
-                        T r = 0;
-                        for(int i = 0; i < N; i++) r += v[i]*v[i]; 
-                        return r;
-                }       
-                
-                void  negate()
-                { for(int i = 0; i < N; i++) v[i] = -v[i]; }
-                
-                
-                T normalize() 
-                { 
-                        T sum(0);
-                        for(int i = 0; i < N; i++) 
-                sum += v[i]*v[i];
-                        sum = T(sqrt(sum));
-            if (sum > GLH_EPSILON)
-                            for(int i = 0; i < N; i++) 
-                    v[i] /= sum;
-                        return sum;
-                }
-                
-                
-                vec<N,T> & set_value( const T * rhs )
-                { for(int i = 0; i < N; i++) v[i] = rhs[i]; return *this; }
-                
-                T & operator [] ( int i )
-                { return v[i]; }
-                
-                const T & operator [] ( int i ) const
-                { return v[i]; }
-
-                vec<N,T> & operator *= ( T d )
-                { for(int i = 0; i < N; i++) v[i] *= d; return *this;}
-                
-                vec<N,T> & operator *= ( const vec<N,T> & u )
-                { for(int i = 0; i < N; i++) v[i] *= u[i]; return *this;}
-                
-                vec<N,T> & operator /= ( T d )
-                { if(d == 0) return *this; for(int i = 0; i < N; i++) v[i] /= d; return *this;}
-                
-                vec<N,T> & operator += ( const vec<N,T> & u )
-                { for(int i = 0; i < N; i++) v[i] += u.v[i]; return *this;}
-                
-                vec<N,T> & operator -= ( const vec<N,T> & u )
-                { for(int i = 0; i < N; i++) v[i] -= u.v[i]; return *this;}
-                
-                
-                vec<N,T> operator - () const
-                { vec<N,T> rv = v; rv.negate(); return rv; }
-                
-                vec<N,T> operator + ( const vec<N,T> &v) const
-                { vec<N,T> rt(*this); return rt += v; }
-                
-                vec<N,T> operator - ( const vec<N,T> &v) const
-                { vec<N,T> rt(*this); return rt -= v; }
-                
-                vec<N,T> operator * ( T d) const
-                { vec<N,T> rt(*this); return rt *= d; }
-                
-                friend bool operator == <> ( const vec<N,T> &v1, const vec<N,T> &v2 );
-                friend bool operator != <> ( const vec<N,T> &v1, const vec<N,T> &v2 );
-                
-                
-        //protected:
-                T v[N];
-        };
-        
-        
-        
-        // vector friend operators
-        
-        template <int N, class T> inline
-                vec<N,T> operator * ( const vec<N,T> & b, T d )
-        {
-                vec<N,T> rt(b);
-                return rt *= d;
-        }
-
-        template <int N, class T> inline
-                vec<N,T> operator * ( T d, const vec<N,T> & b )
-        { return b*d; }
-        
-        template <int N, class T> inline
-                vec<N,T> operator * ( const vec<N,T> & b, const vec<N,T> & d )
-        {
-                vec<N,T> rt(b);
-                return rt *= d;
-        }
-
-        template <int N, class T> inline
-                vec<N,T> operator / ( const vec<N,T> & b, T d )
-        { vec<N,T> rt(b); return rt /= d; }
-        
-        template <int N, class T> inline
-                vec<N,T> operator + ( const vec<N,T> & v1, const vec<N,T> & v2 )
-        { vec<N,T> rt(v1); return rt += v2; }
-        
-        template <int N, class T> inline
-                vec<N,T> operator - ( const vec<N,T> & v1, const vec<N,T> & v2 )
-        { vec<N,T> rt(v1); return rt -= v2; }
-        
-        
-        template <int N, class T> inline
-                bool operator == ( const vec<N,T> & v1, const vec<N,T> & v2 )
-        {
-                for(int i = 0; i < N; i++)
-                        if(v1.v[i] != v2.v[i])
-                                return false;
-                        return true;
-        }
-        
-        template <int N, class T> inline
-                bool operator != ( const vec<N,T> & v1, const vec<N,T> & v2 )
-        { return !(v1 == v2); }
-        
-
-        typedef vec<3,unsigned char> vec3ub;
-        typedef vec<4,unsigned char> vec4ub;
-
-
-
-
-
-        namespace GLH_REAL_NAMESPACE
-        {
-        typedef GLH_REAL real;
-
-        class line;
-        class plane;
-        class matrix4;
-        class quaternion;
-        typedef quaternion rotation; 
-  
-        class vec2 : public vec<2,real>
-        {
-    public:
-                vec2(const real & t = real()) : vec<2,real>(t)
-                {}
-                vec2(const vec<2,real> & t) : vec<2,real>(t)
-                {}
-                vec2(const real * tp) : vec<2,real>(tp)
-                {}
-                
-                vec2(real x, real y )
-                { v[0] = x; v[1] = y; }
-                
-                void get_value(real & x, real & y) const
-                { x = v[0]; y = v[1]; }
-                
-                vec2 & set_value( const real & x, const real & y)
-                { v[0] = x; v[1] = y; return *this; }
-                
-        };
-        
-        
-        class vec3 : public vec<3,real>
-        {
-    public:
-                vec3(const real & t = real()) : vec<3,real>(t)
-                {}
-                vec3(const vec<3,real> & t) : vec<3,real>(t)
-                {}
-                vec3(const real * tp) : vec<3,real>(tp)
-                {}
-                
-                vec3(real x, real y, real z)
-                { v[0] = x; v[1] = y; v[2] = z; }
-                
-                void get_value(real & x, real & y, real & z) const
-                { x = v[0]; y = v[1]; z = v[2]; }
-                
-                vec3 cross( const vec3 &rhs ) const
-                {
-                        vec3 rt;
-                        rt.v[0] = v[1]*rhs.v[2]-v[2]*rhs.v[1];
-                        rt.v[1] = v[2]*rhs.v[0]-v[0]*rhs.v[2];
-                        rt.v[2] = v[0]*rhs.v[1]-v[1]*rhs.v[0];  
-                        return rt;
-                }
-                
-                vec3 & set_value( const real & x, const real & y, const real & z)
-                { v[0] = x; v[1] = y; v[2] = z; return *this; }
-                
-        };
-
-                
-    class vec4 : public vec<4,real>
-    {
-    public:
-        vec4(const real & t = real()) : vec<4,real>(t)
-        {}
-        vec4(const vec<4,real> & t) : vec<4,real>(t)
-        {}
-
-        vec4(const vec<3,real> & t, real fourth)
-
-        { v[0] = t.v[0]; v[1] = t.v[1]; v[2] = t.v[2]; v[3] = fourth; }
-        vec4(const real * tp) : vec<4,real>(tp)
-        {}
-        vec4(real x, real y, real z, real w)
-        { v[0] = x; v[1] = y; v[2] = z; v[3] = w; }
-
-        void get_value(real & x, real & y, real & z, real & w) const
-        { x = v[0]; y = v[1]; z = v[2]; w = v[3]; }
-  
-        vec4 & set_value( const real & x, const real & y, const real & z, const real & w)
-        { v[0] = x; v[1] = y; v[2] = z; v[3] = w; return *this; }
-    };
-
-    inline
-    vec3 homogenize(const vec4 & v)
-    {
-      vec3 rt;
-      assert(v.v[3] != GLH_ZERO);
-      rt.v[0] = v.v[0]/v.v[3];
-      rt.v[1] = v.v[1]/v.v[3];
-      rt.v[2] = v.v[2]/v.v[3];
-      return rt;
-    }
-  
-
-
-    class line
-    {
-    public:
-  
-        line()
-        { set_value(vec3(0,0,0),vec3(0,0,1)); }
-
-        line( const vec3 & p0, const vec3 &p1)
-        { set_value(p0,p1); }
-
-        void set_value( const vec3 &p0, const vec3 &p1)
-        {
-          position = p0;
-          direction = p1-p0;
-          direction.normalize();
-        }
-  
-        bool get_closest_points(const line &line2, 
-                                                  vec3 &pointOnThis,
-                                                  vec3 &pointOnThat)
-        {
-  
-          // quick check to see if parallel -- if so, quit.
-          if(fabs(direction.dot(line2.direction)) == 1.0)
-                  return 0;
-          line l2 = line2;
-  
-          // Algorithm: Brian Jean
-          // 
-          register real u;
-          register real v;
-          vec3 Vr = direction;
-          vec3 Vs = l2.direction;
-          register real Vr_Dot_Vs = Vr.dot(Vs);
-          register real detA = real(1.0 - (Vr_Dot_Vs * Vr_Dot_Vs));
-          vec3 C = l2.position - position;
-          register real C_Dot_Vr =  C.dot(Vr);
-          register real C_Dot_Vs =  C.dot(Vs);
-  
-          u = (C_Dot_Vr - Vr_Dot_Vs * C_Dot_Vs)/detA;
-          v = (C_Dot_Vr * Vr_Dot_Vs - C_Dot_Vs)/detA;
-  
-          pointOnThis = position;
-          pointOnThis += direction * u;
-          pointOnThat = l2.position;
-          pointOnThat += l2.direction * v;
-  
-          return 1;
-        }
-  
-        vec3 get_closest_point(const vec3 &point)
-        {
-          vec3 np = point - position;
-          vec3 rp = direction*direction.dot(np)+position;
-          return rp;
-        }
-  
-        const vec3 & get_position() const {return position;}
-
-        const vec3 & get_direction() const {return direction;}
-  
-    //protected:
-        vec3 position;
-        vec3 direction;
-    };
-  
-  
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-  
-  
-
-
-
-
-
-
-
-
-
-
-
-  // matrix
-
-  
-  class matrix4
-  {
-    
-  public:
-        
-    matrix4() { make_identity(); }
-    
-        matrix4( real r ) 
-        { set_value(r); }
-
-        matrix4( real * m )
-        { set_value(m); }
-    
-    matrix4( real a00, real a01, real a02, real a03,
-               real a10, real a11, real a12, real a13,
-                   real a20, real a21, real a22, real a23,
-                   real a30, real a31, real a32, real a33 )
-        {
-                element(0,0) = a00;
-                element(0,1) = a01;
-                element(0,2) = a02;
-                element(0,3) = a03;
-                
-                element(1,0) = a10;
-                element(1,1) = a11;
-                element(1,2) = a12;
-                element(1,3) = a13;
-                
-                element(2,0) = a20;
-                element(2,1) = a21;
-                element(2,2) = a22;
-                element(2,3) = a23;
-                
-                element(3,0) = a30;
-                element(3,1) = a31;
-                element(3,2) = a32;
-                element(3,3) = a33;
-        }
-            
-    
-    void get_value( real * mp ) const
-        {
-                int c = 0;
-                for(int j=0; j < 4; j++)
-                        for(int i=0; i < 4; i++)
-                                mp[c++] = element(i,j);
-        }
-    
-    
-    const real * get_value() const
-        { return m; }
-    
-        void set_value( real * mp)
-        {
-                int c = 0;
-                for(int j=0; j < 4; j++)
-                        for(int i=0; i < 4; i++)
-                                element(i,j) = mp[c++];
-        }
-    
-        void set_value( real r ) 
-        {
-                for(int i=0; i < 4; i++)
-                        for(int j=0; j < 4; j++)
-                                element(i,j) = r;
-        }
-    
-    void make_identity()
-        {
-                element(0,0) = 1.0;
-                element(0,1) = 0.0;
-                element(0,2) = 0.0; 
-                element(0,3) = 0.0;
-                
-                element(1,0) = 0.0;
-                element(1,1) = 1.0; 
-                element(1,2) = 0.0;
-                element(1,3) = 0.0;
-                
-                element(2,0) = 0.0;
-                element(2,1) = 0.0;
-                element(2,2) = 1.0;
-                element(2,3) = 0.0;
-                
-                element(3,0) = 0.0; 
-                element(3,1) = 0.0; 
-                element(3,2) = 0.0;
-                element(3,3) = 1.0;
-        }
-        
-        
-    static matrix4 identity()
-        {
-                static matrix4 mident (
-                        1.0, 0.0, 0.0, 0.0,
-                        0.0, 1.0, 0.0, 0.0,
-                        0.0, 0.0, 1.0, 0.0,
-                        0.0, 0.0, 0.0, 1.0  );
-                return mident;
-        }
-    
-        
-    void set_scale( real s )
-        {
-                element(0,0) = s;
-                element(1,1) = s;
-                element(2,2) = s;
-        }
-    
-    void set_scale( const vec3 & s )
-        {
-                element(0,0) = s.v[0];
-                element(1,1) = s.v[1];
-                element(2,2) = s.v[2];
-        }
-    
-    
-    void set_translate( const vec3 & t )
-        {
-                element(0,3) = t.v[0];
-                element(1,3) = t.v[1];
-                element(2,3) = t.v[2];
-        }
-    
-        void set_row(int r, const vec4 & t)
-        {
-                element(r,0) = t.v[0];
-                element(r,1) = t.v[1];
-                element(r,2) = t.v[2];
-                element(r,3) = t.v[3];
-        }
-
-        void set_column(int c, const vec4 & t)
-        {
-                element(0,c) = t.v[0];
-                element(1,c) = t.v[1];
-                element(2,c) = t.v[2];
-                element(3,c) = t.v[3];
-        }
-
-    
-        void get_row(int r, vec4 & t) const
-        {
-                t.v[0] = element(r,0);
-                t.v[1] = element(r,1);
-                t.v[2] = element(r,2);
-                t.v[3] = element(r,3);
-        }
-
-        vec4 get_row(int r) const
-        {
-                vec4 v; get_row(r, v);
-                return v;
-        }
-
-        void get_column(int c, vec4 & t) const
-        {
-                t.v[0] = element(0,c);
-                t.v[1] = element(1,c);
-                t.v[2] = element(2,c);
-                t.v[3] = element(3,c);
-        }
-
-        vec4 get_column(int c) const
-        {
-                vec4 v; get_column(c, v);
-                return v;
-        }
-
-    matrix4 inverse() const
-        {
-                matrix4 minv;
-                
-                real r1[8], r2[8], r3[8], r4[8];
-                real *s[4], *tmprow;
-                
-                s[0] = &r1[0];
-                s[1] = &r2[0];
-                s[2] = &r3[0];
-                s[3] = &r4[0];
-                
-                register int i,j,p,jj;
-                for(i=0;i<4;i++)
-                {
-                        for(j=0;j<4;j++)
-                        {
-                                s[i][j] = element(i,j);
-                                if(i==j) s[i][j+4] = 1.0;
-                                else     s[i][j+4] = 0.0;
-                        }
-                }
-                real scp[4];
-                for(i=0;i<4;i++)
-                {
-                        scp[i] = real(fabs(s[i][0]));
-                        for(j=1;j<4;j++)
-                                if(real(fabs(s[i][j])) > scp[i]) scp[i] = real(fabs(s[i][j]));
-                                if(scp[i] == 0.0) return minv; // singular matrix!
-                }
-                
-                int pivot_to;
-                real scp_max;
-                for(i=0;i<4;i++)
-                {
-                        // select pivot row
-                        pivot_to = i;
-                        scp_max = real(fabs(s[i][i]/scp[i]));
-                        // find out which row should be on top
-                        for(p=i+1;p<4;p++)
-                                if(real(fabs(s[p][i]/scp[p])) > scp_max)
-                                { scp_max = real(fabs(s[p][i]/scp[p])); pivot_to = p; }
-                                // Pivot if necessary
-                                if(pivot_to != i)
-                                {
-                                        tmprow = s[i];
-                                        s[i] = s[pivot_to];
-                                        s[pivot_to] = tmprow;
-                                        real tmpscp;
-                                        tmpscp = scp[i];
-                                        scp[i] = scp[pivot_to];
-                                        scp[pivot_to] = tmpscp;
-                                }
-                                
-                                real mji;
-                                // perform gaussian elimination
-                                for(j=i+1;j<4;j++)
-                                {
-                                        mji = s[j][i]/s[i][i];
-                                        s[j][i] = 0.0;
-                                        for(jj=i+1;jj<8;jj++)
-                                                s[j][jj] -= mji*s[i][jj];
-                                }
-                }
-                if(s[3][3] == 0.0) return minv; // singular matrix!
-                
-                //
-                // Now we have an upper triangular matrix.
-                //
-                //  x x x x | y y y y
-                //  0 x x x | y y y y 
-                //  0 0 x x | y y y y
-                //  0 0 0 x | y y y y
-                //
-                //  we'll back substitute to get the inverse
-                //
-                //  1 0 0 0 | z z z z
-                //  0 1 0 0 | z z z z
-                //  0 0 1 0 | z z z z
-                //  0 0 0 1 | z z z z 
-                //
-                
-                real mij;
-                for(i=3;i>0;i--)
-                {
-                        for(j=i-1;j > -1; j--)
-                        {
-                                mij = s[j][i]/s[i][i];
-                                for(jj=j+1;jj<8;jj++)
-                                        s[j][jj] -= mij*s[i][jj];
-                        }
-                }
-                
-                for(i=0;i<4;i++)
-                        for(j=0;j<4;j++)
-                                minv(i,j) = s[i][j+4] / s[i][i];
-                        
-                        return minv;
-        }
-    
-    
-    matrix4 transpose() const
-        {
-                matrix4 mtrans;
-                
-                for(int i=0;i<4;i++)
-                        for(int j=0;j<4;j++)
-                                mtrans(i,j) = element(j,i);             
-                return mtrans;
-        }
-    
-    matrix4 & mult_right( const matrix4 & b )
-        {
-                matrix4 mt(*this);
-                set_value(real(0));
-
-                for(int i=0; i < 4; i++)
-                        for(int j=0; j < 4; j++)
-                                for(int c=0; c < 4; c++)
-                                        element(i,j) += mt(i,c) * b(c,j);
-                return *this;
-        }    
-
-    matrix4 & mult_left( const matrix4 & b )
-        {
-                matrix4 mt(*this);
-                set_value(real(0));
-
-                for(int i=0; i < 4; i++)
-                        for(int j=0; j < 4; j++)
-                                for(int c=0; c < 4; c++)
-                                        element(i,j) += b(i,c) * mt(c,j);
-                return *this;
-        }
-        
-        // dst = M * src
-    void mult_matrix_vec( const vec3 &src, vec3 &dst ) const
-        {
-                real w = (
-                        src.v[0] * element(3,0) +
-                        src.v[1] * element(3,1) + 
-                        src.v[2] * element(3,2) +
-                        element(3,3)          );
-        
-        assert(w != GLH_ZERO);
-
-        dst.v[0]  = (
-                        src.v[0] * element(0,0) +
-                        src.v[1] * element(0,1) +
-                        src.v[2] * element(0,2) +
-                        element(0,3)          ) / w;
-                dst.v[1]  = (
-                        src.v[0] * element(1,0) +
-                        src.v[1] * element(1,1) +
-                        src.v[2] * element(1,2) +
-                        element(1,3)          ) / w;
-                dst.v[2]  = (
-                        src.v[0] * element(2,0) +
-                        src.v[1] * element(2,1) + 
-                        src.v[2] * element(2,2) +
-                        element(2,3)          ) / w;
-        }
-    
-        void mult_matrix_vec( vec3 & src_and_dst) const
-        { mult_matrix_vec(vec3(src_and_dst), src_and_dst); }
-
-
-    // dst = src * M
-    void mult_vec_matrix( const vec3 &src, vec3 &dst ) const
-        {
-                real w = (
-                        src.v[0] * element(0,3) +
-                        src.v[1] * element(1,3) +
-                        src.v[2] * element(2,3) +
-                        element(3,3)          );
-        
-        assert(w != GLH_ZERO);
-
-                dst.v[0]  = (
-                        src.v[0] * element(0,0) +
-                        src.v[1] * element(1,0) + 
-                        src.v[2] * element(2,0) + 
-                        element(3,0)          ) / w;
-                dst.v[1]  = (
-                        src.v[0] * element(0,1) +
-                        src.v[1] * element(1,1) +
-                        src.v[2] * element(2,1) +
-                        element(3,1)          ) / w;
-                dst.v[2]  = (
-                        src.v[0] * element(0,2) +
-                        src.v[1] * element(1,2) +
-                        src.v[2] * element(2,2) +
-                        element(3,2)          ) / w;
-        }
-        
-
-        void mult_vec_matrix( vec3 & src_and_dst) const
-        { mult_vec_matrix(vec3(src_and_dst), src_and_dst); }
-
-        // dst = M * src
-    void mult_matrix_vec( const vec4 &src, vec4 &dst ) const
-        {
-        dst.v[0]  = (
-                        src.v[0] * element(0,0) +
-                        src.v[1] * element(0,1) +
-                        src.v[2] * element(0,2) +
-                        src.v[3] * element(0,3));
-                dst.v[1]  = (
-                        src.v[0] * element(1,0) +
-                        src.v[1] * element(1,1) +
-                        src.v[2] * element(1,2) +
-                        src.v[3] * element(1,3));
-                dst.v[2]  = (
-                        src.v[0] * element(2,0) +
-                        src.v[1] * element(2,1) + 
-                        src.v[2] * element(2,2) +
-                        src.v[3] * element(2,3));
-                dst.v[3] = (
-                        src.v[0] * element(3,0) +
-                        src.v[1] * element(3,1) + 
-                        src.v[2] * element(3,2) +
-                        src.v[3] * element(3,3));
-        }
-    
-        void mult_matrix_vec( vec4 & src_and_dst) const
-        { mult_matrix_vec(vec4(src_and_dst), src_and_dst); }
-
-
-    // dst = src * M
-    void mult_vec_matrix( const vec4 &src, vec4 &dst ) const
-        {
-                dst.v[0]  = (
-                        src.v[0] * element(0,0) +
-                        src.v[1] * element(1,0) + 
-                        src.v[2] * element(2,0) + 
-                        src.v[3] * element(3,0));
-                dst.v[1]  = (
-                        src.v[0] * element(0,1) +
-                        src.v[1] * element(1,1) +
-                        src.v[2] * element(2,1) +
-                        src.v[3] * element(3,1));
-                dst.v[2]  = (
-                        src.v[0] * element(0,2) +
-                        src.v[1] * element(1,2) +
-                        src.v[2] * element(2,2) +
-                        src.v[3] * element(3,2));
-                dst.v[3] = (
-                        src.v[0] * element(0,3) +
-                        src.v[1] * element(1,3) +
-                        src.v[2] * element(2,3) +
-                        src.v[3] * element(3,3));
-        }
-        
-
-        void mult_vec_matrix( vec4 & src_and_dst) const
-        { mult_vec_matrix(vec4(src_and_dst), src_and_dst); }
-
-    
-    // dst = M * src
-    void mult_matrix_dir( const vec3 &src, vec3 &dst ) const
-        {
-                dst.v[0]  = (
-                        src.v[0] * element(0,0) +
-                        src.v[1] * element(0,1) +
-                        src.v[2] * element(0,2) ) ;
-                dst.v[1]  = ( 
-                        src.v[0] * element(1,0) +
-                        src.v[1] * element(1,1) +
-                        src.v[2] * element(1,2) ) ;
-                dst.v[2]  = ( 
-                        src.v[0] * element(2,0) +
-                        src.v[1] * element(2,1) + 
-                        src.v[2] * element(2,2) ) ;
-        }
-        
-
-        void mult_matrix_dir( vec3 & src_and_dst) const
-        { mult_matrix_dir(vec3(src_and_dst), src_and_dst); }
-
-
-        // dst = src * M
-    void mult_dir_matrix( const vec3 &src, vec3 &dst ) const
-        {
-                dst.v[0]  = ( 
-                        src.v[0] * element(0,0) +
-                        src.v[1] * element(1,0) +
-                        src.v[2] * element(2,0) ) ;
-                dst.v[1]  = ( 
-                        src.v[0] * element(0,1) +
-                        src.v[1] * element(1,1) +
-                        src.v[2] * element(2,1) ) ;
-                dst.v[2]  = (
-                        src.v[0] * element(0,2) +
-                        src.v[1] * element(1,2) + 
-                        src.v[2] * element(2,2) ) ;
-        }
-    
-    
-        void mult_dir_matrix( vec3 & src_and_dst) const
-        { mult_dir_matrix(vec3(src_and_dst), src_and_dst); }
-
-
-    real & operator () (int row, int col)
-    { return element(row,col); }
-
-    const real & operator () (int row, int col) const
-    { return element(row,col); }
-
-        real & element (int row, int col)
-    { return m[row | (col<<2)]; }
-
-    const real & element (int row, int col) const
-    { return m[row | (col<<2)]; }
-
-    matrix4 & operator *= ( const matrix4 & mat )
-        {
-                mult_right( mat );
-                return *this;
-        }
-    
-    matrix4 & operator *= ( const real & r )
-        {
-                for (int i = 0; i < 4; ++i)
-        {
-            element(0,i) *= r;
-            element(1,i) *= r;
-            element(2,i) *= r;
-            element(3,i) *= r;
-        }
-                return *this;
-        }
-
-    matrix4 & operator += ( const matrix4 & mat )
-        {
-                for (int i = 0; i < 4; ++i)
-        {
-            element(0,i) += mat.element(0,i);
-            element(1,i) += mat.element(1,i);
-            element(2,i) += mat.element(2,i);
-            element(3,i) += mat.element(3,i);
-        }
-                return *this;
-        }
-
-    friend matrix4 operator * ( const matrix4 & m1,     const matrix4 & m2 );
-    friend bool operator == ( const matrix4 & m1, const matrix4 & m2 );
-    friend bool operator != ( const matrix4 & m1, const matrix4 & m2 );
-    
-  //protected:
-          real m[16];
-  };
-  
-  inline  
-  matrix4 operator * ( const matrix4 & m1, const matrix4 & m2 )
-  {
-          matrix4 product;
-          
-          product = m1;
-          product.mult_right(m2);
-          
-          return product;
-  }
-  
-  inline
-  bool operator ==( const matrix4 &m1, const matrix4 &m2 )
-  {
-          return ( 
-                  m1(0,0) == m2(0,0) &&
-                  m1(0,1) == m2(0,1) &&
-                  m1(0,2) == m2(0,2) &&
-                  m1(0,3) == m2(0,3) &&
-                  m1(1,0) == m2(1,0) &&
-                  m1(1,1) == m2(1,1) &&
-                  m1(1,2) == m2(1,2) &&
-                  m1(1,3) == m2(1,3) &&
-                  m1(2,0) == m2(2,0) &&
-                  m1(2,1) == m2(2,1) &&
-                  m1(2,2) == m2(2,2) &&
-                  m1(2,3) == m2(2,3) &&
-                  m1(3,0) == m2(3,0) &&
-                  m1(3,1) == m2(3,1) &&
-                  m1(3,2) == m2(3,2) &&
-                  m1(3,3) == m2(3,3) );
-  }
-  
-  inline
-  bool operator != ( const matrix4 & m1, const matrix4 & m2 )
-  { return !( m1 == m2 ); }  
-
-
-
-
-
-
-
-
-
-
-
-
-  
-    class quaternion
-    {
-    public:
-    
-    quaternion()
-    {
-        *this = identity();
-    }
-
-    quaternion( const real v[4] )
-    {
-        set_value( v );
-    }
-
-
-    quaternion( real q0, real q1, real q2, real q3 )
-    {
-        set_value( q0, q1, q2, q3 );
-    }
-
-
-    quaternion( const matrix4 & m )
-    {
-        set_value( m );
-    }
-
-
-    quaternion( const vec3 &axis, real radians )
-    {
-        set_value( axis, radians );
-    }
-
-
-    quaternion( const vec3 &rotateFrom, const vec3 &rotateTo )
-    {
-        set_value( rotateFrom, rotateTo );
-    }
-
-    quaternion( const vec3 & from_look, const vec3 & from_up,
-                      const vec3 & to_look, const vec3& to_up)
-    {
-            set_value(from_look, from_up, to_look, to_up);
-    }
-
-    const real * get_value() const
-    {
-        return  &q[0];
-    }
-
-    void get_value( real &q0, real &q1, real &q2, real &q3 ) const
-    {
-        q0 = q[0];
-        q1 = q[1];
-        q2 = q[2];
-        q3 = q[3];
-    }
-
-    quaternion & set_value( real q0, real q1, real q2, real q3 )
-    {
-        q[0] = q0;
-        q[1] = q1;
-        q[2] = q2;
-        q[3] = q3;
-        counter = 0;
-        return *this;
-    }
-
-    void get_value( vec3 &axis, real &radians ) const
-    {
-        radians = real(acos( q[3] ) * GLH_TWO);
-        if ( radians == GLH_ZERO )
-            axis = vec3( 0.0, 0.0, 1.0 );
-        else
-        {
-            axis.v[0] = q[0];
-            axis.v[1] = q[1];
-            axis.v[2] = q[2];
-            axis.normalize();
-        }
-    }
-
-    void get_value( matrix4 & m ) const
-    {
-        real s, xs, ys, zs, wx, wy, wz, xx, xy, xz, yy, yz, zz;
-
-        real norm = q[0] * q[0] + q[1] * q[1] + q[2] * q[2] + q[3] * q[3];
-
-        s = (equivalent(norm,GLH_ZERO)) ? GLH_ZERO : ( GLH_TWO / norm );
-
-        xs = q[0] * s;
-        ys = q[1] * s;
-        zs = q[2] * s;
-
-        wx = q[3] * xs;
-        wy = q[3] * ys;
-        wz = q[3] * zs;
-
-        xx = q[0] * xs;
-        xy = q[0] * ys;
-        xz = q[0] * zs;
-
-        yy = q[1] * ys;
-        yz = q[1] * zs;
-        zz = q[2] * zs;
-
-        m(0,0) = real( GLH_ONE - ( yy + zz ));
-        m(1,0) = real ( xy + wz );
-        m(2,0) = real ( xz - wy );
-
-        m(0,1) = real ( xy - wz );
-        m(1,1) = real ( GLH_ONE - ( xx + zz ));
-        m(2,1) = real ( yz + wx );
-
-        m(0,2) = real ( xz + wy );
-        m(1,2) = real ( yz - wx );
-        m(2,2) = real ( GLH_ONE - ( xx + yy ));
-
-        m(3,0) = m(3,1) = m(3,2) = m(0,3) = m(1,3) = m(2,3) = GLH_ZERO;
-        m(3,3) = GLH_ONE;
-    }
-
-    quaternion & set_value( const real * qp )
-    {
-        memcpy(q,qp,sizeof(real) * 4);
-
-        counter = 0;
-        return *this;
-    }
-
-    quaternion & set_value( const matrix4 & m )
-    {
-        real tr, s;
-        int i, j, k;
-        const int nxt[3] = { 1, 2, 0 };
-
-        tr = m(0,0) + m(1,1) + m(2,2);
-
-        if ( tr > GLH_ZERO )
-        {
-            s = real(sqrt( tr + m(3,3) ));
-            q[3] = real ( s * 0.5 );
-            s = real(0.5) / s;
-
-            q[0] = real ( ( m(1,2) - m(2,1) ) * s );
-            q[1] = real ( ( m(2,0) - m(0,2) ) * s );
-            q[2] = real ( ( m(0,1) - m(1,0) ) * s );
-        }
-        else
-        {
-            i = 0;
-            if ( m(1,1) > m(0,0) )
-              i = 1;
-
-            if ( m(2,2) > m(i,i) )
-              i = 2;
-
-            j = nxt[i];
-            k = nxt[j];
-
-            s = real(sqrt( ( m(i,j) - ( m(j,j) + m(k,k) )) + GLH_ONE ));
-
-            q[i] = real ( s * 0.5 );
-            s = real(0.5 / s);
-
-            q[3] = real ( ( m(j,k) - m(k,j) ) * s );
-            q[j] = real ( ( m(i,j) + m(j,i) ) * s );
-            q[k] = real ( ( m(i,k) + m(k,i) ) * s );
-        }
-
-        counter = 0;
-        return *this;
-    }
-
-    quaternion & set_value( const vec3 &axis, real theta )
-    {
-        real sqnorm = axis.square_norm();
-
-        if (sqnorm <= GLH_EPSILON)
-        {
-            // axis too small.
-            x = y = z = 0.0;
-            w = 1.0;
-        } 
-        else 
-        {
-            theta *= real(0.5);
-            real sin_theta = real(sin(theta));
-
-            if (!equivalent(sqnorm,GLH_ONE)) 
-              sin_theta /= real(sqrt(sqnorm));
-            x = sin_theta * axis.v[0];
-            y = sin_theta * axis.v[1];
-            z = sin_theta * axis.v[2];
-            w = real(cos(theta));
-        }
-        return *this;
-    }
-
-    quaternion & set_value( const vec3 & rotateFrom, const vec3 & rotateTo )
-    {
-        vec3 p1, p2;
-        real alpha;
-
-        p1 = rotateFrom; 
-        p1.normalize();
-        p2 = rotateTo;  
-        p2.normalize();
-
-        alpha = p1.dot(p2);
-
-        if(equivalent(alpha,GLH_ONE))
-        { 
-            *this = identity(); 
-            return *this; 
-        }
-
-        // ensures that the anti-parallel case leads to a positive dot
-        if(equivalent(alpha,-GLH_ONE))
-        {
-            vec3 v;
-
-            if(p1.v[0] != p1.v[1] || p1.v[0] != p1.v[2])
-                v = vec3(p1.v[1], p1.v[2], p1.v[0]);
-            else
-                v = vec3(-p1.v[0], p1.v[1], p1.v[2]);
-
-            v -= p1 * p1.dot(v);
-            v.normalize();
-
-            set_value(v, GLH_PI);
-            return *this;
-        }
-
-        p1 = p1.cross(p2);  
-        p1.normalize();
-        set_value(p1,real(acos(alpha)));
-
-        counter = 0;
-        return *this;
-    }
-
-    quaternion & set_value( const vec3 & from_look, const vec3 & from_up,
-                      const vec3 & to_look, const vec3 & to_up)
-    {
-            quaternion r_look = quaternion(from_look, to_look);
-            
-            vec3 rotated_from_up(from_up);
-            r_look.mult_vec(rotated_from_up);
-            
-            quaternion r_twist = quaternion(rotated_from_up, to_up);
-            
-            *this = r_twist;
-            *this *= r_look;
-            return *this;
-    }
-
-    quaternion & operator *= ( const quaternion & qr )
-    {
-        quaternion ql(*this);
-   
-        w = ql.w * qr.w - ql.x * qr.x - ql.y * qr.y - ql.z * qr.z;
-        x = ql.w * qr.x + ql.x * qr.w + ql.y * qr.z - ql.z * qr.y;
-        y = ql.w * qr.y + ql.y * qr.w + ql.z * qr.x - ql.x * qr.z;
-        z = ql.w * qr.z + ql.z * qr.w + ql.x * qr.y - ql.y * qr.x;
-
-        counter += qr.counter;
-        counter++;
-        counter_normalize();
-        return *this;
-    }
-
-    void normalize()
-    {
-        real rnorm = GLH_ONE / real(sqrt(w * w + x * x + y * y + z * z));
-        if (equivalent(rnorm, GLH_ZERO))
-            return;
-        x *= rnorm;
-        y *= rnorm;
-        z *= rnorm;
-        w *= rnorm;
-        counter = 0;
-    }
-
-    friend bool operator == ( const quaternion & q1, const quaternion & q2 );      
-
-    friend bool operator != ( const quaternion & q1, const quaternion & q2 );
-
-    friend quaternion operator * ( const quaternion & q1, const quaternion & q2 );
-
-    bool equals( const quaternion & r, real tolerance ) const
-    {
-        real t;
-
-        t = (
-                        (q[0]-r.q[0])*(q[0]-r.q[0]) +
-            (q[1]-r.q[1])*(q[1]-r.q[1]) +
-            (q[2]-r.q[2])*(q[2]-r.q[2]) +
-            (q[3]-r.q[3])*(q[3]-r.q[3]) );
-        if(t > GLH_EPSILON) 
-            return false;
-        return 1;
-    }
-
-    quaternion & conjugate()
-    {
-        q[0] *= -GLH_ONE;
-        q[1] *= -GLH_ONE;
-        q[2] *= -GLH_ONE;
-        return *this;
-    }
-
-    quaternion & invert()
-    {
-        return conjugate();
-    }
-
-    quaternion inverse() const
-    {
-        quaternion r = *this;
-        return r.invert();
-    }
-
-    //
-    // Quaternion multiplication with cartesian vector
-    // v' = q*v*q(star)
-    //
-    void mult_vec( const vec3 &src, vec3 &dst ) const
-    {
-        real v_coef = w * w - x * x - y * y - z * z;                     
-        real u_coef = GLH_TWO * (src.v[0] * x + src.v[1] * y + src.v[2] * z);  
-        real c_coef = GLH_TWO * w;                                       
-
-        dst.v[0] = v_coef * src.v[0] + u_coef * x + c_coef * (y * src.v[2] - z * src.v[1]);
-        dst.v[1] = v_coef * src.v[1] + u_coef * y + c_coef * (z * src.v[0] - x * src.v[2]);
-        dst.v[2] = v_coef * src.v[2] + u_coef * z + c_coef * (x * src.v[1] - y * src.v[0]);
-    }
-
-    void mult_vec( vec3 & src_and_dst) const
-    {
-        mult_vec(vec3(src_and_dst), src_and_dst);
-    }
-
-    void scale_angle( real scaleFactor )
-    {
-        vec3 axis;
-        real radians;
-
-        get_value(axis, radians);
-        radians *= scaleFactor;
-        set_value(axis, radians);
-    }
-
-    static quaternion slerp( const quaternion & p, const quaternion & q, real alpha )
-    {
-        quaternion r;
-
-        real cos_omega = p.x * q.x + p.y * q.y + p.z * q.z + p.w * q.w;
-        // if B is on opposite hemisphere from A, use -B instead
-      
-        int bflip;
-        if ( ( bflip = (cos_omega < GLH_ZERO)) )
-            cos_omega = -cos_omega;
-
-        // complementary interpolation parameter
-        real beta = GLH_ONE - alpha;     
-
-        if(cos_omega <= GLH_ONE - GLH_EPSILON)
-            return p;
-
-        real omega = real(acos(cos_omega));
-        real one_over_sin_omega = GLH_ONE / real(sin(omega));
-
-        beta    = real(sin(omega*beta)  * one_over_sin_omega);
-        alpha   = real(sin(omega*alpha) * one_over_sin_omega);
-
-        if (bflip)
-            alpha = -alpha;
-
-        r.x = beta * p.q[0]+ alpha * q.q[0];
-        r.y = beta * p.q[1]+ alpha * q.q[1];
-        r.z = beta * p.q[2]+ alpha * q.q[2];
-        r.w = beta * p.q[3]+ alpha * q.q[3];
-        return r;
-    }
-
-    static quaternion identity()
-    {
-        static quaternion ident( vec3( 0.0, 0.0, 0.0 ), GLH_ONE );
-        return ident;
-    }
-
-    real & operator []( int i )
-    {
-        assert(i < 4);
-        return q[i];
-    }
-
-    const real & operator []( int i ) const
-    {
-        assert(i < 4);
-        return q[i];
-    }
-
-    protected:
-
-        void counter_normalize()
-        {
-            if (counter > GLH_QUATERNION_NORMALIZATION_THRESHOLD)
-                normalize();
-        }
-
-        union 
-        {
-            struct 
-            {
-                real q[4];
-            };
-            struct 
-            {
-                real x;
-                real y;
-                real z;
-                real w;
-            };
-        };
-
-        // renormalization counter
-        unsigned char counter;
-    };
-
-    inline
-    bool operator == ( const quaternion & q1, const quaternion & q2 )
-    {
-        return (equivalent(q1.x, q2.x) &&
-                        equivalent(q1.y, q2.y) &&
-                        equivalent(q1.z, q2.z) &&
-                        equivalent(q1.w, q2.w) );
-    }
-
-    inline
-    bool operator != ( const quaternion & q1, const quaternion & q2 )
-    { 
-        return ! ( q1 == q2 ); 
-    }
-
-    inline
-    quaternion operator * ( const quaternion & q1, const quaternion & q2 )
-    {   
-        quaternion r(q1); 
-        r *= q2; 
-        return r; 
-    }
-  
-      
-    
-
-
-
-
-
-  
-  
-  class plane
-  {
-  public:
-          
-          plane()
-      {
-                  planedistance = 0.0;
-                  planenormal.set_value( 0.0, 0.0, 1.0 );
-      }
-          
-          
-          plane( const vec3 &p0, const vec3 &p1, const vec3 &p2 )
-      {
-                  vec3 v0 = p1 - p0;
-                  vec3 v1 = p2 - p0;
-                  planenormal = v0.cross(v1);  
-                  planenormal.normalize();
-                  planedistance = p0.dot(planenormal);
-      }
-          
-          plane( const vec3 &normal, real distance )
-      {
-                  planedistance = distance;
-                  planenormal = normal;
-                  planenormal.normalize();
-      }
-          
-          plane( const vec3 &normal, const vec3 &point )
-      {
-                  planenormal = normal;
-                  planenormal.normalize();
-                  planedistance = point.dot(planenormal);
-      }
-          
-          void offset( real d )
-      {
-                  planedistance += d;
-      }
-          
-          bool intersect( const line &l, vec3 &intersection ) const
-      {
-                  vec3 pos, dir;
-                  vec3 pn = planenormal;
-                  real pd = planedistance;
-                  
-                  pos = l.get_position();
-                  dir = l.get_direction();
-                  
-                  if(dir.dot(pn) == 0.0) return 0;
-                  pos -= pn*pd;
-                  // now we're talking about a plane passing through the origin
-                  if(pos.dot(pn) < 0.0) pn.negate();
-                  if(dir.dot(pn) > 0.0) dir.negate();
-                  vec3 ppos = pn * pos.dot(pn);
-                  pos = (ppos.length()/dir.dot(-pn))*dir;
-                  intersection = l.get_position();
-                  intersection += pos;
-                  return 1;
-      }
-          void transform( const matrix4 &matrix )
-      {
-                  matrix4 invtr = matrix.inverse();
-                  invtr = invtr.transpose();
-                  
-                  vec3 pntOnplane = planenormal * planedistance;
-                  vec3 newPntOnplane;
-                  vec3 newnormal;
-                  
-                  invtr.mult_dir_matrix(planenormal, newnormal);
-                  matrix.mult_vec_matrix(pntOnplane, newPntOnplane);
-                  
-                  newnormal.normalize();
-                  planenormal = newnormal;
-                  planedistance = newPntOnplane.dot(planenormal);
-      }
-          
-          bool is_in_half_space( const vec3 &point ) const
-      {
-                  
-                  if(( point.dot(planenormal) - planedistance) < 0.0)
-                          return 0;
-                  return 1;
-      }
-          
-          
-          real distance( const vec3 & point ) const 
-      {
-                  return planenormal.dot(point - planenormal*planedistance);
-      }
-          
-          const vec3 &get_normal() const
-      {
-                  return planenormal;
-      }
-          
-          
-          real get_distance_from_origin() const
-      {
-                  return planedistance;
-      }
-          
-          
-          friend bool operator == ( const plane & p1, const plane & p2 );
-          
-          
-          friend bool operator != ( const plane & p1, const plane & p2 );
-          
-  //protected:
-          vec3 planenormal;
-          real planedistance;
-  };
-  
-  inline
-  bool operator == (const plane & p1, const plane & p2 )
-  {
-          return (  p1.planedistance == p2.planedistance && p1.planenormal == p2.planenormal);
-  }
-  
-  inline
-  bool operator != ( const plane & p1, const plane & p2 )
-  { return  ! (p1 == p2); }
-  
-  
-
-  } // "ns_##GLH_REAL"
-
-  // make common typedefs...
-#ifdef GLH_REAL_IS_FLOAT
-  typedef GLH_REAL_NAMESPACE::vec2 vec2f;
-  typedef GLH_REAL_NAMESPACE::vec3 vec3f;
-  typedef GLH_REAL_NAMESPACE::vec4 vec4f;
-  typedef GLH_REAL_NAMESPACE::quaternion quaternionf;
-  typedef GLH_REAL_NAMESPACE::quaternion rotationf;
-  typedef GLH_REAL_NAMESPACE::line linef;
-  typedef GLH_REAL_NAMESPACE::plane planef;
-  typedef GLH_REAL_NAMESPACE::matrix4 matrix4f;
-#endif
-
-  
-
-
-}  // namespace glh
-
-
-
-#endif
-
+/*

+    glh - is a platform-indepenedent C++ OpenGL helper library 

+

+

+    Copyright (c) 2000 Cass Everitt

+        Copyright (c) 2000 NVIDIA Corporation

+    All rights reserved.

+

+    Redistribution and use in source and binary forms, with or

+        without modification, are permitted provided that the following

+        conditions are met:

+

+     * Redistributions of source code must retain the above

+           copyright notice, this list of conditions and the following

+           disclaimer.

+

+     * Redistributions in binary form must reproduce the above

+           copyright notice, this list of conditions and the following

+           disclaimer in the documentation and/or other materials

+           provided with the distribution.

+

+     * The names of contributors to this software may not be used

+           to endorse or promote products derived from this software

+           without specific prior written permission. 

+

+       THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS

+           ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT

+           LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS

+           FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE

+           REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,

+           INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,

+           BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;

+           LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER

+           CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT

+           LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN

+           ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 

+           POSSIBILITY OF SUCH DAMAGE. 

+

+

+    Cass Everitt - cass@r3.nu

+*/

+

+/*

+glh_linear.h

+*/

+

+// Author:  Cass W. Everitt

+

+#ifndef GLH_LINEAR_H

+#define GLH_LINEAR_H

+

+#include <memory.h>

+#include <math.h>

+#include <assert.h>

+

+// only supports float for now...

+#define GLH_REAL_IS_FLOAT

+

+#ifdef GLH_REAL_IS_FLOAT

+# define GLH_REAL float

+# define GLH_REAL_NAMESPACE ns_float

+#endif

+

+#define     GLH_QUATERNION_NORMALIZATION_THRESHOLD  64

+

+#define     GLH_RAD_TO_DEG      GLH_REAL(57.2957795130823208767981548141052)

+#define     GLH_DEG_TO_RAD      GLH_REAL(0.0174532925199432957692369076848861)

+#define     GLH_ZERO            GLH_REAL(0.0)

+#define     GLH_ONE             GLH_REAL(1.0)

+#define     GLH_TWO             GLH_REAL(2.0)

+#define     GLH_EPSILON         GLH_REAL(10e-6)

+#define     GLH_PI              GLH_REAL(3.1415926535897932384626433832795)    

+

+#define     equivalent(a,b)     (((a < b + GLH_EPSILON) && (a > b - GLH_EPSILON)) ? true : false)

+

+namespace glh

+{

+

+        inline GLH_REAL to_degrees(GLH_REAL radians) { return radians*GLH_RAD_TO_DEG; }

+        inline GLH_REAL to_radians(GLH_REAL degrees) { return degrees*GLH_DEG_TO_RAD; }

+

+        // forward declarations for friend template functions.

+        template <int N, class T> class vec;

+        

+        // forward declarations for friend template functions.

+        template <int N, class T>

+                bool operator == ( const vec<N,T> & v1, const vec<N,T> & v2 );

+        

+        // forward declarations for friend template functions.

+        template <int N, class T>

+                bool operator != ( const vec<N,T> & v1, const vec<N,T> & v2 );

+

+        template <int N, class T>       

+        class vec

+        {                               

+    public:

+                int size() const { return N; }

+                

+                vec(const T & t = T()) 

+                { for(int i = 0; i < N; i++) v[i] = t; }

+                vec(const T * tp)

+                { for(int i = 0; i < N; i++) v[i] = tp[i]; }

+                

+                const T * get_value() const

+                { return v; }

+                

+                

+                T dot( const vec<N,T> & rhs ) const

+                { 

+                        T r = 0;

+                        for(int i = 0; i < N; i++) r += v[i]*rhs.v[i];

+                        return r;

+                }

+                

+                T length() const

+                {

+                        T r = 0;

+                        for(int i = 0; i < N; i++) r += v[i]*v[i]; 

+                        return T(sqrt(r));

+                }       

+                

+                T square_norm() const

+                {

+                        T r = 0;

+                        for(int i = 0; i < N; i++) r += v[i]*v[i]; 

+                        return r;

+                }       

+                

+                void  negate()

+                { for(int i = 0; i < N; i++) v[i] = -v[i]; }

+                

+                

+                T normalize() 

+                { 

+                        T sum(0);

+                        for(int i = 0; i < N; i++) 

+                sum += v[i]*v[i];

+                        sum = T(sqrt(sum));

+            if (sum > GLH_EPSILON)

+                            for(int i = 0; i < N; i++) 

+                    v[i] /= sum;

+                        return sum;

+                }

+                

+                

+                vec<N,T> & set_value( const T * rhs )

+                { for(int i = 0; i < N; i++) v[i] = rhs[i]; return *this; }

+                

+                T & operator [] ( int i )

+                { return v[i]; }

+                

+                const T & operator [] ( int i ) const

+                { return v[i]; }

+

+                vec<N,T> & operator *= ( T d )

+                { for(int i = 0; i < N; i++) v[i] *= d; return *this;}

+                

+                vec<N,T> & operator *= ( const vec<N,T> & u )

+                { for(int i = 0; i < N; i++) v[i] *= u[i]; return *this;}

+                

+                vec<N,T> & operator /= ( T d )

+                { if(d == 0) return *this; for(int i = 0; i < N; i++) v[i] /= d; return *this;}

+                

+                vec<N,T> & operator += ( const vec<N,T> & u )

+                { for(int i = 0; i < N; i++) v[i] += u.v[i]; return *this;}

+                

+                vec<N,T> & operator -= ( const vec<N,T> & u )

+                { for(int i = 0; i < N; i++) v[i] -= u.v[i]; return *this;}

+                

+                

+                vec<N,T> operator - () const

+                { vec<N,T> rv = v; rv.negate(); return rv; }

+                

+                vec<N,T> operator + ( const vec<N,T> &v) const

+                { vec<N,T> rt(*this); return rt += v; }

+                

+                vec<N,T> operator - ( const vec<N,T> &v) const

+                { vec<N,T> rt(*this); return rt -= v; }

+                

+                vec<N,T> operator * ( T d) const

+                { vec<N,T> rt(*this); return rt *= d; }

+                

+                friend bool operator == <> ( const vec<N,T> &v1, const vec<N,T> &v2 );

+                friend bool operator != <> ( const vec<N,T> &v1, const vec<N,T> &v2 );

+                

+                

+        //protected:

+                T v[N];

+        };

+        

+        

+        

+        // vector friend operators

+        

+        template <int N, class T> inline

+                vec<N,T> operator * ( const vec<N,T> & b, T d )

+        {

+                vec<N,T> rt(b);

+                return rt *= d;

+        }

+

+        template <int N, class T> inline

+                vec<N,T> operator * ( T d, const vec<N,T> & b )

+        { return b*d; }

+        

+        template <int N, class T> inline

+                vec<N,T> operator * ( const vec<N,T> & b, const vec<N,T> & d )

+        {

+                vec<N,T> rt(b);

+                return rt *= d;

+        }

+

+        template <int N, class T> inline

+                vec<N,T> operator / ( const vec<N,T> & b, T d )

+        { vec<N,T> rt(b); return rt /= d; }

+        

+        template <int N, class T> inline

+                vec<N,T> operator + ( const vec<N,T> & v1, const vec<N,T> & v2 )

+        { vec<N,T> rt(v1); return rt += v2; }

+        

+        template <int N, class T> inline

+                vec<N,T> operator - ( const vec<N,T> & v1, const vec<N,T> & v2 )

+        { vec<N,T> rt(v1); return rt -= v2; }

+        

+        

+        template <int N, class T> inline

+                bool operator == ( const vec<N,T> & v1, const vec<N,T> & v2 )

+        {

+                for(int i = 0; i < N; i++)

+                        if(v1.v[i] != v2.v[i])

+                                return false;

+                        return true;

+        }

+        

+        template <int N, class T> inline

+                bool operator != ( const vec<N,T> & v1, const vec<N,T> & v2 )

+        { return !(v1 == v2); }

+        

+

+        typedef vec<3,unsigned char> vec3ub;

+        typedef vec<4,unsigned char> vec4ub;

+

+

+

+

+

+        namespace GLH_REAL_NAMESPACE

+        {

+        typedef GLH_REAL real;

+

+        class line;

+        class plane;

+        class matrix4;

+        class quaternion;

+        typedef quaternion rotation; 

+  

+        class vec2 : public vec<2,real>

+        {

+    public:

+                vec2(const real & t = real()) : vec<2,real>(t)

+                {}

+                vec2(const vec<2,real> & t) : vec<2,real>(t)

+                {}

+                vec2(const real * tp) : vec<2,real>(tp)

+                {}

+                

+                vec2(real x, real y )

+                { v[0] = x; v[1] = y; }

+                

+                void get_value(real & x, real & y) const

+                { x = v[0]; y = v[1]; }

+                

+                vec2 & set_value( const real & x, const real & y)

+                { v[0] = x; v[1] = y; return *this; }

+                

+        };

+        

+        

+        class vec3 : public vec<3,real>

+        {

+    public:

+                vec3(const real & t = real()) : vec<3,real>(t)

+                {}

+                vec3(const vec<3,real> & t) : vec<3,real>(t)

+                {}

+                vec3(const real * tp) : vec<3,real>(tp)

+                {}

+                

+                vec3(real x, real y, real z)

+                { v[0] = x; v[1] = y; v[2] = z; }

+                

+                void get_value(real & x, real & y, real & z) const

+                { x = v[0]; y = v[1]; z = v[2]; }

+                

+                vec3 cross( const vec3 &rhs ) const

+                {

+                        vec3 rt;

+                        rt.v[0] = v[1]*rhs.v[2]-v[2]*rhs.v[1];

+                        rt.v[1] = v[2]*rhs.v[0]-v[0]*rhs.v[2];

+                        rt.v[2] = v[0]*rhs.v[1]-v[1]*rhs.v[0];  

+                        return rt;

+                }

+                

+                vec3 & set_value( const real & x, const real & y, const real & z)

+                { v[0] = x; v[1] = y; v[2] = z; return *this; }

+                

+        };

+

+                

+    class vec4 : public vec<4,real>

+    {

+    public:

+        vec4(const real & t = real()) : vec<4,real>(t)

+        {}

+        vec4(const vec<4,real> & t) : vec<4,real>(t)

+        {}

+

+        vec4(const vec<3,real> & t, real fourth)

+

+        { v[0] = t.v[0]; v[1] = t.v[1]; v[2] = t.v[2]; v[3] = fourth; }

+        vec4(const real * tp) : vec<4,real>(tp)

+        {}

+        vec4(real x, real y, real z, real w)

+        { v[0] = x; v[1] = y; v[2] = z; v[3] = w; }

+

+        void get_value(real & x, real & y, real & z, real & w) const

+        { x = v[0]; y = v[1]; z = v[2]; w = v[3]; }

+  

+        vec4 & set_value( const real & x, const real & y, const real & z, const real & w)

+        { v[0] = x; v[1] = y; v[2] = z; v[3] = w; return *this; }

+    };

+

+    inline

+    vec3 homogenize(const vec4 & v)

+    {

+      vec3 rt;

+      assert(v.v[3] != GLH_ZERO);

+      rt.v[0] = v.v[0]/v.v[3];

+      rt.v[1] = v.v[1]/v.v[3];

+      rt.v[2] = v.v[2]/v.v[3];

+      return rt;

+    }

+  

+

+

+    class line

+    {

+    public:

+  

+        line()

+        { set_value(vec3(0,0,0),vec3(0,0,1)); }

+

+        line( const vec3 & p0, const vec3 &p1)

+        { set_value(p0,p1); }

+

+        void set_value( const vec3 &p0, const vec3 &p1)

+        {

+          position = p0;

+          direction = p1-p0;

+          direction.normalize();

+        }

+  

+        bool get_closest_points(const line &line2, 

+                                                  vec3 &pointOnThis,

+                                                  vec3 &pointOnThat)

+        {

+  

+          // quick check to see if parallel -- if so, quit.

+          if(fabs(direction.dot(line2.direction)) == 1.0)

+                  return 0;

+          line l2 = line2;

+  

+          // Algorithm: Brian Jean

+          // 

+          register real u;

+          register real v;

+          vec3 Vr = direction;

+          vec3 Vs = l2.direction;

+          register real Vr_Dot_Vs = Vr.dot(Vs);

+          register real detA = real(1.0 - (Vr_Dot_Vs * Vr_Dot_Vs));

+          vec3 C = l2.position - position;

+          register real C_Dot_Vr =  C.dot(Vr);

+          register real C_Dot_Vs =  C.dot(Vs);

+  

+          u = (C_Dot_Vr - Vr_Dot_Vs * C_Dot_Vs)/detA;

+          v = (C_Dot_Vr * Vr_Dot_Vs - C_Dot_Vs)/detA;

+  

+          pointOnThis = position;

+          pointOnThis += direction * u;

+          pointOnThat = l2.position;

+          pointOnThat += l2.direction * v;

+  

+          return 1;

+        }

+  

+        vec3 get_closest_point(const vec3 &point)

+        {

+          vec3 np = point - position;

+          vec3 rp = direction*direction.dot(np)+position;

+          return rp;

+        }

+  

+        const vec3 & get_position() const {return position;}

+

+        const vec3 & get_direction() const {return direction;}

+  

+    //protected:

+        vec3 position;

+        vec3 direction;

+    };

+  

+  

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+  

+  

+

+

+

+

+

+

+

+

+

+

+

+  // matrix

+

+  

+  class matrix4

+  {

+    

+  public:

+        

+    matrix4() { make_identity(); }

+    

+        matrix4( real r ) 

+        { set_value(r); }

+

+        matrix4( real * m )

+        { set_value(m); }

+    

+    matrix4( real a00, real a01, real a02, real a03,

+               real a10, real a11, real a12, real a13,

+                   real a20, real a21, real a22, real a23,

+                   real a30, real a31, real a32, real a33 )

+        {

+                element(0,0) = a00;

+                element(0,1) = a01;

+                element(0,2) = a02;

+                element(0,3) = a03;

+                

+                element(1,0) = a10;

+                element(1,1) = a11;

+                element(1,2) = a12;

+                element(1,3) = a13;

+                

+                element(2,0) = a20;

+                element(2,1) = a21;

+                element(2,2) = a22;

+                element(2,3) = a23;

+                

+                element(3,0) = a30;

+                element(3,1) = a31;

+                element(3,2) = a32;

+                element(3,3) = a33;

+        }

+            

+    

+    void get_value( real * mp ) const

+        {

+                int c = 0;

+                for(int j=0; j < 4; j++)

+                        for(int i=0; i < 4; i++)

+                                mp[c++] = element(i,j);

+        }

+    

+    

+    const real * get_value() const

+        { return m; }

+    

+        void set_value( real * mp)

+        {

+                int c = 0;

+                for(int j=0; j < 4; j++)

+                        for(int i=0; i < 4; i++)

+                                element(i,j) = mp[c++];

+        }

+    

+        void set_value( real r ) 

+        {

+                for(int i=0; i < 4; i++)

+                        for(int j=0; j < 4; j++)

+                                element(i,j) = r;

+        }

+    

+    void make_identity()

+        {

+                element(0,0) = 1.0;

+                element(0,1) = 0.0;

+                element(0,2) = 0.0; 

+                element(0,3) = 0.0;

+                

+                element(1,0) = 0.0;

+                element(1,1) = 1.0; 

+                element(1,2) = 0.0;

+                element(1,3) = 0.0;

+                

+                element(2,0) = 0.0;

+                element(2,1) = 0.0;

+                element(2,2) = 1.0;

+                element(2,3) = 0.0;

+                

+                element(3,0) = 0.0; 

+                element(3,1) = 0.0; 

+                element(3,2) = 0.0;

+                element(3,3) = 1.0;

+        }

+        

+        

+    static matrix4 identity()

+        {

+                static matrix4 mident (

+                        1.0, 0.0, 0.0, 0.0,

+                        0.0, 1.0, 0.0, 0.0,

+                        0.0, 0.0, 1.0, 0.0,

+                        0.0, 0.0, 0.0, 1.0  );

+                return mident;

+        }

+    

+        

+    void set_scale( real s )

+        {

+                element(0,0) = s;

+                element(1,1) = s;

+                element(2,2) = s;

+        }

+    

+    void set_scale( const vec3 & s )

+        {

+                element(0,0) = s.v[0];

+                element(1,1) = s.v[1];

+                element(2,2) = s.v[2];

+        }

+    

+    

+    void set_translate( const vec3 & t )

+        {

+                element(0,3) = t.v[0];

+                element(1,3) = t.v[1];

+                element(2,3) = t.v[2];

+        }

+    

+        void set_row(int r, const vec4 & t)

+        {

+                element(r,0) = t.v[0];

+                element(r,1) = t.v[1];

+                element(r,2) = t.v[2];

+                element(r,3) = t.v[3];

+        }

+

+        void set_column(int c, const vec4 & t)

+        {

+                element(0,c) = t.v[0];

+                element(1,c) = t.v[1];

+                element(2,c) = t.v[2];

+                element(3,c) = t.v[3];

+        }

+

+    

+        void get_row(int r, vec4 & t) const

+        {

+                t.v[0] = element(r,0);

+                t.v[1] = element(r,1);

+                t.v[2] = element(r,2);

+                t.v[3] = element(r,3);

+        }

+

+        vec4 get_row(int r) const

+        {

+                vec4 v; get_row(r, v);

+                return v;

+        }

+

+        void get_column(int c, vec4 & t) const

+        {

+                t.v[0] = element(0,c);

+                t.v[1] = element(1,c);

+                t.v[2] = element(2,c);

+                t.v[3] = element(3,c);

+        }

+

+        vec4 get_column(int c) const

+        {

+                vec4 v; get_column(c, v);

+                return v;

+        }

+

+    matrix4 inverse() const

+        {

+                matrix4 minv;

+                

+                real r1[8], r2[8], r3[8], r4[8];

+                real *s[4], *tmprow;

+                

+                s[0] = &r1[0];

+                s[1] = &r2[0];

+                s[2] = &r3[0];

+                s[3] = &r4[0];

+                

+                register int i,j,p,jj;

+                for(i=0;i<4;i++)

+                {

+                        for(j=0;j<4;j++)

+                        {

+                                s[i][j] = element(i,j);

+                                if(i==j) s[i][j+4] = 1.0;

+                                else     s[i][j+4] = 0.0;

+                        }

+                }

+                real scp[4];

+                for(i=0;i<4;i++)

+                {

+                        scp[i] = real(fabs(s[i][0]));

+                        for(j=1;j<4;j++)

+                                if(real(fabs(s[i][j])) > scp[i]) scp[i] = real(fabs(s[i][j]));

+                                if(scp[i] == 0.0) return minv; // singular matrix!

+                }

+                

+                int pivot_to;

+                real scp_max;

+                for(i=0;i<4;i++)

+                {

+                        // select pivot row

+                        pivot_to = i;

+                        scp_max = real(fabs(s[i][i]/scp[i]));

+                        // find out which row should be on top

+                        for(p=i+1;p<4;p++)

+                                if(real(fabs(s[p][i]/scp[p])) > scp_max)

+                                { scp_max = real(fabs(s[p][i]/scp[p])); pivot_to = p; }

+                                // Pivot if necessary

+                                if(pivot_to != i)

+                                {

+                                        tmprow = s[i];

+                                        s[i] = s[pivot_to];

+                                        s[pivot_to] = tmprow;

+                                        real tmpscp;

+                                        tmpscp = scp[i];

+                                        scp[i] = scp[pivot_to];

+                                        scp[pivot_to] = tmpscp;

+                                }

+                                

+                                real mji;

+                                // perform gaussian elimination

+                                for(j=i+1;j<4;j++)

+                                {

+                                        mji = s[j][i]/s[i][i];

+                                        s[j][i] = 0.0;

+                                        for(jj=i+1;jj<8;jj++)

+                                                s[j][jj] -= mji*s[i][jj];

+                                }

+                }

+                if(s[3][3] == 0.0) return minv; // singular matrix!

+                

+                //

+                // Now we have an upper triangular matrix.

+                //

+                //  x x x x | y y y y

+                //  0 x x x | y y y y 

+                //  0 0 x x | y y y y

+                //  0 0 0 x | y y y y

+                //

+                //  we'll back substitute to get the inverse

+                //

+                //  1 0 0 0 | z z z z

+                //  0 1 0 0 | z z z z

+                //  0 0 1 0 | z z z z

+                //  0 0 0 1 | z z z z 

+                //

+                

+                real mij;

+                for(i=3;i>0;i--)

+                {

+                        for(j=i-1;j > -1; j--)

+                        {

+                                mij = s[j][i]/s[i][i];

+                                for(jj=j+1;jj<8;jj++)

+                                        s[j][jj] -= mij*s[i][jj];

+                        }

+                }

+                

+                for(i=0;i<4;i++)

+                        for(j=0;j<4;j++)

+                                minv(i,j) = s[i][j+4] / s[i][i];

+                        

+                        return minv;

+        }

+    

+    

+    matrix4 transpose() const

+        {

+                matrix4 mtrans;

+                

+                for(int i=0;i<4;i++)

+                        for(int j=0;j<4;j++)

+                                mtrans(i,j) = element(j,i);             

+                return mtrans;

+        }

+    

+    matrix4 & mult_right( const matrix4 & b )

+        {

+                matrix4 mt(*this);

+                set_value(real(0));

+

+                for(int i=0; i < 4; i++)

+                        for(int j=0; j < 4; j++)

+                                for(int c=0; c < 4; c++)

+                                        element(i,j) += mt(i,c) * b(c,j);

+                return *this;

+        }    

+

+    matrix4 & mult_left( const matrix4 & b )

+        {

+                matrix4 mt(*this);

+                set_value(real(0));

+

+                for(int i=0; i < 4; i++)

+                        for(int j=0; j < 4; j++)

+                                for(int c=0; c < 4; c++)

+                                        element(i,j) += b(i,c) * mt(c,j);

+                return *this;

+        }

+        

+        // dst = M * src

+    void mult_matrix_vec( const vec3 &src, vec3 &dst ) const

+        {

+                real w = (

+                        src.v[0] * element(3,0) +

+                        src.v[1] * element(3,1) + 

+                        src.v[2] * element(3,2) +

+                        element(3,3)          );

+        

+        assert(w != GLH_ZERO);

+

+        dst.v[0]  = (

+                        src.v[0] * element(0,0) +

+                        src.v[1] * element(0,1) +

+                        src.v[2] * element(0,2) +

+                        element(0,3)          ) / w;

+                dst.v[1]  = (

+                        src.v[0] * element(1,0) +

+                        src.v[1] * element(1,1) +

+                        src.v[2] * element(1,2) +

+                        element(1,3)          ) / w;

+                dst.v[2]  = (

+                        src.v[0] * element(2,0) +

+                        src.v[1] * element(2,1) + 

+                        src.v[2] * element(2,2) +

+                        element(2,3)          ) / w;

+        }

+    

+        void mult_matrix_vec( vec3 & src_and_dst) const

+        { mult_matrix_vec(vec3(src_and_dst), src_and_dst); }

+

+

+    // dst = src * M

+    void mult_vec_matrix( const vec3 &src, vec3 &dst ) const

+        {

+                real w = (

+                        src.v[0] * element(0,3) +

+                        src.v[1] * element(1,3) +

+                        src.v[2] * element(2,3) +

+                        element(3,3)          );

+        

+        assert(w != GLH_ZERO);

+

+                dst.v[0]  = (

+                        src.v[0] * element(0,0) +

+                        src.v[1] * element(1,0) + 

+                        src.v[2] * element(2,0) + 

+                        element(3,0)          ) / w;

+                dst.v[1]  = (

+                        src.v[0] * element(0,1) +

+                        src.v[1] * element(1,1) +

+                        src.v[2] * element(2,1) +

+                        element(3,1)          ) / w;

+                dst.v[2]  = (

+                        src.v[0] * element(0,2) +

+                        src.v[1] * element(1,2) +

+                        src.v[2] * element(2,2) +

+                        element(3,2)          ) / w;

+        }

+        

+

+        void mult_vec_matrix( vec3 & src_and_dst) const

+        { mult_vec_matrix(vec3(src_and_dst), src_and_dst); }

+

+        // dst = M * src

+    void mult_matrix_vec( const vec4 &src, vec4 &dst ) const

+        {

+        dst.v[0]  = (

+                        src.v[0] * element(0,0) +

+                        src.v[1] * element(0,1) +

+                        src.v[2] * element(0,2) +

+                        src.v[3] * element(0,3));

+                dst.v[1]  = (

+                        src.v[0] * element(1,0) +

+                        src.v[1] * element(1,1) +

+                        src.v[2] * element(1,2) +

+                        src.v[3] * element(1,3));

+                dst.v[2]  = (

+                        src.v[0] * element(2,0) +

+                        src.v[1] * element(2,1) + 

+                        src.v[2] * element(2,2) +

+                        src.v[3] * element(2,3));

+                dst.v[3] = (

+                        src.v[0] * element(3,0) +

+                        src.v[1] * element(3,1) + 

+                        src.v[2] * element(3,2) +

+                        src.v[3] * element(3,3));

+        }

+    

+        void mult_matrix_vec( vec4 & src_and_dst) const

+        { mult_matrix_vec(vec4(src_and_dst), src_and_dst); }

+

+

+    // dst = src * M

+    void mult_vec_matrix( const vec4 &src, vec4 &dst ) const

+        {

+                dst.v[0]  = (

+                        src.v[0] * element(0,0) +

+                        src.v[1] * element(1,0) + 

+                        src.v[2] * element(2,0) + 

+                        src.v[3] * element(3,0));

+                dst.v[1]  = (

+                        src.v[0] * element(0,1) +

+                        src.v[1] * element(1,1) +

+                        src.v[2] * element(2,1) +

+                        src.v[3] * element(3,1));

+                dst.v[2]  = (

+                        src.v[0] * element(0,2) +

+                        src.v[1] * element(1,2) +

+                        src.v[2] * element(2,2) +

+                        src.v[3] * element(3,2));

+                dst.v[3] = (

+                        src.v[0] * element(0,3) +

+                        src.v[1] * element(1,3) +

+                        src.v[2] * element(2,3) +

+                        src.v[3] * element(3,3));

+        }

+        

+

+        void mult_vec_matrix( vec4 & src_and_dst) const

+        { mult_vec_matrix(vec4(src_and_dst), src_and_dst); }

+

+    

+    // dst = M * src

+    void mult_matrix_dir( const vec3 &src, vec3 &dst ) const

+        {

+                dst.v[0]  = (

+                        src.v[0] * element(0,0) +

+                        src.v[1] * element(0,1) +

+                        src.v[2] * element(0,2) ) ;

+                dst.v[1]  = ( 

+                        src.v[0] * element(1,0) +

+                        src.v[1] * element(1,1) +

+                        src.v[2] * element(1,2) ) ;

+                dst.v[2]  = ( 

+                        src.v[0] * element(2,0) +

+                        src.v[1] * element(2,1) + 

+                        src.v[2] * element(2,2) ) ;

+        }

+        

+

+        void mult_matrix_dir( vec3 & src_and_dst) const

+        { mult_matrix_dir(vec3(src_and_dst), src_and_dst); }

+

+

+        // dst = src * M

+    void mult_dir_matrix( const vec3 &src, vec3 &dst ) const

+        {

+                dst.v[0]  = ( 

+                        src.v[0] * element(0,0) +

+                        src.v[1] * element(1,0) +

+                        src.v[2] * element(2,0) ) ;

+                dst.v[1]  = ( 

+                        src.v[0] * element(0,1) +

+                        src.v[1] * element(1,1) +

+                        src.v[2] * element(2,1) ) ;

+                dst.v[2]  = (

+                        src.v[0] * element(0,2) +

+                        src.v[1] * element(1,2) + 

+                        src.v[2] * element(2,2) ) ;

+        }

+    

+    

+        void mult_dir_matrix( vec3 & src_and_dst) const

+        { mult_dir_matrix(vec3(src_and_dst), src_and_dst); }

+

+

+    real & operator () (int row, int col)

+    { return element(row,col); }

+

+    const real & operator () (int row, int col) const

+    { return element(row,col); }

+

+        real & element (int row, int col)

+    { return m[row | (col<<2)]; }

+

+    const real & element (int row, int col) const

+    { return m[row | (col<<2)]; }

+

+    matrix4 & operator *= ( const matrix4 & mat )

+        {

+                mult_right( mat );

+                return *this;

+        }

+    

+    matrix4 & operator *= ( const real & r )

+        {

+                for (int i = 0; i < 4; ++i)

+        {

+            element(0,i) *= r;

+            element(1,i) *= r;

+            element(2,i) *= r;

+            element(3,i) *= r;

+        }

+                return *this;

+        }

+

+    matrix4 & operator += ( const matrix4 & mat )

+        {

+                for (int i = 0; i < 4; ++i)

+        {

+            element(0,i) += mat.element(0,i);

+            element(1,i) += mat.element(1,i);

+            element(2,i) += mat.element(2,i);

+            element(3,i) += mat.element(3,i);

+        }

+                return *this;

+        }

+

+    friend matrix4 operator * ( const matrix4 & m1,     const matrix4 & m2 );

+    friend bool operator == ( const matrix4 & m1, const matrix4 & m2 );

+    friend bool operator != ( const matrix4 & m1, const matrix4 & m2 );

+    

+  //protected:

+          real m[16];

+  };

+  

+  inline  

+  matrix4 operator * ( const matrix4 & m1, const matrix4 & m2 )

+  {

+          matrix4 product;

+          

+          product = m1;

+          product.mult_right(m2);

+          

+          return product;

+  }

+  

+  inline

+  bool operator ==( const matrix4 &m1, const matrix4 &m2 )

+  {

+          return ( 

+                  m1(0,0) == m2(0,0) &&

+                  m1(0,1) == m2(0,1) &&

+                  m1(0,2) == m2(0,2) &&

+                  m1(0,3) == m2(0,3) &&

+                  m1(1,0) == m2(1,0) &&

+                  m1(1,1) == m2(1,1) &&

+                  m1(1,2) == m2(1,2) &&

+                  m1(1,3) == m2(1,3) &&

+                  m1(2,0) == m2(2,0) &&

+                  m1(2,1) == m2(2,1) &&

+                  m1(2,2) == m2(2,2) &&

+                  m1(2,3) == m2(2,3) &&

+                  m1(3,0) == m2(3,0) &&

+                  m1(3,1) == m2(3,1) &&

+                  m1(3,2) == m2(3,2) &&

+                  m1(3,3) == m2(3,3) );

+  }

+  

+  inline

+  bool operator != ( const matrix4 & m1, const matrix4 & m2 )

+  { return !( m1 == m2 ); }  

+

+

+

+

+

+

+

+

+

+

+

+

+  

+    class quaternion

+    {

+    public:

+    

+    quaternion()

+    {

+        *this = identity();

+    }

+

+    quaternion( const real v[4] )

+    {

+        set_value( v );

+    }

+

+

+    quaternion( real q0, real q1, real q2, real q3 )

+    {

+        set_value( q0, q1, q2, q3 );

+    }

+

+

+    quaternion( const matrix4 & m )

+    {

+        set_value( m );

+    }

+

+

+    quaternion( const vec3 &axis, real radians )

+    {

+        set_value( axis, radians );

+    }

+

+

+    quaternion( const vec3 &rotateFrom, const vec3 &rotateTo )

+    {

+        set_value( rotateFrom, rotateTo );

+    }

+

+    quaternion( const vec3 & from_look, const vec3 & from_up,

+                      const vec3 & to_look, const vec3& to_up)

+    {

+            set_value(from_look, from_up, to_look, to_up);

+    }

+

+    const real * get_value() const

+    {

+        return  &q[0];

+    }

+

+    void get_value( real &q0, real &q1, real &q2, real &q3 ) const

+    {

+        q0 = q[0];

+        q1 = q[1];

+        q2 = q[2];

+        q3 = q[3];

+    }

+

+    quaternion & set_value( real q0, real q1, real q2, real q3 )

+    {

+        q[0] = q0;

+        q[1] = q1;

+        q[2] = q2;

+        q[3] = q3;

+        counter = 0;

+        return *this;

+    }

+

+    void get_value( vec3 &axis, real &radians ) const

+    {

+        radians = real(acos( q[3] ) * GLH_TWO);

+        if ( radians == GLH_ZERO )

+            axis = vec3( 0.0, 0.0, 1.0 );

+        else

+        {

+            axis.v[0] = q[0];

+            axis.v[1] = q[1];

+            axis.v[2] = q[2];

+            axis.normalize();

+        }

+    }

+

+    void get_value( matrix4 & m ) const

+    {

+        real s, xs, ys, zs, wx, wy, wz, xx, xy, xz, yy, yz, zz;

+

+        real norm = q[0] * q[0] + q[1] * q[1] + q[2] * q[2] + q[3] * q[3];

+

+        s = (equivalent(norm,GLH_ZERO)) ? GLH_ZERO : ( GLH_TWO / norm );

+

+        xs = q[0] * s;

+        ys = q[1] * s;

+        zs = q[2] * s;

+

+        wx = q[3] * xs;

+        wy = q[3] * ys;

+        wz = q[3] * zs;

+

+        xx = q[0] * xs;

+        xy = q[0] * ys;

+        xz = q[0] * zs;

+

+        yy = q[1] * ys;

+        yz = q[1] * zs;

+        zz = q[2] * zs;

+

+        m(0,0) = real( GLH_ONE - ( yy + zz ));

+        m(1,0) = real ( xy + wz );

+        m(2,0) = real ( xz - wy );

+

+        m(0,1) = real ( xy - wz );

+        m(1,1) = real ( GLH_ONE - ( xx + zz ));

+        m(2,1) = real ( yz + wx );

+

+        m(0,2) = real ( xz + wy );

+        m(1,2) = real ( yz - wx );

+        m(2,2) = real ( GLH_ONE - ( xx + yy ));

+

+        m(3,0) = m(3,1) = m(3,2) = m(0,3) = m(1,3) = m(2,3) = GLH_ZERO;

+        m(3,3) = GLH_ONE;

+    }

+

+    quaternion & set_value( const real * qp )

+    {

+        memcpy(q,qp,sizeof(real) * 4);

+

+        counter = 0;

+        return *this;

+    }

+

+    quaternion & set_value( const matrix4 & m )

+    {

+        real tr, s;

+        int i, j, k;

+        const int nxt[3] = { 1, 2, 0 };

+

+        tr = m(0,0) + m(1,1) + m(2,2);

+

+        if ( tr > GLH_ZERO )

+        {

+            s = real(sqrt( tr + m(3,3) ));

+            q[3] = real ( s * 0.5 );

+            s = real(0.5) / s;

+

+            q[0] = real ( ( m(1,2) - m(2,1) ) * s );

+            q[1] = real ( ( m(2,0) - m(0,2) ) * s );

+            q[2] = real ( ( m(0,1) - m(1,0) ) * s );

+        }

+        else

+        {

+            i = 0;

+            if ( m(1,1) > m(0,0) )

+              i = 1;

+

+            if ( m(2,2) > m(i,i) )

+              i = 2;

+

+            j = nxt[i];

+            k = nxt[j];

+

+            s = real(sqrt( ( m(i,j) - ( m(j,j) + m(k,k) )) + GLH_ONE ));

+

+            q[i] = real ( s * 0.5 );

+            s = real(0.5 / s);

+

+            q[3] = real ( ( m(j,k) - m(k,j) ) * s );

+            q[j] = real ( ( m(i,j) + m(j,i) ) * s );

+            q[k] = real ( ( m(i,k) + m(k,i) ) * s );

+        }

+

+        counter = 0;

+        return *this;

+    }

+

+    quaternion & set_value( const vec3 &axis, real theta )

+    {

+        real sqnorm = axis.square_norm();

+

+        if (sqnorm <= GLH_EPSILON)

+        {

+            // axis too small.

+            x = y = z = 0.0;

+            w = 1.0;

+        } 

+        else 

+        {

+            theta *= real(0.5);

+            real sin_theta = real(sin(theta));

+

+            if (!equivalent(sqnorm,GLH_ONE)) 

+              sin_theta /= real(sqrt(sqnorm));

+            x = sin_theta * axis.v[0];

+            y = sin_theta * axis.v[1];

+            z = sin_theta * axis.v[2];

+            w = real(cos(theta));

+        }

+        return *this;

+    }

+

+    quaternion & set_value( const vec3 & rotateFrom, const vec3 & rotateTo )

+    {

+        vec3 p1, p2;

+        real alpha;

+

+        p1 = rotateFrom; 

+        p1.normalize();

+        p2 = rotateTo;  

+        p2.normalize();

+

+        alpha = p1.dot(p2);

+

+        if(equivalent(alpha,GLH_ONE))

+        { 

+            *this = identity(); 

+            return *this; 

+        }

+

+        // ensures that the anti-parallel case leads to a positive dot

+        if(equivalent(alpha,-GLH_ONE))

+        {

+            vec3 v;

+

+            if(p1.v[0] != p1.v[1] || p1.v[0] != p1.v[2])

+                v = vec3(p1.v[1], p1.v[2], p1.v[0]);

+            else

+                v = vec3(-p1.v[0], p1.v[1], p1.v[2]);

+

+            v -= p1 * p1.dot(v);

+            v.normalize();

+

+            set_value(v, GLH_PI);

+            return *this;

+        }

+

+        p1 = p1.cross(p2);  

+        p1.normalize();

+        set_value(p1,real(acos(alpha)));

+

+        counter = 0;

+        return *this;

+    }

+

+    quaternion & set_value( const vec3 & from_look, const vec3 & from_up,

+                      const vec3 & to_look, const vec3 & to_up)

+    {

+            quaternion r_look = quaternion(from_look, to_look);

+            

+            vec3 rotated_from_up(from_up);

+            r_look.mult_vec(rotated_from_up);

+            

+            quaternion r_twist = quaternion(rotated_from_up, to_up);

+            

+            *this = r_twist;

+            *this *= r_look;

+            return *this;

+    }

+

+    quaternion & operator *= ( const quaternion & qr )

+    {

+        quaternion ql(*this);

+   

+        w = ql.w * qr.w - ql.x * qr.x - ql.y * qr.y - ql.z * qr.z;

+        x = ql.w * qr.x + ql.x * qr.w + ql.y * qr.z - ql.z * qr.y;

+        y = ql.w * qr.y + ql.y * qr.w + ql.z * qr.x - ql.x * qr.z;

+        z = ql.w * qr.z + ql.z * qr.w + ql.x * qr.y - ql.y * qr.x;

+

+        counter += qr.counter;

+        counter++;

+        counter_normalize();

+        return *this;

+    }

+

+    void normalize()

+    {

+        real rnorm = GLH_ONE / real(sqrt(w * w + x * x + y * y + z * z));

+        if (equivalent(rnorm, GLH_ZERO))

+            return;

+        x *= rnorm;

+        y *= rnorm;

+        z *= rnorm;

+        w *= rnorm;

+        counter = 0;

+    }

+

+    friend bool operator == ( const quaternion & q1, const quaternion & q2 );      

+

+    friend bool operator != ( const quaternion & q1, const quaternion & q2 );

+

+    friend quaternion operator * ( const quaternion & q1, const quaternion & q2 );

+

+    bool equals( const quaternion & r, real tolerance ) const

+    {

+        real t;

+

+        t = (

+                        (q[0]-r.q[0])*(q[0]-r.q[0]) +

+            (q[1]-r.q[1])*(q[1]-r.q[1]) +

+            (q[2]-r.q[2])*(q[2]-r.q[2]) +

+            (q[3]-r.q[3])*(q[3]-r.q[3]) );

+        if(t > GLH_EPSILON) 

+            return false;

+        return 1;

+    }

+

+    quaternion & conjugate()

+    {

+        q[0] *= -GLH_ONE;

+        q[1] *= -GLH_ONE;

+        q[2] *= -GLH_ONE;

+        return *this;

+    }

+

+    quaternion & invert()

+    {

+        return conjugate();

+    }

+

+    quaternion inverse() const

+    {

+        quaternion r = *this;

+        return r.invert();

+    }

+

+    //

+    // Quaternion multiplication with cartesian vector

+    // v' = q*v*q(star)

+    //

+    void mult_vec( const vec3 &src, vec3 &dst ) const

+    {

+        real v_coef = w * w - x * x - y * y - z * z;                     

+        real u_coef = GLH_TWO * (src.v[0] * x + src.v[1] * y + src.v[2] * z);  

+        real c_coef = GLH_TWO * w;                                       

+

+        dst.v[0] = v_coef * src.v[0] + u_coef * x + c_coef * (y * src.v[2] - z * src.v[1]);

+        dst.v[1] = v_coef * src.v[1] + u_coef * y + c_coef * (z * src.v[0] - x * src.v[2]);

+        dst.v[2] = v_coef * src.v[2] + u_coef * z + c_coef * (x * src.v[1] - y * src.v[0]);

+    }

+

+    void mult_vec( vec3 & src_and_dst) const

+    {

+        mult_vec(vec3(src_and_dst), src_and_dst);

+    }

+

+    void scale_angle( real scaleFactor )

+    {

+        vec3 axis;

+        real radians;

+

+        get_value(axis, radians);

+        radians *= scaleFactor;

+        set_value(axis, radians);

+    }

+

+    static quaternion slerp( const quaternion & p, const quaternion & q, real alpha )

+    {

+        quaternion r;

+

+        real cos_omega = p.x * q.x + p.y * q.y + p.z * q.z + p.w * q.w;

+        // if B is on opposite hemisphere from A, use -B instead

+      

+        int bflip;

+        if ( ( bflip = (cos_omega < GLH_ZERO)) )

+            cos_omega = -cos_omega;

+

+        // complementary interpolation parameter

+        real beta = GLH_ONE - alpha;     

+

+        if(cos_omega <= GLH_ONE - GLH_EPSILON)

+            return p;

+

+        real omega = real(acos(cos_omega));

+        real one_over_sin_omega = GLH_ONE / real(sin(omega));

+

+        beta    = real(sin(omega*beta)  * one_over_sin_omega);

+        alpha   = real(sin(omega*alpha) * one_over_sin_omega);

+

+        if (bflip)

+            alpha = -alpha;

+

+        r.x = beta * p.q[0]+ alpha * q.q[0];

+        r.y = beta * p.q[1]+ alpha * q.q[1];

+        r.z = beta * p.q[2]+ alpha * q.q[2];

+        r.w = beta * p.q[3]+ alpha * q.q[3];

+        return r;

+    }

+

+    static quaternion identity()

+    {

+        static quaternion ident( vec3( 0.0, 0.0, 0.0 ), GLH_ONE );

+        return ident;

+    }

+

+    real & operator []( int i )

+    {

+        assert(i < 4);

+        return q[i];

+    }

+

+    const real & operator []( int i ) const

+    {

+        assert(i < 4);

+        return q[i];

+    }

+

+    protected:

+

+        void counter_normalize()

+        {

+            if (counter > GLH_QUATERNION_NORMALIZATION_THRESHOLD)

+                normalize();

+        }

+

+        union 

+        {

+            struct 

+            {

+                real q[4];

+            };

+            struct 

+            {

+                real x;

+                real y;

+                real z;

+                real w;

+            };

+        };

+

+        // renormalization counter

+        unsigned char counter;

+    };

+

+    inline

+    bool operator == ( const quaternion & q1, const quaternion & q2 )

+    {

+        return (equivalent(q1.x, q2.x) &&

+                        equivalent(q1.y, q2.y) &&

+                        equivalent(q1.z, q2.z) &&

+                        equivalent(q1.w, q2.w) );

+    }

+

+    inline

+    bool operator != ( const quaternion & q1, const quaternion & q2 )

+    { 

+        return ! ( q1 == q2 ); 

+    }

+

+    inline

+    quaternion operator * ( const quaternion & q1, const quaternion & q2 )

+    {   

+        quaternion r(q1); 

+        r *= q2; 

+        return r; 

+    }

+  

+      

+    

+

+

+

+

+

+  

+  

+  class plane

+  {

+  public:

+          

+          plane()

+      {

+                  planedistance = 0.0;

+                  planenormal.set_value( 0.0, 0.0, 1.0 );

+      }

+          

+          

+          plane( const vec3 &p0, const vec3 &p1, const vec3 &p2 )

+      {

+                  vec3 v0 = p1 - p0;

+                  vec3 v1 = p2 - p0;

+                  planenormal = v0.cross(v1);  

+                  planenormal.normalize();

+                  planedistance = p0.dot(planenormal);

+      }

+          

+          plane( const vec3 &normal, real distance )

+      {

+                  planedistance = distance;

+                  planenormal = normal;

+                  planenormal.normalize();

+      }

+          

+          plane( const vec3 &normal, const vec3 &point )

+      {

+                  planenormal = normal;

+                  planenormal.normalize();

+                  planedistance = point.dot(planenormal);

+      }

+          

+          void offset( real d )

+      {

+                  planedistance += d;

+      }

+          

+          bool intersect( const line &l, vec3 &intersection ) const

+      {

+                  vec3 pos, dir;

+                  vec3 pn = planenormal;

+                  real pd = planedistance;

+                  

+                  pos = l.get_position();

+                  dir = l.get_direction();

+                  

+                  if(dir.dot(pn) == 0.0) return 0;

+                  pos -= pn*pd;

+                  // now we're talking about a plane passing through the origin

+                  if(pos.dot(pn) < 0.0) pn.negate();

+                  if(dir.dot(pn) > 0.0) dir.negate();

+                  vec3 ppos = pn * pos.dot(pn);

+                  pos = (ppos.length()/dir.dot(-pn))*dir;

+                  intersection = l.get_position();

+                  intersection += pos;

+                  return 1;

+      }

+          void transform( const matrix4 &matrix )

+      {

+                  matrix4 invtr = matrix.inverse();

+                  invtr = invtr.transpose();

+                  

+                  vec3 pntOnplane = planenormal * planedistance;

+                  vec3 newPntOnplane;

+                  vec3 newnormal;

+                  

+                  invtr.mult_dir_matrix(planenormal, newnormal);

+                  matrix.mult_vec_matrix(pntOnplane, newPntOnplane);

+                  

+                  newnormal.normalize();

+                  planenormal = newnormal;

+                  planedistance = newPntOnplane.dot(planenormal);

+      }

+          

+          bool is_in_half_space( const vec3 &point ) const

+      {

+                  

+                  if(( point.dot(planenormal) - planedistance) < 0.0)

+                          return 0;

+                  return 1;

+      }

+          

+          

+          real distance( const vec3 & point ) const 

+      {

+                  return planenormal.dot(point - planenormal*planedistance);

+      }

+          

+          const vec3 &get_normal() const

+      {

+                  return planenormal;

+      }

+          

+          

+          real get_distance_from_origin() const

+      {

+                  return planedistance;

+      }

+          

+          

+          friend bool operator == ( const plane & p1, const plane & p2 );

+          

+          

+          friend bool operator != ( const plane & p1, const plane & p2 );

+          

+  //protected:

+          vec3 planenormal;

+          real planedistance;

+  };

+  

+  inline

+  bool operator == (const plane & p1, const plane & p2 )

+  {

+          return (  p1.planedistance == p2.planedistance && p1.planenormal == p2.planenormal);

+  }

+  

+  inline

+  bool operator != ( const plane & p1, const plane & p2 )

+  { return  ! (p1 == p2); }

+  

+  

+

+  } // "ns_##GLH_REAL"

+

+  // make common typedefs...

+#ifdef GLH_REAL_IS_FLOAT

+  typedef GLH_REAL_NAMESPACE::vec2 vec2f;

+  typedef GLH_REAL_NAMESPACE::vec3 vec3f;

+  typedef GLH_REAL_NAMESPACE::vec4 vec4f;

+  typedef GLH_REAL_NAMESPACE::quaternion quaternionf;

+  typedef GLH_REAL_NAMESPACE::quaternion rotationf;

+  typedef GLH_REAL_NAMESPACE::line linef;

+  typedef GLH_REAL_NAMESPACE::plane planef;

+  typedef GLH_REAL_NAMESPACE::matrix4 matrix4f;

+#endif

+

+  

+

+

+}  // namespace glh

+

+

+

+#endif

+

diff --git a/linden/indra/newview/llcloud.cpp b/linden/indra/newview/llcloud.cpp
index 3e9b86a..b325265 100644
--- a/linden/indra/newview/llcloud.cpp
+++ b/linden/indra/newview/llcloud.cpp
@@ -427,7 +427,7 @@ void LLCloudLayer::decompress(LLBitPack &bitpack, LLGroupHeader *group_headerp)
         group_headerp->stride = group_headerp->patch_size;              // offset required to step up one row
         set_group_of_patch_header(group_headerp);
 
-        decode_patch_header(bitpack, &patch_header);
+        decode_patch_header(bitpack, &patch_header, FALSE);
         decode_patch(bitpack, gBuffer);
         decompress_patch(mDensityp, gBuffer, &patch_header); 
 }
diff --git a/linden/indra/newview/llfloaterregioninfo.cpp b/linden/indra/newview/llfloaterregioninfo.cpp
index d4ffe22..4880308 100644
--- a/linden/indra/newview/llfloaterregioninfo.cpp
+++ b/linden/indra/newview/llfloaterregioninfo.cpp
@@ -1244,10 +1244,10 @@ BOOL LLPanelRegionTextureInfo::sendUpdate()
         llinfos << "LLPanelRegionTextureInfo::sendUpdate()" << llendl;

 

         // Make sure user hasn't chosen wacky textures.

-        if (!validateTextureSizes())

-        {

-                return FALSE;

-        }

+        //if (!validateTextureSizes())

+        //{

+        //      return FALSE;

+        //}

 

         LLTextureCtrl* texture_ctrl;

         std::string buffer;

diff --git a/linden/indra/newview/llglsandbox.cpp b/linden/indra/newview/llglsandbox.cpp
index 98c4d06..96fd39b 100644
--- a/linden/indra/newview/llglsandbox.cpp
+++ b/linden/indra/newview/llglsandbox.cpp
@@ -710,7 +710,7 @@ void LLViewerParcelMgr::renderOneSegment(F32 x1, F32 y1, F32 x2, F32 y2, F32 hei
 {
         // HACK: At edge of last region of world, we need to make sure the region
         // resolves correctly so we can get a height value.
-        const F32 BORDER = REGION_WIDTH_METERS - 0.1f;
+        const F32 BORDER = regionp->getWidth() - 0.1f;
 
         F32 clamped_x1 = x1;
         F32 clamped_y1 = y1;
diff --git a/linden/indra/newview/llmapresponders.cpp b/linden/indra/newview/llmapresponders.cpp
index 9d974f2..b6b5c8a 100644
--- a/linden/indra/newview/llmapresponders.cpp
+++ b/linden/indra/newview/llmapresponders.cpp
@@ -120,6 +120,8 @@ void LLMapLayerResponder::result(const LLSD& result)
                 
                         S32 x_regions = map_block["X"];
                         S32 y_regions = map_block["Y"];
+                        S32 size_x_regions = map_block["SizeX"];

+                        S32 size_y_regions = map_block["SizeY"];

                         std::string name = map_block["Name"];
                         S32 access = map_block["Access"];
                         S32 region_flags = map_block["RegionFlags"];
@@ -168,6 +170,8 @@ void LLMapLayerResponder::result(const LLSD& result)
                                 LLWorldMap::getInstance()->mSimInfoMap[handle] = siminfo;
 
                                 siminfo->mHandle = handle;
+                                siminfo->msizeX = size_x_regions;

+                                siminfo->msizeY = size_y_regions;

                                 siminfo->mName.assign( name );
                                 siminfo->mAccess = access;              /*Flawfinder: ignore*/
                                 siminfo->mRegionFlags = region_flags;
diff --git a/linden/indra/newview/llstartup.cpp b/linden/indra/newview/llstartup.cpp
index 5a3a8ee..0b6d75f 100644
--- a/linden/indra/newview/llstartup.cpp
+++ b/linden/indra/newview/llstartup.cpp
@@ -368,6 +368,8 @@ bool idle_startup()
         static U64 first_sim_handle = 0;
         static LLHost first_sim;
         static std::string first_sim_seed_cap;
+        static U32 first_sim_size_x = 256;

+        static U32 first_sim_size_y = 256;
 
         static LLVector3 initial_sun_direction(1.f, 0.f, 0.f);
         static LLVector3 agent_start_position_region(10.f, 10.f, 10.f);         // default for when no space server
@@ -1620,6 +1622,16 @@ bool idle_startup()
                                 first_sim_handle = to_region_handle(region_x, region_y);
                         }
                         
+                        text = LLUserAuth::getInstance()->getResponse("region_size_x");

+                        if(!text.empty()) {

+                                first_sim_size_x = strtoul(text.c_str(), NULL, 10);

+                                LLViewerParcelMgr::getInstance()->init(first_sim_size_x);

+                        }

+

+                        //region Y size is currently unused, major refactoring required. - Patrick Sapinski (2/10/2011)

+                        text = LLUserAuth::getInstance()->getResponse("region_size_y");

+                        if(!text.empty()) first_sim_size_y = strtoul(text.c_str(), NULL, 10);

+                        

                         const std::string look_at_str = LLUserAuth::getInstance()->getResponse("look_at");
                         if (!look_at_str.empty())
                         {
@@ -1915,7 +1927,7 @@ bool idle_startup()
 
                 gAgent.initOriginGlobal(from_region_handle(first_sim_handle));
 
-                LLWorld::getInstance()->addRegion(first_sim_handle, first_sim);
+                LLWorld::getInstance()->addRegion(first_sim_handle, first_sim, first_sim_size_x, first_sim_size_y);
 
                 LLViewerRegion *regionp = LLWorld::getInstance()->getRegionFromHandle(first_sim_handle);
                 LL_INFOS("AppInit") << "Adding initial simulator " << regionp->getOriginGlobal() << LL_ENDL;
diff --git a/linden/indra/newview/llsurface.cpp b/linden/indra/newview/llsurface.cpp
index aaafe0d..66f8076 100644
--- a/linden/indra/newview/llsurface.cpp
+++ b/linden/indra/newview/llsurface.cpp
@@ -303,7 +303,7 @@ void LLSurface::initTextures()
                 mWaterObjp = (LLVOWater *)gObjectList.createObjectViewer(LLViewerObject::LL_VO_WATER, mRegionp);
                 gPipeline.createObject(mWaterObjp);
                 LLVector3d water_pos_global = from_region_handle(mRegionp->getHandle());
-                water_pos_global += LLVector3d(128.0, 128.0, DEFAULT_WATER_HEIGHT);
+                water_pos_global += LLVector3d(mRegionp->getWidth()/2, mRegionp->getWidth()/2, DEFAULT_WATER_HEIGHT);
                 mWaterObjp->setPositionGlobal(water_pos_global);
         }
 }
@@ -356,8 +356,8 @@ void LLSurface::setOriginGlobal(const LLVector3d &origin_global)
         // Hack!
         if (mWaterObjp.notNull() && mWaterObjp->mDrawable.notNull())
         {
-                const F64 x = origin_global.mdV[VX] + 128.0;
-                const F64 y = origin_global.mdV[VY] + 128.0;
+                const F64 x = origin_global.mdV[VX] + (F64)mRegionp->getWidth()/2;
+                const F64 y = origin_global.mdV[VY] + (F64)mRegionp->getWidth()/2;
                 const F64 z = mWaterObjp->getPositionGlobal().mdV[VZ];
 
                 LLVector3d water_origin_global(x, y, z);
@@ -707,14 +707,22 @@ void LLSurface::decompressDCTPatch(LLBitPack &bitpack, LLGroupHeader *gopp, BOOL
 
         while (1)
         {
-                decode_patch_header(bitpack, &ph);
+                decode_patch_header(bitpack, &ph, b_large_patch);
                 if (ph.quant_wbits == END_OF_PATCHES)
                 {
                         break;
                 }
 
-                i = ph.patchids >> 5;
-                j = ph.patchids & 0x1F;
+                if (b_large_patch)

+                {

+                        i = ph.patchids >> 16; //x

+                        j = ph.patchids & 0xFFFF; //y

+                }

+                else

+                {

+                        i = ph.patchids >> 5; //x

+                        j = ph.patchids & 0x1F; //y

+                }
 
                 if ((i >= mPatchesPerEdge) || (j >= mPatchesPerEdge))
                 {
diff --git a/linden/indra/newview/llviewermessage.cpp b/linden/indra/newview/llviewermessage.cpp
index 9ebfd04..4de4768 100755
--- a/linden/indra/newview/llviewermessage.cpp
+++ b/linden/indra/newview/llviewermessage.cpp
@@ -3531,6 +3531,17 @@ void process_teleport_finish(LLMessageSystem* msg, void**)
         U32 teleport_flags;
         msg->getU32Fast(_PREHASH_Info, _PREHASH_TeleportFlags, teleport_flags);
         
+        U32 region_size_x = 256;

+        msg->getU32Fast(_PREHASH_Info, _PREHASH_RegionSizeX, region_size_x);

+        U32 region_size_y = 256;

+        msg->getU32Fast(_PREHASH_Info, _PREHASH_RegionSizeY, region_size_y);

+

+        //and a little hack for Second Life compatibility

+        if (region_size_y == 0 || region_size_x == 0)

+        {

+                region_size_x = 256;

+                region_size_y = 256;

+        }
         
         std::string seedCap;
         msg->getStringFast(_PREHASH_Info, _PREHASH_SeedCapability, seedCap);
@@ -3550,7 +3561,7 @@ void process_teleport_finish(LLMessageSystem* msg, void**)
 
         // Viewer trusts the simulator.
         gMessageSystem->enableCircuit(sim_host, TRUE);
-        LLViewerRegion* regionp =  LLWorld::getInstance()->addRegion(region_handle, sim_host);
+        LLViewerRegion* regionp =  LLWorld::getInstance()->addRegion(region_handle, sim_host, region_size_x, region_size_y);
 
 /*
         // send camera update to new region
@@ -3866,9 +3877,21 @@ void process_crossed_region(LLMessageSystem* msg, void**)
         std::string seedCap;
         msg->getStringFast(_PREHASH_RegionData, _PREHASH_SeedCapability, seedCap);
 
+        U32 region_size_x = 256;

+        msg->getU32(_PREHASH_RegionData, _PREHASH_RegionSizeX, region_size_x);

+        U32 region_size_y = 256;

+        msg->getU32(_PREHASH_RegionData, _PREHASH_RegionSizeY, region_size_y);

+

+        //and a little hack for Second Life compatibility

+        if (region_size_y == 0 || region_size_x == 0)

+        {

+                region_size_x = 256;

+                region_size_y = 256;

+        }
+
         send_complete_agent_movement(sim_host);
 
-        LLViewerRegion* regionp = LLWorld::getInstance()->addRegion(region_handle, sim_host);
+        LLViewerRegion* regionp = LLWorld::getInstance()->addRegion(region_handle, sim_host, region_size_x, region_size_y);
         regionp->setSeedCapability(seedCap);
 
         // Tell the LightShare handler that we have changed regions.
diff --git a/linden/indra/newview/llviewerparcelmgr.cpp b/linden/indra/newview/llviewerparcelmgr.cpp
index 630da54..abe0e9f 100644
--- a/linden/indra/newview/llviewerparcelmgr.cpp
+++ b/linden/indra/newview/llviewerparcelmgr.cpp
@@ -137,6 +137,18 @@ LLViewerParcelMgr::LLViewerParcelMgr()
         mHoverParcel = new LLParcel();
         mCollisionParcel = new LLParcel();
 
+        mBlockedImage = gImageList.getImageFromFile("noentrylines.j2c");
+        mPassImage = gImageList.getImageFromFile("noentrypasslines.j2c");
+
+        init(256);

+}

+

+//moved this stuff out of the constructor and into a function that we can call again after we get the region size.

+//LLViewerParcelMgr needs to be changed so we either get an instance per region, or it handles various region sizes

+//on a single grid properly - Patrick Sapinski (2/10/2011)

+void LLViewerParcelMgr::init(F32 region_size)

+{
+
         mParcelsPerEdge = S32(  REGION_WIDTH_METERS / PARCEL_GRID_STEP_METERS );
         mHighlightSegments = new U8[(mParcelsPerEdge+1)*(mParcelsPerEdge+1)];
         resetSegments(mHighlightSegments);
@@ -144,10 +156,9 @@ LLViewerParcelMgr::LLViewerParcelMgr()
         mCollisionSegments = new U8[(mParcelsPerEdge+1)*(mParcelsPerEdge+1)];
         resetSegments(mCollisionSegments);
 
-        mBlockedImage = gImageList.getImageFromFile("noentrylines.j2c");
-        mPassImage = gImageList.getImageFromFile("noentrypasslines.j2c");
-
-        S32 overlay_size = mParcelsPerEdge * mParcelsPerEdge / PARCEL_OVERLAY_CHUNKS;
+        S32 mParcelOverLayChunks = region_size * region_size / (128 * 128);

+

+        S32 overlay_size = mParcelsPerEdge * mParcelsPerEdge / mParcelOverLayChunks;

         sPackedOverlay = new U8[overlay_size];
 
         mAgentParcelOverlay = new U8[mParcelsPerEdge * mParcelsPerEdge];
@@ -1351,8 +1362,7 @@ void LLViewerParcelMgr::processParcelOverlay(LLMessageSystem *msg, void **user)
                 return;
         }
 
-        S32 parcels_per_edge = LLViewerParcelMgr::getInstance()->mParcelsPerEdge;
-        S32 expected_size = parcels_per_edge * parcels_per_edge / PARCEL_OVERLAY_CHUNKS;
+        S32 expected_size = 1024; //parcels_per_edge * parcels_per_edge / PARCEL_OVERLAY_CHUNKS;

         if (packed_overlay_size != expected_size)
         {
                 llwarns << "Got parcel overlay size " << packed_overlay_size
diff --git a/linden/indra/newview/llviewerparcelmgr.h b/linden/indra/newview/llviewerparcelmgr.h
index 9bf6096..dcdea3c 100644
--- a/linden/indra/newview/llviewerparcelmgr.h
+++ b/linden/indra/newview/llviewerparcelmgr.h
@@ -82,6 +82,8 @@ public:
         LLViewerParcelMgr();
         ~LLViewerParcelMgr();
 
+        void init(F32 region_size);

+

         static void cleanupGlobals();
 
         BOOL    selectionEmpty() const;
diff --git a/linden/indra/newview/llviewerparceloverlay.cpp b/linden/indra/newview/llviewerparceloverlay.cpp
index 0bcd8f3..a31f153 100644
--- a/linden/indra/newview/llviewerparceloverlay.cpp
+++ b/linden/indra/newview/llviewerparceloverlay.cpp
@@ -58,6 +58,7 @@ const U8  OVERLAY_IMG_COMPONENTS = 4;
 LLViewerParcelOverlay::LLViewerParcelOverlay(LLViewerRegion* region, F32 region_width_meters)
 :       mRegion( region ),
         mParcelGridsPerEdge( S32( region_width_meters / PARCEL_GRID_STEP_METERS ) ),
+        mRegionSize(S32(region_width_meters)),

         mDirty( FALSE ),
         mTimeSinceLastUpdate(),
         mOverlayTextureIdx(-1),
@@ -299,7 +300,8 @@ void LLViewerParcelOverlay::uncompressLandOverlay(S32 chunk, U8 *packed_overlay)
 {
         // Unpack the message data into the ownership array
         S32     size    = mParcelGridsPerEdge * mParcelGridsPerEdge;
-        S32 chunk_size = size / PARCEL_OVERLAY_CHUNKS;
+        S32 mParcelOverLayChunks = mRegionSize * mRegionSize / (128 * 128);

+        S32 chunk_size = size / mParcelOverLayChunks;

 
         memcpy(mOwnership + chunk*chunk_size, packed_overlay, chunk_size);              /*Flawfinder: ignore*/
 
diff --git a/linden/indra/newview/llviewerparceloverlay.h b/linden/indra/newview/llviewerparceloverlay.h
index 9bed1dd..d3b5980 100644
--- a/linden/indra/newview/llviewerparceloverlay.h
+++ b/linden/indra/newview/llviewerparceloverlay.h
@@ -98,6 +98,7 @@ private:
         LLViewerRegion* mRegion;
 
         S32                             mParcelGridsPerEdge;
+        S32                             mRegionSize;

 
         LLPointer<LLImageGL> mTexture;
         LLPointer<LLImageRaw> mImageRaw;
diff --git a/linden/indra/newview/llviewerregion.cpp b/linden/indra/newview/llviewerregion.cpp
index a8047b8..4186362 100644
--- a/linden/indra/newview/llviewerregion.cpp
+++ b/linden/indra/newview/llviewerregion.cpp
@@ -65,6 +65,7 @@
 #include "llvoclouds.h"
 #include "llworld.h"
 #include "llspatialpartition.h"
+#include "llviewerparcelmgr.h"
 
 // Viewer object cache version, change if object update
 // format changes. JC
@@ -199,6 +200,8 @@ LLViewerRegion::LLViewerRegion(const U64 &handle,
         if (!gNoRender)
         {
                 mParcelOverlay = new LLViewerParcelOverlay(this, region_width_meters);
+                //Re-init the parcel mgr for this sim

+                LLViewerParcelMgr::getInstance()->init(region_width_meters);
         }
         else
         {
diff --git a/linden/indra/newview/llvlmanager.cpp b/linden/indra/newview/llvlmanager.cpp
index 177093c..68b4d7b 100644
--- a/linden/indra/newview/llvlmanager.cpp
+++ b/linden/indra/newview/llvlmanager.cpp
@@ -57,21 +57,23 @@ LLVLManager::~LLVLManager()
 
 void LLVLManager::addLayerData(LLVLData *vl_datap, const S32 mesg_size)
 {
-        if (LAND_LAYER_CODE == vl_datap->mType)
-        {
-                mLandBits += mesg_size * 8;
-        }
-        else if (WIND_LAYER_CODE == vl_datap->mType)
-        {
-                mWindBits += mesg_size * 8;
-        }
-        else if (CLOUD_LAYER_CODE == vl_datap->mType)
-        {
-                mCloudBits += mesg_size * 8;
-        }
-        else
-        {
-                llerrs << "Unknown layer type!" << (S32)vl_datap->mType << llendl;
+        if (LAND_LAYER_CODE == vl_datap->mType ||

+                AURORA_LAND_LAYER_CODE == vl_datap->mType)

+        {

+                mLandBits += mesg_size * 8;

+        }

+        else if (WIND_LAYER_CODE == vl_datap->mType ||

+                AURORA_WIND_LAYER_CODE == vl_datap->mType)

+        {

+                mWindBits += mesg_size * 8;

+        }

+        else if (CLOUD_LAYER_CODE == vl_datap->mType)

+        {

+                mCloudBits += mesg_size * 8;

+        }

+        else

+        {

+                llerrs << "Unknown layer type!" << (S32)vl_datap->mType << llendl;

         }
 
         mPacketData.put(vl_datap);
@@ -90,18 +92,23 @@ void LLVLManager::unpackData(const S32 num_packets)
                 LLGroupHeader goph;
 
                 decode_patch_group_header(bit_pack, &goph);
-                if (LAND_LAYER_CODE == datap->mType)
-                {
-                        datap->mRegionp->getLand().decompressDCTPatch(bit_pack, &goph, FALSE);
-                }
-                else if (WIND_LAYER_CODE == datap->mType)
-                {
-                        datap->mRegionp->mWind.decompress(bit_pack, &goph);
-
-                }
-                else if (CLOUD_LAYER_CODE == datap->mType)
-                {
-                        datap->mRegionp->mCloudLayer.decompress(bit_pack, &goph);
+                if (LAND_LAYER_CODE == datap->mType)

+                {

+                        datap->mRegionp->getLand().decompressDCTPatch(bit_pack, &goph, FALSE);

+                }

+                else if (AURORA_LAND_LAYER_CODE == datap->mType)

+                {

+                        datap->mRegionp->getLand().decompressDCTPatch(bit_pack, &goph, TRUE);

+                }

+                else if (WIND_LAYER_CODE == datap->mType ||

+                        AURORA_WIND_LAYER_CODE == datap->mType)

+                {

+                        datap->mRegionp->mWind.decompress(bit_pack, &goph);

+                }

+                else if (CLOUD_LAYER_CODE == datap->mType ||

+                        AURORA_CLOUD_LAYER_CODE == datap->mType)

+                {

+                        datap->mRegionp->mCloudLayer.decompress(bit_pack, &goph);

                 }
         }
 
diff --git a/linden/indra/newview/llvowater.cpp b/linden/indra/newview/llvowater.cpp
index c66295a..8af9e4a 100644
--- a/linden/indra/newview/llvowater.cpp
+++ b/linden/indra/newview/llvowater.cpp
@@ -74,7 +74,7 @@ LLVOWater::LLVOWater(const LLUUID &id, const LLPCode pcode, LLViewerRegion *regi
 {
         // Terrain must draw during selection passes so it can block objects behind it.
         mbCanSelect = FALSE;
-        setScale(LLVector3(256.f, 256.f, 0.f)); // Hack for setting scale for bounding boxes/visibility.
+        setScale(LLVector3(mRegionp->getWidth(), mRegionp->getWidth(), 0.f)); // Hack for setting scale for bounding boxes/visibility.

 
         mUseTexture = TRUE;
         mIsEdgePatch = FALSE;
diff --git a/linden/indra/newview/llwind.cpp b/linden/indra/newview/llwind.cpp
index ae98bea..3dcbdf6 100644
--- a/linden/indra/newview/llwind.cpp
+++ b/linden/indra/newview/llwind.cpp
@@ -120,12 +120,12 @@ void LLWind::decompress(LLBitPack &bitpack, LLGroupHeader *group_headerp)
         set_group_of_patch_header(group_headerp);
 
         // X component
-        decode_patch_header(bitpack, &patch_header);
+        decode_patch_header(bitpack, &patch_header, FALSE);
         decode_patch(bitpack, buffer);
         decompress_patch(mVelX, buffer, &patch_header);
 
         // Y component
-        decode_patch_header(bitpack, &patch_header);
+        decode_patch_header(bitpack, &patch_header, FALSE);
         decode_patch(bitpack, buffer);
         decompress_patch(mVelY, buffer, &patch_header);
 
diff --git a/linden/indra/newview/llworld.cpp b/linden/indra/newview/llworld.cpp
index 7866bf8..525195d 100644
--- a/linden/indra/newview/llworld.cpp
+++ b/linden/indra/newview/llworld.cpp
@@ -80,12 +80,12 @@ const S32 WORLD_PATCH_SIZE = 16;
 
 extern LLColor4U MAX_WATER_COLOR;
 
-const U32 LLWorld::mWidth = 256;
+U32 LLWorld::mWidth = 256;
 
 // meters/point, therefore mWidth * mScale = meters per edge
 const F32 LLWorld::mScale = 1.f;
 
-const F32 LLWorld::mWidthInMeters = mWidth * mScale;
+F32 LLWorld::mWidthInMeters = mWidth * mScale;
 
 //
 // Functions
@@ -140,7 +140,7 @@ F32	LLWorld::getRegionMaxHeight() const
         return gHippoLimits->getMaxHeight();
 }
 
-LLViewerRegion* LLWorld::addRegion(const U64 &region_handle, const LLHost &host)
+LLViewerRegion* LLWorld::addRegion(const U64 &region_handle, const LLHost &host, const U32 &region_size_x, const U32 &region_size_y)
 {
         LLMemType mt(LLMemType::MTYPE_REGIONS);
         
@@ -172,9 +172,11 @@ LLViewerRegion* LLWorld::addRegion(const U64 &region_handle, const LLHost &host)
 
         U32 iindex = 0;
         U32 jindex = 0;
+        mWidth = region_size_x;

+        mWidthInMeters = mWidth * mScale;

         from_region_handle(region_handle, &iindex, &jindex);
-        S32 x = (S32)(iindex/mWidth);
-        S32 y = (S32)(jindex/mWidth);
+        S32 x = (S32)(iindex/256);

+        S32 y = (S32)(jindex/256);

         llinfos << "Adding new region (" << x << ":" << y << ")" << llendl;
         llinfos << "Host: " << host << llendl;
 
@@ -902,7 +904,7 @@ void LLWorld::updateWaterObjects()
         }
 
         // Region width in meters.
-        S32 const rwidth = (S32)REGION_WIDTH_U32;
+        S32 const rwidth = (S32)getRegionWidthInMeters();
 
         // The distance we might see into the void
         // when standing on the edge of a region, in meters.
@@ -1256,9 +1258,20 @@ void process_enable_simulator(LLMessageSystem *msg, void **user_data)
         // which simulator should we modify?
         LLHost sim(ip_u32, port);
 
+        U32 region_size_x = 256;

+        msg->getU32Fast(_PREHASH_SimulatorInfo, _PREHASH_RegionSizeX, region_size_x);

+        U32 region_size_y = 256;

+        msg->getU32Fast(_PREHASH_SimulatorInfo, _PREHASH_RegionSizeY, region_size_y);

+

+        if (region_size_y == 0 || region_size_x == 0)

+        {

+                region_size_x = 256;

+                region_size_y = 256;

+        }
+
         // Viewer trusts the simulator.
         msg->enableCircuit(sim, TRUE);
-        LLWorld::getInstance()->addRegion(handle, sim);
+        LLWorld::getInstance()->addRegion(handle, sim, region_size_x, region_size_y);
 
         // give the simulator a message it can use to get ip and port
         llinfos << "simulator_enable() Enabling " << sim << " with code " << msg->getOurCircuitCode() << llendl;
diff --git a/linden/indra/newview/llworld.h b/linden/indra/newview/llworld.h
index 964729d..5dbef6d 100644
--- a/linden/indra/newview/llworld.h
+++ b/linden/indra/newview/llworld.h
@@ -71,7 +71,7 @@ public:
         LLWorld();
         void destroyClass();
 
-        LLViewerRegion* addRegion(const U64 &region_handle, const LLHost &host);
+        LLViewerRegion* addRegion(const U64 &region_handle, const LLHost &host, const U32 &region_size_x, const U32 &region_size_y);
                 // safe to call if already present, does the "right thing" if
                 // hosts are same, or if hosts are different, etc...
         void                    removeRegion(const LLHost &host);
@@ -171,12 +171,12 @@ private:
         region_list_t   mCulledRegionList;
 
         // Number of points on edge
-        static const U32 mWidth;
+        static U32 mWidth;
 
         // meters/point, therefore mWidth * mScale = meters per edge
         static const F32 mScale;
 
-        static const F32 mWidthInMeters;
+        static F32 mWidthInMeters;
 
         F32 mLandFarClip;                                       // Far clip distance for land.
         LLPatchVertexArray              mLandPatch;
diff --git a/linden/indra/newview/llworldmap.cpp b/linden/indra/newview/llworldmap.cpp
index 3ada36f..ead72fa 100644
--- a/linden/indra/newview/llworldmap.cpp
+++ b/linden/indra/newview/llworldmap.cpp
@@ -237,16 +237,27 @@ LLSimInfo* LLWorldMap::simInfoFromPosGlobal(const LLVector3d& pos_global)
         return simInfoFromHandle(handle);
 }
 
-LLSimInfo* LLWorldMap::simInfoFromHandle(const U64 handle)
+LLSimInfo* LLWorldMap::simInfoFromHandle(const U64 findhandle)
 {
-        sim_info_map_t::iterator it = mSimInfoMap.find(handle);
-        if (it != mSimInfoMap.end())
-        {
-                LLSimInfo* sim_info = (*it).second;
-                if (sim_info)
-                {
-                        return sim_info;
-                }
+        std::map<U64, LLSimInfo*>::const_iterator it;

+        for (it = LLWorldMap::getInstance()->mSimInfoMap.begin(); it != LLWorldMap::getInstance()->mSimInfoMap.end(); ++it)

+        {

+                const U64 handle = (*it).first;

+                LLSimInfo* info = (*it).second;

+                if(handle == findhandle)

+                {

+                        return info;

+                }

+                U32 x = 0, y = 0;

+                from_region_handle(findhandle, &x, &y);

+                U32 checkRegionX, checkRegionY;

+                from_region_handle(handle, &checkRegionX, &checkRegionY);

+

+        if(x > checkRegionX && x < (checkRegionX + info->msizeX) &&

+                        y > checkRegionY && y < (checkRegionY + info->msizeY))

+                {

+                        return info;

+                }

         }
         return NULL;
 }
diff --git a/linden/indra/newview/llworldmap.h b/linden/indra/newview/llworldmap.h
index b7089f3..6725f91 100644
--- a/linden/indra/newview/llworldmap.h
+++ b/linden/indra/newview/llworldmap.h
@@ -81,6 +81,8 @@ public:
 
 public:
         U64 mHandle;
+        S32 msizeX;

+        S32 msizeY;

         std::string mName;
 
         F64 mAgentsUpdateTime;
diff --git a/linden/indra/newview/llworldmapview.cpp b/linden/indra/newview/llworldmapview.cpp
index 443ee74..51de598 100644
--- a/linden/indra/newview/llworldmapview.cpp
+++ b/linden/indra/newview/llworldmapview.cpp
@@ -471,8 +471,8 @@ void LLWorldMapView::draw()
                 // When the view isn't panned, 0,0 = center of rectangle
                 F32 bottom =    sPanY + half_height + relative_y;
                 F32 left =              sPanX + half_width + relative_x;
-                F32 top =               bottom + sMapScale ;
-                F32 right =             left + sMapScale ;
+                F32 top =               bottom+ (sMapScale * info->msizeY / REGION_WIDTH_METERS);

+                F32 right =             left + (sMapScale * info->msizeY / REGION_WIDTH_METERS);

 
                 // Switch to world map texture (if available for this region) if either:
                 // 1. Tiles are zoomed out small enough, or
@@ -567,17 +567,21 @@ void LLWorldMapView::draw()
                 center_global.mdV[VX] += 128.0;
                 center_global.mdV[VY] += 128.0;
 
-                S32 draw_size = llround(sMapScale);
+                S32 x_draw_size = llround(sMapScale);

+                S32 y_draw_size = llround(sMapScale);

+                x_draw_size *= info->msizeX / REGION_WIDTH_METERS;

+                y_draw_size *= info->msizeY / REGION_WIDTH_METERS;

+                

                 if (simimage != NULL)
                 {
                         simimage->setBoostLevel(LLViewerImageBoostLevel::BOOST_MAP);
-                        simimage->setKnownDrawSize(llround(draw_size * LLUI::sGLScaleFactor.mV[VX]), llround(draw_size * LLUI::sGLScaleFactor.mV[VY]));
+                        simimage->setKnownDrawSize(llround(x_draw_size * LLUI::sGLScaleFactor.mV[VX]), llround(y_draw_size * LLUI::sGLScaleFactor.mV[VY]));
                 }
 
                 if (overlayimage != NULL)
                 {
                         overlayimage->setBoostLevel(LLViewerImageBoostLevel::BOOST_MAP);
-                        overlayimage->setKnownDrawSize(llround(draw_size * LLUI::sGLScaleFactor.mV[VX]), llround(draw_size * LLUI::sGLScaleFactor.mV[VY]));
+                        overlayimage->setKnownDrawSize(llround(x_draw_size * LLUI::sGLScaleFactor.mV[VX]), llround(y_draw_size * LLUI::sGLScaleFactor.mV[VY]));
                 }
                         
 //              LLTextureView::addDebugImage(simimage);