com/gradoservice/proj4as/proj/ProjTmerc.as

   1 /*******************************************************************************\r
   2 NAME                            TRANSVERSE MERCATOR\r
   3 \r
   4 PURPOSE:        Transforms input longitude and latitude to Easting and\r
   5                 Northing for the Transverse Mercator projection.  The\r
   6                 longitude and latitude must be in radians.  The Easting\r
   7                 and Northing values will be returned in meters.\r
   8 \r
   9 ALGORITHM REFERENCES\r
  10 \r
  11 1.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological\r
  12     Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United\r
  13     State Government Printing Office, Washington D.C., 1987.\r
  14 \r
  15 2.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",\r
  16     U.S. Geological Survey Professional Paper 1453 , United State Government\r
  17     Printing Office, Washington D.C., 1989.\r
  18 *******************************************************************************/\r
  19 \r
  20 package com.gradoservice.proj4as.proj\r
  21 {\r
  22         import com.gradoservice.proj4as.ProjPoint;\r
  23         import com.gradoservice.proj4as.ProjConstants;\r
  24         import com.gradoservice.proj4as.Datum;\r
  25                         \r
  26         public class ProjTmerc extends AbstractProjProjection\r
  27         {\r
  28                 public function ProjTmerc(data:ProjParams)\r
  29                 {\r
  30                         super(data);\r
  31                 }\r
  32                 \r
  33   override public function init():void\r
  34         {\r
  35     this.e0 = ProjConstants.e0fn(this.es);\r
  36     this.e1 = ProjConstants.e1fn(this.es);\r
  37     this.e2 = ProjConstants.e2fn(this.es);\r
  38     this.e3 = ProjConstants.e3fn(this.es);\r
  39     this.ml0 = this.a * ProjConstants.mlfn(this.e0, this.e1, this.e2, this.e3, this.lat0);\r
  40     }\r
  41 \r
  42   /**\r
  43     Transverse Mercator Forward  - long/lat to x/y\r
  44     long/lat in radians\r
  45   */\r
  46   override public function forward(p:ProjPoint):ProjPoint \r
  47         {\r
  48     var lon:Number = p.x;\r
  49     var lat:Number = p.y;\r
  50 \r
  51     var delta_lon:Number = ProjConstants.adjust_lon(lon - this.long0); // Delta longitude\r
  52     var con:Number;    // cone constant\r
  53     var x:Number, y:Number;\r
  54     var sin_phi:Number=Math.sin(lat);\r
  55     var cos_phi:Number=Math.cos(lat);\r
  56 \r
  57     if (this.sphere) {  /* spherical form */\r
  58       var b:Number = cos_phi * Math.sin(delta_lon);\r
  59       if ((Math.abs(Math.abs(b) - 1.0)) < .0000000001)  {\r
  60         trace("tmerc:forward: Point projects into infinity");\r
  61         return null;\r
  62       } else {\r
  63         x = .5 * this.a * this.k0 * Math.log((1.0 + b)/(1.0 - b));\r
  64         con = Math.acos(cos_phi * Math.cos(delta_lon)/Math.sqrt(1.0 - b*b));\r
  65         if (lat < 0) con = - con;\r
  66         y = this.a * this.k0 * (con - this.lat0);\r
  67       }\r
  68     } else {\r
  69       var al:Number  = cos_phi * delta_lon;\r
  70       var als:Number = Math.pow(al,2);\r
  71       var c:Number   = this.ep2 * Math.pow(cos_phi,2);\r
  72       var tq:Number  = Math.tan(lat);\r
  73       var t:Number   = Math.pow(tq,2);\r
  74       con = 1.0 - this.es * Math.pow(sin_phi,2);\r
  75       var n:Number   = this.a / Math.sqrt(con);\r
  76       var ml:Number  = this.a * ProjConstants.mlfn(this.e0, this.e1, this.e2, this.e3, lat);\r
  77 \r
  78       x = this.k0 * n * al * (1.0 + als / 6.0 * (1.0 - t + c + als / 20.0 * (5.0 - 18.0 * t + Math.pow(t,2) + 72.0 * c - 58.0 * this.ep2))) + this.x0;\r
  79       y = this.k0 * (ml - this.ml0 + n * tq * (als * (0.5 + als / 24.0 * (5.0 - t + 9.0 * c + 4.0 * Math.pow(c,2) + als / 30.0 * (61.0 - 58.0 * t + Math.pow(t,2) + 600.0 * c - 330.0 * this.ep2))))) + this.y0;\r
  80 \r
  81     }\r
  82     p.x = x; p.y = y;\r
  83     return p;\r
  84   } // tmercFwd()\r
  85 \r
  86   /**\r
  87     Transverse Mercator Inverse  -  x/y to long/lat\r
  88   */\r
  89    override public function inverse(p:ProjPoint):ProjPoint\r
  90    {\r
  91     var con:Number, phi:Number;  /* temporary angles       */\r
  92     var delta_phi:Number; /* difference between longitudes    */\r
  93     var i:Number;\r
  94     var max_iter:Number = 6;      /* maximun number of iterations */\r
  95     var lat:Number, lon:Number;\r
  96 \r
  97     if (this.sphere) {   /* spherical form */\r
  98       var f:Number = Math.exp(p.x/(this.a * this.k0));\r
  99       var g:Number = .5 * (f - 1/f);\r
 100       var temp:Number = this.lat0 + p.y/(this.a * this.k0);\r
 101       var h:Number = Math.cos(temp);\r
 102       con = Math.sqrt((1.0 - h * h)/(1.0 + g * g));      \r
 103       lat = ProjConstants.asinz(con);\r
 104       if (temp < 0)\r
 105         lat = -lat;\r
 106       if ((g == 0) && (h == 0)) {\r
 107         lon = this.long0;\r
 108       } else {\r
 109         lon = ProjConstants.adjust_lon(Math.atan2(g,h) + this.long0);\r
 110       }\r
 111     } else {    // ellipsoidal form\r
 112       var x:Number= p.x - this.x0;\r
 113       var y:Number = p.y - this.y0;\r
 114 \r
 115       con = (this.ml0 + y / this.k0) / this.a;\r
 116       phi = con;\r
 117       for (i=0;;i++) {\r
 118         delta_phi=((con + this.e1 * Math.sin(2.0*phi) - this.e2 * Math.sin(4.0*phi) + this.e3 * Math.sin(6.0*phi)) / this.e0) - phi;\r
 119         phi += delta_phi;\r
 120         if (Math.abs(delta_phi) <= ProjConstants.EPSLN) break;\r
 121         if (i >= max_iter) {\r
 122           trace("tmerc:inverse: Latitude failed to converge");\r
 123           return null;\r
 124         }\r
 125       } // for()\r
 126       if (Math.abs(phi) < ProjConstants.HALF_PI) {\r
 127         // sincos(phi, &sin_phi, &cos_phi);\r
 128         var sin_phi:Number=Math.sin(phi);\r
 129         var cos_phi:Number=Math.cos(phi);\r
 130         var tan_phi:Number = Math.tan(phi);\r
 131         var c:Number = this.ep2 * Math.pow(cos_phi,2);\r
 132         var cs:Number = Math.pow(c,2);\r
 133         var t:Number = Math.pow(tan_phi,2);\r
 134         var ts :Number= Math.pow(t,2);\r
 135         con = 1.0 - this.es * Math.pow(sin_phi,2);\r
 136         var n:Number = this.a / Math.sqrt(con);\r
 137         var r:Number = n * (1.0 - this.es) / con;\r
 138         var d:Number = x / (n * this.k0);\r
 139         var ds:Number = Math.pow(d,2);\r
 140         lat = phi - (n * tan_phi * ds / r) * (0.5 - ds / 24.0 * (5.0 + 3.0 * t + 10.0 * c - 4.0 * cs - 9.0 * this.ep2 - ds / 30.0 * (61.0 + 90.0 * t + 298.0 * c + 45.0 * ts - 252.0 * this.ep2 - 3.0 * cs)));\r
 141         lon = ProjConstants.adjust_lon(this.long0 + (d * (1.0 - ds / 6.0 * (1.0 + 2.0 * t + c - ds / 20.0 * (5.0 - 2.0 * c + 28.0 * t - 3.0 * cs + 8.0 * this.ep2 + 24.0 * ts))) / cos_phi));\r
 142       } else {\r
 143         lat = ProjConstants.HALF_PI * ProjConstants.sign(y);\r
 144         lon = this.long0;\r
 145       }\r
 146     }\r
 147     p.x = lon;\r
 148     p.y = lat;\r
 149     return p;\r
 150   } // tmercInv()\r
 151                 \r
 152                 \r
 153         }\r
 154 }