# 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 18115 days, 17 hours, 36 minutes ago.