void starblebounce(int j, int k, float uf) // hard spheres collision response // September 2006 // feb 2018... { vector disp; float dist, flt; vector cl; //collision line star s1=stars[j],s2=stars[k]; vector v1=s1.v, v2=s2.v; vector v1icl,v2icl; // collision line initial velocities float v1iclm,v2iclm; // collision line initial velocities magnitudes vector v1fcl,v2fcl; // collision line final velocities float v1fclm,v2fclm; // collision line final velocities magnitudes vector v1ici,v2ici; // collision transverse initial velocities vector v1fci,v2fci; // collision transverse final velocities float m1=s1.m, m2=s2.m; disp = vectdiff(s2.r,s1.r); dist = vabs(disp); if(dist==0.0) return; cl = vecttimes( 1.0/dist , disp ); // collision line unit vector from 1 to 2 v1iclm = vectdot( v1 , cl ); v2iclm = vectdot( v2 , cl ); if (v1iclm-v2iclm<0.0) return; //not colliding v1icl = vecttimes( v1iclm , cl ); v2icl = vecttimes( v2iclm , cl ); v1ici = vectdiff( v1 , v1icl ); v2ici = vectdiff( v2 , v2icl ); v1fclm = ( (m1-m2) * v1iclm + 2.0 * m2 * v2iclm ) / ( m1 + m2 ) ; v2fclm = ( (m2-m1) * v2iclm + 2.0 * m1 * v1iclm ) / ( m1 + m2 ) ; v1fcl = vecttimes( v1fclm , cl ); v2fcl = vecttimes( v2fclm , cl ); if ( uf!=0.0 ) { vector cp; //collision point position float rad1=s1.rad, rad2=s2.rad; vector v1cp,v2cp; // collision-transverse velocities at the contact point vector rv1cp; // contact point relative collision-transverse velocity vector rv1cpuv,rv2cpuv; // '' unit vectors vector w1i=s1.w, w2i=s2.w, w1f, w2f; float cldv1,cldv2; //collision velocity delta for friction normal float i1,i2; // moments of inertia int overfric1=0,overfric2=0; float td1=1.0, td2=1.0; vector cpr1,cpr2; //contact point from ball centers cp = vectsum( vecttimes( 0.5*rad1 - 0.5*rad2 + 0.5*dist , cl ), s1.r ); // i1 = ballmoment(m1,rad1); i2 = ballmoment(m2,rad2); cpr1 = vectdiff(cp,s1.r); cpr2 = vectdiff(cp,s2.r); cldv1 = vabs(vectdiff(v1fcl,v1icl)); cldv2 = vabs(vectdiff(v2fcl,v2icl)); v1cp = vectsum( v1ici, vectcross( w1i, cpr1 ) ); v2cp = vectsum( v2ici, vectcross( w2i, cpr2 ) ); rv1cp = vectdiff( v1cp, v2cp ); //rv2cp = vectdiff( v2cp, v1cp ); rv1cpuv = unitvect(rv1cp); rv2cpuv = vecttimes(-1.0,rv1cpuv); v1fci = vectsum( v1ici, vecttimes( -uf*cldv1, rv1cpuv ) ); v2fci = vectsum( v2ici, vecttimes( -uf*cldv2, rv2cpuv ) ); w1f = vectsum( w1i, vectcross( cpr1, vecttimes( -uf*cldv1*m1/i1, rv1cpuv ) ) ); w2f = vectsum( w2i, vectcross( cpr2, vecttimes( -uf*cldv2*m2/i2, rv2cpuv ) ) ); if( vectdot( vectsum( v1fci, vectcross( w1f, cpr1) ), rv1cpuv) - vectdot( v2cp, rv1cpuv) < 0.0 ) { overfric1=1; td1 = vectdot(rv1cp,rv1cpuv) / ( uf * cldv1 * (1.0 + m1*vabs2(cpr1)/i1) ); } if( vectdot( vectsum( v2fci, vectcross( w2f, cpr2) ), rv2cpuv) - vectdot( v1cp, rv2cpuv) < 0.0 ) { overfric2=1; td2 = vectdot(rv1cp,rv1cpuv) / ( uf * cldv2 * (1.0 + m2*vabs2(cpr2)/i2) ); } if( overfric1 || overfric2 ) { //if(td1>1.0 || td2>1.0){gotoxy(20,1); printf("%4f %4f",td1,td2);}//floating point curiosity; shouldn't happen, but it can if ( td2m1*vabs2(v1ici)+i1*vabs2(w1i)+m2*vabs2(v2ici)+i2*vabs2(w2i) ) { float tlke,a,b; b = 2.0*m1*uf*cldv1*td1*vectdot(v1ici,rv1cpuv) + 2.0*vabs(cpr1)*uf*cldv1*m1*td1*vectdot(w1i,rv1cpuv) 2.0*m2*uf*cldv2*td1*vectdot(v2ici,rv2cpuv) + 2.0*vabs(cpr2)*uf*cldv2*m2*td1*vectdot(w2i,rv2cpuv); a = m1*pow(uf*cldv1*td1,2) + i1*pow(vabs(cpr1)*uf*cldv1*m1/i1*td1,2) + m2*pow(uf*cldv2*td2,2) + i2*pow(vabs(cpr2)*uf*cldv2*m2/i2*td2,2); // td1,td1? // redo // ... tlke = b/a; if(tlke<0.0)tlke=0.0; v1fci = vectsum( v1ici, vecttimes( -uf*cldv1*td1*tlke, rv1cpuv ) ); v2fci = vectsum( v2ici, vecttimes( -uf*cldv2*td1*tlke, rv2cpuv ) ); w1f = vectsum( w1i, vectcross( cpr1, vecttimes( -uf*cldv1*m1/i1*td1*tlke, rv1cpuv ) ) ); w2f = vectsum( w2i, vectcross( cpr2, vecttimes( -uf*cldv2*m2/i2*td1*tlke, rv2cpuv ) ) ); } // stars[j].w=w1f; stars[k].w=w2f; stars[j].v=vectsum(v1fcl,v1fci); stars[k].v=vectsum(v2fcl,v2fci); } else { stars[j].v=vectsum(v1fcl,v1ici); stars[k].v=vectsum(v2fcl,v2ici); } }