class Sphere2D { // Each residue
  int id; // position in the sphereList (based on the original file initialization)
  int relId; // id relative to strand start
  int strandId = 0; // sequence it belongs to (starting on 0)
  String residue = ""; // Residue letter
  Sphere2D fiveP; // five prime backbone adjacent residue
  Sphere2D threeP; // three prime adjacent residue
  Sphere2D pairedWForce; // pair that applies force
  Sphere2D pairedH; // paired Hoogsteen residue DEPRACATE
  Sphere2D pairedS; // paired Sugar residue DEPRACATE
  boolean hasHelixUp = false; // true if this residue is paired and the adjacent previous residue on the same strand (1 lower strandId) is paired to the adjacent next residue (1 higher strandId) of the residue this is paired with.
  boolean hasHelixDown = false;
  // boolean inHairpin = false; // For hairpin optimization.
  double vx=0, vy=0, vz=0; // velocities
  double x, y, z=0; // coordinates
  double forceX=0, forceY=0, forceZ=0;
  float colorNum = 0;
  color baseColor = #FFFFFF;
  String[] bonds; // stores bond type at that location

  Sphere2D(int _id, int _strandId, String _residue, String[] _bonds) {
    this.id = _id;
    this.strandId = _strandId;
    this.residue = _residue;
    this.bonds = _bonds;
  }

  Sphere2D(int _id, int _strandId, String _residue, double xpos, double ypos, double zpos) {
    this.id = _id;
    this.strandId = _strandId;
    this.residue = _residue;
    this.x = xpos;
    this.y = ypos;
    this.z = zpos;
  }

  Sphere2D cloneFunction() {
    Sphere2D result = new Sphere2D(id, strandId, residue, x, y, z);
    result.fiveP = this.fiveP;
    result.threeP = this.threeP;
    result.hasHelixUp = this.hasHelixUp;
    result.hasHelixDown = this.hasHelixDown;
    result.relId = this.relId;
    result.pairedWForce = this.pairedWForce;
    result.pairedH = this.pairedH;
    result.pairedS = this.pairedS;
    result.vx = this.vx;
    result.vy = this.vy;
    result.colorNum = this.colorNum;
    result.baseColor = this.baseColor;
    result.forceX = this.forceX;
    result.forceY = this.forceY;
    result.bonds = stringArrayClone(this.bonds);
    return result;
  }

  boolean isPaired(Sphere2D other) {
    if (other == null) {
      return false;
    }
    return (pairedWForce == other) || (pairedH == other) || (pairedS == other);
  }

  boolean isPaired() {
    return (pairedWForce != null) || (pairedH != null) || (pairedS != null);
  }

  boolean isAdjacent(Sphere2D other) {
    if (other == null) {
      return false;
    }
    return (fiveP == other) || (threeP == other);
  }

  // Distance to another sphere
  double distance(Sphere2D other) {
    double dx = other.x - x;
    double dy = other.y - y;
    //double dz = other.z - z;
    return normFunction(dx, dy, 0);
  }

  // Distance to a point
  double distance(double otherX, double otherY) {
    double dx = otherX - x;
    double dy = otherY - y;
    //double dz = otherZ - z;
    return normFunction(dx, dy, 0);
  }

