1 #include <u.h>
2 #include <libc.h>
3 #include "map.h"
4 
5 int Xstereographic(struct place *place, double *x, double *y);
6 
7 static struct place eastpole;
8 static struct place westpole;
9 static double eastx, easty;
10 static double westx, westy;
11 static double scale;
12 static double pwr;
13 
14 /* conformal map w = ((1+z)^A - (1-z)^A)/((1+z)^A + (1-z)^A),
15    where A<1, maps unit circle onto a convex lune with x= +-1
16    mapping to vertices of angle A*PI at w = +-1 */
17 
18 /* there are cuts from E and W poles to S pole,
19    in absence of a cut routine, error is returned for
20    points outside a polar cap through E and W poles */
21 
22 static int Xlune(struct place *place, double *x, double *y)
23 {
24 	double stereox, stereoy;
25 	double z1x, z1y, z2x, z2y;
26 	double w1x, w1y, w2x, w2y;
27 	double numx, numy, denx, deny;
28 	if(place->nlat.l < eastpole.nlat.l-FUZZ)
29 		return	-1;
30 	Xstereographic(place, &stereox, &stereoy);
31 	stereox *= scale;
32 	stereoy *= scale;
33 	z1x = 1 + stereox;
34 	z1y = stereoy;
35 	z2x = 1 - stereox;
36 	z2y = -stereoy;
37 	cpow(z1x,z1y,&w1x,&w1y,pwr);
38 	cpow(z2x,z2y,&w2x,&w2y,pwr);
39 	numx = w1x - w2x;
40 	numy = w1y - w2y;
41 	denx = w1x + w2x;
42 	deny = w1y + w2y;
43 	cdiv(numx, numy, denx, deny, x, y);
44 	return 1;
45 }	
46 
47 proj
48 lune(double lat, double theta)
49 {
50 	deg2rad(lat, &eastpole.nlat);
51 	deg2rad(-90.,&eastpole.wlon);
52 	deg2rad(lat, &westpole.nlat);
53 	deg2rad(90. ,&westpole.wlon);
54 	Xstereographic(&eastpole, &eastx, &easty);
55 	Xstereographic(&westpole, &westx, &westy);
56 	if(fabs(easty)>FUZZ || fabs(westy)>FUZZ ||
57 	   fabs(eastx+westx)>FUZZ)
58 		abort();
59 	scale = 1/eastx;
60 	pwr = theta/180;
61 	return Xlune;
62 }