aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
Diffstat (limited to 'src')
-rw-r--r--src/PJ_healpix.c206
-rw-r--r--src/projects.h280
2 files changed, 250 insertions, 236 deletions
diff --git a/src/PJ_healpix.c b/src/PJ_healpix.c
index e8d3f38d..e1b48daa 100644
--- a/src/PJ_healpix.c
+++ b/src/PJ_healpix.c
@@ -56,8 +56,8 @@ struct pj_opaque {
};
typedef struct {
- int cn; /* An integer 0--3 indicating the position of the polar cap. */
- double x, y; /* Coordinates of the pole point (point of most extreme latitude on the polar caps). */
+ int cn; /* An integer 0--3 indicating the position of the polar cap. */
+ double x, y; /* Coordinates of the pole point (point of most extreme latitude on the polar caps). */
enum Region {north, south, equatorial} region;
} CapMap;
@@ -160,41 +160,60 @@ static int pnpoly(int nvert, double vert[][2], double testx, double testy) {
int in_image(double x, double y, int proj, int north_square, int south_square) {
if (proj == 0) {
double healpixVertsJit[][2] = {
- {-1.0*M_PI- EPS, M_PI/4.0},
- {-3.0*M_PI/4.0, M_PI/2.0 + EPS},
- {-1.0*M_PI/2.0, M_PI/4.0 + EPS},
- {-1.0*M_PI/4.0, M_PI/2.0 + EPS},
- {0.0, M_PI/4.0 + EPS},
- {M_PI/4.0, M_PI/2.0 + EPS},
- {M_PI/2.0, M_PI/4.0 + EPS},
- {3.0*M_PI/4.0, M_PI/2.0 + EPS},
- {M_PI+ EPS, M_PI/4.0},
- {M_PI+ EPS, -1.0*M_PI/4.0},
- {3.0*M_PI/4.0, -1.0*M_PI/2.0 - EPS},
- {M_PI/2.0, -1.0*M_PI/4.0 - EPS},
- {M_PI/4.0, -1.0*M_PI/2.0 - EPS},
- {0.0, -1.0*M_PI/4.0 - EPS},
- {-1.0*M_PI/4.0, -1.0*M_PI/2.0 - EPS},
- {-1.0*M_PI/2.0, -1.0*M_PI/4.0 - EPS},
- {-3.0*M_PI/4.0, -1.0*M_PI/2.0 - EPS},
- {-1.0*M_PI - EPS, -1.0*M_PI/4.0}
+ {-1.0*M_PI - EPS, M_FORTPI},
+ {-3.0*M_FORTPI, M_HALFPI + EPS},
+ {M_HALFPI, M_FORTPI + EPS},
+ {M_FORTPI, M_HALFPI + EPS},
+ {0.0, M_FORTPI + EPS},
+ {M_FORTPI, M_HALFPI + EPS},
+ {M_HALFPI, M_FORTPI + EPS},
+ {3.0*M_FORTPI, M_HALFPI + EPS},
+ {M_PI + EPS, M_FORTPI},
+ {M_PI + EPS, -M_FORTPI},
+ {3.0*M_FORTPI, -1.0*M_HALFPI - EPS},
+ {M_HALFPI, -1.0*M_FORTPI - EPS},
+ {M_FORTPI, -1.0*M_HALFPI - EPS},
+ {0.0, -1.0*M_FORTPI - EPS},
+ {-1.0*M_FORTPI, -1.0*M_HALFPI - EPS},
+ {-1.0*M_HALFPI, -1.0*M_FORTPI - EPS},
+ {-3.0*M_FORTPI, -1.0*M_HALFPI - EPS},
+ {-1.0*M_PI - EPS, -1.0*M_FORTPI}
};
return pnpoly((int)sizeof(healpixVertsJit)/
sizeof(healpixVertsJit[0]), healpixVertsJit, x, y);
} else {
- double rhealpixVertsJit[][2] = {
- {-1.0*M_PI - EPS, M_PI/4.0 + EPS},
- {-1.0*M_PI + north_square*M_PI/2.0- EPS, M_PI/4.0 + EPS},
- {-1.0*M_PI + north_square*M_PI/2.0- EPS, 3*M_PI/4.0 + EPS},
- {-1.0*M_PI + (north_square + 1.0)*M_PI/2.0 + EPS, 3*M_PI/4.0 + EPS},
- {-1.0*M_PI + (north_square + 1.0)*M_PI/2.0 + EPS, M_PI/4.0 + EPS},
- {M_PI + EPS, M_PI/4.0 + EPS},
- {M_PI + EPS, -1.0*M_PI/4.0 - EPS},
- {-1.0*M_PI + (south_square + 1.0)*M_PI/2.0 + EPS, -1.0*M_PI/4.0 - EPS},
- {-1.0*M_PI + (south_square + 1.0)*M_PI/2.0 + EPS, -3.0*M_PI/4.0 - EPS},
- {-1.0*M_PI + south_square*M_PI/2.0 - EPS, -3.0*M_PI/4.0 - EPS},
- {-1.0*M_PI + south_square*M_PI/2.0 - EPS, -1.0*M_PI/4.0 - EPS},
- {-1.0*M_PI - EPS, -1.0*M_PI/4.0 - EPS}};
+ /**
+ * Assigning each element by index to avoid warnings such as
+ * 'initializer element is not computable at load time'.
+ * Before C99 this was not allowed and to keep as portable as
+ * possible we do it the C89 way here.
+ **/
+ double rhealpixVertsJit[12][2];
+ rhealpixVertsJit[0][0] = -1.0*M_PI - EPS;
+ rhealpixVertsJit[0][1] = M_FORTPI + EPS;
+ rhealpixVertsJit[1][0] = -1.0*M_PI + north_square*M_HALFPI- EPS;
+ rhealpixVertsJit[1][2] = M_FORTPI + EPS;
+ rhealpixVertsJit[2][0] = -1.0*M_PI + north_square*M_HALFPI- EPS;
+ rhealpixVertsJit[2][1] = 3*M_FORTPI + EPS;
+ rhealpixVertsJit[3][0] = -1.0*M_PI + (north_square + 1.0)*M_HALFPI + EPS;
+ rhealpixVertsJit[3][1] = 3*M_FORTPI + EPS;
+ rhealpixVertsJit[4][0] = -1.0*M_PI + (north_square + 1.0)*M_HALFPI + EPS;
+ rhealpixVertsJit[4][1] = M_FORTPI + EPS;
+ rhealpixVertsJit[5][0] = M_PI + EPS;
+ rhealpixVertsJit[5][1] = M_FORTPI + EPS;
+ rhealpixVertsJit[6][0] = M_PI + EPS;
+ rhealpixVertsJit[6][1] = -1.0*M_FORTPI - EPS;
+ rhealpixVertsJit[7][0] = -1.0*M_PI + (south_square + 1.0)*M_HALFPI + EPS;
+ rhealpixVertsJit[7][1] = -1.0*M_FORTPI - EPS;
+ rhealpixVertsJit[8][0] = -1.0*M_PI + (south_square + 1.0)*M_HALFPI + EPS;
+ rhealpixVertsJit[8][1] = -3.0*M_FORTPI - EPS;
+ rhealpixVertsJit[9][0] = -1.0*M_PI + south_square*M_HALFPI - EPS;
+ rhealpixVertsJit[9][1] = -3.0*M_FORTPI - EPS;
+ rhealpixVertsJit[10][0] = -1.0*M_PI + south_square*M_HALFPI - EPS;
+ rhealpixVertsJit[10][1] = -1.0*M_FORTPI - EPS;
+ rhealpixVertsJit[11][0] = -1.0*M_PI - EPS;
+ rhealpixVertsJit[11][1] = -1.0*M_FORTPI - EPS;
+
return pnpoly((int)sizeof(rhealpixVertsJit)/
sizeof(rhealpixVertsJit[0]), rhealpixVertsJit, x, y);
}
@@ -247,9 +266,9 @@ XY healpix_sphere(LP lp) {
if (cn >= 4) {
cn = 3;
}
- lamc = -3*M_PI/4 + (M_PI/2)*cn;
+ lamc = -3*M_FORTPI + (M_HALFPI)*cn;
xy.x = lamc + (lam - lamc)*sigma;
- xy.y = pj_sign(phi)*M_PI/4*(2 - sigma);
+ xy.y = pj_sign(phi)*M_FORTPI*(2 - sigma);
}
return xy;
}
@@ -262,25 +281,25 @@ LP healpix_sphere_inverse(XY xy) {
LP lp;
double x = xy.x;
double y = xy.y;
- double y0 = M_PI/4.0;
+ double y0 = M_FORTPI;
/* Equatorial region. */
if (fabsl(y) <= y0) {
lp.lam = x;
lp.phi = asin(8.0*y/(3.0*M_PI));
- } else if (fabsl(y) < M_PI/2.0) {
+ } else if (fabsl(y) < M_HALFPI) {
double cn = floor(2.0*x/M_PI + 2.0);
double xc, tau;
if (cn >= 4) {
cn = 3;
}
- xc = -3.0*M_PI/4.0 + (M_PI/2.0)*cn;
+ xc = -3.0*M_FORTPI + (M_HALFPI)*cn;
tau = 2.0 - 4.0*fabsl(y)/M_PI;
lp.lam = xc + (x - xc)/tau;
lp.phi = pj_sign(y)*asin(1.0 - pow(tau , 2.0)/3.0);
} else {
lp.lam = -1.0*M_PI;
- lp.phi = pj_sign(y)*M_PI/2.0;
+ lp.phi = pj_sign(y)*M_HALFPI;
}
return (lp);
}
@@ -343,48 +362,48 @@ static CapMap get_cap(double x, double y, int north_square, int south_square,
capmap.x = x;
capmap.y = y;
if (inverse == 0) {
- if (y > M_PI/4.0) {
+ if (y > M_FORTPI) {
capmap.region = north;
- c = M_PI/2.0;
- } else if (y < -1*M_PI/4.0) {
+ c = M_HALFPI;
+ } else if (y < -1*M_FORTPI) {
capmap.region = south;
- c = -1*M_PI/2.0;
+ c = -1*M_HALFPI;
} else {
capmap.region = equatorial;
capmap.cn = 0;
return capmap;
}
/* polar region */
- if (x < -1*M_PI/2.0) {
+ if (x < -1*M_HALFPI) {
capmap.cn = 0;
- capmap.x = (-1*3.0*M_PI/4.0);
+ capmap.x = (-1*3.0*M_FORTPI);
capmap.y = c;
- } else if (x >= -1*M_PI/2.0 && x < 0) {
+ } else if (x >= -1*M_HALFPI && x < 0) {
capmap.cn = 1;
- capmap.x = -1*M_PI/4.0;
+ capmap.x = -1*M_FORTPI;
capmap.y = c;
- } else if (x >= 0 && x < M_PI/2.0) {
+ } else if (x >= 0 && x < M_HALFPI) {
capmap.cn = 2;
- capmap.x = M_PI/4.0;
+ capmap.x = M_FORTPI;
capmap.y = c;
} else {
capmap.cn = 3;
- capmap.x = 3.0*M_PI/4.0;
+ capmap.x = 3.0*M_FORTPI;
capmap.y = c;
}
return capmap;
} else {
double eps;
- if (y > M_PI/4.0) {
+ if (y > M_FORTPI) {
capmap.region = north;
- capmap.x = (-3.0*M_PI/4.0 + north_square*M_PI/2.0);
- capmap.y = M_PI/2.0;
- x = x - north_square*M_PI/2.0;
- } else if (y < -1*M_PI/4.0) {
+ capmap.x = (-3.0*M_FORTPI + north_square*M_HALFPI);
+ capmap.y = M_HALFPI;
+ x = x - north_square*M_HALFPI;
+ } else if (y < -1*M_FORTPI) {
capmap.region = south;
- capmap.x = (-3.0*M_PI/4.0 + south_square*M_PI/2);
- capmap.y = -1*M_PI/2.0;
- x = x - south_square*M_PI/2.0;
+ capmap.x = (-3.0*M_FORTPI + south_square*M_HALFPI);
+ capmap.y = -1*M_HALFPI;
+ x = x - south_square*M_HALFPI;
} else {
capmap.region = equatorial;
capmap.cn = 0;
@@ -394,21 +413,21 @@ static CapMap get_cap(double x, double y, int north_square, int south_square,
x, y moves to when rHEALPix polar square is disassembled. */
eps = 1e-15; /* Kludge. Fuzz to avoid some rounding errors. */
if (capmap.region == north) {
- if (y >= -1*x - M_PI/4.0 - eps && y < x + 5.0*M_PI/4.0 - eps) {
+ if (y >= -1*x - M_FORTPI - eps && y < x + 5.0*M_FORTPI - eps) {
capmap.cn = (north_square + 1) % 4;
- } else if (y > -1*x -1*M_PI/4.0 + eps && y >= x + 5.0*M_PI/4.0 - eps) {
+ } else if (y > -1*x -1*M_FORTPI + eps && y >= x + 5.0*M_FORTPI - eps) {
capmap.cn = (north_square + 2) % 4;
- } else if (y <= -1*x -1*M_PI/4.0 + eps && y > x + 5.0*M_PI/4.0 + eps) {
+ } else if (y <= -1*x -1*M_FORTPI + eps && y > x + 5.0*M_FORTPI + eps) {
capmap.cn = (north_square + 3) % 4;
} else {
capmap.cn = north_square;
}
} else if (capmap.region == south) {
- if (y <= x + M_PI/4.0 + eps && y > -1*x - 5.0*M_PI/4 + eps) {
+ if (y <= x + M_FORTPI + eps && y > -1*x - 5.0*M_FORTPI + eps) {
capmap.cn = (south_square + 1) % 4;
- } else if (y < x + M_PI/4.0 - eps && y <= -1*x - 5.0*M_PI/4.0 + eps) {
+ } else if (y < x + M_FORTPI - eps && y <= -1*x - 5.0*M_FORTPI + eps) {
capmap.cn = (south_square + 2) % 4;
- } else if (y >= x + M_PI/4.0 - eps && y < -1*x - 5.0*M_PI/4.0 - eps) {
+ } else if (y >= x + M_FORTPI - eps && y < -1*x - 5.0*M_FORTPI - eps) {
capmap.cn = (south_square + 3) % 4;
} else {
capmap.cn = south_square;
@@ -433,9 +452,12 @@ static XY combine_caps(double x, double y, int north_square, int south_square,
XY xy;
double v[2];
double a[2];
+ double c[2];
double vector[2];
double v_min_c[2];
double ret_dot[2];
+ double (*tmpRot)[2];
+ int pole = 0;
CapMap capmap = get_cap(x, y, north_square, south_square, inverse);
if (capmap.region == equatorial) {
@@ -443,62 +465,44 @@ static XY combine_caps(double x, double y, int north_square, int south_square,
xy.y = capmap.y;
return xy;
}
- v[0] = x;
- v[1] = y;
+
+ v[0] = x; v[1] = y;
+ c[0] = capmap.x; c[1] = capmap.y;
+
if (inverse == 0) {
/* Rotate (x, y) about its polar cap tip and then translate it to
north_square or south_square. */
- int pole = 0;
- double (*tmpRot)[2];
- double c[2] = {capmap.x, capmap.y};
+ a[0] = -3.0*M_FORTPI + pole*M_HALFPI;
+ a[1] = M_HALFPI;
if (capmap.region == north) {
pole = north_square;
- a[0] = (-3.0*M_PI/4.0 + pole*M_PI/2);
- a[1] = (M_PI/2.0 + pole*0);
tmpRot = rot[get_rotate_index(capmap.cn - pole)];
- vector_sub(v, c, v_min_c);
- dot_product(tmpRot, v_min_c, ret_dot);
- vector_add(ret_dot, a, vector);
} else {
pole = south_square;
- a[0] = (-3.0*M_PI/4.0 + pole*M_PI/2);
- a[1] = (M_PI/-2.0 + pole*0);
tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
- vector_sub(v, c, v_min_c);
- dot_product(tmpRot, v_min_c, ret_dot);
- vector_add(ret_dot, a, vector);
}
- xy.x = vector[0];
- xy.y = vector[1];
- return xy;
} else {
/* Inverse function.
Unrotate (x, y) and then translate it back. */
- int pole = 0;
- double (*tmpRot)[2];
- double c[2] = {capmap.x, capmap.y};
+ a[0] = -3.0*M_FORTPI + capmap.cn*M_HALFPI;
+ a[1] = M_HALFPI;
/* disassemble */
if (capmap.region == north) {
pole = north_square;
- a[0] = (-3.0*M_PI/4.0 + capmap.cn*M_PI/2);
- a[1] = (M_PI/2.0 + capmap.cn*0);
tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
- vector_sub(v, c, v_min_c);
- dot_product(tmpRot, v_min_c, ret_dot);
- vector_add(ret_dot, a, vector);
} else {
pole = south_square;
- a[0] = (-3.0*M_PI/4.0 + capmap.cn*M_PI/2);
- a[1] = (M_PI/-2.0 + capmap.cn*0);
tmpRot = rot[get_rotate_index(capmap.cn - pole)];
- vector_sub(v, c, v_min_c);
- dot_product(tmpRot, v_min_c, ret_dot);
- vector_add(ret_dot, a, vector);
}
- xy.x = vector[0];
- xy.y = vector[1];
- return xy;
}
+
+ vector_sub(v, c, v_min_c);
+ dot_product(tmpRot, v_min_c, ret_dot);
+ vector_add(ret_dot, a, vector);
+
+ xy.x = vector[0];
+ xy.y = vector[1];
+ return xy;
}
@@ -621,9 +625,9 @@ PJ *PROJECTION(healpix) {
P->opaque = Q;
if (P->es) {
- Q->apa = pj_authset(P->es); /* For auth_lat(). */
- Q->qp = pj_qsfn(1.0, P->e, P->one_es); /* For auth_lat(). */
- P->a = P->a*sqrt(0.5*Q->qp); /* Set P->a to authalic radius. */
+ Q->apa = pj_authset(P->es); /* For auth_lat(). */
+ Q->qp = pj_qsfn(1.0, P->e, P->one_es); /* For auth_lat(). */
+ P->a = P->a*sqrt(0.5*Q->qp); /* Set P->a to authalic radius. */
P->ra = 1.0/P->a;
P->fwd = e_healpix_forward;
P->inv = e_healpix_inverse;
diff --git a/src/projects.h b/src/projects.h
index 9aca7d95..ced5a4f0 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -54,15 +54,15 @@ extern "C" {
#endif
#ifndef NULL
-# define NULL 0
+# define NULL 0
#endif
#ifndef FALSE
-# define FALSE 0
+# define FALSE 0
#endif
#ifndef TRUE
-# define TRUE 1
+# define TRUE 1
#endif
#ifndef MAX
@@ -74,11 +74,12 @@ extern "C" {
# define ABS(x) ((x<0) ? (-1*(x)) : x)
#endif
- /* maximum path/filename */
+/* maximum path/filename */
#ifndef MAX_PATH_FILENAME
#define MAX_PATH_FILENAME 1024
#endif
- /* prototype hypot for systems where absent */
+
+/* prototype hypot for systems where absent */
#ifndef _WIN32
extern double hypot(double, double);
#endif
@@ -97,15 +98,23 @@ extern double hypot(double, double);
#define _USE_MATH_DEFINES
#endif
+/* If we still haven't got M_PI*, we rely on our own defines.
+ * For example, this is necessary when compiling with gcc and
+ * the -ansi flag.
+ */
+#ifndef M_PI
+#define M_PI 3.14159265358979310
+#define M_PI_2 1.57079632679489660
+#define M_PI_4 0.78539816339744828
+#endif
+
/* some more useful math constants and aliases */
-#define M_FORTPI M_PI_4 /* pi/4 */
-#define M_HALFPI M_PI_2 /* pi/2 */
-/* M_PI pi */
-#define M_PI_HALFPI 4.71238898038468985769 /* 1.5*pi */
-#define M_TWOPI 6.28318530717958647693 /* 2*pi */
-#define M_TWO_D_PI M_2_PI /* 2/pi */
-#define M_TWOPI_HALFPI 7.85398163397448309616 /* 2.5*pi */
-/* M_SQRT2 sqrt(2) */
+#define M_FORTPI M_PI_4 /* pi/4 */
+#define M_HALFPI M_PI_2 /* pi/2 */
+#define M_PI_HALFPI 4.71238898038468985769 /* 1.5*pi */
+#define M_TWOPI 6.28318530717958647693 /* 2*pi */
+#define M_TWO_D_PI M_2_PI /* 2/pi */
+#define M_TWOPI_HALFPI 7.85398163397448309616 /* 2.5*pi */
/* maximum tag id length for +init and default files */
@@ -133,7 +142,7 @@ struct projFileAPI_t;
/* proj thread context */
typedef struct {
- int last_errno;
+ int last_errno;
int debug_level;
void (*logger)(void *, int, const char *);
void *app_data;
@@ -156,7 +165,7 @@ typedef struct {
#define USE_PROJUV
typedef struct { double u, v; } projUV;
-typedef struct { double r, i; } COMPLEX;
+typedef struct { double r, i; } COMPLEX;
typedef struct { double u, v, w; } projUVW;
#ifndef PJ_LIB__
@@ -175,39 +184,39 @@ typedef union { double f; int i; char *s; } PROJVALUE;
struct PJconsts;
struct PJ_LIST {
- char *id; /* projection keyword */
- struct PJconsts *(*proj)(struct PJconsts*);/* projection entry point */
- char * const *descr; /* description text */
+ char *id; /* projection keyword */
+ struct PJconsts *(*proj)(struct PJconsts*); /* projection entry point */
+ char * const *descr; /* description text */
};
/* Merging this into the PJ_LIST infrastructure is tempting, but may imply ABI breakage. Perhaps at next major version? */
struct PJ_SELFTEST_LIST {
- char *id; /* projection keyword */
- int (* testfunc)(void); /* projection entry point */
+ char *id; /* projection keyword */
+ int (* testfunc)(void); /* projection entry point */
};
struct PJ_ELLPS {
- char *id; /* ellipse keyword name */
- char *major; /* a= value */
- char *ell; /* elliptical parameter */
- char *name; /* comments */
+ char *id; /* ellipse keyword name */
+ char *major; /* a= value */
+ char *ell; /* elliptical parameter */
+ char *name; /* comments */
};
struct PJ_UNITS {
- char *id; /* units keyword */
- char *to_meter; /* multiply by value to get meters */
- char *name; /* comments */
+ char *id; /* units keyword */
+ char *to_meter; /* multiply by value to get meters */
+ char *name; /* comments */
};
struct PJ_DATUMS {
- char *id; /* datum keyword */
- char *defn; /* ie. "to_wgs84=..." */
- char *ellipse_id; /* ie from ellipse table */
- char *comments; /* EPSG code, etc */
+ char *id; /* datum keyword */
+ char *defn; /* ie. "to_wgs84=..." */
+ char *ellipse_id;/* ie from ellipse table */
+ char *comments; /* EPSG code, etc */
};
struct PJ_PRIME_MERIDIANS {
- char *id; /* prime meridian keyword */
- char *defn; /* offset from greenwich in DMS format. */
+ char *id; /* prime meridian keyword */
+ char *defn; /* offset from greenwich in DMS format. */
};
typedef struct {
@@ -218,31 +227,32 @@ typedef struct {
} PJ_Region;
struct DERIVS {
- double x_l, x_p; /* derivatives of x for lambda-phi */
- double y_l, y_p; /* derivatives of y for lambda-phi */
+ double x_l, x_p; /* derivatives of x for lambda-phi */
+ double y_l, y_p; /* derivatives of y for lambda-phi */
};
struct FACTORS {
- struct DERIVS der;
- double h, k; /* meridinal, parallel scales */
- double omega, thetap; /* angular distortion, theta prime */
- double conv; /* convergence */
- double s; /* areal scale factor */
- double a, b; /* max-min scale error */
- int code; /* info as to analytics, see following */
+ struct DERIVS der;
+ double h, k; /* meridinal, parallel scales */
+ double omega, thetap; /* angular distortion, theta prime */
+ double conv; /* convergence */
+ double s; /* areal scale factor */
+ double a, b; /* max-min scale error */
+ int code; /* info as to analytics, see following */
};
-#define IS_ANAL_XL_YL 01 /* derivatives of lon analytic */
-#define IS_ANAL_XP_YP 02 /* derivatives of lat analytic */
-#define IS_ANAL_HK 04 /* h and k analytic */
-#define IS_ANAL_CONV 010 /* convergence analytic */
- /* parameter list struct */
-typedef struct ARG_list {
- struct ARG_list *next;
- char used;
- char param[1]; } paralist;
- /* base projection data structure */
+#define IS_ANAL_XL_YL 01 /* derivatives of lon analytic */
+#define IS_ANAL_XP_YP 02 /* derivatives of lat analytic */
+#define IS_ANAL_HK 04 /* h and k analytic */
+#define IS_ANAL_CONV 010 /* convergence analytic */
+
+/* parameter list struct */
+typedef struct ARG_list {
+ struct ARG_list *next;
+ char used;
+ char param[1]; } paralist;
+/* base projection data structure */
#ifdef PJ_LIB__
/* we need this forward declaration in order to be able to add a
pointer to struct opaque to the typedef struct PJconsts below */
@@ -251,60 +261,60 @@ typedef struct ARG_list {
typedef struct PJconsts {
projCtx_t *ctx;
- XY (*fwd)(LP, struct PJconsts *);
- LP (*inv)(XY, struct PJconsts *);
- XYZ (*fwd3d)(LPZ, struct PJconsts *);
- LPZ (*inv3d)(XYZ, struct PJconsts *);
- void (*spc)(LP, struct PJconsts *, struct FACTORS *);
- void (*pfree)(struct PJconsts *);
- const char *descr;
- paralist *params; /* parameter list */
- int over; /* over-range flag */
- int geoc; /* geocentric latitude flag */
- int is_latlong; /* proj=latlong ... not really a projection at all */
- int is_geocent; /* proj=geocent ... not really a projection at all */
- double
- a, /* major axis or radius if es==0 */
- a_orig, /* major axis before any +proj related adjustment */
- es, /* e ^ 2 */
- es_orig, /* es before any +proj related adjustment */
- e, /* eccentricity */
- ra, /* 1/A */
- one_es, /* 1 - e^2 */
- rone_es, /* 1/one_es */
- lam0, phi0, /* central longitude, latitude */
- x0, y0, /* easting and northing */
- k0, /* general scaling factor */
- to_meter, fr_meter; /* cartesian scaling */
-
- int datum_type; /* PJD_UNKNOWN/3PARAM/7PARAM/GRIDSHIFT/WGS84 */
- double datum_params[7];
- struct _pj_gi **gridlist;
- int gridlist_count;
-
- int has_geoid_vgrids;
- struct _pj_gi **vgridlist_geoid;
- int vgridlist_geoid_count;
- double vto_meter, vfr_meter;
-
- double from_greenwich; /* prime meridian offset (in radians) */
- double long_wrap_center; /* 0.0 for -180 to 180, actually in radians*/
- int is_long_wrap_set;
- char axis[4];
-
- /* New Datum Shift Grid Catalogs */
- char *catalog_name;
- struct _PJ_GridCatalog *catalog;
-
- double datum_date;
-
- struct _pj_gi *last_before_grid;
- PJ_Region last_before_region;
- double last_before_date;
-
- struct _pj_gi *last_after_grid;
- PJ_Region last_after_region;
- double last_after_date;
+ XY (*fwd)(LP, struct PJconsts *);
+ LP (*inv)(XY, struct PJconsts *);
+ XYZ (*fwd3d)(LPZ, struct PJconsts *);
+ LPZ (*inv3d)(XYZ, struct PJconsts *);
+ void (*spc)(LP, struct PJconsts *, struct FACTORS *);
+ void (*pfree)(struct PJconsts *);
+
+ const char *descr;
+ paralist *params; /* parameter list */
+ int over; /* over-range flag */
+ int geoc; /* geocentric latitude flag */
+ int is_latlong; /* proj=latlong ... not really a projection at all */
+ int is_geocent; /* proj=geocent ... not really a projection at all */
+ double a; /* major axis or radius if es==0 */
+ double a_orig; /* major axis before any +proj related adjustment */
+ double es; /* e ^ 2 */
+ double es_orig; /* es before any +proj related adjustment */
+ double e; /* eccentricity */
+ double ra; /* 1/A */
+ double one_es; /* 1 - e^2 */
+ double rone_es; /* 1/one_es */
+ double lam0, phi0; /* central longitude, latitude */
+ double x0, y0; /* easting and northing */
+ double k0; /* general scaling factor */
+ double to_meter, fr_meter; /* cartesian scaling */
+
+ int datum_type; /* PJD_UNKNOWN/3PARAM/7PARAM/GRIDSHIFT/WGS84 */
+ double datum_params[7];
+ struct _pj_gi **gridlist;
+ int gridlist_count;
+
+ int has_geoid_vgrids;
+ struct _pj_gi **vgridlist_geoid;
+ int vgridlist_geoid_count;
+ double vto_meter, vfr_meter;
+
+ double from_greenwich; /* prime meridian offset (in radians) */
+ double long_wrap_center; /* 0.0 for -180 to 180, actually in radians*/
+ int is_long_wrap_set;
+ char axis[4];
+
+ /* New Datum Shift Grid Catalogs */
+ char *catalog_name;
+ struct _PJ_GridCatalog *catalog;
+
+ double datum_date;
+
+ struct _pj_gi *last_before_grid;
+ PJ_Region last_before_region;
+ double last_before_date;
+
+ struct _pj_gi *last_after_grid;
+ PJ_Region last_after_region;
+ double last_after_date;
#ifdef PJ_LIB__
struct pj_opaque *opaque;
@@ -342,16 +352,16 @@ extern struct PJ_PRIME_MERIDIANS pj_prime_meridians[];
#endif
#ifdef PJ_LIB__
- /* repetitive projection code */
+/* repetitive projection code */
#define PROJ_HEAD(id, name) static const char des_##id [] = name
#define ENTRYA(name) \
- C_NAMESPACE_VAR const char * const pj_s_##name = des_##name; \
- C_NAMESPACE PJ *pj_##name(PJ *P) { if (!P) { \
- if( (P = (PJ*) pj_malloc(sizeof(PJ))) != NULL) { \
+ C_NAMESPACE_VAR const char * const pj_s_##name = des_##name; \
+ C_NAMESPACE PJ *pj_##name(PJ *P) { if (!P) { \
+ if( (P = (PJ*) pj_malloc(sizeof(PJ))) != NULL) { \
memset( P, 0, sizeof(PJ) ); \
- P->pfree = freeup; P->fwd = 0; P->inv = 0; \
+ P->pfree = freeup; P->fwd = 0; P->inv = 0; \
P->fwd3d = 0; P->inv3d = 0; \
- P->spc = 0; P->descr = des_##name;
+ P->spc = 0; P->descr = des_##name;
#define ENTRYX } return P; } else {
#define ENTRY0(name) ENTRYA(name) ENTRYX
#define ENTRY1(name, a) ENTRYA(name) P->a = 0; ENTRYX
@@ -413,22 +423,22 @@ typedef struct { float lam, phi; } FLP;
typedef struct { int lam, phi; } ILP;
struct CTABLE {
- char id[MAX_TAB_ID]; /* ascii info */
- LP ll; /* lower left corner coordinates */
- LP del; /* size of cells */
- ILP lim; /* limits of conversion matrix */
- FLP *cvs; /* conversion matrix */
+ char id[MAX_TAB_ID]; /* ascii info */
+ LP ll; /* lower left corner coordinates */
+ LP del; /* size of cells */
+ ILP lim; /* limits of conversion matrix */
+ FLP *cvs; /* conversion matrix */
};
typedef struct _pj_gi {
- char *gridname; /* identifying name of grid, eg "conus" or ntv2_0.gsb */
- char *filename; /* full path to filename */
+ char *gridname; /* identifying name of grid, eg "conus" or ntv2_0.gsb */
+ char *filename; /* full path to filename */
const char *format; /* format of this grid, ie "ctable", "ntv1",
"ntv2" or "missing". */
- int grid_offset; /* offset in file, for delayed loading */
- int must_swap; /* only for NTv2 */
+ int grid_offset; /* offset in file, for delayed loading */
+ int must_swap; /* only for NTv2 */
struct CTABLE *ct;
@@ -438,18 +448,18 @@ typedef struct _pj_gi {
typedef struct {
PJ_Region region;
- int priority; /* higher used before lower */
- double date; /* year.fraction */
- char *definition; /* usually the gridname */
+ int priority; /* higher used before lower */
+ double date; /* year.fraction */
+ char *definition; /* usually the gridname */
PJ_GRIDINFO *gridinfo;
- int available; /* 0=unknown, 1=true, -1=false */
+ int available; /* 0=unknown, 1=true, -1=false */
} PJ_GridCatalogEntry;
typedef struct _PJ_GridCatalog {
char *catalog_name;
- PJ_Region region; /* maximum extent of catalog data */
+ PJ_Region region; /* maximum extent of catalog data */
int entry_count;
PJ_GridCatalogEntry *entries;
@@ -493,18 +503,18 @@ COMPLEX pj_zpolyd1(COMPLEX, COMPLEX *, int, COMPLEX *);
int pj_deriv(LP, double, PJ *, struct DERIVS *);
int pj_factors(LP, PJ *, double, struct FACTORS *);
-struct PW_COEF {/* row coefficient structure */
- int m; /* number of c coefficients (=0 for none) */
- double *c; /* power coefficients */
+struct PW_COEF { /* row coefficient structure */
+ int m; /* number of c coefficients (=0 for none) */
+ double *c; /* power coefficients */
};
/* Approximation structures and procedures */
-typedef struct { /* Chebyshev or Power series structure */
- projUV a, b; /* power series range for evaluation */
- /* or Chebyshev argument shift/scaling */
- struct PW_COEF *cu, *cv;
- int mu, mv; /* maximum cu and cv index (+1 for count) */
- int power; /* != 0 if power series, else Chebyshev */
+typedef struct { /* Chebyshev or Power series structure */
+ projUV a, b; /* power series range for evaluation */
+ /* or Chebyshev argument shift/scaling */
+ struct PW_COEF *cu, *cv;
+ int mu, mv; /* maximum cu and cv index (+1 for count) */
+ int power; /* != 0 if power series, else Chebyshev */
} Tseries;
Tseries *mk_cheby(projUV, projUV, double, projUV *, projUV (*)(projUV), int, int, int);
projUV bpseval(projUV, Tseries *);
generated by cgit v1.2.3 (git 2.25.1) at 2026-01-01 21:18:59 +0000