#!/usr/bin/perl -w
# ballistics.pl, 2006-05-02
# This software does do many little steps to calculate the flight path of
# a projectile. It implements following effects of nature:
# - air drag (constant and curves similar to Achenbach)
# - atmospheric density changes according to height (is hopefully good
#   enough for slinging heights)
# - gravity
# I am not an expert in these things, but the results are very close to
# the data given by Baatz in his article, so I think they are at least
# good enough for analyzing slinging physics. You can try how light clay
# bullets fly or heavy stones. You can try how far the things fly when
# you get in the low area shown by Achenbach. And there are many other
# things to try.
#
# Usage:
# 1) Define your projectile. Necessary is the cw (in English 'cd'?)
#    value, the mass and the surf(ace). Their units are [], [kg] and [m^2].
#    You can give it also a descr(iption) so you know from what
#    object these values came.
# 2) Edit the parameters after line 'PARAMETERS'.
#    speed: Starting speed of projectile [m/s]. Suitable speeds are around
#           30-120. Baatz assumed 75.
#    angle: Starting direction above horizon [degree]. For furthest distance
#           try values between 35 and 45. Baatz used 40. Flat pathes with
#           Apache style could be e.g. 5 or 10.
#    bottom: Calculation stops when heigh is below that value. [m] Value
#            could be 0 for underhand release or -2 for Apach style or
#            -100 if you sling down a cliff.
#    projname: The tag of the projectile to sling.
#    dragmodel: 0: A constant cw is used.
#               1: A simple model of cw changes depending on speed modelled
#                  like plots from Achenbach
# 3) Optionally edit the less important parameters in the state structure
#    found after line 'STATE'.
#    stepsize: Time for each calculation step. [s]
#    grav: Earth gravity. [m/s^2]
#    d_air: Air density at starting point. [kg/m^3] This value is ajusted
#           according to height changes.
# 4) Optionally edit airdrag1() or add other air drag functions.
# 5) Run this perl script.
# 6) Look at the output. The output columns have following meaning:
#    1: Iteration number. Multiplied by stepsize results in the actual time
#       passed.
#    2: Total speed of projectile. [m/s]
#    3: Speed of projectile in horizontal (throwing) direction. [m/s]
#    4: Speed of projectile in vertical (up) direction. [m/s]
#    5: Position of projectile in horizontal (throwing) direction. [m]
#    6: Position of projectile in vertical (up) direction. [m]
#    7: Current air drag. [N]
#    You can plot e.g. with gnuplot column 5 against 6 resulting in
#    the flight pathway:
#    gnuplot> plot 'lead' u 5:6 w l
#    Or you can plot column 1 against 2 showing the development of speed
#    over time.

# Projectile
$projectile{ blackhole } = {
        descr => 'Really dense matter',
		cw => 0, # []
		mass => 0.04, # [kg]
		dens => 1100.3, # [kg/m^3]
		surf => (0.01/10000), # [m^2]
    };

$projectile{ leadball1 } = {
        descr => 'Really smooth lead ball',
		cw => 0.48, # []
		mass => 0.04, # [kg]
		dens => 11.3, # [kg/m^3]
		surf => (2.8/10000), # [m^2]
    };

$projectile{ leadball2 } = {
        descr => 'Smooth lead ball',
		cw => 0.40, # []
		mass => 0.04, # [kg]
		dens => 11.3, # [kg/m^3]
		surf => (2.8/10000), # [m^2] r=0.9454 -> s=2.8079
    };

$projectile{ leadball3 } = {
        descr => 'Smooth heavy lead ball',
		cw => 0.48, # []
		mass => 0.2, # [kg]
		dens => 11.3, # [kg/m^3]
		surf => (3.7/10000), # [m^2]
    };

$projectile{ stoneball1 } = {
        descr => 'Really smooth stone ball',
		cw => 0.48, # []
		mass => 0.04, # [kg]
		dens => 2.7, # [kg/m^3]
		surf => (7.3/10000), # [m^2]
    };

$projectile{ stoneball2 } = {
        descr => 'Smooth stone ball',
		cw => 0.40, # []
		mass => 0.04, # [kg]
		dens => 2.7, # [kg/m^3]
		surf => (7.3/10000), # [m^2]
    };

