Understand shapefiles in different projections (notably OSGB)
[potlatch2.git] / 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 }