/* Rocket.Java
 * Bill True, December 1998
 * Based on the program water.c by Michael Holder.
 */
import java.util.*;
import java.awt.*;
import java.applet.*;
  
// *** This class holds all the physical constants ***
class Phy extends Object {
  static final double ACCEL_GRAVITY =  32.17;    /* ft/s^2 */ 
  static final double AMB_PRESSURE = 14.7;       /* PSI  */
  static final double RHO_WATER = 0.0011229586;  /* slugs/cu in. */
  static final double RHO_AIR = 0.0023783028;    /* slugs/cu. in. */
  static final double GAMMA_AIR = 1.4;
}

// *** This class is the applet ***
public class Rocket extends Applet implements Runnable {
  static final int PREF_WIDTH = 650;
  static final int PREF_HEIGHT = 500;
  Thread simulationThread = null;
  Motor theMotor;

  // sim params
  double time, delta_t, print_interval, graph_time,
         altitude, velocity, acceleration;

  // motor params
  double eng_diameter, nozzle_diameter, discharge_coef, 
         thrust_coef, eng_length, load_fraction,
         gage_pressure, dry_weight, drag_coef;

  // GUI controls
  DoubleReader dtControl, piControl, gtControl, diaControl, 
               nzlControl, dcControl, tcControl, mlControl, 
               lfControl, pControl, wtControl, drgControl;

  Button       goButton;
  GraphCanvas  theGraph;
  Panel        theControlPanel;
  Panel        theKeyPanel;

