// Simulation of dust particle contamination 
// JJB 2006

import java.io.*;
import java.util.*;
import java.awt.*;
import java.applet.Applet;
import java.awt.BorderLayout;
import java.awt.event.*;
import java.awt.image.BufferedImage;
import java.awt.image.RenderedImage;
import java.awt.AWTException;
import java.awt.Robot;

import javax.swing.*;
import javax.swing.event.*;
import javax.media.j3d.*;
import javax.vecmath.*;
import javax.imageio.*;


import com.sun.j3d.utils.image.TextureLoader;
import com.sun.j3d.utils.applet.MainFrame;
import com.sun.j3d.utils.geometry.Sphere;
import com.sun.j3d.utils.geometry.Cone;
import com.sun.j3d.utils.universe.*;
import com.sun.j3d.utils.behaviors.mouse.*;

import org.jdesktop.j3d.loaders.vrml97.*;

import com.sun.image.codec.jpeg.*;

//import org.j3d.geom.*;

// v1.0		12 Jun 2007	:	First labelled version
//							Add gradient for wind magnitude based on highest Z found in geometry
// v1.01	15 Jun 2007 :	Fix for if pDestinationPoint is null
// v1.02	18 Jun 2007 :	The test for having reached the edge of the surface is improved:
//							if the particle is found to have moved further away from its destination
//							point than it was in the previous time step, it is marked as having 
//							reached its destination 
//							A particle will not lift off if it has already been subjected to
//							conditions that gave rise to a detachment fraction larger than the
//							current value
//							Correction to formula for direct placement position
// v1.03	18 Jun 2007 :	Logarithmic wind tunnel velocity profile (only for "show motion")
//							New routine "adjustWind" to calculate this
// v1.04	21 Jun 2007	:	Fix for particles with no destination edge
// v1.05	22 Jun 2007	:	Revert to linear wind profile
//							New "SURFACES" option allows sprinkling on ALL, SIDE or TOP surfaces
//							Wind cone now more "3D" looking
// v1.06	25 Jun 2007	:	New test particle (commented) for compatibility check with mathcad
//							3D image now appears in a separate window
//							New special wind profile command in config file
// v1.07	27 Jun 2007	:	New wind cone sizing and position.
//							Starting view position shows the model at a good size, without needing
//							to zoom in or out.
// v1.08	09 Jul 2007	:	Fix for wind cone appearance and position: was not being recalculated correctly
// v1.09	23 Jul 2007	:	Change to minimum grid line spacing: 10cm
// v1.10	26 Jul 2007	:	Fix for scaling output coordinates written to file: missed GEOMETRTSCALEFACTOR
//				

public class Lander extends Applet {
	final String Version = "Version 1.10";
	final static int GUIWidth = 700;
	final static int GUIHeight = 400;
	Canvas3D Con;
	Dimension CD;
	Label L,LMessage,LDistance,Title;
    Thread T;
	JFrame FImage;
	JButton BStart,BAttach;
	JSlider windSpeed,windDirection;
	JCheckBox CBoxDebug, CBoxMinOne, CBoxEdges, CBoxShowMotion, CBoxMovie;
	JTextField TParticles,TRadius;
	TransformGroup objTrans;
	Transform3D transToWorld;
	FrameBehavior cBeh;
	Alpha frameAlpha;
	Bounds bounds;
	BranchGroup objRoot; 
	BranchGroup landerBG;
	BranchGroup surfaceBG;
	TransformGroup tgWind;
	Transform3D landerTransform;
	Cone windCone;
	Vector3d windTranslation;

	Transform3D inverseTransform;
	PolygonAttributes wireAttributes;
	Appearance wireAppearance;
	Material wireMaterial;
//	PickTool pickTool;

	boolean debug = false;
	boolean bMotion = false;
	boolean bViscousDrag = true;
	boolean bLeadingEdge = true;
	boolean bStagnation = true;
	boolean bShowTracks = true;
	boolean bShowNormals = true;
	boolean bBatch = false;

	float lowestZ = 1.e9f;
	float highestZ = 0.f;
	float lowestX = 1.e9f;
	float highestX = 0.f;
	float lowestY = 1.e9f;
	float highestY = 0.f;

	//String vrmlFile = "Rover_Top_Assy_MrSR-TX-ASM.wrl";
	//String vrmlFile = "Body__52354_solid_wrl2a2b1_2.wrl";
	String configFile = "Lander.cfg";
	String vrmlFile = "Rover_Test.wrl";
	String outputFile = "Simulation.dat";
	String windTriplets = "";
	String S;
	String Surfaces = "ALL";

	PrintStream dataStream;

	int nBlobs;

	Vector blobs = null;
	Vector landerShapes = null;

	int nTriangles = 0;
	Triangle triangle[];

	int currentImage = 0; // Used to label image capture JPEG files

	boolean run = false;
	long timestart;

	boolean sharedDone = false;
	
	final float XPLANE = 500.0f;
	final float YPLANE = 500.0f;
	final float ZPLANE = 500.f;

	float XSURFACE = 5000.0f;
	float YSURFACE = 5000.0f;

	double modelMax = 0;

	final double EPSILON = 0.000001;
	final double BIGEPSILON = 0.2;
	final double TENTH_CENTIMETRE = 0.001;
	final double ONE_CENTIMETRE = 0.01;
	final double TEN_CENTIMETRES = 0.1;

	double GEOMETRYSCALEFACTOR = 0.01;


	String blobTypes[]				= {"Dust Particle"};
	Color3f blobColors[]			= {new Color3f(1.0f,1.0f,1.0f)};

	double nParticles1_2			= 0;				// Number of 1-2 micron particles per metre^2  
	double nParticles2_4			= 0;				// Number of 2-4 micron particles per metre^2
	double nParticles4_8			= 0;				// Number of 4-8 micron particles per metre^2  
	double nParticles8_16			= 0;				// Number of 8-16 micron particles per metre^2
	double nParticles16_32			= 0;				// Number of 16-32 micron particles per metre^2  
	double nParticles32_64			= 0;				// Number of 32-64 micron particles per metre^2
	double nParticles64_128			= 0;				// Number of 64-128 micron particles per metre^2  
	double nParticles128_256		= 0;				// Number of 128-256 micron particles per metre^2

	double nPgen1_2					= 0;				// Number of 1-2 micron particles generated  
	double nPgen2_4					= 0;				// Number of 2-4 micron particle generated
	double nPgen4_8					= 0;				// Number of 4-8 micron particles generated  
	double nPgen8_16				= 0;				// Number of 8-16 micron particle generated
	double nPgen16_32				= 0;				// Number of 16-32 micron particles generated  
	double nPgen32_64				= 0;				// Number of 32-64 micron particles generated
	double nPgen64_128				= 0;				// Number of 64-128 micron particles generated  
	double nPgen128_256				= 0;				// Number of 128-256 micron particles generated

	double cleanroom				= 100;				// default cleanroom level: 100 microns
	double particleSize				= 0;				// if non-zero, specifies Gaussian distribution of this mean value
	double particleSigma			= 20;				// width of Gaussian, in microns
	double particleNumber			= 0;				// used for non-zero particleSize. Number to generate.
	double particleDensity			= 1000.0;			// particle density, kg/m^3	(1g/cm^3 == 1000kg/m^3)	

	double windTimes[];									// array of times when wind changes
	double windSpeeds[];								// array of speeds for wind changes
	double windAngles[];								// array of angles for wind changes

	int nWinds;											// number of wind changes 

	double Density;
	double Viscosity;
	double Gravity;

	double GravityEarth				= 9.80665;			// m/s^2
	double DensityEarth				= 1.293;			// Density of Air at NTP, kg/m^3
	double ViscosityEarth			= 0.000015;			// Kinematic Viscosity of Air at 0C, m^2/sec

	double GravityMars				= 3.726527;			// m/s^2
	double DensityMars				= 0.02563;			// Density of Martian atmosphere at 500 Pascals, kg/m^3
	double ViscosityMars			= 0.000561;			// Kinematic Viscosity of Martian Atmosphere, m^2/sec

	double CC						= 0.926;			// Constant for viscous drag calculation

	double dtime					= 0.0005;			// unit time step on calculations (seconds)
	double DEFAULTSCALE				= 0.001;			// Zoom scale for view in 3D


	int nTypes = blobTypes.length;

	int numShapes					= 0;

	double Time						= 0;
	double wind						= 0.0;
	int currentWind					= 0;				// index of the current wind value
	Vector3f windVector				= new Vector3f(1.f,0.f,0.f); // direction of the wind
	double lastWindChangeTime		= 0;				// time at which the wind last changed
	int numWindChanges				= 0;				// number of changes of wind

	int nExperiments				= 1;				// in batch mode, number of experiments to run
	int iExperiment					= 1;				// in batch mode, number of current experiment

