lune.c (1565B)
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 }