]> git.donarmstrong.com Git - lilypond.git/blobdiff - lily/beam-quanting.cc
Run `make grand-replace'.
[lilypond.git] / lily / beam-quanting.cc
index 178d8da48437e271adb9f0f22f1cf1f56903fd4a..cf933d500bc3d3071d1386fdb172e8391bb51465 100644 (file)
@@ -1,42 +1,57 @@
 /*
   beam-quanting.cc -- implement Beam quanting functions
-  
+
   source file of the GNU LilyPond music typesetter
-  
-  (c)  1997--2003 Han-Wen Nienhuys <hanwen@cs.uu.nl>
-  Jan Nieuwenhuizen <janneke@gnu.org>
-  
 
-  
+  (c) 1997--2008 Han-Wen Nienhuys <hanwen@xs4all.nl>
+  Jan Nieuwenhuizen <janneke@gnu.org>
 */
 
+#include "beam.hh"
 
+#include <algorithm>
+using namespace std;
 
-#include <math.h>
-
-#include "warn.hh"
 #include "grob.hh"
+#include "align-interface.hh"
+#include "international.hh"
+#include "output-def.hh"
+#include "pointer-group-interface.hh"
 #include "staff-symbol-referencer.hh"
-#include "beam.hh"
 #include "stem.hh"
-#include "paper-def.hh"
-#include "group-interface.hh"
-#include "align-interface.hh"
+#include "warn.hh"
+#include "main.hh"
 
-const int INTER_QUANT_PENALTY = 1000;
-const Real SECONDARY_BEAM_DEMERIT  = 10.0;
-const int STEM_LENGTH_DEMERIT_FACTOR = 5;
+Real
+get_detail (SCM alist, SCM sym, Real def)
+{
+  SCM entry = scm_assq (sym, alist);
 
-// possibly ridiculous, but too short stems just won't do
-const int STEM_LENGTH_LIMIT_PENALTY = 5000;
-const int DAMPING_DIRECTION_PENALTY = 800;
-const int MUSICAL_DIRECTION_FACTOR = 400;
-const int IDEAL_SLOPE_FACTOR = 10;
+  if (scm_is_pair (entry))
+    return robust_scm2double (scm_cdr (entry), def);
+  return def;
+}
 
-const Real ROUND_TO_ZERO_SLOPE = 0.05;
-const int ROUND_TO_ZERO_POINTS = 4;
+void
+Beam_quant_parameters::fill (Grob *him)
+{
+  SCM details = him->get_property ("details");
 
-extern bool debug_beam_quanting_flag;
+  /* 
+     TODO: The default values should be copied to define-grobs.scm.
+   */
+  INTER_QUANT_PENALTY = get_detail (details, ly_symbol2scm ("inter-quant-penalty"), 1000.0);
+  SECONDARY_BEAM_DEMERIT = get_detail (details, ly_symbol2scm ("secondary-beam-demerit"), 10.0);
+  STEM_LENGTH_DEMERIT_FACTOR = get_detail (details, ly_symbol2scm ("stem-length-demerit-factor"), 5);
+  REGION_SIZE = get_detail (details, ly_symbol2scm ("region-size"), 2);
+  BEAM_EPS = get_detail (details, ly_symbol2scm ("beam-eps"), 1e-3);
+  STEM_LENGTH_LIMIT_PENALTY = get_detail (details, ly_symbol2scm ("stem-length-limit-penalty"), 5000);
+  DAMPING_DIRECTION_PENALTY = get_detail (details, ly_symbol2scm ("damping-direction-penalty"), 800);
+  HINT_DIRECTION_PENALTY = get_detail (details, ly_symbol2scm ("hint-direction-penalty"), 20);
+  MUSICAL_DIRECTION_FACTOR = get_detail (details, ly_symbol2scm ("musical-direction-factor"), 400);
+  IDEAL_SLOPE_FACTOR = get_detail (details, ly_symbol2scm ("ideal-slope-factor"), 10);
+  ROUND_TO_ZERO_SLOPE = get_detail (details, ly_symbol2scm ("round-to-zero-slope"), 0.02);
+}
 
 static Real
 shrink_extra_weight (Real x, Real fac)
@@ -44,331 +59,336 @@ shrink_extra_weight (Real x, Real fac)
   return fabs (x) * ((x < 0) ? fac : 1.0);
 }
 
-
 struct Quant_score
 {
   Real yl;
   Real yr;
   Real demerits;
 
-#if DEBUG_QUANTING
-  String score_card_;
+#if DEBUG_BEAM_SCORING
+  string score_card_;
 #endif
 };
 
-
 /*
   TODO:
-  
-   - Make all demerits customisable
 
-   - One sensible check per demerit (what's this --hwn)
+  - Make all demerits customisable
 
-   - Add demerits for quants per se, as to forbid a specific quant
-     entirely
+  - One sensible check per demerit (what's this --hwn)
 
+  - Add demerits for quants per se, as to forbid a specific quant
+  entirely
 */
 
-int best_quant_score_idx (Array<Quant_score>  const & qscores)
+int
+best_quant_score_idx (vector<Quant_score> const &qscores)
 {
   Real best = 1e6;
   int best_idx = -1;
-  for (int i = qscores.size (); i--;)
+  for (vsize i = qscores.size (); i--;)
     {
       if (qscores[i].demerits < best)
        {
-         best = qscores [i].demerits ;
+         best = qscores [i].demerits;
          best_idx = i;
        }
     }
 
   return best_idx;
 }
-  
-MAKE_SCHEME_CALLBACK (Beam, quanting, 1);
+
+MAKE_SCHEME_CALLBACK (Beam, quanting, 2);
 SCM
-Beam::quanting (SCM smob)
+Beam::quanting (SCM smob, SCM posns)
 {
   Grob *me = unsmob_grob (smob);
 
-  SCM s = me->get_grob_property ("positions");
-  Real yl = gh_scm2double (gh_car (s));
-  Real yr = gh_scm2double (gh_cdr (s));
+  Beam_quant_parameters parameters;
+  parameters.fill (me);
 
+  Real yl = scm_to_double (scm_car (posns));
+  Real yr = scm_to_double (scm_cdr (posns));
 
   /*
     Calculations are relative to a unit-scaled staff, i.e. the quants are
     divided by the current staff_space.
-    
-   */
+  */
   Real ss = Staff_symbol_referencer::staff_space (me);
-  Real thickness = Beam::get_thickness (me) / ss ;
+  Real thickness = Beam::get_thickness (me) / ss;
   Real slt = Staff_symbol_referencer::line_thickness (me) / ss;
 
-  SCM sdy = me->get_grob_property ("least-squares-dy");
-  Real dy_mus = gh_number_p (sdy) ? gh_scm2double (sdy) : 0.0;
-  
+  Real dy_mus = robust_scm2double (me->get_property ("least-squares-dy"), 0);
   Real straddle = 0.0;
   Real sit = (thickness - slt) / 2;
   Real inter = 0.5;
   Real hang = 1.0 - (thickness - slt) / 2;
   Real quants [] = {straddle, sit, inter, hang };
-  
-  int num_quants = int (sizeof (quants)/sizeof (Real));
-  Array<Real> quantsl;
-  Array<Real> quantsr;
+
+  int num_quants = int (sizeof (quants) / sizeof (Real));
+  vector<Real> quantsl;
+  vector<Real> quantsr;
 
   /*
     going to REGION_SIZE == 2, yields another 0.6 second with
     wtk1-fugue2.
 
-
-    (result indexes between 70 and 575)  ? --hwn. 
+    (result indexes between 70 and 575)  ? --hwn.
 
   */
 
-
-  
   /*
     Do stem computations.  These depend on YL and YR linearly, so we can
     precompute for every stem 2 factors.
-   */
-  Link_array<Grob> stems=
-    Pointer_group_interface__extract_grobs (me, (Grob*)0, "stems");
-  Array<Stem_info> stem_infos;
-  Array<Real> base_lengths;
-  Array<Real> stem_xposns;  
+  */
+  vector<Grob*> stems
+    = extract_grob_array (me, "stems");
+  vector<Stem_info> stem_infos;
+  vector<Real> base_lengths;
+  vector<Real> stem_xposns;
 
-  Drul_array<bool> dirs_found(0,0);
+  Drul_array<bool> dirs_found (0, 0);
   Grob *common[2];
   for (int a = 2; a--;)
-    common[a] = common_refpoint_of_array (stems, me, Axis(a));
+    common[a] = common_refpoint_of_array (stems, me, Axis (a));
 
-  Grob * fvs = first_visible_stem (me);
-  Grob *lvs = last_visible_stem (me);
+  Grob *fvs = first_normal_stem (me);
+  Grob *lvs = last_normal_stem (me);
   Real xl = fvs ? fvs->relative_coordinate (common[X_AXIS], X_AXIS) : 0.0;
   Real xr = fvs ? lvs->relative_coordinate (common[X_AXIS], X_AXIS) : 0.0;
 
   /*
-    We store some info to quickly interpolate.
-
-    Sometimes my head is screwed on backwards.  The stemlength are
-    AFFINE linear in YL and YR. If YL == YR == 0, then we might have
+    We store some info to quickly interpolate.  The stemlength are
+    affine linear in YL and YR. If YL == YR == 0, then we might have
     stem_y != 0.0, when we're cross staff.
-    
-   */
-  for (int i= 0; i < stems.size(); i++)
+
+  */
+  for (vsize i = 0; i < stems.size (); i++)
     {
-      Grob*s = stems[i];
+      Grob *s = stems[i];
 
       Stem_info si (Stem::get_stem_info (s));
       si.scale (1 / ss);
-      stem_infos.push (si);
-      dirs_found[stem_infos.top ().dir_] = true;
+      stem_infos.push_back (si);
+      dirs_found[stem_infos.back ().dir_] = true;
 
-      bool f = to_boolean (s->get_grob_property ("french-beaming"))
-        && s != lvs && s != fvs;
+      bool f = to_boolean (s->get_property ("french-beaming"))
+       && s != lvs && s != fvs;
 
-      base_lengths.push (calc_stem_y (me, s, common, xl, xr,
-                                     Interval (0,0), f) / ss);
-      stem_xposns.push (s->relative_coordinate (common[X_AXIS], X_AXIS));
-    }
+      if (Stem::is_normal_stem (s))
+       {
+         base_lengths.push_back (calc_stem_y (me, s, common, xl, xr, CENTER, 
+                                              Interval (0, 0), f) / ss);
+       }
+      else
+       {
+         base_lengths.push_back (0);
+       }
 
-  bool xstaff= false;
-  if (lvs && fvs)
-    {
-      Grob *commony = fvs->common_refpoint (lvs, Y_AXIS);
-      xstaff = Align_interface::has_interface (commony);
+      stem_xposns.push_back (s->relative_coordinate (common[X_AXIS], X_AXIS));
     }
-  
+  bool xstaff = Align_interface::has_interface (common[Y_AXIS]);
+
   Direction ldir = Direction (stem_infos[0].dir_);
-  Direction rdir = Direction (stem_infos.top ().dir_);
-  bool knee_b = dirs_found[LEFT] && dirs_found[RIGHT];
+  Direction rdir = Direction (stem_infos.back ().dir_);
+  bool is_knee = dirs_found[LEFT] && dirs_found[RIGHT];
 
+  int region_size = (int) parameters.REGION_SIZE;
 
-  int region_size = REGION_SIZE;
   /*
-    Knees are harder, lets try some more possibilities for knees. 
-   */
-  if (knee_b)
+    Knees are harder, lets try some more possibilities for knees.
+  */
+  if (is_knee)
     region_size += 2;
 
   /*
     Asymetry ? should run to <= region_size ?
-   */
-  for (int i = -region_size ; i < region_size; i++)
+  */
+  for (int i = -region_size; i < region_size; i++)
     for (int j = 0; j < num_quants; j++)
       {
-       quantsl.push (i + quants[j] + int (yl));
-       quantsr.push (i + quants[j] + int (yr));
+       quantsl.push_back (i + quants[j] + int (yl));
+       quantsr.push_back (i + quants[j] + int (yr));
       }
 
-  Array<Quant_score> qscores;
-  
-  for (int l =0; l < quantsl.size (); l++)  
-    for (int r =0; r < quantsr.size (); r++)
+  vector<Quant_score> qscores;
+
+  for (vsize l = 0; l < quantsl.size (); l++)
+    for (vsize r = 0; r < quantsr.size (); r++)
       {
        Quant_score qs;
        qs.yl = quantsl[l];
        qs.yr = quantsr[r];
        qs.demerits = 0.0;
-       
-       qscores.push (qs);
+
+       qscores.push_back (qs);
       }
 
   /* This is a longish function, but we don't separate this out into
      neat modular separate subfunctions, as the subfunctions would be
      called for many values of YL, YR. By precomputing various
      parameters outside of the loop, we can save a lot of time. */
-  for (int i = qscores.size (); i--;)
+  for (vsize i = qscores.size (); i--;)
     {
-      Real d =  score_slopes_dy (qscores[i].yl, qscores[i].yr,
-                                dy_mus, yr- yl, 
-                                xr - xl,
-                                xstaff);
+      Real d = score_slopes_dy (qscores[i].yl, qscores[i].yr,
+                               dy_mus, yr- yl,
+                               xr - xl,
+                               xstaff, &parameters);
       qscores[i].demerits += d;
 
-#if DEBUG_QUANTING
-      qscores[i].score_card_ += to_string ("S%.2f",d);
+#if DEBUG_BEAM_SCORING
+      qscores[i].score_card_ += to_string ("S%.2f", d);
 #endif
     }
 
   Real rad = Staff_symbol_referencer::staff_radius (me);
-  int beam_count = get_beam_count (me);
+  Drul_array<int> edge_beam_counts
+    (Stem::beam_multiplicity (stems[0]).length () + 1,
+     Stem::beam_multiplicity (stems.back ()).length () + 1);
+
   Real beam_translation = get_beam_translation (me) / ss;
 
-  Real reasonable_score = (knee_b) ? 200000 : 100;
-  for (int i = qscores.size (); i--;)
+  Real reasonable_score = (is_knee) ? 200000 : 100;
+  for (vsize i = qscores.size (); i--;)
     if (qscores[i].demerits < reasonable_score)
       {
        Real d = score_forbidden_quants (qscores[i].yl, qscores[i].yr,
-                                    rad, slt, thickness, beam_translation,
-                                    beam_count, ldir, rdir); 
+                                        rad, slt, thickness, beam_translation,
+                                        edge_beam_counts, ldir, rdir, &parameters);
        qscores[i].demerits += d;
 
-#if DEBUG_QUANTING
+#if DEBUG_BEAM_SCORING
        qscores[i].score_card_ += to_string (" F %.2f", d);
 #endif
       }
 
-  for (int i = qscores.size (); i--;)
+  for (vsize i = qscores.size (); i--;)
     if (qscores[i].demerits < reasonable_score)
       {
-       Real d=score_stem_lengths (stems, stem_infos,
-                                base_lengths, stem_xposns,
-                                xl, xr,
-                                knee_b,
-                                qscores[i].yl, qscores[i].yr);
-       qscores[i].demerits +=  d;
-
-#if DEBUG_QUANTING
+       Real d = score_stem_lengths (stems, stem_infos,
+                                    base_lengths, stem_xposns,
+                                    xl, xr,
+                                    is_knee,
+                                    qscores[i].yl, qscores[i].yr, &parameters);
+       qscores[i].demerits += d;
+
+#if DEBUG_BEAM_SCORING
        qscores[i].score_card_ += to_string (" L %.2f", d);
 #endif
       }
 
   int best_idx = best_quant_score_idx (qscores);
 
-
-#if DEBUG_QUANTING
-  SCM inspect_quants = me->get_grob_property ("inspect-quants");
-  if (debug_beam_quanting_flag
-      && gh_pair_p (inspect_quants))
+#if DEBUG_BEAM_SCORING
+  SCM inspect_quants = me->get_property ("inspect-quants");
+  if (to_boolean (me->layout ()->lookup_variable (ly_symbol2scm ("debug-beam-scoring")))
+      && scm_is_pair (inspect_quants))
     {
       Drul_array<Real> ins = ly_scm2interval (inspect_quants);
 
-      int i = 0;
-
       Real mindist = 1e6;
-      for (; i < qscores.size(); i ++)
+      for (vsize i = 0; i < qscores.size (); i++)
        {
-         Real d =fabs (qscores[i].yl- ins[LEFT]) + fabs (qscores[i].yr - ins[RIGHT]);
+         Real d = fabs (qscores[i].yl- ins[LEFT]) + fabs (qscores[i].yr - ins[RIGHT]);
          if (d < mindist)
            {
              best_idx = i;
-             mindist= d;
+             mindist = d;
            }
        }
       if (mindist > 1e5)
-       programming_error ("Could not find quant.");
+       programming_error ("cannot find quant");
     }
 #endif
+
+  Interval final_positions;
+  if (best_idx < 0)
+    {
+      warning (_ ("no feasible beam position"));
+      final_positions = Interval (0, 0);
+    }
+  else
+    {
+      final_positions = Drul_array<Real> (qscores[best_idx].yl,
+                                         qscores[best_idx].yr);
+    }
   
-  me->set_grob_property ("positions",
-                        ly_interval2scm (Drul_array<Real> (qscores[best_idx].yl,
-                                         qscores[best_idx].yr)));
-#if DEBUG_QUANTING
-  if (debug_beam_quanting_flag)
+#if DEBUG_BEAM_SCORING
+  if (best_idx >= 0
+      && to_boolean (me->layout ()->lookup_variable (ly_symbol2scm ("debug-beam-scoring"))))
     {
       qscores[best_idx].score_card_ += to_string ("i%d", best_idx);
-      
+
       // debug quanting
-      me->set_grob_property ("quant-score",
-                            scm_makfrom0str (qscores[best_idx].score_card_.to_str0 ()));
+      me->set_property ("quant-score",
+                       ly_string2scm (qscores[best_idx].score_card_));
     }
 #endif
 
-  return SCM_UNSPECIFIED;
+  return ly_interval2scm (final_positions);
 }
 
 Real
-Beam::score_stem_lengths (Link_array<Grob> const &stems,
-                         Array<Stem_info> const &stem_infos,
-                         Array<Real> const &base_stem_ys,
-                         Array<Real> const &stem_xs,
-                         Real xl, Real xr, 
-                         bool knee, 
-                         Real yl, Real yr)
+Beam::score_stem_lengths (vector<Grob*> const &stems,
+                         vector<Stem_info> const &stem_infos,
+                         vector<Real> const &base_stem_ys,
+                         vector<Real> const &stem_xs,
+                         Real xl, Real xr,
+                         bool knee,
+                         Real yl, Real yr,
+
+                         Beam_quant_parameters const *parameters)
 {
-  Real limit_penalty = STEM_LENGTH_LIMIT_PENALTY;
+  Real limit_penalty = parameters->STEM_LENGTH_LIMIT_PENALTY;
   Drul_array<Real> score (0, 0);
   Drul_array<int> count (0, 0);
-  
-  for (int i=0; i < stems.size (); i++)
+
+  for (vsize i = 0; i < stems.size (); i++)
     {
-      Grobs = stems[i];
-      if (Stem::invisible_b (s))
+      Grob *s = stems[i];
+      if (!Stem::is_normal_stem (s))
        continue;
 
       Real x = stem_xs[i];
-      Real dx = xr-xl;
-      Real beam_y = dx ? yr *(x - xl)/dx + yl * ( xr - x)/dx : (yr + yl)/2;
+      Real dx = xr - xl;
+      Real beam_y = dx ? yr * (x - xl) / dx + yl * (xr - x) / dx : (yr + yl) / 2;
       Real current_y = beam_y + base_stem_ys[i];
-      Real length_pen = STEM_LENGTH_DEMERIT_FACTOR;
-      
+      Real length_pen = parameters->STEM_LENGTH_DEMERIT_FACTOR;
+
       Stem_info info = stem_infos[i];
       Direction d = info.dir_;
 
-      score[d] += limit_penalty * (0 >? (d * (info.shortest_y_ - current_y)));
-      
+      score[d] += limit_penalty * max (0.0, (d * (info.shortest_y_ - current_y)));
+
       Real ideal_diff = d * (current_y - info.ideal_y_);
       Real ideal_score = shrink_extra_weight (ideal_diff, 1.5);
-      
+
       /* We introduce a power, to make the scoring strictly
          convex. Otherwise a symmetric knee beam (up/down/up/down)
          does not have an optimum in the middle. */
       if (knee)
        ideal_score = pow (ideal_score, 1.1);
-      
+
       score[d] += length_pen * ideal_score;
 
-      count[d] ++;
+      count[d]++;
     }
 
   Direction d = DOWN;
   do
-    { 
-      score[d] /= (count[d] >? 1);
-    }
+    score[d] /= max (count[d], 1);
   while (flip (&d) != DOWN);
 
-  return score[LEFT]+score[RIGHT];
+  return score[LEFT] + score[RIGHT];
 }
 
 Real
 Beam::score_slopes_dy (Real yl, Real yr,
                       Real dy_mus, Real dy_damp,
                       Real dx,
-                      bool xstaff)
+                      bool xstaff,
+
+                      Beam_quant_parameters const *parameters)
 {
   Real dy = yr - yl;
   Real dem = 0.0;
@@ -379,35 +399,35 @@ Beam::score_slopes_dy (Real yl, Real yr,
 
     TODO: find a way to incorporate the complexity of the beam in this
     penalty.
-   */
-  if (fabs (dy/dx) > ROUND_TO_ZERO_SLOPE
-      && sign (dy_damp) != sign (dy))
+  */
+  if (sign (dy_damp) != sign (dy))
     {
-      dem += DAMPING_DIRECTION_PENALTY;
+      if (!dy)
+       {
+         if (fabs (dy_damp / dx) > parameters->ROUND_TO_ZERO_SLOPE)
+           dem += parameters->DAMPING_DIRECTION_PENALTY;
+         else
+           dem += parameters->HINT_DIRECTION_PENALTY;
+       }
+      else
+       dem += parameters->DAMPING_DIRECTION_PENALTY;
     }
+  
+  dem += parameters->MUSICAL_DIRECTION_FACTOR
+    * max (0.0, (fabs (dy) - fabs (dy_mus)));
 
-   dem += MUSICAL_DIRECTION_FACTOR * (0 >? (fabs (dy) - fabs (dy_mus)));
-
-
-   Real slope_penalty = IDEAL_SLOPE_FACTOR;
+  Real slope_penalty = parameters->IDEAL_SLOPE_FACTOR;
 
-   /* Xstaff beams tend to use extreme slopes to get short stems. We
-      put in a penalty here. */
-   if (xstaff)
-     slope_penalty *= 10;
+  /* Xstaff beams tend to use extreme slopes to get short stems. We
+     put in a penalty here. */
+  if (xstaff)
+    slope_penalty *= 10;
 
-   /* Huh, why would a too steep beam be better than a too flat one ? */
-   dem += shrink_extra_weight (fabs (dy_damp) - fabs (dy), 1.5)
-     * slope_penalty;
+  /* Huh, why would a too steep beam be better than a too flat one ? */
+  dem += shrink_extra_weight (fabs (dy_damp) - fabs (dy), 1.5)
+    * slope_penalty;
 
-   /*
-     almost zero slopes look like errors in horizontal beams. 
-    */
-   if (fabs (dy) > 1e-3
-       && fabs (dy / dx) < ROUND_TO_ZERO_SLOPE)
-     dem += ROUND_TO_ZERO_POINTS;
-   
-   return dem;
+  return dem;
 }
 
 static Real
@@ -416,124 +436,106 @@ my_modf (Real x)
   return x - floor (x);
 }
 
-
 /*
   TODO: The fixed value SECONDARY_BEAM_DEMERIT is probably flawed:
-   because for 32nd and 64th beams the forbidden quants are relatively
-   more important than stem lengths.
+  because for 32nd and 64th beams the forbidden quants are relatively
+  more important than stem lengths.
 */
 Real
 Beam::score_forbidden_quants (Real yl, Real yr,
                              Real radius,
                              Real slt,
                              Real thickness, Real beam_translation,
-                             int beam_count,
-                             Direction ldir, Direction rdir)
+                             Drul_array<int> beam_counts,
+                             Direction ldir, Direction rdir,
+
+                             Beam_quant_parameters const *parameters)
 {
   Real dy = yr - yl;
-  Drul_array<Real> y(yl,yr);
-  Drul_array<Direction> dirs(ldir,rdir);
-  
-  Real extra_demerit = SECONDARY_BEAM_DEMERIT / beam_count;
+  Drul_array<Real> y (yl, yr);
+  Drul_array<Direction> dirs (ldir, rdir);
 
-  /*
-    Inside the staff, inter quants are forbidden.
-   */
-  Real dem = 0.0;
-  Direction d = LEFT;
-  do
-    {
-      if (fabs (y[d]) <= (radius + 0.5) && fabs (my_modf (y[d]) - 0.5) < 1e-3)
-       dem += INTER_QUANT_PENALTY;
-    }
-  while ((flip (&d))!= LEFT); 
+  Real extra_demerit = parameters->SECONDARY_BEAM_DEMERIT / (max (beam_counts[LEFT], beam_counts[RIGHT]));
 
+  Direction d = LEFT;
+  Real dem = 0.0;
+  Real eps = parameters->BEAM_EPS;
 
-  for (int j = 1; j <= beam_count; j++)
+  do
     {
-      do
+      for (int j = 1; j <= beam_counts[d]; j++)
        {
-         /*
-           see if the outer staffline falls in a beam-gap
-           
-           This test is too weak; we should really check all lines.
-          */
          Direction stem_dir = dirs[d];
-         Real gap1 =  y[d] - stem_dir * ((j-1) * beam_translation + thickness / 2 - slt/2 );
-         Real gap2 = y[d] - stem_dir * (j * beam_translation - thickness / 2 + slt/2);
+
+         /*
+           The 2.2 factor is to provide a little leniency for
+           borderline cases. If we do 2.0, then the upper outer line
+           will be in the gap of the (2, sit) quant, leading to a
+           false demerit.
+         */
+         Real gap1 = y[d] - stem_dir * ((j - 1) * beam_translation + thickness / 2 - slt / 2.2);
+         Real gap2 = y[d] - stem_dir * (j * beam_translation - thickness / 2 + slt / 2.2);
 
          Interval gap;
          gap.add_point (gap1);
          gap.add_point (gap2);
 
-         if (gap.contains (radius))
-           dem += extra_demerit;
+         for (Real k = -radius;
+              k <= radius + eps; k += 1.0)
+           if (gap.contains (k))
+             {
+               Real dist = min (fabs (gap[UP] - k), fabs (gap[DOWN] - k));
+
+               /*
+                 this parameter is tuned to grace-stem-length.ly
+               */
+               Real fixed_demerit = 0.4;
+
+               dem += extra_demerit
+                 * (fixed_demerit
+                    + (1 - fixed_demerit) * (dist / gap.length ()) * 2);
+             }
        }
-      while ((flip (&d))!= LEFT); 
     }
+  while ((flip (&d)) != LEFT);
 
-
-  
-  // todo: use beam_count of outer stems.
-  if (beam_count >= 2)
+  if (max (beam_counts[LEFT], beam_counts[RIGHT]) >= 2)
     {
       Real straddle = 0.0;
       Real sit = (thickness - slt) / 2;
       Real inter = 0.5;
       Real hang = 1.0 - (thickness - slt) / 2;
 
-      Real eps = 1e-3;
-
-      // hmm, without Interval/Drul_array, you get ~ 4x same code...
-      if (fabs (y[LEFT] - dirs[LEFT] * beam_translation) < radius + inter)
-       {
-         if (dirs[LEFT] == UP && dy <= eps
-             && fabs (my_modf (y[LEFT]) - sit) < eps)
-           dem += extra_demerit;
-         
-         if (dirs[LEFT] == DOWN && dy >= eps
-             && fabs (my_modf (y[LEFT]) - hang) < eps)
-           dem += extra_demerit;
-       }
-
-      if (fabs (y[RIGHT] - dirs[RIGHT] * beam_translation) < radius + inter)
-       {
-         if (dirs[RIGHT] == UP && dy >= eps
-             && fabs (my_modf (y[RIGHT]) - sit) < eps)
-           dem += extra_demerit;
-         
-         if (dirs[RIGHT] == DOWN && dy <= eps
-             && fabs (my_modf (y[RIGHT]) - hang) < eps)
-           dem += extra_demerit;
-       }
-      
-      if (beam_count >= 3)
+      Direction d = LEFT;
+      do
        {
-         if (fabs (y[LEFT] - 2 * dirs[LEFT] * beam_translation) < radius + inter)
+         if (beam_counts[d] >= 2
+             && fabs (y[d] - dirs[d] * beam_translation) < radius + inter)
            {
-             if (dirs[LEFT] == UP && dy <= eps
-                 && fabs (my_modf (y[LEFT]) - straddle) < eps)
+             if (dirs[d] == UP && dy <= eps
+                 && fabs (my_modf (y[d]) - sit) < eps)
                dem += extra_demerit;
-             
-             if (dirs[LEFT] == DOWN && dy >= eps
-                 && fabs (my_modf (y[LEFT]) - straddle) < eps)
+
+             if (dirs[d] == DOWN && dy >= eps
+                 && fabs (my_modf (y[d]) - hang) < eps)
                dem += extra_demerit;
            }
-         
-         if (fabs (y[RIGHT] - 2 * dirs[RIGHT] * beam_translation) < radius + inter)
+
+         if (beam_counts[d] >= 3
+             && fabs (y[d] - 2 * dirs[d] * beam_translation) < radius + inter)
            {
-             if (dirs[RIGHT] == UP && dy >= eps
-                 && fabs (my_modf (y[RIGHT]) - straddle) < eps)
+             if (dirs[d] == UP && dy <= eps
+                 && fabs (my_modf (y[d]) - straddle) < eps)
                dem += extra_demerit;
-             
-             if (dirs[RIGHT] == DOWN && dy <= eps
-                 && fabs (my_modf (y[RIGHT]) - straddle) < eps)
+
+             if (dirs[d] == DOWN && dy >= eps
+                 && fabs (my_modf (y[d]) - straddle) < eps)
                dem += extra_demerit;
            }
        }
+      while (flip (&d) != LEFT);
     }
-  
+
   return dem;
 }
 
-