/*
 * Copyright (c) Ronald A Chernich 2004. The copyright holder grants the right
 * to use and copy this work provided this copyright message is retained.
 *
 * This program was written to be compatible with Java 1.1.x for maximum
 * compatibility when executed within a browser based environment.
 */
import java.text.*;


/**
 * Calculates properties of a given harmonic cam shape for a four-stroke,
 * poppet valve engine.
 * The original concept for this program appeared in SIC magazine. See:
 *<p>
 * Jenkins, Rodrick: "Cam Design for Four-stroke cycle Engines other than
 * Radials", Strictly Internal Combustion Magazine, ISSN 10446567,
 * Postal Instant Press, Kent WA, Volume 3, Issue #18, page 3:
 *<p>
 * The methods calcNoseRad() and calcTables() are based on this article,
 * but caution: this dog has fleas; Jenkins got the lift right for a
 * "mushroom" follower, but the approach and results for deriving velocity
 * and acceleration are highly questionable.
 *<p>
 * Currently (v1.1) this method is replaced by a version based on equations
 * in an engineering text and at least produces the same results as that text
 * which the original would not.
 */

public class CamCalc
{
  public static final String VERSION = "1.2";

  private static final double G_IMP = 32.2;
  private static final double G_MET = 9.8;

  private double m_camAngle;
  private double m_valveLift;
  private double m_flankRad;
  private double m_baseRad;
  private double m_noseRad;
  private double m_lMax;
  private int m_rpm;

  private double m_noseLiftStart;

  private double [] m_lift;
  private double [] m_xDisp;
  private double [] m_velocity;
  private double [] m_accel;

  private int m_inc = 2;
  private boolean m_isMetric = false;

  /**
   * Prepare table of results for params
   * @param ca Cam Angle (degrees)
   * @param lift Cam lift (metric or imperial)
   * @param lift Cam flank raduis (metric or imperial)
   * @param lift Cam base circle radius (metric or imperial)
   * @param rpm Engine speed (Revolutions per minute)
   * @exception Exception if a gross param error detected.
   */
  public CamCalc (double ca, double lift, double fr, double br, int rpm)
    throws Exception
  {
    m_camAngle = ca;
    m_valveLift = lift;
    m_flankRad = fr;
    m_baseRad = br;
    m_rpm = rpm;

    if (m_camAngle > 180.0) {
      throw new Exception("Cam angle must be <= 180 degrees");
    }
    if (m_camAngle == 180.0) {
      m_camAngle = 179.99;
    }

    calcNoseRad();
  }

  /**
   * Sets the increment used when creating the displayable output table.
   * @param inc Degrees per table row increment
   */
  public void setIncrement (int val)
  {
    if (val > 0) {
      m_inc = val;
    }
  }

  /**
   * Sets the rounding mode, conversions and G figures. Note that if metric
   * units are required, thie method must be called before any output data
   * is retrieved.
   * @param isMetric True for metric formatting and calculations
   * @exception RuntimeException if trying to set after horse bolted.
   */
  public void setMetric (boolean isMetric)
    throws RuntimeException
  {
    if (m_lift != null) {
      throw new RuntimeException("Cannot change units after fetching data");
    }
    m_isMetric = isMetric;
  }

  /**
   * Formats the input data (together with calculated nose radius and
   * minimum tappet width).
   */
  public String getInDataSummary()
  {
    if (m_lift == null) {
      calcTables2();
      //calcTables();
    }
    return m_isMetric ? inMetric() : inImp();
  }

  /**
   * Creates a summary table of characteristics from calculated data
   */
  public String getOutDataSummary()
  {
    if (m_lift == null) {
      calcTables2();
      //calcTables();
    }
    return m_isMetric ? outMetric() : outImp();
  }

