1*59cc4ca5SDavid du Colombier #include <u.h>
2*59cc4ca5SDavid du Colombier #include <libc.h>
37dd7cddfSDavid du Colombier #include "map.h"
47dd7cddfSDavid du Colombier
57dd7cddfSDavid du Colombier int Xstereographic(struct place *place, double *x, double *y);
67dd7cddfSDavid du Colombier
77dd7cddfSDavid du Colombier static struct place eastpole;
87dd7cddfSDavid du Colombier static struct place westpole;
97dd7cddfSDavid du Colombier static double eastx, easty;
107dd7cddfSDavid du Colombier static double westx, westy;
117dd7cddfSDavid du Colombier static double scale;
127dd7cddfSDavid du Colombier static double pwr;
137dd7cddfSDavid du Colombier
147dd7cddfSDavid du Colombier /* conformal map w = ((1+z)^A - (1-z)^A)/((1+z)^A + (1-z)^A),
157dd7cddfSDavid du Colombier where A<1, maps unit circle onto a convex lune with x= +-1
167dd7cddfSDavid du Colombier mapping to vertices of angle A*PI at w = +-1 */
177dd7cddfSDavid du Colombier
187dd7cddfSDavid du Colombier /* there are cuts from E and W poles to S pole,
197dd7cddfSDavid du Colombier in absence of a cut routine, error is returned for
207dd7cddfSDavid du Colombier points outside a polar cap through E and W poles */
217dd7cddfSDavid du Colombier
Xlune(struct place * place,double * x,double * y)227dd7cddfSDavid du Colombier static Xlune(struct place *place, double *x, double *y)
237dd7cddfSDavid du Colombier {
247dd7cddfSDavid du Colombier double stereox, stereoy;
257dd7cddfSDavid du Colombier double z1x, z1y, z2x, z2y;
267dd7cddfSDavid du Colombier double w1x, w1y, w2x, w2y;
277dd7cddfSDavid du Colombier double numx, numy, denx, deny;
287dd7cddfSDavid du Colombier if(place->nlat.l < eastpole.nlat.l-FUZZ)
297dd7cddfSDavid du Colombier return -1;
307dd7cddfSDavid du Colombier Xstereographic(place, &stereox, &stereoy);
317dd7cddfSDavid du Colombier stereox *= scale;
327dd7cddfSDavid du Colombier stereoy *= scale;
337dd7cddfSDavid du Colombier z1x = 1 + stereox;
347dd7cddfSDavid du Colombier z1y = stereoy;
357dd7cddfSDavid du Colombier z2x = 1 - stereox;
367dd7cddfSDavid du Colombier z2y = -stereoy;
377dd7cddfSDavid du Colombier cpow(z1x,z1y,&w1x,&w1y,pwr);
387dd7cddfSDavid du Colombier cpow(z2x,z2y,&w2x,&w2y,pwr);
397dd7cddfSDavid du Colombier numx = w1x - w2x;
407dd7cddfSDavid du Colombier numy = w1y - w2y;
417dd7cddfSDavid du Colombier denx = w1x + w2x;
427dd7cddfSDavid du Colombier deny = w1y + w2y;
437dd7cddfSDavid du Colombier cdiv(numx, numy, denx, deny, x, y);
447dd7cddfSDavid du Colombier return 1;
457dd7cddfSDavid du Colombier }
467dd7cddfSDavid du Colombier
477dd7cddfSDavid du Colombier proj
lune(double lat,double theta)487dd7cddfSDavid du Colombier lune(double lat, double theta)
497dd7cddfSDavid du Colombier {
507dd7cddfSDavid du Colombier deg2rad(lat, &eastpole.nlat);
517dd7cddfSDavid du Colombier deg2rad(-90.,&eastpole.wlon);
527dd7cddfSDavid du Colombier deg2rad(lat, &westpole.nlat);
537dd7cddfSDavid du Colombier deg2rad(90. ,&westpole.wlon);
547dd7cddfSDavid du Colombier Xstereographic(&eastpole, &eastx, &easty);
557dd7cddfSDavid du Colombier Xstereographic(&westpole, &westx, &westy);
567dd7cddfSDavid du Colombier if(fabs(easty)>FUZZ || fabs(westy)>FUZZ ||
577dd7cddfSDavid du Colombier fabs(eastx+westx)>FUZZ)
587dd7cddfSDavid du Colombier abort();
597dd7cddfSDavid du Colombier scale = 1/eastx;
607dd7cddfSDavid du Colombier pwr = theta/180;
617dd7cddfSDavid du Colombier return Xlune;
627dd7cddfSDavid du Colombier }
63