/* 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]
   }
}