  /**
   * @return Table generated for metric values
   */
  private String outMetric ()
  {
    DecimalFormat decForm = new DecimalFormat();
    StringBuffer sb = new StringBuffer();
    int idx = 0;
    sb.append("Angle\tLift\tVel\tAccel\tAccel\n");
    sb.append("Degrees\tmm\tm/s\tm/s/s\tg\n");

    int angle = 360 - (int)(m_camAngle / 2);
    while (idx < m_camAngle) {
      sb.append(angle).append('\t');
      decForm.applyPattern("0.00");
      sb.append(decForm.format(m_lift[idx]));
      sb.append('\t').append(decForm.format(m_velocity[idx]));
      decForm.applyPattern("###0");
      sb.append('\t').append(decForm.format(m_accel[idx]));
      sb.append('\t').append(decForm.format(m_accel[idx] / G_MET));
      sb.append('\n');
      idx += m_inc;
      angle += m_inc;
      angle %= 360;
    }
    return sb.toString();
  }

  /**
   * @return Table generated for Imperial values
   */
  private String outImp ()
  {
    DecimalFormat decForm = new DecimalFormat();
    StringBuffer sb = new StringBuffer();
    int idx = 0;
    sb.append("Angle\tLift\tVel\tAccel\tAccel\n");
    sb.append("Degrees\tthous\tft/s\tft/s/s\tg\n");

    int angle = 360 - (int)(m_camAngle / 2);
    while (idx <= m_camAngle) {
      sb.append(angle).append('\t');
      decForm.applyPattern("##0");
      sb.append(decForm.format(m_lift[idx] * 1000));
      decForm.applyPattern("0.00");
      sb.append('\t').append(decForm.format(m_velocity[idx]));
      decForm.applyPattern("###0");
      sb.append('\t').append(decForm.format(m_accel[idx]));
      sb.append('\t').append(decForm.format(m_accel[idx] / G_IMP));
      sb.append('\n');
      idx += m_inc;
      angle += m_inc;
      angle %= 360;
    }
    return sb.toString();
  }

  /**
   * @return Parameter for Imperial values
   */
  public String inImp ()
  {
    DecimalFormat decForm = new DecimalFormat("###.###");
    int [] gHolder = new int[2];
    getMinMaxG(gHolder);
    StringBuffer sb = new StringBuffer();
    sb.append("        Cam Open Angle (degrees): ").append(m_camAngle);
    sb.append("\n              Valve Lift (thous): ");
    sb.append((int)((m_valveLift * 1000) + 0.5));
    sb.append("\n           Flank Radius (inches): ").append(m_flankRad);
    sb.append("\n            Base Radius (inches): ").append(m_baseRad);
    sb.append("\n              Engine Speed (RPM): ").append(m_rpm);
    sb.append("\n \n         Nose Radius (inches): ");
    sb.append(decForm.format(m_noseRad));
    sb.append("\nTang follower width min (inches): ");
    sb.append(decForm.format(m_lMax * 2));
    sb.append("\n            Max acceleration (g): ").append(gHolder[1]);
    sb.append("\n            Min acceleration (g): ").append(gHolder[0]);
    sb.append("\n        Lift Area (inch degrees): ");
    sb.append(decForm.format(getLiftArea()));
    decForm.applyPattern("##0.0#");
    sb.append("\n Nose Transition delay (degrees): ");
    sb.append(decForm.format(m_noseLiftStart));
    sb.append('\n');
    return sb.toString();
  }

  /**
   * @return Parameter for metric values
   */
  public String inMetric ()
  {
    DecimalFormat decForm = new DecimalFormat();
    int [] gHolder = new int[2];
    getMinMaxG(gHolder);
    StringBuffer sb = new StringBuffer();
    sb.append("\n       Cam Open Angle (degrees): ").append(m_camAngle);
    sb.append("\n                Valve Lift (mm): ");
    sb.append(decForm.format(m_valveLift));
    sb.append("\n              Flank Radius (mm): ").append(m_flankRad);
    sb.append("\n               Base Radius (mm): ").append(m_baseRad);
    sb.append("\n             Engine Speed (RPM): ").append(m_rpm);
    sb.append("\n \n               Nose Radius (mm): ");
    sb.append(decForm.format(m_noseRad));
    sb.append("\n   Tang follower width min (mm): ");
    sb.append(decForm.format(m_lMax * 2));
    sb.append("\n           Max acceleration (g): ").append(gHolder[1]);
    sb.append("\n           Min acceleration (g): ").append(gHolder[0]);
    sb.append("\n         Lift Area (mm degrees): ");
    sb.append(decForm.format(getLiftArea()));
    decForm.applyPattern("##0.0#");
    sb.append("\nNose Transition delay (degrees): ");
    sb.append(decForm.format(m_noseLiftStart));
    sb.append('\n');
    return sb.toString();
  }

