Skip to content

Instantly share code, notes, and snippets.

@rouault rouault/testproj.c Secret
Created Jun 23, 2015

Embed
What would you like to do?
Regular code derived from PJ_tmer.c + pj_mlfn.c
// 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
You can’t perform that action at this time.