/*

Automag.c version 1.0  08/17/94

This program calculates the behavior of an automag paintball gun using
a simple model of the gun. See file "automag.txt" for an explanation
of the program.

I have no proprietary feelings about this code. If you can improve it,
please feel free to do so, offer it to anyone you want, I don't need to
be notified or credited, or anything.
                                                - Paul Reiser -
*/


#define VAX                     /* Compiler Specification               */

#ifdef LSC                      /* For use with Lightspeed C compiler   */
#include "math.h"
#include "unix.h"
#include "stdio.h"
#include "ctype.h"
#include "storage.h"
#define MAC
#endif

#ifdef THC                      /* For use with Think C compiler        */
#include "math.h"
#include "unix.h"
#include "stdio.h"
#include "ctype.h"
#include "Memory.h"
#define MAC
#endif

#ifdef VAX                      /* For use with VAX C compiler  */
#include <math.h>
#include <stdio.h>
#endif

/*      General Globals and defines     */

#define pi      3.1415926535

typedef struct  /* Thermodynamic variables      */

{ double Ps;            /* Gas pressure of source               (psi)   */
  double Pr;            /* Gas pressure after reg. valve        (psi)   */
  double Pa;            /* Gas pressure in air chamber          (psi)   */
  double Pb;            /* Gas pressure in barrel               (psi)   */
  double Px;            /* Atmospheric pressure                 (psi)   */

  double Vr;            /* regulator volume                     (in^2)  */
  double Va;            /* air chamber volume                   (in^2)  */
  
  double T;             /* Temperature                          (deg F) */
  double m;             /* Atomic weight of gas                         */
  
  double GC;            /* Gas consumption                      (oz)    */
}VARpv;

typedef struct  /* Regulator valve variables                            */

{ double S;     /* Surface area of regulator piston             (in^2)  */
  double x0;    /* Nominal R-valve to top of R-seal             (in)    */
  double Cl;    /* Leakage flow conductance                             */
  double Z;     /* Flow conductance/length                              */
  double T;     /* Position of regulator nut                    (turns) */
  double K;     /* spring constant of the reg piston spring     (lb/in) */
  double H;     /* Regulator nut turns per inch                         */
  double x;     /* Position of regulator valve                  (in)    */
}VARr;

typedef struct  /* on-off valve variables                               */

{ double C;     /* Leakage conductance                                  */
  double Z;     /* Flow conductance/Length                              */
  double x;     /* On-off pin position                          (in)    */
}VARo;

typedef struct  /* Power tube variables                                 */

{ double M;     /* Mass of bolt                                 (oz)    */
  double S;     /* Area of bolt                                 (in^2)  */
  double xc;    /* Position of bolt when it latches             (in)    */
  double Cl;    /* Leakage conductance                                  */
  double Z;     /* flow conductance/length                              */
  double K;     /* Spring constant of power tube spring         (lb/in) */
  double R;     /* Friction coefficient                      (lb-sec/in)*/
  double x0;    /* Relaxed length of power tube spring          (in)    */
  double xi;    /* Length of p.t. spring when bolt is latched   (in)    */
  double xm;    /* Maximum position of bolt                     (in)    */
  double x;     /* Position of bolt                             (in)    */
  double v;     /* Velocity of bolt                             (fps)   */
}VARp;

typedef struct  /* barrel & ball variables                              */

{ double M;     /* Mass of ball                                         */
  double D;     /* Diameter of barrel and ball                          */
  double L;     /* Length of barrel                             (in)    */
  double x;     /* Position of ball                             (in)    */
  double v;     /* Velocity of ball                             (fps)   */
}VARb;

typedef struct  /* trigger variables                                    */

{ double B;     /* Sear angle                                   (deg)   */
  double Bo;    /* Sear angle at which on-off valve turns on    (deg)   */
  double Bb;    /* Sear angle at which bolt is released         (deg)   */
  double xs;    /* Sear pivot to on-off pin distance            (in)    */
}VARt;

/*      Functions       */

/*int   tstep();*/

double Qsr();           /* Rate gas enters the regulator        (m/s)   */
double Qra();           /* Rate gas enters the air chamber      (m/s)   */
double Qab();           /* Rate gas enters the barrel           (m/s)   */
double Qbx();           /* Rate gas leaves the barrel           (m/s)   */

