com/gradoservice/proj4as/proj/ProjOmerc.as

   1 /*******************************************************************************\r
   2 NAME                       OBLIQUE MERCATOR (HOTINE) \r
   3 \r
   4 PURPOSE:        Transforms input longitude and latitude to Easting and\r
   5                 Northing for the Oblique 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 PROGRAMMER              DATE\r
  10 ----------              ----\r
  11 T. Mittan               Mar, 1993\r
  12 \r
  13 ALGORITHM REFERENCES\r
  14 \r
  15 1.  Snyder, John P., "Map Projections--A Working Manual", U.S. Geological\r
  16     Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United\r
  17     State Government Printing Office, Washington D.C., 1987.\r
  18 \r
  19 2.  Snyder, John P. and Voxland, Philip M., "An Album of Map Projections",\r
  20     U.S. Geological Survey Professional Paper 1453 , United State Government\r
  21     Printing Office, Washington D.C., 1989.\r
  22 *******************************************************************************/\r
  23 \r
  24 package com.gradoservice.proj4as.proj\r
  25 {\r
  26         import com.gradoservice.proj4as.ProjPoint;\r
  27         import com.gradoservice.proj4as.ProjConstants;\r
  28         import com.gradoservice.proj4as.Datum;\r
  29                 \r
  30         public class ProjOmerc extends AbstractProjProjection\r
  31         {\r
  32                 private var d:Number;           \r
  33                 private var f:Number;           \r
  34                 private var h:Number;           \r
  35                 private var j:Number;           \r
  36                 private var l:Number;   \r
  37                 private var p:Number;   \r
  38                 private var u:Number;           \r
  39                 private var dlon:Number;                \r
  40                 private var al:Number;\r
  41                 private var bl:Number;          \r
  42                 private var el:Number;          \r
  43                 private var com:Number;         \r
  44                 private var cos_p20:Number;             \r
  45                 private var cosaz:Number;               \r
  46                 private var cosgam:Number;              \r
  47                 private var at1:Number;         \r
  48                 private var gam:Number;         \r
  49                 private var gama:Number;\r
  50                 private var lon1:Number;\r
  51                 private var lon2:Number;\r
  52                 private var sin_p20:Number;                                             \r
  53                 private var sinaz:Number;       \r
  54                 private var singam:Number;              \r
  55                 private var ts:Number;          \r
  56                 private var ts1:Number;         \r
  57                 private var ts2:Number;         \r
  58                 \r
  59                 public function ProjOmerc(data:ProjParams)\r
  60                 {\r
  61                         super(data);\r
  62                 }\r
  63                 \r
  64                 \r
  65                   override public function init():void\r
  66                   {\r
  67                     if (!this.mode) this.mode=0;\r
  68                     if (!this.lon1)   {this.lon1=0;this.mode=1;}\r
  69                     if (!this.lon2)   this.lon2=0;\r
  70                     if (!this.lat2)    this.lat2=0;\r
  71                 \r
  72                     /* Place parameters in static storage for common use\r
  73                       -------------------------------------------------*/\r
  74                     var temp:Number = this.b/ this.a;\r
  75                     var es:Number = 1.0 - Math.pow(temp,2);\r
  76                     var e:Number = Math.sqrt(es);\r
  77                 \r
  78                     this.sin_p20=Math.sin(this.lat0);\r
  79                     this.cos_p20=Math.cos(this.lat0);\r
  80                 \r
  81                     this.con = 1.0 - this.es * this.sin_p20 * this.sin_p20;\r
  82                     this.com = Math.sqrt(1.0 - es);\r
  83                     this.bl = Math.sqrt(1.0 + this.es * Math.pow(this.cos_p20,4.0)/(1.0 - es));\r
  84                     this.al = this.a * this.bl * this.k0 * this.com / this.con;\r
  85                     if (Math.abs(this.lat0) < ProjConstants.EPSLN) {\r
  86                        this.ts = 1.0;\r
  87                        this.d = 1.0;\r
  88                        this.el = 1.0;\r
  89                     } else {\r
  90                        this.ts = ProjConstants.tsfnz(this.e,this.lat0,this.sin_p20);\r
  91                        this.con = Math.sqrt(this.con);\r
  92                        this.d = this.bl * this.com / (this.cos_p20 * this.con);\r
  93                        if ((this.d * this.d - 1.0) > 0.0) {\r
  94                           if (this.lat0 >= 0.0) {\r
  95                              this.f = this.d + Math.sqrt(this.d * this.d - 1.0);\r
  96                           } else {\r
  97                              this.f = this.d - Math.sqrt(this.d * this.d - 1.0);\r
  98                           }\r
  99                        } else {\r
 100                          this.f = this.d;\r
 101                        }\r
 102                        this.el = this.f * Math.pow(this.ts,this.bl);\r
 103                     }\r
 104                 \r
 105                     //this.longc=52.60353916666667;\r
 106                 \r
 107                     if (this.mode != 0) {\r
 108                        this.g = .5 * (this.f - 1.0/this.f);\r
 109                        this.gama = ProjConstants.asinz(Math.sin(this.alpha) / this.d);\r
 110                        this.longc= this.longc - ProjConstants.asinz(this.g * Math.tan(this.gama))/this.bl;\r
 111                 \r
 112                        /* Report parameters common to format B\r
 113                        -------------------------------------*/\r
 114                        //genrpt(azimuth * R2D,"Azimuth of Central Line:    ");\r
 115                        //cenlon(lon_origin);\r
 116                       // cenlat(lat_origin);\r
 117                 \r
 118                        this.con = Math.abs(this.lat0);\r
 119                        if ((this.con > ProjConstants.EPSLN) && (Math.abs(this.con - ProjConstants.HALF_PI) > ProjConstants.EPSLN)) {\r
 120                             this.singam=Math.sin(this.gama);\r
 121                             this.cosgam=Math.cos(this.gama);\r
 122                 \r
 123                             this.sinaz=Math.sin(this.alpha);\r
 124                             this.cosaz=Math.cos(this.alpha);\r
 125                 \r
 126                             if (this.lat0>= 0) {\r
 127                                this.u =  (this.al / this.bl) * Math.atan(Math.sqrt(this.d*this.d - 1.0)/this.cosaz);\r
 128                             } else {\r
 129                                this.u =  -(this.al / this.bl) *Math.atan(Math.sqrt(this.d*this.d - 1.0)/this.cosaz);\r
 130                             }\r
 131                           } else {\r
 132                             trace("omerc:Init:DataError");\r
 133                           }\r
 134                        } else {\r
 135                        this.sinphi =Math. sin(this.at1);\r
 136                        this.ts1 = ProjConstants.tsfnz(this.e,this.lat1,this.sinphi);\r
 137                        this.sinphi = Math.sin(this.lat2);\r
 138                        this.ts2 = ProjConstants.tsfnz(this.e,this.lat2,this.sinphi);\r
 139                        this.h = Math.pow(this.ts1,this.bl);\r
 140                        this.l = Math.pow(this.ts2,this.bl);\r
 141                        this.f = this.el/this.h;\r
 142                        this.g = .5 * (this.f - 1.0/this.f);\r
 143                        this.j = (this.el * this.el - this.l * this.h)/(this.el * this.el + this.l * this.h);\r
 144                        this.p = (this.l - this.h) / (this.l + this.h);\r
 145                        this.dlon = this.lon1 - this.lon2;\r
 146                        if (this.dlon < -ProjConstants.PI) this.lon2 = this.lon2 - 2.0 * ProjConstants.PI;\r
 147                        if (this.dlon > ProjConstants.PI) this.lon2 = this.lon2 + 2.0 * ProjConstants.PI;\r
 148                        this.dlon = this.lon1 - this.lon2;\r
 149                        this.longc = .5 * (this.lon1 + this.lon2) -Math.atan(this.j * Math.tan(.5 * this.bl * this.dlon)/this.p)/this.bl;\r
 150                        this.dlon  = ProjConstants.adjust_lon(this.lon1 - this.longc);\r
 151                        this.gama = Math.atan(Math.sin(this.bl * this.dlon)/this.g);\r
 152                        this.alpha = ProjConstants.asinz(this.d * Math.sin(this.gama));\r
 153                 \r
 154                        /* Report parameters common to format A\r
 155                        -------------------------------------*/\r
 156                 \r
 157                        if (Math.abs(this.lat1 - this.lat2) <= ProjConstants.EPSLN) {\r
 158                           trace("omercInitDataError");\r
 159                           //return(202);\r
 160                        } else {\r
 161                           this.con = Math.abs(this.lat1);\r
 162                        }\r
 163                        if ((this.con <= ProjConstants.EPSLN) || (Math.abs(this.con - ProjConstants.HALF_PI) <= ProjConstants.EPSLN)) {\r
 164                            trace("omercInitDataError");\r
 165                                 //return(202);\r
 166                        } else {\r
 167                          if (Math.abs(Math.abs(this.lat0) - ProjConstants.HALF_PI) <= ProjConstants.EPSLN) {\r
 168                             trace("omercInitDataError");\r
 169                             //return(202);\r
 170                          }\r
 171                        }\r
 172                 \r
 173                        this.singam=Math.sin(this.gam);\r
 174                        this.cosgam=Math.cos(this.gam);\r
 175                 \r
 176                        this.sinaz=Math.sin(this.alpha);\r
 177                        this.cosaz=Math.cos(this.alpha);  \r
 178                 \r
 179                 \r
 180                        if (this.lat0 >= 0) {\r
 181                           this.u =  (this.al/this.bl) * Math.atan(Math.sqrt(this.d * this.d - 1.0)/this.cosaz);\r
 182                        } else {\r
 183                           this.u = -(this.al/this.bl) * Math.atan(Math.sqrt(this.d * this.d - 1.0)/this.cosaz);\r
 184                        }\r
 185                      }\r
 186                   }\r
 187                 \r
 188                 \r
 189                   /* Oblique Mercator forward equations--mapping lat,long to x,y\r
 190                     ----------------------------------------------------------*/\r
 191                   override public function forward(p:ProjPoint):ProjPoint\r
 192                   {\r
 193                     var theta:Number;           /* angle                                        */\r
 194                     var sin_phi:Number, cos_phi:Number;/* sin and cos value                             */\r
 195                     var b:Number;               /* temporary values                             */\r
 196                     var c:Number, t:Number, tq:Number;  /* temporary values                             */\r
 197                     var con:Number, n:Number, ml:Number;        /* cone constant, small m                       */\r
 198                     var q:Number,us:Number,vl:Number;\r
 199                     var ul:Number,vs:Number;\r
 200                     var s:Number;\r
 201                     var dlon:Number;\r
 202                     var ts1:Number;\r
 203                 \r
 204                     var lon:Number=p.x;\r
 205                     var lat:Number=p.y;\r
 206                     /* Forward equations\r
 207                       -----------------*/\r
 208                     sin_phi = Math.sin(lat);\r
 209                     dlon = ProjConstants.adjust_lon(lon - this.longc);\r
 210                     vl = Math.sin(this.bl * dlon);\r
 211                     if (Math.abs(Math.abs(lat) - ProjConstants.HALF_PI) > ProjConstants.EPSLN) {\r
 212                        ts1 = ProjConstants.tsfnz(this.e,lat,sin_phi);\r
 213                        q = this.el / (Math.pow(ts1,this.bl));\r
 214                        s = .5 * (q - 1.0 / q);\r
 215                        t = .5 * (q + 1.0/ q);\r
 216                        ul = (s * this.singam - vl * this.cosgam) / t;\r
 217                        con = Math.cos(this.bl * dlon);\r
 218                        if (Math.abs(con) < .0000001) {\r
 219                           us = this.al * this.bl * dlon;\r
 220                        } else {\r
 221                           us = this.al * Math.atan((s * this.cosgam + vl * this.singam) / con)/this.bl;\r
 222                           if (con < 0) us = us + ProjConstants.PI * this.al / this.bl;\r
 223                        }\r
 224                     } else {\r
 225                        if (lat >= 0) {\r
 226                           ul = this.singam;\r
 227                        } else {\r
 228                           ul = -this.singam;\r
 229                        }\r
 230                        us = this.al * lat / this.bl;\r
 231                     }\r
 232                     if (Math.abs(Math.abs(ul) - 1.0) <= ProjConstants.EPSLN) {\r
 233                        //alert("Point projects into infinity","omer-for");\r
 234                        trace("omercFwdInfinity");\r
 235                        //return(205);\r
 236                     }\r
 237                     vs = .5 * this.al * Math.log((1.0 - ul)/(1.0 + ul)) / this.bl;\r
 238                     us = us - this.u;\r
 239                     var x:Number = this.x0 + vs * this.cosaz + us * this.sinaz;\r
 240                     var y:Number = this.y0 + us * this.cosaz - vs * this.sinaz;\r
 241                 \r
 242                     p.x=x;\r
 243                     p.y=y;\r
 244                     return p;\r
 245                   }\r
 246                 \r
 247                  override public function inverse(p:ProjPoint):ProjPoint\r
 248                  {\r
 249                     var delta_lon:Number;       /* Delta longitude (Given longitude - center    */\r
 250                     var theta:Number;           /* angle                                        */\r
 251                     var delta_theta:Number;     /* adjusted longitude                           */\r
 252                     var sin_phi:Number, cos_phi:Number;/* sin and cos value                             */\r
 253                     var b:Number;               /* temporary values                             */\r
 254                     var c:Number, t:Number, tq:Number;  /* temporary values                             */\r
 255                     var con:Number, n:Number, ml:Number;        /* cone constant, small m                       */\r
 256                     var vs:Number,us:Number,q:Number,s:Number,ts1:Number;\r
 257                     var vl:Number,ul:Number,bs:Number;\r
 258                     var dlon:Number;\r
 259                     var  flag:Number;\r
 260                 \r
 261                     /* Inverse equations\r
 262                       -----------------*/\r
 263                     p.x -= this.x0;\r
 264                     p.y -= this.y0;\r
 265                     flag = 0;\r
 266                     vs = p.x * this.cosaz - p.y * this.sinaz;\r
 267                     us = p.y * this.cosaz + p.x * this.sinaz;\r
 268                     us = us + this.u;\r
 269                     q = Math.exp(-this.bl * vs / this.al);\r
 270                     s = .5 * (q - 1.0/q);\r
 271                     t = .5 * (q + 1.0/q);\r
 272                     vl = Math.sin(this.bl * us / this.al);\r
 273                     ul = (vl * this.cosgam + s * this.singam)/t;\r
 274                     if (Math.abs(Math.abs(ul) - 1.0) <= ProjConstants.EPSLN)\r
 275                        {\r
 276                        var lon:Number = this.longc;\r
 277                        if (ul >= 0.0) {\r
 278                           var lat:Number = ProjConstants.HALF_PI;\r
 279                        } else {\r
 280                          lat = -ProjConstants.HALF_PI;\r
 281                        }\r
 282                     } else {\r
 283                        con = 1.0 / this.bl;\r
 284                        ts1 =Math.pow((this.el / Math.sqrt((1.0 + ul) / (1.0 - ul))),con);\r
 285                        lat = ProjConstants.phi2z(this.e,ts1);\r
 286                        //if (flag != 0)\r
 287                           //return(flag);\r
 288                        //~ con = Math.cos(this.bl * us /al);\r
 289                        theta = this.longc - Math.atan2((s * this.cosgam - vl * this.singam) , con)/this.bl;\r
 290                        lon = ProjConstants.adjust_lon(theta);\r
 291                     }\r
 292                     p.x=lon;\r
 293                     p.y=lat;\r
 294                     return p;\r
 295                   }\r
 296                 \r
 297                 \r
 298         }\r
 299 }