// Bouncing Ball Gravity Simulator
// Takes into account: gravity, energy loss due to collisions, 
// air friction for laminar flow, and drag for turbulent flow.
// 3/24/2008
// You may reproduce this code as long as you give me credit for the original:
// Burak Kanber
// Technical Director, Niqos Inc
// The Cooper Union

// Please remember that this is essentially a scale model. Some of the values have been
// adjusted contrary to real world values for your viewing pleasure. If you want this 
// to be really realistic, you'll have to adjust the scale_factor below to something larger,
// say .0005, not .0001. 
// Remember, that little red ball on your screen is 10 meters (33 feet!) in diameter.
// The values have been adjusted so that the user can play with the ball and have it bounce
// around and such. With more realistic numbers relating to scale, the ball won't fly very 
// far horizontally before drag takes a hold of it, and that's not much fun now is it?

// Global parameters
// Acceleration due to gravity on earth: 9.81 m/s^2
// Change it to 1.6 to simulate moon gravity.
// Change it to 26 to simulate Jupiter gravity.
// Change it to 0.0 to simulate no gravity at all. 
float ay = 9.8;

// Velocities in x and y direction (set value explicitly for given initial velocity)
// Velocity in m/s
float vx = random(-25,25);
float vy = random(-25,25);

// Geometry of screen, in meters
int width = 500;
int height = 500;

// Ball diameter, in meters
int diameter = 10;

// Ball mass, in kg
float mass = 3;


// Environmental variables
//
// rho, the density of air (or whatever medium you're in)
// ambient atmospheric earth (SATP), rho=1.168 kg/m^3
// for STP, rho=1.292 kg/m^3
float rho=1.168;
// characteristic area-- for a sphere is the cross-sectional area in m^2
// you can change this if you want to experiment with different values, 
// but it's probably just easier to change the drag coeffecient
float characteristic_area = (PI/4)*diameter*diameter;
// coeffecient of drag. For a rough sphere, Cd = .4. Dimensionless
float drag_coeff = .4;
// Dynamic viscosity of air. For SATP, 3.54(10^-5) Ns/m^2
float viscosity = 3.54*pow(10,-5);
// Reynolds number. Dimensionless, used to determine whether flow is 
// laminar or turbulent. Re = rho*velocity*diameter/viscosity.
float reynolds = 0;
// Cutoff of turbulent/laminar flow. In general, for boundary layers, 
// Re > 10^6 displays turbulent flow.
float turbulence = 1000000;
// Drag acceleration, just defined here to be used later.
// Ad = .5*rho*v^2*Cd*A/m
float drag_accel_x=0;
float drag_accel_y=0;
// Drag force scale factor. dimensionless constant used to adjust for scale.
// I _THINK_ if you want to simulate real life, change this to .07.
// I use .0004 because its more fun.
// If you set the diam to 10, the mass to 3, and the scale_coeff to .0006, 
// the simulation roughly emulates a 1 (not 10) meter ball of 3kg.
float scale_coeff=.0004;
// Air Friction Coefficient, dimensionless. Used for laminar flow case.
float fric_coeff=.1;
// Energy kept after a collision. Dimensionless. 
// Also known as coeffecient of restitution. Must be negative, or 
// the ball won't change direction. If made larger than 1, the ball will 
// bounce back higher than it dropped....
float v_decay=-.727;
float h_decay=-.727;
// Incremental time in seconds. Smaller values give you slower animation, 
// but more accurate (though not necessarily to the human eye) motion.
// 0.14 "feels" about right. Take into consideration the dimensions and the
// frame rate and speed of your computer when setting this. I haven't done
// the necessary calculations to figure out the best value for this.
float t = 0.14;
// X and Y Position
float y = 100.0;
float x = 200.0;
// Do Not Edit Below This Line
PFont font;
String lam_turb="";
int initMousePositionX=0;
int initMousePositionY=0;
boolean dragging=false;

