//******************************************************************************
// QuickGun.java:	Applet
//
//	Java version of Milora's QuickGun code to calculate 2 stage light gas gun
//	dynamics.
//
//	Developed by:  L.R. Baylor		7-Apr-1997 
//
//******************************************************************************
import java.applet.*;
import java.awt.*;
import java.net.URL;
import QuickGunFrame;
import QGmain;
import graph.*;

import IDD_QuickGun;
//==============================================================================
// Main Class for applet QuickGun
//
//==============================================================================
public class QuickGun extends Applet implements Runnable
{
	// THREAD SUPPORT:
	//		m_QuickGun	is the Thread object for the applet
	//--------------------------------------------------------------------------
	Thread	 m_QuickGun = null;
	IDD_QuickGun dlg;


	// ANIMATION SUPPORT:
	//		m_Graphics		used for storing the applet's Graphics context
	//		m_Images[]		the array of Image objects for the animation
	//		m_nCurrImage	the index of the next image to be displayed
	//		m_ImgWidth		width of each image
	//		m_ImgHeight		height of each image
	//		m_fAllLoaded	indicates whether all images have been loaded
	//		NUM_IMAGES 		number of images used in the animation
	//--------------------------------------------------------------------------
	private Graphics m_Graphics;
	private Image	 m_Images[];
	private int 	 m_nCurrImage;
	private int 	 m_nImgWidth  = 0;
	private int 	 m_nImgHeight = 0;
	private boolean  m_fAllLoaded = false;
	private final int NUM_IMAGES = 18;


	private int ticket = 1;

	// QGUN variables
	private float apel, pp, u1, x1, tim, agap, zeta, mu, ml, rhs, mno, mpist;
	private float t, tr, p, pr, accel, ke, sp, xp, densdwn, cdwn, dmpist;
	private float gr, wr, g, w, m, d, a, l, volpt, v, vr, vb, lb, pb, db, lp, mp, ap, mstar;
	private float eseod, topen, cd, ldwn0, sigc, ti, pro, houtpist;
	private float ka, tp, up, sigpel, sigpnew, pstar, sigsta, p1, c;
	private float apnew, apelav, ustar, delchar, term1, term2, achar, s, lambda, tpmean, tmag;
	private float sdot, tfmean, ufluid, amean, afluid,  tl, tup, xdot, sigstaf, frac;
	private float dt, ta, wspeed, upnew, dup, k3, tchar, f2, f1, k1, h, hout, work, hdm;
	private float udwn, ldwn, dens1, densb, machfn, cstar, u2, t0,  x2, pistl, fcoeff, gap;
	private float tr0, massr, mdotv, pdwn, massdwn, tdwn, avalve, c2, c1;
	private float valvefn, pratio, rcrit, time, dm, ui, mass,  dens, astar, uchar;
	private float testx, xmin, umin, nend, ttrans1, ttrans2, trel, p0;
	private float vdwn0,  massdwn0, pdwn0, pr0, dens0, densr, densr0, massr0, mass0, test, ttrans;

	final static int nMaxsize = 2001;

	private float tstar[] = new float[nMaxsize];
	private float sigstar[] = new float[nMaxsize];
	private float fu1[] = new float[nMaxsize];
	private float ftime[] = new float[nMaxsize];
	private float fpbreech[] = new float[nMaxsize];
	private float ftbreech[] = new float[nMaxsize];
	private float fx1[] = new float[nMaxsize];
	private float fxpel[] = new float[nMaxsize];
	private float fupel[] = new float[nMaxsize];
	private float ftp[] = new float[nMaxsize];
	private float fpbase[] = new float[nMaxsize];

	private long ln, li, lj, lk, lna, lnpoint, lthis, lcount;

	private double deps, ddel;

	// STANDALONE APPLICATION SUPPORT:
	//		m_fStandAlone will be set to true if applet is run standalone
	//--------------------------------------------------------------------------
	boolean m_fStandAlone = false;

	// STANDALONE APPLICATION SUPPORT
	// 	The main() method acts as the applet's entry point when it is run
	// as a standalone application. It is ignored if the applet is run from
	// within an HTML page.
	//--------------------------------------------------------------------------
	public static void main(String args[])
	{
		// Create Toplevel Window to contain applet QuickGun
		//----------------------------------------------------------------------
		QuickGunFrame frame = new QuickGunFrame("QuickGun");

		// Must show Frame before we size it so insets() will return valid values
		//----------------------------------------------------------------------
		frame.show();
        frame.hide();
		frame.resize(frame.insets().left + frame.insets().right  + 500,
					 frame.insets().top  + frame.insets().bottom + 400);

		// The following code starts the applet running within the frame window.
		// It also calls GetParameters() to retrieve parameter values from the
		// command line, and sets m_fStandAlone to true to prevent init() from
		// trying to get them from the HTML page.
		//----------------------------------------------------------------------
		QuickGun applet_QuickGun = new QuickGun();

		frame.add("Center", applet_QuickGun);
		applet_QuickGun.m_fStandAlone = true;
		applet_QuickGun.init();
		applet_QuickGun.start();
        frame.show();
	}

	// QuickGun Class Constructor
	//--------------------------------------------------------------------------
	public QuickGun()
	{
		// TODO: Add constructor code here
	}

	// APPLET INFO SUPPORT:
	//		The getAppletInfo() method returns a string describing the applet's
	// author, copyright date, or miscellaneous information.
    //--------------------------------------------------------------------------
	public String getAppletInfo()
	{
		return "Name: QuickGun\r\n" +
		       "Author: L.R. Baylor\r\n" +
		       "Affiliation: Oak Ridge National Laboratory";
	}


