xref: /plan9-contrib/sys/src/cmd/map/libmap/lambert.c (revision 219b2ee8daee37f4aad58d63f21287faa8e4ffdc)
1 #include "map.h"
2 
3 static struct coord stdp0, stdp1;
4 static double k;
5 
6 static int
7 Xlambert(struct place *place, double *x, double *y)
8 {
9 	double r;
10 	if(place->nlat.l < -80.*RAD)
11 		return(-1);
12 	if(place->nlat.l > 89.*RAD)
13 		r = 0;	/* slovenly */
14 	else
15 		r = stdp0.c*exp(0.5*k*log(
16 		   (1+stdp0.s)*(1-place->nlat.s)/((1-stdp0.s)*(1+place->nlat.s))));
17 	if(stdp1.l<0.)
18 		r = -r;
19 	*x = - r*sin(k * place->wlon.l);
20 	*y = - r*cos(k * place->wlon.l);
21 	return(1);
22 }
23 
24 proj
25 lambert(double par0, double par1)
26 {
27 	double temp;
28 	if(fabs(par0)>fabs(par1)){
29 		temp = par0;
30 		par0 = par1;
31 		par1 = temp;
32 	}
33 	deg2rad(par0, &stdp0);
34 	deg2rad(par1, &stdp1);
35 	if(fabs(par1+par0)<.1)
36 		return(mercator());
37 	if(fabs(par1-par0)<.1)
38 		return(perspective(-1.));
39 	if(fabs(par0)>89.5||fabs(par1)>89.5)
40 		return(0);
41 	k = 2*log(stdp1.c/stdp0.c)/log(
42 		(1+stdp0.s)*(1-stdp1.s)/((1-stdp0.s)*(1+stdp1.s)));
43 	return(Xlambert);
44 }
45