qcsrc/common/weapons/calculations.qc

   1 #include "calculations.qh"
   2
   3 // =============================
   4 //  Explosion Force Calculation
   5 // =============================
   6
   7 float explosion_calcpush_getmultiplier(vector explosion_v, vector target_v)
   8 {
   9         float a;
  10         a  = explosion_v * (explosion_v - target_v);
  11
  12         if(a <= 0)
  13                 // target is too fast to be hittable by this
  14                 return 0;
  15
  16         a /= (explosion_v * explosion_v);
  17                 // we know we can divide by this, or above a would be == 0
  18
  19         return a;
  20 }
  21
  22 #if 0
  23 vector explosion_calcpush(vector explosion_v, float explosion_m, vector target_v, float target_m, float elasticity)
  24 {
  25         // solution of the equations:
  26         //    v'                = v + a vp             // central hit
  27         //    m*v'   + mp*vp'   = m*v + mp*vp          // conservation of momentum
  28         //    m*v'^2 + mp*vp'^2 = m*v^2 + mp*vp^2      // conservation of energy (ELASTIC hit)
  29         // -> a = 0                                    // case 1: did not hit
  30         // -> a = 2*mp*(vp^2 - vp.v) / ((m+mp) * vp^2) // case 2: did hit
  31         //                                             // non-elastic hits are somewhere between these two
  32
  33         // this would be physically correct, but we don't do that
  34         return explosion_v * explosion_calcpush_getmultiplier(explosion_v, target_v) * (
  35                 (1 + elasticity) * (
  36                         explosion_m
  37                 ) / (
  38                         target_m + explosion_m
  39                 )
  40         ); // note: this factor is at least 0, at most 2
  41 }
  42 #endif
  43
  44 // simplified formula, tuned so that if the target has velocity 0, we get exactly the original force
  45 vector damage_explosion_calcpush(vector explosion_f, vector target_v, float speedfactor)
  46 {
  47         // if below 1, the formulas make no sense (and would cause superjumps)
  48         if(speedfactor < 1)
  49                 return explosion_f;
  50
  51 #if 0
  52         float m;
  53         // find m so that
  54         //   speedfactor * (1 + e) * m / (1 + m) == 1
  55         m = 1 / ((1 + 0) * speedfactor - 1);
  56         vector v;
  57         v = explosion_calcpush(explosion_f * speedfactor, m, target_v, 1, 0);
  58         // the factor we then get is:
  59         //   1
  60         LOG_INFOF("MASS: %f\nv: %v -> %v\nENERGY BEFORE == %f + %f = %f\nENERGY AFTER >= %f\n",
  61                 m,
  62                 target_v, target_v + v,
  63                 target_v * target_v, m * explosion_f * speedfactor * explosion_f * speedfactor, target_v * target_v + m * explosion_f * speedfactor * explosion_f * speedfactor,
  64                 (target_v + v) * (target_v + v));
  65         return v;
  66 #endif
  67         return explosion_f * explosion_calcpush_getmultiplier(explosion_f * speedfactor, target_v);
  68 }
  69
  70
  71 // =========================
  72 //  Shot Spread Calculation
  73 // =========================
  74
  75 vector cliptoplane(vector v, vector p)
  76 {
  77         return v - (v * p) * p;
  78 }
  79
  80 vector solve_cubic_pq(float p, float q)
  81 {
  82         float D, u, v, a;
  83         D = q*q/4.0 + p*p*p/27.0;
  84         if(D < 0)
  85         {
  86                 // irreducibilis
  87                 a = 1.0/3.0 * acos(-q/2.0 * sqrt(-27.0/(p*p*p)));
  88                 u = sqrt(-4.0/3.0 * p);
  89                 // a in range 0..pi/3
  90                 // cos(a)
  91                 // cos(a + 2pi/3)
  92                 // cos(a + 4pi/3)
  93                 return
  94                         u *
  95                         (
  96                                 '1 0 0' * cos(a + 2.0/3.0*M_PI)
  97                                 +
  98                                 '0 1 0' * cos(a + 4.0/3.0*M_PI)
  99                                 +
 100                                 '0 0 1' * cos(a)
 101                         );
 102         }
 103         else if(D == 0)
 104         {
 105                 // simple
 106                 if(p == 0)
 107                         return '0 0 0';
 108                 u = 3*q/p;
 109                 v = -u/2;
 110                 if(u >= v)
 111                         return '1 1 0' * v + '0 0 1' * u;
 112                 else
 113                         return '0 1 1' * v + '1 0 0' * u;
 114         }
 115         else
 116         {
 117                 // cardano
 118                 u = cbrt(-q/2.0 + sqrt(D));
 119                 v = cbrt(-q/2.0 - sqrt(D));
 120                 return '1 1 1' * (u + v);
 121         }
 122 }
 123 vector solve_cubic_abcd(float a, float b, float c, float d)
 124 {
 125         // y = 3*a*x + b
 126         // x = (y - b) / 3a
 127         float p, q;
 128         vector v;
 129         p = (9*a*c - 3*b*b);
 130         q = (27*a*a*d - 9*a*b*c + 2*b*b*b);
 131         v = solve_cubic_pq(p, q);
 132         v = (v -  b * '1 1 1') * (1.0 / (3.0 * a));
 133         if(a < 0)
 134                 v += '1 0 -1' * (v.z - v.x); // swap x, z
 135         return v;
 136 }
 137
 138 vector findperpendicular(vector v)
 139 {
 140         vector p;
 141         p.x = v.z;
 142         p.y = -v.x;
 143         p.z = v.y;
 144         return normalize(cliptoplane(p, v));
 145 }
 146
 147 #ifdef SVQC
 148         int W_GunAlign(entity this, int preferred_align)
 149         {
 150                 if(this.m_gunalign)
 151                         return this.m_gunalign; // no adjustment needed
 152
 153                 entity own = this.owner;
 154
 155                 if(preferred_align < 1 || preferred_align > 4)
 156                         preferred_align = 3; // default
 157
 158                 for(int j = 4; j > 1; --j) // > 1 as 1 is just center again
 159                 {
 160                         int taken = 0;
 161                         for(int slot = 0; slot < MAX_WEAPONSLOTS; ++slot)
 162                         {
 163                                 .entity weaponentity = weaponentities[slot];
 164                                 if(own.(weaponentity).m_gunalign == j) // we know it can't be ours thanks to the above check
 165                                         taken |= BIT(j);
 166                                 if(own.(weaponentity).m_gunalign == preferred_align)
 167                                         taken |= BIT(preferred_align);
 168                         }
 169
 170                         if(!(taken & BIT(preferred_align)))
 171                                 return preferred_align; // prefer the recommended
 172                         if(!(taken & BIT(j)))
 173                                 return j; // or fall back if it's not available
 174                 }
 175
 176                 return preferred_align; // return it anyway
 177         }
 178 #else
 179         int W_GunAlign(entity this, int preferred_align)
 180         {
 181                 return this.m_gunalign > 0 ? this.m_gunalign : preferred_align;
 182         }
 183 #endif
 184
 185 #if 0
 186 int W_GetGunAlignment(entity player)
 187 {
 188         int gunalign = STAT(GUNALIGN, player);
 189         if(gunalign < 1 || gunalign > 4)
 190                 gunalign = 3; // default value
 191         --gunalign;
 192
 193         return gunalign;
 194 }
 195 #endif
 196
 197 vector W_CalculateSpread(vector forward, float spread, float spreadfactor, float spreadstyle)
 198 {
 199         float sigma;
 200         vector v1 = '0 0 0', v2;
 201         float dx, dy, r;
 202         spread *= spreadfactor; //g_weaponspreadfactor;
 203         if(spread <= 0)
 204                 return forward;
 205
 206         switch(spreadstyle)
 207         {
 208                 case 0:
 209                 {
 210                         // this is the baseline for the spread value!
 211                         // standard deviation: sqrt(2/5)
 212                         // density function: sqrt(1-r^2)
 213                         return forward + randomvec() * spread;
 214                 }
 215                 case 1:
 216                 {
 217                         // same thing, basically
 218                         return normalize(forward + cliptoplane(randomvec() * spread, forward));
 219                 }
 220                 case 2:
 221                 {
 222                         // circle spread... has at sigma=1 a standard deviation of sqrt(1/2)
 223                         sigma = spread * 0.89442719099991587855; // match baseline stddev
 224                         v1 = findperpendicular(forward);
 225                         v2 = cross(forward, v1);
 226                         // random point on unit circle
 227                         dx = random() * 2 * M_PI;
 228                         dy = sin(dx);
 229                         dx = cos(dx);
 230                         // radius in our dist function
 231                         r = random();
 232                         r = sqrt(r);
 233                         return normalize(forward + (v1 * dx + v2 * dy) * r * sigma);
 234                 }
 235                 case 3: // gauss 3d
 236                 {
 237                         sigma = spread * 0.44721359549996; // match baseline stddev
 238                         // note: 2D gaussian has sqrt(2) times the stddev of 1D, so this factor is right
 239                         v1 = forward;
 240                         v1_x += gsl_ran_gaussian(sigma);
 241                         v1_y += gsl_ran_gaussian(sigma);
 242                         v1_z += gsl_ran_gaussian(sigma);
 243                         return v1;
 244                 }
 245                 case 4: // gauss 2d
 246                 {
 247                         sigma = spread * 0.44721359549996; // match baseline stddev
 248                         // note: 2D gaussian has sqrt(2) times the stddev of 1D, so this factor is right
 249                         v1_x = gsl_ran_gaussian(sigma);
 250                         v1_y = gsl_ran_gaussian(sigma);
 251                         v1_z = gsl_ran_gaussian(sigma);
 252                         return normalize(forward + cliptoplane(v1, forward));
 253                 }
 254                 case 5: // 1-r
 255                 {
 256                         sigma = spread * 1.154700538379252; // match baseline stddev
 257                         v1 = findperpendicular(forward);
 258                         v2 = cross(forward, v1);
 259                         // random point on unit circle
 260                         dx = random() * 2 * M_PI;
 261                         dy = sin(dx);
 262                         dx = cos(dx);
 263                         // radius in our dist function
 264                         r = random();
 265                         r = solve_cubic_abcd(-2, 3, 0, -r) * '0 1 0';
 266                         return normalize(forward + (v1 * dx + v2 * dy) * r * sigma);
 267                 }
 268                 case 6: // 1-r^2
 269                 {
 270                         sigma = spread * 1.095445115010332; // match baseline stddev
 271                         v1 = findperpendicular(forward);
 272                         v2 = cross(forward, v1);
 273                         // random point on unit circle
 274                         dx = random() * 2 * M_PI;
 275                         dy = sin(dx);
 276                         dx = cos(dx);
 277                         // radius in our dist function
 278                         r = random();
 279                         r = sqrt(1 - r);
 280                         r = sqrt(1 - r);
 281                         return normalize(forward + (v1 * dx + v2 * dy) * r * sigma);
 282                 }
 283                 case 7: // (1-r) (2-r)
 284                 {
 285                         sigma = spread * 1.224744871391589; // match baseline stddev
 286                         v1 = findperpendicular(forward);
 287                         v2 = cross(forward, v1);
 288                         // random point on unit circle
 289                         dx = random() * 2 * M_PI;
 290                         dy = sin(dx);
 291                         dx = cos(dx);
 292                         // radius in our dist function
 293                         r = random();
 294                         r = 1 - sqrt(r);
 295                         r = 1 - sqrt(r);
 296                         return normalize(forward + (v1 * dx + v2 * dy) * r * sigma);
 297                 }
 298                 default:
 299                         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)!");
 300         }
 301
 302         return '0 0 0';
 303         /*
 304          * how to derive falloff functions:
 305          * rho(r) := (2-r) * (1-r);
 306          * a : 0;
 307          * b : 1;
 308          * rhor(r) := r * rho(r);
 309          * cr(t) := integrate(rhor(r), r, a, t);
 310          * scr(t) := integrate(rhor(r) * r^2, r, a, t);
 311          * variance : scr(b) / cr(b);
 312          * solve(cr(r) = rand * cr(b), r), programmmode:false;
 313          * sqrt(0.4 / variance), numer;
 314          */
 315 }