Declaring variables mid-function changes results

Here's an alternate version of the example.

' Added Method BRL_Rotation:TQuat(), the original one, for comparison.
' Added extra displayed information.
' Reduced to single numerical example.

' Every rotation has two representations as a quaternion.
' Notice that the "good" quaternion is not the BRL quaternion. It is the other one.

' Also notice that the bad angles, although wrong, are closely related to the correct values.
' Pitch has become -Pitch and Yaw has become 180 - Yaw.


Strict

Local vec:TVec3
Local mat:TMat4

'For Local n=1 To 20
'	vec=vec3(Rnd(-90,90),Rnd(-180,180),Rnd(-180,180))
	vec = vec3( -60, 30, -120 )
	mat=TMat4.FromRotation(vec.toquat())

	Print 
'	Print n
	Print "     Original angles: " + vec.ToString()
	Print "        Becomes quat: " + vec.ToQuat().ToString()

	Print "    Mat to quat good: " + mat.rotation().ToString()
	Print "     Mat to quat bad: " + mat.buggyrotation().ToString()
	Print "Restored angles good: " + mat.rotation().ToEuler().ToString()
	Print " Restored angles bad: " + mat.buggyrotation().ToEuler().ToString()
	Print " BRL restored angles: " + mat.BRL_Rotation().ToEuler().ToString()	
'Next

Global Infinity#=1 Shl 23

