Geodesic slicing using geodesic library from thing:10725

geodesic.scad

include <platonic.scad>

//////////////////////////////////////////////////////////////////////////////////////////////////
// Matrix Math Helpers
function remove(list, removeI, i=0) = (i >= len(list))
    ? []
    : concat(i==removeI?[]:[list[i]], remove(list, removeI, i+1));
    
function det(m) = let(r=[for(i=[0:1:len(m)-1]) i]) det_help(m, 0, r);
function det_help(m, i, r) = len(r) == 0 ? 1 :
 m[len(m)-len(r)][r[i]]*det_help(m,0,remove(r,i)) - (i+1<len(r)?
det_help(m, i+1, r) : 0);

function matrix_invert(m) = let(r=[for(i=[0:len(m)-1]) i]) [for(i=r)
[for(j=r)
 ((i+j)%2==0 ? 1:-1) * matrix_minor(m,0,remove(r,j),remove(r,i))]] / det(m);
function matrix_minor(m,k,ri, rj) = let(len_r=len(ri)) len_r == 0 ? 1 :
 m[ri[0]][rj[k]]*matrix_minor(m,0,remove(ri,0),remove(rj,k)) -
(k+1<len_r?matrix_minor(m,k+1,ri,rj) : 0);

//////////////////////////////////////////////////////////////////////////////////////////////////
// Polyhedron Building Code

function sp_index_of_or_add(list, value) = let(
    searchRes = search([value], list)
) [
    searchRes == [[]] ? len(list) : searchRes[0],
    searchRes == [[]] ? concat(list, [value]) : list
];

function sp_indexes_of_or_add(list, values, indexes=[], i=0) = i >= len(values)
    ? [indexes, list]
    : let(res=sp_index_of_or_add(list, values[i])) sp_indexes_of_or_add(res[1], values, concat(indexes, [res[0]]), i+1);

function polyhedron_from_faces(inputFaces, outputPoints=[], outputFaces=[], i=0) = (i >= len(inputFaces))
    ? [outputPoints, outputFaces]
    : let(
        inputFace=inputFaces[i],
        updatedPointsRes=sp_indexes_of_or_add(outputPoints, inputFace)
    ) polyhedron_from_faces(
        inputFaces,
        updatedPointsRes[1],
        concat(outputFaces, [updatedPointsRes[0]]),
        i + 1
    );

//////////////////////////////////////////////////////////////////////////////////////////////////
// Geodesicing

function mult_mat4_vec3(mat4, vec3) = let(
    vec4 = [vec3[0], vec3[1], vec3[2], 1],
    result = mat4 * vec4
) [result[0], result[1], result[2]];

function normalize_point(p, radius=1) = sph_to_cart(sphu_from_cart(p, radius));
function normalize_points(pts, radius=1) = [for (p=pts) normalize_point(p, radius)];

function point_between(a, b, f) = a + (b - a) * f;

function divide_line(a, b, n) = [
    for (i=[0:n+1]) point_between(a, b, i / (n+1))
];

function flatten(input, output, i=0) = 
    i >= len(input)
        ? []
        : concat(input[i], flatten(input, output, i+1));
    
function divide_face(
    points,
    face,
    divisionCount
) = (divisionCount == 0) ? [[points[face[0]], points[face[1]], points[face[2]]]] : let(
    leftPoints=divide_line(points[face[0]], points[face[1]], divisionCount),
    rightPoints=divide_line(points[face[0]], points[face[2]], divisionCount)
) flatten([
    for (i=[1 : divisionCount+1]) lattice_line(
        i==1 ? [leftPoints[0]]
             : divide_line(leftPoints[i-1], rightPoints[i-1], i-2),
        divide_line(leftPoints[i], rightPoints[i], i-1)
    )
]);

function lattice_line(
    shortPoints,
    longPoints
) = concat(
    [[longPoints[1], shortPoints[0], longPoints[0]]],
    len(longPoints)>2?[for (i=[1:len(shortPoints)-1]) [longPoints[i], shortPoints[i], shortPoints[i-1]]]:[],
    len(longPoints)>2?[for (i=[1:len(longPoints)-2]) [longPoints[i+1], shortPoints[i], longPoints[i] ]]:[],
    []
);
    
function geodesic_icos_faces(geodesicFactor=3, faceIndexes=[0:19], radius=1) = let(
    icos=icosahedron(radius),
    points=icos[0],
    faces=icos[1],
    raw=flatten([for(faceIndex=faceIndexes) divide_face(points, faces[faceIndex], geodesicFactor-1)]),
    geo=let(r=polyhedron_from_faces(raw)) [ normalize_points(r[0], radius), r[1] ]
) geo;
    
function select_poly_faces(poly, faceIndexes) = 
    polyhedron_from_faces([ for (i=faceIndexes) face_points(poly[0], poly[1][i]) ]);

function normal_poly_width_radius(poly, newRadius) = [
    normalize_points(poly[0], newRadius),
    poly[1]
];

function face_center(points, face) = [
    (points[face[0]][0]+points[face[1]][0]+points[face[2]][0])/3,
    (points[face[0]][1]+points[face[1]][1]+points[face[2]][1])/3,
    (points[face[0]][2]+points[face[1]][2]+points[face[2]][2])/3
];
    
function face_points(points, face) = [
    for (p=face) points[p]
];
    
function cart_rad(c) = sqrt(c[0]*c[0]+c[1]*c[1]+c[2]*c[2]);

