com/gradoservice/proj4as/proj/ProjAea.as

   1 package com.gradoservice.proj4as.proj
   2 {
   3         import com.gradoservice.proj4as.ProjConstants;
   4         import com.gradoservice.proj4as.ProjPoint;
   5         import com.gradoservice.proj4as.Datum;
   6
   7         public class ProjAea extends AbstractProjProjection
   8         {
   9                 public function ProjAea(data:ProjParams)
  10                 {
  11                         super(data);
  12                 }
  13
  14   override public function init():void
  15         {
  16     if (Math.abs(this.lat1 + this.lat2) <ProjConstants.EPSLN) {
  17        trace("aeaInitEqualLatitudes");
  18        return;
  19     }
  20     this.temp = this.b / this.a;
  21     this.es = 1.0 - Math.pow(this.temp,2);
  22     this.e3 = Math.sqrt(this.es);
  23
  24     this.sin_po=Math.sin(this.lat1);
  25     this.cos_po=Math.cos(this.lat1);
  26     this.t1=this.sin_po
  27     this.con = this.sin_po;
  28     this.ms1 = ProjConstants.msfnz(this.e3,this.sin_po,this.cos_po);
  29     this.qs1 = ProjConstants.qsfnz(this.e3,this.sin_po,this.cos_po);
  30
  31     this.sin_po=Math.sin(this.lat2);
  32     this.cos_po=Math.cos(this.lat2);
  33     this.t2=this.sin_po;
  34     this.ms2 = ProjConstants.msfnz(this.e3,this.sin_po,this.cos_po);
  35     this.qs2 = ProjConstants.qsfnz(this.e3,this.sin_po,this.cos_po);
  36
  37     this.sin_po=Math.sin(this.lat0);
  38     this.cos_po=Math.cos(this.lat0);
  39     this.t3=this.sin_po;
  40     this.qs0 = ProjConstants.qsfnz(this.e3,this.sin_po,this.cos_po);
  41
  42     if (Math.abs(this.lat1 - this.lat2) > ProjConstants.EPSLN) {
  43       this.ns0 = (this.ms1 * this.ms1 - this.ms2 *this.ms2)/ (this.qs2 - this.qs1);
  44     } else {
  45       this.ns0 = this.con;
  46     }
  47     this.c = this.ms1 * this.ms1 + this.ns0 * this.qs1;
  48     this.rh = this.a * Math.sqrt(this.c - this.ns0 * this.qs0)/this.ns0;
  49   }
  50
  51 /* Albers Conical Equal Area forward equations--mapping lat,long to x,y
  52   -------------------------------------------------------------------*/
  53   override public function forward(p:ProjPoint):ProjPoint
  54   {
  55     var lon:Number=p.x;
  56     var lat:Number=p.y;
  57
  58     this.sin_phi=Math.sin(lat);
  59     this.cos_phi=Math.cos(lat);
  60
  61     var qs:Number = ProjConstants.qsfnz(this.e3,this.sin_phi,this.cos_phi);
  62     var rh1:Number =this.a * Math.sqrt(this.c - this.ns0 * qs)/this.ns0;
  63     var theta:Number = this.ns0 * ProjConstants.adjust_lon(lon - this.long0);
  64     var x:Number = rh1 * Math.sin(theta) + this.x0;
  65     var y:Number = this.rh - rh1 * Math.cos(theta) + this.y0;
  66
  67     p.x = x;
  68     p.y = y;
  69     return p;
  70   }
  71
  72
  73   override public function inverse(p:ProjPoint):ProjPoint
  74   {
  75     var rh1:Number;
  76     var qs:Number;
  77     var con:Number;
  78     var theta:Number
  79     var lon:Number;
  80     var lat:Number;
  81
  82     p.x -= this.x0;
  83     p.y = this.rh - p.y + this.y0;
  84     if (this.ns0 >= 0) {
  85       rh1 = Math.sqrt(p.x *p.x + p.y * p.y);
  86       con = 1.0;
  87     } else {
  88       rh1 = -Math.sqrt(p.x * p.x + p.y *p.y);
  89       con = -1.0;
  90     }
  91     theta = 0.0;
  92     if (rh1 != 0.0) {
  93       theta = Math.atan2(con * p.x, con * p.y);
  94     }
  95     con = rh1 * this.ns0 / this.a;
  96     qs = (this.c - con * con) / this.ns0;
  97     if (this.e3 >= 1e-10) {
  98       con = 1 - .5 * (1.0 -this.es) * Math.log((1.0 - this.e3) / (1.0 + this.e3))/this.e3;
  99       if (Math.abs(Math.abs(con) - Math.abs(qs)) > .0000000001 ) {
 100           lat = this.phi1z(this.e3,qs);
 101       } else {
 102           if (qs >= 0) {
 103              lat = .5 * Math.PI;
 104           } else {
 105              lat = -.5 * Math.PI;
 106           }
 107       }
 108     } else {
 109       lat = this.phi1z(e3,qs);
 110     }
 111
 112     lon = ProjConstants.adjust_lon(theta/this.ns0 + this.long0);
 113     p.x = lon;
 114     p.y = lat;
 115     return p;
 116   }
 117
 118 /* Function to compute phi1, the latitude for the inverse of the
 119    Albers Conical Equal-Area projection.
 120 -------------------------------------------*/
 121    private function phi1z (eccent:Number,qs:Number):Number
 122    {
 123     var con:Number;
 124     var com:Number
 125     var dphi:Number;
 126     var phi:Number = ProjConstants.asinz(.5 * qs);
 127     if (eccent < ProjConstants.EPSLN) return phi;
 128
 129     var eccnts:Number = eccent * eccent;
 130     for (var i:int = 1; i <= 25; i++) {
 131         sinphi = Math.sin(phi);
 132         cosphi = Math.cos(phi);
 133         con = eccent * sinphi;
 134         com = 1.0 - con * con;
 135         dphi = .5 * com * com / cosphi * (qs / (1.0 - eccnts) - sinphi / com + .5 / eccent * Math.log((1.0 - con) / (1.0 + con)));
 136         phi = phi + dphi;
 137         if (Math.abs(dphi) <= 1e-7) return phi;
 138     }
 139     trace("aea:phi1z:Convergence error");
 140     return 0;
 141   }
 142
 143
 144         }
 145 }