/*
rifled.cpp   v1.05Linux

A C program to calculate external rifle ballistics.
Inspired by an article titled 'Ballistics on the Home Computer'
by William C. Davis Jr. published in American Rifleman, June 1983.

Features:
   -Tabular format attempts to answer basic rifle sighting questions.
   -Empirical look-up table method is used for accurate velocity results.
   -Temperature and altitude effects are calculated.

Disclaimer:
No warranty is expressed or implied. Use at your own risk.
(c)copyright 1995 by Ross A. Winters, all rights reserved.

Revision History
v1.00     Original Release. Compiled in Turbo C++ v3.0 for 8088/8086 processor
		  w/floating point emulation, RAW, 8/26/95
v1.02     Add max. range and range increment. RAW 11/1/95
v1.03     Add comma delimited file output, RAW 11/15/95
v1.04     Edited for Linux/CGI interactive web page. RAW 9/17/99
          usage: "rifle MZVEL WEIGHT HEIGHT COEF TGTRANGE ALTITUDE TEMP MXRANGE DELTAX >FILE.OUT"
*/

#include <stdio.h>
#include <stdlib.h> 
#include <math.h>

/* Input variable ranges, some are 'best guess' */
#define V_MAX 4500
#define V_MIN 400
#define WT_MIN 1
#define WT_MAX 500
#define H_MIN 0.001
#define H_MAX 6.0
#define R0_MIN 50
#define R0_MAX 1000
#define BC_MIN 0.001
#define BC_MAX 5
#define AL_MIN -500
#define AL_MAX 25000
#define TF_MIN -60
#define TF_MAX 120
#define X_MIN 0
#define X_MAX 1000
#define DX_MIN 1
#define DX_MAX 100

/* return 2 range error types */
#define RANGE_ERR_INPUT 1  /* Input data range errors */
#define RANGE_ERR_FLIGHT -1 /* Sometimes velocity range errors occur during flight. */

double v,g,h,c1,r0,al,tf; /* ballistic related input parameters */
double mx,dx; /* table max range and delta x input parameters */
double vb,as,bs,cs,at,bt,ct; /* ballistic calculation variables */
double sv,tv,su,x,f,af,df,rf,ro,u,en,tu,t,d,ym,hm,e,e0; /* ballistic calculation variables */

void zero_angle();
int ballistics(),do_table1(double v1),do_table2(double v2);

/* functions */
int do_table1(double v1)
{
if(v1>=400 && v1<500){vb=450.0;as=33958.1;bs=-51.756;cs=.0473157;at=35.269;bt=-0.115879;ct=2.34824E-4;return 0;};
if(v1>=500 && v1<600){vb=550.0;as=29218.7;bs=-43.640;cs=.035161;at=25.733;bt=-0.079629;ct=1.36695E-4;return 0;};
if(v1>=600 && v1<700){vb=650.0;as=25192.0;bs=-37.095;cs=.0311391;at=18.9904;bt=-0.057305;ct=9.20586E-5;return 0;};
if(v1>=700 && v1<800){vb=750.0;as=21792.0;bs=-30.896;cs=.0313196;at=14.1144;bt=-0.041349;ct=6.9381E-5;return 0;};
if(v1>=800 && v1<900){vb=850.0;as=19020.1;bs=-24.4588;cs=.0330303;at=10.6366;bt=-0.0288527;ct=5.58719E-5;return 0;};
if(v1>=900 && v1<1000){vb=950.0;as=16906.5;bs=-17.8388;cs=.0325134;at=8.27886;bt=-0.01884;ct=4.41677E-5;return 0;};
if(v1>=1000 && v1<1050){vb=1025.0;as=15747.2;bs=-13.2423;cs=.0276581;at=7.10211;bt=-0.0129298;ct=3.32035E-5;return 0;};
if(v1>=1050 && v1<1075){vb=1062.5;as=15288.1;bs=-11.3183;cs=.0233460;at=6.66207;bt=-0.0106524;ct=2.68757E-5;return 0;};
if(v1>=1075 && v1<1100){vb=1087.0;as=15019.2;bs=-10.2383;cs=.0199585;at=6.41183;bt=-9.41688E-3;ct=2.2611E-5;return 0;};
if(v1>=1100 && v1<1110){vb=1105.0;as=14846.0;bs=-9.58136;cs=.0164795;at=6.25378;bt=-8.66849E-3;ct=2.0396E-5;return 0;};
if(v1>=1110 && v1<1120){vb=1115.0;as=14751.9;bs=-9.23136;cs=.0165939;at=6.16907;bt=-8.28369E-3;ct=1.83284E-5;return 0;};
if(v1>=1120 && v1<1130){vb=1125.0;as=14661.2;bs=-8.93002;cs=.0147629;at=6.08802;bt=-7.93184E-3;ct=1.70991E-5;return 0;};
if(v1>=1130 && v1<1150){vb=1140.0;as=14530.6;bs=-8.48200;cs=.0116909;at=5.97274;bt=-7.45342E-3;ct=1.50953E-5;return 0;};
if(v1>=1150 && v1<1250){vb=1200.0;as=14062.5;bs=-7.25680;cs=8.50001E-3;at=5.57228;bt=-6.058E-3;ct=9.6174E-6;return 0;};
if(v1>=1250 && v1<1500){vb=1375.0;as=12977.5;bs=-5.54640;cs=3.06207E-3;at=4.72557;bt=-4.05656E-3;ct=3.72586E-6;return 0;};
if(v1>=1500 && v1<2000){vb=1750.0;as=11189.1;bs=-4.31589;cs=9.60342E-4;at=3.56969;bt=-2.49074E-3;ct=1.27366E-6;return 0;};
if(v1>=2000 && v1<2500){vb=2250.0;as=9222.06;bs=-3.66916;cs=4.54010E-4;at=2.57801;bt=-1.63915E-3;ct=5.68494E-7;return 0;};
if(v1>=2500 && v1<3000){vb=2750.0;as=7492.73;bs=-3.26733;cs=3.71407E-4;at=1.88267;bt=-1.1924E-3;ct=3.52828E-7;return 0;};
if(v1>=3000 && v1<3500){vb=3250.0;as=5950.45;bs=-2.90704;cs=3.48885E-4;at=1.36656;bt=-8.97071E-4;ct=2.45908E-7;return 0;};
if(v1>=3500 && v1<4000){vb=3750.0;as=4582.38;bs=-2.57414;cs=3.10524E-4;at=.974457;bt=-6.88003E-4;ct=1.74925E-7;return 0;};
if(v1>=4000 && v1<=4500){vb=4250.0;as=3369.09;bs=-2.29123;cs=2.57709E-4;at=.670381;bt=-5.40083E-4;ct=1.24357E-7;return 0;};
}