module aimAt(p) {
    rotate([0, 0, atan2(p.y,p.x)])
        rotate([0, atan2(norm([p.x,p.y]), p.z), 0])
        children();
}

module vertexLabels(poly) {
    points=poly[0];
    
    for (i = [0:len(points)-1]) { 
        pt = points[i];
        translate(pt) aimAt(pt) {
            text(str(i), size=12, halign="center", valign="center", $fn=1);
            translate([0,-10,0]) rotate([0,90,0]) cylinder(r=1,h=15,center=true);
        }
    }
}

module faceLabels(poly) {
    points=poly[0];
    faces=poly[1];
    radius=sph_from_cart(points[0])[2];
    
    for (i = [0:len(faces)-1]) { 
        face=faces[i];
        
        placeOnFace(points, face) {
            translate([0,-10,0]) rotate([0,90,0]) cylinder(r=1,h=15,center=true);
            text(str(i), size=15, halign="center", valign="center", $fn=1);
        }
    }
}

module simplePolyhedron(poly) {
    polyhedron(poly[0], poly[1], 1);
}

module placeOnFace(points, face, invert = false) {
    placeOnFaceAtPoint(
        points,
        face,
        face_center(points, face),
        invert
    ) children();
}

function placeOnFaceAtPointMatrix(points, face, center, pointXAt) = let(
    expFaceP = face_points(points, face),
    expFace = [
        expFaceP[2],
        expFaceP[0],
        expFaceP[1],
    ],
    
    v1 = expFace[1] - expFace[0],
    v2 = expFace[2] - expFace[0],
            
    kPrime = -normalize_point(cross(v1, v2)),
    jPrime = normalize_point(expFace[0] - center),
    iPrime = normalize_point(cross(jPrime, kPrime)),
    
    //echo(cart_rad(cross(v1, v2)), cart_rad(center));
    
    transformMatrix = [
        [iPrime[0], jPrime[0], kPrime[0], center[0]],
        [iPrime[1], jPrime[1], kPrime[1], center[1]],
        [iPrime[2], jPrime[2], kPrime[2], center[2]],
        [0, 0, 0,  1]
    ]
) transformMatrix;

module placeOnFaceAtPoint(points, face, center, invert = false) {
    transformMatrix = placeOnFaceAtPointMatrix(points, face, center);
        
    // The mirror here handles the case where the vertices are in an order that would send Z towards, rather than away from, the origin
    if (invert) {
        mirror([0,0,cart_rad(transformMatrix*[0,0,1,1]) < cart_rad(transformMatrix*[0,0,0,1])?1:0])
            multmatrix(matrix_invert(transformMatrix)) {
            children();
            // # cylinder(r=2,h=50);
            // rotate([90,0,0]) cylinder(r=2,h=200,center=true);
        }
    } else {
        multmatrix(transformMatrix) 
            mirror([0,0,cart_rad(transformMatrix*[0,0,1,1]) < cart_rad(transformMatrix*[0,0,0,1])?1:0]) {
            children();
            // # cylinder(r=2,h=50);
            // rotate([90,0,0]) cylinder(r=2,h=200,center=true);
        }
    }
}

module placeOnFaces(poly) {
    points=poly[0];
    faces=poly[1];
    
    for (face = faces) {
        placeOnFace(points, face) children(); 
    }
}

simplePolyhedron();

geodesic_example.scad

include <geodesic.scad>

geo=geodesic_icos_faces(geodesicFactor=2, faceIndexes=[0:19], radius=150);
points=geo[0];
faces=geo[1];

polyhedron(points, faces, 1);

color("blue") vertexLabels(geo);
color("red") faceLabels(geo);

maths_geodesic.scad

/*
	License: This code is placed in the public Domain
	Contributed By: Willliam A Adams
	August 2011
*/

// A couple of useful constants
Cpi = 3.14159;
Cphi = 1.61803399;
Cepsilon = 0.00000001;

 
// Function: clean
//
// Parameters:
//	n - A number that might be very close to zero
// Description:  
//	There are times when you want a very small number to 
// 	just be zero, instead of being that very small number.
//	This function will compare the number to an arbitrarily small 
//	number.  If it is smaller than the 'epsilon', then zero will be 
// 	returned.  Otherwise, the original number will be returned.
//

function clean(n) = (n < 0) ? ((n < -Cepsilon) ? n : 0) : 
	(n < Cepsilon) ? 0 : n; 

// Function: safediv
//
// Parameters
//	n - The numerator
//	d - The denominator
//
// Description:
//	Since division by zero is generally not a desirable thing, safediv
//	will return '0' whenever there is a division by zero.  Although this will
//	mask some erroneous division by zero errors, it is often the case
//	that you actually want this behavior.  So, it makes it convenient.
function safediv(n,d) = (d==0) ? 0 : n/d;


//==================================
// Degrees
//==================================

function DEGREES(radians) = (180/Cpi) * radians;

function RADIANS(degrees) = Cpi/180 * degrees;

function deg(deg, min=0, sec=0) = [deg, min, sec];

function deg_to_dec(d) = d[0] + d[1]/60 + d[2]/60/60;


//==================================
//  Spherical coordinates
//==================================

// create an instance of a spherical coordinate
// long - rotation around z -axis
// lat - latitude, starting at 0 == 'north pole'
// rad - distance from center
function sph(long, lat, rad=1) = [long, lat, rad];

// Convert spherical to cartesian
function sph_to_cart(s) = [
	clean(s[2]*sin(s[1])*cos(s[0])),  
	clean(s[2]*sin(s[1])*sin(s[0])),
	clean(s[2]*cos(s[1]))
	];

