/* check_mat.c  verify matricies from demo_mat.c */

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

#ifndef M_PI
#define M_PI 3.14159
#endif

static void prmat(double m[16]);
static void matmul(double a[16], double b[16], double c[16]); /* a*b into c */
static void matcpy(double a[16], double c[16]); /* a into c */
static void vecmul(double a[16], double p[4], double pp[4]); /* a*p into pp */
static double mmatrix[16]= {1.0, 0.0, 0.0, 0.0,
                            0.0, 1.0, 0.0, 0.0,
                            0.0, 0.0, 1.0, 0.0,
                            0.0, 0.0, 0.0, 1.0};
static double pmatrix[16];
static double rmatrix[16]= {1.0, 0.0, 0.0, 0.0,
                            0.0, 1.0, 0.0, 0.0,
                            0.0, 0.0, 1.0, 0.0,
                            0.0, 0.0, 0.0, 1.0};
static double tmatrix[16]= {1.0, 0.0, 0.0, 0.0,
                            0.0, 1.0, 0.0, 0.0,
                            0.0, 0.0, 1.0, 0.0,
                            0.0, 0.0, 0.0, 1.0};
static double c[16];
static double p[4] = {0.5, 1.5, 2.5, 1.0};
static double pll[4];
static double pur[4];
static double pllf[4];
static double purf[4];
static double pp[4];
static double size = 2.0;
static double xmin, xmax, ymin, ymax, near, far;
static double w=300.0;
static double h=400.0;
static double x, y, z;