  Label        burnLbl, impulseLbl, altLbl, timeLbl;
  Panel        theResultPanel;

public void init() {
  String paramStr;

  // initialize simulation parameters from web page
  try { paramStr = getParameter("TIMESTEP");
        delta_t = Double.valueOf(paramStr).doubleValue();
  } catch (Exception e) { delta_t = 0.001; }
   
  try { paramStr = getParameter("PRINT_INTERVAL");
        print_interval = Double.valueOf(paramStr).doubleValue();
  } catch (Exception e) { print_interval = 0.1; }

  try { paramStr = getParameter("GRAPH_TIME");
        graph_time = Double.valueOf(paramStr).doubleValue();
  } catch (Exception e) { graph_time = 1.0; }
  
  try { paramStr = getParameter("MOTOR_DIA");
        eng_diameter = Double.valueOf(paramStr).doubleValue();
  } catch (Exception e) { eng_diameter = 1.0; }

  try { paramStr = getParameter("NOZZLE_DIA");
        nozzle_diameter = Double.valueOf(paramStr).doubleValue();
  } catch (Exception e) { nozzle_diameter = 0.5; }

  try { paramStr = getParameter("DISCHARGE_COEF");
        discharge_coef = Double.valueOf(paramStr).doubleValue();
  } catch (Exception e) { discharge_coef = 1.0; }

  try { paramStr = getParameter("THRUST_COEF");
        thrust_coef = Double.valueOf(paramStr).doubleValue();
  } catch (Exception e) { thrust_coef = 1.0; }

  try { paramStr = getParameter("MOTOR_LEN");
        eng_length = Double.valueOf(paramStr).doubleValue();
  } catch (Exception e) { eng_length = 10.0; }

  try { paramStr = getParameter("LOAD_FRAC");
        load_fraction = Double.valueOf(paramStr).doubleValue();
  } catch (Exception e) { load_fraction = 0.5; }

  try { paramStr = getParameter("GAUGE_PRESSURE");
        gage_pressure = Double.valueOf(paramStr).doubleValue();
  } catch (Exception e) { gage_pressure = 100.0; }

  try { paramStr = getParameter("DRY_WEIGHT");
        dry_weight = Double.valueOf(paramStr).doubleValue();
  } catch (Exception e) { dry_weight = 10.0; }

  try { paramStr = getParameter("DRAG_COEF");
        drag_coef = Double.valueOf(paramStr).doubleValue();
  } catch (Exception e) { drag_coef = 0.0; }
  
  // Set window size and add all the components
  resize(PREF_WIDTH,PREF_HEIGHT);              
  theGraph = new GraphCanvas();
  goButton = new Button("Launch!");

  dtControl  = new DoubleReader("Delta time (sec.)",     0.0, 99999.0, delta_t);
  piControl  = new DoubleReader("Print Interval (sec.)", 0.0, 99999.0, print_interval);
  gtControl  = new DoubleReader("Graph Time (sec.)",     0.0, 99999.0, graph_time);
  diaControl = new DoubleReader("Motor Diameter (in.)",  0.0, 99999.0, eng_diameter);
  nzlControl = new DoubleReader("Nozzle Diameter (in.)", 0.0, 99999.0, nozzle_diameter);
  dcControl  = new DoubleReader("Discharge Coefficient", 0.0, 99999.0, discharge_coef);
  tcControl  = new DoubleReader("Thrust Coefficient",    0.0, 99999.0, thrust_coef);
  mlControl  = new DoubleReader("Motor Length (in.)",    0.0, 99999.0, eng_length);
  lfControl  = new DoubleReader("Load Fraction (water)", 0.0, 1.0,     load_fraction);
  pControl   = new DoubleReader("Pressure (gauge psi.)", 0.0, 99999.0, gage_pressure);
  wtControl  = new DoubleReader("Dry Weight (oz.)",      0.0, 99999.0, dry_weight);
  drgControl = new DoubleReader("Drag Coefficient",      0.0, 99999.0, drag_coef);
  theControlPanel = new Panel();
  theControlPanel.setLayout(new GridLayout(4,3)); 
  theControlPanel.add(dtControl);  
  theControlPanel.add(piControl);
  theControlPanel.add(gtControl);
  theControlPanel.add(diaControl);
  theControlPanel.add(nzlControl);
  theControlPanel.add(dcControl);
  theControlPanel.add(tcControl);
  theControlPanel.add(mlControl);
  theControlPanel.add(lfControl);
  theControlPanel.add(pControl);
  theControlPanel.add(wtControl);
  theControlPanel.add(drgControl);

  theKeyPanel = new Panel();
  theKeyPanel.setLayout(new GridLayout(6,1));
  Label tmpLabel;
  tmpLabel = new Label("_____ Pressure");
  tmpLabel.setForeground(Color.blue);
  theKeyPanel.add(tmpLabel);
  tmpLabel = new Label("_____ Mass");
  tmpLabel.setForeground(Color.green);
  theKeyPanel.add(tmpLabel);
  tmpLabel = new Label("_____ Acceleration");
  tmpLabel.setForeground(Color.orange);
  theKeyPanel.add(tmpLabel);
  tmpLabel = new Label("_____ Velocity");
  tmpLabel.setForeground(Color.red);
  theKeyPanel.add(tmpLabel);
  tmpLabel = new Label("_____ Altitude");
  tmpLabel.setForeground(Color.black);
  theKeyPanel.add(tmpLabel);
  theKeyPanel.add(goButton);

  burnLbl    = new Label("        ");
  impulseLbl = new Label("        ");
  altLbl     = new Label("        ");
  timeLbl    = new Label("        ");
  theResultPanel = new Panel();
  theResultPanel.setLayout(new GridLayout(8,1));
  tmpLabel = new Label("Burnout Time (sec.)");
  theResultPanel.add(tmpLabel);
  theResultPanel.add(burnLbl);
  tmpLabel = new Label("Total Impulse (lb-sec)");
  theResultPanel.add(tmpLabel);
  theResultPanel.add(impulseLbl);
  tmpLabel = new Label("Max Altitude (ft.)");
  theResultPanel.add(tmpLabel);
  theResultPanel.add(altLbl);
  tmpLabel = new Label("Total Flight Time (sec.)");
  theResultPanel.add(tmpLabel);
  theResultPanel.add(timeLbl);

  setLayout(new FlowLayout());
  add(theControlPanel);
  add(theResultPanel);
  add(theGraph);
  add(theKeyPanel);
}

public boolean action(Event e, Object arg) {
  System.out.println("action");

  if (e.target == goButton) {
    System.out.println("goButton pressed");
    simulationThread = new Thread(this);
    simulationThread.start();
    return true;
  } else
    return false;
}

public void run() {
  final int THRUST_PHASE   = 0;
  final int COAST_PHASE    = 1;
  final int DESCENT_PHASE  = 2;
  final int IMPACT_PHASE   = 3;

  Graphics g;
  String paramStr;

  System.out.println("run");

  if (simulationThread != null) {

    // get current values from all the controls  
    try { print_interval = piControl.getValue(); } 
    catch (Exception e) { ; }
    System.out.println("read print_interval: " + print_interval);
 
    try { delta_t = dtControl.getValue(); } 
    catch (Exception e) { ; }
    System.out.println("read delta_t: " + delta_t);

    try { graph_time = gtControl.getValue(); } 
    catch (Exception e) { ; }
    System.out.println("read graph_time: " + graph_time);

    try { eng_diameter = diaControl.getValue(); } 
    catch (Exception e) { ; }
    System.out.println("read eng_diameter: " + eng_diameter);
   
    try { nozzle_diameter = nzlControl.getValue(); } 
    catch (Exception e) { ; }
    System.out.println("read nozzle_diameter: " + nozzle_diameter);

    try { discharge_coef = dcControl.getValue(); } 
    catch (Exception e) { ; }
    System.out.println("read discharge_coef: " + discharge_coef);

    try { thrust_coef = tcControl.getValue(); } 
    catch (Exception e) { ; }
    System.out.println("read thrust_coef: " + thrust_coef);

    try { eng_length = mlControl.getValue(); } 
    catch (Exception e) { ; }
    System.out.println("read eng_length: " + eng_length);

    try { load_fraction = lfControl.getValue(); } 
    catch (Exception e) { ; }
    System.out.println("read load_fraction: " + load_fraction);

    try { gage_pressure = pControl.getValue(); } 
    catch (Exception e) { ; }
    System.out.println("read gage_pressure: " + gage_pressure);

    try { dry_weight = wtControl.getValue(); } 
    catch (Exception e) { ; }
    System.out.println("read dry_weight: " + dry_weight);

    try { drag_coef = drgControl.getValue(); } 
    catch (Exception e) { ; }
    System.out.println("read drag_coef: " + drag_coef);

    // These initial values come only from the web page,
    // they're loaded again for each run.

    try { paramStr = getParameter("T_ZERO");
          time = Double.valueOf(paramStr).doubleValue();
    } catch (Exception e) { time = 0.0; }
    System.out.println("time =" + Double.toString(time));  

    try { paramStr = getParameter("V_ZERO");
          velocity = Double.valueOf(paramStr).doubleValue();
    } catch (Exception e) { velocity = 0.0; }
    System.out.println("velocity =" + Double.toString(velocity));  

    try { paramStr = getParameter("H_ZERO");
          altitude = Double.valueOf(paramStr).doubleValue();
    } catch (Exception e) { altitude = 0.0; }
    System.out.println("altitude =" + Double.toString(altitude));  


    // create and initialize the motor
    dry_weight /= 16.0;  // convert oz to lbs
    theMotor = new Motor(drag_coef, thrust_coef, discharge_coef,
                         eng_diameter, eng_length, nozzle_diameter,
                         gage_pressure, load_fraction, dry_weight);
   
    // calculate burnout time
    double burnTime = theMotor.solve_t();
    showStatus("burnout time: " + burnTime);
    System.out.println("burnout time: " + burnTime);
    burnLbl.setText(formatDouble(burnTime, 8));
   
    // start the drawing
    g = theGraph.getGraphics();
    theGraph.drawAxis(g);

    // set up for each line (parameter) on the graph
    GraphLine pressLine = new GraphLine(theGraph, Color.blue,
                               time, time + graph_time,
                              -theMotor.getAve_pressure(), 
                               theMotor.getAve_pressure(),
                               time, 0.0);

    GraphLine massLine  = new GraphLine(theGraph, Color.green,
                               time, time + graph_time,
                              -theMotor.getMass(), theMotor.getMass(),
                               time, 0.0);

    double estMaxAccel = 1.25 * (theMotor.getThrust()/theMotor.getMass());
    GraphLine acclLine = new GraphLine(theGraph, Color.orange, 
                               time, time + graph_time,
                               -estMaxAccel, estMaxAccel,
                               time, 0.0);

    double estMaxVel = estMaxAccel * burnTime;
    GraphLine velLine = new GraphLine(theGraph, Color.red,
                               time, time + graph_time,
                              -estMaxVel, estMaxVel,
                               time, 0.0);

    double estMaxAlt = 2 * estMaxVel * burnTime;
    GraphLine altLine = new GraphLine(theGraph, Color.black,
                               time, time + graph_time,
                              -estMaxAlt, estMaxAlt,
                               time, 0.0);

    // set up time loop
    int print_iterations = (int) (print_interval/delta_t);
       
    // run the time loop
    double total_impulse = 0.0;
    double maxAlt = 0.0;
    long iteration_count = 0;
    int print_count = 0;

    System.out.println("time, ave_pressure, thrust, mass, drag, "
                     + "acceleration, velocity, altitude");

    int phase = THRUST_PHASE;
    while (phase 	!= IMPACT_PHASE) {
      theMotor.step_pressure(delta_t);
      acceleration = (theMotor.getThrust() - theMotor.getDrag(velocity))
                     /theMotor.getMass() - Phy.ACCEL_GRAVITY;
      velocity += acceleration * delta_t;
      altitude += velocity * delta_t;

      print_count--;
      if (print_count <= 0) {
        print_count = print_iterations;

        // print simulation state
        System.out.print(time);
        System.out.print(", " + theMotor.getAve_pressure());
        System.out.print(", " + theMotor.getThrust());
        System.out.print(", " + theMotor.getMass());
        System.out.print(", " + theMotor.getDrag(velocity));
        System.out.print(", " + acceleration);
        System.out.print(", " + velocity);
        System.out.println(", " + altitude);

        // plot simulation state
        pressLine.extend(time, theMotor.getAve_pressure());
        massLine.extend(time, theMotor.getMass());
        acclLine.extend(time, acceleration);
        velLine.extend(time, velocity);
        altLine.extend(time, altitude);
      }

      time += delta_t;
      iteration_count++;

      switch (phase) {
        case THRUST_PHASE : { 
          if (theMotor.getThrust() > 0.0001) {
            total_impulse += theMotor.getThrust() * delta_t;
          } 
          else {
            showStatus("total impulse: " + total_impulse);
            impulseLbl.setText(formatDouble(total_impulse, 8));
            phase = COAST_PHASE;
          }
          break;
        }

        case COAST_PHASE : { 
          if (velocity < 0.0) {
            maxAlt = altitude;
            showStatus("max altitude: " + maxAlt);
            altLbl.setText(formatDouble(maxAlt, 8));
            phase = DESCENT_PHASE;
          }
          break;
        }

        case DESCENT_PHASE : { 
          if (altitude <= 0.0) {
            showStatus("total flight time: " + time);
            timeLbl.setText(formatDouble(time, 8));
            phase = IMPACT_PHASE;
          }
          break;
       }

      } // end switch
           
    } // end while
    
    System.out.println("flight time: " + time);
    System.out.println("max altitude: " + maxAlt);
    System.out.println("total impulse: " + total_impulse); 
    System.out.println("simulation steps: " + iteration_count);
    simulationThread = null;
  }
}


public static String formatDouble(double num, int width){
  String hashStr = new String("###################################");
  String tmpStr;
  int    dotLoc;

  tmpStr = Double.toString(num);
  dotLoc = tmpStr.indexOf(".");
  if (dotLoc > width+1)
    tmpStr = hashStr.substring(0, width);
  else if (dotLoc == width-1)
    tmpStr = tmpStr.substring(0, width-1);
  else if (tmpStr.length() > width)
    tmpStr = tmpStr.substring(0, width);

  return tmpStr;
}

}

