Examples with RationalNumber cbit.vcell.matrix.RationalNumber used on opensource projects

Example 6 with RationalNumber

use of cbit.vcell.matrix.RationalNumber in project vcell by virtualcell.

the class SVDTest method getRationalMatrixFast.

/**
 * Insert the method's description here.
 * Creation date: (5/13/2003 11:06:29 AM)
 * @return Jama.Matrix
 */
public static RationalMatrixFast getRationalMatrixFast(SimpleMatrix matrix) {
    RationalMatrixFast rMatrix = new RationalMatrixFast(matrix.getNumRows(), matrix.getNumCols());
    for (int i = 0; i < rMatrix.getNumRows(); i++) {
        for (int j = 0; j < rMatrix.getNumCols(); j++) {
            RationalNumber r = RationalNumber.getApproximateFraction(matrix.get_elem(i, j));
            rMatrix.set_elem(i, j, r.getNumBigInteger().longValue(), r.getDenBigInteger().longValue());
        }
    }
    return rMatrix;
}

Example 7 with RationalNumber

use of cbit.vcell.matrix.RationalNumber in project vcell by virtualcell.

the class ModelUnitConverter method convertExprWithUnitFactors.

private static void convertExprWithUnitFactors(SymbolTable oldSymbolTable, SymbolTable newSymbolTable, VCUnitDefinition oldExprUnit, VCUnitDefinition newExprUnit, Expression expr) throws ExpressionException {
    if (expr == null) {
        return;
    }
    String[] symbols = expr.getSymbols();
    if (symbols != null) {
        for (String s : symbols) {
            SymbolTableEntry boundSTE = expr.getSymbolBinding(s);
            SymbolTableEntry newSTE = newSymbolTable.getEntry(s);
            SymbolTableEntry oldSTE = oldSymbolTable.getEntry(s);
            if (boundSTE == null) {
                System.err.println("symbol '" + s + "' in expression " + expr.infix() + " was not bound to the model");
                continue;
            }
            if (oldSTE == null) {
                System.err.println("symbol '" + s + "' in expression " + expr.infix() + " was not found in the new symbol table");
                continue;
            }
            if (newSTE == null) {
                System.err.println("symbol '" + s + "' in expression " + expr.infix() + " was not bound to the old symbol table");
                continue;
            }
            if (newSTE != boundSTE) {
                System.err.println("symbol '" + s + "' in expression " + expr.infix() + " is not bound properly (binding doesn't match new symbol table)");
                continue;
            }
            if (oldSTE.getUnitDefinition() == null) {
                System.err.println("symbol '" + s + "' in expression " + expr.infix() + " is has a null unit in old model, can't convert");
                continue;
            }
            if (newSTE.getUnitDefinition() == null) {
                System.err.println("symbol '" + s + "' in expression " + expr.infix() + " is has a null unit in new model, can't convert");
                continue;
            }
            if (oldSTE.getUnitDefinition().isTBD()) {
                continue;
            }
            if (newSTE.getUnitDefinition().isTBD()) {
                continue;
            }
            VCUnitDefinition oldToNewConversionUnit = oldSTE.getUnitDefinition().divideBy(newSTE.getUnitDefinition());
            RationalNumber conversionFactor = oldToNewConversionUnit.getDimensionlessScale();
            if (conversionFactor.doubleValue() != 1.0) {
                Expression spcSymbol = new Expression(newSTE, newSTE.getNameScope());
                expr.substituteInPlace(spcSymbol, Expression.mult(spcSymbol, new Expression(conversionFactor)));
            }
        }
    }
    if (oldExprUnit == null || oldExprUnit.isTBD() || newExprUnit == null || newExprUnit.isTBD()) {
        Expression flattened = expr.flatten();
        Expression origExp = new Expression(expr);
        expr.substituteInPlace(origExp, flattened);
        return;
    }
    VCUnitDefinition oldToNewConversionUnit = newExprUnit.divideBy(oldExprUnit);
    RationalNumber conversionFactor = oldToNewConversionUnit.getDimensionlessScale();
    if (conversionFactor.doubleValue() != 1.0) {
        Expression oldExp = new Expression(expr);
        expr.substituteInPlace(oldExp, Expression.mult(oldExp, new Expression(conversionFactor)));
    }
    Expression flattened = expr.flatten();
    Expression origExp = new Expression(expr);
    expr.substituteInPlace(origExp, flattened);
}

