// unit-circle // // Demo of distributing particles on the unit circle // real{} initPosns = load("-"); 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 // ============================ begin setting up gridding of domain vec2 xDom = [-1,1]; vec2 yDom = [-1,1]; real xSamples = floor((xDom - xDom)/2); // Lamont verify logic floor vs ceil here real ySamples = floor((yDom - yDom)/2); vec4 qWinDim = [xDom,xDom,yDom,yDom]; // i.e. [XMIN, XMAX, YMIN, YMAX] (required for the query function) vec2 qGridDim = [xSamples,ySamples]; // how many grid cells you want in each direction for the uniform grid (required for the query function) vec2 qCellDim = [(xDom - xDom)/xSamples,(yDom - yDom)/ySamples]; // the length in each direction for a cell (required for the query function) // ============================ end setting up gridding of domain 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 - 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.001, pos + 0.01); } outPos = pos; } } // can add statements in here to do some global computation global{ iter+=1; real var = variance{P.energy | P in Particle.all}; stdev = sqrt(var); } initially {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: tially {Particle(ii, initPosns{ii*2}, initPosns{ii*2+1}) | ii in 0 .. numOfParticles-1 };