Warning: Attempt to read property "date" on null in /usr/local/www/websvn.planix.org/blame.php on line 247

Warning: Attempt to read property "msg" on null in /usr/local/www/websvn.planix.org/blame.php on line 247
WebSVN – planix.SVN – Blame – /os/branches/feature_posix/sys/src/cmd/map/libmap/albers.c – Rev 2

Subversion Repositories planix.SVN

Rev

Go to most recent revision | Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 1 
#include <u.h>
2 
#include <libc.h>
3 
#include "map.h"
4 
 
5 
/* For Albers formulas see Deetz and Adams "Elements of Map Projection", */
6 
/* USGS Special Publication No. 68, GPO 1921 */
7 
 
8 
static double r0sq, r1sq, d2, n, den, sinb1, sinb2;
9 
static struct coord plat1, plat2;
10 
static southpole;
11 
 
12 
static double num(double s)
13 
{
14 
	if(d2==0)
15 
		return(1);
16 
	s = d2*s*s;
17 
	return(1+s*(2./3+s*(3./5+s*(4./7+s*5./9))));
18 
}
19 
 
20 
/* Albers projection for a spheroid, good only when N pole is fixed */
21 
 
22 
static int
23 
Xspalbers(struct place *place, double *x, double *y)
24 
{
25 
	double r = sqrt(r0sq-2*(1-d2)*place->nlat.s*num(place->nlat.s)/n);
26 
	double t = n*place->wlon.l;
27 
	*y = r*cos(t);
28 
	*x = -r*sin(t);
29 
	if(!southpole)
30 
		*y = -*y;
31 
	else
32 
		*x = -*x;
33 
	return(1);
34 
}
35 
 
36 
/* lat1, lat2: std parallels; e2: squared eccentricity */
37 
 
38 
static proj albinit(double lat1, double lat2, double e2)
39 
{
40 
	double r1;
41 
	double t;
42 
	for(;;) {
43 
		if(lat1 < -90)
44 
			lat1 = -180 - lat1;
45 
		if(lat2 > 90)
46 
			lat2 = 180 - lat2;
47 
		if(lat1 <= lat2)
48 
			break;
49 
		t = lat1; lat1 = lat2; lat2 = t;
50 
	}
51 
	if(lat2-lat1 < 1) {
52 
		if(lat1 > 89)
53 
			return(azequalarea());
54 
		return(0);
55 
	}
56 
	if(fabs(lat2+lat1) < 1)
57 
		return(cylequalarea(lat1));
58 
	d2 = e2;
59 
	den = num(1.);
60 
	deg2rad(lat1,&plat1);
61 
	deg2rad(lat2,&plat2);
62 
	sinb1 = plat1.s*num(plat1.s)/den;
63 
	sinb2 = plat2.s*num(plat2.s)/den;
64 
	n = (plat1.c*plat1.c/(1-e2*plat1.s*plat1.s) -
65 
	    plat2.c*plat2.c/(1-e2*plat2.s*plat2.s)) /
66 
	    (2*(1-e2)*den*(sinb2-sinb1));
67 
	r1 = plat1.c/(n*sqrt(1-e2*plat1.s*plat1.s));
68 
	r1sq = r1*r1;
69 
	r0sq = r1sq + 2*(1-e2)*den*sinb1/n;
70 
	southpole = lat1<0 && plat2.c>plat1.c;
71 
	return(Xspalbers);
72 
}
73 
 
74 
proj
75 
sp_albers(double lat1, double lat2)
76 
{
77 
	return(albinit(lat1,lat2,EC2));
78 
}
79 
 
80 
proj
81 
albers(double lat1, double lat2)
82 
{
83 
	return(albinit(lat1,lat2,0.));
84 
}
85 
 
86 
static double scale = 1;
87 
static double twist = 0;
88 
 
89 
void
90 
albscale(double x, double y, double lat, double lon)
91 
{
92 
	struct place place;
93 
	double alat, alon, x1,y1;
94 
	scale = 1;
95 
	twist = 0;
96 
	invalb(x,y,&alat,&alon);
97 
	twist = lon - alon;
98 
	deg2rad(lat,&place.nlat);
99 
	deg2rad(lon,&place.wlon);
100 
	Xspalbers(&place,&x1,&y1);
101 
	scale = sqrt((x1*x1+y1*y1)/(x*x+y*y));
102 
}
103 
 
104 
void
105 
invalb(double x, double y, double *lat, double *lon)
106 
{
107 
	int i;
108 
	double sinb_den, sinp;
109 
	x *= scale;
110 
	y *= scale;
111 
	*lon = atan2(-x,fabs(y))/(RAD*n) + twist;
112 
	sinb_den = (r0sq - x*x - y*y)*n/(2*(1-d2));
113 
	sinp = sinb_den;
114 
	for(i=0; i<5; i++)
115 
		sinp = sinb_den/num(sinp);
116 
	*lat = asin(sinp)/RAD;
117 
}