#!/usr/bin/perl
# STS-107 Columbia debris trajectory simulation
# Copyright (C) 2003  Ian Kluft
#    derived from open source works of Robert L. McCoy, US Army Ballistics Research Lab
#    
# ported to Perl by Ian Kluft from mctraj.bas program for bullet trajectories
# original author unknown
# * PROGRAM [MCTRAJ.BAS], MAR. 1987. [TANDY 1000 HX, 07/90, REV. 02/93]
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, write to the Free Software
#    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

use strict;
use Getopt::Long;

# POINT MASS TRAJECTORIES FOR SMALL ARMS
# THE PROGRAM REQUIRES AN INPUT TABLE OF
# DRAG COEFFICIENT (CD) VERSUS MACH NUMBER (M).
# ADDITIONAL REQUIRED INPUTS ARE:
# MUZZLE VELOCITY (FT/SEC); BALLISTIC COEFFICIENT (LB/IN 2);
# HEIGHT OF SIGHT LINE ABOVE BORE CENTERLINE (INCHES);
# GUN ELEVATION ANGLE (MINUTES); RATIO AIR DENSITY TO
# STANDARD DENSITY; AIR TEMPERATURE (DEG F); RANGE PRINT
# INTERVAL (YARDS/METERS); RANGE TO TERMINATE TRAJECTORY (YARDS/METERS);
# RANGE WIND SPEED (MPH--POSITIVE IF WIND BLOWS FROM GUN TO
# TARGET); AND CROSSWIND SPEED (MPH--POSITIVE IF WIND BLOWS
# FROM LEFT TO RIGHT ACROSS FIRING LINE).
# THE PROGRAM ALSO PROVIDES AN OPTION TO ADJUST THE TRAJECTORY
# TO PASS THROUGH A SPECIFIED POINT IN SPACE, DENOTED BY
# RMATCH (YARDS/METERS), AND HMATCH (INCHES).  THE GUN ELEVATION ANGLE
# IS ADJUSTED UNTIL THE TRAJECTORY PASSES THROUGH THE POINT
# (RMATCH,HMATCH).  IF NO ADJUSTMENT IS DESIRED, INPUT ZEROS
# FOR RMATCH AND HMATCH, AND THE TRAJECTORY WILL BE RUN WITH
# THE INPUT ELEVATION ANGLE.
# THE PROGRAM SOLVES THE TRAJECTORY BY NUMERICAL INTEGRATION
# USING THE HEUN METHOD, WHICH IS AN ITERATIVELY APPLIED
# SECOND ORDER PREDICTOR-CORRECTOR TECHNIQUE.
#
# THE TRAJECTORY OUTPUT IS RANGE (YARDS/METERS); HEIGHT RELATIVE TO
# LINE OF SIGHT (INCHES); DEFLECTION (INCHES); TOTAL VELOCITY
# (FT/SEC); TIME OF FLIGHT (SECONDS); RANGE COMPONENT OF
# VELOCITY (VX--FT/SEC); VERTICAL COMPONENT OF VELOCITY
# (VY--FT/SEC); AND HORIZONTAL COMPONENT OF VELOCITY (VZ--FT/SEC).

# variable conversion notes:
# G  - $gravity   acceleration due to grapvity 32.174 ft/sec^2
# E1 - deleted    constant near-zero floating point number
# U1 - $fmt       string for BASIC's PRINT USING, converted to printf format
# K  - deleted    arbitrary name of drag function
# M  - $cd[$x][0] mach number in drag coefficient table
# D  - $cd[$x][1] drag coefficient (at a mach number) in drag coefficient table
# YM - deleted    flag for displaying range in yard or meters (use metric)
# SI - deleted    flag for computation by US Army or ICAO methods (use ICAO)
# K2 - projectile string for projectile identification
# V0 - velocity0  originally muzzle velocity, ft/sec
# C  - b_coeff    ballistic coefficient
# H0 - height0    initial height above ground
# P0 - angle0     initial projectile angle (minutes)
# R0 - air_density ratio of air density to sea level
# TDF- temp       temperature (F)
# W1 - wind_range downrange wind velocity vector
# W3 - wind_cross crosswind velocity vector
# R8 - deleted    range to point that projectile must pass through
# H8 - deleted    height of point that projectile must pass through
# P2 - print_step range to print each line
# R3 - deleted    range to end printing (simulation needs to run to impact)
# PR2- range_dist downrange distance, also index for table output

