Dymaxion.frink

View or download Dymaxion.frink in plain text format


/* 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[5]] / sqrt[10] ]
   class var gt   = garc / 2.0
   class var gdve = sqrt[3 + sqrt[5]] / sqrt[5 + sqrt[5]]
   class var gel  = sqrt[8] / sqrt[5 + sqrt[5]]
   class var gs1  = sqrt[5 + 2 * sqrt[5] ] / sqrt[15]
   class var sqrt3 = sqrt[3]
   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]
   }
}


View or download Dymaxion.frink in plain text format


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 17683 days, 10 hours, 56 minutes ago.