pd.if.org Git - isea/blob - pfuncs.c

   1 #ifndef ISEA_STATIC
   2 #define ISEA_STATIC
   3 #endif
   4
   5 enum isea_poly { ISEA_NONE, ISEA_ICOSAHEDRON = 20 };
   6 enum isea_topology { ISEA_HEXAGON=6, ISEA_TRIANGLE=3, ISEA_DIAMOND=4 };
   7 enum isea_address_form { ISEA_GEO, ISEA_Q2DI, ISEA_SEQNUM, ISEA_INTERLEAVE,
   8         ISEA_PLANE, ISEA_Q2DD, ISEA_PROJTRI, ISEA_VERTEX2DD, ISEA_HEX
   9 };
  10
  11 struct isea_dgg {
  12         int     polyhedron; /* ignored, icosahedron */
  13         double  o_lat, o_lon, o_az; /* orientation, radians */
  14         int     pole; /* true if standard snyder */
  15         int     topology; /* ignored, hexagon */
  16         int     aperture; /* valid values depend on partitioning method */
  17         int     resolution;
  18         double  radius; /* radius of the earth in meters, ignored 1.0 */
  19         int     output; /* an isea_address_form */
  20         int     triangle; /* triangle of last transformed point */
  21         int     quad; /* quad of last transformed point */
  22         unsigned long serial;
  23 };
  24
  25 struct isea_pt {
  26         double x, y;
  27 };
  28
  29 struct isea_geo {
  30         double lon, lat;
  31 };
  32
  33 struct isea_address {
  34         int     type; /* enum isea_address_form */
  35         int     number;
  36         double  x,y; /* or i,j or lon,lat depending on type */
  37 };
  38
  39 /* ENDINC */
  40
  41 enum snyder_polyhedron {
  42         SNYDER_POLY_HEXAGON, SNYDER_POLY_PENTAGON,
  43         SNYDER_POLY_TETRAHEDRON, SNYDER_POLY_CUBE,
  44         SNYDER_POLY_OCTAHEDRON, SNYDER_POLY_DODECAHEDRON,
  45         SNYDER_POLY_ICOSAHEDRON
  46 };
  47
  48 struct snyder_constants {
  49         double          g, G, theta, ea_w, ea_a, ea_b, g_w, g_a, g_b;
  50 };
  51
  52 /* TODO put these in radians to avoid a later conversion */
  53 ISEA_STATIC
  54 struct snyder_constants constants[] = {
  55         {23.80018260, 62.15458023, 60.0, 3.75, 1.033, 0.968, 5.09, 1.195, 1.0},
  56         {20.07675127, 55.69063953, 54.0, 2.65, 1.030, 0.983, 3.59, 1.141, 1.027},
  57         {0.0},
  58         {0.0},
  59         {0.0},
  60         {0.0},
  61         {37.37736814, 36.0, 30.0, 17.27, 1.163, 0.860, 13.14, 1.584, 1.0},
  62 };
  63
  64 #define E 52.62263186
  65 #define F 10.81231696
  66
  67 #define DEG60 1.04719755119659774614
  68 #define DEG120 2.09439510239319549229
  69 #define DEG72 1.25663706143591729537
  70 #define DEG90 1.57079632679489661922
  71 #define DEG144 2.51327412287183459075
  72 #define DEG36 0.62831853071795864768
  73 #define DEG108 1.88495559215387594306
  74 #define DEG180 M_PI
  75 /* sqrt(5)/M_PI */
  76 #define ISEA_SCALE 0.8301572857837594396028083
  77
  78 /* 26.565051177 degrees */
  79 #define V_LAT 0.46364760899944494524
  80
  81 #define RAD2DEG (180.0/M_PI)
  82 #define DEG2RAD (M_PI/180.0)
  83
  84 ISEA_STATIC
  85 struct isea_geo vertex[] = {
  86         {0.0, DEG90},
  87         {DEG180, V_LAT},
  88         {-DEG108, V_LAT},
  89         {-DEG36, V_LAT},
  90         {DEG36, V_LAT},
  91         {DEG108, V_LAT},
  92         {-DEG144, -V_LAT},
  93         {-DEG72, -V_LAT},
  94         {0.0, -V_LAT},
  95         {DEG72, -V_LAT},
  96         {DEG144, -V_LAT},
  97         {0.0, -DEG90}
  98 };
  99
 100 /* TODO make an isea_pt array of the vertices as well */
 101
 102 static int      tri_v1[] = {0, 0, 0, 0, 0, 0, 6, 7, 8, 9, 10, 2, 3, 4, 5, 1, 11, 11, 11, 11, 11};
 103
 104 /* 52.62263186 */
 105 #define E_RAD 0.91843818702186776133
 106
 107 /* 10.81231696 */
 108 #define F_RAD 0.18871053072122403508
 109
 110 /* triangle Centers */
 111 struct isea_geo icostriangles[] = {
 112         {0.0, 0.0},
 113         {-DEG144, E_RAD},
 114         {-DEG72, E_RAD},
 115         {0.0, E_RAD},
 116         {DEG72, E_RAD},
 117         {DEG144, E_RAD},
 118         {-DEG144, F_RAD},
 119         {-DEG72, F_RAD},
 120         {0.0, F_RAD},
 121         {DEG72, F_RAD},
 122         {DEG144, F_RAD},
 123         {-DEG108, -F_RAD},
 124         {-DEG36, -F_RAD},
 125         {DEG36, -F_RAD},
 126         {DEG108, -F_RAD},
 127         {DEG180, -F_RAD},
 128         {-DEG108, -E_RAD},
 129         {-DEG36, -E_RAD},
 130         {DEG36, -E_RAD},
 131         {DEG108, -E_RAD},
 132         {DEG180, -E_RAD},
 133 };
 134
 135 static double
 136 az_adjustment(int triangle)
 137 {
 138         double          adj;
 139
 140         struct isea_geo v;
 141         struct isea_geo c;
 142
 143         v = vertex[tri_v1[triangle]];
 144         c = icostriangles[triangle];
 145
 146         /* TODO looks like the adjustment is always either 0 or 180 */
 147         /* at least if you pick your vertex carefully */
 148         adj = atan2(cos(v.lat) * sin(v.lon - c.lon),
 149                     cos(c.lat) * sin(v.lat)
 150                     - sin(c.lat) * cos(v.lat) * cos(v.lon - c.lon));
 151         return adj;
 152 }
 153
 154 /* R tan(g) sin(60) */
 155 #define TABLE_G 0.6615845383
 156
 157 /* H = 0.25 R tan g = */
 158 #define TABLE_H 0.1909830056
 159
 160 #define RPRIME 0.91038328153090290025
 161
 162 ISEA_STATIC
 163 struct isea_pt
 164 isea_triangle_xy(int triangle)
 165 {
 166         struct isea_pt  c;
 167         double Rprime = 0.91038328153090290025;
 168
 169         triangle = (triangle - 1) % 20;
 170
 171         c.x = TABLE_G * ((triangle % 5) - 2) * 2.0;
 172         if (triangle > 9) {
 173                 c.x += TABLE_G;
 174         }
 175         switch (triangle / 5) {
 176         case 0:
 177                 c.y = 5.0 * TABLE_H;
 178                 break;
 179         case 1:
 180                 c.y = TABLE_H;
 181                 break;
 182         case 2:
 183                 c.y = -TABLE_H;
 184                 break;
 185         case 3:
 186                 c.y = -5.0 * TABLE_H;
 187                 break;
 188         default:
 189                 /* should be impossible */
 190                 exit(EXIT_FAILURE);
 191         };
 192         c.x *= Rprime;
 193         c.y *= Rprime;
 194
 195         return c;
 196 }
 197
 198 /* snyder eq 14 */
 199 static double
 200 sph_azimuth(double f_lon, double f_lat, double t_lon, double t_lat)
 201 {
 202         double          az;
 203
 204         az = atan2(cos(t_lat) * sin(t_lon - f_lon),
 205                    cos(f_lat) * sin(t_lat)
 206                    - sin(f_lat) * cos(t_lat) * cos(t_lon - f_lon)
 207                 );
 208         return az;
 209 }
 210
 211 /* coord needs to be in radians */
 212 ISEA_STATIC
 213 int
 214 isea_snyder_forward(struct isea_geo * ll, struct isea_pt * out)
 215 {
 216         int             i;
 217
 218         /*
 219          * spherical distance from center of polygon face to any of its
 220          * vertexes on the globe
 221          */
 222         double          g;
 223
 224         /*
 225          * spherical angle between radius vector to center and adjacent edge
 226          * of spherical polygon on the globe
 227          */
 228         double          G;
 229
 230         /*
 231          * plane angle between radius vector to center and adjacent edge of
 232          * plane polygon
 233          */
 234         double          theta;
 235
 236         /* additional variables from snyder */
 237         double          q, Rprime, H, Ag, Azprime, Az, dprime, f, rho,
 238                         x, y;
 239
 240         /* variables used to store intermediate results */
 241         double          cot_theta, tan_g, az_offset;
 242
 243         /* how many multiples of 60 degrees we adjust the azimuth */
 244         int             Az_adjust_multiples;
 245
 246         struct snyder_constants c;
 247
 248         /*
 249          * TODO by locality of reference, start by trying the same triangle
 250          * as last time
 251          */
 252
 253         /* TODO put these constants in as radians to begin with */
 254         c = constants[SNYDER_POLY_ICOSAHEDRON];
 255         theta = c.theta * DEG2RAD;
 256         g = c.g * DEG2RAD;
 257         G = c.G * DEG2RAD;
 258
 259         for (i = 1; i <= 20; i++) {
 260                 double          z;
 261                 struct isea_geo center;
 262
 263                 center = icostriangles[i];
 264
 265                 /* step 1 */
 266 #if 0
 267                 z = sph_distance(center.lon, center.lat, ll->lon, ll->lat);
 268 #else
 269                 z = acos(sin(center.lat) * sin(ll->lat)
 270                          + cos(center.lat) * cos(ll->lat) * cos(ll->lon - center.lon));
 271 #endif
 272
 273                 /* not on this triangle */
 274                 if (z > g + 0.000005) { /* TODO DBL_EPSILON */
 275                         continue;
 276                 }
 277                 Az = sph_azimuth(ll->lon, ll->lat, center.lon, center.lat);
 278
 279                 Az = atan2(cos(ll->lat) * sin(ll->lon - center.lon),
 280                            cos(center.lat) * sin(ll->lat)
 281                 - sin(center.lat) * cos(ll->lat) * cos(ll->lon - center.lon)
 282                         );
 283
 284                 /* step 2 */
 285
 286                 /* This calculates "some" vertex coordinate */
 287                 az_offset = az_adjustment(i);
 288
 289                 Az -= az_offset;
 290
 291                 /* TODO I don't know why we do this.  It's not in snyder */
 292                 /* maybe because we should have picked a better vertex */
 293                 if (Az < 0.0) {
 294                         Az += 2.0 * M_PI;
 295                 }
 296                 /*
 297                  * adjust Az for the point to fall within the range of 0 to
 298                  * 2(90 - theta) or 60 degrees for the hexagon, by
 299                  * and therefore 120 degrees for the triangle
 300                  * of the icosahedron
 301                  * subtracting or adding multiples of 60 degrees to Az and
 302                  * recording the amount of adjustment
 303                  */
 304
 305                 Az_adjust_multiples = 0;
 306                 while (Az < 0.0) {
 307                         Az += DEG120;
 308                         Az_adjust_multiples--;
 309                 }
 310                 while (Az > DEG120 + DBL_EPSILON) {
 311                         Az -= DEG120;
 312                         Az_adjust_multiples++;
 313                 }
 314
 315                 /* step 3 */
 316                 cot_theta = 1.0 / tan(theta);
 317                 tan_g = tan(g); /* TODO this is a constant */
 318
 319                 /* Calculate q from eq 9. */
 320                 /* TODO cot_theta is cot(30) */
 321                 q = atan2(tan_g, cos(Az) + sin(Az) * cot_theta);
 322
 323                 /* not in this triangle */
 324                 if (z > q + 0.000005) {
 325                         continue;
 326                 }
 327                 /* step 4 */
 328
 329                 /* Apply equations 5-8 and 10-12 in order */
 330
 331                 /* eq 5 */
 332                 /* Rprime = 0.9449322893 * R; */
 333                 /* R' in the paper is for the truncated */
 334                 Rprime = 0.91038328153090290025;
 335
 336                 /* eq 6 */
 337                 H = acos(sin(Az) * sin(G) * cos(g) - cos(Az) * cos(G));
 338
 339                 /* eq 7 */
 340                 /* Ag = (Az + G + H - DEG180) * M_PI * R * R / DEG180; */
 341                 Ag = Az + G + H - DEG180;
 342
 343                 /* eq 8 */
 344                 Azprime = atan2(2.0 * Ag, Rprime * Rprime * tan_g * tan_g - 2.0 * Ag * cot_theta);
 345
 346                 /* eq 10 */
 347                 /* cot(theta) = 1.73205080756887729355 */
 348                 dprime = Rprime * tan_g / (cos(Azprime) + sin(Azprime) * cot_theta);
 349
 350                 /* eq 11 */
 351                 f = dprime / (2.0 * Rprime * sin(q / 2.0));
 352
 353                 /* eq 12 */
 354                 rho = 2.0 * Rprime * f * sin(z / 2.0);
 355
 356                 /*
 357                  * add back the same 60 degree multiple adjustment from step
 358                  * 2 to Azprime
 359                  */
 360
 361                 Azprime += DEG120 * Az_adjust_multiples;
 362
 363                 /* calculate rectangular coordinates */
 364
 365                 x = rho * sin(Azprime);
 366                 y = rho * cos(Azprime);
 367
 368                 /*
 369                  * TODO
 370                  * translate coordinates to the origin for the particular
 371                  * hexagon on the flattened polyhedral map plot
 372                  */
 373
 374                 out->x = x;
 375                 out->y = y;
 376
 377                 return i;
 378         }
 379
 380         /*
 381          * should be impossible, this implies that the coordinate is not on
 382          * any triangle
 383          */
 384
 385         fprintf(stderr, "impossible transform: %f %f is not on any triangle\n",
 386                 ll->lon * RAD2DEG, ll->lat * RAD2DEG);
 387
 388         exit(EXIT_FAILURE);
 389
 390         /* not reached */
 391         return 0;               /* supresses a warning */
 392 }
 393
 394 /*
 395  * return the new coordinates of any point in orginal coordinate system.
 396  * Define a point (newNPold) in orginal coordinate system as the North Pole in
 397  * new coordinate system, and the great circle connect the original and new
 398  * North Pole as the lon0 longitude in new coordinate system, given any point
 399  * in orginal coordinate system, this function return the new coordinates.
 400  */
 401
 402 #define PRECISION 0.0000000000005
 403
 404 /*
 405  * TODO rename to isea_coordtrans and take an struct pointer to the result
 406  * struct TODO remove function entirely.  we use the synder one. the only
 407  * difference is a 180 degree longitude rotation
 408  */
 409 ISEA_STATIC
 410 struct isea_geo
 411 coordtrans(struct isea_geo * newNPold, struct isea_geo * ptold, double lon0)
 412 {
 413         double          cosptnewlat, cosptnewlon;
 414         struct isea_geo ptnew;
 415
 416         cosptnewlat = sin(newNPold->lat) * sin(ptold->lat) +
 417                 cos(newNPold->lat) * cos(ptold->lat) * cos(newNPold->lon - ptold->lon);
 418
 419         if (cosptnewlat > 1.)
 420                 cosptnewlat = 1.0;
 421         if (cosptnewlat < -1.)
 422                 cosptnewlat = -1.0;
 423
 424         ptnew.lat = acos(cosptnewlat);
 425
 426         if (fabs(ptnew.lat - 0.0) < PRECISION * 100000)
 427                 ptnew.lon = 0.0;
 428         else if (fabs(ptnew.lat - M_PI) < PRECISION * 100000)
 429                 ptnew.lon = 0.0;
 430         else {
 431                 cosptnewlon = (sin(ptold->lat) * cos(newNPold->lat) - cos(ptold->lat) *
 432                                sin(newNPold->lat) * cos(newNPold->lon - ptold->lon)) / sin(ptnew.lat);
 433                 if (cosptnewlon > 1.)
 434                         cosptnewlon = 1.0;
 435                 if (cosptnewlon < -1.)
 436                         cosptnewlon = -1.0;
 437
 438                 ptnew.lon = acos(cosptnewlon);
 439
 440                 if ((ptold->lon - newNPold->lon) >= 0 && (ptold->lon - newNPold->lon) < 180)
 441                         ptnew.lon = -ptnew.lon + lon0;
 442
 443                 else
 444                         ptnew.lon = ptnew.lon + lon0;
 445
 446                 if (ptnew.lon > M_PI)
 447                         ptnew.lon -= 2 * M_PI;
 448                 if (ptnew.lon < -M_PI)
 449                         ptnew.lon += 2 * M_PI;
 450         }
 451         ptnew.lat = M_PI / 2 - ptnew.lat;
 452         return ptnew;
 453 }
 454
 455 /* formula from Snyder, Map Projections: A working manual, p31 */
 456 /*
 457  * old north pole at np in new coordinates
 458  * could be simplified a bit with fewer intermediates
 459  *
 460  * TODO take a result pointer
 461  */
 462 ISEA_STATIC
 463 struct isea_geo
 464 snyder_ctran(struct isea_geo * np, struct isea_geo * pt)
 465 {
 466         struct isea_geo npt;
 467         double          alpha, phi, lambda, lambda0, beta, lambdap, phip;
 468         double          sin_phip;
 469         double          lp_b;   /* lambda prime minus beta */
 470         double          cos_p, sin_a;
 471
 472         phi = pt->lat;
 473         lambda = pt->lon;
 474         alpha = np->lat;
 475         beta = np->lon;
 476         lambda0 = beta;
 477
 478         cos_p = cos(phi);
 479         sin_a = sin(alpha);
 480
 481         /* mpawm 5-7 */
 482         sin_phip = sin_a * sin(phi) - cos(alpha) * cos_p * cos(lambda - lambda0);
 483
 484         /* mpawm 5-8b */
 485
 486         /* use the two argument form so we end up in the right quadrant */
 487         lp_b = atan2(cos_p * sin(lambda - lambda0),
 488            (sin_a * cos_p * cos(lambda - lambda0) + cos(alpha) * sin(phi)));
 489
 490         lambdap = lp_b + beta;
 491
 492         /* normalize longitude */
 493         /* TODO can we just do a modulus ? */
 494         lambdap = fmod(lambdap, 2 * M_PI);
 495         while (lambdap > M_PI)
 496                 lambdap -= 2 * M_PI;
 497         while (lambdap < -M_PI)
 498                 lambdap += 2 * M_PI;
 499
 500         phip = asin(sin_phip);
 501
 502         npt.lat = phip;
 503         npt.lon = lambdap;
 504
 505         return npt;
 506 }
 507
 508 ISEA_STATIC
 509 struct isea_geo
 510 isea_ctran(struct isea_geo * np, struct isea_geo * pt, double lon0)
 511 {
 512         struct isea_geo npt;
 513
 514         np->lon += M_PI;
 515         npt = snyder_ctran(np, pt);
 516         np->lon -= M_PI;
 517
 518         npt.lon -= (M_PI - lon0 + np->lon);
 519
 520         /*
 521          * snyder is down tri 3, isea is along side of tri1 from vertex 0 to
 522          * vertex 1 these are 180 degrees apart
 523          */
 524         npt.lon += M_PI;
 525         /* normalize longitude */
 526         npt.lon = fmod(npt.lon, 2 * M_PI);
 527         while (npt.lon > M_PI)
 528                 npt.lon -= 2 * M_PI;
 529         while (npt.lon < -M_PI)
 530                 npt.lon += 2 * M_PI;
 531
 532         return npt;
 533 }
 534
 535 /* in radians */
 536 #define ISEA_STD_LAT 1.01722196792335072101
 537 #define ISEA_STD_LON .19634954084936207740
 538
 539 ISEA_STATIC
 540 struct isea_dgg isea_standard_dgg = {
 541         20, 58.28252559, 11.25, 0.0,
 542         6, 4, 0, 6, 1.0
 543 };
 544
 545 /* fuller's at 5.2454 west, 2.3009 N, adjacent at 7.46658 deg */
 546
 547 ISEA_STATIC
 548 int
 549 isea_grid_init(struct isea_dgg * g)
 550 {
 551         if (!g)
 552                 return 0;
 553
 554         g->polyhedron = 20;
 555         g->o_lat = ISEA_STD_LAT;
 556         g->o_lon = ISEA_STD_LON;
 557         g->o_az = 0.0;
 558         g->aperture = 4;
 559         g->resolution = 6;
 560         g->radius = 1.0;
 561         g->topology = 6;
 562
 563         return 1;
 564 }
 565
 566 ISEA_STATIC
 567 int
 568 isea_orient_isea(struct isea_dgg * g)
 569 {
 570         if (!g)
 571                 return 0;
 572         g->o_lat = ISEA_STD_LAT;
 573         g->o_lon = ISEA_STD_LON;
 574         g->o_az = 0.0;
 575         return 1;
 576 }
 577
 578 ISEA_STATIC
 579 int
 580 isea_orient_pole(struct isea_dgg * g)
 581 {
 582         if (!g)
 583                 return 0;
 584         g->o_lat = M_PI / 2.0;
 585         g->o_lon = 0.0;
 586         g->o_az = 0;
 587         return 1;
 588 }
 589
 590 ISEA_STATIC
 591 int
 592 isea_orient_dymax(struct isea_dgg * g)
 593 {
 594         if (!g)
 595                 return 0;
 596         g->o_lat = 2.300882 * M_PI / 180.0;
 597         g->o_lon = -5.24539 * M_PI / 180.0;
 598         g->o_az = 7.46658 * M_PI / 180.0;
 599         return 1;
 600 }
 601
 602 ISEA_STATIC
 603 int
 604 isea_transform(struct isea_dgg * g, struct isea_geo * in,
 605                struct isea_pt * out)
 606 {
 607         struct isea_geo i, pole;
 608         int             tri;
 609
 610         pole.lat = g->o_lat;
 611         pole.lon = g->o_lon;
 612
 613         i = isea_ctran(&pole, in, g->o_az);
 614
 615         tri = isea_snyder_forward(&i, out);
 616         out->x *= g->radius;
 617         out->y *= g->radius;
 618         g->triangle = tri;
 619
 620         return tri;
 621 }
 622
 623 #define DOWNTRI(tri) (((tri - 1) / 5) % 2 == 1)
 624
 625 /* get which triangle a point is in, assumes radius = 1.0 */
 626 ISEA_STATIC
 627 int isea_xy_triangle(struct isea_pt *p, double radius) {
 628         double dist, test;
 629         int tri, closest;
 630         struct isea_pt tc, pt;
 631         pt.x = p->x / radius;
 632         pt.y = p->y / radius;
 633
 634         tc = isea_triangle_xy(1);
 635         dist = (pt.x - tc.x) * (pt.x - tc.x) + (pt.y - tc.y) * (pt.y - tc.y);
 636         closest = 1;
 637
 638         /* TODO just calculate it directly, no need to actually do all this
 639          * work
 640          */
 641         for (tri = 2; tri <= 20; tri++) {
 642                 tc = isea_triangle_xy(tri);
 643                 test = (pt.x - tc.x) * (pt.x - tc.x) + (pt.y - tc.y) * (pt.y - tc.y);
 644                 if (test < dist) {
 645                         dist = test;
 646                         closest = tri;
 647                 }
 648         }
 649         return closest;
 650 }
 651
 652 ISEA_STATIC
 653 int isea_projtri(struct isea_pt *xy, double radius) {
 654         int tri;
 655
 656         tri = isea_xy_triangle(xy, radius);
 657
 658         xy->x /= radius;
 659         xy->y /= radius;
 660         xy->x += 0.5;
 661         xy->y += 2.0 * .14433756729740644112;
 662
 663         return tri;
 664 }
 665
 666 ISEA_STATIC
 667 void
 668 isea_rotate(struct isea_pt * pt, double degrees)
 669 {
 670         double          rad;
 671
 672         double          x, y;
 673
 674         rad = -degrees * M_PI / 180.0;
 675         while (rad >= 2.0 * M_PI) rad -= 2.0 * M_PI;
 676         while (rad <= -2.0 * M_PI) rad += 2.0 * M_PI;
 677
 678         x = pt->x * cos(rad) + pt->y * sin(rad);
 679         y = -pt->x * sin(rad) + pt->y * cos(rad);
 680
 681         pt->x = x;
 682         pt->y = y;
 683 }
 684
 685 ISEA_STATIC
 686 int isea_tri_plane(int tri, struct isea_pt *pt, double radius) {
 687         struct isea_pt tc; /* center of triangle */
 688
 689         if (DOWNTRI(tri)) {
 690                 isea_rotate(pt, 180.0);
 691         }
 692         tc = isea_triangle_xy(tri);
 693         tc.x *= radius;
 694         tc.y *= radius;
 695         pt->x += tc.x;
 696         pt->y += tc.y;
 697
 698         return tri;
 699 }
 700
 701 ISEA_STATIC
 702 int isea_projection(struct isea_dgg *g, double xin, double yin, double *xout,
 703                 double *yout) {
 704         return 1;
 705 }
 706
 707 /* convert projected triangle coords to quad xy coords, return quad number */
 708 ISEA_STATIC
 709 int
 710 isea_ptdd(int tri, struct isea_pt *pt) {
 711         int             downtri, quad;
 712
 713         downtri = (((tri - 1) / 5) % 2 == 1);
 714         quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;
 715
 716         isea_rotate(pt, downtri ? 240.0 : 60.0);
 717         if (downtri) {
 718                 pt->x += 0.5;
 719                 /* pt->y += cos(30.0 * M_PI / 180.0); */
 720                 pt->y += .86602540378443864672;
 721         }
 722         return quad;
 723 }
 724
 725 ISEA_STATIC
 726 int
 727 isea_dddi_ap3odd(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di)
 728 {
 729         struct isea_pt  v;
 730         double          hexwidth;
 731         double          sidelength;     /* in hexes */
 732         int             d, i;
 733         int             maxcoord;
 734         struct hex      h;
 735
 736         /* This is the number of hexes from apex to base of a triangle */
 737         sidelength = (pow(2.0, g->resolution) + 1.0) / 2.0;
 738
 739         /* apex to base is cos(30deg) */
 740         hexwidth = cos(M_PI / 6.0) / sidelength;
 741
 742         /* TODO I think sidelength is always x.5, so
 743          * (int)sidelength * 2 + 1 might be just as good
 744          */
 745         maxcoord = (int) (sidelength * 2.0 + 0.5);
 746
 747         v = *pt;
 748         hexbin2(0, hexwidth, v.x, v.y, &h.x, &h.y);
 749         h.iso = 0;
 750         hex_iso(&h);
 751
 752         d = h.x - h.z;
 753         i = h.x + h.y + h.y;
 754
 755         /*
 756          * you want to test for max coords for the next quad in the same
 757          * "row" first to get the case where both are max
 758          */
 759         if (quad <= 5) {
 760                 if (d == 0 && i == maxcoord) {
 761                         /* north pole */
 762                         quad = 0;
 763                         d = 0;
 764                         i = 0;
 765                 } else if (i == maxcoord) {
 766                         /* upper right in next quad */
 767                         quad += 1;
 768                         if (quad == 6)
 769                                 quad = 1;
 770                         i = maxcoord - d;
 771                         d = 0;
 772                 } else if (d == maxcoord) {
 773                         /* lower right in quad to lower right */
 774                         quad += 5;
 775                         d = 0;
 776                 }
 777         } else if (quad >= 6) {
 778                 if (i == 0 && d == maxcoord) {
 779                         /* south pole */
 780                         quad = 11;
 781                         d = 0;
 782                         i = 0;
 783                 } else if (d == maxcoord) {
 784                         /* lower right in next quad */
 785                         quad += 1;
 786                         if (quad == 11)
 787                                 quad = 6;
 788                         d = maxcoord - i;
 789                         i = 0;
 790                 } else if (i == maxcoord) {
 791                         /* upper right in quad to upper right */
 792                         quad = (quad - 4) % 5;
 793                         i = 0;
 794                 }
 795         }
 796
 797         di->x = d;
 798         di->y = i;
 799
 800         g->quad = quad;
 801         return quad;
 802 }
 803
 804 ISEA_STATIC
 805 int
 806 isea_dddi(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di) {
 807         struct isea_pt  v;
 808         double          hexwidth;
 809         int             sidelength;     /* in hexes */
 810         struct hex      h;
 811
 812         if (g->aperture == 3 && g->resolution % 2 != 0) {
 813                 return isea_dddi_ap3odd(g, quad, pt, di);
 814         }
 815         /* todo might want to do this as an iterated loop */
 816         if (g->aperture >0) {
 817                 sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
 818         } else {
 819                 sidelength = g->resolution;
 820         }
 821
 822         hexwidth = 1.0 / sidelength;
 823
 824         v = *pt;
 825         isea_rotate(&v, -30.0);
 826         hexbin2(0, hexwidth, v.x, v.y, &h.x, &h.y);
 827         h.iso = 0;
 828         hex_iso(&h);
 829
 830         /* we may actually be on another quad */
 831         if (quad <= 5) {
 832                 if (h.x == 0 && h.z == -sidelength) {
 833                         /* north pole */
 834                         quad = 0;
 835                         h.z = 0;
 836                         h.y = 0;
 837                         h.x = 0;
 838                 } else if (h.z == -sidelength) {
 839                         quad = quad + 1;
 840                         if (quad == 6)
 841                                 quad = 1;
 842                         h.y = sidelength - h.x;
 843                         h.z = h.x - sidelength;
 844                         h.x = 0;
 845                 } else if (h.x == sidelength) {
 846                         quad += 5;
 847                         h.y = -h.z;
 848                         h.x = 0;
 849                 }
 850         } else if (quad >= 6) {
 851                 if (h.z == 0 && h.x == sidelength) {
 852                         /* south pole */
 853                         quad = 11;
 854                         h.x = 0;
 855                         h.y = 0;
 856                         h.z = 0;
 857                 } else if (h.x == sidelength) {
 858                         quad = quad + 1;
 859                         if (quad == 11)
 860                                 quad = 6;
 861                         h.x = h.y + sidelength;
 862                         h.y = 0;
 863                         h.z = -h.x;
 864                 } else if (h.y == -sidelength) {
 865                         quad -= 4;
 866                         h.y = 0;
 867                         h.z = -h.x;
 868                 }
 869         }
 870         di->x = h.x;
 871         di->y = -h.z;
 872
 873         g->quad = quad;
 874         return quad;
 875 }
 876
 877 ISEA_STATIC
 878 int isea_ptdi(struct isea_dgg *g, int tri, struct isea_pt *pt,
 879                 struct isea_pt *di) {
 880         struct isea_pt  v;
 881         int             quad;
 882
 883         v = *pt;
 884         quad = isea_ptdd(tri, &v);
 885         quad = isea_dddi(g, quad, &v, di);
 886         return quad;
 887 }
 888
 889 /* q2di to seqnum */
 890 ISEA_STATIC
 891 int isea_disn(struct isea_dgg *g, int quad, struct isea_pt *di) {
 892         int             sidelength;
 893         int             sn, height;
 894         int             hexes;
 895
 896         if (quad == 0) {
 897                 g->serial = 1;
 898                 return g->serial;
 899         }
 900         /* hexes in a quad */
 901         hexes = (int) (pow(g->aperture, g->resolution) + 0.5);
 902         if (quad == 11) {
 903                 g->serial = 1 + 10 * hexes + 1;
 904                 return g->serial;
 905         }
 906         if (g->aperture == 3 && g->resolution % 2 == 1) {
 907                 height = (int) (pow(g->aperture, (g->resolution - 1) / 2.0));
 908                 sn = ((int) di->x) * height;
 909                 sn += ((int) di->y) / height;
 910                 sn += (quad - 1) * hexes;
 911                 sn += 2;
 912         } else {
 913                 sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
 914                 sn = (quad - 1) * hexes + sidelength * di->x + di->y + 2;
 915         }
 916
 917         g->serial = sn;
 918         return sn;
 919 }
 920
 921 /* TODO just encode the quad in the d or i coordinate
 922  * quad is 0-11, which can be four bits.
 923  * d' = d << 4 + q, d = d' >> 4, q = d' & 0xf
 924  */
 925 /* convert a q2di to global hex coord */
 926 ISEA_STATIC
 927 int isea_hex(struct isea_dgg *g, int tri,
 928                 struct isea_pt *pt, struct isea_pt *hex) {
 929         struct isea_pt v;
 930         int sidelength;
 931         int d, i, x, y, quad;
 932         double oddsl;
 933
 934         quad = isea_ptdi(g, tri, pt, &v);
 935
 936         hex->x = ((int)v.x << 4) + quad;
 937         hex->y = v.y;
 938
 939         return 1;
 940
 941         d = v.x;
 942         i = v.y;
 943
 944         /* Aperture 3 odd resolutions */
 945         if (g->aperture == 3 && g->resolution % 2 != 0) {
 946                 int offset = (int)(pow(3.0, g->resolution - 1) + 0.5);
 947
 948                 oddsl = (pow(2.0, g->resolution) + 1.0) / 2.0;
 949
 950                 d += offset * ((g->quad-1) % 5);
 951                 i += offset * ((g->quad-1) % 5);
 952
 953                 if (quad == 0) {
 954                         d = 0;
 955                         i = offset;
 956                 } else if (quad == 11) {
 957                         d = 2 * offset;
 958                         i = 0;
 959                 } else if (quad > 5) {
 960                         d += offset;
 961                 }
 962
 963                 x = (2*d - i) /3;
 964                 y = (2*i - d) /3;
 965
 966                 hex->x = x + offset / 3;
 967                 hex->y = y + 2 * offset / 3;
 968                 return 1;
 969         }
 970
 971         /* aperture 3 even resolutions and aperture 4 */
 972         sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
 973         if (g->quad == 0) {
 974                 hex->x = 0;
 975                 hex->y = sidelength;
 976         } else if (g->quad == 11) {
 977                 hex->x = sidelength * 2;
 978                 hex->y = 0;
 979         } else {
 980                 hex->x = d + sidelength * ((g->quad-1) % 5);
 981                 if (g->quad > 5) hex->x += sidelength;
 982                 hex->y = i + sidelength * ((g->quad-1) % 5);
 983         }
 984
 985         return 1;
 986 }
 987
 988 ISEA_STATIC
 989 struct isea_pt
 990 isea_forward(struct isea_dgg *g, struct isea_geo *in)
 991 {
 992         int             tri, downtri, quad;
 993         struct isea_pt  out, coord;
 994
 995         tri = isea_transform(g, in, &out);
 996
 997         downtri = (((tri - 1) / 5) % 2 == 1);
 998         quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;
 999
1000         if (g->output == ISEA_PLANE) {
1001                 isea_tri_plane(tri, &out, g->radius);
1002                 return out;
1003         }
1004
1005         /* convert to isea standard triangle size */
1006         out.x = out.x / g->radius * ISEA_SCALE;
1007         out.y = out.y / g->radius * ISEA_SCALE;
1008         out.x += 0.5;
1009         out.y += 2.0 * .14433756729740644112;
1010
1011         switch (g->output) {
1012         case ISEA_PROJTRI:
1013                 /* nothing to do, already in projected triangle */
1014                 break;
1015         case ISEA_VERTEX2DD:
1016                 g->quad = isea_ptdd(tri, &out);
1017                 break;
1018         case ISEA_Q2DD:
1019                 /* Same as above, we just don't print as much */
1020                 g->quad = isea_ptdd(tri, &out);
1021                 break;
1022         case ISEA_Q2DI:
1023                 g->quad = isea_ptdi(g, tri, &out, &coord);
1024                 return coord;
1025                 break;
1026         case ISEA_SEQNUM:
1027                 isea_ptdi(g, tri, &out, &coord);
1028                 /* disn will set g->serial */
1029                 isea_disn(g, g->quad, &coord);
1030                 return coord;
1031                 break;
1032         case ISEA_HEX:
1033                 isea_hex(g, tri, &out, &coord);
1034                 return coord;
1035                 break;
1036         }
1037
1038         return out;
1039 }