Example 8 with RationalNumber

use of cbit.vcell.matrix.RationalNumber in project vcell by virtualcell.

the class ASTPowerTerm method toSymbol.

public String toSymbol(RationalNumber power, UnitTextFormat format) {
    if (jjtGetChild(0) instanceof ASTIdNode || jjtGetChild(0) instanceof ASTIntegerBaseNode) {
        if (jjtGetNumChildren() == 1) {
            if (power.equals(RationalNumber.ONE)) {
                return jjtGetChild(0).toSymbol(power, format);
            } else {
                switch(format) {
                    case plain:
                        {
                            return jjtGetChild(0).toSymbol(power, format) + power.infix();
                        }
                    case unicode:
                        {
                            String powerString = ASTRationalNumberExponent.getUnicodeExponent(power);
                            return jjtGetChild(0).toSymbol(power, format) + powerString;
                        }
                    default:
                        {
                            throw new RuntimeException("format " + format.name() + " not supported by unit symbol");
                        }
                }
            }
        } else {
            return jjtGetChild(0).toSymbol(power, format) + jjtGetChild(1).toSymbol(power, format);
        }
    } else if (jjtGetNumChildren() == 1) {
        return jjtGetChild(0).toSymbol(power, format);
    } else if (jjtGetNumChildren() == 2) {
        if (jjtGetChild(1) instanceof ASTNegative) {
            ASTNegative negNode = (ASTNegative) jjtGetChild(1);
            RationalNumber exponent = negNode.getRationalNumber();
            return jjtGetChild(0).toSymbol(power.mult(exponent), format);
        } else if (jjtGetChild(1) instanceof ASTRationalNumberExponent) {
            ASTRationalNumberExponent rationalNumberExponent = (ASTRationalNumberExponent) jjtGetChild(1);
            RationalNumber exponent = rationalNumberExponent.value;
            return jjtGetChild(0).toSymbol(power.mult(exponent), format);
        } else {
            throw new RuntimeException("unexpected second child " + jjtGetChild(1).getClass().getName());
        }
    } else {
        throw new RuntimeException("unexpected unit format");
    }
}

Example 9 with RationalNumber

use of cbit.vcell.matrix.RationalNumber in project vcell by virtualcell.

the class StructureAnalyzer method refreshTotalDependancies.

/**
 * This method was created by a SmartGuide.
 * @param b cbit.vcell.math.Matrix
 * @param vars java.lang.String[]
 */
