]> pd.if.org Git - isea/blob - pfuncs.c
Initial checkin.
[isea] / 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 }