#include <stdlib.h>
#include <conio.h>
#include <stdio.h>
#include <malloc.h>
#include <math.h>
#include <assert.h>
#include <string.h>

/*******************************************************************
 *  Model de W.J. Rheingans de vol de les fletxes en l'atmosfera   *
 *                                                                 *
 *  Implementat per Josep Gregori, Desembre 2002                   *
 *  Fonts :                                                        *
 *   W.J. Rheingans, "Exterior Ballistics of Bows and Arrows",     *
 *       http://home.att.net/sajackson/Archery_Physics.html        *
 *   Smallball.c, (C)1996 Ole-Hj.Kristensen (ohk@edeber.nta.no)    *
 *                                                                 *
 *******************************************************************/

// Constants
static const double g=9.81;
static const double PI=3.141592653589793;
// Conversio d'unitats
static const double inch_m=39.3701;
static const double m_inch=0.0254;
static const double yard_m=1.09361;
static const double m_yard=0.9144;
static const double feet_m=3.28084;
static const double m_feet=0.3048;
static const double kg_pound=0.454;
static const double pound_kg=2.2026431;
static const double grains_kg=15418.501;
static const double deg_rad=57.2957795130;
static const double rad_deg=0.017453292519943;

// Semafors de proces
int pause  = 1;
int print  = 0;
int file   = 0;
int zeroit = 0;

// Output stream.
FILE *outstream = stdout;
int nlns = 0, lpp = 60;

// Nombre de passos en la integracio de la trajectoria
const int STEPSperSECOND = 200;
int NSTEP;

// Variables que defineixen un punt de la trajectoria
class vars
{ public :
  float t;
  float vct[4];

  vars() { t=x()=vx()=y()=vy()=0.; }
  vars(const float *z) { init(z); }
  vars(const float tm,const float *z) { init(tm,z); }
  vars(const vars &z) { init(z.t,z.vct); }

  vars &operator=(const vars &z) { init(z.t,z.vct); return *this; }

  float &x()  { return vct[1]; }
  float &vx() { return vct[0]; }
  float &y()  { return vct[3]; }
  float &vy() { return vct[2]; }

  void init(const float *z)  { t=0.; set_vals(z); }
  void init(const float tm,const float *z)  { t=tm; set_vals(z); }

  void get_vals(float *z);
  void set_vals(const float *z);

  float V()    { return sqrt(vx()*vx()+vy()*vy()); }
  float alfa() { return atan(vy()/vx())*deg_rad; }
};

void vars::set_vals(const float *z)
{ for(int i=0;i<4;i++) vct[i]=z[i]; }

void vars::get_vals(float *z)
{ for(int i=0;i<4;i++) z[i]=vct[i]; }


/*****************************
 *   LES CLASSES PROJECTIL   *
 *****************************/

// Projectile generic
class projectile
{ public :
   float point;   // factor
   float mass;    // Kg
   vars xy;       // m, m/s
   vars dxdy;     // m/s m/s2

   projectile(float p,float weight);

   void start(float V,float alfa,float x,float y);
   void start(float V,float alfa) { start(V,alfa,0.,0.); }
   virtual void step(float dh,vars &nxtp);
   virtual void step(float dh,const vars &dd,vars &nxtp);
   virtual void DragCnt()=0;
   virtual void derivate()=0;
   virtual void derivate(vars &dd)=0;
   virtual void derivate(const vars &z,vars &dd)=0;
};

projectile::projectile(float p,float weight)
{ point=p;
  mass=weight/grains_kg/g;
}

void projectile::start(float V,float alfa,float x,float y)
{ xy.t=0.;
  xy.x()=x;
  xy.y()=y;
  xy.vx()=V*cos(alfa*rad_deg)*m_feet;
  xy.vy()=V*sin(alfa*rad_deg)*m_feet;
  derivate();
}

void projectile::step(float dh,vars &nxtp)
{ for(int i=0;i<4;i++)
     nxtp.vct[i] = xy.vct[i] + dh * dxdy.vct[i];
}

void projectile::step(float dh,const vars &dd,vars &nxtp)
{ for(int i=0;i<4;i++)
     nxtp.vct[i] = xy.vct[i] + dh * dd.vct[i];
}


