### icosa medusa

turtling is fun... :-)
spherical tiling based on icosahedron (reminds me of medusa's hairdo - hence the name)
and again using reinder's occlusion methods from "Cubic space division #2"

```// created by florian berger (flockaroo) - 2018

// using reinder's occlusion magic from  "Cubic space division #2"

Canvas.setpenopacity(1.);

// Global code will be evaluated once.
const turtle = new Turtle();
const polygonList = [];

const PI2 = Math.PI*2.0;

function mcos(x) {
return Math.cos(x);
}

function msin(x) {
return Math.sin(x);
}

function cos2(x) {
return [Math.cos(x[0]),Math.cos(x[1])];
}

function sin2(x) {
return [Math.sin(x[0]),Math.sin(x[1])];
}

function SC(x) {
return [Math.sin(x),Math.cos(x)];
}

return [a[0]+b[0],a[1]+b[1]];
}

return [a[0]+b[0],a[1]+b[1],a[2]+b[2]];
}

function sub3(a,b) {
return [a[0]-b[0],a[1]-b[1],a[2]-b[2]];
}

function dot3(a,b) {
return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];
}

function scale2(a,b) {
return [a[0]*b,a[1]*b];
}

function scale3(a,b) {
return [a[0]*b,a[1]*b,a[2]*b];
}

function scale4(a,b) {
return [a[0]*b,a[1]*b,a[2]*b,a[3]*b];
}

function mul4(a,b) {
return [a[0]*b[0],a[1]*b[1],a[2]*b[2],a[3]*b[3]];
}

function mymix(a,b,f) {
return a*(1.0-f)+b*f;
}

function mymix22(a,b,f) {
return [a[0]*(1.0-f[0])+b[0]*f[0],a[1]*(1.0-f[1])+b[1]*f[1]];
}

function mix3(a,b,f) {
}

function length2(a) {
return Math.sqrt(a[0]*a[0]+a[1]*a[1]);
}

function length3(a) {
return Math.sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]);
}

function normalize3(a) {
return scale3(a,1.0/length3(a));
}

function cross(a,b) {
return [
a[1]*b[2]-b[1]*a[2],
a[2]*b[0]-b[2]*a[0],
a[0]*b[1]-b[0]*a[1]
];
}

const G=(.5+Math.sqrt(5./4.));
//const PI2=(3.141592653*2.);
const PI=3.141592653;

function fract(a) {return a-Math.floor(a);}
function floor2(a) { return [Math.floor(a[0]),Math.floor(a[1])];}
function fract2(a) { return [fract(a[0]),fract(a[1])];}

// noise funcs by Morgan McGuire https://www.shadertoy.com/view/4dS3Wd

function hash(n) { return fract(Math.sin(n) * 1.0e4); }
function hash2(p) { return fract(1.0e4 * Math.sin(17.0 * p[0] + p[1] * 0.1) * (0.1 + Math.abs(Math.sin(p[1] * 13.0 + p[0])))); }

function noise(x) {
var i = Math.floor(x);
var f = fract(x);
var u = f * f * (3.0 - 2.0 * f);
return mymix(hash(i), hash(i + 1.0), u);
}

function noise2(x) {
var i = floor2(x);
var f = fract2(x);
return mymix(r1, r2, f[1]);
}

function getRand01Sph(pos)
{
var res = [1024,1024];
//var texc=((pos.xy*123.+pos.z)*res+.5)/res;
var texc=[ pos[0]*123.+pos[2]+.5/res[0],
pos[1]*123.+pos[2]+.5/res[1] ];
var n=1.0-noise2(scale2(texc,256.));
return [0,0,0,0];
return [n,n,n,n];
}

//const vec4 p0 = vec4( 1, G, -G ,0 )/length(vec2(1,G));

const pI = scale4([ 1, G, -G ,0 ],1.0/length2([1,G]));

/*vec3 icosaPosRaw[12] = vec3[] (
-p0.xwz,  p0.xwy, -p0.xwy,  p0.xwz,
p0.wyx, -p0.wzx,  p0.wzx, -p0.wyx,
p0.yxw,  p0.zxw, -p0.zxw, -p0.yxw
);*/

const icosaPosRaw = [
[-pI[0],-pI[3],-pI[2]/*.xwz*/],
[ pI[0], pI[3], pI[1]/*.xwy*/],
[-pI[0],-pI[3],-pI[1]/*.xwy*/],
[ pI[0], pI[3], pI[2]/*.xwz*/],

[ pI[3], pI[1], pI[0]/*.wyx*/],
[-pI[3],-pI[2],-pI[0]/*.wzx*/],
[ pI[3], pI[2], pI[0]/*.wzx*/],
[-pI[3],-pI[1],-pI[0]/*.wyx*/],

[ pI[1], pI[0], pI[3]/*.yxw*/],
[ pI[2], pI[0], pI[3]/*.zxw*/],
[-pI[2],-pI[0],-pI[3]/*.zxw*/],
[-pI[1],-pI[0],-pI[3]/*.yxw*/]
];

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

// get icosahedron triangle
function getIcosaTri(idx)
{
var i1 = posIdx[(idx%20)*3+0];
var i2 = posIdx[(idx%20)*3+1];
var i3 = posIdx[(idx%20)*3+2];

var p1=icosaPosRaw[i1];
var p2=icosaPosRaw[i2];
var p3=icosaPosRaw[i3];
return [p1,p2,p3];
}

// subdivide 1 triangle into 4 triangles and give back closest triangle
function getTriSubDiv(idx, p1, p2, p3)
{

if     (idx==0) { p1=p1; p2=p4; p3=p6; }
else if(idx==1) { p1=p6; p2=p5; p3=p3; }
else if(idx==2) { p1=p6; p2=p4; p3=p5; }
else if(idx==3) { p1=p4; p2=p2; p3=p5; }
return [p1, p2, p3];
}

var triStripIndex = [0,1,2,1,3,2];

function mixSq3(a,b,f) { return mymix3(a,b,Math.cos(f*PI)*.5+.5); }

function geomTangentCurve(pos1, pos2, tan1, tan2, r1, r2,
rSegNum, tSegNum, vIdx)
{
var l = length3(sub3(pos1,pos2));
l*=.4;
var i=Math.floor(vIdx/3/2)%tSegNum;
//{  // converted some loops into proper vertex index values
var fact, fact2;
fact=Math.max(0.,(i)/(tSegNum)); // force >=0 because of sqrt below
fact=-Math.cos(fact*PI)*.5+.5; // more homogeneous steps
fact2=Math.max(0.,(i+1)/(tSegNum)); // force >=0 because of sqrt below
fact2=-Math.cos(fact2*PI)*.5+.5; // more homogeneous steps

var ta = mix3(tan1,tan2,fact);
var tn = mix3(tan1,tan2,fact2);

var dph=PI*2./(rSegNum);
//vec3 b1=normalize(vec3(ta.x,-ta.y,0));
var b1=normalize3(cross(ta,p1));
var b2=normalize3(cross(ta,b1));
//vec3 b3=normalize(vec3(tn.x,-tn.y,0));
var b3=normalize3(cross(tn,p2));
var b4=normalize3(cross(tn,b3));
var r_1 = mymix(r1,r2,fact);
var r_2 = mymix(r1,r2,fact2);
var j=Math.floor(vIdx/3/2/tSegNum)%rSegNum;
//{
var ph  = (j)*dph;
var ph2 = ph+dph;
var v = [v1,v2,v3,v4];
var pos = v[triStripIndex[vIdx%6]];
var normal = normalize3(cross(sub3(v[1],v[0]),sub3(v[2],v[0])));
//}
//}
return [pos,normal]
}

function calcAngle(v1, v2)
{
return Math.acos(dot3(v1,v2)/length3(v1)/length3(v2));
}

// distance to 2 torus segments in a triangle
// each torus segment spans from the middle of one side to the middle of another side
function geomTruchet(p1, p2, p3, dz, rSegNum, tSegNum, trNum,
{

var d = 10000.0;
// random rotation of torus-start-edges
if      (rnd>.75) { var d=p1; p1=p2; p2=d; }
else if (rnd>.50) { var d=p1; p1=p3; p3=d; }
else if (rnd>.25) { var d=p2; p2=p3; p3=d; }

var p4 = scale3(p1,(1.-dz));
var p5 = scale3(p2,(1.-dz));
var p6 = scale3(p3,(1.-dz));

// FIXME: why is this necessary - very seldom actually!?
var xchg=false;
if(dot3(cross(sub3(p2,p1),sub3(p3,p1)),p1)>0.0) {
var dummy;
dummy=p2; p2=p3; p3=dummy;
dummy=p5; p5=p6; p6=dummy;
xchg=true;
}

var lp1 = length3(p1);
var lp4 = length3(p4);

var r,r1,r2,fact,ang,fullAng;
var n = normalize3(cross(sub3(p2,p1),sub3(p3,p1)));

// torus segments:
// actually i have to fade from one torus into another
// because not all triangles are equilateral
var m;
//    std::vector <vec3> p;
var v1,v2,v3,v4,v5,v6;
var tubeNum=rSegNum*tSegNum*2*3;
var i=Math.floor(idx/(tubeNum))%trNum;
{
if(i==0) { v1=p1; v2=p2; v3=p3; v4=p4; v5=p5; v6=p6; }
if(i==1) { v1=p2; v2=p3; v3=p1; v4=p5; v5=p6; v6=p4; }
if(i==2) { v1=p3; v2=p1; v3=p2; v4=p6; v5=p4; v6=p5; }
//if(dot(cross(v2-v1,v3-v1),v1)>0.0) { vec3 dummy=v2; v2=v3; v3=dummy; }
//if(dot(cross(v5-v4,v6-v4),v4)>0.0) { vec3 dummy=v5; v5=v6; v6=dummy; }

fullAng = calcAngle(sub3(v3,v1),sub3(v2,v1));
//ang = calcAngle(pos2-v1,v2-v1);
var dang=fullAng/(tSegNum);
//if (fullAng<.001) break;

//float r1, r2;
//r1=length(v2-v1)*.5f; r1=length(v3-v1)*.5f;
var pos1, pos2, pos3;
// FIXME: why is this necessary - very seldom actually!? - see above
if(xchg) {
}
if(rnd2>.25)
{
}
var tan1 = normalize3(cross(sub3(v2,v1),v1));
var tan2 = normalize3(cross(sub3(v3,v1),v1));
idx%tubeNum);

return pn;
}
}

// final shape
function geom_medusa(rNum, tNum, subdiv, idx)
{
var p1,p2,p3;

var icosaFaceNum = 20;
var subDivNum = 4;

var trNum = 3; // tubes per truchet segemnt
var truchetNum=rNum*tNum*2*3*trNum; // 2 triangles * 3 vertices * trNum tubes

//for(int i1=0;i1<icosaFaceNum;i1++)
var idiv=truchetNum; for(var i=0;i<subdiv;i++) idiv*=subDivNum;
var pi = getIcosaTri(Math.floor(idx/idiv));
var p_subDivNum_i = 1;
for(var i=0;i<subdiv;i++)
{
idiv=Math.floor(idiv/subDivNum);
var isub = Math.floor(idx/idiv)%subDivNum;
pi = getTriSubDiv(isub,pi[0],pi[1],pi[2]);
p_subDivNum_i*=subDivNum;
}
var pn = geomTruchet(pi[0],pi[1],pi[2],0.12/(1+subdiv),rNum,tNum,trNum,-1.,idx%truchetNum);

return pn;
}

function medusaTri(idx)
{
var pos = [0,0,0];
var normal = [0,0,0];
var pn;
pn=geom_medusa(8,6,0,idx);
//pn=geomTangentCurve([-2,0,0], [0,0,2], [1,0,0], [0,0,-1], .1, .1, 10, 10, idx);
pos=pn[0];
//pn = getIcosaTri(Math.floor(idx/3));
//pos=pn[idx%3];
return pos;
}

function rotX(ph,v) {
return [ v[0],v[1]*mcos(ph)+v[2]*msin(ph), v[2]*mcos(ph)-v[1]*msin(ph) ];
}

function rotY(ph,v) {
return [ v[0]*mcos(ph)+v[2]*msin(ph), v[1], v[2]*mcos(ph)-v[0]*msin(ph) ];
}

function project(p)
{
p[2]+=180;
return [p[0]/p[2]*180.,p[1]/p[2]*180.,p[2]];
}

{
var z = p0[2]+p1[2]+p2[2]+p3[2];
var idx=0;

// hmm, why is the one below not working... !?
//    if(quads[i+8]>z) { idx=i; break; }
//}
quads.splice(idx, 0, p0[0], p0[1], p1[0], p1[1], p2[0], p2[1], p3[0], p3[1], z);
}

function walk(i) {
var num = 8*6*3*20;
if(i==0){
for(let j=0;j<num;j++) {
var p0=medusaTri(j*6);
var p1=medusaTri(j*6+1);
var p2=medusaTri(j*6+2);
var p3=medusaTri(j*6+4);
p0=scale3(p0,70.0);
p1=scale3(p1,70.0);
p2=scale3(p2,70.0);
p3=scale3(p3,70.0);
p0=rotX(1.4,p0);
p1=rotX(1.4,p1);
p2=rotX(1.4,p2);
p3=rotX(1.4,p3);
p0=rotY(.1,p0);
p1=rotY(.1,p1);
p2=rotY(.1,p2);
p3=rotY(.1,p3);
p0=project(p0);
p1=project(p1);
p2=project(p2);
p3=project(p3);

if(cross(sub3(p1,p0),sub3(p2,p0))[2]<0.0)
{
}
}
}

const p = new Polygon();
p.cp.push([p0[0], p0[1]]);
p.cp.push([p1[0], p1[1]]);
p.cp.push([p2[0], p2[1]]);
p.cp.push([p3[0], p3[1]]);
drawPolygon(turtle, p);
/*turtle.penup();
turtle.goto(p0);
turtle.pendown();
turtle.goto(p1);
turtle.goto(p2);
turtle.goto(p0);*/
/*turtle.goto(p3);
turtle.goto(p2);*/
return i <= num/1.8;
}

////////////////////////////
// reinder's occlusion code parts from "Cubic space division #2"
////////////////////////////

function drawPolygon(turtle, p) {
let vis = true;
for (let j=0; j<polygonList.length; j++) {
if(!p.boolean(polygonList[j])) {
vis = false;
break;
}
}
if (vis) {
p.draw(turtle, 0);
polygonList.push(p);
}
}

// polygon functions
function LineSegment(p1, p2) {
this.p1 = p1;
this.p2 = p2;
}
function Polygon() {
this.cp = []; // clip path: array of [x,y] pairs
this.dp = []; // 2d line to draw: array of linesegments
}
for (let i=s, l=this.cp.length; i<l; i++) {
this.dp.push(new LineSegment(this.cp[i], this.cp[(i+1)%l]));
}
}
Polygon.prototype.createPoly = function(x,y,c,r,a) {
this.cp = [];
for (let i=0; i<c; i++) {
this.cp.push( [x + Math.sin(i*Math.PI*2/c+a) * r, y + Math.cos(i*Math.PI*2/c+a) * r] );
}
}
Polygon.prototype.draw = function(t, inp=0) {
if (this.dp.length ==0) {
return;
}
for (let i=0, l=this.dp.length; i<l; i++) {
const d = this.dp[i];
if (!vec2_equal(d.p1, t.pos())) {
t.penup();
t.goto([d.p1[0]+inp*(Math.random()-.5), d.p1[1]+inp*(Math.random()-.5)]);
t.pendown();
}
t.goto([d.p2[0]+inp*(Math.random()-.5), d.p2[1]+inp*(Math.random()-.5)]);
}
}
Polygon.prototype.inside = function(p) {
// find number of i ntersection points from p to far away
const p1 = [0.1, -1000];
let int = 0;
for (let i=0, l=this.cp.length; i<l; i++) {
if (vec2_find_segment_intersect(p, p1, this.cp[i], this.cp[(i+1)%l])) {
int ++;
}
}
return int & 1;
}
Polygon.prototype.boolean = function(p, diff = true) {
// very naive polygon diff algorithm - made this up myself
const ndp = [];
for (let i=0, l=this.dp.length; i<l; i++) {
const ls = this.dp[i];

// find all intersections with clip path
const int = [];
for (let j=0, cl=p.cp.length; j<cl; j++) {
const pint = vec2_find_segment_intersect(ls.p1,ls.p2,p.cp[j],p.cp[(j+1)%cl]);
if (pint) {
int.push(pint);
}
}
if (int.length == 0) { // 0 intersections, inside or outside?
if (diff == !p.inside(ls.p1)) {
ndp.push(ls);
}
} else {
int.push(ls.p1);
int.push(ls.p2);
// order intersection points on line ls.p1 to ls.p2
const cmp = [ls.p2[0]-ls.p1[0], ls.p2[1]-ls.p1[1]];
int.sort( (a,b) => {
const db = vec2_dot([b[0]-ls.p1[0], b[1]-ls.p1[1]], cmp);
const da = vec2_dot([a[0]-ls.p1[0], a[1]-ls.p1[1]], cmp);
return da - db;
});
for (let j=0; j<int.length-1; j++) {
if (!vec2_equal(int[j], int[j+1])) {
if (diff == !p.inside([(int[j][0]+int[j+1][0])/2,(int[j][1]+int[j+1][1])/2])) {
ndp.push(new LineSegment(int[j], int[j+1]));
}
}
}
}
}
this.dp = ndp;
return this.dp.length > 0;
}

// vec functions
const vec2_equal = (a,b) => vec2_dist_sqr(a,b) < 0.01;
const vec2_dot = (a, b) => a[0]*b[0]+a[1]*b[1];
const vec2_dist_sqr = (a, b) => (a[0]-b[0])*(a[0]-b[0]) + (a[1]-b[1])*(a[1]-b[1]);
//port of http://paulbourke.net/geometry/pointlineplane/Helpers.cs
function vec2_find_segment_intersect(l1p1, l1p2, l2p1, l2p2) {
const d = (l2p2[1] - l2p1[1]) * (l1p2[0] - l1p1[0]) - (l2p2[0] - l2p1[0]) * (l1p2[1] - l1p1[1]);
const n_a = (l2p2[0] - l2p1[0]) * (l1p1[1] - l2p1[1]) - (l2p2[1] - l2p1[1]) * (l1p1[0] - l2p1[0]);
const n_b = (l1p2[0] - l1p1[0]) * (l1p1[1] - l2p1[1]) - (l1p2[1] - l1p1[1]) * (l1p1[0] - l2p1[0]);
if (d == 0) {
return false;
}
const ua = n_a / d;
const ub = n_b / d;
if (ua >= 0 && ua <= 1 && ub >= 0 && ub <= 1) {
return [l1p1[0] + (ua * (l1p2[0] - l1p1[0])), l1p1[1] + (ua * (l1p2[1] - l1p1[1])) ];
}
return false;
}
```