// Convert from cartesian to spherical
function sph_from_cart(c) = sph(
	atan2(c[1],c[0]), 
	atan2(sqrt(c[0]*c[0]+c[1]*c[1]), c[2]), 
	sqrt(c[0]*c[0]+c[1]*c[1]+c[2]*c[2])
	);

function sphu_from_cart(c, rad=1) = sph(
	atan2(c[1],c[0]), 
	atan2(sqrt(c[0]*c[0]+c[1]*c[1]), c[2]), 
	rad
	);

// compute the chord distance between two points on a sphere
function sph_dist(c1, c2) = sqrt(
	c1[2]*c1[2] + c2[2]*c2[2] - 
	2*c1[2]*c2[2]*
	((cos(c1[1])*cos(c2[1])) + cos(c1[0]-c2[0])*sin(c1[1])*sin(c2[1]))   
	);


//==========================================
//	Geodesic calculations
// 
//  Reference: Geodesic Math and How to Use It
//  By: Hugh Kenner
//  Second Paperback Edition (2003), p.74-75
//  http://www.amazon.com/Geodesic-Math-How-Hugh-Kenner/dp/0520239318
//
//  The book was used for reference, so if you want to check the math, 
//  you can plug in various numbers to various routines and see if you get
//  the same numbers in the book.
//
//  In general, there are enough routines here to implement the various
//  pieces necessary to make geodesic objects.
//==========================================

function poly_sum_interior_angles(sides) = (sides-2)*180;
function poly_single_interior_angle(pq) = poly_sum_interior_angles(pq[0])/pq[0];



// Given a set of coordinates, return the frequency
// Simply calculated by adding up the values of the coordinates
function geo_freq(xyz) = xyz[0]+xyz[1]+xyz[2];

// Convert between the 2D coordinates of vertices on the face triangle
// to the 3D vertices needed to calculate spherical coordinates
function geo_tri2_tri3(xyf) = [xyf[1], xyf[0]-xyf[1], xyf[2]-xyf[0]];

// Given coordinates for a vertex on the octahedron face
// return the spherical coordinates for the vertex
// class 1, method 1
function octa_class1(c) = sph(
	atan(safediv(c[0], c[1])),
	atan(sqrt(c[0]*c[0]+c[1]*c[1])/c[2]),
	1 
	);

function octa_class2(c) = sph(
	atan(c[0]/c[1]),
	atan( sqrt( 2*(c[0]*c[0]+c[1]*c[1])) /c[2]),
	1 
	);

function icosa_class1(c) = octa_class1(
	[
		c[0]*sin(72),  
		c[1]+c[0]*cos(72),  
		geo_freq(c)/2+c[2]/Cphi
	]);

function icosa_class2(c) = sph(
	atan([c0]/c[1]), 
	atan(sqrt(c[0]*c[0]+c[1]*c[1]))/cos(36)*c[2],
	1
	);
 
function tetra_class1(c) = octa_class1(
	[
		sqrt(3*c[0]),  
		2*c[1]-c[0],  
		(3*c[2]-c[0]-c[1])/sqrt(2)
	]);

function class1_icosa_chord_factor(v1, v2, freq) = sph_dist( 
		icosa_class1(geo_tri2_tri3( [v1[0], v1[1], freq])),
		icosa_class1(geo_tri2_tri3( [v2[0], v2[1], freq]))
	);

nice_rendering.scad

// Icosahedron

// base coordinates
// source:  http://dmccooey.com/polyhedra/Icosahedron.txt
// generated by  http://kitwallace.co.uk/3d/solid-to-scad.xq
Name = "Icosahedron";
include <geodesic.scad>

poly=geodesic_icos_faces(geodesicFactor=2, faceIndexes=[0:19], radius=150);

points = poly[0];
faces = poly[1];
edges = [];
// --------------------------------- 


// cut holes out of shell
eps=0.02;
radius=150;
shell_ratio=0.3;
prism_base_ratio =0.6;
prism_height_ratio=0.8;
prism_scale=0.5;
nfaces = [];
scale=1;

//insert

spoints = normalize(centre_points(points),radius);
sfaces = lhs_faces(faces,spoints);
cfaces =  select_nsided_faces(sfaces,nfaces);

module shape() {
   difference() {
      polyhedron(spoints,sfaces);
      scale(1-shell_ratio) polyhedron(spoints,sfaces);
     face_prisms_in(cfaces,spoints,prism_base_ratio,prism_scale,prism_height_ratio);
   }
};


scale(scale)
 place_on_largest_face(sfaces,spoints)
   shape();
 


// ruler(10);

// functions for the construction of polyhedra
// chris wallace
// see http://kitwallace.tumblr.com/tagged/polyhedra for info


//  functions for creating the matrices for transforming a single point

function m_translate(v) = [ [1, 0, 0, 0],
                            [0, 1, 0, 0],
                            [0, 0, 1, 0],
                            [v.x, v.y, v.z, 1  ] ];
                            
function m_rotate(v) =  [ [1,  0,         0,        0],
                          [0,  cos(v.x),  sin(v.x), 0],
                          [0, -sin(v.x),  cos(v.x), 0],
                          [0,  0,         0,        1] ]
                      * [ [ cos(v.y), 0,  -sin(v.y), 0],
                          [0,         1,  0,        0],
                          [ sin(v.y), 0,  cos(v.y), 0],
                          [0,         0,  0,        1] ]
                      * [ [ cos(v.z),  sin(v.z), 0, 0],
                          [-sin(v.z),  cos(v.z), 0, 0],
                          [ 0,         0,        1, 0],
                          [ 0,         0,        0, 1] ];
                            
