-- To Do: better integration formula!!!

BEGIN
-- computes particle accelerations and returns also the max approx relative ang. velocity (rad)
part: Particle;
rpq, ap: Sp.Vector;
spq, cp, cpq, vpq: REAL;
maxav ← 0.0;
FOR i: NAT IN [0..sys.np) DO
part ← sys.part[i];
ap ← [0.0, 0.0, 0.0];
FOR j: NAT IN [0..sys.np) DO
IF i#j THEN
{rpq ← Sp.Sub [part.pos, sys.part[j].pos]; -- relative displacement
spq ← Sp.LengthSq [rpq];
cpq ← sys.part[j].charge/(spq*Real.SqRt[spq]);
ap ← Sp.Sub [ap, Sp.Mul [cpq, rpq]];
IF i < j THEN
{vpq ← Sp.LengthSq [Sp.Sub [part.vel, sys.part[j].vel]]; -- relative velocity
maxav ← MAX [maxav, vpq/spq]}
}
ENDLOOP;
cp ← sys.K * sys.part[i].charge / sys.part[i].mass;
sys.part[i].acc ← Sp.Mul [cp, ap];
ENDLOOP;
maxav ← Real.SqRt[maxav]
END;

BEGIN
sys: REF System = NARROW [obj.parms.data];
part: Particle;
pos: Sp.Vector;
ddt, tt, maxav: REAL;
tt ← sys.dt;
WHILE tt > 0.0 DO
maxav ← ComputeAcc[sys];
ddt ← MIN[tt, sys.mang/maxav];
FOR i: NAT IN [0..sys.np) DO
part ← sys.part[i];
pos ← Sp.Add [part.pos, Sp.Mul [ddt, part.vel]];
sys.part[i].pos ← Sp.Add [pos, Sp.Mul [0.5*ddt*ddt, part.acc]];
sys.part[i].vel ← Sp.Add [part.vel, Sp.Mul [ddt, part.acc]]
ENDLOOP;
tt ← tt-ddt;
ENDLOOP;
RETURN [NIL]
END;

BEGIN
sys: REF System = NARROW [parms.data];
xp, yp, zp: REAL;
FOR i: NAT IN [0..sys.np) DO
[xp, yp, zp] ← sys.part[i].pos;
Draw.DrawDot [context, sf, [xp, yp, 1.0], 2, IF sys.part[i].charge < 0.0 THEN black ELSE red]
ENDLOOP
END;

BEGIN
-- called with access = none
sys: REF System ← curSys;
sysOb: Draw.Object ← curSysOb;
Process.SetPriority[Process.priorityBackground];
sys.quit ← FALSE;
WHILE NOT sys.quit AND Draw.Alive[dr] DO
[] ← Draw.Modify [obj: sysOb, Action: TimeStep, erase: FALSE, repaint: TRUE]
ENDLOOP
END;

BEGIN
-- called with access = none for immediate effect
curSys.quit ← TRUE
END;

BEGIN
-- called with access = none
N: NAT;
sys: REF System;
c: Draw.Change;
Bug.out.PutF ["\n N. Particles = "]; N ← Bug.in.GetInt[];
sys ← NEW [System [N]];
Bug.out.PutF ["\n dt = "]; sys.dt ← Bug.in.GetReal [];
Bug.out.PutF ["\n max ang vel = "]; sys.mang ← Bug.in.GetReal [];
Bug.out.PutF ["\n grav.constant = "]; sys.K ← Bug.in.GetReal [];
ThrowSys [sys];
c ← Draw.GetWriteRights [dr];
BEGIN ENABLE UNWIND => {Draw.Release[dr, c]};
IF curSys # NIL THEN
{curSys.quit ← TRUE;
Draw.Remove [curSysOb]};
curSysOb ← Draw.Make[SysPainter, [data: sys]];
curSys ← sys;
Draw.Add [dr, curSysOb, 3];
Draw.DoPaintAll [dr];
Draw.Release [dr, c]
END
END;

BEGIN
xs, ys, zs, m, mt: REAL ← 0.0;
slow, dij: REAL;
ch: REAL ← 1.0;
pos, vel, dst: Sp.Vector;
Ep, Ec: REAL ← 0.0;
FOR i: NAT IN [0..sys.np) DO
sys.part[i].mass ← 1.0; sys.part[i].charge ← ch; ch ← -ch
ENDLOOP;
sys.part[0].mass ← sys.np-1;
sys.part[0].pos ← [0.0, 0.0, 0.0];
sys.part[0].vel ← [0.0, 1.0, 0.0];
sys.part[1].pos ← [0.5, 0.0, 0.0];
sys.part[1].vel ← [0.0, -0.7, 0.0];
FOR i: NAT IN [2..sys.np) DO
sys.part[i].pos ← Sp.Random [posBox];
sys.part[i].vel ← Sp.Random [velBox];
ENDLOOP;
FOR i: NAT IN [0..sys.np) DO
pos ← sys.part[i].pos;
vel ← sys.part[i].vel;
m ← sys.part[i].mass; ch ← sys.part[i].charge;
Ec ← Ec + m*(vel.x*vel.x + vel.y*vel.y + vel.z*vel.z);
FOR j: NAT IN (i..sys.np) DO
dst ← Sp.Sub [sys.part[j].pos, pos];
dij ← SqRt[dst.x*dst.x+dst.y*dst.y+dst.z*dst.z];
Ep ← Ep + ch*sys.part[j].charge/dij
ENDLOOP;
mt ← mt + m;
xs ← xs + m*vel.x; ys ← ys + m*vel.y; zs ← zs+ m*vel.z;
ENDLOOP;
xs ← xs/mt; ys ← ys/mt; zs ← zs/mt;
Ec ← Ec - mt*(xs*xs + ys*ys + zs*zs); Ec ← Ec/2;
Ep ← - Ep*sys.K;
slow ← IF Ep < 0 THEN SqRt[-Ep/(2.0*Ec)] ELSE 1.0;
FOR i: NAT IN [0..sys.np) DO
vel ← sys.part[i].vel;
sys.part[i].vel ← [slow*(vel.x-xs), slow*(vel.y-ys), slow*(vel.z-zs)]
ENDLOOP
END;

BEGIN
-- called with write access
ThrowSys [curSys];
Draw.DoPaintAll [dr]
END;