void setup() {
  font=loadFont("Arial-Black-48.vlw");
  textFont(font,16);
  smooth();
  size(width,height+50);
}
void drawbg() {
  // Called by draw(), takes care of the ball and the BG
  fill(255);
  rect(0,0,width, height);
  fill(0);
  rect(0,height,width,50);
  fill(255);
  textFont(font,16);
  text("Velocity (m/s): " + nf(sqrt(vx*vx + vy*vy),0,1) ,20,height+20);
  text(lam_turb,20,height+40);
  text("Burak Kanber",width-150,height+20);
  text("www.niqos.com",width-150,height+40);
  fill(255,0,0);
  ellipse(x,y,ceil(diameter),ceil(diameter));
  if (dragging) {
    // This handles drawing the aiming line when you shoot the ball
    strokeWeight(2);
    line(initMousePositionX,initMousePositionY,mouseX,mouseY);
    strokeWeight(1);
  }
}

void draw() {
  // Draw the background and the ball.
  drawbg();
  if (!dragging) {
    // Kinematic equation for vertical motion, from
    // Xf = Xi + Vi*t + .5*a*t^2
    y = y + (vy*t) + (.5*ay*(t*t));
    // Same equation, for horizontal motion, without gravity term
    x = x + (vx*t); 
    // We need to define the air friction losses explicitly here.
    float fric_y=0;
    float fric_x=0;
    // Determine reynolds number 
    reynolds = rho*sqrt(vx*vx + vy*vy)*diameter / viscosity;
    // Determine if flow is laminar or turbulent.
    if (reynolds > turbulence) {
      lam_turb="Turbulent";
      // Calculate drag force
      drag_accel_x = scale_coeff*vx*vx*rho*drag_coeff*characteristic_area/mass;

      if (vx != 0) {
        vx += -(vx/abs(vx))*drag_accel_x*t;
      }
      drag_accel_y = scale_coeff*vy*vy*rho*drag_coeff*characteristic_area/mass;
      if (y < height-diameter/2) {
        if (vy == 0) {
          vy = vy + (ay*t);
        } 
        else {
          vy = vy + (ay*t) - (vy/abs(vy))*drag_accel_y*t;
        }
      }
    } 
    else {
      lam_turb="Laminar";
      // To get the "direction" of the ball, we do vy/abs(vy)
      // This contains a possible divide-by-zero, (which is why we
      // separated this stuff in the first place), so we check for that here.
      // This isn't proper from a Fluid Mechanics standpoint, as real air 
      // friction isn't linear. I'll make this real air friction eventually.
      if (!(vy==0)) { 
        fric_y= (-vy/abs(vy))*fric_coeff; 
      }
      if (!(vx==0)) { 
        fric_x= (-vx/abs(vx))*fric_coeff; 
      }
      // Y- velocity gets affected by both gravity and friction.
      // First part from Vf = Vi + a*t (first derivative of kinematic eqn above).
      if (y < height-diameter/2) {
        vy = vy + (ay*t) + fric_y*t/mass;
      }
      // X Velocity only bothered by friction
      vx = vx + fric_x*t/mass;
    }


    // Test to see if ball has flown past boundaries.
    // If so, force it back in, otherwise it'll get stuck because of the
    // energy decay.
    // The decay terms are negative, so they reverse the ball's velocity.
    // They're smaller than 1 so they lessen the ball's velocity as well.
    // In real life, it may not be a linear relationship like this, but
    // it's a good approximation.
    if (y > height - (diameter/2)) { 
      y=height - (diameter/2); 
      vy = vy*v_decay;
    }

    if (x > width - (diameter/2)) { 
      x=width - (diameter/2); 
      vx=vx*h_decay;
    }
    if (x < diameter/2) { 
      x=diameter/2;
      vx=vx*h_decay; 
    }
  }
}
void mousePressed() {
  // So the drawbg() and draw() functions know whet to do
  dragging = true;
  // Remember the initial mouse position.
  initMousePositionX=mouseX;
  initMousePositionY=mouseY;
  // Set the ball here
  x=mouseX;
  y=mouseY;
  // Get rid of the velocities. This step isn't necessary, but you want to
  // control everything you can.
  vx=0;
  vy=0;
}

void mouseReleased() {
  // Start motion again:
  dragging=false;
  // Subtract the initial mouse position from the current one to give us 
  // the velocity in each direction. The bigger the difference, 
  // the faster she flies.
  vx = initMousePositionX - mouseX;
  vy = initMousePositionY - mouseY;
}