tlune.c - plan9port - [fork] Plan 9 from user space
  HTML git clone git://src.adamsgaard.dk/plan9port
   DIR Log
   DIR Files
   DIR Refs
   DIR README
   DIR LICENSE
       ---
       tlune.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 }