function vec3(v) = [v.x, v.y, v.z];
function transform(v, m)  = vec3([v.x, v.y, v.z, 1] * m);
                            
function matrix_to(p0, p) = 
                       m_rotate([0, atan2(sqrt(pow(p[0], 2) + pow(p[1], 2)), p[2]), 0]) 
                     * m_rotate([0, 0, atan2(p[1], p[0])]) 
                     * m_translate(p0);

function matrix_from(p0, p) = 
                      m_translate(-p0)
                      * m_rotate([0, 0, -atan2(p[1], p[0])]) 
                      * m_rotate([0, -atan2(sqrt(pow(p[0], 2) + pow(p[1], 2)), p[2]), 0]); 

function transform_points(list, matrix, i = 0) = 
    i < len(list) 
       ? concat([ transform(list[i], matrix) ], transform_points(list, matrix, i + 1))
       : [];


//  convert from point indexes to point coordinates

function as_points(indexes,points,i=0) =
     i < len(indexes) 
        ?  concat([points[indexes[i]]], as_points(indexes,points,i+1))
        : [];

//  basic vector functions
function normal_r(face) =
     cross(face[1]-face[0],face[2]-face[0]);

function normal(face) =
     - normal_r(face) / norm(normal_r(face));

function centre(points) = 
      vsum(points) / len(points);

// sum a list of vectors
function vsum(points,i=0) =  
      i < len(points)
        ?  (points[i] + vsum(points,i+1))
        :  [0,0,0];

function ssum(list,i=0) =  
      i < len(list)
        ?  (list[i] + ssum(list,i+1))
        :  0;


// add a vector to a list of vectors
function vadd(points,v,i=0) =
      i < len(points)
        ?  concat([points[i] + v], vadd(points,v,i+1))
        :  [];

function reverse_r(v,n) =
      n == 0 
        ? [v[0]]
        : concat([v[n]],reverse_r(v,n-1));

function reverse(v) = reverse_r(v, len(v)-1);

function sum_norm(points,i=0) =
    i < len(points)
       ?  norm(points[i]) + sum_norm(points,i+1)
       : 0 ;

function average_radius(points) =
       sum_norm(points) / len(points);


// select one dimension of a list of vectors
function slice(v,k,i=0) =
   i <len(v)
      ?  concat([v[i][k]], slice(v,k,i+1))
      : [];

function max(v, max=-9999999999999999,i=0) =
     i < len(v) 
        ?  v[i] > max 
            ?  max(v, v[i], i+1 )
            :  max(v, max, i+1 ) 
        : max;

function min(v, min=9999999999999999,i=0) =
     i < len(v) 
        ?  v[i] < min 
            ?  min(v, v[i], i+1 )
            :  min(v, min, i+1 ) 
        : min;

function project(pts,i=0) =
     i < len(pts)
        ? concat([[pts[i][0],pts[i][1]]], project(pts,i+1))
        : [];
        
function contains(n, list, i=0) =
     i < len(list) 
        ?  n == list[i]
           ?  true
           :  contains(n,list,i+1)
        : false;

// normalize the points to have origin at 0,0,0 
function centre_points(points) = 
     vadd(points, - centre(points));

//scale to average radius = radius
function normalize(points,radius) =
    points * radius /average_radius(points);
 
function select_nsided_faces(faces,nsides,i=0) =
  len(nsides) == 0
     ?  faces
     :  i < len(faces)
         ?  contains(len(faces[i]), nsides)
             ? concat([faces[i]],  select_nsided_faces(faces,nsides,i+1))
             : select_nsided_faces(faces,nsides,i+1)
         : [];
         
function longest_edge(face,max=-1,i=0) =
       i < len(face)
          ?  norm(face[i] - face[(i+1)% len(face)]) > max
             ?  longest_edge(face, norm(face[i] - face[(i+1)% len(face)]),i+1)
             :  longest_edge(face, max,i+1)
          : max ;

function point_edges(point,edges,i=0) =
    i < len(edges) 
       ? point == edges[i][0] || point == edges[i][1]
         ? concat([edges[i]], point_edges(point,edges,i+1))
         : point_edges(point,edges,i+1)
       : [];

function select_nedged_points(points,edges,nedges,i=0) =
     i < len(points) 
         ?  len(point_edges(i,edges)) == nedges
             ? concat([i],  select_nedged_points(points,edges,nedges,i+1))
             : select_nedged_points(points,edges,nedges,i+1)
         : [];

function triangle(a,b) = norm(cross(a,b))/2;

function face_area_centre(face,centre,i=0) =
    i < len(face)
       ?  triangle(
                face[i] - centre,
                face[(i+1) % len(face)] - centre)
          + face_area_centre(face,centre,i+1)
       : 0 ;

function face_area(face) = face_area_centre(face,centre(face));

function face_areas(faces,points,i=0) =
   i < len(faces)
      ? concat([[i,  face_area(as_points(faces[i],points))]] ,
               face_areas(faces,points,i+1))
      : [] ;
 
function max_area(areas, max=[-1,-1], i=0) =
   i <len(areas)
      ? areas[i][1] > max[1]
         ?  max_area(areas,areas[i],i+1)
         :  max_area(areas,max,i+1)
      : max;