int UpDate_r();         /* Update regulator valve                       */
int UpDate_o();
int UpDate_p();
int UpDate_b();
int UpDate_t();

/*----------------------------------------------------------------------*/

main()

{ int   i;

  VARr  Zr;
  VARo  Zo;
  VARp  Zp;
  VARb  Zb;
  VARt  Zt;
  VARpv PV;

  double t;     /* time                                 (sec)   */
  double dt;    /* time increment                               */
  int n;        /* number of iterations                         */

/* Thermodynamic variables                                      */

  PV.Px = 14.6959;      /* atmospheric pressure         (psi)   */
  PV.Ps = 1000.0;       /* source pressure              (psi)   */
  PV.Pr =  449.0;       /* regulator pressure           (psi)   */
  PV.Pa =  449.0;       /* air chamber pressure         (psi)   */
  PV.Pb =  PV.Px;       /* barrel pressure              (psi)   */

  PV.Vr = 0.05;         /* Volume of regulator chamber  (in^3)  */
  PV.Va = 0.463;                /* Volume of air chamber        (in^3)  */
  
  PV.T  = 70.0;         /* Temperature                  (deg F) */
  PV.m  = 44.0;         /* Atomic weight of gas                 */
  
  PV.GC = 0.0;          /* Gas consumption              (oz)    */

/* Regulator valve variables                                    */

  Zr.S  = 0.19635;      /* Surface area of regulator piston     (in^2)  */
  Zr.x0 = 0.0625;       /* Nominal R-valve to top of R-seal     (in)    */
  Zr.Cl = 0.0;          /* Leakage flow conductance                     */
  Zr.Z  = 200.0;        /* Flow conductance/length                      */
  Zr.T  = 0.0;          /* Position of regulator nut            (turns) */
  Zr.K  = 1367.5;       /* spring const. of reg piston spring   (lb/in) */
  Zr.H  = 32.0;         /* Regulator nut turns per inch                 */
  Zr.x  = Zr.x0;        /* Position of regulator valve          (in)    */ 

/* on-off valve variables                                               */

  Zo.C  = 0.0;          /* Leakage conductance                          */
  Zo.Z  = 200.0;        /* Flow conductance/Length                      */
  Zo.x  = 0.0;          /* on-off pin position                          */

/* Power tube variables                                                 */

  Zp.M  = 1.848;        /* Mass of bolt                         (oz)    */
  Zp.S  = 0.181;        /* Area of bolt                         (in^2)  */
  Zp.xc = -0.1;         /* Position of bolt when it latches     (in)    */
  Zp.Cl = 0.0;          /* Leakage conductance                          */
  Zp.Z  = 360.0;        /* flow conductance/length                      */
  Zp.K  = 11.3;         /* Spring constant of power tube spring (lb/in) */
  Zp.R  = 1.0;          /* Friction coefficient              (lb-sec/in)*/
  Zp.x0 = 1.89;         /* Relaxed length of power tube spring  (in)    */
  Zp.xi = 1.65;         /* Length of pts. when bolt is latched  (in)    */
  Zp.xm = 1.18;         /* Maximum position of bolt             (in)    */
  Zp.x  = Zp.xc;        /* Position of bolt                     (in)    */
  Zp.v  = 0.0;          /* Velocity of bolt                     (fps)   */

/* barrel & ball variables                                              */

  Zb.M  = 0.113;        /* Mass of ball                         (oz)    */
  Zb.D  = 0.68;         /* Diameter of barrel and ball          (in)    */
  Zb.L  = 11.0;         /* Length of barrel                     (in)    */
  Zb.x  = 0.00;         /* Position of ball                     (in)    */
  Zb.v  = 0.00;         /* Velocity of ball                     (fps)   */

/* trigger variables                                                    */

  Zt.B  = 0.0;  /* Sear angle                                   (deg)   */
  Zt.Bo = 2.0;  /* Sear angle at which on-off valve turns on    (deg)   */
  Zt.Bb = 4.0;  /* Sear angle at which bolt is released (deg)   */
  Zt.xs = 0.8;  /* Sear pivot to on-off pin distance            (in)    */

  t     = 0.0;  /* beginning time                               (sec)   */
  dt    = 0.0002;/* time step                                   (sec)   */
  n     = 75;   /* number of steps                                      */
  
/* Step through time and print out some variables                       */

  for(i=0;i<n;i++)
  { 
    printf("\nt=%6.4lf Pr=%6.2lf Pa=%6.2lf ",t,PV.Pr,PV.Pa);
    printf("GC=%e Zb.v=%6.4lf Zp.x=%6.4lf Zt.B=%3.1lf",PV.GC,Zb.v,Zp.x,Zt.B);
    tstep(&t,&dt,&PV,&Zr,&Zo,&Zp,&Zb,&Zt);
  }

  printf("\n\nDONE");

}