static StructureAnalyzer.Dependency[] refreshTotalDependancies(RationalMatrix nullSpaceMatrix, SpeciesContextMapping[] speciesContextMappings, MathMapping_4_8 argMathMapping, boolean bFast) throws Exception {
    // System.out.println("StructureAnalyzer.refreshTotalDependancies()");
    SimulationContext argSimContext = argMathMapping.getSimulationContext();
    if (nullSpaceMatrix == null) {
        // System.out.println("the matrix has full rank, there are no dependencies");
        return new StructureAnalyzer.Dependency[0];
    }
    if (speciesContextMappings.length != nullSpaceMatrix.getNumCols()) {
        // throw new Exception("varName array not same dimension as b matrix");
        System.out.println("varName array not same dimension as b matrix");
        nullSpaceMatrix.show();
        for (int i = 0; i < speciesContextMappings.length; i++) {
            System.out.println("scm[" + i + "] is " + speciesContextMappings[i].getSpeciesContext().getName());
        }
    }
    // System.out.println("there are "+nullSpaceMatrix.rows+" dependencies");
    Vector<Dependency> dependencyList = new Vector<Dependency>();
    for (int i = 0; i < nullSpaceMatrix.getNumRows(); i++) {
        // 
        // find first variable
        // 
        Expression exp = null;
        Expression constantExp = null;
        String constantName = null;
        SpeciesContextMapping firstSCM = null;
        boolean bFirst = true;
        for (int j = 0; j < nullSpaceMatrix.getNumCols(); j++) {
            RationalNumber coeff = nullSpaceMatrix.get_elem(i, j);
            if (coeff.doubleValue() != 0.0) {
                if (bFirst) {
                    if (coeff.doubleValue() != 1.0) {
                        System.out.println("i=" + i + " j=" + j);
                        nullSpaceMatrix.show();
                        throw new Exception("expecting a coefficient of 1.0, instead coeff = " + coeff.infix());
                    }
                    firstSCM = speciesContextMappings[j];
                    // 
                    // first term of dependancy expression   ("K_CalciumCyt")
                    // 
                    SpeciesContext firstSC = firstSCM.getSpeciesContext();
                    constantName = "K_" + firstSC.getName() + "_total";
                    exp = new Expression(constantName);
                    // 
                    // first term of K expression
                    // 
                    StructureMapping firstSM = argSimContext.getGeometryContext().getStructureMapping(firstSC.getStructure());
                    SpeciesContextSpec firstSCS = argSimContext.getReactionContext().getSpeciesContextSpec(firstSC);
                    SymbolTableEntry scSTE = null;
                    if (bFast) {
                        scSTE = firstSCS.getSpeciesContext();
                    } else {
                        scSTE = firstSCS.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
                    }
                    constantExp = Expression.mult(new Expression(coeff), firstSM.getNormalizedConcentrationCorrection(argSimContext, argMathMapping), new Expression(scSTE, argMathMapping.getNameScope()));
                    bFirst = false;
                } else {
                    // 
                    // add term to dependancy expression     (" - 2*IP3Cyt ")
                    // 
                    SpeciesContextMapping scm = speciesContextMappings[j];
                    SpeciesContext sc = scm.getSpeciesContext();
                    StructureMapping sm = argSimContext.getGeometryContext().getStructureMapping(sc.getStructure());
                    SpeciesContextSpec scs = argSimContext.getReactionContext().getSpeciesContextSpec(sc);
                    exp = Expression.add(exp, Expression.negate(Expression.mult(new Expression(coeff), sm.getNormalizedConcentrationCorrection(argSimContext, argMathMapping), new Expression(sc, argMathMapping.getNameScope()))));
                    // 
                    // add term to K expression
                    // 
                    SymbolTableEntry scSTE = null;
                    if (bFast) {
                        scSTE = sc;
                    } else {
                        scSTE = scs.getParameterFromRole(SpeciesContextSpec.ROLE_InitialConcentration);
                    }
                    constantExp = Expression.add(constantExp, Expression.mult(new Expression(coeff), sm.getNormalizedConcentrationCorrection(argSimContext, argMathMapping), new Expression(scSTE, argMathMapping.getNameScope())));
                }
            }
        }
        if (firstSCM != null) {
            // 
            // store totalMass parameter (e.g. K_xyz_total = xyz_init + wzy_init)
            // 
            // MathMapping.MathMappingParameter totalMassParameter = mathMapping.addMathMappingParameter(constantName,constantExp.flatten(),MathMapping.PARAMETER_ROLE_TOTALMASS,cbit.vcell.units.VCUnitDefinition.UNIT_uM);
            // 
            // store dependency parameter (e.g. xyz = K_xyz_total - wzy)
            // 
            StructureMapping sm = argSimContext.getGeometryContext().getStructureMapping(firstSCM.getSpeciesContext().getStructure());
            exp = Expression.mult(exp, Expression.invert(sm.getNormalizedConcentrationCorrection(argSimContext, argMathMapping)));
            exp = exp.flatten();
            // exp.bindExpression(mathMapping_temp);
            // firstSCM.setDependencyExpression(exp);
            Dependency dependency = new Dependency();
            dependency.dependencyExpression = exp;
            dependency.speciesContextMapping = firstSCM;
            dependency.invariantSymbolName = constantName;
            dependency.conservedMoietyExpression = constantExp.flatten();
            dependencyList.add(dependency);
        }
    }
    return (StructureAnalyzer.Dependency[]) BeanUtils.getArray(dependencyList, StructureAnalyzer.Dependency.class);
}