// Classe fletxa amb tot el que determina el seu comportament balistic
class arrow : public projectile
{ private :
   float feather;  // factor
   float diam;     // mm
   float length;   // cm
   float flength;  // mm
   float fheight;  // mm
   float Rf;

   void DragCnt();

  public :

   arrow(float p,float weight,float f,float dm,float ln,float fln,float fhg);

   void derivate();
   void derivate(vars &dd);
   void derivate(const vars &z,vars &dd);

   float Rfv() { return Rf; }
};

arrow::arrow(float p,float weight,float f,float dm,float ln,float fln,
             float fhg) : projectile(p,weight)
{ feather=f; diam=dm; length=ln; flength=fln; fheight=fhg;
  DragCnt();
}

void arrow::DragCnt()
{ Rf = 7574.58 * ( 0.0000013*point*diam*diam*1.E-6 +    // mm*mm
                   0.000000035*length*diam*1.E-5 +      // cm*mm
                   0.000000077*6*0.5*feather*flength*fheight*1.E-6 ) / mass;
}

void arrow::derivate()
{ float   V = xy.V();
  dxdy.t    = xy.t;
  dxdy.x()  = xy.vx();
  dxdy.vx() = -Rf*xy.vx()*V;
  dxdy.y()  = xy.vy();
  dxdy.vy() = -Rf*xy.vy()*V-g;
}

void arrow::derivate(vars &dd)
{ float   V = xy.V();
  dd.t    = xy.t;
  dd.x()  = xy.vx();
  dd.vx() = -Rf*xy.vx()*V;
  dd.y()  = xy.vy();
  dd.vy() = -Rf*xy.vy()*V-g;
}

void arrow::derivate(const vars &z,vars &dd)
{ float V = z.V();
  dd.t    = z.t;
  dd.x()  = z.vx();
  dd.vx() = -Rf*z.vx()*V;
  dd.y()  = z.vy();
  dd.vy() = -Rf*z.vy()*V-g;
}


/************************************************
 *	 Runge-Kutta de quart ordre per calcular el  *
 *	 seguent punt de la trajectoria.             *
 ************************************************/
void RungeKutta(projectile &p,float h)
{ int   i;
  float xmid, halfstep, h6;
  vars  dy_mid,dy_trial,y_trial;

  halfstep = h * 0.5;
  h6 = h / 6.0;

  // t+1/2*h
  xmid = p.xy.t + halfstep;

  // k1 es p.xy.dxdy
  // Calcula y+1/2*h*k1
  p.step(halfstep,y_trial);
  // t+1/2*h
  y_trial.t = xmid;
  // Calcula k2 = f(t+1/2*h,y+1/2*h*k1)
  p.derivate(y_trial,dy_trial);
  // Calcula y+1/2*h*k2
  p.step(halfstep,dy_trial,y_trial);
  // Calcula k3 = f(t+1/2*h,y+1/2*h*k2)
  p.derivate(y_trial,dy_mid);

  // Calcula y+h*k3
  p.step(h,dy_mid,y_trial);
  // Per reutilitzar l'espai de dy_tryal.
  for(i=0;i<4;i++) dy_mid.vct[i] += dy_trial.vct[i];
  // t+h
  y_trial.t = p.xy.t+h;
  // Calcula k4 = f(t+h,y+h*k3)
  p.derivate(y_trial,dy_trial);

  // Calcula y(t+h). Al retorn p.xy conte el seguent punt.
  for(i=0;i<4;i++)
     p.xy.vct[i] += h6*(p.dxdy.vct[i]+dy_trial.vct[i]+2.*dy_mid.vct[i]);
  p.xy.t += h;
}


/**********************************************
 *  CLASSE TRAJECTORIA COM A VECTOR DE PUNTS  *
 **********************************************/

class trace
{ public :
   int n;
   vars *xy;
   projectile &p;
   float mxt;

   trace(int nstep,projectile &pr,float maxt);
   ~trace();

   void calc_traj();
};

trace::trace(int nstep,projectile &pr,float maxt) : p(pr)
{
  xy = new vars[nstep+1];
  n = nstep;
  mxt = maxt;
}