/*----------------------------------------------------------------------*/

int tstep(pt,pdt,pPV,pZr,pZo,pZp,pZb,pZt)
#define t  (*pt)
#define dt (*pdt)
#define PV (*pPV)
#define Zr (*pZr)
#define Zo (*pZo)
#define Zp (*pZp)
#define Zb (*pZb)
#define Zt (*pZt)

  double t;
  double dt;
  VARpv PV;
  VARr  Zr;
  VARo  Zo;
  VARp  Zp;
  VARb  Zb;
  VARt  Zt;

/* step forward in time dt seconds              */

{ int i;
  double hQsr,hQra,hQab,hQbx;
  double Vb;
  double Vbo;
  double TK;

  TK = 273.16 + (PV.T-32.0)/1.8;        /* Temperature in kelvin        */
  
/* Update source pressure       */

  { double c[4];
    static int ok;
    c[0] = -47.922400;
    c[1] =  0.51399206;
    c[2] = -1.67373965e-03;
    c[3] =  1.89693767e-06;
    PV.Ps = 0.0;
    for(i=0;i<=3;i++)PV.Ps = PV.Ps + c[i]*pow(TK,i);
    PV.Ps = exp(PV.Ps);
    if(ok == 0)
    { printf("\nPV.Ps = %6.2lf",PV.Ps);
      ok = 1;
    }
  }

/* Update moving parts  */

  UpDate_t(&Zt,&Zb,&Zp,&PV,t,dt);/* trigger variables           */
  UpDate_o(&PV,&Zo,&Zt,t,dt);   /* on-off valve                 */
  UpDate_p(&PV,&Zp,&Zt,t,dt);   /* power tube                   */
  UpDate_r(&PV,&Zr,t,dt);       /* regulator valve and piston   */

  Vbo = 0.25*pi*Zb.D*Zb.D*(Zb.x+0.1);   /* save old barrel volume       */
  UpDate_b(&PV,&Zb,&Zp,t,dt);           /* ball & barrel                */

/*  Update Pressures    */

  hQsr = Qsr(PV.Ps,PV.Pr,&Zr);
  hQra = Qra(PV.Pr,PV.Pa,&Zo,&Zt);
  hQab = Qab(PV.Pa,PV.Pb,&Zp);
  hQbx = Qbx(PV.Pb,PV.Px,&Zb);

  PV.Pr = PV.Pr + (hQsr - hQra)*dt/PV.Vr;
  
  PV.Pa = PV.Pa + (hQra - hQab)*dt/PV.Va;
  
  if(Zb.x < Zb.L)
  { Vb = 0.25*pi*Zb.D*Zb.D*(Zb.x+0.1);  /* Vb = new barrel volume       */
    PV.Pb = (PV.Pb*Vbo + (hQab - hQbx)*dt)/Vb;
  }
  else PV.Pb = PV.Px;
  
/* Update gas consumption       */

  PV.GC = PV.GC + 4.85e-4*hQsr*dt*PV.m/TK;

/* Update time                  */

  t = t+dt;

  return;

}

#undef t
#undef dt
#undef PV
#undef Zr
#undef Zo
#undef Zp
#undef Zb
#undef Zt

/*----------------------------------------------------------------------*/

double Qsr(Ps,Pr,pZr)
#define Zr (*pZr)
  double Ps;    /* source pressure                      (psi)   */
  double Pr;    /* regulator pressure                   (psi)   */
  VARr   Zr;    /* regulator variables                          */

/* calculate the rate of flow from source to regulator          */