int do_table2(double v2)
{
if(v2<=43041.0 && v2>36664.2){vb=350.0;as=39663.0;bs=-63.768;cs=.04758391;at=49.669;bt=-0.184500;ct=4.82396E-4;return 0;};
if(v2<=36664.2 && v2>31488.6){vb=450.0;as=33958.1;bs=-51.756;cs=.0473157;at=35.269;bt=-0.115879;ct=2.34824E-4;return 0;};
if(v2<=31488.6 && v2>27124.6){vb=550.0;as=29218.7;bs=-43.640;cs=.035161;at=25.733;bt=-0.079629;ct=1.36695E-4;return 0;};
if(v2<=27124.6 && v2>23415.1){vb=650.0;as=25192.0;bs=-37.095;cs=.0311391;at=18.9904;bt=-0.057305;ct=9.20586E-5;return 0;};
if(v2<=23415.1 && v2>20325.5){vb=750.0;as=21792.0;bs=-30.896;cs=.0313196;at=14.1144;bt=-0.041349;ct=6.9381E-5;return 0;};
if(v2<=20325.5 && v2>17879.9){vb=850.0;as=19020.1;bs=-24.4588;cs=.0330303;at=10.6366;bt=-0.0288527;ct=5.58719E-5;return 0;};
if(v2<=17879.9 && v2>16095.6){vb=950.0;as=16906.5;bs=-17.8388;cs=.0325134;at=8.27886;bt=-0.01884;ct=4.41677E-5;return 0;};
if(v2<=16095.6 && v2>15433.3){vb=1025.0;as=15747.2;bs=-13.2423;cs=.0276581;at=7.10211;bt=-0.0129298;ct=3.32035E-5;return 0;};
if(v2<=15433.3 && v2>15150.3){vb=1062.5;as=15288.1;bs=-11.3183;cs=.0233460;at=6.66207;bt=-0.0106524;ct=2.68757E-5;return 0;};
if(v2<=15150.3 && v2>14894.3){vb=1087.0;as=15019.2;bs=-10.2383;cs=.0199585;at=6.41183;bt=-9.41688E-3;ct=2.2611E-5;return 0;};
if(v2<=14894.3 && v2>14798.5){vb=1105.0;as=14846.0;bs=-9.58136;cs=.0164795;at=6.25378;bt=-8.66849E-3;ct=2.0396E-5;return 0;};
if(v2<=14798.5 && v2>14706.2){vb=1115.0;as=14751.9;bs=-9.23136;cs=.0165939;at=6.16907;bt=-8.28369E-3;ct=1.83284E-5;return 0;};
if(v2<=14706.2 && v2>14616.9){vb=1125.0;as=14661.2;bs=-8.93002;cs=.0147629;at=6.08802;bt=-7.93184E-3;ct=1.70991E-5;return 0;};
if(v2<=14616.9 && v2>14447.0){vb=1140.0;as=14530.6;bs=-8.48200;cs=.0116909;at=5.97274;bt=-7.45342E-3;ct=1.50953E-5;return 0;};
if(v2<=14447.0 && v2>13720.5){vb=1200.0;as=14062.5;bs=-7.25680;cs=8.50001E-3;at=5.57228;bt=-6.058E-3;ct=9.6174E-6;return 0;};
if(v2<=13720.5 && v2>12330.3){vb=1375.0;as=12977.5;bs=-5.54640;cs=3.06207E-3;at=4.72557;bt=-4.05656E-3;ct=3.72586E-6;return 0;};
if(v2<=12330.3 && v2>10168.1){vb=1750.0;as=11189.1;bs=-4.31589;cs=9.60342E-4;at=3.56969;bt=-2.49074E-3;ct=1.27366E-6;return 0;};
if(v2<=10168.1 && v2>8332.83){vb=2250.0;as=9222.06;bs=-3.66916;cs=4.54010E-4;at=2.57801;bt=-1.63915E-3;ct=5.68494E-7;return 0;};
if(v2<=8332.83 && v2>6699.05){vb=2750.0;as=7492.73;bs=-3.26733;cs=3.71407E-4;at=1.88267;bt=-1.1924E-3;ct=3.52828E-7;return 0;};
if(v2<=6699.05 && v2>5245.45){vb=3250.0;as=5950.45;bs=-2.90704;cs=3.48885E-4;at=1.36656;bt=-8.97071E-4;ct=2.45908E-7;return 0;};
if(v2<=5245.45 && v2>3958.11){vb=3750.0;as=4582.38;bs=-2.57414;cs=3.10524E-4;at=.974457;bt=-6.88003E-4;ct=1.74925E-7;return 0;};
if(v2<=3958.11 && v2>2812.29){vb=4250.0;as=3369.09;bs=-2.29123;cs=2.57709E-4;at=.670381;bt=-5.40083E-4;ct=1.24357E-7;return 0;};
return -1;
}