    public Lander(String[] args) {

		vrmlFile = null;


		// Set some defaults before reading configuration file
		Density = DensityMars;
		Viscosity = ViscosityMars;
		Gravity = GravityMars;
		String environment = "Mars";
		nWinds = 0; // will be set later, if not in config file

		timestart = System.currentTimeMillis();

		if (args.length > 0)
		{
			configFile = args[0];
		}

		System.out.println("\n\nSimulation: " + Version + "\n\n");
		System.out.println("Config file: " + configFile);

		// Read configuration file
		try
		{
			BufferedReader F = new BufferedReader(new InputStreamReader(new FileInputStream(configFile)));
			while ((S = F.readLine()) != null)
			{
				//System.out.println("Read:" + S);
				StringTokenizer T = new StringTokenizer(S, " :\t");
				while (T.hasMoreTokens())
				{
					String token = new String(T.nextToken());
					if (token.equals("#")) break; // comment
					if (token.equals("VRML"))
					{
						if(vrmlFile == null) vrmlFile = T.nextToken();
						System.out.println("VRML model file: " + vrmlFile);
						break;
					}
					if (token.equals("OUTPUT"))
					{
						outputFile = T.nextToken();
						System.out.println("Simulation output: " + outputFile);
						break;
					}
					if (token.equals("SURFACES"))
					{
						Surfaces = T.nextToken();
						System.out.println("Surfaces: " + Surfaces);
						break;
					}
					if (token.equals("DEBUG"))
					{
						if (T.nextToken().equals("ON"))
						{
							debug = true;
						}
						else
						{
							debug = false;
						}
						break;
					}
					if (token.equals("BATCH"))
					{
						if (T.nextToken().equals("ON"))
						{
							bBatch = true;
						}
						else
						{
							bBatch = false;
						}
						break;
					}
					if (token.equals("MOTION"))
					{
						if (T.nextToken().equals("ON"))
						{
							bMotion = true;
						}
						else
						{
							bMotion = false;
						}
						break;
					}
					if (token.equals("NORMALS"))
					{
						if (T.nextToken().equals("ON"))
						{
							bShowNormals = true;
						}
						else
						{
							bShowNormals = false;
						}
						break;
					}
					if (token.equals("VISCOUS"))
					{
						if (T.nextToken().equals("ON"))
						{
							bViscousDrag = true;
						}
						else
						{
							bViscousDrag = false;
						}
						break;
					}
					if (token.equals("LEADING"))
					{
						if (T.nextToken().equals("ON"))
						{
							bLeadingEdge = true;
							System.out.println("Leading Edge Calculation Determines if Particles will Detach");
						}
						else
						{
							bLeadingEdge = false;
							System.out.println("All Particles will Detach");
						}
						break;
					}
					if (token.equals("STAGNATION"))
					{
						if (T.nextToken().equals("ON"))
						{
							bStagnation = true;
						}
						else
						{
							bStagnation = false;
						}
						break;
					}
					if (token.equals("TRACKS"))
					{
						if (T.nextToken().equals("ON"))
						{
							bShowTracks = true;
						}
						else
						{
							bShowTracks = false;
						}
						break;
					}
					if (token.equals("ENVIRONMENT"))
					{
						environment = T.nextToken();
						System.out.println("Environment: "+environment);
						if(environment.equals("Mars")) {
							Density = DensityMars;
							Viscosity = ViscosityMars;
							Gravity = GravityMars;
						} else if(environment.equals("Earth")) {
							Density = DensityEarth;
							Viscosity = ViscosityEarth;
							Gravity = GravityEarth;
						}

						break;
					}
					if (token.equals("CLEANROOM"))
					{
						cleanroom = Double.valueOf(T.nextToken()).doubleValue();
						break;
					}
					if (token.equals("PARTICLEDENSITY"))
					{
						particleDensity = Double.valueOf(T.nextToken()).doubleValue();
						System.out.println("Particles have density: " + particleDensity + " kg/m^3");
						break;
					}
					if (token.equals("PARTICLESIZE"))
					{
						particleSize = Double.valueOf(T.nextToken()).doubleValue();
						if (particleSize != 0) System.out.println("Particles will be generated Gaussian with mean: " + particleSize + " microns");
						break;
					}
					if (token.equals("PARTICLESIGMA"))
					{
						particleSigma = Double.valueOf(T.nextToken()).doubleValue();
						if (particleSize != 0) System.out.println("Particle size Gaussian distribution sigma: " + particleSigma + " microns");
						break;
					}
					if (token.equals("PARTICLENUMBER"))
					{
						particleNumber = Integer.valueOf(T.nextToken()).intValue();
						System.out.println("If PARTICLESIZE is specified, " + particleNumber + " particles will be generated");
						break;
					}
					if (token.equals("EXPERIMENTS"))
					{
						nExperiments = Integer.valueOf(T.nextToken()).intValue();
						System.out.println("In batch mode there will be " + nExperiments + " experiments");
						break;
					}
					if (token.equals("WIND"))
					{
						// wind parameters are triplets specified as (T,S,D)
						// where T=time, S=speed, D=direction angle
						// e.g. (0,0.1,0)(2,0.2,0)(4,0.4,0)
						windTriplets = new String(S);
						StringTokenizer B = new StringTokenizer(S,"(");
						nWinds = B.countTokens() - 1;
						windTimes = new double[nWinds];
						windSpeeds = new double[nWinds];
						windAngles = new double[nWinds];
						int iwind = 0;
						System.out.println("Wind values specified: "+nWinds);
						String triplet = B.nextToken();
						while(B.hasMoreTokens()) {
							triplet = B.nextToken();
							StringTokenizer t = new StringTokenizer(triplet," ,)");
							windTimes[iwind] = Double.valueOf(t.nextToken()).doubleValue();
							windSpeeds[iwind] = Double.valueOf(t.nextToken()).doubleValue();
							windAngles[iwind] = Double.valueOf(t.nextToken()).doubleValue();
							System.out.println(" Wind: "+windTimes[iwind]+" "+windSpeeds[iwind]+" "+windAngles[iwind]);;
							iwind++;
						}
						if (windTimes[0] < 0)
						{
							// wind profile ... first number indicates delta time between changes
							// second is wind speed (+ve implies constant, -ve implies increment)
							// third is direction (+ve constant, -ve increment)
						}
						else if (windTimes[0] != 0)
						{
							System.out.println("Error: first wind triplet must specify wind at T=0");
							System.exit(0);
						}
						break;
					}
					if (token.equals("SCALE"))
					{
						GEOMETRYSCALEFACTOR = Double.valueOf(T.nextToken()).doubleValue();
						System.out.println("VRML positional data will be scaled by " + GEOMETRYSCALEFACTOR);
						break;
					}

				}
			}
		}
		catch (Exception e)
		{
			System.out.println("Error or no configiration file");
			System.exit(0);
		}

		if (args.length > 1)
		{
			vrmlFile = args[0];
		}

		System.out.println("VRML file: " + vrmlFile);

		if(nWinds == 0) {
			nWinds = 1;
			windTimes = new double[nWinds];
			windSpeeds = new double[nWinds];
			windAngles = new double[nWinds];
			windTimes[0] = 0;
			windSpeeds[0] = 0;
			windAngles[0] = 0;
		}
		currentWind = 0;
		wind = windSpeeds[0];
		double angle = windAngles[0];
		windVector = new Vector3f((float)Math.cos(angle / 360. * 2.0 * Math.PI), (float)Math.sin(angle / 360. * 2.0 * Math.PI), 0.0f);

		if (particleSize == 0)
		{
			System.out.println("Calculating contaminant particle densities for specified cleanroom of " + cleanroom + " microns");

			nParticles1_2 = nPerSquareMetre(1.5, cleanroom);
			nParticles2_4 = nPerSquareMetre(3, cleanroom);
			nParticles4_8 = nPerSquareMetre(6, cleanroom);
			nParticles8_16 = nPerSquareMetre(12, cleanroom);
			nParticles16_32 = nPerSquareMetre(24, cleanroom);
			nParticles32_64 = nPerSquareMetre(48, cleanroom);
			nParticles64_128 = nPerSquareMetre(96, cleanroom);
			nParticles128_256 = nPerSquareMetre(192, cleanroom);

			System.out.println("Particles 1-2 microns:		" + nParticles1_2);
			System.out.println("Particles 2-4 microns:		" + nParticles2_4);
			System.out.println("Particles 4-8 microns:		" + nParticles4_8);
			System.out.println("Particles 8-16 microns:		" + nParticles8_16);
			System.out.println("Particles 16-32 microns:	" + nParticles16_32);
			System.out.println("Particles 32-64 microns:	" + nParticles32_64);
			System.out.println("Particles 64-128 microns:	" + nParticles64_128);
			System.out.println("Particles 128-256 microns:	" + nParticles128_256);

			System.out.println("Check for 100 microns:		" + nPerSquareMetre(100, cleanroom));
		}
		else
		{
			System.out.println(particleNumber+" particles per square meter, of size " + particleSize + " microns, will be generated");
		}


		setLayout(new BorderLayout());

		JPanel pMessages = new JPanel();
		pMessages.setLayout(new GridLayout(1,2));

		LMessage = new Label("                                                                                         ");
		pMessages.add(LMessage);

		LDistance = new Label("");
		pMessages.add(LDistance);

		add("North",pMessages);

		JPanel pControls = new JPanel();
		pControls.setLayout(new GridLayout(3,1));

		JPanel pView = new JPanel();
		pView.setLayout(new GridLayout(6,2));

		JButton BSide = new JButton("Top View");
		BSide.addActionListener(new eventHandler());
		pView.add(BSide);

		JButton BEnd  = new JButton("Side View");
		BEnd.addActionListener(new eventHandler());
		pView.add(BEnd);

		JButton BZoomIn = new JButton("Zoom In");
		BZoomIn.addActionListener(new eventHandler());
		pView.add(BZoomIn);

		JButton BZoomOut = new JButton("Zoom Out");
		BZoomOut.addActionListener(new eventHandler());
		pView.add(BZoomOut);

		JButton BLeft = new JButton("Left");
		BLeft.addActionListener(new eventHandler());
		pView.add(BLeft);

		JButton BRight = new JButton("Right");
		BRight.addActionListener(new eventHandler());
		pView.add(BRight);

		JButton BNew = new JButton("New");
		BNew.addActionListener(new eventHandler());
		pView.add(BNew);

		if (bBatch)
		{
			BStart = new JButton("Stop");
		}
		else
		{
			BStart = new JButton("Start");
		}
		BStart.addActionListener(new eventHandler());
		pView.add(BStart);

		BAttach = new JButton("Hide Lander");
		BAttach.addActionListener(new eventHandler());
		pView.add(BAttach);

		JButton BSnapShot = new JButton("Save View as JPEG");
		BSnapShot.addActionListener(new eventHandler());
		pView.add(BSnapShot);

		JButton BWriteData = new JButton("Write Data File");
		BWriteData.addActionListener(new eventHandler());
		pView.add(BWriteData);

		pControls.add(pView);

		JPanel pOptions = new JPanel(new GridLayout(1,4));

		CBoxDebug = new JCheckBox("Debug", debug);
		pOptions.add(CBoxDebug);

		CBoxShowMotion = new JCheckBox("Show Particle Motion",bMotion);
		pOptions.add(CBoxShowMotion);

		CBoxEdges = new JCheckBox("Show Tracks", bShowTracks);
		pOptions.add(CBoxEdges);

		CBoxMovie = new JCheckBox("Save Image Every Step", false);
		pOptions.add(CBoxMovie);

		pControls.add(pOptions);

		JPanel pParameters = new JPanel();
		pParameters.setLayout(new GridLayout(1,2));

		windSpeed = new JSlider(JSlider.HORIZONTAL,0,5000,0);
		windSpeed.addChangeListener(new SliderListener());

		windSpeed.setBorder(BorderFactory.createTitledBorder("Wind Speed (cm/s)"));
		windSpeed.setMajorTickSpacing(500);
		windSpeed.setMinorTickSpacing(100);
		windSpeed.setPaintTicks(true);
		windSpeed.setPaintLabels(true);

		pParameters.add(windSpeed);

		windDirection = new JSlider(JSlider.HORIZONTAL,0,360,0);
		windDirection.addChangeListener(new SliderListener());

		windDirection.setBorder(BorderFactory.createTitledBorder("Wind Direction (degrees)"));
		windDirection.setMajorTickSpacing(60);
		windDirection.setMinorTickSpacing(10);
		windDirection.setPaintTicks(true);
		windDirection.setPaintLabels(true);

		pParameters.add(windDirection);

		pControls.add(pParameters);

		add("Center",pControls);


		// Now the Scene Graph and Canvas etc.

		GraphicsConfiguration config = SimpleUniverse
			.getPreferredConfiguration();


		BranchGroup scene = createSceneGraph();

		Con = new Canvas3D(config);

		SimpleUniverse u = new SimpleUniverse(Con);

		//add("Center", Con);

		FImage = new JFrame("3D Lander Simulation");
		JPanel PC = new JPanel(new BorderLayout());
		PC.add(Con, BorderLayout.CENTER);
		FImage.getContentPane().add(PC);
		FImage.setSize(800, 800);
		FImage.setLocation(new Point(this.getX() + GUIWidth + 10, this.getY()));
		//FImage.pack();
		FImage.setVisible(true);
		FImage.addWindowListener(new WindowAdapter()
		{
			public void windowClosing(WindowEvent e) { System.exit(0); }
		});


		ViewingPlatform vp = u.getViewingPlatform();
        vp.setNominalViewingTransform();

		u.addBranchGraph(scene);

		setupTriangles();

		System.out.println("Lowest Z coordinate found "+lowestZ);
		System.out.println("Highest Z coordinate found " + highestZ);

		// Scale the view so that the model is visible

		Transform3D t1 = new Transform3D();
		objTrans.getTransform(t1);
		double scale = t1.getScale();
		scale *= 200/(highestZ-lowestZ);
		t1.setScale(scale);
		objTrans.setTransform(t1);

		createSurface();

		createWindCone();

		initializeSim();

		if (bBatch) run = true;

    }

    public static void main(String[] args) {
			new MainFrame(new Lander(args), GUIWidth, GUIHeight);
    }

	public double nPerSquareMetre(double d, double dc)
	{
		// Implements Partha's calculation of the number of contaminant particles per square metre
		// of size d microns, given a cleanroom specified to dc microns
		if (dc == 0) return 0;
		double constant = 0.926;
		double ml10_2 = 1/(Math.log(10)*Math.log(10));
		if (d < dc) constant = 0.33;

		double n = -constant*Math.log(d)*Math.log(d) * ml10_2;
		n = Math.pow(10, n);
		n *= Math.pow(10, constant * Math.log(dc)* Math.log(dc) * ml10_2);

		// convert to number per square metre

		//n *= 0.09290304;
		return n;
	}

	   public void takeSnapShot() {
			System.out.print("Snapshot...");
			String pref = "00000";
			currentImage++;
			if(currentImage > 9) pref = "0000";
			if(currentImage > 99) pref = "000";
			if(currentImage > 999) pref = "00";
			if(currentImage > 9999) pref = "0";
			if (currentImage > 99999) pref = "";
			Point location = Con.getLocationOnScreen();
			System.out.println("Window " + location.x + " " + location.y + " Con " + Con.getX() + " " + Con.getY());

			int xpos = location.x;// + Con.getX();
			int ypos = location.y;// + Con.getY();

		   try {
			//BufferedImage screencapture = new Robot().createScreenCapture(
			//	  new Rectangle(Toolkit.getDefaultToolkit().getScreenSize()));

			BufferedImage screencapture = new Robot().createScreenCapture(
			   new Rectangle(xpos, ypos,
							  Con.getWidth(), Con.getHeight()));
			   // Save as JPEG
			File file = new File("Lander"+pref+currentImage+".png");
			ImageIO.write(screencapture, "png", file);
		}
		catch (Exception e) { System.err.println(e); }
		System.out.println("Done");

	   }

	public void createWindCone()
	{

		tgWind = new TransformGroup();
		tgWind.setCapability(TransformGroup.ALLOW_TRANSFORM_WRITE);
		tgWind.setCapability(TransformGroup.ALLOW_TRANSFORM_READ);

		BranchGroup windBG = new BranchGroup();
		Transform3D twind = new Transform3D();
		windTranslation = new Vector3d();
		windTranslation.x = -(3*modelMax);
		windTranslation.y = 0;
		windTranslation.z = (highestZ);
		Matrix3d twist = new Matrix3d();
		twist.rotZ(-Math.PI * 0.5);

		twind.setRotation(twist);
		twind.setTranslation(windTranslation);

		tgWind.setTransform(twind);

		Color3f c = new Color3f(0.9f, 0.9f, 0.1f);

		Appearance windAppearance = new Appearance();
		ColoringAttributes windColoring = new ColoringAttributes();
		windColoring.setColor(c);
		windAppearance.setColoringAttributes(windColoring);

		Material mat = new Material();
		mat.setAmbientColor(c);
		mat.setSpecularColor(new Color3f(1.0f, 1.0f, 1.0f));
		mat.setDiffuseColor(c);
		mat.setShininess(12.0f);
		mat.setLightingEnable(true);

		windAppearance.setMaterial(mat);

		float radius = (float)(modelMax * 0.05);
		float height = (float)(modelMax * 0.5);

		System.out.println("Creating wind cone of radius " + radius + " height " + height + " at " + windTranslation);

		windCone = new Cone(radius, height, Cone.GENERATE_NORMALS, windAppearance);

		tgWind.addChild(windCone);
		windBG.addChild(tgWind);
		objTrans.addChild(windBG);

	}

	public void createSurface()
	{
		System.out.println("Creating surface at Z=" + lowestZ);

		modelMax = (highestX - lowestX);
		if (highestY - lowestY > modelMax) modelMax = highestY - lowestY;
		if (highestZ - lowestZ > modelMax) modelMax = highestZ - lowestZ;

		System.out.println("Model extent: " + modelMax);

		XSURFACE = YSURFACE = (float)(int)(modelMax * 10);

		System.out.println("Surface size " + 2 * XSURFACE + "x" + 2 * YSURFACE + " m^2");

		/*

		Color3f surfaceColor = new Color3f(0.0f, 0.0f, 0.0f);
		Box surfaceBox = new Box(-XSURFACE / 2, XSURFACE / 2, -YSURFACE / 2, YSURFACE / 2, lowestZ - 0.1f, lowestZ, surfaceColor);
		Shape3D shape = surfaceBox.getShape3D();
		shape.setCapability(Shape3D.ALLOW_GEOMETRY_READ);
		shape.setCapability(Shape3D.ALLOW_APPEARANCE_READ);
		shape.setCapability(Shape3D.ALLOW_APPEARANCE_WRITE);
		surfaceBG = surfaceBox.getBranchGroup();
		objTrans.addChild(surfaceBG);

		*/

		// Create the surface grid
		double gridSpacing = 10.*(double)((int)(XSURFACE/20)/10);
		if (gridSpacing <= 0) gridSpacing = 0.1;

		System.out.println("Grid line spacing: " + gridSpacing + "m");
		int nlinesX = (int)(XSURFACE / gridSpacing) + 1;
		int nlinesY = (int)(YSURFACE / gridSpacing) + 1;
		LineArray surfaceGridLines = new LineArray(2 * nlinesX + 2 * nlinesY, LineArray.COORDINATES);

		int icoord = 0;
		for (int i = 0; i < nlinesX; i++)
		{
			Point3d coord = new Point3d(-XSURFACE / 2 + i * gridSpacing, -YSURFACE / 2, lowestZ);
			surfaceGridLines.setCoordinate(icoord, coord);
			icoord++;
			coord = new Point3d(-XSURFACE / 2 + i * gridSpacing, YSURFACE / 2, lowestZ);
			surfaceGridLines.setCoordinate(icoord, coord);
			icoord++;
		}
		for (int i = 0; i < nlinesY; i++)
		{
			Point3d coord = new Point3d(-XSURFACE / 2, -YSURFACE / 2 + i * gridSpacing, lowestZ);
			surfaceGridLines.setCoordinate(icoord, coord);
			icoord++;
			coord = new Point3d(XSURFACE / 2, -YSURFACE / 2 + i * gridSpacing, lowestZ);
			surfaceGridLines.setCoordinate(icoord, coord);
			icoord++;
		}

		Appearance surfaceGridAppearance = new Appearance();
		LineAttributes surfaceLineAttributes = new LineAttributes(1.0f, LineAttributes.PATTERN_SOLID, false);
		surfaceGridAppearance.setLineAttributes(surfaceLineAttributes);

		ColoringAttributes gridColoring = new ColoringAttributes();
		gridColoring.setColor(new Color3f(0.5f, 0.5f, 0.5f));
		surfaceGridAppearance.setColoringAttributes(gridColoring);

		Shape3D surfaceGrid = new Shape3D(surfaceGridLines, surfaceGridAppearance);
		BranchGroup gridBG = new BranchGroup();
		gridBG.addChild(surfaceGrid);

		objTrans.addChild(gridBG);
	}
    
