import java.lang.Math.*;

class Ball {
	final double pi = Math.PI;
//	physical constants_______________________________________
	final double grav = 9.81;	//  acceleration due to gravity in m/s²
	final double pAtm = 101000;	//  atmospheric pressure in Pa
	final double gamma = 1.4;	//  Cp/Cv ratio for adiabatic compression
//	parameters describing the ball; default values supplied by applet
	double mass;		//  mass of ball in kg 
	double radius;		//  radius of undistorted ball in m
	double kappa;		//  modulus of elasticity of ball in N/m
	double decay;		//  decay time for oscillation in s
	double press0;		//  gauge pressure in undistorted ball, in Pa
//	the next two parameters set the accuracy of simulation___
	int N;				//  no. of slices = size of arrays
	double dt;			//  time increment
//	derived constants__________________________________________
	double vol0;		//  volume of undistorted ball (m^3)
	double strs0;		//  surface stress due to initial pressure (N/m)
	double damp;		//  damping factor for surface waves
	double dens;		//  surface density of casing of ball (kg/m²)
	double ds0;			//  width of undistorted ring elements (m)
//	Arrays_________________________________________________
	double r0[];		//undistorted radius of ring N
	double r1[];		//mean of undistorted radii of rings N, N+1
	double r[];	 		//actual radius of ring N
	double h[];	 		//height of ring N above surface
	double vel_r[];  	//radial velocity of element of ring N
	double vel_v[]; 	//vertical velocity of element
	double kvf[];		// inertia factor	
//	Variables_________NB up is +ve__________________________
	double press;		// actual pressure in ball
	double drop;		// drop height
	double h_av;		// average height (centre of mass)
	double vel_av;  	// average vertical velocity of ball
	double fr = 0; 		// radius of "footprint" 
	double frmax = 0;	// max. radius of footprint
	double contact = 0;	// contact time
	double hBounce = 0;	// height of bounce, calculated at liftoff
	double effic = 0;	// height of bounce / drop
	int stage = 0;		// flag for switch
	
//	Constructor*****************************************
	Ball (BouncyApplet o) {
		N = o.nRings;
		press0 = o.press * 1e3;		// kPa -> Pa
		drop = o.drop;				// m
		mass = o.mass * 1e-3;		// g -> kg
		radius = o.radius * 1e-2;	// cm -> m
		decay = o.decay * 1e-3;		// ms -> s
		kappa = o.kappa * 1e3;		// kN/m -> N/m
		dt = o.dt * 1e-6;			// µs -> s
		strs0 = radius*press0/2;
		dens = mass/(4*pi*radius*radius);
		ds0 = 2*radius*Math.sin(pi/N/2.);
		damp = (1-Math.exp(-N*N*dt/decay/4/pi/pi));
		h_av = 0.0443*Math.sqrt(drop);  // 10ms before impact
		vel_av = -Math.sqrt(2*grav*(drop - h_av));

		r0 = new double[N];
		r1 = new double[N];
		r = new double[N+1];
		h = new double[N+1];
		vel_r = new double[N];
		vel_v = new double[N];
		kvf = new double[N];
		
		for (int m=0; m<N; m++) {
			r0[m]= r[m] = radius * Math.sin(pi*(m+0.5)/N);
			r1[m]= radius * Math.sin(pi*(m+1.)/N);
			h[m] = h_av + radius*(1 - Math.cos(pi*(m+0.5)/N));
			vel_r[m] = 0;
			vel_v[m] = vel_av;
			kvf[m] = 1/(dens*r0[m]*ds0);
		}
		r[N] = -r[N-1];		// this allows neat calculations which 
		h[N] = h[N-1];		// include the top ring with the rest.
		
		vol0 = calcVol();	// used to calculate pressure changes
	}

//	methods ***********************************************
	double calcVol() {
		double vol = r[0]*r[0]*pi*(h[1]-h[0])/2.;
		for (int m=1; m<N; m++) {
			vol += r[m]*r[m]*pi*(h[m+1]-h[m-1])/2.;
		}
		return vol;
	}

/*	This method checks whether the ball is touching the floor. 
 *  Rings that have fallen "through" the floor are put back on 
 *  the floor and their velocities, if downward, set to zero.
 *  The radius of the "footprint" and the contact time are also
 *  found; the rebound height is calculated at lift-off.
 */	
	int onFloor() {
		int m1 = -1;			// this value indicates no contact
		for (int m=0; m<N/2; m++) {
			if (h[m]<= 0.) {
				m1 = m;
				h[m] = 0.;
				if (vel_v[m] < 0) vel_v[m] = 0.;
			}
		}
		if (m1 >= 0) {
			fr = r[m1];			// radius of "footprint"
	//	  if (h[m1+2] > 2.*h[m1+1])
	//		  fr += (r[m1+1]-r[m1])/(h[m1+2]/h[m1+1]-1);
		}
		else fr = 0;
		
		if (fr > 0) {
			contact += dt;
			if (fr > frmax)frmax = fr;
		}
		switch (stage){
			case 0:				// before impact
				if (fr > 0) stage = 1; 
				break;
			case 1:				// in contact with the floor
				if (fr == 0) {  // lifting off
					stage = 2; 
					hBounce = h_av + (vel_av*vel_av)/2/9.81;
					effic = hBounce/drop;
					hBounce = 0.0443*Math.sqrt(hBounce);
				}
				break;
			case 2:				// after lift-off
				if (h_av > hBounce)
					stage = 3;	//stop iterating
		}
		return stage;
	}
	
