
<<-- COGParticles.mesa: particles by inverse square law>>
<<-- last modified by Stolfi - October 12, 1982 12:30 am>>
-- To Do: better integration formula!!!

-- To Run:   run COGAll; run COGParticles


DIRECTORY
  IO USING [PutF, GetChar, GetInt, GetReal],
  COGDebug USING [in, out],
  Rope USING [ROPE],
  Real USING [SqRt],
  GraphicsColor USING [red, black],
  Process USING [SetPriority, priorityBackground],
  COGCart USING [Point],
  COGSpace USING [Point, Vector, Box, Random, Add, Sub, Mul, LengthSq],
  COGDrawing;
COGParticles: CEDAR PROGRAM
  IMPORTS
    IO, COGDrawing, COGSpace, COGDebug, Real, Process =
BEGIN
  OPEN
    Draw: COGDrawing, Cart: COGCart, Sp: COGSpace, Bug: COGDebug,
    Real, Rope, GraphicsColor, IO;
dr: Draw.Drawing = Draw.MakeDrawing["Particles", [min: [-1.0, -1.0], max: [1.0, 1.0]]];
posBox: Sp.Box = [[-0.5, -0.5, -0.5], [0.5, 0.5, 0.5]];
velBox: Sp.Box = [[-0.1, -0.1, -0.1], [0.1, 0.1, 0.1]];
Particle: TYPE = RECORD [pos, vel, acc: Sp.Vector, mass, charge: REAL];
System: TYPE = RECORD
  [K: REAL, -- grav. constant
   dt: REAL, -- time step for display
   mang: REAL, -- maximum rel. angular displacement during integration step 
   quit: BOOL, -- set to FALSE by Start button, to TRUE by Stop button
   part: SEQUENCE np: NAT OF Particle];

ComputeAcc: PROC [sys: REF System] RETURNS [maxav: REAL] = TRUSTED
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; 
TimeStep: Draw.ObjectAction = TRUSTED
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; 

SysPainter: Draw.PainterProc = TRUSTED
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; 

curSys: REF System;
curSysOb: Draw.Object;

Start: Draw.MenuProc = TRUSTED
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;
Stop: Draw.MenuProc =
BEGIN
  -- called with access = none for immediate effect
  curSys.quit _ TRUE
END;

NewSys: Draw.MenuProc =
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;
ThrowSys: PROC [sys: REF System] = TRUSTED
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;

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

Bug.out.PutF["\nHello! say something: "]; [] _ Bug.in.GetChar[];
Draw.AddMenuAction [dr, "Start", Start, NIL, none];
Draw.AddMenuAction [dr, "Stop", Stop, NIL, none];
Draw.AddMenuAction [dr, "NewSys", NewSys, NIL, none];
Draw.AddMenuAction [dr, "Reset", Reset, NIL, write];
Bug.out.PutF["\nEnjoy -- bye bye!"];
END...