  // Rotate by angle (in radians) about center of circle rx, ry
  void rotateSphere(double angle, double rx, double ry) {
    double x2 = x - rx;
    double y2 = y - ry;
    double ca = Math.cos(angle);
    double sa = Math.sin(angle);
    double x2b = x2 * ca - y2 * sa;
    double y2b = x2 * sa + y2 * ca; 
    x = x2b + rx;
    y = y2b + ry;
  }

  
  // Compute force on this sphere by every other outside force and every other sphere, update velocity and position accordingly.
  void update() {
    // Stay within boundaries:
    if (ds.padding > 0) { // slider at 0 turns off all padding
      if (x < ds.padding) {
        forceX += ff.padWeight * ( ds.padding - x );
      } else if (x > (width - ds.padding)) {
        forceX -= ff.padWeight * (x - (width - ds.padding));
      }
      if (y < ds.padding) {
        forceY += ff.padWeight * (ds.padding - y);
      } else if (y > (height - ds.padding)) {
        forceY -= ff.padWeight * (y - (height - ds.padding));
      }
    }
    
    // Attract to center of mass
    forceX -= ff.cmWeight * (cmx - (width/2.0));
    forceY -= ff.cmWeight * (cmy - (height/2.0));
    
    
    // Straighten
    if (fiveP != null && fiveP.fiveP != null) {
      double dx = fiveP.x - fiveP.fiveP.x; // Vector from the base 2 positions upstream to base upstream
      double dy = fiveP.y - fiveP.fiveP.y;
      //double dz = fiveP.z - fiveP.fiveP.z;
      double len = normFunction(dx, dy, 0);
      dx /= len; // normalize
      dy /= len;
      //dz /= len;
      
      double dx2 = x - fiveP.x; // Vector from upstream base to this base
      double dy2 = y - fiveP.y;
      //double dz2 = z - fiveP.z;
      double len2 = normFunction(dx2, dy2, 0);
      dx2 /= len2;
      dy2 /= len2;
      //dz2 /= len2;
      
      PVector p55 = new PVector((float)dx, (float)dy, 0);
      PVector p5 = new PVector((float)dx2, (float)dy2, 0);
      double angTerm = angleBetweenFunct(p5, p55); // bending angle: ideally 0

      double px = fiveP.x + len2 * dx; // Continuation of line made by 5'5' to 5': Ideal point of this base
      double py = fiveP.y + len2 * dy;
      //double pz = fiveP.z + len2 * dz;
      double ddx = px - x; // Vector from this base's position to ideal position. Pull is proportional to deviation from ideal
      double ddy = py - y;
      //double ddz = pz - z;
      
      double straightFactor = ff.straightWeight * angTerm * angTerm; // Proportional to square of angle
      
      forceX += ddx * straightFactor;
      forceY += ddy * straightFactor;
      //forceZ += ddz * straightFactor;
      
      /**
      // debugging
      double sForce = Math.sqrt(Math.pow(ddx * straightFactor, 2) + Math.pow(ddy * straightFactor, 2));
      sForceTot += sForce;
      if (sForce > sMax) { sMax = sForce; }
      ++sCount;
      */
    }
    
    // Straighten
    if (threeP != null && threeP.threeP != null) {
      double dx = threeP.x - threeP.threeP.x; // Vector from the base 2 positions upstream to base upstream
      double dy = threeP.y - threeP.threeP.y;
      //double dz = threeP.z - threeP.threeP.z;
      double len = normFunction(dx, dy, 0);
      dx /= len; // normalize
      dy /= len;
      //dz /= len;
      
      double dx2 = x - threeP.x; // Vector from upstream base to this base
      double dy2 = y - threeP.y;
      //double dz2 = z - fiveP.z;
      double len2 = normFunction(dx2, dy2, 0);
      dx2 /= len2;
      dy2 /= len2;
      //dz2 /= len2;
      
      PVector p55 = new PVector((float)dx, (float)dy, 0);
      PVector p5 = new PVector((float)dx2, (float)dy2, 0);
      double angTerm = angleBetweenFunct(p5, p55); // bending angle: ideally 0

      double px = threeP.x + len2 * dx; // Continuation of line made by 5'5' to 5': Ideal point of this base
      double py = threeP.y + len2 * dy;
      //double pz = threeP.z + len2 * dz;
      double ddx = px - x; // Vector from this base's position to ideal position. Pull is proportional to deviation from ideal
      double ddy = py - y;
      //double ddz = pz - z;
      
      double straightFactor = ff.straightWeight * angTerm * angTerm; // Proportional to square of angle
      
      forceX += ddx * straightFactor;
      forceY += ddy * straightFactor;
      //forceZ += ddz * straightFactor;
      
      /**
      // debugging
      double sForce = Math.sqrt(Math.pow(ddx * straightFactor, 2) + Math.pow(ddy * straightFactor, 2));
      sForceTot += sForce;
      if (sForce > sMax) { sMax = sForce; }
      ++sCount;
      */
    }
    
    /*
    // Temperature slider adds random shake
    if (ff.jitterWeight > 0) {
      PVector randVec = PVector.random3D(); // random direction
      PVector randVec2 = PVector.random2D(); // random magnitude
      double jitter = ff.jitterWeight * randVec2.x * randVec2.y;
      forceX += jitter * randVec.x;
      forceY += jitter * randVec.y;
      forceZ += jitter * randVec.z;
    }
    */
    
    double ax = forceX / ff.mass;
    vx = (ff.damping * vx) + ax; // Old: vx = ff.damping * (vx + ax);
    
    double ay = forceY / ff.mass;
    vy = (ff.damping * vy) + ay;

    double vel = normFunction(vx, vy, 0);
    
    // debugging
    //double vtemp = Math.sqrt(Math.pow(vx, 2) + Math.pow(vy, 2));
    //if (vtemp > vMax) { vMax = vtemp; }
    
    // Limit velocities
    if (vel > ff.velMax) {
      double velReduce = ff.velMax / vel;
      vx *= velReduce;
      vy *= velReduce;
    }
    
    //vTot += Math.sqrt(Math.pow(vx, 2) + Math.pow(vy, 2)); // debugging

    // Apply velocity to residue
    x += vx;
    y += vy;
    //z += vz;
    
    forceX = 0;
    forceY = 0;
    //forceZ = 0;
  }

