Commit dac9fecd authored by Calvin Metcalf's avatar Calvin Metcalf
Browse files

somerc

parent 505f0e77

src/projCode/somerc.js

+89 −110
Original line number Diff line number Diff line
@@ -27,21 +27,16 @@ proj4.Proj.somerc = { 
    this.R = this.k0 * semiMajorAxis * Math.sqrt(1 - e2) / (1 - e2 * Math.pow(sinPhy0, 2.0));

    this.alpha = Math.sqrt(1 + e2 / (1 - e2) * Math.pow(Math.cos(phy0), 4.0));

    this.b0 = Math.asin(sinPhy0 / this.alpha);

    this.K = Math.log(Math.tan(Math.PI / 4.0 + this.b0 / 2.0))

            - this.alpha

            * Math.log(Math.tan(Math.PI / 4.0 + phy0 / 2.0))

            + this.alpha

            * e / 2

            * Math.log((1 + e * sinPhy0)

            / (1 - e * sinPhy0));

    var k1 = Math.log(Math.tan(Math.PI / 4.0 + this.b0 / 2.0));

    var k2 = Math.log(Math.tan(Math.PI / 4.0 + phy0 / 2.0));

    var k3 = Math.log((1 + e * sinPhy0) / (1 - e * sinPhy0));

    this.K = k1 - this.alpha * k2 + this.alpha * e / 2 * k3;

  },





  forward: function(p) {

    var Sa1 = Math.log(Math.tan(Math.PI / 4.0 - p.y / 2.0));

    var Sa2 = this.e / 2.0

            * Math.log((1 + this.e * Math.sin(p.y))

            / (1 - this.e * Math.sin(p.y)));

    var Sa2 = this.e / 2.0 * Math.log((1 + this.e * Math.sin(p.y)) / (1 - this.e * Math.sin(p.y)));

    var S = -this.alpha * (Sa1 + Sa2) + this.K;



    // spheric latitude
@@ -51,16 +46,11 @@ proj4.Proj.somerc = { 
    var I = this.alpha * (p.x - this.lambda0);



    // psoeudo equatorial rotation

    var rotI = Math.atan(Math.sin(I)

            / (Math.sin(this.b0) * Math.tan(b) +

               Math.cos(this.b0) * Math.cos(I)));

    var rotI = Math.atan(Math.sin(I) / (Math.sin(this.b0) * Math.tan(b) + Math.cos(this.b0) * Math.cos(I)));



    var rotB = Math.asin(Math.cos(this.b0) * Math.sin(b) -

                         Math.sin(this.b0) * Math.cos(b) * Math.cos(I));

    var rotB = Math.asin(Math.cos(this.b0) * Math.sin(b) - Math.sin(this.b0) * Math.cos(b) * Math.cos(I));



    p.y = this.R / 2.0

            * Math.log((1 + Math.sin(rotB)) / (1 - Math.sin(rotB)))

            + this.y0;

    p.y = this.R / 2.0 * Math.log((1 + Math.sin(rotB)) / (1 - Math.sin(rotB))) + this.y0;

    p.x = this.R * rotI + this.x0;

    return p;

  },
@@ -72,11 +62,8 @@ proj4.Proj.somerc = { 
    var rotI = Y / this.R;

    var rotB = 2 * (Math.atan(Math.exp(X / this.R)) - Math.PI / 4.0);



    var b = Math.asin(Math.cos(this.b0) * Math.sin(rotB)

            + Math.sin(this.b0) * Math.cos(rotB) * Math.cos(rotI));

    var I = Math.atan(Math.sin(rotI)

            / (Math.cos(this.b0) * Math.cos(rotI) - Math.sin(this.b0)

            * Math.tan(rotB)));

    var b = Math.asin(Math.cos(this.b0) * Math.sin(rotB) + Math.sin(this.b0) * Math.cos(rotB) * Math.cos(rotI));

    var I = Math.atan(Math.sin(rotI) / (Math.cos(this.b0) * Math.cos(rotI) - Math.sin(this.b0) * Math.tan(rotB)));



    var lambda = this.lambda0 + I / this.alpha;
@@ -84,21 +71,13 @@ proj4.Proj.somerc = { 
    var phy = b;

    var prevPhy = -1000.0;

    var iteration = 0;

    while (Math.abs(phy - prevPhy) > 0.0000001)

    {

      if (++iteration > 20)

      {

    while (Math.abs(phy - prevPhy) > 0.0000001) {

      if (++iteration > 20) {

        proj4.reportError("omercFwdInfinity");

        return;

      }

      //S = Math.log(Math.tan(Math.PI / 4.0 + phy / 2.0));

      S = 1.0

              / this.alpha

              * (Math.log(Math.tan(Math.PI / 4.0 + b / 2.0)) - this.K)

              + this.e

              * Math.log(Math.tan(Math.PI / 4.0

              + Math.asin(this.e * Math.sin(phy))

              / 2.0));

      S = 1.0 / this.alpha * (Math.log(Math.tan(Math.PI / 4.0 + b / 2.0)) - this.K) + this.e * Math.log(Math.tan(Math.PI / 4.0 + Math.asin(this.e * Math.sin(phy)) / 2.0));

      prevPhy = phy;

      phy = 2.0 * Math.atan(Math.exp(S)) - Math.PI / 2.0;

    }