  /**
   * Calculates the Cam nose radius given base circle, lift and flanks
   * (assumes rising and falling flanks are the same)
   */
  private void calcNoseRad ()
    throws Exception
  {
    double p3 = Math.cos(Math.toRadians(m_camAngle / 2.0));
    double ar1 = -(2.0 * m_baseRad * m_flankRad)
                 + Math.pow((m_baseRad + m_valveLift), 2.0)
                 + Math.pow(m_baseRad, 2.0)
                 + (2 * (m_baseRad + m_valveLift) *
                        (m_flankRad - m_baseRad) * p3);
    double br1 = ((m_baseRad + m_valveLift - m_flankRad) -
                  ((m_baseRad - m_flankRad) * p3)) * 2.0;
    m_noseRad = ar1 / br1;
    if (m_noseRad < 0) {
      String msg = "Flank radius too small to reach stated lift" +
                   "\n(nose radius = " + m_noseRad + ")";
      throw new Exception(msg);
    }
  }

  /**
   * Builds tables of data for each degree of cam rotation over the cam
   * angle using the formulae presented in:
   *   Bevan, Thomas (ed.): "The Theory of Machines", p288 et seq. 
   */
  private void calcTables2 ()
  {
    m_lift = new double[180];
    m_accel = new double[180];
    m_velocity = new double[180];

    // Assume the cam will be rotating at 1/2 engine speed
    double w = 2.0 * Math.PI * ((m_rpm / 2) / 60);
    double alpha = m_camAngle / 2;
    double op = m_flankRad - m_baseRad;
    double uop = op / (m_isMetric ? 1000 : 12);  // meters or feet
    double oq = m_baseRad + m_valveLift - m_noseRad;
    double uoq = oq / (m_isMetric ? 1000 : 12);  // meters or feet

    // Find transition angles to and from nose radius lift from flank
    // radius lift relative to the commencement of cam lift.
    double pq2 = Math.pow(op, 2) + Math.pow(oq, 2) + 
                (2 * op * oq * Math.cos(Math.toRadians(alpha)));
    double pq = Math.pow(pq2, 0.5);
    double sinPhi = (oq / pq) * Math.sin(Math.toRadians(alpha));
    double t1 = Math.toDegrees(Math.asin(sinPhi));
    double t2 = m_camAngle - t1;

    m_noseLiftStart = t1;

    for (int idx = 0; idx <= m_camAngle; idx++) {
      double theta = idx;
      if (theta <= t1) {
        // lift on the rising flank
        double delta = Math.toRadians(theta);
        m_lift[idx] = op * ( 1 - Math.cos(delta));
        m_velocity[idx] = w * uop * Math.sin(delta);
        m_accel[idx] = Math.pow(w, 2) * uop * Math.cos(delta);
        // find offset of tangential contact point
        double m = op * Math.sin(Math.toRadians(theta - t1));
        m_lMax = Math.max(m_lMax, m);
      }
      else if ((theta > t1) && (theta < t2)) {
        // lift on the nose
        double delta = Math.toRadians(alpha - theta);
        m_lift[idx] = m_noseRad - m_baseRad + (oq * Math.cos(delta));
        m_velocity[idx] = w * uoq * Math.sin(delta);
        m_accel[idx] = -Math.pow(w, 2) * uoq * Math.cos(delta);
        // find offset of tangential contact point
        double m = oq * Math.sin(delta);
        m_lMax = Math.max(m_lMax, m);
      }
      else if (theta >= t2) {
        // lift on the falling flank
        double delta = Math.toRadians(m_camAngle - theta);
        m_lift[idx] = op * ( 1 - Math.cos(delta));
        m_velocity[idx] = -w * uop * Math.sin(delta);
        m_accel[idx] = Math.pow(w, 2) * uop * Math.cos(delta);
        // no need to find offset of tangential contact point
      }
      else {
        // no lift
        m_lift[idx] = 0.0;
      }
    }
  }