  /* UNUSED
  double[] directionTo(Spring2D other) {
    double[] result = new double[3];
    result[0] = other.x - x;
    result[1] = other.y - y;
    result[2] = other.z - z;
    double d = distance(other);
    if (d > 0.0) {
      result[0] /= d;
      result[1] /= d;
      result[2] /= d;
    }
    return result;
  }
  */

  // Draw bond and backbone lines
  void displayLines(boolean[] drawn) { // TODO: MOVE somewhere else?
    // Backbone
    drawer.Stroke(ds.bbCol);
    drawer.StrokeWeight(ds.bbWeight / ds.zoomFactor);
    if (threeP != null && !drawn[threeP.id]) { // If threeP has not been drawn yet.
      //       line((float)x, (float)y, (float)threeP.x, (float)threeP.y);
      drawer.Line((float)x, (float)y, (float)threeP.x, (float)threeP.y);
    }
    if (fiveP != null && !drawn[fiveP.id]) { // If threeP has not been drawn yet.
      //       line((float)x, (float)y, (float)fiveP.x, (float)fiveP.y);
      drawer.Line((float)x, (float)y, (float)fiveP.x, (float)fiveP.y);
    }
    
    // Base pairing
    if (pairedWForce != null && "cWW".equals(bonds[pairedWForce.id]) && !drawn[pairedWForce.id]) { // If pairedW has not been drawn yet. // ADDing: Check if its cWW!
      drawer.Stroke(ds.bpCol);
      drawer.StrokeWeight(ds.bpWeight / ds.zoomFactor);
      //       line((float)x, (float)y, (float)pairedWForce.x, (float)pairedWForce.y);
      drawer.Line((float)x, (float)y, (float)pairedWForce.x, (float)pairedWForce.y);
    }
  }