	public BranchGroup createSceneGraph() {

		// Create the root of the branch graph
		objRoot = new BranchGroup();
		//objRoot.setCapability(BranchGroup.ALLOW_CHILDREN_WRITE);
		//objRoot.setCapability(BranchGroup.ALLOW_CHILDREN_EXTEND);

		objTrans = new TransformGroup();
		objTrans.setCapability(TransformGroup.ALLOW_TRANSFORM_WRITE);
		objTrans.setCapability(TransformGroup.ALLOW_TRANSFORM_READ);
		objTrans.setCapability(TransformGroup.ALLOW_CHILDREN_READ);
		objTrans.setCapability(TransformGroup.ALLOW_CHILDREN_WRITE);
		objTrans.setCapability(TransformGroup.ALLOW_CHILDREN_EXTEND);
		objTrans.setCapability(TransformGroup.ALLOW_LOCAL_TO_VWORLD_READ);

		objRoot.addChild(objTrans);

		landerTransform = new Transform3D();
		landerTransform.rotX(-0.5*Math.PI);
		landerTransform.setScale(DEFAULTSCALE);

		objTrans.setTransform(landerTransform);

		inverseTransform = new Transform3D();
		inverseTransform.rotX(0.5*Math.PI);
		inverseTransform.setScale(1/DEFAULTSCALE);


		wireAttributes = new PolygonAttributes();
		wireAttributes.setPolygonMode(PolygonAttributes.POLYGON_LINE);
		wireAttributes.setCullFace(PolygonAttributes.CULL_NONE);

		wireMaterial = new Material();
		Color3f wireColor = new Color3f(0.0f,0.0f,1.0f);
		Color3f specColor = new Color3f(0.5f, 0.5f, 0.5f);
		wireMaterial.setAmbientColor(wireColor);
		wireMaterial.setSpecularColor(specColor);
		wireMaterial.setDiffuseColor(wireColor);
		wireMaterial.setLightingEnable(true);

		wireAppearance = new Appearance();
		wireAppearance.setPolygonAttributes(wireAttributes);
		wireAppearance.setMaterial(wireMaterial);


		System.out.println("Loading VRML file: "+vrmlFile);
		VrmlLoader loader = new VrmlLoader();

		landerShapes = new Vector(10000);

		VrmlScene scene = null;
		try {
			scene = (VrmlScene) loader.load(vrmlFile);
		} catch (Exception e) {
			System.err.println("Error: "+e);
			System.exit(0);
		}

		landerBG = scene.getSceneGroup();

		recursiveGetShapes(landerBG);

		System.out.println("TOTAL: "+numShapes+" shapes read");


		landerBG.setCapability(BranchGroup.ALLOW_BOUNDS_READ);
		landerBG.setCapability(BranchGroup.ALLOW_CHILDREN_READ);
		landerBG.setCapability(BranchGroup.ALLOW_DETACH);

		objTrans.addChild(landerBG);

		BoundingSphere bl = (BoundingSphere) landerBG.getBounds();
	
		System.out.println("Lander bounds "+bl);

		//bounds = new BoundingSphere(new Point3d(0.0, 0.0, 0.0), 10000.0);


		Point3d bc = new Point3d();
		bl.getCenter(bc);
		double radius = bl.getRadius();
		bounds = new BoundingSphere(bc, 100.0*radius);


		BranchGroup frameBG = new BranchGroup();
		frameBG.setCapability(BranchGroup.ALLOW_DETACH);


		frameAlpha = new Alpha(-1, Alpha.INCREASING_ENABLE |
				     Alpha.DECREASING_ENABLE,
				     0, 0,
				     2000, 1000, 200,
				     2000, 1000, 200);

		// Create the frame behavior
		cBeh = new FrameBehavior(frameAlpha,objTrans);
		cBeh.setSchedulingBounds(bounds);
		objRoot.addChild(cBeh);

        // Set up the background
        Color3f bgColor = new Color3f(0.1f, 0.1f, 0.1f);
        Background bgNode = new Background(bgColor);
        bgNode.setApplicationBounds(bounds);
        objRoot.addChild(bgNode);

		// Set up the lighting
		BranchGroup lightsBG = new BranchGroup();
		lightsBG.setCapability(BranchGroup.ALLOW_DETACH);

		DirectionalLight lgt7 = new DirectionalLight();
		lgt7.setColor(new Color3f( 1.0f, 1.0f, 1.0f ));
		lgt7.setDirection(0.0f,-1.0f,-1.0f);
		lgt7.setInfluencingBounds(bounds);
		lgt7.setEnable(true);

		lightsBG.addChild(lgt7);

		AmbientLight lgt1 = new AmbientLight(new Color3f( 0.2f, 0.2f, 0.2f ));
		lgt1.setInfluencingBounds(bounds);
		lgt1.setEnable(true);

		lightsBG.addChild(lgt1);

		DirectionalLight lgt2 = new DirectionalLight();
		lgt2.setColor(new Color3f( 1.0f, 1.0f, 1.0f ));
		lgt2.setDirection(1.0f,1.0f,1.0f);
		lgt2.setInfluencingBounds(bounds);
		lgt2.setEnable(true);

		lightsBG.addChild(lgt2);

		DirectionalLight lgt3 = new DirectionalLight();
		lgt3.setColor(new Color3f(1.0f, 1.0f, 1.0f));
		lgt3.setDirection(-1.0f, -1.0f, 1.0f);
		lgt3.setInfluencingBounds(bounds);
		lgt3.setEnable(true);

		//lightsBG.addChild(lgt3);		

		objTrans.addChild(lightsBG);

		// Create the rotate behavior node
	    MouseRotate behavior = new MouseRotate();
		behavior.setTransformGroup(objTrans);
		objRoot.addChild(behavior);
		behavior.setSchedulingBounds(bounds);

        // Have Java 3D perform optimizations on this scene graph.
        //objRoot.compile()

		return objRoot;
    }

	void recursiveGetShapes(Object value) {
		if(value instanceof SceneGraphObject) {
			SceneGraphObject sg = (SceneGraphObject) value;

			if(sg instanceof Shape3D) {
				Shape3D shape = (Shape3D) sg;

				shape.setCapability(Shape3D.ALLOW_GEOMETRY_READ);
				shape.setCapability(Shape3D.ALLOW_GEOMETRY_WRITE);
				shape.setCapability(Shape3D.ALLOW_APPEARANCE_READ);
				shape.setCapability(Shape3D.ALLOW_APPEARANCE_WRITE);
				shape.setCapability(Shape3D.ALLOW_APPEARANCE_OVERRIDE_WRITE);
				shape.setCapability(Shape3D.ALLOW_LOCAL_TO_VWORLD_READ);

				int numGeometries = shape.numGeometries();
				for (int i = 0; i < numGeometries; i++)
				{
					Geometry g = shape.getGeometry(i);
					g.setCapability(Geometry.ALLOW_INTERSECT);
					if (g instanceof GeometryArray)
					{
						GeometryArray ga = (GeometryArray)g;
						ga.setCapability(GeometryArray.ALLOW_COUNT_READ);
						ga.setCapability(GeometryArray.ALLOW_FORMAT_READ);
						ga.setCapability(GeometryArray.ALLOW_COORDINATE_READ);
						ga.setCapability(GeometryArray.ALLOW_NORMAL_READ);
					}
				}

				//Appearance a = shape.getAppearance();
				//a.setCapability(Appearance.ALLOW_POLYGON_ATTRIBUTES_READ);
				//a.setCapability(Appearance.ALLOW_POLYGON_ATTRIBUTES_WRITE);
				//a.setCapability(Appearance.ALLOW_MATERIAL_READ);
				//a.setCapability(Appearance.ALLOW_MATERIAL_WRITE);

				//a.setPolygonAttributes(wireAttributes);
				//a.setMaterial(wireMaterial);

				shape.setAppearanceOverrideEnable(true);
				shape.setAppearance(wireAppearance);

				landerShapes.add(shape);

				numShapes++;

			} else if(sg instanceof Group) {
				Group g = (Group) sg;
				//g.setCapability(Group.ENABLE_COLLISION_REPORTING);
				g.setCapability(Group.ALLOW_LOCAL_TO_VWORLD_READ);
				g.setCapability(Group.ALLOW_CHILDREN_WRITE);
				g.setCapability(Group.ALLOW_CHILDREN_EXTEND);
				// recurse on children
				java.util.Enumeration enumChildren = g.getAllChildren( );
				while( enumChildren.hasMoreElements())	
					recursiveGetShapes( enumChildren.nextElement());
			} else if(sg instanceof Link) {
				Link l = (Link) sg;
				SharedGroup shg = l.getSharedGroup();
				System.out.println("Shared Group "+shg);

				BranchGroup bg = new BranchGroup();

				for(int ic=0;ic<shg.numChildren();ic++) {
					Node child = shg.getChild(ic);
					if (child instanceof Light) continue;
					System.out.println("Cloning child "+child+" to new BranchGroup");
					Node cloned = child.cloneTree(true);
					bg.addChild(cloned);
				}

				Group parentOfLink = (Group) l.getParent();

				//if(parentOfLink == null) parentOfLink = objTrans;
				parentOfLink.addChild(bg);

				java.util.Enumeration enumChildren = bg.getAllChildren( );
				while( enumChildren.hasMoreElements())	
					recursiveGetShapes( enumChildren.nextElement());


			}
			else if (sg instanceof Leaf)
			{
				System.out.println("Untreated leaf " + numShapes + " " + sg);
			}
			else
			{
				System.out.println("Untreated scenegraph object "+sg);
			}
		}
	}

	void destroyParticles() {
		if(blobs != null) {

			for(int i=0;i<blobs.size();i++) {

				Blob b = (Blob) blobs.elementAt(i);

				BranchGroup blobBG = b.blobBG;

				blobBG.detach();

				blobBG = null;

			}
			System.gc();		
		}	
	}

	void destroyStagnationNormals()
	{
		java.util.Enumeration enumChildren = objTrans.getAllChildren();
		while (enumChildren.hasMoreElements())
		{
			Node n = (Node)enumChildren.nextElement();
			if (n instanceof BranchGroup)
			{
				String user = (String) n.getUserData();
				if (user == null) continue;
				if (!user.equals("Stag") && !user.equals("Norm")) continue;
				BranchGroup bn = (BranchGroup)n; 
				bn.detach();
				bn = null;
			}
		}
		System.gc();
	}

	void setupTriangles()
	{

		Point3d coords[];
		Vector3f normals[];

		// Looks at all triangles in the geometry, calculates area and saves each

		nTriangles = 0;

		for(int is=0;is<landerShapes.size();is++) {
			Shape3D shape = (Shape3D) landerShapes.elementAt(is);
			TriangleArray ta = shapeGetTriangles(shape);
			if(ta == null) {
				System.out.println("Shape is null!");
				System.exit(0);
				continue;
			}
			nTriangles += ta.getVertexCount()/3;
		}

		triangle = new Triangle[nTriangles];
		
		nTriangles = 0;

		for(int is=0;is<landerShapes.size();is++) {

			Shape3D shape = (Shape3D) landerShapes.elementAt(is);

			Transform3D lt = new Transform3D();

			try {
				shape.getLocalToVworld(lt);
			} catch (Exception e) {
				System.out.println("Error "+e);
				System.exit(0);
				continue;
			}

			TriangleArray ta = shapeGetTriangles(shape);

			int nvertices = ta.getVertexCount();

			coords = new Point3d[nvertices];
			normals = new Vector3f[nvertices];

			for(int iv=0;iv<nvertices;iv++) {
				coords[iv] = new Point3d();
				normals[iv] = new Vector3f();
			}
			ta.getCoordinates(0,coords);
			ta.getNormals(0,normals);

			int nn = (nvertices/3)-1; // number of neighbouring triangles for each triangle
			int n[] = new int[nn];
			int firstTriangle = nTriangles;
			int lastTriangle = nTriangles + (nvertices/3);

			for(int iv=0;iv<nvertices;iv+=3) {

				Point3d v1 = new Point3d(coords[iv].x, coords[iv].y, coords[iv].z);
				Point3d v2 = new Point3d(coords[iv+1].x, coords[iv+1].y, coords[iv+1].z);
				Point3d v3 = new Point3d(coords[iv+2].x, coords[iv+2].y, coords[iv+2].z);
													
				lt.transform(v1);
				inverseTransform.transform(v1);
				lt.transform(v2);
				inverseTransform.transform(v2);
				lt.transform(v3);
				inverseTransform.transform(v3);

				if(v3.z < lowestZ) lowestZ = (float) v3.z;
				if(v2.z < lowestZ) lowestZ = (float) v2.z;
				if(v1.z < lowestZ) lowestZ = (float) v1.z;
				if (v3.z > highestZ) highestZ = (float)v3.z;
				if (v2.z > highestZ) highestZ = (float)v2.z;
				if (v1.z > highestZ) highestZ = (float)v1.z;

				if (v3.x < lowestX) lowestX = (float)v3.x;
				if (v2.x < lowestX) lowestX = (float)v2.x;
				if (v1.x < lowestX) lowestX = (float)v1.x;
				if (v3.x > highestX) highestX = (float)v3.x;
				if (v2.x > highestX) highestX = (float)v2.x;
				if (v1.x > highestX) highestX = (float)v1.x;

				if (v3.y < lowestY) lowestY = (float)v3.y;
				if (v2.y < lowestY) lowestY = (float)v2.y;
				if (v1.y < lowestY) lowestY = (float)v1.y;
				if (v3.y > highestY) highestY = (float)v3.y;
				if (v2.y > highestY) highestY = (float)v2.y;
				if (v1.y > highestY) highestY = (float)v1.y;


				Vector3f normal = new Vector3f(normals[iv]);

				lt.transform(normal);
				inverseTransform.transform(normal);

				triangle[nTriangles] = new Triangle(v1,v2,v3,normal);

				// Store the neighbouring triangles, i.e. the indices of the triangles in
				// the same shape

				int in=0;
				for(int it=firstTriangle;it<lastTriangle;it++) {
					if(it == nTriangles) continue;
					n[in++] = it;
				} 

				triangle[nTriangles].setNeighbours(nn,n);

				nTriangles++;
			}

			// For the triangles in this shape, find the groups which are coplanar,
			// and for each group, calculate the centre of mass of the group,
			// which is used as the stagnation point when wind is blowing
			// perpendicular to the group

			int ngroups = 1;
			int ninshape = nvertices / 3;
			int group[] = new int[ninshape];
			for(int i=0;i<ninshape;i++) group[i] = -1;
			group[0] = 1;
			if (CBoxDebug.isSelected()) System.out.println("Shape with " + ninshape + " triangles");
			for (int it1 = firstTriangle; it1 < lastTriangle-1; it1++)
			{
				if(group[it1-firstTriangle] == -1) {
					ngroups++;
					group[it1-firstTriangle] = ngroups;
				}
				for (int it2 = it1 + 1; it2 < lastTriangle; it2++)
				{
					boolean coplanar = CoplanarTriangles(triangle[it1], triangle[it2]);
					if (CBoxDebug.isSelected()) System.out.println("Comparing " + triangle[it1].print() + " with " + triangle[it2].print() + " Coplanar=" + coplanar);
					if(coplanar) {
						group[it2-firstTriangle] = ngroups;
					}
				}
			}
			if(group[ninshape-1] == -1) {
				ngroups++;
				group[ninshape-1] = ngroups;
			}
			if (CBoxDebug.isSelected())
			{
				System.out.println("Groups " + ngroups);
				for (int i = 0; i < ninshape; i++)
				{
					System.out.println("  Triangle " + i + " Group " + group[i]);
				}
			}

			// For each group, calculate the centre of mass
			for (int ig = 0; ig < ngroups; ig++)
			{
				double sumarea = 0;
				double wX = 0;
				double wY = 0;
				double wZ = 0;
				for (int it = firstTriangle; it < lastTriangle; it++)
				{
					if (group[it - firstTriangle] != (ig+1)) continue;
					double area = triangle[it].area;
					Point3d centroid = triangle[it].centroid;
					sumarea += area;
					wX += centroid.x * area;
					wY += centroid.y * area;
					wZ += centroid.z * area;
				}
				Point3d CofM = new Point3d(wX/sumarea,wY/sumarea,wZ/sumarea);
				if (CBoxDebug.isSelected()) System.out.println("Group " + (ig + 1) + " Centre of Mass " + CofM);
				// Store this position in each Triangle
				for (int it = firstTriangle; it < lastTriangle; it++)
				{
					if (group[it - firstTriangle] != (ig + 1)) continue;
					triangle[it].setStagnation(CofM);
				}

			}

		}
	}

