com/gradoservice/proj4as/proj/ProjSoMerc.as

   1 /*******************************************************************************\r
   2 NAME                       SWISS OBLIQUE MERCATOR\r
   3 \r
   4 PURPOSE:        Swiss projection.\r
   5 WARNING:  X and Y are inverted (weird) in the swiss coordinate system. Not\r
   6    here, since we want X to be horizontal and Y vertical.\r
   7 \r
   8 ALGORITHM REFERENCES\r
   9 1. "Formules et constantes pour le Calcul pour la\r
  10  projection cylindrique conforme Г  axe oblique et pour la transformation entre\r
  11  des systГЁmes de rГ©fГ©rence".\r
  12  http://www.swisstopo.admin.ch/internet/swisstopo/fr/home/topics/survey/sys/refsys/switzerland.parsysrelated1.31216.downloadList.77004.DownloadFile.tmp/swissprojectionfr.pdf\r
  13 \r
  14 *******************************************************************************/\r
  15 \r
  16 \r
  17 package com.gradoservice.proj4as.proj\r
  18 {\r
  19         import com.gradoservice.proj4as.ProjPoint;\r
  20         import com.gradoservice.proj4as.ProjConstants;\r
  21         import com.gradoservice.proj4as.Datum;\r
  22                 \r
  23         public class ProjSoMerc extends AbstractProjProjection\r
  24         {\r
  25                 private var lambda0:Number;\r
  26                 private var R:Number;\r
  27                 private var b0:Number;\r
  28                 private var K:Number;\r
  29                 \r
  30                 \r
  31                 public function ProjSoMerc(data:ProjParams)\r
  32                 {\r
  33                         super(data);\r
  34                 }\r
  35                 \r
  36                 \r
  37                   override public function init():void\r
  38                   {\r
  39                     var phy0:Number = this.lat0;\r
  40                     this.lambda0 = this.long0;\r
  41                     var sinPhy0:Number = Math.sin(phy0);\r
  42                     var semiMajorAxis:Number = this.a;\r
  43                     var invF:Number = this.rf;\r
  44                     var flattening:Number = 1 / invF;\r
  45                     var e2:Number = 2 * flattening - Math.pow(flattening, 2);\r
  46                     var e:Number = this.e = Math.sqrt(e2);\r
  47                     this.R = semiMajorAxis * Math.sqrt(1 - e2) / (1 - e2 * Math.pow(sinPhy0, 2.0));\r
  48                     this.alpha = Math.sqrt(1 + e2 / (1 - e2) * Math.pow(Math.cos(phy0), 4.0));\r
  49                     this.b0 = Math.asin(sinPhy0 / this.alpha);\r
  50                     this.K = Math.log(Math.tan(Math.PI / 4.0 + this.b0 / 2.0))\r
  51                             - this.alpha\r
  52                             * Math.log(Math.tan(Math.PI / 4.0 + phy0 / 2.0))\r
  53                             + this.alpha\r
  54                             * e / 2\r
  55                             * Math.log((1 + e * sinPhy0)\r
  56                             / (1 - e * sinPhy0));\r
  57                   }\r
  58                 \r
  59                 \r
  60                   override public function forward(p:ProjPoint):ProjPoint\r
  61                   {\r
  62                     var Sa1:Number = Math.log(Math.tan(Math.PI / 4.0 - p.y / 2.0));\r
  63                     var Sa2:Number = this.e / 2.0\r
  64                             * Math.log((1 + this.e * Math.sin(p.y))\r
  65                             / (1 - this.e * Math.sin(p.y)));\r
  66                     var S:Number = -this.alpha * (Sa1 + Sa2) + this.K;\r
  67                 \r
  68                         // spheric latitude\r
  69                     var b:Number = 2.0 * (Math.atan(Math.exp(S)) - Math.PI / 4.0);\r
  70                 \r
  71                         // spheric longitude\r
  72                     var I:Number = this.alpha * (p.x - this.lambda0);\r
  73                 \r
  74                         // psoeudo equatorial rotation\r
  75                     var rotI:Number = Math.atan(Math.sin(I)\r
  76                             / (Math.sin(this.b0) * Math.tan(b) +\r
  77                                Math.cos(this.b0) * Math.cos(I)));\r
  78                 \r
  79                     var rotB:Number = Math.asin(Math.cos(this.b0) * Math.sin(b) -\r
  80                                          Math.sin(this.b0) * Math.cos(b) * Math.cos(I));\r
  81                 \r
  82                     p.y = this.R / 2.0\r
  83                             * Math.log((1 + Math.sin(rotB)) / (1 - Math.sin(rotB)))\r
  84                             + this.y0;\r
  85                     p.x = this.R * rotI + this.x0;\r
  86                     return p;\r
  87                   }\r
  88                 \r
  89                  override public function inverse(p:ProjPoint):ProjPoint\r
  90                  {\r
  91                     var Y:Number = p.x - this.x0;\r
  92                     var X:Number = p.y - this.y0;\r
  93                 \r
  94                     var rotI:Number = Y / this.R;\r
  95                     var rotB:Number = 2 * (Math.atan(Math.exp(X / this.R)) - Math.PI / 4.0);\r
  96                 \r
  97                     var b:Number = Math.asin(Math.cos(this.b0) * Math.sin(rotB)\r
  98                             + Math.sin(this.b0) * Math.cos(rotB) * Math.cos(rotI));\r
  99                     var I:Number = Math.atan(Math.sin(rotI)\r
 100                             / (Math.cos(this.b0) * Math.cos(rotI) - Math.sin(this.b0)\r
 101                             * Math.tan(rotB)));\r
 102                 \r
 103                     var lambda:Number = this.lambda0 + I / this.alpha;\r
 104                 \r
 105                     var S:Number = 0.0;\r
 106                     var phy:Number = b;\r
 107                     var prevPhy:Number = -1000.0;\r
 108                     var iteration:Number = 0;\r
 109                     while (Math.abs(phy - prevPhy) > 0.0000001)\r
 110                     {\r
 111                       if (++iteration > 20)\r
 112                       {\r
 113                         trace("omercFwdInfinity");\r
 114                         return null;\r
 115                       }\r
 116                       //S = Math.log(Math.tan(Math.PI / 4.0 + phy / 2.0));\r
 117                       S = 1.0\r
 118                               / this.alpha\r
 119                               * (Math.log(Math.tan(Math.PI / 4.0 + b / 2.0)) - this.K)\r
 120                               + this.e\r
 121                               * Math.log(Math.tan(Math.PI / 4.0\r
 122                               + Math.asin(this.e * Math.sin(phy))\r
 123                               / 2.0));\r
 124                       prevPhy = phy;\r
 125                       phy = 2.0 * Math.atan(Math.exp(S)) - Math.PI / 2.0;\r
 126                     }\r
 127                 \r
 128                     p.x = lambda;\r
 129                     p.y = phy;\r
 130                     return p;\r
 131                   }\r
 132                 \r
 133                 \r
 134         }\r
 135 }