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