int main(int argc, char * argv[])
{
  printf("output from check_mat.c \n");
  
  xmin=-size;
  xmax= size;
  ymin=-size*h/w;
  ymax= size*h/w;
  near= 1.0*size;
  far = 100.0*size;
  
  printf("xmin=%f, xmax=%f \n", xmin, xmax);
  printf("ymin=%f, ymax=%f \n", ymin, ymax);
  printf("near=%f, far =%f \n", near, far);

  pll[0]=xmin;        pll[1]=ymin;        pll[2]=-near; pll[3]=1.0;
  pur[0]=xmax;        pur[1]=ymax;        pur[2]=-near; pur[3]=1.0;
  pllf[0]=100.0*xmax; pllf[1]=100.0*ymin; pllf[2]=-far; pllf[3]=1.0;
  purf[0]=100.0*xmin; purf[1]=100.0*ymax; purf[2]=-far; purf[3]=1.0;
  
  pmatrix[4*0+0]=2.0*near/(xmax-xmin);
  pmatrix[4*0+1]=0.0;
  pmatrix[4*0+2]=(xmax+xmin)/(xmax-xmin);
  pmatrix[4*0+3]=0.0;
  pmatrix[4*1+0]=0.0;
  pmatrix[4*1+1]=2.0*near/(ymax-ymin);
  pmatrix[4*1+2]=(ymax+ymin)/(ymax-ymin);
  pmatrix[4*1+3]=0.0;
  pmatrix[4*2+0]=0.0;
  pmatrix[4*2+1]=0.0;
  pmatrix[4*2+2]=-(far+near)/(far-near);
  pmatrix[4*2+3]=-2.0*far*near/(far-near);
  pmatrix[4*3+0]=0.0;
  pmatrix[4*3+1]=0.0;
  pmatrix[4*3+2]=-1.0;
  pmatrix[4*3+3]=0.0;
  
  printf("\n perspective pmatrix \n");
  prmat(pmatrix);



  printf("\n mmatrix identity \n");
  prmat(mmatrix);
  
  vecmul(mmatrix, p, pp);
  printf("\n mmatrix * p=%f, %f, %f, %f \n", p[0], p[1], p[2], p[3]);
  printf("   equals pp=%f, %f, %f, %f \n", pp[0], pp[1], pp[2], pp[3]);

  rmatrix[4*1+1]= cos(M_PI/6.0); /* 30 degrees */
  rmatrix[4*2+2]= cos(M_PI/6.0); /* 30 degrees */
  rmatrix[4*1+2]=-sin(M_PI/6.0); /* 30 degrees */
  rmatrix[4*2+1]= sin(M_PI/6.0); /* 30 degrees */
  matmul(mmatrix, rmatrix, c);
  matcpy(c, mmatrix);
  printf("\n mmatrix 30 deg rotated about X \n");
  prmat(mmatrix);
  
  vecmul(mmatrix, p, pp);
  printf("\n mmatrix * p=%f, %f, %f, %f \n", p[0], p[1], p[2], p[3]);
  printf("   equals pp=%f, %f, %f, %f \n", pp[0], pp[1], pp[2], pp[3]);

  rmatrix[4*1+1]= 1.0;
  rmatrix[4*2+2]= 1.0;
  rmatrix[4*1+2]= 0.0;
  rmatrix[4*2+1]= 0.0;
  rmatrix[4*0+0]= cos(M_PI/6.0); /* 30 degrees */
  rmatrix[4*2+2]= cos(M_PI/6.0); /* 30 degrees */
  rmatrix[4*0+2]= sin(M_PI/6.0); /* 30 degrees */
  rmatrix[4*2+0]=-sin(M_PI/6.0); /* 30 degrees */
  matmul(mmatrix, rmatrix, c);
  matcpy(c, mmatrix);
  printf("\n mmatrix rotated 30 deg about X and Y \n");
  prmat(mmatrix);
  
  vecmul(mmatrix, p, pp);
  printf("\n mmatrix * p=%f, %f, %f, %f \n", p[0], p[1], p[2], p[3]);
  printf("   equals pp=%f, %f, %f, %f \n", pp[0], pp[1], pp[2], pp[3]);

  rmatrix[4*0+0]= 1.0;
  rmatrix[4*2+2]= 1.0;
  rmatrix[4*0+2]= 0.0;
  rmatrix[4*2+0]= 0.0;
  rmatrix[4*0+0]= cos(M_PI/6.0); /* 30 degrees */
  rmatrix[4*1+1]= cos(M_PI/6.0); /* 30 degrees */
  rmatrix[4*0+1]=-sin(M_PI/6.0); /* 30 degrees */
  rmatrix[4*1+0]= sin(M_PI/6.0); /* 30 degrees */
  matmul(mmatrix, rmatrix, c);
  matcpy(c, mmatrix);
  printf("\n mmatrix rotated 30 deg about X,Y, and Z \n");
  prmat(mmatrix);
  
  vecmul(mmatrix, p, pp);
  printf("\n mmatrix * p=%f, %f, %f, %f \n", p[0], p[1], p[2], p[3]);
  printf("   equals pp=%f, %f, %f, %f \n", pp[0], pp[1], pp[2], pp[3]);

  tmatrix[4*0+3] = 1.0;
  tmatrix[4*1+3] = 2.0;
  tmatrix[4*2+3] = -20.0;
  printf("\n tmatrix  \n");
  prmat(tmatrix);
  matmul(mmatrix, tmatrix, c);
  matcpy(c, mmatrix);
  printf("\n mmatrix rotated and translated \n");
  prmat(mmatrix);
  
  vecmul(mmatrix, p, pp);
  printf("\n mmatrix * p=%f, %f, %f, %f \n", p[0], p[1], p[2], p[3]);
  printf("   equals pp=%f, %f, %f, %f \n", pp[0], pp[1], pp[2], pp[3]);
  
  matmul(pmatrix, mmatrix, c);
  vecmul(c, p, pp);
  printf("\n pmatrix * mmatrix * p=%f, %f, %f, %f \n", p[0], p[1], p[2], p[3]);
  printf("   equals pp=%f, %f, %f, %f \n", pp[0], pp[1], pp[2], pp[3]);
  printf("   normal pp=%f, %f, %f, %f \n", pp[0]/pp[3], pp[1]/pp[3], pp[2]/pp[3], pp[3]/pp[3]);
  printf("  x/z,y/z pp=%f, %f \n", pp[0]/pp[2], pp[1]/pp[2]);

  x = w*(1.0+pp[0]/pp[3])/2.0;
  y = h*(1.0+pp[1]/pp[3])/2.0;
  z =   (1.0+pp[2]/pp[3])/2.0;
  printf("\n  x scr, y scr=%f, %f  at relative z=%f)\n", x, y, z);
  printf("width, height =%f, %f \n", w, h);
  printf("from x=0, y=0 in lower left \n");  

  printf("\n now check corners \n");
  vecmul(pmatrix, pll, pp);
  printf("\n pmatrix * pll=%f, %f, %f, %f \n", pll[0], pll[1], pll[2], pll[3]);
  printf("pp=%f, %f, %f, %f,  x/z=%f,  y/z=%f \n",
         pp[0], pp[1], pp[2], pp[3], pp[0]/pp[2], pp[1]/pp[2]);
  x = w*(1.0+pp[0]/pp[3])/2.0;
  y = h*(1.0+pp[1]/pp[3])/2.0;
  printf("  x scr, y scr=%f, %f  at z=%f \n", x, y, (1.0+pp[2]/pp[3])/2.0);

  vecmul(pmatrix, pur, pp);
  printf("\n pmatrix * pur=%f, %f, %f, %f \n", pur[0], pur[1], pur[2], pur[3]);
  printf("pp=%f, %f, %f, %f,  x/z=%f,  y/z=%f \n",
         pp[0], pp[1], pp[2], pp[3], pp[0]/pp[2], pp[1]/pp[2]);
  x = w*(1.0+pp[0]/pp[3])/2.0;
  y = h*(1.0+pp[1]/pp[3])/2.0;
  printf("  x scr, y scr=%f, %f  at z=%f \n", x, y, (1.0+pp[2]/pp[3])/2.0);

  vecmul(pmatrix, pllf, pp);
  printf("\n pmatrix * pll=%f, %f, %f, %f \n", pllf[0], pllf[1], pllf[2], pllf[3]);
  printf("pp=%f, %f, %f, %f,  x/z=%f,  y/z=%f \n",
         pp[0], pp[1], pp[2], pp[3], pp[0]/pp[2], pp[1]/pp[2]);
  x = w*(1.0+pp[0]/pp[3])/2.0;
  y = h*(1.0+pp[1]/pp[3])/2.0;
  printf("  x scr, y scr=%f, %f  at z=%f \n", x, y, (1.0+pp[2]/pp[3])/2.0);

  vecmul(pmatrix, purf, pp);
  printf("\n pmatrix * purf=%f, %f, %f, %f \n", purf[0], purf[1], purf[2], purf[3]);
  printf("pp=%f, %f, %f, %f,  x/z=%f,  y/z=%f \n",
         pp[0], pp[1], pp[2], pp[3], pp[0]/pp[2], pp[1]/pp[2]);
  x = w*(1.0+pp[0]/pp[3])/2.0;
  y = h*(1.0+pp[1]/pp[3])/2.0;
  printf("  x scr, y scr=%f, %f  at z=%f \n", x, y, (1.0+pp[2]/pp[3])/2.0);

  return 0;
}

