sys/src/cmd/map/libmap/lune.c

   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 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 }