  /**
   * Builds tables of data for each degree of cam rotation over the
   * "action" sector. Variable names are the same as used in Jenkins'
   * BASIC program for debug purposes (uggh).
   */
  private void calcTables ()
  {
    m_lift = new double[180];
    m_accel = new double[180];
    m_velocity = new double[180];

    double q = Math.PI/180.0;
    double p3 = m_camAngle / 2.0;
    double r0 = m_baseRad;
    double r1 = m_noseRad;
    double rr = m_flankRad;
    double rf = rr;

    // Calculate the transition point:
    //  p1 Start of flank rise rel max lift (ie zero and useless)
    //  p2 From rising flank to nose radius
    //  p3 Half cam action angle (look at the way he derives it!)
    //  p4 Nose radius to falling flank
    //  p5 Falling flank back to base circle
    // p1 is a bit useless. Aside: this is typical of the confused thinking
    // imposed by (uggh) BASIC, especially early BASIC interpreters. The
    // unhelpful identifiers cause  the developer to loos track of what's
    // going on and entropy (in the form of chaos) invariably follows...
    double h = m_baseRad + m_valveLift - m_noseRad;
    double p1 = 0;
    double ap2 = (Math.pow((rr-r0), 2) + Math.pow((rr-r1), 2) - Math.pow(h, 2))/
                 (2 * (rr-r0) * (rr-r1));
    double p2 = (Math.PI/2 - Math.atan(ap2/Math.sqrt(-ap2*ap2+1)))/q;
    double ap4 = (Math.pow(h, 2) + Math.pow((rf-r1), 2) - Math.pow((rf-r0), 2))/
                 (2 * h * (rf-r1));
    double bp4 = Math.PI/2 - Math.atan(ap4/Math.sqrt(-ap4*ap4+1));
    double p4 = p3 + (bp4 / q);
    double p5 = 2 * p3;
//System.err.println(p1+" "+p2+" "+p3+" "+p4+" "+p5);

    m_noseLiftStart = p2;

    // For 180 degrees of arc, apply appropriate calculation to find
    // lift applied to a tangential follower. Retain the maximum deviation
    // from the reference to where the tangential follower contacts the cam.
    // This will be the minimum width/2 for a fully tangential cam follower.
    for (int idx = 0; idx <= m_camAngle; idx++) {
      double p = idx;
      double r = 0;
      if ((p >= p1) && (p <= p2)) {
        r = rr - (rr - r0) * Math.cos(Math.toRadians(p-p1));
        double m = (rr - r0) * Math.sin(Math.toRadians(p-p1));
        m_lMax = Math.max(m_lMax, m);
      }
      else if ((p > p2) && (p <= p4)) {
        r = r1 + h * Math.cos(Math.toRadians(p3-p));
        double m = h * Math.sin(Math.toRadians(p3-p));
        m_lMax = Math.max(m_lMax, m);
      }
      else if ((p >= p4) && (p <= p5)) {
        r = rf - (rf - r0) * Math.cos(Math.toRadians(p5-p));
        double m = (rf - r0) * Math.sin(Math.toRadians(p5-p));
        m_lMax = Math.max(m_lMax, m);
      }
      else { // ((p > p5) && (p < p1))
        r = r0;
      }

      m_lift[idx] = r - r0;
    }

    // FIXUP: I don't get it. Using the "averaging" method, "t"
    // *should* be time for one rev, not time for 5 revs!
    // Plus, by omitting parenthesis, Jenkins leaves the result in the
    // hands of the BASIC interpreter's commutability.  Dangerous.
    // Alltogether, this dog don't hunt...
    int rpm = m_rpm / 2;
    double t = 1.0 / rpm / 6.0;
    int lim = (int)m_camAngle - 2;
    for (int idx = 1; idx < lim; idx++) {
      double dv = ((m_lift[idx+1] - m_lift[idx-1])/2.0)/t;
      m_velocity[idx] = dv / (m_isMetric ? 1000.0 : 12.0);
    }

    for (int idx = 2; idx < lim-1; idx++) {
      double dv = ((m_velocity[idx+1] - m_velocity[idx-1])/2.0)/t;
      m_accel[idx] = dv;
      // FIXUP: orig prog placed the G figure in the acceleration matrix
      // and applied G again at table print time!  Go figure...
      //m_accel[idx] = dv / (m_isMetric ? 9.8 : 32.2);
    }
  }