function bbox(v) = [
   [min(slice(spoints,0)), max(slice(spoints,0))],
   [min(slice(spoints,1)), max(slice(spoints,1))],
   [min(slice(spoints,2)), max(slice(spoints,2))]
];

// check that all faces have a lhs orientation
function cosine_between(u, v) =(u * v) / (norm(u) * norm(v));

function lhs_faces(faces,points,i=0) =
     i < len(faces) 
        ?  cosine_between(normal(as_points(faces[i],points)),
                         centre(as_points(faces[i],points))) < 0
            ?  concat([reverse(faces[i])],lhs_faces(faces,points,i+1))
            :  concat([faces[i]],lhs_faces(faces,points,i+1))
        : [] ;


function fs(p) = f(p[0],p[1],p[2]);

function modulate_point(p) =
    spherical_to_xyz(fs(xyz_to_spherical(p)));

function modulate_points(points,i=0) =
   i < len(points)
      ? concat([modulate_point(points[i])],modulate_points(points,i+1))
      : [];

function xyz_to_spherical(p) =
    [ norm(p), acos(p.z/ norm(p)), atan2(p.x,p.y)] ;

function spherical_to_xyz_full(r,theta,phi) =
    [ r * sin(theta) * cos(phi),
      r * sin(theta) * sin(phi),
      r * cos(theta)];

function spherical_to_xyz(s) =
     spherical_to_xyz_full(s[0],s[1],s[2]);

function select_large_faces(faces,points, min,i=0) =
  i < len(faces)
     ?  face_area(as_points(faces[i],points)) > min
       ? concat([faces[i]],  select_large_faces(faces,points,min,i+1))
       :select_large_faces(faces,points,min,i+1)     
     : [];
        
function lower(char) =
    contains(char,"abcdefghijklmnopqrstuvwxyz") ;

function char_layer(char) =
    lower(char) 
         ? str(char,"_")
         : char;

module write_char(font,char) {
	linear_extrude(height=1,convexity=10) 
      import(file=str("write/",font,".dxf"),layer=char_layer(char));
};

module write_centred_char(font,char) {
	linear_extrude(height=1,convexity=10) 
      translate([-2.5,-4,0])
          import(file=str("write/",font,".dxf"),layer=char_layer(char));
};
module engrave_face_word(faces,points,word,font,ratio,thickness) {
    for (i=[0:len(faces) - 1]) 
      if (i <len(word))
        assign (f = as_points(faces[i],points)) 
        assign (n = normal(f), c = centre(f))
        assign(s = longest_edge(f) / 20* ratio)
           orient_to(c,n)  
                translate([0,0,-thickness+eps])
                     scale([s,s,thickness])
                          write_centred_char(font,word[i]);
}
                          
module orient_to(centre, normal) {   
      translate(centre)
      rotate([0, 0, atan2(normal[1], normal[0])]) //rotation
      rotate([0, atan2(sqrt(pow(normal[0], 2)+pow(normal[1], 2)),normal[2]), 0])
      children();
}

module orient_from(centre, normal) {   
      rotate([0, -atan2(sqrt(pow(normal[0], 2)+pow(normal[1], 2)),normal[2]), 0])
      rotate([0, 0, -atan2(normal[1], normal[0])]) //rotation
      translate(-centre)
      children();
}

module place_on_largest_face(faces,points) {
  assign (largest = max_area(face_areas(faces,points)))
  assign (lpoints = as_points(faces[largest[0]],points))
  assign (n = normal(lpoints),c = centre(lpoints))
  orient_from(c,-n)
  children();
}
              
module make_edge(edge, points, r) {
    assign(p0 = points[edge[0]], p1 = points[edge[1]])
    assign(v = p1 -p0 )
     orient_to(p0,v)
       cylinder(r=r, h=norm(v)); 
}

module make_edges(points, edges, r) {
   for (i =[0:len(edges)-1])
      make_edge(edges[i],points, r);
}

module make_vertices(points,r) { 
   for (i = [0:len(points)-1])
      translate(points[i]) sphere(r); 
}

module face_prism (face,prism_base_ratio,prism_scale,prism_height_ratio) {
    assign (n = normal(face), c= centre(face))
    assign (m = matrix_from(c,n))
    assign (tpts =  prism_base_ratio * transform_points(face,m))
    assign (max_length = longest_edge(face))
    assign (xy = project(tpts)) 
      linear_extrude(height=prism_height_ratio * max_length, scale=prism_scale) 
          polygon(points=xy);
}

module face_prisms_in(faces,points,prism_base_ratio,prism_scale,prism_height_ratio) {
    for (i=[0:len(faces) - 1]) 
       assign (f = as_points(faces[i],points)) 
       assign (n = normal(f), c = centre(f))
       orient_to(c,n) 
          translate([0,0,eps]) 
               mirror() rotate([0,180,0]) 
                   face_prism(f,prism_base_ratio,prism_scale,prism_height_ratio);
}

module face_prisms_out(faces,points,prism_base_ratio,prism_scale,prism_height_ratio) {
    for (i=[0:len(faces) - 1]) 
       assign (f = as_points(faces[i],points)) 
       assign (n = normal(f), c = centre(f))
       orient_to(c,n) 
          translate([0,0,-eps]) 
               face_prism(f,prism_base_ratio,prism_scale,prism_height_ratio);
}