// *** This class implements the motor ***
class Motor extends Object {

  // simulation coefficients (these were global variables)
  double discharge_coef;      /* Adjusts nozzle velocity  */
  double thrust_coef;         /* Adjust thrust efficiency */
  double drag_coef;
  
  // motor state variables(these were in the motor_t c-structure)
  double pressure;        /* p.s.i. abs  */
  double total_volume;    /* cu. in.     */
  double air_volume;      /* cu. in.     */
  double nozzle_area;     /* sq. in.     */
  double dry_mass;        /* slugs       */
  double ave_pressure;    /* p.s.i. abs  */
  double ave_volume;      /* cu. in.     */

  // rocket parameters (was a variable in "main")
  double proj_area;

  Motor(double drg_coef, double thrst_coef, double disc_coef,
        double eng_dia, double eng_len, double noz_dia,
        double gag_pres, double ld_frac, double dry_wt) {

    System.out.println("Motor (constructor)");
    // set state variables and calculate derived quantities
    discharge_coef = disc_coef;
    thrust_coef = thrst_coef;
    drag_coef = drg_coef;
    
    pressure = gag_pres + Phy.AMB_PRESSURE;
    ave_pressure = pressure;
    total_volume = eng_len * Math.PI * eng_dia * eng_dia / 4.0;
    air_volume = total_volume * (1.0 - ld_frac);
    nozzle_area = Math.PI * noz_dia * noz_dia / 4.0;
    dry_mass = dry_wt / Phy.ACCEL_GRAVITY;

    proj_area = Math.PI * eng_dia * eng_dia / 576.0;  // eng_dia is in., area is sq ft.

    System.out.println("air volume: " + air_volume);
    System.out.println("total volume: " + total_volume);
    System.out.println("nozzle area: " + nozzle_area);
    System.out.println("dry mass: " + dry_mass);
  }