	void initializeSim()
	{
		// Called at the start of the simulation or whenever the "New" button is pressed
		Time = 0;
		lastWindChangeTime = 0;
		numWindChanges = 0;
		currentImage = 0;
		setupParticles();

		currentWind = 0;

		wind = windSpeeds[currentWind];
		double angle = windAngles[currentWind];

		windSpeed.setValue((int)(wind * 100)); // Set the slider to the wind speed
		windDirection.setValue((int)(angle)); // Set the slider to the wind angle
	}


	void setupParticles() {

		Point3d coords[] = null;
		Vector3f normals[] = null;

		destroyParticles();
		destroyStagnationNormals();

		blobs = new Vector(1000); 

		nBlobs = 0;

		nPgen1_2 = 0;
		nPgen2_4 = 0;
		nPgen4_8 = 0;
		nPgen8_16 = 0;
		nPgen16_32 = 0;
		nPgen32_64 = 0;
		nPgen64_128 = 0;
		nPgen128_256 = 0;

		for(int type=0;type<nTypes;type++) {
				double x = 0;
				double y = 0;
				double z = 0;
				double Vx = 0.0;
				double Vy = 0.0;
				double Vz = 0.0;

				System.out.println("Type "+type+" over "+nTriangles+" triangles");

				double smallestArea = 100000.0;
				double totalArea = 0;

				for(int it=0;it<nTriangles;it++) {

					if (it % 100 == 0) System.out.print(".");

					Triangle T = triangle[it];

					// Ignore "triangles" from geometry that have zero area 

					if(T.area <= 0) continue;

					// Note that geometry is in arbitrary units, but surface density is per metre squared.
					double area = T.area * GEOMETRYSCALEFACTOR * GEOMETRYSCALEFACTOR;
					totalArea += area;

					if (area < smallestArea) smallestArea = area;

					// Show the stagnation points for this triangle if requested

					if (bStagnation)
					{
						Box box = new Box(-5, 5, -5, 5, -5, 5, new Color3f(1.f, 0.f, 1.f));

						Shape3D shape = box.getShape3D();

						Transform3D boxT = new Transform3D();
						boxT.set(new Vector3d(T.stagnation.x, T.stagnation.y, T.stagnation.z));

						TransformGroup boxTG = new TransformGroup(boxT);
						boxTG.setCapability(TransformGroup.ALLOW_TRANSFORM_WRITE);
						boxTG.setCapability(TransformGroup.ALLOW_TRANSFORM_READ);
						boxTG.setCapability(TransformGroup.ALLOW_BOUNDS_READ);
						boxTG.setCapability(TransformGroup.ALLOW_BOUNDS_WRITE);

						boxTG.addChild(box.getBranchGroup());

						BranchGroup boxBG = new BranchGroup();
						boxBG.setCapability(BranchGroup.ALLOW_DETACH);
						boxBG.setCapability(BranchGroup.ALLOW_CHILDREN_EXTEND);

						boxBG.addChild(boxTG);
						boxBG.setUserData(new String("Stag"));

						objTrans.addChild(boxBG);
					}

					if (bShowNormals) {
						
						LineArray normalArray = new LineArray(2, LineArray.COORDINATES);
						normalArray.setCapability(LineArray.ALLOW_COORDINATE_WRITE);
						double xx = (T.xA+T.xB+T.xC)/3;
						double yy = (T.yA+T.yB+T.yC)/3;
						double zz = (T.zA+T.zB+T.zC)/3;
						Point3d coord = new Point3d(xx,yy,zz);
						normalArray.setCoordinate(0, coord);
						xx += 100*(double) T.normal.x;
						yy += 100*(double) T.normal.y;
						zz += 100*(double) T.normal.z;
						coord = new Point3d(xx,yy,zz);
						normalArray.setCoordinate(1, coord);

						Appearance normalAppearance = new Appearance();
						LineAttributes normalAttributes = new LineAttributes(1.0f, LineAttributes.PATTERN_SOLID, false);
						normalAppearance.setLineAttributes(normalAttributes);

						ColoringAttributes normalColoring = new ColoringAttributes();
						normalColoring.setColor(new Color3f(0.0f, 1.0f, 0.f));
						normalAppearance.setColoringAttributes(normalColoring);

						Shape3D normalLine = new Shape3D(normalArray, normalAppearance);
						BranchGroup normalBG = new BranchGroup();
						normalBG.setCapability(BranchGroup.ALLOW_DETACH);
						normalBG.setUserData(new String("Norm"));

						normalBG.addChild(normalLine);

						objTrans.addChild(normalBG);
					}

					int ngenerated = 0;
					int n1_2 = 0;
					int n2_4 = 0;
					int n4_8 = 0;
					int n8_16 = 0;
					int n16_32 = 0;
					int n32_64 = 0;
					int n64_128 = 0;
					int n128_256 = 0;

					if (particleSize == 0)
					{
						// we now generate contaminant particles according the surface area, and the
						// densities calculated according to the surface cleanliness level

						n1_2 = (int)(nParticles1_2 * area);
						n2_4 = (int)(nParticles2_4 * area);
						n4_8 = (int)(nParticles4_8 * area);
						n8_16 = (int)(nParticles8_16 * area);
						n16_32 = (int)(nParticles16_32 * area);
						n32_64 = (int)(nParticles32_64 * area);
						n64_128 = (int)(nParticles64_128 * area);
						n128_256 = (int)(nParticles128_256 * area);

						ngenerated = n1_2 + n2_4 + n4_8 + n8_16 + n16_32 + n32_64 + n64_128 + n128_256;

						if (ngenerated == 0) continue; // nothing to do if no particles

						nPgen1_2 += n1_2;
						nPgen2_4 += n2_4;
						nPgen4_8 += n4_8;
						nPgen8_16 += n8_16;
						nPgen16_32 += n16_32;
						nPgen32_64 += n32_64;
						nPgen64_128 += n64_128;
						nPgen128_256 += n128_256;
					}
					else
					{
						// a fixed size of particle will be generated
						ngenerated = (int) (particleNumber * area);
						if (Surfaces.equals("TOP"))
						{
							if (T.normal.angle(new Vector3f(0.f, 0.f, 1.f)) > EPSILON) ngenerated = 0;
						}
						else if (Surfaces.equals("SIDE"))
						{
							if (T.normal.angle(new Vector3f(-1.f, 0.f, 0.f)) > EPSILON) ngenerated = 0;
						} 

					}

					if(CBoxDebug.isSelected())System.out.println("Triangle "+it+" Area "+area+" Will generate "+ngenerated+" normal "+T.normal);

					for(int ig=0;ig<ngenerated;ig++) {

						// calculate a size for this particle

						double rr = 1e-6; // base radius of particle

						if (particleSize == 0)
						{

							if (ig < n1_2)
							{
								rr *= (1 + Math.random()); // between 1 and 2 microns
							}
							else if (ig < n1_2 + n2_4)
							{
								rr *= 2 * (1 + Math.random()); // between 2 and 4 microns
							}
							else if (ig < n1_2 + n2_4 + n4_8)
							{
								rr *= 4 * (1 + Math.random()); // between 4 and 8 microns
							}
							else if (ig < n1_2 + n2_4 + n4_8 + n8_16)
							{
								rr *= 8 * (1 + Math.random()); // between 8 and 16 microns
							}
							else if (ig < n1_2 + n2_4 + n4_8 + n8_16 + n16_32)
							{
								rr *= 16 * (1 + Math.random()); // between 16 and 32 microns
							}
							else if (ig < n1_2 + n2_4 + n4_8 + n8_16 + n16_32 + n32_64)
							{
								rr *= 32 * (1 + Math.random()); // between 32 and 64 microns
							}
							else if (ig < n1_2 + n2_4 + n4_8 + n8_16 + n16_32 + n32_64 + n64_128)
							{
								rr *= 64 * (1 + Math.random()); // between 64 and 128 microns
							}
							else
							{
								rr *= 128 * (1 + Math.random()); // between 128 and 256 microns
							}
						}
						else
						{
							// user has asked for a specific particle size
							double g = gaussian(particleSize, particleSigma);
							if(g > 0) rr *= g;
						}

						// Calculate a random point position inside the triangle
						// using barycentric coordinates s,t,u

						double s = Math.random();
						double t = Math.random();
						if(s+t > 1) {
							s = 1-s;
							t = 1-t;
						}
						double u = 1 - s - t;

						x = T.xA*s + T.xB*t + T.xC*u;
						y = T.yA*s + T.yB*t + T.yC*u;
						z = T.zA*s + T.zB*t + T.zC*u;

						Point3d p = new Point3d(x,y,z);

						// mass = 4/3. Pi. R^3 * Density

						double Mass = Math.PI*4.0/3.0*rr*rr*rr*particleDensity;

						Blob b = new Blob("Surface particle "+nBlobs,type,p.x,p.y,p.z,Vx,Vy,Vz,Mass,it,rr,blobColors[type]);

						blobs.addElement(b);

						b.calculateEdges();

						nBlobs++;

						BranchGroup blobBG = b.blobBG;

						objTrans.addChild(blobBG);

					}

				}
				System.out.println("Smallest area found in geometry = " + smallestArea+ "m^2");
				System.out.println("Total Surface area of geometry = " + totalArea + "m^2");

		}

// test particles

		// a test particle
		/*
		double rr = 0.5*128. * 1.0e-6; // microns
		double Mass = Math.PI * 4.0 / 3.0 * rr * rr * rr * particleDensity;
		System.out.println("Test particle radius " + rr + " Mass " + Mass);
		double xp = 0;
		double yp = 0;
		double zp = 0.095;
		double effectiveWind = adjustWind(zp, wind);
		Blob b = new Blob("Test Particle 1",0, xp, yp, zp, 
			effectiveWind*windVector.x, effectiveWind*windVector.y, effectiveWind*windVector.z, 
			Mass, -1, rr, blobColors[0]);
		b.started = true;
		b.Xfree = xp;
		b.Yfree = yp;
		b.Zfree = zp;
		b.Tfree = Time;
		blobs.addElement(b);
		nBlobs++;
		BranchGroup blobBG = b.blobBG;
		objTrans.addChild(blobBG);
		 * 
		 * 
		 */
		
		


		nBlobs = blobs.size();

		System.out.println("Generated "+nBlobs+" particles");
		if (particleSize == 0)
		{
			System.out.println("  1-2 microns:		" + nPgen1_2);
			System.out.println("  2-4 microns:		" + nPgen2_4);
			System.out.println("  4-8 microns:		" + nPgen4_8);
			System.out.println("  8-16 microns:		" + nPgen8_16);
			System.out.println("  16-32 microns:	" + nPgen16_32);
			System.out.println("  32-64 microns:	" + nPgen32_64);
			System.out.println("  64-128 microns:	" + nPgen64_128);
			System.out.println("  128-256 microns:	" + nPgen128_256);
		}

	}

	public double gaussian(double mean, double sigma)
	{
		// Generates a random number in a gaussian distribution with
		// mean "mean" and width "sigma"
		double v1, v2, s;
		do
		{
			v1 = 2 * Math.random() - 1;   // between -1.0 and 1.0
			v2 = 2 * Math.random() - 1;   // between -1.0 and 1.0
			s = v1 * v1 + v2 * v2;
		} while (s >= 1 || s == 0);
		double multiplier = Math.sqrt(-2 * Math.log(s) / s);
		return v1 * multiplier * sigma + mean;

		// could use v2 as well, but I can't be arsed

	}