module face_prisms_through(faces,points,prism_base_ratio,prism_scale,prism_height_ratio) {
    for (i=[0:len(faces) - 1]) 
       assign (f = as_points(faces[i],points)) 
       assign (n = normal(f), c = centre(f))
       orient_to(c,n) 
          translate([0,0,prism_height_ratio*longest_edge(f)/2]) 
               mirror() rotate([0,180,0]) 
                   face_prism(f,prism_base_ratio,prism_scale,prism_height_ratio);
}
module ruler(n) {
   for (i=[0:n-1]) 
       translate([(i-n/2 +0.5)* 10,0,0]) cube([9.8,5,2], center=true);
}

module ground(size=50) {
   translate([0,0,-size]) cube(2*size,center=true);
}

module cross_section(size=50) {
   translate([0,0,-size]) cube(2*size);
}

module ring(radius,thickness,height) {
   difference() {
      cylinder(h=height,r=radius);
      translate([0,0,-eps]) cylinder(h=height+2*eps,r=radius - thickness);
   }
}

platonic.scad

include <maths_geodesic.scad>


// Information about platonic solids 
// This information is useful in constructing the various solids
// can be found here: http://en.wikipedia.org/wiki/Platonic_solid
// V - vertices
// E - edges
// F - faces
// number, V, E, F, schlafli symbol, dihedral angle, element, name
//tetrahedron = [1, 4, 6, 4, [3,3], 70.5333, "fire", "tetrahedron"];
//hexahedron = [2, 8, 12, 6, [4,3], 90, "earth", "cube"];
//octahedron = [3, 6, 12, 8, [3,4], 109.467, "air", "air"];
//dodecahedron = [4, 20, 30, 12, [5,3], 116.565, "ether", "universe"];
//icosahedron = [5, 12, 30, 20, [3,5], 138.190, "water", "water"];

// Schlafli representation for the platonic solids
// Given this representation, we have enough information
// to derive a number of other attributes of the solids
tetra_sch = [3,3];
hexa_sch = [4,3];
octa_sch = [3,4];
dodeca_sch = [5,3];
icosa_sch = [3,5];

// Given the schlafli representation, calculate
// the number of edges, vertices, and faces for the solid
function plat_edges(pq) = (2*pq[0]*pq[1])/
	((2*pq[0])-(pq[0]*pq[1])+(2*pq[1]));
function plat_vertices(pq) = (2*plat_edges(pq))/pq[1];
function plat_faces(pq) = (2*plat_edges(pq))/pq[0];


// Calculate angular deficiency of each vertex in a platonic solid 
// p - sides
// q - number of edges per vertex
//function angular_defect(pq) = 360 - (poly_single_interior_angle(pq)*pq[1]);
function plat_deficiency(pq) = DEGREES(2*Cpi - pq[1]*Cpi*(1-2/pq[0]));

function plat_dihedral(pq) = 2 * asin( cos(180/pq[1])/sin(180/pq[0]));

function plat_circumradius(pq, a) = 
	(a/2)*
	tan(Cpi/pq[1])*
	tan(plat_dihedral(pq)/2);

function plat_midradius(pq, a) = 
	(a/2)*
	cot(Cpi/pq[0])*
	tan(plat_dihedral(pq)/2);

function plat_inradius(pq,a) = 
	a/(2*tan(DEGREES(Cpi/pq[0])))*
	sqrt((1-cos(plat_dihedral(pq)))/(1+cos(plat_dihedral(pq))));

//================================================
//	Tetrahedron
//================================================
tetra_cart = [
	[+1, +1, +1],
	[-1, -1, +1],
	[-1, +1, -1],
	[+1, -1, -1]
];

function tetra_unit(rad=1) = [
	sph_to_cart(sphu_from_cart(tetra_cart[0], rad)), 
	sph_to_cart(sphu_from_cart(tetra_cart[1], rad)),
	sph_to_cart(sphu_from_cart(tetra_cart[2], rad)),
	sph_to_cart(sphu_from_cart(tetra_cart[3], rad)),
	];


tetrafaces = [
	[0, 3, 1],
	[0,1,2],
	[2,1,3],
	[0,2,3]
];

tetra_edges = [
	[0,1],
	[0,2],
	[0,3], 
	[1,2], 
	[1,3], 
	[2,3],	
	];

function tetrahedron(rad=1) = [tetra_unit(rad), tetrafaces, tetra_edges];


//================================================
//	Hexahedron - Cube 
//================================================
// vertices for a unit cube with sides of length 1
hexa_cart = [
	[0.5, 0.5, 0.5], 
	[-0.5, 0.5, 0.5], 
	[-0.5, -0.5, 0.5], 
	[0.5, -0.5, 0.5],
	[0.5, 0.5, -0.5], 
	[-0.5, 0.5, -0.5], 
	[-0.5, -0.5, -0.5], 
	[0.5, -0.5, -0.5],
];

function hexa_unit(rad=1) = [
	sph_to_cart(sphu_from_cart(hexa_cart[0], rad)), 
	sph_to_cart(sphu_from_cart(hexa_cart[1], rad)),
	sph_to_cart(sphu_from_cart(hexa_cart[2], rad)),
	sph_to_cart(sphu_from_cart(hexa_cart[3], rad)),
	sph_to_cart(sphu_from_cart(hexa_cart[4], rad)), 
	sph_to_cart(sphu_from_cart(hexa_cart[5], rad)), 
	sph_to_cart(sphu_from_cart(hexa_cart[6], rad)),
	sph_to_cart(sphu_from_cart(hexa_cart[7], rad)),
	];


