qcsrc/common/weapons/calculations.qc

   1 // =============================
   2 //  Explosion Force Calculation
   3 // =============================
   4 float explosion_calcpush_getmultiplier(vector explosion_v, vector target_v)
   5 {
   6         float a;
   7         a  = explosion_v * (explosion_v - target_v);
   8
   9         if(a <= 0)
  10                 // target is too fast to be hittable by this
  11                 return 0;
  12
  13         a /= (explosion_v * explosion_v);
  14                 // we know we can divide by this, or above a would be == 0
  15
  16         return a;
  17 }
  18
  19 #if 0
  20 vector explosion_calcpush(vector explosion_v, float explosion_m, vector target_v, float target_m, float elasticity)
  21 {
  22         // solution of the equations:
  23         //    v'                = v + a vp             // central hit
  24         //    m*v'   + mp*vp'   = m*v + mp*vp          // conservation of momentum
  25         //    m*v'^2 + mp*vp'^2 = m*v^2 + mp*vp^2      // conservation of energy (ELASTIC hit)
  26         // -> a = 0                                    // case 1: did not hit
  27         // -> a = 2*mp*(vp^2 - vp.v) / ((m+mp) * vp^2) // case 2: did hit
  28         //                                             // non-elastic hits are somewhere between these two
  29
  30         // this would be physically correct, but we don't do that
  31         return explosion_v * explosion_calcpush_getmultiplier(explosion_v, target_v) * (
  32                 (1 + elasticity) * (
  33                         explosion_m
  34                 ) / (
  35                         target_m + explosion_m
  36                 )
  37         ); // note: this factor is at least 0, at most 2
  38 }
  39 #endif
  40
  41 // simplified formula, tuned so that if the target has velocity 0, we get exactly the original force
  42 vector damage_explosion_calcpush(vector explosion_f, vector target_v, float speedfactor)
  43 {
  44         // if below 1, the formulas make no sense (and would cause superjumps)
  45         if(speedfactor < 1)
  46                 return explosion_f;
  47
  48 #if 0
  49         float m;
  50         // find m so that
  51         //   speedfactor * (1 + e) * m / (1 + m) == 1
  52         m = 1 / ((1 + 0) * speedfactor - 1);
  53         vector v;
  54         v = explosion_calcpush(explosion_f * speedfactor, m, target_v, 1, 0);
  55         // the factor we then get is:
  56         //   1
  57         print(sprintf("MASS: %f\nv: %v -> %v\nENERGY BEFORE == %f + %f = %f\nENERGY AFTER >= %f\n",
  58                 m,
  59                 target_v, target_v + v,
  60                 target_v * target_v, m * explosion_f * speedfactor * explosion_f * speedfactor, target_v * target_v + m * explosion_f * speedfactor * explosion_f * speedfactor,
  61                 (target_v + v) * (target_v + v)));
  62         return v;
  63 #endif
  64         return explosion_f * explosion_calcpush_getmultiplier(explosion_f * speedfactor, target_v);
  65 }
  66
  67
  68 // =========================
  69 //  Shot Spread Calculation
  70 // =========================
  71
  72 vector cliptoplane(vector v, vector p)
  73 {
  74         return v - (v * p) * p;
  75 }
  76
  77 vector solve_cubic_pq(float p, float q)
  78 {
  79         float D, u, v, a;
  80         D = q*q/4.0 + p*p*p/27.0;
  81         if(D < 0)
  82         {
  83                 // irreducibilis
  84                 a = 1.0/3.0 * acos(-q/2.0 * sqrt(-27.0/(p*p*p)));
  85                 u = sqrt(-4.0/3.0 * p);
  86                 // a in range 0..pi/3
  87                 // cos(a)
  88                 // cos(a + 2pi/3)
  89                 // cos(a + 4pi/3)
  90                 return
  91                         u *
  92                         (
  93                                 '1 0 0' * cos(a + 2.0/3.0*M_PI)
  94                                 +
  95                                 '0 1 0' * cos(a + 4.0/3.0*M_PI)
  96                                 +
  97                                 '0 0 1' * cos(a)
  98                         );
  99         }
 100         else if(D == 0)
 101         {
 102                 // simple
 103                 if(p == 0)
 104                         return '0 0 0';
 105                 u = 3*q/p;
 106                 v = -u/2;
 107                 if(u >= v)
 108                         return '1 1 0' * v + '0 0 1' * u;
 109                 else
 110                         return '0 1 1' * v + '1 0 0' * u;
 111         }
 112         else
 113         {
 114                 // cardano
 115                 u = cbrt(-q/2.0 + sqrt(D));
 116                 v = cbrt(-q/2.0 - sqrt(D));
 117                 return '1 1 1' * (u + v);
 118         }
 119 }
 120 vector solve_cubic_abcd(float a, float b, float c, float d)
 121 {
 122         // y = 3*a*x + b
 123         // x = (y - b) / 3a
 124         float p, q;
 125         vector v;
 126         p = (9*a*c - 3*b*b);
 127         q = (27*a*a*d - 9*a*b*c + 2*b*b*b);
 128         v = solve_cubic_pq(p, q);
 129         v = (v -  b * '1 1 1') * (1.0 / (3.0 * a));
 130         if(a < 0)
 131                 v += '1 0 -1' * (v_z - v_x); // swap x, z
 132         return v;
 133 }
 134
 135 vector findperpendicular(vector v)
 136 {
 137         vector p;
 138         p_x = v_z;
 139         p_y = -v_x;
 140         p_z = v_y;
 141         return normalize(cliptoplane(p, v));
 142 }
 143
 144 vector W_CalculateSpread(vector forward, float spread, float spreadfactor, float spreadstyle)
 145 {
 146         float sigma;
 147         vector v1 = '0 0 0', v2;
 148         float dx, dy, r;
 149         float sstyle;
 150         spread *= spreadfactor; //g_weaponspreadfactor;
 151         if(spread <= 0)
 152                 return forward;
 153         sstyle = spreadstyle; //autocvar_g_projectiles_spread_style;
 154
 155         if(sstyle == 0)
 156         {
 157                 // this is the baseline for the spread value!
 158                 // standard deviation: sqrt(2/5)
 159                 // density function: sqrt(1-r^2)
 160                 return forward + randomvec() * spread;
 161         }
 162         else if(sstyle == 1)
 163         {
 164                 // same thing, basically
 165                 return normalize(forward + cliptoplane(randomvec() * spread, forward));
 166         }
 167         else if(sstyle == 2)
 168         {
 169                 // circle spread... has at sigma=1 a standard deviation of sqrt(1/2)
 170                 sigma = spread * 0.89442719099991587855; // match baseline stddev
 171                 v1 = findperpendicular(forward);
 172                 v2 = cross(forward, v1);
 173                 // random point on unit circle
 174                 dx = random() * 2 * M_PI;
 175                 dy = sin(dx);
 176                 dx = cos(dx);
 177                 // radius in our dist function
 178                 r = random();
 179                 r = sqrt(r);
 180                 return normalize(forward + (v1 * dx + v2 * dy) * r * sigma);
 181         }
 182         else if(sstyle == 3) // gauss 3d
 183         {
 184                 sigma = spread * 0.44721359549996; // match baseline stddev
 185                 // note: 2D gaussian has sqrt(2) times the stddev of 1D, so this factor is right
 186                 v1 = forward;
 187                 v1_x += gsl_ran_gaussian(sigma);
 188                 v1_y += gsl_ran_gaussian(sigma);
 189                 v1_z += gsl_ran_gaussian(sigma);
 190                 return v1;
 191         }
 192         else if(sstyle == 4) // gauss 2d
 193         {
 194                 sigma = spread * 0.44721359549996; // match baseline stddev
 195                 // note: 2D gaussian has sqrt(2) times the stddev of 1D, so this factor is right
 196                 v1_x = gsl_ran_gaussian(sigma);
 197                 v1_y = gsl_ran_gaussian(sigma);
 198                 v1_z = gsl_ran_gaussian(sigma);
 199                 return normalize(forward + cliptoplane(v1, forward));
 200         }
 201         else if(sstyle == 5) // 1-r
 202         {
 203                 sigma = spread * 1.154700538379252; // match baseline stddev
 204                 v1 = findperpendicular(forward);
 205                 v2 = cross(forward, v1);
 206                 // random point on unit circle
 207                 dx = random() * 2 * M_PI;
 208                 dy = sin(dx);
 209                 dx = cos(dx);
 210                 // radius in our dist function
 211                 r = random();
 212                 r = solve_cubic_abcd(-2, 3, 0, -r) * '0 1 0';
 213                 return normalize(forward + (v1 * dx + v2 * dy) * r * sigma);
 214         }
 215         else if(sstyle == 6) // 1-r^2
 216         {
 217                 sigma = spread * 1.095445115010332; // match baseline stddev
 218                 v1 = findperpendicular(forward);
 219                 v2 = cross(forward, v1);
 220                 // random point on unit circle
 221                 dx = random() * 2 * M_PI;
 222                 dy = sin(dx);
 223                 dx = cos(dx);
 224                 // radius in our dist function
 225                 r = random();
 226                 r = sqrt(1 - r);
 227                 r = sqrt(1 - r);
 228                 return normalize(forward + (v1 * dx + v2 * dy) * r * sigma);
 229         }
 230         else if(sstyle == 7) // (1-r) (2-r)
 231         {
 232                 sigma = spread * 1.224744871391589; // match baseline stddev
 233                 v1 = findperpendicular(forward);
 234                 v2 = cross(forward, v1);
 235                 // random point on unit circle
 236                 dx = random() * 2 * M_PI;
 237                 dy = sin(dx);
 238                 dx = cos(dx);
 239                 // radius in our dist function
 240                 r = random();
 241                 r = 1 - sqrt(r);
 242                 r = 1 - sqrt(r);
 243                 return normalize(forward + (v1 * dx + v2 * dy) * r * sigma);
 244         }
 245         else
 246                 error("g_projectiles_spread_style must be 0 (sphere), 1 (flattened sphere), 2 (circle), 3 (gauss 3D), 4 (gauss plane), 5 (linear falloff), 6 (quadratic falloff), 7 (stronger falloff)!");
 247         return '0 0 0';
 248         /*
 249          * how to derive falloff functions:
 250          * rho(r) := (2-r) * (1-r);
 251          * a : 0;
 252          * b : 1;
 253          * rhor(r) := r * rho(r);
 254          * cr(t) := integrate(rhor(r), r, a, t);
 255          * scr(t) := integrate(rhor(r) * r^2, r, a, t);
 256          * variance : scr(b) / cr(b);
 257          * solve(cr(r) = rand * cr(b), r), programmmode:false;
 258          * sqrt(0.4 / variance), numer;
 259          */
 260 }