Skip to content

Instantly share code, notes, and snippets.

@Gro-Tsen

Gro-Tsen/index.html

Last active April 30, 2020 14:38
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save Gro-Tsen/3d1c1fb2bbe3e4166e6beae2a9a957f1 to your computer and use it in GitHub Desktop.
Save Gro-Tsen/3d1c1fb2bbe3e4166e6beae2a9a957f1 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
index.html
spherical-interactions.html
Raw
spherical-interactions.html
<!DOCTYPE html>
<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
<head>
<meta charset="utf-8" />
<title>Spherical Interactions</title>
</head>
<body>
<h1>Spherical Interactions</h1>
<canvas id="canvas" width="1000" height="750">
</canvas>
<p>Energy: <span id="energy">?</span></p>
<p>Angular momentum: (<span id="angmom-0">?</span>, <span id="angmom-1">?</span>, <span id="angmom-2">?</span>)</p>
<script type="text/javascript">
// <![CDATA[
"use strict";
function gaussian() {
// Generate a Gaussian random variable by Box-Muller.
"use strict";
var u0 = Math.random();
var u1 = Math.random();
return Math.sqrt(-2*Math.log(u0))*Math.cos(2*Math.PI*u1);
}
function dotprod(vec1, vec2) {
// Compute and return dot product of two vectors.
"use strict";
var v = 0;
for ( var i=0 ; i<vec1.length ; i++ )
v += vec1[i]*vec2[i];
return v;
}
function renormalize(vec) {
// Alter vector to make it of unit length.
"use strict";
var norm = Math.sqrt(dotprod(vec, vec));
for ( var i=0 ; i<vec.length ; i++ )
vec[i] /= norm;
}
var nbpoints = 12; // Number of points to simulate
var colortab = [ // Color of the points (should have nbpoints entries)
"rgb(255,0,0)", "rgb(255,128,0)", "rgb(255,255,0)", "rgb(128,255,0)",
"rgb(0,255,0)", "rgb(0,255,128)", "rgb(0,255,255)", "rgb(0,128,255)",
"rgb(0,0,255)", "rgb(128,0,255)", "rgb(255,0,255)", "rgb(255,0,128)"
];
// Force constant:
var coulomb_constant = -0.2; // Positive is attractive, negative repulsive.
function points_acc_func(pos, vel) {
// Compute acceleration of points in function of position and velocity.
"use strict";
var acc = new Array(nbpoints);
for ( var i=0 ; i<nbpoints ; i++ ) {
acc[i] = [0, 0, 0];
for ( var j=0 ; j<nbpoints ; j++ )
if ( j != i )
for ( var k=0 ; k<3 ; k++ ) {
var cosrho = dotprod(pos[j],pos[i]);
/// Euclidean elastic force:
// acc[i][k] += coulomb_constant*(pos[j][k] - pos[i][k]);
/// Test force:
// acc[i][k] += coulomb_constant*(pos[j][k] - pos[i][k]) * cosrho;
/// Combined force of previous two:
// acc[i][k] += coulomb_constant*(pos[j][k] - pos[i][k]) * (1+cosrho)*0.5;
/// Spherical Coulomb force (V = −k/tan(ρ)):
// acc[i][k] += coulomb_constant*(pos[j][k] - pos[i][k]) / Math.pow(1-cosrho*cosrho, 1.5);
/// Whatever this is:
// acc[i][k] += coulomb_constant*(pos[j][k] - pos[i][k]) * cosrho / Math.pow(1-cosrho*cosrho, 1.5);
/// Euclidean Coulomb force:
acc[i][k] += coulomb_constant*(pos[j][k] - pos[i][k]) / Math.pow(2*(1-cosrho), 1.5);
}
/// Damping (this breaks conservation of energy (duh)):
// for ( var k=0 ; k<3 ; k++ )
// acc[i][k] -= vel[i][k]*0.01;
/// Make points stay on the sphere (essentially D'Alembert's principle):
var ortho = dotprod(pos[i], acc[i]) + dotprod(vel[i], vel[i]);
ortho /= dotprod(pos[i], pos[i]);
for ( var k=0 ; k<3 ; k++ )
acc[i][k] -= ortho*pos[i][k];
}
return acc;
}
// State variables:
var time = 0; // Absolute time
var points_pos; // Euclidean position of points (vecs of length 3)
var points_vel; // Euclidean velocities of points (vecs of length 3)
var keep_past_length = 10;
var past_pos; // Past velocities
function angular_momentum() {
// Compute angular momentum (a conserved quantity).
"use strict";
var angmom = [ 0,0,0 ];
for ( var i=0 ; i<nbpoints ; i++ ) {
angmom[0] += points_pos[i][1] * points_vel[i][2] - points_pos[i][2] * points_vel[i][1];
angmom[1] += points_pos[i][2] * points_vel[i][0] - points_pos[i][0] * points_vel[i][2];
angmom[2] += points_pos[i][0] * points_vel[i][1] - points_pos[i][1] * points_vel[i][0];
}
return angmom;
}
function energy() {
// Compute total energy (a conserved quantity).
"use strict";
var energy = 0;
for ( var i=0 ; i<nbpoints ; i++ ) {
// Kinetic part:
energy += dotprod(points_vel[i], points_vel[i])*0.5;
// Potential part:
for ( var j=0 ; j<i ; j++ ) {
var cosrho = dotprod(points_pos[j],points_pos[i]);
/// Euclidean elastic force:
// energy -= coulomb_constant * cosrho;
/// Test force:
// energy -= coulomb_constant * 0.5*cosrho*cosrho;
/// Combined force of previous two:
// energy -= coulomb_constant * cosrho * (1+0.5*cosrho)*0.5;
/// Spherical Coulomb force (V = −k/tan(ρ)):
// energy -= coulomb_constant * cosrho/Math.sqrt(1-cosrho*cosrho);
/// Whatever this is:
// energy -= coulomb_constant / Math.sqrt(1-cosrho*cosrho);
/// Euclidean Coulomb force:
energy -= coulomb_constant / Math.sqrt(2*(1-cosrho));
}
}
return energy;
}
function init_points() {
// Start with uniformly distributed random points, zero velocities.
"use strict";
points_pos = new Array(nbpoints);
points_vel = new Array(nbpoints);
for ( var i=0 ; i<nbpoints ; i++ ) {
points_pos[i] = [gaussian(), gaussian(), gaussian()];
renormalize(points_pos[i])
points_vel[i] = [0,0,0];
// points_vel[i] = [gaussian()*0.2, gaussian()*0.2, gaussian()*0.2];
var ortho = dotprod(points_pos[i], points_vel[i]);
for ( var k=0 ; k<3 ; k++ )
points_vel[i][k] -= ortho*points_pos[i][k];
}
past_pos = new Array();
past_pos.push(points_pos);
}
init_points();
function runge_kutta_evolve(step) {
// Fourth-order Runge-Kutta
"use strict";
var points_pos_1 = points_pos;
var points_vel_1 = points_vel;
var points_acc_1 = points_acc_func(points_pos_1, points_vel_1);
var points_pos_2 = new Array(nbpoints);
var points_vel_2 = new Array(nbpoints);
for ( var i=0 ; i<nbpoints ; i++ ) {
points_pos_2[i] = new Array(3);
points_vel_2[i] = new Array(3);
for ( var k=0 ; k<3 ; k++ ) {
points_pos_2[i][k] = points_pos[i][k] + points_vel_1[i][k]*step/2;
points_vel_2[i][k] = points_vel[i][k] + points_acc_1[i][k]*step/2;
}
}
var points_acc_2 = points_acc_func(points_pos_2, points_vel_2);
var points_pos_3 = new Array(nbpoints);
var points_vel_3 = new Array(nbpoints);
for ( var i=0 ; i<nbpoints ; i++ ) {
points_pos_3[i] = new Array(3);
points_vel_3[i] = new Array(3);
for ( var k=0 ; k<3 ; k++ ) {
points_pos_3[i][k] = points_pos[i][k] + points_vel_2[i][k]*step/2;
points_vel_3[i][k] = points_vel[i][k] + points_acc_2[i][k]*step/2;
}
}
var points_acc_3 = points_acc_func(points_pos_3, points_vel_3);
var points_pos_4 = new Array(nbpoints);
var points_vel_4 = new Array(nbpoints);
for ( var i=0 ; i<nbpoints ; i++ ) {
points_pos_4[i] = new Array(3);
points_vel_4[i] = new Array(3);
for ( var k=0 ; k<3 ; k++ ) {
points_pos_4[i][k] = points_pos[i][k] + points_vel_3[i][k]*step;
points_vel_4[i][k] = points_vel[i][k] + points_acc_3[i][k]*step;
}
}
var points_acc_4 = points_acc_func(points_pos_4, points_vel_4);
var new_points_pos = new Array(nbpoints);
var new_points_vel = new Array(nbpoints);
for ( var i=0 ; i<nbpoints ; i++ ) {
new_points_pos[i] = new Array(3);
new_points_vel[i] = new Array(3);
for ( var k=0 ; k<3 ; k++ ) {
new_points_pos[i][k] = points_pos[i][k] + (points_vel_1[i][k]+points_vel_2[i][k]*2+points_vel_3[i][k]*2+points_vel_4[i][k])*step/6;
new_points_vel[i][k] = points_vel[i][k] + (points_acc_1[i][k]+points_acc_2[i][k]*2+points_acc_3[i][k]*2+points_acc_4[i][k])*step/6;
}
}
// Make sure points are still on the sphere and moving tangentially to it:
for ( var i=0 ; i<nbpoints ; i++ ) {
renormalize(new_points_pos[i]);
var ortho = dotprod(new_points_pos[i], new_points_vel[i]);
for ( var k=0 ; k<3 ; k++ )
new_points_vel[i][k] -= ortho*new_points_pos[i][k];
}
time += step;
points_pos = new_points_pos;
points_vel = new_points_vel;
}
var canvas;
var ctx;
canvas = document.getElementById("canvas");
if ( typeof(canvas.getContext) != "function" ) {
alert("Your browser does not support the HTML5 <canvas> element.\n"
+ "This page will not function.");
throw new Error("Canvas unsupported");
}
ctx = canvas.getContext("2d");
ctx.lineCap = "round";
ctx.lineJoin = "round";
var facets = [
// The six views of the sphere, each with its orthonormal frame:
{ cen: [125,375], frame: [[0,0,1], [0,1,0], [-1,0,0]] },
{ cen: [375,125], frame: [[1,0,0], [0,0,-1], [0,1,0]] },
{ cen: [375,375], frame: [[1,0,0], [0,1,0], [0,0,1]] },
{ cen: [375,625], frame: [[1,0,0], [0,0,1], [0,-1,0]] },
{ cen: [625,375], frame: [[0,0,-1], [0,1,0], [1,0,0]] },
{ cen: [875,375], frame: [[-1,0,0], [0,1,0], [0,0,-1]] },
];
function draw_points() {
// Draw everything.
"use strict";
ctx.clearRect(0, 0, 1000, 750);
for ( var f=0 ; f<facets.length ; f++ ) {
ctx.fillStyle = "rgb(128,128,128)";
// Draw grey cricles:
ctx.beginPath();
ctx.arc(facets[f].cen[0], facets[f].cen[1], 120, 0, 2*Math.PI, false);
ctx.fill();
// Draw dotted "tail":
for ( var t=0 ; t<past_pos.length ; t++ ) {
for ( var i=0 ; i<nbpoints ; i++ ) {
var prj = new Array(3);
for ( var k=0 ; k<3 ; k++ )
prj[k] = dotprod(past_pos[t][i], facets[f].frame[k]);
if ( prj[2] >= 0 ) {
ctx.fillStyle = colortab[i];
ctx.fillRect(facets[f].cen[0] + prj[0]*120, facets[f].cen[1] - prj[1]*120, 1, 1);
}
}
}
// Draw the points themselves:
for ( var i=0 ; i<nbpoints ; i++ ) {
var prj = new Array(3);
for ( var k=0 ; k<3 ; k++ )
prj[k] = dotprod(points_pos[i], facets[f].frame[k]);
if ( prj[2] >= 0 ) {
var x = prj[0];
var y = prj[1];
var z = prj[2];
// Draw each point as a little ellipse (projection of
// a small tangent circle):
ctx.save();
ctx.translate(facets[f].cen[0] + x*120, facets[f].cen[1] - y*120);
ctx.rotate(Math.atan2(-y, x));
ctx.scale(z, 1);
ctx.fillStyle = colortab[i];
ctx.beginPath();
ctx.arc(x, y, 3, 0, 2*Math.PI, false);
ctx.fill();
ctx.restore();
}
}
}
// Compute and display energy and angular momentum:
var engy = energy();
var elt = document.getElementById("energy");
while ( elt.lastChild )
elt.removeChild(elt.lastChild);
elt.appendChild(document.createTextNode(engy.toFixed(3)));
var angmom = angular_momentum();
for ( var k=0 ; k<3 ; k++ ) {
var elt = document.getElementById("angmom-"+k);
while ( elt.lastChild )
elt.removeChild(elt.lastChild);
elt.appendChild(document.createTextNode(angmom[k].toFixed(3)));
}
}
draw_points();
function update() {
"use strict";
var newtime = time + 0.05;
var prevstep = 0.05;
while ( time < newtime ) {
// Try to find appropriate step size (note: this uses
// conservation of energy, so simplify or rewrite all this if
// damping is in place):
var stepsize = newtime - time;
if ( stepsize > 2*prevstep )
stepsize = 2*prevstep;
var pre_time = time;
var pre_points_pos = points_pos;
var pre_points_vel = points_vel;
var pre_engy = energy();
while ( 1 ) {
runge_kutta_evolve(stepsize);
var post_engy = energy();
// console.log("time="+time+", stepsize="+stepsize+", relerror="+Math.abs(post_engy-pre_engy)/stepsize);
if ( Math.abs(post_engy-pre_engy)/stepsize < 1.e-3
|| stepsize < 0.001 )
break; // Don't decrease step size further
// console.log("halving stepsize");
stepsize /= 2; // Try smaller step size
// Yuck :-\
time = pre_time;
points_pos = pre_points_pos;
points_vel = pre_points_vel;
}
prevstep = stepsize;
}
// Record previous position:
past_pos.push(points_pos);
if ( past_pos.length >= keep_past_length )
past_pos.shift();
// Draw, rinse, repeat:
draw_points();
window.setTimeout(update, 50);
}
window.setTimeout(update, 50);
// ]]>
</script>
</body>
</html>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment