//debugging
/**
 double vTot;
 double sForceTot;
 double sCount;
 double sMax;
 double vMax;
 */

class SphereList extends ArrayList<Sphere2D> {
  ArrayList<Integer> displayOrder = new ArrayList<Integer>(); // Order that bases are drawn

  HashMap<Integer, Integer> relaxedIds = new HashMap<Integer, Integer>(); // Set of bases to ignore tension forces on (HashMap used because Set is not supported in Processing)

  public SphereList() {
    super();
    relaxedIds = new HashMap<Integer, Integer>();
    displayOrder = new ArrayList<Integer>();
  }

  SphereList cloneFunction() {
    SphereList spheres = new SphereList();
    for (int i = 0; i < size(); ++i) {
      spheres.add(get(i).cloneFunction());
    }
    for (int i = 0; i < spheres.size(); ++i) {
      if (spheres.get(i).fiveP != null) { // still points to original array
        spheres.get(i).fiveP = spheres.get(spheres.get(i).fiveP.id);
      }
      if (spheres.get(i).threeP != null) { // still points to original array
        spheres.get(i).threeP = spheres.get(spheres.get(i).threeP.id);
      }
      if (spheres.get(i).pairedWForce != null) { // still points to original array
        spheres.get(i).pairedWForce = spheres.get(spheres.get(i).pairedWForce.id);
      }
      if (get(i).pairedH != null) { // still points to original array
        get(i).pairedH = get(get(i).pairedH.id);
      }
      if (get(i).pairedS != null) { // still points to original array
        get(i).pairedS = get(get(i).pairedS.id);
      }
    }
    spheres.relaxedIds = (HashMap<Integer, Integer>)relaxedIds.clone();
    spheres.displayOrder = (ArrayList<Integer>)displayOrder.clone();
    return spheres;
  }

  /** Returns sphere with this id. */
  public Sphere2D get(int id) {
    return (Sphere2D)(super.get(id));
  }

  /** Returns number of spheres. */
  public int size() {
    return super.size();
  }

  // Compute forces, update velocities and positions
  void update() {
    /**
     // debugging
     double bbDiffTot = 0;
     int bbCount = 0;
     double bpDiffTot = 0;
     int bpCount = 0;
     double minPosViol = 100;
     double tensionTot = 0;
     int TCount = 0;
     double tMax = 0;
     double eForceTot = 0;
     int eCount = 0;
     double eMax = 0;
     double vdTot = 0;
     int vdCount = 0;
     double vdMax = 0;
     sForceTot = 0;
     sMax = 0;
     sCount = 0;
     vTot = 0;
     vMax = 0;
     */

    double eCutoff = Math.sqrt(ff.backboneDist) * ff.eCutoffScale;
    if (size() > 1800) {
      eCutoff /= 2;
    }

    // Interactions between residues
    for (int i = 0, l = size(); i < l; ++i) {
      Sphere2D s1 = get(i);

      for (int j = i+1; j < l; ++j) {
        Sphere2D s2 = get(j);
        double d = s1.distance(s2);
        if (d < 1) {
          d = 1; // avoid singularity
        }

        double d2 = d * d;
        double dx = s1.x - s2.x;
        double dy = s1.y - s2.y;

        double dx0 = dx/d;
        double dy0 = dy/d;

        double pushX = (ff.backboneDist * dx0) / d2;
        double pushY = (ff.backboneDist * dy0) / d2;


        // Electrostatic repulsion
        if ((d < eCutoff) && (s1.pairedWForce != s2) && (!s1.isAdjacent(s2))) {
          s1.forceX += ff.eWeight * pushX;
          s1.forceY += ff.eWeight * pushY;

          s2.forceX -= ff.eWeight * pushX;
          s2.forceY -= ff.eWeight * pushY;

          // debugging
          /**
           double eForce = Math.sqrt(Math.pow(ff.eWeight * pushX, 2) + Math.pow(ff.eWeight * pushY, 2));
           eForceTot += eForce;
           if (eForce > eMax) { eMax = eForce; }
           ++eCount;
           */
        }

        // Van Der Waals repulsion
        s1.forceX += ff.vdWeight * pushX;
        s1.forceY += ff.vdWeight * pushY;

        s2.forceX -= ff.vdWeight * pushX;
        s2.forceY -= ff.vdWeight * pushY;

        // debugging
        /**
         double vdForce = Math.sqrt(Math.pow(ff.vdWeight * pushX, 2) + Math.pow(ff.vdWeight * pushY, 2));
         vdTot += vdForce;
         if (vdForce > vdMax) { vdMax = vdForce; }
         ++vdCount;
         */

        // Tension force
        double viol = -1;
        if ( s1.pairedWForce == s2 && !(relaxedIds.containsKey(s1.id) || relaxedIds.containsKey(s2.id))) { //  && !(s1.inHairpin && s2.inHairpin)
          viol = d - ff.bpDist;
          //bpDiffTot += viol; // debuging
          //++bpCount; // debuging
        } else if (s1.isAdjacent(s2)) {
          viol = d - ff.backboneDist;
          //bbDiffTot += viol; // debuging
          //++bbCount; // debuging
        }

        //if (viol > 0 && viol < minPosViol) { minPosViol = viol; } // debugging

        // ?TODO: Add a soft (perhaps scaled somehow) outwards restoring force for when violation is negative?
        if (viol > 0.1) {
          s1.forceX -= ff.stretchWeight * viol * dx0;
          s1.forceY -= ff.stretchWeight * viol * dy0;

          s2.forceX += ff.stretchWeight * viol * dx0;
          s2.forceY += ff.stretchWeight * viol * dy0;

          // debugging
          /**
           double TForce = Math.sqrt(Math.pow(ff.stretchWeight * viol * dx0, 2) + Math.pow(ff.stretchWeight * viol * dy0, 2));
           tensionTot += TForce;
           if (TForce > tMax) { tMax = TForce; }
           ++TCount;
           */
        }
      }
    }

    if (attachedIds.size() == 0 || simulateAllMode) { // simulate all residues
      for (int i = 0; i < size(); ++i) {
        get(i).update();
      }
    } else { // only simulate selected residues
      for (Iterator<Integer> it = attachedIds.keySet().iterator(); it.hasNext(); ) {
        get(it.next()).update();
      }
    }

    // Force debugging
    /**
     if (debugPrints) {
     // Check for 0s in denominators
     double bbAvg = bbCount > 0 ? bbDiffTot / bbCount : 0;
     double bpAvg = bpCount > 0 ? bpDiffTot / bpCount : 0;
     double tAvg = TCount > 0 ? tensionTot / TCount : 0;
     double vAvg = size() > 0 ? vTot / size() : 0;
     double eAvg = eCount > 0 ? eForceTot / eCount : 0;
     double sAvg = sCount > 0 ? sForceTot / sCount : 0;
     double vdAvg = vdCount > 0 ? vdTot / vdCount : 0;
     println("\n\nBackbone: " + ff.backboneDist);
     println("Avg backbone viol: " + bbAvg);
     println("Basepair: " + ff.bpDist);
     println("Avg basepair viol: " + bpAvg);
     println("minPosViol: " + minPosViol);
     println("TensionAvg: " + tAvg + ",  tMax: " + tMax);
     
     println("\nE-Avg: " + eAvg + ",  eMax: " + eMax + ",  eCount: " + eCount + ",  eCutoff: " + eCutoff);
     
     println("\nvd-Avg: " + vdAvg + ",  vdMax: " + vdMax);
     
     println("\nStraight Avg: " + sAvg + ",  sMax: " + sMax);
     
     println("\nV-Avg: " + vAvg + ",  vMax: " + vMax);
     }
     */
  }

  // Make helices rigid by fixing residue coordinates
  void fixHelices() {
    SphereList spheres = this;
    for (int i = 0; i < spheres.size(); ++i) {
      Sphere2D res = spheres.get(i);
      if (res.pairedWForce != null && (!res.hasHelixUp) && res.hasHelixDown && res.id < res.pairedWForce.id) { // found beginning of helix strand. 
        Sphere2D resEnd = res;
        int j;
        for (j = i+1; j < spheres.size(); ++j) { // from here till the end of the helix
          resEnd = spheres.get(j);
          if (!resEnd.hasHelixDown) { // If this is the end of the helix
            break;
          }
        }
        
        //if (rigidHelices) ...

        int len = j - i + 1; // length of helix
        double stacks = len - 1; // how many bonds between residues of the same backbone in one side of the helix
        double dx = resEnd.x - res.x; // Distance between the first and last residues in the helix (vector from start to end)
        double dy = resEnd.y - res.y;
        double dx2 = resEnd.pairedWForce.x - res.pairedWForce.x; // vector running along the other half of the helix (same direction as lower half)
        double dy2 = resEnd.pairedWForce.y - res.pairedWForce.y;
        double dxInterval = dx / stacks; // Step intervals between the first and last residues in the helix
        double dyInterval = dy / stacks;
        double dx2Interval = dx2 / stacks;
        double dy2Interval = dy2 / stacks;
        double distanceSum = 0;
        int count = 0;
        for (int k = 0; k < len; ++k) {
          // SET NEW COORDINATES FOR EACH SPHERE IN BOTH SIDES OF HELIX
          Sphere2D resNew = spheres.get(i + k);
          if ( !relaxedIds.containsKey(resNew.id) && !relaxedIds.containsKey(resNew.pairedWForce.id) ) {
            resNew.x = res.x + (k * dxInterval);
            resNew.y = res.y + (k * dyInterval);

            resNew.pairedWForce.x = res.pairedWForce.x + (k * dx2Interval);
            resNew.pairedWForce.y = res.pairedWForce.y + (k * dy2Interval);
            ++count;
            distanceSum += resNew.pairedWForce.distance(resNew);
          }
        }

        double distanceAvg = 0;
        if (count > 0) {
          distanceAvg = distanceSum / count;
        }

        if (distanceAvg < ff.bpDist) { 
          distanceAvg = ff.bpDist;
        }
        double dxMid = (dx+dx2) / 2.0; // Distance between midpoints of starts and ends of helices
        double dyMid = (dy+dy2) / 2.0;
        PVector midLine = new PVector((float)dxMid, (float)dyMid);
        float angle = midLine.heading(); // Gives angle from 0 to PI running clockwise starting on positive x-axis to negative x-axis,
        // and 0 to -PI running counter clockwise

        double xDifference = (distanceAvg * sin(angle)) / 2.0; // Right angle correction
        double yDifference = (-distanceAvg * cos(angle)) / 2.0;

        double xStartMidpoint = (res.x + res.pairedWForce.x) / 2.0; // Midpoint between the starts of the helix
        double yStartMidpoint = (res.y + res.pairedWForce.y) / 2.0;

        double xMidInterval = dxMid / stacks; // Step interval along midline
        double yMidInterval = dyMid / stacks;
        double xNew, yNew;

        double mul = 0;
        for (int k = 0; k < len; ++k) { // Determine which side of the midline the strand is on
          Sphere2D resNew = spheres.get(i + k);
          double sign = (midLine.x) * (resNew.y - yStartMidpoint) - (midLine.y) * (resNew.x - xStartMidpoint);
          if (sign > 0) {
            mul -= 1;
          } else if (sign < 0) {
            mul += 1;
          }
        }
        if (mul >= 0) {
          mul = 1;
        } else {
          mul = -1;
        }

        // Pull residues to midline
        for (int k = 0; k < len; ++k) {
          Sphere2D resNew = spheres.get(i + k);

          if ( !relaxedIds.containsKey(resNew.id) && !relaxedIds.containsKey(resNew.pairedWForce.id) ) {
            xNew = xStartMidpoint + (k * xMidInterval);
            yNew = yStartMidpoint + (k * yMidInterval);

            if (simulateAllMode || (attachedIds.size() == 0 || attachedIds.containsKey(resNew.id))) {
              resNew.x = xNew + (xDifference * mul);
              resNew.y = yNew + (yDifference * mul);
            }

            if (simulateAllMode || (attachedIds.size() == 0 || attachedIds.containsKey(resNew.pairedWForce.id))) {
              resNew.pairedWForce.x = xNew - (xDifference * mul);
              resNew.pairedWForce.y = yNew - (yDifference * mul);
            }
          }
        }
        
        if (rigidHairpins && !rigidLoops) {
          // Circularize hairpins
          res = resEnd;
          while (res.threeP != null) {
            res = res.threeP;
            if (res.isPaired()) {
              break;
            }
          }

          if (resEnd.pairedWForce.id == res.id && resEnd.id != res.id) {
            //println("\nFound hairpin from residue " + resEnd.id + " - " + res.id);
            //Ensure last bond of the helix is within a range of the basepair distance
            //double dist = res.distance(resEnd);
            double sideL =  res.distance(resEnd); // ff.backboneDist
            //double distDiff = dist - ff.bpDist;
            //distDiff = distDiff < 0 ? distDiff*-1 : distDiff;
            //println("distDiff: " + distDiff);

            //if (distDiff < 10) {

            // ***Circularize Hairpin***
            // println("Circularizing!");

            Sphere2D resStart = resEnd;
            resEnd = res;

            //resStart.inHairpin = true;
            //resEnd.inHairpin = true;

            // Treat it like a regular polygon of n-sides
            int n = (resEnd.id - resStart.id) + 1;
            
            // Vector from end to start
            PVector v1 = new PVector((float)(resStart.x - resEnd.x), (float)(resStart.y - resEnd.y));
            // Vector from fiveP to start
            PVector v2 = new PVector((float)(resStart.x - resStart.fiveP.x), (float)(resStart.y - resStart.fiveP.y));

            double angChange = (2 * PI) / (double) n; // Exterior angle
            //            println("angChange: " + angChange);
            if (v1.y * v2.x > 0) {
              angChange *= -1;
            }

            float ang = v1.heading();
            res = resStart;
            double targetX = res.x;
            double targetY = res.y;

            for (int k = 0; k < n - 2; ++k) {
              ang += angChange;

              double xDiff = sideL * cos(ang);
              double yDiff = sideL * sin(ang);

              targetX += xDiff;
              targetY += yDiff;

              res = res.threeP;

              res.forceX += ff.hairpinWeight * (targetX - res.x);
              res.forceY += ff.hairpinWeight * (targetY - res.y);

              // res.threeP.forceX +=  ff.hairpinWeight * ((res.x + xDiff) - res.threeP.x);  Force correction
              //res.threeP.forceY +=  ff.hairpinWeight * ((res.y + yDiff) - res.threeP.y);

              //res.threeP.x = res.x + xDiff;  Perfect correction
              //res.threeP.y = res.y + yDiff;

              //res = res.threeP;  Place for the other
              //res.inHairpin = true;  Unimplemented idea
            }
          }
        }
      } // if res.pairedForce ...
    } // for int i = 0; i < spheres.size ...
  } // fixHelices
  
  // Recursive algorithm: find helix, search border for next hairpin, if found call again, if not call circularize. Treat previous pairs as residues on the circumference

  //ArrayList<Integer> searchEnds = new ArrayList<Integer>();
  
  void fixLoops() {
    //println("\n\n\n*****STARTING FIX LOOPS*****");
    SphereList spheres = this;
    //searchEnds = new ArrayList<Integer>();
    for (int i = 0; i < spheres.size(); ++i) {
      Sphere2D res = spheres.get(i);
      if (res.pairedWForce != null && res.pairedWForce.strandId == res.strandId && (!res.hasHelixUp) && res.hasHelixDown && res.id < res.pairedWForce.id) { // found beginning of an intra-strand helix
        //println("\n\nCalling findLoop: " + res.id + " , from fixLoops");
        findLoop(res);
        i = res.pairedWForce.id;
      }
    }
  }
  
  void findLoop(Sphere2D helixStart) {
    SphereList spheres = this;
    Sphere2D loopStart = helixStart; // loopStart: residue on the end of the 5' to 3' helix strand
    for (int i = helixStart.id+1; i < spheres.size(); ++i) { // from here till the end of the helix
      loopStart = spheres.get(i);
      if (!loopStart.hasHelixDown) {
        break;
      }
    }
    Sphere2D loopEnd = loopStart.pairedWForce;
    //println("loopStart: " + loopStart.id);
    
    ArrayList<Sphere2D> loopResidues = new ArrayList<Sphere2D>();
    int endSearch = loopEnd.id;
    //searchEnds.add(endSearch); // For ensuring we aren't messing with pseudoknotted helixes
    Sphere2D searchRes = loopStart;
    while (searchRes.id < (endSearch-1)) {
      searchRes = searchRes.threeP;
      loopResidues.add(searchRes);
      if (searchRes.pairedWForce != null && searchRes.id < searchRes.pairedWForce.id) {
        // Interstrand helix or pseudoknot, abort circularization of this loop
        if (searchRes.pairedWForce.strandId != searchRes.strandId || searchRes.pairedWForce.id > endSearch) {
          return;
        }
        /*
        if (searchRes.pairedWForce.strandId != searchRes.strandId || searchRes.pairedWForce.id > searchEnds.get(searchEnds.size() - 1)) {
          searchEnds.remove(searchEnds.size() - 1);
          return;
        }
        */
        
        /*
        for (int i = 0, int l = searchEnds.size(); i < l; ++i) {
          if (searchRes.pairedWForce.id > searchEnds.get(i))
        }
        */
        
        if ((!searchRes.hasHelixUp) && searchRes.hasHelixDown) {
          //println("found nested helix inside start: " + loopStart.id + "  at searchRes " + searchRes.id);
          findLoop(searchRes);
        }
        //println("loop at searchRes: " + searchRes.id);
        searchRes = searchRes.pairedWForce;
        //println("Also adding the inside of helix: " + searchRes.id);
        loopResidues.add(searchRes);
      }
    }
    
    // ***Circularize*** //
    //println("\ncircularizing loop from loopStart " + loopStart.id + " to " + loopEnd.id);
    double dist = loopStart.distance(loopEnd);
    
    // Treat it like a regular polygon of n-sides
    int n = loopResidues.size() + 2;
    double sideL =  ff.backboneDist; //  ff.bpDist;
    
    // Vector from end to start
    PVector v1 = new PVector((float)(loopStart.x - loopEnd.x), (float)(loopStart.y - loopEnd.y));
    // Vector from fiveP to start
    PVector v2 = new PVector((float)(loopStart.x - loopStart.fiveP.x), (float)(loopStart.y - loopStart.fiveP.y));
    
    double angChange = (2 * PI) / (double) n; // Exterior angle
    // println("angChange: " + angChange);
    if (v1.y * v2.x > 0) {
      angChange *= -1;
    }
    
    float ang = v1.heading();
    Sphere2D res = loopStart;
    double targetX = res.x;
    double targetY = res.y;
  
    for (int k = 0; k < n - 2; ++k) {
      ang += angChange;
  
      double xDiff = sideL * cos(ang);
      double yDiff = sideL * sin(ang);
  
      targetX += xDiff;
      targetY += yDiff;
      
      if (loopResidues.get(k).id < res.id) {
        println("\n****ERROR!!!!****\n");
      }
      res = loopResidues.get(k);
      
      //fill(#FF0D0D);
      //ellipse((float)targetX, (float)targetY, 15, 15);
      
      res.forceX += ff.loopWeight * (targetX - res.x);
      res.forceY += ff.loopWeight * (targetY - res.y);
  
      // res.threeP.forceX +=  ff.hairpinWeight * ((res.x + xDiff) - res.threeP.x);  Force correction
      //res.threeP.forceY +=  ff.hairpinWeight * ((res.y + yDiff) - res.threeP.y);
  
      //res.threeP.x = res.x + xDiff;  Perfect correction
      //res.threeP.y = res.y + yDiff;
      
      // Push away and straighten nested stems
      if (res.hasHelixDown) {
        
        float newAng = ang + (float)angChange;
        if (angChange < 0) {
          newAng += PI / 2.0;
        } else {
          newAng -= PI / 2.0;
        }
        
        Sphere2D outerHelixEnd = res.threeP;
        while (outerHelixEnd.hasHelixDown) {
          outerHelixEnd = outerHelixEnd.threeP;
        }
        //println("\n outerHelixEnd: " + outerHelixEnd.id);
        
        int outerHelixLength = outerHelixEnd.id - res.id;
        double helixDiffX = outerHelixLength * sideL * cos(newAng);
        double helixDiffY = outerHelixLength * sideL * sin(newAng);
        
        double outer1TargetX = targetX + helixDiffX;
        double outer1TargetY = targetY + helixDiffY;
        
        Sphere2D outerHelixEnd2 = outerHelixEnd.pairedWForce;
        
        if (angChange < 0) {
          newAng -= PI / 2.0;
        } else {
          newAng += PI / 2.0;
        }
        
        double outer2TargetX = outer1TargetX + (dist * cos(newAng));
        double outer2TargetY = outer1TargetY + (dist * sin(newAng));
        
        /*
        fill(#0DF9FF);
        ellipse((float)outer1TargetX, (float)outer1TargetY, 15, 15);
        ellipse((float)outer2TargetX, (float)outer2TargetY, 15, 15);
        */
        
        outerHelixEnd.forceX += ff.branchStraightenWeight * (outer1TargetX - outerHelixEnd.x);
        outerHelixEnd.forceY += ff.branchStraightenWeight * (outer1TargetY - outerHelixEnd.y);
        
        outerHelixEnd2.forceX += ff.branchStraightenWeight * (outer2TargetX - outerHelixEnd2.x);
        outerHelixEnd2.forceY += ff.branchStraightenWeight * (outer2TargetY - outerHelixEnd2.y);
      }
    }
  }


  // Set all velocities to zero
  void stopVelocities() {
    for (int i = 0; i < size(); ++i) {
      Sphere2D sphere = get(i);
      sphere.vx = 0.0;
      sphere.vy = 0.0;
    }
  }


  /** Returns coordinates of sphere that is closest to given coordinates */
  double[] findClosest(double x, double y) {
    int closestId = 0;
    double dBest = get(0).distance(x, y);
    for (int i = 1; i < size(); ++i) {
      double d = get(i).distance(x, y);
      if (d < dBest) {
        dBest = d;
        closestId = i;
      }
    }
    double[] output = {closestId, dBest};
    return output;
  }

  PVector computeCenterOfMass() {
    cmx = 0.0;
    cmy = 0.0;
    for (int i = 0; i < size(); ++i) {
      cmx += get(i).x;
      cmy += get(i).y;
    }
    cmx /= size(); // center of mass
    cmy /= size();
    return new PVector(cmx, cmy);
  }

  void display() {
    boolean[] drawn = new boolean[displayOrder.size()];
    float displayDiameter = (float) (ff.radius * ds.radiusShowMul);
    float offsetX = ds.textOffsX * displayDiameter;
    float offsetY = ds.textOffsY * displayDiameter;
    drawer.TextSize(ds.textReduction * displayDiameter);
    textAlign(LEFT);
    for (int id : displayOrder) {
      get(id).displayLines(drawn);
      get(id).displayBody(displayDiameter, offsetX, offsetY);
      drawn[id] = true;
    }
  }

  void displayLabels() {
    float radius2 = (float) (ff.radius * ds.radiusShowMul);
    drawer.TextSize(constrain(ds.textReduction * radius2, 6, 65));
    textAlign(CENTER, CENTER);
    drawer.Stroke(#9B9B9B);
    drawer.StrokeWeight(constrain(sliders[0].value * ds.labelWeight, .7, ds.labelWeight));
    drawer.Fill(ds.textCol);

    if (!capturingScreen) {
      for (Iterator<Integer> it = relaxedIds.keySet().iterator(); it.hasNext(); ) {
        PVector pv = findPositionForPrint( it.next(), radius2, radius2 );
        if (pv != null) {
          drawer.Text("R", pv.x, pv.y);
        }
      }
    }

    drawer.Fill(ds.textCol);
    int[] strandStarts = currentState.sim.starts;
    for (int i = 0; i < strandStarts.length; ++i) {
      int start = strandStarts[i];
      int end;
      if (i < strandStarts.length-1) {
        end = strandStarts[i+1];
      } else {
        end = size();
      }
      PVector pv = findPositionForPrint(start, (radius2 + 10), radius2 );
      if (pv != null) {
        String textString = str(1);
        float ntX = (float)get(start).x;
        float ntY = (float)get(start).y;
        float dx = pv.x - ntX;
        float dy = pv.y - ntY;
        float magn = sqrt( dx*dx + dy*dy );
        float normX = dx/magn;
        float normY = dy/magn;
        drawer.Fill(ds.textCol);
        drawer.Line( pv.x-(normX*8), pv.y-(normY*8), ntX + (normX * radius2/2.0), ntY + (normY * radius2/2.0));
        drawer.Text(textString, pv.x, pv.y);
      }
      for (int j = start+9; j < end; j += 10) {
        pv = findPositionForPrint(j, (radius2 + 10), radius2 );
        if (pv != null) {
          String textString = str((j-start)+1);
          float ntX = (float)get(j).x;
          float ntY = (float)get(j).y;
          float dx = pv.x - ntX;
          float dy = pv.y - ntY;
          float magn = sqrt( dx*dx + dy*dy );
          float normX = dx/magn;
          float normY = dy/magn;
          drawer.Fill(ds.textCol);
          drawer.Line( pv.x-(normX*8), pv.y-(normY*8), ntX + (normX * radius2/2.0), ntY + (normY * radius2/2.0));
          drawer.Text(textString, pv.x, pv.y);
        }
      }
    }
  }

  void displayInfo() {
    float radius2 = (float)(ff.radius * ds.radiusShowMul);
    drawer.TextSize(constrain(ds.textReduction * radius2, 6, 65));
    textAlign(CENTER, CENTER);
    drawer.Fill(ds.textCol);

    int[] strandStarts = currentState.sim.starts;
    for (int i = 0; i < strandStarts.length; ++i) {
      int start = strandStarts[i];
      int end;
      if (i < strandStarts.length-1) {
        end = strandStarts[i+1];
      } else {
        end = size();
      }
      int id = int((start + end) / 2);
      PVector pv = findPositionForPrint( id, radius2 * 1.5, radius2 );
      if (pv != null) {
        String textString = "Strand " + (get(id).strandId + 1);
        drawer.Text(textString, pv.x, pv.y);
      }

      // 5' 3' labels
      if (end - start >= 4) {
        pv = findPositionForPrint( start+1, radius2 * 1.5, radius2 );
        if (pv != null) {
          drawer.Text("5'", pv.x, pv.y);
        }
        if ((end - start) % 10 == 0) {
          --end;
        }
        pv = findPositionForPrint( end-1, radius2 * 1.5, radius2 );
        if (pv != null) {
          drawer.Text("3'", pv.x, pv.y);
        }
      }
    }
  }

  PVector findPositionForPrint(int i0, double delta, double targetBuffer) {
    double x = get(i0).x;
    double y = get(i0).y;
    for (int ix = -1; ix <= 1; ++ix) {
      for (int iy = 1; iy >= -1; --iy) {
        if (ix == 0 && iy == 0) { 
          continue;
        }
        double dist = delta;
        if (ix != 0 && iy != 0) { 
          dist *= .707;
        }

        double x2 = x + (ix * dist);
        double y2 = y + (iy * dist);
        boolean violation = false;
        for (int ii = 0; ii < size(); ++ii) {
          if (ii == i0) { 
            continue;
          }
          Sphere2D sphere = get(ii);
          double d = normFunction(x2-sphere.x, y2-sphere.y, 0);
          if (d < targetBuffer) {
            violation = true;
            break;
          }
        }
        if (!violation) {
          return new PVector((float)x2, (float)y2);
        }
      }
    }
    return null;
  }
}