static void prmat(double m[16])
{
  printf("m[ 0.. 3] %f, %f, %f, %f \n",
          m[0], m[1], m[2], m[3]);
  printf("m[ 4.. 7] %f, %f, %f, %f \n",
          m[4], m[5], m[6], m[7]);
  printf("m[ 8..11] %f, %f, %f, %f \n",
          m[8], m[9], m[10], m[11]);
  printf("m[12..15] %f, %f, %f, %f \n",
          m[12], m[13], m[14], m[15]);
}

static void matmul(double a[16], double b[16], double c[16])
{
  int i, j, k;
  
  for(i=0; i<4; i++)
  {
    for(j=0; j<4; j++)
    {
      c[4*i+j]=0.0;
      for(k=0; k<4; k++)
      {
        c[4*i+j] = c[4*i+j] + a[4*i+k] * b[4*k+j];
      }
    }
  }
}

static void matcpy(double a[16], double c[16])
{
  int i, j;
  
  for(i=0; i<4; i++)
  {
    for(j=0; j<4; j++)
    {
      c[4*i+j]=a[4*i+j];
    }
  }
}


static void vecmul(double a[16], double p[4], double pp[4])
{
  int i, j;
  
  for(i=0; i<4; i++)
  {
    pp[i] = 0.0;
    for(j=0; j<4; j++)
    {
      pp[i]= pp[i] + a[4*i+j] * p[j];
    }
  }
}

