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

   1 #include <u.h>
   2 #include <libc.h>
   3 #include "map.h"
   4
   5 #define ORTHRAD 1000
   6 static double viewpt;
   7
   8 static int
   9 Xperspective(struct place *place, double *x, double *y)
  10 {
  11         double r;
  12         if(viewpt<=1+FUZZ && fabs(place->nlat.s<=viewpt+.01))
  13                 return(-1);
  14         r = place->nlat.c*(viewpt - 1.)/(viewpt - place->nlat.s);
  15         *x = - r*place->wlon.s;
  16         *y = - r*place->wlon.c;
  17         if(r>4.)
  18                 return(-1);
  19         if(fabs(viewpt)>1 && place->nlat.s<1/viewpt ||
  20            fabs(viewpt)<=1 && place->nlat.s<viewpt)
  21                         return 0;
  22         return(1);
  23 }
  24
  25 proj
  26 perspective(double radius)
  27 {
  28         viewpt = radius;
  29         if(viewpt >= ORTHRAD)
  30                 return(Xorthographic);
  31         if(fabs(viewpt-1.)<.0001)
  32                 return(0);
  33         return(Xperspective);
  34 }
  35
  36         /* called from various conformal projections,
  37            but not from stereographic itself */
  38 int
  39 Xstereographic(struct place *place, double *x, double *y)
  40 {
  41         double v = viewpt;
  42         int retval;
  43         viewpt = -1;
  44         retval = Xperspective(place, x, y);
  45         viewpt = v;
  46         return retval;
  47 }
  48
  49 proj
  50 stereographic(void)
  51 {
  52         viewpt = -1.;
  53         return(Xperspective);
  54 }
  55
  56 proj
  57 gnomonic(void)
  58 {
  59         viewpt = 0.;
  60         return(Xperspective);
  61 }
  62
  63 int
  64 plimb(double *lat, double *lon, double res)
  65 {
  66         static first = 1;
  67         if(viewpt >= ORTHRAD)
  68                 return olimb(lat, lon, res);
  69         if(first) {
  70                 first = 0;
  71                 *lon = -180;
  72                 if(fabs(viewpt) < .01)
  73                         *lat = 0;
  74                 else if(fabs(viewpt)<=1)
  75                         *lat = asin(viewpt)/RAD;
  76                 else
  77                         *lat = asin(1/viewpt)/RAD;
  78         } else
  79                 *lon += res;
  80         if(*lon <= 180)
  81                 return 1;
  82         first = 1;
  83         return -1;
  84 }