# TODO:
# convert air density ratio into a table from winds aloft archive
our $air_density = 1;
# convert temperature into a table from winds aloft archive
our $temp = -40;
# convert wind direction into a table from winds aloft archive
our $wind_range = 0;
our $wind_cross = 0;

# define program constants
our $gravity = 32.174;

our @cd = (
	# GS ballistic model for a small sphere
	# data from  http://internet.cybermesa.com/~jbm/downloads/text/mcgs.txt
	[ 0.00, 0.4662 ], [ 0.05, 0.4689 ], [ 0.10, 0.4717 ], [ 0.15, 0.4745 ],
	[ 0.20, 0.4772 ], [ 0.25, 0.4800 ], [ 0.30, 0.4827 ], [ 0.35, 0.4852 ],
	[ 0.40, 0.4882 ], [ 0.45, 0.4920 ], [ 0.50, 0.4970 ], [ 0.55, 0.5080 ],
	[ 0.60, 0.5260 ], [ 0.65, 0.5590 ], [ 0.70, 0.5920 ], [ 0.75, 0.6258 ],
	[ 0.80, 0.6610 ], [ 0.85, 0.6985 ], [ 0.90, 0.7370 ], [ 0.95, 0.7757 ],
	[ 1.00, 0.8140 ], [ 1.05, 0.8512 ], [ 1.10, 0.8870 ], [ 1.15, 0.9210 ],
	[ 1.20, 0.9510 ], [ 1.25, 0.9740 ], [ 1.30, 0.9910 ], [ 1.35, 0.9990 ],
	[ 1.40, 1.0030 ], [ 1.45, 1.0060 ], [ 1.50, 1.0080 ], [ 1.55, 1.0090 ],
	[ 1.60, 1.0090 ], [ 1.65, 1.0090 ], [ 1.70, 1.0090 ], [ 1.75, 1.0080 ],
	[ 1.80, 1.0070 ], [ 1.85, 1.0060 ], [ 1.90, 1.0040 ], [ 1.95, 1.0025 ],
	[ 2.00, 1.0010 ], [ 2.05, 0.9990 ], [ 2.10, 0.9970 ], [ 2.15, 0.9956 ],
	[ 2.20, 0.9940 ], [ 2.25, 0.9916 ], [ 2.30, 0.9890 ], [ 2.35, 0.9869 ],
	[ 2.40, 0.9850 ], [ 2.45, 0.9830 ], [ 2.50, 0.9810 ], [ 2.55, 0.9790 ],
	[ 2.60, 0.9770 ], [ 2.65, 0.9750 ], [ 2.70, 0.9730 ], [ 2.75, 0.9710 ],
	[ 2.80, 0.9690 ], [ 2.85, 0.9670 ], [ 2.90, 0.9650 ], [ 2.95, 0.9630 ],
	[ 3.00, 0.9610 ], [ 3.05, 0.9589 ], [ 3.10, 0.9570 ], [ 3.15, 0.9555 ],
	[ 3.20, 0.9540 ], [ 3.25, 0.9520 ], [ 3.30, 0.9500 ], [ 3.35, 0.9485 ],
	[ 3.40, 0.9470 ], [ 3.45, 0.9450 ], [ 3.50, 0.9430 ], [ 3.55, 0.9414 ],
	[ 3.60, 0.9400 ], [ 3.65, 0.9385 ], [ 3.70, 0.9370 ], [ 3.75, 0.9355 ],
	[ 3.80, 0.9340 ], [ 3.85, 0.9325 ], [ 3.90, 0.9310 ], [ 3.95, 0.9295 ],
	[ 4.00, 0.9280 ],

	# make a wild guess on increased drag at hypersonic speeds
	# these two set the rate for extrapolation above Mach 6
	[ 5.0, 0.94042491 ],
	#[ 5.0, 1.0 ],
	#[ 6.0, 1.2 ], 
);

