-
-
Save rouault/946104d0b98e8e8cc564 to your computer and use it in GitHub Desktop.
Regular code derived from PJ_tmer.c + pj_mlfn.c
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
// Compile with gcc -O2 testproj.c -o testproj -lm | |
#include <math.h> | |
#include <stdio.h> | |
#ifndef M_PI | |
# define M_PI 3.14159265358979323846 /* pi */ | |
#endif | |
#define HALFPI (M_PI/2) | |
#define EPS10 1.e-10 | |
#define FC1 1. | |
#define FC2 .5 | |
#define FC3 .16666666666666666666 | |
#define FC4 .08333333333333333333 | |
#define FC5 .05 | |
#define FC6 .03333333333333333333 | |
#define FC7 .02380952380952380952 | |
#define FC8 .01785714285714285714 | |
#define EN_SIZE 5 | |
#define C00 1. | |
#define C02 .25 | |
#define C04 .046875 | |
#define C06 .01953125 | |
#define C08 .01068115234375 | |
#define C22 .75 | |
#define C44 .46875 | |
#define C46 .01302083333333333333 | |
#define C48 .00712076822916666666 | |
#define C66 .36458333333333333333 | |
#define C68 .00569661458333333333 | |
#define C88 .3076171875 | |
#define EPS 1e-11 | |
#define MAX_ITER 10 | |
double | |
pj_mlfn(double phi, double sphi, double cphi, const double *en) { | |
cphi *= sphi; | |
sphi *= sphi; | |
return(en[0] * phi - cphi * (en[1] + sphi*(en[2] | |
+ sphi*(en[3] + sphi*en[4])))); | |
} | |
double | |
pj_inv_mlfn(double arg, double es, double _1_div_1_m_es, const double *en/*, | |
double* sinphinew, double* cosphinew*/) { | |
double s, t, phi, k = _1_div_1_m_es; | |
int i; | |
phi = arg; | |
for (i = MAX_ITER; i ; --i) { /* rarely goes over 2 iterations */ | |
double cosphi = cos(phi); | |
s = sin(phi); | |
t = 1. - es * s * s; | |
t = (pj_mlfn(phi, s, cosphi, en) - arg) * (t * sqrt(t)) * k; | |
if (fabs(t) < EPS) | |
{ | |
//sin(A+B)=sin A cos B + cos A sin B | |
//sin(A-B)=sin A cos B - cos A sin B | |
//cos(A+B)=cos A cos B - sin A sin B | |
//cos(A-B)=cos A cos B + sin A sin B | |
//*sinphinew = s - cosphi * t; | |
//*cosphinew = cosphi + s * t; | |
phi -= t; | |
return phi; | |
} | |
phi -= t; | |
} | |
//pj_ctx_set_errno( ctx, -17 ); | |
return phi; | |
} | |
void compute_mlfn_coeffs(double* en, double es) | |
{ | |
double t; | |
en[0] = C00 - es * (C02 + es * (C04 + es * (C06 + es * C08))); | |
en[1] = es * (C22 - es * (C04 + es * (C06 + es * C08))); | |
en[2] = (t = es * es) * (C44 - es * (C46 + es * C48)); | |
en[3] = (t *= es) * (C66 - es * C68); | |
en[4] = t * es * C88; | |
} | |
typedef struct { double x, y; } XY; | |
typedef struct { double lam, phi; } LP; | |
int main(int argc, char* argv[]) | |
{ | |
double y0 = 0; | |
double x0 = 500000.; | |
double lam0 = -3; | |
double k0 = 0.9996; | |
double phi0 = 0.; | |
double a = 6378137; | |
double ra = 1./a; | |
double invf = 298.257223563; | |
double f = 1 / invf; | |
double es = 2 *f - f*f; | |
double en[EN_SIZE]; | |
double ml0, esp; | |
int i; | |
double acc = 0; | |
double inv_k0 = 1. / k0; | |
double _1_div_1_m_es = 1. / (1.-es); | |
compute_mlfn_coeffs(en, es); | |
ml0 = pj_mlfn(phi0, sin(phi0), cos(phi0), en); | |
esp = es / (1. - es); | |
#define NITER (100*1000*1000) | |
for(i=0;i<NITER;i++) | |
{ | |
LP lp; | |
XY xy; | |
double n, con, cosphi, d, ds, sinphi, t; | |
xy.x = (1000000.*i) / NITER; | |
xy.y = (9000000.*i) / NITER; | |
xy.x *= ra; | |
xy.y *= ra; | |
lp.phi = pj_inv_mlfn(ml0 + xy.y * inv_k0, es, _1_div_1_m_es, en/*, | |
&sinphi, &cosphi*/); | |
if (fabs(lp.phi) >= HALFPI) { | |
lp.phi = xy.y < 0. ? -HALFPI : HALFPI; | |
lp.lam = 0.; | |
} else { | |
double invcosphi; | |
sinphi = sin(lp.phi); | |
cosphi = cos(lp.phi); | |
invcosphi = 1.0 / cosphi; | |
t = fabs(cosphi) > 1e-10 ? sinphi * invcosphi : 0.; | |
n = esp * cosphi * cosphi; | |
d = xy.x * sqrt(con = 1. - es * sinphi * sinphi) * inv_k0; | |
con *= t; | |
t *= t; | |
ds = d * d; | |
lp.phi -= (con * ds * _1_div_1_m_es) * FC2 * (1. - | |
ds * FC4 * (5. + t * (3. - 9. * n) + n * (1. - 4 * n) - | |
ds * FC6 * (61. + t * (90. - 252. * n + | |
45. * t) + 46. * n | |
- ds * FC8 * (1385. + t * (3633. + t * (4095. + 1574. * t)) ) | |
))); | |
lp.lam = d*(FC1 - | |
ds*FC3*( 1. + 2.*t + n - | |
ds*FC5*(5. + t*(28. + 24.*t + 8.*n) + 6.*n | |
- ds * FC7 * (61. + t * (662. + t * (1320. + 720. * t)) ) | |
))) * invcosphi; | |
} | |
//if( (i % (NITER/100)) == 0 ) | |
// printf("%.16g %.16g --> %.16g %.16g\n", xy.x * a, xy.y * a, | |
// lp.lam / M_PI * 180, lp.phi / M_PI * 180); | |
acc += lp.phi + lp.lam; | |
} | |
printf("Average of all longitude and latitudes: %.16g\n", acc / NITER); | |
return 0; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment