Commit df349606 authored by Michael Adair's avatar Michael Adair
Browse files

closes #66: change implementation to match PROJ.4 code for spherical and...

closes #66: change implementation to match PROJ.4 code for spherical and non-spherical forms for sinu projection
parent a4641a0d
Loading
Loading
Loading
Loading
+118 −78
Original line number Diff line number Diff line
@@ -1215,7 +1215,61 @@ Proj4js.common = {
  {
    var temp= e*sinphi;
    return a/Math.sqrt(1.0 - temp*temp);
  },
  
  //code from the PROJ.4 pj_mlfn.c file;  this may be useful for other projections
  pj_enfn: function(es) {
    var en = new Array();
    en[0] = this.C00 - es * (this.C02 + es * (this.C04 + es * (this.C06 + es * this.C08)));
    en[1] = es * (this.C22 - es * (this.C04 + es * (this.C06 + es * this.C08)));
    var t = es * es;
    en[2] = t * (this.C44 - es * (this.C46 + es * this.C48));
    t *= es;
    en[3] = t * (this.C66 - es * this.C68);
    en[4] = t * es * this.C88;
    return en;
  },
  
  pj_mlfn: function(phi, sphi, cphi, en) {
    cphi *= sphi;
    sphi *= sphi;
    return(en[0] * phi - cphi * (en[1] + sphi*(en[2]+ sphi*(en[3] + sphi*en[4]))));
  },
  
  pj_inv_mlfn: function(arg, es, en) {
    var k = 1./(1.-es);
    var phi = arg;
    for (var i = Proj4js.common.MAX_ITER; i ; --i) { /* rarely goes over 2 iterations */
      var s = Math.sin(phi);
      t = 1. - es * s * s;
      //t = this.pj_mlfn(phi, s, Math.cos(phi), en) - arg;
      //phi -= t * (t * Math.sqrt(t)) * k;
      t = (this.pj_mlfn(phi, s, Math.cos(phi), en) - arg) * (t * Math.sqrt(t)) * k;
      phi -= t;
      if (Math.abs(t) < Proj4js.common.EPSLN)
        return phi;
    }
    Proj4js.reportError("cass:pj_inv_mlfn: Convergence error");
    return phi;
  },

/* meridinal distance for ellipsoid and inverse
**	8th degree - accurate to < 1e-5 meters when used in conjuction
**		with typical major axis values.
**	Inverse determines phi to EPS (1e-11) radians, about 1e-6 seconds.
*/
  C00: 1.0,
  C02: .25,
  C04: .046875,
  C06: .01953125,
  C08: .01068115234375,
  C22: .75,
  C44: .46875,
  C46: .01302083333333333333,
  C48: .00712076822916666666,
  C66: .36458333333333333333,
  C68: .00569661458333333333,
  C88: .3076171875  

};

