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