trace::~trace()
{
  if(xy) delete(xy);
  xy=0;
}

void trace::calc_traj()
{ float h;

  if( !xy ) return;
  // Primer punt
  xy[0] = p.xy;
  // Mida de pas
  h = mxt/n;
  // Calcul de la trajectoria en n punts
  for(int k=0;k<n;k++)
  { // Calcula seguent punt
    RungeKutta(p,h);
    // el guarda en la trajectoria
    xy[k+1] = p.xy;
    // el deriva per calcular el seguent.
    p.derivate();
  }
}


/***********************************
 *  FUNCIONS d'INTERPOLACIO LINEAL  *
 ***********************************/

float linint(float xo,float yo,float x1,float y1,float x)
{
  return yo + (x-xo)*(y1-yo)/(x1-xo);
}

void linint(int j,const vars xy[],float t,vars &vout)
{ for(int i=0;i<4;i++)
     vout.vct[i] = linint(xy[j-1].t,xy[j-1].vct[i],xy[j].t,xy[j].vct[i],t);
  vout.t = t;
}


/*********************************
 *   ENTRADA i SORTIDA de DADES  *
 *********************************/

/* ------------------------------------------------------------- */
/* Params for the simulation */
struct calcparam
{ // flight pars
  float     v0;         /* Initial velocity : fps */
  float     grounde;    /* Ground elevation */
  float     maxtime;    /* Run simulation for maxtime, seconds */
  float     maxrange;   /* Display results up to maxrange, m */
  float     angle;      /* Shot elevation on ground, degres */
  float     zero;       /* Where to zero the sights, meters */
  // arrow description
  float     point;      // Head form factor. 1 for ogival,2.3 paralel,7 blunt.
  float     weight;     // Arrow weight : grains
  float     diam;       // Arrow diameter : mm
  float     length;     // Arrow length : cm
  float     feather;    // Feather surface factor. 1 plastic, 2 feathers.
  float     flength;    // Feather length : mm.
  float     fheight;    // Feather height : mm.
};

/* ------------------------------------------------------------------ */
/* Display parameters */
void printpar(calcparam &p,arrow &a)
{ int k;

  fprintf(outstream,"\n\n");
  if( print ) { fprintf(outstream,"\n\n"); nlns += 2; }
  fprintf(outstream,"W.J.Rheingans model, integration by Runge-Kutta on %d steps per"
                    " second\n\n",STEPSperSECOND);
  fprintf(outstream,"Arrow :\n");
  fprintf(outstream,"   Head form factor       : %.2f \n",p.point);
  fprintf(outstream,"   Weight                 : %.2f grains \n",p.weight);
  fprintf(outstream,"   Diameter               : %.2f mm \n",p.diam);
  fprintf(outstream,"   Length                 : %.2f cm \n",p.length);
  fprintf(outstream,"   Feather surface factor : %.2f \n",p.feather);
  fprintf(outstream,"   Feather length         : %.2f mm \n",p.flength);
  fprintf(outstream,"   Feather height         : %.2f mm \n",p.fheight);
  fprintf(outstream,"   Retardation coef.      : %.5E \n",a.Rfv());
  fprintf(outstream,"\nFligth parameters : \n");
  fprintf(outstream,"   Initial velocity       : %.2f fps \n",p.v0);
  fprintf(outstream,"   Ground elevation       : %.3f deg \n",p.grounde);
  fprintf(outstream,"   Maximum range          : %.2f m\n",p.maxrange);
  fprintf(outstream,"   Zero range             : %.2f m \n", p.zero);
  fprintf(outstream,"   Elevation              : %.3f deg \n",p.angle);
  fprintf(outstream,"   Max time               : %.3f sec \n", p.maxtime);

  nlns += 21;
}

