Package: proj
Version: 4.4.9d-2
Severity: wishlist
Also sent to the [EMAIL PROTECTED]
Patch in attachment can be applied to proj-4.5.0. To apply the patch:
tar xzf proj-4.5.0.tgz
cd proj-4.5.0
patch -p1 <../proj-4.5.0-geod_api.patch
Following functions were added:
GEODESIC_T *GEOD_init(int, char **, GEODESIC_T *geodesic);
GEODESIC_T *GEOD_init_plus(const char *args, GEODESIC_T *geodesic);
They work the same way as proj_init and proj_init_plus, except that they
accept additional parameter GEODESIC_T *geodesic.
If geodesic is 0, memory for GEODESIC_T will be allocated and pointer to
allocated memory will be returned.
Functions geodetic_for, geod_inv, geod_pre are changed to accept pointer
to the GEODESIC_T. All global variables are moved to this structure.
Functions geod_for, geod_inv, geod_pre become members of the proj library.
I also include changes to Makefile.in, but may be it is better to
regenerate it again (or even exclude from distribution as it is
generated file).
Following 2 functions are not included into the patch, they just make
usage of the geod functions convinient.
void
geod_fwd(GEODESIC_T *GEODESIC, double lon1, double lat1, double S,
double az, double *lon2, double *lat2)
{
GEODESIC->p1.u = lat1;
GEODESIC->p1.v = lon1;
al12 = az;
GEODESIC->DIST = S; //* GEODESIC->FR_METER;
geod_pre(GEODESIC);
geod_for(GEODESIC);
*lat2 = GEODESIC->p2.u;
*lon2 = GEODESIC->p2.v;
}
int
geod_rev(GEODESIC_T *GEODESIC, double lon1, double lat1, double lon2,
double lat2, double *S, double *az)
{
GEODESIC->p1.u = lat1;
GEODESIC->p1.v = lon1;
GEODESIC->p2.u = lat2;
GEODESIC->p2.v = lon2;
int ret = geod_inv(GEODESIC);
if(al12 < 0)
{
al12 += TWOPI;
}
*az = al12;
*S = GEODESIC->DIST; // * GEODESIC->TO_METER;
return ret;
}
Hope it will be useful.
-- System Information:
Debian Release: testing/unstable
APT prefers testing
APT policy: (500, 'testing')
Architecture: i386 (i686)
Shell: /bin/sh linked to /bin/bash
Kernel: Linux 2.6.17-2-k7
Locale: LANG=en_US.UTF-8, LC_CTYPE=en_US.UTF-8 (charmap=UTF-8)
Versions of packages proj depends on:
ii libc6 2.3.6.ds1-7 GNU C Library: Shared libraries
proj recommends no packages.
-- no debconf information
diff -Naur proj-4.5.0-orig/src/geod.c proj-4.5.0/src/geod.c
--- proj-4.5.0-orig/src/geod.c 2006-12-04 11:18:34.000000000 +0100
+++ proj-4.5.0/src/geod.c 2006-12-05 12:35:35.000000000 +0100
@@ -23,7 +23,12 @@
pline[50], /* work string */
*usage =
"%s\nusage: %s [ -afFIptTwW [args] ] [ +opts[=arg] ] [ files ]\n";
- static void
+
+
+static GEODESIC_T Geodesic;
+static GEODESIC_T *GEODESIC = &Geodesic;
+
+static void
printLL(double p, double l) {
if (oform) {
(void)printf(oform, p * RAD_TO_DEG); TAB;
@@ -33,28 +38,31 @@
(void)fputs(rtodms(pline, l, 'E', 'W'),stdout);
}
}
- static void
+
+
+static void
do_arc(void) {
double az;
- printLL(phi2, lam2); putchar('\n');
- for (az = al12; n_alpha--; ) {
- al12 = az = adjlon(az + del_alpha);
- geod_pre();
- geod_for();
- printLL(phi2, lam2); putchar('\n');
+ printLL(GEODESIC->p2.u, GEODESIC->p2.v); putchar('\n');
+ for (az = GEODESIC->ALPHA12; GEODESIC->n_alpha--; ) {
+ GEODESIC->ALPHA12 = az = adjlon(az + GEODESIC->del_alpha);
+ geod_pre(GEODESIC);
+ geod_for(GEODESIC);
+ printLL(GEODESIC->p2.u, GEODESIC->p2.v); putchar('\n');
}
}
- static void /* generate intermediate geodesic coordinates */
+
+static void /* generate intermediate geodesic coordinates */
do_geod(void) {
double phil, laml, del_S;
- phil = phi2;
- laml = lam2;
- printLL(phi1, lam1); putchar('\n');
- for ( geod_S = del_S = geod_S / n_S; --n_S; geod_S += del_S) {
- geod_for();
- printLL(phi2, lam2); putchar('\n');
+ phil = GEODESIC->p2.u;
+ laml = GEODESIC->p2.v;
+ printLL(GEODESIC->p1.u, GEODESIC->p1.v); putchar('\n');
+ for ( GEODESIC->DIST = del_S = GEODESIC->DIST / GEODESIC->n_S; --(GEODESIC->n_S); GEODESIC->DIST += del_S) {
+ geod_for(GEODESIC);
+ printLL(GEODESIC->p2.u, GEODESIC->p2.v); putchar('\n');
}
printLL(phil, laml); putchar('\n');
}
@@ -76,51 +84,51 @@
fputs(line, stdout);
continue;
}
- phi1 = dmstor(s, &s);
- lam1 = dmstor(s, &s);
+ GEODESIC->p1.u = dmstor(s, &s);
+ GEODESIC->p1.v = dmstor(s, &s);
if (inverse) {
- phi2 = dmstor(s, &s);
- lam2 = dmstor(s, &s);
- geod_inv();
+ GEODESIC->p2.u = dmstor(s, &s);
+ GEODESIC->p2.v = dmstor(s, &s);
+ geod_inv(GEODESIC);
} else {
- al12 = dmstor(s, &s);
- geod_S = strtod(s, &s) * to_meter;
- geod_pre();
- geod_for();
+ GEODESIC->ALPHA12 = dmstor(s, &s);
+ GEODESIC->DIST = strtod(s, &s) * GEODESIC->TO_METER;
+ geod_pre(GEODESIC);
+ geod_for(GEODESIC);
}
if (!*s && (s > line)) --s; /* assumed we gobbled \n */
if (pos_azi) {
- if (al12 < 0.) al12 += TWOPI;
- if (al21 < 0.) al21 += TWOPI;
+ if (GEODESIC->ALPHA12 < 0.) GEODESIC->ALPHA12 += TWOPI;
+ if (GEODESIC->ALPHA21 < 0.) GEODESIC->ALPHA21 += TWOPI;
}
if (fullout) {
- printLL(phi1, lam1); TAB;
- printLL(phi2, lam2); TAB;
+ printLL(GEODESIC->p1.u, GEODESIC->p1.v); TAB;
+ printLL(GEODESIC->p2.u, GEODESIC->p2.v); TAB;
if (oform) {
- (void)printf(oform, al12 * RAD_TO_DEG); TAB;
- (void)printf(oform, al21 * RAD_TO_DEG); TAB;
- (void)printf(osform, geod_S * fr_meter);
+ (void)printf(oform, GEODESIC->ALPHA12 * RAD_TO_DEG); TAB;
+ (void)printf(oform, GEODESIC->ALPHA21 * RAD_TO_DEG); TAB;
+ (void)printf(osform, GEODESIC->DIST * GEODESIC->FR_METER);
} else {
- (void)fputs(rtodms(pline, al12, 0, 0), stdout); TAB;
- (void)fputs(rtodms(pline, al21, 0, 0), stdout); TAB;
- (void)printf(osform, geod_S * fr_meter);
+ (void)fputs(rtodms(pline, GEODESIC->ALPHA12, 0, 0), stdout); TAB;
+ (void)fputs(rtodms(pline, GEODESIC->ALPHA21, 0, 0), stdout); TAB;
+ (void)printf(osform, GEODESIC->DIST * GEODESIC->FR_METER);
}
} else if (inverse)
if (oform) {
- (void)printf(oform, al12 * RAD_TO_DEG); TAB;
- (void)printf(oform, al21 * RAD_TO_DEG); TAB;
- (void)printf(osform, geod_S * fr_meter);
+ (void)printf(oform, GEODESIC->ALPHA12 * RAD_TO_DEG); TAB;
+ (void)printf(oform, GEODESIC->ALPHA21 * RAD_TO_DEG); TAB;
+ (void)printf(osform, GEODESIC->DIST * GEODESIC->FR_METER);
} else {
- (void)fputs(rtodms(pline, al12, 0, 0), stdout); TAB;
- (void)fputs(rtodms(pline, al21, 0, 0), stdout); TAB;
- (void)printf(osform, geod_S * fr_meter);
+ (void)fputs(rtodms(pline, GEODESIC->ALPHA12, 0, 0), stdout); TAB;
+ (void)fputs(rtodms(pline, GEODESIC->ALPHA21, 0, 0), stdout); TAB;
+ (void)printf(osform, GEODESIC->DIST * GEODESIC->FR_METER);
}
else {
- printLL(phi2, lam2); TAB;
+ printLL(GEODESIC->p2.u, GEODESIC->p2.v); TAB;
if (oform)
- (void)printf(oform, al21 * RAD_TO_DEG);
+ (void)printf(oform, GEODESIC->ALPHA21 * RAD_TO_DEG);
else
- (void)fputs(rtodms(pline, al21, 0, 0), stdout);
+ (void)fputs(rtodms(pline, GEODESIC->ALPHA21, 0, 0), stdout);
}
(void)fputs(s, stdout);
}
@@ -209,12 +217,12 @@
eargv[eargc++] = *argv;
}
/* done with parameter and control input */
- geod_set(pargc, pargv); /* setup projection */
- if ((n_alpha || n_S) && eargc)
+ GEOD_init(pargc, pargv, GEODESIC); /* setup projection */
+ if ((GEODESIC->n_alpha || GEODESIC->n_S) && eargc)
emess(1,"files specified for arc/geodesic mode");
- if (n_alpha)
+ if (GEODESIC->n_alpha)
do_arc();
- else if (n_S)
+ else if (GEODESIC->n_S)
do_geod();
else { /* process input file list */
if (eargc == 0) /* if no specific files force sysin */
diff -Naur proj-4.5.0-orig/src/geodesic.h proj-4.5.0/src/geodesic.h
--- proj-4.5.0-orig/src/geodesic.h 2006-12-04 11:18:34.000000000 +0100
+++ proj-4.5.0/src/geodesic.h 2006-12-05 12:35:35.000000000 +0100
@@ -1,3 +1,8 @@
+#ifndef __GEODESIC_H__
+#define __GEODESIC_H__
+
+#include "projects.h"
+
#ifndef lint
static char GEODESIC_H_ID[] = "@(#)geodesic.h 4.3 95/08/19 GIE REL";
#endif
@@ -12,40 +17,37 @@
# define GEOD_EXTERN
#endif
-GEOD_EXTERN struct geodesic {
+typedef struct geodesic {
double A;
- double LAM1, PHI1, ALPHA12;
- double LAM2, PHI2, ALPHA21;
+
+ projUV p1, p2;
+
+ double ALPHA12;
+ double ALPHA21;
+
double DIST;
double ONEF, FLAT, FLAT2, FLAT4, FLAT64;
int ELLIPSE;
-} GEODESIC;
+ double FR_METER, TO_METER, del_alpha;
+ int n_alpha, n_S;
+
+
+ double th1,costh1,sinth1,sina12,cosa12,M,N,c1,c2,D,P,s1;
+ int merid, signS;
+} GEODESIC_T;
-# define geod_a GEODESIC.A
-# define lam1 GEODESIC.LAM1
-# define phi1 GEODESIC.PHI1
-# define al12 GEODESIC.ALPHA12
-# define lam2 GEODESIC.LAM2
-# define phi2 GEODESIC.PHI2
-# define al21 GEODESIC.ALPHA21
-# define geod_S GEODESIC.DIST
-# define geod_f GEODESIC.FLAT
-# define onef GEODESIC.ONEF
-# define f2 GEODESIC.FLAT2
-# define f4 GEODESIC.FLAT4
-# define ff2 GEODESIC.FLAT4
-# define f64 GEODESIC.FLAT64
-# define ellipse GEODESIC.ELLIPSE
-
-
-GEOD_EXTERN int n_alpha, n_S;
-GEOD_EXTERN double to_meter, fr_meter, del_alpha;
-void geod_set(int, char **);
-void geod_for(void);
-void geod_pre(void);
-void geod_inv(void);
+ GEODESIC_T *GEOD_init(int, char **, GEODESIC_T *g);
+ GEODESIC_T *GEOD_init_plus(const char *args, GEODESIC_T *g);
+ void geod_for(GEODESIC_T *g);
+ void geod_pre(GEODESIC_T *g);
+ int geod_inv(GEODESIC_T *g);
+
+
#ifdef __cplusplus
}
#endif
+
+#endif // __GEODESIC_H__
+
diff -Naur proj-4.5.0-orig/src/geod_for.c proj-4.5.0/src/geod_for.c
--- proj-4.5.0-orig/src/geod_for.c 2006-12-04 11:18:34.000000000 +0100
+++ proj-4.5.0/src/geod_for.c 2006-12-05 12:35:35.000000000 +0100
@@ -4,103 +4,110 @@
# include "projects.h"
# include "geodesic.h"
# define MERI_TOL 1e-9
- static double
-th1,costh1,sinth1,sina12,cosa12,M,N,c1,c2,D,P,s1;
- static int
-merid, signS;
- void
-geod_pre(void) {
- al12 = adjlon(al12); /* reduce to +- 0-PI */
- signS = fabs(al12) > HALFPI ? 1 : 0;
- th1 = ellipse ? atan(onef * tan(phi1)) : phi1;
- costh1 = cos(th1);
- sinth1 = sin(th1);
- if ((merid = fabs(sina12 = sin(al12)) < MERI_TOL)) {
- sina12 = 0.;
- cosa12 = fabs(al12) < HALFPI ? 1. : -1.;
- M = 0.;
+
+
+// input: al12, ELLIPSE, ONEF, phi1, FLAT, FLAT4
+// output: al12 (ajusting) and!!! double s1, D, P, c1, M, N, cosa12, sina12, sinth1, costh1, th1, int signS, merid
+
+void
+geod_pre(GEODESIC_T *GEODESIC) {
+ GEODESIC->ALPHA12 = adjlon(GEODESIC->ALPHA12); /* reduce to +- 0-PI */
+ GEODESIC->signS = fabs(GEODESIC->ALPHA12) > HALFPI ? 1 : 0;
+ GEODESIC->th1 = GEODESIC->ELLIPSE ? atan(GEODESIC->ONEF * tan(GEODESIC->p1.u)) : GEODESIC->p1.u;
+ GEODESIC->costh1 = cos(GEODESIC->th1);
+ GEODESIC->sinth1 = sin(GEODESIC->th1);
+ if ((GEODESIC->merid = fabs(GEODESIC->sina12 = sin(GEODESIC->ALPHA12)) < MERI_TOL)) {
+ GEODESIC->sina12 = 0.;
+ GEODESIC->cosa12 = fabs(GEODESIC->ALPHA12) < HALFPI ? 1. : -1.;
+ GEODESIC->M = 0.;
} else {
- cosa12 = cos(al12);
- M = costh1 * sina12;
+ GEODESIC->cosa12 = cos(GEODESIC->ALPHA12);
+ GEODESIC->M = GEODESIC->costh1 * GEODESIC->sina12;
}
- N = costh1 * cosa12;
- if (ellipse) {
- if (merid) {
- c1 = 0.;
- c2 = f4;
- D = 1. - c2;
- D *= D;
- P = c2 / D;
+ GEODESIC->N = GEODESIC->costh1 * GEODESIC->cosa12;
+ if (GEODESIC->ELLIPSE) {
+ if (GEODESIC->merid) {
+ GEODESIC->c1 = 0.;
+ GEODESIC->c2 = GEODESIC->FLAT4;
+ GEODESIC->D = 1. - GEODESIC->c2;
+ GEODESIC->D *= GEODESIC->D;
+ GEODESIC->P = GEODESIC->c2 / GEODESIC->D;
} else {
- c1 = geod_f * M;
- c2 = f4 * (1. - M * M);
- D = (1. - c2)*(1. - c2 - c1 * M);
- P = (1. + .5 * c1 * M) * c2 / D;
+ GEODESIC->c1 = GEODESIC->FLAT * GEODESIC->M;
+ GEODESIC->c2 = GEODESIC->FLAT4 * (1. - GEODESIC->M * GEODESIC->M);
+ GEODESIC->D = (1. - GEODESIC->c2)*(1. - GEODESIC->c2 - GEODESIC->c1 * GEODESIC->M);
+ GEODESIC->P = (1. + .5 * GEODESIC->c1 * GEODESIC->M) * GEODESIC->c2 / GEODESIC->D;
}
}
- if (merid) s1 = HALFPI - th1;
+ if (GEODESIC->merid) GEODESIC->s1 = HALFPI - GEODESIC->th1;
else {
- s1 = (fabs(M) >= 1.) ? 0. : acos(M);
- s1 = sinth1 / sin(s1);
- s1 = (fabs(s1) >= 1.) ? 0. : acos(s1);
+ GEODESIC->s1 = (fabs(GEODESIC->M) >= 1.) ? 0. : acos(GEODESIC->M);
+ GEODESIC->s1 = GEODESIC->sinth1 / sin(GEODESIC->s1);
+ GEODESIC->s1 = (fabs(GEODESIC->s1) >= 1.) ? 0. : acos(GEODESIC->s1);
}
}
- void
-geod_for(void) {
- double d,sind,u,V,X,ds,cosds,sinds,ss,de;
-
- if (ellipse) {
- d = geod_S / (D * geod_a);
- if (signS) d = -d;
- u = 2. * (s1 - d);
+
+// input: ELLIPSE, DIST, A and!!! D, signS, s1
+// output:
+
+void
+geod_for(GEODESIC_T *GEODESIC) {
+ double d,sind,u,V,X,ds,cosds,sinds,ss = 0,de;
+
+ if (GEODESIC->ELLIPSE) {
+ d = GEODESIC->DIST / (GEODESIC->D * GEODESIC->A);
+ if (GEODESIC->signS) d = -d;
+ u = 2. * (GEODESIC->s1 - d);
V = cos(u + d);
- X = c2 * c2 * (sind = sin(d)) * cos(d) * (2. * V * V - 1.);
- ds = d + X - 2. * P * V * (1. - 2. * P * cos(u)) * sind;
- ss = s1 + s1 - ds;
+ X = GEODESIC->c2 * GEODESIC->c2 * (sind = sin(d)) * cos(d) * (2. * V * V - 1.);
+ ds = d + X - 2. * GEODESIC->P * V * (1. - 2. * GEODESIC->P * cos(u)) * sind;
+ ss = GEODESIC->s1 + GEODESIC->s1 - ds;
} else {
- ds = geod_S / geod_a;
- if (signS) ds = - ds;
+ ds = GEODESIC->DIST / GEODESIC->A;
+ if (GEODESIC->signS) ds = - ds;
}
cosds = cos(ds);
sinds = sin(ds);
- if (signS) sinds = - sinds;
- al21 = N * cosds - sinth1 * sinds;
- if (merid) {
- phi2 = atan( tan(HALFPI + s1 - ds) / onef);
- if (al21 > 0.) {
- al21 = PI;
- if (signS)
+ if (GEODESIC->signS) sinds = - sinds;
+ GEODESIC->ALPHA21 = GEODESIC->N * cosds - GEODESIC->sinth1 * sinds;
+ if (GEODESIC->merid) {
+ GEODESIC->p2.u = atan( tan(HALFPI + GEODESIC->s1 - ds) / GEODESIC->ONEF);
+ if (GEODESIC->ALPHA21 > 0.) {
+ GEODESIC->ALPHA21 = PI;
+ if (GEODESIC->signS)
de = PI;
else {
- phi2 = - phi2;
+ GEODESIC->p2.u = - GEODESIC->p2.u;
de = 0.;
}
} else {
- al21 = 0.;
- if (signS) {
- phi2 = - phi2;
+ GEODESIC->ALPHA21 = 0.;
+ if (GEODESIC->signS) {
+ GEODESIC->p2.u = - GEODESIC->p2.u;
de = 0;
} else
de = PI;
}
} else {
- al21 = atan(M / al21);
- if (al21 > 0)
- al21 += PI;
- if (al12 < 0.)
- al21 -= PI;
- al21 = adjlon(al21);
- phi2 = atan(-(sinth1 * cosds + N * sinds) * sin(al21) /
- (ellipse ? onef * M : M));
- de = atan2(sinds * sina12 ,
- (costh1 * cosds - sinth1 * sinds * cosa12));
- if (ellipse)
- if (signS)
- de += c1 * ((1. - c2) * ds +
- c2 * sinds * cos(ss));
+ GEODESIC->ALPHA21 = atan(GEODESIC->M / GEODESIC->ALPHA21);
+ if (GEODESIC->ALPHA21 > 0)
+ GEODESIC->ALPHA21 += PI;
+ if (GEODESIC->ALPHA12 < 0.)
+ GEODESIC->ALPHA21 -= PI;
+ GEODESIC->ALPHA21 = adjlon(GEODESIC->ALPHA21);
+ GEODESIC->p2.u = atan(-(GEODESIC->sinth1 * cosds + GEODESIC->N * sinds) * sin(GEODESIC->ALPHA21) /
+ (GEODESIC->ELLIPSE ? GEODESIC->ONEF * GEODESIC->M : GEODESIC->M));
+ de = atan2(sinds * GEODESIC->sina12 ,
+ (GEODESIC->costh1 * cosds - GEODESIC->sinth1 * sinds * GEODESIC->cosa12));
+ if (GEODESIC->ELLIPSE)
+ {
+ if (GEODESIC->signS)
+ de += GEODESIC->c1 * ((1. - GEODESIC->c2) * ds +
+ GEODESIC->c2 * sinds * cos(ss));
else
- de -= c1 * ((1. - c2) * ds -
- c2 * sinds * cos(ss));
+ de -= GEODESIC->c1 * ((1. - GEODESIC->c2) * ds -
+ GEODESIC->c2 * sinds * cos(ss));
+ }
}
- lam2 = adjlon( lam1 + de );
+ GEODESIC->p2.v = adjlon( GEODESIC->p1.v + de );
}
diff -Naur proj-4.5.0-orig/src/geod_inv.c proj-4.5.0/src/geod_inv.c
--- proj-4.5.0-orig/src/geod_inv.c 2006-12-04 11:18:34.000000000 +0100
+++ proj-4.5.0/src/geod_inv.c 2006-12-05 12:35:35.000000000 +0100
@@ -4,24 +4,26 @@
# include "projects.h"
# include "geodesic.h"
# define DTOL 1e-12
- void
-geod_inv(void) {
+
+int
+geod_inv(GEODESIC_T *GEODESIC)
+{
double th1,th2,thm,dthm,dlamm,dlam,sindlamm,costhm,sinthm,cosdthm,
sindthm,L,E,cosd,d,X,Y,T,sind,tandlammp,u,v,D,A,B;
- if (ellipse) {
- th1 = atan(onef * tan(phi1));
- th2 = atan(onef * tan(phi2));
+ if (GEODESIC->ELLIPSE) {
+ th1 = atan(GEODESIC->ONEF * tan(GEODESIC->p1.u));
+ th2 = atan(GEODESIC->ONEF * tan(GEODESIC->p2.u));
} else {
- th1 = phi1;
- th2 = phi2;
+ th1 = GEODESIC->p1.u;
+ th2 = GEODESIC->p2.u;
}
thm = .5 * (th1 + th2);
dthm = .5 * (th2 - th1);
- dlamm = .5 * ( dlam = adjlon(lam2 - lam1) );
+ dlamm = .5 * ( dlam = adjlon(GEODESIC->p2.v - GEODESIC->p1.v) );
if (fabs(dlam) < DTOL && fabs(dthm) < DTOL) {
- al12 = al21 = geod_S = 0.;
- return;
+ GEODESIC->ALPHA12 = GEODESIC->ALPHA21 = GEODESIC->DIST = 0.;
+ return -1;
}
sindlamm = sin(dlamm);
costhm = cos(thm); sinthm = sin(thm);
@@ -29,7 +31,7 @@
L = sindthm * sindthm + (cosdthm * cosdthm - sinthm * sinthm)
* sindlamm * sindlamm;
d = acos(cosd = 1 - L - L);
- if (ellipse) {
+ if (GEODESIC->ELLIPSE) {
E = cosd + cosd;
sind = sin( d );
Y = sinthm * cosdthm;
@@ -42,18 +44,19 @@
D = 4. * T * T;
A = D * E;
B = D + D;
- geod_S = geod_a * sind * (T - f4 * (T * X - Y) +
- f64 * (X * (A + (T - .5 * (A - E)) * X) -
- Y * (B + E * Y) + D * X * Y));
+ GEODESIC->DIST = GEODESIC->A * sind * (T - GEODESIC->FLAT4 * (T * X - Y) +
+ GEODESIC->FLAT64 * (X * (A + (T - .5 * (A - E)) * X) -
+ Y * (B + E * Y) + D * X * Y));
tandlammp = tan(.5 * (dlam - .25 * (Y + Y - E * (4. - X)) *
- (f2 * T + f64 * (32. * T - (20. * T - A)
- * X - (B + 4.) * Y)) * tan(dlam)));
+ (GEODESIC->FLAT2 * T + GEODESIC->FLAT64 * (32. * T - (20. * T - A)
+ * X - (B + 4.) * Y)) * tan(dlam)));
} else {
- geod_S = geod_a * d;
+ GEODESIC->DIST = GEODESIC->A * d;
tandlammp = tan(dlamm);
}
u = atan2(sindthm , (tandlammp * costhm));
v = atan2(cosdthm , (tandlammp * sinthm));
- al12 = adjlon(TWOPI + v - u);
- al21 = adjlon(TWOPI - v - u);
+ GEODESIC->ALPHA12 = adjlon(TWOPI + v - u);
+ GEODESIC->ALPHA21 = adjlon(TWOPI - v - u);
+ return 0;
}
diff -Naur proj-4.5.0-orig/src/geod_set.c proj-4.5.0/src/geod_set.c
--- proj-4.5.0-orig/src/geod_set.c 2006-12-04 11:18:34.000000000 +0100
+++ proj-4.5.0/src/geod_set.c 2006-12-05 12:35:35.000000000 +0100
@@ -4,17 +4,27 @@
#define _IN_GEOD_SET
+#include <stdlib.h>
#include <string.h>
#include "projects.h"
#include "geodesic.h"
#include "emess.h"
- void
-geod_set(int argc, char **argv) {
- paralist *start = 0, *curr;
+
+GEODESIC_T *
+GEOD_init(int argc, char **argv, GEODESIC_T *GEODESIC)
+{
+ paralist *start = 0, *curr = 0;
double es;
char *name;
int i;
+
+ if(0 == GEODESIC)
+ {
+ GEODESIC = malloc(sizeof(GEODESIC_T));
+ }
+ memset(GEODESIC, 0, sizeof(GEODESIC_T));
+
/* put arguments into internal linked list */
if (argc <= 0)
emess(1, "no arguments in initialization list");
@@ -24,49 +34,50 @@
else
start = curr = pj_mkparam(argv[i]);
/* set elliptical parameters */
- if (pj_ell_set(start, &geod_a, &es)) emess(1,"ellipse setup failure");
+ if (pj_ell_set(start, &GEODESIC->A, &es)) emess(1,"ellipse setup failure");
/* set units */
- if (name = pj_param(start, "sunits").s) {
+ if ((name = pj_param(start, "sunits").s)) {
char *s;
struct PJ_UNITS *unit_list = pj_get_units_ref();
for (i = 0; (s = unit_list[i].id) && strcmp(name, s) ; ++i) ;
if (!s)
emess(1,"%s unknown unit conversion id", name);
- fr_meter = 1. / (to_meter = atof(unit_list[i].to_meter));
+ GEODESIC->FR_METER = 1. / (GEODESIC->TO_METER = atof(unit_list[i].to_meter));
} else
- to_meter = fr_meter = 1.;
- if (ellipse = es != 0.) {
- onef = sqrt(1. - es);
- geod_f = 1 - onef;
- f2 = geod_f/2;
- f4 = geod_f/4;
- f64 = geod_f*geod_f/64;
+ GEODESIC->TO_METER = GEODESIC->FR_METER = 1.;
+ if ((GEODESIC->ELLIPSE = (es != 0.))) {
+ GEODESIC->ONEF = sqrt(1. - es);
+ GEODESIC->FLAT = 1 - GEODESIC->ONEF;
+ GEODESIC->FLAT2 = GEODESIC->FLAT/2;
+ GEODESIC->FLAT4 = GEODESIC->FLAT/4;
+ GEODESIC->FLAT64 = GEODESIC->FLAT*GEODESIC->FLAT/64;
} else {
- onef = 1.;
- geod_f = f2 = f4 = f64 = 0.;
+ GEODESIC->ONEF = 1.;
+ GEODESIC->FLAT = GEODESIC->FLAT2 = GEODESIC->FLAT4 = GEODESIC->FLAT64 = 0.;
}
/* check if line or arc mode */
if (pj_param(start, "tlat_1").i) {
double del_S;
#undef f
- phi1 = pj_param(start, "rlat_1").f;
- lam1 = pj_param(start, "rlon_1").f;
+
+ GEODESIC->p1.u = pj_param(start, "rlat_1").f;
+ GEODESIC->p1.v = pj_param(start, "rlon_1").f;
if (pj_param(start, "tlat_2").i) {
- phi2 = pj_param(start, "rlat_2").f;
- lam2 = pj_param(start, "rlon_2").f;
- geod_inv();
- geod_pre();
- } else if (geod_S = pj_param(start, "dS").f) {
- al12 = pj_param(start, "rA").f;
- geod_pre();
- geod_for();
+ GEODESIC->p2.u = pj_param(start, "rlat_2").f;
+ GEODESIC->p2.v = pj_param(start, "rlon_2").f;
+ geod_inv(GEODESIC);
+ geod_pre(GEODESIC);
+ } else if ((GEODESIC->DIST = pj_param(start, "dS").f)) {
+ GEODESIC->ALPHA12 = pj_param(start, "rA").f;
+ geod_pre(GEODESIC);
+ geod_for(GEODESIC);
} else emess(1,"incomplete geodesic/arc info");
- if ((n_alpha = pj_param(start, "in_A").i) > 0) {
- if (!(del_alpha = pj_param(start, "rdel_A").f))
+ if ((GEODESIC->n_alpha = pj_param(start, "in_A").i) > 0) {
+ if (!(GEODESIC->del_alpha = pj_param(start, "rdel_A").f))
emess(1,"del azimuth == 0");
- } else if (del_S = fabs(pj_param(start, "ddel_S").f)) {
- n_S = geod_S / del_S + .5;
- } else if ((n_S = pj_param(start, "in_S").i) <= 0)
+ } else if ((del_S = fabs(pj_param(start, "ddel_S").f))) {
+ GEODESIC->n_S = GEODESIC->DIST / del_S + .5;
+ } else if ((GEODESIC->n_S = pj_param(start, "in_S").i) <= 0)
emess(1,"no interval divisor selected");
}
/* free up linked list */
@@ -74,4 +85,53 @@
curr = start->next;
pj_dalloc(start);
}
+ return GEODESIC;
}
+
+GEODESIC_T *GEOD_init_plus(const char *definition, GEODESIC_T *geod)
+{
+#define MAX_ARG 200
+ char *argv[MAX_ARG];
+ char *defn_copy;
+ int argc = 0, i;
+
+ /* make a copy that we can manipulate */
+ defn_copy = strdup(definition);
+
+ /* split into arguments based on '+' and trim white space */
+
+ for( i = 0; defn_copy[i] != '\0'; i++ )
+ {
+ switch( defn_copy[i] )
+ {
+ case '+':
+ if( i == 0 || defn_copy[i-1] == '\0' )
+ {
+ if( argc+1 == MAX_ARG )
+ {
+ //pj_errno = -44;
+ return NULL;
+ }
+
+ argv[argc++] = defn_copy + i + 1;
+ }
+ break;
+
+ case ' ':
+ case '\t':
+ case '\n':
+ defn_copy[i] = '\0';
+ break;
+
+ default:
+ /* do nothing */;
+ }
+ }
+
+ /* perform actual initialization */
+ GEODESIC_T *ret_geod = GEOD_init(argc, argv, geod);
+
+ free( defn_copy );
+ return ret_geod;
+}
+
diff -Naur proj-4.5.0-orig/src/Makefile.am proj-4.5.0/src/Makefile.am
--- proj-4.5.0-orig/src/Makefile.am 2006-12-04 11:18:34.000000000 +0100
+++ proj-4.5.0/src/Makefile.am 2006-12-05 12:35:35.000000000 +0100
@@ -2,7 +2,7 @@
INCLUDES = -DPROJ_LIB=\"$(pkgdatadir)\" @JNI_INCLUDE@
-include_HEADERS = projects.h nad_list.h proj_api.h org_proj4_Projections.h
+include_HEADERS = projects.h nad_list.h proj_api.h org_proj4_Projections.h geodesic.h
EXTRA_DIST = makefile.vc proj.def
@@ -10,7 +10,7 @@
cs2cs_SOURCES = cs2cs.c gen_cheb.c p_series.c
nad2nad_SOURCES = nad2nad.c
nad2bin_SOURCES = nad2bin.c
-geod_SOURCES = geod.c geod_set.c geod_for.c geod_inv.c geodesic.h
+geod_SOURCES = geod.c
proj_LDADD = libproj.la
cs2cs_LDADD = libproj.la
@@ -59,7 +59,8 @@
nad_cvt.c nad_init.c nad_intr.c emess.c emess.h \
pj_apply_gridshift.c pj_datums.c pj_datum_set.c pj_transform.c \
geocent.c geocent.h pj_utils.c pj_gridinfo.c pj_gridlist.c \
- jniproj.c
+ jniproj.c \
+ geod_set.c geod_for.c geod_inv.c geodesic.h
install-exec-local:
diff -Naur proj-4.5.0-orig/src/Makefile.in proj-4.5.0/src/Makefile.in
--- proj-4.5.0-orig/src/Makefile.in 2006-12-04 11:18:34.000000000 +0100
+++ proj-4.5.0/src/Makefile.in 2006-12-05 12:35:35.000000000 +0100
@@ -90,7 +90,7 @@
nad_cvt.lo nad_init.lo nad_intr.lo emess.lo \
pj_apply_gridshift.lo pj_datums.lo pj_datum_set.lo \
pj_transform.lo geocent.lo pj_utils.lo pj_gridinfo.lo \
- pj_gridlist.lo jniproj.lo
+ pj_gridlist.lo jniproj.lo geod_set.lo geod_for.lo geod_inv.lo
libproj_la_OBJECTS = $(am_libproj_la_OBJECTS)
binPROGRAMS_INSTALL = $(INSTALL_PROGRAM)
PROGRAMS = $(bin_PROGRAMS)
@@ -98,8 +98,7 @@
p_series.$(OBJEXT)
cs2cs_OBJECTS = $(am_cs2cs_OBJECTS)
cs2cs_DEPENDENCIES = libproj.la
-am_geod_OBJECTS = geod.$(OBJEXT) geod_set.$(OBJEXT) geod_for.$(OBJEXT) \
- geod_inv.$(OBJEXT)
+am_geod_OBJECTS = geod.$(OBJEXT)
geod_OBJECTS = $(am_geod_OBJECTS)
geod_DEPENDENCIES = libproj.la
am_nad2bin_OBJECTS = nad2bin.$(OBJEXT)
@@ -160,6 +159,7 @@
EXEEXT = @EXEEXT@
F77 = @F77@
FFLAGS = @FFLAGS@
+GREP = @GREP@
INSTALL_DATA = @INSTALL_DATA@
INSTALL_PROGRAM = @INSTALL_PROGRAM@
INSTALL_SCRIPT = @INSTALL_SCRIPT@
@@ -188,12 +188,9 @@
SHELL = @SHELL@
STRIP = @STRIP@
VERSION = @VERSION@
-ac_ct_AR = @ac_ct_AR@
ac_ct_CC = @ac_ct_CC@
ac_ct_CXX = @ac_ct_CXX@
ac_ct_F77 = @ac_ct_F77@
-ac_ct_RANLIB = @ac_ct_RANLIB@
-ac_ct_STRIP = @ac_ct_STRIP@
am__fastdepCC_FALSE = @am__fastdepCC_FALSE@
am__fastdepCC_TRUE = @am__fastdepCC_TRUE@
am__fastdepCXX_FALSE = @am__fastdepCXX_FALSE@
@@ -210,35 +207,42 @@
build_os = @build_os@
build_vendor = @build_vendor@
datadir = @datadir@
+datarootdir = @datarootdir@
+docdir = @docdir@
+dvidir = @dvidir@
exec_prefix = @exec_prefix@
host = @host@
host_alias = @host_alias@
host_cpu = @host_cpu@
host_os = @host_os@
host_vendor = @host_vendor@
+htmldir = @htmldir@
includedir = @includedir@
infodir = @infodir@
install_sh = @install_sh@
libdir = @libdir@
libexecdir = @libexecdir@
+localedir = @localedir@
localstatedir = @localstatedir@
mandir = @mandir@
mkdir_p = @mkdir_p@
oldincludedir = @oldincludedir@
+pdfdir = @pdfdir@
prefix = @prefix@
program_transform_name = @program_transform_name@
+psdir = @psdir@
sbindir = @sbindir@
sharedstatedir = @sharedstatedir@
sysconfdir = @sysconfdir@
target_alias = @target_alias@
INCLUDES = -DPROJ_LIB=\"$(pkgdatadir)\" @JNI_INCLUDE@
-include_HEADERS = projects.h nad_list.h proj_api.h org_proj4_Projections.h
+include_HEADERS = projects.h nad_list.h proj_api.h org_proj4_Projections.h geodesic.h
EXTRA_DIST = makefile.vc proj.def
proj_SOURCES = proj.c gen_cheb.c p_series.c
cs2cs_SOURCES = cs2cs.c gen_cheb.c p_series.c
nad2nad_SOURCES = nad2nad.c
nad2bin_SOURCES = nad2bin.c
-geod_SOURCES = geod.c geod_set.c geod_for.c geod_inv.c geodesic.h
+geod_SOURCES = geod.c
proj_LDADD = libproj.la
cs2cs_LDADD = libproj.la
nad2nad_LDADD = libproj.la
@@ -283,7 +287,8 @@
nad_cvt.c nad_init.c nad_intr.c emess.c emess.h \
pj_apply_gridshift.c pj_datums.c pj_datum_set.c pj_transform.c \
geocent.c geocent.h pj_utils.c pj_gridinfo.c pj_gridlist.c \
- jniproj.c
+ jniproj.c \
+ geod_set.c geod_for.c geod_inv.c geodesic.h
all: proj_config.h
$(MAKE) $(AM_MAKEFLAGS) all-am
@@ -515,9 +520,9 @@
@AMDEP_TRUE@@am__include@ @[EMAIL PROTECTED]/$(DEPDIR)/[EMAIL PROTECTED]@
@AMDEP_TRUE@@am__include@ @[EMAIL PROTECTED]/$(DEPDIR)/[EMAIL PROTECTED]@
@AMDEP_TRUE@@am__include@ @[EMAIL PROTECTED]/$(DEPDIR)/[EMAIL PROTECTED]@
[EMAIL PROTECTED]@@am__include@ @[EMAIL PROTECTED]/$(DEPDIR)/[EMAIL PROTECTED]@
[EMAIL PROTECTED]@@am__include@ @[EMAIL PROTECTED]/$(DEPDIR)/[EMAIL PROTECTED]@
[EMAIL PROTECTED]@@am__include@ @[EMAIL PROTECTED]/$(DEPDIR)/[EMAIL PROTECTED]@
[EMAIL PROTECTED]@@am__include@ @[EMAIL PROTECTED]/$(DEPDIR)/[EMAIL PROTECTED]@
[EMAIL PROTECTED]@@am__include@ @[EMAIL PROTECTED]/$(DEPDIR)/[EMAIL PROTECTED]@
[EMAIL PROTECTED]@@am__include@ @[EMAIL PROTECTED]/$(DEPDIR)/[EMAIL PROTECTED]@
@AMDEP_TRUE@@am__include@ @[EMAIL PROTECTED]/$(DEPDIR)/[EMAIL PROTECTED]@
@AMDEP_TRUE@@am__include@ @[EMAIL PROTECTED]/$(DEPDIR)/[EMAIL PROTECTED]@
@AMDEP_TRUE@@am__include@ @[EMAIL PROTECTED]/$(DEPDIR)/[EMAIL PROTECTED]@
