From: Nathan Wagner Date: Fri, 27 May 2011 07:01:51 +0000 (+0000) Subject: Initial checkin. X-Git-Tag: 0.1^0 X-Git-Url: https://pd.if.org/git/?p=isea;a=commitdiff_plain Initial checkin. --- 7ed53098e43d4fb96c47345f18750fd588a04240 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); +}