  double getMass() { 
    return dry_mass + (total_volume-ave_volume) * Phy.RHO_WATER; 
  }

  double getDrag(double vel) {
    double d;
    d = proj_area * drag_coef * Phy.RHO_AIR * vel * vel / 2.0;
    return ( (vel < 0) ? -d : d); 
  }

  double getAve_pressure() {
    return ave_pressure;
  }

  double getThrust() {
    return discharge_coef * thrust_coef * 2.0 * nozzle_area * 
           (ave_pressure - Phy.AMB_PRESSURE);
  }

  double ex_velocity(double pressure, double area) {
    return discharge_coef * area * 
           Math.sqrt(2*(pressure-Phy.AMB_PRESSURE)/Phy.RHO_WATER);
  }

  double adiabatic (double p1, double v1, double v2){
    return p1 * Math.pow(v1/v2,Phy.GAMMA_AIR);
  }


  double solve_t() {
    final int ITERATIONS = 100000;

    int i;
    double v;
    double sum;

    double delta_v;
    double sigma;
    double last_t;
    double next_t;

    delta_v=(total_volume-air_volume)/(ITERATIONS-1);
    sum=0.0;
    last_t=delta_v/ex_velocity(pressure,nozzle_area);

/* printf("dV=%10.4e  last_t=%10.4e\n",delta_v,last_t); */

    for (i=1; i<ITERATIONS; i++) {
	v=air_volume+i*delta_v;
	next_t=delta_v / ex_velocity(adiabatic(pressure,air_volume,v),nozzle_area);
	sum=sum+(last_t+next_t)/2.0;
/* 	printf("V=%10.4e  next_t=%10.4e\n",v,next_t); */
	last_t=next_t;
    }
    return sum;
  }

