GENERIC MODULE ElectricEvaluator (LocalAttributes, LocalOps, LocalEvaluator);
  (*
    WHERE
      /\ IN LocalAttributes
          /\ TYPE T
      /\ INTERFACE LocalOps = LocalOps(LocalAttributes)
      /\ INTERFACE LocalEvaluator = LocalEvaluator(LocalAttributes)  
  *)

REVEAL
  T = LocalEvaluator.T BRANDED OBJECT
      numCharges: CARDINAL;
      totCharge: LONGREAL;
      q: REF ARRAY OF LocalAttributes.T;
      w: REF ARRAY OF REAL;
      radius: LONGREAL;
    OVERRIDES
      eval := TEval;
      print := TPrint;
    END;

PROCEDURE New(
    READONLY q: ARRAY OF LocalAttributes.T; (* Charge positions *)
    READONLY w: ARRAY OF REAL;              (* Charge magnitudes *)
    radius: LONGREAL;
  ): T =
  VAR wt := NEW(REF ARRAY OF REAL, NUMBER(q));
      qt := NEW(REF ARRAY OF LocalAttributes.T, NUMBER(q));
      n: CARDINAL := 0;
      totCharge: LONGREAL := 0.0d0;
      s: LONGREAL;
      dsdqn: LocalAttributes.T;
  BEGIN
    WITH
      w = wt^,
      q = qt^,
      tolSqr = radius*radius
    DO
      <* ASSERT NUMBER(q) = NUMBER(w) *>
      <* ASSERT radius > 0.0d0 *>
      FOR i := 1 TO LAST(q) DO
        WITH
          wn = w[n],
          qn = q[n]
        DO
          totCharge := totCharge + FLOAT(wn, LONGREAL);
          IF n < i THEN wn := w[i]; qn := q[i] END;
          INC(n);
          FOR j := 0 TO n-2 DO
            WITH wj = w[j], qj = q[j] DO
              LocalAttributes.DistSqr(qn, qj, s, dsdqn);
              IF s < tolSqr THEN
                CombineCharges(wj, qj, wn, qn);
                DEC(n); EXIT
              END
            END
          END
        END
      END;
      IF n < NUMBER(qt^) THEN
        WITH 
          wtNew = NEW(REF ARRAY OF REAL, n), 
          qtNew = NEW(REF ARRAY OF LocalAttributes.T, n)
        DO
          wtNew^ := SUBARRAY(wt^, 0, n);
          qtNew^ := SUBARRAY(qt^, 0, n);
          wt := wtNew; qt := qtNew
        END
      END;
      RETURN NEW (T,
        numCharges := n,
        totCharge := totCharge,
        radius := radius,
        w := wt,
        q := qt
      }
    END
  END New;
  
PROCEDURE CombineCharges(
    VAR wa: REAL;
    VAR pa: LocalAttributes.T;
    READONLY wb: REAL;
    READONLY pb: LocalAttributes.T;
  ) =
  (*
    Combines "wa,pa" and "wb,pb" into a single charge, stored into 
    the former's location.
    Should alco compute a new radius, or better and ellipsoid... *)
  BEGIN
    WITH 
      wwa = FLOAT(wa, LONGREAL),
      wwb = FLOAT(wb, LONGREAL),
      wwt = wwa+wwb
      sa = wwa/wwt, sb = wwb/wwt 
    DO
      pa := LocalAttributes.Mix(sa, pa, sb, pb);
      wa := FLOAT(wwt)
    END
  END CombineCharges;
  
PROCEDURE TEval(
    self: T;
    READONLY p: LocalAttributes.T;       (* The local attributes of a texel *)
    VAR (*OUT*) E: LONGREAL;             (* Energy term for this texel *)
    VAR (*OUT*) dEdp: LocalAttributes.T; (* Partial derivatives "dE/dp" *)
  ) =
  BEGIN
    E := 0.0d0;
    dEdp := LocalAttributes.Zeros;
    WITH
      q = self.q^,
      w = self.w^,
      radius = self.radius
    DO
      FOR j := 0 TO a.nCharges-1 DO
        WITH qj = q[j], wj = FLOAT(w[j], LONGREAL) DO
          AddSingleChargeEnergy(p, qj, wj, radius, E, dEdp);
        END
      END
    END;
  END TEval;

PROCEDURE AddSingleChargeEnergy(
    READONLY p: LocalAttributes.T;      (* Position of test charge *)
    READONLY q: LocalAttributes.T;      (* Position of fixed diffuse charge *)
    wq: LONGREAL;                       (* Strength of "q" *)
    rq: LONGREAL;                       (* Nominal radius of "q" *)
    VAR (*I/O*) E: LONGREAL;            (* Accumulated energy *)
    VAR (*I/O*) dEdp: LocalAttributes.T;(* Accumulated partials "dE/dp" *)
  ) = 
  (*
    Adds to "E" and "dEdp" the electric potential of a point-like negative 
    test charge at point "p" due to a fuzzy charge of magnitude "w" AND 
    extend "radius" fixed at "q".
  *)
  VAR s: LONGREAL;               (* Distance squared between "p" and "q" *)
      dsdp: LocalAttributes.T;   (* Partials "ds/dp[*]" *)
  BEGIN
    (* Compute ss, dsdp: *)
    s := rq*rq; 
    dsdp := LocalAttributes.Zeros;
    LocalAttributes.DistSqr(p, q, s, dsdp)
    
    (* Now compute energy and derivatives: *)
    WITH
      r = Math.sqrt(s),
      e = - (wq / r - wq / rq),
      deds = 0.5d0 * e / s
    DO
      E := E + e;
      EVAL LocalAttributes.Modify(dEdp, deds, dsdp) 
    END
  END AddSingleChargeEnergy;

PROCEDURE TPrint(self: T; wr: Wr.T) =

  PROCEDURE FR(x: REAL): TEXT =
    BEGIN
      RETURN Fmt.Real(x, Fmt.Style.Sci, 4)
    END FR;

  <* FATAL Wr.Failure, Thread.Alerted *>
  VAR wr: Wr.T;
  BEGIN
    WITH
      fName = name & ".plot",
      radius = FLOAT(self.radius),
      q = self.pt^,
      w = self.wt^
    DO
      Wr.PutText(wr, "# num charges = " & Fmt.Int(self.numCharges) & "\n");
      Wr.PutText(wr, "# tot charge = " & Fmt.Real(self.totCharge) & "\n");
      Wr.PutText(wr, "# num attributes = " & Fmt.Int(LocalAttributes.N) & "\n");
      Wr.PutText(wr, "# attribute description = " & LocalAttributes.Description & "\n");
      FOR j := 0 TO self.numCharges-1 DO
        WITH qr = q[j], wj = w[j] DO 
          Wr.PutText(wr, FR(radius)); Wr.PutChar(wr, ' ');
          Wr.PutText(wr, FR(wj));  Wr.PutChar(wr, ' ');
          Wr.PutChar(wr, ' ');
          LocalAttributes.Print(wr, qj);
          Wr.PutChar(wr, '\n');
        END
      END;
      Wr.Flush(wr);
    END
  END TPrint;


BEGIN
END ElectricEvaluator.