vcs.maemo.org Git - navit-package/blob - navit/transform.c

   1 #include <assert.h>
   2 #include <stdio.h>
   3 #include <math.h>
   4 #include <limits.h>
   5 #include <glib.h>
   6 #include <string.h>
   7 #include "config.h"
   8 #include "coord.h"
   9 #include "debug.h"
  10 #include "map.h"
  11 #include "transform.h"
  12 #include "projection.h"
  13 #include "point.h"
  14
  15 struct transformation {
  16         long scale;             /* Scale factor */
  17         int angle;              /* Rotation angle */
  18         double cos_val,sin_val; /* cos and sin of rotation angle */
  19         enum projection pro;
  20         struct map_selection *map_sel;
  21         struct map_selection *screen_sel;
  22         struct point screen_center;
  23         struct coord map_center;        /* Center of source rectangle */
  24 };
  25
  26 struct transformation *
  27 transform_new(void)
  28 {
  29         struct transformation *this_;
  30
  31         this_=g_new0(struct transformation, 1);
  32
  33         return this_;
  34 }
  35
  36 static const double gar2geo_units = 360.0/(1<<24);
  37 static const double geo2gar_units = 1/(360.0/(1<<24));
  38
  39 void
  40 transform_to_geo(enum projection pro, struct coord *c, struct coord_geo *g)
  41 {
  42         switch (pro) {
  43         case projection_mg:
  44                 g->lng=c->x/6371000.0/M_PI*180;
  45                 g->lat=atan(exp(c->y/6371000.0))/M_PI*360-90;
  46                 break;
  47         case projection_garmin:
  48                 g->lng=c->x*gar2geo_units;
  49                 g->lat=c->y*gar2geo_units;
  50                 break;
  51         default:
  52                 break;
  53         }
  54 }
  55
  56 void
  57 transform_from_geo(enum projection pro, struct coord_geo *g, struct coord *c)
  58 {
  59         switch (pro) {
  60         case projection_mg:
  61                 c->x=g->lng*6371000.0*M_PI/180;
  62                 c->y=log(tan(M_PI_4+g->lat*M_PI/360))*6371000.0;
  63                 break;
  64         case projection_garmin:
  65                 c->x=g->lng*geo2gar_units;
  66                 c->y=g->lat*geo2gar_units;
  67                 break;
  68         default:
  69                 break;
  70         }
  71 }
  72
  73 void
  74 transform_from_to(struct coord *cfrom, enum projection from, struct coord *cto, enum projection to)
  75 {
  76         struct coord_geo g;
  77         transform_to_geo(from, cfrom, &g);
  78         transform_from_geo(to, &g, cto);
  79 }
  80
  81 void
  82 transform_geo_to_cart(struct coord_geo *geo, double a, double b, struct coord_geo_cart *cart)
  83 {
  84         double n,ee=1-b*b/(a*a);
  85         n = a/sqrt(1-ee*sin(geo->lat)*sin(geo->lat));
  86         cart->x=n*cos(geo->lat)*cos(geo->lng);
  87         cart->y=n*cos(geo->lat)*sin(geo->lng);
  88         cart->z=n*(1-ee)*sin(geo->lat);
  89 }
  90
  91 void
  92 transform_cart_to_geo(struct coord_geo_cart *cart, double a, double b, struct coord_geo *geo)
  93 {
  94         double lat,lati,n,ee=1-b*b/(a*a), lng = atan(cart->y/cart->x);
  95
  96         lat = atan(cart->z / sqrt((cart->x * cart->x) + (cart->y * cart->y)));
  97         do
  98         {
  99                 lati = lat;
 100
 101                 n = a / sqrt(1-ee*sin(lat)*sin(lat));
 102                 lat = atan((cart->z + ee * n * sin(lat)) / sqrt(cart->x * cart->x + cart->y * cart->y));
 103         }
 104         while (fabs(lat - lati) >= 0.000000000000001);
 105
 106         geo->lng=lng/M_PI*180;
 107         geo->lat=lat/M_PI*180;
 108 }
 109
 110
 111 void
 112 transform_datum(struct coord_geo *from, enum map_datum from_datum, struct coord_geo *to, enum map_datum to_datum)
 113 {
 114 }
 115
 116 int
 117 transform(struct transformation *t, enum projection pro, struct coord *c, struct point *p, int count, int unique)
 118 {
 119         struct coord c1;
 120         int xcn, ycn;
 121         struct coord_geo g;
 122 #ifdef AVOID_FLOAT
 123         int xc,yc;
 124 #else
 125         double xc,yc;
 126 #endif
 127         int i,j = 0;
 128         for (i=0; i < count; i++) {
 129                 if (pro == t->pro) {
 130                         xc=c[i].x;
 131                         yc=c[i].y;
 132                 } else {
 133                         transform_to_geo(pro, &c[i], &g);
 134                         transform_from_geo(t->pro, &g, &c1);
 135                         xc=c1.x;
 136                         yc=c1.y;
 137                 }
 138 //              dbg(2,"0x%x, 0x%x - 0x%x,0x%x contains 0x%x,0x%x\n", t->r.lu.x, t->r.lu.y, t->r.rl.x, t->r.rl.y, c->x, c->y);
 139 //              ret=coord_rect_contains(&t->r, c);
 140                 xc-=t->map_center.x;
 141                 yc-=t->map_center.y;
 142                 yc=-yc;
 143                 if (t->angle) {
 144                         xcn=xc*t->cos_val+yc*t->sin_val;
 145                         ycn=-xc*t->sin_val+yc*t->cos_val;
 146                         xc=xcn;
 147                         yc=ycn;
 148                 }
 149                 xc=xc*16;
 150                 yc=yc*16;
 151 #ifndef AVOID_FLOAT
 152                 if (t->scale!=1) {
 153                         xc=xc/(double)(t->scale);
 154                         yc=yc/(double)(t->scale);
 155                 }
 156 #else
 157                 if (t->scale!=1) {
 158                         xc=xc/t->scale;
 159                         yc=yc/t->scale;
 160                 }
 161 #endif
 162                 xc+=t->screen_center.x;
 163                 yc+=t->screen_center.y;
 164                 if (xc < -0x8000)
 165                         xc=-0x8000;
 166                 if (xc > 0x7fff) {
 167                         xc=0x7fff;
 168                 }
 169                 if (yc < -0x8000)
 170                         yc=-0x8000;
 171                 if (yc > 0x7fff)
 172                         yc=0x7fff;
 173                 if (j == 0 || !unique || p[j-1].x != xc || p[j-1].y != yc) {
 174                         p[j].x=xc;
 175                         p[j].y=yc;
 176                         j++;
 177                 }
 178         }
 179         return j;
 180 }
 181
 182 void
 183 transform_reverse(struct transformation *t, struct point *p, struct coord *c)
 184 {
 185         int xc,yc;
 186         xc=p->x;
 187         yc=p->y;
 188         xc-=t->screen_center.x;
 189         yc-=t->screen_center.y;
 190         xc=xc*t->scale/16;
 191         yc=-yc*t->scale/16;
 192         if (t->angle) {
 193                 int xcn, ycn;
 194                 xcn=xc*t->cos_val+yc*t->sin_val;
 195                 ycn=-xc*t->sin_val+yc*t->cos_val;
 196                 xc=xcn;
 197                 yc=ycn;
 198         }
 199         c->x=t->map_center.x+xc;
 200         c->y=t->map_center.y+yc;
 201 }
 202
 203 enum projection
 204 transform_get_projection(struct transformation *this_)
 205 {
 206         return this_->pro;
 207 }
 208
 209 void
 210 transform_set_projection(struct transformation *this_, enum projection pro)
 211 {
 212         this_->pro=pro;
 213 }
 214
 215 static int
 216 min4(int v1,int v2, int v3, int v4)
 217 {
 218         int res=v1;
 219         if (v2 < res)
 220                 res=v2;
 221         if (v3 < res)
 222                 res=v3;
 223         if (v4 < res)
 224                 res=v4;
 225         return res;
 226 }
 227
 228 static int
 229 max4(int v1,int v2, int v3, int v4)
 230 {
 231         int res=v1;
 232         if (v2 > res)
 233                 res=v2;
 234         if (v3 > res)
 235                 res=v3;
 236         if (v4 > res)
 237                 res=v4;
 238         return res;
 239 }
 240
 241 struct map_selection *
 242 transform_get_selection(struct transformation *this_, enum projection pro, int order)
 243 {
 244
 245         struct map_selection *ret,*curri,*curro;
 246         struct coord_geo g;
 247         int i;
 248
 249         ret=map_selection_dup(this_->map_sel);
 250         curri=this_->map_sel;
 251         curro=ret;
 252         while (curri) {
 253                 if (this_->pro != pro) {
 254                         transform_to_geo(this_->pro, &curri->u.c_rect.lu, &g);
 255                         transform_from_geo(pro, &g, &curro->u.c_rect.lu);
 256                         dbg(1,"%f,%f", g.lat, g.lng);
 257                         transform_to_geo(this_->pro, &curri->u.c_rect.rl, &g);
 258                         transform_from_geo(pro, &g, &curro->u.c_rect.rl);
 259                         dbg(1,": - %f,%f\n", g.lat, g.lng);
 260                 }
 261                 dbg(1,"transform rect for %d is %d,%d - %d,%d\n", pro, curro->u.c_rect.lu.x, curro->u.c_rect.lu.y, curro->u.c_rect.rl.x, curro->u.c_rect.rl.y);
 262                 for (i = 0 ; i < layer_end ; i++)
 263                         curro->order[i]+=order;
 264                 curri=curri->next;
 265                 curro=curro->next;
 266         }
 267         return ret;
 268 }
 269
 270 struct coord *
 271 transform_center(struct transformation *this_)
 272 {
 273         return &this_->map_center;
 274 }
 275
 276 void
 277 transform_set_angle(struct transformation *t,int angle)
 278 {
 279         t->angle=angle;
 280         t->cos_val=cos(M_PI*t->angle/180);
 281         t->sin_val=sin(M_PI*t->angle/180);
 282 }
 283
 284 int
 285 transform_get_angle(struct transformation *this_,int angle)
 286 {
 287         return this_->angle;
 288 }
 289
 290 void
 291 transform_set_screen_selection(struct transformation *t, struct map_selection *sel)
 292 {
 293         map_selection_destroy(t->screen_sel);
 294         t->screen_sel=map_selection_dup(sel);
 295         if (sel) {
 296                 t->screen_center.x=(sel->u.p_rect.rl.x-sel->u.p_rect.lu.x)/2;
 297                 t->screen_center.y=(sel->u.p_rect.rl.y-sel->u.p_rect.lu.y)/2;
 298         }
 299 }
 300
 301 #if 0
 302 void
 303 transform_set_size(struct transformation *t, int width, int height)
 304 {
 305         t->width=width;
 306         t->height=height;
 307 }
 308 #endif
 309
 310 void
 311 transform_get_size(struct transformation *t, int *width, int *height)
 312 {
 313         struct point_rect *r;
 314         if (t->screen_sel) {
 315                 r=&t->screen_sel->u.p_rect;
 316                 *width=r->rl.x-r->lu.x;
 317                 *height=r->rl.y-r->lu.y;
 318         }
 319 }
 320
 321 void
 322 transform_setup(struct transformation *t, struct pcoord *c, int scale, int angle)
 323 {
 324         t->pro=c->pro;
 325         t->map_center.x=c->x;
 326         t->map_center.y=c->y;
 327         t->scale=scale;
 328         transform_set_angle(t, angle);
 329 }
 330
 331 #if 0
 332
 333 void
 334 transform_setup_source_rect_limit(struct transformation *t, struct coord *center, int limit)
 335 {
 336         t->center=*center;
 337         t->scale=1;
 338         t->angle=0;
 339         t->r.lu.x=center->x-limit;
 340         t->r.rl.x=center->x+limit;
 341         t->r.rl.y=center->y-limit;
 342         t->r.lu.y=center->y+limit;
 343 }
 344 #endif
 345
 346 void
 347 transform_setup_source_rect(struct transformation *t)
 348 {
 349         int i;
 350         struct coord screen[4];
 351         struct point screen_pnt[4];
 352         struct point_rect *pr;
 353         struct map_selection *ms,*msm,*next,**msm_last;
 354         ms=t->map_sel;
 355         while (ms) {
 356                 next=ms->next;
 357                 g_free(ms);
 358                 ms=next;
 359         }
 360         t->map_sel=NULL;
 361         msm_last=&t->map_sel;
 362         ms=t->screen_sel;
 363         while (ms) {
 364                 msm=g_new0(struct map_selection, 1);
 365                 *msm=*ms;
 366                 pr=&ms->u.p_rect;
 367                 screen_pnt[0].x=pr->lu.x;
 368                 screen_pnt[0].y=pr->lu.y;
 369                 screen_pnt[1].x=pr->rl.x;
 370                 screen_pnt[1].y=pr->lu.y;
 371                 screen_pnt[2].x=pr->lu.x;
 372                 screen_pnt[2].y=pr->rl.y;
 373                 screen_pnt[3].x=pr->rl.x;
 374                 screen_pnt[3].y=pr->rl.y;
 375                 for (i = 0 ; i < 4 ; i++) {
 376                         transform_reverse(t, &screen_pnt[i], &screen[i]);
 377                 }
 378                 msm->u.c_rect.lu.x=min4(screen[0].x,screen[1].x,screen[2].x,screen[3].x);
 379                 msm->u.c_rect.rl.x=max4(screen[0].x,screen[1].x,screen[2].x,screen[3].x);
 380                 msm->u.c_rect.rl.y=min4(screen[0].y,screen[1].y,screen[2].y,screen[3].y);
 381                 msm->u.c_rect.lu.y=max4(screen[0].y,screen[1].y,screen[2].y,screen[3].y);
 382                 *msm_last=msm;
 383                 msm_last=&msm->next;
 384                 ms=ms->next;
 385         }
 386 }
 387
 388 long
 389 transform_get_scale(struct transformation *t)
 390 {
 391         return t->scale;
 392 }
 393
 394 void
 395 transform_set_scale(struct transformation *t, long scale)
 396 {
 397         t->scale=scale;
 398 }
 399
 400
 401 int
 402 transform_get_order(struct transformation *t)
 403 {
 404         int scale=t->scale;
 405         int order=0;
 406         while (scale > 1) {
 407                 order++;
 408                 scale>>=1;
 409         }
 410         order=18-order;
 411         if (order < 0)
 412                 order=0;
 413         return order;
 414 }
 415
 416
 417 void
 418 transform_geo_text(struct coord_geo *g, char *buffer)
 419 {
 420         double lng=g->lng;
 421         double lat=g->lat;
 422         char lng_c='E';
 423         char lat_c='N';
 424
 425         if (lng < 0) {
 426                 lng=-lng;
 427                 lng_c='W';
 428         }
 429         if (lat < 0) {
 430                 lat=-lat;
 431                 lat_c='S';
 432         }
 433
 434         sprintf(buffer,"%02.0f%07.4f%c %03.0f%07.4f%c", floor(lat), fmod(lat*60,60), lat_c, floor(lng), fmod(lng*60,60), lng_c);
 435
 436 }
 437
 438 #define TWOPI (M_PI*2)
 439 #define GC2RAD(c) ((c) * TWOPI/(1<<24))
 440 #define minf(a,b) ((a) < (b) ? (a) : (b))
 441
 442 static double
 443 transform_distance_garmin(struct coord *c1, struct coord *c2)
 444 {
 445 #ifdef USE_HALVESINE
 446         static const int earth_radius = 6371*1000; //m change accordingly
 447 // static const int earth_radius  = 3960; //miles
 448
 449 //Point 1 cords
 450         float lat1  = GC2RAD(c1->y);
 451         float long1 = GC2RAD(c1->x);
 452
 453 //Point 2 cords
 454         float lat2  = GC2RAD(c2->y);
 455         float long2 = GC2RAD(c2->x);
 456
 457 //Haversine Formula
 458         float dlong = long2-long1;
 459         float dlat  = lat2-lat1;
 460
 461         float sinlat  = sinf(dlat/2);
 462         float sinlong = sinf(dlong/2);
 463
 464         float a=(sinlat*sinlat)+cosf(lat1)*cosf(lat2)*(sinlong*sinlong);
 465         float c=2*asinf(minf(1,sqrt(a)));
 466 #ifdef AVOID_FLOAT
 467         return round(earth_radius*c);
 468 #else
 469         return earth_radius*c;
 470 #endif
 471 #else
 472 #define GMETER 2.3887499999999999
 473         double dx,dy;
 474         dx=c1->x-c2->x;
 475         dy=c1->y-c2->y;
 476         return sqrt(dx*dx+dy*dy)*GMETER;
 477 #undef GMETER
 478 #endif
 479 }
 480
 481 double
 482 transform_scale(int y)
 483 {
 484         struct coord c;
 485         struct coord_geo g;
 486         c.x=0;
 487         c.y=y;
 488         transform_to_geo(projection_mg, &c, &g);
 489         return 1/cos(g.lat/180*M_PI);
 490 }
 491
 492 #ifdef AVOID_FLOAT
 493 static int
 494 tab_sqrt[]={14142,13379,12806,12364,12018,11741,11517,11333,11180,11051,10943,10850,10770,10701,10640,10587,10540,10499,10462,10429,10400,10373,10349,10327,10307,10289,10273,10257,10243,10231,10219,10208};
 495 #endif
 496
 497 double
 498 transform_distance(enum projection pro, struct coord *c1, struct coord *c2)
 499 {
 500         if (pro == projection_mg) {
 501 #ifndef AVOID_FLOAT
 502         double dx,dy,scale=transform_scale((c1->y+c2->y)/2);
 503         dx=c1->x-c2->x;
 504         dy=c1->y-c2->y;
 505         return sqrt(dx*dx+dy*dy)/scale;
 506 #else
 507         int dx,dy,f,scale=15539;
 508         dx=c1->x-c2->x;
 509         dy=c1->y-c2->y;
 510         if (dx < 0)
 511                 dx=-dx;
 512         if (dy < 0)
 513                 dy=-dy;
 514         while (dx > 20000 || dy > 20000) {
 515                 dx/=10;
 516                 dy/=10;
 517                 scale/=10;
 518         }
 519         if (! dy)
 520                 return dx*10000/scale;
 521         if (! dx)
 522                 return dy*10000/scale;
 523         if (dx > dy) {
 524                 f=dx*8/dy-8;
 525                 if (f >= 32)
 526                         return dx*10000/scale;
 527                 return dx*tab_sqrt[f]/scale;
 528         } else {
 529                 f=dy*8/dx-8;
 530                 if (f >= 32)
 531                         return dy*10000/scale;
 532                 return dy*tab_sqrt[f]/scale;
 533         }
 534 #endif
 535         } else if (pro == projection_garmin) {
 536                 return transform_distance_garmin(c1, c2);
 537         } else {
 538                 printf("Unknown projection: %d\n", pro);
 539                 return 0;
 540         }
 541 }
 542
 543 int
 544 transform_distance_sq(struct coord *c1, struct coord *c2)
 545 {
 546         int dx=c1->x-c2->x;
 547         int dy=c1->y-c2->y;
 548
 549         if (dx > 32767 || dy > 32767 || dx < -32767 || dy < -32767)
 550                 return INT_MAX;
 551         else
 552                 return dx*dx+dy*dy;
 553 }
 554
 555 int
 556 transform_distance_line_sq(struct coord *l0, struct coord *l1, struct coord *ref, struct coord *lpnt)
 557 {
 558         int vx,vy,wx,wy;
 559         int c1,c2;
 560         int climit=1000000;
 561         struct coord l;
 562
 563         vx=l1->x-l0->x;
 564         vy=l1->y-l0->y;
 565         wx=ref->x-l0->x;
 566         wy=ref->y-l0->y;
 567
 568         c1=vx*wx+vy*wy;
 569         if ( c1 <= 0 ) {
 570                 if (lpnt)
 571                         *lpnt=*l0;
 572                 return transform_distance_sq(l0, ref);
 573         }
 574         c2=vx*vx+vy*vy;
 575         if ( c2 <= c1 ) {
 576                 if (lpnt)
 577                         *lpnt=*l1;
 578                 return transform_distance_sq(l1, ref);
 579         }
 580         while (c1 > climit || c2 > climit) {
 581                 c1/=256;
 582                 c2/=256;
 583         }
 584         l.x=l0->x+vx*c1/c2;
 585         l.y=l0->y+vy*c1/c2;
 586         if (lpnt)
 587                 *lpnt=l;
 588         return transform_distance_sq(&l, ref);
 589 }
 590
 591 int
 592 transform_distance_polyline_sq(struct coord *c, int count, struct coord *ref, struct coord *lpnt, int *pos)
 593 {
 594         int i,dist,distn;
 595         struct coord lp;
 596         if (count < 2)
 597                 return INT_MAX;
 598         if (pos)
 599                 *pos=0;
 600         dist=transform_distance_line_sq(&c[0], &c[1], ref, lpnt);
 601         for (i=2 ; i < count ; i++) {
 602                 distn=transform_distance_line_sq(&c[i-1], &c[i], ref, &lp);
 603                 if (distn < dist) {
 604                         dist=distn;
 605                         if (lpnt)
 606                                 *lpnt=lp;
 607                         if (pos)
 608                                 *pos=i-1;
 609                 }
 610         }
 611         return dist;
 612 }
 613
 614
 615 void
 616 transform_print_deg(double deg)
 617 {
 618         printf("%2.0f:%2.0f:%2.4f", floor(deg), fmod(deg*60,60), fmod(deg*3600,60));
 619 }
 620
 621 #ifdef AVOID_FLOAT
 622 static int tab_atan[]={0,262,524,787,1051,1317,1584,1853,2126,2401,2679,2962,3249,3541,3839,4142,4452,4770,5095,5430,5774,6128,6494,6873,7265,7673,8098,8541,9004,9490,10000,10538};
 623
 624 static int
 625 atan2_int_lookup(int val)
 626 {
 627         int len=sizeof(tab_atan)/sizeof(int);
 628         int i=len/2;
 629         int p=i-1;
 630         for (;;) {
 631                 i>>=1;
 632                 if (val < tab_atan[p])
 633                         p-=i;
 634                 else
 635                         if (val < tab_atan[p+1])
 636                                 return p+(p>>1);
 637                         else
 638                                 p+=i;
 639         }
 640 }
 641
 642 static int
 643 atan2_int(int dx, int dy)
 644 {
 645         int f,mul=1,add=0,ret;
 646         if (! dx) {
 647                 return dy < 0 ? 180 : 0;
 648         }
 649         if (! dy) {
 650                 return dx < 0 ? -90 : 90;
 651         }
 652         if (dx < 0) {
 653                 dx=-dx;
 654                 mul=-1;
 655         }
 656         if (dy < 0) {
 657                 dy=-dy;
 658                 add=180*mul;
 659                 mul*=-1;
 660         }
 661         while (dx > 20000 || dy > 20000) {
 662                 dx/=10;
 663                 dy/=10;
 664         }
 665         if (dx > dy) {
 666                 ret=90-atan2_int_lookup(dy*10000/dx);
 667         } else {
 668                 ret=atan2_int_lookup(dx*10000/dy);
 669         }
 670         return ret*mul+add;
 671 }
 672 #endif
 673
 674 int
 675 transform_get_angle_delta(struct coord *c1, struct coord *c2, int dir)
 676 {
 677         int dx=c2->x-c1->x;
 678         int dy=c2->y-c1->y;
 679 #ifndef AVOID_FLOAT
 680         double angle;
 681         angle=atan2(dx,dy);
 682         angle*=180/M_PI;
 683 #else
 684         int angle;
 685         angle=atan2_int(dx,dy);
 686 #endif
 687         if (dir == -1)
 688                 angle=angle-180;
 689         if (angle < 0)
 690                 angle+=360;
 691         return angle;
 692 }
 693
 694 int
 695 transform_within_border(struct transformation *this_, struct point *p, int border)
 696 {
 697         struct map_selection *ms=this_->screen_sel;
 698         while (ms) {
 699                 struct point_rect *r=&ms->u.p_rect;
 700                 if (p->x >= r->lu.x+border && p->x <= r->rl.x-border &&
 701                     p->y >= r->lu.y+border && p->y <= r->rl.y-border)
 702                         return 1;
 703                 ms=ms->next;
 704         }
 705         return 0;
 706 }
 707
 708 /*
 709 Note: there are many mathematically equivalent ways to express these formulas. As usual, not all of them are computationally equivalent.
 710
 711 L = latitude in radians (positive north)
 712 Lo = longitude in radians (positive east)
 713 E = easting (meters)
 714 N = northing (meters)
 715
 716 For the sphere
 717
 718 E = r Lo
 719 N = r ln [ tan (pi/4 + L/2) ]
 720
 721 where
 722
 723 r = radius of the sphere (meters)
 724 ln() is the natural logarithm
 725
 726 For the ellipsoid
 727
 728 E = a Lo
 729 N = a * ln ( tan (pi/4 + L/2) * ( (1 - e * sin (L)) / (1 + e * sin (L))) ** (e/2)  )
 730
 731
 732                                                e
 733                                                -
 734                    pi     L     1 - e sin(L)   2
 735     =  a ln( tan( ---- + ---) (--------------)   )
 736                    4      2     1 + e sin(L)
 737
 738
 739 where
 740
 741 a = the length of the semi-major axis of the ellipsoid (meters)
 742 e = the first eccentricity of the ellipsoid
 743
 744
 745 */
 746
 747