/* ------------------------------------------------------------------ */
/* Enter parameters */
void getpar(calcparam &calc)
{ float c,v2,dx;

  printf("\nArrow parameters :");
  printf("\n Head form factor :");
  scanf("%f",&calc.point);
  printf("\n Weight (grains) :");
  scanf("%f",&calc.weight);
  printf("\n Diameter (mm) :");
  scanf("%f",&calc.diam);
  printf("\n Length (cm) :");
  scanf("%f",&calc.length);
  printf("\n Feather surface factor :");
  scanf("%f",&calc.feather);
  printf("\n Feather lengh (mm) :");
  scanf("%f",&calc.flength);
  printf("\n Feather height (mm) :");
  scanf("%f",&calc.fheight);

  printf("\nFlight parameters :");
  printf("\n Enter initial velocity (fps): ");
  scanf("%f",&calc.v0);
  printf("\n Enter ground elevation angle (deg): ");
  scanf("%f",&calc.grounde);
  printf("\n Enter max range (m): ");
  scanf("%f",&calc.maxrange);
  printf("\n Enter zero range (m): ");
  scanf("%f",&calc.zero);
  if( zeroit )
    { v2 = calc.v0*m_feet;  v2 *= v2;
      dx = calc.zero * cos(calc.grounde*rad_deg);
      calc.angle = 0.5*asin(g*calc.zero/v2)*deg_rad;
    }
  else
    { printf("\n Enter elevation angle (deg): ");
      scanf("%f",&calc.angle);
    }
  printf("\n Enter max time (sec): ");
  scanf("%f",&calc.maxtime);
  if(calc.maxrange < 0.000001 && calc.maxtime > 0.000001)
    calc.maxrange = calc.v0 * calc.maxtime;
  else if(calc.maxrange > 0.000001 && calc.maxtime < 0.000001)
    calc.maxtime = calc.maxrange / calc.v0;
}


/* --------------------------------------------------------------- */
/* Display result */

void printpoint(calcparam &p,vars &vout)
{ float ycor = vout.y();
  if(fabs(p.grounde)>0.0001) ycor -= vout.x()*tan(p.grounde/180.*PI);
  fprintf(outstream, "%8.3f %8.1f %8.1f %8.2f %8.1f %8.3f %8.1f\n",
          vout.t, vout.V()*feet_m, vout.vx()*feet_m, vout.x(),
          vout.vy()*feet_m, ycor, vout.alfa());
}

void printtable(calcparam &p,trace &tj)
{ int   i, j;
  float t,gry,ycor,dx,D;
  vars vout;

  // Titols
  fprintf(outstream,"\n %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s\n",
          "  D","time ", "Vel  ", " Vx   ", " X   ", " Vy   ", " Y   ", " Alfa ");
  fprintf(outstream," %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s\n",
          " meter","sec  ", "fps  ", " fps  ", " meter", " fps  ", " meter", "degres");
  // Primer punt t=0
  fprintf(outstream,"%8.2f %8.3f %8.2f %8.2f %8.2f %8.2f %8.3f %8.2f\n",
          0.,tj.xy[0].t, tj.xy[0].V()*feet_m, tj.xy[0].vx()*feet_m, tj.xy[0].x(),
          tj.xy[0].vy()*feet_m, tj.xy[0].y(), tj.xy[0].alfa());
  nlns += 4;

  // Trajectoria a un punt cada metre.
  for(i = 1, j = 0; i <= p.maxrange; i += 1)
  { dx = i * cos(p.grounde*rad_deg);
    while (j < NSTEP && tj.xy[j].x() < dx) j++;
    // Interpolar pel punt precis
    t = linint(tj.xy[j-1].x(),tj.xy[j-1].t,tj.xy[j].x(),tj.xy[j].t,(float)dx);
    linint(j,tj.xy,t,vout);
    ycor = vout.y();
    D = vout.x();
    if(fabs(p.grounde)>0.0001)
      { gry = vout.x()*tan(p.grounde*rad_deg);
        ycor -= gry;
        D = sqrt(vout.x()*vout.x()+gry*gry);
      }
    // Formatejar i treure dades
    fprintf(outstream,"%8.2f %8.3f %8.2f %8.2f %8.2f %8.2f %8.3f %8.2f\n",
            D,vout.t, vout.V()*feet_m, vout.vx()*feet_m, vout.x(),
            vout.vy()*feet_m, ycor, vout.alfa());
    // Control de pagina
    nlns++;
    // A cada pantalla parar i tornar a treure titols, si cal.
    if( ! print && nlns >= 22 && pause)
      { printf("\n Press Enter to continue ");
        getch();
        printf("\n");
        fprintf(outstream,"\n %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s\n",
                "  D","time ", "Vel  ", " Vx   ", " X   ", " Vy   ", " Y   ", " Alfa ");
        fprintf(outstream," %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s\n",
                " meter","sec  ", "fps  ", " fps  ", " meter", " fps  ", " meter", "degres");
        nlns = 3;
      }
    else if( print && nlns >= lpp )
      { fprintf(outstream,"\n\n\n");
        fprintf(outstream,"\n %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s\n",
                "  D","time ", "Vel  ", " Vx   ", " X   ", " Vy   ", " Y   ", " Alfa ");
        fprintf(outstream," %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s %8.8s\n",
                " meter","sec  ", "fps  ", " fps  ", " meter", " fps  ", " meter", "degres");
        nlns = 6;
      }
  }
  fprintf(outstream,"#%c", 'L' - '@'); // Final de sortida.
}

