import java.*; // converted from original Fortran-90 code...yikes! public class BB { double sigma, delta, theta; public static int npts = 0; double PI=3.14159265; double FT2METERS = 0.3048; // convert feet to meters; double G = 9.8; // acceleration of gravity double LBS2KG = 0.45359237; // weight in pounds to mass in kg. double RHOZERO = 1.2250; // density of air at sealevel, kg/cu.m /* BASEBALL CONSTANTS - by the rules, the circumference of a baseball must be no less than 9 inches and no more than 9.25 inches. The weight of a baseball must be no less than 5.00 ounces and no more than 5.25 ounces. I assume my ideal ball lies exactly in the middle. In SI units: */ double DIAM=(9.125/(12*PI))*FT2METERS;// diameter of a baseball (m) double MASS=(5.125/16)*LBS2KG; // mass of a baseball (kg) double SREF = 0.25*PI*DIAM*DIAM; // frontal area (sq.m) /** * @param rz new Rho-Zero value for trajectory */ public void setRhoZero(double rz) { RHOZERO = rz; } /* FUNCTION Accel(time, position, velocity) RESULT(acceleration) * not PURE * --------------------------------------------------------------------------- * PURPOSE - Compute the acceleration (vector) for a spherical projectile * moving through a viscous medium. Assume Mach number is small enough * that wave drag may be neglected. Ignore added mass term. * NOTE - position has units of meters, but first argument to SimpleAtmosphere * is in kilometers. Be sure to remember to multiply by 0.001 */ public double[] Accel(double time, double[] position, double[] velocity) { double[] acceleration = new double[2]; double[] VERTICAL = {0., 1.}; double cd, density, dragMagnitude, q, reynolds, vmag, vsq; double[] drag = new double[2]; double[] unitVelocity = new double[2]; vsq = Math.pow(velocity[0],2) + Math.pow(velocity[1],2) ; vmag = Math.sqrt(vsq); unitVelocity[0]=velocity[0]/vmag; unitVelocity[1]=velocity[1]/vmag; SimpleAtmosphere(0.001*position[1]) ; //first arg is altitude in kilometers density=sigma*RHOZERO; q=0.5*density*vsq; reynolds=density*vmag*DIAM/MetricViscosity(theta); cd=CdSphere(reynolds); dragMagnitude=cd*q*SREF; drag[0]=-dragMagnitude*unitVelocity[0]; drag[1]=-dragMagnitude*unitVelocity[1]; acceleration[0]=drag[0]/MASS - G*VERTICAL[0]; acceleration[1]=drag[1]/MASS - G*VERTICAL[1]; return acceleration; } /* FUNCTION CdSphere(r) RESULT(d) * --------------------------------------------------------------------------- * PURPOSE - Compute the drag coefficient of a sphere as a function of * Reynolds number. Assumes Mach number is small. * Taken from Chow, "Computational Aerodynamics" REAL,INTENT(IN):: r ! Reynolds number REAL:: d ! drag coefficient based on cross-section area !---------------------------------------------------------------------------- */ public double CdSphere(double r) { double d; if (r <= 0) { d=0; } else if (r <= 1.0) { d=24.0/r; } else if (r <= 400.0) { d=24*Math.pow(r,-0.646); } else if (r <= 3.e5) { d=0.5; } else if (r <= 2.e6) { d=3.66E-4*Math.pow(r,0.4275); } else { d=0.18; } return d; } /* SUBROUTINE CorrectFinalPosition(initialAltitude, a1, a2) ! --------------------------------------------------------------------------- ! PURPOSE - Correct the final point of a trajectory so that its y-coordinate ! is exactly equal to initial altitude. Assume that the altitudes of a1 ! and a2 straddle initialAltitude. */ public TrajectoryPoint CorrectFinalPosition(double initialAltitude, TrajectoryPoint a1, TrajectoryPoint a2) { double fraction; fraction=(initialAltitude-a1.position[1])/(a2.position[1]-a1.position[1]); TrajectoryPoint a3 = new TrajectoryPoint(); a3.time =a1.time + fraction*(a2.time - a1.time); a3.position[0]=a1.position[0] + fraction*(a2.position[0]-a1.position[0]); a3.position[1]=a1.position[1] + fraction*(a2.position[1]-a1.position[1]); a3.velocity[0]=a1.velocity[0] + fraction*(a2.velocity[0]-a1.velocity[0]); a3.velocity[1]=a1.velocity[1] + fraction*(a2.velocity[1]-a1.velocity[1]); a3.acceleration[0]=a1.acceleration[0] + fraction*(a2.acceleration[0]-a1.acceleration[0]); a3.acceleration[1]=a1.acceleration[1] + fraction*(a2.acceleration[1]-a1.acceleration[1]); return a3; } /* !+ FUNCTION MetricViscosity(theta) RESULT(visc) ! ------------------------------------------------------------------------- ! PURPOSE - Compute viscosity using Sutherland's formula. ! Returns viscosity in kg/(meter-sec) */ public double MetricViscosity(double theta) { double BETAVISC = 1.458E-6; // viscosity term, N sec/(sq.m sqrt(deg K) double SUTH = 110.4; // Sutherland's constant, deg K double TZERO = 288.15; // temperature at sealevel, deg K double temp = TZERO * theta; double visc = BETAVISC * Math.sqrt(temp*temp*temp) / (temp+SUTH); return visc; } /* !+ SUBROUTINE SimpleAtmosphere(alt,sigma,delta,theta) ! ------------------------------------------------------------------------- ! PURPOSE - Compute the characteristics of the lower atmosphere. ! NOTES-Correct to 20 km. Only approximate above there */ public void SimpleAtmosphere(double alt) { double h; // h=geopotential altitude double REARTH = 6369.0; // radius of the Earth (km) double GMR = 34.163195; // gas constant h = alt * REARTH / (alt+REARTH); //convert geometric to geopotential alt if (h < 11.0) { theta=1.0+(-6.5/288.15)*h; // Troposphere delta=Math.pow(theta,(GMR/6.5)); } else { theta=216.65/288.15; // Stratosphere delta=0.2233611*Math.exp(-GMR*(h-11.0)/216.65); } sigma=delta/theta; return; } /* !+ SUBROUTINE BaseBallKutta(p1, h, p2) ! --------------------------------------------------------------------------- ! PURPOSE - Advance one time-like step in a trajectory. This is a system ! of four first order ordinary differential equations. Use fourth-order ! Runge-Kutta equation to advance one time step. */ public TrajectoryPoint BaseBallKutta(TrajectoryPoint p1, double h) { double HALF=.5; double SIXTH = 1./6.; double [] dx1, dx2, dx3, dx4; double [] dv1, dv2, dv3, dv4; dx1 = new double[2]; dx2 = new double[2]; dx3 = new double[2]; dx4 = new double[2]; dv1 = new double[2]; dv2 = new double[2]; dv3 = new double[2]; dv4 = new double[2]; double[] position = new double[2]; double[] velocity = new double[2]; double t=p1.time; double[] x=p1.position; double[] v=p1.velocity; double[] a=Accel(t, x, v); dx1[0] = h*v[0]; dx1[1] = h*v[1]; dv1[0] = h*a[0]; dv1[1] = h*a[1]; position[0] = x[0] + HALF*dx1[0]; position[1] = x[1] + HALF*dx1[1]; velocity[0] = v[0] + HALF*dv1[0]; velocity[1] = v[1] + HALF*dv1[1]; a = Accel(t+HALF*h, position, velocity); //a=Accel(t+HALF*h, x+HALF*dx1, v+HALF*dv1); dx2[0] = h*(v[0] + HALF*dv1[0]); dx2[1] = h*(v[1] + HALF*dv1[1]); dv2[0] = h*a[0]; dv2[1] = h*a[1]; position[0] = x[0] + HALF*dx2[0]; position[1] = x[1] + HALF*dx2[1]; velocity[0] = v[0] + HALF*dv2[0]; velocity[1] = v[1] + HALF*dv2[1]; a = Accel(t+HALF*h, position, velocity); //a=Accel(t+HALF*h, x+HALF*dx2, v+HALF*dv2); dx3[0] = h*(v[0] + HALF*dv2[0]); dx3[1] = h*(v[1] + HALF*dv2[1]); dv3[0] = h*a[0] ; dv3[1] = h*a[1] ; position[0] = x[0] + dx3[0]; position[1] = x[1] + dx3[1]; velocity[0] = v[0] + dv3[0]; velocity[1] = v[1] + dv3[1]; a = Accel(t+h, position, velocity); //a=Accel(t+h, x+dx3, v+dv3); dx4[0] = h*(v[0] + dv3[0]); dx4[1] = h*(v[1] + dv3[1]); dv4[0] = h*a[0]; dv4[1] = h*a[1]; TrajectoryPoint p2 = new TrajectoryPoint(); p2.time = t + h; p2.position[0] = p1.position[0] + SIXTH*(dx1[0]+dx2[0]+dx2[0]+dx3[0]+dx3[0]+dx4[0]); p2.position[1] = p1.position[1] + SIXTH*(dx1[1]+dx2[1]+dx2[1]+dx3[1]+dx3[1]+dx4[1]); p2.velocity[0] = p1.velocity[0] + SIXTH*(dv1[0]+dv2[0]+dv2[0]+dv3[0]+dv3[0]+dv4[0]); p2.velocity[1] = p1.velocity[1] + SIXTH*(dv1[1]+dv2[1]+dv2[1]+dv3[1]+dv3[1]+dv4[1]); p2.acceleration=Accel(p2.time, p2.position, p2.velocity); return p2; } /* !+ SUBROUTINE Trajectory(initialAltitude, initialVelocity, initialTheta, & dt, normalized, npts, history) ! --------------------------------------------------------------------------- ! PURPOSE - Compute a trajectory, performing numerical solution of a set of ! ordinary differential equations with a fixed time step. Halt the ! calculation when the altitude is less than the initial altitude and ! correct the final point to have the same altitude as the initial altitude. */ public TrajectoryPoint[] Trajectory(double initialAltitude, double initialVelocity, double initialTheta, double dt, boolean normalized, int maxhistory) { double t; double[] position = new double[2]; double[] velocity = new double[2]; double[] acceleration = new double[2]; t=0.0; position[0]=0.0; position[1]=initialAltitude; velocity[0]=initialVelocity*Math.cos(initialTheta*PI/180); velocity[1]=initialVelocity*Math.sin(initialTheta*PI/180); acceleration=Accel(t, position, velocity) ; TrajectoryPoint[] history = new TrajectoryPoint[maxhistory]; history[0] = new TrajectoryPoint(); history[0].time=t; history[0].position=position; history[0].velocity=velocity; history[0].acceleration=acceleration; int k; for (k=1; k