// enumerate the faces of a cube
// vertex order is clockwise winding
hexafaces = [
	[0,3,2,1],	// top
	[0,1,5,4],
	[1,2,6,5],
	[2,3,7,6],
	[3,0,4,7],
	[4,5,6,7],	// bottom
];

hexa_edges = [
	[0,1],
	[0,3], 
	[0,4], 
	[1,2],
	[1,5],
	[2,3],
	[2,6], 
	[3,7], 
	[4,5], 	
	[4,7], 
	[5,4],
	[5,6], 
	[6,7], 
	];


function hexahedron(rad=1) =[hexa_unit(rad), hexafaces, hexa_edges];


//================================================
//	Octahedron 
//================================================

octa_cart = [
	[+1, 0, 0],  // + x axis
	[-1, 0, 0],	// - x axis
	[0, +1, 0],	// + y axis
	[0, -1, 0],	// - y axis
	[0, 0, +1],	// + z axis
	[0, 0, -1] 	// - z axis
];

function octa_unit(rad=1) = [
	sph_to_cart(sphu_from_cart(octa_cart[0], rad)), 
	sph_to_cart(sphu_from_cart(octa_cart[1], rad)),
	sph_to_cart(sphu_from_cart(octa_cart[2], rad)),
	sph_to_cart(sphu_from_cart(octa_cart[3], rad)),
	sph_to_cart(sphu_from_cart(octa_cart[4], rad)), 
	sph_to_cart(sphu_from_cart(octa_cart[5], rad)), 
	];

octafaces = [
	[4,2,0],
	[4,0,3],
	[4,3,1],
	[4,1,2],
	[5,0,2],
	[5,3,0],
	[5,1,3],
	[5,2,1]
	];

octa_edges = [
	[0,2], 
	[0,3],
	[0,4],
	[0,5],
	[1,2],
	[1,3],
	[1,4],
	[1,5],
	[2,4], 
	[2,5],
	[3,4],
	[3,5],
	];

function octahedron(rad=1) = [octa_unit(rad), octafaces, octa_edges];

//================================================
//	Dodecahedron
//================================================
// (+-1, +-1, +-1)
// (0, +-1/Cphi, +-Cphi)
// (+-1/Cphi, +-Cphi, 0)
// (+-Cphi, 0, +-1/Cphi)

dodeca_cart = [
	[+1, +1, +1],			// 0, 0
	[+1, -1, +1],			// 0, 1
	[-1, -1, +1],			// 0, 2
	[-1, +1, +1],			// 0, 3

	[+1, +1, -1],			// 1, 4
	[-1, +1, -1],			// 1, 5
	[-1, -1, -1],			// 1, 6
	[+1, -1, -1],			// 1, 7

	[0, +1/Cphi, +Cphi],		// 2, 8
	[0, -1/Cphi, +Cphi],		// 2, 9
	[0, -1/Cphi, -Cphi],		// 2, 10
	[0, +1/Cphi, -Cphi],		// 2, 11

	[-1/Cphi, +Cphi, 0],		// 3, 12
	[+1/Cphi, +Cphi, 0],		// 3, 13
	[+1/Cphi, -Cphi, 0],		// 3, 14
	[-1/Cphi, -Cphi, 0],		// 3, 15

	[-Cphi, 0, +1/Cphi],		// 4, 16
	[+Cphi, 0, +1/Cphi],		// 4, 17
	[+Cphi, 0, -1/Cphi],		// 4, 18
	[-Cphi, 0, -1/Cphi],		// 4, 19
];

function dodeca_unit(rad=1) = [
	sph_to_cart(sphu_from_cart(dodeca_cart[0], rad)), 
	sph_to_cart(sphu_from_cart(dodeca_cart[1], rad)),
	sph_to_cart(sphu_from_cart(dodeca_cart[2], rad)),
	sph_to_cart(sphu_from_cart(dodeca_cart[3], rad)),
	sph_to_cart(sphu_from_cart(dodeca_cart[4], rad)), 
	sph_to_cart(sphu_from_cart(dodeca_cart[5], rad)), 
	sph_to_cart(sphu_from_cart(dodeca_cart[6], rad)),
	sph_to_cart(sphu_from_cart(dodeca_cart[7], rad)),
	sph_to_cart(sphu_from_cart(dodeca_cart[8], rad)),
	sph_to_cart(sphu_from_cart(dodeca_cart[9], rad)), 
	sph_to_cart(sphu_from_cart(dodeca_cart[10], rad)), 
	sph_to_cart(sphu_from_cart(dodeca_cart[11], rad)),
	sph_to_cart(sphu_from_cart(dodeca_cart[12], rad)),
	sph_to_cart(sphu_from_cart(dodeca_cart[13], rad)),
	sph_to_cart(sphu_from_cart(dodeca_cart[14], rad)), 
	sph_to_cart(sphu_from_cart(dodeca_cart[15], rad)), 
	sph_to_cart(sphu_from_cart(dodeca_cart[16], rad)),
	sph_to_cart(sphu_from_cart(dodeca_cart[17], rad)),
	sph_to_cart(sphu_from_cart(dodeca_cart[18], rad)),
	sph_to_cart(sphu_from_cart(dodeca_cart[19], rad)), 
	];



