/*
 *	clipquat.h
 *
 *	クォータニオンライブラリを定義します。
 *
 *	CLiP - Common Library for P/ECE
 *	Copyright (C) 2006 Naoyuki Sawa
 *
 *	* Mon Sep 11 01:20:26 JST 2006 Naoyuki Sawa
 *	- 作成開始。
 *	* Thu Sep 14 13:34:52 JST 2006 Naoyuki Sawa
 *	- 1st リリース。
 */
#ifndef __CLIP_QUAT_H__
#define __CLIP_QUAT_H__

/*****************************************************************************
 *	quaternion
 *****************************************************************************/

/* クォータニオン構造体の定義です。 */
typedef struct _quaternion {
	fixed x;	/* 虚数成分(I) */
	fixed y;	/* 虚数成分(J) */
	fixed z;	/* 虚数成分(K) */
	fixed w;	/* 実数成分(1) */
} quaternion;

extern const quaternion QUATERNION_0;	/* ゼロクォータニオン(和の単位元) */
extern const quaternion QUATERNION_1;	/* 単位クォータニオン(積の単位元) */

/* クォータニオンの初期化関数は用意していません。
 * GCC拡張機能を利用して、次のように書けます。
 *	quaternion q;
 *	...
 *	q = (quaternion){x_value,y_value,z_value,w_value};
 *	...
 */

/* クォータニオンの和を求めます。
 * [in,out]
 *	戻り値 = a + b
 */
quaternion qadd(quaternion a, quaternion b);

/* クォータニオンの和の逆元を求めます。(負元)
 * [in,out]
 *	戻り値 = -a
 */
quaternion qneg(quaternion a);

/* クォータニオンの積を求めます。(クォータニオン積)
 * [in,out]
 *	戻り値 = a * b
 * [note]
 *	* 回転を表すクォータニオンを合成するときは、左から掛け合わせてください。
 *	- 位置を表すクォータニオンpを、回転を表すクォータニオンq1で回転し、その後、回転を表すクォータニオンq2で回転するには、
 *		quaternion p  = (quaternion){ x, y, z, 0 };
 *		quaternion q1 = qrot(Ax1, Ay1, Az1, d1);
 *		quaternion q2 = qrot(Ax2, Ay2, Az2, d2);
 *		p = vrotq(q1, p);
 *		p = vrotq(q2, p);
 *	  回転を表すクォータニオンq1とq2をあらかじめ合成してqを求め、それをもって、位置を表すクォータニオンpを回転するには、
 *		quaternion p  = (quaternion){ x, y, z, 0 };
 *		quaternion q1 = qrot(Ax1, Ay1, Az1, d1);
 *		quaternion q2 = qrot(Ax2, Ay2, Az2, d2);
 *		quaternion q  = qmul(q2, q1);  // qmul(q1, q2) ではないことに注意!!
 *		p = vrotq(q, p);
 *	  どちらの方法でも、回転後の位置を表すクォータニオンpは同じ結果となります。
 */
quaternion qmul(quaternion a, quaternion b);

/* クォータニオンの積の逆元を求めます。(逆クォータニオン)
 * [in,out]
 *	戻り値 = a^-1
 */
quaternion qinv(quaternion a);

/* クォータニオンの各成分を同じ値で乗算します。
 * [in,out]
 *	戻り値 = a * b
 */
quaternion qfmul(quaternion a, fixed b);
quaternion qimul(quaternion a, int   b);

/* クォータニオンの各成分を同じ値で除算します。
 * [in,out]
 *	戻り値 = a / b
 */
quaternion qfdiv(quaternion a, fixed b);
quaternion qidiv(quaternion a, int   b);

/* 共役クォータニオンを求めます。
 * [in,out]
 *	戻り値 = a*
 */
quaternion qconj(quaternion a);

/* クォータニオンのノルムを求めます。
 * [in,out]
 *	戻り値 = ||a||
 */
fixed qnorm(quaternion a);

/* クォータニオンのノルムの2乗を求めます。
 * [in,out]
 *	戻り値 = ||a||^2
 */
fixed qnorm2(quaternion a);

/* クォータニオンを正規化します。
 * [in,out]
 *	戻り値 = a / ||a||
 */
quaternion qnormal(quaternion a);

/* クォータニオンの内積を求めます。
 * [in,out]
 *	戻り値 = a・b
 * [note]
 * 	* クォータニオンの内積の定義は、以下の通りです。
 *		a・b = ||a||||b||cosθ
 *	- 「3D-CGプログラマーのためのクォータニオン入門」(工学社)には、クォータニオンの内積の定義が載っていないようですが、
 *	  「クォータニオンによる座標変換(http://www.purose.net/~y_kim/QuaternionTransform.pdf)」p.14 に記述されていました。
 *	* クォータニオンa,bが回転を表すクォータニオンである場合、それらはθ/2に基づいて作成されています。(qrot()の実装参照)
 *	  従って、クォータニオンの内積が返すcosθも、実際にはcosθ/2となります。
 *	  回転を表すクォータニオンa,bの成す角度を求めるには、以下のように計算します。
 *		θ = acos(a・b / ||a||||b||) * 2
 *	  a,bは回転を表すクォータニオンなので、||a||=||b||=1 と仮定すれば、
 *		θ = acos(a・b) * 2
 *	- ところが、球面線形補間を行うときは、θ = acos(a・b) * 2 としてはいけません。
 *	  θが実質θ/2であるまま、気にせずに進めると、正しい回転の補間結果になります。
 *	  qslerp()の実装のコメントを参照してください。(2006/09/14現在、まだ、理解できていません)
 */
