3
// Copyright 1997-1998 by David K. McAllister
5
// I used code Copyright 1997 by Jonathan P. Leech
6
// as an example in implenting this.
8
// This file implements the dynamics of particle actions.
14
#define SQRT2PI 2.506628274631000502415765284811045253006
15
#define ONEOVERSQRT2PI (1. / SQRT2PI)
17
// To offset [0 .. 1] vectors to [-.5 .. .5]
18
static pVector vHalf(0.5, 0.5, 0.5);
20
static inline pVector RandVec()
22
return pVector(drand48(), drand48(), drand48());
25
// Return a random number with a normal distribution.
26
static inline float NRand(float sigma = 1.0f)
28
#define ONE_OVER_SIGMA_EXP (1.0f / 0.7975f)
30
if(sigma == 0) return 0;
37
while(drand48() > expf(-fsqr(y - 1.0f)*0.5f));
40
return y * sigma * ONE_OVER_SIGMA_EXP;
42
return -y * sigma * ONE_OVER_SIGMA_EXP;
45
void PAAvoid::Execute(ParticleGroup *group)
47
float magdt = magnitude * dt;
53
if(look_ahead < P_MAXFLOAT)
55
for(int i = 0; i < group->p_count; i++)
57
Particle &m = group->list[i];
59
// p2 stores the plane normal (the a,b,c of the plane eqn).
60
// Old and new distances: dist(p,plane) = n * p + d
61
// radius1 stores -n*p, which is d.
62
float dist = m.pos * position.p2 + position.radius1;
66
float vm = m.vel.length();
67
pVector Vn = m.vel / vm;
68
// float dot = Vn * position.p2;
70
pVector tmp = (position.p2 * (magdt / (dist*dist+epsilon))) + Vn;
71
m.vel = tmp * (vm / tmp.length());
77
for(int i = 0; i < group->p_count; i++)
79
Particle &m = group->list[i];
81
// p2 stores the plane normal (the a,b,c of the plane eqn).
82
// Old and new distances: dist(p,plane) = n * p + d
83
// radius1 stores -n*p, which is d.
84
float dist = m.pos * position.p2 + position.radius1;
86
float vm = m.vel.length();
87
pVector Vn = m.vel / vm;
88
// float dot = Vn * position.p2;
90
pVector tmp = (position.p2 * (magdt / (dist*dist+epsilon))) + Vn;
91
m.vel = tmp * (vm / tmp.length());
98
// Compute the inverse matrix of the plane basis.
99
pVector &u = position.u;
100
pVector &v = position.v;
102
// The normalized bases are needed inside the loop.
103
pVector un = u / position.radius1Sqr;
104
pVector vn = v / position.radius2Sqr;
107
float wx = u.y*v.z-u.z*v.y;
108
float wy = u.z*v.x-u.x*v.z;
109
float wz = u.x*v.y-u.y*v.x;
111
float det = 1/(wz*u.x*v.y-wz*u.y*v.x-u.z*wx*v.y-u.x*v.z*wy+v.z*wx*u.y+u.z*v.x*wy);
113
pVector s1((v.y*wz-v.z*wy), (v.z*wx-v.x*wz), (v.x*wy-v.y*wx));
115
pVector s2((u.y*wz-u.z*wy), (u.z*wx-u.x*wz), (u.x*wy-u.y*wx));
118
// See which particles bounce.
119
for(int i = 0; i < group->p_count; i++)
121
Particle &m = group->list[i];
123
// See if particle's current and next positions cross plane.
124
// If not, couldn't bounce, so keep going.
125
pVector pnext(m.pos + m.vel * dt * look_ahead);
127
// p2 stores the plane normal (the a,b,c of the plane eqn).
128
// Old and new distances: dist(p,plane) = n * p + d
129
// radius1 stores -n*p, which is d.
130
float distold = m.pos * position.p2 + position.radius1;
131
float distnew = pnext * position.p2 + position.radius1;
133
// Opposite signs if product < 0
134
// There is no faster way to do this.
135
if(distold * distnew >= 0)
138
float nv = position.p2 * m.vel;
139
float t = -distold / nv;
141
// Actual intersection point p(t) = pos + vel t
142
pVector phit(m.pos + m.vel * t);
144
// Offset from origin in plane, p - origin
145
pVector offset(phit - position.p1);
147
// Dot product with basis vectors of old frame
148
// in terms of new frame gives position in uv frame.
149
float upos = offset * s1;
150
float vpos = offset * s2;
152
// Did it cross plane outside triangle?
153
if(upos < 0 || vpos < 0 || upos > 1 || vpos > 1)
156
// A hit! A most palpable hit!
157
// Compute distance to the three edges.
158
pVector uofs = (un * (un * offset)) - offset;
159
float udistSqr = uofs.length2();
160
pVector vofs = (vn * (vn * offset)) - offset;
161
float vdistSqr = vofs.length2();
163
pVector foffset((u + v) - offset);
164
pVector fofs = (un * (un * foffset)) - foffset;
165
float fdistSqr = fofs.length2();
166
pVector gofs = (un * (un * foffset)) - foffset;
167
float gdistSqr = gofs.length2();
170
if(udistSqr <= vdistSqr && udistSqr <= fdistSqr
171
&& udistSqr <= gdistSqr) S = uofs;
172
else if(vdistSqr <= fdistSqr && vdistSqr <= gdistSqr) S = vofs;
173
else if(fdistSqr <= gdistSqr) S = fofs;
178
// We now have a vector to safety.
179
float vm = m.vel.length();
180
pVector Vn = m.vel / vm;
183
pVector tmp = (S * (magdt / (t*t+epsilon))) + Vn;
184
m.vel = tmp * (vm / tmp.length());
190
// Compute the inverse matrix of the plane basis.
191
pVector &u = position.u;
192
pVector &v = position.v;
194
// The normalized bases are needed inside the loop.
195
pVector un = u / position.radius1Sqr;
196
pVector vn = v / position.radius2Sqr;
198
// f is the third (non-basis) triangle edge.
204
float wx = u.y*v.z-u.z*v.y;
205
float wy = u.z*v.x-u.x*v.z;
206
float wz = u.x*v.y-u.y*v.x;
208
float det = 1/(wz*u.x*v.y-wz*u.y*v.x-u.z*wx*v.y-u.x*v.z*wy+v.z*wx*u.y+u.z*v.x*wy);
210
pVector s1((v.y*wz-v.z*wy), (v.z*wx-v.x*wz), (v.x*wy-v.y*wx));
212
pVector s2((u.y*wz-u.z*wy), (u.z*wx-u.x*wz), (u.x*wy-u.y*wx));
215
// See which particles bounce.
216
for(int i = 0; i < group->p_count; i++)
218
Particle &m = group->list[i];
220
// See if particle's current and next positions cross plane.
221
// If not, couldn't bounce, so keep going.
222
pVector pnext(m.pos + m.vel * dt * look_ahead);
224
// p2 stores the plane normal (the a,b,c of the plane eqn).
225
// Old and new distances: dist(p,plane) = n * p + d
226
// radius1 stores -n*p, which is d.
227
float distold = m.pos * position.p2 + position.radius1;
228
float distnew = pnext * position.p2 + position.radius1;
230
// Opposite signs if product < 0
231
// Is there a faster way to do this?
232
if(distold * distnew >= 0)
235
float nv = position.p2 * m.vel;
236
float t = -distold / nv;
238
// Actual intersection point p(t) = pos + vel t
239
pVector phit(m.pos + m.vel * t);
241
// Offset from origin in plane, p - origin
242
pVector offset(phit - position.p1);
244
// Dot product with basis vectors of old frame
245
// in terms of new frame gives position in uv frame.
246
float upos = offset * s1;
247
float vpos = offset * s2;
249
// Did it cross plane outside triangle?
250
if(upos < 0 || vpos < 0 || (upos + vpos) > 1)
253
// A hit! A most palpable hit!
254
// Compute distance to the three edges.
255
pVector uofs = (un * (un * offset)) - offset;
256
float udistSqr = uofs.length2();
257
pVector vofs = (vn * (vn * offset)) - offset;
258
float vdistSqr = vofs.length2();
259
pVector foffset(offset - u);
260
pVector fofs = (fn * (fn * foffset)) - foffset;
261
float fdistSqr = fofs.length2();
263
if(udistSqr <= vdistSqr && udistSqr <= fdistSqr) S = uofs;
264
else if(vdistSqr <= fdistSqr) S = vofs;
269
// We now have a vector to safety.
270
float vm = m.vel.length();
271
pVector Vn = m.vel / vm;
274
pVector tmp = (S * (magdt / (t*t+epsilon))) + Vn;
275
m.vel = tmp * (vm / tmp.length());
281
float r1Sqr = fsqr(position.radius1);
282
float r2Sqr = fsqr(position.radius2);
284
// See which particles bounce.
285
for(int i = 0; i < group->p_count; i++)
287
Particle &m = group->list[i];
289
// See if particle's current and next positions cross plane.
290
// If not, couldn't bounce, so keep going.
291
pVector pnext(m.pos + m.vel * dt * look_ahead);
293
// p2 stores the plane normal (the a,b,c of the plane eqn).
294
// Old and new distances: dist(p,plane) = n * p + d
295
// radius1 stores -n*p, which is d. radius1Sqr stores d.
296
float distold = m.pos * position.p2 + position.radius1Sqr;
297
float distnew = pnext * position.p2 + position.radius1Sqr;
299
// Opposite signs if product < 0
300
// Is there a faster way to do this?
301
if(distold * distnew >= 0)
304
// Find position at the crossing point by parameterizing
305
// p(t) = pos + vel * t
306
// Solve dist(p(t),plane) = 0 e.g.
307
// n * p(t) + D = 0 ->
308
// n * p + t (n * v) + D = 0 ->
309
// t = -(n * p + D) / (n * v)
310
// Could factor n*v into distnew = distold + n*v and save a bit.
311
// Safe since n*v != 0 assured by quick rejection test.
312
// This calc is indep. of dt because we have established that it
313
// will hit before dt. We just want to know when.
314
float nv = position.p2 * m.vel;
315
float t = -distold / nv;
317
// Actual intersection point p(t) = pos + vel t
318
pVector phit(m.pos + m.vel * t);
320
// Offset from origin in plane, phit - origin
321
pVector offset(phit - position.p1);
323
float rad = offset.length2();
325
if(rad > r1Sqr || rad < r2Sqr)
328
// A hit! A most palpable hit!
332
// We now have a vector to safety.
333
float vm = m.vel.length();
334
pVector Vn = m.vel / vm;
337
pVector tmp = (S * (magdt / (t*t+epsilon))) + Vn;
338
m.vel = tmp * (vm / tmp.length());
344
float rSqr = position.radius1 * position.radius1;
346
// See which particles are aimed toward the sphere.
347
for(int i = 0; i < group->p_count; i++)
349
Particle &m = group->list[i];
351
// First do a ray-sphere intersection test and
352
// see if it's soon enough.
353
// Can I do this faster without t?
354
float vm = m.vel.length();
355
pVector Vn = m.vel / vm;
357
pVector L = position.p1 - m.pos;
360
float disc = rSqr - (L * L) + v * v;
362
continue; // I'm not heading toward it.
364
// Compute length for second rejection test.
365
float t = v - sqrtf(disc);
366
if(t < 0 || t > (vm * look_ahead))
369
// Get a vector to safety.
375
pVector tmp = (S * (magdt / (t*t+epsilon))) + Vn;
376
m.vel = tmp * (vm / tmp.length());
385
void PABounce::Execute(ParticleGroup *group)
387
switch(position.type)
391
// Compute the inverse matrix of the plane basis.
392
pVector &u = position.u;
393
pVector &v = position.v;
396
float wx = u.y*v.z-u.z*v.y;
397
float wy = u.z*v.x-u.x*v.z;
398
float wz = u.x*v.y-u.y*v.x;
400
float det = 1/(wz*u.x*v.y-wz*u.y*v.x-u.z*wx*v.y-u.x*v.z*wy+v.z*wx*u.y+u.z*v.x*wy);
402
pVector s1((v.y*wz-v.z*wy), (v.z*wx-v.x*wz), (v.x*wy-v.y*wx));
404
pVector s2((u.y*wz-u.z*wy), (u.z*wx-u.x*wz), (u.x*wy-u.y*wx));
407
// See which particles bounce.
408
for(int i = 0; i < group->p_count; i++)
410
Particle &m = group->list[i];
412
// See if particle's current and next positions cross plane.
413
// If not, couldn't bounce, so keep going.
414
pVector pnext(m.pos + m.vel * dt);
416
// p2 stores the plane normal (the a,b,c of the plane eqn).
417
// Old and new distances: dist(p,plane) = n * p + d
418
// radius1 stores -n*p, which is d.
419
float distold = m.pos * position.p2 + position.radius1;
420
float distnew = pnext * position.p2 + position.radius1;
422
// Opposite signs if product < 0
423
// Is there a faster way to do this?
424
if(distold * distnew >= 0)
427
// Find position at the crossing point by parameterizing
428
// p(t) = pos + vel * t
429
// Solve dist(p(t),plane) = 0 e.g.
430
// n * p(t) + D = 0 ->
431
// n * p + t (n * v) + D = 0 ->
432
// t = -(n * p + D) / (n * v)
433
// Could factor n*v into distnew = distold + n*v and save a bit.
434
// Safe since n*v != 0 assured by quick rejection test.
435
// This calc is indep. of dt because we have established that it
436
// will hit before dt. We just want to know when.
437
float nv = position.p2 * m.vel;
438
float t = -distold / nv;
440
// Actual intersection point p(t) = pos + vel t
441
pVector phit(m.pos + m.vel * t);
443
// Offset from origin in plane, p - origin
444
pVector offset(phit - position.p1);
446
// Dot product with basis vectors of old frame
447
// in terms of new frame gives position in uv frame.
448
float upos = offset * s1;
449
float vpos = offset * s2;
451
// Did it cross plane outside triangle?
452
if(upos < 0 || vpos < 0 || (upos + vpos) > 1)
455
// A hit! A most palpable hit!
457
// Compute tangential and normal components of velocity
458
pVector vn(position.p2 * nv); // Normal Vn = (V.N)N
459
pVector vt(m.vel - vn); // Tangent Vt = V - Vn
461
// Compute new velocity heading out:
462
// Don't apply friction if tangential velocity < cutoff
463
if(vt.length2() <= cutoffSqr)
464
m.vel = vt - vn * resilience;
466
m.vel = vt * oneMinusFriction - vn * resilience;
472
float r1Sqr = fsqr(position.radius1);
473
float r2Sqr = fsqr(position.radius2);
475
// See which particles bounce.
476
for(int i = 0; i < group->p_count; i++)
478
Particle &m = group->list[i];
480
// See if particle's current and next positions cross plane.
481
// If not, couldn't bounce, so keep going.
482
pVector pnext(m.pos + m.vel * dt);
484
// p2 stores the plane normal (the a,b,c of the plane eqn).
485
// Old and new distances: dist(p,plane) = n * p + d
486
// radius1 stores -n*p, which is d. radius1Sqr stores d.
487
float distold = m.pos * position.p2 + position.radius1Sqr;
488
float distnew = pnext * position.p2 + position.radius1Sqr;
490
// Opposite signs if product < 0
491
// Is there a faster way to do this?
492
if(distold * distnew >= 0)
495
// Find position at the crossing point by parameterizing
496
// p(t) = pos + vel * t
497
// Solve dist(p(t),plane) = 0 e.g.
498
// n * p(t) + D = 0 ->
499
// n * p + t (n * v) + D = 0 ->
500
// t = -(n * p + D) / (n * v)
501
// Could factor n*v into distnew = distold + n*v and save a bit.
502
// Safe since n*v != 0 assured by quick rejection test.
503
// This calc is indep. of dt because we have established that it
504
// will hit before dt. We just want to know when.
505
float nv = position.p2 * m.vel;
506
float t = -distold / nv;
508
// Actual intersection point p(t) = pos + vel t
509
pVector phit(m.pos + m.vel * t);
511
// Offset from origin in plane, phit - origin
512
pVector offset(phit - position.p1);
514
float rad = offset.length2();
516
if(rad > r1Sqr || rad < r2Sqr)
519
// A hit! A most palpable hit!
521
// Compute tangential and normal components of velocity
522
pVector vn(position.p2 * nv); // Normal Vn = (V.N)N
523
pVector vt(m.vel - vn); // Tangent Vt = V - Vn
525
// Compute new velocity heading out:
526
// Don't apply friction if tangential velocity < cutoff
527
if(vt.length2() <= cutoffSqr)
528
m.vel = vt - vn * resilience;
530
m.vel = vt * oneMinusFriction - vn * resilience;
536
// See which particles bounce.
537
for(int i = 0; i < group->p_count; i++)
539
Particle &m = group->list[i];
541
// See if particle's current and next positions cross plane.
542
// If not, couldn't bounce, so keep going.
543
pVector pnext(m.pos + m.vel * dt);
545
// p2 stores the plane normal (the a,b,c of the plane eqn).
546
// Old and new distances: dist(p,plane) = n * p + d
547
// radius1 stores -n*p, which is d.
548
float distold = m.pos * position.p2 + position.radius1;
549
float distnew = pnext * position.p2 + position.radius1;
551
// Opposite signs if product < 0
552
if(distold * distnew >= 0)
555
// Compute tangential and normal components of velocity
556
float nmag = m.vel * position.p2;
557
pVector vn(position.p2 * nmag); // Normal Vn = (V.N)N
558
pVector vt(m.vel - vn); // Tangent Vt = V - Vn
560
// Compute new velocity heading out:
561
// Don't apply friction if tangential velocity < cutoff
562
if(vt.length2() <= cutoffSqr)
563
m.vel = vt - vn * resilience;
565
m.vel = vt * oneMinusFriction - vn * resilience;
571
// Compute the inverse matrix of the plane basis.
572
pVector &u = position.u;
573
pVector &v = position.v;
576
float wx = u.y*v.z-u.z*v.y;
577
float wy = u.z*v.x-u.x*v.z;
578
float wz = u.x*v.y-u.y*v.x;
580
float det = 1/(wz*u.x*v.y-wz*u.y*v.x-u.z*wx*v.y-u.x*v.z*wy+v.z*wx*u.y+u.z*v.x*wy);
582
pVector s1((v.y*wz-v.z*wy), (v.z*wx-v.x*wz), (v.x*wy-v.y*wx));
584
pVector s2((u.y*wz-u.z*wy), (u.z*wx-u.x*wz), (u.x*wy-u.y*wx));
587
// See which particles bounce.
588
for(int i = 0; i < group->p_count; i++)
590
Particle &m = group->list[i];
592
// See if particle's current and next positions cross plane.
593
// If not, couldn't bounce, so keep going.
594
pVector pnext(m.pos + m.vel * dt);
596
// p2 stores the plane normal (the a,b,c of the plane eqn).
597
// Old and new distances: dist(p,plane) = n * p + d
598
// radius1 stores -n*p, which is d.
599
float distold = m.pos * position.p2 + position.radius1;
600
float distnew = pnext * position.p2 + position.radius1;
602
// Opposite signs if product < 0
603
if(distold * distnew >= 0)
606
// Find position at the crossing point by parameterizing
607
// p(t) = pos + vel * t
608
// Solve dist(p(t),plane) = 0 e.g.
609
// n * p(t) + D = 0 ->
610
// n * p + t (n * v) + D = 0 ->
611
// t = -(n * p + D) / (n * v)
612
float t = -distold / (position.p2 * m.vel);
614
// Actual intersection point p(t) = pos + vel t
615
pVector phit(m.pos + m.vel * t);
617
// Offset from origin in plane, p - origin
618
pVector offset(phit - position.p1);
620
// Dot product with basis vectors of old frame
621
// in terms of new frame gives position in uv frame.
622
float upos = offset * s1;
623
float vpos = offset * s2;
625
// Crossed plane outside bounce region if !(0<=[uv]pos<=1)
626
if(upos < 0 || upos > 1 || vpos < 0 || vpos > 1)
629
// A hit! A most palpable hit!
631
// Compute tangential and normal components of velocity
632
float nmag = m.vel * position.p2;
633
pVector vn(position.p2 * nmag); // Normal Vn = (V.N)N
634
pVector vt(m.vel - vn); // Tangent Vt = V - Vn
636
// Compute new velocity heading out:
637
// Don't apply friction if tangential velocity < cutoff
638
if(vt.length2() <= cutoffSqr)
639
m.vel = vt - vn * resilience;
641
m.vel = vt * oneMinusFriction - vn * resilience;
647
// Sphere that particles bounce off
648
// The particles are always forced out of the sphere.
649
for(int i = 0; i < group->p_count; i++)
651
Particle &m = group->list[i];
653
// See if particle's next position is inside domain.
655
pVector pnext(m.pos + m.vel * dt);
657
if(position.Within(pnext))
659
// See if we were inside on previous timestep.
660
bool pinside = position.Within(m.pos);
662
// Normal to surface. This works for a sphere. Isn't
663
// computed quite right, should extrapolate particle
664
// position to surface.
665
pVector n(m.pos - position.p1);
668
// Compute tangential and normal components of velocity
669
float nmag = m.vel * n;
671
pVector vn(n * nmag); // Normal Vn = (V.N)N
672
pVector vt = m.vel - vn; // Tangent Vt = V - Vn
676
// Previous position was inside. If normal component of
677
// velocity points in, reverse it. This effectively
678
// repels particles which would otherwise be trapped
685
// Previous position was outside -> particle will cross
686
// surface boundary. Reverse normal component of velocity,
687
// and apply friction (if Vt >= cutoff) and resilience.
689
// Compute new velocity heading out:
690
// Don't apply friction if tangential velocity < cutoff
691
if(vt.length2() <= cutoffSqr)
692
m.vel = vt - vn * resilience;
694
m.vel = vt * oneMinusFriction - vn * resilience;
704
// Set the secondary position of each particle to be its position.
705
void PACallActionList::Execute(ParticleGroup *group)
707
pCallActionList(action_list_num);
710
// Set the secondary position of each particle to be its position.
711
void PACopyVertexB::Execute(ParticleGroup *group)
717
for(i = 0; i < group->p_count; i++)
719
Particle &m = group->list[i];
726
for(i = 0; i < group->p_count; i++)
728
Particle &m = group->list[i];
735
void PADamping::Execute(ParticleGroup *group)
737
// This is important if dt is != 1.
739
pVector scale(one - ((one - damping) * dt));
741
for(int i = 0; i < group->p_count; i++)
743
Particle &m = group->list[i];
744
float vSqr = m.vel.length2();
746
if(vSqr >= vlowSqr && vSqr <= vhighSqr)
755
// Exert force on each particle away from explosion center
756
void PAExplosion::Execute(ParticleGroup *group)
758
float radius = velocity * age;
759
float magdt = magnitude * dt;
760
float oneOverSigma = 1.0f / stdev;
761
float inexp = -0.5f*fsqr(oneOverSigma);
762
float outexp = ONEOVERSQRT2PI * oneOverSigma;
764
for(int i = 0; i < group->p_count; i++)
766
Particle &m = group->list[i];
768
// Figure direction to particle.
769
pVector dir(m.pos - center);
770
float distSqr = dir.length2();
771
float dist = sqrtf(distSqr);
772
float DistFromWaveSqr = fsqr(radius - dist);
774
float Gd = exp(DistFromWaveSqr * inexp) * outexp;
776
m.vel += dir * (Gd * magdt / (dist * (distSqr + epsilon)));
782
// Follow the next particle in the list
783
void PAFollow::Execute(ParticleGroup *group)
785
float magdt = magnitude * dt;
786
float max_radiusSqr = max_radius * max_radius;
788
if(max_radiusSqr < P_MAXFLOAT)
790
for(int i = 0; i < group->p_count - 1; i++)
792
Particle &m = group->list[i];
794
// Accelerate toward the particle after me in the list.
795
pVector tohim(group->list[i+1].pos - m.pos); // tohim = p1 - p0
796
float tohimlenSqr = tohim.length2();
798
if(tohimlenSqr < max_radiusSqr)
800
// Compute force exerted between the two bodies
801
m.vel += tohim * (magdt / (sqrtf(tohimlenSqr) * (tohimlenSqr + epsilon)));
807
for(int i = 0; i < group->p_count - 1; i++)
809
Particle &m = group->list[i];
811
// Accelerate toward the particle after me in the list.
812
pVector tohim(group->list[i+1].pos - m.pos); // tohim = p1 - p0
813
float tohimlenSqr = tohim.length2();
815
// Compute force exerted between the two bodies
816
m.vel += tohim * (magdt / (sqrtf(tohimlenSqr) * (tohimlenSqr + epsilon)));
821
// Inter-particle gravitation
822
void PAGravitate::Execute(ParticleGroup *group)
824
float magdt = magnitude * dt;
825
float max_radiusSqr = max_radius * max_radius;
827
if(max_radiusSqr < P_MAXFLOAT)
829
for(int i = 0; i < group->p_count; i++)
831
Particle &m = group->list[i];
833
// Add interactions with other particles
834
for(int j = i + 1; j < group->p_count; j++)
836
Particle &mj = group->list[j];
838
pVector tohim(mj.pos - m.pos); // tohim = p1 - p0
839
float tohimlenSqr = tohim.length2();
841
if(tohimlenSqr < max_radiusSqr)
843
// Compute force exerted between the two bodies
844
pVector acc(tohim * (magdt / (sqrtf(tohimlenSqr) * (tohimlenSqr + epsilon))));
854
for(int i = 0; i < group->p_count; i++)
856
Particle &m = group->list[i];
858
// Add interactions with other particles
859
for(int j = i + 1; j < group->p_count; j++)
861
Particle &mj = group->list[j];
863
pVector tohim(mj.pos - m.pos); // tohim = p1 - p0
864
float tohimlenSqr = tohim.length2();
866
// Compute force exerted between the two bodies
867
pVector acc(tohim * (magdt / (sqrtf(tohimlenSqr) * (tohimlenSqr + epsilon))));
876
// Acceleration in a constant direction
877
void PAGravity::Execute(ParticleGroup *group)
879
pVector ddir(direction * dt);
881
for(int i = 0; i < group->p_count; i++)
883
// Step velocity with acceleration
884
group->list[i].vel += ddir;
888
// Accelerate particles along a line
889
void PAJet::Execute(ParticleGroup *group)
891
float magdt = magnitude * dt;
892
float max_radiusSqr = max_radius * max_radius;
894
if(max_radiusSqr < P_MAXFLOAT)
896
for(int i = 0; i < group->p_count; i++)
898
Particle &m = group->list[i];
900
// Figure direction to particle.
901
pVector dir(m.pos - center);
903
// Distance to jet (force drops as 1/r^2)
904
// Soften by epsilon to avoid tight encounters to infinity
905
float rSqr = dir.length2();
907
if(rSqr < max_radiusSqr)
912
// Step velocity with acceleration
913
m.vel += accel * (magdt / (rSqr + epsilon));
919
for(int i = 0; i < group->p_count; i++)
921
Particle &m = group->list[i];
923
// Figure direction to particle.
924
pVector dir(m.pos - center);
926
// Distance to jet (force drops as 1/r^2)
927
// Soften by epsilon to avoid tight encounters to infinity
928
float rSqr = dir.length2();
933
// Step velocity with acceleration
934
m.vel += accel * (magdt / (rSqr + epsilon));
939
// Get rid of older particles
940
void PAKillOld::Execute(ParticleGroup *group)
942
// Must traverse list in reverse order so Remove will work
943
for(int i = group->p_count-1; i >= 0; i--)
945
Particle &m = group->list[i];
947
if(!((m.age < age_limit) ^ kill_less_than))
952
// Match velocity to near neighbors
953
void PAMatchVelocity::Execute(ParticleGroup *group)
955
float magdt = magnitude * dt;
956
float max_radiusSqr = max_radius * max_radius;
958
if(max_radiusSqr < P_MAXFLOAT)
960
for(int i = 0; i < group->p_count; i++)
962
Particle &m = group->list[i];
964
// Add interactions with other particles
965
for(int j = i + 1; j < group->p_count; j++)
967
Particle &mj = group->list[j];
969
pVector tohim(mj.pos - m.pos); // tohim = p1 - p0
970
float tohimlenSqr = tohim.length2();
972
if(tohimlenSqr < max_radiusSqr)
974
// Compute force exerted between the two bodies
975
pVector acc(mj.vel * (magdt / (tohimlenSqr + epsilon)));
985
for(int i = 0; i < group->p_count; i++)
987
Particle &m = group->list[i];
989
// Add interactions with other particles
990
for(int j = i + 1; j < group->p_count; j++)
992
Particle &mj = group->list[j];
994
pVector tohim(mj.pos - m.pos); // tohim = p1 - p0
995
float tohimlenSqr = tohim.length2();
997
// Compute force exerted between the two bodies
998
pVector acc(mj.vel * (magdt / (tohimlenSqr + epsilon)));
1007
void PAMove::Execute(ParticleGroup *group)
1009
// Step particle positions forward by dt, and age the particles.
1010
for(int i = 0; i < group->p_count; i++)
1012
Particle &m = group->list[i];
1015
m.pos += m.vel * dt;
1019
// Accelerate particles towards a line
1020
void PAOrbitLine::Execute(ParticleGroup *group)
1022
float magdt = magnitude * dt;
1023
float max_radiusSqr = max_radius * max_radius;
1025
if(max_radiusSqr < P_MAXFLOAT)
1027
for(int i = 0; i < group->p_count; i++)
1029
Particle &m = group->list[i];
1031
// Figure direction to particle from base of line.
1032
pVector f(m.pos - p);
1034
pVector w(axis * (f * axis));
1036
// Direction from particle to nearest point on line.
1037
pVector into = w - f;
1039
// Distance to line (force drops as 1/r^2, normalize by 1/r)
1040
// Soften by epsilon to avoid tight encounters to infinity
1041
float rSqr = into.length2();
1043
if(rSqr < max_radiusSqr)
1044
// Step velocity with acceleration
1045
m.vel += into * (magdt / (sqrtf(rSqr) + (rSqr + epsilon)));
1050
// Removed because it causes pipeline stalls.
1051
for(int i = 0; i < group->p_count; i++)
1053
Particle &m = group->list[i];
1055
// Figure direction to particle from base of line.
1056
pVector f(m.pos - p);
1058
pVector w(axis * (f * axis));
1060
// Direction from particle to nearest point on line.
1061
pVector into = w - f;
1063
// Distance to line (force drops as 1/r^2, normalize by 1/r)
1064
// Soften by epsilon to avoid tight encounters to infinity
1065
float rSqr = into.length2();
1067
// Step velocity with acceleration
1068
m.vel += into * (magdt / (sqrtf(rSqr) + (rSqr + epsilon)));
1073
// Accelerate particles towards a point
1074
void PAOrbitPoint::Execute(ParticleGroup *group)
1076
float magdt = magnitude * dt;
1077
float max_radiusSqr = max_radius * max_radius;
1079
if(max_radiusSqr < P_MAXFLOAT)
1081
for(int i = 0; i < group->p_count; i++)
1083
Particle &m = group->list[i];
1085
// Figure direction to particle.
1086
pVector dir(center - m.pos);
1088
// Distance to gravity well (force drops as 1/r^2, normalize by 1/r)
1089
// Soften by epsilon to avoid tight encounters to infinity
1090
float rSqr = dir.length2();
1092
// Step velocity with acceleration
1093
if(rSqr < max_radiusSqr)
1094
m.vel += dir * (magdt / (sqrtf(rSqr) + (rSqr + epsilon)));
1099
// Avoids pipeline stalls.
1100
for(int i = 0; i < group->p_count; i++)
1102
Particle &m = group->list[i];
1104
// Figure direction to particle.
1105
pVector dir(center - m.pos);
1107
// Distance to gravity well (force drops as 1/r^2, normalize by 1/r)
1108
// Soften by epsilon to avoid tight encounters to infinity
1109
float rSqr = dir.length2();
1111
// Step velocity with acceleration
1112
m.vel += dir * (magdt / (sqrtf(rSqr) + (rSqr + epsilon)));
1117
// Accelerate in random direction each time step
1118
void PARandomAccel::Execute(ParticleGroup *group)
1120
for(int i = 0; i < group->p_count; i++)
1122
Particle &m = group->list[i];
1124
pVector acceleration;
1125
gen_acc.Generate(acceleration);
1127
// dt will affect this by making a higher probability of
1128
// being near the original velocity after unit time. Smaller
1129
// dt approach a normal distribution instead of a square wave.
1130
m.vel += acceleration * dt;
1134
// Immediately displace position randomly
1135
void PARandomDisplace::Execute(ParticleGroup *group)
1137
for(int i = 0; i < group->p_count; i++)
1139
Particle &m = group->list[i];
1141
pVector displacement;
1142
gen_disp.Generate(displacement);
1144
// dt will affect this by making a higher probability of
1145
// being near the original position after unit time. Smaller
1146
// dt approach a normal distribution instead of a square wave.
1147
m.pos += displacement * dt;
1151
// Immediately assign a random velocity
1152
void PARandomVelocity::Execute(ParticleGroup *group)
1154
for(int i = 0; i < group->p_count; i++)
1156
Particle &m = group->list[i];
1159
gen_vel.Generate(velocity);
1161
// Shouldn't multiply by dt because velocities are
1162
// invariant of dt. How should dt affect this?
1168
// Produce coefficients of a velocity function v(t)=at^2 + bt + c
1169
// satisfying initial x(0)=x0,v(0)=v0 and desired x(t)=xf,v(t)=vf,
1170
// where x = x(0) + integrate(v(T),0,t)
1171
static inline void _pconstrain(float x0, float v0, float xf, float vf,
1172
float t, float *a, float *b, float *c)
1175
*b = 2 * (-t*vf - 2*t*v0 + 3*xf - 3*x0) / (t * t);
1176
*a = 3 * (t*vf + t*v0 - 2*xf + 2*x0) / (t * t * t);
1180
// Over time, restore particles to initial positions
1181
// Put all particles on the surface of a statue, explode the statue,
1182
// and then suck the particles back to the original position. Cool!
1183
void PARestore::Execute(ParticleGroup *group)
1187
for(int i = 0; i < group->p_count; i++)
1189
Particle &m = group->list[i];
1191
// Already constrained, keep it there.
1193
m.vel = pVector(0,0,0);
1198
float t = time_left;
1199
float dtSqr = dt * dt;
1200
float tSqrInv2dt = dt * 2.0f / (t * t);
1201
float tCubInv3dtSqr = dtSqr * 3.0f / (t * t * t);
1203
for(int i = 0; i < group->p_count; i++)
1206
Particle &m = group->list[i];
1208
// Solve for a desired-behavior velocity function in each axis
1209
// _pconstrain(m.pos.x, m.vel.x, m.posB.x, 0., timeLeft, &a, &b, &c);
1211
// Figure new velocity at next timestep
1212
// m.vel.x = a * dtSqr + b * dt + c;
1214
float b = (-2*t*m.vel.x + 3*m.posB.x - 3*m.pos.x) * tSqrInv2dt;
1215
float a = (t*m.vel.x - m.posB.x - m.posB.x + m.pos.x + m.pos.x) * tCubInv3dtSqr;
1217
// Figure new velocity at next timestep
1220
b = (-2*t*m.vel.y + 3*m.posB.y - 3*m.pos.y) * tSqrInv2dt;
1221
a = (t*m.vel.y - m.posB.y - m.posB.y + m.pos.y + m.pos.y) * tCubInv3dtSqr;
1223
// Figure new velocity at next timestep
1226
b = (-2*t*m.vel.z + 3*m.posB.z - 3*m.pos.z) * tSqrInv2dt;
1227
a = (t*m.vel.z - m.posB.z - m.posB.z + m.pos.z + m.pos.z) * tCubInv3dtSqr;
1229
// Figure new velocity at next timestep
1232
Particle &m = group->list[i];
1234
// XXX Optimize this.
1235
// Solve for a desired-behavior velocity function in each axis
1236
float a, b, c; // Coefficients of velocity function needed
1238
_pconstrain(m.pos.x, m.vel.x, m.posB.x, 0.,
1239
timeLeft, &a, &b, &c);
1241
// Figure new velocity at next timestep
1242
m.vel.x = a * dtSqr + b * dt + c;
1244
_pconstrain(m.pos.y, m.vel.y, m.posB.y, 0.,
1245
timeLeft, &a, &b, &c);
1247
// Figure new velocity at next timestep
1248
m.vel.y = a * dtSqr + b * dt + c;
1250
_pconstrain(m.pos.z, m.vel.z, m.posB.z, 0.,
1251
timeLeft, &a, &b, &c);
1253
// Figure new velocity at next timestep
1254
m.vel.z = a * dtSqr + b * dt + c;
1263
// Kill particles with positions on wrong side of the specified domain
1264
void PASink::Execute(ParticleGroup *group)
1266
// Must traverse list in reverse order so Remove will work
1267
for(int i = group->p_count-1; i >= 0; i--)
1269
Particle &m = group->list[i];
1271
// Remove if inside/outside flag matches object's flag
1272
if(!(position.Within(m.pos) ^ kill_inside))
1277
// Kill particles with velocities on wrong side of the specified domain
1278
void PASinkVelocity::Execute(ParticleGroup *group)
1280
// Must traverse list in reverse order so Remove will work
1281
for(int i = group->p_count-1; i >= 0; i--)
1283
Particle &m = group->list[i];
1285
// Remove if inside/outside flag matches object's flag
1286
if(!(velocity.Within(m.vel) ^ kill_inside))
1291
// Randomly add particles to the system
1292
void PASource::Execute(ParticleGroup *group)
1294
int rate = int(floor(particle_rate * dt));
1296
// Dither the fraction particle in time.
1297
if(drand48() < particle_rate * dt - float(rate))
1300
// Don't emit more than it can hold.
1301
if(group->p_count + rate > group->max_particles)
1302
rate = group->max_particles - group->p_count;
1304
pVector pos, posB, vel, col, siz;
1308
for(int i = 0; i < rate; i++)
1310
position.Generate(pos);
1312
velocity.Generate(vel);
1313
color.Generate(col);
1314
float ag = age + NRand(age_sigma);
1316
group->Add(pos, pos, siz, vel, col, alpha, ag);
1321
for(int i = 0; i < rate; i++)
1323
position.Generate(pos);
1324
positionB.Generate(posB);
1326
velocity.Generate(vel);
1327
color.Generate(col);
1328
float ag = age + NRand(age_sigma);
1330
group->Add(pos, posB, siz, vel, col, alpha, ag);
1335
void PASpeedLimit::Execute(ParticleGroup *group)
1337
float min_sqr = min_speed*min_speed;
1338
float max_sqr = max_speed*max_speed;
1340
for(int i = 0; i < group->p_count; i++)
1342
Particle &m = group->list[i];
1343
float sSqr = m.vel.length2();
1344
if(sSqr<min_sqr && sSqr)
1346
float s = sqrtf(sSqr);
1347
m.vel *= (min_speed/s);
1349
else if(sSqr>max_sqr)
1351
float s = sqrtf(sSqr);
1352
m.vel *= (max_speed/s);
1357
// Change color of all particles toward the specified color
1358
void PATargetColor::Execute(ParticleGroup *group)
1360
float scaleFac = scale * dt;
1362
for(int i = 0; i < group->p_count; i++)
1364
Particle &m = group->list[i];
1365
m.color += (color - m.color) * scaleFac;
1366
m.alpha += (alpha - m.alpha) * scaleFac;
1370
// Change sizes of all particles toward the specified size
1371
void PATargetSize::Execute(ParticleGroup *group)
1373
float scaleFac_x = scale.x * dt;
1374
float scaleFac_y = scale.y * dt;
1375
float scaleFac_z = scale.z * dt;
1377
for(int i = 0; i < group->p_count; i++)
1379
Particle &m = group->list[i];
1380
pVector dif(size - m.size);
1381
dif.x *= scaleFac_x;
1382
dif.y *= scaleFac_y;
1383
dif.z *= scaleFac_z;
1388
// Change velocity of all particles toward the specified velocity
1389
void PATargetVelocity::Execute(ParticleGroup *group)
1391
float scaleFac = scale * dt;
1393
for(int i = 0; i < group->p_count; i++)
1395
Particle &m = group->list[i];
1396
m.vel += (velocity - m.vel) * scaleFac;
1400
// Immediately displace position using vortex
1401
// Vortex tip at center, around axis, with magnitude
1402
// and tightness exponent
1403
void PAVortex::Execute(ParticleGroup *group)
1405
float magdt = magnitude * dt;
1406
float max_radiusSqr = max_radius * max_radius;
1408
if(max_radiusSqr < P_MAXFLOAT)
1410
for(int i = 0; i < group->p_count; i++)
1412
Particle &m = group->list[i];
1414
// Vector from tip of vortex
1415
pVector offset(m.pos - center);
1417
// Compute distance from particle to tip of vortex.
1418
float rSqr = offset.length2();
1420
// Don't do anything to particle if too close or too far.
1421
if(rSqr > max_radiusSqr)
1424
float r = sqrtf(rSqr);
1426
// Compute normalized offset vector.
1427
pVector offnorm(offset / r);
1429
// Construct orthogonal vector frame in which to rotate
1430
// transformed point around origin
1431
float axisProj = offnorm * axis; // offnorm . axis
1433
// Components of offset perpendicular and parallel to axis
1434
pVector w(axis * axisProj); // parallel component
1435
pVector u(offnorm - w); // perpendicular component
1437
// Perpendicular component completing frame:
1438
pVector v(axis ^ u);
1440
// Figure amount of rotation
1441
// Resultant is (cos theta) u + (sin theta) v
1442
float theta = magdt / (rSqr + epsilon);
1443
float s = sinf(theta);
1444
float c = cosf(theta);
1446
offset = (u * c + v * s + w) * r;
1448
// Translate back to object space
1449
m.pos = offset + center;
1454
for(int i = 0; i < group->p_count; i++)
1456
Particle &m = group->list[i];
1458
// Vector from tip of vortex
1459
pVector offset(m.pos - center);
1461
// Compute distance from particle to tip of vortex.
1462
float rSqr = offset.length2();
1464
float r = sqrtf(rSqr);
1466
// Compute normalized offset vector.
1467
pVector offnorm(offset / r);
1469
// Construct orthogonal vector frame in which to rotate
1470
// transformed point around origin
1471
float axisProj = offnorm * axis; // offnorm . axis
1473
// Components of offset perpendicular and parallel to axis
1474
pVector w(axis * axisProj); // parallel component
1475
pVector u(offnorm - w); // perpendicular component
1477
// Perpendicular component completing frame:
1478
pVector v(axis ^ u);
1480
// Figure amount of rotation
1481
// Resultant is (cos theta) u + (sin theta) v
1482
float theta = magdt / (rSqr + epsilon);
1483
float s = sinf(theta);
1484
float c = cosf(theta);
1486
offset = (u * c + v * s + w) * r;
1488
// Translate back to object space
1489
m.pos = offset + center;
1494
////////////////////////////////////////////////////////////////////////////////
1495
// Stuff for the pDomain.
1497
pDomain::pDomain(PDomainEnum dtype, float a0, float a1,
1498
float a2, float a3, float a4, float a5,
1499
float a6, float a7, float a8)
1505
p1 = pVector(a0, a1, a2);
1509
p1 = pVector(a0, a1, a2);
1510
pVector tmp(a3, a4, a5);
1511
// p2 is vector from p1 to other endpoint.
1516
// p1 is the min corner. p2 is the max corner.
1519
p1.x = a0; p2.x = a3;
1523
p1.x = a3; p2.x = a0;
1527
p1.y = a1; p2.y = a4;
1531
p1.y = a4; p2.y = a1;
1535
p1.z = a2; p2.z = a5;
1539
p1.z = a5; p2.z = a2;
1544
p1 = pVector(a0, a1, a2);
1545
pVector tp2 = pVector(a3, a4, a5);
1546
pVector tp3 = pVector(a6, a7, a8);
1551
// The rest of this is needed for bouncing.
1552
radius1Sqr = u.length();
1553
pVector tu = u / radius1Sqr;
1554
radius2Sqr = v.length();
1555
pVector tv = v / radius2Sqr;
1557
p2 = tu ^ tv; // This is the non-unit normal.
1558
p2.normalize(); // Must normalize it.
1560
// radius1 stores the d of the plane eqn.
1561
radius1 = -(p1 * p2);
1566
p1 = pVector(a0, a1, a2);
1567
u = pVector(a3, a4, a5);
1568
v = pVector(a6, a7, a8);
1570
// The rest of this is needed for bouncing.
1571
radius1Sqr = u.length();
1572
pVector tu = u / radius1Sqr;
1573
radius2Sqr = v.length();
1574
pVector tv = v / radius2Sqr;
1576
p2 = tu ^ tv; // This is the non-unit normal.
1577
p2.normalize(); // Must normalize it.
1579
// radius1 stores the d of the plane eqn.
1580
radius1 = -(p1 * p2);
1585
p1 = pVector(a0, a1, a2);
1586
p2 = pVector(a3, a4, a5);
1587
p2.normalize(); // Must normalize it.
1589
// radius1 stores the d of the plane eqn.
1590
radius1 = -(p1 * p2);
1594
p1 = pVector(a0, a1, a2);
1597
radius1 = a3; radius2 = a4;
1601
radius1 = a4; radius2 = a3;
1603
radius1Sqr = radius1 * radius1;
1604
radius2Sqr = radius2 * radius2;
1609
// p2 is a vector from p1 to the other end of cylinder.
1610
// p1 is apex of cone.
1612
p1 = pVector(a0, a1, a2);
1613
pVector tmp(a3, a4, a5);
1618
radius1 = a6; radius2 = a7;
1622
radius1 = a7; radius2 = a6;
1624
radius1Sqr = fsqr(radius1);
1626
// Given an arbitrary nonzero vector n, make two orthonormal
1627
// vectors u and v forming a frame [u,v,n.normalize()].
1629
float p2l2 = n.length2(); // Optimize this.
1632
// radius2Sqr stores 1 / (p2.p2)
1633
// XXX Used to have an actual if.
1634
radius2Sqr = p2l2 ? 1.0f / p2l2 : 0.0f;
1636
// Find a vector orthogonal to n.
1637
pVector basis(1.0f, 0.0f, 0.0f);
1638
if(fabs(basis * n) > 0.999)
1639
basis = pVector(0.0f, 1.0f, 0.0f);
1641
// Project away N component, normalize and cross to get
1642
// second orthonormal vector.
1643
u = basis - n * (basis * n);
1650
p1 = pVector(a0, a1, a2);
1652
float tmp = 1./radius1;
1653
radius2Sqr = -0.5f*fsqr(tmp);
1654
radius2 = ONEOVERSQRT2PI * tmp;
1659
p1 = pVector(a0, a1, a2); // Center point
1660
p2 = pVector(a3, a4, a5); // Normal (not used in Within and Generate)
1665
radius1 = a6; radius2 = a7;
1669
radius1 = a7; radius2 = a6;
1672
// Find a vector orthogonal to n.
1673
pVector basis(1.0f, 0.0f, 0.0f);
1674
if(fabs(basis * p2) > 0.999)
1675
basis = pVector(0.0f, 1.0f, 0.0f);
1677
// Project away N component, normalize and cross to get
1678
// second orthonormal vector.
1679
u = basis - p2 * (basis * p2);
1682
radius1Sqr = -(p1 * p2); // D of the plane eqn.
1688
// Determines if pos is inside the domain
1689
bool pDomain::Within(const pVector &pos) const
1694
return !((pos.x < p1.x) || (pos.x > p2.x) ||
1695
(pos.y < p1.y) || (pos.y > p2.y) ||
1696
(pos.z < p1.z) || (pos.z > p2.z));
1698
// Distance from plane = n * p + d
1699
// Inside is the positive half-space.
1700
return pos * p2 >= -radius1;
1703
pVector rvec(pos - p1);
1704
float rSqr = rvec.length2();
1705
return rSqr <= radius1Sqr && rSqr >= radius2Sqr;
1710
// This is painful and slow. Might be better to do quick
1711
// accept/reject tests.
1712
// Let p2 = vector from base to tip of the cylinder
1713
// x = vector from base to test point
1715
// dist = ------ = projected distance of x along the axis
1716
// p2. p2 ranging from 0 (base) to 1 (tip)
1718
// rad = x - dist * p2 = projected vector of x along the base
1719
// p1 is the apex of the cone.
1721
pVector x(pos - p1);
1723
// Check axial distance
1724
// radius2Sqr stores 1 / (p2.p2)
1725
float dist = (p2 * x) * radius2Sqr;
1726
if(dist < 0.0f || dist > 1.0f)
1729
// Check radial distance; scale radius along axis for cones
1730
pVector xrad = x - p2 * dist; // Radial component of x
1731
float rSqr = xrad.length2();
1734
return (rSqr <= fsqr(dist * radius1) &&
1735
rSqr >= fsqr(dist * radius2));
1737
return (rSqr <= radius1Sqr && rSqr >= fsqr(radius2));
1741
pVector x(pos - p1);
1742
// return exp(-0.5 * xSq * Sqr(oneOverSigma)) * ONEOVERSQRT2PI * oneOverSigma;
1743
float Gx = expf(x.length2() * radius2Sqr) * radius2;
1744
return (drand48() < Gx);
1752
return false; // XXX Is there something better?
1756
// Generate a random point uniformly distrbuted within the domain
1757
void pDomain::Generate(pVector &pos) const
1765
pos = p1 + p2 * drand48();
1768
// Scale and translate [0,1] random to fit box
1769
pos.x = p1.x + (p2.x - p1.x) * drand48();
1770
pos.y = p1.y + (p2.y - p1.y) * drand48();
1771
pos.z = p1.z + (p2.z - p1.z) * drand48();
1775
float r1 = drand48();
1776
float r2 = drand48();
1778
pos = p1 + u * r1 + v * r2;
1780
pos = p1 + u * (1.0f-r1) + v * (1.0f-r2);
1784
pos = p1 + u * drand48() + v * drand48();
1786
case PDPlane: // How do I sensibly make a point on an infinite plane?
1790
// Place on [-1..1] sphere
1791
pos = RandVec() - vHalf;
1794
// Scale unit sphere pos by [0..r] and translate
1795
// (should distribute as r^2 law)
1796
if(radius1 == radius2)
1797
pos = p1 + pos * radius1;
1799
pos = p1 + pos * (radius2 + drand48() * (radius1 - radius2));
1804
// For a cone, p2 is the apex of the cone.
1805
float dist = drand48(); // Distance between base and tip
1806
float theta = drand48() * 2.0f * float(M_PI); // Angle around axis
1807
// Distance from axis
1808
float r = radius2 + drand48() * (radius1 - radius2);
1810
float x = r * cosf(theta); // Weighting of each frame vector
1811
float y = r * sinf(theta);
1813
// Scale radius along axis for cones
1820
pos = p1 + p2 * dist + u * x + v * y;
1824
pos.x = p1.x + NRand(radius1);
1825
pos.y = p1.y + NRand(radius1);
1826
pos.z = p1.z + NRand(radius1);
1831
float theta = drand48() * 2.0f * float(M_PI); // Angle around normal
1832
// Distance from center
1833
float r = radius2 + drand48() * (radius1 - radius2);
1835
float x = r * cosf(theta); // Weighting of each frame vector
1836
float y = r * sinf(theta);
1838
pos = p1 + u * x + v * y;
1842
pos = pVector(0,0,0);
added
removed
src/thirdparty/particles/actions.cpp
