source/blender/blenkernel/intern/particle_system.c

   1 /* particle_system.c
   2  *
   3  *
   4  * $Id$
   5  *
   6  * ***** BEGIN GPL LICENSE BLOCK *****
   7  *
   8  * This program is free software; you can redistribute it and/or
   9  * modify it under the terms of the GNU General Public License
  10  * as published by the Free Software Foundation; either version 2
  11  * of the License, or (at your option) any later version.
  12  *
  13  * This program is distributed in the hope that it will be useful,
  14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  * GNU General Public License for more details.
  17  *
  18  * You should have received a copy of the GNU General Public License
  19  * along with this program; if not, write to the Free Software Foundation,
  20  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  21  *
  22  * The Original Code is Copyright (C) 2007 by Janne Karhu.
  23  * All rights reserved.
  24  *
  25  * The Original Code is: all of this file.
  26  *
  27  * Contributor(s): Raul Fernandez Hernandez (Farsthary), Stephen Swhitehorn.
  28  *
  29  * ***** END GPL LICENSE BLOCK *****
  30  */
  31
  32 #include <stddef.h>
  33 #include "BLI_storage.h" /* _LARGEFILE_SOURCE */
  34
  35 #include <stdlib.h>
  36 #include <math.h>
  37 #include <string.h>
  38
  39 #include "MEM_guardedalloc.h"
  40
  41 #include "DNA_anim_types.h"
  42 #include "DNA_boid_types.h"
  43 #include "DNA_particle_types.h"
  44 #include "DNA_mesh_types.h"
  45 #include "DNA_meshdata_types.h"
  46 #include "DNA_modifier_types.h"
  47 #include "DNA_object_force.h"
  48 #include "DNA_object_types.h"
  49 #include "DNA_material_types.h"
  50 #include "DNA_curve_types.h"
  51 #include "DNA_group_types.h"
  52 #include "DNA_scene_types.h"
  53 #include "DNA_texture_types.h"
  54 #include "DNA_ipo_types.h" // XXX old animation system stuff... to be removed!
  55 #include "DNA_listBase.h"
  56
  57 #include "BLI_edgehash.h"
  58 #include "BLI_rand.h"
  59 #include "BLI_jitter.h"
  60 #include "BLI_math.h"
  61 #include "BLI_blenlib.h"
  62 #include "BLI_kdtree.h"
  63 #include "BLI_kdopbvh.h"
  64 #include "BLI_listbase.h"
  65 #include "BLI_threads.h"
  66 #include "BLI_storage.h" /* For _LARGEFILE64_SOURCE;  zlib needs this on some systems */
  67 #include "BLI_utildefines.h"
  68
  69 #include "BKE_main.h"
  70 #include "BKE_animsys.h"
  71 #include "BKE_boids.h"
  72 #include "BKE_cdderivedmesh.h"
  73 #include "BKE_collision.h"
  74 #include "BKE_displist.h"
  75 #include "BKE_effect.h"
  76 #include "BKE_particle.h"
  77 #include "BKE_global.h"
  78
  79 #include "BKE_DerivedMesh.h"
  80 #include "BKE_object.h"
  81 #include "BKE_material.h"
  82 #include "BKE_cloth.h"
  83 #include "BKE_depsgraph.h"
  84 #include "BKE_lattice.h"
  85 #include "BKE_pointcache.h"
  86 #include "BKE_mesh.h"
  87 #include "BKE_modifier.h"
  88 #include "BKE_scene.h"
  89 #include "BKE_bvhutils.h"
  90
  91 #include "PIL_time.h"
  92
  93 #include "RE_shader_ext.h"
  94
  95 /* fluid sim particle import */
  96 #ifndef DISABLE_ELBEEM
  97 #include "DNA_object_fluidsim.h"
  98 #include "LBM_fluidsim.h"
  99 #include <zlib.h>
 100 #include <string.h>
 101
 102 #ifdef WIN32
 103 #ifndef snprintf
 104 #define snprintf _snprintf
 105 #endif
 106 #endif
 107
 108 #endif // DISABLE_ELBEEM
 109
 110 /************************************************/
 111 /*                      Reacting to system events                       */
 112 /************************************************/
 113
 114 static int particles_are_dynamic(ParticleSystem *psys) {
 115         if(psys->pointcache->flag & PTCACHE_BAKED)
 116                 return 0;
 117
 118         if(psys->part->type == PART_HAIR)
 119                 return psys->flag & PSYS_HAIR_DYNAMICS;
 120         else
 121                 return ELEM3(psys->part->phystype, PART_PHYS_NEWTON, PART_PHYS_BOIDS, PART_PHYS_FLUID);
 122 }
 123 int psys_get_current_display_percentage(ParticleSystem *psys)
 124 {
 125         ParticleSettings *part=psys->part;
 126
 127         if((psys->renderdata && !particles_are_dynamic(psys)) /* non-dynamic particles can be rendered fully */
 128                 || (part->child_nbr && part->childtype) /* display percentage applies to children */
 129                 || (psys->pointcache->flag & PTCACHE_BAKING)) /* baking is always done with full amount */
 130                 return 100;
 131
 132         return psys->part->disp;
 133 }
 134
 135 static int tot_particles(ParticleSystem *psys, PTCacheID *pid)
 136 {
 137         if(pid && psys->pointcache->flag & PTCACHE_EXTERNAL)
 138                 return pid->cache->totpoint;
 139         else if(psys->part->distr == PART_DISTR_GRID && psys->part->from != PART_FROM_VERT)
 140                 return psys->part->grid_res * psys->part->grid_res * psys->part->grid_res;
 141         else
 142                 return psys->part->totpart;
 143 }
 144
 145 void psys_reset(ParticleSystem *psys, int mode)
 146 {
 147         PARTICLE_P;
 148
 149         if(ELEM(mode, PSYS_RESET_ALL, PSYS_RESET_DEPSGRAPH)) {
 150                 if(mode == PSYS_RESET_ALL || !(psys->flag & PSYS_EDITED)) {
 151                         /* don't free if not absolutely necessary */
 152                         if(psys->totpart != tot_particles(psys, NULL)) {
 153                                 psys_free_particles(psys);
 154                                 psys->totpart= 0;
 155                         }
 156
 157                         psys->totkeyed= 0;
 158                         psys->flag &= ~(PSYS_HAIR_DONE|PSYS_KEYED);
 159
 160                         if(psys->edit && psys->free_edit) {
 161                                 psys->free_edit(psys->edit);
 162                                 psys->edit = NULL;
 163                                 psys->free_edit = NULL;
 164                         }
 165                 }
 166         }
 167         else if(mode == PSYS_RESET_CACHE_MISS) {
 168                 /* set all particles to be skipped */
 169                 LOOP_PARTICLES
 170                         pa->flag |= PARS_NO_DISP;
 171         }
 172
 173         /* reset children */
 174         if(psys->child) {
 175                 MEM_freeN(psys->child);
 176                 psys->child= 0;
 177         }
 178
 179         psys->totchild= 0;
 180
 181         /* reset path cache */
 182         psys_free_path_cache(psys, psys->edit);
 183
 184         /* reset point cache */
 185         BKE_ptcache_invalidate(psys->pointcache);
 186
 187         if(psys->fluid_springs) {
 188                 MEM_freeN(psys->fluid_springs);
 189                 psys->fluid_springs = NULL;
 190         }
 191
 192         psys->tot_fluidsprings = psys->alloc_fluidsprings = 0;
 193 }
 194
 195 static void realloc_particles(ParticleSimulationData *sim, int new_totpart)
 196 {
 197         ParticleSystem *psys = sim->psys;
 198         ParticleSettings *part = psys->part;
 199         ParticleData *newpars = NULL;
 200         BoidParticle *newboids = NULL;
 201         PARTICLE_P;
 202         int totpart, totsaved = 0;
 203
 204         if(new_totpart<0) {
 205                 if(part->distr==PART_DISTR_GRID  && part->from != PART_FROM_VERT) {
 206                         totpart= part->grid_res;
 207                         totpart*=totpart*totpart;
 208                 }
 209                 else
 210                         totpart=part->totpart;
 211         }
 212         else
 213                 totpart=new_totpart;
 214
 215         if(totpart != psys->totpart) {
 216                 if(psys->edit && psys->free_edit) {
 217                         psys->free_edit(psys->edit);
 218                         psys->edit = NULL;
 219                         psys->free_edit = NULL;
 220                 }
 221
 222                 if(totpart) {
 223                         newpars= MEM_callocN(totpart*sizeof(ParticleData), "particles");
 224                         if(psys->part->phystype == PART_PHYS_BOIDS)
 225                                 newboids= MEM_callocN(totpart*sizeof(BoidParticle), "boid particles");
 226                 }
 227
 228                 if(psys->particles) {
 229                         totsaved=MIN2(psys->totpart,totpart);
 230                         /*save old pars*/
 231                         if(totsaved) {
 232                                 memcpy(newpars,psys->particles,totsaved*sizeof(ParticleData));
 233
 234                                 if(psys->particles->boid)
 235                                         memcpy(newboids, psys->particles->boid, totsaved*sizeof(BoidParticle));
 236                         }
 237
 238                         if(psys->particles->keys)
 239                                 MEM_freeN(psys->particles->keys);
 240
 241                         if(psys->particles->boid)
 242                                 MEM_freeN(psys->particles->boid);
 243
 244                         for(p=0, pa=newpars; p<totsaved; p++, pa++) {
 245                                 if(pa->keys) {
 246                                         pa->keys= NULL;
 247                                         pa->totkey= 0;
 248                                 }
 249                         }
 250
 251                         for(p=totsaved, pa=psys->particles+totsaved; p<psys->totpart; p++, pa++)
 252                                 if(pa->hair) MEM_freeN(pa->hair);
 253
 254                         MEM_freeN(psys->particles);
 255                         psys_free_pdd(psys);
 256                 }
 257
 258                 psys->particles=newpars;
 259                 psys->totpart=totpart;
 260
 261                 if(newboids) {
 262                         LOOP_PARTICLES
 263                                 pa->boid = newboids++;
 264                 }
 265         }
 266
 267         if(psys->child) {
 268                 MEM_freeN(psys->child);
 269                 psys->child=0;
 270                 psys->totchild=0;
 271         }
 272 }
 273
 274 static int get_psys_child_number(struct Scene *scene, ParticleSystem *psys)
 275 {
 276         int nbr;
 277
 278         if(!psys->part->childtype)
 279                 return 0;
 280
 281         if(psys->renderdata)
 282                 nbr= psys->part->ren_child_nbr;
 283         else
 284                 nbr= psys->part->child_nbr;
 285
 286         return get_render_child_particle_number(&scene->r, nbr);
 287 }
 288
 289 static int get_psys_tot_child(struct Scene *scene, ParticleSystem *psys)
 290 {
 291         return psys->totpart*get_psys_child_number(scene, psys);
 292 }
 293
 294 static void alloc_child_particles(ParticleSystem *psys, int tot)
 295 {
 296         if(psys->child){
 297                 /* only re-allocate if we have to */
 298                 if(psys->part->childtype && psys->totchild == tot) {
 299                         memset(psys->child, 0, tot*sizeof(ChildParticle));
 300                         return;
 301                 }
 302
 303                 MEM_freeN(psys->child);
 304                 psys->child=0;
 305                 psys->totchild=0;
 306         }
 307
 308         if(psys->part->childtype) {
 309                 psys->totchild= tot;
 310                 if(psys->totchild)
 311                         psys->child= MEM_callocN(psys->totchild*sizeof(ChildParticle), "child_particles");
 312         }
 313 }
 314
 315 void psys_calc_dmcache(Object *ob, DerivedMesh *dm, ParticleSystem *psys)
 316 {
 317         /* use for building derived mesh mapping info:
 318
 319            node: the allocated links - total derived mesh element count
 320            nodearray: the array of nodes aligned with the base mesh's elements, so
 321                                   each original elements can reference its derived elements
 322         */
 323         Mesh *me= (Mesh*)ob->data;
 324         PARTICLE_P;
 325
 326         /* CACHE LOCATIONS */
 327         if(!dm->deformedOnly) {
 328                 /* Will use later to speed up subsurf/derivedmesh */
 329                 LinkNode *node, *nodedmelem, **nodearray;
 330                 int totdmelem, totelem, i, *origindex;
 331
 332                 if(psys->part->from == PART_FROM_VERT) {
 333                         totdmelem= dm->getNumVerts(dm);
 334                         totelem= me->totvert;
 335                         origindex= dm->getVertDataArray(dm, CD_ORIGINDEX);
 336                 }
 337                 else { /* FROM_FACE/FROM_VOLUME */
 338                         totdmelem= dm->getNumFaces(dm);
 339                         totelem= me->totface;
 340                         origindex= dm->getFaceDataArray(dm, CD_ORIGINDEX);
 341                 }
 342
 343                 nodedmelem= MEM_callocN(sizeof(LinkNode)*totdmelem, "psys node elems");
 344                 nodearray= MEM_callocN(sizeof(LinkNode *)*totelem, "psys node array");
 345
 346                 for(i=0, node=nodedmelem; i<totdmelem; i++, origindex++, node++) {
 347                         node->link= SET_INT_IN_POINTER(i);
 348
 349                         if(*origindex != -1) {
 350                                 if(nodearray[*origindex]) {
 351                                         /* prepend */
 352                                         node->next = nodearray[*origindex];
 353                                         nodearray[*origindex]= node;
 354                                 }
 355                                 else
 356                                         nodearray[*origindex]= node;
 357                         }
 358                 }
 359
 360                 /* cache the verts/faces! */
 361                 LOOP_PARTICLES {
 362                         if(pa->num < 0) {
 363                                 pa->num_dmcache = -1;
 364                                 continue;
 365                         }
 366
 367                         if(psys->part->from == PART_FROM_VERT) {
 368                                 if(nodearray[pa->num])
 369                                         pa->num_dmcache= GET_INT_FROM_POINTER(nodearray[pa->num]->link);
 370                         }
 371                         else { /* FROM_FACE/FROM_VOLUME */
 372                                 /* Note that sometimes the pa->num is over the nodearray size, this is bad, maybe there is a better place to fix this,
 373                                  * but for now passing NULL is OK. every face will be searched for the particle so its slower - Campbell */
 374                                 pa->num_dmcache= psys_particle_dm_face_lookup(ob, dm, pa->num, pa->fuv, pa->num < totelem ? nodearray[pa->num] : NULL);
 375                         }
 376                 }
 377
 378                 MEM_freeN(nodearray);
 379                 MEM_freeN(nodedmelem);
 380         }
 381         else {
 382                 /* TODO PARTICLE, make the following line unnecessary, each function
 383                  * should know to use the num or num_dmcache, set the num_dmcache to
 384                  * an invalid value, just incase */
 385
 386                 LOOP_PARTICLES
 387                         pa->num_dmcache = -1;
 388         }
 389 }
 390
 391 static void distribute_particles_in_grid(DerivedMesh *dm, ParticleSystem *psys)
 392 {
 393         ParticleData *pa=0;
 394         float min[3], max[3], delta[3], d;
 395         MVert *mv, *mvert = dm->getVertDataArray(dm,0);
 396         int totvert=dm->getNumVerts(dm), from=psys->part->from;
 397         int i, j, k, p, res=psys->part->grid_res, size[3], axis;
 398
 399         mv=mvert;
 400
 401         /* find bounding box of dm */
 402         VECCOPY(min,mv->co);
 403         VECCOPY(max,mv->co);
 404         mv++;
 405
 406         for(i=1; i<totvert; i++, mv++){
 407                 min[0]=MIN2(min[0],mv->co[0]);
 408                 min[1]=MIN2(min[1],mv->co[1]);
 409                 min[2]=MIN2(min[2],mv->co[2]);
 410
 411                 max[0]=MAX2(max[0],mv->co[0]);
 412                 max[1]=MAX2(max[1],mv->co[1]);
 413                 max[2]=MAX2(max[2],mv->co[2]);
 414         }
 415
 416         VECSUB(delta,max,min);
 417
 418         /* determine major axis */
 419         axis = (delta[0]>=delta[1])?0:((delta[1]>=delta[2])?1:2);
 420
 421         d = delta[axis]/(float)res;
 422
 423         size[axis]=res;
 424         size[(axis+1)%3]=(int)ceil(delta[(axis+1)%3]/d);
 425         size[(axis+2)%3]=(int)ceil(delta[(axis+2)%3]/d);
 426
 427         /* float errors grrr.. */
 428         size[(axis+1)%3] = MIN2(size[(axis+1)%3],res);
 429         size[(axis+2)%3] = MIN2(size[(axis+2)%3],res);
 430
 431         min[0]+=d/2.0f;
 432         min[1]+=d/2.0f;
 433         min[2]+=d/2.0f;
 434
 435         for(i=0,p=0,pa=psys->particles; i<res; i++){
 436                 for(j=0; j<res; j++){
 437                         for(k=0; k<res; k++,p++,pa++){
 438                                 pa->fuv[0]=min[0]+(float)i*d;
 439                                 pa->fuv[1]=min[1]+(float)j*d;
 440                                 pa->fuv[2]=min[2]+(float)k*d;
 441                                 pa->flag |= PARS_UNEXIST;
 442                                 pa->hair_index=0; /* abused in volume calculation */
 443                         }
 444                 }
 445         }
 446
 447         /* enable particles near verts/edges/faces/inside surface */
 448         if(from==PART_FROM_VERT){
 449                 float vec[3];
 450
 451                 pa=psys->particles;
 452
 453                 min[0]-=d/2.0f;
 454                 min[1]-=d/2.0f;
 455                 min[2]-=d/2.0f;
 456
 457                 for(i=0,mv=mvert; i<totvert; i++,mv++){
 458                         sub_v3_v3v3(vec,mv->co,min);
 459                         vec[0]/=delta[0];
 460                         vec[1]/=delta[1];
 461                         vec[2]/=delta[2];
 462                         (pa     +((int)(vec[0]*(size[0]-1))*res
 463                                 +(int)(vec[1]*(size[1]-1)))*res
 464                                 +(int)(vec[2]*(size[2]-1)))->flag &= ~PARS_UNEXIST;
 465                 }
 466         }
 467         else if(ELEM(from,PART_FROM_FACE,PART_FROM_VOLUME)){
 468                 float co1[3], co2[3];
 469
 470                 MFace *mface= NULL, *mface_array;
 471                 float v1[3], v2[3], v3[3], v4[4], lambda;
 472                 int a, a1, a2, a0mul, a1mul, a2mul, totface;
 473                 int amax= from==PART_FROM_FACE ? 3 : 1;
 474
 475                 totface=dm->getNumFaces(dm);
 476                 mface_array= dm->getFaceDataArray(dm,CD_MFACE);
 477
 478                 for(a=0; a<amax; a++){
 479                         if(a==0){ a0mul=res*res; a1mul=res; a2mul=1; }
 480                         else if(a==1){ a0mul=res; a1mul=1; a2mul=res*res; }
 481                         else{ a0mul=1; a1mul=res*res; a2mul=res; }
 482
 483                         for(a1=0; a1<size[(a+1)%3]; a1++){
 484                                 for(a2=0; a2<size[(a+2)%3]; a2++){
 485                                         mface= mface_array;
 486
 487                                         pa=psys->particles + a1*a1mul + a2*a2mul;
 488                                         VECCOPY(co1,pa->fuv);
 489                                         co1[a]-=d/2.0f;
 490                                         VECCOPY(co2,co1);
 491                                         co2[a]+=delta[a] + 0.001f*d;
 492                                         co1[a]-=0.001f*d;
 493
 494                                         /* lets intersect the faces */
 495                                         for(i=0; i<totface; i++,mface++){
 496                                                 VECCOPY(v1,mvert[mface->v1].co);
 497                                                 VECCOPY(v2,mvert[mface->v2].co);
 498                                                 VECCOPY(v3,mvert[mface->v3].co);
 499
 500                                                 if(isect_axial_line_tri_v3(a,co1, co2, v2, v3, v1, &lambda)){
 501                                                         if(from==PART_FROM_FACE)
 502                                                                 (pa+(int)(lambda*size[a])*a0mul)->flag &= ~PARS_UNEXIST;
 503                                                         else /* store number of intersections */
 504                                                                 (pa+(int)(lambda*size[a])*a0mul)->hair_index++;
 505                                                 }
 506
 507                                                 if(mface->v4){
 508                                                         VECCOPY(v4,mvert[mface->v4].co);
 509
 510                                                         if(isect_axial_line_tri_v3(a,co1, co2, v4, v1, v3, &lambda)){
 511                                                                 if(from==PART_FROM_FACE)
 512                                                                         (pa+(int)(lambda*size[a])*a0mul)->flag &= ~PARS_UNEXIST;
 513                                                                 else
 514                                                                         (pa+(int)(lambda*size[a])*a0mul)->hair_index++;
 515                                                         }
 516                                                 }
 517                                         }
 518
 519                                         if(from==PART_FROM_VOLUME){
 520                                                 int in=pa->hair_index%2;
 521                                                 if(in) pa->hair_index++;
 522                                                 for(i=0; i<size[0]; i++){
 523                                                         if(in || (pa+i*a0mul)->hair_index%2)
 524                                                                 (pa+i*a0mul)->flag &= ~PARS_UNEXIST;
 525                                                         /* odd intersections == in->out / out->in */
 526                                                         /* even intersections -> in stays same */
 527                                                         in=(in + (pa+i*a0mul)->hair_index) % 2;
 528                                                 }
 529                                         }
 530                                 }
 531                         }
 532                 }
 533         }
 534
 535         if(psys->part->flag & PART_GRID_INVERT){
 536                 for(i=0; i<size[0]; i++){
 537                         for(j=0; j<size[1]; j++){
 538                                 pa=psys->particles + res*(i*res + j);
 539                                 for(k=0; k<size[2]; k++, pa++){
 540                                         pa->flag ^= PARS_UNEXIST;
 541                                 }
 542                         }
 543                 }
 544         }
 545 }
 546
 547 /* modified copy from rayshade.c */
 548 static void hammersley_create(float *out, int n, int seed, float amount)
 549 {
 550         RNG *rng;
 551         double p, t, offs[2];
 552         int k, kk;
 553
 554         rng = rng_new(31415926 + n + seed);
 555         offs[0]= rng_getDouble(rng) + amount;
 556         offs[1]= rng_getDouble(rng) + amount;
 557         rng_free(rng);
 558
 559         for (k = 0; k < n; k++) {
 560                 t = 0;
 561                 for (p = 0.5, kk = k; kk; p *= 0.5, kk >>= 1)
 562                         if (kk & 1) /* kk mod 2 = 1 */
 563                                 t += p;
 564
 565                 out[2*k + 0]= fmod((double)k/(double)n + offs[0], 1.0);
 566                 out[2*k + 1]= fmod(t + offs[1], 1.0);
 567         }
 568 }
 569
 570 /* modified copy from effect.c */
 571 static void init_mv_jit(float *jit, int num, int seed2, float amount)
 572 {
 573         RNG *rng;
 574         float *jit2, x, rad1, rad2, rad3;
 575         int i, num2;
 576
 577         if(num==0) return;
 578
 579         rad1= (float)(1.0/sqrt((float)num));
 580         rad2= (float)(1.0/((float)num));
 581         rad3= (float)sqrt((float)num)/((float)num);
 582
 583         rng = rng_new(31415926 + num + seed2);
 584         x= 0;
 585                 num2 = 2 * num;
 586         for(i=0; i<num2; i+=2) {
 587
 588                 jit[i]= x + amount*rad1*(0.5f - rng_getFloat(rng));
 589                 jit[i+1]= i/(2.0f*num) + amount*rad1*(0.5f - rng_getFloat(rng));
 590
 591                 jit[i]-= (float)floor(jit[i]);
 592                 jit[i+1]-= (float)floor(jit[i+1]);
 593
 594                 x+= rad3;
 595                 x -= (float)floor(x);
 596         }
 597
 598         jit2= MEM_mallocN(12 + 2*sizeof(float)*num, "initjit");
 599
 600         for (i=0 ; i<4 ; i++) {
 601                 BLI_jitterate1(jit, jit2, num, rad1);
 602                 BLI_jitterate1(jit, jit2, num, rad1);
 603                 BLI_jitterate2(jit, jit2, num, rad2);
 604         }
 605         MEM_freeN(jit2);
 606         rng_free(rng);
 607 }
 608
 609 static void psys_uv_to_w(float u, float v, int quad, float *w)
 610 {
 611         float vert[4][3], co[3];
 612
 613         if(!quad) {
 614                 if(u+v > 1.0f)
 615                         v= 1.0f-v;
 616                 else
 617                         u= 1.0f-u;
 618         }
 619
 620         vert[0][0]= 0.0f; vert[0][1]= 0.0f; vert[0][2]= 0.0f;
 621         vert[1][0]= 1.0f; vert[1][1]= 0.0f; vert[1][2]= 0.0f;
 622         vert[2][0]= 1.0f; vert[2][1]= 1.0f; vert[2][2]= 0.0f;
 623
 624         co[0]= u;
 625         co[1]= v;
 626         co[2]= 0.0f;
 627
 628         if(quad) {
 629                 vert[3][0]= 0.0f; vert[3][1]= 1.0f; vert[3][2]= 0.0f;
 630                 interp_weights_poly_v3( w,vert, 4, co);
 631         }
 632         else {
 633                 interp_weights_poly_v3( w,vert, 3, co);
 634                 w[3]= 0.0f;
 635         }
 636 }
 637
 638 /* Find the index in "sum" array before "value" is crossed. */
 639 static int binary_search_distribution(float *sum, int n, float value)
 640 {
 641         int mid, low=0, high=n;
 642
 643         if(value == 0.f)
 644                 return 0;
 645
 646         while(low <= high) {
 647                 mid= (low + high)/2;
 648
 649                 if(sum[mid] < value && value <= sum[mid+1])
 650                         return mid;
 651
 652                 if(sum[mid] >= value)
 653                         high= mid - 1;
 654                 else if(sum[mid] < value)
 655                         low= mid + 1;
 656                 else
 657                         return mid;
 658         }
 659
 660         return low;
 661 }
 662
 663 /* the max number if calls to rng_* funcs within psys_thread_distribute_particle
 664  * be sure to keep up to date if this changes */
 665 #define PSYS_RND_DIST_SKIP 2
 666
 667 /* note: this function must be thread safe, for from == PART_FROM_CHILD */
 668 #define ONLY_WORKING_WITH_PA_VERTS 0
 669 static void psys_thread_distribute_particle(ParticleThread *thread, ParticleData *pa, ChildParticle *cpa, int p)
 670 {
 671         ParticleThreadContext *ctx= thread->ctx;
 672         Object *ob= ctx->sim.ob;
 673         DerivedMesh *dm= ctx->dm;
 674         ParticleData *tpa;
 675 /*      ParticleSettings *part= ctx->sim.psys->part; */
 676         float *v1, *v2, *v3, *v4, nor[3], orco1[3], co1[3], co2[3], nor1[3];
 677         float cur_d, min_d, randu, randv;
 678         int from= ctx->from;
 679         int cfrom= ctx->cfrom;
 680         int distr= ctx->distr;
 681         int i, intersect, tot;
 682         int rng_skip_tot= PSYS_RND_DIST_SKIP; /* count how many rng_* calls wont need skipping */
 683
 684         if(from == PART_FROM_VERT) {
 685                 /* TODO_PARTICLE - use original index */
 686                 pa->num= ctx->index[p];
 687                 pa->fuv[0] = 1.0f;
 688                 pa->fuv[1] = pa->fuv[2] = pa->fuv[3] = 0.0;
 689
 690 #if ONLY_WORKING_WITH_PA_VERTS
 691                 if(ctx->tree){
 692                         KDTreeNearest ptn[3];
 693                         int w, maxw;
 694
 695                         psys_particle_on_dm(ctx->dm,from,pa->num,pa->num_dmcache,pa->fuv,pa->foffset,co1,0,0,0,orco1,0);
 696                         transform_mesh_orco_verts((Mesh*)ob->data, &orco1, 1, 1);
 697                         maxw = BLI_kdtree_find_n_nearest(ctx->tree,3,orco1,NULL,ptn);
 698
 699                         for(w=0; w<maxw; w++){
 700                                 pa->verts[w]=ptn->num;
 701                         }
 702                 }
 703 #endif
 704         }
 705         else if(from == PART_FROM_FACE || from == PART_FROM_VOLUME) {
 706                 MFace *mface;
 707
 708                 pa->num = i = ctx->index[p];
 709                 mface = dm->getFaceData(dm,i,CD_MFACE);
 710
 711                 switch(distr){
 712                 case PART_DISTR_JIT:
 713                         ctx->jitoff[i] = fmod(ctx->jitoff[i],(float)ctx->jitlevel);
 714                         psys_uv_to_w(ctx->jit[2*(int)ctx->jitoff[i]], ctx->jit[2*(int)ctx->jitoff[i]+1], mface->v4, pa->fuv);
 715                         ctx->jitoff[i]++;
 716                         break;
 717                 case PART_DISTR_RAND:
 718                         randu= rng_getFloat(thread->rng);
 719                         randv= rng_getFloat(thread->rng);
 720                         rng_skip_tot -= 2;
 721
 722                         psys_uv_to_w(randu, randv, mface->v4, pa->fuv);
 723                         break;
 724                 }
 725                 pa->foffset= 0.0f;
 726
 727                 /* experimental */
 728                 if(from==PART_FROM_VOLUME){
 729                         MVert *mvert=dm->getVertDataArray(dm,CD_MVERT);
 730
 731                         tot=dm->getNumFaces(dm);
 732
 733                         psys_interpolate_face(mvert,mface,0,0,pa->fuv,co1,nor,0,0,0,0);
 734
 735                         normalize_v3(nor);
 736                         mul_v3_fl(nor,-100.0);
 737
 738                         VECADD(co2,co1,nor);
 739
 740                         min_d=2.0;
 741                         intersect=0;
 742
 743                         for(i=0,mface=dm->getFaceDataArray(dm,CD_MFACE); i<tot; i++,mface++){
 744                                 if(i==pa->num) continue;
 745
 746                                 v1=mvert[mface->v1].co;
 747                                 v2=mvert[mface->v2].co;
 748                                 v3=mvert[mface->v3].co;
 749
 750                                 if(isect_line_tri_v3(co1, co2, v2, v3, v1, &cur_d, 0)){
 751                                         if(cur_d<min_d){
 752                                                 min_d=cur_d;
 753                                                 pa->foffset=cur_d*50.0f; /* to the middle of volume */
 754                                                 intersect=1;
 755                                         }
 756                                 }
 757                                 if(mface->v4){
 758                                         v4=mvert[mface->v4].co;
 759
 760                                         if(isect_line_tri_v3(co1, co2, v4, v1, v3, &cur_d, 0)){
 761                                                 if(cur_d<min_d){
 762                                                         min_d=cur_d;
 763                                                         pa->foffset=cur_d*50.0f; /* to the middle of volume */
 764                                                         intersect=1;
 765                                                 }
 766                                         }
 767                                 }
 768                         }
 769                         if(intersect==0)
 770                                 pa->foffset=0.0;
 771                         else switch(distr){
 772                                 case PART_DISTR_JIT:
 773                                         pa->foffset*= ctx->jit[p%(2*ctx->jitlevel)];
 774                                         break;
 775                                 case PART_DISTR_RAND:
 776                                         pa->foffset*=BLI_frand();
 777                                         break;
 778                         }
 779                 }
 780         }
 781         else if(from == PART_FROM_PARTICLE) {
 782                 tpa=ctx->tpars+ctx->index[p];
 783                 pa->num=ctx->index[p];
 784                 pa->fuv[0]=tpa->fuv[0];
 785                 pa->fuv[1]=tpa->fuv[1];
 786                 /* abusing foffset a little for timing in near reaction */
 787                 pa->foffset=ctx->weight[ctx->index[p]];
 788                 ctx->weight[ctx->index[p]]+=ctx->maxweight;
 789         }
 790         else if(from == PART_FROM_CHILD) {
 791                 MFace *mf;
 792
 793                 if(ctx->index[p] < 0) {
 794                         cpa->num=0;
 795                         cpa->fuv[0]=cpa->fuv[1]=cpa->fuv[2]=cpa->fuv[3]=0.0f;
 796                         cpa->pa[0]=cpa->pa[1]=cpa->pa[2]=cpa->pa[3]=0;
 797                         return;
 798                 }
 799
 800                 mf= dm->getFaceData(dm, ctx->index[p], CD_MFACE);
 801
 802                 randu= rng_getFloat(thread->rng);
 803                 randv= rng_getFloat(thread->rng);
 804                 rng_skip_tot -= 2;
 805
 806                 psys_uv_to_w(randu, randv, mf->v4, cpa->fuv);
 807
 808                 cpa->num = ctx->index[p];
 809
 810                 if(ctx->tree){
 811                         KDTreeNearest ptn[10];
 812                         int w,maxw;//, do_seams;
 813                         float maxd,mind,/*dd,*/totw=0.0;
 814                         int parent[10];
 815                         float pweight[10];
 816
 817                         psys_particle_on_dm(dm,cfrom,cpa->num,DMCACHE_ISCHILD,cpa->fuv,cpa->foffset,co1,nor1,NULL,NULL,orco1,NULL);
 818                         transform_mesh_orco_verts((Mesh*)ob->data, &orco1, 1, 1);
 819                         maxw = BLI_kdtree_find_n_nearest(ctx->tree,4,orco1,NULL,ptn);
 820
 821                         maxd=ptn[maxw-1].dist;
 822                         mind=ptn[0].dist;
 823                         /*dd=maxd-mind;*/ /*UNUSED*/
 824
 825                         /* the weights here could be done better */
 826                         for(w=0; w<maxw; w++){
 827                                 parent[w]=ptn[w].index;
 828                                 pweight[w]=(float)pow(2.0,(double)(-6.0f*ptn[w].dist/maxd));
 829                         }
 830                         for(;w<10; w++){
 831                                 parent[w]=-1;
 832                                 pweight[w]=0.0f;
 833                         }
 834
 835                         for(w=0,i=0; w<maxw && i<4; w++){
 836                                 if(parent[w]>=0){
 837                                         cpa->pa[i]=parent[w];
 838                                         cpa->w[i]=pweight[w];
 839                                         totw+=pweight[w];
 840                                         i++;
 841                                 }
 842                         }
 843                         for(;i<4; i++){
 844                                 cpa->pa[i]=-1;
 845                                 cpa->w[i]=0.0f;
 846                         }
 847
 848                         if(totw>0.0f) for(w=0; w<4; w++)
 849                                 cpa->w[w]/=totw;
 850
 851                         cpa->parent=cpa->pa[0];
 852                 }
 853         }
 854
 855         if(rng_skip_tot > 0) /* should never be below zero */
 856                 rng_skip(thread->rng, rng_skip_tot);
 857 }
 858
 859 static void *exec_distribution(void *data)
 860 {
 861         ParticleThread *thread= (ParticleThread*)data;
 862         ParticleSystem *psys= thread->ctx->sim.psys;
 863         ParticleData *pa;
 864         ChildParticle *cpa;
 865         int p, totpart;
 866
 867         if(thread->ctx->from == PART_FROM_CHILD) {
 868                 totpart= psys->totchild;
 869                 cpa= psys->child;
 870
 871                 for(p=0; p<totpart; p++, cpa++) {
 872                         if(thread->ctx->skip) /* simplification skip */
 873                                 rng_skip(thread->rng, PSYS_RND_DIST_SKIP * thread->ctx->skip[p]);
 874
 875                         if((p+thread->num) % thread->tot == 0)
 876                                 psys_thread_distribute_particle(thread, NULL, cpa, p);
 877                         else /* thread skip */
 878                                 rng_skip(thread->rng, PSYS_RND_DIST_SKIP);
 879                 }
 880         }
 881         else {
 882                 totpart= psys->totpart;
 883                 pa= psys->particles + thread->num;
 884                 for(p=thread->num; p<totpart; p+=thread->tot, pa+=thread->tot)
 885                         psys_thread_distribute_particle(thread, pa, NULL, p);
 886         }
 887
 888         return 0;
 889 }
 890
 891 /* not thread safe, but qsort doesn't take userdata argument */
 892 static int *COMPARE_ORIG_INDEX = NULL;
 893 static int compare_orig_index(const void *p1, const void *p2)
 894 {
 895         int index1 = COMPARE_ORIG_INDEX[*(const int*)p1];
 896         int index2 = COMPARE_ORIG_INDEX[*(const int*)p2];
 897
 898         if(index1 < index2)
 899                 return -1;
 900         else if(index1 == index2) {
 901                 /* this pointer comparison appears to make qsort stable for glibc,
 902                  * and apparently on solaris too, makes the renders reproducable */
 903                 if(p1 < p2)
 904                         return -1;
 905                 else if(p1 == p2)
 906                         return 0;
 907                 else
 908                         return 1;
 909         }
 910         else
 911                 return 1;
 912 }
 913
 914 /* creates a distribution of coordinates on a DerivedMesh       */
 915 /*                                                                                                                      */
 916 /* 1. lets check from what we are emitting                                      */
 917 /* 2. now we know that we have something to emit from so        */
 918 /*        let's calculate some weights                                                  */
 919 /* 2.1 from even distribution                                                           */
 920 /* 2.2 and from vertex groups                                                           */
 921 /* 3. next we determine the indexes of emitting thing that      */
 922 /*        the particles will have                                                               */
 923 /* 4. let's do jitter if we need it                                                     */
 924 /* 5. now we're ready to set the indexes & distributions to     */
 925 /*        the particles                                                                                 */
 926 /* 6. and we're done!                                                                           */
 927
 928 /* This is to denote functionality that does not yet work with mesh - only derived mesh */
 929 static int psys_threads_init_distribution(ParticleThread *threads, Scene *scene, DerivedMesh *finaldm, int from)
 930 {
 931         ParticleThreadContext *ctx= threads[0].ctx;
 932         Object *ob= ctx->sim.ob;
 933         ParticleSystem *psys= ctx->sim.psys;
 934         Object *tob;
 935         ParticleData *pa=0, *tpars= 0;
 936         ParticleSettings *part;
 937         ParticleSystem *tpsys;
 938         ParticleSeam *seams= 0;
 939         ChildParticle *cpa=0;
 940         KDTree *tree=0;
 941         DerivedMesh *dm= NULL;
 942         float *jit= NULL;
 943         int i, seed, p=0, totthread= threads[0].tot;
 944         int /*no_distr=0,*/ cfrom=0;
 945         int tot=0, totpart, *index=0, children=0, totseam=0;
 946         //int *vertpart=0;
 947         int jitlevel= 1, distr;
 948         float *weight=0,*sum=0,*jitoff=0;
 949         float cur, maxweight=0.0, tweight, totweight, co[3], nor[3], orco[3], ornor[3];
 950
 951         if(ob==0 || psys==0 || psys->part==0)
 952                 return 0;
 953
 954         part=psys->part;
 955         totpart=psys->totpart;
 956         if(totpart==0)
 957                 return 0;
 958
 959         if (!finaldm->deformedOnly && !finaldm->getFaceDataArray(finaldm, CD_ORIGINDEX)) {
 960                 printf("Can't create particles with the current modifier stack, disable destructive modifiers\n");
 961 // XXX          error("Can't paint with the current modifier stack, disable destructive modifiers");
 962                 return 0;
 963         }
 964
 965         BLI_srandom(31415926 + psys->seed);
 966
 967         if(from==PART_FROM_CHILD){
 968                 distr=PART_DISTR_RAND;
 969                 BLI_srandom(31415926 + psys->seed + psys->child_seed);
 970
 971                 if(part->from!=PART_FROM_PARTICLE && part->childtype==PART_CHILD_FACES){
 972                         dm= finaldm;
 973                         children=1;
 974
 975                         tree=BLI_kdtree_new(totpart);
 976
 977                         for(p=0,pa=psys->particles; p<totpart; p++,pa++){
 978                                 psys_particle_on_dm(dm,part->from,pa->num,pa->num_dmcache,pa->fuv,pa->foffset,co,nor,0,0,orco,ornor);
 979                                 transform_mesh_orco_verts((Mesh*)ob->data, &orco, 1, 1);
 980                                 BLI_kdtree_insert(tree, p, orco, ornor);
 981                         }
 982
 983                         BLI_kdtree_balance(tree);
 984
 985                         totpart=get_psys_tot_child(scene, psys);
 986                         cfrom=from=PART_FROM_FACE;
 987                 }
 988                 else{
 989                         /* no need to figure out distribution */
 990                         int child_nbr= get_psys_child_number(scene, psys);
 991
 992                         totpart= get_psys_tot_child(scene, psys);
 993                         alloc_child_particles(psys, totpart);
 994                         cpa=psys->child;
 995                         for(i=0; i<child_nbr; i++){
 996                                 for(p=0; p<psys->totpart; p++,cpa++){
 997                                         float length=2.0;
 998                                         cpa->parent=p;
 999
1000                                         /* create even spherical distribution inside unit sphere */
1001                                         while(length>=1.0f){
1002                                                 cpa->fuv[0]=2.0f*BLI_frand()-1.0f;
1003                                                 cpa->fuv[1]=2.0f*BLI_frand()-1.0f;
1004                                                 cpa->fuv[2]=2.0f*BLI_frand()-1.0f;
1005                                                 length=len_v3(cpa->fuv);
1006                                         }
1007
1008                                         cpa->num=-1;
1009                                 }
1010                         }
1011                         /* dmcache must be updated for parent particles if children from faces is used */
1012                         psys_calc_dmcache(ob, finaldm, psys);
1013
1014                         return 0;
1015                 }
1016         }
1017         else{
1018                 dm= CDDM_from_mesh((Mesh*)ob->data, ob);
1019
1020                 /* special handling of grid distribution */
1021                 if(part->distr==PART_DISTR_GRID && from != PART_FROM_VERT){
1022                         distribute_particles_in_grid(dm,psys);
1023                         dm->release(dm);
1024                         return 0;
1025                 }
1026
1027                 /* we need orco for consistent distributions */
1028                 DM_add_vert_layer(dm, CD_ORCO, CD_ASSIGN, get_mesh_orco_verts(ob));
1029
1030                 distr=part->distr;
1031
1032                 if(from==PART_FROM_VERT){
1033                         MVert *mv= dm->getVertDataArray(dm, CD_MVERT);
1034                         float (*orcodata)[3]= dm->getVertDataArray(dm, CD_ORCO);
1035                         int totvert = dm->getNumVerts(dm);
1036
1037                         tree=BLI_kdtree_new(totvert);
1038
1039                         for(p=0; p<totvert; p++){
1040                                 if(orcodata) {
1041                                         VECCOPY(co,orcodata[p])
1042                                         transform_mesh_orco_verts((Mesh*)ob->data, &co, 1, 1);
1043                                 }
1044                                 else
1045                                         VECCOPY(co,mv[p].co)
1046                                 BLI_kdtree_insert(tree,p,co,NULL);
1047                         }
1048
1049                         BLI_kdtree_balance(tree);
1050                 }
1051         }
1052
1053         /* 1. */
1054         switch(from){
1055                 case PART_FROM_VERT:
1056                         tot = dm->getNumVerts(dm);
1057                         break;
1058                 case PART_FROM_VOLUME:
1059                 case PART_FROM_FACE:
1060                         tot = dm->getNumFaces(dm);
1061                         break;
1062                 case PART_FROM_PARTICLE:
1063                         if(psys->target_ob)
1064                                 tob=psys->target_ob;
1065                         else
1066                                 tob=ob;
1067
1068                         if((tpsys=BLI_findlink(&tob->particlesystem,psys->target_psys-1))){
1069                                 tpars=tpsys->particles;
1070                                 tot=tpsys->totpart;
1071                         }
1072                         break;
1073         }
1074
1075         if(tot==0){
1076                 /*no_distr=1;*/ /*UNUSED*/
1077                 if(children){
1078                         if(G.f & G_DEBUG)
1079                                 fprintf(stderr,"Particle child distribution error: Nothing to emit from!\n");
1080                         if(psys->child) {
1081                                 for(p=0,cpa=psys->child; p<totpart; p++,cpa++){
1082                                         cpa->fuv[0]=cpa->fuv[1]=cpa->fuv[2]=cpa->fuv[3]= 0.0;
1083                                         cpa->foffset= 0.0f;
1084                                         cpa->parent=0;
1085                                         cpa->pa[0]=cpa->pa[1]=cpa->pa[2]=cpa->pa[3]=0;
1086                                         cpa->num= -1;
1087                                 }
1088                         }
1089                 }
1090                 else {
1091                         if(G.f & G_DEBUG)
1092                                 fprintf(stderr,"Particle distribution error: Nothing to emit from!\n");
1093                         for(p=0,pa=psys->particles; p<totpart; p++,pa++){
1094                                 pa->fuv[0]=pa->fuv[1]=pa->fuv[2]= pa->fuv[3]= 0.0;
1095                                 pa->foffset= 0.0f;
1096                                 pa->num= -1;
1097                         }
1098                 }
1099
1100                 if(dm != finaldm) dm->release(dm);
1101                 return 0;
1102         }
1103
1104         /* 2. */
1105
1106         weight=MEM_callocN(sizeof(float)*tot, "particle_distribution_weights");
1107         index=MEM_callocN(sizeof(int)*totpart, "particle_distribution_indexes");
1108         sum=MEM_callocN(sizeof(float)*(tot+1), "particle_distribution_sum");
1109         jitoff=MEM_callocN(sizeof(float)*tot, "particle_distribution_jitoff");
1110
1111         /* 2.1 */
1112         if((part->flag&PART_EDISTR || children) && ELEM(from,PART_FROM_PARTICLE,PART_FROM_VERT)==0){
1113                 MVert *v1, *v2, *v3, *v4;
1114                 float totarea=0.0, co1[3], co2[3], co3[3], co4[3];
1115                 float (*orcodata)[3];
1116
1117                 orcodata= dm->getVertDataArray(dm, CD_ORCO);
1118
1119                 for(i=0; i<tot; i++){
1120                         MFace *mf=dm->getFaceData(dm,i,CD_MFACE);
1121
1122                         if(orcodata) {
1123                                 VECCOPY(co1, orcodata[mf->v1]);
1124                                 VECCOPY(co2, orcodata[mf->v2]);
1125                                 VECCOPY(co3, orcodata[mf->v3]);
1126                                 transform_mesh_orco_verts((Mesh*)ob->data, &co1, 1, 1);
1127                                 transform_mesh_orco_verts((Mesh*)ob->data, &co2, 1, 1);
1128                                 transform_mesh_orco_verts((Mesh*)ob->data, &co3, 1, 1);
1129                         }
1130                         else {
1131                                 v1= (MVert*)dm->getVertData(dm,mf->v1,CD_MVERT);
1132                                 v2= (MVert*)dm->getVertData(dm,mf->v2,CD_MVERT);
1133                                 v3= (MVert*)dm->getVertData(dm,mf->v3,CD_MVERT);
1134                                 VECCOPY(co1, v1->co);
1135                                 VECCOPY(co2, v2->co);
1136                                 VECCOPY(co3, v3->co);
1137                         }
1138
1139                         if (mf->v4){
1140                                 if(orcodata) {
1141                                         VECCOPY(co4, orcodata[mf->v4]);
1142                                         transform_mesh_orco_verts((Mesh*)ob->data, &co4, 1, 1);
1143                                 }
1144                                 else {
1145                                         v4= (MVert*)dm->getVertData(dm,mf->v4,CD_MVERT);
1146                                         VECCOPY(co4, v4->co);
1147                                 }
1148                                 cur= area_quad_v3(co1, co2, co3, co4);
1149                         }
1150                         else
1151                                 cur= area_tri_v3(co1, co2, co3);
1152
1153                         if(cur>maxweight)
1154                                 maxweight=cur;
1155
1156                         weight[i]= cur;
1157                         totarea+=cur;
1158                 }
1159
1160                 for(i=0; i<tot; i++)
1161                         weight[i] /= totarea;
1162
1163                 maxweight /= totarea;
1164         }
1165         else if(from==PART_FROM_PARTICLE){
1166                 float val=(float)tot/(float)totpart;
1167                 for(i=0; i<tot; i++)
1168                         weight[i]=val;
1169                 maxweight=val;
1170         }
1171         else{
1172                 float min=1.0f/(float)(MIN2(tot,totpart));
1173                 for(i=0; i<tot; i++)
1174                         weight[i]=min;
1175                 maxweight=min;
1176         }
1177
1178         /* 2.2 */
1179         if(ELEM3(from,PART_FROM_VERT,PART_FROM_FACE,PART_FROM_VOLUME)){
1180                 float *vweight= psys_cache_vgroup(dm,psys,PSYS_VG_DENSITY);
1181
1182                 if(vweight){
1183                         if(from==PART_FROM_VERT) {
1184                                 for(i=0;i<tot; i++)
1185                                         weight[i]*=vweight[i];
1186                         }
1187                         else { /* PART_FROM_FACE / PART_FROM_VOLUME */
1188                                 for(i=0;i<tot; i++){
1189                                         MFace *mf=dm->getFaceData(dm,i,CD_MFACE);
1190                                         tweight = vweight[mf->v1] + vweight[mf->v2] + vweight[mf->v3];
1191
1192                                         if(mf->v4) {
1193                                                 tweight += vweight[mf->v4];
1194                                                 tweight /= 4.0;
1195                                         }
1196                                         else {
1197                                                 tweight /= 3.0;
1198                                         }
1199
1200                                         weight[i]*=tweight;
1201                                 }
1202                         }
1203                         MEM_freeN(vweight);
1204                 }
1205         }
1206
1207         /* 3. */
1208         totweight= 0.0f;
1209         for(i=0;i<tot; i++)
1210                 totweight += weight[i];
1211
1212         if(totweight > 0.0f)
1213                 totweight= 1.0f/totweight;
1214
1215         sum[0]= 0.0f;
1216         for(i=0;i<tot; i++)
1217                 sum[i+1]= sum[i]+weight[i]*totweight;
1218
1219         if((part->flag&PART_TRAND) || (part->simplify_flag&PART_SIMPLIFY_ENABLE)) {
1220                 float pos;
1221
1222                 for(p=0; p<totpart; p++) {
1223                         /* In theory sys[tot] should be 1.0, but due to float errors this is not necessarily always true, so scale pos accordingly. */
1224                         pos= BLI_frand() * sum[tot];
1225                         index[p]= binary_search_distribution(sum, tot, pos);
1226                         index[p]= MIN2(tot-1, index[p]);
1227                         jitoff[index[p]]= pos;
1228                 }
1229         }
1230         else {
1231                 double step, pos;
1232
1233                 step= (totpart <= 1)? 0.5: 1.0/(totpart-1);
1234                 pos= 1e-16f; /* tiny offset to avoid zero weight face */
1235                 i= 0;
1236
1237                 for(p=0; p<totpart; p++, pos+=step) {
1238                         while((i < tot) && (pos > sum[i+1]))
1239                                 i++;
1240
1241                         index[p]= MIN2(tot-1, i);
1242
1243                         /* avoid zero weight face */
1244                         if(p == totpart-1 && weight[index[p]] == 0.0f)
1245                                 index[p]= index[p-1];
1246
1247                         jitoff[index[p]]= pos;
1248                 }
1249         }
1250
1251         MEM_freeN(sum);
1252
1253         /* for hair, sort by origindex, allows optimizations in rendering */
1254         /* however with virtual parents the children need to be in random order */
1255         if(part->type == PART_HAIR && !(part->childtype==PART_CHILD_FACES && part->parents!=0.0)) {
1256                 if(from != PART_FROM_PARTICLE) {
1257                         COMPARE_ORIG_INDEX = NULL;
1258
1259                         if(from == PART_FROM_VERT) {
1260                                 if(dm->numVertData)
1261                                         COMPARE_ORIG_INDEX= dm->getVertDataArray(dm, CD_ORIGINDEX);
1262                         }
1263                         else {
1264                                 if(dm->numFaceData)
1265                                         COMPARE_ORIG_INDEX= dm->getFaceDataArray(dm, CD_ORIGINDEX);
1266                         }
1267
1268                         if(COMPARE_ORIG_INDEX) {
1269                                 qsort(index, totpart, sizeof(int), compare_orig_index);
1270                                 COMPARE_ORIG_INDEX = NULL;
1271                         }
1272                 }
1273         }
1274
1275         /* weights are no longer used except for FROM_PARTICLE, which needs them zeroed for indexing */
1276         if(from==PART_FROM_PARTICLE){
1277                 for(i=0; i<tot; i++)
1278                         weight[i]=0.0f;
1279         }
1280
1281         /* 4. */
1282         if(distr==PART_DISTR_JIT && ELEM(from,PART_FROM_FACE,PART_FROM_VOLUME)) {
1283                 jitlevel= part->userjit;
1284
1285                 if(jitlevel == 0) {
1286                         jitlevel= totpart/tot;
1287                         if(part->flag & PART_EDISTR) jitlevel*= 2;      /* looks better in general, not very scietific */
1288                         if(jitlevel<3) jitlevel= 3;
1289                 }
1290
1291                 jit= MEM_callocN((2+ jitlevel*2)*sizeof(float), "jit");
1292
1293                 /* for small amounts of particles we use regular jitter since it looks
1294                  * a bit better, for larger amounts we switch to hammersley sequence
1295                  * because it is much faster */
1296                 if(jitlevel < 25)
1297                         init_mv_jit(jit, jitlevel, psys->seed, part->jitfac);
1298                 else
1299                         hammersley_create(jit, jitlevel+1, psys->seed, part->jitfac);
1300                 BLI_array_randomize(jit, 2*sizeof(float), jitlevel, psys->seed); /* for custom jit or even distribution */
1301         }
1302
1303         /* 5. */
1304         ctx->tree= tree;
1305         ctx->seams= seams;
1306         ctx->totseam= totseam;
1307         ctx->sim.psys= psys;
1308         ctx->index= index;
1309         ctx->jit= jit;
1310         ctx->jitlevel= jitlevel;
1311         ctx->jitoff= jitoff;
1312         ctx->weight= weight;
1313         ctx->maxweight= maxweight;
1314         ctx->from= (children)? PART_FROM_CHILD: from;
1315         ctx->cfrom= cfrom;
1316         ctx->distr= distr;
1317         ctx->dm= dm;
1318         ctx->tpars= tpars;
1319
1320         if(children) {
1321                 totpart= psys_render_simplify_distribution(ctx, totpart);
1322                 alloc_child_particles(psys, totpart);
1323         }
1324
1325         if(!children || psys->totchild < 10000)
1326                 totthread= 1;
1327
1328         seed= 31415926 + ctx->sim.psys->seed;
1329         for(i=0; i<totthread; i++) {
1330                 threads[i].rng= rng_new(seed);
1331                 threads[i].tot= totthread;
1332         }
1333
1334         return 1;
1335 }
1336
1337 static void distribute_particles_on_dm(ParticleSimulationData *sim, int from)
1338 {
1339         DerivedMesh *finaldm = sim->psmd->dm;
1340         ListBase threads;
1341         ParticleThread *pthreads;
1342         ParticleThreadContext *ctx;
1343         int i, totthread;
1344
1345         pthreads= psys_threads_create(sim);
1346
1347         if(!psys_threads_init_distribution(pthreads, sim->scene, finaldm, from)) {
1348                 psys_threads_free(pthreads);
1349                 return;
1350         }
1351
1352         totthread= pthreads[0].tot;
1353         if(totthread > 1) {
1354                 BLI_init_threads(&threads, exec_distribution, totthread);
1355
1356                 for(i=0; i<totthread; i++)
1357                         BLI_insert_thread(&threads, &pthreads[i]);
1358
1359                 BLI_end_threads(&threads);
1360         }
1361         else
1362                 exec_distribution(&pthreads[0]);
1363
1364         psys_calc_dmcache(sim->ob, finaldm, sim->psys);
1365
1366         ctx= pthreads[0].ctx;
1367         if(ctx->dm != finaldm)
1368                 ctx->dm->release(ctx->dm);
1369
1370         psys_threads_free(pthreads);
1371 }
1372
1373 /* ready for future use, to emit particles without geometry */
1374 static void distribute_particles_on_shape(ParticleSimulationData *sim, int UNUSED(from))
1375 {
1376         ParticleSystem *psys = sim->psys;
1377         PARTICLE_P;
1378
1379         fprintf(stderr,"Shape emission not yet possible!\n");
1380
1381         LOOP_PARTICLES {
1382                 pa->fuv[0]=pa->fuv[1]=pa->fuv[2]=pa->fuv[3]= 0.0;
1383                 pa->foffset= 0.0f;
1384                 pa->num= -1;
1385         }
1386 }
1387 static void distribute_particles(ParticleSimulationData *sim, int from)
1388 {
1389         PARTICLE_PSMD;
1390         int distr_error=0;
1391
1392         if(psmd){
1393                 if(psmd->dm)
1394                         distribute_particles_on_dm(sim, from);
1395                 else
1396                         distr_error=1;
1397         }
1398         else
1399                 distribute_particles_on_shape(sim, from);
1400
1401         if(distr_error){
1402                 ParticleSystem *psys = sim->psys;
1403                 PARTICLE_P;
1404
1405                 fprintf(stderr,"Particle distribution error!\n");
1406
1407                 LOOP_PARTICLES {
1408                         pa->fuv[0]=pa->fuv[1]=pa->fuv[2]=pa->fuv[3]= 0.0;
1409                         pa->foffset= 0.0f;
1410                         pa->num= -1;
1411                 }
1412         }
1413 }
1414
1415 /* threaded child particle distribution and path caching */
1416 ParticleThread *psys_threads_create(ParticleSimulationData *sim)
1417 {
1418         ParticleThread *threads;
1419         ParticleThreadContext *ctx;
1420         int i, totthread;
1421
1422         if(sim->scene->r.mode & R_FIXED_THREADS)
1423                 totthread= sim->scene->r.threads;
1424         else
1425                 totthread= BLI_system_thread_count();
1426
1427         threads= MEM_callocN(sizeof(ParticleThread)*totthread, "ParticleThread");
1428         ctx= MEM_callocN(sizeof(ParticleThreadContext), "ParticleThreadContext");
1429
1430         ctx->sim = *sim;
1431         ctx->dm= ctx->sim.psmd->dm;
1432         ctx->ma= give_current_material(sim->ob, sim->psys->part->omat);
1433
1434         memset(threads, 0, sizeof(ParticleThread)*totthread);
1435
1436         for(i=0; i<totthread; i++) {
1437                 threads[i].ctx= ctx;
1438                 threads[i].num= i;
1439                 threads[i].tot= totthread;
1440         }
1441
1442         return threads;
1443 }
1444
1445 void psys_threads_free(ParticleThread *threads)
1446 {
1447         ParticleThreadContext *ctx= threads[0].ctx;
1448         int i, totthread= threads[0].tot;
1449
1450         /* path caching */
1451         if(ctx->vg_length)
1452                 MEM_freeN(ctx->vg_length);
1453         if(ctx->vg_clump)
1454                 MEM_freeN(ctx->vg_clump);
1455         if(ctx->vg_kink)
1456                 MEM_freeN(ctx->vg_kink);
1457         if(ctx->vg_rough1)
1458                 MEM_freeN(ctx->vg_rough1);
1459         if(ctx->vg_rough2)
1460                 MEM_freeN(ctx->vg_rough2);
1461         if(ctx->vg_roughe)
1462                 MEM_freeN(ctx->vg_roughe);
1463
1464         if(ctx->sim.psys->lattice){
1465                 end_latt_deform(ctx->sim.psys->lattice);
1466                 ctx->sim.psys->lattice= NULL;
1467         }
1468
1469         /* distribution */
1470         if(ctx->jit) MEM_freeN(ctx->jit);
1471         if(ctx->jitoff) MEM_freeN(ctx->jitoff);
1472         if(ctx->weight) MEM_freeN(ctx->weight);
1473         if(ctx->index) MEM_freeN(ctx->index);
1474         if(ctx->skip) MEM_freeN(ctx->skip);
1475         if(ctx->seams) MEM_freeN(ctx->seams);
1476         //if(ctx->vertpart) MEM_freeN(ctx->vertpart);
1477         BLI_kdtree_free(ctx->tree);
1478
1479         /* threads */
1480         for(i=0; i<totthread; i++) {
1481                 if(threads[i].rng)
1482                         rng_free(threads[i].rng);
1483                 if(threads[i].rng_path)
1484                         rng_free(threads[i].rng_path);
1485         }
1486
1487         MEM_freeN(ctx);
1488         MEM_freeN(threads);
1489 }
1490
1491 /* set particle parameters that don't change during particle's life */
1492 void initialize_particle(ParticleSimulationData *sim, ParticleData *pa, int p)
1493 {
1494         ParticleSettings *part = sim->psys->part;
1495         ParticleTexture ptex;
1496         Material *ma=0;
1497         //IpoCurve *icu=0; // XXX old animation system
1498         int totpart;
1499
1500         totpart=sim->psys->totpart;
1501
1502         ptex.life=ptex.size=ptex.exist=ptex.length=1.0;
1503         ptex.time=(float)p/(float)totpart;
1504
1505         BLI_srandom(sim->psys->seed + p + 125);
1506
1507         if(part->from!=PART_FROM_PARTICLE && part->type!=PART_FLUID){
1508                 ma=give_current_material(sim->ob,part->omat);
1509
1510                 /* TODO: needs some work to make most blendtypes generally usefull */
1511                 psys_get_texture(sim,ma,pa,&ptex,MAP_PA_INIT);
1512         }
1513
1514         if(part->type==PART_HAIR)
1515                 pa->time= 0.0f;
1516         //else if(part->type==PART_REACTOR && (part->flag&PART_REACT_STA_END)==0)
1517         //      pa->time= 300000.0f;    /* max frame */
1518         else{
1519                 //icu=find_ipocurve(psys->part->ipo,PART_EMIT_TIME);
1520                 //if(icu){
1521                 //      calc_icu(icu,100*ptex.time);
1522                 //      ptex.time=icu->curval;
1523                 //}
1524
1525                 pa->time= part->sta + (part->end - part->sta)*ptex.time;
1526         }
1527
1528         if(part->type!=PART_HAIR && part->distr!=PART_DISTR_GRID && part->from != PART_FROM_VERT){
1529                 if(ptex.exist < BLI_frand())
1530                         pa->flag |= PARS_UNEXIST;
1531                 else
1532                         pa->flag &= ~PARS_UNEXIST;
1533         }
1534
1535         pa->hair_index=0;
1536         /* we can't reset to -1 anymore since we've figured out correct index in distribute_particles */
1537         /* usage other than straight after distribute has to handle this index by itself - jahka*/
1538         //pa->num_dmcache = DMCACHE_NOTFOUND; /* assume we dont have a derived mesh face */
1539 }
1540 static void initialize_all_particles(ParticleSimulationData *sim)
1541 {
1542         //IpoCurve *icu=0; // XXX old animation system
1543         ParticleSystem *psys = sim->psys;
1544         PARTICLE_P;
1545
1546         LOOP_PARTICLES
1547                 initialize_particle(sim, pa, p);
1548
1549         if(psys->part->type != PART_FLUID) {
1550 #if 0 // XXX old animation system
1551                 icu=find_ipocurve(psys->part->ipo,PART_EMIT_FREQ);
1552                 if(icu){
1553                         float time=psys->part->sta, end=psys->part->end;
1554                         float v1, v2, a=0.0f, t1,t2, d;
1555
1556                         p=0;
1557                         pa=psys->particles;
1558
1559
1560                         calc_icu(icu,time);
1561                         v1=icu->curval;
1562                         if(v1<0.0f) v1=0.0f;
1563
1564                         calc_icu(icu,time+1.0f);
1565                         v2=icu->curval;
1566                         if(v2<0.0f) v2=0.0f;
1567
1568                         for(p=0, pa=psys->particles; p<totpart && time<end; p++, pa++){
1569                                 while(a+0.5f*(v1+v2) < (float)(p+1) && time<end){
1570                                         a+=0.5f*(v1+v2);
1571                                         v1=v2;
1572                                         time++;
1573                                         calc_icu(icu,time+1.0f);
1574                                         v2=icu->curval;
1575                                 }
1576                                 if(time<end){
1577                                         if(v1==v2){
1578                                                 pa->time=time+((float)(p+1)-a)/v1;
1579                                         }
1580                                         else{
1581                                                 d=(float)sqrt(v1*v1-2.0f*(v2-v1)*(a-(float)(p+1)));
1582                                                 t1=(-v1+d)/(v2-v1);
1583                                                 t2=(-v1-d)/(v2-v1);
1584
1585                                                 /* the root between 0-1 is the correct one */
1586                                                 if(t1>0.0f && t1<=1.0f)
1587                                                         pa->time=time+t1;
1588                                                 else
1589                                                         pa->time=time+t2;
1590                                         }
1591                                 }
1592
1593                                 pa->dietime = pa->time+pa->lifetime;
1594                                 pa->flag &= ~PARS_UNEXIST;
1595                         }
1596                         for(; p<totpart; p++, pa++){
1597                                 pa->flag |= PARS_UNEXIST;
1598                         }
1599                 }
1600 #endif // XXX old animation system
1601         }
1602 }
1603 /* sets particle to the emitter surface with initial velocity & rotation */
1604 void reset_particle(ParticleSimulationData *sim, ParticleData *pa, float dtime, float cfra)
1605 {
1606         Object *ob = sim->ob;
1607         ParticleSystem *psys = sim->psys;
1608         ParticleSettings *part;
1609         ParticleTexture ptex;
1610         ParticleKey state;
1611         //IpoCurve *icu=0; // XXX old animation system
1612         float fac, phasefac, nor[3]={0,0,0},loc[3],vel[3]={0.0,0.0,0.0},rot[4],q2[4];
1613         float r_vel[3],r_ave[3],r_rot[4],vec[3],p_vel[3]={0.0,0.0,0.0};
1614         float x_vec[3]={1.0,0.0,0.0}, utan[3]={0.0,1.0,0.0}, vtan[3]={0.0,0.0,1.0}, rot_vec[3]={0.0,0.0,0.0};
1615         float q_phase[4], r_phase;
1616         int p = pa - psys->particles;
1617         part=psys->part;
1618
1619         ptex.ivel=1.0;
1620         ptex.life=1.0;
1621
1622         /* we need to get every random even if they're not used so that they don't effect eachother */
1623         r_vel[0] = 2.0f * (PSYS_FRAND(p + 10) - 0.5f);
1624         r_vel[1] = 2.0f * (PSYS_FRAND(p + 11) - 0.5f);
1625         r_vel[2] = 2.0f * (PSYS_FRAND(p + 12) - 0.5f);
1626
1627         r_ave[0] = 2.0f * (PSYS_FRAND(p + 13) - 0.5f);
1628         r_ave[1] = 2.0f * (PSYS_FRAND(p + 14) - 0.5f);
1629         r_ave[2] = 2.0f * (PSYS_FRAND(p + 15) - 0.5f);
1630
1631         r_rot[0] = 2.0f * (PSYS_FRAND(p + 16) - 0.5f);
1632         r_rot[1] = 2.0f * (PSYS_FRAND(p + 17) - 0.5f);
1633         r_rot[2] = 2.0f * (PSYS_FRAND(p + 18) - 0.5f);
1634         r_rot[3] = 2.0f * (PSYS_FRAND(p + 19) - 0.5f);
1635         normalize_qt(r_rot);
1636
1637         r_phase = PSYS_FRAND(p + 20);
1638
1639         if(part->from==PART_FROM_PARTICLE){
1640                 float speed;
1641                 ParticleSimulationData tsim= {0};
1642                 tsim.scene= sim->scene;
1643                 tsim.ob= psys->target_ob ? psys->target_ob : ob;
1644                 tsim.psys = BLI_findlink(&tsim.ob->particlesystem, sim->psys->target_psys-1);
1645
1646                 state.time = pa->time;
1647                 if(pa->num == -1)
1648                         memset(&state, 0, sizeof(state));
1649                 else
1650                         psys_get_particle_state(&tsim, pa->num, &state, 1);
1651                 psys_get_from_key(&state, loc, nor, rot, 0);
1652
1653                 mul_qt_v3(rot, vtan);
1654                 mul_qt_v3(rot, utan);
1655
1656                 speed= normalize_v3_v3(p_vel, state.vel);
1657                 mul_v3_fl(p_vel, dot_v3v3(r_vel, p_vel));
1658                 VECSUB(p_vel, r_vel, p_vel);
1659                 normalize_v3(p_vel);
1660                 mul_v3_fl(p_vel, speed);
1661
1662                 VECCOPY(pa->fuv, loc); /* abusing pa->fuv (not used for "from particle") for storing emit location */
1663         }
1664         else{
1665                 /* get precise emitter matrix if particle is born */
1666                 if(part->type!=PART_HAIR && pa->time < cfra && pa->time >= sim->psys->cfra) {
1667                         /* we have to force RECALC_ANIM here since where_is_objec_time only does drivers */
1668                         BKE_animsys_evaluate_animdata(&sim->ob->id, sim->ob->adt, pa->time, ADT_RECALC_ANIM);
1669                         where_is_object_time(sim->scene, sim->ob, pa->time);
1670                 }
1671
1672                 /* get birth location from object               */
1673                 if(part->tanfac!=0.0)
1674                         psys_particle_on_emitter(sim->psmd, part->from,pa->num, pa->num_dmcache, pa->fuv,pa->foffset,loc,nor,utan,vtan,0,0);
1675                 else
1676                         psys_particle_on_emitter(sim->psmd, part->from,pa->num, pa->num_dmcache, pa->fuv,pa->foffset,loc,nor,0,0,0,0);
1677
1678                 /* get possible textural influence */
1679                 psys_get_texture(sim, give_current_material(sim->ob,part->omat), pa, &ptex, MAP_PA_IVEL|MAP_PA_LIFE);
1680
1681                 //if(vg_vel && pa->num != -1)
1682                 //      ptex.ivel*=psys_particle_value_from_verts(sim->psmd->dm,part->from,pa,vg_vel);
1683
1684                 /* particles live in global space so    */
1685                 /* let's convert:                                               */
1686                 /* -location                                                    */
1687                 mul_m4_v3(ob->obmat,loc);
1688
1689                 /* -normal                                                              */
1690                 mul_mat3_m4_v3(ob->obmat,nor);
1691                 normalize_v3(nor);
1692
1693                 /* -tangent                                                             */
1694                 if(part->tanfac!=0.0){
1695                         //float phase=vg_rot?2.0f*(psys_particle_value_from_verts(sim->psmd->dm,part->from,pa,vg_rot)-0.5f):0.0f;
1696                         float phase=0.0f;
1697                         mul_v3_fl(vtan,-(float)cos(M_PI*(part->tanphase+phase)));
1698                         fac=-(float)sin(M_PI*(part->tanphase+phase));
1699                         VECADDFAC(vtan,vtan,utan,fac);
1700
1701                         mul_mat3_m4_v3(ob->obmat,vtan);
1702
1703                         VECCOPY(utan,nor);
1704                         mul_v3_fl(utan,dot_v3v3(vtan,nor));
1705                         VECSUB(vtan,vtan,utan);
1706
1707                         normalize_v3(vtan);
1708                 }
1709
1710
1711                 /* -velocity                                                    */
1712                 if(part->randfac!=0.0){
1713                         mul_mat3_m4_v3(ob->obmat,r_vel);
1714                         normalize_v3(r_vel);
1715                 }
1716
1717                 /* -angular velocity                                    */
1718                 if(part->avemode==PART_AVE_RAND){
1719                         mul_mat3_m4_v3(ob->obmat,r_ave);
1720                         normalize_v3(r_ave);
1721                 }
1722
1723                 /* -rotation                                                    */
1724                 if(part->randrotfac != 0.0f){
1725                         mat4_to_quat(rot,ob->obmat);
1726                         mul_qt_qtqt(r_rot,r_rot,rot);
1727                 }
1728         }
1729
1730         if(part->phystype==PART_PHYS_BOIDS && pa->boid) {
1731                 BoidParticle *bpa = pa->boid;
1732                 float dvec[3], q[4], mat[3][3];
1733
1734                 VECCOPY(pa->state.co,loc);
1735
1736                 /* boids don't get any initial velocity  */
1737                 pa->state.vel[0]=pa->state.vel[1]=pa->state.vel[2]=0.0f;
1738
1739                 /* boids store direction in ave */
1740                 if(fabs(nor[2])==1.0f) {
1741                         sub_v3_v3v3(pa->state.ave, loc, ob->obmat[3]);
1742                         normalize_v3(pa->state.ave);
1743                 }
1744                 else {
1745                         VECCOPY(pa->state.ave, nor);
1746                 }
1747                 /* and gravity in r_ve */
1748                 bpa->gravity[0] = bpa->gravity[1] = 0.0f;
1749                 bpa->gravity[2] = -1.0f;
1750                 if((sim->scene->physics_settings.flag & PHYS_GLOBAL_GRAVITY)
1751                         && sim->scene->physics_settings.gravity[2]!=0.0f)
1752                         bpa->gravity[2] = sim->scene->physics_settings.gravity[2];
1753
1754                 /* calculate rotation matrix */
1755                 project_v3_v3v3(dvec, r_vel, pa->state.ave);
1756                 sub_v3_v3v3(mat[0], pa->state.ave, dvec);
1757                 normalize_v3(mat[0]);
1758                 negate_v3_v3(mat[2], r_vel);
1759                 normalize_v3(mat[2]);
1760                 cross_v3_v3v3(mat[1], mat[2], mat[0]);
1761
1762                 /* apply rotation */
1763                 mat3_to_quat_is_ok( q,mat);
1764                 copy_qt_qt(pa->state.rot, q);
1765
1766                 bpa->data.health = part->boids->health;
1767                 bpa->data.mode = eBoidMode_InAir;
1768                 bpa->data.state_id = ((BoidState*)part->boids->states.first)->id;
1769                 bpa->data.acc[0]=bpa->data.acc[1]=bpa->data.acc[2]=0.0f;
1770         }
1771         else {
1772                 /* conversion done so now we apply new: */
1773                 /* -velocity from:                                              */
1774
1775                 /*              *reactions                                              */
1776                 if(dtime>0.0f){
1777                         VECSUB(vel,pa->state.vel,pa->prev_state.vel);
1778                 }
1779
1780                 /*              *emitter velocity                               */
1781                 if(dtime!=0.0 && part->obfac!=0.0){
1782                         VECSUB(vel,loc,pa->state.co);
1783                         mul_v3_fl(vel,part->obfac/dtime);
1784                 }
1785
1786                 /*              *emitter normal                                 */
1787                 if(part->normfac!=0.0)
1788                         VECADDFAC(vel,vel,nor,part->normfac);
1789
1790                 /*              *emitter tangent                                */
1791                 if(sim->psmd && part->tanfac!=0.0)
1792                         VECADDFAC(vel,vel,vtan,part->tanfac);
1793                         //VECADDFAC(vel,vel,vtan,part->tanfac*(vg_tan?psys_particle_value_from_verts(sim->psmd->dm,part->from,pa,vg_tan):1.0f));
1794
1795                 /*              *emitter object orientation             */
1796                 if(part->ob_vel[0]!=0.0) {
1797                         normalize_v3_v3(vec, ob->obmat[0]);
1798                         VECADDFAC(vel, vel, vec, part->ob_vel[0]);
1799                 }
1800                 if(part->ob_vel[1]!=0.0) {
1801                         normalize_v3_v3(vec, ob->obmat[1]);
1802                         VECADDFAC(vel, vel, vec, part->ob_vel[1]);
1803                 }
1804                 if(part->ob_vel[2]!=0.0) {
1805                         normalize_v3_v3(vec, ob->obmat[2]);
1806                         VECADDFAC(vel, vel, vec, part->ob_vel[2]);
1807                 }
1808
1809                 /*              *texture                                                */
1810                 /* TODO */
1811
1812                 /*              *random                                                 */
1813                 if(part->randfac!=0.0)
1814                         VECADDFAC(vel,vel,r_vel,part->randfac);
1815
1816                 /*              *particle                                               */
1817                 if(part->partfac!=0.0)
1818                         VECADDFAC(vel,vel,p_vel,part->partfac);
1819
1820                 //icu=find_ipocurve(psys->part->ipo,PART_EMIT_VEL);
1821                 //if(icu){
1822                 //      calc_icu(icu,100*((pa->time-part->sta)/(part->end-part->sta)));
1823                 //      ptex.ivel*=icu->curval;
1824                 //}
1825
1826                 mul_v3_fl(vel,ptex.ivel);
1827
1828                 VECCOPY(pa->state.vel,vel);
1829
1830                 /* -location from emitter                               */
1831                 VECCOPY(pa->state.co,loc);
1832
1833                 /* -rotation                                                    */
1834                 pa->state.rot[0]=1.0;
1835                 pa->state.rot[1]=pa->state.rot[2]=pa->state.rot[3]=0.0;
1836
1837                 if(part->rotmode){
1838                         /* create vector into which rotation is aligned */
1839                         switch(part->rotmode){
1840                                 case PART_ROT_NOR:
1841                                         copy_v3_v3(rot_vec, nor);
1842                                         break;
1843                                 case PART_ROT_VEL:
1844                                         copy_v3_v3(rot_vec, vel);
1845                                         break;
1846                                 case PART_ROT_GLOB_X:
1847                                 case PART_ROT_GLOB_Y:
1848                                 case PART_ROT_GLOB_Z:
1849                                         rot_vec[part->rotmode - PART_ROT_GLOB_X] = 1.0f;
1850                                         break;
1851                                 case PART_ROT_OB_X:
1852                                 case PART_ROT_OB_Y:
1853                                 case PART_ROT_OB_Z:
1854                                         copy_v3_v3(rot_vec, ob->obmat[part->rotmode - PART_ROT_OB_X]);
1855                                         break;
1856                         }
1857
1858                         /* create rotation quat */
1859                         negate_v3(rot_vec);
1860                         vec_to_quat( q2,rot_vec, OB_POSX, OB_POSZ);
1861
1862                         /* randomize rotation quat */
1863                         if(part->randrotfac!=0.0f)
1864                                 interp_qt_qtqt(rot, q2, r_rot, part->randrotfac);
1865                         else
1866                                 copy_qt_qt(rot,q2);
1867
1868                         /* rotation phase */
1869                         phasefac = part->phasefac;
1870                         if(part->randphasefac != 0.0f)
1871                                 phasefac += part->randphasefac * r_phase;
1872                         axis_angle_to_quat( q_phase,x_vec, phasefac*(float)M_PI);
1873
1874                         /* combine base rotation & phase */
1875                         mul_qt_qtqt(pa->state.rot, rot, q_phase);
1876                 }
1877
1878                 /* -angular velocity                                    */
1879
1880                 pa->state.ave[0] = pa->state.ave[1] = pa->state.ave[2] = 0.0;
1881
1882                 if(part->avemode){
1883                         switch(part->avemode){
1884                                 case PART_AVE_SPIN:
1885                                         VECCOPY(pa->state.ave,vel);
1886                                         break;
1887                                 case PART_AVE_RAND:
1888                                         VECCOPY(pa->state.ave,r_ave);
1889                                         break;
1890                         }
1891                         normalize_v3(pa->state.ave);
1892                         mul_v3_fl(pa->state.ave,part->avefac);
1893
1894                         //icu=find_ipocurve(psys->part->ipo,PART_EMIT_AVE);
1895                         //if(icu){
1896                         //      calc_icu(icu,100*((pa->time-part->sta)/(part->end-part->sta)));
1897                         //      mul_v3_fl(pa->state.ave,icu->curval);
1898                         //}
1899                 }
1900         }
1901
1902
1903         if(part->type == PART_HAIR){
1904                 pa->lifetime = 100.0f;
1905         }
1906         else{
1907                 pa->lifetime = part->lifetime*ptex.life;
1908
1909                 if(part->randlife != 0.0)
1910                         pa->lifetime *= 1.0f - part->randlife * PSYS_FRAND(p + 21);
1911         }
1912
1913         pa->dietime = pa->time + pa->lifetime;
1914
1915         if(sim->psys->pointcache && sim->psys->pointcache->flag & PTCACHE_BAKED &&
1916                 sim->psys->pointcache->mem_cache.first) {
1917                 float dietime = psys_get_dietime_from_cache(sim->psys->pointcache, p);
1918                 pa->dietime = MIN2(pa->dietime, dietime);
1919         }
1920
1921         if(pa->time > cfra)
1922                 pa->alive = PARS_UNBORN;
1923         else if(pa->dietime <= cfra)
1924                 pa->alive = PARS_DEAD;
1925         else
1926                 pa->alive = PARS_ALIVE;
1927
1928         pa->state.time = cfra;
1929 }
1930 static void reset_all_particles(ParticleSimulationData *sim, float dtime, float cfra, int from)
1931 {
1932         ParticleData *pa;
1933         int p, totpart=sim->psys->totpart;
1934         //float *vg_vel=psys_cache_vgroup(sim->psmd->dm,sim->psys,PSYS_VG_VEL);
1935         //float *vg_tan=psys_cache_vgroup(sim->psmd->dm,sim->psys,PSYS_VG_TAN);
1936         //float *vg_rot=psys_cache_vgroup(sim->psmd->dm,sim->psys,PSYS_VG_ROT);
1937
1938         for(p=from, pa=sim->psys->particles+from; p<totpart; p++, pa++)
1939                 reset_particle(sim, pa, dtime, cfra);
1940
1941         //if(vg_vel)
1942         //      MEM_freeN(vg_vel);
1943 }
1944 /************************************************/
1945 /*                      Particle targets                                        */
1946 /************************************************/
1947 ParticleSystem *psys_get_target_system(Object *ob, ParticleTarget *pt)
1948 {
1949         ParticleSystem *psys = NULL;
1950
1951         if(pt->ob == NULL || pt->ob == ob)
1952                 psys = BLI_findlink(&ob->particlesystem, pt->psys-1);
1953         else
1954                 psys = BLI_findlink(&pt->ob->particlesystem, pt->psys-1);
1955
1956         if(psys)
1957                 pt->flag |= PTARGET_VALID;
1958         else
1959                 pt->flag &= ~PTARGET_VALID;
1960
1961         return psys;
1962 }
1963 /************************************************/
1964 /*                      Keyed particles                                         */
1965 /************************************************/
1966 /* Counts valid keyed targets */
1967 void psys_count_keyed_targets(ParticleSimulationData *sim)
1968 {
1969         ParticleSystem *psys = sim->psys, *kpsys;
1970         ParticleTarget *pt = psys->targets.first;
1971         int keys_valid = 1;
1972         psys->totkeyed = 0;
1973
1974         for(; pt; pt=pt->next) {
1975                 kpsys = psys_get_target_system(sim->ob, pt);
1976
1977                 if(kpsys && kpsys->totpart) {
1978                         psys->totkeyed += keys_valid;
1979                         if(psys->flag & PSYS_KEYED_TIMING && pt->duration != 0.0f)
1980                                 psys->totkeyed += 1;
1981                 }
1982                 else {
1983                         keys_valid = 0;
1984                 }
1985         }
1986
1987         psys->totkeyed *= psys->flag & PSYS_KEYED_TIMING ? 1 : psys->part->keyed_loops;
1988 }
1989
1990 static void set_keyed_keys(ParticleSimulationData *sim)
1991 {
1992         ParticleSystem *psys = sim->psys;
1993         ParticleSimulationData ksim= {0};
1994         ParticleTarget *pt;
1995         PARTICLE_P;
1996         ParticleKey *key;
1997         int totpart = psys->totpart, k, totkeys = psys->totkeyed;
1998
1999         ksim.scene= sim->scene;
2000
2001         /* no proper targets so let's clear and bail out */
2002         if(psys->totkeyed==0) {
2003                 free_keyed_keys(psys);
2004                 psys->flag &= ~PSYS_KEYED;
2005                 return;
2006         }
2007
2008         if(totpart && psys->particles->totkey != totkeys) {
2009                 free_keyed_keys(psys);
2010
2011                 key = MEM_callocN(totpart*totkeys*sizeof(ParticleKey), "Keyed keys");
2012
2013                 LOOP_PARTICLES {
2014                         pa->keys = key;
2015                         pa->totkey = totkeys;
2016                         key += totkeys;
2017                 }
2018         }
2019
2020         psys->flag &= ~PSYS_KEYED;
2021
2022
2023         pt = psys->targets.first;
2024         for(k=0; k<totkeys; k++) {
2025                 ksim.ob = pt->ob ? pt->ob : sim->ob;
2026                 ksim.psys = BLI_findlink(&ksim.ob->particlesystem, pt->psys - 1);
2027
2028                 LOOP_PARTICLES {
2029                         key = pa->keys + k;
2030                         key->time = -1.0; /* use current time */
2031
2032                         psys_get_particle_state(&ksim, p%ksim.psys->totpart, key, 1);
2033
2034                         if(psys->flag & PSYS_KEYED_TIMING){
2035                                 key->time = pa->time + pt->time;
2036                                 if(pt->duration != 0.0f && k+1 < totkeys) {
2037                                         copy_particle_key(key+1, key, 1);
2038                                         (key+1)->time = pa->time + pt->time + pt->duration;
2039                                 }
2040                         }
2041                         else if(totkeys > 1)
2042                                 key->time = pa->time + (float)k / (float)(totkeys - 1) * pa->lifetime;
2043                         else
2044                                 key->time = pa->time;
2045                 }
2046
2047                 if(psys->flag & PSYS_KEYED_TIMING && pt->duration!=0.0f)
2048                         k++;
2049
2050                 pt = (pt->next && pt->next->flag & PTARGET_VALID)? pt->next : psys->targets.first;
2051         }
2052
2053         psys->flag |= PSYS_KEYED;
2054 }
2055 /************************************************/
2056 /*                      Reactors                                                        */
2057 /************************************************/
2058 //static void push_reaction(ParticleSimulationData *sim, int pa_num, int event, ParticleKey *state)
2059 //{
2060 //      Object *rob;
2061 //      ParticleSystem *rpsys;
2062 //      ParticleSettings *rpart;
2063 //      ParticleData *pa;
2064 //      ListBase *lb=&sim->psys->effectors;
2065 //      ParticleEffectorCache *ec;
2066 //      ParticleReactEvent *re;
2067 //
2068 //      if(lb->first) for(ec = lb->first; ec; ec= ec->next){
2069 //              if(ec->type & PSYS_EC_REACTOR){
2070 //                      /* all validity checks already done in add_to_effectors */
2071 //                      rob=ec->ob;
2072 //                      rpsys=BLI_findlink(&rob->particlesystem,ec->psys_nbr);
2073 //                      rpart=rpsys->part;
2074 //                      if(rpsys->part->reactevent==event){
2075 //                              pa=sim->psys->particles+pa_num;
2076 //                              re= MEM_callocN(sizeof(ParticleReactEvent), "react event");
2077 //                              re->event=event;
2078 //                              re->pa_num = pa_num;
2079 //                              re->ob = sim->ob;
2080 //                              re->psys = sim->psys;
2081 //                              re->size = pa->size;
2082 //                              copy_particle_key(&re->state,state,1);
2083 //
2084 //                              switch(event){
2085 //                                      case PART_EVENT_DEATH:
2086 //                                              re->time=pa->dietime;
2087 //                                              break;
2088 //                                      case PART_EVENT_COLLIDE:
2089 //                                              re->time=state->time;
2090 //                                              break;
2091 //                                      case PART_EVENT_NEAR:
2092 //                                              re->time=state->time;
2093 //                                              break;
2094 //                              }
2095 //
2096 //                              BLI_addtail(&rpsys->reactevents, re);
2097 //                      }
2098 //              }
2099 //      }
2100 //}
2101 //static void react_to_events(ParticleSystem *psys, int pa_num)
2102 //{
2103 //      ParticleSettings *part=psys->part;
2104 //      ParticleData *pa=psys->particles+pa_num;
2105 //      ParticleReactEvent *re=psys->reactevents.first;
2106 //      int birth=0;
2107 //      float dist=0.0f;
2108 //
2109 //      for(re=psys->reactevents.first; re; re=re->next){
2110 //              birth=0;
2111 //              if(part->from==PART_FROM_PARTICLE){
2112 //                      if(pa->num==re->pa_num && pa->alive==PARS_UNBORN){
2113 //                              if(re->event==PART_EVENT_NEAR){
2114 //                                      ParticleData *tpa = re->psys->particles+re->pa_num;
2115 //                                      float pa_time=tpa->time + pa->foffset*tpa->lifetime;
2116 //                                      if(re->time >= pa_time){
2117 //                                              pa->time=pa_time;
2118 //                                              pa->dietime=pa->time+pa->lifetime;
2119 //                                      }
2120 //                              }
2121 //                              else{
2122 //                                      pa->time=re->time;
2123 //                                      pa->dietime=pa->time+pa->lifetime;
2124 //                              }
2125 //                      }
2126 //              }
2127 //              else{
2128 //                      dist=len_v3v3(pa->state.co, re->state.co);
2129 //                      if(dist <= re->size){
2130 //                              if(pa->alive==PARS_UNBORN){
2131 //                                      pa->time=re->time;
2132 //                                      pa->dietime=pa->time+pa->lifetime;
2133 //                                      birth=1;
2134 //                              }
2135 //                              if(birth || part->flag&PART_REACT_MULTIPLE){
2136 //                                      float vec[3];
2137 //                                      VECSUB(vec,pa->state.co, re->state.co);
2138 //                                      if(birth==0)
2139 //                                              mul_v3_fl(vec,(float)pow(1.0f-dist/re->size,part->reactshape));
2140 //                                      VECADDFAC(pa->state.vel,pa->state.vel,vec,part->reactfac);
2141 //                                      VECADDFAC(pa->state.vel,pa->state.vel,re->state.vel,part->partfac);
2142 //                              }
2143 //                              if(birth)
2144 //                                      mul_v3_fl(pa->state.vel,(float)pow(1.0f-dist/re->size,part->reactshape));
2145 //                      }
2146 //              }
2147 //      }
2148 //}
2149 //void psys_get_reactor_target(ParticleSimulationData *sim, Object **target_ob, ParticleSystem **target_psys)
2150 //{
2151 //      Object *tob;
2152 //
2153 //      tob = sim->psys->target_ob ? sim->psys->target_ob : sim->ob;
2154 //
2155 //      *target_psys = BLI_findlink(&tob->particlesystem, sim->psys->target_psys-1);
2156 //      if(*target_psys)
2157 //              *target_ob=tob;
2158 //      else
2159 //              *target_ob=0;
2160 //}
2161 /************************************************/
2162 /*                      Point Cache                                                     */
2163 /************************************************/
2164 void psys_make_temp_pointcache(Object *ob, ParticleSystem *psys)
2165 {
2166         PointCache *cache = psys->pointcache;
2167
2168         if(cache->flag & PTCACHE_DISK_CACHE && cache->mem_cache.first == NULL) {
2169                 PTCacheID pid;
2170                 BKE_ptcache_id_from_particles(&pid, ob, psys);
2171                 cache->flag &= ~PTCACHE_DISK_CACHE;
2172                 BKE_ptcache_disk_to_mem(&pid);
2173                 cache->flag |= PTCACHE_DISK_CACHE;
2174         }
2175 }
2176 static void psys_clear_temp_pointcache(ParticleSystem *psys)
2177 {
2178         if(psys->pointcache->flag & PTCACHE_DISK_CACHE)
2179                 BKE_ptcache_free_mem(&psys->pointcache->mem_cache);
2180 }
2181 void psys_get_pointcache_start_end(Scene *scene, ParticleSystem *psys, int *sfra, int *efra)
2182 {
2183         ParticleSettings *part = psys->part;
2184
2185         *sfra = MAX2(1, (int)part->sta);
2186         *efra = MIN2((int)(part->end + part->lifetime + 1.0), scene->r.efra);
2187 }
2188
2189 /************************************************/
2190 /*                      Effectors                                                       */
2191 /************************************************/
2192 void psys_update_particle_tree(ParticleSystem *psys, float cfra)
2193 {
2194         if(psys) {
2195                 PARTICLE_P;
2196
2197                 if(!psys->tree || psys->tree_frame != cfra) {
2198
2199                         BLI_kdtree_free(psys->tree);
2200
2201                         psys->tree = BLI_kdtree_new(psys->totpart);
2202
2203                         LOOP_SHOWN_PARTICLES {
2204                                 if(pa->alive == PARS_ALIVE) {
2205                                         if(pa->state.time == cfra)
2206                                                 BLI_kdtree_insert(psys->tree, p, pa->prev_state.co, NULL);
2207                                         else
2208                                                 BLI_kdtree_insert(psys->tree, p, pa->state.co, NULL);
2209                                 }
2210                         }
2211                         BLI_kdtree_balance(psys->tree);
2212
2213                         psys->tree_frame = psys->cfra;
2214                 }
2215         }
2216 }
2217
2218 static void psys_update_effectors(ParticleSimulationData *sim)
2219 {
2220         pdEndEffectors(&sim->psys->effectors);
2221         sim->psys->effectors = pdInitEffectors(sim->scene, sim->ob, sim->psys, sim->psys->part->effector_weights);
2222         precalc_guides(sim, sim->psys->effectors);
2223 }
2224
2225 /*********************************************************************************************************
2226                     SPH fluid physics
2227
2228  In theory, there could be unlimited implementation of SPH simulators
2229
2230  This code uses in some parts adapted algorithms from the pseudo code as outlined in the Research paper:
2231
2232  Titled: Particle-based Viscoelastic Fluid Simulation.
2233  Authors: Simon Clavet, Philippe Beaudoin and Pierre Poulin
2234  Website: http://www.iro.umontreal.ca/labs/infographie/papers/Clavet-2005-PVFS/
2235
2236  Presented at Siggraph, (2005)
2237
2238 ***********************************************************************************************************/
2239 #define PSYS_FLUID_SPRINGS_INITIAL_SIZE 256
2240 ParticleSpring *add_fluid_spring(ParticleSystem *psys, ParticleSpring *spring)
2241 {
2242         /* Are more refs required? */
2243         if(psys->alloc_fluidsprings == 0 || psys->fluid_springs == NULL) {
2244                 psys->alloc_fluidsprings = PSYS_FLUID_SPRINGS_INITIAL_SIZE;
2245                 psys->fluid_springs = (ParticleSpring*)MEM_callocN(psys->alloc_fluidsprings * sizeof(ParticleSpring), "Particle Fluid Springs");
2246         }
2247         else if(psys->tot_fluidsprings == psys->alloc_fluidsprings) {
2248                 /* Double the number of refs allocated */
2249                 psys->alloc_fluidsprings *= 2;
2250                 psys->fluid_springs = (ParticleSpring*)MEM_reallocN(psys->fluid_springs, psys->alloc_fluidsprings * sizeof(ParticleSpring));
2251         }
2252
2253         memcpy(psys->fluid_springs + psys->tot_fluidsprings, spring, sizeof(ParticleSpring));
2254         psys->tot_fluidsprings++;
2255
2256         return psys->fluid_springs + psys->tot_fluidsprings - 1;
2257 }
2258
2259 void  delete_fluid_spring(ParticleSystem *psys, int j)
2260 {
2261         if (j != psys->tot_fluidsprings - 1)
2262                 psys->fluid_springs[j] = psys->fluid_springs[psys->tot_fluidsprings - 1];
2263
2264         psys->tot_fluidsprings--;
2265
2266         if (psys->tot_fluidsprings < psys->alloc_fluidsprings/2 && psys->alloc_fluidsprings > PSYS_FLUID_SPRINGS_INITIAL_SIZE){
2267                 psys->alloc_fluidsprings /= 2;
2268                 psys->fluid_springs = (ParticleSpring*)MEM_reallocN(psys->fluid_springs,  psys->alloc_fluidsprings * sizeof(ParticleSpring));
2269         }
2270 }
2271
2272 EdgeHash *build_fluid_springhash(ParticleSystem *psys)
2273 {
2274         EdgeHash *springhash = NULL;
2275         ParticleSpring *spring;
2276         int i = 0;
2277
2278         springhash = BLI_edgehash_new();
2279
2280         for(i=0, spring=psys->fluid_springs; i<psys->tot_fluidsprings; i++, spring++)
2281                 BLI_edgehash_insert(springhash, spring->particle_index[0], spring->particle_index[1], SET_INT_IN_POINTER(i+1));
2282
2283         return springhash;
2284 }
2285 static void particle_fluidsim(ParticleSystem *psys, int own_psys, ParticleData *pa, float dtime, float mass, float *gravity, EdgeHash *springhash)
2286 {
2287         SPHFluidSettings *fluid = psys->part->fluid;
2288         KDTreeNearest *ptn = NULL;
2289         ParticleData *npa;
2290         ParticleSpring *spring = NULL;
2291
2292         float temp[3];
2293         float q, q1, u, I, D, rij, d, Lij;
2294         float pressure_near, pressure;
2295         float p=0, pnear=0;
2296
2297         float omega = fluid->viscosity_omega;
2298         float beta = fluid->viscosity_beta;
2299         float massfactor = 1.0f/mass;
2300         float spring_k = fluid->spring_k;
2301         float h = fluid->radius;
2302         float L = fluid->rest_length * fluid->radius;
2303
2304         int n, neighbours = BLI_kdtree_range_search(psys->tree, h, pa->prev_state.co, NULL, &ptn);
2305         int spring_index = 0, index = own_psys ? pa - psys->particles : -1;
2306
2307         /* pressure and near pressure */
2308         for(n=own_psys?1:0; n<neighbours; n++) {
2309                 sub_v3_v3(ptn[n].co, pa->prev_state.co);
2310                 mul_v3_fl(ptn[n].co, 1.f/ptn[n].dist);
2311                 q = ptn[n].dist/h;
2312
2313                 if(q < 1.f) {
2314                         q1 = 1.f - q;
2315
2316                         p += q1*q1;
2317                         pnear += q1*q1*q1;
2318                 }
2319         }
2320
2321         p *= mass;
2322         pnear *= mass;
2323         pressure =  fluid->stiffness_k * (p - fluid->rest_density);
2324         pressure_near = fluid->stiffness_knear * pnear;
2325
2326         /* main calculations */
2327         for(n=own_psys?1:0; n<neighbours; n++) {
2328                 npa = psys->particles + ptn[n].index;
2329
2330                 rij = ptn[n].dist;
2331                 q = rij/h;
2332                 q1 = 1.f-q;
2333
2334                 /* Double Density Relaxation - Algorithm 2 (can't be thread safe!)*/
2335                 D =  dtime * dtime * (pressure + pressure_near*q1)*q1 * 0.5f;
2336                 madd_v3_v3fl(pa->state.co, ptn[n].co, -D * massfactor);
2337                 if(own_psys)
2338                         madd_v3_v3fl(npa->state.co, ptn[n].co, D * massfactor);
2339
2340                 if(index < ptn[n].index) {
2341                         /* Viscosity - Algorithm 5 */
2342                         if(omega > 0.f  || beta > 0.f) {
2343                                 sub_v3_v3v3(temp, pa->state.vel, npa->state.vel);
2344                                 u = dot_v3v3(ptn[n].co, temp);
2345
2346                                 if (u > 0){
2347                                         I = dtime * (q1 * (omega * u + beta * u*u)) * 0.5f;
2348                                         madd_v3_v3fl(pa->state.vel, ptn[n].co, -I * massfactor);
2349
2350                                         if(own_psys)
2351                                                 madd_v3_v3fl(npa->state.vel, ptn[n].co, I * massfactor);
2352                                 }
2353                         }
2354
2355                         if(spring_k > 0.f) {
2356                                 /* Viscoelastic spring force - Algorithm 4*/
2357                                 if (fluid->flag & SPH_VISCOELASTIC_SPRINGS && springhash){
2358                                         spring_index = GET_INT_FROM_POINTER(BLI_edgehash_lookup(springhash, index, ptn[n].index));
2359
2360                                         if(spring_index) {
2361                                                 spring = psys->fluid_springs + spring_index - 1;
2362                                         }
2363                                         else {
2364                                                 ParticleSpring temp_spring;
2365                                                 temp_spring.particle_index[0] = index;
2366                                                 temp_spring.particle_index[1] = ptn[n].index;
2367                                                 temp_spring.rest_length = (fluid->flag & SPH_CURRENT_REST_LENGTH) ? rij : L;
2368                                                 temp_spring.delete_flag = 0;
2369
2370                                                 spring = add_fluid_spring(psys, &temp_spring);
2371                                         }
2372
2373                                         Lij = spring->rest_length;
2374                                         d = fluid->yield_ratio * Lij;
2375
2376                                         if (rij > Lij + d) // Stretch, 25 is just a multiplier for plasticity_constant value to counter default dtime of 1/25
2377                                                 spring->rest_length += dtime * 25.f * fluid->plasticity_constant * (rij - Lij - d);
2378                                         else if(rij < Lij - d) // Compress
2379                                                 spring->rest_length -= dtime * 25.f * fluid->plasticity_constant * (Lij - d - rij);
2380                                 }
2381                                 else { /* PART_SPRING_HOOKES - Hooke's spring force */
2382                                         /* L is a factor of radius */
2383                                         D = 0.5 * dtime * dtime * 10.f * fluid->spring_k * (1.f - L/h) * (L - rij);
2384
2385                                         madd_v3_v3fl(pa->state.co, ptn[n].co, -D * massfactor);
2386                                         if(own_psys)
2387                                                 madd_v3_v3fl(npa->state.co, ptn[n].co, D * massfactor);
2388                                 }
2389                         }
2390                 }
2391         }
2392
2393         /* Artificial buoyancy force in negative gravity direction  */
2394         if (fluid->buoyancy >= 0.f && gravity) {
2395                 float B = -dtime * dtime * fluid->buoyancy * (p - fluid->rest_density) * 0.5f;
2396                 madd_v3_v3fl(pa->state.co, gravity, -B * massfactor);
2397         }
2398
2399         if(ptn)
2400                 MEM_freeN(ptn);
2401 }
2402
2403 static void apply_particle_fluidsim(Object *ob, ParticleSystem *psys, ParticleData *pa, float dtime, float *gravity, EdgeHash *springhash){
2404         ParticleTarget *pt;
2405
2406         particle_fluidsim(psys, 1, pa, dtime, psys->part->mass, gravity, springhash);
2407
2408         /*----check other SPH systems (Multifluids) , each fluid has its own parameters---*/
2409         for(pt=psys->targets.first; pt; pt=pt->next) {
2410                 ParticleSystem *epsys = psys_get_target_system(ob, pt);
2411
2412                 if(epsys)
2413                         particle_fluidsim(epsys, 0, pa, dtime, psys->part->mass, gravity, NULL);
2414         }
2415         /*----------------------------------------------------------------*/
2416 }
2417
2418 static void apply_fluid_springs(ParticleSystem *psys, float timestep){
2419         SPHFluidSettings *fluid = psys->part->fluid;
2420         ParticleData *pa1, *pa2;
2421         ParticleSpring *spring = psys->fluid_springs;
2422
2423         float h = fluid->radius;
2424         float massfactor = 1.0f/psys->part->mass;
2425         float D, Rij[3], rij, Lij;
2426         int i;
2427
2428         if((fluid->flag & SPH_VISCOELASTIC_SPRINGS)==0 || fluid->spring_k == 0.f)
2429                 return;
2430
2431         /* Loop through the springs */
2432         for(i=0; i<psys->tot_fluidsprings; i++, spring++) {
2433                 Lij = spring->rest_length;
2434
2435                 if (Lij > h) {
2436                         spring->delete_flag = 1;
2437                 }
2438                 else {
2439                         pa1 = psys->particles + spring->particle_index[0];
2440                         pa2 = psys->particles + spring->particle_index[1];
2441
2442                         sub_v3_v3v3(Rij, pa2->prev_state.co, pa1->prev_state.co);
2443                         rij = normalize_v3(Rij);
2444
2445                         /* Calculate displacement and apply value */
2446                         D =  0.5f * timestep * timestep * 10.f * fluid->spring_k * (1.f - Lij/h) * (Lij - rij);
2447
2448                         madd_v3_v3fl(pa1->state.co, Rij, -D * pa1->state.time * pa1->state.time * massfactor);
2449                         madd_v3_v3fl(pa2->state.co, Rij, D * pa2->state.time * pa2->state.time * massfactor);
2450                 }
2451         }
2452
2453         /* Loop through springs backwaqrds - for efficient delete function */
2454         for (i=psys->tot_fluidsprings-1; i >= 0; i--) {
2455                 if(psys->fluid_springs[i].delete_flag)
2456                         delete_fluid_spring(psys, i);
2457         }
2458 }
2459
2460 /************************************************/
2461 /*                      Newtonian physics                                       */
2462 /************************************************/
2463 /* gathers all forces that effect particles and calculates a new state for the particle */
2464 static void apply_particle_forces(ParticleSimulationData *sim, int p, float dfra, float cfra)
2465 {
2466         ParticleSettings *part = sim->psys->part;
2467         ParticleData *pa = sim->psys->particles + p;
2468         EffectedPoint epoint;
2469         ParticleKey states[5], tkey;
2470         float timestep = psys_get_timestep(sim);
2471         float force[3],impulse[3],dx[4][3],dv[4][3],oldpos[3];
2472         float dtime=dfra*timestep, time, pa_mass=part->mass, fac /*, fra=sim->psys->cfra*/;
2473         int i, steps=1;
2474
2475         /* maintain angular velocity */
2476         VECCOPY(pa->state.ave,pa->prev_state.ave);
2477         VECCOPY(oldpos,pa->state.co);
2478
2479         if(part->flag & PART_SIZEMASS)
2480                 pa_mass*=pa->size;
2481
2482         switch(part->integrator){
2483                 case PART_INT_EULER:
2484                         steps=1;
2485                         break;
2486                 case PART_INT_MIDPOINT:
2487                         steps=2;
2488                         break;
2489                 case PART_INT_RK4:
2490                         steps=4;
2491                         break;
2492                 case PART_INT_VERLET:
2493                         steps=1;
2494                         break;
2495         }
2496
2497         copy_particle_key(states,&pa->state,1);
2498
2499         for(i=0; i<steps; i++){
2500                 force[0]=force[1]=force[2]=0.0;
2501                 impulse[0]=impulse[1]=impulse[2]=0.0;
2502                 /* add effectors */
2503                 pd_point_from_particle(sim, pa, states+i, &epoint);
2504                 if(part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)
2505                         pdDoEffectors(sim->psys->effectors, sim->colliders, part->effector_weights, &epoint, force, impulse);
2506
2507                 /* calculate air-particle interaction */
2508                 if(part->dragfac!=0.0f){
2509                         fac=-part->dragfac*pa->size*pa->size*len_v3(states[i].vel);
2510                         VECADDFAC(force,force,states[i].vel,fac);
2511                 }
2512
2513                 /* brownian force */
2514                 if(part->brownfac!=0.0){
2515                         force[0]+=(BLI_frand()-0.5f)*part->brownfac;
2516                         force[1]+=(BLI_frand()-0.5f)*part->brownfac;
2517                         force[2]+=(BLI_frand()-0.5f)*part->brownfac;
2518                 }
2519
2520                 /* force to acceleration*/
2521                 mul_v3_fl(force,1.0f/pa_mass);
2522
2523                 /* add global acceleration (gravitation) */
2524                 if(psys_uses_gravity(sim)
2525                         /* normal gravity is too strong for hair so it's disabled by default */
2526                         && (part->type != PART_HAIR || part->effector_weights->flag & EFF_WEIGHT_DO_HAIR)) {
2527                         float gravity[3];
2528                         VECCOPY(gravity, sim->scene->physics_settings.gravity);
2529                         mul_v3_fl(gravity, part->effector_weights->global_gravity);
2530                         VECADD(force,force,gravity);
2531                 }
2532
2533                 /* calculate next state */
2534                 VECADD(states[i].vel,states[i].vel,impulse);
2535
2536                 switch(part->integrator){
2537                         case PART_INT_EULER:
2538                                 VECADDFAC(pa->state.co,states->co,states->vel,dtime);
2539                                 VECADDFAC(pa->state.vel,states->vel,force,dtime);
2540                                 break;
2541                         case PART_INT_MIDPOINT:
2542                                 if(i==0){
2543                                         VECADDFAC(states[1].co,states->co,states->vel,dtime*0.5f);
2544                                         VECADDFAC(states[1].vel,states->vel,force,dtime*0.5f);
2545                                         /*fra=sim->psys->cfra+0.5f*dfra;*/
2546                                 }
2547                                 else{
2548                                         VECADDFAC(pa->state.co,states->co,states[1].vel,dtime);
2549                                         VECADDFAC(pa->state.vel,states->vel,force,dtime);
2550                                 }
2551                                 break;
2552                         case PART_INT_RK4:
2553                                 switch(i){
2554                                         case 0:
2555                                                 VECCOPY(dx[0],states->vel);
2556                                                 mul_v3_fl(dx[0],dtime);
2557                                                 VECCOPY(dv[0],force);
2558                                                 mul_v3_fl(dv[0],dtime);
2559
2560                                                 VECADDFAC(states[1].co,states->co,dx[0],0.5f);
2561                                                 VECADDFAC(states[1].vel,states->vel,dv[0],0.5f);
2562                                                 /*fra=sim->psys->cfra+0.5f*dfra;*/
2563                                                 break;
2564                                         case 1:
2565                                                 VECADDFAC(dx[1],states->vel,dv[0],0.5f);
2566                                                 mul_v3_fl(dx[1],dtime);
2567                                                 VECCOPY(dv[1],force);
2568                                                 mul_v3_fl(dv[1],dtime);
2569
2570                                                 VECADDFAC(states[2].co,states->co,dx[1],0.5f);
2571                                                 VECADDFAC(states[2].vel,states->vel,dv[1],0.5f);
2572                                                 break;
2573                                         case 2:
2574                                                 VECADDFAC(dx[2],states->vel,dv[1],0.5f);
2575                                                 mul_v3_fl(dx[2],dtime);
2576                                                 VECCOPY(dv[2],force);
2577                                                 mul_v3_fl(dv[2],dtime);
2578
2579                                                 VECADD(states[3].co,states->co,dx[2]);
2580                                                 VECADD(states[3].vel,states->vel,dv[2]);
2581                                                 /*fra=cfra;*/
2582                                                 break;
2583                                         case 3:
2584                                                 VECADD(dx[3],states->vel,dv[2]);
2585                                                 mul_v3_fl(dx[3],dtime);
2586                                                 VECCOPY(dv[3],force);
2587                                                 mul_v3_fl(dv[3],dtime);
2588
2589                                                 VECADDFAC(pa->state.co,states->co,dx[0],1.0f/6.0f);
2590                                                 VECADDFAC(pa->state.co,pa->state.co,dx[1],1.0f/3.0f);
2591                                                 VECADDFAC(pa->state.co,pa->state.co,dx[2],1.0f/3.0f);
2592                                                 VECADDFAC(pa->state.co,pa->state.co,dx[3],1.0f/6.0f);
2593
2594                                                 VECADDFAC(pa->state.vel,states->vel,dv[0],1.0f/6.0f);
2595                                                 VECADDFAC(pa->state.vel,pa->state.vel,dv[1],1.0f/3.0f);
2596                                                 VECADDFAC(pa->state.vel,pa->state.vel,dv[2],1.0f/3.0f);
2597                                                 VECADDFAC(pa->state.vel,pa->state.vel,dv[3],1.0f/6.0f);
2598                                 }
2599                                 break;
2600                         case PART_INT_VERLET:   /* Verlet integration */
2601                                 VECADDFAC(pa->state.vel,pa->state.vel,force,dtime);
2602                                 VECADDFAC(pa->state.co,pa->state.co,pa->state.vel,dtime);
2603
2604                                 VECSUB(pa->state.vel,pa->state.co,oldpos);
2605                                 mul_v3_fl(pa->state.vel,1.0f/dtime);
2606                                 break;
2607                 }
2608         }
2609
2610         /* damp affects final velocity */
2611         if(part->dampfac!=0.0)
2612                 mul_v3_fl(pa->state.vel,1.0f-part->dampfac);
2613
2614         VECCOPY(pa->state.ave, states->ave);
2615
2616         /* finally we do guides */
2617         time=(cfra-pa->time)/pa->lifetime;
2618         CLAMP(time,0.0,1.0);
2619
2620         VECCOPY(tkey.co,pa->state.co);
2621         VECCOPY(tkey.vel,pa->state.vel);
2622         tkey.time=pa->state.time;
2623
2624         if(part->type != PART_HAIR) {
2625                 if(do_guides(sim->psys->effectors, &tkey, p, time)) {
2626                         VECCOPY(pa->state.co,tkey.co);
2627                         /* guides don't produce valid velocity */
2628                         VECSUB(pa->state.vel,tkey.co,pa->prev_state.co);
2629                         mul_v3_fl(pa->state.vel,1.0f/dtime);
2630                         pa->state.time=tkey.time;
2631                 }
2632         }
2633 }
2634 static void rotate_particle(ParticleSettings *part, ParticleData *pa, float dfra, float timestep)
2635 {
2636         float rotfac, rot1[4], rot2[4]={1.0,0.0,0.0,0.0}, dtime=dfra*timestep;
2637
2638         if((part->flag & PART_ROT_DYN)==0){
2639                 if(part->avemode==PART_AVE_SPIN){
2640                         float angle;
2641                         float len1 = len_v3(pa->prev_state.vel);
2642                         float len2 = len_v3(pa->state.vel);
2643
2644                         if(len1==0.0f || len2==0.0f)
2645                                 pa->state.ave[0]=pa->state.ave[1]=pa->state.ave[2]=0.0f;
2646                         else{
2647                                 cross_v3_v3v3(pa->state.ave,pa->prev_state.vel,pa->state.vel);
2648                                 normalize_v3(pa->state.ave);
2649                                 angle=dot_v3v3(pa->prev_state.vel,pa->state.vel)/(len1*len2);
2650                                 mul_v3_fl(pa->state.ave,saacos(angle)/dtime);
2651                         }
2652
2653                         axis_angle_to_quat(rot2,pa->state.vel,dtime*part->avefac);
2654                 }
2655         }
2656
2657         rotfac=len_v3(pa->state.ave);
2658         if(rotfac==0.0){ /* unit_qt(in VecRotToQuat) doesn't give unit quat [1,0,0,0]?? */
2659                 rot1[0]=1.0;
2660                 rot1[1]=rot1[2]=rot1[3]=0;
2661         }
2662         else{
2663                 axis_angle_to_quat(rot1,pa->state.ave,rotfac*dtime);
2664         }
2665         mul_qt_qtqt(pa->state.rot,rot1,pa->prev_state.rot);
2666         mul_qt_qtqt(pa->state.rot,rot2,pa->state.rot);
2667
2668         /* keep rotation quat in good health */
2669         normalize_qt(pa->state.rot);
2670 }
2671
2672 /* convert from triangle barycentric weights to quad mean value weights */
2673 static void intersect_dm_quad_weights(float *v1, float *v2, float *v3, float *v4, float *w)
2674 {
2675         float co[3], vert[4][3];
2676
2677         VECCOPY(vert[0], v1);
2678         VECCOPY(vert[1], v2);
2679         VECCOPY(vert[2], v3);
2680         VECCOPY(vert[3], v4);
2681
2682         co[0]= v1[0]*w[0] + v2[0]*w[1] + v3[0]*w[2] + v4[0]*w[3];
2683         co[1]= v1[1]*w[0] + v2[1]*w[1] + v3[1]*w[2] + v4[1]*w[3];
2684         co[2]= v1[2]*w[0] + v2[2]*w[1] + v3[2]*w[2] + v4[2]*w[3];
2685
2686         interp_weights_poly_v3( w,vert, 4, co);
2687 }
2688
2689 /* check intersection with a derivedmesh */
2690 int psys_intersect_dm(Scene *scene, Object *ob, DerivedMesh *dm, float *vert_cos, float *co1, float* co2, float *min_d, int *min_face, float *min_w,
2691                                                   float *face_minmax, float *pa_minmax, float radius, float *ipoint)
2692 {
2693         MFace *mface=0;
2694         MVert *mvert=0;
2695         int i, totface, intersect=0;
2696         float cur_d, cur_uv[2], v1[3], v2[3], v3[3], v4[3], min[3], max[3], p_min[3],p_max[3];
2697         float cur_ipoint[3];
2698
2699         if(dm==0){
2700                 psys_disable_all(ob);
2701
2702                 dm=mesh_get_derived_final(scene, ob, 0);
2703                 if(dm==0)
2704                         dm=mesh_get_derived_deform(scene, ob, 0);
2705
2706                 psys_enable_all(ob);
2707
2708                 if(dm==0)
2709                         return 0;
2710         }
2711
2712
2713
2714         if(pa_minmax==0){
2715                 INIT_MINMAX(p_min,p_max);
2716                 DO_MINMAX(co1,p_min,p_max);
2717                 DO_MINMAX(co2,p_min,p_max);
2718         }
2719         else{
2720                 VECCOPY(p_min,pa_minmax);
2721                 VECCOPY(p_max,pa_minmax+3);
2722         }
2723
2724         totface=dm->getNumFaces(dm);
2725         mface=dm->getFaceDataArray(dm,CD_MFACE);
2726         mvert=dm->getVertDataArray(dm,CD_MVERT);
2727
2728         /* lets intersect the faces */
2729         for(i=0; i<totface; i++,mface++){
2730                 if(vert_cos){
2731                         VECCOPY(v1,vert_cos+3*mface->v1);
2732                         VECCOPY(v2,vert_cos+3*mface->v2);
2733                         VECCOPY(v3,vert_cos+3*mface->v3);
2734                         if(mface->v4)
2735                                 VECCOPY(v4,vert_cos+3*mface->v4)
2736                 }
2737                 else{
2738                         VECCOPY(v1,mvert[mface->v1].co);
2739                         VECCOPY(v2,mvert[mface->v2].co);
2740                         VECCOPY(v3,mvert[mface->v3].co);
2741                         if(mface->v4)
2742                                 VECCOPY(v4,mvert[mface->v4].co)
2743                 }
2744
2745                 if(face_minmax==0){
2746                         INIT_MINMAX(min,max);
2747                         DO_MINMAX(v1,min,max);
2748                         DO_MINMAX(v2,min,max);
2749                         DO_MINMAX(v3,min,max);
2750                         if(mface->v4)
2751                                 DO_MINMAX(v4,min,max)
2752                         if(isect_aabb_aabb_v3(min,max,p_min,p_max)==0)
2753                                 continue;
2754                 }
2755                 else{
2756                         VECCOPY(min, face_minmax+6*i);
2757                         VECCOPY(max, face_minmax+6*i+3);
2758                         if(isect_aabb_aabb_v3(min,max,p_min,p_max)==0)
2759                                 continue;
2760                 }
2761
2762                 if(radius>0.0f){
2763                         if(isect_sweeping_sphere_tri_v3(co1, co2, radius, v2, v3, v1, &cur_d, cur_ipoint)){
2764                                 if(cur_d<*min_d){
2765                                         *min_d=cur_d;
2766                                         VECCOPY(ipoint,cur_ipoint);
2767                                         *min_face=i;
2768                                         intersect=1;
2769                                 }
2770                         }
2771                         if(mface->v4){
2772                                 if(isect_sweeping_sphere_tri_v3(co1, co2, radius, v4, v1, v3, &cur_d, cur_ipoint)){
2773                                         if(cur_d<*min_d){
2774                                                 *min_d=cur_d;
2775                                                 VECCOPY(ipoint,cur_ipoint);
2776                                                 *min_face=i;
2777                                                 intersect=1;
2778                                         }
2779                                 }
2780                         }
2781                 }
2782                 else{
2783                         if(isect_line_tri_v3(co1, co2, v1, v2, v3, &cur_d, cur_uv)){
2784                                 if(cur_d<*min_d){
2785                                         *min_d=cur_d;
2786                                         min_w[0]= 1.0 - cur_uv[0] - cur_uv[1];
2787                                         min_w[1]= cur_uv[0];
2788                                         min_w[2]= cur_uv[1];
2789                                         min_w[3]= 0.0f;
2790                                         if(mface->v4)
2791                                                 intersect_dm_quad_weights(v1, v2, v3, v4, min_w);
2792                                         *min_face=i;
2793                                         intersect=1;
2794                                 }
2795                         }
2796                         if(mface->v4){
2797                                 if(isect_line_tri_v3(co1, co2, v1, v3, v4, &cur_d, cur_uv)){
2798                                         if(cur_d<*min_d){
2799                                                 *min_d=cur_d;
2800                                                 min_w[0]= 1.0 - cur_uv[0] - cur_uv[1];
2801                                                 min_w[1]= 0.0f;
2802                                                 min_w[2]= cur_uv[0];
2803                                                 min_w[3]= cur_uv[1];
2804                                                 intersect_dm_quad_weights(v1, v2, v3, v4, min_w);
2805                                                 *min_face=i;
2806                                                 intersect=1;
2807                                         }
2808                                 }
2809                         }
2810                 }
2811         }
2812         return intersect;
2813 }
2814
2815 void particle_intersect_face(void *userdata, int index, const BVHTreeRay *ray, BVHTreeRayHit *hit)
2816 {
2817         ParticleCollision *col = (ParticleCollision *) userdata;
2818         MFace *face = col->md->mfaces + index;
2819         MVert *x = col->md->x;
2820         MVert *v = col->md->current_v;
2821         float vel[3], co1[3], co2[3], uv[2], ipoint[3], temp[3], t;
2822         float x0[3], x1[3], x2[3], x3[3];
2823         float *t0=x0, *t1=x1, *t2=x2, *t3=(face->v4 ? x3 : NULL);
2824
2825         /* move collision face to start of timestep */
2826         madd_v3_v3v3fl(t0, x[face->v1].co, v[face->v1].co, col->cfra);
2827         madd_v3_v3v3fl(t1, x[face->v2].co, v[face->v2].co, col->cfra);
2828         madd_v3_v3v3fl(t2, x[face->v3].co, v[face->v3].co, col->cfra);
2829         if(t3)
2830                 madd_v3_v3v3fl(t3, x[face->v4].co, v[face->v4].co, col->cfra);
2831
2832         /* calculate average velocity of face */
2833         copy_v3_v3(vel, v[ face->v1 ].co);
2834         add_v3_v3(vel, v[ face->v2 ].co);
2835         add_v3_v3(vel, v[ face->v3 ].co);
2836         mul_v3_fl(vel, 0.33334f*col->dfra);
2837
2838         /* substract face velocity, in other words convert to
2839            a coordinate system where only the particle moves */
2840         madd_v3_v3v3fl(co1, col->co1, vel, -col->f);
2841         sub_v3_v3v3(co2, col->co2, vel);
2842
2843         do
2844         {
2845                 if(ray->radius == 0.0f) {
2846                         if(isect_line_tri_v3(co1, co2, t0, t1, t2, &t, uv)) {
2847                                 if(t >= 0.0f && t < hit->dist/col->ray_len) {
2848                                         hit->dist = col->ray_len * t;
2849                                         hit->index = index;
2850
2851                                         /* calculate normal that's facing the particle */
2852                                         normal_tri_v3( col->nor,t0, t1, t2);
2853                                         VECSUB(temp, co2, co1);
2854                                         if(dot_v3v3(col->nor, temp) > 0.0f)
2855                                                 negate_v3(col->nor);
2856
2857                                         VECCOPY(col->vel,vel);
2858
2859                                         col->hit_ob = col->ob;
2860                                         col->hit_md = col->md;
2861                                 }
2862                         }
2863                 }
2864                 else {
2865                         if(isect_sweeping_sphere_tri_v3(co1, co2, ray->radius, t0, t1, t2, &t, ipoint)) {
2866                                 if(t >=0.0f && t < hit->dist/col->ray_len) {
2867                                         hit->dist = col->ray_len * t;
2868                                         hit->index = index;
2869
2870                                         interp_v3_v3v3(temp, co1, co2, t);
2871
2872                                         VECSUB(col->nor, temp, ipoint);
2873                                         normalize_v3(col->nor);
2874
2875                                         VECCOPY(col->vel,vel);
2876
2877                                         col->hit_ob = col->ob;
2878                                         col->hit_md = col->md;
2879                                 }
2880                         }
2881                 }
2882
2883                 t1 = t2;
2884                 t2 = t3;
2885                 t3 = NULL;
2886
2887         } while(t2);
2888 }
2889 /* Particle - Mesh collision code
2890  * Features:
2891  * - point and swept sphere to mesh surface collisions
2892  * - moving colliders (but not yet rotating or deforming colliders)
2893  * - friction & damping
2894  * - angular momentum <-> linear momentum
2895  * - high accuracy by re-applying particle acceleration after collision
2896  * - behaves relatively well even if limit of 10 collisions per simulation step is exceeded
2897  * Main parts:
2898  * 1. check for all possible deflectors for closest intersection on particle path
2899  * 2. if deflection was found calculate new coordinates or kill the particle
2900  */
2901 static void deflect_particle(ParticleSimulationData *sim, int p, float dfra, float cfra){
2902         Object *ground_ob = NULL;
2903         ParticleSettings *part = sim->psys->part;
2904         ParticleData *pa = sim->psys->particles + p;
2905         ParticleCollision col;
2906         ColliderCache *coll;
2907         BVHTreeRayHit hit;
2908         float ray_dir[3], acc[3];
2909         float radius = ((part->flag & PART_SIZE_DEFL)?pa->size:0.0f), boid_z = 0.0f;
2910         float timestep = psys_get_timestep(sim) * dfra;
2911         float inv_timestep = 1.0f/timestep;
2912         int deflections=0, max_deflections=10;
2913
2914         /* get acceleration (from gravity, forcefields etc. to be re-applied after collision) */
2915         sub_v3_v3v3(acc, pa->state.vel, pa->prev_state.vel);
2916         mul_v3_fl(acc, inv_timestep);
2917
2918         /* set values for first iteration */
2919         copy_v3_v3(col.co1, pa->prev_state.co);
2920         copy_v3_v3(col.co2, pa->state.co);
2921         copy_v3_v3(col.ve1, pa->prev_state.vel);
2922         copy_v3_v3(col.ve2, pa->state.vel);
2923         col.f = 0.0f;