# graph3D.frink

```/*    This is a simple but rather interesting program that graphs equations in    3 dimensions.    You enter equations in terms of x, y, and z something like one of the    following:     y = sin[x] + cos[z]     x^2 + y^2 + z^2 = 81      (This is a sphere)     y cos[x] = x sin[z]    This version of the program can also graph INEQUALITIES, which have    less-than or greater-than symbols instead of just equals.    Inequalities are important for graphing infinitely-thin objects and making    them printable and sliceable.  For example, you might need to convert:    x^2 + y^2 - z^2 = 81      into    abs[x^2 + y^2 - z^2 - 81] <= 2    to give the walls some thickness and give your slicer a chance to print it    successfully.  Also increasing the value of "res" below will make the    voxels in the .obj file larger.    Here is an egg.  Modify the 1.5 and 1.1 to change its aspect ratio:    (x^2 + y^2 + z^2)^2 = 6 (1.5 z^3 + (1.5 - 1.1) z (y^2 + x^2))    Here is a Klein bottle:    (x^2 + y^2 + z^2 + 2y - 1)((x^2 + y^2 + z^2 - 2y - 1)^2 - 8 z^2) + 16 x z (x^2 + y^2 + z^2 - 2y - 1) = 0    Here is a heart:    320 (-x^2 z^3 - 9 y^2 z^3 / 80 + (x^2 + 9 y^2 / 4 + z^2 - 1)^3) <= 0        This uses a recursive method to subdivide and test cuboids.    You can also use logical relations like AND and OR to combine multiple    shapes. */ lasteq = "" xmin = -10 xmax =  10 ymin = -10 ymax =  10 zmin = -10 zmax =  10 // Change the doublings to vary the number of voxels.  This is the number // of doublings, so if the number is 10 we have 2^10=1024 doublings for // a resolution of 1024x1024x1024. // Be warned that increasing the doublings by 1 makes 8 times as many voxels! doublings = 9 r = 2^doublings   // Number of voxels on each axis res = 254/in      // Resolution at which to render the .obj file // If there are arguments to the program, graph them, otherwise prompt. while func = (length[ARGS] > 0 ? ARGS@0 : input["Enter equation: ", lasteq]) {    hasInequality = false    certEq = undef    lasteq = certFunc = func    // If there's an inequality, let's make a test equation to see if we can    // fill in an entire cuboid using the "CERTAINLY" comparators.    if func =~ %r/([<>]|!=)/    {       hasInequality = true       certFunc =~ %s/<=/ CLE /g  // Replace <= with certainly less than or equals       certFunc =~ %s/>=/ CGE /g  // Replace >= with certainly greater than or equals       certFunc =~ %s/</ CLT /g   // Replace <  with certainly less than       certFunc =~ %s/>/ CGT /g   // Replace >  with certainly greater than       certFunc =~ %s/!=/ CNE /g  // Replace =  with certainly not equals       certFunc =~ %s/=/ CEQ /g   // Replace =  with certainly equals       certEq = parseToExpression[certFunc]    }    // These replacements turn normal comparator and equality tests into    // "POSSIBLY EQUALS" tests.    func =~ %s/<=/ PLE /g  // Replace <= with possibly less than or equals    func =~ %s/>=/ PGE /g  // Replace >= with possibly greater than or equals    func =~ %s/</ PLT /g   // Replace <  with possibly less than    func =~ %s/>/ PGT /g   // Replace >  with possibly greater than    func =~ %s/!=/ PNE /g  // Replace =  with possibly not equals    func =~ %s/=/ PEQ /g   // Replace =  with possibly equals    eq = parseToExpression[func]    println[inputForm[eq]]    if certEq != undef       println[inputForm[certEq]]    // Scale factors on each axis    sx = r / (xmax-xmin)    sy = r / (ymax-ymin)    sz = r / (zmax-zmin)        v = callJava["frink.graphics.VoxelArray", "construct", [xmin sx, xmax sx, ymin sy, ymax sy, zmin sz, zmax sz, false]]    // Perform the recursive testing of the volume    testCube[xmin, xmax, ymin, ymax, zmin, zmax, v, eq, certEq, doublings, sx, sy, sz]    v.projectX[undef].show["X"]    v.projectY[undef].show["Y"]    v.projectZ[undef].show["Z"]    filename = "graph3D.obj"    print["Writing \$filename..."]    w = new Writer[filename]    w.println[v.toObjFormat["graph3D", 1 / (res mm)]]    w.close[]    println["done."]    if length[ARGS] > 0       exit[] } // Recursive function to test an interval containing the specified bounds. // If no possible solution exists, the recursion halts.  If the entire cube // is filled with a "certainly" equation, it is filled and recursion halts. // If only a possible solution exists, this breaks it down into 8 sub-cuboids // and tests each of them recursively. // level is the maximum number of levels to split, so the total // resolution of the final graph will be 2^level. testCube[x1, x2, y1, y2, z1, z2, v, eq, certEq, level, sx, sy, sz] := {    nextLevel = level - 1    x = new interval[x1, x2]    y = new interval[y1, y2]    z = new interval[z1, z2]        // Test the cuboid.  If it possibly contains solutions, recursively    // subdivide.    res = eval[eq]        if res or res==undef    {       if (nextLevel >= 0)       {          // Do we have inequalities and a CERTAINLY test?          if (certEq != undef) AND (eval[certEq] == true)          {             // If the entire cuboid is a solution, then fill the rectangle             // and stop further recursion on this cuboid.             v.fillCube[x1 sx, x2 sx, y1 sy, y2 sy, z1 sz, z2 sz, true]             return          }          // Further subdivide the cuboid into 8 octants and recursively          // test them all          cx = (x1 + x2)/2          cy = (y1 + y2)/2          cz = (z1 + z2)/2                    testCube[x1, cx, y1, cy, z1, cz, v, eq, certEq, nextLevel, sx, sy, sz]          testCube[cx, x2, y1, cy, z1, cz, v, eq, certEq, nextLevel, sx, sy, sz]          testCube[x1, cx, cy, y2, z1, cz, v, eq, certEq, nextLevel, sx, sy, sz]          testCube[cx, x2, cy, y2, z1, cz, v, eq, certEq, nextLevel, sx, sy, sz]          testCube[x1, cx, y1, cy, cz, z2, v, eq, certEq, nextLevel, sx, sy, sz]          testCube[cx, x2, y1, cy, cz, z2, v, eq, certEq, nextLevel, sx, sy, sz]          testCube[x1, cx, cy, y2, cz, z2, v, eq, certEq, nextLevel, sx, sy, sz]          testCube[cx, x2, cy, y2, cz, z2, v, eq, certEq, nextLevel, sx, sy, sz]                 } else            if (res)             // Valid point at lowest level; fill it               v.fillCube[x1 sx, x2 sx, y1 sy, y2 sy, z1 sz, z2 sz,true]    } } ```