	// The init() method is called by the AWT when an applet is first loaded or
	// reloaded.  Override this method to perform whatever initialization your
	// applet needs, such as initializing data structures, loading images or
	// fonts, creating frame windows, setting the layout manager, or adding UI
	// components.
    //--------------------------------------------------------------------------
	public void init()
	{
		LayoutManager m_layout;

        // If you use a ResourceWizard-generated "control creator" class to
        // arrange controls in your applet, you may want to call its
        // CreateControls() method from within this method. Remove the following
        // call to resize() before adding the call to CreateControls();
        // CreateControls() does its own resizing.
        //----------------------------------------------------------------------

		dlg = new IDD_QuickGun( this );
		dlg.CreateControls();

		dlg.IDC_COMBO1.addItem("H2");
		dlg.IDC_COMBO1.addItem("He");
		dlg.IDC_COMBO1.addItem("N2");
		dlg.IDC_ResVol.setText("0.5");
		dlg.IDC_ResPress.setText("32.0");
		dlg.IDC_ResTemp.setText("300");
		dlg.IDC_FillGasPress.setText("2.5");
		dlg.IDC_FillGasTemp.setText("300");
		dlg.IDC_PelRelPress.setText("6.5");
		dlg.IDC_PelLen.setText("2.7");
		dlg.IDC_PelDia.setText("2.7");
		dlg.IDC_PelDens.setText("0.2");
		dlg.IDC_PistMass.setText("10");
		dlg.IDC_PistLen.setText("3.5");
		dlg.IDC_PistClear.setText("0.001");
		dlg.IDC_PistFricCoeff.setText("0.1");
		dlg.IDC_PistDiam.setText("2.5");
		dlg.IDC_PumpTubeLen.setText(".35");
		dlg.IDC_BreechVol.setText("0.25");
		dlg.IDC_BreechMach.setText("0.3");
		dlg.IDC_ValveEseOd.setText("2.5");
		dlg.IDC_ValveRespTime.setText("1.0");
		dlg.IDC_ValveDischCoeff.setText("0.70");
		dlg.IDC_ValvePistDeadSpace.setText("1.0");

		dlg.IDC_COMBO2.addItem("Vpel (km/s)");
		dlg.IDC_COMBO2.addItem("Xpel (m)");
		dlg.IDC_COMBO2.addItem("Pbase (bar)");
		dlg.IDC_COMBO2.addItem("Pbreech (bar)");
		dlg.IDC_COMBO2.addItem("Tbreech (1000K)");
		dlg.IDC_COMBO2.addItem("Vpiston (km/s)");
		dlg.IDC_COMBO3.addItem("Xpel (m)");
		dlg.IDC_COMBO3.addItem("Pbase (bar)");
		dlg.IDC_COMBO3.addItem("Pbreech (bar)");
		dlg.IDC_COMBO3.addItem("Tbreech (1000K)");
		dlg.IDC_COMBO3.addItem("Vpiston (km/s)");
		dlg.IDC_COMBO3.addItem("Vpel (km/s)");


}

	// Place additional applet clean up code here.  destroy() is called when
	// when you applet is terminating and being unloaded.
	//-------------------------------------------------------------------------
	public void destroy()
	{
		// TODO: Place applet cleanup code here
	}

	// ANIMATION SUPPORT
	//		The imageUpdate() method is called repeatedly by the AWT while
	// images are being constructed.  The flags parameter provides information
	// about the status of images under construction.  The AppletWizard's
	// initial implementation of this method checks whether the ALLBITS flag is
	// set.  When the ALLBITS is set, this method knows that an image is
	// completely loaded. When all the images are completely loaded, this
	// method sets the m_fAllLoaded variable to true so that animation can begin.
	//--------------------------------------------------------------------------
	public boolean imageUpdate(Image img, int flags, int x, int y, int w, int h)
	{
		// Nothing to do if images are all loaded
		//----------------------------------------------------------------------
		if (m_fAllLoaded)
			return false;

		// Want all bits to be available before painting
		//----------------------------------------------------------------------
		if ((flags & ALLBITS) == 0)
			return true;

		// All bits are available, so increment loaded count of fully
		// loaded images, starting animation if all images are loaded
		//----------------------------------------------------------------------
		if (++m_nCurrImage == NUM_IMAGES)
		{
			m_nCurrImage = 0;
			m_fAllLoaded = true;
		}				

		return false;
	}

    // ANIMATION SUPPORT:
    //		Draws the next image, if all images are currently loaded
    //--------------------------------------------------------------------------
	private void displayImage(Graphics g)
	{
		if (!m_fAllLoaded)
			return;

		// Draw Image in center of applet
		//----------------------------------------------------------------------
		g.drawImage(m_Images[m_nCurrImage], 
				   (size().width - m_nImgWidth)   / 2,
				   (size().height - m_nImgHeight) / 2, null);
	}

