com/gradoservice/proj4as/proj/ProjStere.as

   1 package com.gradoservice.proj4as.proj\r
   2 {\r
   3         import com.gradoservice.proj4as.ProjPoint;\r
   4         import com.gradoservice.proj4as.ProjConstants;\r
   5         import com.gradoservice.proj4as.Datum;\r
   6                 \r
   7         public class ProjStere extends AbstractProjProjection\r
   8         {\r
   9                 private const TOL:Number=1.e-8;\r
  10                 private const NITER:Number=8;\r
  11                 private const CONV:Number=1.e-10;\r
  12                 private const S_POLE:Number=0;\r
  13                 private const N_POLE:Number=1;\r
  14                 private const OBLIQ:Number=2;\r
  15                 private const EQUIT:Number=3;\r
  16                 \r
  17                 private var akm1:Number;\r
  18                 private var cosph0:Number;\r
  19                 private var cosX1:Number;\r
  20                 private var sinX1:Number;               \r
  21                 private var phi0:Number;\r
  22                 private var phits:Number;\r
  23                 private var sinph0:Number;\r
  24                 \r
  25                 public function ProjStere(data:ProjParams)\r
  26                 {\r
  27                         super(data);\r
  28                 }\r
  29                 \r
  30                 private function ssfn_(phit:Number, sinphi:Number, eccen:Number):Number \r
  31                 {\r
  32                         sinphi *= eccen;\r
  33                         return (Math.tan (.5 * (ProjConstants.HALF_PI + phit)) * Math.pow((1. - sinphi) / (1. + sinphi), .5 * eccen));\r
  34                 }\r
  35                 \r
  36                 \r
  37                   override public function init():void\r
  38                   {\r
  39                         this.phits = this.lat_ts ? this.lat_ts : ProjConstants.HALF_PI;\r
  40                     var t:Number = Math.abs(this.lat0);\r
  41                         if ((Math.abs(t) - ProjConstants.HALF_PI) < ProjConstants.EPSLN) {\r
  42                                 this.mode = this.lat0 < 0. ? this.S_POLE : this.N_POLE;\r
  43                         } else {\r
  44                                 this.mode = t > ProjConstants.EPSLN ? this.OBLIQ : this.EQUIT;\r
  45                     }\r
  46                         this.phits = Math.abs(this.phits);\r
  47                         if (this.es) {\r
  48                                 var X:Number;\r
  49                 \r
  50                                 switch (this.mode) {\r
  51                                 case this.N_POLE:\r
  52                                 case this.S_POLE:\r
  53                                         if (Math.abs(this.phits - ProjConstants.HALF_PI) < ProjConstants.EPSLN) {\r
  54                                                 this.akm1 = 2. * this.k0 / Math.sqrt(Math.pow(1+this.e,1+this.e)*Math.pow(1-this.e,1-this.e));\r
  55                                         } else {\r
  56                           t = Math.sin(this.phits);\r
  57                                                 this.akm1 = Math.cos(this.phits) / ProjConstants.tsfnz(this.e, this.phits, t);\r
  58                                                 t *= this.e;\r
  59                                                 this.akm1 /= Math.sqrt(1. - t * t);\r
  60                                         }\r
  61                                         break;\r
  62                                 case this.EQUIT:\r
  63                                         this.akm1 = 2. * this.k0;\r
  64                                         break;\r
  65                                 case this.OBLIQ:\r
  66                                         t = Math.sin(this.lat0);\r
  67                                         X = 2. * Math.atan(this.ssfn_(this.lat0, t, this.e)) - ProjConstants.HALF_PI;\r
  68                                         t *= this.e;\r
  69                                         this.akm1 = 2. * this.k0 * Math.cos(this.lat0) / Math.sqrt(1. - t * t);\r
  70                                         this.sinX1 = Math.sin(X);\r
  71                                         this.cosX1 = Math.cos(X);\r
  72                                         break;\r
  73                                 }\r
  74                         } else {\r
  75                                 switch (this.mode) {\r
  76                                 case this.OBLIQ:\r
  77                                         this.sinph0 = Math.sin(this.lat0);\r
  78                                         this.cosph0 = Math.cos(this.lat0);\r
  79                                 case this.EQUIT:\r
  80                                         this.akm1 = 2. * this.k0;\r
  81                                         break;\r
  82                                 case this.S_POLE:\r
  83                                 case this.N_POLE:\r
  84                                         this.akm1 = Math.abs(this.phits - ProjConstants.HALF_PI) >= ProjConstants.EPSLN ?\r
  85                                            Math.cos(this.phits) / Math.tan(ProjConstants.FORTPI - .5 * this.phits) :\r
  86                                            2. * this.k0 ;\r
  87                                         break;\r
  88                                 }\r
  89                         }\r
  90                   } \r
  91                 \r
  92                 // Stereographic forward equations--mapping lat,long to x,y\r
  93                   override public function forward(p:ProjPoint):ProjPoint\r
  94                   {\r
  95                     var lon:Number = p.x;\r
  96                     var lat:Number = p.y;\r
  97                     var x:Number, y:Number, A:Number, X:Number;\r
  98                     var sinX:Number,cosX:Number;\r
  99                     \r
 100                     if (this.sphere) {\r
 101                         var  sinphi:Number, cosphi:Number, coslam:Number, sinlam:Number;\r
 102                 \r
 103                         sinphi = Math.sin(lat);\r
 104                         cosphi = Math.cos(lat);\r
 105                         coslam = Math.cos(lon);\r
 106                         sinlam = Math.sin(lon);\r
 107                         switch (this.mode) {\r
 108                         case this.EQUIT:\r
 109                                 y = 1. + cosphi * coslam;\r
 110                                 if (y <= ProjConstants.EPSLN) {\r
 111                           trace('ERROR');\r
 112                           //F_ERROR;\r
 113                         }\r
 114                         y = this.akm1 / y;\r
 115                                 x = y * cosphi * sinlam;\r
 116                         y *= sinphi;\r
 117                                 break;\r
 118                         case this.OBLIQ:\r
 119                                 y = 1. + this.sinph0 * sinphi + this.cosph0 * cosphi * coslam;\r
 120                                 if (y <= ProjConstants.EPSLN) {\r
 121                           trace('ERROR');\r
 122                           //F_ERROR;\r
 123                         }\r
 124                         y = this.akm1 / y;\r
 125                                 x = y * cosphi * sinlam;\r
 126                                 y *= this.cosph0 * sinphi - this.sinph0 * cosphi * coslam;\r
 127                                 break;\r
 128                         case this.N_POLE:\r
 129                                 coslam = -coslam;\r
 130                                 lat = -lat;\r
 131                         //Note  no break here so it conitnues through S_POLE\r
 132                         case this.S_POLE:\r
 133                                 if (Math.abs(lat - ProjConstants.HALF_PI) < this.TOL) {\r
 134                           trace('ERROR');\r
 135                           //F_ERROR;\r
 136                         }\r
 137                         y = this.akm1 * Math.tan(ProjConstants.FORTPI + .5 * lat)\r
 138                                 x = sinlam * y;\r
 139                                 y *= coslam;\r
 140                                 break;\r
 141                         }\r
 142                     } else {\r
 143                         coslam = Math.cos(lon);\r
 144                         sinlam = Math.sin(lon);\r
 145                         sinphi = Math.sin(lat);\r
 146                         if (this.mode == this.OBLIQ || this.mode == this.EQUIT) {\r
 147                         X = 2. * Math.atan(this.ssfn_(lat, sinphi, this.e));\r
 148                                 sinX = Math.sin(X - ProjConstants.HALF_PI);\r
 149                                 cosX = Math.cos(X);\r
 150                         }\r
 151                         switch (this.mode) {\r
 152                         case this.OBLIQ:\r
 153                                 A = this.akm1 / (this.cosX1 * (1. + this.sinX1 * sinX + this.cosX1 * cosX * coslam));\r
 154                                 y = A * (this.cosX1 * sinX - this.sinX1 * cosX * coslam);\r
 155                                 x = A * cosX;\r
 156                                 break;\r
 157                         case this.EQUIT:\r
 158                                 A = 2. * this.akm1 / (1. + cosX * coslam);\r
 159                                 y = A * sinX;\r
 160                                 x = A * cosX;\r
 161                                 break;\r
 162                         case this.S_POLE:\r
 163                                 lat = -lat;\r
 164                                 coslam = - coslam;\r
 165                                 sinphi = -sinphi;\r
 166                         case this.N_POLE:\r
 167                                 x = this.akm1 * ProjConstants.tsfnz(this.e, lat, sinphi);\r
 168                                 y = - x * coslam;\r
 169                                 break;\r
 170                         }\r
 171                         x = x * sinlam;\r
 172                     }\r
 173                     p.x = x*this.a + this.x0;\r
 174                     p.y = y*this.a + this.y0;\r
 175                     return p;\r
 176                   }\r
 177                 \r
 178                 \r
 179                 //* Stereographic inverse equations--mapping x,y to lat/long\r
 180                  override public function inverse(p:ProjPoint):ProjPoint\r
 181                  {\r
 182                     var x:Number = (p.x - this.x0)/this.a;   /* descale and de-offset */\r
 183                     var y:Number = (p.y - this.y0)/this.a;\r
 184                     var lon:Number, lat:Number;\r
 185                 \r
 186                     var cosphi:Number, sinphi:Number, tp:Number=0.0, phi_l:Number=0.0, rho:Number, halfe:Number=0.0, pi2:Number=0.0;\r
 187                     var i:int;\r
 188                 \r
 189                     if (this.sphere) {\r
 190                         var  c:Number, rh:Number, sinc:Number, cosc:Number;\r
 191                 \r
 192                       rh = Math.sqrt(x*x + y*y);\r
 193                       c = 2. * Math.atan(rh / this.akm1);\r
 194                         sinc = Math.sin(c);\r
 195                         cosc = Math.cos(c);\r
 196                         lon = 0.;\r
 197                         switch (this.mode) {\r
 198                         case this.EQUIT:\r
 199                                 if (Math.abs(rh) <= ProjConstants.EPSLN) {\r
 200                                         lat = 0.;\r
 201                                 } else {\r
 202                                         lat = Math.asin(y * sinc / rh);\r
 203                         }\r
 204                                 if (cosc != 0. || x != 0.) lon = Math.atan2(x * sinc, cosc * rh);\r
 205                                 break;\r
 206                         case this.OBLIQ:\r
 207                                 if (Math.abs(rh) <= ProjConstants.EPSLN) {\r
 208                                         lat = this.phi0;\r
 209                                 } else {\r
 210                                         lat = Math.asin(cosc * sinph0 + y * sinc * cosph0 / rh);\r
 211                         }\r
 212                         c = cosc - sinph0 * Math.sin(lat);\r
 213                                 if (c != 0. || x != 0.) {\r
 214                                         lon = Math.atan2(x * sinc * cosph0, c * rh);\r
 215                         }\r
 216                                 break;\r
 217                         case this.N_POLE:\r
 218                                 y = -y;\r
 219                         case this.S_POLE:\r
 220                                 if (Math.abs(rh) <= ProjConstants.EPSLN) {\r
 221                                         lat = this.phi0;\r
 222                                 } else {\r
 223                                         lat = Math.asin(this.mode == this.S_POLE ? -cosc : cosc);\r
 224                         }\r
 225                                 lon = (x == 0. && y == 0.) ? 0. : Math.atan2(x, y);\r
 226                                 break;\r
 227                         }\r
 228                     } else {\r
 229                         rho = Math.sqrt(x*x + y*y);\r
 230                         switch (this.mode) {\r
 231                         case this.OBLIQ:\r
 232                         case this.EQUIT:\r
 233                         tp = 2. * Math.atan2(rho * this.cosX1 , this.akm1);\r
 234                                 cosphi = Math.cos(tp);\r
 235                                 sinphi = Math.sin(tp);\r
 236                         if( rho == 0.0 ) {\r
 237                                   phi_l = Math.asin(cosphi * this.sinX1);\r
 238                         } else {\r
 239                                   phi_l = Math.asin(cosphi * this.sinX1 + (y * sinphi * this.cosX1 / rho));\r
 240                         }\r
 241                 \r
 242                                 tp = Math.tan(.5 * (ProjConstants.HALF_PI + phi_l));\r
 243                                 x *= sinphi;\r
 244                                 y = rho * this.cosX1 * cosphi - y * this.sinX1* sinphi;\r
 245                                 pi2 = ProjConstants.HALF_PI;\r
 246                                 halfe = .5 * this.e;\r
 247                                 break;\r
 248                         case this.N_POLE:\r
 249                                 y = -y;\r
 250                         case this.S_POLE:\r
 251                         tp = - rho / this.akm1\r
 252                                 phi_l = ProjConstants.HALF_PI - 2. * Math.atan(tp);\r
 253                                 pi2 = -ProjConstants.HALF_PI;\r
 254                                 halfe = -.5 * this.e;\r
 255                                 break;\r
 256                         }\r
 257                         for (i = this.NITER; i--; phi_l = lat) \r
 258                                 { //check this\r
 259                                 sinphi = this.e * Math.sin(phi_l);\r
 260                                 lat = 2. * Math.atan(tp * Math.pow((1.+sinphi)/(1.-sinphi), halfe)) - pi2;\r
 261                                 if (Math.abs(phi_l - lat) < this.CONV) \r
 262                                         {\r
 263                                         if (this.mode == this.S_POLE) lat = -lat;\r
 264                                         lon = (x == 0. && y == 0.) ? 0. : Math.atan2(x, y);\r
 265                                         p.x = lon;\r
 266                                         p.y = lat\r
 267                                         return p;\r
 268                                         }\r
 269                                 }\r
 270                         }\r
 271                     return null;\r
 272                         }\r
 273                 \r
 274 \r
 275                 \r
 276                 \r
 277         }\r
 278 }