fixed qdot(quaternion a, quaternion b);

/* 回転を表すクォータニオンを作成します。
 * [in]
 *	x,y,z		回転軸ベクトル。
 *	d		回転角度。
 * [out]
 *	戻り値		回転を表すクォータニオン。
 * [note]
 *	* 回転軸ベクトルは、正規化されていなくても構いません。
 */
quaternion qrot(fixed x, fixed y, fixed z, fixed d);

/* クォータニオンの球面線形補間を行います。
 * [in]
 *	a		t=0のときの回転を表すクォータニオン。
 *	b		t=1のときの回転を表すクォータニオン。
 * [out]
 *	戻り値		球面線形補間されたクォータニオン。
 * [note]
 *	* 回転を表すクォータニオンは、正規化されている必要があります。
 *	- クォータニオンに対して演算を繰り返すと、誤差が蓄積します。
 *	  適切なタイミングでqnormal()を呼び出し、正規化してください。
 *	  (演算結果に対して、毎回、正規化を行う必要はありません。)
 */
quaternion qslerp(quaternion a, quaternion b, fixed t);

/* 回転を表すクォータニオンで、ベクトルを回転します。
 * [in]
 *	v		回転されるベクトル。
 *	q		回転を表すクォータニオン。
 * [out]
 *	戻り値		回転されたベクトル。
 * [note]
 *	* 回転を表すクォータニオンは、正規化されている必要があります。
 *	- クォータニオンに対して演算を繰り返すと、誤差が蓄積します。
 *	  適切なタイミングでqnormal()を呼び出し、正規化してください。
 *	  (演算結果に対して、毎回、正規化を行う必要はありません。)
 *	* 回転を表すクォータニオンは、位置を表すクォータニオンの左から掛けられるのですが、
 *	  自己代入演算マクロの利便性を考えて、ベクトル、クォータニオンの引数順としました。
 *	  利便性を優先して妥協するならば、vxform()の引数順序も同様にすべきだったのですが、
 *	  vxform()の引数順序はもう変更が難しいです。(長期間安定して使われてきた、アセンブラ化されている等の理由)
 *	  従って、vxform()の仕様は変更しないことにします。(2006/09/14 現在)
 */
vector vrotq(vector v, quaternion q);

/* 回転を表すクォータニオンで、行列を変換します。
 * [in]
 *	m		変換される行列。
 *	q		回転を表すクォータニオン。
 * [out]
 *	戻り値		変換された行列。
 * [note]
 *	* 回転を表すクォータニオンは、正規化されている必要があります。
 *	- クォータニオンに対して演算を繰り返すと、誤差が蓄積します。
 *	  適切なタイミングでqnormal()を呼び出し、正規化してください。
 *	  (演算結果に対して、毎回、正規化を行う必要はありません。)
 *	* 回転を表すクォータニオンは、回転される既存の行列の左から掛けられるのですが、
 *	  自己代入演算マクロの利便性を考えて、行列、クォータニオンの引数順としました。
 */
matrix mrotq(matrix m, quaternion q);

/*****************************************************************************
 *	自己代入演算
 *****************************************************************************/

/* quaternion */
#define    qadd_(a, b)	do { (a) =    qadd((a), (b)); } while(0)
#define    qneg_(a   )	do { (a) =    qneg((a)     ); } while(0)
#define    qmul_(a, b)	do { (a) =    qmul((a), (b)); } while(0)
#define    qinv_(a   )	do { (a) =    qinv((a)     ); } while(0)
#define   qfmul_(a, b)	do { (a) =   qfmul((a), (b)); } while(0)
#define   qimul_(a, b)	do { (a) =   qimul((a), (b)); } while(0)
#define   qfdiv_(a, b)	do { (a) =   qfdiv((a), (b)); } while(0)
#define   qidiv_(a, b)	do { (a) =   qidiv((a), (b)); } while(0)
#define   qconj_(a   )	do { (a) =   qconj((a)     ); } while(0)
#define qnormal_(a   )	do { (a) = qnormal((a)     ); } while(0)

/* vector */
#define vrotq_(v, q)	do { (v) = vrotq((v), (q)); } while(0)

/* matrix */
#define mrotq_(m, q)	do { (m) = mrotq((m), (q)); } while(0)

#endif /*__CLIP_QUAT_H__*/