	public double LineLineDistance(Point3d p1, Point3d p2, Point3d p3, Point3d p4)
	{
		// Converted to Java from Dan Sunday's http://geometryalgorithms.com/Archive/algorithm_0106/algorithm_0106.htm

		// This calculates the distance between two 3D line segments in space

		Vector3d u = new Vector3d(p2.x - p1.x, p2.y - p1.y, p2.z - p1.z);
		Vector3d v = new Vector3d(p4.x - p3.x, p4.y - p3.y, p4.z - p3.z);
		Vector3d w = new Vector3d(p1.x - p3.x, p1.y - p3.y, p1.z - p3.z);

		double a = u.dot(u);
		double b = u.dot(v);
		double c = v.dot(v);
		double d = u.dot(w);
		double e = v.dot(w);
		double D = a * c - b * b;

		double sc, sN, sD = D;
		double tc, tN, tD = D;

		if (D < EPSILON)
		{
			sN = 0;
			sD = 1;
			tN = e;
			tD = c;
		}
		else
		{
			sN = (b * e - c * d);
			tN = (a * e - b * d);
			if (sN < 0)
			{
				sN = 0;
				tN = e;
				tD = c;
			}
			else if (sN > sD)
			{
				sN = sD;
				tN = e + b;
				tD = c;
			}
		}

		if (tN < 0)
		{
			tN = 0;
			if (-d < 0)
			{
				sN = 0;
			}
			else if (-d > a)
			{
				sN = sD;
			}
			else
			{
				sN = -d;
				sD = a;
			}
		}else if(tN > tD) {
			tN = tD;
			if((-d+b) < 0) {
				sN = 0;
			} else if((-d+b) > a) {
				sN = sD;
			} else {
				sN = (-d+b);
				sD = a;
			}
		}

		if(Math.abs(sN) < EPSILON) {
			sc = 0;
		} else {
			sc = sN/sD;
		}
		if(Math.abs(tN) < EPSILON) {
			tc = 0;
		} else {
			tc = tN/tD;
		}

		Vector3d dP = new Vector3d(w.x + sc * u.x - tc * v.x, w.y + sc * u.y - tc * v.y, w.z + sc * u.z - tc * v.z);

		//System.out.println("sc, tc " + sc + " " + tc+" w"+w+" u"+u+" v"+v);

		// psave is the position on segment p3-p4 of the closest approach of segment p1-p2
		psave = new Point3d(p3.x+tc * v.x, p3.y+tc * v.y,p3.z+tc * v.z);

		//System.out.println("Vector " + dP);
		//System.out.println("psave " + psave);
		//System.out.println("dP length " + dP.length());

		return dP.length();

	}

	public Point3d psave;

	public double OldLineLineDistance(Point3d p1, Point3d p2, Point3d p3, Point3d p4) {
		// From an algorithm by Paul Bourke at http://astronomy.swin.edu.au/~pbourke/geometry/lineline3d/ 

		// Calculate the line segment PaPb that is the shortest route between   
		// two lines P1P2 and P3P4. Calculate also the values of mua and mub where
		// Pa = P1 + mua (P2 - P1)
		// Pb = P3 + mub (P4 - P3)

		// return the length of the line segment. Negative value -1 implies lines don't cross.

		// This algorithm assumes lines are infinitely long, and that p1,p2 and p3,p4 are just
		// points on those lines, rather than the endpoints of the lines.
		// See Dan Sunday's algorithm, converted to Java, for what we need in this simulation

		Point3d p13,p43,p21,pa,pb;
		double d1343,d4321,d1321,d4343,d2121;
		double numer,denom;
		double mua,mub;

		p13 = new Point3d(p1.x - p3.x, p1.y - p3.y, p1.z - p3.z);
		p43 = new Point3d(p4.x - p3.x, p4.y - p3.y, p4.z - p3.z);

		if (Math.abs(p43.x)  < EPSILON && Math.abs(p43.y)  < EPSILON && Math.abs(p43.z)  < EPSILON) return -1;

		p21 = new Point3d(p2.x - p1.x, p2.y - p1.y, p2.z - p1.z);

		if (Math.abs(p21.x)  < EPSILON && Math.abs(p21.y)  < EPSILON && Math.abs(p21.z)  < EPSILON) return -1;

		d1343 = p13.x * p43.x + p13.y * p43.y + p13.z * p43.z;
		d4321 = p43.x * p21.x + p43.y * p21.y + p43.z * p21.z;
		d1321 = p13.x * p21.x + p13.y * p21.y + p13.z * p21.z;
		d4343 = p43.x * p43.x + p43.y * p43.y + p43.z * p43.z;
		d2121 = p21.x * p21.x + p21.y * p21.y + p21.z * p21.z;

		denom = d2121 * d4343 - d4321 * d4321;
		
		if (Math.abs(denom) < EPSILON) return -1;
		
		numer = d1343 * d4321 - d1321 * d4343;

		mua = numer / denom;
		mub = (d1343 + d4321 * mua) / d4343;

		pa = new Point3d(p1.x + mua * p21.x, p1.y + mua * p21.y, p1.z + mua * p21.z);
		pb = new Point3d(p3.x + mub * p43.x, p3.y + mub * p43.y, p3.z + mub * p43.z);

		// save point a
		psave = new Point3d(pa);

		System.out.println("Intersection " + pa+ " " + pb);
		return pa.distance(pb);
}

	public Point3d GetIntersection(Point3d p1, Point3d p2, Triangle T) {
		// For the line defined by points p1 and p2, look at each side of
		// the triangle T and see if the line intersects it. If the line
		// intersects then the position of the intersection is returned
		// If the line intersects two sides, then the position of the
		// intersection closest to p2 is returned. If no side is intersected,
		// the routine returns null

		double d;
		double distance = 100000;
		Point3d pIntersect = null;

		d = LineLineDistance(p1, p2, T.A, T.B);
		// psave is the intersection point calculated in LineLineDistance
		if(d >= 0 && d < BIGEPSILON) 
		{
			distance = p2.distance(psave);
			pIntersect = new Point3d(psave);
			if(CBoxDebug.isSelected())System.out.println("Distance of " + p1 + " " + p2 + " from triangle side " + T.A + " " + T.B + " is " + d);
		} 
		d = LineLineDistance(p1, p2, T.B, T.C);
		if (d >= 0 && d < BIGEPSILON)
		{
			if (p2.distance(psave) < distance)
			{
				pIntersect = new Point3d(psave);
				distance = p2.distance(psave);
			}
			if(CBoxDebug.isSelected())System.out.println("Distance of " + p1 + " " + p2 + " from triangle side " + T.B + " " + T.C + " is " + d);
		}
		d = LineLineDistance(p1, p2, T.C, T.A);
		if (d >= 0 && d < BIGEPSILON)
		{
			if (p2.distance(psave) < distance)
			{
				pIntersect = new Point3d(psave);
				distance = p2.distance(psave);
			}
			if(CBoxDebug.isSelected())System.out.println("Distance of " + p1 + " " + p2 + " from triangle side " + T.C + " " + T.A + " is " + d);
		}
		if (CBoxDebug.isSelected() && pIntersect == null) System.out.println("GetIntersection: no intersect!");

		return pIntersect;
	}

	public boolean SameSide(Point3d p1, Point3d p2, Point3d A, Point3d B) {
		// Checks if points p1 and p2 are both on the same side of line AB

		Vector3d BA = new Vector3d(B.x-A.x, B.y-A.y, B.z-A.z);

		Vector3d p1A = new Vector3d(p1.x-A.x, p1.y-A.y, p1.z-A.z);

		Vector3d p2A = new Vector3d(p2.x-A.x, p2.y-A.y, p2.z-A.z);

		p1A.cross(BA,p1A);
		p2A.cross(BA,p2A);

		if(p1A.dot(p2A) >= 0) return true;
		return false;
	}

	public boolean CoplanarTriangles(Triangle T1, Triangle T2)
	{
		// returns true if Triangle T1 is coplanar with triangle T2

		// The method used is to compare first the normals of each triangle. If the
		// normals are different, then the triangles cannot be coplanar.
		// Then a vector between a vertex on one triangle to a vertex on
		// the other triangle is created: if this is perpendicular to the
		// normal, then the triangles are coplanar

		if (T1.normal.angle(T2.normal) != 0) return false;
		Vector3f v = new Vector3f((float)(T1.A.x - T2.A.x), (float)(T1.A.y - T2.A.y), (float)(T1.A.z - T2.A.z));
		if (Math.abs(Math.abs(v.angle(T1.normal)) - Math.PI * 0.5) > EPSILON) return false;

		return true;
	}

	public boolean PointInTriangle(Point3d p, Triangle T) {
		// Checks if point p is inside triangle T
		if(!SameSide(p,T.A,T.B,T.C)) return false;
		if(!SameSide(p,T.B,T.A,T.C)) return false;
		if(!SameSide(p,T.C,T.A,T.B)) return false;
		return true;
	}


   public String getAppletInfo () {
      return "Lander by J.J.B.";
   }




	double adjustWind(double windSpeed, double Z)
	{
		// For a nominal wind speed of "windSpeed" above the boundary
		// layer, calculate the effective wind at a lower position


		// From Greeley/Iversen paper u/ustar = ln((z+z0)/z0)/0.4
		// 0.4(1) is von Karman's constant
		// See also http://www.sciencedirect.com/science?_ob=ArticleURL&_udi=B6V93-49Y9C9M-9&_user=1010281&_coverDate=04%2F02%2F2004&_rdoc=1&_fmt=&_orig=search&_sort=d&view=c&_acct=C000050264&_version=1&_urlVersion=0&_userid=1010281&md5=ec596b56e4d9d980a21ee1f65ff76fd7#sec1

		//double zeff = Z - lowestZ;
		//double zzero = 1.97*(highestZ-lowestZ);

		//double effectiveWind = wind * Math.log((zeff + zzero) / zzero) / 0.41;

		//System.out.println(" zeff, zzero, wind, eff " + zeff + " " + zzero + " " + wind + " " + effectiveWind);

		// Linear drop from highestZ to lowestZ
		double effectiveWind = wind * (Z - lowestZ) / (highestZ - lowestZ);
		if (effectiveWind > wind) effectiveWind = wind;


		return effectiveWind;
	}

	TriangleArray shapeGetTriangles(Shape3D shape) {
		// Finds the TriangleArray geometry in the Shape3D, else returns null
		int numGeometries = shape.numGeometries();
		//System.out.println("Shape3D "+shape+" "+numGeometries+" geometries");
		for(int i=0;i<numGeometries;i++) {
			Geometry g = shape.getGeometry(i);
			if(g instanceof GeometryArray) {
				if(g instanceof TriangleArray) {
					// This is what we deal with
					return (TriangleArray) g;
				} else {
					System.out.print("Untreated ");
					if(g instanceof QuadArray) System.out.print("QuadArray");
					if(g instanceof PointArray) System.out.print("PointArray");
					if(g instanceof LineArray) System.out.print("LineArray");
					if(g instanceof IndexedGeometryArray) System.out.print("IndexedGeometryArray");
					if(g instanceof LineStripArray) System.out.print("LineStripArray");
					if(g instanceof TriangleFanArray) System.out.print("TriangleFanArray");
					if(g instanceof TriangleStripArray) System.out.print("TriangleStripArray");
					System.out.println("");
					System.exit(0);

				}
			} else {
				System.out.println("Untreated geometry");
			}
		}
		return null;
	}

	public static double erf(double z)
	{
		// From MyMath.java
		double t = 1.0 / (1.0 + 0.5 * Math.abs(z));

		// use Horner's method
		double ans = 1 - t * Math.exp(-z * z - 1.26551223 +
											t * (1.00002368 +
											t * (0.37409196 +
											t * (0.09678418 +
											t * (-0.18628806 +
											t * (0.27886807 +
											t * (-1.13520398 +
											t * (1.48851587 +
											t * (-0.82215223 +
											t * (0.17087277))))))))));
		if (z >= 0) return ans;
		else return -ans;
	}


	class Triangle {
		// a triangle in the model.
		Vector3f normal;
		double xA,yA,zA;
		double xB,yB,zB;
		double xC,yC,zC;

		Point3d A,B,C;

		Point3d centroid;

		Point3d stagnation; // Position of the stagnation point for the surface on which this triangle lies

		double area;

		int neighbour[];
		int nn;

		Triangle(Point3d v1, Point3d v2, Point3d v3, Vector3f _normal) {

			xA = v1.x;
			yA = v1.y;
			zA = v1.z;
			xB = v2.x;
			yB = v2.y;
			zB = v2.z;
			xC = v3.x;
			yC = v3.y;
			zC = v3.z;

			A = new Point3d(v1);
			B = new Point3d(v2);
			C = new Point3d(v3);

			// Triangle defined by two vectors, AB, AC
			// and area is half cross product

			normal = null;
			if (_normal != null)
			{
				normal = new Vector3f(_normal);
				normal.normalize();
			}

			double xAB = xB-xA;
			double yAB = yB-yA;
			double zAB = zB-zA;

			double xAC = xC-xA;
			double yAC = yC-yA;
			double zAC = zC-zA;

			double crossx = yAB*zAC - zAB*yAC;
			double crossy = -(xAB*zAC - zAB*xAC);
			double crossz = xAB*yAC - yAB*xAC;
			
			// calculate centroid (centre of mass) of the triangle

			centroid = new Point3d((xA + xB + xC) / 3, (yA + yB + yC) / 3, (zA + zB + zC) / 3);

			// calculate area of the triangle

			area = 0.5* Math.sqrt( crossx*crossx + crossy*crossy + crossz*crossz );

			if(CBoxDebug.isSelected()) {
				System.out.println("New triangle Area "+area);
				System.out.println("Centroid " + centroid);
				System.out.println("V1: "+xA+" "+yA+" "+zA);
				System.out.println("V2: "+xB+" "+yB+" "+zB);
				System.out.println("V3: "+xC+" "+yC+" "+zC+"\n");
			}

		}

		String print()
		{
			return new String("(" + (float)xA + "," + (float)yA + "," + (float)zA + "),("
				+ (float)xB + "," + (float)yB + "," + (float)zB + "),("
				+ (float)xC + "," + (float)yC + "," + (float)zC + ")");
		}

		void setNeighbours(int _nn, int _n[]) {
			nn = _nn;
			neighbour = new int[nn];
			for(int i=0;i<nn;i++) {
				neighbour[i] = _n[i];
			}
		}

		void setStagnation(Point3d CofM)
		{
			stagnation = new Point3d(CofM);
		}
		
	}


   class Blob {

	   // current position, velocity, mass, size
		double Xpos;	// centimetres
		double Ypos;
		double Zpos;

	   double Xfree;	// position when particle leaves surface, centimetres
	   double Yfree;
	   double Zfree;
	   double Tfree;	// time at which particle leaves surface, seconds