	// QuickGun Paint Handler
	//--------------------------------------------------------------------------
	public void paint(Graphics g)
	{
		// ANIMATION SUPPORT:
		//		The following code displays a status message until all the
		// images are loaded. Then it calls displayImage to display the current
		// image.
		//----------------------------------------------------------------------
//		if (m_fAllLoaded)
//		{
//			Rectangle r = g.getClipRect();
			
//			g.clearRect(r.x, r.y, r.width, r.height);
//			displayImage(g);
//		}
//		else
//			g.drawString("Loading images...", 10, 20);

		// TODO: Place additional applet Paint code here
	}
	//--------------------------------------------------------------------------
	public boolean action(Event event, Object obj)
	{
		Float fFloat;

		// if this action sprung from Run button...
		Object oTarget = event.target;
		if (oTarget instanceof Button)
		{

			String sElement = dlg.IDC_COMBO1.getSelectedItem();
			if (sElement == "He") {
				String sgas = sElement;
				g =1.667f;	/*	ratio of specific heats	*/
				w = 4.0f;	/*	molecular weight	*/
				wr = 4.0f;
				gr =1.667f;	
			}
			else if (sElement == "H2") {
				String sgas = sElement;
				g =1.4f;
				w = 2.0f;
				gr =1.4f;
				wr = 2.0f;
			}
			else if (sElement == "N2") {
				String sgas = sElement;
				g =1.4f;
				w = 28f;
				gr =1.4f;
				wr = 28f;
			}


			String sInput = dlg.IDC_ResVol.getText();
			fFloat = Float.valueOf(sInput);
			vr = fFloat.floatValue()/1000.f;	//units of cc

			sInput = dlg.IDC_ResPress.getText();
			fFloat = Float.valueOf(sInput);
			pr0 = fFloat.floatValue()*6890.0f*14.7f; //units of psi

			sInput = dlg.IDC_ResTemp.getText();
			fFloat = Float.valueOf(sInput);
			tr0 = fFloat.floatValue();

			sInput = dlg.IDC_FillGasPress.getText();
			fFloat = Float.valueOf(sInput);
			p0 = fFloat.floatValue()*6890.0f*14.7f;  //units of psi

			sInput = dlg.IDC_FillGasTemp.getText();
			fFloat = Float.valueOf(sInput);
			t0 = fFloat.floatValue();

			sInput = dlg.IDC_PelRelPress.getText();
			fFloat = Float.valueOf(sInput);
			pb = fFloat.floatValue()*6890.0f*14.7f; //units of psi

			sInput = dlg.IDC_PelLen.getText();
			fFloat = Float.valueOf(sInput);
			lp = fFloat.floatValue()/1000.0f; //units of m

			sInput = dlg.IDC_PelDia.getText();
			fFloat = Float.valueOf(sInput);
			db = fFloat.floatValue()/1000.0f; //units of m

			sInput = dlg.IDC_PelDens.getText();
			fFloat = Float.valueOf(sInput);
			sp = fFloat.floatValue();		//units of g/cc
			mp=3140.f*.25f*sp*lp*db*db;			/*	pellet mass	*/
			ap=.25f*3.14f*db*db;

			sInput = dlg.IDC_PistMass.getText();
			fFloat = Float.valueOf(sInput);
			m = fFloat.floatValue()/1000.0f;	//units of kg

			sInput = dlg.IDC_PistLen.getText();
			fFloat = Float.valueOf(sInput);
			pistl = fFloat.floatValue()/100.0f; //units of m

			sInput = dlg.IDC_PistClear.getText();
			fFloat = Float.valueOf(sInput);
			gap = fFloat.floatValue()/100.0f;  //units of m

			sInput = dlg.IDC_PistFricCoeff.getText();
			fFloat = Float.valueOf(sInput);
			fcoeff = fFloat.floatValue();	  //units of m

			sInput = dlg.IDC_PistDiam.getText();
			fFloat = Float.valueOf(sInput);
			d = fFloat.floatValue()/100.0f;  //units of m
			a=(3.14f/4)*d*d;

			sInput = dlg.IDC_PumpTubeLen.getText();
			fFloat = Float.valueOf(sInput);
			l = fFloat.floatValue();  //units of m
			volpt=1000.f*a*l;		/*	pump tube volume	*/

			sInput = dlg.IDC_BreechVol.getText();
			fFloat = Float.valueOf(sInput);
			vb = fFloat.floatValue();  //units of m3
			vb=.000001f*vb;
			lb=vb/a;

			sInput = dlg.IDC_BreechMach.getText();
			fFloat = Float.valueOf(sInput);
			mstar = fFloat.floatValue();

			sInput = dlg.IDC_ValveEseOd.getText();
			fFloat = Float.valueOf(sInput);
			eseod = fFloat.floatValue()/100.0f;

			sInput = dlg.IDC_ValveRespTime.getText();
			fFloat = Float.valueOf(sInput);
			topen = fFloat.floatValue()*.001f;

			sInput = dlg.IDC_ValveDischCoeff.getText();
			fFloat = Float.valueOf(sInput);
			cd = fFloat.floatValue();

			sInput = dlg.IDC_ValvePistDeadSpace.getText();
			fFloat = Float.valueOf(sInput);
			ldwn0 = fFloat.floatValue();	// units of cm



			//Run the number crunching part of the code
			//QuickGunFrame.setCursor(3);
			qgun_calc();
			//QuickGunFrame.setCursor(0);
			
			return true;
		}
		return false;
	}

	//		The start() method is called when the page containing the applet
	// first appears on the screen. The AppletWizard's initial implementation
	// of this method starts execution of the applet's thread.
	//--------------------------------------------------------------------------
	public void start()
	{
		if (m_QuickGun == null)
		{
			m_QuickGun = new Thread(this);
			m_QuickGun.start();
		}
		// TODO: Place additional applet start code here
	}
	
	//		The stop() method is called when the page containing the applet is
	// no longer on the screen. The AppletWizard's initial implementation of
	// this method stops execution of the applet's thread.
	//--------------------------------------------------------------------------
	public void stop()
	{
		if (m_QuickGun != null)
		{
			m_QuickGun.stop();
			m_QuickGun = null;
		}

		// TODO: Place additional applet stop code here
	}

