PROCEDURE AccumulateTerm(
          READONLY u, v, p, q: LR3.T;
          VAR eDu, eDv, eDp, eDq: LR3.T;
        ) =
        (*
          Adds to "e" the energy for one hinge with endpoints
          "u" and "v", and wingtips "p" and "q".  
          Also adds to "eDu", "eDv", "eDp", "eDq" the 
          corresponding gradients.
        *)
        BEGIN


          IF NOT LR3.Equal(u,v) THEN
            WITH
              r = LR3.Sub(v, u),
              a = LR3.Sub(q, u), 
              b = LR3.Sub(v, p)
            DO
              VAR eTerm: LONGREAL;
                  eDa, eDb, eDr: LR3.T;
              BEGIN
                Compute_e_from_a_b_r(a, b, r, eTerm, eDa, eDb, eDr);
                e := e + eTerm;
                eDp := LR3.Sub(eDp, eDb);
                eDq := LR3.Add(eDq, eDa);
                eDu := LR3.Sub(eDu, LR3.Add(eDr, eDa));
                eDv := LR3.Add(eDv, LR3.Add(eDb, eDr))
              END
            END
          END
        END AccumulateTerm;
              
      PROCEDURE Compute_e_from_a_b_r(
          READONLY a, b, r: LR3.T;
          VAR eTerm: LONGREAL;
          VAR eDa, eDb, eDr: LR3.T;
        ) =
        (*
          Computes the hinge energy "eTerm" given the spine vector "r"
          and two vectors "a" and "b" that define the exterior dihedral
          angle.  Also sets the term's derivatives "eDa", "eDb", "eDr".
          
          Specifically, the energy is the length of the hinge times
          some concave function of the angle between the projections
          of "a" and "b" on a plane perpendicular to "r".
          The energy is minimum when the angle is zero.
        *)
        BEGIN
          WITH
            h = LR3.Norm(r)
          DO
            IF r = 0.0d0 THEN
              eTerm := 0.0d0; 
              eDa := LR3.T{0.0d0, 0.0d0, 0.0d0};
              eDb := LR3.T{0.0d0, 0.0d0, 0.0d0};
              eDr := LR3.T{0.0d0, 0.0d0, 0.0d0};
            ELSE
              VAR z: LONGREAL; (* Energy density per unit length *)
                  zDa, zDb, zDt: LR3.T;
              BEGIN
                WITH
                  f = 1.0d0/h,
                  t = LR3.Scale(f, r)  
                DO
                  Compute_z_from_a_b_t(a, b, t, z, zDa, zDb, zDt);
                  eTerm := h * z;
                  eDa := LR3.Scale(h, zDa);
                  eDb := LR3.Scale(h, zDb);
                  eDr := LR3.Mix(z - LR3.Dot(eDz, t), t, 1.0d0, eDz)
                END
              END
            END
         END
        END Compute_e_from_a_b_r;
        
      PROCEDURE Compute_z_from_a_b_t(
          READONLY a, b, t: LR3.T;
          VAR z: LONGREAL;
          VAR zDa, zDb, zDt: LR3.T;
        ) =
        (*
          Computes the linear energy density "z" of the hinge (and its
          partial derivatives) given the direction vector "t" of the
          spine and two vectors "a", "b" that define the external
          dihedral angle.
        *)
        BEGIN
          WITH
            at = LR3.Dot(a, t),
            an = LR3.Sub(a, LR3.Scale(at, t)),
            bt = LR3.Dot(b, t),
            bn = LR3.Sub(b, LR3.Scale(bt, t))
          DO
            VAR zDan, zDbn: LR3.T;
            BEGIN
              Compute_z_from_an_bn(an, bn, z, zDan, zDbn);
              WITH
                zDap = LR3.Dot(zDan, t),
                zDbp = LR3.Dot(zDbn, t),
              DO
                zDa := LR3.Mix(1.0d0, zDan, -zDap, t);
                zDb := LR3.Mix(1.0d0, zDbn, -zDbp, t);
                zDt := LR3.Add(
                  LR3.Mix(-at, zDan, -zDap, a),
                  LR3.Mix(-bt, zDbn, -zDap, b)
                )
              END
            END
          END
        END Compute_z_from_a_b_t;
        
      PROCEDURE Compute_z_from_an_bn(
          READONLY an, bn: LR3.T;
          VAR z: LONGREAL;
          VAR zDan, zDbn: LR3.T;
        ) =
        (*
          Computes the energy density "z" given two vectors "an" and "bn",
          orthogonal to the spine, that define the external dihedral angle. *)
        BEGIN
          WITH
            ar = LR3.Norm(an),
            br = LR3.Norm(bn),
            ab = LR3.Dot(an, bn),
            den = ar*br,
            cos = ab / den
          DO
            es := 1.0d0 - cos;
            WITH
              esDcos = -1.0d0,
              esDab = esDcos/den,
              esDden = - esDcos*cos/den,
              esDar = esDden * br,
              esDbr = esDden * ar
            DO
              esDan := LR3.Mix(esDab, bn, esDar/ar, an);
              esDbn := LR3.Mix(esDab, an, esDbr/br, bn);
            END
          END
        END Compute_es_from_an_bn;