# function for interpolation in drag table
# input mach number
# output drag coefficient
sub interpolate_drag
{
	my ( $mach ) = @_;
	my ( $i, $pos, $cd );

	$i = 0;
	while (( $mach > $cd[$i][0]) and $i < $#cd ) {
		$i++;
	}
	if ( !defined $cd[$i+1]) {
		$i--;
	}
	#$pos = ($cd[$i+1][1] - $cd[$i][1]) / ($cd[$i+1][0] - $cd[$i][0]);
	$pos = ($mach - $cd[$i][0]) / ($cd[$i+1][0] - $cd[$i][0]);
	#$cd = $cd[$i][1] + $pos * ( $mach - $cd[$i][0]);
	$cd = $cd[$i][1] + $pos * ( $cd[$i+1][1] - $cd[$i][1]);
	#print "debug: M=$mach -> CD=$cd\n";

	return $cd;
}

#
# main
#

# defaults for command-line arguments
our ( $projectile, $velocity0, $height0, $angle0, $print_step, $b_coeff );
$projectile = "space shuttle tile";
$velocity0 = 24370.14;    # Mach 22.41
$height0 =  226748;       #  226748 ft altitude
$angle0 = -9.966;         # angle = -9.966 minutes (15.34 ft descent per mile)
$print_step = 1;          # print step = 1 mile
$b_coeff = .3744153/36.0; # ballistic coefficient lbs/in^2

# process command-line arguments
if ( ! GetOptions (
	"projectile:s" => \$projectile,
	"velocity:s" => \$velocity0,
	"height:s" => \$height0,
	"angle:s" => \$angle0,
	"bc:s" => \$b_coeff,
	))
{
	die "usage: $0 --projectile=projectile-id\n";
}

# Initialize computational parameters
my $UC = 5280;  # range displayed in miles
my $DINT = 1;
my $D3 = $DINT * $UC;

# ICAO standard settings
my $RH1 = -.00002926;
my $RH2 = -.0000000001;
my $TK1 = -.000006858;
my $TK2 = -.00000000002776;
# Note: PIR = -(PI/8)*(RHO0/144)
my $PIR = -.000208551;
my $VV1 = 49.0223;

# initialize input values for trajectory calculation
my $J = 0;
my $P1 = $wind_range;
my $P3 = $wind_cross;
my $R4 = $D3 * $print_step;
my $R6 = 0;
my $PR4 = $print_step;
my $PR6 = 0;
my $wind_range = (22 * $wind_range) / 15;
my $wind_cross = (22 * $wind_cross) / 15;
my $done = 0;

# print table heading
our $fmt = "%6s %9.1f     %7.1f   %7.1f  %7.3f  %7.1f  %7.1f %6.1f\n";
printf "%6s %9s     %7s   %7s  %7s  %7s  %7s %6s\n",
	"RANGE", "HEIGHT", "DEFL.", "MACH", "TIME", "VX", "VY", "VZ";
printf "%6s %9s     %7s   %7s  %7s  %7s  %7s %6s\n",
	"(MILES)", "(FT)", "(FT)", "(FPS)", "(SEC)", "(FPS)", "(FPS)", "(FPS)";

my $V3 = $velocity0 * cos($angle0 / 3437.74677);
my $V4 = $velocity0 * sin($angle0 / 3437.74677);
my $V5 = 0;
my $R1 = 0;
my $PR1 = 0;
my $H1 = $height0;
my $D1 = 0;
my $T1 = 0;
my $PX3 = $PR4;
my @E;
my @H;

