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