// 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.2; // actual particle radius 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[1] - xDom[0])/RR); // Lamont verify logic floor vs ceil here real ySamples = floor((yDom[1] - yDom[0])/RR); vec4 qWinDim = [xDom[0],xDom[1],yDom[0],yDom[1]]; // 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[1] - xDom[0])/xSamples,(yDom[1] - yDom[0])/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; real energy = ii; // should eventually do convergence test based on variance of energies update { // print positions print(atan2(pos[1],pos[0]), " "); if (ii == numOfParticles-1) { print("\n"); } // compute energy and forces on us energy = 0; vec2 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); // oscillating if moved close to where we were two steps ago real osci = |pos - posOld2|/|pos - posOld1|; if (osci < 0.1) { // altering step size slightly disrupts oscillation hh *= 0.9; } } outPos = pos; if (iter >= iterMax) { stabilize; } } } // can add statements in here to do some global computation global{ iter+=1; } 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 };