  /**
   * Find min/max acceleration (G)
   */
  private void getMinMaxG (int [] holder)
  {
    double vMin = 0;
    double vMax = 0;
    int lim = (int)m_camAngle;
    for (int idx = 0; idx < lim; idx++) {
      vMin = Math.min(vMin, m_accel[idx]);
      vMax = Math.max(vMax, m_accel[idx]);
    }
    holder[0] = (int)((vMin / (m_isMetric ? G_MET : G_IMP)) + .5);
    holder[1] = (int)((vMax / (m_isMetric ? G_MET : G_IMP)) + .5);
  }

  /**
   * @return Approximation of the area under cam lift curve (unit * degrees)
   */
  private double getLiftArea ()
  {
    double res = 0;
    double dy = 0;
    int lim = (int)m_camAngle - 1;
    for (int idx = 0; idx < lim; idx++) {
      double y1 = m_lift[idx];
      double y2 = m_lift[idx + 1];
      res += (dy + Math.abs((y2 - y1)/2));
      dy = y2;
    }
    return res;
  }

  /**
   * Command line driver for Cam Calculator class.  See "Usage" message
   * for argument list. The cmd line driver does not yet do metric.
   */
  public static void main (String [] args)
  {
    if (args.length < 5) {
      System.err.println("Usage: java CamCalc <N> <L> <RR> <R0> <RT> [-ax -m]");
      System.err.println("  where N = Cam angle (degrees)");
      System.err.println("        L = Valve Lift (inches)");
      System.err.println("       RR = Flank Radius (inches)");
      System.err.println("       R0 = Base Radius (inches)");
      System.err.println("       RT = Engine Speed (RPM)");
      System.err.println("      -ax = Table data every 'x' degrees");
      System.err.println("      -m  = Metric Mode [default is Imperial]");
      System.exit(1);
    }

    CamCalc camCalc = null;
    try {
      double ca = Double.parseDouble(args[0]);
      double lf = Double.parseDouble(args[1]);
      double fr = Double.parseDouble(args[2]);
      double br = Double.parseDouble(args[3]);
      int rpm = Integer.parseInt(args[4]);
      camCalc = new CamCalc(ca, lf, fr, br, rpm);
    }
    catch (Exception ex) {
      System.err.println(ex.getMessage());
      System.exit(1);
    }

    int idx = 5;
    while (idx < args.length) {
      if (args[idx].startsWith("-a")) {
        camCalc.setIncrement(Integer.parseInt(args[idx].substring(2)));
      }
      else if (args[idx].equals("-m")) {
        camCalc.setMetric(true);
      }
      else {
        System.err.println("Invalid argument '" + args[idx] + "'");
        System.exit(1);
      }
      ++idx;
    }
    System.out.println(camCalc.getInDataSummary());
    System.out.println(camCalc.getOutDataSummary());
    System.exit(0);
  }

} // endof CamCalc