int ballistics() /* Ballistics Calculation Routine */
{
int y=0;

   af=1.00-3.59596E-5*al+4.7741E-10*pow(al,2);
   rf=518.67/(459.67+tf);
   ro=af*rf;
   x=0;
   do_table1(v);
   sv=as+bs*(v-vb)+cs*pow(v-vb,2);
   tv=at+bt*(v-vb)+ct*pow(v-vb,2);
   zero_angle(); /* find the aiming angle */
   do {
	  su=sv+3*x/(c1/ro);
	  if(do_table2(su)) return RANGE_ERR_FLIGHT; /* trap velocity errors */
	  u=vb+(-bs-sqrt(pow(bs,2)-4*cs*(as-su)))/(2*cs); /* velocity, speed */
	  en=pow(u,2)*g/450400; /* energy */
	  f=14.0069+6.59285*(u/v-.65)-1.94051*pow((u/v-.65),2);
	  tu=at+bt*(u-vb)+ct*pow(u-vb,2);
	  t=(c1/ro)*(tu-tv); /* time of flight */
	  d=12*f*pow(t,2); /* total drop */
	  ym=12*4.05*pow(t,2); /* estimated elev. max, originally k=4.05 */
	  df=176*(t-3*x/v); /* 10mph x-wind deflection */
	  e=(x*e0-d)+h; /* elevation, y trajectory */
	  hm=ym+0.4*h;
	  if (x==0) {hm=h;};

	  printf ("%5g\t%+.2f\t% .0f\t% .0f\t%.3f\t%.2f\t%+.2f\n",x,e,u,en,t,d,hm,df);

	  x=x+dx;
      }
   while (x<=mx);
return 0;
}

void zero_angle()
/* figure how much drop occurs, then compensate for the aiming angle */
{
	  su=sv+3*r0/(c1/ro);
	  do_table2(su);
	  u=vb+(-bs-sqrt(pow(bs,2)-4*cs*(as-su)))/(2*cs);
	  f=14.0069+6.59285*(u/v-.65)-1.94051*pow((u/v-.65),2);
	  tu=at+bt*(u-vb)+ct*pow(u-vb,2);
	  t=(c1/ro)*(tu-tv);
	  d=12*f*pow(t,2);
	  if (r0==0){e0=0;}else{e0=(d-h)/r0;}
}

/* the main */
int main(int argc, char *argv[])
{
v=atof(argv[1]);
g=atof(argv[2]);
h=atof(argv[3]);
c1=atof(argv[4]);
r0=atof(argv[5]);
al=atof(argv[6]);
tf=atof(argv[7]);
mx=atof(argv[8]);
dx=atof(argv[9]);
 
if  (
     v<V_MIN || v>V_MAX ||
     g<WT_MIN || g>WT_MAX ||
     h<H_MIN || h>H_MAX ||
     r0<R0_MIN || r0>R0_MAX ||
     c1<BC_MIN || c1>BC_MAX ||
     al<AL_MIN || al>AL_MAX ||
     tf<TF_MIN || tf>TF_MAX ||
     x<X_MIN || x>X_MAX ||
     dx<DX_MIN || dx>DX_MAX
     )
      return RANGE_ERR_INPUT; 
      else 
      {h=-h;return ballistics();}
}