pipeline: prefer to skip gmqcc build by downloading from beta.xonotic.org and update...
[xonotic/xonotic-data.pk3dir.git] / 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",
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 u * vec3(
94                         cos(a + 2.0/3.0*M_PI),
95                         cos(a + 4.0/3.0*M_PI),
96                         cos(a)
97                 );
98         }
99         else if(D == 0)
100         {
101                 // simple
102                 if(p == 0)
103                         return '0 0 0';
104                 u = 3*q/p;
105                 v = -u/2;
106                 return (u >= v) ? vec3(v, v, u) : vec3(u, v, v);
107         }
108         else
109         {
110                 // cardano
111                 //u = cbrt(-q/2.0 + sqrt(D));
112                 //v = cbrt(-q/2.0 - sqrt(D));
113                 a = cbrt(-q/2.0 + sqrt(D)) + cbrt(-q/2.0 - sqrt(D));
114                 return vec3(a, a, a);
115         }
116 }
117 vector solve_cubic_abcd(float a, float b, float c, float d)
118 {
119         // y = 3*a*x + b
120         // x = (y - b) / 3a
121         float p, q;
122         vector v;
123         p = (9*a*c - 3*b*b);
124         q = (27*a*a*d - 9*a*b*c + 2*b*b*b);
125         v = solve_cubic_pq(p, q);
126         v = (v -  b * '1 1 1') * (1.0 / (3.0 * a));
127         if(a < 0)
128                 v += '1 0 -1' * (v.z - v.x); // swap x, z
129         return v;
130 }
131
132 vector findperpendicular(vector v)
133 {
134         return normalize(cliptoplane(vec3(v.z, -v.x, v.y), v));
135 }
136
137 #ifdef SVQC
138         int W_GunAlign(entity this, int preferred_align)
139         {
140                 if(this.m_gunalign)
141                         return this.m_gunalign; // no adjustment needed
142
143                 entity own = this.owner;
144
145                 if(preferred_align < 1 || preferred_align > 4)
146                         preferred_align = 3; // default
147
148                 for(int j = 4; j > 1; --j) // > 1 as 1 is just center again
149                 {
150                         int taken = 0;
151                         for(int slot = 0; slot < MAX_WEAPONSLOTS; ++slot)
152                         {
153                                 .entity weaponentity = weaponentities[slot];
154                                 if(own.(weaponentity).m_gunalign == j) // we know it can't be ours thanks to the above check
155                                         taken |= BIT(j);
156                                 if(own.(weaponentity).m_gunalign == preferred_align)
157                                         taken |= BIT(preferred_align);
158                         }
159
160                         if(!(taken & BIT(preferred_align)))
161                                 return preferred_align; // prefer the recommended
162                         if(!(taken & BIT(j)))
163                                 return j; // or fall back if it's not available
164                 }
165
166                 return preferred_align; // return it anyway
167         }
168 #else
169         int W_GunAlign(entity this, int preferred_align)
170         {
171                 return this.m_gunalign > 0 ? this.m_gunalign : preferred_align;
172         }
173 #endif
174
175 #if 0
176 int W_GetGunAlignment(entity player)
177 {
178         int gunalign = STAT(GUNALIGN, player);
179         if(gunalign < 1 || gunalign > 4)
180                 gunalign = 3; // default value
181         --gunalign;
182
183         return gunalign;
184 }
185 #endif
186
187 vector W_CalculateSpread(vector forward, float spread, float spreadfactor, float spreadstyle)
188 {
189         float sigma;
190         vector v1 = '0 0 0', v2;
191         float dx, dy, r;
192         spread *= spreadfactor; //autocvar_g_weaponspreadfactor;
193         if(spread <= 0)
194                 return forward;
195
196         switch(spreadstyle)
197         {
198                 case 0:
199                 {
200                         // this is the baseline for the spread value!
201                         // standard deviation: sqrt(2/5)
202                         // density function: sqrt(1-r^2)
203                         return forward + randomvec() * spread;
204                 }
205                 case 1:
206                 {
207                         // same thing, basically
208                         return normalize(forward + cliptoplane(randomvec() * spread, forward));
209                 }
210                 case 2:
211                 {
212                         // circle spread... has at sigma=1 a standard deviation of sqrt(1/2)
213                         sigma = spread * 0.89442719099991587855; // match baseline stddev
214                         v1 = findperpendicular(forward);
215                         v2 = cross(forward, v1);
216                         // random point on unit circle
217                         dx = random() * 2 * M_PI;
218                         dy = sin(dx);
219                         dx = cos(dx);
220                         // radius in our dist function
221                         r = random();
222                         r = sqrt(r);
223                         return normalize(forward + (v1 * dx + v2 * dy) * r * sigma);
224                 }
225                 case 3: // gauss 3d
226                 {
227                         sigma = spread * 0.44721359549996; // match baseline stddev
228                         // note: 2D gaussian has sqrt(2) times the stddev of 1D, so this factor is right
229                         v1 = forward;
230                         v1_x += gsl_ran_gaussian(sigma);
231                         v1_y += gsl_ran_gaussian(sigma);
232                         v1_z += gsl_ran_gaussian(sigma);
233                         return v1;
234                 }
235                 case 4: // gauss 2d
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_x = gsl_ran_gaussian(sigma);
240                         v1_y = gsl_ran_gaussian(sigma);
241                         v1_z = gsl_ran_gaussian(sigma);
242                         return normalize(forward + cliptoplane(v1, forward));
243                 }
244                 case 5: // 1-r
245                 {
246                         sigma = spread * 1.154700538379252; // match baseline stddev
247                         v1 = findperpendicular(forward);
248                         v2 = cross(forward, v1);
249                         // random point on unit circle
250                         dx = random() * 2 * M_PI;
251                         dy = sin(dx);
252                         dx = cos(dx);
253                         // radius in our dist function
254                         r = random();
255                         r = solve_cubic_abcd(-2, 3, 0, -r) * '0 1 0';
256                         return normalize(forward + (v1 * dx + v2 * dy) * r * sigma);
257                 }
258                 case 6: // 1-r^2
259                 {
260                         sigma = spread * 1.095445115010332; // match baseline stddev
261                         v1 = findperpendicular(forward);
262                         v2 = cross(forward, v1);
263                         // random point on unit circle
264                         dx = random() * 2 * M_PI;
265                         dy = sin(dx);
266                         dx = cos(dx);
267                         // radius in our dist function
268                         r = random();
269                         r = sqrt(1 - r);
270                         r = sqrt(1 - r);
271                         return normalize(forward + (v1 * dx + v2 * dy) * r * sigma);
272                 }
273                 case 7: // (1-r) (2-r)
274                 {
275                         sigma = spread * 1.224744871391589; // match baseline stddev
276                         v1 = findperpendicular(forward);
277                         v2 = cross(forward, v1);
278                         // random point on unit circle
279                         dx = random() * 2 * M_PI;
280                         dy = sin(dx);
281                         dx = cos(dx);
282                         // radius in our dist function
283                         r = random();
284                         r = 1 - sqrt(r);
285                         r = 1 - sqrt(r);
286                         return normalize(forward + (v1 * dx + v2 * dy) * r * sigma);
287                 }
288                 default:
289                         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)!");
290         }
291
292         return '0 0 0';
293         /*
294          * how to derive falloff functions:
295          * rho(r) := (2-r) * (1-r);
296          * a : 0;
297          * b : 1;
298          * rhor(r) := r * rho(r);
299          * cr(t) := integrate(rhor(r), r, a, t);
300          * scr(t) := integrate(rhor(r) * r^2, r, a, t);
301          * variance : scr(b) / cr(b);
302          * solve(cr(r) = rand * cr(b), r), programmmode:false;
303          * sqrt(0.4 / variance), numer;
304          */
305 }