  void step_pressure(double delta_t) {
    int i;
    double new_press = 0.0;
    double new_vol = 0.0;
    double delta_v;
    double ave_press;

    /* solve for dP and dV on adiabat */

    ave_press = pressure;
    
    if (air_volume < total_volume)
    {
	for (i=0; i<3; i++) {
        delta_v = delta_t*ex_velocity(ave_press,nozzle_area);
        new_vol = air_volume+delta_v;
        new_press = adiabatic(pressure,air_volume,new_vol);
        ave_press=(pressure+new_press)/2;
      }

	ave_volume = (air_volume+new_vol)/2.0;
	ave_pressure = ave_press;
	pressure = new_press;
	air_volume = new_vol;
    } 
    else {
	pressure = Phy.AMB_PRESSURE;
      ave_pressure = Phy.AMB_PRESSURE;
    }
  }

} 

// *** This class is the canvas used for the graph ***
class GraphCanvas extends Canvas {

  public Dimension minimumSize() {
    return new Dimension(100, 100);
  }

  public Dimension preferredSize() {
    return new Dimension(300, 300);
  }

  public void paint (Graphics g) {
    drawAxis(g);
  }

  public void drawAxis (Graphics g) {
    int w = this.size().width;
    int h = this.size().height;
    g.clearRect(0, 0, w, h);
    g.setColor(Color.black);
    g.drawLine(w/10, h/10, w/10, 9*h/10); // y-axis
    g.drawLine(w/10, h/2, 9*w/10, h/2);   // x-axis
  }

