#version 2 // unit-circle // // Demo of distributing particles on the unit circle // //real[] initPosns = load_sequence("-"); real[] initPosns = load_sequence("100.nrrd"); int numOfParticles = length(initPosns)/2; input int iterMax = 300; // maximum number of steps to run input real rr = 0.5; // actual particle radius real stdev = ∞; real threshold = 0.08; // covnergence threshold real RR = rr+0.0; // neighbor query radius (MUST be >= rr; should get same results for any RR >= rr) real hhInit = 10.0; // initial integration step size; can err too big, will be trimmed down during iterations int iter = 1; // which iteration we're on real stepMax = 2*π/numOfParticles; // limit on distance to travel per iter strand Particle (int ii, real posx, real posy) { vec2 pos = normalize([posx,posy]); real hh = hhInit; vec2 posOld1 = pos; // remember last TWO positions vec2 posOld2 = pos; output vec2 outPos = pos; int id = ii; real energy = ii; // HEY: can do convergence test based on variance of energies vec2 force = [0,0]; // or can test convergence based on sum of |force| update { if( iter > iterMax) { stabilize; } // compute energy and forces on us energy = 0; force = [0,0]; foreach (Particle p_j in sphere(RR)) { vec2 r_ij = (pos - p_j.pos)/rr; vec2 d_ij = normalize(r_ij); if (|r_ij| < 1) { energy += (1 - |r_ij|)^4; force += - (-4*(1 - |r_ij|)^3) * d_ij; } } force /= rr; // smaller particles make larger forces // update position based on force posOld2 = posOld1; // shuffle saved positions down posOld1 = pos; if (energy > 0.0) { // we have neighbors tensor[2,2] pten = identity[2] - normalize(pos)⊗normalize(pos); force = pten•force; // project force onto tangent surface vec2 step = hh*force; if (|step| > stepMax) { // decrease hh by factor by which step was too big hh *= stepMax/|step|; // and find smaller step step = hh*force; } // take step and re-find implicit surface pos = normalize(pos + step); real travel = |pos - posOld1| + |posOld1 - posOld2|; if (travel > 0) { // if we've moved in the past two steps, but we've moved back // to where we were two steps ago, we're oscillating ==> // okay = 0. Two steps in the same direction ==> okay = 1. real okay = |pos - posOld2|/travel; // slow down if oscillating, speed up a little if not hh *= lerp(0.8, 1.001, 0, okay, 1); } } if (id == 0 && iter == 140) { new Particle(2,pos[0] + 0.001, pos[1] + 0.01); } outPos = pos; } } // can add statements in here to do some global computation update { // real var = variance{P.energy | P in Particle.all}; real avg = mean{P.energy | P in Particle.all}; print("iter ", iter, ": avg ", avg, "\n"); iter+=1; real var = sum{(avg-P.energy)^2 | P in Particle.all}; stdev = sqrt(var); } create_collection {Particle(ii, initPosns[ii*2], initPosns[ii*2+1]) | ii in 0 .. numOfParticles-1 }
Click to toggle
does not end with </html> tag
does not end with </body> tag
The output has ended thus: lection {Particle(ii, initPosns[ii*2], initPosns[ii*2+1]) | ii in 0 .. numOfParticles-1 }