Example 10 with RationalNumber

use of cbit.vcell.matrix.RationalNumber in project vcell by virtualcell.

the class RayCaster method rayCastXYZ_rational.

// 
// compute ray intersection in X-Z plane with ray cast from y = -Infinity to +Infinity
// hitList is indexed as (i+numX*k)
// 
public static RayCastResultsRational rayCastXYZ_rational(SurfaceCollection surfaceCollection, RationalNumber[] samplesX, RationalNumber[] samplesY, RationalNumber[] samplesZ) {
    int numX = samplesX.length;
    int numY = samplesY.length;
    int numZ = samplesZ.length;
    long t1 = System.currentTimeMillis();
    HitListRational[] hitListsXY = new HitListRational[numX * numY];
    for (int j = 0; j < numY; j++) {
        for (int i = 0; i < numX; i++) {
            // samplesX[i], samplesY[j]
            hitListsXY[i + numX * j] = new HitListRational();
        }
    }
    HitListRational[] hitListsXZ = new HitListRational[numX * numZ];
    for (int k = 0; k < numZ; k++) {
        for (int i = 0; i < numX; i++) {
            // samplesX[i], samplesZ[k]
            hitListsXZ[i + numX * k] = new HitListRational();
        }
    }
    HitListRational[] hitListsYZ = new HitListRational[numY * numZ];
    for (int k = 0; k < numZ; k++) {
        for (int j = 0; j < numY; j++) {
            // samplesY[j], samplesZ[k]
            hitListsYZ[j + numY * k] = new HitListRational();
        }
    }
    // 
    // here we expect either triangles or nonplanar quads which will be broken up into two triangles.
    // then process using barycentric-like coordinates.
    // 
    int hitCount = 0;
    cbit.vcell.geometry.surface.Triangle[] triangles = new cbit.vcell.geometry.surface.Triangle[2];
    for (int j = 0; j < surfaceCollection.getSurfaceCount(); j++) {
        Surface surface = surfaceCollection.getSurfaces(j);
        verifySurface(surface);
        for (int k = 0; k < surface.getPolygonCount(); k++) {
            Polygon polygon = surface.getPolygons(k);
            Node[] nodes = polygon.getNodes();
            double p0x = nodes[0].getX();
            double p0y = nodes[0].getY();
            double p0z = nodes[0].getZ();
            double p1x = nodes[1].getX();
            double p1y = nodes[1].getY();
            double p1z = nodes[1].getZ();
            double p2x = nodes[2].getX();
            double p2y = nodes[2].getY();
            double p2z = nodes[2].getZ();
            double minX = Math.min(p0x, Math.min(p1x, p2x));
            double maxX = Math.max(p0x, Math.max(p1x, p2x));
            double minY = Math.min(p0y, Math.min(p1y, p2y));
            double maxY = Math.max(p0y, Math.max(p1y, p2y));
            double minZ = Math.min(p0z, Math.min(p1z, p2z));
            double maxZ = Math.max(p0z, Math.max(p1z, p2z));
            if (polygon.getNodeCount() == 4) {
                double p3x = nodes[3].getX();
                double p3y = nodes[3].getY();
                double p3z = nodes[3].getZ();
                minX = Math.min(minX, p3x);
                maxX = Math.max(maxX, p3x);
                minY = Math.min(minY, p3y);
                maxY = Math.max(maxY, p3y);
                minZ = Math.min(minZ, p3z);
                maxZ = Math.max(maxZ, p3z);
            }
            {
                double epsilon = 1e-8;
                minX -= epsilon;
                maxX += epsilon;
                minY -= epsilon;
                maxY += epsilon;
                minZ -= epsilon;
                maxZ += epsilon;
            }
            // precompute ray indices that are within the bounding box of the quad.
            int startI = numX;
            int endI = -1;
            for (int ii = 0; ii < numX; ii++) {
                double sampleX = samplesX[ii].doubleValue();
                if (sampleX >= minX && sampleX <= maxX) {
                    startI = Math.min(startI, ii);
                    endI = Math.max(endI, ii);
                }
            }
            int startJ = numY;
            int endJ = -1;
            for (int jj = 0; jj < numY; jj++) {
                double sampleY = samplesY[jj].doubleValue();
                if (sampleY >= minY && sampleY <= maxY) {
                    startJ = Math.min(startJ, jj);
                    endJ = Math.max(endJ, jj);
                }
            }
            int startK = numZ;
            int endK = -1;
            for (int kk = 0; kk < numZ; kk++) {
                double sampleZ = samplesZ[kk].doubleValue();
                if (sampleZ >= minZ && sampleZ <= maxZ) {
                    startK = Math.min(startK, kk);
                    endK = Math.max(endK, kk);
                }
            }
            // 
            // convert to quads to triangles if necessary
            // 
            int numTriangles;
            if (polygon.getNodeCount() == 3) {
                numTriangles = 1;
                triangles[0] = new cbit.vcell.geometry.surface.Triangle(polygon.getNodes(0), polygon.getNodes(1), polygon.getNodes(2));
            } else if (polygon.getNodeCount() == 4) {
                numTriangles = 2;
                triangles[0] = new cbit.vcell.geometry.surface.Triangle(polygon.getNodes(0), polygon.getNodes(1), polygon.getNodes(2));
                triangles[1] = new cbit.vcell.geometry.surface.Triangle(polygon.getNodes(0), polygon.getNodes(2), polygon.getNodes(3));
            } else {
                throw new RuntimeException("polygons with " + polygon.getNodeCount() + " edges are supported");
            }
            for (int triIndex = 0; triIndex < numTriangles; triIndex++) {
                cbit.vcell.geometry.surface.Triangle triangle = triangles[triIndex];
                Node node0 = triangle.getNodes()[0];
                RationalNumber ax = RationalNumber.getApproximateFraction(node0.getX());
                RationalNumber ay = RationalNumber.getApproximateFraction(node0.getY());
                RationalNumber az = RationalNumber.getApproximateFraction(node0.getZ());
                Node node1 = triangle.getNodes()[1];
                RationalNumber bx = RationalNumber.getApproximateFraction(node1.getX());
                RationalNumber by = RationalNumber.getApproximateFraction(node1.getY());
                RationalNumber bz = RationalNumber.getApproximateFraction(node1.getZ());
                Node node2 = triangle.getNodes()[2];
                RationalNumber cx = RationalNumber.getApproximateFraction(node2.getX());
                RationalNumber cy = RationalNumber.getApproximateFraction(node2.getY());
                RationalNumber cz = RationalNumber.getApproximateFraction(node2.getZ());
                // bx-ax;
                RationalNumber v1x = bx.sub(ax);
                // by-ay;
                RationalNumber v1y = by.sub(ay);
                // bz-az;
                RationalNumber v1z = bz.sub(az);
                // cx-ax;
                RationalNumber v0x = cx.sub(ax);
                // cy-ay;
                RationalNumber v0y = cy.sub(ay);
                // cz-az;
                RationalNumber v0z = cz.sub(az);
                // v1y*v0z - v1z*v0y;
                RationalNumber nx = v1y.mult(v0z).sub(v1z.mult(v0y));
                // -(v1x*v0z - v1z*v0x);
                RationalNumber ny = v1x.mult(v0z).sub(v1z.mult(v0x)).minus();
                // v1x*v0y - v1y*v0x;
                RationalNumber nz = v1x.mult(v0y).sub(v1y.mult(v0x));
                // here we probably don't have to normalize the vector ... so don't have to do this.
                // Math.sqrt(nx*nx+ny*ny+nz*nz);
                RationalNumber normalLength = RationalNumber.getApproximateFraction(Math.sqrt(nx.mult(nx).add(ny.mult(ny)).add(nz.mult(nz)).doubleValue()));
                nx = nx.div(normalLength);
                ny = ny.div(normalLength);
                nz = nz.div(normalLength);
                // SAMPLE IN XY
                {
                    // v1x*v1x + v1y*v1y;
                    RationalNumber dot11_xy = v1x.mult(v1x).add(v1y.mult(v1y));
                    // v0x*v1x + v0y*v1y;
                    RationalNumber dot01_xy = v0x.mult(v1x).add(v0y.mult(v1y));
                    // v0x*v0x + v0y*v0y;
                    RationalNumber dot00_xy = v0x.mult(v0x).add(v0y.mult(v0y));
                    // get normal in z direction (A,B,C) ... z component of (B-A)X(C-A)
                    // boolean entering = (nz < 0);
                    // 1.0 / (dot00_xy*dot11_xy - dot01_xy*dot01_xy);
                    RationalNumber invDenom = dot00_xy.mult(dot11_xy).sub(dot01_xy.mult(dot01_xy)).inverse();
                    for (int rayIndexI = startI; rayIndexI <= endI; rayIndexI++) {
                        for (int rayIndexJ = startJ; rayIndexJ <= endJ; rayIndexJ++) {
                            HitListRational hitList = hitListsXY[rayIndexI + rayIndexJ * numX];
                            // samplesX[rayIndexI]-ax;
                            RationalNumber vpx = samplesX[rayIndexI].sub(ax);
                            // samplesY[rayIndexJ]-ay;
                            RationalNumber vpy = samplesY[rayIndexJ].sub(ay);
                            // v0x*vpx + v0y*vpy;
                            RationalNumber dot0p_xy = v0x.mult(vpx).add(v0y.mult(vpy));
                            // v1x*vpx + v1y*vpy;
                            RationalNumber dot1p_xy = v1x.mult(vpx).add(v1y.mult(vpy));
                            // (dot11_xy * dot0p_xy - dot01_xy * dot1p_xy) * invDenom;
                            RationalNumber u = (dot11_xy.mult(dot0p_xy).sub(dot01_xy.mult(dot1p_xy))).mult(invDenom);
                            // (dot00_xy * dot1p_xy - dot01_xy * dot0p_xy) * invDenom;
                            RationalNumber v = (dot00_xy.mult(dot1p_xy).sub(dot01_xy.mult(dot0p_xy))).mult(invDenom);
                            // if ((u > -epsilon) && (v > -epsilon) && (u + v < 1+epsilon)){
                            if (u.ge(RationalNumber.ZERO) && v.ge(RationalNumber.ZERO) && u.add(v).le(RationalNumber.ONE)) {
                                // az + (u*v0z) + (v*v1z);
                                RationalNumber v2z = az.add(u.mult(v0z)).add(v.mult(v1z));
                                if (nz.isZero()) {
                                    System.err.println("normal is zero for triangle index " + triIndex + " in XY plane: vertices [(" + ax + "," + ay + "," + az + "),(" + bx + "," + by + "," + bz + "),(" + cx + "," + cy + "," + cz + ")] in plane YZ");
                                }
                                RationalNumber centroidZ = az.add(bz).add(cz).div(new RationalNumber(3));
                                String debugMessage = "u,v=(" + u + "," + v + "),nodes=(" + node0.getGlobalIndex() + "," + node1.getGlobalIndex() + "," + node2.getGlobalIndex() + ")";
                                hitList.addHitEvent(new HitEventRational(surface, polygon, nz, v2z, centroidZ, debugMessage));
                                // System.out.println("hit xy @ x="+samplesX[rayIndexI]+", y="+samplesY[rayIndexJ]+" ... ray hit at z="+v2z);
                                hitCount++;
                            }
                        }
                    }
                }
                // SAMPLE IN XZ
                {
                    // v1x*v1x + v1z*v1z;
                    RationalNumber dot11_xz = v1x.mult(v1x).add(v1z.mult(v1z));
                    // v0x*v1x + v0z*v1z;
                    RationalNumber dot01_xz = v0x.mult(v1x).add(v0z.mult(v1z));
                    // v0x*v0x + v0z*v0z;
                    RationalNumber dot00_xz = v0x.mult(v0x).add(v0z.mult(v0z));
                    // get normal in y direction (A,B,C) ... y component of (B-A)X(C-A)
                    // boolean entering = (ny < 0);
                    // 1.0 / (dot00_xz*dot11_xz - dot01_xz*dot01_xz);
                    RationalNumber invDenom = dot00_xz.mult(dot11_xz).sub(dot01_xz.mult(dot01_xz)).inverse();
                    for (int rayIndexI = startI; rayIndexI <= endI; rayIndexI++) {
                        for (int rayIndexK = startK; rayIndexK <= endK; rayIndexK++) {
                            HitListRational hitList = hitListsXZ[rayIndexI + rayIndexK * numX];
                            // samplesX[rayIndexI]-ax;
                            RationalNumber vpx = samplesX[rayIndexI].sub(ax);
                            // samplesZ[rayIndexK]-az;
                            RationalNumber vpz = samplesZ[rayIndexK].sub(az);
                            // v0x*vpx + v0z*vpz;
                            RationalNumber dot0p_xz = v0x.mult(vpx).add(v0z.mult(vpz));
                            // v1x*vpx + v1z*vpz;
                            RationalNumber dot1p_xz = v1x.mult(vpx).add(v1z.mult(vpz));
                            // (dot11_xz * dot0p_xz - dot01_xz * dot1p_xz) * invDenom;
                            RationalNumber u = dot11_xz.mult(dot0p_xz).sub(dot01_xz.mult(dot1p_xz)).mult(invDenom);
                            // (dot00_xz * dot1p_xz - dot01_xz * dot0p_xz) * invDenom;
                            RationalNumber v = dot00_xz.mult(dot1p_xz).sub(dot01_xz.mult(dot0p_xz)).mult(invDenom);
                            // if ((u > -epsilon) && (v > -epsilon) && (u + v < 1+epsilon)){
                            if (u.ge(RationalNumber.ZERO) && v.ge(RationalNumber.ZERO) && u.add(v).le(RationalNumber.ONE)) {
                                // ay + (u*v0y) + (v*v1y);
                                RationalNumber v2y = ay.add(u.mult(v0y)).add(v.mult(v1y));
                                if (ny.isZero()) {
                                    System.err.println("normal is zero for triangle index " + triIndex + " in XZ plane (" + rayIndexI + "," + rayIndexK + "): vertices [(" + ax + "," + ay + "," + az + "),(" + bx + "," + by + "," + bz + "),(" + cx + "," + cy + "," + cz + ")] in plane YZ");
                                }
                                RationalNumber centroidY = ay.add(by).add(cy).div(new RationalNumber(3));
                                String debugMessage = "u,v=(" + u + "," + v + "),nodes=(" + node0.getGlobalIndex() + "," + node1.getGlobalIndex() + "," + node2.getGlobalIndex() + ")";
                                hitList.addHitEvent(new HitEventRational(surface, polygon, ny, v2y, centroidY, debugMessage));
                                // System.out.println("hit xz @ x="+samplesX[rayIndexI]+", z="+samplesZ[rayIndexK]+" ... ray hit at y="+v2y);
                                hitCount++;
                            }
                        }
                    }
                }
                // SAMPLE IN YZ
                {
                    // v1y*v1y + v1z*v1z;
                    RationalNumber dot11_yz = v1y.mult(v1y).add(v1z.mult(v1z));
                    // v0y*v1y + v0z*v1z;
                    RationalNumber dot01_yz = v0y.mult(v1y).add(v0z.mult(v1z));
                    // v0y*v0y + v0z*v0z;
                    RationalNumber dot00_yz = v0y.mult(v0y).add(v0z.mult(v0z));
                    // get normal in y direction (A,B,C) ... y component of (B-A)X(C-A)
                    // boolean entering = (nx < 0);
                    // 1.0 / (dot00_yz*dot11_yz - dot01_yz*dot01_yz);
                    RationalNumber invDenom = dot00_yz.mult(dot11_yz).sub(dot01_yz.mult(dot01_yz)).inverse();
                    for (int rayIndexJ = startJ; rayIndexJ <= endJ; rayIndexJ++) {
                        for (int rayIndexK = startK; rayIndexK <= endK; rayIndexK++) {
                            HitListRational hitList = hitListsYZ[rayIndexJ + rayIndexK * numY];
                            // samplesY[rayIndexJ]-ay;
                            RationalNumber vpy = samplesY[rayIndexJ].sub(ay);
                            // samplesZ[rayIndexK]-az;
                            RationalNumber vpz = samplesZ[rayIndexK].sub(az);
                            // v0y*vpy + v0z*vpz;
                            RationalNumber dot0p_yz = v0y.mult(vpy).add(v0z.mult(vpz));
                            // v1y*vpy + v1z*vpz;
                            RationalNumber dot1p_yz = v1y.mult(vpy).add(v1z.mult(vpz));
                            // (dot11_yz * dot0p_yz - dot01_yz * dot1p_yz) * invDenom;
                            RationalNumber u = dot11_yz.mult(dot0p_yz).sub(dot01_yz.mult(dot1p_yz)).mult(invDenom);
                            // (dot00_yz * dot1p_yz - dot01_yz * dot0p_yz) * invDenom;
                            RationalNumber v = dot00_yz.mult(dot1p_yz).sub(dot01_yz.mult(dot0p_yz)).mult(invDenom);
                            // if ((u > -epsilon) && (v > -epsilon) && (u + v < 1+epsilon)){
                            if (u.ge(RationalNumber.ZERO) && v.ge(RationalNumber.ZERO) && u.add(v).le(RationalNumber.ONE)) {
                                // ax + (u*v0x) + (v*v1x);
                                RationalNumber v2x = ax.add(u.mult(v0x)).add(v.mult(v1x));
                                if (ny.isZero()) {
                                    System.err.println("normal is zero for triangle index " + triIndex + " in YZ plane: vertices [(" + ax + "," + ay + "," + az + "),(" + bx + "," + by + "," + bz + "),(" + cx + "," + cy + "," + cz + ")] in plane YZ");
                                }
                                RationalNumber centroidX = ax.add(bx).add(cx).div(new RationalNumber(3));
                                String debugMessage = "u,v=(" + u + "," + v + "),nodes=(" + node0.getGlobalIndex() + "," + node1.getGlobalIndex() + "," + node2.getGlobalIndex() + ")";
                                hitList.addHitEvent(new HitEventRational(surface, polygon, nx, v2x, centroidX, debugMessage));
                                // System.out.println("hit yz @ y="+samplesY[rayIndexJ]+", z="+samplesZ[rayIndexK]+" ... ray hit at x="+v2x);
                                hitCount++;
                            }
                        }
                    }
                }
            }
        }
    }
    long t2 = System.currentTimeMillis();
    System.out.println("***********************(RationalNumber UV) ray cast XYZ (" + hitCount + " hits), " + (t2 - t1) + "ms *********************");
    return new RayCastResultsRational(hitListsXY, hitListsXZ, hitListsYZ, numX, numY, numZ);
}