git.hungrycats.org Git - xscreensaver/blob - hacks/glx/dymaxionmap-coords.c

   1 /* http://www.rwgrayprojects.com/rbfnotes/maps/graymap6.html
   2    Slightly modified by jwz for xscreensaver
   3  */
   4
   5
   6 /**************************************************************/
   7 /*                                                            */
   8 /* This C program is copyrighted by  Robert W. Gray and may   */
   9 /* not be used in ANY for-profit project without written      */
  10 /* permission.                                                */
  11 /*                                                            */
  12 /**************************************************************/
  13
  14 /* (Note: Robert Gray has kindly given me his permission to include
  15    this code in xscreensaver. -- Jamie Zawinski, Apr 2018.)
  16  */
  17
  18
  19 /**************************************************************/
  20 /*                                                            */
  21 /* This C program contains the Dymaxion map coordinate        */
  22 /* transformation routines for converting longitude/latitude  */
  23 /* points to (X, Y) points on the Dymaxion map.               */
  24 /*                                                            */
  25 /* This version uses the exact transformation equations.      */
  26 /**************************************************************/
  27
  28 #include <stdio.h>
  29 #include <stdlib.h>
  30 #include <math.h>
  31
  32 #ifdef HAVE_CONFIG_H
  33 # include "config.h"
  34 #endif
  35
  36 #include "dymaxionmap-coords.h"
  37
  38 /************************************************************************/
  39 /* NOTE: in C, array indexing starts with element zero (0).  I choose   */
  40 /*       to start my array indexing with elemennt one (1) so all arrays */
  41 /*       are defined one element longer than they need to be.           */
  42 /************************************************************************/
  43
  44 /************************************************************************/
  45 /* global variables accessable to all procedures                        */
  46 /************************************************************************/
  47
  48 static double v_x[13], v_y[13], v_z[13];
  49 static double center_x[21], center_y[21], center_z[21];
  50 static double garc, gt, gdve, gel;
  51
  52 /********************************************/
  53 /*      function pre-definitions            */
  54 /********************************************/
  55
  56 static double radians(double degrees);
  57 static void rotate(double angle, double *x, double *y);
  58 static void r2(int axis, double alpha, double *x, double *y, double *z);
  59 static void init_stuff(void);
  60 /*static void convert_s_t_p(double lng, double lat, double *x, double *y);*/
  61 static void s_to_c(double theta, double phi, double *x, double *y, double *z);
  62 static void c_to_s(double *theta, double *phi, double x, double y, double z);
  63 static void s_tri_info(double x, double y, double z,
  64                  int *tri, int *lcd);
  65 static void dymax_point(int tri, int lcd,
  66                  double x, double y, double z,
  67                  double *dx, double *dy);
  68 static void conv_ll_t_sc(double lng, double lat, double *theta, double *phi);
  69
  70
  71 /****************************************/
  72 /*      function definitions            */
  73 /****************************************/
  74
  75
  76 void
  77 /* convert_s_t_p */
  78 dymaxion_convert
  79 (double lng, double lat, double *x, double *y)
  80 {
  81   /***********************************************************/
  82   /* This is the main control procedure.                     */
  83   /***********************************************************/
  84
  85   double theta, phi;
  86   double hx, hy, hz;
  87   double px = 0, py = 0;
  88   int tri, hlcd;
  89
  90   static int initted = 0;
  91   if (! initted) {
  92     init_stuff();
  93     initted = 1;
  94   }
  95
  96   /* Convert the given (long.,lat.) coordinate into spherical */
  97   /* polar coordinates (r, theta, phi) with radius=1.         */
  98   /* Angles are given in radians, NOT degrees.                */
  99
 100   conv_ll_t_sc(lng, lat, &theta, &phi);
 101
 102   /* convert the spherical polar coordinates into cartesian   */
 103   /* (x, y, z) coordinates.                                   */
 104
 105   s_to_c(theta, phi, &hx, &hy, &hz);
 106
 107   /* determine which of the 20 spherical icosahedron triangles */
 108   /* the given point is in and the LCD triangle.               */
 109
 110   s_tri_info(hx, hy, hz, &tri, &hlcd);
 111
 112   /* Determine the corresponding Fuller map plane (x, y) point */
 113
 114   dymax_point(tri, hlcd, hx, hy, hz, &px, &py);
 115   *x = px;
 116   *y = py;
 117
 118 } /* end convert_s_t_p */
 119
 120
 121 static void conv_ll_t_sc(double lng, double lat, double *theta, double *phi)
 122 {
 123   /* convert (long., lat.) point into spherical polar coordinates */
 124   /* with r=radius=1.  Angles are given in radians.               */
 125
 126   double h_theta, h_phi;
 127
 128   h_theta = 90.0 - lat ;
 129   h_phi = lng;
 130   if (lng < 0.0) {h_phi = lng + 360.0;}
 131   *theta = radians(h_theta);
 132   *phi = radians(h_phi);
 133
 134 } /* end conv_ll_t_sc */
 135
 136
 137 static double radians(double degrees)
 138 {
 139     /* convert angles in degrees into angles in radians */
 140
 141     double pi2, c1;
 142
 143     pi2 = 2 * 3.14159265358979323846;
 144     c1 = pi2 / 360;
 145     return(c1 * degrees);
 146
 147 } /* end of radians function */
 148
 149
 150 static void init_stuff()
 151 {
 152    /* initializes the global variables which includes the */
 153    /* vertix coordinates and mid-face coordinates.        */
 154
 155    double /* i, */ hold_x, hold_y, hold_z, magn;
 156    /* double theta, phi; */
 157
 158    /* Cartesian coordinates for the 12 vertices of icosahedron */
 159
 160    v_x[1] =    0.420152426708710003;
 161    v_y[1] =    0.078145249402782959;
 162    v_z[1] =    0.904082550615019298;
 163    v_x[2] =    0.995009439436241649 ;
 164    v_y[2] =   -0.091347795276427931 ;
 165    v_z[2] =    0.040147175877166645 ;
 166    v_x[3] =    0.518836730327364437 ;
 167    v_y[3] =    0.835420380378235850 ;
 168    v_z[3] =    0.181331837557262454 ;
 169    v_x[4] =   -0.414682225320335218 ;
 170    v_y[4] =    0.655962405434800777 ;
 171    v_z[4] =    0.630675807891475371 ;
 172    v_x[5] =   -0.515455959944041808 ;
 173    v_y[5] =   -0.381716898287133011 ;
 174    v_z[5] =    0.767200992517747538 ;
 175    v_x[6] =    0.355781402532944713 ;
 176    v_y[6] =   -0.843580002466178147 ;
 177    v_z[6] =    0.402234226602925571 ;
 178    v_x[7] =    0.414682225320335218 ;
 179    v_y[7] =   -0.655962405434800777 ;
 180    v_z[7] =   -0.630675807891475371 ;
 181    v_x[8] =    0.515455959944041808 ;
 182    v_y[8] =    0.381716898287133011 ;
 183    v_z[8] =   -0.767200992517747538 ;
 184    v_x[9] =   -0.355781402532944713 ;
 185    v_y[9] =    0.843580002466178147 ;
 186    v_z[9] =   -0.402234226602925571 ;
 187    v_x[10] =   -0.995009439436241649 ;
 188    v_y[10] =    0.091347795276427931 ;
 189    v_z[10] =   -0.040147175877166645 ;
 190    v_x[11] =   -0.518836730327364437  ;
 191    v_y[11] =   -0.835420380378235850  ;
 192    v_z[11] =   -0.181331837557262454  ;
 193    v_x[12] =   -0.420152426708710003 ;
 194    v_y[12] =   -0.078145249402782959 ;
 195    v_z[12] =   -0.904082550615019298 ;
 196
 197    /* now calculate mid face coordinates             */
 198
 199    hold_x = (v_x[1] + v_x[2] + v_x[3]) / 3.0 ;
 200    hold_y = (v_y[1] + v_y[2] + v_y[3]) / 3.0 ;
 201    hold_z = (v_z[1] + v_z[2] + v_z[3]) / 3.0 ;
 202    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 203    center_x[1] = hold_x / magn;
 204    center_y[1] = hold_y / magn;
 205    center_z[1] = hold_z / magn;
 206
 207    hold_x = (v_x[1] + v_x[3] + v_x[4]) / 3.0 ;
 208    hold_y = (v_y[1] + v_y[3] + v_y[4]) / 3.0 ;
 209    hold_z = (v_z[1] + v_z[3] + v_z[4]) / 3.0 ;
 210    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 211    center_x[2] = hold_x / magn;
 212    center_y[2] = hold_y / magn;
 213    center_z[2] = hold_z / magn;
 214
 215    hold_x = (v_x[1] + v_x[4] + v_x[5]) / 3.0 ;
 216    hold_y = (v_y[1] + v_y[4] + v_y[5]) / 3.0 ;
 217    hold_z = (v_z[1] + v_z[4] + v_z[5]) / 3.0 ;
 218    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 219    center_x[3] = hold_x / magn;
 220    center_y[3] = hold_y / magn;
 221    center_z[3] = hold_z / magn;
 222
 223    hold_x = (v_x[1] + v_x[5] + v_x[6]) / 3.0 ;
 224    hold_y = (v_y[1] + v_y[5] + v_y[6]) / 3.0 ;
 225    hold_z = (v_z[1] + v_z[5] + v_z[6]) / 3.0 ;
 226    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 227    center_x[4] = hold_x / magn;
 228    center_y[4] = hold_y / magn;
 229    center_z[4] = hold_z / magn;
 230
 231    hold_x = (v_x[1] + v_x[2] + v_x[6]) / 3.0 ;
 232    hold_y = (v_y[1] + v_y[2] + v_y[6]) / 3.0 ;
 233    hold_z = (v_z[1] + v_z[2] + v_z[6]) / 3.0 ;
 234    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 235    center_x[5] = hold_x / magn;
 236    center_y[5] = hold_y / magn;
 237    center_z[5] = hold_z / magn;
 238
 239    hold_x = (v_x[2] + v_x[3] + v_x[8]) / 3.0 ;
 240    hold_y = (v_y[2] + v_y[3] + v_y[8]) / 3.0 ;
 241    hold_z = (v_z[2] + v_z[3] + v_z[8]) / 3.0 ;
 242    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 243    center_x[6] = hold_x / magn;
 244    center_y[6] = hold_y / magn;
 245    center_z[6] = hold_z / magn;
 246
 247    hold_x = (v_x[8] + v_x[3] + v_x[9]) / 3.0 ;
 248    hold_y = (v_y[8] + v_y[3] + v_y[9]) / 3.0 ;
 249    hold_z = (v_z[8] + v_z[3] + v_z[9]) / 3.0 ;
 250    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 251    center_x[7] = hold_x / magn;
 252    center_y[7] = hold_y / magn;
 253    center_z[7] = hold_z / magn;
 254
 255    hold_x = (v_x[9] + v_x[3] + v_x[4]) / 3.0 ;
 256    hold_y = (v_y[9] + v_y[3] + v_y[4]) / 3.0 ;
 257    hold_z = (v_z[9] + v_z[3] + v_z[4]) / 3.0 ;
 258    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 259    center_x[8] = hold_x / magn;
 260    center_y[8] = hold_y / magn;
 261    center_z[8] = hold_z / magn;
 262
 263    hold_x = (v_x[10] + v_x[9] + v_x[4]) / 3.0 ;
 264    hold_y = (v_y[10] + v_y[9] + v_y[4]) / 3.0 ;
 265    hold_z = (v_z[10] + v_z[9] + v_z[4]) / 3.0 ;
 266    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 267    center_x[9] = hold_x / magn;
 268    center_y[9] = hold_y / magn;
 269    center_z[9] = hold_z / magn;
 270
 271    hold_x = (v_x[5] + v_x[10] + v_x[4]) / 3.0 ;
 272    hold_y = (v_y[5] + v_y[10] + v_y[4]) / 3.0 ;
 273    hold_z = (v_z[5] + v_z[10] + v_z[4]) / 3.0 ;
 274    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 275    center_x[10] = hold_x / magn;
 276    center_y[10] = hold_y / magn;
 277    center_z[10] = hold_z / magn;
 278
 279    hold_x = (v_x[5] + v_x[11] + v_x[10]) / 3.0 ;
 280    hold_y = (v_y[5] + v_y[11] + v_y[10]) / 3.0 ;
 281    hold_z = (v_z[5] + v_z[11] + v_z[10]) / 3.0 ;
 282    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 283    center_x[11] = hold_x / magn;
 284    center_y[11] = hold_y / magn;
 285    center_z[11] = hold_z / magn;
 286
 287    hold_x = (v_x[5] + v_x[6] + v_x[11]) / 3.0 ;
 288    hold_y = (v_y[5] + v_y[6] + v_y[11]) / 3.0 ;
 289    hold_z = (v_z[5] + v_z[6] + v_z[11]) / 3.0 ;
 290    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 291    center_x[12] = hold_x / magn;
 292    center_y[12] = hold_y / magn;
 293    center_z[12] = hold_z / magn;
 294
 295    hold_x = (v_x[11] + v_x[6] + v_x[7]) / 3.0 ;
 296    hold_y = (v_y[11] + v_y[6] + v_y[7]) / 3.0 ;
 297    hold_z = (v_z[11] + v_z[6] + v_z[7]) / 3.0 ;
 298    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 299    center_x[13] = hold_x / magn;
 300    center_y[13] = hold_y / magn;
 301    center_z[13] = hold_z / magn;
 302
 303    hold_x = (v_x[7] + v_x[6] + v_x[2]) / 3.0 ;
 304    hold_y = (v_y[7] + v_y[6] + v_y[2]) / 3.0 ;
 305    hold_z = (v_z[7] + v_z[6] + v_z[2]) / 3.0 ;
 306    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 307    center_x[14] = hold_x / magn;
 308    center_y[14] = hold_y / magn;
 309    center_z[14] = hold_z / magn;
 310
 311    hold_x = (v_x[8] + v_x[7] + v_x[2]) / 3.0 ;
 312    hold_y = (v_y[8] + v_y[7] + v_y[2]) / 3.0 ;
 313    hold_z = (v_z[8] + v_z[7] + v_z[2]) / 3.0 ;
 314    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 315    center_x[15] = hold_x / magn;
 316    center_y[15] = hold_y / magn;
 317    center_z[15] = hold_z / magn;
 318
 319    hold_x = (v_x[12] + v_x[9] + v_x[8]) / 3.0 ;
 320    hold_y = (v_y[12] + v_y[9] + v_y[8]) / 3.0 ;
 321    hold_z = (v_z[12] + v_z[9] + v_z[8]) / 3.0 ;
 322    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 323    center_x[16] = hold_x / magn;
 324    center_y[16] = hold_y / magn;
 325    center_z[16] = hold_z / magn;
 326
 327    hold_x = (v_x[12] + v_x[9] + v_x[10]) / 3.0 ;
 328    hold_y = (v_y[12] + v_y[9] + v_y[10]) / 3.0 ;
 329    hold_z = (v_z[12] + v_z[9] + v_z[10]) / 3.0 ;
 330    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 331    center_x[17] = hold_x / magn;
 332    center_y[17] = hold_y / magn;
 333    center_z[17] = hold_z / magn;
 334
 335    hold_x = (v_x[12] + v_x[11] + v_x[10]) / 3.0 ;
 336    hold_y = (v_y[12] + v_y[11] + v_y[10]) / 3.0 ;
 337    hold_z = (v_z[12] + v_z[11] + v_z[10]) / 3.0 ;
 338    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 339    center_x[18] = hold_x / magn;
 340    center_y[18] = hold_y / magn;
 341    center_z[18] = hold_z / magn;
 342
 343    hold_x = (v_x[12] + v_x[11] + v_x[7]) / 3.0 ;
 344    hold_y = (v_y[12] + v_y[11] + v_y[7]) / 3.0 ;
 345    hold_z = (v_z[12] + v_z[11] + v_z[7]) / 3.0 ;
 346    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 347    center_x[19] = hold_x / magn;
 348    center_y[19] = hold_y / magn;
 349    center_z[19] = hold_z / magn;
 350
 351    hold_x = (v_x[12] + v_x[8] + v_x[7]) / 3.0 ;
 352    hold_y = (v_y[12] + v_y[8] + v_y[7]) / 3.0 ;
 353    hold_z = (v_z[12] + v_z[8] + v_z[7]) / 3.0 ;
 354    magn = sqrt(hold_x * hold_x + hold_y * hold_y + hold_z * hold_z);
 355    center_x[20] = hold_x / magn;
 356    center_y[20] = hold_y / magn;
 357    center_z[20] = hold_z / magn;
 358
 359    garc = 2.0 * asin( sqrt( 5 - sqrt(5)) / sqrt(10) );
 360    gt = garc / 2.0;
 361
 362    gdve = sqrt( 3 + sqrt(5) ) / sqrt( 5 + sqrt(5) );
 363    gel = sqrt(8) / sqrt(5 + sqrt(5));
 364
 365 } /* end of int_stuff procedure */
 366
 367
 368 static void s_to_c(double theta, double phi, double *x, double *y, double *z)
 369 {
 370     /* Covert spherical polar coordinates to cartesian coordinates. */
 371     /* The angles are given in radians.                             */
 372
 373     *x = sin(theta) * cos(phi);
 374     *y = sin(theta) * sin(phi);
 375     *z = cos(theta);
 376
 377  } /* end s_to_c */
 378
 379
 380 static void c_to_s(double *lng, double *lat, double x, double y, double z)
 381    {
 382     /* convert cartesian coordinates into spherical polar coordinates. */
 383     /* The angles are given in radians.                                */
 384
 385     double a;
 386
 387     if (x>0.0 && y>0.0){a = radians(0.0);}
 388     if (x<0.0 && y>0.0){a = radians(180.0);}
 389     if (x<0.0 && y<0.0){a = radians(180.0);}
 390     if (x>0.0 && y<0.0){a = radians(360.0);}
 391     *lat = acos(z);
 392     if (x==0.0 && y>0.0){*lng = radians(90.0);}
 393     if (x==0.0 && y<0.0){*lng = radians(270.0);}
 394     if (x>0.0 && y==0.0){*lng = radians(0.0);}
 395     if (x<0.0 && y==0.0){*lng = radians(180.0);}
 396     if (x!=0.0 && y!=0.0){*lng = atan(y/x) + a;}
 397
 398 } /* end c_to_s */
 399
 400
 401 void s_tri_info(double x, double y, double z,
 402                  int *tri, int *lcd)
 403 {
 404   /* Determine which triangle and LCD triangle the point is in. */
 405
 406   double  h_dist1, h_dist2, h_dist3, h1, h2, h3;
 407   int i, h_tri, h_lcd ;
 408   int v1 = 0, v2 = 0, v3 = 0;
 409
 410   h_tri = 0;
 411   h_dist1 = 9999.0;
 412
 413   /* Which triangle face center is the closest to the given point */
 414   /* is the triangle in which the given point is in.              */
 415
 416   for (i = 1; i <=20; i = i + 1)
 417    {
 418      h1 = center_x[i] - x;
 419      h2 = center_y[i] - y;
 420      h3 = center_z[i] - z;
 421      h_dist2 = sqrt(h1 * h1 + h2 * h2 + h3 * h3);
 422      if (h_dist2 < h_dist1)
 423       {
 424         h_tri = i;
 425         h_dist1 = h_dist2;
 426       } /* end the if statement */
 427    }  /* end the for statement */
 428
 429    *tri = h_tri;
 430
 431    /* Now the LCD triangle is determined. */
 432
 433    switch (h_tri)
 434    {
 435     case 1:  v1 =  1; v2 =  3; v3 =  2; break;
 436     case 2:  v1 =  1; v2 =  4; v3 =  3; break;
 437     case 3:  v1 =  1; v2 =  5; v3 =  4; break;
 438     case 4:  v1 =  1; v2 =  6; v3 =  5; break;
 439     case 5:  v1 =  1; v2 =  2; v3 =  6; break;
 440     case 6:  v1 =  2; v2 =  3; v3 =  8; break;
 441     case 7:  v1 =  3; v2 =  9; v3 =  8; break;
 442     case 8:  v1 =  3; v2 =  4; v3 =  9; break;
 443     case 9:  v1 =  4; v2 = 10; v3 =  9; break;
 444     case 10: v1 =  4; v2 =  5; v3 = 10; break;
 445     case 11: v1 =  5; v2 = 11; v3 = 10; break;
 446     case 12: v1 =  5; v2 =  6; v3 = 11; break;
 447     case 13: v1 =  6; v2 =  7; v3 = 11; break;
 448     case 14: v1 =  2; v2 =  7; v3 =  6; break;
 449     case 15: v1 =  2; v2 =  8; v3 =  7; break;
 450     case 16: v1 =  8; v2 =  9; v3 = 12; break;
 451     case 17: v1 =  9; v2 = 10; v3 = 12; break;
 452     case 18: v1 = 10; v2 = 11; v3 = 12; break;
 453     case 19: v1 = 11; v2 =  7; v3 = 12; break;
 454     case 20: v1 =  8; v2 = 12; v3 =  7; break;
 455    } /* end of switch statement */
 456
 457    h1 = x - v_x[v1];
 458    h2 = y - v_y[v1];
 459    h3 = z - v_z[v1];
 460    h_dist1 = sqrt(h1 * h1 + h2 * h2 + h3 * h3);
 461
 462    h1 = x - v_x[v2];
 463    h2 = y - v_y[v2];
 464    h3 = z - v_z[v2];
 465    h_dist2 = sqrt(h1 * h1 + h2 * h2 + h3 * h3);
 466
 467    h1 = x - v_x[v3];
 468    h2 = y - v_y[v3];
 469    h3 = z - v_z[v3];
 470    h_dist3 = sqrt(h1 * h1 + h2 * h2 + h3 * h3);
 471
 472    if ( (h_dist1 <= h_dist2) && (h_dist2 <= h_dist3) ) {h_lcd = 1; }
 473    if ( (h_dist1 <= h_dist3) && (h_dist3 <= h_dist2) ) {h_lcd = 6; }
 474    if ( (h_dist2 <= h_dist1) && (h_dist1 <= h_dist3) ) {h_lcd = 2; }
 475    if ( (h_dist2 <= h_dist3) && (h_dist3 <= h_dist1) ) {h_lcd = 3; }
 476    if ( (h_dist3 <= h_dist1) && (h_dist1 <= h_dist2) ) {h_lcd = 5; }
 477    if ( (h_dist3 <= h_dist2) && (h_dist2 <= h_dist1) ) {h_lcd = 4; }
 478
 479    *lcd = h_lcd;
 480
 481 } /* end s_tri_info */
 482
 483
 484 static void dymax_point(int tri, int lcd,
 485                  double x, double y, double z,
 486                  double *px, double *py)
 487 {
 488   int axis, v1 = 0;
 489   double hlng, hlat, h0x, h0y, h0z, h1x, h1y, h1z;
 490
 491   double gs;
 492   double gx, gy, gz, ga1,ga2,ga3,ga1p,ga2p,ga3p,gxp,gyp/*,gzp*/;
 493
 494
 495   /* In order to rotate the given point into the template spherical */
 496   /* triangle, we need the spherical polar coordinates of the center */
 497   /* of the face and one of the face vertices. So set up which vertex */
 498   /* to use.                                                          */
 499
 500    switch (tri)
 501    {
 502     case 1:  v1 =  1;  break;
 503     case 2:  v1 =  1;  break;
 504     case 3:  v1 =  1;  break;
 505     case 4:  v1 =  1;  break;
 506     case 5:  v1 =  1;  break;
 507     case 6:  v1 =  2;  break;
 508     case 7:  v1 =  3;  break;
 509     case 8:  v1 =  3;  break;
 510     case 9:  v1 =  4;  break;
 511     case 10: v1 =  4;  break;
 512     case 11: v1 =  5;  break;
 513     case 12: v1 =  5;  break;
 514     case 13: v1 =  6;  break;
 515     case 14: v1 =  2;  break;
 516     case 15: v1 =  2;  break;
 517     case 16: v1 =  8;  break;
 518     case 17: v1 =  9;  break;
 519     case 18: v1 = 10;  break;
 520     case 19: v1 = 11;  break;
 521     case 20: v1 =  8;  break;
 522    } /* end of switch statement */
 523
 524    h0x = x;
 525    h0y = y;
 526    h0z = z;
 527
 528    h1x = v_x[v1];
 529    h1y = v_y[v1];
 530    h1z = v_z[v1];
 531
 532    c_to_s(&hlng, &hlat, center_x[tri], center_y[tri], center_z[tri]);
 533
 534    axis = 3;
 535    r2(axis,hlng,&h0x,&h0y,&h0z);
 536    r2(axis,hlng,&h1x,&h1y,&h1z);
 537
 538    axis = 2;
 539    r2(axis,hlat,&h0x,&h0y,&h0z);
 540    r2(axis,hlat,&h1x,&h1y,&h1z);
 541
 542    c_to_s(&hlng,&hlat,h1x,h1y,h1z);
 543    hlng = hlng - radians(90.0);
 544
 545    axis = 3;
 546    r2(axis,hlng,&h0x,&h0y,&h0z);
 547
 548    /* exact transformation equations */
 549
 550    gz = sqrt(1 - h0x * h0x - h0y * h0y);
 551    gs = sqrt( 5 + 2 * sqrt(5) ) / ( gz * sqrt(15) );
 552
 553    gxp = h0x * gs ;
 554    gyp = h0y * gs ;
 555
 556    ga1p = 2.0 * gyp / sqrt(3.0) + (gel / 3.0) ;
 557    ga2p = gxp - (gyp / sqrt(3)) +  (gel / 3.0) ;
 558    ga3p = (gel / 3.0) - gxp - (gyp / sqrt(3));
 559
 560    ga1 = gt + atan( (ga1p - 0.5 * gel) / gdve);
 561    ga2 = gt + atan( (ga2p - 0.5 * gel) / gdve);
 562    ga3 = gt + atan( (ga3p - 0.5 * gel) / gdve);
 563
 564    gx = 0.5 * (ga2 - ga3) ;
 565
 566    gy = (1.0 / (2.0 * sqrt(3)) ) * (2 * ga1 - ga2 - ga3);
 567
 568    /* Re-scale so plane triangle edge length is 1. */
 569
 570    x = gx / garc;
 571    y = gy / garc;
 572
 573   /* rotate and translate to correct position          */
 574
 575   switch (tri)
 576    {
 577      case  1: rotate(240.0,&x, &y);
 578               *px = x + 2.0; *py = y + 7.0 / (2.0 * sqrt(3.0)) ; break;
 579      case  2: rotate(300.0, &x, &y); *px = x + 2.0;
 580               *py = y + 5.0 / (2.0 * sqrt(3.0)) ; break;
 581      case  3: rotate(0.0, &x, &y);
 582              *px = x + 2.5; *py = y + 2.0 / sqrt(3.0); break;
 583      case  4: rotate(60.0, &x, &y);
 584               *px = x + 3.0; *py = y + 5.0 / (2.0 * sqrt(3.0)) ; break;
 585      case  5: rotate(180.0, &x, &y);
 586               *px = x + 2.5; *py = y + 4.0 * sqrt(3.0) / 3.0; break;
 587      case  6: rotate(300.0, &x, &y);
 588               *px = x + 1.5; *py = y + 4.0 * sqrt(3.0) / 3.0; break;
 589      case  7: rotate(300.0, &x, &y);
 590               *px = x + 1.0; *py = y + 5.0 / (2.0 * sqrt(3.0)) ; break;
 591      case  8: rotate(0.0, &x, &y);
 592               *px = x + 1.5; *py = y + 2.0 / sqrt(3.0); break;
 593      case  9: if (lcd > 2)
 594               {
 595               rotate(300.0, &x, &y);
 596               *px = x + 1.5; *py = y + 1.0 / sqrt(3.0);
 597               }
 598               else
 599               {
 600               rotate(0.0, &x, &y);
 601               *px = x + 2.0; *py = y + 1.0 / (2.0 * sqrt(3.0));
 602               }
 603               break;
 604
 605      case 10: rotate(60.0, &x, &y);
 606               *px = x + 2.5; *py = y + 1.0 / sqrt(3.0); break;
 607      case 11: rotate(60.0, &x, &y);
 608               *px = x + 3.5; *py = y + 1.0 / sqrt(3.0); break;
 609      case 12: rotate(120.0, &x, &y);
 610               *px = x + 3.5; *py = y + 2.0 / sqrt(3.0); break;
 611      case 13: rotate(60.0, &x, &y);
 612               *px = x + 4.0; *py = y + 5.0 / (2.0 * sqrt(3.0)); break;
 613      case 14: rotate(0.0, &x, &y);
 614               *px = x + 4.0; *py = y + 7.0 / (2.0 * sqrt(3.0)) ; break;
 615      case 15: rotate(0.0, &x, &y);
 616               *px = x + 5.0; *py = y + 7.0 / (2.0 * sqrt(3.0)) ; break;
 617      case 16: if (lcd < 4)
 618               {
 619                 rotate(60.0, &x, &y);
 620                 *px = x + 0.5; *py = y + 1.0 / sqrt(3.0);
 621                }
 622                else
 623                {
 624                 rotate(0.0, &x, &y);
 625                 *px = x + 5.5; *py = y + 2.0 / sqrt(3.0);
 626                }
 627                break;
 628      case 17: rotate(0.0, &x, &y);
 629               *px = x + 1.0; *py = y + 1.0 / (2.0 * sqrt(3.0)); break;
 630      case 18: rotate(120.0, &x, &y);
 631               *px = x + 4.0; *py = y + 1.0 / (2.0 * sqrt(3.0)); break;
 632      case 19: rotate(120.0, &x, &y);
 633               *px = x + 4.5; *py = y + 2.0 / sqrt(3.0); break;
 634      case 20: rotate(300.0, &x, &y);
 635               *px = x + 5.0; *py = y + 5.0 / (2.0 * sqrt(3.0)); break;
 636
 637    } /* end switch statement */
 638
 639 } /* end of dymax_point */
 640
 641
 642 static void rotate(double angle, double *x, double *y)
 643 {
 644   /* Rotate the point to correct orientation in XY-plane. */
 645
 646   double ha, hx, hy ;
 647
 648   ha = radians(angle);
 649   hx = *x;
 650   hy = *y;
 651   *x = hx * cos(ha) - hy * sin(ha);
 652   *y = hx * sin(ha) + hy * cos(ha);
 653
 654 } /* end rotate procedure */
 655
 656
 657 static void r2(int axis, double alpha, double *x, double *y, double *z)
 658 {
 659   /* Rotate a 3-D point about the specified axis.         */
 660
 661   double a, b, c;
 662
 663   a = *x;
 664   b = *y;
 665   c = *z;
 666   if (axis == 1)
 667   {
 668    *y = b * cos(alpha) + c * sin(alpha);
 669    *z = c * cos(alpha) - b * sin(alpha);
 670   }
 671
 672   if (axis == 2)
 673   {
 674    *x = a * cos(alpha) - c * sin(alpha);
 675    *z = a * sin(alpha) + c * cos(alpha);
 676   }
 677
 678   if (axis == 3)
 679   {
 680    *x = a * cos(alpha) + b * sin(alpha);
 681    *y = b * cos(alpha) - a * sin(alpha);
 682   }
 683
 684 } /* end of r2 */
 685