/*
{ double S;     / Surface area of regulator piston              (in^2)  /
  double x0;    / Nominal R-valve to top of R-seal              (in)    /
  double Cl;    / Leakage flow conductance                              /
  double Z;     / Flow conductance/length                               /
  double T;     / Position of regulator nut                     (turns) /
  double K;     / spring constant of the reg piston spring      (lb/in) /
  double H;     / Regulator nut turns per inch                          /
  double x;     / Position of regulator valve                   (in)    /
}VARr;
*/

{ double Q,L;

  if(Zr.x < 0.0)L=0; else L=Zr.x;

  Q  = (Zr.Cl+Zr.Z*L)*(Ps-Pr);

  return(Q);
}

#undef Zr

/*----------------------------------------------------------------------*/

double Qra(Pr,Pa,pZo,pZt)
#define Zo (*pZo)
#define Zt (*pZt)
  double Pr;    /* regulator pressure                           (psi)   */
  double Pa;    /* air chamber pressure                         (psi)   */
  VARo   Zo;    /* on-off valve variables                               */
  VARt   Zt;    /* trigger variables                                    */

/* calculate the rate of flow from regulator to air chamber             */

/*
{ double C;     / Leakage conductance                                   /
  double Z;     / Flow conductance/Length                               /
  double x;     / On-off pin position                           (in)    /
}VARo;

{ double B;     / Sear angle                                    (deg)   /
  double Bo;    / Sear angle at which on-off valve turns on     (deg)   /
  double Bb;    / Sear angle at which bolt is released          (deg)   /
  double xs;    / Sear pivot to on-off pin distance             (in)    /
}VARt;
*/

{ double Q,L;

  if(Zo.x < 0.0)L=0.0; else L=Zo.x;
  
  Q = (Zo.C + Zo.Z*L)*(Pr-Pa);

  return(Q);
}

#undef Zo
#undef Zt

/*----------------------------------------------------------------------*/

double Qab(Pa,Pb,pZp)
#define Zp (*pZp)
  double Pa;    /* air chamber pressure                 (psi)   */
  double Pb;    /* barrel pressure                      (psi)   */
  VARp   Zp;    /* power tube exit variables                    */
  
/* calculate the rate of flow from air chamber to barrel        */

/*
{ double M;     / Mass of bolt                                  (oz)    /
  double S;     / Area of bolt                                  (in^2)  /
  double xc;    / Position of bolt when it latches              (in)    /
  double Cl;    / Leakage conductance                                   /
  double Z;     / flow conductance/length                               /
  double K;     / Spring constant of power tube spring          (lb/in) /
  double R;     / Friction coefficient                       (lb-sec/in)/
  double x0;    / Relaxed length of power tube spring           (in)    /
  double xi;    / Length of p.t. spring when bolt is latched    (in)    /
  double xm     / Maximum position of bolt                      (in)    /
  double x;     / Position of bolt                              (in)    /
  double v;     / Velocity of bolt                              (fps)   /
}VARp;

*/

{ double Q,L;

  if(Zp.x < 0.0)L=0; else L=Zp.x;
  
  Q = (Zp.Cl+Zp.Z*L)*(Pa-Pb);

  return(Q);

}

#undef Zp

/*----------------------------------------------------------------------*/

double Qbx(Pb,Px,pZb)
#define Zb (*pZb)
  double Pb;    /* barrel pressure                      (psi)   */
  double Px;    /* Atmospheric pressure                 (psi)   */
  VARb   Zb;    /* barrel vent variables                        */

/* calculate the rate of flow from barrel to atmosphere */

/*
{ double M;     / Mass of ball                                          /
  double D;     / Diameter of barrel and ball                           /
  double L;     / Length of barrel                              (in)    /
  double x;     / Position of ball                              (in)    /
  double v;     / Velocity of ball                              (fps)   /
}VARb;
*/

{ double Q;

  Q = 0.0;
  return(Q);

}

#undef Zb

/*----------------------------------------------------------------------*/

UpDate_r(pPV,pZr,t,dt)
#define   PV (*pPV)
#define   Zr (*pZr)
  VARpv   PV;           /* Pressures and Volumes        */
  VARr    Zr;           /* regulator valve variables    */
  double  t;
  double  dt;

/*
  Update the position of the regulator valve. Its assumed
  that the top of the valve is attached to the regulator piston.
*/