		double Vx;		// metres/sec
		double Vy;
		double Vz;

		double Mass;	// kilograms
		double radius;	// radius of particle

	   double maxDF;	// maximum value of the detachment fraction seen by this blob

		double S;		// Ratio of the particle's density to that of Air or the planetary atmosphere
		double T;		// T = (1/18)*d^2*S*C/Viscosity

	   double FoverTau; // 6*Pi*R^3/(Wa.N.Z0.sqrt(2.R.alpha.Z0))
						// i.e. this is F/tau

	   double detachmentFraction;

		double size;	// size of plotted point in display (not physically relevant)

		LineArray windLineArray;
	   LineArray normalArray;
		Shape3D windLine;
	   BranchGroup windLineBG;

		int itriangle; // index of the Triangle that this Blob is on

		// transform in which the coordinates are set
		Transform3D    blobMat;

	   double windAngle; // Angle the wind makes with the surface the particle is on
	   Vector3f rollDirection; // the direction in which the particle will roll
	   Point3d pLeadingEdge; // the position of the leading edge intersection into the current wind
	   Point3d pDestinationPoint; // the position of the trailing edge intersection away from the current wind

	   double dStag; // The ratio of the distance of the blob to the stagnation point, and the distance of
					 // the destination point from the stagnation point

	   double dDestination; // The distance the particle mist travel to its destination edge on the surface

		int type;
		Color3f color;
		String name;

		boolean landed;
		boolean started;

		BranchGroup    blobBG;
		TransformGroup blobTG;

		Blob (String _name, int _type, double X, double Y, double Z, 
			double Vx0, double Vy0, double Vz0, double M0, int _itriangle, double _radius, Color3f _color) {
			Xpos = X;
			Ypos = Y;
			Zpos = Z;
			Vx = Vx0;
			Vy = Vy0;
			Vz = Vz0;
			Mass = M0;
			type = _type;
			name = new String(blobTypes[type]);
			radius = _radius;
			color = _color;
			itriangle = _itriangle;

			name = new String(_name);

			landed = false;
			started = false;
			pLeadingEdge = null;
			pDestinationPoint = null;
			rollDirection = null;

			S = particleDensity / Density;
			T = 4 * radius * radius * S * CC / (18 * Viscosity);

			//System.out.println("S = "+S+" CC = " + CC+ " Viscosity=" +Viscosity);

			double WA = 0.012;  // Newtons/metre. For wood or Lycopodium
			double Nlyco = 300; // A constant for Lycopodium
			double Zzero = 4.0e-10; // inter atomic distance for Lycopodium
			double Alpha = 20;	// Another constant

			FoverTau = 6 * Math.PI * radius * radius * radius;
			FoverTau /= (WA * Nlyco * Zzero * Math.sqrt(2 * radius * Alpha * Zzero));

			maxDF = 0.0;

			//System.out.println("FoverTau = " + FoverTau);

			size = Math.abs((highestZ-lowestZ) / 200) ; // for display purposes only

			Box box = new Box(-size/2,size/2,-size/2,size/2,-size/2,size/2,null);

			Shape3D shape = box.getShape3D();

			shape.setUserData(new String("B"));

			blobMat = new Transform3D();
			blobMat.set(new Vector3d(Xpos,Ypos,Zpos));

			blobTG = new TransformGroup(blobMat);
			blobTG.setCapability(TransformGroup.ALLOW_TRANSFORM_WRITE);
			blobTG.setCapability(TransformGroup.ALLOW_TRANSFORM_READ );
			blobTG.setCapability(TransformGroup.ALLOW_BOUNDS_READ);
			blobTG.setCapability(TransformGroup.ALLOW_BOUNDS_WRITE);

			blobTG.addChild(box.getBranchGroup());

			blobBG = new BranchGroup();
			blobBG.setCapability(BranchGroup.ALLOW_DETACH);
			blobBG.setCapability(BranchGroup.ALLOW_CHILDREN_EXTEND);

			blobBG.addChild(blobTG);

			if(CBoxDebug.isSelected()) System.out.println("Created "+name+" type "+blobTypes[_type]+
				" at ("+Xpos+","+Ypos+","+Zpos+
				") (cm.) Radius "+radius+" (m) Mass "+Mass+" (kg) S "+S+" T "+T+" (s)");

			if (CBoxEdges.isSelected())
			{
				windLineArray = new LineArray(2, LineArray.COORDINATES);
				windLineArray.setCapability(LineArray.ALLOW_COORDINATE_WRITE);
				Point3d coord = new Point3d(Xpos, Ypos, Zpos);
				windLineArray.setCoordinate(0, coord);
				coord = new Point3d(Xpos, Ypos, Zpos);
				windLineArray.setCoordinate(1, coord);

				Appearance windLineAppearance = new Appearance();
				LineAttributes windLineAttributes = new LineAttributes(1.0f, LineAttributes.PATTERN_SOLID, false);
				windLineAppearance.setLineAttributes(windLineAttributes);

				ColoringAttributes windLineColoring = new ColoringAttributes();
				windLineColoring.setColor(new Color3f(1.0f, 0.f, 0.f));
				windLineAppearance.setColoringAttributes(windLineColoring);

				windLine = new Shape3D(windLineArray, windLineAppearance);
				windLineBG = new BranchGroup();
				windLineBG.setCapability(BranchGroup.ALLOW_DETACH);
				windLineBG.addChild(windLine);

				blobBG.addChild(windLineBG);
			}


		}

	   void calculateEdges() {

		   // Check the current wind direction compared with the normal to the
		   // surface the particle is sitting on.

		   Point3d p = new Point3d(Xpos, Ypos, Zpos); // current position

		   // If the particle is already moving, then we do not calculate the leading/trailing edge,
		   // as that was already done, and the assumption is that the wind does not change once the
		   // particle starts moving

		   if (started) return;

		   // We do not calculate for landed particles

		   if (landed) return;

		   pLeadingEdge = null;
		   pDestinationPoint = null;

		   if (CBoxEdges.isSelected())
		   {
			   // the windlinearray is just a visual line for the display that connects the particle and
			   // the computed trailing edge of the lander volume

			   windLineArray.setCoordinate(0, p);
			   windLineArray.setCoordinate(1, p);
		   }

		   // If the particle is not on a triangle (not on a surface), we do not compute leading/trailing edge
		   if (itriangle == -1) return;

		   Vector3f normal = triangle[itriangle].normal;

		   windAngle = windVector.angle(normal);

		   // Is the wind vector bearing on the particle? 


		   if (windAngle < 0.5 * Math.PI - EPSILON) return; // surface is away from the wind, and not parallel to it


		   if(CBoxDebug.isSelected()) System.out.println("Surface for this particle is at " + windAngle + " radians to the wind vector");

		   Triangle TN = triangle[itriangle];

		   if (windAngle > 0.5 * Math.PI + EPSILON)
		   {
			   // wind bears on surface and is not parallel
			   // This is where we calculate the roll direction based on the stagnation point for the
			   // wind on the triangle
			   //System.out.println("Surface is in wind");


			   // for acute angles, roll direction is computed from stagnation point, otherwise it is
			   // computed from the wind vector directly

			   if (windAngle > 0.75 * Math.PI)
			   {
				   //System.out.println("Compute using stagnation point");
				   rollDirection = new Vector3f((float)(Xpos - TN.stagnation.x), (float)(Ypos - TN.stagnation.y), (float)(Zpos - TN.stagnation.z));
				   rollDirection.normalize();
			   }
			   else
			   {
				   //System.out.println("Compute using wind component directly");
					rollDirection = new Vector3f();
				   rollDirection.cross(windVector, normal);
				   rollDirection.cross(normal, rollDirection);
			   }


			   // Create a long vector in the roll direction and check where it intersects the furthest edge
			   // of the surface we are on

			   Point3d infinity = new Point3d(Xpos + 10000.f * rollDirection.x, Ypos + 10000.f * rollDirection.y, Zpos + 10000.f * rollDirection.z);

			   if (CBoxEdges.isSelected()) windLineArray.setCoordinate(1, infinity);

			   double shortest = p.distance(infinity);

			   Triangle T = TN;

			   if(CBoxDebug.isSelected()) System.out.println("\n\nChecking " + name + " for furthest wind edge (" + Xpos + " " + Ypos + " " + Zpos + ")");

			   for (int ineighbour = -1; ineighbour < T.nn; ineighbour++)
			   {
				   if(CBoxDebug.isSelected()) System.out.println("\nNeighbour " + ineighbour);
				   if (ineighbour >= 0) TN = triangle[T.neighbour[ineighbour]];

				   Point3d pIntersect = GetIntersection(p, infinity, TN);

				   if (pIntersect != null)
				   {
					   if (infinity.distance(pIntersect) < shortest)
					   {
						   shortest = infinity.distance(pIntersect);
						   pDestinationPoint = new Point3d(pIntersect);
						   if (CBoxEdges.isSelected()) windLineArray.setCoordinate(1, pDestinationPoint);
						   if(CBoxDebug.isSelected()) System.out.println("Trailing Edge " + pDestinationPoint + " Shortest " + shortest);
					   }
				   }
			   }

			   dStag = 0;

			   if (CBoxDebug.isSelected()) System.out.println("Destination point: " + pDestinationPoint);
			   if (pDestinationPoint == null)
			   {
				   System.err.println("Destination point is NULL!");
				   return;
			   }

			   // Compute dStag, the ratio of the distance of the point from the stagnation point, to the distance
			   // of the destination point from the stagnation point. This will be used to set the rolling velocity
			   // when the particle moves
			   TN = triangle[itriangle];
			   dStag = TN.stagnation.distance(p) / TN.stagnation.distance(pDestinationPoint);

			   //System.out.println("dStag = " + dStag);

			   double effectiveWind = adjustWind(wind,Zpos);

			   double U = effectiveWind;

			   dDestination = p.distance(pDestinationPoint);

			   double dEdge = 0;

			   if (windAngle > 0.75 * Math.PI)
			   {
					// Compute the distance of the particle from the destination edge
					dEdge = dDestination;

				   // calculate the effective wind speed
				   U *= dStag;
				   if (Math.abs(U) < EPSILON) U = EPSILON;

			   } else {
				   // we need to calculate the distance of the particle/blob to the leading edge of the surface in the wind
				   // Create a long vector that points from the particle into the wind,
				   // and check where it intersects with each of the object's triangles,
				   // picking the intersection point closest to the wind.

				   if (Math.abs(U) < EPSILON) U = EPSILON;

				   infinity = new Point3d(Xpos - 10000.f * rollDirection.x, Ypos - 10000.f * rollDirection.y, Zpos - 10000.f * rollDirection.z);
				   //if(CBoxEdges.isSelected()) windLineArray.setCoordinate(1, infinity);

				   shortest = p.distance(infinity);

				   T = TN;

				   if (CBoxDebug.isSelected()) System.out.println("\n\nChecking " + name + " for nearest wind edge (" + Xpos + " " + Ypos + " " + Zpos + ")");

				   for (int ineighbour = -1; ineighbour < T.nn; ineighbour++)
				   {
					   if (CBoxDebug.isSelected()) System.out.println("\nNeighbour " + ineighbour);
					   if (ineighbour >= 0) TN = triangle[T.neighbour[ineighbour]];

					   Point3d pIntersect = GetIntersection(p, infinity, TN);

					   if (pIntersect != null)
					   {
						   if (infinity.distance(pIntersect) < shortest)
						   {
							   shortest = infinity.distance(pIntersect);
							   pLeadingEdge = new Point3d(pIntersect);
							   //if(CBoxEdges.isSelected()) windLineArray.setCoordinate(1, pLeadingEdge);
							   if (CBoxDebug.isSelected()) System.out.println("Leading Edge " + pLeadingEdge + " Shortest " + shortest);
						   }
					   }
				   }

				   if (pLeadingEdge == null)
				   {
					   System.out.println("Particle in Acute wind has no leading edge! " + p + " Infinity = " + infinity);
					   return;
				   }
				   // Compute the distance of the particle from the leading edge
				   dEdge = p.distance(pLeadingEdge);
			   }

			   // calculate the boundary layer thickness (note dEdge is scaled from the VRML model coordinates)
			   double delta = Math.sqrt(Viscosity * dEdge * GEOMETRYSCALEFACTOR / U);

			   //System.out.println("Acute wind particle, delta = " + delta);

			   double tau = Density * Viscosity * U / delta;
			   double F = FoverTau * tau;

			   detachmentFraction = 0;
			   if (F >= 3)
			   {
				   double sigmaF = 70. / 3.; // F50 / 3

				   double a = (2.5 - F) / sigmaF / 1.4142135623;

				   detachmentFraction = -erf(a);

				   // The particle will not lift off if it has already been subjected to
				   // conditions that created a larger detachmentFraction than the current value

				   if (detachmentFraction > maxDF)
				   {
					   maxDF = detachmentFraction;
				   }
				   else
				   {
					   detachmentFraction = 0;
				   }
			   }
			   //System.out.println("Detachment Frac " + detachmentFraction);

			   if (CBoxDebug.isSelected()) System.out.println("Effective Wind = " + effectiveWind + " particle radius = " + radius + " tau=" + tau + " F=" + F + " detachment fraction = " + detachmentFraction);

			   // Use the detachment fraction to determine if the particle moves
			   if (Math.random() < detachmentFraction || !bLeadingEdge)
			   {
				   if (CBoxDebug.isSelected()) System.out.println("Particle will detach");
				   // We create a non-null point that will only be used as a check
				   pLeadingEdge = new Point3d(1, 1, 1);
			   }
			   else
			   {
				   if (CBoxDebug.isSelected()) System.out.println("Particle will NOT detach");
				   pLeadingEdge = null;
			   }
			   return;
		   }



		   // special case for when wind is parallel to surface 
		   // we need to calculate the perpendicular distance of the point from the leading edge of the wind

		   rollDirection = new Vector3f();

		   rollDirection.cross(windVector, normal);
		   rollDirection.cross(normal, rollDirection);

		   if (windAngle < 0.5 * Math.PI + EPSILON)
		   {
			   // Create a long vector that points from the particle into the wind,
			   // and check where it intersects with each of the object's triangles,
			   // picking the intersection point closest to the wind.

			   Point3d infinity = new Point3d(Xpos - 10000.f * windVector.x, Ypos - 10000.f * windVector.y, Zpos - 10000.f * windVector.z);
			   //if(CBoxEdges.isSelected()) windLineArray.setCoordinate(1, infinity);

			   double shortest = p.distance(infinity);

			   Triangle T = TN;

			   if(CBoxDebug.isSelected()) System.out.println("\n\nChecking " + name + " for nearest wind edge (" + Xpos + " " + Ypos + " " + Zpos + ")");

			   for (int ineighbour = -1; ineighbour < T.nn; ineighbour++)
			   {
				   if(CBoxDebug.isSelected()) System.out.println("\nNeighbour " + ineighbour);
				   if (ineighbour >= 0) TN = triangle[T.neighbour[ineighbour]];

				   Point3d pIntersect = GetIntersection(p, infinity, TN);

				   if (pIntersect != null)
				   {
					   if (infinity.distance(pIntersect) < shortest)
					   {
						   shortest = infinity.distance(pIntersect);
						   pLeadingEdge = new Point3d(pIntersect);
						   //if(CBoxEdges.isSelected()) windLineArray.setCoordinate(1, pLeadingEdge);
						   if(CBoxDebug.isSelected()) System.out.println("Leading Edge " + pLeadingEdge + " Shortest " + shortest);
					   }
				   }
			   }

			   if (pLeadingEdge == null)
			   {
				   System.out.println("Particle in wind has no leading edge! " + p + " Infinity = " + infinity);
				   return;
			   }

			   // We have the distance to the leading edge

			   double U = adjustWind(wind, Zpos);
			   if (Math.abs(U) < EPSILON) U = EPSILON;


			   double X = GEOMETRYSCALEFACTOR * p.distance(pLeadingEdge); // distance of particle to leading edge in metres
			   double Re = X * U / Viscosity;

			   if (Math.abs(Re) < EPSILON) Re = EPSILON;

			   double Cf = 0.0583 / Math.pow(Re, 0.2);
			   double tau = 0.5 * Density * U * U * Cf;
			   double F = FoverTau * tau;

			   detachmentFraction = 0;
			   if (F >= 3)
			   {
				   double sigmaF = 70. / 3.; // F50 / 3

				   double a = (2.5 - F) / sigmaF / 1.4142135623;

				   detachmentFraction = -erf(a);
				   // The particle will not lift off if it has already been subjected to
				   // conditions that created a larger detachmentFraction than the current value

				   if (detachmentFraction > maxDF)
				   {
					   maxDF = detachmentFraction;
				   }
				   else
				   {
					   detachmentFraction = 0;
				   }

			   }

			   if (CBoxDebug.isSelected()) System.out.println("Effective Wind = " + U + " particle radius = " + radius + " X= " + X + " tau=" + tau + " Re=" + Re + " Cf=" + Cf + " F=" + F + " detachment fraction = " + detachmentFraction);

				// Use the detachment fraction to determine if the particle moves

			   if (Math.random() < detachmentFraction || !bLeadingEdge)
			   {
				   if (CBoxDebug.isSelected()) System.out.println("Particle will detach");
			   }
			   else
			   {
				   if (CBoxDebug.isSelected()) System.out.println("Particle will NOT detach");
				   pLeadingEdge = null;
			   }

			   // Now we do the same thing, but for the trailing edge. This is where the particle will roll to.

			   infinity = new Point3d(Xpos + 10000.f * windVector.x, Ypos + 10000.f * windVector.y, Zpos + 10000.f * windVector.z);
			   if(CBoxEdges.isSelected()) windLineArray.setCoordinate(1, infinity);

			   shortest = p.distance(infinity);

			   TN = triangle[itriangle];
			   T = TN;

			   if(CBoxDebug.isSelected()) System.out.println("\n\nChecking " + name + " for furthest wind edge (" + Xpos + " " + Ypos + " " + Zpos + ")");

			   for (int ineighbour = -1; ineighbour < T.nn; ineighbour++)
			   {
				   if(CBoxDebug.isSelected()) System.out.println("\nNeighbour " + ineighbour);
				   if (ineighbour >= 0) TN = triangle[T.neighbour[ineighbour]];

				   Point3d pIntersect = GetIntersection(p, infinity, TN);

				   if (pIntersect != null)
				   {
					   if (infinity.distance(pIntersect) < shortest)
					   {
						   shortest = infinity.distance(pIntersect);
						   pDestinationPoint = new Point3d(pIntersect);
						   if(CBoxEdges.isSelected()) windLineArray.setCoordinate(1, pDestinationPoint);
						   if(CBoxDebug.isSelected()) System.out.println("Trailing Edge " + pDestinationPoint + " Shortest " + shortest);
					   }
				   }
			   }
			   if (pDestinationPoint == null)
			   {
				   System.out.println("Particle in wind has no destination point! " + p + " Infinity = " + infinity);
				   return;
			   }

			   dDestination = p.distance(pDestinationPoint);

		   }
	   }

