// NOW ALBERS var EradiusKM = 6378.2064; var Ec = 0.082271854; // eccentricity var Ecsq = Ec * Ec; var Eradius = 1.0; var dtor = Math.PI/180.0; function calc_q(lat) { var s,es; s = Math.sin(lat); es = s * Ec; // echo "lat=lat,es=es,s=s,Ecsq=Ecsq\n"; return (1.0 - Ecsq) * ((s / (1 - es * es))- (1 / (2 * Ec)) * Math.log((1 - es) / (1 + es))); } function calc_msq(lat) { var s,es,c; s = Math.sin(lat); es = s * Ec; c = Math.cos(lat); return c * c/ (1 - es * es); } var cone_const,albersProjectX,albersProjectY; var albersTargetY,albersTargetX,albersCx,albersCy; function albersSetup(p) /* * Pre-compute a bunch of variables which are used by the * albers_project() * * southlat Southern latitude for Albers projection * northlat Northern latitude for Albers projection * originlon Longitude for origin of projected map * originlat Latitude for origin of projected map - * often (northlat + southlat) / 2 */ { // global Eradius,cone_const,middlelon,bigC,r0; var q1,m1sq,q2,m2sq,q0; middlelon = p.originlon; q1 = calc_q(p.southlat); // echo "q1=q1\n"; m1sq = calc_msq(p.southlat); q2 = calc_q(p.northlat); // echo "q2=q2\n"; m2sq = calc_msq(p.northlat); q0 = calc_q(p.originlat); cone_const = (m1sq - m2sq) / (q2 - q1); bigC = m1sq + cone_const * q1; r0 = Eradius * Math.sqrt(bigC - cone_const * q0) / cone_const; } function albersProject(lond,latd) /* * Project lon/lat (in radians) according to albers_setup and * return the results via xp, yp. Units of *xp and *yp are same * as the units used by CONST_Eradius. */ { var lon,lat,q,theta,r,xp,yp,rs; // echo "dtor=dtor,cone_const=cone_const\n"; lon = dtor * lond; lat = dtor * latd; q = calc_q(lat); theta = cone_const * (lon - middlelon); r = Eradius * Math.sqrt(bigC - cone_const * q) / cone_const; xp = r * Math.sin(theta); yp = r0 - r * Math.cos(theta); rs = new Object(); rs . x = xp; rs . y = yp; return rs; } function albersProjectP(p) { return albersProject(p.x,p.y); } function albersProjectX(lond) { var rs,dl; // alert("lond = "+lond); rs = albersProject(lond,albersCy); // alert("rs + "+rs.x+" "+rs.y+" target "+albersTargetX); dl = (rs.x)-albersTargetX; //alert(" lond = "+lond+" lat = "+albersCy+" DL X ="+ dl); return dl; } function albersProjectY(latd) { var rs; rs = albersProject(albersCx,latd); dl = (rs.y)-albersTargetY; // echo " lond = albersCx, lat = latd DL Y = dl\n"; return dl; } //general solver; newton-raphson; if an end of the interval // is reached and the derivative is negative at this point, the end is // returned, and result kind = 'pastend' // the other result kinds are 'failure' and 'success' function newtonIterate(n) { var f,dx,nx,ny,atmin,atmax; f = n.theFunction; // alert(f); if (n.cur_iter == 0) { n.cur_x = n.initial_x; n.cur_y = f(n.cur_x); n.cur_slope = 0.0; } n.cur_iter = n.cur_iter + 1; // if (n.verbose) // echo 'cur_x = '.n.cur_x.' cur_y = '.n.cur_y.' cur_slope = '.n.cur_slope.' iter '.n.cur_iter."\n"; // alert(n.cur_y); if (Math.abs(n.cur_y) < n.tolerance) { n.resultKind = 'success'; return true; } if (n.cur_iter == 1) { if (n.initial_jump == 0.0) dx = n.max_jump; else dx = n.initial_jump; } else if (Math.abs(n.cur_slope) < 0.00001) dx = n.max_jump; else { dx = -(n.cur_y/n.cur_slope); if (dx < -(n.max_jump)) dx = -(n.max_jump); else if (dx > n.max_jump) dx = n.max_jump; } // echo "dx = dx\n"; nx = n.cur_x + dx; if (nx > n.max_x) { nx = n.max_x; atmax = true; } else if (nx < n.min_x) { nx = n.min_x; atmin = true; } ny = f(nx); // echo "ny = ny\n"; n.cur_slope = (ny - n.cur_y)/dx; // if at end, and cur_slope indicates the solution beyond the extreme, // check the derivative using min_jump before concluding pastend if (atmax && (n.cur_slope < 0)) { // echo "ATMAX XXX\n"; dx = n.min_jump; n.cur_y = f(n.max_x-dx); n.cur_slope = (ny - n.cur_y)/dx; if (n.cur_slope < 0) { n.cur_x = n.max_x; n.resultKind = 'pastEnd'; return true; } } else if (atmin && (n.cur_slope > 0)) { // echo "ATMIN XXX\n"; dx = n.min_jump; n.cur_y = ny; ny = f(n.min_x+dx); n.cur_slope = (ny - n.cur_y)/dx; if (n.cur_slope > 0) { n.cur_x = n.min_x; n.resultKind = 'pastEnd'; return true; } } n.cur_x = nx; n.cur_y = ny; return false; } function newtonSolve(n,iv) { var f,done; f = n.theFunction; n.cur_iter = 0; n.initial_x = iv; done = false; n.resultKind = 'failure'; while (!done && (n.cur_iter < n.max_iter)) done = newtonIterate(n); return n.cur_x; } function mkAlbersNS() { var n; n = new Object(); n.tolerance = 0.0000001; n.max_x = 200.0; n.min_x = -200.0; n.max_jump = 5; n.max_iter = 100; n.verbose = false; n.initial_jump = 1; n.cur_iter = 0; n.cur_slope = 0.0; return n; } function albersInverse(p) { var nx,ny,vx,vy,rs,nit, done; nx = mkAlbersNS(); nx.theFunction = albersProjectX; albersCx = -100; nx.initial_x = albersCx;//the middle of the US, more or less ny = mkAlbersNS(); albersCy = 40; ny.initial_x = albersCy; ny.theFunction = albersProjectY; albersTargetX = p.x; albersTargetY = p.y; done= false; // n.resultKind = 'failure'; while (!done && (nx.cur_iter < 100)) { // echo "\n\n****X****\n"; newtonIterate(nx); albersCx = nx.cur_x; // echo "\n\n****Y****\n"; newtonIterate(ny); albersCy = ny.cur_x; vx = nx.cur_y; vy = ny.cur_y; dst = Math.max(Math.abs(nx.cur_y),Math.abs(ny.cur_y)); // echo "X albersCx Y albersCy VX vx VY vy \n"; // echo " dst = dst\n\n"; done = (dst<0.00001); } nit = nx.cur_iter; // alert("iterations: "+nit);//remove // echo "Albers inverse: nit iterations\n
"; rs = new Object(); rs.x = nx.cur_x; rs.y = ny.cur_x; return rs; } AlbersProject0 = new Object(); AlbersProject0 . southlat = dtor*29.5; AlbersProject0 . northlat = dtor * 45.5; AlbersProject0 . originlon = -96.0 * dtor; AlbersProject0 . originlat = dtor * 23.0; //AlbersProject0 . xform = mk_scale(10); //AlbersProject.setup(); //alert("hoooo"); albersSetup(AlbersProject0);