/*
{ double S;     / Surface area of regulator piston              (in^2)  /
  double x0;    / Nominal R-valve to top of R-seal              (in)    /
  double Cl;    / Leakage flow conductance                              /
  double Z;     / Flow conductance/length                               /
  double T;     / Position of regulator nut                     (turns) /
  double K;     / spring constant of the reg piston spring      (lb/in) /
  double H;     / Regulator nut turns per inch                          /
  double x;     / Position of regulator valve                   (in)    /
}VARr;
*/

{ Zr.x = Zr.x0 + Zr.T/Zr.H - (PV.Pr-PV.Px)*Zr.S/Zr.K;
  if(Zr.x > Zr.x0)      Zr.x = Zr.x0;
  if(Zr.x < 0.0)        Zr.x = 0.0;
  return(0);
}

#undef Zr
#undef PV

/*----------------------------------------------------------------------*/

int UpDate_o(pPV,pZo,pZt,t,dt)
#define PV (*pPV)
#define Zo (*pZo)
#define Zt (*pZt)
  VARpv   PV;           /* Pressures and Volumes        */
  VARo    Zo;           /* on-off valve variables       */
  VARt    Zt;           /* trigger variables            */
  double  t;
  double  dt;

/* Update the on-off valve                      */

/*
{ double C;     / Leakage conductance                                   /
  double Z;     / Flow conductance/Length                               /
  double x;     / On-off pin position                           (in)    /
}VARo;

{ double B;     / Sear angle                                    (deg)   /
  double Bo;    / Sear angle at which on-off valve turns on     (deg)   /
  double Bb;    / Sear angle at which bolt is released          (deg)   /
  double xs;    / Sear pivot to on-off pin distance             (in)    /
}VARt;
*/

{ Zo.x = -Zt.xs*(Zt.B-Zt.Bo)*pi/180.;   /* position of on-off pin       */
  return(0);
}

#undef Zo
#undef Zt
#undef PV

/*----------------------------------------------------------------------*/

int UpDate_p(pPV,pZp,pZt,t,dt)
#define PV (*pPV)
#define Zp (*pZp)
#define Zt (*pZt)
  VARpv  PV;            /* Pressures and Volumes        */
  VARp   Zp;            /* power tube variables         */
  VARt   Zt;            /* trigger variables            */
  double t;
  double dt;

/*
  Update the position and velocity of the bolt. If the bolt is inside
  the air chamber, both the air chamber pressure and the power tube spring
  exert a force on the bolt. When its outside the air chamber, only the
  spring exerts a force. If the bolt position ever gets below Zp.xbl and the
  sear angle is smaller than Zp.Bb, the bolt is stopped.
  
{ double M;     / Mass of bolt                                  (oz)    /
  double S;     / Area of bolt                                  (in^2)  /
  double xc;    / Position of bolt when it latches              (in)    /
  double Cl;    / Leakage conductance                                   /
  double Z;     / flow conductance/length                               /
  double K;     / Spring constant of power tube spring          (lb/in) /
  double R;     / Friction coefficient                       (lb-sec/in)/
  double x0;    / Relaxed length of power tube spring           (in)    /
  double xi;    / Length of p.t. spring when bolt is latched    (in)    /
  double xm     / Maximum position of bolt                      (in)    /
  double x;     / Position of bolt                              (in)    /
  double v;     / Velocity of bolt                              (fps)   /
}VARp;


{ double B;     / Sear angle                                    (deg)   /
  double Bo;    / Sear angle at which on-off valve turns on     (deg)   /
  double Bb;    / Sear angle at which bolt is released          (deg)   /
  double xs;    / Sear pivot to on-off pin distance             (in)    /
}VARt;

*/

{ if((Zp.x <= Zp.xc) && (Zt.B < Zt.Bb)) /* if Bolt is latched, stop it  */
  {
    Zp.x = Zp.xc;
    Zp.v = 0.0;
  }
  else
  { double vv,F,A;
  
/*  Velocity (in/sec)            */
    vv = Zp.v*12.0;
    
/*  Force (lb)  */
    F = -Zp.K*(Zp.x+Zp.x0-Zp.xi-Zp.xc) + (PV.Pa-PV.Pb)*Zp.S;

    if(Zp.x <= 0.0)F = F - Zp.R*Zp.v;
    
/*  Acceleration (in/sec^2)     */
    A =  514.784*F/Zp.M;

    Zp.x = Zp.x + (vv + 0.5*A*dt)*dt;
    Zp.v = (vv + A*dt)/12.0;
  }

/* if bolt is at maximum position, stop it      */

  if(Zp.x > Zp.xm)
  { Zp.x = Zp.xm;
    Zp.v = 0.0;
  }

/* if bolt is at minimum position, stop it      */

  if(Zp.x < Zp.xc)
  { Zp.x = Zp.xc;
    Zp.v = 0.0;
  }
  
  return;
}