while ( !$done ) {
	# begin trajectory calculation
	my $C3 = ($PIR * $air_density) / $b_coeff;
	my $B1 = sqrt(($V3 - $wind_range) ** 2 + $V4 ** 2 + ($V5 - $wind_cross) ** 2);
	my $T0 = ($temp + 459.67) * exp(($TK1 + $TK2 * $H1) * $H1) - 459.67;
	my $V1 = $VV1 * sqrt($T0 + 459.67);

	my $C1 = interpolate_drag( $B1 / $V1 );

	my $C4 = ($C3 * $C1 * $B1 * exp(($RH1 + $RH2 * $H1) * $H1)) / $V3;
	my $A1 = $C4 * ($V3 - $wind_range);
	my $A2 = $C4 * $V4 - $gravity / $V3;
	my $A3 = $C4 * ($V5 - $wind_cross);
	my $range_dist;
	my $H2;

	do {
		# apply Euler predictor formula
		my $R2 = $R1 + $D3;
		$range_dist = $PR1 + $DINT;
		my $V6 = $V3 + $A1 * $D3;
		my $V7 = $V4 + $A2 * $D3;
		my $V8 = $V5 + $A3 * $D3;
		my $B2 = sqrt(($V6 - $wind_range) ** 2 + $V7 ** 2 + ($V8 - $wind_cross) ** 2);
		my $A4;
		my $A5;
		my $A6;
		my $E2;

		# apply iterative Heun corrector formula
		do {
			my $U = $B2;
			$T0 = ($temp + 459.67) * exp(($TK1 + $TK2 * $H1) * $H1) - 459.67;
			$V1 = $VV1 * sqrt($T0 + 459.67);

			my $C2 = interpolate_drag( $B2 / $V1 );

			my $C5 = ($C3 * $C2 * $B2 * exp(($RH1 + $RH2 * $H1) * $H1)) / $V6;
			$A4 = $C5 * ($V6 - $wind_range);
			$A5 = $C5 * $V7 - $gravity / $V6;
			$A6 = $C5 * ($V8 - $wind_cross);
			$V6 = $V3 + .5 * ($A1 + $A4) * $D3;
			$V7 = $V4 + .5 * ($A2 + $A5) * $D3;
			$V8 = $V5 + .5 * ($A3 + $A6) * $D3;
			$B2 = sqrt(($V6 - $wind_range) ** 2 + $V7 ** 2
				+ ($V8 - $wind_cross) ** 2);
			$E2 = abs(($B2 - $U) / $B2);
		} while $E2 > .00001;


		# compute values at R2
		$H2 = $H1 + (($V4 + $V7) / ($V3 + $V6)) * $D3;
		my $D2 = $D1 + (($V5 + $V8) / ($V3 + $V6)) * $D3;
		my $T2 = $T1 + (2 * $D3) / ($V3 + $V6);

		# reset conditions at R1 to new conditions at R2
		$R1 = $R2;
		$PR1 = $range_dist;
		$H1 = $H2;
		$D1 = $D2;
		$T1 = $T2;
		$V3 = $V6;
		$V4 = $V7;
		$V5 = $V8;
		$A1 = $A4;
		$A2 = $A5;
		$A3 = $A6;

		# check print conditions
		if ( $range_dist >= $PX3 ) {
			printf $fmt, $range_dist, $H2, $D2, $V1, $T2,
				$V6, $V7, $V8;
		}

		# INCREMENT PRINT RANGE.
		$PX3 = $PX3 + $PR4;
	} while ( $H2 >= 0 );

	# INTERATION TO ADJUST ELEVATION FOR MATCH RANGE AND HEIGHT.
	if ( $H2 < .00001 ) {
		# final trajectory
		$done = 1;
	} else {
		$J = $J + 1;
		$E[$J] = $angle0;
		$H[$J] = $H2;
		if ( $J >= 2 ) {
			# adjust elevation angle
			$E[$J + 1] = $E[$J] - $H[$J] * ($E[$J - 1]
				- $E[$J]) / ($H[$J - 1] - $H[$J]);
		} else {
			$E[$J + 1] = $angle0 + 2;
		}
		$angle0 = $E[$J + 1];
	}
}