# Dymaxion.frink

``` /* Conversion of Lat/Long coordinates into x,y coordinates in the Dymaxion    (Buckminster Fuller) icosahedral projection.     See:  http://www.rwgrayprojects.com/rbfnotes/maps/graymap1.html */ class Dymaxion {    // Vertices of icosahedron (array of points)    class var v = undef    // Triangles of icosahedron (array of Triangle)    class var t = undef     class var garc = 2.0 * asin[ sqrt[5 - sqrt] / sqrt ]    class var gt   = garc / 2.0    class var gdve = sqrt[3 + sqrt] / sqrt[5 + sqrt]    class var gel  = sqrt / sqrt[5 + sqrt]    class var gs1  = sqrt[5 + 2 * sqrt ] / sqrt    class var sqrt3 = sqrt    class var trimap = [[0,0,0], [1,3,2], [1,4,3], [1,5,4], [1,6,5], [1,2,6],                        [2,3,8], [3,9,8], [3,4,9], [4,10,9], [4,5,10],                        [5,11,10], [5,6,11], [6,7,11], [2,7,6], [2,8,7],                        [8,9,12], [9,10,12], [10,11,12], [11,7,12], [8,12,7]]    // Rotation and translation for each triangle    class var rAndT = [[0 degrees,   0,  0],                [240 degrees, 2,  7/2/sqrt3],                [300 degrees, 2,  5/2/sqrt3],                [0 degrees,   5/2, 2/sqrt3],                [60 degrees, 3, 5/2/sqrt3],                [180 degrees, 5/2, 4/3 sqrt3],                [300 degrees, 3/2, 4/3 sqrt3],                [300 degrees, 1, 5/2/sqrt3],                [0 degrees, 3/2, 2/sqrt3],                [300 degrees, 3/2, 1/sqrt3],  // 9 can be split                [60 degrees, 5/2, 1/sqrt3],                [60 degrees, 7/2, 1/sqrt3],                [120 degrees, 7/2, 2/sqrt3],                [60 degrees, 4, 5/2/sqrt3],                [0 degrees, 4, 7/2/sqrt3],                [0 degrees, 5, 7/2/sqrt3],                [60 degrees, 1/2, 1/sqrt3],  // 16 can be split                [0 degrees, 1, 1/2/sqrt3],                [120 degrees, 4, 1/2/sqrt3],                [120 degrees, 9/2, 2/sqrt3],                [300 degrees, 5, 5/2/sqrt3]]           class var initialized = false    /** Convert [lat, long] to x,y coordinates in the Dymaxion projection.        This is the main function that the user will use.  If splitTriangles is        true (the default,) this will make a canonical Fuller projection that        splits region 16 (containing Antarctica) and region 9 (containing        Indonesia) into sub-triangles so that they're contiguous.  However, this        is not usually necessary when mapping an image onto the icosahedron,        and is usually undesirable when building a 3-dimensional model, for        which cases splitTriangles should be set to false.            returns [x, y, triangle, lcd]    */    class latLongToXY[lat, long, splitTriangles=true] :=    {       if ! initialized          init[]              [theta, phi] = latLongToSpherical[lat, long]       point = sphericalToCartesian[theta, phi]              /* determine which of the 20 spherical icosahedron triangles */       /* the given point is in and the LCD triangle.               */       [triangle, lcd] = getTriangleInfo[point]       /* Determine the corresponding Fuller map plane (x, y) point */       [x,y] = dymaxionPoint[point, triangle, lcd, splitTriangles]              return [x,y,triangle,lcd]    }        /** Initialize class-level variables */    class init[] :=    {       if initialized          return              v = new array       /* Cartesian coordinates for the 12 vertices of the icosahedron */       v@1  = new Point3D[ 0.420152426708710003,  0.078145249402782959,  0.904082550615019298]       v@2  = new Point3D[ 0.995009439436241649, -0.091347795276427931,  0.040147175877166645]       v@3  = new Point3D[ 0.518836730327364437,  0.835420380378235850,  0.181331837557262454]       v@4  = new Point3D[-0.414682225320335218,  0.655962405434800777,  0.630675807891475371]       v@5  = new Point3D[-0.515455959944041808, -0.381716898287133011,  0.767200992517747538]       v@6  = new Point3D[ 0.355781402532944713, -0.843580002466178147,  0.402234226602925571]       v@7  = new Point3D[ 0.414682225320335218, -0.655962405434800777, -0.630675807891475371]       v@8  = new Point3D[ 0.515455959944041808,  0.381716898287133011, -0.767200992517747538]       v@9  = new Point3D[-0.355781402532944713,  0.843580002466178147, -0.402234226602925571]       v@10 = new Point3D[-0.995009439436241649,  0.091347795276427931, -0.040147175877166645]       v@11 = new Point3D[-0.518836730327364437, -0.835420380378235850, -0.181331837557262454]       v@12 = new Point3D[-0.420152426708710003, -0.078145249402782959, -0.904082550615019298]       t = new array       // Create triangles for each side       for i = 1 to 20       {          [a,b,c] = trimap@i          t@i = new Triangle[v@a,v@b,v@c]       }       initialized = true    }    /** Convert lat/long to spherical coordinates [theta, phi] with radius 1 */    class latLongToSpherical[lat, long] :=    {       theta = 90 degrees - lat ;       phi = long mod circle       return[theta, phi]    }    /** Converts spherical coordinates [theta, phi] to 3-D cartesian        coordinates.        returns:  Point3D */    class sphericalToCartesian[theta, phi] :=    {       st = sin[theta]       return new Point3D[st * cos[phi],                          st * sin[phi],                          cos[theta]]    }    /** Convert a point in 3-D cartesian coordinates to lat,long coordinates.        returns[lat, long]    */    class cartesianToLatLong[p is Point3D] :=    {       return [acos[p.z], atan[p.y, p.x]]    }        /** Determine which triangle and LCD triangle a point is in.        Returns the index of the triangle        and its LCD triangle:        [triangle, lcd]    */    class getTriangleInfo[p is Point3D] :=    {       if ! initialized          init[]              tri = 0       h_dist1 = 9999.0       var v1       var v2       var v3       /* Which triangle face center is the closest to the given point */       /* is the triangle in which the given point is in.              */       for i = 1 to 20       {          h_dist2 = p.getDistanceSquared[t@i.center]          if (h_dist2 < h_dist1)          {             tri = i;             h_dist1 = h_dist2;          }       }       [v1,v2,v3] = trimap@tri       h_dist1 = p.getDistanceSquared[v@v1]       h_dist2 = p.getDistanceSquared[v@v2]       h_dist3 = p.getDistanceSquared[v@v3]       if ( (h_dist1 <= h_dist2) && (h_dist2 <= h_dist3) )          lcd = 1       if ( (h_dist1 <= h_dist3) && (h_dist3 <= h_dist2) )          lcd = 6       if ( (h_dist2 <= h_dist1) && (h_dist1 <= h_dist3) )          lcd = 2       if ( (h_dist2 <= h_dist3) && (h_dist3 <= h_dist1) )          lcd = 3       if ( (h_dist3 <= h_dist1) && (h_dist1 <= h_dist2) )          lcd = 5;       if ( (h_dist3 <= h_dist2) && (h_dist2 <= h_dist1) )          lcd = 4;       return [tri, lcd]    }    /** Gets the [x,y] coordinates of the point in the Dymaxion projection.        returns [x,y]    */    class dymaxionPoint[p is Point3D, tri, lcd, splitTriangles] :=    {       v1 = trimap@tri@0       p1 = v@v1       [lat, long] = cartesianToLatLong[t@tri.center]       axis = 3       p = p.rotate[axis,long]       p1 = p1.rotate[axis,long]       axis = 2       p  =  p.rotate[axis,lat]       p1 = p1.rotate[axis,lat]       [lat, long] = cartesianToLatLong[p1]       long = long - 90 degrees       axis = 3       p = p.rotate[axis,long]       /* exact transformation equations */       gz = sqrt[1 - (p.x)^2 - (p.y)^2]       gs = gs1 / gz       gxp = p.x * gs       gyp = p.y * gs       ga1p = 2.0 * gyp / sqrt3 + (gel / 3.0)        ga2p = gxp - (gyp / sqrt3) +  (gel / 3.0)        ga3p = (gel / 3.0) - gxp - (gyp / sqrt3)       ga1 = gt + atan[(ga1p - 0.5 * gel) / gdve]       ga2 = gt + atan[(ga2p - 0.5 * gel) / gdve]       ga3 = gt + atan[(ga3p - 0.5 * gel) / gdve]       gx = 0.5 * (ga2 - ga3)        gy = (1.0 / (2.0 * sqrt3) ) * (2 * ga1 - ga2 - ga3)       /* Re-scale so plane triangle edge length is 1. */       x = gx / garc       y = gy / garc       if splitTriangles       {          if tri== 9             if (lcd > 2)                return rotateAndTranslate[300 degrees, 3/2, 1/sqrt3, x, y]             else                return rotateAndTranslate[0 degrees, 2, 1/2/sqrt3, x, y]                       if tri==16             if (lcd < 4)                return rotateAndTranslate[60 degrees, 1/2, 1/sqrt3, x, y]             else                return rotateAndTranslate[0 degrees, 11/2, 2/sqrt3, x, y]       }       [angle, xoff, yoff] = rAndT@tri       return rotateAndTranslate[angle, xoff, yoff, x, y]    }    /** Rotates around the specified angle in the plane.        returns [rx, ry] */    class rotateAndTranslate[angle, xoff, yoff, x, y] :=    {       ca = cos[angle]       sa = sin[angle]       return [x * ca - y * sa + xoff,               x * sa + y * ca + yoff]    } } class Point3D {    var x    var y    var z    new[xx,yy,zz] :=    {       x = xx       y = yy       z = zz    }    /* Gets the distance to another point. */    getDistance[other is Point3D] := sqrt[(x - other.x)^2 + (y - other.y)^2 + (z - other.z)^2]    /** Gets the square of the distance.  This is sufficient to compare whether        one point is closer than another and eliminates the need to perform        a square root. */    getDistanceSquared[other is Point3D] := (x - other.x)^2 + (y - other.y)^2 + (z - other.z)^2        /* Rotate a 3-D point around the specified axis.       axis:  1=x-axis, 2=y-axis, 3=z=axis  */    rotate[axis, angle] :=    {       ca = cos[angle]       sa = sin[angle]       if (axis == 1)   // x-axis          return new Point3D[x,                             y * ca + z * sa,                             z * ca - y * sa]       if (axis == 2)  // y-axis          return new Point3D[x * ca - z * sa,                             y,                             x * sa + z * ca]       if (axis == 3)  // z-axis          return new Point3D[x * ca + y * sa,                             y * ca - x * sa,                             z]    } } /** A triangle defined by three Point3D coordinates with a center point. */ class Triangle {    var p1    var p2    var p3    var center    // Point3D    /** Construct a Triangle and calculate its center point. */    new[p1a, p2a, p3a] :=    {       p1 = p1a       p2 = p2a       p3 = p3a       hx = (p1.x + p2.x + p3.x) / 3       hy = (p1.y + p2.y + p3.y) / 3       hz = (p1.z + p2.z + p3.z) / 3       mag = sqrt[hx^2 + hy^2 + hz^2]       center = new Point3D[hx/mag, hy/mag, hz/mag]    } } ```

This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 18491 days, 5 hours, 1 minutes ago.