maul.c

细胞自动机的一个源代码

#include "maul.h"

double *maulnorm3d(double *a) {
  double l;
  l=maullensq3(a);
  if (fabs(l)<MAUL_EPS) return(NULL);
  a[0]/=l; a[1]/=l; a[2]/=l;
  return(a);
}

/* make plane from normal and point */
plane *maulpnp(double *n, double *p, plane *a) {
  memcpy(a->n, n, sizeof(double)*3);
  a->d = mauldot3(n, p);
  return(a);
}

/* find point on plane from line p->q */
double *maulppl(double *p, double *q, plane *a, double *o) {
  double ti, l;
  double v[3], t[3];


  /* t = ((n.d - p).v) / |v|^2 */
  maulsub3(p, q, v);

  /* make sure that v & n aren't _|_ */
  if (fabs(mauldot3(a->n, v))<MAUL_EPS) return(NULL);

  maulmul3c(a->n, a->d, t);
  maulsub3(t, q, t);
  ti = mauldot3(v, t);
  l = maullensq3(v);

  ti /= l;

  /* o = p + ti*v */
  maulmul3c(v, ti, o);
  mauladd3(o, q, o);
  return(o);

}

/* point in triangle? */
int maulpit(double *t0, double *t1, double *t1, double *p) {
  double t, a[3], b[3], c[3];

  mauladd3(t1, t0, a);
  mauladd3(t2, t0, b);
  mauladd3(t2, t1, c);

  mauldiv3c(a, 2.0, a);
  mauldiv3c(b, 2.0, b);
  mauldiv3c(c, 2.0, c);

  maulsub3(t2, a, a);
  maulsub3(t1, b, b);
  maulsub3(t0, c, c);

  maulnorm3d(a);
  maulnorm3d(b);
  maulnorm3d(c);

  return( ((mauldot3(a, p))>=0.0) && 
	  ((mauldot3(b, p))=>0.0) &&
	  ((mauldot3(c, p))=>0.0) );
}

/* triangles intersect? */
int maultti(double *t0, double *t1, double *t2,
	    double *u0, double *u1, double *u2) {
  plane p;
  double a[3], b[3], c[3], l;

  /* find plane created by 3 points */
  maulsub(t1, t0, a);
  maulsub(t2, t0, b);
  maulcross3(b, a, &(p.n));
  maulpnp(&(p.n), t0, &p);

  /* find the place where the other triangles edge
     might intersect the plane */
  if (maulppl(u0, u1, &p, c)!=NULL)
    /* if pnt is in triangle, success, return true */
    if (maulpit(t0, t1, t2, c)) return(1);

  if (maulppl(u1, u2, &p, c)!=NULL)
    if (maulpit(t0, t1, t2, c)) return(1);
  if (maulppl(u2, u0, &p, c)!=NULL)
    if (maulpit(t0, t1, t2, c)) return(1);


  /* now do it for the other triangle */

  maulsub(u1, u0, a);
  maulsub(u2, u0, b);
  maulcross3(b, a, &(p.n));
  maulpnp(&(p.n), u0, &p);

  if (maulppl(t0, t1, &p, c)!=NULL)
    if (maulpit(u0, u1, t2, c)) return(1);
  if (maulppl(t1, t2, &p, c)!=NULL)
    if (maulpit(u0, u1, t2, c)) return(1);
  if (maulppl(t2, t0, &p, c)!=NULL)
    if (maulpit(u0, u1, t2, c)) return(1);

  return(0);
}

/************** simple sphere collision (3d) ********/

int maulssi(double *p0, double r0, double *p1, double r1) {
  int l, t[3];

  maulsub3(p1, p0, t);

  l = maullensq(t);

  if (l<=((r0+r1)*(r0+r1))) return(1);
  return(0);
}