/* vektor3fklasse.cc

 Rechnen mit 3D Vektoren mit float (32-Bit-Fliesszahlen)
Anwendung: einfuegen in C++ Programmen
*/

#ifndef VEKTOR3FKLASSE_CC
#define VEKTOR3FKLASSE_CC
#include <math.h>
#ifndef PI
#define PI 3.14159265358979323846
#define ZWEIPI (2.0*PI)
#define PIHALBE (0.5*PI)
#define GRAD (PI/180.0)
#endif

inline float absf(float x) {return (x<0) ? -x : x;}

class Vektor3f
{
public:
 float x,y,z;
 Vektor3f() {x=y=z=0;}
 Vektor3f(float a,float b=0,float c=0) {x=a; y=b; z=c;}
 float get(int i) {if(i==0) return x; if(i==1) return y; if(i==2) return z;
		    return sqrtf(x*x+y*y+z*z); //Betrag
		   }
 Vektor3f operator=(float a) {x=a; y=z=0; return *this;}
 Vektor3f operator=(Vektor3f v) {x=v.x; y=v.y; z=v.z; return *this;}
 Vektor3f operator+=(Vektor3f v) {x+=v.x; y+=v.y; z+=v.z; return *this;}
 Vektor3f operator-=(Vektor3f v) {x-=v.x; y-=v.y; z-=v.z; return *this;}
 Vektor3f operator*=(float a) {x*=a; y*=a; z*=a; return *this;}
 Vektor3f operator/=(float a) {x/=a; y/=a; z/=a; return *this;}
 Vektor3f einheitsvektor()
	{float wu=sqrtf(x*x+y*y+z*z); x/=wu; y/=wu; z/=wu; return *this;}
 float& operator[](int i) {return i==0?x:(i==1?y:z);}
 friend Vektor3f operator+(Vektor3f,Vektor3f);
 friend Vektor3f operator-(Vektor3f,Vektor3f);
 friend Vektor3f operator*(Vektor3f,float);
 friend Vektor3f operator*(float,Vektor3f);
 friend Vektor3f operator*(Vektor3f,Vektor3f);
 friend Vektor3f operator/(Vektor3f,float);
 friend Vektor3f operator/(Vektor3f,Vektor3f);//provi.
 friend Vektor3f operator-(Vektor3f);
 void drehenxy(float sw,float cw);
 void drehenyx(float sw,float cw);
 void drehenxy(float w) {drehenxy(sinf(w),cosf(w));}
 void drehenyx(float w) {drehenxy(-sinf(w),cosf(w));}
 void drehenzx(float sw,float cw);
 void drehenxz(float sw,float cw);
 void drehenzx(float w) {drehenzx(sinf(w),cosf(w));}
 void drehenxz(float w) {drehenzx(-sinf(w),cosf(w));}
 void drehenyz(float sw,float cw);
 void drehenzy(float sw,float cw);
 void drehenyz(float w) {drehenyz(sinf(w),cosf(w));}
 void drehenzy(float w) {drehenyz(-sinf(w),cosf(w));}
};

Vektor3f operator+(Vektor3f v1,Vektor3f v2)
{
 Vektor3f v; v.x=v1.x+v2.x; v.y=v1.y+v2.y; v.z=v1.z+v2.z;
 return v;
}
Vektor3f operator-(Vektor3f v1,Vektor3f v2)
{
 Vektor3f v; v.x=v1.x-v2.x; v.y=v1.y-v2.y; v.z=v1.z-v2.z;
 return v;
}
Vektor3f operator*(Vektor3f v1,float c)
{
 Vektor3f v; v.x=v1.x*c; v.y=v1.y*c; v.z=v1.z*c;
 return v;
}
Vektor3f operator*(float c,Vektor3f v1)
{
 Vektor3f v; v.x=v1.x*c; v.y=v1.y*c; v.z=v1.z*c;
 return v;
}
Vektor3f operator/(Vektor3f v1,float c)
{
 Vektor3f v; v.x=v1.x/c; v.y=v1.y/c; v.z=v1.z/c;
 return v;
}
Vektor3f operator-(Vektor3f v1)
{
 Vektor3f v; v.x= -v1.x; v.y= -v1.y; v.z= -v1.z;
 return v;
}

float abs(Vektor3f v)				//Betrag eines Vektors
{
 return sqrtf(v.x*v.x+v.y*v.y+v.z*v.z);
}
float absquadrat(Vektor3f v)
{
 return v.x*v.x+v.y*v.y+v.z*v.z;
}

float Vektor3fprodukt(Vektor3f v1,Vektor3f v2)	// v1*v2
{
 return v1.x*v2.x+v1.y*v2.y+v1.z*v2.z;
}

Vektor3f operator*(Vektor3f v1,Vektor3f v2)
{
 Vektor3f v; v.x=v1.x*v2.x; v.y=v1.y*v2.y; v.z=v1.z*v2.z;
 return v;
}
Vektor3f operator/(Vektor3f v1,Vektor3f v2) //provi.
{
 Vektor3f v; v.x=v1.x/v2.x; v.y=v1.y/v2.y; v.z=v1.z/v2.z;
 return v;
}

float xwinkel(Vektor3f v)		//Winkel zwischen Vektor und x-Achse
{return atan2(sqrtf(v.z*v.z+v.y*v.y),v.x);}	//kommt im Bereich +-PI heraus

Vektor3f vektorkreuzprodukt(Vektor3f b,Vektor3f a)	//Kreuzprodukt
{
 Vektor3f v;
 v.x =  b.y*a.z-b.z*a.y;
 v.y = -b.x*a.z+b.z*a.x;
 v.z =  b.x*a.y-b.y*a.x;
 return v;
}