	// THREAD SUPPORT
	//		The run() method is called when the applet's thread is started. If
	// your applet performs any ongoing activities without waiting for user
	// input, the code for implementing that behavior typically goes here. For
	// example, for an applet that performs animation, the run() method controls
	// the display of images.
	//--------------------------------------------------------------------------
	public void run()
	{
		repaint();

		m_Graphics = getGraphics();
		m_nCurrImage   = 0;
		m_Images   = new Image[NUM_IMAGES];

		// Load in all the images
		//----------------------------------------------------------------------
//		String strImage;

		// For each image in the animation, this method first constructs a
		// string containing the path to the image file; then it begins loading
		// the image into the m_Images array.  Note that the call to getImage
		// will return before the image is completely loaded.
		//----------------------------------------------------------------------
//		for (int i = 1; i <= NUM_IMAGES; i++)
//		{
			// Build path to next image
			//------------------------------------------------------------------
//			strImage = "images/img00" + ((i < 10) ? "0" : "") + i + ".gif";
//			if (m_fStandAlone)
//				m_Images[i-1] = Toolkit.getDefaultToolkit().getImage(strImage);
//			else
//				m_Images[i-1] = getImage(getDocumentBase(), strImage);

			// Get width and height of one image.
			// Assuming all images are same width and height
			//------------------------------------------------------------------
//			if (m_nImgWidth == 0)
//			{
//				try
//				{
					// The getWidth() and getHeight() methods of the Image class
					// return -1 if the dimensions are not yet known. The
					// following code keeps calling getWidth() and getHeight()
					// until they return actual values.
					// NOTE: This is only executed once in this loop, since we
					//       are assuming all images are the same width and
					//       height.  However, since we do not want to duplicate
					//       the above image load code, the code resides in the
					//       loop.
					//----------------------------------------------------------
//					while ((m_nImgWidth = m_Images[i-1].getWidth(null)) < 0)
//						Thread.sleep(1);

//					while ((m_nImgHeight = m_Images[i-1].getHeight(null)) < 0)
//						Thread.sleep(1);						
//				}
//				catch (InterruptedException e)
//				{
					// TODO: Place exception-handling code here in case an
					//       InterruptedException is thrown by Thread.sleep(),
					//		 meaning that another thread has interrupted this one
//				}
//			}

			// Force image to fully load
			//------------------------------------------------------------------
//			m_Graphics.drawImage(m_Images[i-1], -1000, -1000, this);
//		}

		// Wait until all images are fully loaded
		//----------------------------------------------------------------------
//		while (!m_fAllLoaded)
//		{
//			try
//			{
//				Thread.sleep(10);
//			}
//			catch (InterruptedException e)
//			{
				// TODO: Place exception-handling code here in case an
				//       InterruptedException is thrown by Thread.sleep(),
				//		 meaning that another thread has interrupted this one
//			}
//		}
		
//		repaint();

//		while (true)
//		{
//			try
//			{
				// Draw next image in animation
				//--------------------------------------------------------------
//				displayImage(m_Graphics);
//				m_nCurrImage++;
//				if (m_nCurrImage == NUM_IMAGES)
//					m_nCurrImage = 0;

				// TODO:  Add additional thread-specific code here
//				Thread.sleep(50);
//			}
//			catch (InterruptedException e)
//			{
				// TODO: Place exception-handling code here in case an
				//       InterruptedException is thrown by Thread.sleep(),
				//		 meaning that another thread has interrupted this one
//				stop();
//			}
//		}
	}

	// MOUSE SUPPORT:
	//		The mouseDown() method is called if the mouse button is pressed
	// while the mouse cursor is over the applet's portion of the screen.
	//--------------------------------------------------------------------------
	public boolean mouseDown(Event evt, int x, int y)
	{
		// TODO: Place applet mouseDown code here
		return true;
	}

	// MOUSE SUPPORT:
	//		The mouseUp() method is called if the mouse button is released
	// while the mouse cursor is over the applet's portion of the screen.
	//--------------------------------------------------------------------------
	public boolean mouseUp(Event evt, int x, int y)
	{
		// TODO: Place applet mouseUp code here
		return true;
	}