Function Vec3:TVec3( x#=0.0,y#=0.0,z#=0.0 )
	Local t:TVec3=New TVec3
	t.x=x
	t.y=y
	t.z=z
	Return t
EndFunction

Function Quat:TQuat(x#,y#,z#,w#)
	Local q:TQuat=New TQuat
	q.x=x
	q.y=y
	q.z=z
	q.w=w
	Return q
EndFunction

Type TVec3 {value}

	Field x#,y#,z# {attribute}
	
	Method Copy:TVec3()
		Return Vec3( x,y,z )
	End Method
	
	Method ToString$()
		Return "Vec3("+x+","+y+","+z+")"
	End Method
	
	Method ToField$()
		Return x+","+y+","+z
	End Method
	
	Method FromField:TVec3( t$ )
		Local bits$[]=t.Split( "," )
		If bits.length<>3 Throw "Format error"
		x=bits[0].ToFloat()
		y=bits[1].ToFloat()
		z=bits[2].ToFloat()
		Return Self
	End Method

	Method Pointer:Float Ptr()
		Return Varptr x
	EndMethod

	Method Length#()
		Return Sqr( x*x+y*y+z*z )
	End Method
	
	Method Dot#( v:TVec3 )
		Return x*v.x+y*v.y+z*v.z
	End Method
	
	Method Inverse:TVec3()
		Return Vec3( -x,-y,-z )
	End Method
	
	Method Reciprocal:TVec3()
		Return Vec3( 1/x,1/y,1/z )
	End Method
	
	Method Normalize:TVec3()
		Local t#=Sqr( x*x+y*y+z*z )
		If t=0.0 Return vec3(0.0,0.0,0.0)
		Return Vec3( x/t,y/t,z/t )
	End Method
	
	Method Scale:TVec3( scale# )
		Return Vec3( x*scale,y*scale,z*scale )
	End Method 
	
	Method DistanceTo#( v:TVec3 )
		Local dx#=x-v.x,dy#=y-v.y,dz#=z-v.z
		Return Sqr( dx*dx+dy*dy+dz*dz )
	End Method
	
	Method Plus:TVec3( v:TVec3 )
		Return Vec3( x+v.x,y+v.y,z+v.z )
	End Method
	
	Method Minus:TVec3( v:TVec3 )
		Return Vec3( x-v.x,y-v.y,z-v.z )
	End Method
	
	Method Times:TVec3( v:TVec3 )
		Return Vec3( x*v.x,y*v.y,z*v.z )
	End Method
	
	Method DividedBy:TVec3( v:TVec3 )
		Return Vec3( x/v.x,y/v.y,z/v.z )
	End Method
	
	Method Cross:TVec3( v:TVec3 )
		Return Vec3( y*v.z-z*v.y,z*v.x-x*v.z,x*v.y-y*v.x )
	End Method
	
	Method Blend:TVec3( v:TVec3,Alpha# )
		Local beta#=1-Alpha
		Return Vec3( x*beta+v.x*Alpha,y*beta+v.y*Alpha,z*beta+v.z*Alpha )
	End Method
	
	Method ToQuat:TQuat()
		Local c1#=Cos(-z/2)
		Local s1#=Sin(-z/2)
		Local c2#=Cos(-x/2)
		Local s2#=Sin(-x/2)
		Local c3#=Cos( y/2)
		Local s3#=Sin( y/2)
		Local c1_c2#=c1*c2
		Local s1_s2#=s1*s2
		Local t:TQuat=New TQuat
		t.x=c1*s2*c3 - s1*c2*s3;
		t.y=c1_c2*s3 + s1_s2*c3;
		t.z=s1*c2*c3 + c1*s2*s3;
		t.w=c1_c2*c3 - s1_s2*s3;
		Return t
	EndMethod
	
	Function FromQuat:TVec3(quat:TQuat)
		Return quat.toeuler()
	EndFunction
	
EndType

Type TQuat

	Field x#,y#,z#,w#=1.0 {attribute}
	
	Method ToString$()
		Return "Quat("+x+","+y+","+z+","+w+")"
	End Method
	
	Method Inverse:TQuat()
		Return Quat( -x,-y,-z,w )
	End Method

	Method ToAngleAxis( angle# Var,axis:TVec3 Var )
		Local t#=ACos(w)
		angle=t*2
		If angle>0
			t=1/Sin(t)
			axis=Vec3( x*t,y*t,z*t )
		Else
			axis=Vec3(0,0,0)
		EndIf
	End Method
	
	Method Yaw#()
		Return ATan2( 2*y*w-2*z*x,1-2*y*y-2*x*x )
	End Method
	
	Method Pitch#()
		Return -ASin( 2*z*y+2*x*w )
	End Method
	
	Method Roll#()
		Return -ATan2( 2*z*w-2*y*x,1-2*z*z-2*x*x )
	End Method
	
	Method ToEuler:TVec3()
		Local x2#=x*x,y2#=y*y,z2#=z*z
		Local euler:TVec3=New TVec3
		euler.y#=ATan2( 2*y*w-2*z*x,1-2*y2-2*x2 )
		euler.x#=-ASin( 2*z*y+2*x*w )
		euler.z#=-ATan2( 2*z*w-2*y*x,1-2*z2-2*x2 )
		Return euler
	End Method
	
	Method Times:TQuat( q:TQuat )
		Local result:TQuat=New TQuat
		result.x#=w*q.x + x*q.w + y*q.z - z*q.y
		result.y#=w*q.y + y*q.w + z*q.x - x*q.z
		result.z#=w*q.z + z*q.w + x*q.y - y*q.x
		result.w#=w*q.w - x*q.x - y*q.y - z*q.z
		Return result
	EndMethod
	
	Method Slerp:TQuat( q:TQuat,a# )
		Local b#=1-a,f
		Local d#=x*q.x+y*q.y+z*q.z+w*q.w
		If d<0
			d=-d
			f=True
		EndIf
		If d<1
			Local om#=acos_( d )
			Local si#=sin_( om )
			a=sin_( a*om )/si
			b=sin_( b*om )/si
		EndIf
		If f a=-a
		Return Quat( x*b + q.x*a , y*b + q.y*a , z*b + q.z*a , w*b + q.w*a )
	EndMethod
	
	Method Normalize:TQuat()
		Local q:TQuat=New TQuat
		Local m#=Sqr(x*x+y*y+z*z+w*w)
		q.x=x/m
		q.y=y/m
		q.z=z/m
		q.w=w/m
		Return q
	EndMethod
	
	Function FromAngleAxis:TQuat( angle#,axis:TVec3 )
		angle:/2
		Local s#=Sin(angle)
		Local t:TQuat=New TQuat
		t.x=s*axis.x
		t.y=s*axis.y
		t.z=s*axis.z
		t.w=Cos(angle)
		Return t
	EndFunction
	
	Function FromEuler:TQuat(euler:TVec3)
		Return euler.toquat()
	EndFunction	

EndType

Type TMat4

	Field ix#,iy#,iz#,iw# {attribute}
	Field jx#,jy#,jz#,jw# {attribute}
	Field kx#,ky#,kz#,kw# {attribute}
	Field tx#,ty#,tz#,tw# {attribute}

	Method Copy:TMat4()
		Local t:TMat4=New TMat4
		MemCopy t,Self,SizeOf Self
		Return t
	End Method
	
	Method ToString$()
		Local t$="Mat4{~n"
		t:+ix+","+iy+","+iz+","+iw+"~n"
		t:+jx+","+jy+","+jz+","+jw+"~n"
		t:+kx+","+ky+","+kz+","+kw+"~n"
		t:+tx+","+ty+","+tz+","+tw+"~n"
		Return t+"}~n"
	End Method
	
	Method Pointer:Float Ptr()
		Return Varptr ix
	EndMethod
	
	Method Translation:TVec3()
		Return Vec3(tx,ty,tz)
	End Method

	Method Rotation:TQuat()		
		Local mult:Float
		Local x#,y#,z#,w#
		Local fourWSquaredMinus1:Float = +ix+jy+kz
		Local fourXSquaredMinus1:Float = +ix-jy-kz
		Local fourYSquaredMinus1:Float = +jy-ix-kz
		Local fourZSquaredMinus1:Float = +kz-ix-jy
		Local biggestIndex=0
		Local fourBiggestSquaredMinus1:Float=fourWSquaredMinus1
		If fourXSquaredMinus1 > fourBiggestSquaredMinus1
			fourBiggestSquaredMinus1 = fourXSquaredMinus1 
			biggestIndex=1
		EndIf
		If fourYSquaredMinus1 > fourBiggestSquaredMinus1
			fourBiggestSquaredMinus1 = fourYSquaredMinus1 
			biggestIndex=2
		EndIf
		If fourZSquaredMinus1 > fourBiggestSquaredMinus1
			fourBiggestSquaredMinus1 = fourZSquaredMinus1 
			biggestIndex=3
		EndIf
		Local biggestValue:Float=Sqr(fourBiggestSquaredMinus1+1.0)*0.5
		mult=0.25/biggestValue
		Select BiggestIndex
			Case 0
				w#=BiggestValue
				x#=(jz-ky)*mult
				y#=-(kx-iz)*mult
				z#=-(iy-jx)*mult
			Case 1
				x#=-BiggestValue
				w#=-( jz - ky )*mult
				y#=( iy + jx )*mult
				z#=( kx + iz )*mult
			Case 2
				y#=BiggestValue
				w#=-(kx-iz)*mult
				x#=-(iy+jx)*mult
				z#=(jz+ky)*mult
			Case 3
				z#=-BiggestValue
				w#=(iy-jx)*mult
				x#=(kx+iz)*mult
				y#=-(jz+ky)*mult
		EndSelect
		Return quat(-x,y,z,w)
	EndMethod

	Method BuggyRotation:TQuat()		
		Local mult:Float
		Local fourWSquaredMinus1:Float = +ix+jy+kz
		Local fourXSquaredMinus1:Float = +ix-jy-kz
		Local fourYSquaredMinus1:Float = +jy-ix-kz
		Local fourZSquaredMinus1:Float = +kz-ix-jy
		Local biggestIndex=0
		Local fourBiggestSquaredMinus1:Float=fourWSquaredMinus1
		If fourXSquaredMinus1 > fourBiggestSquaredMinus1
			fourBiggestSquaredMinus1 = fourXSquaredMinus1 
			biggestIndex=1
		EndIf
		If fourYSquaredMinus1 > fourBiggestSquaredMinus1
			fourBiggestSquaredMinus1 = fourYSquaredMinus1 
			biggestIndex=2
		EndIf
		If fourZSquaredMinus1 > fourBiggestSquaredMinus1
			fourBiggestSquaredMinus1 = fourZSquaredMinus1 
			biggestIndex=3
		EndIf
		Local biggestValue:Float=Sqr(fourBiggestSquaredMinus1+1.0)*0.5
		Local x#,y#,z#,w#
		mult=0.25/biggestValue
		Select BiggestIndex
			Case 0
				w#=BiggestValue
				x#=(jz-ky)*mult
				y#=-(kx-iz)*mult
				z#=-(iy-jx)*mult
			Case 1
				x#=-BiggestValue
				w#=-( jz - ky )*mult
				y#=( iy + jx )*mult
				z#=( kx + iz )*mult
			Case 2
				y#=BiggestValue
				w#=-(kx-iz)*mult
				x#=-(iy+jx)*mult
				z#=(jz+ky)*mult
			Case 3
				z#=-BiggestValue
				w#=(iy-jx)*mult
				x#=(kx+iz)*mult
				y#=-(jz+ky)*mult
		EndSelect
		Return quat(-x,y,z,w)
	EndMethod
	
	Method BRL_Rotation:TQuat()
		Local sc:TVec3=Scale()
		Local iv:TVec3=Vec3(ix,iy,iz).Scale(1/sc.x)
		Local jv:TVec3=Vec3(jx,jy,jz).Scale(1/sc.y)
		Local kv:TVec3=Vec3(kx,ky,kz).Scale(1/sc.z)
		Local x#,y#,z#,w#
		Local t#=iv.x+jv.y+kv.z+1
		If t>0
			t=Sqr(t)*2
			x=(kv.y-jv.z)/t
			y=(iv.z-kv.x)/t
			z=(jv.x-iv.y)/t
			w=t/4
		Else If iv.x>jv.y And iv.x>kv.z
			t=Sqr(iv.x-jv.y-kv.z+1)*2
			x=t/4
			y=(jv.x+iv.y)/t
			z=(iv.z+kv.x)/t
			w=(kv.y-jv.z)/t
		Else If jv.y>kv.z
			t=Sqr(jv.y-kv.z-iv.x+1)*2
			x=(jv.x+iv.y)/t
			y=t/4
			z=(kv.y+jv.z)/t
			w=(iv.z-kv.x)/t
		Else
			t=Sqr(kv.z-jv.y-iv.x+1)*2
			x=(iv.z+kv.x)/t
			y=(kv.y+jv.z)/t
			z=t/4
			w=(jv.x-iv.y)/t
		EndIf
		Return Quat( x,y,z,w )
	End Method
	
	Method Scale:TVec3()
		Return Vec3( Vec3(ix,iy,iz).Length(),Vec3(jx,jy,jz).Length(),Vec3(kx,ky,kz).Length() )
	End Method
	
	Method Times:TMat4( m:TMat4 )
		Local t:TMat4=New TMat4
		t.ix= ix*m.ix + jx*m.iy + kx*m.iz + tx*m.iw
		t.iy= iy*m.ix + jy*m.iy + ky*m.iz + ty*m.iw
		t.iz= iz*m.ix + jz*m.iy + kz*m.iz + tz*m.iw
		t.iw= iw*m.ix + jw*m.iy + kw*m.iz + tw*m.iw
		t.jx= ix*m.jx + jx*m.jy + kx*m.jz + tx*m.jw
		t.jy= iy*m.jx + jy*m.jy + ky*m.jz + ty*m.jw
		t.jz= iz*m.jx + jz*m.jy + kz*m.jz + tz*m.jw
		t.jw= iw*m.jx + jw*m.jy + kw*m.jz + tw*m.jw
		t.kx= ix*m.kx + jx*m.ky + kx*m.kz + tx*m.kw
		t.ky= iy*m.kx + jy*m.ky + ky*m.kz + ty*m.kw
		t.kz= iz*m.kx + jz*m.ky + kz*m.kz + tz*m.kw
		t.kw= iw*m.kx + jw*m.ky + kw*m.kz + tw*m.kw
		t.tx= ix*m.tx + jx*m.ty + kx*m.tz + tx*m.tw
		t.ty= iy*m.tx + jy*m.ty + ky*m.tz + ty*m.tw
		t.tz= iz*m.tx + jz*m.ty + kz*m.tz + tz*m.tw
		t.tw= iw*m.tx + jw*m.ty + kw*m.tz + tw*m.tw
		Return t
	End Method
	
	Method TimesPoint:TVec3( v:TVec3 )
		Return Vec3(..
			ix*v.x + jx*v.y + kx*v.z + tx,..
			iy*v.x + jy*v.y + ky*v.z + ty,..
			iz*v.x + jz*v.y + kz*v.z + tz )
	End Method
	
	Method TimesVector:TVec3( v:TVec3 )
		Return Vec3(..
			ix*v.x + jx*v.y + kx*v.z,..
			iy*v.x + jy*v.y + ky*v.z,..
			iz*v.x + jz*v.y + kz*v.z )
	End Method

	Method TimesNormal:TVec3( v:TVec3 )
		Local m:TMat4=Inverse()
		Return Vec3(..
			m.ix*v.x + m.iy*v.y + m.iz*v.z,..
			m.jx*v.x + m.jy*v.y + m.jz*v.z,..
			m.kx*v.x + m.ky*v.y + m.kz*v.z )
	End Method

	Method Determinant#()
		Assert Abs(iw)<=.001 And Abs(jw)<=.001 And Abs(kw<=.001) And tw>=1-.001
		Return ix*(jy*kz-jz*ky) - iy*(jx*kz-jz*kx) + iz*(jx*ky-jy*kx)
	End Method
	
	Method I:TVec3()
		Return Vec3(ix,iy,iz)
	End Method
	
	Method J:TVec3()
		Return Vec3(jx,jy,jz)
	End Method
		
	Method K:TVec3()
		Return Vec3(kx,ky,kz)
	End Method
		
	Method Cofactor:TMat4()
		Local t:TMat4=New TMat4
		t.ix= (jy*kz-jz*ky) ; t.iy=-(jx*kz-jz*kx) ; 	t.iz= (jx*ky-jy*kx)
		t.jx=-(iy*kz-iz*ky) ; t.jy= (ix*kz-iz*kx) ; 	t.jz=-(ix*ky-iy*kx)
		t.kx= (iy*jz-iz*jy) ; t.ky=-(ix*jz-iz*jx) ; 	t.kz= (ix*jy-iy*jx)
		Return t
	End Method
	
	Method Inverse:TMat4()
		Local c#=1.0/Determinant()
		Local t:TMat4=New TMat4
		t.ix= c * ( jy*kz - jz*ky )
		t.iy=-c * ( iy*kz - iz*ky )
		t.iz= c * ( iy*jz - iz*jy )
		t.jx=-c * ( jx*kz - jz*kx )
		t.jy= c * ( ix*kz - iz*kx )
		t.jz=-c * ( ix*jz - iz*jx )
		t.kx= c * ( jx*ky - jy*kx )
		t.ky=-c * ( ix*ky - iy*kx )
		t.kz= c * ( ix*jy - iy*jx )
		t.tx=-( tx*t.ix + ty*t.jx + tz*t.kx )
		t.ty=-( tx*t.iy + ty*t.jy + tz*t.ky )
		t.tz=-( tx*t.iz + ty*t.jz + tz*t.kz )
		t.tw=1
		Return t
	End Method
	
	Method Transpose:TMat4()
		Local t:TMat4=New TMat4
		t.ix=ix ; t.iy=jx ; t.iz=kx ; t.iw=tx
		t.jx=iy ; t.jy=jy ; t.jz=ky ; t.jw=ty
		t.kx=iz ; t.ky=jz ; t.kz=kz ; t.kw=tz
		t.tx=iw ; t.ty=jw ; t.tz=kw ; t.tw=tw
		Return t
	End Method
	
	Function FromDirection:TMat4( dir:TVec3,up:TVec3=Null )
		If Not up up=Vec3(0,1,0)
		Local j:TVec3=up.Normalize()
		Local k:TVec3=dir.Normalize()
		Local i:TVec3=j.Cross(k).Normalize()
		j=k.Cross(i).Normalize()
		Local t:TMat4=New TMat4
		t.ix=i.x ; t.iy=i.y ; t.iz=i.z
		t.jx=j.x ; t.jy=j.y ; t.jz=j.z
		t.kx=k.x ; t.ky=k.y ; t.kz=k.z
		t.tw=1
		Return t
	EndFunction
	
	Function FromTranslation:TMat4( v:TVec3 )
		Local t:TMat4=New TMat4
		t.ix=1 ; t.jy=1 ;	t.kz=1
		t.tx=v.x ; t.ty=v.y ; t.tz=v.z
		t.tw=1
		Return t
	EndFunction
	
	Function FromRotation:TMat4( q:TQuat )
		Local t:TMat4=New TMat4
		Local xx#=q.x*q.x,yy#=q.y*q.y,zz#=q.z*q.z
		Local xy#=q.x*q.y,xz#=q.x*q.z,yz#=q.y*q.z
		Local wx#=q.w*q.x,wy#=q.w*q.y,wz#=q.w*q.z
		t.ix=1-2*(yy+zz) ; t.iy=  2*(xy-wz) ; t.iz=  2*(xz+wy)
		t.jx=  2*(xy+wz) ; t.jy=1.0-2*(xx+zz) ; t.jz=  2*(yz-wx)
		t.kx=  2*(xz-wy) ; t.ky=  2*(yz+wx) ; t.kz=1.0-2*(xx+yy)
		t.tw=1
		Return t
	EndFunction
	
	Function FromScale:TMat4( v:TVec3 )
		Local t:TMat4=New TMat4
		t.ix=v.x
		t.jy=v.y
		t.kz=v.z
		t.tw=1
		Return t
	EndFunction
	
	Function FromUniformScale:TMat4( sc# )
		Local t:TMat4=New TMat4
		t.ix=sc
		t.jy=sc
		t.kz=sc
		t.tw=1
		Return t
	EndFunction
	
	Function FromTransRotScale:TMat4( trans:TVec3,rot:TQuat,scale:TVec3 )
		Local t:TMat4
		If trans
			t=FromTranslation(trans)
			If rot t=t.Times(FromRotation(rot))
			If scale t=t.Times(FromScale(scale))
		Else If rot
			t=FromRotation(rot)
			If scale t=t.Times(FromScale(scale))
		Else If scale
			t=FromScale(scale)
		Else
			t=Identity()
		EndIf
		Return t
	EndFunction
	
	Function FromYaw:TMat4( yaw# )	'untested!
		Local t:TMat4=New TMat4
		Local s#=Sin(yaw),c#=Cos(yaw)
		t.ix=c  ; t.iz=s ; t.jy=1
		t.kx=-s ; t.kz=c ; t.tw=1
		Return t
	EndFunction
	
	Function FromPitch:TMat4( pitch# )	'untested!
		Local t:TMat4=New TMat4
		Local s#=Sin(pitch),c#=Cos(pitch)
		t.ix=1  ; t.jy=c ; t.jz=s
		t.ky=-s ; t.kz=c ; t.tw=1
		Return t
	EndFunction
	
	Function FromRoll:TMat4( roll# )	'untested!
		Local t:TMat4=New TMat4
		Local s#=Sin(roll),c#=Cos(roll)
		t.ix=c ; t.iy=s ; t.jx=-s
		t.jy=c ; t.kz=1 ; t.tw=1
		Return t
	EndFunction
	
	Function FromYawPitchRoll:TMat4( yaw#,pitch#,roll# )
		Local t:TMat4
		If yaw
			t=FromYaw(yaw)
			If pitch t=t.Times(FromPitch(pitch))
			If roll t=t.Times(FromRoll(roll))
		Else If pitch
			t=FromPitch(pitch)
			If roll t=t.Times(FromRoll(roll))
		Else If roll
			t=FromRoll(roll)
		Else
			t=Identity()
		EndIf
		Return t
	EndFunction
	
	Function FromOrtho:TMat4( l#,r#,b#,t#,n#,f# )
		Local q:TMat4=New TMat4
		q.ix=2/(r-l)
		q.tx=-(r+l)/(r-l)
		q.jy=2/(t-b)
		q.ty=-(t+b)/(t-b)
		q.kz=2/(f-n)
		q.tz=-(f+n)/(f-n)
		q.tw=1
		Return q
	EndFunction
	
	Function FromFrustum:TMat4( near_left#,near_right#,near_bottom#,near_top#,near#,far# )
		Local t:TMat4=New TMat4
		Local near2#=near*2
		Local w#=near_right-near_left
		Local h#=near_top-near_bottom
		Local d#=far-near
		t.ix=near2/w
		t.jy=near2/h
		t.kx=(near_right+near_left)/w
		t.ky=(near_top+near_bottom)/h
		t.kz=(far+near)/d
		t.kw=1
		t.tz=-(far*near2)/d
		Return t
	EndFunction
	
	Function FromInfiniteFrustum:TMat4( near_left#,near_right#,near_bottom#,near_top#,near# )
		Local t:TMat4=New TMat4
		Local near2#=near*2
		Local w#=near_right-near_left
		Local h#=near_top-near_bottom
		Local e#=.001
		t.ix=near2/w
		t.jy=near2/w
		t.kx=(near_right+near_left)/w
		t.ky=(near_top+near_bottom)/h
		t.kz=1-e
		t.kw=1
		t.tz=near*(e-1)
		Return t
	EndFunction

	Function Identity:TMat4()
		Local t:TMat4=New TMat4
		t.ix=1 ; t.jy=1 ; t.kz=1 ; t.tw=1
		Return t
	EndFunction
	
	Method Normalize:TMat4()
		Local mat:TMat4=New TMat4
		Local vec:TVec3
		vec=i().normalize()
		mat.ix=vec.x
		mat.iy=vec.y
		mat.iz=vec.z
		vec=j().normalize()
		mat.jx=vec.x
		mat.jy=vec.y
		mat.jz=vec.z
		vec=k().normalize()
		mat.kx=vec.x
		mat.ky=vec.y
		mat.kz=vec.z
		mat.tx=tx
		mat.ty=ty
		mat.tz=tz
		mat.tw=1.0
		Return mat
	EndMethod
	
EndType

And the output is

     Original angles: Vec3(-60.0000000,30.0000000,-120.000000)
        Becomes quat: Quat(0.0473671779,0.530330062,0.789149106,0.306186199)
    Mat to quat good: Quat(-0.0473671854,-0.530330062,-0.789149106,-0.306186199)
     Mat to quat bad: Quat(0.789149106,-0.306186199,0.0473671854,-0.530330062)
Restored angles good: Vec3(-59.9999924,29.9999924,-119.999992)
 Restored angles bad: Vec3(59.9999924,150.000000,-119.999992)
 BRL restored angles: Vec3(-59.9999733,29.9999962,-119.999977)