  // Draw base
  // TODO: contains code that should be part of DisplaySettings such as color palettes etc
  void displayBody(float displayDiam, float offsetX, float offsetY) {
    drawer.Fill(#AAAAAA); // must be called before NoStroke
    if (ds.outlineMode) {
      drawer.Stroke(#000000);
      drawer.StrokeWeight( (sqrt((float)(ff.radius * ds.radiusShowMul)) / 3.0) / ds.zoomFactor );
    } else {
      drawer.NoStroke();
    }
    
    // Light colors
    if (ds.colorMode == DisplaySettings.COLOR_MODE_LIGHT) {
      color c = color(#AAAAAA); // grey
      int id = (strandId) % 9;
      switch (id) {
      case 0: 
        c = color(#BFF9FF); // light blue 
        break;
      case 1:
        c = color(#FFE5B2);  // light orange 
        break;
      case 2: 
        c = color(#B2FFB9); // green 
        break;
      case 3: 
        c = color(#FFCBF6); // pink 
        break;
      case 4: 
        c = color(#FFFF98); // gold 
        break;
      case 5: 
        c = color(#FF9A98); // light red 
        break;
      case 6: 
        c = color(#DA8BFF);  // purple
        break;
      case 7: 
        c = color(#FFFFFF); // white
        break;
      case 8: 
        c = color(#7EB9FF); // Blue
        break;
      default: {
        if (debugPrints) { println("Internal error in color pallete"); }
        c = color(#FFFFFF);
        break;
        }
      }
      //       fill(c);
      drawer.Fill(c);

    // Dark pastels
    } else if (ds.colorMode == DisplaySettings.COLOR_MODE_FEW) {
      // RGB values from Show Me the Numbers from Stephen Few:
      color c = color(#AAAAAA); // grey
      int id = (strandId) % 8;
      switch (id) {
      case 0: 
        c = color(#65B2E8); // light blue #65B2E8
        break;
      case 1:
        c = color(#FAA43A);  // light orange #FAA43A
        break;
      case 2: 
        c = color(#60BD68); // green #60BD68
        break;
      case 3: 
        c = color(#F17CB0); // salmon #F17CB0
        break;
      case 4: 
        c = color(#F2DE44); // gold #F2DE44
        break;
      case 5: 
        c = color(#F15854); // light red #F15854
        break;
      case 6: 
        c = color(#B276B2);  // purple #B276B2
        break;
      case 7: 
        c = color(#C0C0C0); // light grey #C0C0C0
        break;
      default: {
        if (debugPrints) { println("Internal error in color pallete"); }
        c = color(#FFFFFF);
        break;
        }
      }
      //       fill(c);
      drawer.Fill(c);

    // Color by nucleotide type
    } else if (ds.colorMode == DisplaySettings.COLOR_MODE_RESIDUE) {
      String temp = residue.toUpperCase();
      if ("G".equals(temp)) {
        // fill(#64E4FF); //  #5DA5DA
        drawer.Fill(#64E4FF); //  #5DA5DA
      } else if ("A".equals(temp)) {
        // fill(#FF7471); // #60BD68
        drawer.Fill(#FF7471); // #60BD68
      } else if ("U".equals(temp)) {
        // fill(#64FF72); //  #98FFAA
        drawer.Fill(#64FF72); //  #98FFAA
      } else if ("T".equals(temp)) {
        // fill(#17FF8D);
        drawer.Fill(#17FF8D);
      } else if ("C".equals(temp)) {
        //  fill(#FFF68B); // #F17CB0
        drawer.Fill(#FFF68B); // #F17CB0
      } else {
        // fill(#FFFFFF);
        drawer.Fill(#FFFFFF);
      }
    
    // Rainbow, gap at dark purple
    } else if (ds.colorMode == DisplaySettings.COLOR_MODE_RAINBOW) {
      double frac = ( ds.colorMax * (double)(id) / ((double)currentState.sim.lenTot + (ds.colorMax/8.0)));
      if (frac > ((ds.colorMax*5.0) / 8.0)) {
        frac += ds.colorMax / 8.0;
      }
      // fill(color((int)frac, ds.colorMax, ds.colorMax));
      drawer.Fill(color((int)frac, ds.colorMax, ds.colorMax));
        
    // Bright neons
    } else if (ds.colorMode == DisplaySettings.COLOR_MODE_BYG) {
      color c = color(#AAAAAA); // grey
      int id = (strandId) % 9;
      switch (id) {
      case 0: 
        c = color(#00FFFF); // cyan
        break;
      case 1:
        c = color(#FF8400); // orange
        break;
      case 2:
        c = color(#00FF00); // green
        break;
      case 3:
        c = color(#FE00FF); // magenta
        break;
      case 4:
        c = color(#FFF63E); // yellow
        break;
      case 5: 
        c = color(#FF0000); // red
        break;
      case 6:
        c = color(#0092FF); // blue
        break;
      case 7:
        c = color(#D400FF);  // purple
        break;
      case 8:
        c = color(#FFFFFF); // white
        break;
      default: 
      { 
        if (debugPrints) { println("Internal error in color pallete"); }
        c = color(#FFFFFF);
        break;
      }
    }
    // fill(c);
    drawer.Fill(c);
  
    } else if (ds.colorMode == DisplaySettings.COLOR_MODE_GREY) {
      // colorMode(RGB);
      color c = color(#DBDBDB);
      // fill(c);
      drawer.Fill(c);
      
    } else if (ds.colorMode == DisplaySettings.COLOR_MODE_WHITE) {
      drawer.Fill(#FFFFFF);
      if (ds.outlineMode) {
        drawer.Stroke(#000000);
      } else {
        drawer.NoStroke();
      }	
      float saveWeight = drawer.strokeWidth;
      // drawer.StrokeWeight( (sqrt((float)(ff.radius * ds.radiusShowMul)) / 3.0) / ds.zoomFactor );
      // fill(#FFFFFF);
      drawer.StrokeWeight(saveWeight / ds.zoomFactor);
    // TODO: add: hairpin loop, internal loop, bulge, psuedoknot
    } else if (ds.colorMode == DisplaySettings.COLOR_MODE_STRUCTURE) {
      color c = #BFF9FF; // Default light blue
      // Helix
      if (pairedWForce != null) {
        c = #FFF63E; // Yellow
      }
      // fill(c);
      drawer.Fill(c);
    
    // User input
    } else if (ds.colorMode == DisplaySettings.COLOR_MODE_CUSTOM) {
      // fill(baseColor);
      drawer.Fill(baseColor);
    } else if (ds.colorMode == DisplaySettings.COLOR_MODE_RNADNA) {
      color c = color(#AAAAAA); // grey
      if (residue.length() > 0){
       if (residue.toUpperCase().equals(residue)) { // upper case: RNA
         c = color(#FF8400); // orange
       } else { // lower case: DNA
         c = color(#00FFFF); // cyan
       }
      }
      drawer.Fill(c);
    }
    
    // If it is selected, highlight and enlarge
    if (attachedIds.containsKey(id)) {
      displayDiam *= ds.highlightMul;
      drawer.Stroke(ds.highlightCol);
    }
    
    // Draw residue
    drawer.Ellipse((float)x, (float)y, displayDiam, displayDiam);
       //  ellipse((float)x, (float)y, displayDiam, displayDiam);

    // Draw letter
    if (ds.textMode) {
//             fill(ds.textCol);
      drawer.Fill(ds.textCol);
//             text(residue, (float)x + offsetX, (float)y + offsetY);
      drawer.Text(residue, (float)x + offsetX, (float)y + offsetY);
    }
  }
}