Discus Trajectory

After a hiatus of a few decades, I'm thinking about again throwing the discus as a hobby:)

This page has some quick calculations of trajectory & distance for throwing a discus (in a vacuum).

I decided to do some Physics calculations to consider the effects of angle & velocity on distance:

Meter units plot for a few selected velocities (C# code)

Feet units plot for a few (different) selected velocities (C# code)

Note: C# plotting/drawing code not included at this time.

Some trajectory plots (generated via R code, below)

Pretty-printed R code screenshot for above plot

R Code for above plot (as text):

### Discus calculations and plotting

plotDiscusTraj = function() {

degToRad = pi / 180

metersToFeet = 100 / 2.54 / 12

f = function( a, b, c, x ) { return( a * x^2 + b * x + c ) }

g = 32.2

h = 1.2 * metersToFeet

theta = 40 * degToRad

v0 = c( 55, 65, 75, 85 ) ### Feet per sec

### y(t) coeff

a_t = -g / 2

b_t = v0 * sin( theta )

c_t = h

### y(x) coeff

a_x = -g / ( 2 * v0^2 * cos( theta )^2 )

b_x = rep( tan( theta ), length( v0 ) )

c_x = rep( h, length( v0 ) )

### Flight duration --> distance calculations

tf = ( -b_t - sqrt( b_t^2 - 4 * a_t * c_t ) ) / ( 2 * a_t )

dist = v0 * cos( theta ) * tf

vertexX = -b_x / ( 2 * a_x ) ### Parabola x-coord = -b/(2a)

vertexY = f( a_x, b_x, c_x, vertexX );

minX = 0; maxX = max( dist ) * 1.05 ### Axis limits

minY = 0; maxY = max( vertexY ) * 1.05

plot( 1, type="n", main="Trajectory & Distance by Velocity", xaxs="i", yaxs="i",

xlab="dist in feet", ylab="elev in feet", xlim=c(minX, maxX), ylim=c(minY, maxY), cex.lab=1.1 )

### Plot trajectory for each initial velocity

for ( i in 1:length( a_x ) ) {

x = seq( 0, dist[ i ], length.out=200 )

y = f( a_x[ i ], b_x[ i ], c_x[ i ], x )

lines( x, y )

text( vertexX[ i ], vertexY[ i ] - 2, sprintf( "%.0f f/s", v0[ i ] ) );

}

Pretty-printed R code screenshot for above plot

Here is the above code as text (no indentation):

public class Trajectory {

public enum Units { Mks, Fps };

public static readonly float degToRad = ( float ) Math.PI / 180;

public static readonly float radToDeg = 1 / degToRad;

public static readonly float metersToFeet = ( float ) ( 100 / 2.54 / 12 );

public static readonly float feetToMeters = 1 / metersToFeet;

public static readonly float gMks = 9.807f;

public static readonly float gFps = gMks * metersToFeet;

public readonly float g;

public readonly float h;

public readonly float v;

public readonly float dist;

public readonly float theta;

public readonly Units units;

public readonly int nSamples;

public readonly float a, b, c;

private readonly List<PointF> trajPoints = new List<PointF>();

private readonly PointF vertex;

// This static factory method returns an immutable instance of this class

// Internal fields are stored assuming the units specified in this invocation

// thetaDeg param is assumed to be in degrees and is converted to radians before being stored in instance.theta

public static Trajectory getInstance( float v, float h, float thetaDeg, Units units = Units.Fps, int nSamples = 200 ) {

return new Trajectory( v, h, thetaDeg, units, nSamples );

}

private Trajectory( float v, float h, float thetaDeg, Units units, int nSamples = 200 ) {

this.v = v;

this.h = h;

this.theta = thetaDeg * degToRad;

this.units = units;

this.nSamples = nSamples;

g = units == Units.Mks ? gMks : gFps;

// Calculate distance quadratic coefficients

float cosVal = ( float ) Math.Cos( theta );

a = -g / ( 2 * v * v * cosVal * cosVal );

b = ( float ) Math.Tan( theta );

c = h;

// Calculate max/throw distance

// Calculate x for which y = 0 in quadratic: y = ax^2 + bx + c using Quadratic Equation

dist = ( -b - ( float ) Math.Sqrt( b * b - 4 * a * c ) ) / ( 2 * a );

// Calculate vertex, the x-coord of vertex = -b/(2a)

float vertexX = -b / ( 2 * a );

float vertexY = a * vertexX * vertexX + b * vertexX + c;

vertex = new PointF( vertexX, vertexY );

// Calculate trajectory points

for ( int i = 0; i < nSamples; ++i ) {

float x = ( float ) i / ( float ) ( nSamples - 1 ) * ( float ) dist;

float y = ( float ) ( a * x * x + b * x + c );

trajPoints.Add( new PointF( x, y ) );

}

public PointF Vertex { get { return new PointF( vertex.X, vertex.Y ); } }

public PointF[] TrajPoints { get { return trajPoints.ToArray(); } } // Return array of PointF copies

private static string ptfToStr( PointF pt ) {

return "( " + pt.X.ToString( "0.00" ) + ", " + pt.Y.ToString( "0.00" ) + " )";

}

public string summStr() {

string unitStr = units == Units.Fps ? "f" : "m";

return "v0: " + v.ToString( "0." ) + " " + unitStr + "/s" +

" dist: " + dist.ToString( "0.0" ) + " " + unitStr +

" angle: " + theta * Trajectory.radToDeg + " deg";

}

public override string ToString() {

string distUnits = units == Units.Fps ? "feet" : "meters";

StringBuilder sb = new StringBuilder();

sb.AppendLine( "Trajectory object dump:" );

sb.AppendLine( " h: " + h + " " + distUnits );

sb.AppendLine( " v: " + v + " " + distUnits + "/sec" );

sb.AppendLine( " dist: " + dist + " " + distUnits );

sb.AppendLine( " vertex: " + vertex + " (" + distUnits + ")" );

sb.AppendLine( " theta: " + theta / degToRad );

sb.AppendLine( " nSamples: " + nSamples );

sb.AppendLine( " units: " + units.ToString() );

sb.AppendLine( " g: " + g + " " + distUnits + "/sec^2" );

sb.AppendLine( " a: " + a );

sb.AppendLine( " b: " + b );

sb.AppendLine( " c: " + c );

sb.AppendLine( " trajPoints:" );

for ( int i = 0; i < trajPoints.Count; ++i )

sb.AppendLine( " trajPoints[" + i + "]: " + ptfToStr( trajPoints[ i ] )+ " " + distUnits );

return sb.ToString();

}