Home My Page Projects Code Snippets Project Openings diderot

# SCM Repository

[diderot] View of /tests/lamont-tests/unit-circle.diderot
 [diderot] / tests / lamont-tests / unit-circle.diderot

# View of /tests/lamont-tests/unit-circle.diderot

Tue Sep 27 20:54:47 2016 UTC (2 years, 11 months ago) by glk
File size: 3429 byte(s)
`initial result of svn export --username anonsvn --password=anonsvn https://svn.smlnj-gforge.cs.uchicago.edu/svn/diderot/branches/vis15/src/tests/`
```#version 2

// unit-circle
//
// Demo of distributing particles on the unit circle
//
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 }

```