### SH Shapes

Spherical Harmonic Shapes

```// You can find the Turtle API reference here: https://turtletoy.net/syntax
Canvas.setpenopacity(0.1);

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

const PI = Math.PI;
const PI2 = Math.PI*2;
const sqrt2 = Math.sqrt(2);

const RANDOM = 0.1;

const WIREFRAME = 1

// Change COEF to combine the linear combination of
let COEF = [
0.1,
0.1,0.1,0.1,
1,0.1,0.1,0.1,1,
0.01,0.1,0.1,0.01,0.01,0.01,0.1,
0.1,0.001,0.001,0.001,0.001,0.001,0.001,0.001,1,
0.01,0.01,0.01,0.01,0.001,0.01,0.01,0.01,0.001,0.001,0.001,
1,0.001,0.001,0.001,0.0001,0.001,0.0001,0.001,0.0001,0.001,0.0001,0.001,1,
//0.0001,0.0001,0.0001,0.0001,1,0.0001,0.0001,0.0001,0.0001,0.0001,0.0001,0.0001,0.0001,0.0001,0.0001
];

const LEVEL = Math.sqrt(COEF.length);

let scale = [1.0,1.0,1.0]
let rotation = quat(0.0,3.1,3.1);
let position = [0,0,1.5];

// Begin SH Math
function factorial(x)
{
if(x==0){return 1;}
else{return x*factorial(x-1);}
}

function P(l,m,x)
{
let pmm = 1.0;
if(m>0)
{
let somx2 = Math.sqrt((1-x)*(1+x));
let fact = 1.0;
for(let i=1;i<=m;i++)
{
pmm *= (-fact)*somx2;
fact += 2.0;
}

}
if(l==m) return pmm;
let pmmp1 = x*(2*m+1.0)*pmm;
if(l==m+1){return pmmp1;}
let p11 = 0.0;

for(let i = m+2; i<=l;i++)
{
pll = ((2.0*11-1.0)*x*pmmp1 - (i+m-1)*pmm) / (i-m);
pmm = pmmp1;
pmmp1 = pll;
}
return pll;
}

function K(l,m)
{
let temp = ((2.0*l+1.0)*factorial(l-m))/(4.0*PI*factorial(l+m));
return Math.sqrt(temp);
}

function SH(l,m,theta,phi)
{
// theta in range [0...PI]
// phi in range [0...2PI]
if(m==0){return K(1,0)*P(l,m,Math.cos(theta));}
else if(m>0){return sqrt2*K(l,m)*Math.cos(m*phi)*P(l,m,Math.cos(theta));}
else
{
return sqrt2*K(l,-m)*Math.sin(-m*phi)*P(l,-m,Math.cos(theta));
}
}

function SHBasisIndex(l,m)
{
return l*(l+1) + m;
}

if(RANDOM)
{
for(let i=0;i<COEF.length;i++)
{
COEF[i] = COEF[i]*Math.random()*0.5;
}
}

function GetSHScale(theta,phi)
{
let sum = 0;
for(var l=0;l<LEVEL;l++)
{
for(var m=-l; m<=l;m++)
{
let index = SHBasisIndex(l,m);
sum += COEF[index]*SH(l,m,theta,phi);
}
}

return Math.abs(sum);
}

// End SH Math

// Quaternion from Eular angles
// https://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles
function quat(yaw,pitch,roll)
{
let cy = Math.cos(yaw*0.5);
let sy = Math.sin(yaw*0.5);
let cp = Math.cos(pitch*0.5);
let sp = Math.sin(pitch*0.5);
let cr = Math.cos(roll*0.5);
let sr = Math.sin(roll*0.5);

let w = cy * cp * cr + sy * sp * sr;
let x = cy * cp * sr - sy * sp * cr;
let y = sy * cp * sr + cy * sp * cr;
let z = sy * cp * cr - cy * sp * sr;

return [x,y,z,w];
}

// ref https://github.com/toji/gl-matrix/blob/master/src/mat4.js
// fromRotationTranslationScale
// matrix is stored in the row pattern
function make_transform(o,pos,quat,scale)
{
// scale temp vars
let sx = scale;
let sy = scale;
let sz = scale;

// translation temp vars
let tx = pos;
let ty = pos;
let tz = pos;

//Quat calc temp variables
let x = quat;
let y = quat;
let z = quat;
let w = quat;

let xx = x*x*2.0;
let xy = x*y*2.0;
let xz = x*z*2.0;
let wx = x*w*2.0;

let yy = y*y*2.0;
let yz = y*z*2.0;
let wy = y*w*2.0;

let zz = z*z*2.0;
let wz = z*w*2.0;

o = (1-yy-zz)*sx;
o = (xy-wz)*sy;
o = (xz+wy)*sz;
o = tx;
o = (xy+wz)*sx;
o = (1-xx-zz)*sy;
o = (yz-wx)*sz;
o = ty;
o = (xz-wy)*sx;
o = (yz+wx)*sy;
o = (1-xx-yy)*sz;
o = tz;
o = 0;
o = 0;
o = 0;
o = 1;

return o;
}

function transform_vector(v,m)
{
let x = v, y = v, z = v, w = v;
let o = [0, 0, 0, 0];
o = m * x + m * y + m * z + m * w;
o = m * x + m * y + m * z + m * w;
o = m * x + m * y + m * z + m * w;
o = m * x + m * y + m * z + m * w;
return o;
}

let focus = 200;
function project(v){
let x = v;
let y = v;
let z = v;
let s = focus/z;
return [x*s,y*s];
}

//transfrom matrix
let transform = new Float32Array(16);

let Mesh = function(verts,edges,tsf)
{
let out = {}
out.verts = verts;
out.edges = edges;
out.transform = tsf
return out;
}

// fill in the SH shape vertice and edge list
function make_sh_shape(nu,nv)
{
let o = {};
o.edges = [];
o.verts = [];
let du = 1.0/nu;
let dv = 1.0/nv;
nu +=1;
nv +=1;
let u = 0,v = 0;
let i = 0,j = 0;

if(!WIREFRAME)
{
o.verts.push([0,0,0,1.0]);
}

// push vertice
for(let i=0;i<nu;i++)
{
u = i*du;
for(let j=0;j<nv;j++)
{
v = j*dv;

let cu = Math.cos(u*PI2);
let su = Math.sin(u*PI2);
let cv = Math.cos(v*PI);
let sv = Math.sin(v*PI);

let x = sv*cu;
let y = sv*su;
let z = cv;

let sh_scale = GetSHScale(v*PI,u*PI2);

o.verts.push([x*sh_scale,y*sh_scale,z*sh_scale,1.0]);

//o.verts.push([u,v,1.0,1.0]);
}
}

for(let i=0;i<nu;i++)
{
u = i*du;
for(let j=0;j<nv;j++)
{
v = j*dv;
let id = i*nv+j;
let inext = (i+1)%nu;
let jnext = (j+1)%nv;
let id_down = i*nv+jnext;
let id_right = inext*nv+j;
if(WIREFRAME)
{
o.edges.push(id,id_down);
o.edges.push(id,id_right);
}
else
{
o.edges.push(0,id);
}

}

}
return o;
}

function draw_mesh(tctx,mesh)
{
let edges = mesh.edges;
let verts = mesh.verts;
let tsf = mesh.transform;

for (let i = 0, len = edges.length/2; i < len; i++) {

let i0 = edges[2*i];
let i1 = edges[2*i+1];

// perform transform
let v0 = transform_vector(verts[i0],tsf);
let v1 = transform_vector(verts[i1],tsf);

// projection
v0 = project(v0);
v1 = project(v1);

// draw edge
tctx.penup();
tctx.goto([v0,v0]);
tctx.pendown();
tctx.goto([v1,v1]);

}
}

function draw_mesh_substep(tctx,mesh,step)
{
let edges = mesh.edges;
let len = edges.length/2;
if(step>=len-2)
{
return false;
}

let verts = mesh.verts;
let tsf = mesh.transform;

let i0 = edges[2*step];
let i1 = edges[2*step+1];

// perform transform
let v0 = transform_vector(verts[i0],tsf);
let v1 = transform_vector(verts[i1],tsf);

// projection
v0 = project(v0);
v1 = project(v1);

// draw edge
tctx.penup();
tctx.goto([v0,v0]);
tctx.pendown();
tctx.goto([v1,v1]);
return true;
}

make_transform(transform,position,rotation,scale);

let btl = make_sh_shape(200,100);
let bottle = Mesh(btl.verts,btl.edges,transform);

// The walk function will be called until it returns false.
function walk(i) {

//draw_mesh(turtle,bottle);
let stop = draw_mesh_substep(turtle,bottle,i);
return stop;
}
```