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 20145 days, 21 hours, 26 minutes ago.