--- /dev/null
+CFLAGS=-Wall -I.. -I.
+LDFLAGS=-lm
+TESTS= t/hextest.t t/hexbin.t t/projtri.t t/q2dd.t t/q2di.t t/q2di_r3.t \
+ t/seqnum.t t/globehex.t t/plane.t
+
+all: test transform PJ_isea.c
+
+isea.o: isea.c pfuncs.c hfuncs.c
+
+transform: transform.o isea.o
+
+hextest: hextest.o hex.o tap.o
+
+t/%.t: t/%.o t/tap.o isea.o
+ $(CC) -o $@ $+ $(LDFLAGS)
+
+test: $(TESTS)
+ @prove --exec ''
+
+snyder.o: snyder.c
+ cc -c $(CFLAGS) -o $@ $+
+
+PJ_isea.c: pjhead.c hfuncs.c pfuncs.c pjtail.c
+ cat $+ > $@
+
+clean:
+ rm -f isea *.o gshhs.pdf *.mp transform PJ_isea.c
+ rm -f t/*.o t/*.t
+
+gshhs_isea.mp: transform gshhs_c.asc
+ ./transform -O plane < gshhs_c.asc |./mkmp > $@
+
+gshhs_isear.mp: transform gshhs_c.asc
+ ./transform -p 58.28252559 -l -168.75 -z 0 -O plane < gshhs_c.asc |./mkmp > $@
+
+gshhs_isear2.mp: transform gshhs_c.asc
+ ./transform -p 58.28252559 -l -168.75 -z 144 -O plane < gshhs_c.asc |./mkmp > $@
+
+gshhs_isear3.mp: transform gshhs_c.asc
+ ./transform -p 58.28252559 -l -168.75 -z 72 -O plane < gshhs_c.asc |./mkmp > $@
+
+gshhs_isear4.mp: transform gshhs_c.asc
+ ./transform -p 58.28252559 -l -168.75 -z -144 -O plane < gshhs_c.asc |./mkmp > $@
+
+gshhs_isear5.mp: transform gshhs_c.asc
+ ./transform -p 58.28252559 -l -168.75 -z -72 -O plane < gshhs_c.asc |./mkmp > $@
+
+gshhs_sym1.mp: transform gshhs_c.asc
+ ./transform -p 48.28252559 -l 106 -z -100 -O plane < gshhs_c.asc |./mkmp > $@
+
+gshhs_sym2.mp: transform gshhs_c.asc
+ ./transform -p 58.28252559 -l 124 -z -77 -O plane < gshhs_c.asc |./mkmp > $@
+
+gshhs_sym3.mp: transform gshhs_c.asc
+ ./transform -p 58.28252559 -l 124 -z -92 -O plane < gshhs_c.asc |./mkmp > $@
+
+gshhs_snyder.mp: transform gshhs_c.asc
+ ./transform -P -O plane < gshhs_c.asc |./mkmp > $@
+
+gshhs_world.mp: transform gshhs_c.asc
+ ./transform -p 90 -l 0 -z -22.5 -O plane < gshhs_c.asc |./mkmp > $@
+
+gshhs_dymax.mp: transform gshhs_c.asc
+ ./transform -p 2.300882 -l -5.24539 -z 7.46658 -O plane < gshhs_c.asc |./mkmp > $@
+
+gshhs_eurasia.mp: transform gshhs_c.asc
+ #./transform -l -145 -p 49 -z 160 -O plane < gshhs_c.asc |./mkmp > $@
+ ./transform -l -155 -p 55 -z -174 -O plane < gshhs_c.asc |./mkmp > $@
+
+gshhs_eurquad.mp: transform gshhs_c.asc
+ #./transform -l -40 -p 35 -z -118 -O plane < gshhs_c.asc |./mkmp > $@
+ ./transform -l 90 -p 55 -z 35 -O plane < gshhs_c.asc |./mkmp > $@
+
+gshhs_dymax2.mp:
+ ./transform -p 50.103201 -l -143.47849 -z 156.08388 -O plane < gshhs_c.asc |./mkmp > $@
+
+MP=gshhs_eurasia.mp gshhs_isea.mp gshhs_dymax.mp \
+ gshhs_snyder.mp gshhs_eurquad.mp gshhs_world.mp \
+ gshhs_isear.mp gshhs_isear2.mp gshhs_isear3.mp gshhs_isear4.mp \
+ gshhs_isear5.mp gshhs_dymax2.mp\
+
+gshhs.pdf: $(MP) gshhs.tex
+ for i in $(MP); do mpost $$i; done;
+ pdflatex gshhs.tex
--- /dev/null
+Adding isea to proj4
+
+make PJ_isea.c
+cp PJ_isea.c path/to/proj-4.7.0
+cd path/to/proj-4.7.0
+patch -p1 < ../isea/proj4.patch
+
+now build proj.4 as usual for your system
+
+-- postgis
+
+There is too much numerical instability to try to project directly with postgis.
+The interruptions make it almost impossible to reconstitute the input polygons.
+
+Might be able to do it by computing segments, and then ordering the segments
+in crossing order.
+
+
+ISEA Params
+
+Aperture will set the scaling property
+
+Resolution is either the resolution level for the given aperture,
+or the total number of hexes along a triangle side if aperture is zero.
+
+I think the ISEA plane projection is offsetting to avoid negative Y coords.
--- /dev/null
+#include <stdlib.h>
+#include <math.h>
+#include <stdio.h>
+
+#include "hfuncs.c"
--- /dev/null
+#ifndef HEX_H_
+#define HEX_H_ 1
+
+struct hex {
+ int iso;
+ int x, y, z;
+};
+
+int hexbin2(int horizontal, double width, double x, double y,
+ int *i, int *j);
+
+int hexbin(int horizontal, double width, double x, double y,
+ int *i, int *j);
+int hex_xy(struct hex *h);
+int hex_iso(struct hex *h);
+
+int hex_distance(struct hex *a, struct hex *b);
+int hex_bearing(struct hex *a, struct hex *b);
+int hex_at(struct hex *h, int distance, struct hex **list);
+int hex_within(struct hex *h, int distance, struct hex **list);
+int hex_dir(struct hex *h, int direction);
+
+#endif
--- /dev/null
+#ifndef ISEA_STATIC
+#define ISEA_STATIC
+#endif
+
+struct hex {
+ int iso;
+ int x, y, z;
+};
+
+/* y *must* be positive down as the xy /iso conversion assumes this */
+ISEA_STATIC
+int hex_xy(struct hex *h) {
+ if (!h->iso) return 1;
+ if (h->x >= 0) {
+ h->y = -h->y - (h->x+1)/2;
+ } else {
+ /* need to round toward -inf, not toward zero, so x-1 */
+ h->y = -h->y - h->x/2;
+ }
+ h->iso = 0;
+
+ return 1;
+}
+
+ISEA_STATIC
+int hex_iso(struct hex *h) {
+ if (h->iso) return 1;
+
+ if (h->x >= 0) {
+ h->y = (-h->y - (h->x+1)/2);
+ } else {
+ /* need to round toward -inf, not toward zero, so x-1 */
+ h->y = (-h->y - (h->x)/2);
+ }
+
+ h->z = -h->x - h->y;
+ h->iso = 1;
+ return 1;
+}
+
+ISEA_STATIC
+int hexbin2(int horizontal, double width, double x, double y,
+ int *i, int *j) {
+ double z, rx, ry, rz;
+ double abs_dx, abs_dy, abs_dz;
+ int ix, iy, iz, s;
+ struct hex h;
+
+ x = x / cos(30 * M_PI / 180.0); /* rotated X coord */
+ y = y - x / 2.0; /* adjustment for rotated X */
+
+ /* adjust for actual hexwidth */
+ x /= width;
+ y /= width;
+
+ z = -x - y;
+
+ ix = rx = floor(x + 0.5);
+ iy = ry = floor(y + 0.5);
+ iz = rz = floor(z + 0.5);
+
+ s = ix + iy + iz;
+
+ if (s) {
+ abs_dx = fabs(rx - x);
+ abs_dy = fabs(ry - y);
+ abs_dz = fabs(rz - z);
+
+ if (abs_dx >= abs_dy && abs_dx >= abs_dz) {
+ ix -= s;
+ } else if (abs_dy >= abs_dx && abs_dy >= abs_dz) {
+ iy -= s;
+ } else {
+ iz -= s;
+ }
+ }
+ h.x = ix;
+ h.y = iy;
+ h.z = iz;
+ h.iso = 1;
+
+ hex_xy(&h);
+ *i = h.x;
+ *j = h.y;
+ return ix * 100 + iy;
+}
--- /dev/null
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <float.h>
+
+#include "hfuncs.c"
+#include "pfuncs.c"
--- /dev/null
+#ifndef ISEA_H_
+#define ISEA_H_ 1
+
+#include "hex.h"
+
+enum isea_poly { ISEA_NONE, ISEA_ICOSAHEDRON = 20 };
+enum isea_topology { ISEA_HEXAGON=6, ISEA_TRIANGLE=3, ISEA_DIAMOND=4 };
+enum isea_address_form { ISEA_GEO, ISEA_Q2DI, ISEA_SEQNUM, ISEA_INTERLEAVE,
+ ISEA_PLANE, ISEA_Q2DD, ISEA_PROJTRI, ISEA_VERTEX2DD, ISEA_HEX
+};
+
+/* TODO calculate this */
+#define ISEA_SCALE 0.8301572857837594396028083
+
+struct isea_dgg {
+ int polyhedron; /* ignored, icosahedron */
+ double o_lat, o_lon, o_az; /* orientation, radians */
+ int pole; /* true if standard snyder */
+ int topology; /* ignored, hexagon */
+ int aperture; /* valid values depend on partitioning method */
+ int resolution;
+ double radius; /* radius of the earth in meters, ignored 1.0 */
+ int output; /* an isea_address_form */
+ int triangle; /* triangle of last transformed point */
+ int quad; /* quad of last transformed point */
+ unsigned long serial;
+};
+
+struct isea_pt {
+ double x, y;
+};
+
+struct isea_geo {
+ double lon, lat;
+};
+
+extern struct isea_dgg isea_standard_dgg;
+
+struct isea_address {
+ int type; /* enum isea_address_form */
+ int number;
+ double x,y; /* or i,j or lon,lat depending on type */
+};
+
+#if 0
+struct isea_address {
+ unsigned validity; /* which are good */
+ struct coord geo;
+ struct point plane;
+ unsigned long seqnum;
+ unsigned quad;
+ unsigned triangle;
+ struct point quad_xy;
+ struct point quad_ij;
+ struct point triangle_xy;
+ unsigned vertex;
+};
+#endif
+
+/* set the orientation */
+int isea_orient(struct isea_dgg *dgg, double lon, double lat, double az);
+int isea_orient_isea(struct isea_dgg *d);
+int isea_orient_pole(struct isea_dgg *d);
+int isea_orient_dymax(struct isea_dgg *d);
+
+/* set a clip region */
+int isea_region(struct isea_dgg *dgg, double lon, double lat);
+
+/* ensure validity */
+int isea_valid(struct isea_dgg *dgg);
+
+/* set the resolution by picking the nearest area */
+int isea_resolution_area(struct isea_dgg *dgg, double area);
+
+/* the characteristic length scale */
+double isea_cell_cls(struct isea_dgg *dgg);
+
+/* intercell distance as projected */
+double isea_cell_width(struct isea_dgg *dgg);
+
+double isea_cell_area(struct isea_dgg *dgg);
+
+/* total number of cells */
+double isea_cell_count(struct isea_dgg *dgg);
+
+/* set the resolution by picking the nearest intercell distance */
+int isea_resolution_distance(struct isea_dgg *dgg, double distance);
+
+/* grid size in cells */
+int isea_grid_size(struct isea_dgg *dgg);
+
+int isea_grid_init(struct isea_dgg *d);
+
+/* const? */
+struct isea_dgg *isea_grid_standard(void);
+
+/* coordinate projection */
+
+/* returns triangle, out may be null */
+int isea_transform(struct isea_dgg *dgg, struct isea_geo *in,
+ struct isea_pt *out);
+
+struct isea_pt isea_forward(struct isea_dgg *g, struct isea_geo *in);
+
+/* returns seqnum */
+int isea_readdress(struct isea_dgg *dgg, int input, struct isea_pt *in,
+ int output, struct isea_pt *out);
+
+struct isea_pt isea_triangle_xy(int triangle);
+int isea_xy_triangle(struct isea_pt *p, double radius);
+
+/* grid changes */
+/*
+ * the idea is that if you have a cartesian coordinate, and just want
+ * to change the resolution, it's just a re-bin, not a new
+ * projection.
+ */
+int isea_tranform_grid(void);
+
+/* grid generation */
+
+/* address conversion */
+void isea_rotate(struct isea_pt *pt, double degrees);
+int isea_dddi(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di);
+int isea_ptdd(int triangle, struct isea_pt *pt);
+int isea_ptdi(struct isea_dgg *g, int tri, struct isea_pt *pt, struct isea_pt *di);
+int isea_hex(struct isea_dgg *g, int tri, struct isea_pt *pt, struct isea_pt *di);
+int isea_disn(struct isea_dgg *g, int quad, struct isea_pt *h);
+
+/* binning point values */
+
+/* presence/absence binning */
+int isea_snyder_forward(struct isea_geo *ll, struct isea_pt *out);
+
+#endif
--- /dev/null
+#ifndef ISEA_STATIC
+#define ISEA_STATIC
+#endif
+
+enum isea_poly { ISEA_NONE, ISEA_ICOSAHEDRON = 20 };
+enum isea_topology { ISEA_HEXAGON=6, ISEA_TRIANGLE=3, ISEA_DIAMOND=4 };
+enum isea_address_form { ISEA_GEO, ISEA_Q2DI, ISEA_SEQNUM, ISEA_INTERLEAVE,
+ ISEA_PLANE, ISEA_Q2DD, ISEA_PROJTRI, ISEA_VERTEX2DD, ISEA_HEX
+};
+
+struct isea_dgg {
+ int polyhedron; /* ignored, icosahedron */
+ double o_lat, o_lon, o_az; /* orientation, radians */
+ int pole; /* true if standard snyder */
+ int topology; /* ignored, hexagon */
+ int aperture; /* valid values depend on partitioning method */
+ int resolution;
+ double radius; /* radius of the earth in meters, ignored 1.0 */
+ int output; /* an isea_address_form */
+ int triangle; /* triangle of last transformed point */
+ int quad; /* quad of last transformed point */
+ unsigned long serial;
+};
+
+struct isea_pt {
+ double x, y;
+};
+
+struct isea_geo {
+ double lon, lat;
+};
+
+struct isea_address {
+ int type; /* enum isea_address_form */
+ int number;
+ double x,y; /* or i,j or lon,lat depending on type */
+};
+
+/* ENDINC */
+
+enum snyder_polyhedron {
+ SNYDER_POLY_HEXAGON, SNYDER_POLY_PENTAGON,
+ SNYDER_POLY_TETRAHEDRON, SNYDER_POLY_CUBE,
+ SNYDER_POLY_OCTAHEDRON, SNYDER_POLY_DODECAHEDRON,
+ SNYDER_POLY_ICOSAHEDRON
+};
+
+struct snyder_constants {
+ double g, G, theta, ea_w, ea_a, ea_b, g_w, g_a, g_b;
+};
+
+/* TODO put these in radians to avoid a later conversion */
+ISEA_STATIC
+struct snyder_constants constants[] = {
+ {23.80018260, 62.15458023, 60.0, 3.75, 1.033, 0.968, 5.09, 1.195, 1.0},
+ {20.07675127, 55.69063953, 54.0, 2.65, 1.030, 0.983, 3.59, 1.141, 1.027},
+ {0.0},
+ {0.0},
+ {0.0},
+ {0.0},
+ {37.37736814, 36.0, 30.0, 17.27, 1.163, 0.860, 13.14, 1.584, 1.0},
+};
+
+#define E 52.62263186
+#define F 10.81231696
+
+#define DEG60 1.04719755119659774614
+#define DEG120 2.09439510239319549229
+#define DEG72 1.25663706143591729537
+#define DEG90 1.57079632679489661922
+#define DEG144 2.51327412287183459075
+#define DEG36 0.62831853071795864768
+#define DEG108 1.88495559215387594306
+#define DEG180 M_PI
+/* sqrt(5)/M_PI */
+#define ISEA_SCALE 0.8301572857837594396028083
+
+/* 26.565051177 degrees */
+#define V_LAT 0.46364760899944494524
+
+#define RAD2DEG (180.0/M_PI)
+#define DEG2RAD (M_PI/180.0)
+
+ISEA_STATIC
+struct isea_geo vertex[] = {
+ {0.0, DEG90},
+ {DEG180, V_LAT},
+ {-DEG108, V_LAT},
+ {-DEG36, V_LAT},
+ {DEG36, V_LAT},
+ {DEG108, V_LAT},
+ {-DEG144, -V_LAT},
+ {-DEG72, -V_LAT},
+ {0.0, -V_LAT},
+ {DEG72, -V_LAT},
+ {DEG144, -V_LAT},
+ {0.0, -DEG90}
+};
+
+/* TODO make an isea_pt array of the vertices as well */
+
+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};
+
+/* 52.62263186 */
+#define E_RAD 0.91843818702186776133
+
+/* 10.81231696 */
+#define F_RAD 0.18871053072122403508
+
+/* triangle Centers */
+struct isea_geo icostriangles[] = {
+ {0.0, 0.0},
+ {-DEG144, E_RAD},
+ {-DEG72, E_RAD},
+ {0.0, E_RAD},
+ {DEG72, E_RAD},
+ {DEG144, E_RAD},
+ {-DEG144, F_RAD},
+ {-DEG72, F_RAD},
+ {0.0, F_RAD},
+ {DEG72, F_RAD},
+ {DEG144, F_RAD},
+ {-DEG108, -F_RAD},
+ {-DEG36, -F_RAD},
+ {DEG36, -F_RAD},
+ {DEG108, -F_RAD},
+ {DEG180, -F_RAD},
+ {-DEG108, -E_RAD},
+ {-DEG36, -E_RAD},
+ {DEG36, -E_RAD},
+ {DEG108, -E_RAD},
+ {DEG180, -E_RAD},
+};
+
+static double
+az_adjustment(int triangle)
+{
+ double adj;
+
+ struct isea_geo v;
+ struct isea_geo c;
+
+ v = vertex[tri_v1[triangle]];
+ c = icostriangles[triangle];
+
+ /* TODO looks like the adjustment is always either 0 or 180 */
+ /* at least if you pick your vertex carefully */
+ adj = atan2(cos(v.lat) * sin(v.lon - c.lon),
+ cos(c.lat) * sin(v.lat)
+ - sin(c.lat) * cos(v.lat) * cos(v.lon - c.lon));
+ return adj;
+}
+
+/* R tan(g) sin(60) */
+#define TABLE_G 0.6615845383
+
+/* H = 0.25 R tan g = */
+#define TABLE_H 0.1909830056
+
+#define RPRIME 0.91038328153090290025
+
+ISEA_STATIC
+struct isea_pt
+isea_triangle_xy(int triangle)
+{
+ struct isea_pt c;
+ double Rprime = 0.91038328153090290025;
+
+ triangle = (triangle - 1) % 20;
+
+ c.x = TABLE_G * ((triangle % 5) - 2) * 2.0;
+ if (triangle > 9) {
+ c.x += TABLE_G;
+ }
+ switch (triangle / 5) {
+ case 0:
+ c.y = 5.0 * TABLE_H;
+ break;
+ case 1:
+ c.y = TABLE_H;
+ break;
+ case 2:
+ c.y = -TABLE_H;
+ break;
+ case 3:
+ c.y = -5.0 * TABLE_H;
+ break;
+ default:
+ /* should be impossible */
+ exit(EXIT_FAILURE);
+ };
+ c.x *= Rprime;
+ c.y *= Rprime;
+
+ return c;
+}
+
+/* snyder eq 14 */
+static double
+sph_azimuth(double f_lon, double f_lat, double t_lon, double t_lat)
+{
+ double az;
+
+ az = atan2(cos(t_lat) * sin(t_lon - f_lon),
+ cos(f_lat) * sin(t_lat)
+ - sin(f_lat) * cos(t_lat) * cos(t_lon - f_lon)
+ );
+ return az;
+}
+
+/* coord needs to be in radians */
+ISEA_STATIC
+int
+isea_snyder_forward(struct isea_geo * ll, struct isea_pt * out)
+{
+ int i;
+
+ /*
+ * spherical distance from center of polygon face to any of its
+ * vertexes on the globe
+ */
+ double g;
+
+ /*
+ * spherical angle between radius vector to center and adjacent edge
+ * of spherical polygon on the globe
+ */
+ double G;
+
+ /*
+ * plane angle between radius vector to center and adjacent edge of
+ * plane polygon
+ */
+ double theta;
+
+ /* additional variables from snyder */
+ double q, Rprime, H, Ag, Azprime, Az, dprime, f, rho,
+ x, y;
+
+ /* variables used to store intermediate results */
+ double cot_theta, tan_g, az_offset;
+
+ /* how many multiples of 60 degrees we adjust the azimuth */
+ int Az_adjust_multiples;
+
+ struct snyder_constants c;
+
+ /*
+ * TODO by locality of reference, start by trying the same triangle
+ * as last time
+ */
+
+ /* TODO put these constants in as radians to begin with */
+ c = constants[SNYDER_POLY_ICOSAHEDRON];
+ theta = c.theta * DEG2RAD;
+ g = c.g * DEG2RAD;
+ G = c.G * DEG2RAD;
+
+ for (i = 1; i <= 20; i++) {
+ double z;
+ struct isea_geo center;
+
+ center = icostriangles[i];
+
+ /* step 1 */
+#if 0
+ z = sph_distance(center.lon, center.lat, ll->lon, ll->lat);
+#else
+ z = acos(sin(center.lat) * sin(ll->lat)
+ + cos(center.lat) * cos(ll->lat) * cos(ll->lon - center.lon));
+#endif
+
+ /* not on this triangle */
+ if (z > g + 0.000005) { /* TODO DBL_EPSILON */
+ continue;
+ }
+ Az = sph_azimuth(ll->lon, ll->lat, center.lon, center.lat);
+
+ Az = atan2(cos(ll->lat) * sin(ll->lon - center.lon),
+ cos(center.lat) * sin(ll->lat)
+ - sin(center.lat) * cos(ll->lat) * cos(ll->lon - center.lon)
+ );
+
+ /* step 2 */
+
+ /* This calculates "some" vertex coordinate */
+ az_offset = az_adjustment(i);
+
+ Az -= az_offset;
+
+ /* TODO I don't know why we do this. It's not in snyder */
+ /* maybe because we should have picked a better vertex */
+ if (Az < 0.0) {
+ Az += 2.0 * M_PI;
+ }
+ /*
+ * adjust Az for the point to fall within the range of 0 to
+ * 2(90 - theta) or 60 degrees for the hexagon, by
+ * and therefore 120 degrees for the triangle
+ * of the icosahedron
+ * subtracting or adding multiples of 60 degrees to Az and
+ * recording the amount of adjustment
+ */
+
+ Az_adjust_multiples = 0;
+ while (Az < 0.0) {
+ Az += DEG120;
+ Az_adjust_multiples--;
+ }
+ while (Az > DEG120 + DBL_EPSILON) {
+ Az -= DEG120;
+ Az_adjust_multiples++;
+ }
+
+ /* step 3 */
+ cot_theta = 1.0 / tan(theta);
+ tan_g = tan(g); /* TODO this is a constant */
+
+ /* Calculate q from eq 9. */
+ /* TODO cot_theta is cot(30) */
+ q = atan2(tan_g, cos(Az) + sin(Az) * cot_theta);
+
+ /* not in this triangle */
+ if (z > q + 0.000005) {
+ continue;
+ }
+ /* step 4 */
+
+ /* Apply equations 5-8 and 10-12 in order */
+
+ /* eq 5 */
+ /* Rprime = 0.9449322893 * R; */
+ /* R' in the paper is for the truncated */
+ Rprime = 0.91038328153090290025;
+
+ /* eq 6 */
+ H = acos(sin(Az) * sin(G) * cos(g) - cos(Az) * cos(G));
+
+ /* eq 7 */
+ /* Ag = (Az + G + H - DEG180) * M_PI * R * R / DEG180; */
+ Ag = Az + G + H - DEG180;
+
+ /* eq 8 */
+ Azprime = atan2(2.0 * Ag, Rprime * Rprime * tan_g * tan_g - 2.0 * Ag * cot_theta);
+
+ /* eq 10 */
+ /* cot(theta) = 1.73205080756887729355 */
+ dprime = Rprime * tan_g / (cos(Azprime) + sin(Azprime) * cot_theta);
+
+ /* eq 11 */
+ f = dprime / (2.0 * Rprime * sin(q / 2.0));
+
+ /* eq 12 */
+ rho = 2.0 * Rprime * f * sin(z / 2.0);
+
+ /*
+ * add back the same 60 degree multiple adjustment from step
+ * 2 to Azprime
+ */
+
+ Azprime += DEG120 * Az_adjust_multiples;
+
+ /* calculate rectangular coordinates */
+
+ x = rho * sin(Azprime);
+ y = rho * cos(Azprime);
+
+ /*
+ * TODO
+ * translate coordinates to the origin for the particular
+ * hexagon on the flattened polyhedral map plot
+ */
+
+ out->x = x;
+ out->y = y;
+
+ return i;
+ }
+
+ /*
+ * should be impossible, this implies that the coordinate is not on
+ * any triangle
+ */
+
+ fprintf(stderr, "impossible transform: %f %f is not on any triangle\n",
+ ll->lon * RAD2DEG, ll->lat * RAD2DEG);
+
+ exit(EXIT_FAILURE);
+
+ /* not reached */
+ return 0; /* supresses a warning */
+}
+
+/*
+ * return the new coordinates of any point in orginal coordinate system.
+ * Define a point (newNPold) in orginal coordinate system as the North Pole in
+ * new coordinate system, and the great circle connect the original and new
+ * North Pole as the lon0 longitude in new coordinate system, given any point
+ * in orginal coordinate system, this function return the new coordinates.
+ */
+
+#define PRECISION 0.0000000000005
+
+/*
+ * TODO rename to isea_coordtrans and take an struct pointer to the result
+ * struct TODO remove function entirely. we use the synder one. the only
+ * difference is a 180 degree longitude rotation
+ */
+ISEA_STATIC
+struct isea_geo
+coordtrans(struct isea_geo * newNPold, struct isea_geo * ptold, double lon0)
+{
+ double cosptnewlat, cosptnewlon;
+ struct isea_geo ptnew;
+
+ cosptnewlat = sin(newNPold->lat) * sin(ptold->lat) +
+ cos(newNPold->lat) * cos(ptold->lat) * cos(newNPold->lon - ptold->lon);
+
+ if (cosptnewlat > 1.)
+ cosptnewlat = 1.0;
+ if (cosptnewlat < -1.)
+ cosptnewlat = -1.0;
+
+ ptnew.lat = acos(cosptnewlat);
+
+ if (fabs(ptnew.lat - 0.0) < PRECISION * 100000)
+ ptnew.lon = 0.0;
+ else if (fabs(ptnew.lat - M_PI) < PRECISION * 100000)
+ ptnew.lon = 0.0;
+ else {
+ cosptnewlon = (sin(ptold->lat) * cos(newNPold->lat) - cos(ptold->lat) *
+ sin(newNPold->lat) * cos(newNPold->lon - ptold->lon)) / sin(ptnew.lat);
+ if (cosptnewlon > 1.)
+ cosptnewlon = 1.0;
+ if (cosptnewlon < -1.)
+ cosptnewlon = -1.0;
+
+ ptnew.lon = acos(cosptnewlon);
+
+ if ((ptold->lon - newNPold->lon) >= 0 && (ptold->lon - newNPold->lon) < 180)
+ ptnew.lon = -ptnew.lon + lon0;
+
+ else
+ ptnew.lon = ptnew.lon + lon0;
+
+ if (ptnew.lon > M_PI)
+ ptnew.lon -= 2 * M_PI;
+ if (ptnew.lon < -M_PI)
+ ptnew.lon += 2 * M_PI;
+ }
+ ptnew.lat = M_PI / 2 - ptnew.lat;
+ return ptnew;
+}
+
+/* formula from Snyder, Map Projections: A working manual, p31 */
+/*
+ * old north pole at np in new coordinates
+ * could be simplified a bit with fewer intermediates
+ *
+ * TODO take a result pointer
+ */
+ISEA_STATIC
+struct isea_geo
+snyder_ctran(struct isea_geo * np, struct isea_geo * pt)
+{
+ struct isea_geo npt;
+ double alpha, phi, lambda, lambda0, beta, lambdap, phip;
+ double sin_phip;
+ double lp_b; /* lambda prime minus beta */
+ double cos_p, sin_a;
+
+ phi = pt->lat;
+ lambda = pt->lon;
+ alpha = np->lat;
+ beta = np->lon;
+ lambda0 = beta;
+
+ cos_p = cos(phi);
+ sin_a = sin(alpha);
+
+ /* mpawm 5-7 */
+ sin_phip = sin_a * sin(phi) - cos(alpha) * cos_p * cos(lambda - lambda0);
+
+ /* mpawm 5-8b */
+
+ /* use the two argument form so we end up in the right quadrant */
+ lp_b = atan2(cos_p * sin(lambda - lambda0),
+ (sin_a * cos_p * cos(lambda - lambda0) + cos(alpha) * sin(phi)));
+
+ lambdap = lp_b + beta;
+
+ /* normalize longitude */
+ /* TODO can we just do a modulus ? */
+ lambdap = fmod(lambdap, 2 * M_PI);
+ while (lambdap > M_PI)
+ lambdap -= 2 * M_PI;
+ while (lambdap < -M_PI)
+ lambdap += 2 * M_PI;
+
+ phip = asin(sin_phip);
+
+ npt.lat = phip;
+ npt.lon = lambdap;
+
+ return npt;
+}
+
+ISEA_STATIC
+struct isea_geo
+isea_ctran(struct isea_geo * np, struct isea_geo * pt, double lon0)
+{
+ struct isea_geo npt;
+
+ np->lon += M_PI;
+ npt = snyder_ctran(np, pt);
+ np->lon -= M_PI;
+
+ npt.lon -= (M_PI - lon0 + np->lon);
+
+ /*
+ * snyder is down tri 3, isea is along side of tri1 from vertex 0 to
+ * vertex 1 these are 180 degrees apart
+ */
+ npt.lon += M_PI;
+ /* normalize longitude */
+ npt.lon = fmod(npt.lon, 2 * M_PI);
+ while (npt.lon > M_PI)
+ npt.lon -= 2 * M_PI;
+ while (npt.lon < -M_PI)
+ npt.lon += 2 * M_PI;
+
+ return npt;
+}
+
+/* in radians */
+#define ISEA_STD_LAT 1.01722196792335072101
+#define ISEA_STD_LON .19634954084936207740
+
+ISEA_STATIC
+struct isea_dgg isea_standard_dgg = {
+ 20, 58.28252559, 11.25, 0.0,
+ 6, 4, 0, 6, 1.0
+};
+
+/* fuller's at 5.2454 west, 2.3009 N, adjacent at 7.46658 deg */
+
+ISEA_STATIC
+int
+isea_grid_init(struct isea_dgg * g)
+{
+ if (!g)
+ return 0;
+
+ g->polyhedron = 20;
+ g->o_lat = ISEA_STD_LAT;
+ g->o_lon = ISEA_STD_LON;
+ g->o_az = 0.0;
+ g->aperture = 4;
+ g->resolution = 6;
+ g->radius = 1.0;
+ g->topology = 6;
+
+ return 1;
+}
+
+ISEA_STATIC
+int
+isea_orient_isea(struct isea_dgg * g)
+{
+ if (!g)
+ return 0;
+ g->o_lat = ISEA_STD_LAT;
+ g->o_lon = ISEA_STD_LON;
+ g->o_az = 0.0;
+ return 1;
+}
+
+ISEA_STATIC
+int
+isea_orient_pole(struct isea_dgg * g)
+{
+ if (!g)
+ return 0;
+ g->o_lat = M_PI / 2.0;
+ g->o_lon = 0.0;
+ g->o_az = 0;
+ return 1;
+}
+
+ISEA_STATIC
+int
+isea_orient_dymax(struct isea_dgg * g)
+{
+ if (!g)
+ return 0;
+ g->o_lat = 2.300882 * M_PI / 180.0;
+ g->o_lon = -5.24539 * M_PI / 180.0;
+ g->o_az = 7.46658 * M_PI / 180.0;
+ return 1;
+}
+
+ISEA_STATIC
+int
+isea_transform(struct isea_dgg * g, struct isea_geo * in,
+ struct isea_pt * out)
+{
+ struct isea_geo i, pole;
+ int tri;
+
+ pole.lat = g->o_lat;
+ pole.lon = g->o_lon;
+
+ i = isea_ctran(&pole, in, g->o_az);
+
+ tri = isea_snyder_forward(&i, out);
+ out->x *= g->radius;
+ out->y *= g->radius;
+ g->triangle = tri;
+
+ return tri;
+}
+
+#define DOWNTRI(tri) (((tri - 1) / 5) % 2 == 1)
+
+/* get which triangle a point is in, assumes radius = 1.0 */
+ISEA_STATIC
+int isea_xy_triangle(struct isea_pt *p, double radius) {
+ double dist, test;
+ int tri, closest;
+ struct isea_pt tc, pt;
+ pt.x = p->x / radius;
+ pt.y = p->y / radius;
+
+ tc = isea_triangle_xy(1);
+ dist = (pt.x - tc.x) * (pt.x - tc.x) + (pt.y - tc.y) * (pt.y - tc.y);
+ closest = 1;
+
+ /* TODO just calculate it directly, no need to actually do all this
+ * work
+ */
+ for (tri = 2; tri <= 20; tri++) {
+ tc = isea_triangle_xy(tri);
+ test = (pt.x - tc.x) * (pt.x - tc.x) + (pt.y - tc.y) * (pt.y - tc.y);
+ if (test < dist) {
+ dist = test;
+ closest = tri;
+ }
+ }
+ return closest;
+}
+
+ISEA_STATIC
+int isea_projtri(struct isea_pt *xy, double radius) {
+ int tri;
+
+ tri = isea_xy_triangle(xy, radius);
+
+ xy->x /= radius;
+ xy->y /= radius;
+ xy->x += 0.5;
+ xy->y += 2.0 * .14433756729740644112;
+
+ return tri;
+}
+
+ISEA_STATIC
+void
+isea_rotate(struct isea_pt * pt, double degrees)
+{
+ double rad;
+
+ double x, y;
+
+ rad = -degrees * M_PI / 180.0;
+ while (rad >= 2.0 * M_PI) rad -= 2.0 * M_PI;
+ while (rad <= -2.0 * M_PI) rad += 2.0 * M_PI;
+
+ x = pt->x * cos(rad) + pt->y * sin(rad);
+ y = -pt->x * sin(rad) + pt->y * cos(rad);
+
+ pt->x = x;
+ pt->y = y;
+}
+
+ISEA_STATIC
+int isea_tri_plane(int tri, struct isea_pt *pt, double radius) {
+ struct isea_pt tc; /* center of triangle */
+
+ if (DOWNTRI(tri)) {
+ isea_rotate(pt, 180.0);
+ }
+ tc = isea_triangle_xy(tri);
+ tc.x *= radius;
+ tc.y *= radius;
+ pt->x += tc.x;
+ pt->y += tc.y;
+
+ return tri;
+}
+
+ISEA_STATIC
+int isea_projection(struct isea_dgg *g, double xin, double yin, double *xout,
+ double *yout) {
+ return 1;
+}
+
+/* convert projected triangle coords to quad xy coords, return quad number */
+ISEA_STATIC
+int
+isea_ptdd(int tri, struct isea_pt *pt) {
+ int downtri, quad;
+
+ downtri = (((tri - 1) / 5) % 2 == 1);
+ quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;
+
+ isea_rotate(pt, downtri ? 240.0 : 60.0);
+ if (downtri) {
+ pt->x += 0.5;
+ /* pt->y += cos(30.0 * M_PI / 180.0); */
+ pt->y += .86602540378443864672;
+ }
+ return quad;
+}
+
+ISEA_STATIC
+int
+isea_dddi_ap3odd(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di)
+{
+ struct isea_pt v;
+ double hexwidth;
+ double sidelength; /* in hexes */
+ int d, i;
+ int maxcoord;
+ struct hex h;
+
+ /* This is the number of hexes from apex to base of a triangle */
+ sidelength = (pow(2.0, g->resolution) + 1.0) / 2.0;
+
+ /* apex to base is cos(30deg) */
+ hexwidth = cos(M_PI / 6.0) / sidelength;
+
+ /* TODO I think sidelength is always x.5, so
+ * (int)sidelength * 2 + 1 might be just as good
+ */
+ maxcoord = (int) (sidelength * 2.0 + 0.5);
+
+ v = *pt;
+ hexbin2(0, hexwidth, v.x, v.y, &h.x, &h.y);
+ h.iso = 0;
+ hex_iso(&h);
+
+ d = h.x - h.z;
+ i = h.x + h.y + h.y;
+
+ /*
+ * you want to test for max coords for the next quad in the same
+ * "row" first to get the case where both are max
+ */
+ if (quad <= 5) {
+ if (d == 0 && i == maxcoord) {
+ /* north pole */
+ quad = 0;
+ d = 0;
+ i = 0;
+ } else if (i == maxcoord) {
+ /* upper right in next quad */
+ quad += 1;
+ if (quad == 6)
+ quad = 1;
+ i = maxcoord - d;
+ d = 0;
+ } else if (d == maxcoord) {
+ /* lower right in quad to lower right */
+ quad += 5;
+ d = 0;
+ }
+ } else if (quad >= 6) {
+ if (i == 0 && d == maxcoord) {
+ /* south pole */
+ quad = 11;
+ d = 0;
+ i = 0;
+ } else if (d == maxcoord) {
+ /* lower right in next quad */
+ quad += 1;
+ if (quad == 11)
+ quad = 6;
+ d = maxcoord - i;
+ i = 0;
+ } else if (i == maxcoord) {
+ /* upper right in quad to upper right */
+ quad = (quad - 4) % 5;
+ i = 0;
+ }
+ }
+
+ di->x = d;
+ di->y = i;
+
+ g->quad = quad;
+ return quad;
+}
+
+ISEA_STATIC
+int
+isea_dddi(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di) {
+ struct isea_pt v;
+ double hexwidth;
+ int sidelength; /* in hexes */
+ struct hex h;
+
+ if (g->aperture == 3 && g->resolution % 2 != 0) {
+ return isea_dddi_ap3odd(g, quad, pt, di);
+ }
+ /* todo might want to do this as an iterated loop */
+ if (g->aperture >0) {
+ sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
+ } else {
+ sidelength = g->resolution;
+ }
+
+ hexwidth = 1.0 / sidelength;
+
+ v = *pt;
+ isea_rotate(&v, -30.0);
+ hexbin2(0, hexwidth, v.x, v.y, &h.x, &h.y);
+ h.iso = 0;
+ hex_iso(&h);
+
+ /* we may actually be on another quad */
+ if (quad <= 5) {
+ if (h.x == 0 && h.z == -sidelength) {
+ /* north pole */
+ quad = 0;
+ h.z = 0;
+ h.y = 0;
+ h.x = 0;
+ } else if (h.z == -sidelength) {
+ quad = quad + 1;
+ if (quad == 6)
+ quad = 1;
+ h.y = sidelength - h.x;
+ h.z = h.x - sidelength;
+ h.x = 0;
+ } else if (h.x == sidelength) {
+ quad += 5;
+ h.y = -h.z;
+ h.x = 0;
+ }
+ } else if (quad >= 6) {
+ if (h.z == 0 && h.x == sidelength) {
+ /* south pole */
+ quad = 11;
+ h.x = 0;
+ h.y = 0;
+ h.z = 0;
+ } else if (h.x == sidelength) {
+ quad = quad + 1;
+ if (quad == 11)
+ quad = 6;
+ h.x = h.y + sidelength;
+ h.y = 0;
+ h.z = -h.x;
+ } else if (h.y == -sidelength) {
+ quad -= 4;
+ h.y = 0;
+ h.z = -h.x;
+ }
+ }
+ di->x = h.x;
+ di->y = -h.z;
+
+ g->quad = quad;
+ return quad;
+}
+
+ISEA_STATIC
+int isea_ptdi(struct isea_dgg *g, int tri, struct isea_pt *pt,
+ struct isea_pt *di) {
+ struct isea_pt v;
+ int quad;
+
+ v = *pt;
+ quad = isea_ptdd(tri, &v);
+ quad = isea_dddi(g, quad, &v, di);
+ return quad;
+}
+
+/* q2di to seqnum */
+ISEA_STATIC
+int isea_disn(struct isea_dgg *g, int quad, struct isea_pt *di) {
+ int sidelength;
+ int sn, height;
+ int hexes;
+
+ if (quad == 0) {
+ g->serial = 1;
+ return g->serial;
+ }
+ /* hexes in a quad */
+ hexes = (int) (pow(g->aperture, g->resolution) + 0.5);
+ if (quad == 11) {
+ g->serial = 1 + 10 * hexes + 1;
+ return g->serial;
+ }
+ if (g->aperture == 3 && g->resolution % 2 == 1) {
+ height = (int) (pow(g->aperture, (g->resolution - 1) / 2.0));
+ sn = ((int) di->x) * height;
+ sn += ((int) di->y) / height;
+ sn += (quad - 1) * hexes;
+ sn += 2;
+ } else {
+ sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
+ sn = (quad - 1) * hexes + sidelength * di->x + di->y + 2;
+ }
+
+ g->serial = sn;
+ return sn;
+}
+
+/* TODO just encode the quad in the d or i coordinate
+ * quad is 0-11, which can be four bits.
+ * d' = d << 4 + q, d = d' >> 4, q = d' & 0xf
+ */
+/* convert a q2di to global hex coord */
+ISEA_STATIC
+int isea_hex(struct isea_dgg *g, int tri,
+ struct isea_pt *pt, struct isea_pt *hex) {
+ struct isea_pt v;
+ int sidelength;
+ int d, i, x, y, quad;
+ double oddsl;
+
+ quad = isea_ptdi(g, tri, pt, &v);
+
+ hex->x = ((int)v.x << 4) + quad;
+ hex->y = v.y;
+
+ return 1;
+
+ d = v.x;
+ i = v.y;
+
+ /* Aperture 3 odd resolutions */
+ if (g->aperture == 3 && g->resolution % 2 != 0) {
+ int offset = (int)(pow(3.0, g->resolution - 1) + 0.5);
+
+ oddsl = (pow(2.0, g->resolution) + 1.0) / 2.0;
+
+ d += offset * ((g->quad-1) % 5);
+ i += offset * ((g->quad-1) % 5);
+
+ if (quad == 0) {
+ d = 0;
+ i = offset;
+ } else if (quad == 11) {
+ d = 2 * offset;
+ i = 0;
+ } else if (quad > 5) {
+ d += offset;
+ }
+
+ x = (2*d - i) /3;
+ y = (2*i - d) /3;
+
+ hex->x = x + offset / 3;
+ hex->y = y + 2 * offset / 3;
+ return 1;
+ }
+
+ /* aperture 3 even resolutions and aperture 4 */
+ sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
+ if (g->quad == 0) {
+ hex->x = 0;
+ hex->y = sidelength;
+ } else if (g->quad == 11) {
+ hex->x = sidelength * 2;
+ hex->y = 0;
+ } else {
+ hex->x = d + sidelength * ((g->quad-1) % 5);
+ if (g->quad > 5) hex->x += sidelength;
+ hex->y = i + sidelength * ((g->quad-1) % 5);
+ }
+
+ return 1;
+}
+
+ISEA_STATIC
+struct isea_pt
+isea_forward(struct isea_dgg *g, struct isea_geo *in)
+{
+ int tri, downtri, quad;
+ struct isea_pt out, coord;
+
+ tri = isea_transform(g, in, &out);
+
+ downtri = (((tri - 1) / 5) % 2 == 1);
+ quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;
+
+ if (g->output == ISEA_PLANE) {
+ isea_tri_plane(tri, &out, g->radius);
+ return out;
+ }
+
+ /* convert to isea standard triangle size */
+ out.x = out.x / g->radius * ISEA_SCALE;
+ out.y = out.y / g->radius * ISEA_SCALE;
+ out.x += 0.5;
+ out.y += 2.0 * .14433756729740644112;
+
+ switch (g->output) {
+ case ISEA_PROJTRI:
+ /* nothing to do, already in projected triangle */
+ break;
+ case ISEA_VERTEX2DD:
+ g->quad = isea_ptdd(tri, &out);
+ break;
+ case ISEA_Q2DD:
+ /* Same as above, we just don't print as much */
+ g->quad = isea_ptdd(tri, &out);
+ break;
+ case ISEA_Q2DI:
+ g->quad = isea_ptdi(g, tri, &out, &coord);
+ return coord;
+ break;
+ case ISEA_SEQNUM:
+ isea_ptdi(g, tri, &out, &coord);
+ /* disn will set g->serial */
+ isea_disn(g, g->quad, &coord);
+ return coord;
+ break;
+ case ISEA_HEX:
+ isea_hex(g, tri, &out, &coord);
+ return coord;
+ break;
+ }
+
+ return out;
+}
--- /dev/null
+/*
+ * options
+ *
+ * orientation
+ * -l longitude
+ * -z azimuth
+ * -p latitude
+ *
+ * -I isea standard orientation
+ * -D dymaxion orientation
+ * -P vertex at pole
+ *
+ * -a aperture
+ * -r resolution
+ * -t h t d hex tri diamond
+ *
+ * addressing
+ * -O q2dd seqnum interleave plane q2dd projtri vertex2dd
+ * -I
+ *
+ * operation
+ * -g generate grid
+ * -x transform
+ * -b bin point values (average)
+ * -e bin point presence
+ * -G spacing resolution generate graticule
+ */
+
+#if 0
+void usage(void) {
+ fprintf(stderr, "%s", "isea [-l <lon>] [-z <azimuth>] [-p <latitude>] [-I] [-D] [-P] [-a <aperture>] [-r <resolution>] [-t <topology>] [-A <addressing>] [-g] [-b] [-e] [-x]\n");
+}
+
+void die(char *s) {
+ fprintf(stderr, "%s\n", s);
+ exit(EXIT_FAILURE);
+}
+
+static int verbose = 0;
+static int precision = 10;
+
+void options(int ac, char **av, struct isea_dgg *g) {
+ int i = 1;
+ char *optarg;
+
+ for (i=1;i<ac;i++) {
+ if (av[i][0] != '-') {
+ break;
+ }
+ switch (av[i][1]) {
+ case 'v':
+ verbose = 1; break;
+ case 'I':
+ isea_orient_isea(g); break;
+ case 'P':
+ isea_orient_pole(g); break;
+ case 'D':
+ isea_orient_dymax(g); break;
+ case 'R':
+ g->radius = strtod(av[++i],NULL); break;
+ case 'l':
+ g->o_lon = strtod(av[++i],NULL); break;
+ case 'p':
+ g->o_lat = strtod(av[++i],NULL); break;
+ case 'z':
+ g->o_az = strtod(av[++i],NULL); break;
+ case 'd':
+ precision = atoi(av[++i]); break;
+ case 'a':
+ g->aperture = atoi(av[++i]); break;
+ case 'r':
+ g->resolution = atoi(av[++i]); break;
+ case 't':
+ switch (av[++i][0]) {
+ case 'h':
+ g->topology = ISEA_HEXAGON;
+ break;
+ case 't':
+ g->topology = ISEA_TRIANGLE;
+ break;
+ case 'd':
+ g->topology = ISEA_DIAMOND;
+ break;
+ default:
+ usage();
+ die("Invalid topology");
+ break;
+ };
+ break;
+ case 'O':
+ if (av[i][2]) {
+ optarg = &av[i][2];
+ } else {
+ optarg = av[++i];
+ }
+ if (!strcmp(optarg, "projtri")) {
+ g->output = ISEA_PROJTRI;
+ }
+ else if (!strcmp(optarg, "vertex2dd")) {
+ g->output = ISEA_VERTEX2DD;
+ }
+ else if (!strcmp(optarg, "q2dd")) {
+ g->output = ISEA_Q2DD;
+ }
+ else if (!strcmp(optarg, "q2di")) {
+ g->output = ISEA_Q2DI;
+ }
+ else if (!strcmp(optarg, "seqnum")) {
+ g->output = ISEA_SEQNUM;
+ }
+ else if (!strcmp(optarg, "plane")) {
+ g->output = ISEA_PLANE;
+ } else {
+ fprintf(stderr, "unknown output format: %s\n",optarg);
+ die("exiting");
+ }
+ break;
+ default:
+ usage();
+ die("Unknown option");
+ }
+ }
+}
+
+#endif
+
+#define PROJ_PARMS__ \
+ struct isea_dgg dgg;
+
+#define PJ_LIB__
+#include <projects.h>
+
+PROJ_HEAD(isea, "Icosahedral Snyder Equal Area") "\n\tSph";
+
+FORWARD(s_forward);
+ struct isea_pt out;
+ struct isea_geo in;
+
+ in.lon = lp.lam;
+ in.lat = lp.phi;
+
+ out = isea_forward(&P->dgg, &in);
+
+ xy.x = out.x;
+ xy.y = out.y;
+
+ return xy;
+}
+FREEUP; if (P) pj_dalloc(P); }
+
+ENTRY0(isea)
+ char *opt;
+
+ P->fwd = s_forward;
+ isea_grid_init(&P->dgg);
+
+ P->dgg.output = ISEA_PLANE;
+/* P->dgg.radius = P->a; /* otherwise defaults to 1 */
+ /* calling library will scale, I think */
+
+ opt = pj_param(P->params, "sorient").s;
+ if (opt) {
+ if (!strcmp(opt, "isea")) {
+ isea_orient_isea(&P->dgg);
+ } else if (!strcmp(opt, "pole")) {
+ isea_orient_pole(&P->dgg);
+ } else {
+ E_ERROR(-34);
+ }
+ }
+
+ if (pj_param(P->params, "tazi").i) {
+ P->dgg.o_az = pj_param(P->params, "razi").f;
+ }
+
+ if (pj_param(P->params, "tlon_0").i) {
+ P->dgg.o_lon = pj_param(P->params, "rlon_0").f;
+ }
+
+ if (pj_param(P->params, "tlat_0").i) {
+ P->dgg.o_lat = pj_param(P->params, "rlat_0").f;
+ }
+
+ if (pj_param(P->params, "taperture").i) {
+ P->dgg.aperture = pj_param(P->params, "iaperture").i;
+ }
+
+ if (pj_param(P->params, "tresolution").i) {
+ P->dgg.resolution = pj_param(P->params, "iresolution").i;
+ }
+
+ opt = pj_param(P->params, "smode").s;
+ if (opt) {
+ if (!strcmp(opt, "plane")) {
+ P->dgg.output = ISEA_PLANE;
+ } else if (!strcmp(opt, "di")) {
+ P->dgg.output = ISEA_Q2DI;
+ }
+ else if (!strcmp(opt, "dd")) {
+ P->dgg.output = ISEA_Q2DD;
+ }
+ else if (!strcmp(opt, "hex")) {
+ P->dgg.output = ISEA_HEX;
+ }
+ else {
+ E_ERROR(-34);
+ }
+ }
+
+ if (pj_param(P->params, "trescale").i) {
+ P->dgg.radius = ISEA_SCALE;
+ }
+
+ if (pj_param(P->params, "tresolution").i) {
+ P->dgg.resolution = pj_param(P->params, "iresolution").i;
+ } else {
+ P->dgg.resolution = 4;
+ }
+
+ if (pj_param(P->params, "taperture").i) {
+ P->dgg.aperture = pj_param(P->params, "iaperture").i;
+ } else {
+ P->dgg.aperture = 3;
+ }
+
+ENDENTRY(P)
+
+#if 0
+int
+main(int ac, char *av[])
+{
+ struct isea_dgg dgg;
+ struct isea_geo in;
+ struct isea_pt out, coord;
+ int tri,calctri;
+ int graticule = 0;
+ int quad = 0;
+
+ char line[1024];
+ char *rest, *nl;
+
+ isea_grid_init(&dgg);
+ dgg.output = ISEA_PLANE;
+ dgg.radius = ISEA_SCALE; /* 0.8301572857837594396028083; */
+
+ options(ac, av, &dgg);
+
+ if (0) {
+ printf("isea lat = %lf, lon = %lf, az = %lf\n",
+ dgg.o_lat * 180.0 / M_PI, dgg.o_lon * 180.0 / M_PI,
+ dgg.o_az * 180.0 / M_PI);
+ }
+
+ while (fgets(line, sizeof line, stdin)) {
+ line[sizeof line-1] = 0;
+ if (line[0] == '#') {
+ continue;
+ }
+ if (!strncmp(line, "ENDPOLY", 7)) {
+ printf("%s", line);
+ continue;
+ }
+ in.lon = strtod(line, &rest);
+ in.lat = strtod(rest, &rest);
+
+ nl = strchr(rest, '\n');
+ if (nl) *nl = 0;
+
+ in.lon *= M_PI / 180.0;
+ in.lat *= M_PI / 180.0;
+
+ out = isea_forward(&dgg, &in);
+ tri = dgg.triangle;
+ quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;
+ coord = out;
+
+ if (dgg.output == ISEA_PLANE) {
+ calctri = isea_xy_triangle(&coord,
+ dgg.output == ISEA_PLANE ? dgg.radius : 1.0);
+ if (tri != calctri) {
+ fprintf(stderr,
+ "got %d (%f, %f), expected %d\n",
+ calctri, out.x, out.y, tri);
+ }
+ }
+
+ in.lon *= 180.0 / M_PI;
+ in.lat *= 180.0 / M_PI;
+
+ switch (dgg.output) {
+ case ISEA_PLANE:
+ printf("%d %.*f %.*f%s\n", tri, precision, out.x, precision, out.y, rest);
+ break;
+ case ISEA_PROJTRI:
+ printf("%d %.*f %.*f%s\n", tri-1, precision, out.x,
+ precision, out.y, rest);
+ break;
+ case ISEA_VERTEX2DD:
+ printf("%d %d %.*f %.*f%s\n", quad, tri-1, precision,
+ out.x, precision, out.y, rest);
+ break;
+ case ISEA_Q2DD:
+ /* Same as above, we just don't print as much */
+ printf("%d %.*f %.*f%s\n", quad, precision,
+ out.x, precision, out.y, rest);
+ break;
+ case ISEA_Q2DI:
+ printf("%d %.0lf %.0lf%s\n", dgg.quad, coord.x, coord.y, rest);
+ break;
+ case ISEA_SEQNUM:
+ quad = isea_ptdi(&dgg, tri, &out, &coord);
+ printf("%d%s\n", isea_disn(&dgg, quad, &coord), rest);
+ break;
+ }
+ fflush(stdout);
+ }
+
+ if (graticule) {
+
+ }
+
+ exit(EXIT_SUCCESS);
+}
+#endif
--- /dev/null
+/*
+ * This code was entirely written by Nathan Wagner
+ * and is in the public domain.
+ */
+
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <float.h>
+
+/*
+ * Proj 4 provides its own entry points into
+ * the code, so none of the library functions
+ * need to be global
+ */
+#define ISEA_STATIC static
--- /dev/null
+/*
+ * Proj 4 integration code follows
+ */
+
+#define PROJ_PARMS__ \
+ struct isea_dgg dgg;
+
+#define PJ_LIB__
+#include <projects.h>
+
+PROJ_HEAD(isea, "Icosahedral Snyder Equal Area") "\n\tSph";
+
+FORWARD(s_forward);
+ struct isea_pt out;
+ struct isea_geo in;
+
+ in.lon = lp.lam;
+ in.lat = lp.phi;
+
+ out = isea_forward(&P->dgg, &in);
+
+ xy.x = out.x;
+ xy.y = out.y;
+
+ return xy;
+}
+FREEUP; if (P) pj_dalloc(P); }
+
+ENTRY0(isea)
+ char *opt;
+
+ P->fwd = s_forward;
+ isea_grid_init(&P->dgg);
+
+ P->dgg.output = ISEA_PLANE;
+/* P->dgg.radius = P->a; /* otherwise defaults to 1 */
+ /* calling library will scale, I think */
+
+ opt = pj_param(P->ctx,P->params, "sorient").s;
+ if (opt) {
+ if (!strcmp(opt, "isea")) {
+ isea_orient_isea(&P->dgg);
+ } else if (!strcmp(opt, "pole")) {
+ isea_orient_pole(&P->dgg);
+ } else {
+ E_ERROR(-34);
+ }
+ }
+
+ if (pj_param(P->ctx,P->params, "tazi").i) {
+ P->dgg.o_az = pj_param(P->ctx,P->params, "razi").f;
+ }
+
+ if (pj_param(P->ctx,P->params, "tlon_0").i) {
+ P->dgg.o_lon = pj_param(P->ctx,P->params, "rlon_0").f;
+ }
+
+ if (pj_param(P->ctx,P->params, "tlat_0").i) {
+ P->dgg.o_lat = pj_param(P->ctx,P->params, "rlat_0").f;
+ }
+
+ if (pj_param(P->ctx,P->params, "taperture").i) {
+ P->dgg.aperture = pj_param(P->ctx,P->params, "iaperture").i;
+ }
+
+ if (pj_param(P->ctx,P->params, "tresolution").i) {
+ P->dgg.resolution = pj_param(P->ctx,P->params, "iresolution").i;
+ }
+
+ opt = pj_param(P->ctx,P->params, "smode").s;
+ if (opt) {
+ if (!strcmp(opt, "plane")) {
+ P->dgg.output = ISEA_PLANE;
+ } else if (!strcmp(opt, "di")) {
+ P->dgg.output = ISEA_Q2DI;
+ }
+ else if (!strcmp(opt, "dd")) {
+ P->dgg.output = ISEA_Q2DD;
+ }
+ else if (!strcmp(opt, "hex")) {
+ P->dgg.output = ISEA_HEX;
+ }
+ else {
+ /* TODO verify error code. Possibly eliminate magic */
+ E_ERROR(-34);
+ }
+ }
+
+ if (pj_param(P->ctx,P->params, "trescale").i) {
+ P->dgg.radius = ISEA_SCALE;
+ }
+
+ if (pj_param(P->ctx,P->params, "tresolution").i) {
+ P->dgg.resolution = pj_param(P->ctx,P->params, "iresolution").i;
+ } else {
+ P->dgg.resolution = 4;
+ }
+
+ if (pj_param(P->ctx,P->params, "taperture").i) {
+ P->dgg.aperture = pj_param(P->ctx,P->params, "iaperture").i;
+ } else {
+ P->dgg.aperture = 3;
+ }
+
+ENDENTRY(P)
--- /dev/null
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <float.h>
+
+#define ISEA_STATIC static
+#include "hfuncs.c"
+#include "pfuncs.c"
+
+
--- /dev/null
+diff --git a/src/Makefile.am b/src/Makefile.am
+index 676ae91..a8ccaca 100644
+--- a/src/Makefile.am
++++ b/src/Makefile.am
+@@ -25,7 +25,7 @@ libproj_la_LDFLAGS = -no-undefined -version-info 6:6:6
+
+ libproj_la_SOURCES = \
+ projects.h pj_list.h \
+- PJ_aeqd.c PJ_gnom.c PJ_laea.c PJ_mod_ster.c \
++ PJ_aeqd.c PJ_gnom.c PJ_laea.c PJ_isea.c PJ_mod_ster.c \
+ PJ_nsper.c PJ_nzmg.c PJ_ortho.c PJ_stere.c PJ_sterea.c \
+ PJ_aea.c PJ_bipc.c PJ_bonne.c PJ_eqdc.c \
+ PJ_imw_p.c PJ_krovak.c PJ_lcc.c PJ_poly.c \
+diff --git a/src/Makefile.in b/src/Makefile.in
+index 0f4d323..8bf2686 100644
+--- a/src/Makefile.in
++++ b/src/Makefile.in
+@@ -57,7 +57,7 @@ am__installdirs = "$(DESTDIR)$(libdir)" "$(DESTDIR)$(bindir)" \
+ libLTLIBRARIES_INSTALL = $(INSTALL)
+ LTLIBRARIES = $(lib_LTLIBRARIES)
+ libproj_la_LIBADD =
+-am_libproj_la_OBJECTS = PJ_aeqd.lo PJ_gnom.lo PJ_laea.lo \
++am_libproj_la_OBJECTS = PJ_aeqd.lo PJ_gnom.lo PJ_laea.lo PJ_isea.lo \
+ PJ_mod_ster.lo PJ_nsper.lo PJ_nzmg.lo PJ_ortho.lo PJ_stere.lo \
+ PJ_sterea.lo PJ_aea.lo PJ_bipc.lo PJ_bonne.lo PJ_eqdc.lo \
+ PJ_imw_p.lo PJ_krovak.lo PJ_lcc.lo PJ_poly.lo PJ_rpoly.lo \
+@@ -259,7 +259,7 @@ lib_LTLIBRARIES = libproj.la
+ libproj_la_LDFLAGS = -no-undefined -version-info 6:6:6
+ libproj_la_SOURCES = \
+ projects.h pj_list.h \
+- PJ_aeqd.c PJ_gnom.c PJ_laea.c PJ_mod_ster.c \
++ PJ_aeqd.c PJ_gnom.c PJ_laea.c PJ_isea.c PJ_mod_ster.c \
+ PJ_nsper.c PJ_nzmg.c PJ_ortho.c PJ_stere.c PJ_sterea.c \
+ PJ_aea.c PJ_bipc.c PJ_bonne.c PJ_eqdc.c \
+ PJ_imw_p.c PJ_krovak.c PJ_lcc.c PJ_poly.c \
+@@ -464,6 +464,7 @@ distclean-compile:
+ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_krovak.Plo@am__quote@
+ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_labrd.Plo@am__quote@
+ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_laea.Plo@am__quote@
++@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_isea.Plo@am__quote@
+ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_lagrng.Plo@am__quote@
+ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_larr.Plo@am__quote@
+ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/PJ_lask.Plo@am__quote@
+diff --git a/src/makefile.vc b/src/makefile.vc
+index 54394d4..057a194 100644
+--- a/src/makefile.vc
++++ b/src/makefile.vc
+@@ -25,7 +25,8 @@ misc = \
+ PJ_chamb.obj PJ_hammer.obj PJ_lagrng.obj PJ_larr.obj \
+ PJ_lask.obj PJ_nocol.obj PJ_ob_tran.obj PJ_oea.obj \
+ PJ_tpeqd.obj PJ_vandg.obj PJ_vandg2.obj PJ_vandg4.obj \
+- PJ_wag7.obj pj_latlong.obj PJ_krovak.obj pj_geocent.obj
++ PJ_wag7.obj pj_latlong.obj PJ_krovak.obj pj_geocent.obj \
++ PJ_isea.obj
+
+ pseudo = \
+ PJ_boggs.obj PJ_collg.obj PJ_crast.obj PJ_denoy.obj \
+diff --git a/src/pj_list.h b/src/pj_list.h
+index b965cc7..e421c3a 100644
+--- a/src/pj_list.h
++++ b/src/pj_list.h
+@@ -47,6 +47,7 @@ PROJ_HEAD(gs50, "Mod. Stererographics of 50 U.S.")
+ PROJ_HEAD(hammer, "Hammer & Eckert-Greifendorff")
+ PROJ_HEAD(hatano, "Hatano Asymmetrical Equal Area")
+ PROJ_HEAD(imw_p, "Internation Map of the World Polyconic")
++PROJ_HEAD(isea, "Icosahedral Snyder Equal Area")
+ PROJ_HEAD(kav5, "Kavraisky V")
+ PROJ_HEAD(kav7, "Kavraisky VII")
+ PROJ_HEAD(krovak, "Krovak")
--- /dev/null
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+
+#include "isea.h"
+#include "hex.h"
+
+void die(char *s) {
+ fprintf(stderr, "%s\n", s);
+ exit(EXIT_FAILURE);
+}
+
+static int verbose = 0;
+static int precision = 10;
+
+int
+main(int ac, char *av[])
+{
+ struct isea_dgg dgg;
+ struct isea_geo in;
+ struct isea_pt out, coord;
+ int tri,calctri;
+ int graticule = 0;
+ int quad = 0;
+
+ char line[1024];
+ char *rest, *nl;
+
+ isea_grid_init(&dgg);
+ dgg.output = ISEA_PLANE;
+ dgg.radius = 0.8301572857837594396028083;
+
+ options(ac, av, &dgg);
+
+ if (0) {
+ printf("isea lat = %lf, lon = %lf, az = %lf\n",
+ dgg.o_lat * 180.0 / M_PI, dgg.o_lon * 180.0 / M_PI,
+ dgg.o_az * 180.0 / M_PI);
+ }
+
+ while (fgets(line, sizeof line, stdin)) {
+ if (line[0] == '#') {
+ continue;
+ }
+ if (!strncmp(line, "ENDPOLY", 7)) {
+ printf("%s", line);
+ continue;
+ }
+ in.lon = strtod(line, &rest);
+ in.lat = strtod(rest, &rest);
+
+ nl = strchr(rest, '\n');
+ if (nl) *nl = 0;
+
+ in.lon *= M_PI / 180.0;
+ in.lat *= M_PI / 180.0;
+
+ out = isea_forward(&dgg, &in);
+ tri = dgg.triangle;
+ quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;
+ coord = out;
+
+ if (dgg.output == ISEA_PLANE) {
+ calctri = isea_xy_triangle(&coord,
+ dgg.output == ISEA_PLANE ? dgg.radius : 1.0);
+ if (tri != calctri) {
+ fprintf(stderr,
+ "got %d (%f, %f), expected %d\n",
+ calctri, out.x, out.y, tri);
+ }
+ }
+
+ in.lon *= 180.0 / M_PI;
+ in.lat *= 180.0 / M_PI;
+
+ switch (dgg.output) {
+ case ISEA_PLANE:
+ printf("%d %.*f %.*f%s\n", tri, precision, out.x, precision, out.y, rest);
+ break;
+ case ISEA_PROJTRI:
+ printf("%d %.*f %.*f%s\n", tri-1, precision, out.x,
+ precision, out.y, rest);
+ break;
+ case ISEA_VERTEX2DD:
+ printf("%d %d %.*f %.*f%s\n", quad, tri-1, precision,
+ out.x, precision, out.y, rest);
+ break;
+ case ISEA_Q2DD:
+ /* Same as above, we just don't print as much */
+ printf("%d %.*f %.*f%s\n", quad, precision,
+ out.x, precision, out.y, rest);
+ break;
+ case ISEA_Q2DI:
+ printf("%d %.0lf %.0lf%s\n", dgg.quad, coord.x, coord.y, rest);
+ break;
+ case ISEA_SEQNUM:
+ quad = isea_ptdi(&dgg, tri, &out, &coord);
+ printf("%d%s\n", isea_disn(&dgg, quad, &coord), rest);
+ break;
+ }
+ fflush(stdout);
+ }
+
+ if (graticule) {
+
+ }
+
+ exit(EXIT_SUCCESS);
+}
--- /dev/null
+5.3700177,2.0396587,49900 Corvallis
+5.3578962,2.0451813,42300 Aloha
+5.3589320,2.0461297,41700 Tigard
+5.3696974,2.0417109,41400 Albany
+5.3593468,2.0468240,35700 LakeOswego
+5.3644930,2.0429646,32600 Keizer
+5.3612037,2.0414984,26800 MacMinnville
+5.3605206,2.0477630,26100 OregonCity
+5.3998953,2.0361897,23300 GrantsPass
+5.3596352,2.0460631,23100 Tualatin
+5.3600231,2.0473897,22500 WestLinn
+5.3590763,2.0476873,20700 Milwaukie
+5.3593678,2.0628116,20600 CitrusPark
+5.3626963,2.0448257,20400 Woodburn
+5.3888368,2.0368818,20300 Roseburg
+
--- /dev/null
+struct points cities[] = {
+{-123.28, 44.57, "49900 Corvallis"},
+{-122.87, 45.49, "42300 Aloha"},
+{-122.77, 45.43, "41700 Tigard"},
+{-123.09, 44.62, "41400 Albany"},
+{-122.70, 45.41, "35700 LakeOswego"},
+{-123.02, 45.00, "32600 Keizer"},
+{-123.19, 45.21, "26800 MacMinnville"},
+{-122.60, 45.34, "26100 OregonCity"},
+{-123.32, 42.44, "23300 GrantsPass"},
+{-122.77, 45.38, "23100 Tualatin"},
+{-122.64, 45.37, "22500 WestLinn"},
+{-122.62, 45.44, "20700 Milwaukie"},
+{-121.17, 45.60, "20600 CitrusPark"},
+{-122.86, 45.15, "20400 Woodburn"},
+{-123.36, 43.22, "20300 Roseburg"},
+{0.0, 0.0, NULL}
+};
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include "tap.h"
+#include "isea.h"
+
+#include "testing.h"
+
+#include "hexes.pt"
+
+struct dgg {
+ int x, y;
+ char *data;
+} dgg[] = {
+{0, 0, "vertex 1"},
+{1, 0, "vertex 2 r4 1 0 0"},
+{1, 1, "step 1 r4 1 0 1"},
+{1, 2, "step 2 r4 1 0 2"},
+{1, 3, "step 3 r4 1 0 3"},
+{49, 6, "triangle center 1 r4 1 3 6"},
+{2, 0, "vertex 3"},
+{2, 1, "step 3.1"},
+{3, 0, "vertex 4"},
+{4, 0, "vertex 5"},
+{5, 0, "vertex 6"},
+{6, 0, "vertex 7"},
+{7, 0, "vertex 8"},
+{8, 0, "vertex 9"},
+{9, 0, "vertex 10"},
+{10, 0, "vertex 11"},
+{11, 0, "vertex 12"},
+{0, 0, NULL}
+};
+
+struct dggap3 {
+ int x, y;
+ char *data;
+} dggap3[] = {
+{0, 0, "vertex 1"},
+{1, 0, "vertex 2 r4 1 0 0"},
+{1, 1, "step 1 r4 1 0 1"},
+{1, 2, "step 2 r4 1 0 2"},
+{1, 3, "step 3 r4 1 0 3"},
+{49, 6, "triangle center 1 r4 1 3 6"},
+{6, 9, "vertex 3"},
+{6, 9, "step 3.1"},
+{9, 12, "vertex 4"},
+{12, 15, "vertex 5"},
+{15, 18, "vertex 6"},
+{9, 3, "vertex 7"},
+{12, 6, "vertex 8"},
+{15, 9, "vertex 9"},
+{18, 12, "vertex 10"},
+{21, 15, "vertex 11"},
+{15, 0, "vertex 12"},
+{0, 0, NULL}
+};
+
+int main(void) {
+ int i, j;
+
+ struct isea_dgg g;
+ struct isea_pt xy;
+ struct isea_geo ll;
+
+ isea_grid_init(&g);
+ isea_orient_pole(&g);
+ g.output = ISEA_HEX;
+ g.radius = ISEA_SCALE;
+ g.aperture = 3;
+ g.resolution = 4;
+
+ i=0;
+ while (input[i].data) i++;
+
+ plan_tests(i * 2 * 1);
+
+ for (j=0; j<i; j++) {
+ ll.lon = input[j].lon * M_PI / 180.0;
+ ll.lat = input[j].lat * M_PI / 180.0;
+ xy = isea_forward(&g, &ll);
+ ok(xy.x == dgg[j].x, "hex X %d == %d (%f, %f %s)",
+ (int)xy.x, dgg[j].x,
+ input[j].lon, input[j].lat, input[j].data
+ );
+ ok(xy.y == dgg[j].y, "hex Y %d == %d", (int)xy.y, dgg[j].y);
+ }
+
+ if (0) {
+ g.resolution = 3;
+ for (j=0; j<i; j++) {
+ if (j == 7 || j == 2 || j == 3) continue;
+ ll.lon = input[j].lon * M_PI / 180.0;
+ ll.lat = input[j].lat * M_PI / 180.0;
+ xy = isea_forward(&g, &ll);
+ skip_start(j == 7 || j == 2 || j == 3, 2, "skipping ap3 input %d as it falls on a hex side and is numerically unstable", j);
+ ok(xy.x == dggap3[j].x, "hex ap3 X %d == %d (%f, %f %s)",
+ (int)xy.x, dggap3[j].x,
+ input[j].lon, input[j].lat, input[j].data
+ );
+ ok(xy.y == dggap3[j].y, "hex ap3 Y %d == %d", (int)xy.y,
+ dggap3[j].y);
+ skip_end;
+ }
+ }
+
+ return exit_status();
+
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "tap.h"
+#include "hex.h"
+
+void iso(struct hex *h) {
+ hex_iso(h);
+ printf("iso(%d,%d,%d)", h->x, h->y, h->z);
+}
+void xy(struct hex *h) {
+ hex_xy(h);
+ printf("xy(%d,%d)", h->x, h->y);
+}
+
+int main(void) {
+ int x, y;
+ double p, q;
+
+ plan_tests(5);
+
+ hexbin2(0, 1.0, 0.0, 0.0, &x, &y);
+ ok(x == 0 && y == 0, "0.0 0.0 -> 0 0");
+
+ hexbin2(0, 0.1111111, 0.288675, 0.500000, &x, &y);
+ ok(x == 3, "center bin X %d == 3", x);
+ ok(y == -5, "center bin Y %d == -5", y);
+
+ hexbin2(0, 0.1111111, 0.866025, 0.500000, &x, &y);
+ ok(x == 9, "vertex bin X %d == 9", x);
+ ok(y == -5, "vertex bin Y %d == -5", y);
+
+ return exit_status();
+
+ for (p = -4.0; p < 4.0; p+= 0.4) {
+ for (q = -4.0; q < 4.0; q += 0.4) {
+ hexbin2(0, 1.0, p, q, &x, &y);
+ printf("(%f, %f) = %d %d\n", p, q, x, y);
+ }
+ }
+
+}
+
--- /dev/null
+struct hexes {
+ double lon, lat;
+ char *data;
+} input[] = {
+{0.0, 90.0, "vertex 1 r4 0 0 0"},
+{180.0, 26.565051177, "vertex 2 r4 1 0 0"},
+{180.0, 33.61337882400000000000, "step 1 r4 1 0 1"},
+{180.0, 40.66170647100000000000, "step 2 r4 1 0 2"},
+{180.0, 47.71003411800000000000, "step 3 r4 1 0 3"},
+{-144.0, 52.62263186, "triangle center 1 r4 1 3 6"},
+{-108.0, 26.565051177, "vertex 3 r4 2 0 0"},
+{-108.0, 33.61337882400000000000, "step 3.1 @2 0 1"},
+{-36.0, 26.565051177, "vertex 4 r4 3 0 0"},
+{36.0, 26.565051177, "vertex 5 r4 4 0 0"},
+{108.0, 26.565051177, "vertex 6 r4 5 0 0"},
+{-144.0, -26.565051177, "vertex 7 r4 6 0 0"},
+{-72.0, -26.565051177, "vertex 8 r4 7 0 0"},
+{0.0, -26.565051177, "vertex 9 r4 8 0 0"},
+{72.0, -26.565051177, "vertex 10 r4 9 0 0"},
+{144.0, -26.565051177, "vertex 11 r4 10 0 0"},
+{0.0, -90.0, "vertex 12 r4 11 0 0"},
+{0.0, 0.0, NULL}
+};
--- /dev/null
+struct dgg {
+ double a, b;
+ char *data;
+} dgg[] = {
+{2.4038229, 1.7320508, "-170 0"},
+{2.4807823, 1.7320508, ""},
+{2.5576709, 1.7320508, ""},
+{2.6347720, 1.7320508, ""},
+{2.7123722, 1.7320508, ""},
+{2.7907653, 1.7320508, ""},
+{2.8702566, 1.7320508, ""},
+{2.9511682, 1.7320508, ""},
+{0.5262954, 0.8812070, ""},
+{0.5915988, 0.9189100, ""},
+{0.6576797, 0.9570618, ""},
+{0.7244121, 0.9955898, ""},
+{0.7916692, 1.0344207, ""},
+{0.8593228, 1.0734805, ""},
+{0.9272443, 1.1126950, ""},
+{0.9953042, 1.1519894, ""},
+{1.0559646, 1.1870117, ""},
+{1.1159695, 1.2216555, ""},
+{1.1757536, 1.2561719, ""},
+{1.2352032, 1.2904951, ""},
+{2.2944727, 1.3247144, ""},
+{2.3541038, 1.3591424, ""},
+{2.4140125, 1.3937308, ""},
+{2.4740850, 1.4284136, ""},
+{2.5387350, 1.4657393, ""},
+{2.6067418, 1.5050031, ""},
+{2.6745455, 1.5441496, ""},
+{2.7420169, 1.5831042, ""},
+{2.8090276, 1.6217929, ""},
+{2.8754501, 1.6601419, ""},
+{2.9411579, 1.6980783, ""},
+{3.0077369, 1.7320508, ""},
+{3.0894860, 1.7320508, ""},
+{3.1696460, 1.7320508, ""},
+{3.2485491, 1.7320508, ""},
+{3.3265084, 1.7320508, "Origin"},
+{3.4038229, 1.7320508, ""},
+{3.4807823, 1.7320508, ""},
+{3.5576709, 1.7320508, ""},
+{3.6347720, 1.7320508, ""},
+{3.7123722, 1.7320508, ""},
+{3.7907653, 1.7320508, ""},
+{3.8702566, 1.7320508, ""},
+{3.9511682, 1.7320508, ""},
+{4.0262954, 1.7168692, ""},
+{4.0915988, 1.6791662, ""},
+{4.1576797, 1.6410144, ""},
+{4.2244121, 1.6024864, ""},
+{4.2916692, 1.5636555, ""},
+{4.3593228, 1.5245957, ""},
+{4.4272443, 1.4853812, ""},
+{4.4953042, 1.4460868, ""},
+{4.5559646, 1.4110645, ""},
+{4.6159695, 1.3764207, ""},
+{4.6757536, 1.3419043, ""},
+{4.7352032, 1.3075811, ""},
+{3.5000000, 0.0513527, ""},
+{3.5000000, 0.1202087, ""},
+{3.5000000, 0.1893854, ""},
+{3.5000000, 0.2587510, ""},
+{3.5000000, 0.3334024, ""},
+{3.5000000, 0.4119300, ""},
+{3.5000000, 0.4902229, ""},
+{3.5000000, 0.5681322, ""},
+{3.5000000, 0.6455095, ""},
+{3.5000000, 0.7222076, ""},
+{3.5000000, 0.7980804, ""},
+{2.0077369, 1.7320508, ""},
+{2.0894860, 1.7320508, ""},
+{2.1696460, 1.7320508, ""},
+{2.2485491, 1.7320508, ""},
+{2.3265084, 1.7320508, "equator eastern most"},
+{4.7500000, 0.4330127, "south pole"},
+{1.2761529, 0.5048887, "Along zero meridian"},
+{1.3025780, 0.5767523, ""},
+{1.3298119, 0.6489373, ""},
+{1.3585553, 0.7219970, ""},
+{1.3896844, 0.7967593, ""},
+{1.4019639, 0.8614600, ""},
+{3.4116293, 0.9270811, ""},
+{3.4000858, 1.0044103, ""},
+{3.3908933, 1.0805036, ""},
+{3.3823070, 1.1567335, ""},
+{3.3741720, 1.2327090, ""},
+{3.3662149, 1.3080960, ""},
+{3.3578317, 1.3825880, ""},
+{3.3475739, 1.4562424, ""},
+{3.3318447, 1.5308112, ""},
+{3.3292870, 1.5973969, ""},
+{3.3278348, 1.6644314, ""},
+{3.3265084, 1.7320508, "Origin"},
+{2.3278348, 1.7996702, ""},
+{2.3292870, 1.8667047, ""},
+{2.3318447, 1.9332904, ""},
+{2.3475739, 2.0078592, ""},
+{2.3578317, 2.0815136, ""},
+{2.3662149, 2.1560056, ""},
+{2.3741720, 2.2313926, ""},
+{2.3823070, 2.3073681, ""},
+{2.3908933, 2.3835981, ""},
+{1.5698877, 2.4423555, ""},
+{1.5086905, 2.4910171, ""},
+{0.5450643, 2.5108918, ""},
+{0.4951716, 2.4679071, ""},
+{0.4459901, 2.4035673, ""},
+{0.3970902, 2.3421449, ""},
+{0.3481931, 2.2824671, ""},
+{0.2991700, 2.2236506, ""},
+{4.0000000, 2.5980762, "North Pole"},
+{4.7500000, 0.4330127, "South Pole"},
+{4.7991700, 0.3744256, "Along 180E meridian"},
+{4.8481931, 0.3156091, ""},
+{4.8970902, 0.2559313, ""},
+{4.9459901, 0.1945089, ""},
+{4.9951716, 0.1301692, ""},
+{5.0450643, 0.0871844, ""},
+{0.5086905, 0.9730845, ""},
+{0.5698877, 1.0217461, ""},
+{2.3908933, 1.0805036, ""},
+{2.3823070, 1.1567335, ""},
+{2.3741720, 1.2327090, ""},
+{2.3662149, 1.3080960, ""},
+{2.3578317, 1.3825880, ""},
+{2.3475739, 1.4562424, ""},
+{2.3318447, 1.5308112, ""},
+{2.3292870, 1.5973969, ""},
+{2.3278348, 1.6644314, ""},
+{2.3265084, 1.7320508, "equator at 180E"},
+{2.3278348, 1.7996702, ""},
+{2.3292870, 1.8667047, ""},
+{2.3318447, 1.9332904, ""},
+{2.3475739, 2.0078592, ""},
+{2.3578317, 2.0815136, ""},
+{2.3662149, 2.1560056, ""},
+{2.3741720, 2.2313926, ""},
+{2.3823070, 2.3073681, ""},
+{2.3908933, 2.3835981, ""},
+{5.8301981, 0.0173358, ""},
+{5.9029388, 0.0460034, ""},
+{4.4470283, 2.5154572, ""},
+{4.3848560, 2.5371732, ""},
+{4.3045453, 2.5475957, ""},
+{4.2269021, 2.5592331, ""},
+{4.1507711, 2.5717403, ""},
+{4.0753229, 2.5847873, ""},
+{4.0000000, 2.5980762, "North Pole"},
+{0,0,0}
+};
--- /dev/null
+2.4038229 1.7320508 -170 0
+2.4807823 1.7320508
+2.5576709 1.7320508
+2.6347720 1.7320508
+2.7123722 1.7320508
+2.7907653 1.7320508
+2.8702566 1.7320508
+2.9511682 1.7320508
+0.5262954 0.8812070
+0.5915988 0.9189100
+0.6576797 0.9570618
+0.7244121 0.9955898
+0.7916692 1.0344207
+0.8593228 1.0734805
+0.9272443 1.1126950
+0.9953042 1.1519894
+1.0559646 1.1870117
+1.1159695 1.2216555
+1.1757536 1.2561719
+1.2352032 1.2904951
+2.2944727 1.3247144
+2.3541038 1.3591424
+2.4140125 1.3937308
+2.4740850 1.4284136
+2.5387350 1.4657393
+2.6067418 1.5050031
+2.6745455 1.5441496
+2.7420169 1.5831042
+2.8090276 1.6217929
+2.8754501 1.6601419
+2.9411579 1.6980783
+3.0077369 1.7320508
+3.0894860 1.7320508
+3.1696460 1.7320508
+3.2485491 1.7320508
+3.3265084 1.7320508 Origin
+3.4038229 1.7320508
+3.4807823 1.7320508
+3.5576709 1.7320508
+3.6347720 1.7320508
+3.7123722 1.7320508
+3.7907653 1.7320508
+3.8702566 1.7320508
+3.9511682 1.7320508
+4.0262954 1.7168692
+4.0915988 1.6791662
+4.1576797 1.6410144
+4.2244121 1.6024864
+4.2916692 1.5636555
+4.3593228 1.5245957
+4.4272443 1.4853812
+4.4953042 1.4460868
+4.5559646 1.4110645
+4.6159695 1.3764207
+4.6757536 1.3419043
+4.7352032 1.3075811
+3.5000000 0.0513527
+3.5000000 0.1202087
+3.5000000 0.1893854
+3.5000000 0.2587510
+3.5000000 0.3334024
+3.5000000 0.4119300
+3.5000000 0.4902229
+3.5000000 0.5681322
+3.5000000 0.6455095
+3.5000000 0.7222076
+3.5000000 0.7980804
+2.0077369 1.7320508
+2.0894860 1.7320508
+2.1696460 1.7320508
+2.2485491 1.7320508
+2.3265084 1.7320508 equator eastern most
+4.7500000 0.4330127 south pole
+1.2761529 0.5048887 Along zero meridian
+1.3025780 0.5767523
+1.3298119 0.6489373
+1.3585553 0.7219970
+1.3896844 0.7967593
+1.4019639 0.8614600
+3.4116293 0.9270811
+3.4000858 1.0044103
+3.3908933 1.0805036
+3.3823070 1.1567335
+3.3741720 1.2327090
+3.3662149 1.3080960
+3.3578317 1.3825880
+3.3475739 1.4562424
+3.3318447 1.5308112
+3.3292870 1.5973969
+3.3278348 1.6644314
+3.3265084 1.7320508 Origin
+2.3278348 1.7996702
+2.3292870 1.8667047
+2.3318447 1.9332904
+2.3475739 2.0078592
+2.3578317 2.0815136
+2.3662149 2.1560056
+2.3741720 2.2313926
+2.3823070 2.3073681
+2.3908933 2.3835981
+1.5698877 2.4423555
+1.5086905 2.4910171
+0.5450643 2.5108918
+0.4951716 2.4679071
+0.4459901 2.4035673
+0.3970902 2.3421449
+0.3481931 2.2824671
+0.2991700 2.2236506
+4.0000000 2.5980762 North Pole
+4.7500000 0.4330127 South Pole
+4.7991700 0.3744256 Along 180E meridian
+4.8481931 0.3156091
+4.8970902 0.2559313
+4.9459901 0.1945089
+4.9951716 0.1301692
+5.0450643 0.0871844
+0.5086905 0.9730845
+0.5698877 1.0217461
+2.3908933 1.0805036
+2.3823070 1.1567335
+2.3741720 1.2327090
+2.3662149 1.3080960
+2.3578317 1.3825880
+2.3475739 1.4562424
+2.3318447 1.5308112
+2.3292870 1.5973969
+2.3278348 1.6644314
+2.3265084 1.7320508 equator at 180E
+2.3278348 1.7996702
+2.3292870 1.8667047
+2.3318447 1.9332904
+2.3475739 2.0078592
+2.3578317 2.0815136
+2.3662149 2.1560056
+2.3741720 2.2313926
+2.3823070 2.3073681
+2.3908933 2.3835981
+5.8301981 0.0173358
+5.9029388 0.0460034
+4.4470283 2.5154572
+4.3848560 2.5371732
+4.3045453 2.5475957
+4.2269021 2.5592331
+4.1507711 2.5717403
+4.0753229 2.5847873
+4.0000000 2.5980762 North Pole
--- /dev/null
+struct input {
+ double a, b;
+ char *data;
+} input[] = {
+{-175, 0, "-170 0"},
+{-170, 0, ""},
+{-165, 0, ""},
+{-160, 0, ""},
+{-155, 0, ""},
+{-150, 0, ""},
+{-145, 0, ""},
+{-140, 0, ""},
+{-135, 0, ""},
+{-130, 0, ""},
+{-125, 0, ""},
+{-120, 0, ""},
+{-115, 0, ""},
+{-110, 0, ""},
+{-105, 0, ""},
+{-100, 0, ""},
+{-95, 0, ""},
+{-90, 0, ""},
+{-85, 0, ""},
+{-80, 0, ""},
+{-75, 0, ""},
+{-70, 0, ""},
+{-65, 0, ""},
+{-60, 0, ""},
+{-55, 0, ""},
+{-50, 0, ""},
+{-45, 0, ""},
+{-40, 0, ""},
+{-35, 0, ""},
+{-30, 0, ""},
+{-25, 0, ""},
+{-20, 0, ""},
+{-15, 0, ""},
+{-10, 0, ""},
+{-5, 0, ""},
+{0, 0, "Origin"},
+{5, 0, ""},
+{10, 0, ""},
+{15, 0, ""},
+{20, 0, ""},
+{25, 0, ""},
+{30, 0, ""},
+{35, 0, ""},
+{40, 0, ""},
+{45, 0, ""},
+{50, 0, ""},
+{55, 0, ""},
+{60, 0, ""},
+{65, 0, ""},
+{70, 0, ""},
+{75, 0, ""},
+{80, 0, ""},
+{85, 0, ""},
+{90, 0, ""},
+{95, 0, ""},
+{100, 0, ""},
+{105, 0, ""},
+{110, 0, ""},
+{115, 0, ""},
+{120, 0, ""},
+{125, 0, ""},
+{130, 0, ""},
+{135, 0, ""},
+{140, 0, ""},
+{145, 0, ""},
+{150, 0, ""},
+{155, 0, ""},
+{160, 0, ""},
+{165, 0, ""},
+{170, 0, ""},
+{175, 0, ""},
+{180, 0, "equator eastern most"},
+{0, -90, "south pole"},
+{0, -85, "Along zero meridian"},
+{0, -80, ""},
+{0, -75, ""},
+{0, -70, ""},
+{0, -65, ""},
+{0, -60, ""},
+{0, -55, ""},
+{0, -50, ""},
+{0, -45, ""},
+{0, -40, ""},
+{0, -35, ""},
+{0, -30, ""},
+{0, -25, ""},
+{0, -20, ""},
+{0, -15, ""},
+{0, -10, ""},
+{0, -5, ""},
+{0, 0, "Origin"},
+{0, 5, ""},
+{0, 10, ""},
+{0, 15, ""},
+{0, 20, ""},
+{0, 25, ""},
+{0, 30, ""},
+{0, 35, ""},
+{0, 40, ""},
+{0, 45, ""},
+{0, 50, ""},
+{0, 55, ""},
+{0, 60, ""},
+{0, 65, ""},
+{0, 70, ""},
+{0, 75, ""},
+{0, 80, ""},
+{0, 85, ""},
+{0, 90, "North Pole"},
+{180, -90, "South Pole"},
+{180, -85, "Along 180E meridian"},
+{180, -80, ""},
+{180, -75, ""},
+{180, -70, ""},
+{180, -65, ""},
+{180, -60, ""},
+{180, -55, ""},
+{180, -50, ""},
+{180, -45, ""},
+{180, -40, ""},
+{180, -35, ""},
+{180, -30, ""},
+{180, -25, ""},
+{180, -20, ""},
+{180, -15, ""},
+{180, -10, ""},
+{180, -5, ""},
+{180, 0, "equator at 180E"},
+{180, 5, ""},
+{180, 10, ""},
+{180, 15, ""},
+{180, 20, ""},
+{180, 25, ""},
+{180, 30, ""},
+{180, 35, ""},
+{180, 40, ""},
+{180, 45, ""},
+{180, 50, ""},
+{180, 55, ""},
+{180, 60, ""},
+{180, 65, ""},
+{180, 70, ""},
+{180, 75, ""},
+{180, 80, ""},
+{180, 85, ""},
+{180, 90, "North Pole"},
+{0,0,0}
+};
--- /dev/null
+-175 0 -175 0
+-170 0 -170 0
+-165 0 -165 0
+-160 0
+-155 0
+-150 0
+-145 0
+-140 0
+-135 0
+-130 0
+-125 0
+-120 0
+-115 0
+-110 0
+-105 0
+-100 0
+-95 0
+-90 0
+-85 0
+-80 0
+-75 0
+-70 0
+-65 0
+-60 0
+-55 0
+-50 0
+-45 0
+-40 0
+-35 0
+-30 0
+-25 0
+-20 0
+-15 0
+-10 0
+-5 0
+0 0 Origin
+5 0
+10 0
+15 0
+20 0
+25 0
+30 0
+35 0
+40 0
+45 0
+50 0
+55 0
+60 0
+65 0
+70 0
+75 0
+80 0
+85 0
+90 0
+95 0
+100 0
+105 0
+110 0
+115 0
+120 0
+125 0
+130 0
+135 0
+140 0
+145 0
+150 0
+155 0
+160 0
+165 0
+170 0
+175 0
+180 0 equator eastern most
+0 -90 south pole
+0 -85 Along zero meridian
+0 -80
+0 -75
+0 -70
+0 -65
+0 -60
+0 -55
+0 -50
+0 -45
+0 -40
+0 -35
+0 -30
+0 -25
+0 -20
+0 -15
+0 -10
+0 -5
+0 0 Origin
+0 5
+0 10
+0 15
+0 20
+0 25
+0 30
+0 35
+0 40
+0 45
+0 50
+0 55
+0 60
+0 65
+0 70
+0 75
+0 80
+0 85
+0 90 North Pole
+180 -90 South Pole
+180 -85 Along 180E meridian
+180 -80
+180 -75
+180 -70
+180 -65
+180 -60
+180 -55
+180 -50
+180 -45
+180 -40
+180 -35
+180 -30
+180 -25
+180 -20
+180 -15
+180 -10
+180 -5
+180 0 equator at 180E
+180 5
+180 10
+180 15
+180 20
+180 25
+180 30
+180 35
+180 40
+180 45
+180 50
+180 55
+180 60
+180 65
+180 70
+180 75
+180 80
+180 85
+180 90 North Pole
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "tap.h"
+#include "hex.h"
+
+void iso(struct hex *h) {
+ hex_iso(h);
+ printf("iso(%d,%d,%d)", h->x, h->y, h->z);
+}
+void xy(struct hex *h) {
+ hex_xy(h);
+ printf("xy(%d,%d)", h->x, h->y);
+}
+
+int main(void) {
+ int x, y;
+
+ struct hex h;
+
+ plan_tests(121);
+
+ for (x = -5; x<6; x++) {
+ for (y = -5 ; y<6;y++) {
+ h.x = x;
+ h.y = y;
+ h.iso = 0;
+ printf("# ");
+ xy(&h); printf(" = "); iso(&h); printf(" = ");
+ xy(&h); printf("\n");
+ ok(h.x == x && h.y == y, "xy -> iso -> xy");
+ }
+ }
+
+ return exit_status();
+
+}
+
--- /dev/null
+#!/usr/bin/perl
+
+my $fields = shift;
+$fields = 3 unless $fields;
+$" = ', ';
+my @data = ();
+
+print "struct data {\n\tdouble a, b;\n\tchar *data;\n} dgg[] = {\n";
+while (<>) {
+ chomp;
+ my @f = split(/\s+/, $_, $fields);
+ #my $s = pop @f;
+ push @f, '';
+ push @data, qq|{@f[0,1], "$f[2]"}|;
+}
+print join(",\n", @data);
+print "\n};\n";
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include "tap.h"
+#include "isea.h"
+
+#include "testing.h"
+
+void iso(struct hex *h) {
+ hex_iso(h);
+ printf("iso(%d,%d,%d)", h->x, h->y, h->z);
+}
+void xy(struct hex *h) {
+ hex_xy(h);
+ printf("xy(%d,%d)", h->x, h->y);
+}
+
+#include "hexes2.pt"
+#include "hexes2.dgg"
+
+int main(void) {
+ int i, j;
+
+ char isea[64], sahr[64];
+ struct isea_dgg g;
+ struct isea_pt xy;
+ struct isea_geo ll;
+
+ isea_grid_init(&g);
+ g.output = ISEA_PLANE;
+ g.radius = ISEA_SCALE;
+
+ i=0;
+ while (input[i].data) i++;
+
+ plan_tests(i * 2);
+
+ for (j=0; j<i; j++) {
+ skip(2, "don't have known output");
+ continue;
+ /* convert to radians */
+ ll.lon = input[j].a * M_PI / 180.0;
+ ll.lat = input[j].b * M_PI / 180.0;
+ xy = isea_forward(&g, &ll);
+ sprintf(isea, "%.7f", xy.x);
+ sprintf(sahr, "%.7f", dgg[j].a);
+ ok(!strcmp(isea,sahr), "plane X (%f, %f %s) %s != %s / %s",
+ input[j].a, input[j].b, input[j].data,
+ isea, sahr,
+ xy.x, dgg[j].a, dgg[j].data);
+ sprintf(isea, "%.7f", xy.y);
+ sprintf(sahr, "%.7f", dgg[j].b);
+ ok(!strcmp(isea,sahr), "plane Y (%f, %f %s) %f != %f",
+ input[j].a, input[j].b, input[j].data,
+ xy.y, dgg[j].b);
+ }
+
+ return exit_status();
+
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include "tap.h"
+#include "isea.h"
+
+#include "testing.h"
+
+#include "cities.pt"
+
+void iso(struct hex *h) {
+ hex_iso(h);
+ printf("iso(%d,%d,%d)", h->x, h->y, h->z);
+}
+void xy(struct hex *h) {
+ hex_xy(h);
+ printf("xy(%d,%d)", h->x, h->y);
+}
+
+struct dgg {
+ int triangle;
+ double x, y;
+ char *data;
+} dgg[] = {
+{5, 0.5485949, 0.1666408, "49900 Corvallis"},
+{5, 0.5498729, 0.1533820, "42300 Aloha"},
+{5, 0.5485337, 0.1538048, "41700 Tigard"},
+{5, 0.5469778, 0.1653373, "41400 Albany"},
+{5, 0.5477250, 0.1538168, "35700 LakeOswego"},
+{5, 0.5484942, 0.1602033, "32600 Keizer"},
+{5, 0.5514087, 0.1580878, "26800 MacMinnville"},
+{5, 0.5463249, 0.1543639, "26100 OregonCity"},
+{5, 0.5366604, 0.1942501, "23300 GrantsPass"},
+{5, 0.5482398, 0.1544470, "23100 Tualatin"},
+{5, 0.5468969, 0.1541197, "22500 WestLinn"},
+{5, 0.5471126, 0.1531509, "20700 Milwaukie"},
+{5, 0.5338688, 0.1458412, "20600 CitrusPark"},
+{5, 0.5477809, 0.1577168, "20400 Woodburn"},
+{5, 0.5415903, 0.1843270, "20300 Roseburg"},
+{0, 0.0, 0.0, NULL}
+};
+
+int main(void) {
+ int i, j;
+
+ char isea[64], sahr[64];
+ struct isea_dgg g;
+ struct isea_pt xy;
+ struct isea_geo ll;
+
+ isea_grid_init(&g);
+ g.output = ISEA_PROJTRI;
+ g.radius = ISEA_SCALE;
+
+ i=0;
+ while (cities[i].data) i++;
+
+ plan_tests(i * 3);
+
+ for (j=0; j<i; j++) {
+ ll.lon = cities[j].lon * M_PI / 180.0;
+ ll.lat = cities[j].lat * M_PI / 180.0;
+ xy = isea_forward(&g, &ll);
+ sprintf(isea, "%.7f", xy.x);
+ sprintf(sahr, "%.7f", dgg[j].x);
+ ok(!strcmp(isea,sahr), "projtri X (%f, %f %s) %s != %s",
+ cities[j].lon, cities[j].lat, cities[j].data,
+ isea, sahr,
+ xy.x, dgg[j].x);
+ ok(!strcmp(isea,sahr), "projtri Y (%f, %f %s) %f != %f",
+ cities[j].lon, cities[j].lat, cities[j].data,
+ xy.y, dgg[j].y);
+ ok(g.triangle-1 == dgg[j].triangle,
+ "projtri T (%f, %f %s) %d != %d",
+ cities[j].lon, cities[j].lat, cities[j].data,
+ g.triangle-1, dgg[j].triangle);
+ }
+
+ return exit_status();
+
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include "tap.h"
+#include "isea.h"
+
+#include "testing.h"
+
+#include "cities.pt"
+
+struct dgg {
+ int quad;
+ double x, y;
+ char *data;
+} dgg[] = {
+{1, 0.3700177, 0.3076079, "49900 Corvallis"},
+{1, 0.3578962, 0.3131305, "42300 Aloha"},
+{1, 0.3589320, 0.3140789, "41700 Tigard"},
+{1, 0.3696974, 0.3096601, "41400 Albany"},
+{1, 0.3593468, 0.3147732, "35700 LakeOswego"},
+{1, 0.3644930, 0.3109138, "32600 Keizer"},
+{1, 0.3612037, 0.3094476, "26800 MacMinnville"},
+{1, 0.3605206, 0.3157122, "26100 OregonCity"},
+{1, 0.3998953, 0.3041388, "23300 GrantsPass"},
+{1, 0.3596352, 0.3140123, "23100 Tualatin"},
+{1, 0.3600231, 0.3153389, "22500 WestLinn"},
+{1, 0.3590763, 0.3156365, "20700 Milwaukie"},
+{1, 0.3593678, 0.3307608, "20600 CitrusPark"},
+{1, 0.3626963, 0.3127749, "20400 Woodburn"},
+{1, 0.3888368, 0.3048309, "20300 Roseburg"},
+{1, 0.0, 0.0, NULL}
+};
+
+int main(void) {
+ int i, j;
+
+ char isea[64], sahr[64];
+ struct isea_dgg g;
+ struct isea_pt xy;
+ struct isea_geo ll;
+
+ isea_grid_init(&g);
+ g.output = ISEA_Q2DD;
+ g.radius = ISEA_SCALE;
+
+ i=0;
+ while (cities[i].data) i++;
+
+ plan_tests(i * 3);
+
+ for (j=0; j<i; j++) {
+ ll.lon = cities[j].lon * M_PI / 180.0;
+ ll.lat = cities[j].lat * M_PI / 180.0;
+ xy = isea_forward(&g, &ll);
+ sprintf(isea, "%.7f", xy.x);
+ sprintf(sahr, "%.7f", dgg[j].x);
+ ok(!strcmp(isea,sahr), "q2dd X (%f, %f %s) %s == %s",
+ cities[j].lon, cities[j].lat, cities[j].data,
+ isea, sahr,
+ xy.x, dgg[j].x);
+ ok(!strcmp(isea,sahr), "q2dd Y (%f, %f %s) %f == %f",
+ cities[j].lon, cities[j].lat, cities[j].data,
+ xy.y, dgg[j].y);
+ ok(g.quad == dgg[j].quad,
+ "q2dd Q (%f, %f %s) %d == %d",
+ cities[j].lon, cities[j].lat, cities[j].data,
+ g.quad+1, dgg[j].quad);
+ }
+
+ return exit_status();
+
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include "tap.h"
+#include "isea.h"
+
+#include "testing.h"
+
+#include "hexes.pt"
+
+struct dgg {
+ int quad, x, y;
+ char *data;
+} dgg[] = {
+{0, 0, 0, "vertex 1"},
+{1, 0, 0, "vertex 2 r4 1 0 0"},
+{1, 0, 1, "step 1 r4 1 0 1"},
+{1, 0, 2, "step 2 r4 1 0 2"},
+{1, 0, 3, "step 3 r4 1 0 3"},
+{1, 3, 6, "triangle center 1 r4 1 3 6"},
+{2, 0, 0, "vertex 3"},
+{2, 0, 1, "step 3.1"},
+{3, 0, 0, "vertex 4"},
+{4, 0, 0, "vertex 5"},
+{5, 0, 0, "vertex 6"},
+{6, 0, 0, "vertex 7"},
+{7, 0, 0, "vertex 8"},
+{8, 0, 0, "vertex 9"},
+{9, 0, 0, "vertex 10"},
+{10, 0, 0, "vertex 11"},
+{11, 0, 0, "vertex 12"},
+{0, 0, 0, NULL}
+};
+
+int main(void) {
+ int i, j;
+
+ struct isea_dgg g;
+ struct isea_pt xy;
+ struct isea_geo ll;
+
+ isea_grid_init(&g);
+ isea_orient_pole(&g);
+ g.output = ISEA_Q2DI;
+ g.radius = ISEA_SCALE;
+ g.aperture = 3;
+ g.resolution = 4;
+
+ i=0;
+ while (input[i].data) i++;
+
+ plan_tests(i * 3);
+
+ for (j=0; j<i; j++) {
+ ll.lon = input[j].lon * M_PI / 180.0;
+ ll.lat = input[j].lat * M_PI / 180.0;
+ xy = isea_forward(&g, &ll);
+ ok(xy.x == dgg[j].x, "q2di X %d == %d (%f, %f %s)",
+ (int)xy.x, dgg[j].x,
+ input[j].lon, input[j].lat, input[j].data
+ );
+ ok(xy.y == dgg[j].y, "q2di Y %d == %d", (int)xy.y, dgg[j].y);
+ ok(g.quad == dgg[j].quad,
+ "q2di Q %d == %d",
+ g.quad, dgg[j].quad);
+ }
+
+ return exit_status();
+
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include "tap.h"
+#include "isea.h"
+
+#include "testing.h"
+
+#include "hexes.pt"
+
+struct dgg {
+ int quad, x, y;
+ char *data;
+} dgg[] = {
+{0, 0, 0, "vertex 1"},
+{1, 0, 0, "vertex 2 r4 1 0 0"},
+{1, 1, 2, "step 1 r4 1 0 1"},
+{1, 0, 3, "step 2 r4 1 0 2"},
+{1, 0, 3, "step 3 r4 1 0 3"},
+{1, 3, 6, "triangle center 1 r4 1 3 6"},
+{2, 0, 0, "vertex 3"},
+{2, 1, 2, "step 3.1"},
+{3, 0, 0, "vertex 4"},
+{4, 0, 0, "vertex 5"},
+{5, 0, 0, "vertex 6"},
+{6, 0, 0, "vertex 7"},
+{7, 0, 0, "vertex 8"},
+{8, 0, 0, "vertex 9"},
+{9, 0, 0, "vertex 10"},
+{10, 0, 0, "vertex 11"},
+{11, 0, 0, "vertex 12"},
+{0, 0, 0, NULL}
+};
+
+int main(void) {
+ int i, j;
+
+ struct isea_dgg g;
+ struct isea_pt xy;
+ struct isea_geo ll;
+
+ isea_grid_init(&g);
+ isea_orient_pole(&g);
+ g.output = ISEA_Q2DI;
+ g.radius = ISEA_SCALE;
+ g.aperture = 3;
+ g.resolution = 3;
+
+ i=0;
+ while (input[i].data) i++;
+
+ plan_tests(i * 3-3);
+
+ for (j=0; j<i; j++) {
+ if (j==7) continue; /* exactly on line,
+ could bin to either 2,1,2
+ or 1, 7, 8, depending on
+ exact rounding
+ */
+ ll.lon = input[j].lon * M_PI / 180.0;
+ ll.lat = input[j].lat * M_PI / 180.0;
+ xy = isea_forward(&g, &ll);
+ ok(xy.x == dgg[j].x, "q2di X %d == %d (%f, %f %s)",
+ (int)xy.x, dgg[j].x,
+ input[j].lon, input[j].lat, input[j].data
+ );
+ ok(xy.y == dgg[j].y, "q2di Y %d == %d", (int)xy.y, dgg[j].y);
+ ok(g.quad == dgg[j].quad,
+ "q2di Q %d == %d",
+ g.quad, dgg[j].quad);
+ }
+
+ return exit_status();
+
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include "tap.h"
+#include "isea.h"
+
+#include "testing.h"
+
+#include "hexes.pt"
+
+struct dgg {
+ int seq_a3r4, seq_a3r3;
+ char *data;
+} dgg[] = {
+{1, 1, "vertex 1"},
+{2, 2, "vertex 2 r4 1 0 0"},
+{3, 5, "step 1 r4 1 0 1"},
+{4, 3, "step 2 r4 1 0 2"},
+{5, 3, "step 3 r4 1 0 3"},
+{35, 13, "triangle center 1 r4 1 3 6"},
+{83, 29, "vertex 3"},
+{84, 32, "step 3.1"},
+{164, 56, "vertex 4"},
+{245, 83, "vertex 5"},
+{326, 110, "vertex 6"},
+{407, 137, "vertex 7"},
+{488, 164, "vertex 8"},
+{569, 191, "vertex 9"},
+{650, 218, "vertex 10"},
+{731, 245, "vertex 11"},
+{812, 272, "vertex 12"},
+{0, 0, NULL}
+};
+
+int main(void) {
+ int i, j;
+
+ struct isea_dgg g;
+ struct isea_pt xy;
+ struct isea_geo ll;
+
+ isea_grid_init(&g);
+ isea_orient_pole(&g);
+ g.output = ISEA_SEQNUM;
+ g.radius = ISEA_SCALE;
+ g.aperture = 3;
+ g.resolution = 4;
+
+ i=0;
+ while (input[i].data) i++;
+
+ plan_tests(i * 2 - 1);
+
+ for (j=0; j<i; j++) {
+ ll.lon = input[j].lon * M_PI / 180.0;
+ ll.lat = input[j].lat * M_PI / 180.0;
+ xy = isea_forward(&g, &ll);
+ ok(g.serial == dgg[j].seq_a3r4,
+ "seqnum a3r4 %d == %d (%f, %f %s)",
+ g.serial, dgg[j].seq_a3r4,
+ input[j].lon, input[j].lat, input[j].data
+ );
+ }
+
+ g.resolution = 3;
+ for (j=0; j<i; j++) {
+ if (j==7) continue; /* falls on hex edge, so
+ translated hex depends on rounding
+ or precision */
+ ll.lon = input[j].lon * M_PI / 180.0;
+ ll.lat = input[j].lat * M_PI / 180.0;
+ xy = isea_forward(&g, &ll);
+ ok(g.serial == dgg[j].seq_a3r3,
+ "seqnum a3r3 %d == %d (%f, %f %s)",
+ g.serial, dgg[j].seq_a3r3,
+ input[j].lon, input[j].lat, input[j].data
+ );
+ }
+
+ return exit_status();
+
+}
--- /dev/null
+.Dd December 20, 2004
+.Os
+.Dt TAP 3
+.Sh NAME
+.Nm tap
+.Nd write tests that implement the Test Anything Protocol
+.Sh SYNOPSIS
+.In tap.h
+.Sh DESCRIPTION
+The
+.Nm
+library provides functions for writing test scripts that produce output
+consistent with the Test Anything Protocol. A test harness that parses
+this protocol can run these tests and produce useful reports indicating
+their success or failure.
+.Ss PRINTF STRINGS
+In the descriptions that follow, for any function that takes as the
+last two parameters
+.Dq Fa char * , Fa ...
+it can be assumed that the
+.Fa char *
+is a
+.Fn printf
+-like format string, and the optional arguments are values to be placed
+in that string.
+.Ss TEST PLANS
+.Bl -tag -width indent
+.It Xo
+.Ft int
+.Fn plan_tests "unsigned int"
+.Xc
+.It Xo
+.Ft int
+.Fn plan_no_plan "void"
+.Xc
+.It Xo
+.Ft int
+.Fn plan_skip_all "char *" "..."
+.Xc
+.El
+.Pp
+You must first specify a test plan. This indicates how many tests you
+intend to run, and allows the test harness to notice if any tests were
+missed, or if the test program exited prematurely.
+.Pp
+To do this, use
+.Fn plan_tests ,
+which always returns 0. The function will cause your program to exit
+prematurely if you specify 0 tests.
+.Pp
+In some situations you may not know how many tests you will be running, or
+you are developing your test program, and do not want to update the
+.Fn plan_tests
+parameter every time you make a change. For those situations use
+.Fn plan_no_plan .
+It returns 0, and indicates to the test harness that an indeterminate number
+of tests will be run.
+.Pp
+Both
+.Fn plan_tests
+and
+.Fn plan_no_plan
+will cause your test program to exit prematurely with a diagnostic
+message if they are called more than once.
+.Pp
+If your test program detects at run time that some required functionality
+is missing (for example, it relies on a database connection which is not
+present, or a particular configuration option that has not been included
+in the running kernel) use
+.Fn plan_skip_all ,
+passing as parameters a string to display indicating the reason for skipping
+the tests.
+.Ss SIMPLE TESTS
+.Bl -tag -width indent
+.It Xo
+.Ft unsigned int
+.Fn ok "expression" "char *" "..."
+.Xc
+.It Xo
+.Ft unsigned int
+.Fn ok1 "expression"
+.Xc
+.It Xo
+.Ft unsigned int
+.Fn pass "char *" "..."
+.Xc
+.It Xo
+.Ft unsigned int
+.Fn fail "char *" "..."
+.Xc
+.El
+.Pp
+Tests are implemented as expressions checked by calls to the
+.Fn ok
+and
+.Fn ok1
+macros. In both cases
+.Fa expression
+should evaluate to true if the test succeeded.
+.Pp
+.Fn ok
+allows you to specify a name, or comment, describing the test which will
+be included in the output.
+.Fn ok1
+is for those times when the expression to be tested is self
+explanatory and does not need an associated comment. In those cases
+the test expression becomes the comment.
+.Pp
+These four calls are equivalent:
+.Bd -literal -offset indent
+int i = 5;
+
+ok(i == 5, "i equals 5"); /* Overly verbose */
+ok(i == 5, "i equals %d", i); /* Just to demonstrate printf-like
+ behaviour of the test name */
+ok(i == 5, "i == 5"); /* Needless repetition */
+ok1(i == 5); /* Just right */
+.Ed
+.Pp
+It is good practice to ensure that the test name describes the meaning
+behind the test rather than what you are testing. Viz
+.Bd -literal -offset indent
+ok(db != NULL, "db is not NULL"); /* Not bad, but */
+ok(db != NULL, "Database conn. succeeded"); /* this is better */
+.Ed
+.Pp
+.Fn ok
+and
+.Fn ok1
+return 1 if the expression evaluated to true, and 0 if it evaluated to
+false. This lets you chain calls from
+.Fn ok
+to
+.Fn diag
+to only produce diagnostic output if the test failed. For example, this
+code will include diagnostic information about why the database connection
+failed, but only if the test failed.
+.Bd -literal -offset indent
+ok(db != NULL, "Database conn. succeeded") ||
+ diag("Database error code: %d", dberrno);
+.Ed
+.Pp
+You also have
+.Fn pass
+and
+.Fn fail .
+From the Test::More documentation:
+.Bd -literal -offset indent
+Sometimes you just want to say that the tests have passed.
+Usually the case is you've got some complicated condition
+that is difficult to wedge into an ok(). In this case,
+you can simply use pass() (to declare the test ok) or fail
+(for not ok).
+
+Use these very, very, very sparingly.
+.Ed
+.Pp
+These are synonyms for ok(1, ...) and ok(0, ...).
+.Ss SKIPPING TESTS
+.Bl -tag -width indent
+.It Xo
+.Ft int
+.Fn skip "unsigned int" "char *" "..."
+.Xc
+.It Xo
+.Fn skip_start "expression" "unsigned int" "char *" "..."
+.Xc
+.It Xo
+.Sy skip_end
+.Xc
+.El
+.Pp
+Sets of tests can be skipped. Ordinarily you would do this because
+the test can't be run in this particular testing environment.
+.Pp
+For example, suppose some tests should be run as root. If the test is
+not being run as root then the tests should be skipped. In this
+implementation, skipped tests are flagged as being ok, with a special
+message indicating that they were skipped. It is your responsibility
+to ensure that the number of tests skipped (the first parameter to
+.Fn skip )
+is correct for the number of tests to skip.
+.Pp
+One way of implementing this is with a
+.Dq do { } while(0);
+loop, or an
+.Dq if( ) { } else { }
+construct, to ensure that there are no additional side effects from the
+skipped tests.
+.Bd -literal -offset indent
+if(getuid() != 0) {
+ skip(1, "because test only works as root");
+} else {
+ ok(do_something_as_root() == 0, "Did something as root");
+}
+.Ed
+.Pp
+Two macros are provided to assist with this. The previous example could
+be re-written as follows.
+.Bd -literal -offset indent
+skip_start(getuid() != 0, 1, "because test only works as root");
+
+ok(do_something_as_root() == 0, "Did something as root");
+
+skip_end; /* It's a macro, no parentheses */
+.Ed
+.Ss MARKING TESTS AS Dq TODO
+.Bl -tag -width indent
+.It Xo
+.Ft void
+.Fn todo_start "char *" "..."
+.Xc
+.It Xo
+.Ft void
+.Fn todo_end "void"
+.Xc
+.El
+.Pp
+Sets of tests can be flagged as being
+.Dq TODO .
+These are tests that you expect to fail, probably because you haven't
+fixed a bug, or finished a new feature yet. These tests will still be
+run, but with additional output that indicates that they are expected
+to fail. Should a test start to succeed unexpectedly, tools like
+.Xr prove 1
+will indicate this, and you can move the test out of the todo
+block. This is much more useful than simply commenting out (or
+.Dq #ifdef 0 ... #endif )
+the tests.
+.Bd -literal -offset indent
+todo_start("dwim() not returning true yet");
+
+ok(dwim(), "Did what the user wanted");
+
+todo_end();
+.Ed
+.Pp
+Should
+.Fn dwim
+ever start succeeding you will know about it as soon as you run the
+tests. Note that
+.Em unlike
+the
+.Fn skip_*
+family, additional code between
+.Fn todo_start
+and
+.Fn todo_end
+.Em is
+executed.
+.Ss SKIP vs. TODO
+From the Test::More documentation;
+.Bd -literal -offset indent
+If it's something the user might not be able to do, use SKIP.
+This includes optional modules that aren't installed, running
+under an OS that doesn't have some feature (like fork() or
+symlinks), or maybe you need an Internet connection and one
+isn't available.
+
+If it's something the programmer hasn't done yet, use TODO.
+This is for any code you haven't written yet, or bugs you have
+yet to fix, but want to put tests in your testing script
+(always a good idea).
+.Ed
+.Ss DIAGNOSTIC OUTPUT
+.Bl -tag -width indent
+.It Xo
+.Fr unsigned int
+.Fn diag "char *" "..."
+.Xc
+.El
+.Pp
+If your tests need to produce diagnostic output, use
+.Fn diag .
+It ensures that the output will not be considered by the TAP test harness.
+.Fn diag
+adds the necessary trailing
+.Dq \en
+for you.
+.Bd -literal -offset indent
+diag("Expected return code 0, got return code %d", rcode);
+.Ed
+.Pp
+.Fn diag
+always returns 0.
+.Ss EXIT STATUS
+.Bl -tag -width indent
+.It Xo
+.Fr int
+.Fn exit_status void
+.Xc
+.El
+.Pp
+For maximum compatability your test program should return a particular
+exit code. This is calculated by
+.Fn exit_status
+so it is sufficient to always return from
+.Fn main
+with either
+.Dq return exit_status();
+or
+.Dq exit(exit_status());
+as appropriate.
+.Sh EXAMPLES
+The
+.Pa tests
+directory in the source distribution contains numerous tests of
+.Nm
+functionality, written using
+.Nm .
+Examine them for examples of how to construct test suites.
+.Sh COMPATABILITY
+.Nm
+strives to be compatible with the Perl Test::More and Test::Harness
+modules. The test suite verifies that
+.Nm
+is bug-for-bug compatible with their behaviour. This is why some
+functions which would more naturally return nothing return constant
+values.
+.Pp
+If the
+.Lb libpthread
+is found at compile time,
+.Nm
+.Em should
+be thread safe. Indications to the contrary (and test cases that expose
+incorrect behaviour) are very welcome.
+.Sh SEE ALSO
+.Xr Test::More 1 ,
+.Xr Test::Harness 1 ,
+.Xr prove 1
+.Sh STANDARDS
+.Nm
+requires a
+.St -isoC-99
+compiler. Some of the
+.Nm
+functionality is implemented as variadic macros, and that functionality
+was not formally codified until C99. Patches to use
+.Nm
+with earlier compilers that have their own implementation of variadic
+macros will be gratefully received.
+.Sh HISTORY
+.Nm
+was written to help improve the quality and coverage of the FreeBSD
+regression test suite, and released in the hope that others find it
+a useful tool to help improve the quality of their code.
+.Sh AUTHORS
+.An "Nik Clayton" Aq nik@ngo.org.uk ,
+.Aq nik@FreeBSD.org
+.Pp
+.Nm
+would not exist without the efforts of
+.An "Michael G Schwern" Aq schqern@pobox.com ,
+.An "Andy Lester" Aq andy@petdance.com ,
+and the countless others who have worked on the Perl QA programme.
+.Sh BUGS
+Ideally, running the tests would have no side effects on the behaviour
+of the application you are testing. However, it is not always possible
+to avoid them. The following side effects of using
+.Nm
+are known.
+.Bl -bullet -offset indent
+.It
+stdout is set to unbuffered mode after calling any of the
+.Fn plan_*
+functions.
+.El
--- /dev/null
+/*-
+ * Copyright (c) 2004 Nik Clayton
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+#define _GNU_SOURCE
+#include <ctype.h>
+#include <stdarg.h>
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "tap.h"
+
+static int no_plan = 0;
+static int skip_all = 0;
+static int have_plan = 0;
+static unsigned int test_count = 0; /* Number of tests that have been run */
+static unsigned int e_tests = 0; /* Expected number of tests to run */
+static unsigned int failures = 0; /* Number of tests that failed */
+static char *todo_msg = NULL;
+static char *todo_msg_fixed = "libtap malloc issue";
+static int todo = 0;
+static int test_died = 0;
+
+/* Encapsulate the pthread code in a conditional. In the absence of
+ libpthread the code does nothing */
+#ifdef HAVE_LIBPTHREAD
+#include <pthread.h>
+static pthread_mutex_t M = PTHREAD_MUTEX_INITIALIZER;
+# define LOCK pthread_mutex_lock(&M);
+# define UNLOCK pthread_mutex_unlock(&M);
+#else
+# define LOCK
+# define UNLOCK
+#endif
+
+static void _expected_tests(unsigned int);
+static void _tap_init(void);
+static void _cleanup(void);
+
+/*
+ * Generate a test result.
+ *
+ * ok -- boolean, indicates whether or not the test passed.
+ * test_name -- the name of the test, may be NULL
+ * test_comment -- a comment to print afterwards, may be NULL
+ */
+unsigned int
+_gen_result(int ok, const char *func, char *file, unsigned int line,
+ char *test_name, ...)
+{
+ va_list ap;
+ char *local_test_name = NULL;
+ char *c;
+ int name_is_digits;
+
+ LOCK;
+
+ test_count++;
+
+ /* Start by taking the test name and performing any printf()
+ expansions on it */
+ if(test_name != NULL) {
+ va_start(ap, test_name);
+ vasprintf(&local_test_name, test_name, ap);
+ va_end(ap);
+
+ /* Make sure the test name contains more than digits
+ and spaces. Emit an error message and exit if it
+ does */
+ if(local_test_name) {
+ name_is_digits = 1;
+ for(c = local_test_name; *c != '\0'; c++) {
+ if(!isdigit(*c) && !isspace(*c)) {
+ name_is_digits = 0;
+ break;
+ }
+ }
+
+ if(name_is_digits) {
+ diag(" You named your test '%s'. You shouldn't use numbers for your test names.", local_test_name);
+ diag(" Very confusing.");
+ }
+ }
+ }
+
+ if(!ok) {
+ printf("not ");
+ failures++;
+ }
+
+ printf("ok %d", test_count);
+
+ if(test_name != NULL) {
+ printf(" - ");
+
+ /* Print the test name, escaping any '#' characters it
+ might contain */
+ if(local_test_name != NULL) {
+ flockfile(stdout);
+ for(c = local_test_name; *c != '\0'; c++) {
+ if(*c == '#')
+ fputc('\\', stdout);
+ fputc((int)*c, stdout);
+ }
+ funlockfile(stdout);
+ } else { /* vasprintf() failed, use a fixed message */
+ printf("%s", todo_msg_fixed);
+ }
+ }
+
+ /* If we're in a todo_start() block then flag the test as being
+ TODO. todo_msg should contain the message to print at this
+ point. If it's NULL then asprintf() failed, and we should
+ use the fixed message.
+
+ This is not counted as a failure, so decrement the counter if
+ the test failed. */
+ if(todo) {
+ printf(" # TODO %s", todo_msg ? todo_msg : todo_msg_fixed);
+ if(!ok)
+ failures--;
+ }
+
+ printf("\n");
+
+ if(!ok)
+ diag(" Failed %stest (%s:%s() at line %d)",
+ todo ? "(TODO) " : "", file, func, line);
+
+ free(local_test_name);
+
+ UNLOCK;
+
+ /* We only care (when testing) that ok is positive, but here we
+ specifically only want to return 1 or 0 */
+ return ok ? 1 : 0;
+}
+
+/*
+ * Initialise the TAP library. Will only do so once, however many times it's
+ * called.
+ */
+void
+_tap_init(void)
+{
+ static int run_once = 0;
+
+ LOCK;
+
+ if(!run_once) {
+ atexit(_cleanup);
+
+ /* stdout needs to be unbuffered so that the output appears
+ in the same place relative to stderr output as it does
+ with Test::Harness */
+ setbuf(stdout, 0);
+ run_once = 1;
+ }
+
+ UNLOCK;
+}
+
+/*
+ * Note that there's no plan.
+ */
+int
+plan_no_plan(void)
+{
+
+ LOCK;
+
+ _tap_init();
+
+ if(have_plan != 0) {
+ fprintf(stderr, "You tried to plan twice!\n");
+ test_died = 1;
+ UNLOCK;
+ exit(255);
+ }
+
+ have_plan = 1;
+ no_plan = 1;
+
+ UNLOCK;
+
+ return 0;
+}
+
+/*
+ * Note that the plan is to skip all tests
+ */
+int
+plan_skip_all(char *reason)
+{
+
+ LOCK;
+
+ _tap_init();
+
+ skip_all = 1;
+
+ printf("1..0");
+
+ if(reason != NULL)
+ printf(" # Skip %s", reason);
+
+ printf("\n");
+
+ UNLOCK;
+
+ exit(0);
+}
+
+/*
+ * Note the number of tests that will be run.
+ */
+int
+plan_tests(unsigned int tests)
+{
+
+ LOCK;
+
+ _tap_init();
+
+ if(have_plan != 0) {
+ fprintf(stderr, "You tried to plan twice!\n");
+ test_died = 1;
+ UNLOCK;
+ exit(255);
+ }
+
+ if(tests == 0) {
+ fprintf(stderr, "You said to run 0 tests! You've got to run something.\n");
+ test_died = 1;
+ UNLOCK;
+ exit(255);
+ }
+
+ have_plan = 1;
+
+ _expected_tests(tests);
+
+ UNLOCK;
+
+ return 0;
+}
+
+unsigned int
+diag(char *fmt, ...)
+{
+ va_list ap;
+
+ LOCK;
+
+ fputs("# ", stderr);
+
+ va_start(ap, fmt);
+ vfprintf(stderr, fmt, ap);
+ va_end(ap);
+
+ fputs("\n", stderr);
+
+ UNLOCK;
+
+ return 0;
+}
+
+void
+_expected_tests(unsigned int tests)
+{
+
+ LOCK;
+
+ printf("1..%d\n", tests);
+ e_tests = tests;
+
+ UNLOCK;
+}
+
+int
+skip(unsigned int n, char *fmt, ...)
+{
+ va_list ap;
+ char *skip_msg;
+
+ LOCK;
+
+ va_start(ap, fmt);
+ asprintf(&skip_msg, fmt, ap);
+ va_end(ap);
+
+ while(n-- > 0) {
+ test_count++;
+ printf("ok %d # skip %s\n", test_count,
+ skip_msg != NULL ?
+ skip_msg : "libtap():malloc() failed");
+ }
+
+ free(skip_msg);
+
+ UNLOCK;
+
+ return 1;
+}
+
+void
+todo_start(char *fmt, ...)
+{
+ va_list ap;
+
+ LOCK;
+
+ va_start(ap, fmt);
+ vasprintf(&todo_msg, fmt, ap);
+ va_end(ap);
+
+ todo = 1;
+
+ UNLOCK;
+}
+
+void
+todo_end(void)
+{
+
+ LOCK;
+
+ todo = 0;
+ free(todo_msg);
+
+ UNLOCK;
+}
+
+int
+exit_status(void)
+{
+ int r;
+
+ LOCK;
+
+ /* If there's no plan, just return the number of failures */
+ if(no_plan || !have_plan) {
+ UNLOCK;
+ return failures;
+ }
+
+ /* Ran too many tests? Return the number of tests that were run
+ that shouldn't have been */
+ if(e_tests < test_count) {
+ r = test_count - e_tests;
+ UNLOCK;
+ return r;
+ }
+
+ /* Return the number of tests that failed + the number of tests
+ that weren't run */
+ r = failures + e_tests - test_count;
+ UNLOCK;
+
+ return r;
+}
+
+/*
+ * Cleanup at the end of the run, produce any final output that might be
+ * required.
+ */
+void
+_cleanup(void)
+{
+
+ LOCK;
+
+ /* If plan_no_plan() wasn't called, and we don't have a plan,
+ and we're not skipping everything, then something happened
+ before we could produce any output */
+ if(!no_plan && !have_plan && !skip_all) {
+ diag("Looks like your test died before it could output anything.");
+ UNLOCK;
+ return;
+ }
+
+ if(test_died) {
+ diag("Looks like your test died just after %d.", test_count);
+ UNLOCK;
+ return;
+ }
+
+
+ /* No plan provided, but now we know how many tests were run, and can
+ print the header at the end */
+ if(!skip_all && (no_plan || !have_plan)) {
+ printf("1..%d\n", test_count);
+ }
+
+ if((have_plan && !no_plan) && e_tests < test_count) {
+ diag("Looks like you planned %d tests but ran %d extra.",
+ e_tests, test_count - e_tests);
+ UNLOCK;
+ return;
+ }
+
+ if((have_plan || !no_plan) && e_tests > test_count) {
+ diag("Looks like you planned %d tests but only ran %d.",
+ e_tests, test_count);
+ UNLOCK;
+ return;
+ }
+
+ if(failures)
+ diag("Looks like you failed %d tests of %d.",
+ failures, test_count);
+
+ UNLOCK;
+}
--- /dev/null
+/*-
+ * Copyright (c) 2004 Nik Clayton
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+ * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+ * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+ * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+ * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+ * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+ * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+ * SUCH DAMAGE.
+ */
+
+/* '## __VA_ARGS__' is a gcc'ism. C99 doesn't allow the token pasting
+ and requires the caller to add the final comma if they've ommitted
+ the optional arguments */
+#ifdef __GNUC__
+# define ok(e, test, ...) ((e) ? \
+ _gen_result(1, __func__, __FILE__, __LINE__, \
+ test, ## __VA_ARGS__) : \
+ _gen_result(0, __func__, __FILE__, __LINE__, \
+ test, ## __VA_ARGS__))
+
+# define ok1(e) ((e) ? \
+ _gen_result(1, __func__, __FILE__, __LINE__, "%s", #e) : \
+ _gen_result(0, __func__, __FILE__, __LINE__, "%s", #e))
+
+# define pass(test, ...) ok(1, test, ## __VA_ARGS__);
+# define fail(test, ...) ok(0, test, ## __VA_ARGS__);
+
+# define skip_start(test, n, fmt, ...) \
+ do { \
+ if((test)) { \
+ skip(n, fmt, ## __VA_ARGS__); \
+ continue; \
+ }
+#elif __STDC_VERSION__ >= 199901L /* __GNUC__ */
+# define ok(e, ...) ((e) ? \
+ _gen_result(1, __func__, __FILE__, __LINE__, \
+ __VA_ARGS__) : \
+ _gen_result(0, __func__, __FILE__, __LINE__, \
+ __VA_ARGS__))
+
+# define ok1(e) ((e) ? \
+ _gen_result(1, __func__, __FILE__, __LINE__, "%s", #e) : \
+ _gen_result(0, __func__, __FILE__, __LINE__, "%s", #e))
+
+# define pass(...) ok(1, __VA_ARGS__);
+# define fail(...) ok(0, __VA_ARGS__);
+
+# define skip_start(test, n, ...) \
+ do { \
+ if((test)) { \
+ skip(n, __VA_ARGS__); \
+ continue; \
+ }
+#else /* __STDC_VERSION__ */
+# error "Needs gcc or C99 compiler for variadic macros."
+#endif /* __STDC_VERSION__ */
+
+# define skip_end } while(0);
+
+unsigned int _gen_result(int, const char *, char *, unsigned int, char *, ...);
+
+int plan_no_plan(void);
+int plan_skip_all(char *);
+int plan_tests(unsigned int);
+
+unsigned int diag(char *, ...);
+
+int skip(unsigned int, char *, ...);
+
+void todo_start(char *, ...);
+void todo_end(void);
+
+int exit_status(void);
--- /dev/null
+struct points {
+ double lon, lat;
+ char *data;
+};
--- /dev/null
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+
+#include "isea.h"
+#include "hex.h"
+
+/*
+ * options
+ *
+ * orientation
+ * -l longitude
+ * -z azimuth
+ * -p latitude
+ *
+ * -I isea standard orientation
+ * -D dymaxion orientation
+ * -P vertex at pole
+ *
+ * -a aperture
+ * -r resolution
+ * -t h t d hex tri diamond
+ *
+ * addressing
+ * -O q2dd seqnum interleave plane q2dd projtri vertex2dd
+ * -I
+ *
+ * operation
+ * -g generate grid
+ * -x transform
+ * -b bin point values (average)
+ * -e bin point presence
+ * -G spacing resolution generate graticule
+ */
+
+void usage(void) {
+ fprintf(stderr, "%s", "isea [-l <lon>] [-z <azimuth>] [-p <latitude>] [-I] [-D] [-P] [-a <aperture>] [-r <resolution>] [-t <topology>] [-A <addressing>] [-g] [-b] [-e] [-x]\n");
+}
+
+void die(char *s) {
+ fprintf(stderr, "%s\n", s);
+ exit(EXIT_FAILURE);
+}
+
+static int verbose = 0;
+static int precision = 10;
+
+void options(int ac, char **av, struct isea_dgg *g) {
+ int i = 1;
+ char *optarg;
+
+ for (i=1;i<ac;i++) {
+ if (av[i][0] != '-') {
+ break;
+ }
+ switch (av[i][1]) {
+ case 'v':
+ verbose = 1; break;
+ case 'I':
+ isea_orient_isea(g); break;
+ case 'P':
+ isea_orient_pole(g); break;
+ case 'D':
+ isea_orient_dymax(g); break;
+ case 'R':
+ g->radius = strtod(av[++i],NULL); break;
+ case 'l':
+ g->o_lon = strtod(av[++i],NULL) * M_PI / 180.0; break;
+ case 'p':
+ g->o_lat = strtod(av[++i],NULL) * M_PI / 180.0; break;
+ case 'z':
+ g->o_az = strtod(av[++i],NULL) * M_PI / 180.0; break;
+ case 'd':
+ precision = atoi(av[++i]); break;
+ case 'a':
+ g->aperture = atoi(av[++i]); break;
+ case 'r':
+ g->resolution = atoi(av[++i]); break;
+ case 't':
+ switch (av[++i][0]) {
+ case 'h':
+ g->topology = ISEA_HEXAGON;
+ break;
+ case 't':
+ g->topology = ISEA_TRIANGLE;
+ break;
+ case 'd':
+ g->topology = ISEA_DIAMOND;
+ break;
+ default:
+ usage();
+ die("Invalid topology");
+ break;
+ };
+ break;
+ case 'O':
+ if (av[i][2]) {
+ optarg = &av[i][2];
+ } else {
+ optarg = av[++i];
+ }
+ if (!strcmp(optarg, "projtri")) {
+ g->output = ISEA_PROJTRI;
+ }
+ else if (!strcmp(optarg, "vertex2dd")) {
+ g->output = ISEA_VERTEX2DD;
+ }
+ else if (!strcmp(optarg, "q2dd")) {
+ g->output = ISEA_Q2DD;
+ }
+ else if (!strcmp(optarg, "q2di")) {
+ g->output = ISEA_Q2DI;
+ }
+ else if (!strcmp(optarg, "seqnum")) {
+ g->output = ISEA_SEQNUM;
+ }
+ else if (!strcmp(optarg, "plane")) {
+ g->output = ISEA_PLANE;
+ } else {
+ fprintf(stderr, "unknown output format: %s\n",optarg);
+ die("exiting");
+ }
+ break;
+ default:
+ usage();
+ die("Unknown option");
+ }
+ }
+}
+
+int
+main(int ac, char *av[])
+{
+ struct isea_dgg dgg;
+ struct isea_geo in;
+ struct isea_pt out, coord;
+ int tri,calctri;
+ int graticule = 0;
+ int quad = 0;
+
+ char line[1024];
+ char *rest, *nl;
+
+ isea_grid_init(&dgg);
+ dgg.output = ISEA_PLANE;
+ dgg.radius = ISEA_SCALE; /* 0.8301572857837594396028083; */
+
+ options(ac, av, &dgg);
+
+ if (verbose) {
+ fprintf(stderr, "isea lat = %lf, lon = %lf, az = %lf, output = %s\n",
+ dgg.o_lat * 180.0 / M_PI, dgg.o_lon * 180.0 / M_PI,
+ dgg.o_az * 180.0 / M_PI, dgg.output == ISEA_PLANE ? "plane" : "not plane");
+ }
+
+ while (fgets(line, sizeof line, stdin)) {
+ line[sizeof line-1] = 0;
+ if (line[0] == '#') {
+ continue;
+ }
+ if (!strncmp(line, "ENDPOLY", 7)) {
+ printf("%s", line);
+ continue;
+ }
+ in.lon = strtod(line, &rest);
+ in.lat = strtod(rest, &rest);
+
+ nl = strchr(rest, '\n');
+ if (nl) *nl = 0;
+
+ in.lon *= M_PI / 180.0;
+ in.lat *= M_PI / 180.0;
+
+ out = isea_forward(&dgg, &in);
+ tri = dgg.triangle;
+ quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;
+ coord = out;
+
+ if (dgg.output == ISEA_PLANE) {
+ calctri = isea_xy_triangle(&coord,
+ dgg.output == ISEA_PLANE ? dgg.radius : 1.0);
+#if 0
+ if (tri != calctri) {
+ fprintf(stderr,
+ "got %d (%f, %f), expected %d\n",
+ calctri, out.x, out.y, tri);
+ }
+#endif
+ }
+
+ in.lon *= 180.0 / M_PI;
+ in.lat *= 180.0 / M_PI;
+
+ switch (dgg.output) {
+ case ISEA_PLANE:
+ printf("%.*f %.*f%s\n", precision, out.x, precision, out.y, rest);
+ break;
+ case ISEA_PROJTRI:
+ printf("%d %.*f %.*f%s\n", tri-1, precision, out.x,
+ precision, out.y, rest);
+ break;
+ case ISEA_VERTEX2DD:
+ printf("%d %d %.*f %.*f%s\n", quad, tri-1, precision,
+ out.x, precision, out.y, rest);
+ break;
+ case ISEA_Q2DD:
+ /* Same as above, we just don't print as much */
+ printf("%d %.*f %.*f%s\n", quad, precision,
+ out.x, precision, out.y, rest);
+ break;
+ case ISEA_Q2DI:
+ printf("%d %.0lf %.0lf%s\n", dgg.quad, coord.x, coord.y, rest);
+ break;
+ case ISEA_SEQNUM:
+ quad = isea_ptdi(&dgg, tri, &out, &coord);
+ printf("%d%s\n", isea_disn(&dgg, quad, &coord), rest);
+ break;
+ }
+ fflush(stdout);
+ }
+
+ if (graticule) {
+
+ }
+
+ exit(EXIT_SUCCESS);
+}
projects / isea / commitdiff