Vektor3f sin(Vektor3f a) {Vektor3f v(sinf(a.x),sinf(a.y),sinf(a.z)); return v;}
Vektor3f cos(Vektor3f a) {Vektor3f v(cosf(a.x),cosf(a.y),cosf(a.z)); return v;}
Vektor3f tanf(Vektor3f a) {Vektor3f v(tanf(a.x),tanf(a.y),tanf(a.z)); return v;}
Vektor3f asin(Vektor3f a) {Vektor3f v(asinf(a.x),asinf(a.y),asinf(a.z));return v;}
Vektor3f acos(Vektor3f a) {Vektor3f v(acosf(a.x),acosf(a.y),acosf(a.z));return v;}
Vektor3f atanf(Vektor3f a) {Vektor3f v(atanf(a.x),atanf(a.y),atanf(a.z));return v;}
Vektor3f sqrt(Vektor3f a) {Vektor3f v(sqrtf(a.x),sqrtf(a.y),sqrtf(a.z));return v;}
Vektor3f log(Vektor3f a) {Vektor3f v(log(a.x),log(a.y),log(a.z));return v;}
Vektor3f log10(Vektor3f a)
	{Vektor3f v(log10(a.x),log10(a.y),log10(a.z));return v;}
Vektor3f exp(Vektor3f a) {Vektor3f v(exp(a.x),exp(a.y),exp(a.z));return v;}
Vektor3f pow(float b,Vektor3f a)
	{Vektor3f v(pow(b,a.x),pow(b,a.y),pow(b,a.z));return v;}
Vektor3f pow(Vektor3f b,Vektor3f a)
	{Vektor3f v(pow(b.x,a.x),pow(b.y,a.y),pow(b.z,a.z));return v;}

float vektorwinkel(Vektor3f v1,Vektor3f v2)	//Winkel zwischen 2 Vektoren
{
 float a=abs(vektorkreuzprodukt(v1,v2))/(abs(v1)*abs(v2));
 return asinf(a);
}
float vektorwinkel2(Vektor3f a,Vektor3f b)   //bessere Variante
{
 float cosa=(a.x*b.x+a.y*b.y+a.z*b.z)/(abs(a)*abs(b));
 return acosf(cosa);
}
float sinvektorwinkel(Vektor3f v1,Vektor3f v2)	//sin(Winkel) zw. 2 Vektoren
{
 float sina=abs(vektorkreuzprodukt(v1,v2))/(abs(v1)*abs(v2));
 return sina;
}
float cosvektorwinkel(Vektor3f a,Vektor3f b)   //cos(Winkel)
{
 return (a.x*b.x+a.y*b.y+a.z*b.z)/(abs(a)*abs(b));
}

void Vektor3f::drehenxy(float sina,float cosa)
{					//Vektor in xy-Ebene
 float h=x*cosa-y*sina;		//im gegenuhrzeigersinn drehen
 y=x*sina+y*cosa; x=h;
}
void Vektor3f::drehenyx(float sina,float cosa)
{					//Vektor in xy-Ebene
 float h=x*cosa+y*sina;		//im uhrzeigersinn drehen
 y=y*cosa-x*sina; x=h;
}
void Vektor3f::drehenzx(float sina,float cosa)
{					//Vektor in zx-Ebene drehen
 float h=z*cosa-x*sina;
 x=z*sina+x*cosa; z=h;
}
void Vektor3f::drehenxz(float sina,float cosa)
{
 float h=z*cosa+x*sina;
 x=x*cosa-z*sina; z=h;
}
void Vektor3f::drehenyz(float sina,float cosa)
{					//Vektor in yz-Ebene drehen
 float h=y*cosa-z*sina;
 z=y*sina+z*cosa; y=h;
}
void Vektor3f::drehenzy(float sina,float cosa)
{
 float h=y*cosa+z*sina;
 z=z*cosa-y*sina; y=h;
}

/** weiss noch nicht wie machen:
float vektorwinkelv(Vektor3f v1,Vektor3f v2)	//Winkel zwischen 2 Vektoren
{						//mit Vorzeichen
}
vektor vektordrehen(Vektor3f v,float w,Vektor3f axe) //v um Winkel w drehen
{						//um die Drehachse axe
						//im Uhrzeigersinn
 Vektor3f vneu;
 return vneu;
}
**/

float sinwinkel_zwischen_vektor_und_flaeche
	(Vektor3f v,Vektor3f p2,Vektor3f p3)
{
 float sina,cosa,wu,wu2;
 // 1.Drehung so dass p2 auf xz-Ebene zu liegen kommt:
 wu=sqrtf(wu2=p2.x*p2.x+p2.y*p2.y);
 sina=p2.y/wu; cosa=p2.x/wu;
 v.drehenyx(sina,cosa);
 p2.x=wu; p2.y=0;
 p3.drehenyx(sina,cosa);
 // 2.Drehung so dass p2 auf x-Achse zu liegen kommt:
 wu=sqrtf(wu2+p2.z*p2.z);
 sina=p2.z/wu; cosa=p2.x/wu;
 v.drehenzx(sina,cosa);
 p3.drehenzx(sina,cosa);
 // 3.Drehung so dass p3 auf xy-Flaeche zu liegen kommt:
 wu=sqrtf(p3.y*p3.y+p3.z*p3.z);
 sina=p3.z/wu; cosa=p3.y/wu;
 v.drehenzy(sina,cosa);
 // Der gesuchte Winkel ist jetzt der Winkel von v auf xy-Ebene:
 return absf(v.z)/abs(v);
}

float winkel_zwischen_vektor_und_flaeche
	(Vektor3f v,Vektor3f p1,Vektor3f p2,Vektor3f p3)
{
 return asinf(sinwinkel_zwischen_vektor_und_flaeche(v-p1,p2-p1,p3-p1));
}
#endif
