Download or view 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]
}
}
Download or view 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 19966 days, 23 hours, 4 minutes ago.