1 #include <extUtil.h>
2 
3 /* angle between two vectors in rad, en is reference */
v_angle_ref(double * v0,double * v1,double * en)4 double v_angle_ref( double *v0, double *v1, double *en )
5 {
6   double vn0[3], vn1[3], v0v1[3];
7   double fi;
8 
9   v_norm(  v0, vn0 );
10   v_norm(  v1, vn1 );
11 
12   fi=v_sprod( vn0, vn1 );
13   if(fi>1.) fi=1.;
14   if(fi<-1.) fi=-1.;
15   fi=acos(fi);
16 
17   /* fi > PI if (v0 x v1)*en < 0. (check if the crosprod is in dir of en, else its more than pi) */
18   v_prod( v0, v1, v0v1);
19   if(v_sprod( v0v1, en) < 0.) fi=2.*PI-fi;
20   /* printf("a:%lf v1:%lf %lf %lf v2:%lf %lf %lf\n", fi*180./pi, v0[0], v0[1], v0[2], v1[0], v1[1], v1[2] ); */
21   return(fi);
22 }
23