  public void drawLine (Graphics g, Color col,
                        double x1Val, double y1Val, double x2Val, double y2Val,
                        double xScale, double yScale) {
    int x1, y1, x2, y2;
    int w = this.size().width;
    int h = this.size().height;

    g.setColor(col);

    x1 = w/10 + (int) (x1Val * xScale);
    y1 = h/2 - (int) (y1Val * yScale);
    x2 = w/10 + (int) (x2Val * xScale);
    y2 = h/2 - (int) (y2Val * yScale);

    g.drawLine(x1, y1, x2, y2);
  }

  public double calcXScale (double minVal, double maxVal) {
    return ( (8*this.size().width/10) / (maxVal - minVal) ); 
  }  

  public double calcYScale (double minVal, double maxVal) {
    return ( (8*this.size().height/10) / (maxVal - minVal) );
  }  

}

// *** This class implements the extendable drawing lines ***
class GraphLine extends Object {
  GraphCanvas grphCanvas;
  Graphics gr;
  Color clr;
  double xScale, yScale;
  double xOld, yOld;

  // constructor
  public GraphLine (GraphCanvas gc, Color c, 
                    double xMin, double xMax, double yMin, double yMax,
                    double xInit, double yInit) {
    grphCanvas = gc;
    gr = grphCanvas.getGraphics();
    clr = c;
    xScale = grphCanvas.calcXScale(xMin, xMax);
    yScale = grphCanvas.calcYScale(yMin, yMax);
    xOld = xInit;
    yOld = yInit;
  }

  public void extend (double x, double y) {
    grphCanvas.drawLine(gr, clr, xOld, yOld, x, y, xScale, yScale);
    xOld = x;
    yOld = y;
  }
}

// *** This class implements a control for displaying ***
// *** and editing a real number                      *** 
class DoubleReader extends Panel {
  TextField txtIn;
  double min, max;
  double val;  // stores last successful conversion

  // constructor
  public DoubleReader (String labelTxt, double minVal, 
                       double maxVal, double initialVal) 
  {
    setLayout(new FlowLayout(FlowLayout.LEFT, 1, 2));
    add (txtIn = new TextField(Double.toString(initialVal), 6));
    add (new Label(labelTxt));
    val = initialVal;
    min = minVal;
    max = maxVal;
  }

  public double getValue() {
    double tmpVal;

    try {  
      tmpVal = Double.valueOf(txtIn.getText()).doubleValue();
    } 
    catch (Exception e) {  
      txtIn.setText("ERROR");
      tmpVal = val;
    }
    
    if ((tmpVal < min) | (tmpVal > max)) {
      txtIn.setText("RANGE");
      tmpVal = val;
    }
     
    val = tmpVal;
    return val;
  }

}