@@ -3879,7 +3933,18 @@ Proj4js.Proj.sinu = {
	init: function() {
		/* Place parameters in static storage for common use
		  -------------------------------------------------*/
		this.R = 6370997.0; //Radius of earth
		  

		if (!this.sphere) {
		  this.en = Proj4js.common.pj_enfn(this.es);
    } else {
      this.n = 1.;
      this.m = 0.;
      this.es = 0;
      this.C_y = Math.sqrt((this.m + 1.) / this.n);
      this.C_x = this.C_y/(this.m + 1.);
    }
		  
	},

	/* Sinusoidal forward equations--mapping lat,long to x,y
@@ -3890,9 +3955,28 @@ Proj4js.Proj.sinu = {
		var lat=p.y;	
		/* Forward equations
		-----------------*/
		delta_lon = Proj4js.common.adjust_lon(lon - this.long0);
		x = this.R * delta_lon * Math.cos(lat) + this.x0;
		y = this.R * lat + this.y0;
		lon = Proj4js.common.adjust_lon(lon - this.long0);
		
		if (this.sphere) {
      if (!this.m) {
        lat = this.n != 1. ? Math.asin(this.n * Math.sin(lat)): lat;
      } else {
        k = this.n * Math.sin(lat);
        for (var i = Proj4js.common.MAX_ITER; i ; --i) {
          lat -= V = (this.m * lat + Math.sin(lat) - k) / (this.m + Math.cos(lat));
          if (Math.abs(V) < Proj4js.common.EPSLN) break;
        }
      }
      x = this.a * this.C_x * lon * (this.m + Math.cos(lat));
      y = this.a * this.C_y * lat;

		} else {
		  
		  var s = Math.sin(lat);
		  var c = Math.cos(lat);
      y = this.a * Proj4js.common.pj_mlfn(lat, s, c, this.en);
      x = this.a * lon * c / Math.sqrt(1. - this.es * s * s);
		}

		p.x=x;
		p.y=y;	
@@ -3906,16 +3990,27 @@ Proj4js.Proj.sinu = {
		  -----------------*/
		p.x -= this.x0;
		p.y -= this.y0;
		lat = p.y / this.R;
		if (Math.abs(lat) > Proj4js.common.HALF_PI) {
		    Proj4js.reportError("sinu:Inv:DataError");
		}
		temp = Math.abs(lat) - Proj4js.common.HALF_PI;
		if (Math.abs(temp) > Proj4js.common.EPSLN) {
			temp = this.long0+ p.x / (this.R *Math.cos(lat));
			lon = Proj4js.common.adjust_lon(temp);
		lat = p.y / this.a;
		
		if (this.sphere) {
		  
      p.y /= this.C_y;
      lat = this.m ? Math.asin((this.m * p.y + Math.sin(p.y)) / this.n) :
        ( this.n != 1. ? Math.asin(sin(p.y) / this.n) : p.y );
      lon = p.x / (this.C_x * (this.m + Math.cos(p.y)));
		  
		} else {
		  lat = Proj4js.common.pj_inv_mlfn(p.y/this.a, this.es, this.en)
		  var s = Math.abs(lat);
      if (s < Proj4js.common.HALF_PI) {
        s = Math.sin(lat);
        temp = this.long0 + p.x * Math.sqrt(1. - this.es * s * s) /(this.a * Math.cos(lat));
        //temp = this.long0 + p.x / (this.a * Math.cos(lat));
        lon = Proj4js.common.adjust_lon(temp);
      } else if ((s - Proj4js.common.EPSLN) < Proj4js.common.HALF_PI) {
        lon = this.long0;
      }
		  
		}
		  
		p.x=lon;
@@ -4234,8 +4329,8 @@ ALGORITHM REFERENCES
Proj4js.Proj.cass = {
  init : function() {
    if (!this.sphere) {
      this.en = this.pj_enfn(this.es)
      this.m0 = this.pj_mlfn(this.lat0, Math.sin(this.lat0), Math.cos(this.lat0), this.en);
      this.en = Proj4js.common.pj_enfn(this.es)
      this.m0 = Proj4js.common.pj_mlfn(this.lat0, Math.sin(this.lat0), Math.cos(this.lat0), this.en);
    }
  },

@@ -4264,7 +4359,7 @@ Proj4js.Proj.cass = {
        //ellipsoid
      this.n = Math.sin(phi);
      this.c = Math.cos(phi);
      y = this.pj_mlfn(phi, this.n, this.c, this.en);
      y = Proj4js.common.pj_mlfn(phi, this.n, this.c, this.en);
      this.n = 1./Math.sqrt(1. - this.es * this.n * this.n);
      this.tn = Math.tan(phi); 
      this.t = this.tn * this.tn;
@@ -4295,7 +4390,7 @@ Proj4js.Proj.cass = {
      lam = Math.atan2(Math.tan(x), Math.cos(this.dd));
    } else {
      /* ellipsoid */
      var ph1 = this.pj_inv_mlfn(this.m0 + y, this.es, this.en);
      var ph1 = Proj4js.common.pj_inv_mlfn(this.m0 + y, this.es, this.en);
      this.tn = Math.tan(ph1); 
      this.t = this.tn * this.tn;
      this.n = Math.sin(ph1);
@@ -4310,62 +4405,7 @@ Proj4js.Proj.cass = {
    p.x = Proj4js.common.adjust_lon(this.long0+lam);
    p.y = phi;
    return p;
  },//lamazInv()


  //code from the PROJ.4 pj_mlfn.c file;  this may be useful for other projections
  pj_enfn: function(es) {
    var en = new Array();
    en[0] = this.C00 - es * (this.C02 + es * (this.C04 + es * (this.C06 + es * this.C08)));
    en[1] = es * (this.C22 - es * (this.C04 + es * (this.C06 + es * this.C08)));
    var t = es * es;
    en[2] = t * (this.C44 - es * (this.C46 + es * this.C48));
    t *= es;
    en[3] = t * (this.C66 - es * this.C68);
    en[4] = t * es * this.C88;
    return en;
  },
  
  pj_mlfn: function(phi, sphi, cphi, en) {
    cphi *= sphi;
    sphi *= sphi;
    return(en[0] * phi - cphi * (en[1] + sphi*(en[2]+ sphi*(en[3] + sphi*en[4]))));
  },
  
  pj_inv_mlfn: function(arg, es, en) {
    var k = 1./(1.-es);
    var phi = arg;
    for (var i = Proj4js.common.MAX_ITER; i ; --i) { /* rarely goes over 2 iterations */
      var s = Math.sin(phi);
      t = 1. - es * s * s;
      //t = this.pj_mlfn(phi, s, Math.cos(phi), en) - arg;
      //phi -= t * (t * Math.sqrt(t)) * k;
      t = (this.pj_mlfn(phi, s, Math.cos(phi), en) - arg) * (t * Math.sqrt(t)) * k;
      phi -= t;
      if (Math.abs(t) < Proj4js.common.EPSLN)
        return phi;
    }
    Proj4js.reportError("cass:pj_inv_mlfn: Convergence error");
    return phi;
  },

/* meridinal distance for ellipsoid and inverse
**	8th degree - accurate to < 1e-5 meters when used in conjuction
**		with typical major axis values.
**	Inverse determines phi to EPS (1e-11) radians, about 1e-6 seconds.
*/
  C00: 1.0,
  C02: .25,
  C04: .046875,
  C06: .01953125,
  C08: .01068115234375,
  C22: .75,
  C44: .46875,
  C46: .01302083333333333333,
  C48: .00712076822916666666,
  C66: .36458333333333333333,
  C68: .00569661458333333333,
  C88: .3076171875
  }//cassInv()

}
/* ======================================================================
@@ -5062,9 +5102,9 @@ Proj4js.Proj.laea = {
          y = cosz * rh;
          break;
        case this.OBLIQ:
          phi = (Math.abs(rh) <= Proj4js.common.EPSLN) ? this.phi0 : Math.asin(cosz * sinph0 + y * sinz * cosph0 / rh);
          x *= sinz * cosph0;
          y = (cosz - Math.sin(phi) * sinph0) * rh;
          phi = (Math.abs(rh) <= Proj4js.common.EPSLN) ? this.phi0 : Math.asin(cosz * this.sinph0 + y * sinz * this.cosph0 / rh);
          x *= sinz * this.cosph0;
          y = (cosz - Math.sin(phi) * this.sinph0) * rh;
          break;
        case this.N_POLE:
          y = -y;
+12 −10

File changed.

Preview size limit exceeded, changes collapsed.

+55 −1
Original line number Diff line number Diff line
@@ -1200,7 +1200,61 @@ Proj4js.common = {
  {
    var temp= e*sinphi;
    return a/Math.sqrt(1.0 - temp*temp);
  },
  
  //code from the PROJ.4 pj_mlfn.c file;  this may be useful for other projections
  pj_enfn: function(es) {
    var en = new Array();
    en[0] = this.C00 - es * (this.C02 + es * (this.C04 + es * (this.C06 + es * this.C08)));
    en[1] = es * (this.C22 - es * (this.C04 + es * (this.C06 + es * this.C08)));
    var t = es * es;
    en[2] = t * (this.C44 - es * (this.C46 + es * this.C48));
    t *= es;
    en[3] = t * (this.C66 - es * this.C68);
    en[4] = t * es * this.C88;
    return en;
  },
  
  pj_mlfn: function(phi, sphi, cphi, en) {
    cphi *= sphi;
    sphi *= sphi;
    return(en[0] * phi - cphi * (en[1] + sphi*(en[2]+ sphi*(en[3] + sphi*en[4]))));
  },
  
  pj_inv_mlfn: function(arg, es, en) {
    var k = 1./(1.-es);
    var phi = arg;
    for (var i = Proj4js.common.MAX_ITER; i ; --i) { /* rarely goes over 2 iterations */
      var s = Math.sin(phi);
      t = 1. - es * s * s;
      //t = this.pj_mlfn(phi, s, Math.cos(phi), en) - arg;
      //phi -= t * (t * Math.sqrt(t)) * k;
      t = (this.pj_mlfn(phi, s, Math.cos(phi), en) - arg) * (t * Math.sqrt(t)) * k;
      phi -= t;
      if (Math.abs(t) < Proj4js.common.EPSLN)
        return phi;
    }
    Proj4js.reportError("cass:pj_inv_mlfn: Convergence error");
    return phi;
  },

/* meridinal distance for ellipsoid and inverse
**	8th degree - accurate to < 1e-5 meters when used in conjuction
**		with typical major axis values.
**	Inverse determines phi to EPS (1e-11) radians, about 1e-6 seconds.
*/
  C00: 1.0,
  C02: .25,
  C04: .046875,
  C06: .01953125,
  C08: .01068115234375,
  C22: .75,
  C44: .46875,
  C46: .01302083333333333333,
  C48: .00712076822916666666,
  C66: .36458333333333333333,
  C68: .00569661458333333333,
  C88: .3076171875  

};

+5 −60
Original line number Diff line number Diff line
@@ -27,8 +27,8 @@ ALGORITHM REFERENCES
Proj4js.Proj.cass = {
  init : function() {
    if (!this.sphere) {
      this.en = this.pj_enfn(this.es)
      this.m0 = this.pj_mlfn(this.lat0, Math.sin(this.lat0), Math.cos(this.lat0), this.en);
      this.en = Proj4js.common.pj_enfn(this.es)
      this.m0 = Proj4js.common.pj_mlfn(this.lat0, Math.sin(this.lat0), Math.cos(this.lat0), this.en);
    }
  },

@@ -57,7 +57,7 @@ Proj4js.Proj.cass = {
        //ellipsoid
      this.n = Math.sin(phi);
      this.c = Math.cos(phi);
      y = this.pj_mlfn(phi, this.n, this.c, this.en);
      y = Proj4js.common.pj_mlfn(phi, this.n, this.c, this.en);
      this.n = 1./Math.sqrt(1. - this.es * this.n * this.n);
      this.tn = Math.tan(phi); 
      this.t = this.tn * this.tn;
@@ -88,7 +88,7 @@ Proj4js.Proj.cass = {
      lam = Math.atan2(Math.tan(x), Math.cos(this.dd));
    } else {
      /* ellipsoid */
      var ph1 = this.pj_inv_mlfn(this.m0 + y, this.es, this.en);
      var ph1 = Proj4js.common.pj_inv_mlfn(this.m0 + y, this.es, this.en);
      this.tn = Math.tan(ph1); 
      this.t = this.tn * this.tn;
      this.n = Math.sin(ph1);
@@ -103,61 +103,6 @@ Proj4js.Proj.cass = {
    p.x = Proj4js.common.adjust_lon(this.long0+lam);
    p.y = phi;
    return p;
  },//lamazInv()


  //code from the PROJ.4 pj_mlfn.c file;  this may be useful for other projections
  pj_enfn: function(es) {
    var en = new Array();
    en[0] = this.C00 - es * (this.C02 + es * (this.C04 + es * (this.C06 + es * this.C08)));
    en[1] = es * (this.C22 - es * (this.C04 + es * (this.C06 + es * this.C08)));
    var t = es * es;
    en[2] = t * (this.C44 - es * (this.C46 + es * this.C48));
    t *= es;
    en[3] = t * (this.C66 - es * this.C68);
    en[4] = t * es * this.C88;
    return en;
  },
  
  pj_mlfn: function(phi, sphi, cphi, en) {
    cphi *= sphi;
    sphi *= sphi;
    return(en[0] * phi - cphi * (en[1] + sphi*(en[2]+ sphi*(en[3] + sphi*en[4]))));
  },
  
  pj_inv_mlfn: function(arg, es, en) {
    var k = 1./(1.-es);
    var phi = arg;
    for (var i = Proj4js.common.MAX_ITER; i ; --i) { /* rarely goes over 2 iterations */
      var s = Math.sin(phi);
      t = 1. - es * s * s;
      //t = this.pj_mlfn(phi, s, Math.cos(phi), en) - arg;
      //phi -= t * (t * Math.sqrt(t)) * k;
      t = (this.pj_mlfn(phi, s, Math.cos(phi), en) - arg) * (t * Math.sqrt(t)) * k;
      phi -= t;
      if (Math.abs(t) < Proj4js.common.EPSLN)
        return phi;
    }
    Proj4js.reportError("cass:pj_inv_mlfn: Convergence error");
    return phi;
  },

/* meridinal distance for ellipsoid and inverse
**	8th degree - accurate to < 1e-5 meters when used in conjuction
**		with typical major axis values.
**	Inverse determines phi to EPS (1e-11) radians, about 1e-6 seconds.
*/
  C00: 1.0,
  C02: .25,
  C04: .046875,
  C06: .01953125,
  C08: .01068115234375,
  C22: .75,
  C44: .46875,
  C46: .01302083333333333333,
  C48: .00712076822916666666,
  C66: .36458333333333333333,
  C68: .00569661458333333333,
  C88: .3076171875
  }//cassInv()

}
+55 −14
Original line number Diff line number Diff line
@@ -31,7 +31,18 @@ Proj4js.Proj.sinu = {
	init: function() {
		/* Place parameters in static storage for common use
		  -------------------------------------------------*/
		this.R = 6370997.0; //Radius of earth
		  

		if (!this.sphere) {
		  this.en = Proj4js.common.pj_enfn(this.es);
    } else {
      this.n = 1.;
      this.m = 0.;
      this.es = 0;
      this.C_y = Math.sqrt((this.m + 1.) / this.n);
      this.C_x = this.C_y/(this.m + 1.);
    }
		  
	},

	/* Sinusoidal forward equations--mapping lat,long to x,y
@@ -42,9 +53,28 @@ Proj4js.Proj.sinu = {
		var lat=p.y;	
		/* Forward equations
		-----------------*/
		delta_lon = Proj4js.common.adjust_lon(lon - this.long0);
		x = this.R * delta_lon * Math.cos(lat) + this.x0;
		y = this.R * lat + this.y0;
		lon = Proj4js.common.adjust_lon(lon - this.long0);
		
		if (this.sphere) {
      if (!this.m) {
        lat = this.n != 1. ? Math.asin(this.n * Math.sin(lat)): lat;
      } else {
        k = this.n * Math.sin(lat);
        for (var i = Proj4js.common.MAX_ITER; i ; --i) {
          lat -= V = (this.m * lat + Math.sin(lat) - k) / (this.m + Math.cos(lat));
          if (Math.abs(V) < Proj4js.common.EPSLN) break;
        }
      }
      x = this.a * this.C_x * lon * (this.m + Math.cos(lat));
      y = this.a * this.C_y * lat;

		} else {
		  
		  var s = Math.sin(lat);
		  var c = Math.cos(lat);
      y = this.a * Proj4js.common.pj_mlfn(lat, s, c, this.en);
      x = this.a * lon * c / Math.sqrt(1. - this.es * s * s);
		}

		p.x=x;
		p.y=y;	
@@ -58,16 +88,27 @@ Proj4js.Proj.sinu = {
		  -----------------*/
		p.x -= this.x0;
		p.y -= this.y0;
		lat = p.y / this.R;
		if (Math.abs(lat) > Proj4js.common.HALF_PI) {
		    Proj4js.reportError("sinu:Inv:DataError");
		}
		temp = Math.abs(lat) - Proj4js.common.HALF_PI;
		if (Math.abs(temp) > Proj4js.common.EPSLN) {
			temp = this.long0+ p.x / (this.R *Math.cos(lat));
			lon = Proj4js.common.adjust_lon(temp);
		lat = p.y / this.a;
		
		if (this.sphere) {
		  
      p.y /= this.C_y;
      lat = this.m ? Math.asin((this.m * p.y + Math.sin(p.y)) / this.n) :
        ( this.n != 1. ? Math.asin(sin(p.y) / this.n) : p.y );
      lon = p.x / (this.C_x * (this.m + Math.cos(p.y)));
		  
		} else {
		  lat = Proj4js.common.pj_inv_mlfn(p.y/this.a, this.es, this.en)
		  var s = Math.abs(lat);
      if (s < Proj4js.common.HALF_PI) {
        s = Math.sin(lat);
        temp = this.long0 + p.x * Math.sqrt(1. - this.es * s * s) /(this.a * Math.cos(lat));
        //temp = this.long0 + p.x / (this.a * Math.cos(lat));
        lon = Proj4js.common.adjust_lon(temp);
      } else if ((s - Proj4js.common.EPSLN) < Proj4js.common.HALF_PI) {
        lon = this.long0;
      }
		  
		}
		  
		p.x=lon;
Loading