$projectile{ clayball1 } = {
        descr => 'Really smooth clay ball',
		cw => 0.48, # []
		mass => 0.04, # [kg]
		dens => 1.8, # [kg/m^3]
		surf => (9.6/10000), # [m^2]
    };

$projectile{ clayball2 } = {
        descr => 'Smooth clay ball',
		cw => 0.40, # []
		mass => 0.04, # [kg]
		dens => 1.8, # [kg/m^3]
		surf => (9.6/10000), # [m^2]
    };

$projectile{ clayball3 } = {
        descr => 'Smooth light clay ball',
		cw => 0.40, # []
		mass => 0.027, # [kg]
		dens => 1.8, # [kg/m^3]
		surf => (7.3/10000), # [m^2]
    };

$projectile{ marble } = {
        descr => 'marble',
		cw => 0.3, # []
		mass => 0.05, # [kg]
		dens => 2.7, # [kg/m^3]
		surf => (5.0/10000), # [m^2]
    };

# PARAMETERS
$speed=75; # [m/s]
$angle=40; # [degree]
$bottom=0; # [m]
$projname='leadball2';
$dragmodel=0; # 0: constant cw, 1: modified cw similar to Achenbach

# STATE
$state = {
		stepsize => 0.001, # [s]
		grav => 9.81, # [m/s^2]
		d_air => 1.14, # [kg/m^3]
		vx => cos(0.5*3.1415927*($angle/90))*$speed,
		vy => sin(0.5*3.1415927*($angle/90))*$speed,
		vges => $speed,
		drag => 0,
		px => 0,
		py => 0,
		projectile => $projectile{$projname},
		dragmodel => $dragmodel,
	};

# Air drag modeled to Achenbach
sub airdrag1 {
	my ($cw, $v) = @_;

	#$v = ($v/60)*50;

	if( $v<50 ) {
		# High cw at low speed
		return $cw*(0.9+$v/250);
	} else {
		return $cw*(0.1+($v-50)/170);
	}
}

# Air drag for given parameters in N.
# $cw [], $d_air [kg/m^3], $v [m/s], $q [m^2], $height [m]
# $opt has currently no effect, but could be in the future.
sub airdrag {
	my ($cw, $d_air, $v, $q, $height, $opt) = @_;

	if( $opt==1 ) {
		$cw=airdrag1($cw,$v);
	}

	return $cw * 0.5 * $d_air * exp(-0.000106*$height) * ($v**2) * $q;
}

# Modify position according to speed 
sub modpos {
	my $state = shift;
	$state->{px} += $state->{vx} * $state->{stepsize};
	$state->{py} += $state->{vy} * $state->{stepsize};
}

sub sgn {
	my $Wert = shift;
	return $Wert <=> 0;
}

# Modify speed
sub modspeed {
	my $state = shift;

	# direction of movement
	$dir = atan2($state->{vy},$state->{vx});

	# Total speed
	$vges = sqrt($state->{vx}**2+$state->{vy}**2);
	$state->{vges} = $vges;

	# Total air drag
	$drag = airdrag( $state->{projectile}->{cw}, $state->{d_air}, $vges,
		$state->{projectile}->{surf}, $state->{py}, $state->{dragmodel} );
	$state->{drag} = $drag;

	# Separate drags
	$dragx = cos($dir)*$drag;
	$dragy = sin($dir)*$drag;

	$state->{vx} += ((-1)*$dragx/$state->{projectile}->{mass}) *
		$state->{stepsize};

	$state->{vy} += (((-1)*$dragy/$state->{projectile}->{mass}) - $state->{grav}) *
		$state->{stepsize};
}

# main

for( $i=0; $state->{py}>=$bottom; $i++ ) {
	printf( "%d %f %f %f %f %f %f\n", $i, $state->{vges}, $state->{vx}, $state->{vy}, $state->{px}, $state->{py}, $state->{drag} );
	modpos( $state );
	modspeed( $state );
}
printf( "%d %f %f %f %f %f %f\n", $i, $state->{vges}, $state->{vx}, $state->{vy}, $state->{px}, $state->{py}, $state->{drag} );

