// unit-circle // // Demo of distributing particles on the unit circle // vec2{} initPosns = load("-"); int numOfParticles = length(initPosns); input int iterMax = 500; // maximum number of steps to run input real rad = 0.2; // particle radius real hh0 = 1000; // initial step size, can be big int iter = 1; // which iteration we're on real stepMax = 2*rad; // limit on distance to travel per iter // inter-particle energy, and its derivative function real φ(real r) = (1 - r)^4; function real φ'(real r) = -4*(1 - r)^3; //function real φ(real r) = (1/r)*(1-r)^3; //function real φ'(real r) = 3 - 1/(r^2) - 2*r; function real enr(vec2 x) = φ(|x|/rad); function vec2 frc(vec2 x) = φ'(|x|/rad) * (1/rad) * x/|x|; // chain rule strand Particle (int ii, vec2 pos0) { vec2 pos = pos0/|pos0|; real hh = hh0; output vec2 outPos = pos; vec2 dpos = [0,0]; int id = ii; real energyLast = 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 { // compute energy and forces on us if (!( pos[0] <= 0 || pos[0] > 0 )) { print(ii"HEY pos = ", pos, "\n"); } energyLast = 0; force = [0,0]; int ncount = 0; foreach (Particle p_j in sphere(rad)) { ncount += 1; vec2 r_ij = p_j.pos - pos; energyLast += enr(r_ij); force += frc(r_ij); } vec2 norm = normalize(pos); if (ncount < 2) { if (iter % 3 == 0) { vec2 npos = pos + 0.2*rad*[norm[1],-norm[0]]; new Particle(0, npos); } } // update position based on force force = (identity[2] - norm⊗norm)•force; // project force onto tangent surface dpos = hh*force; if (|dpos| > stepMax) { // decrease hh by factor by which step was too big hh *= stepMax/|dpos|; // and find smaller step dpos = hh*force; // == stepMax*normalize(force); } vec2 posLast = pos; // take step and re-apply constraint pos = normalize(pos + dpos); real energy = 0; foreach (Particle p_j in sphere(rad)) { energy += enr(p_j.pos - pos); } if (energy - energyLast > 0.5*(pos - posLast)•(-force)) { hh *= 0.5; // backtrack pos = posLast; } else { hh *= 1.1; } outPos = pos; } } global { /* real mvmt = mean { |P.dpos|/rad | P in Particle.all}; if (mvmt < 0.0001) { print("stabilized at iter ", iter, "\n"); stabilize; } */ iter+=1; if (iter > iterMax) { print("numActive() = ", numActive(), "\n"); stabilize; // HEY really want "die" here: didn't converge } // real var = variance{P.energy | P in Particle.all}; // real var = sum{(avg-P.energy)^2 | P in Particle.all}; } initially {Particle(ii, initPosns[ii]) | ii in 1 .. numOfParticles };
Click to toggle
does not end with </html> tag
does not end with </body> tag
The output has ended thus: n Particle.all}; } initially {Particle(ii, initPosns[ii]) | ii in 1 .. numOfParticles };