// These are the pentagon faces
// but CGAL has a problem rendering if things are 
// not EXACTLY coplanar
// so use the triangle faces instead
//dodeca_faces=[ 
//	[1,9,8,0,17],
//	[9,1,14,15,2],
//	[9,2,16,3,8],
//	[8,3,12,13,0],
//	[0,13,4,18,17],
//	[1,17,18,7,14],
//	[15,14,7,10,6],
//	[2,15,6,19,16],
//	[16,19,5,12,3],
//	[12,5,11,4,13],
//	[18,4,11,10,7],
//	[19,6,10,11,5]
//	];
dodeca_faces = [
	[1,9,8], 
	[1,8,0],
	[1,0,17],
	
	[9,1,14],
	[9,14,15],
	[9,15,2],
	
	[9,2,16],
	[9,16,3],
	[9,3,8],
	
	[8,3,12],
	[8,12,13],
	[8,13,0],
	
	[0,13,4],
	[0,4,18],
	[0,18,17],
	
	[1,17,18],
	[1,18,7],
	[1,7,14],
	
	[15,14,7],
	[15,7,10],
	[15,10,6],
	
	[2,15,6],
	[2,6,19],
	[2,19,16],
	
	[16,19,5],
	[16,5,12],
	[16,12,3],
	
	[12,5,11],
	[12,11,4],
	[12,4,13],
	
	[18,4,11],
	[18,11,10],
	[18,10,7],
	
	[19,6,10],
	[19,10,11],
	[19,11,5]
	];

dodeca_edges=[
	[0,8],
	[0,13],
	[0,17],

	[1,9],
	[1,14],
	[1,17],

	[2,9],
	[2,15],
	[2,16],

	[3,8],
	[3,12],
	[3,16],

	[4,11],
	[4,13],
	[4,18],

	[5,11],
	[5,12],
	[5,19],

	[6,10],
	[6,15],
	[6,19],

	[7,10],
	[7,14],
	[7,18],

	[8,9],
	[10,11],
	[12,13],
	[14,15],
	[16,19],
	[17,18],
	];

function dodecahedron(rad=1) = [dodeca_unit(rad), dodeca_faces, dodeca_edges];

//================================================
//	Icosahedron
//================================================
//
// (0, +-1, +-Cphi)
// (+-Cphi, 0, +-1)
// (+-1, +-Cphi, 0)

icosa_cart = [
	[0, +1, +Cphi],	// 0
	[0, +1, -Cphi],	// 1
	[0, -1, -Cphi],	// 2
	[0, -1, +Cphi],	// 3

	[+Cphi, 0, +1],	// 4
	[+Cphi, 0, -1],	// 5
	[-Cphi, 0, -1],	// 6
	[-Cphi, 0, +1],	// 7

	[+1, +Cphi, 0],	// 8
	[+1, -Cphi, 0],	// 9
	[-1, -Cphi, 0],	// 10
	[-1, +Cphi, 0]	// 11
	];

function icosa_sph(rad=1) = [
	sphu_from_cart(icosa_cart[0], rad), 
	sphu_from_cart(icosa_cart[1], rad),
	sphu_from_cart(icosa_cart[2], rad),
	sphu_from_cart(icosa_cart[3], rad),
	sphu_from_cart(icosa_cart[4], rad), 
	sphu_from_cart(icosa_cart[5], rad), 
	sphu_from_cart(icosa_cart[6], rad),
	sphu_from_cart(icosa_cart[7], rad),
	sphu_from_cart(icosa_cart[8], rad),
	sphu_from_cart(icosa_cart[9], rad), 
	sphu_from_cart(icosa_cart[10], rad), 
	sphu_from_cart(icosa_cart[11], rad),
	];

function icosa_unit(rad=1) = [
	sph_to_cart(sphu_from_cart(icosa_cart[0], rad)), 
	sph_to_cart(sphu_from_cart(icosa_cart[1], rad)),
	sph_to_cart(sphu_from_cart(icosa_cart[2], rad)),
	sph_to_cart(sphu_from_cart(icosa_cart[3], rad)),
	sph_to_cart(sphu_from_cart(icosa_cart[4], rad)), 
	sph_to_cart(sphu_from_cart(icosa_cart[5], rad)), 
	sph_to_cart(sphu_from_cart(icosa_cart[6], rad)),
	sph_to_cart(sphu_from_cart(icosa_cart[7], rad)),
	sph_to_cart(sphu_from_cart(icosa_cart[8], rad)),
	sph_to_cart(sphu_from_cart(icosa_cart[9], rad)), 
	sph_to_cart(sphu_from_cart(icosa_cart[10], rad)), 
	sph_to_cart(sphu_from_cart(icosa_cart[11], rad)),
	];

icosa_faces = [ 
	[3,0,4],
	[3,4,9],
	[3,9,10],
	[3,10,7],
	[3,7,0],
	[0,8,4],
	[0,7,11],
	[0,11,8],
	[4,8,5],
	[4,5,9],
	[7,10,6],
	[7,6,11],
	[9,5,2],
	[9,2,10],
	[2,6,10],
	[1,5,8],
	[1,8,11],
	[1,11,6],
	[5,1,2],
	[2,1,6]
	];

icosa_edges = [
	[0,3],
	[0,4],
	[0,7],
	[0,8],
	[0,11],
	[1,5],
	[1,8],
	[1,11],
	[1,6],
	[1,2],
	[2,5],
	[2,6],
	[2,9],
	[2,10],
	[3,4],
	[3,9],
	[3,10],
	[3,7],
	[4,5],
	[4,8],
	[4,9],
	[5,8],
	[5,9],
	[6,7],		
	[6,10],
	[6,11],
	[7,10],
	[7,11],
	[8,11],
	[9,10],
	];

function icosahedron(rad=1) = [icosa_unit(rad), icosa_faces, icosa_edges];