Subversion Repositories planix.SVN

Rev

Rev 2 | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
2 - 1 
#include <u.h>
2 
#include <libc.h>
3 
#include "map.h"
4 
 
5 
static struct place gywhem, gyehem;
6 
static struct coord gytwist;
7 
static double gyconst, gykc, gyside;
8 
 
9 
 
10 
static void
11 
dosquare(double z1, double z2, double *x, double *y)
12 
{
13 
	double w1,w2;
14 
	w1 = z1 -1;
15 
	if(fabs(w1*w1+z2*z2)>.000001) {
16 
		cdiv(z1+1,z2,w1,z2,&w1,&w2);
17 
		w1 *= gyconst;
18 
		w2 *= gyconst;
19 
		if(w1<0)
20 
			w1 = 0;
21 
		elco2(w1,w2,gykc,1.,1.,x,y);
22 
	} else {
23 
		*x = gyside;
24 
		*y = 0;
25 
	}
26 
}
27 
 
28 
int
29 
Xguyou(struct place *place, double *x, double *y)
30 
{
31 
	int ew;		/*which hemisphere*/
32 
	double z1,z2;
33 
	struct place pl;
34 
	ew = place->wlon.l<0;
35 
	copyplace(place,&pl);
36 
	norm(&pl,ew?&gyehem:&gywhem,&gytwist);
37 
	Xstereographic(&pl,&z1,&z2);
38 
	dosquare(z1/2,z2/2,x,y);
39 
	if(!ew)
40 
		*x -= gyside;
41 
	return(1);
42 
}
43 
 
44 
proj
45 
guyou(void)
46 
{
47 
	double junk;
48 
	gykc = 1/(3+2*sqrt(2.));
49 
	gyconst = -(1+sqrt(2.));
50 
	elco2(-gyconst,0.,gykc,1.,1.,&gyside,&junk);
51 
	gyside *= 2;
52 
	latlon(0.,90.,&gywhem);
53 
	latlon(0.,-90.,&gyehem);
54 
	deg2rad(0.,&gytwist);
55 
	return(Xguyou);
56 
}
57 
 
58 
int
59 
guycut(struct place *g, struct place *og, double *cutlon)
60 
{
61 
	int c;
62 
	c = picut(g,og,cutlon);
63 
	if(c!=1)
64 
		return(c);
65 
	*cutlon = 0.;
66 
	if(g->nlat.c<.7071||og->nlat.c<.7071)
67 
		return(ckcut(g,og,0.));
68 
	return(1);
69 
}
70 
 
71 
static int
72 
Xsquare(struct place *place, double *x, double *y)
73 
{
74 
	double z1,z2;
75 
	double r, theta;
76 
	struct place p;
77 
	copyplace(place,&p);
78 
	if(place->nlat.l<0) {
79 
		p.nlat.l = -p.nlat.l;
80 
		p.nlat.s = -p.nlat.s;
81 
	}
82 
	if(p.nlat.l<FUZZ && fabs(p.wlon.l)>PI-FUZZ){
83 
		*y = -gyside/2;
84 
		*x = p.wlon.l>0?0:gyside;
85 
		return(1);
86 
	}
87 
	Xstereographic(&p,&z1,&z2);
88 
	r = sqrt(sqrt(hypot(z1,z2)/2));
89 
	theta = atan2(z1,-z2)/4;
90 
	dosquare(r*sin(theta),-r*cos(theta),x,y);
91 
	if(place->nlat.l<0)
92 
		*y = -gyside - *y;
93 
	return(1);
94 
}
95 
 
96 
proj
97 
square(void)
98 
{
99 
	guyou();
100 
	return(Xsquare);
101 
}