scalar lambda = 0; if (mag(lambdaDenominator) < VSMALL) { //edge case if (mag(lambdaNominator) < VSMALL) { // we're very close to the face & moving ~parallel, return lambda explicitly if (lambdaDenominator > 0 && lambdaNominator > 0) { if (lambdaDenominator < lambdaNominator) lambda = 1.1; //needs to be > 1.0, particle stays in cell else // travelling along the face lambda = 1.0; // assume we can move to the end position (to machine precision) } else if (lambdaDenominator < 0 && lambdaNominator < 0) { if (lambdaDenominator > lambdaNominator) lambda = 1.1; //needs to be > 1.0, particle stays in cell else // travelling along the face lambda = 1.0; // assume we can move to the end position (to machine precision) } else { //moving away from face, needs to be < 0 lambda = -1.0; } } else { lamdba = 1.1; // not passing through this face } } else { lambda = lambdaNominator/lambdaDenominator; } return lambda;