// cross_product.c  A X B = C   C dot A = 0,  c dot B = 0   
// cos(angle between A B) = (A dot B) / (|A|*|B|)
//                        = just A dot B if normalized length = 1
//
//  C = |A|*|B|*sin(angle between A B) normal to plane A B
//
//          | i  j  k|  c1 = a2*b3 - a3*b2
//  C = det |a1 a2 a3|  c2 = a3*b1 - a1*b3
//          |b1 b2 b3|  c3 = a1*b2 - a2*b1
//
// sub  Q-P is V  v1=q1-p1  v2=q2-p2  v3=q3-p3

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

static void cross(double A[], double B[], double C[])
{                     // zero based subscripts
  C[0] = A[1]*B[2] - A[2]*B[1];
  C[1] = A[2]*B[0] - A[0]*B[2];
  C[2] = A[0]*B[1] - A[1]*B[0];
} // end cross

static void sub(double P[], double Q[], double V[])
{                     // zero based subscripts
  V[0] = Q[0]-P[0];
  V[1] = Q[1]-P[1];
  V[2] = Q[2]-P[2];
} // end sub

static double len(double C[])
{
  return sqrt(C[0]*C[0]+C[1]*C[1]+C[2]*C[2]);
} // end len

static void normv(double V[], double C[])
{
  double L = len(V);
  C[0] = V[0]/L;
  C[1] = V[1]/L;
  C[2] = V[2]/L;
} // end normv

static double dot(double A[], double B[])
{
  return A[0]*B[0]+A[1]*B[1]+A[2]*B[2];
} // end dot

int main()
{
  double vi[3] = {1.0, 0.0, 0.0};
  double vj[3] = {0.0, 1.0, 0.0};
  double vk[3] = {0.0, 0.0, 1.0};
  double C[3];
  double vij45[3] = {0.7071, 0.7071, 0.0};
  double vjk45[3] = {0.0, 0.7071, 0.7071};
  double P[3] = {1.0, 2.0, 3.0};
  double Q1[3] = {4.0, 3.0, 2.0};
  double Q2[3] = {6.0, 7.0, 5.0};
  double V1[3];
  double V2[3];
  double Pi = 3.141592653589793238462643;
  double angle;

  printf("cross_product.c running \n");
  printf("unit vectors \n");
  printf("vi= |%f, %f, %f| \n", vi[0], vi[1], vi[2]);
  printf("vj= |%f, %f, %f| \n", vj[0], vj[1], vj[2]);
  printf("vk= |%f, %f, %f} \n", vk[0], vk[1], vk[2]);

  cross(vi, vj, C);
  printf("C= vi X vj = |%f, %f, %f| \n", C[0], C[1], C[2]);
  printf("vi dot C = %f \n", dot(vi, C));
  printf("vj dot C = %f \n", dot(vj, C));
  cross(vi, vk, C);
  printf("C= vi X vk = |%f, %f, %f| \n", C[0], C[1], C[2]);
  printf("vi dot C = %f \n", dot(vi, C));
  printf("vk dot C = %f \n", dot(vk, C));
  cross(vj, vk, C);
  printf("C= vj X vk = |%f, %f, %f| \n", C[0], C[1], C[2]);
  printf("vj dot C = %f \n", dot(vj, C));
  printf("vk dot C = %f \n", dot(vk, C));
  cross(vk, vj, C);
  printf("C= vk X vj = |%f, %f, %f| \n", C[0], C[1], C[2]);
  printf("vk dot C = %f \n", dot(vk, C));
  printf("vj dot C = %f \n", dot(vj, C));
  printf(" \n");

  printf("vij45= |%f, %f, %f| \n", vij45[0], vij45[1], vij45[2]);
  printf("vjk45= |%f, %f, %f| \n", vjk45[0], vjk45[1], vjk45[2]);
  cross(vij45, vjk45, C);
  printf("vij45 X vjk45 = |%f, %f, %f| \n", C[0], C[1], C[2]);
  printf("vij45 dot C = %f \n", dot(vij45, C));
  printf("vjk45 dot C = %f \n", dot(vjk45, C));
  printf(" \n");


  printf("P= |%f, %f, %f| \n", P[0], P[1], P[2]);
  printf("Q1= |%f, %f, %f| \n", Q1[0], Q1[1], Q1[2]);
  printf("Q2= |%f, %f, %f} \n", Q2[0], Q2[1], Q2[2]);
  sub(P, Q1, V1);
  printf("V1=Q1-P = |%f, %f, %f| \n", V1[0], V1[1], V1[2]);
  printf("len(V1) = %f \n", len(V1));
  sub(P, Q2, V2);
  printf("V2=Q2-P = |%f, %f, %f| \n", V2[0], V2[1], V2[2]);
  printf("len(V2) = %f \n", len(V2));
  angle = acos(dot(V1,V2)/(len(V1)*len(V2)));
  angle = 180.0*angle/Pi;
  printf("angle V1 to V2 = %f degrees \n", angle);
  cross(V1, V2, C);
  printf("C= V1 X V2 = |%f, %f, %f| \n", C[0], C[1], C[2]);
  printf("V1 dot C = %f \n", dot(V1, C));
  printf("V2 dot C = %f \n", dot(V2, C));
  printf("len(V1 X V2) = %f \n", len(C));
  normv(C, V1);
  printf("normv(V1 X V2) = |%f, %f, %f| \n", V1[0], V1[1], V1[2]);
  printf("len(V1 X V2) = %f \n", len(V1));

  printf(" \n");
  printf("cross_product.c finished \n");
  return 0;
} // end cross_product.c