		void update() {

			if (landed) return;

			Point3d p = new Point3d(Xpos, Ypos, Zpos); // current position

			if (!started)
			{

				// Is the wind vector bearing on the particle? If not, return

				if (windAngle < 0.5 * Math.PI - EPSILON) return;

				// Wind is bearing ... check the regime

				double effectiveWind = adjustWind(wind, Zpos);

				if (windAngle > 0.5 * Math.PI + EPSILON)
				{
					// wind is at an acute angle
					if (windAngle > 0.75 * Math.PI)
					{
						// Check if particle will detach: signified by non null pLeadingEdge

						if (pLeadingEdge == null) return;

						if (itriangle == -1) return; // should never happen!
						// Use distance from stagnation point to set the particle's velocity
						Vx = dStag * effectiveWind * rollDirection.x;
						Vy = dStag * effectiveWind * rollDirection.y;
						Vz = dStag * effectiveWind * rollDirection.z;

						//System.out.println("Particle will roll to edge with dStag " + dStag + " " + Vx + " " + Vy + " " + Vz);

						started = true;


					}
					else
					{
						// Check if particle will detach: signified by non null pLeadingEdge

						if (pLeadingEdge == null) return;

						// Just use the component of the wind's vector to set the particle's velocity
						Vx = effectiveWind * rollDirection.x;
						Vy = effectiveWind * rollDirection.y;
						Vz = effectiveWind * rollDirection.z;

						//System.out.println("Particle will roll to edge " + Vx + " " + Vy + " " + Vz);

						started = true;
					}


				}
				else
				{

					// wind is parallel 

					if (pLeadingEdge == null) return;

					// pLeadingEdge is not null if particle should detach

					// Particle will start rolling to trailing edge
					Vx = effectiveWind * rollDirection.x;
					Vy = effectiveWind * rollDirection.y;
					Vz = effectiveWind * rollDirection.z;

					started = true;
				}
			}
			if (started && !landed && itriangle != -1)
			{

				// The particle is rolling. Instead of simulating it rolling to the edge, we
				// can place it directly at the edge, in essentially zero time.

				double particleSpeed = Math.sqrt(Vx * Vx + Vy * Vy + Vz * Vz);
				double stepDistance = particleSpeed * dtime;

				//if (pDestinationPoint != null)
				//{
				//	System.out.println("Distance " + p.distance(pDestinationPoint) + " tolerance " + stepDistance);
				//}

				if (pDestinationPoint == null)
				{
					// This should not happen: needs to be investigated
					System.err.println("Destination point is null, setting to current position");
					Xfree = Xpos;
					Yfree = Ypos;
					Zfree = Zpos;
					Tfree = Time;
					if (CBoxDebug.isSelected()) System.out.println("Particle " + name + " leaves surface, Vx,Vy,Vz "
						 + Vx + " " + Vy + " " + Vz +
						 " X,Y,Z " + Xpos + " " + Ypos + " " + Zpos + " (cm) time " + Time);
					itriangle = -1;
				}
				else if (!CBoxShowMotion.isSelected())
				{
					// Place particle directly at edge
					Xfree = Xpos = pDestinationPoint.x;
					Yfree = Ypos = pDestinationPoint.y;
					Zfree = Zpos = pDestinationPoint.z;
					Tfree = Time;
					if (CBoxDebug.isSelected()) System.out.println("Particle " + name + " leaves surface, Vx,Vy,Vz "
						 + Vx + " " + Vy + " " + Vz +
						 " X,Y,Z " + Xpos + " " + Ypos + " " + Zpos + " (cm) time " + Time);
					itriangle = -1;
				}
				else if (p.distance(pDestinationPoint) > dDestination)
				{
					// The particle is further from its destination than it was in the previous step: we
					// take it that it is now released as it has reached its destination
					if (CBoxDebug.isSelected()) System.out.println("Particle " + name + " leaves surface, Vx,Vy,Vz "
						 + Vx + " " + Vy + " " + Vz +
						 " X,Y,Z " + Xpos + " " + Ypos + " " + Zpos + " (cm) time " + Time);
					Xfree = Xpos = pDestinationPoint.x;
					Yfree = Ypos = pDestinationPoint.y;
					Zfree = Zpos = pDestinationPoint.z;
					Tfree = Time;
					itriangle = -1;
				}
				else
				{
					// we compute the distance the particle must go before reaching its destination on the edge
					dDestination = p.distance(pDestinationPoint);
				}
			}

			if (started && !landed && !CBoxShowMotion.isSelected())
				{
					// If we do not simulate the free fall motion of the released particle,
					// calculate the intercept with the surface directly, and place the particle there
					//System.out.println("Particle will be directly moved to end position on surface");

					//System.out.println("wind = "+(wind*windVector.x)+" T= "+T+" alpha = "+T/(wind*windVector.x));
					if (bViscousDrag)
					{  
						double effectiveWind = adjustWind(wind, Zpos);
						//System.out.println("Xpos,Ypos,Zpos,lowestZ,T,Gravity,wind,effectivewind,Gravity*T " +Xpos+" "+Ypos+" "+ Zpos + " " + lowestZ + " " + T + " " + Gravity+" "+wind+" "+effectiveWind+" "+(Gravity*T));
						double dx = (Zpos - lowestZ) * (effectiveWind * windVector.x) / (2 * Gravity * T);
						double dy = (Zpos - lowestZ) * (effectiveWind * windVector.y) / (2 * Gravity * T);
						Xpos += dx / GEOMETRYSCALEFACTOR;
						Ypos += dy / GEOMETRYSCALEFACTOR;
						Zpos = lowestZ;
						//System.out.println("End Xpos,Ypos,Zpos " + Xpos + " " + Ypos + " " + Zpos);

					}
					else
					{
						// Calculate the final velocity in z 
						double v = Math.sqrt(Vz * Vz + 2 * Gravity * (Zpos - lowestZ) * GEOMETRYSCALEFACTOR);

						//System.out.println(" v = " + v);
						// and thus the time taken to fall
						double t = (v - Vz) / Gravity;
						//System.out.println(" t = " + t);
						// which gives the final X and Y positions (in fact accelx and accely are zero here)

						double effectiveWind = adjustWind(wind, Zpos);

						Xpos += effectiveWind * windVector.x * t / GEOMETRYSCALEFACTOR;
						Ypos += effectiveWind * windVector.y * t / GEOMETRYSCALEFACTOR;
						Zpos = lowestZ;

					}
					landed = true;
			}

			if(Zpos <= lowestZ) {
				if (false)
				{
					System.out.println("\nParticle " + name + " lands, Vx,Vy,Vz "
					+ Vx + " " + Vy + " " + Vz +
					" X,Y,Z " + Xpos + " " + Ypos + " " + Zpos + " (cm) time " + Time);
					System.out.println(" Fell " + GEOMETRYSCALEFACTOR * (Zfree - Zpos) + " metres in " + (Time - Tfree) + " seconds");
				}
				blobMat.set(new Vector3d(Xpos, Ypos, Zpos));
				blobTG.setTransform(blobMat);
				landed = true;
				return;
			}

			// To arrive here, particle must a) have started moving, and b) not landed


			if (itriangle == -1) {
				// free particles

				double effectiveWind = adjustWind(wind, Zpos);

				if (bViscousDrag)
				{ // we add terms for the viscous drag 
					Vx = effectiveWind * windVector.x;
					Vy = effectiveWind * windVector.y;
					Vz = -Gravity*T;

				}
				else
				{
					Vx = effectiveWind * windVector.x;
					Vy = effectiveWind * windVector.y;
					Vz += -Gravity*dtime;
				}

			}

			if (CBoxDebug.isSelected()) System.out.println("Particle " + name + " in motion, wind "+wind+" Vx,Vy,Vz "
				 + Vx + " " + Vy + " " + Vz +
				 " X,Y,Z " + Xpos + " " + Ypos + " " + Zpos + " (cm) time " + Time + " Triangle "+itriangle+" Destination "+pDestinationPoint);
			// Calculate the change in position (metres)

			double dx = (Vx*dtime);// 
			double dy = (Vy*dtime);// 
			double dz = (Vz*dtime);//

			// Update the position 

			Xpos += dx / GEOMETRYSCALEFACTOR;
			Ypos += dy / GEOMETRYSCALEFACTOR;
			Zpos += dz / GEOMETRYSCALEFACTOR;

			blobMat.set(new Vector3d(Xpos, Ypos, Zpos));
			blobTG.setTransform(blobMat);

		}
		
	}




	class FrameBehavior extends Behavior {

	    WakeupOnElapsedFrames w = new WakeupOnElapsedFrames(0);

		public void initialize() {
			wakeupOn(w);
		}

