e2ir.c on line 40

module geometry;
import std.stdio,std.algorithm,std.math;

double dist(double[] p1, double[] p2){
    if((p1.length!=p2.length)||(p1.length!=3)){
        throw new Exception("Distance meant for 3d points only.");
    }
    return sqrt((p2[0]-p1[0])^^2+(p2[1]-p1[1])^^2+(p2[2]-p1[2])^^2);
}
double[] cross(double[] v1, double[] v2){
    return [v1[1]*v2[2]-v1[2]*v2[1],v1[2]*v2[0]-v1[0]*v2[2],v1[0]*v2[1]-v1[1]*v2[0]];
}
double[] unit(double[] v){
    auto output = v.dup;
    output[] /= len(output);
    return output;
}
double len(double[] v){
    return sqrt(v[0]^^2+v[1]^^2+v[2]^^2);
}

double[] makeNormal(double[][] ps){
    //Return a vector normal to the plane containing the 3 points.
    double[] s1;
    double[] s2;
    s1[] += ps[0][]-ps[1][];
    s2[] += ps[1][]-ps[2][];
    return unit(cross(s1,s2));
}

class tri{
    double[][] ps;
    double[][] es;
    double[] center;
    double[] normal;
    double radius;
    double zRadius;
    this(double[][] points,double[] normal){
        this.ps = points;
        this.es = [this.ps[1][]-this.ps[0][],this.ps[2][]-this.ps[1][],this.ps[0][]-this.ps[2][]];
        this.normal = normal;
        this.center[] += this.ps[0][]+this.ps[1][]+this.ps[2][];
        this.center[] /= 3;
        this.zRadius = max(abs(this.center[2]-ps[0][2]),abs(this.center[2]-ps[1][2]),abs(this.center[2]-ps[2][2]));
        this.radius = max(dist(center,ps[0]),dist(center,ps[1]),dist(center,ps[2]));
    }
    this(double[][] points){
        this(points,makeNormal(points));
    }
}