	void calcChange() {
//		calculate stresses, forces and motions in each ring element
		double dr, dr1;
		double dh, dh1, ds, ds1 = ds0*Math.cos(pi/2/N);
		double strsP, strsE;
		double forceN, forceNr, force_r, forceNv, force_v;
		double forceSr, forceSv;
		int m;
	
		press = (press0+pAtm)*Math.pow((vol0/calcVol()),gamma) - pAtm;
		vel_av =  h_av = 0;
		
		forceSr =  forceSv = 0;
		dr1 = r[0]*2;
		dh1 = 0;
 		for (m=0; m<N; m++) {
//			calculate spacings and stresses
			dr = r[m+1]-r[m];
			dh = h[m+1]-h[m];
			ds = Math.sqrt(dr*dr+ dh*dh);
			if (ds >2*ds0) stage = 3;
			strsP = kappa*(ds/ds0 -1.) + strs0;	// polar stress (N-S)
			strsE = kappa*(r[m]/r0[m]-1.) + strs0; // equatorial stress
			forceN = strsP*r1[m];	// force up on element m, due to m+1
//			calculate radial components of force and motion
			forceNr = forceN*dr/ds;
			force_r = - strsE*ds1	   //inward force due to eq. stress
					  + press*r[m]*(dh+dh1)/2. //outward force due to pressure
					  + forceNr
					  + forceSr;		//force on element m due to m-1
			forceSr = - forceNr;		//store to use next time around
			vel_r[m] += force_r*kvf[m]*dt;
			r[m] += vel_r[m]*dt;
//			calculate vertical components of force and motion
			forceNv = forceN*dh/ds;
			force_v = - press*r[m]*(dr+dr1)/2.  //downward if r[m+1] > r[m-1]
					  + forceNv + forceSv;
			forceSv = -forceNv;			 // store, again
			vel_v[m] += (force_v*kvf[m]- grav)*dt;
			vel_av += vel_v[m]*r0[m];
			h[m] += vel_v[m]*dt;
			h_av += h[m]*r0[m];
			dh1 = dh; dr1 = dr;
		}
		r[N] = -r[N-1];		// update for future use
		h[N] = h[N-1];
	
		vel_av *= dens*ds0*2.*pi/mass;
		h_av = h_av*dens*ds0*2.*pi/mass - radius;
	}
		
	void damping(){
		double forceNr, forceNv, forceSr, forceSv;
		forceSr = forceSv = 0;
		for (int m=0; m<N-1; m++) {
			forceNr = (vel_r[m+1] -vel_r[m])*r1[m];
			forceNv = (vel_v[m+1] -vel_v[m])*r1[m];
			vel_r[m] += damp*(forceNr + forceSr)/r0[m];
			vel_v[m] += damp*(forceNv + forceSv)/r0[m];
			forceSr = -forceNr; forceSv = -forceNv;
		}
		vel_r[N-1] += damp*forceSr/r0[N-1];
		vel_v[N-1] += damp*forceSv/r0[N-1];
	}
}