thoming.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
       ---
       thoming.c (2021B)
       ---
            1 #include <u.h>
            2 #include <libc.h>
            3 #include "map.h"
            4 
            5 static struct coord p0;                /* standard parallel */
            6 
            7 int first;
            8 
            9 static double
           10 trigclamp(double x)
           11 {
           12         return x>1? 1: x<-1? -1: x;
           13 }
           14 
           15 static struct coord az;                /* azimuth of p0 seen from place */
           16 static struct coord rad;        /* angular dist from place to p0 */
           17 
           18 static int
           19 azimuth(struct place *place)
           20 {
           21         if(place->nlat.c < FUZZ) {
           22                 az.l = PI/2 + place->nlat.l - place->wlon.l;
           23                 sincos(&az);
           24                 rad.l = fabs(place->nlat.l - p0.l);
           25                 if(rad.l > PI)
           26                         rad.l = 2*PI - rad.l;
           27                 sincos(&rad);
           28                 return 1;
           29         }
           30         rad.c = trigclamp(p0.s*place->nlat.s +        /* law of cosines */
           31                 p0.c*place->nlat.c*place->wlon.c);
           32         rad.s = sqrt(1 - rad.c*rad.c);
           33         if(fabs(rad.s) < .001) {
           34                 az.s = 0;
           35                 az.c = 1;
           36         } else {
           37                 az.s = trigclamp(p0.c*place->wlon.s/rad.s); /* sines */
           38                 az.c = trigclamp((p0.s - rad.c*place->nlat.s)
           39                                 /(rad.s*place->nlat.c));
           40         }
           41         rad.l = atan2(rad.s, rad.c);
           42         return 1;
           43 }
           44 
           45 static int
           46 Xmecca(struct place *place, double *x, double *y)
           47 {
           48         if(!azimuth(place))
           49                 return 0;
           50         *x = -place->wlon.l;
           51         *y = fabs(az.s)<.02? -az.c*rad.s/p0.c: *x*az.c/az.s;
           52         return fabs(*y)>2? -1:
           53                rad.c<0? 0:
           54                1;
           55 }
           56 
           57 proj
           58 mecca(double par)
           59 {
           60         first = 1;
           61         if(fabs(par)>80.)
           62                 return(0);
           63         deg2rad(par,&p0);
           64         return(Xmecca);
           65 }
           66 
           67 static int
           68 Xhoming(struct place *place, double *x, double *y)
           69 {
           70         if(!azimuth(place))
           71                 return 0;
           72         *x = -rad.l*az.s;
           73         *y = -rad.l*az.c;
           74         return place->wlon.c<0? 0: 1;
           75 }
           76 
           77 proj
           78 homing(double par)
           79 {
           80         first = 1;
           81         if(fabs(par)>80.)
           82                 return(0);
           83         deg2rad(par,&p0);
           84         return(Xhoming);
           85 }
           86 
           87 int
           88 hlimb(double *lat, double *lon, double res)
           89 {
           90         if(first) {
           91                 *lon = -90;
           92                 *lat = -90;
           93                 first = 0;
           94                 return 0;
           95         }
           96         *lat += res;
           97         if(*lat <= 90)
           98                 return 1;
           99         if(*lon == 90)
          100                 return -1;
          101         *lon = 90;
          102         *lat = -90;
          103         return 0;
          104 }
          105 
          106 int
          107 mlimb(double *lat, double *lon, double res)
          108 {
          109         int ret = !first;
          110         if(fabs(p0.s) < .01)
          111                 return -1;
          112         if(first) {
          113                 *lon = -180;
          114                 first = 0;
          115         } else
          116                 *lon += res;
          117         if(*lon > 180)
          118                 return -1;
          119         *lat = atan(-cos(*lon*RAD)/p0.s*p0.c)/RAD;
          120         return ret;
          121 }