	    public void processStimulus(Enumeration criteria) {

			if(!run) {
				wakeupOn(w);
				return;
			}

			int nlanded = 0;
			int nstarted = 0;
			int nlander = 0;

			Time += dtime;

			// check if we need to change the wind

			boolean windChanged = false;

			// check first for special wind commands

			if (windTimes[0] < 0)
			{
				if(Time - lastWindChangeTime >= Math.abs(windTimes[0])) {
					lastWindChangeTime = Time;
					if(windSpeeds[0] >= 0) {
						wind = windSpeeds[0];
					} else {
						wind = numWindChanges*Math.abs(windSpeeds[0]);
					}
					double angle = windAngles[0];
					if(windAngles[0] < 0) {
						angle = numWindChanges*Math.abs(windAngles[0]) % 360;
					}
					numWindChanges++;

					windVector = new Vector3f((float)Math.cos(angle / 360. * 2.0 * Math.PI), (float)Math.sin(angle / 360. * 2.0 * Math.PI), 0.0f);
					LMessage.setText("Time= "+Time+" Wind set to " + wind + " metres/sec, angle " + angle + " degrees");
					System.out.println("Wind set to " + wind + " metres/sec, angle " + angle + " degrees");

					windSpeed.setValue((int)(wind * 100)); // Set the slider to the wind speed
					windDirection.setValue((int)(angle)); // Set the slider to the wind angle

					windChanged = true;

				}
			}
			else if (currentWind < nWinds - 1)
			{
				if (Time > windTimes[currentWind + 1])
				{
					currentWind++;

					wind = windSpeeds[currentWind];
					double angle = windAngles[currentWind];
					windVector = new Vector3f((float)Math.cos(angle / 360. * 2.0 * Math.PI), (float)Math.sin(angle / 360. * 2.0 * Math.PI), 0.0f);
					LMessage.setText("Wind set to " + wind + " metres/sec, angle " + angle + " degrees");
					System.out.println("Wind set to " + wind + " metres/sec, angle " + angle + " degrees");

					windSpeed.setValue((int)(wind * 100)); // Set the slider to the wind speed
					windDirection.setValue((int)(angle)); // Set the slider to the wind angle

					windChanged = true;
				}
			}

			debug = CBoxDebug.isSelected();
		
			for(int i=0;i<blobs.size();i++) {

				Blob b = (Blob) blobs.elementAt(i);

				if(b.landed) {
					nlanded++;
				} else {
					// if the wind just changed then recalculate edges
					if (windChanged) b.calculateEdges();
					b.update();
					if (b.started)
					{
						nstarted++;
					}
					else
					{
						nlander++;
					}
				}

				if (i % 1000 == 0) System.out.print("!");

			}

			if (CBoxMovie.isSelected())
			{
				takeSnapShot();
			}

			boolean bFinish = false;
			// Simulation finishes if no particles in motion, and no more wind changes
			if(nstarted == 0 && currentWind == nWinds-1 && bBatch) bFinish = true;
			
			if(bFinish) {			
				LMessage.setText("T = "+(float)Time+" s All Particles Stationary");
				System.out.println("All particles stationary");
				run = false;
				if (bBatch)
				{
					System.out.println("Simulation complete");
					writeResults();
					if (iExperiment < nExperiments)
					{
						iExperiment++;
						System.out.println("Starting experiment " + iExperiment);
						initializeSim();
						run = true;
					}
					else
					{
						System.out.println("Processing Complete");
						double elapsed = (System.currentTimeMillis() - timestart) * 0.001;
						System.out.println("Time elapsed for simulation = " + elapsed + "secs, i.e. " + elapsed / (double)nExperiments + " secs/experiment");
						System.exit(0);
					}
				}
				// all particles are now stationary
			} else {
				LMessage.setText("T = "+(float)Time+" s "+nlander+" on Lander "+(nstarted)+" in Motion, "+nlanded+" on Surface");
			}
			
			wakeupOn(w);
		}
	    public FrameBehavior(Alpha A, TransformGroup T) {
		}
	}

	public class SliderListener implements ChangeListener {
		public void stateChanged(ChangeEvent e) {
			JSlider source = (JSlider)e.getSource();
			if (!source.getValueIsAdjusting()) {
				if(source.equals(windSpeed)) {
					wind = 0.01*(double)source.getValue(); // slider is in cm/sec
					LMessage.setText("Wind speed set to "+wind+" metres/sec");
					System.out.println("Wind speed set to " + wind + " metres/sec");
				} else if(source.equals(windDirection)) {
					double angle = (double)source.getValue();
					windVector = new Vector3f((float)Math.cos(angle/360.*2.0*Math.PI),(float)Math.sin(angle/360.*2.0*Math.PI),0.0f);
					LMessage.setText("Wind direction set to "+angle+" degrees");
					System.out.println("Wind direction set to " + angle + " degrees");

					Transform3D t1 = new Transform3D();
					tgWind.getTransform(t1);
					double scale = t1.getScale();
					t1.rotZ(angle/360.*2.0*Math.PI - Math.PI*0.5);
					//t1.setScale(scale);
					Vector3d translation = new Vector3d();
					//translation.x = -800;
					//translation.y = 0;
					//translation.z = 0;
					t1.setTranslation(windTranslation);
					tgWind.setTransform(t1);

				}

				// recalculate all leading/trailing edges for particles
				System.out.print("Calculate new edges ...");
				for (int i = 0; i < blobs.size(); i++)
				{
					Blob b = (Blob)blobs.elementAt(i);
					b.calculateEdges();
				}
				System.out.println(" Done");


			} else if(source.equals(windDirection)) {

					double angle = (double)source.getValue();
					windVector = new Vector3f((float)Math.cos(angle/360.*2.0*Math.PI),(float)Math.sin(angle/360.*2.0*Math.PI),0.0f);
					LMessage.setText("Wind direction set to "+angle+" degrees");

					Transform3D t1 = new Transform3D();
					tgWind.getTransform(t1);
					double scale = t1.getScale();
					t1.rotZ(angle/360.*2.0*Math.PI - Math.PI*0.5);
					//t1.setScale(scale);
					Vector3d translation = new Vector3d();
					translation.x = -800;
					translation.y = 0;
					translation.z = 0;
					t1.setTranslation(windTranslation);
					tgWind.setTransform(t1);
			}
		}
	}


	public class eventHandler implements ActionListener{
										 

		public void actionPerformed(ActionEvent event) {

			String command = event.getActionCommand();
			System.out.println(command);

			if(command.equals("Hide Lander")) {
				landerBG.detach();
				BAttach.setText("Show Lander");
				return;
			}
			if(command.equals("Show Lander")) {
				objTrans.addChild(landerBG);
				BAttach.setText("Hide Lander");
				return;
			}

			if(command.equals("New")) {
				run = false;
				BStart.setText("Start");

				landerBG.detach();
				objTrans.addChild(landerBG);
				BAttach.setText("Hide Lander");

				cBeh.setEnable(false);

				initializeSim();				

				cBeh.setEnable(true);

				return;
			}
			if(command.equals("Start")) {
				run = true;
				BStart.setText("Stop");
				return;
			}
			if(command.equals("Stop")) {
				run = false;
				BStart.setText("Start");
				return;
			}
			if(command.equals("Save Current View")) {
				takeSnapShot();
				return;
			}

			if (command.equals("Write Data File"))
			{
				writeResults();
				return;
			}

			if(command.equals("Top View")) {
				Vector3d translation = new Vector3d();
				Transform3D t1 = new Transform3D();
				translation.z = 0;
				t1.setTranslation(translation);
				t1.setScale(DEFAULTSCALE);
				objTrans.setTransform(t1);
				return;
			}
			if(command.equals("Side View")) {
				Vector3d translation = new Vector3d();
				Transform3D t1 = new Transform3D();
				t1.rotX(-0.5 * Math.PI);
				t1.setScale(DEFAULTSCALE);
				translation.z = 0;	
				t1.setTranslation(translation);
				objTrans.setTransform(t1);
				//		Matrix3d twist = new Matrix3d();
				//twist.rotZ(-Math.PI * 0.5);

				//twind.setRotation(twist);

				return;
			}
			if(command.equals("Zoom In")) {
				Transform3D t1 = new Transform3D();
				objTrans.getTransform(t1);
				double scale = t1.getScale();
				LMessage.setText("Old zoom scale "+scale);
				scale *= 1.2;
				t1.setScale(scale);
				objTrans.setTransform(t1);
				return;
			}
			if(command.equals("Zoom Out")) {
				Transform3D t1 = new Transform3D();
				objTrans.getTransform(t1);
				double scale = t1.getScale();
				LMessage.setText("Old zoom scale "+scale);
				scale /= 1.2;
				t1.setScale(scale);
				objTrans.setTransform(t1);
				return;
			}
			if (command.equals("Right"))
			{
				Vector3d translation = new Vector3d();
				Transform3D t1 = new Transform3D();
				objTrans.getTransform(t1);
				t1.get(translation);
				translation.x += 0.1;
				t1.setTranslation(translation);
				objTrans.setTransform(t1);
				return;
			}
			if (command.equals("Left"))
			{
				Vector3d translation = new Vector3d();
				Transform3D t1 = new Transform3D();
				objTrans.getTransform(t1);
				t1.get(translation);
				translation.x -= 0.1;
				t1.setTranslation(translation);
				objTrans.setTransform(t1);
				return;
			}




		}

	}


	void writeResults()
	{
		int nblobs = blobs.size();
		if (iExperiment == 1)
		{
			try
			{
				dataStream = new PrintStream(new FileOutputStream(new File(outputFile)));
			}
			catch (Exception e)
			{
				System.err.println("Whoops: cannot open " + outputFile + " Error:" + e);
				return;
			}

			Date date = new Date();
			dataStream.println("# Date:   " + date);
			dataStream.println("# Config: " + configFile);
			dataStream.println("# VRML:   " + vrmlFile);
			dataStream.println("# "+windTriplets);
			dataStream.println("# Codes are: 1=OnLander, 2=InMotion, 3=LandedOnSurface");
			dataStream.println("# Code\tX\tY\tZ\tRadius\tXfree\tYfree\tZfree");
			dataStream.println(nblobs*nExperiments);
		}
		System.out.println("Results are being written to " + outputFile);
		for (int i = 0; i < nblobs; i++)
		{

			Blob b = (Blob)blobs.elementAt(i);

			int code = 1;
			if (b.started) code = 2;
			if (b.landed) code = 3;

			dataStream.println(code + "\t" + (b.Xpos * GEOMETRYSCALEFACTOR) + 
				"\t" + (b.Ypos * GEOMETRYSCALEFACTOR) + 
				"\t" + (b.Zpos * GEOMETRYSCALEFACTOR) + 
				"\t" + b.radius + 
				"\t" + (b.Xfree * GEOMETRYSCALEFACTOR) + 
				"\t" + (b.Yfree * GEOMETRYSCALEFACTOR) + 
				"\t" + (b.Zfree * GEOMETRYSCALEFACTOR));

		}
		if (iExperiment == nExperiments)
		{
			dataStream.flush();
			dataStream.close();
			System.out.println("Results written");
		}
	}


}





class MouseRotate extends MouseBehavior {
    double x_angle, y_angle;
    double x_factor = .03;
    double y_factor = .03;

  private MouseBehaviorCallback callback = null;


  /**
   * Creates a rotate behavior given the transform group.
   * @param transformGroup The transformGroup to operate on.
   */
  public MouseRotate(TransformGroup transformGroup) {
    super(transformGroup);
  }

  /**
   * Creates a default mouse rotate behavior.
   **/
  public MouseRotate() {
      super(0);
  }

  /**
   * Creates a rotate behavior.
   * Note that this behavior still needs a transform
   * group to work on (use setTransformGroup(tg)) and
   * the transform group must add this behavior.
   * @param flags interesting flags (wakeup conditions).
   */
  public MouseRotate(int flags) {
      super(flags);
   }

  public void initialize() {
    super.initialize();
    x_angle = 0;
    y_angle = 0;
    if ((flags & INVERT_INPUT) == INVERT_INPUT) {
       invert = true;
       x_factor *= -1;
       y_factor *= -1;
    }
  }

  /**
   * Return the x-axis movement multipler.
   **/

  public double getXFactor() {
    return x_factor;
  }
  
  /**
   * Return the y-axis movement multipler.
   **/

  public double getYFactor() {
    return y_factor;
  }
  
  /**
   * Set the x-axis amd y-axis movement multipler with factor.
   **/

  public void setFactor( double factor) {
    x_factor = y_factor = factor;
    
  }
  
  /**
   * Set the x-axis amd y-axis movement multipler with xFactor and yFactor
   * respectively.
   **/

  public void setFactor( double xFactor, double yFactor) {
    x_factor = xFactor;
    y_factor = yFactor;    
  }

  public void processStimulus (Enumeration criteria) {
      WakeupCriterion wakeup;
      AWTEvent[] event;
      int id;
      int dx, dy;

      while (criteria.hasMoreElements()) {
         wakeup = (WakeupCriterion) criteria.nextElement();
         if (wakeup instanceof WakeupOnAWTEvent) {
            event = ((WakeupOnAWTEvent)wakeup).getAWTEvent();
            for (int i=0; i<event.length; i++) { 
	      processMouseEvent((MouseEvent) event[i]);

	      if (((buttonPress)&&((flags & MANUAL_WAKEUP) == 0)) ||
		  ((wakeUp)&&((flags & MANUAL_WAKEUP) != 0))){
		
		id = event[i].getID();
		if ((id == MouseEvent.MOUSE_DRAGGED) && 
		    !((MouseEvent)event[i]).isMetaDown() && 
		    !((MouseEvent)event[i]).isAltDown()){
		  
                  x = ((MouseEvent)event[i]).getX();
                  y = ((MouseEvent)event[i]).getY();

                  dx = x - x_last;
                  dy = y - y_last;

		  if (!reset){	    
		    x_angle = dy * y_factor;
		    y_angle = dx * x_factor;
		    
		    transformX.rotX(x_angle);
		    transformY.rotY(y_angle);

			//System.out.println(x_angle + " " + y_angle);
	    
		    transformGroup.getTransform(currXform);
		    
		    //Vector3d translation = new Vector3d();
		    //Matrix3f rotation = new Matrix3f();
		    Matrix4d mat = new Matrix4d();
		    
		    // Remember old matrix
		    currXform.get(mat);
		    
		    // Translate to origin
		    currXform.setTranslation(new Vector3d(0.0,0.0,0.0));
		    if (invert) {
			currXform.mul(currXform, transformX);
			currXform.mul(currXform, transformY);
		    } else {
			currXform.mul(transformX, currXform);
			currXform.mul(transformY, currXform);
		    }
		    
		    // Set old translation back
		    Vector3d translation = new 
		      Vector3d(mat.m03, mat.m13, mat.m23);
		    currXform.setTranslation(translation);
		    
		    // Update xform
		    transformGroup.setTransform(currXform);

		    transformChanged( currXform );

                    if (callback!=null)
                        callback.transformChanged( MouseBehaviorCallback.TRANSLATE,
                                               currXform );


		  }
		  else {
		    reset = false;
		  }

                  x_last = x;
                  y_last = y;
               }
               else if (id == MouseEvent.MOUSE_PRESSED) {
                  x_last = ((MouseEvent)event[i]).getX();
                  y_last = ((MouseEvent)event[i]).getY();
               }
	      }
	    }
         }
      }

      wakeupOn (mouseCriterion);
      
   }

  /**
    * Users can overload this method  which is called every time
    * the Behavior updates the transform
    *
    * Default implementation does nothing
    */
  public void transformChanged( Transform3D transform ) {
  }

  /**
    * The transformChanged method in the callback class will
    * be called every time the transform is updated
    */
  public void setupCallback( MouseBehaviorCallback callback ) {
      this.callback = callback;
  }
}





           