	//--------------------------------------------------------------------------
	public void qgun_calc()
	{
	/*
	Initialyze reservoir and conditions downstream of reservoir and in 
	front of piston.  u is velocity of piston(initially at rest) ,x1 is piston's
	pdwn, tdwn are behind piston.
	*/
	String soutstr;

	int jcount, jj, i, j, k, icount;
	int kcount, kk, n, npoint = 1, na;
	float del, eps;
	boolean b_low_p_bypass = false;
 //   double ddata[] = new double[2*nMaxsize];

	dlg.IDC_DEBUGTXT.setText("Starting Calculation");


	sigc=(2/(g-1))*((float)Math.sqrt(g*8314.4f/w));
	ti=1.5f*(float)Math.sqrt(2*l*m/(.25f*3.14f*eseod*eseod*pr0))+topen;
	/*ti is characteristic piston transit time*/
	agap=3.14f*d*gap/2;
	/*	calculate the mach number for leakage past piston assuming 
	a friction factor f=.005	*/

	mu=1.f;
	ml=0;
	zeta=4.f*.005f*pistl/gap;
	eps=100.f;
	icount=0;

	while(Math.abs(eps/zeta) > .02f)	{
		icount++;
		if (icount>1000)	 {
//	    printf("\n\aSorry: calculation failed to converge in mach number loop");
			dlg.IDC_DEBUGTXT.setText("Sorry: calculation failed to converge in mach number loop");
			return;
	    }
		mno=(mu+ml)/2;
		rhs=(1.f+mno*mno)/(g*mno*mno);
		rhs=rhs+(.5f*(g+1.f)/g)*(float)Math.log((.5f*(g+1.f)*mno*mno)/(1.f+.5f*(g-1.f)*mno*mno));
		eps=zeta-rhs;
		if (eps > 0)
	     	mu=mno;
		else  
	    	ml=mno;
	}
	mpist=mno;


	n=2000;   /* number of time steps for the pellet dynamics calculation.*/
	u1=0.f;
	x1=0.f;
	p=p0;
	t=t0;
	dens0=p0*w/(8314.f*t0);
	dens=dens0;
	densr=pr0*wr/(8314.f*tr0);
	massr0=densr*vr;
	massr=massr0;
	tr=tr0;
	pr=pr0;
	pdwn0=2.f*p0;
	pdwn=pdwn0;
	mass0=dens0*a*(l+lb);
	mass=mass0;
	ldwn0=.01f*ldwn0;  /*ldwn0 is initial distance between piston and valve.*/
	vdwn0=ldwn0*a;
	massdwn0=pdwn0*vdwn0*wr/(8314.f*tr0);
	massdwn=massdwn0;
	tdwn=tr;
	/*
	cd is the valve discharge coefficient. c1, c2 are coefficients. 
	rcrit is critical pressure ratio below which valve is choked.
	Valve description is from Handbook of Fluid Mechanics page 21-30.
	*/
	c1=(float)Math.sqrt(2.f*gr*wr/(8314.f*(gr-1.))); /*gas constant is given by 8314/w*/
	rcrit=(float)Math.pow((2.f/(gr+1.)),(gr/(gr-1.)));
	c2=(float)Math.sqrt((gr*wr/8314.f)*Math.pow((2./(gr+1.)),((gr+1.)/(gr-1.))));

	/* 
	Begin fast loop. For pressures less than release, ie. no movement of the pellet.
	This section just calculates piston dynamics and thermodynamic
	conditions ahead of the piston.  2000 time steps max are allowed.  The loop is
	exited when the downstream pressure exceeds the pellet release value.  If this 
	condition does not occur, a warning message is printed.  
	*/ 

	for (i= 1; i <= 2000; ++i) {

		h=ti/2000.f; /*h is delta time.*/ 
		tim=h*i; 	/* time */
		avalve=(tim<= topen) ? 3.14f*.25f*eseod*eseod*tim/topen : 3.14f*.25f*eseod*eseod;
		accel=pdwn*a/m-p*a/m-fcoeff*(p+pdwn)*d*pistl/m;
		/* accel is piston acceleration.*/
		u2=u1+accel*h;		/*	update velocity	*/
		
		/* Use trapezoidal rule to get new piston position x2. */
		x2=x1+(u2+u1)*(h/2);
		
		/* Calculate valve flow and conditions behind piston(downstream of valve)*/  
		ldwn=ldwn0+x2;		/* distance of piston from valve */
		pratio=(pdwn/pr>=1.f) ? 1.f : pdwn/pr; /* downstream to upstream pressure ratio*/
		valvefn=(pratio < rcrit) ? 1.f:
			(c1/c2)*(float)Math.pow((pratio),(1.f/gr))*(float)Math.sqrt(1.f-(float)Math.pow((pratio),((gr-1.f)/gr)));
		/* Add enthalpy from reservoir at tr to volume behind piston and allow
		piston to expand isentropically; ie. d(me)/dt = h*dm/dt - p*dv/dt   */
	
		mdotv = avalve*cd*c2*pr*valvefn/(float)Math.sqrt(tr);	/* valve mass flow rate	*/
		udwn=massdwn*tdwn;
		hdm=mdotv*h*tr*gr;
		work=pdwn*u2*h*a*(gr-1.f)*wr/8314.f;
		massdwn=massdwn+mdotv*h;
		tdwn=(udwn+hdm-work)/massdwn;
		pdwn=massdwn*8314.f*tdwn/(wr*a*ldwn);
		densdwn=massdwn/(a*ldwn);
		cdwn=(float)Math.sqrt(g*tdwn*8314.f/w);
	
		/* Calculate new reservoir conditions (adiabatic expansion). */
		massr=massr-mdotv*h;
		pr=pr0*(float)Math.pow((massr/massr0),gr);
		tr=tr0*(float)Math.pow((pr/pr0),((gr-1.f)/gr));
	
 		x1=x2;  /* Reset x1 and u1 */
		u1=u2;
	
		/* compress gas in front of piston with mass addition due to leakage
		around piston	 */
	
		c = (float)Math.sqrt(g*t*8314.f/w);
		if(pdwn>=p)		{
    		dmpist = -agap*cdwn*densdwn*mpist*h;
    		houtpist = g*tdwn*dmpist*(1.f+.5f*(g-1.f)*mpist*mpist);
		}
		else	{
    		dmpist = agap*c*dens*mpist*h;
    		houtpist = g*t*dmpist*(1.f+.5f*(g-1.f)*mpist*mpist);
		}
		ui=mass*t;
		work=p*u2*h*a*(g-1.f)*w/8314.f;
		mass=mass-dmpist;
		t=(ui-houtpist+work)/mass;
		p=mass*8314.f*t/(w*a*(l-x2+lb));
		dens=p*w/(8314.f*t);
		ke=.5f*m*u1*u1;		/*	piston kinetic energy	*/
		ttrans1=tim;
	
		/* Find time that pellet releases  */ 
		trel=ttrans1;
		if (p > pb)  break;  /* pellet breaks free, break to slow loop  */
		if ((l-x1) < (.001f*l))  u1=-u1;  /* reflect piston if it approaches head  */
		if (u1 < 0)	 {
    		System.out.println("WARNING: pellet failed to break free");
			//printf("\nWARNING: pellet failed to break free");
    		//printf("\n");
    		return;
		}

	}	/* 	foot of fast loop		*/

	dlg.IDC_DEBUGTXT.setText("Out of fast loop");


	/*
	Out of fast loop. Begin slow loop. Pellet starts to move at pb.
	massO is the initial propellant mass in the pump tube.  Mass loss is
	accounted for as gas flows through barrel at assumed mach no mstar. 
	*/

	/* ti is a characteristic time for the piston to travel from present position
	to the end and  75% of the way back to present position at constant speed.
	*/

	ti=1.75f*(l-x2)/u1;
	densb=dens;  /* densb is density in pump tube at release */
	machfn=1.f/(1.f+.5f*(g-1.f)*mstar*mstar);
	nend=0;
	umin=Math.abs(u1);
	xmin=Math.abs(l-x2);

	jcount = 25;
	jj = 0;
	for(j=0; j<= n-1; ++j) {	/* begin slow loop. Calculate properties of gas in pump tube*/
		npoint=j+1;
		h=ti/n;  /* h is delta time */
		time=tim+j*h;
		accel=pdwn*a/m-p*a/m-fcoeff*(p+pdwn)*d*pistl/m;
		u2=u1+accel*h;
	
		/*  Find piston minimum speed  */
		test=Math.abs(u2);
		if (test < umin) umin = test;
	
		/* Use trapezoidal rule to get piston position x  */
		x2=x1+(u1+u2)*(h/2);
	
		/* Find distance of closest approach */
		testx=Math.abs(l-x2);
		if (testx < xmin) xmin=testx;
	
		/* 
		dens and cstar are conditions at throat. Mass loss is cal-
		culated from these values using isentropic thermodynamic relations.
		Calculate conditions in front of piston allowing mass loss(dm)
		and enthalpy loss(including kinetic energy) through throat.
		*/
		dens=(float)Math.pow((machfn),(1.f/(g-1.f)))*p*w/(8314.f*t);
		cstar=(float)Math.sqrt(machfn)*(float)Math.sqrt(g*t*8314.f/w);
	
		/* 
		Take mass and enthalpy out of volume downstream of piston.
		Remove total enthalpy  at t and allow piston to isentroically compress
		gas ahead of it; ie. d(me)/dt = -(h+.5v^2)*dm/dt + p*dv/dt
		*/
		if(pdwn>=p)		{
    		dmpist = -agap*cdwn*densdwn*mpist*h;
    		houtpist = g*tdwn*dmpist*(1.f+.5f*(g-1.f)*mpist*mpist);
		}
		else	{
			dmpist = agap*c*dens*mpist*h;
			houtpist = g*t*dmpist*(1.f+.5f*(g-1.f)*mpist*mpist);
		}
		dm=h*dens*ap*cstar*mstar;
		ui=mass*t;
		hout=g*t*dm*(1.f+.5f*(g-1.f)*mstar*mstar)+houtpist;
		work=p*u2*h*a*(g-1.f)*w/8314.f;
		mass=mass-dm-dmpist;
		t=(ui-hout+work)/mass;
		p=mass*8314.f*t/(w*a*(l-x2+lb));
		c=(float)Math.sqrt(g*t*8314.f/w);
		/* Calculate valve flow and conditions in front of piston */
		ldwn=ldwn0+x2;
		pratio=(pdwn/pr>=1.f) ? 1.f : pdwn/pr; /* downstream to upstream pressure ratio*/
		avalve=(time<= topen) ? 3.14f*.25f*eseod*eseod*time/topen : 3.14f*.25f*eseod*eseod;
		valvefn=(pratio < rcrit) ? 1.f:
			(c1/c2)*(float)Math.pow((pratio),(1.f/gr))*(float)Math.sqrt(1.-(float)Math.pow((pratio),((gr-1.f)/gr)));
		mdotv = avalve*cd*c2*pr*valvefn/(float)Math.sqrt(tr);
	
		/* Add enthalpy from reservoir at tr to volume behind piston and allow
		piston to expand isentropically; ie. d(me)/dt = h*dm/dt - P*dV/dt.
		*/
		udwn=massdwn*tdwn;	/*	internal energy	*/
		hdm=mdotv*h*tr*gr;		/* enthalpy added  */
		work=pdwn*u2*h*a*(gr-1.f)*wr/8314.f;	/* work done by gas pushing piston	*/
		massdwn=massdwn+mdotv*h;
		tdwn=(udwn+hdm-work)/massdwn;
		pdwn=massdwn*8314.f*tdwn/(wr*a*ldwn);
		densdwn=massdwn/(a*ldwn);
		cdwn=(float)Math.sqrt(g*tdwn*8314.f/w);
		/* Update reservoir conditions. */
		massr=massr-mdotv*h;
		pr=pr0*(float)Math.pow((massr/massr0),gr);
		tr=tr0*(float)Math.pow((pr/pr0),((gr-1.f)/gr));
		x1=x2;	/*	update piston position	*/
		tstar[j]=j*h;		/*	time array	*/
		u1=u2;	/*	update speed	*/
		/* Find piston transit time in fast loop */
		if (u1 > 0) ttrans2 = j*h;	
		if (l-x1 < .001f*l)  u1=-u1;	/*	reflect piston	*/
		/* Sigstar corresponds to the throat conditions. sigma=2*a/(gamma-1)
		t and p are pump tube values(ie. reservoir).
		*/
		t=t0*(float)Math.pow((p/p0),((g-1.)/g));
		sigstar[j]=(2/(g-1.f))*cstar;
		ke=.5f*m*u1*u1;
		if (jcount == 25) {
			//lthis=(long)(j/25.0);
			//if (lthis-j/25.0 == 0) {

			//	fprintf(fp, "\n%3.1f\t %1.3f\t %1.3f\t %3.3f\t %3.2f", tstar[j]*1000000,
			//	 x1, abs(u1)/1000, p/(6890.*10000.0), t/1000.0);
			//save data into arrays
			ftime[jj] = tstar[j]*1000;
			fu1[jj] = Math.abs(u1)/1000;
			fpbreech[jj] = p/100000.0f;
			ftbreech[jj] = t/1000.0f;
			fx1[jj] = x1;
			jcount = 0;
			jj = jj+1;
		}
		jcount = jcount + 1;

	}			/*  foot of slow loop 	*/

	//fprintf(fp, "\n%3.1f\t %1.3f\t %1.3f\t %3.3f\t %3.2f", 0.0, x1, abs(u1)/1000,
	// p/(6890.*10000.0), t/1000.0);
	ftime[jj] = 0.0f;
	fu1[jj] = Math.abs(u1)/1000;
	fpbreech[jj] = p/(100000.0f);
	ftbreech[jj] = t/1000.0f;
	fx1[jj] = x1;

	ttrans=ttrans1+ttrans2;

	//fprintf(fp,"\n%s \t%f\t%s","piston transit time", 1000*ttrans, "ms");
	//fprintf(fp,"\n%s \t%f\t%s","piston impact speed", umin, "m/s");
	//fprintf(fp,"\n%s \t%f\t%s","closest approach", 100*xmin, "cm");
	//fprintf(fp,"\n%s \t%f\t%s","pellet release time", 1000*trel, "ms");

	dlg.IDC_DEBUGTXT.setText("Begin pellet dynamics section");

	/* Begin pellet dynamics calculation section */

	up=0.f;		/*	initial pellet speed	*/
	xp=0.f;		/*	initial pellet position	*/
	ka=2.0f*g/(g-1.f);
	k1=.5f*(g-1.f)*mstar+1.f;

	/* Multiply pb,and densb by conversions to change to throat value. */
	p1=pb*(float)Math.pow(machfn,(g/(g-1.f)));
	k3= p1*ap/mp;
	dens1=densb*(float)Math.pow((machfn),(1.0f/(g-1.f)));
	sigsta=sigstar[0];
	tmag=5;
	dt=tmag*h;
	/* Initialize pressure behind pellet to p1. */
	pp=p1;
	f1=k1*sigsta/sigstar[0];
	f2=up/sigstar[0];
	sigpel=(f1-f2)*sigstar[0];

	//fprintf(fp, "\n%s \t%s,\t%s,\t%s","time(micro s)", "Xpel(m)",
	//"Upel(km/s)", "Pbase(10ksi)");

	kcount = 5;
	kk = 0;
	for (k=1; k <= npoint/tmag; ++k) {		/*	start pellet acceleration loop	*/
		if (pp <= 300000) b_low_p_bypass = true;/*	quit if pressure is too low	 */
		if (pp <= 300000) break;	/*	quit if pressure is too low	 */
		dup=pp*dt*k3/p1;
		/* 
		dup is the incremental pellet velocity change
		pp is pressure behind the pellet.
		*/
		upnew=up+dup;		/*	update pellet speed		*/
		xp=xp+.5f*(up+upnew)*dt;		/*	update pellet position	*/
		up=upnew;
		tp=k*dt;	/*	real time	*/
		apel=.5f*(g-1.f)*sigpel;	/*	sound speed at pellet	*/
	
		del=1.f;		/*	initialize convergence parameter	*/		
		while (Math.abs(del)>0.1f) {		/*	apel iteration loop	 */

			/* 
			perform iteration to find time that pressure pulse started
			from breech along p characteristic and value of sigstar at that time
			(tpmean)
			*/
			tup=tp;		/* upper limit on time	 */
			tl=0;		/*	lower limit on time	 */
			eps=.5f;		/* initialize convergence parameter	 */
			icount=0;
		
			while(Math.abs(eps)>.01){
				icount++;
				if (icount>1000)	 {
    				System.out.println("Sorry: calculation failed to converge in characteristic loop");
    				System.out.println("Try changing input parameters slightly(eg. first stage pressure");
					//printf("\n\aSorry: calculation failed to converge in characteristic loop");
    				//printf("\nTry changing input parameters slightly(eg. first stage pressure)");
    				//printf("\n");
    				return;
    			}
        		tpmean=.5f*(tup+tl);
        		na=(int)(tmag*tpmean/dt);
        		ta=na*dt/tmag;
        		frac=(tpmean-ta)*tmag/dt;
        		sigsta=sigstar[na]+frac*(sigstar[na+1]-sigstar[na]);
        		astar=.5f*(g-1.f)*sigsta;
        		ustar=mstar*astar*up/(mstar*astar+up);
        		xdot=xp/(tp-tpmean);
        		wspeed=.5f*(astar+ustar+apel+up);
        		eps=(xdot-wspeed)/wspeed;
        		if (eps > 0)
        			tup=tpmean;
        		else  
        			tl = tpmean;
    		}
     
			/* have found sigsta for the next time step.
		
			Calculate average speed and sound speed on characteristic
			*/
			uchar=(up+ustar)/2.f;
			achar=(apel+astar)/2.f;
			tup=tp;
			tl=0.f;
		
			/*  perform iteration to find time tfmean that a fluid particle started
			from breech and reached pellet at tp at average speed.  This is
			needed to calculate entropy generation */
		
			eps=.5f;		/* initialize convergence parameter	*/
			icount =0;
		
			while(Math.abs(eps)>.01){
				icount++;
				if (icount>1000)	 {
    				System.out.println("Sorry: calculation failed to converge in ufluid loop");
    				System.out.println("Try changing input parameters slightly(eg. first stage pressure");
    			//printf("\n\aSorry: calculation failed to converge in ufluid loop");
    			//printf("\nTry changing input parameters slightly(eg. first stage pressure)");
    			//printf("\n");
    				return;
    			}		
				tfmean=.5f*(tup+tl);
        		if (tfmean < .0000001) break;
        		na=(int)(tmag*tfmean/dt);
        		ta=na*dt/tmag;
        		frac=(tfmean-ta)*tmag/dt;
        		sigstaf=sigstar[na]+frac*(sigstar[na+1]-sigstar[na]);
        		afluid=.5f*(g-1.f)*sigstaf;
        		ufluid=.5f*(up+afluid*mstar*up/(mstar*afluid+up));
        	
        		/* ufluid is avg. particle speed(avg. of breach and pellet speeds) */
        	
        		xdot=xp/(tp-tfmean);
        		eps=(xdot-ufluid)/ufluid;
        		if (eps > 0)  
        			tup=tfmean;
        		else  
        			tl = tfmean;
    		}
    	
    		/* Calculate entropy change following fluid particle.
			Lambda is friction coefficient defined by f=-(lambda/2d)u^2 where
			f = (1/rho)*grad p. Use avg. lambda from Handbook of Fluid Mechanics. 
			amean is avg. sound speed along particle path between tfmean and tp
			sdot defined as ds/dt = - f*u/temp = -f*u*gamma*R/a^2
			*/
		
			lambda=.02f;
			amean=(apel+afluid)/2.f;
			sdot=.5f*lambda*g*(8314.f/w)*(float)Math.pow(ufluid,3)/(db*amean*amean);
			s=sdot*(tp-tfmean);
		
			/* Calculate change in P Riemann variable(sigma + u) along characteristic.
			This is due to the effect of friction. 
			*/
			term1=-lambda*(tp-tpmean)*uchar*uchar*(1+(g-1)*uchar/achar)/(2.f*db);
			term2=achar*s*w/(g*8314.f);
		
			/* delchar is change in sigma + u on P characteristic connecting tstar & tp.
			signew is value so computed of sigma at the pellet. pp is pressure at
			pellet resulting from new sigma and entropy increase.
			*/
		
			delchar=term1+term2;
			sigpnew=sigsta+ustar+delchar-up;
			apnew=.5f*(g-1.f)*sigpnew;
			apelav=(apel+apnew)/2.f;
			del=(apel-apnew)/apelav;
			apel=apnew;
		
			/* Loop back to WHILE to get convergence of apel and apnew */

		}
		if (b_low_p_bypass == false) {	/*	quit if pressure is too low	 */
			pstar=p1*(float)Math.pow((sigsta/sigstar[0]),(2.f*g/(g-1.f)));
			pp=pstar*(float)Math.pow((sigpnew/sigsta),(2.f*g/(g-1.f)))*(float)Math.exp(-s*w/8314.f);
		}
	
		sigpel=sigpnew;
		if (kcount == 5) {
		//lthis=(long)(k/5.0);
		//if (lthis-k/5.0 == 0) {		/*   print only every 5 time steps	*/
			/*fprintf(fp, "\n%3.1f\t %1.4f\t %1.3f\t %1.3e", tp*1000000, xp, 
			up/1000.0, pp*14.5/1.0e9);*/
			ftp[kk] = tp*1000;
			fupel[kk] = up/1000.0f;
			fpbase[kk] = pp/1.0e5f;
			fxpel[kk] = xp;
			kcount = 0;
			kk = kk+1;
		}
		kcount = kcount+1;

	}		/*	foot of pellet acceleration loop	*/

	dlg.IDC_DEBUGTXT.setText("End of pellet acceleration loop");



	//	output to the screen	
	soutstr = "Piston Transit Time (ms):  " + 1000.0f*ttrans;
	dlg.IDC_TXT1.setText(soutstr);
	soutstr = "Piston Impact Speed (m/s):  " + umin;
	dlg.IDC_TXT2.setText(soutstr);
	soutstr = "Piston Closest Approach (cm):  " + 100.0f*xmin;
	dlg.IDC_TXT3.setText(soutstr);
	soutstr = "Final Pellet Speed (m/s):  " + up;
	dlg.IDC_TXT4.setText(soutstr);

	//Save some data and plot
	QGmain.setnddata(kk-1);

	//Get signals to plot and pass to QCmain
	String sElement = dlg.IDC_COMBO2.getSelectedItem();
	if (sElement.startsWith("Vpe")) {
		QGmain.loadddata1(sElement, ftp, fupel);
	}
	if (sElement.startsWith("Xpe")) {
		QGmain.loadddata1(sElement, ftp, fxpel);
	}
	if (sElement.startsWith("Pba")) {
		QGmain.loadddata1(sElement, ftp, fpbase);
	}
	if (sElement.startsWith("Pbr")) {
		QGmain.loadddata1(sElement, ftime, fpbreech);
	}
	if (sElement.startsWith("Tbr")) {
		QGmain.loadddata1(sElement, ftime, ftbreech);
	}
	if (sElement.startsWith("Vpi")) {
		QGmain.loadddata1(sElement, ftime, fx1);
	}

	sElement = dlg.IDC_COMBO3.getSelectedItem();
	if (sElement.startsWith("Vpe")) {
		QGmain.loadddata2(sElement, ftp, fupel);
	}
	if (sElement.startsWith("Xpe")) {
		QGmain.loadddata2(sElement, ftp, fxpel);
	}
	if (sElement.startsWith("Pba")) {
		QGmain.loadddata2(sElement, ftp, fpbase);
	}
	if (sElement.startsWith("Pbr")) {
		QGmain.loadddata2(sElement, ftime, fpbreech);
	}
	if (sElement.startsWith("Tbr")) {
		QGmain.loadddata2(sElement, ftime, ftbreech);
	}
	if (sElement.startsWith("Vpi")) {
		QGmain.loadddata2(sElement, ftime, fx1);
	}




	//Set flag to plot data
	QGmain.setTicket(ticket);

	//printf("\npiston impact speed(m/s) %f", umin);
	//printf("\npiston closest approach(cm) %f", 100*xmin);
	//printf("\npiston transit time(ms) %f", 1000*ttrans);
	//printf("\nfinal pellet velocity(m/s) %f", up);
	//printf("\npellet release time(ms) %f", 1000*trel);

	return;

	}

}