/* ------------------------------------------------------------------ */
/* Crude calculation of elevation, given a desired zero range */
void zeroin(calcparam &p, trace &tj)
{ int    i;
  float  t,Dx,Dy,e;
  vars   vout;

  Dx=p.zero*cos(p.grounde*rad_deg);
  for (i = 1; i < NSTEP && tj.xy[i].x() < Dx; i++);
  // Interpolar pel punt precis
  t = linint(tj.xy[i-1].x(),tj.xy[i-1].t,tj.xy[i].x(),tj.xy[i].t,Dx);
  linint(i,tj.xy,t,vout);
  // Corregir l'angle d'elevacio
  Dy=p.zero*sin(p.grounde*rad_deg);
  e = (vout.y()-Dy)*cos(p.grounde*rad_deg);
  p.angle += atan(-e/p.zero)*deg_rad;
}


/************************************************************************
 *  Initialize the system and compute tables. The integration of the    *
 *  differential equations is done with a 4th order Runge-Kutta method, *
 *  because it is robust and fairly accurate.                           *
 *  The adventurous may use a Burlisch-Stoer or predictor-corrector     *
 *  method if they wish, it should be able to give the same accuracy    *
 *  with less computer time.                                            *
 ************************************************************************/

/* -------------------------------------------------------------- */
void compute(calcparam &p)
{ int k;

  NSTEP = p.maxtime * STEPSperSECOND;    /* Number of steps in calculation */

  arrow a(p.point,p.weight,p.feather,p.diam,p.length,p.flength,p.fheight);

  trace tj(NSTEP,a,p.maxtime);
  if( ! tj.xy ) return;

  a.start(p.v0,p.angle+p.grounde,0.,0.);
  tj.calc_traj();

  if( ! zeroit )
  { printpar(p,a);
    printtable(p,tj);
  }

  /* Zero it twice, as the approximation is crude, but converges */
  if( zeroit )
   { zeroin(p,tj);
     a.start(p.v0,p.angle+p.grounde,0.,0.);
     tj.calc_traj();

     zeroin(p,tj);
     a.start(p.v0,p.angle+p.grounde,0.,0.);
     tj.calc_traj();

     printpar(p,a);
     printtable(p,tj);
   }
}

/* ------------------------------------------------------------- */
main(int argc, char **argv)
{ float c;
  unsigned char  cs;

  for(int i=1; i<argc && argv[i]; i++)
  { switch (cs = *argv[i] == '-' ? *(argv[i]+1) : *argv[i])
     { case 'h': case 'H':
        { printf(" -p option turns on printer\n");
          printf(" -f filename send output to file\n");
          printf(" -z zeroes the shot at the given range\n");
          return 1;
        }
       case 'p': case 'P':
        { pause = 0;
          print = 1;
          outstream = stdprn;
          break; 
        }
       case 'f': case 'F':
        { pause = 0;
          file = 1;
          outstream = fopen(argv[i]+2,"w");
          if( ! outstream ) outstream = stdout;
          break;
        }
       case 'z': case 'Z': zeroit = 1; break;
       default : ;
     }
  }

  // Obtenir parametres.
  calcparam calc;
  getpar(calc);
  // Calcular trajectoria i lelevacio correcte
  compute(calc);

  if( file ) fclose(outstream);
  exit(0);
}