#undef Zp
#undef Zt
#undef PV

/*----------------------------------------------------------------------*/

int UpDate_b(pPV,pZb,pZp,t,dt)
#define PV (*pPV)
#define Zb (*pZb)
#define Zp (*pZp)
  VARpv  PV;            /* Pressures and Volumes        */
  VARb   Zb;            /* barrel & ball variables      */
  VARp   Zp;            /* power tube variables         */
  double t;
  float  dt;

/* Update the position and velocity of the ball                 */

/*
{ double M;     / Mass of ball                                          /
  double D;     / Diameter of barrel and ball                           /
  double L;     / Length of barrel                              (in)    /
  double x;     / Position of ball                              (in)    /
  double v;     / Velocity of ball                              (fps)   /
}VARb;

{ double M;     / Mass of bolt                                  (oz)    /
  double S;     / Area of bolt                                  (in^2)  /
  double xc;    / Position of bolt when it latches              (in)    /
  double Cl;    / Leakage conductance                                   /
  double Z;     / flow conductance/length                               /
  double K;     / Spring constant of power tube spring          (lb/in) /
  double R;     / Friction coefficient                       (lb-sec/in)/
  double x0;    / Relaxed length of power tube spring           (in)    /
  double xi;    / Length of p.t. spring when bolt is latched    (in)    /
  double xm     / Maximum position of bolt                      (in)    /
  double x;     / Position of bolt                              (in)    /
  double v;     / Velocity of bolt                              (fps)   /
}VARp;

*/

{ if(Zb.x < Zb.L)
  { double A,vv;

/* Velocity (in/sec)    */  
      vv = Zb.v*12.0;

/* Acceleration (in/sec^2)      */
      A = 6177.4*(PV.Pb-PV.Px)*0.25*pi*Zb.D*Zb.D/Zb.M;

      Zb.x = Zb.x + (vv + 0.5*A*dt)*dt;
      Zb.v = (vv + A*dt)/12.0;
  }

  if(Zb.x < 0.0)
  {  Zb.x = 0.0;
     Zb.v = 0.0;
  }

  return(0);
}

#undef Zb
#undef Zp
#undef PV

/*-----------------------------------------------------------------------*/

int UpDate_t(pZt,pZb,pZp,pPV,t,dt)
#define  Zt (*pZt)
#define  Zb (*pZb)
#define  Zp (*pZp)
#define  PV (*pPV)
  VARt   Zt;            /* Trigger variables            */
  VARb   Zb;            /* ball and barrel variables    */
  VARp   Zp;            /* power tube variables         */
  VARpv  PV;            /* Pressures and volumes        */
  double t;
  float  dt;

/* Update the position of the trigger                   */
/* put a ball in the chamber if so desired              */

/*
{ double B;     / Sear angle                                    (deg)   /
  double Bo;    / Sear angle at which on-off valve turns on     (deg)   /
  double Bb;    / Sear angle at which bolt is released          (deg)   /
  double xs;    / Sear pivot to on-off pin distance             (in)    /
}VARt;

{ double M;     / Mass of ball                                          /
  double D;     / Diameter of barrel and ball                           /
  double L;     / Length of barrel                              (in)    /
  double x;     / Position of ball                              (in)    /
  double v;     / Velocity of ball                              (fps)   /
}VARb;

*/

{ 
  if(t >= 0.005)Zt.B=Zt.Bb; else Zt.B=0.0; /* pull trigger at t=0.005 sec*/
  return(0);
}

#undef Zt
#undef Zb
#undef Zp
#undef PV

/*-----------------------------------------------------------------------*/

int Err0r(n)
  int (n);

{
        printf("\nError %d",n);
        exit();
        return(0);
}