From 7ed53098e43d4fb96c47345f18750fd588a04240 Mon Sep 17 00:00:00 2001 From: Nathan Wagner Date: Fri, 27 May 2011 07:01:51 +0000 Subject: [PATCH 1/1] Initial checkin. --- Makefile | 84 ++++ README | 26 ++ hex.c | 5 + hex.h | 23 + hfuncs.c | 86 ++++ isea.c | 7 + isea.h | 135 ++++++ pfuncs.c | 1039 ++++++++++++++++++++++++++++++++++++++++++++ pjexamples.c | 324 ++++++++++++++ pjhead.c | 16 + pjtail.c | 105 +++++ proj.c | 10 + proj4.patch | 69 +++ t/cities.c | 110 +++++ t/cities.in | 16 + t/cities.pt | 18 + t/globehex.c | 111 +++++ t/hexbin.c | 43 ++ t/hexes.pt | 23 + t/hexes2.dgg | 152 +++++++ t/hexes2.dgg.plane | 146 +++++++ t/hexes2.pt | 152 +++++++ t/hexes2.txt | 146 +++++++ t/hextest.c | 38 ++ t/mkc | 17 + t/plane.c | 62 +++ t/projtri.c | 83 ++++ t/q2dd.c | 74 ++++ t/q2di.c | 72 +++ t/q2di_r3.c | 77 ++++ t/seqnum.c | 84 ++++ t/tap.3 | 368 ++++++++++++++++ t/tap.c | 436 +++++++++++++++++++ t/tap.h | 89 ++++ t/testing.h | 4 + transform.c | 228 ++++++++++ 36 files changed, 4478 insertions(+) create mode 100644 Makefile create mode 100644 README create mode 100644 hex.c create mode 100644 hex.h create mode 100644 hfuncs.c create mode 100644 isea.c create mode 100644 isea.h create mode 100644 pfuncs.c create mode 100644 pjexamples.c create mode 100644 pjhead.c create mode 100644 pjtail.c create mode 100644 proj.c create mode 100644 proj4.patch create mode 100644 t/cities.c create mode 100644 t/cities.in create mode 100644 t/cities.pt create mode 100644 t/globehex.c create mode 100644 t/hexbin.c create mode 100644 t/hexes.pt create mode 100644 t/hexes2.dgg create mode 100644 t/hexes2.dgg.plane create mode 100644 t/hexes2.pt create mode 100644 t/hexes2.txt create mode 100644 t/hextest.c create mode 100755 t/mkc create mode 100644 t/plane.c create mode 100644 t/projtri.c create mode 100644 t/q2dd.c create mode 100644 t/q2di.c create mode 100644 t/q2di_r3.c create mode 100644 t/seqnum.c create mode 100644 t/tap.3 create mode 100644 t/tap.c create mode 100644 t/tap.h create mode 100644 t/testing.h create mode 100644 transform.c diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..5264e60 --- /dev/null +++ b/Makefile @@ -0,0 +1,84 @@ +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 diff --git a/README b/README new file mode 100644 index 0000000..ad500e4 --- /dev/null +++ b/README @@ -0,0 +1,26 @@ +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. diff --git a/hex.c b/hex.c new file mode 100644 index 0000000..677b768 --- /dev/null +++ b/hex.c @@ -0,0 +1,5 @@ +#include +#include +#include + +#include "hfuncs.c" diff --git a/hex.h b/hex.h new file mode 100644 index 0000000..5940918 --- /dev/null +++ b/hex.h @@ -0,0 +1,23 @@ +#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 diff --git a/hfuncs.c b/hfuncs.c new file mode 100644 index 0000000..931df7f --- /dev/null +++ b/hfuncs.c @@ -0,0 +1,86 @@ +#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; +} diff --git a/isea.c b/isea.c new file mode 100644 index 0000000..510fcde --- /dev/null +++ b/isea.c @@ -0,0 +1,7 @@ +#include +#include +#include +#include + +#include "hfuncs.c" +#include "pfuncs.c" diff --git a/isea.h b/isea.h new file mode 100644 index 0000000..8f96cf7 --- /dev/null +++ b/isea.h @@ -0,0 +1,135 @@ +#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 diff --git a/pfuncs.c b/pfuncs.c new file mode 100644 index 0000000..af581da --- /dev/null +++ b/pfuncs.c @@ -0,0 +1,1039 @@ +#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; +} diff --git a/pjexamples.c b/pjexamples.c new file mode 100644 index 0000000..1f62694 --- /dev/null +++ b/pjexamples.c @@ -0,0 +1,324 @@ +/* + * 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 ] [-z ] [-p ] [-I] [-D] [-P] [-a ] [-r ] [-t ] [-A ] [-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;iradius = 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 + +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 diff --git a/pjhead.c b/pjhead.c new file mode 100644 index 0000000..484f73d --- /dev/null +++ b/pjhead.c @@ -0,0 +1,16 @@ +/* + * This code was entirely written by Nathan Wagner + * and is in the public domain. + */ + +#include +#include +#include +#include + +/* + * Proj 4 provides its own entry points into + * the code, so none of the library functions + * need to be global + */ +#define ISEA_STATIC static diff --git a/pjtail.c b/pjtail.c new file mode 100644 index 0000000..19123d4 --- /dev/null +++ b/pjtail.c @@ -0,0 +1,105 @@ +/* + * Proj 4 integration code follows + */ + +#define PROJ_PARMS__ \ + struct isea_dgg dgg; + +#define PJ_LIB__ +#include + +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) diff --git a/proj.c b/proj.c new file mode 100644 index 0000000..d69cffe --- /dev/null +++ b/proj.c @@ -0,0 +1,10 @@ +#include +#include +#include +#include + +#define ISEA_STATIC static +#include "hfuncs.c" +#include "pfuncs.c" + + diff --git a/proj4.patch b/proj4.patch new file mode 100644 index 0000000..7fa41bb --- /dev/null +++ b/proj4.patch @@ -0,0 +1,69 @@ +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") diff --git a/t/cities.c b/t/cities.c new file mode 100644 index 0000000..b8e3d3e --- /dev/null +++ b/t/cities.c @@ -0,0 +1,110 @@ +#include +#include +#include +#include + +#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); +} diff --git a/t/cities.in b/t/cities.in new file mode 100644 index 0000000..4b8c67c --- /dev/null +++ b/t/cities.in @@ -0,0 +1,16 @@ +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 + diff --git a/t/cities.pt b/t/cities.pt new file mode 100644 index 0000000..6e214c6 --- /dev/null +++ b/t/cities.pt @@ -0,0 +1,18 @@ +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} +}; diff --git a/t/globehex.c b/t/globehex.c new file mode 100644 index 0000000..4da098f --- /dev/null +++ b/t/globehex.c @@ -0,0 +1,111 @@ +#include +#include +#include +#include + +#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 +#include + +#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); + } + } + +} + diff --git a/t/hexes.pt b/t/hexes.pt new file mode 100644 index 0000000..9c04aec --- /dev/null +++ b/t/hexes.pt @@ -0,0 +1,23 @@ +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} +}; diff --git a/t/hexes2.dgg b/t/hexes2.dgg new file mode 100644 index 0000000..2ce4446 --- /dev/null +++ b/t/hexes2.dgg @@ -0,0 +1,152 @@ +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} +}; diff --git a/t/hexes2.dgg.plane b/t/hexes2.dgg.plane new file mode 100644 index 0000000..35a4e65 --- /dev/null +++ b/t/hexes2.dgg.plane @@ -0,0 +1,146 @@ +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 diff --git a/t/hexes2.pt b/t/hexes2.pt new file mode 100644 index 0000000..aa764bb --- /dev/null +++ b/t/hexes2.pt @@ -0,0 +1,152 @@ +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} +}; diff --git a/t/hexes2.txt b/t/hexes2.txt new file mode 100644 index 0000000..8b40c1c --- /dev/null +++ b/t/hexes2.txt @@ -0,0 +1,146 @@ +-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 diff --git a/t/hextest.c b/t/hextest.c new file mode 100644 index 0000000..bc8249b --- /dev/null +++ b/t/hextest.c @@ -0,0 +1,38 @@ +#include +#include + +#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(); + +} + diff --git a/t/mkc b/t/mkc new file mode 100755 index 0000000..2d08162 --- /dev/null +++ b/t/mkc @@ -0,0 +1,17 @@ +#!/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"; diff --git a/t/plane.c b/t/plane.c new file mode 100644 index 0000000..66306f9 --- /dev/null +++ b/t/plane.c @@ -0,0 +1,62 @@ +#include +#include +#include +#include + +#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 +#include +#include +#include + +#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 +#include +#include +#include + +#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 +#include +#include +#include + +#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 +#include +#include +#include + +#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 +#include +#include +#include + +#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 +#include +#include +#include + +#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 +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; +} diff --git a/t/tap.h b/t/tap.h new file mode 100644 index 0000000..e64f9d5 --- /dev/null +++ b/t/tap.h @@ -0,0 +1,89 @@ +/*- + * 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); diff --git a/t/testing.h b/t/testing.h new file mode 100644 index 0000000..d1c007e --- /dev/null +++ b/t/testing.h @@ -0,0 +1,4 @@ +struct points { + double lon, lat; + char *data; +}; diff --git a/transform.c b/transform.c new file mode 100644 index 0000000..087eb42 --- /dev/null +++ b/transform.c @@ -0,0 +1,228 @@ +#include +#include +#include +#include + +#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 ] [-z ] [-p ] [-I] [-D] [-P] [-a ] [-r ] [-t ] [-A ] [-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;iradius = 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); +} -- 2.40.0