// reflection in e is 2* projection on e - identity, 2P_e - I
// w = 2*(v, e)e  - v;
// Rv=  2*(w, b)/(b,b) b  - w;


	for (i=0; i<3; i++)	bb[i] = e1[i]+e2[i]; // compute bisector and associated quantities
						nsbb=normsq3(bb);
	if (nsbb > tol) { 	// do nothing if bb is numerically zero, identity rotation anyway
						twoovernormsqbb = 2.0/nsbb; 

// rotate view direction unit vector vd

						twov1dotbb = dot3(vd, bb)*twoovernormsqbb;	//  reflect in bb
	for (i=0; i<3; i++) w[i] = twov1dotbb*bb[i]-vd[i];
						twowdote1 = 2.0*dot3(w, e1); //  then in e1
	for (i=0; i<3; i++) vd[i] = twowdote1*e1[i]-w[i];

// then do the same for v2
	
						twov2dotbb = dot3(i2, bb)*twoovernormsqbb;
	for (i=0; i<3; i++) w[i] = twov2dotbb*bb[i]-i2[i];
						twowdote1 = 2.0*dot3(w, e1);
	for (i=0; i<3; i++) i2[i] = twowdote1*e1[i]-w[i];

// That's it!