1 28994509 2004-04-21 devnull #include <u.h>
2 28994509 2004-04-21 devnull #include <libc.h>
3 28994509 2004-04-21 devnull #include <stdio.h>
4 28994509 2004-04-21 devnull #include "map.h"
5 28994509 2004-04-21 devnull #include "iplot.h"
7 28994509 2004-04-21 devnull #define NVERT 20 /* max number of vertices in a -v polygon */
8 28994509 2004-04-21 devnull #define HALFWIDTH 8192 /* output scaled to fit in -HALFWIDTH,HALFWIDTH */
9 28994509 2004-04-21 devnull #define LONGLINES (HALFWIDTH*4) /* permissible segment lengths */
10 28994509 2004-04-21 devnull #define SHORTLINES (HALFWIDTH/8)
11 28994509 2004-04-21 devnull #define SCALERATIO 10 /* of abs to rel data (see map(5)) */
12 28994509 2004-04-21 devnull #define RESOL 2. /* coarsest resolution for tracing grid (degrees) */
13 28994509 2004-04-21 devnull #define TWO_THRD 0.66666666666666667
15 28994509 2004-04-21 devnull int normproj(double, double, double *, double *);
16 28994509 2004-04-21 devnull int posproj(double, double, double *, double *);
17 28994509 2004-04-21 devnull int picut(struct place *, struct place *, double *);
18 28994509 2004-04-21 devnull double reduce(double);
19 28994509 2004-04-21 devnull short getshort(FILE *);
20 28994509 2004-04-21 devnull char *mapindex(char *);
21 28994509 2004-04-21 devnull proj projection;
24 c8b6342d 2005-01-13 devnull static char *mapdir = "#9/map"; /* default map directory */
25 28994509 2004-04-21 devnull struct file {
26 28994509 2004-04-21 devnull char *name;
27 28994509 2004-04-21 devnull char *color;
28 28994509 2004-04-21 devnull char *style;
30 28994509 2004-04-21 devnull static struct file dfltfile = {
31 28994509 2004-04-21 devnull "world", BLACK, SOLID /* default map */
33 28994509 2004-04-21 devnull static struct file *file = &dfltfile; /* list of map files */
34 28994509 2004-04-21 devnull static int nfile = 1; /* length of list */
35 28994509 2004-04-21 devnull static char *currcolor = BLACK; /* current color */
36 28994509 2004-04-21 devnull static char *gridcolor = BLACK;
37 28994509 2004-04-21 devnull static char *bordcolor = BLACK;
39 28994509 2004-04-21 devnull extern struct index index[];
40 28994509 2004-04-21 devnull int halfwidth = HALFWIDTH;
42 28994509 2004-04-21 devnull static int (*cut)(struct place *, struct place *, double *);
43 28994509 2004-04-21 devnull static int (*limb)(double*, double*, double);
44 28994509 2004-04-21 devnull static void dolimb(void);
45 28994509 2004-04-21 devnull static int onlimb;
46 28994509 2004-04-21 devnull static int poles;
47 28994509 2004-04-21 devnull static double orientation[3] = { 90., 0., 0. }; /* -o option */
48 28994509 2004-04-21 devnull static int oriented; /* nonzero if -o option occurred */
49 28994509 2004-04-21 devnull static int upright; /* 1 if orientation[0]==90, -1 if -90, else 0*/
50 28994509 2004-04-21 devnull static int delta = 1; /* -d setting */
51 28994509 2004-04-21 devnull static double limits[4] = { /* -l parameters */
52 28994509 2004-04-21 devnull -90., 90., -180., 180.
54 28994509 2004-04-21 devnull static double klimits[4] = { /* -k parameters */
55 28994509 2004-04-21 devnull -90., 90., -180., 180.
57 28994509 2004-04-21 devnull static int limcase;
58 28994509 2004-04-21 devnull static double rlimits[4]; /* limits expressed in radians */
59 28994509 2004-04-21 devnull static double lolat, hilat, lolon, hilon;
60 28994509 2004-04-21 devnull static double window[4] = { /* option -w */
61 28994509 2004-04-21 devnull -90., 90., -180., 180.
63 28994509 2004-04-21 devnull static int windowed; /* nozero if option -w */
64 28994509 2004-04-21 devnull static struct vert { double x, y; } v[NVERT+2]; /*clipping polygon*/
65 28994509 2004-04-21 devnull static struct edge { double a, b, c; } e[NVERT]; /* coeffs for linear inequality */
66 28994509 2004-04-21 devnull static int nvert; /* number of vertices in clipping polygon */
68 28994509 2004-04-21 devnull static double rwindow[4]; /* window, expressed in radians */
69 28994509 2004-04-21 devnull static double params[2]; /* projection params */
70 28994509 2004-04-21 devnull /* bounds on output values before scaling; found by coarse survey */
71 28994509 2004-04-21 devnull static double xmin = 100.;
72 28994509 2004-04-21 devnull static double xmax = -100.;
73 28994509 2004-04-21 devnull static double ymin = 100.;
74 28994509 2004-04-21 devnull static double ymax = -100.;
75 28994509 2004-04-21 devnull static double xcent, ycent;
76 28994509 2004-04-21 devnull static double xoff, yoff;
77 28994509 2004-04-21 devnull double xrange, yrange;
78 28994509 2004-04-21 devnull static int left = -HALFWIDTH;
79 28994509 2004-04-21 devnull static int right = HALFWIDTH;
80 28994509 2004-04-21 devnull static int bottom = -HALFWIDTH;
81 28994509 2004-04-21 devnull static int top = HALFWIDTH;
82 28994509 2004-04-21 devnull static int longlines = SHORTLINES; /* drop longer segments */
83 28994509 2004-04-21 devnull static int shortlines = SHORTLINES;
84 28994509 2004-04-21 devnull static int bflag = 1; /* 0 for option -b */
85 28994509 2004-04-21 devnull static int s1flag = 0; /* 1 for option -s1 */
86 28994509 2004-04-21 devnull static int s2flag = 0; /* 1 for option -s2 */
87 28994509 2004-04-21 devnull static int rflag = 0; /* 1 for option -r */
88 28994509 2004-04-21 devnull static int kflag = 0; /* 1 if option -k occurred */
89 28994509 2004-04-21 devnull static int xflag = 0; /* 1 for option -x */
90 28994509 2004-04-21 devnull int vflag = 1; /* -1 if option -v occurred */
91 28994509 2004-04-21 devnull static double position[3]; /* option -p */
92 28994509 2004-04-21 devnull static double center[3] = {0., 0., 0.}; /* option -c */
93 28994509 2004-04-21 devnull static struct coord crot; /* option -c */
94 28994509 2004-04-21 devnull static double grid[3] = { 10., 10., RESOL }; /* option -g */
95 28994509 2004-04-21 devnull static double dlat, dlon; /* resolution for tracing grid in lat and lon */
96 28994509 2004-04-21 devnull static double scaling; /* to compute final integer output */
97 28994509 2004-04-21 devnull static struct file *track; /* options -t and -u */
98 28994509 2004-04-21 devnull static int ntrack; /* number of tracks present */
99 28994509 2004-04-21 devnull static char *symbolfile; /* option -y */
101 28994509 2004-04-21 devnull void clamp(double *px, double v);
102 28994509 2004-04-21 devnull void clipinit(void);
103 28994509 2004-04-21 devnull double diddle(struct place *, double, double);
104 28994509 2004-04-21 devnull double diddle(struct place *, double, double);
105 28994509 2004-04-21 devnull void dobounds(double, double, double, double, int);
106 28994509 2004-04-21 devnull void dogrid(double, double, double, double);
107 28994509 2004-04-21 devnull int duple(struct place *, double);
108 b8f742db 2005-01-11 devnull #define fmax map_fmax
109 b8f742db 2005-01-11 devnull #define fmin map_fmin
110 28994509 2004-04-21 devnull double fmax(double, double);
111 28994509 2004-04-21 devnull double fmin(double, double);
112 28994509 2004-04-21 devnull void getdata(char *);
113 28994509 2004-04-21 devnull int gridpt(double, double, int);
114 28994509 2004-04-21 devnull int inpoly(double, double);
115 28994509 2004-04-21 devnull int inwindow(struct place *);
116 28994509 2004-04-21 devnull void pathnames(void);
117 28994509 2004-04-21 devnull int pnorm(double);
118 28994509 2004-04-21 devnull void radbds(double *w, double *rw);
119 28994509 2004-04-21 devnull void revlon(struct place *, double);
120 28994509 2004-04-21 devnull void satellite(struct file *);
121 28994509 2004-04-21 devnull int seeable(double, double);
122 28994509 2004-04-21 devnull void windlim(void);
123 28994509 2004-04-21 devnull void realcut(void);
126 fa325e9b 2020-01-10 cross option(char *s)
129 28994509 2004-04-21 devnull if(s[0]=='-' && (s[1]<'0'||s[1]>'9'))
130 28994509 2004-04-21 devnull return(s[1]!='.'&&s[1]!=0);
132 28994509 2004-04-21 devnull return(0);
136 28994509 2004-04-21 devnull conv(int k, struct coord *g)
138 28994509 2004-04-21 devnull g->l = (0.0001/SCALERATIO)*k;
139 28994509 2004-04-21 devnull sincos(g);
143 28994509 2004-04-21 devnull main(int argc, char *argv[])
145 28994509 2004-04-21 devnull int i,k;
146 28994509 2004-04-21 devnull char *s, *t, *style;
147 28994509 2004-04-21 devnull double x, y;
148 28994509 2004-04-21 devnull double lat, lon;
149 28994509 2004-04-21 devnull double *wlim;
150 28994509 2004-04-21 devnull double dd;
151 28994509 2004-04-21 devnull if(sizeof(short)!=2)
152 28994509 2004-04-21 devnull abort(); /* getshort() won't work */
153 c8b6342d 2005-01-13 devnull mapdir = unsharp(mapdir);
154 28994509 2004-04-21 devnull s = getenv("MAP");
156 28994509 2004-04-21 devnull file[0].name = s;
157 28994509 2004-04-21 devnull s = getenv("MAPDIR");
159 28994509 2004-04-21 devnull mapdir = s;
160 fa325e9b 2020-01-10 cross if(argc<=1)
161 28994509 2004-04-21 devnull error("usage: map projection params options");
162 28994509 2004-04-21 devnull for(k=0;index[k].name;k++) {
163 28994509 2004-04-21 devnull s = index[k].name;
164 28994509 2004-04-21 devnull t = argv[1];
165 28994509 2004-04-21 devnull while(*s == *t){
166 28994509 2004-04-21 devnull if(*s==0) goto found;
171 28994509 2004-04-21 devnull fprintf(stderr,"projections:\n");
172 28994509 2004-04-21 devnull for(i=0;index[i].name;i++) {
173 28994509 2004-04-21 devnull fprintf(stderr,"%s",index[i].name);
174 28994509 2004-04-21 devnull for(k=0; k<index[i].npar; k++)
175 28994509 2004-04-21 devnull fprintf(stderr," p%d", k);
176 28994509 2004-04-21 devnull fprintf(stderr,"\n");
178 28994509 2004-04-21 devnull exits("error");
180 28994509 2004-04-21 devnull argv += 2;
181 28994509 2004-04-21 devnull argc -= 2;
182 28994509 2004-04-21 devnull cut = index[k].cut;
183 28994509 2004-04-21 devnull limb = index[k].limb;
184 28994509 2004-04-21 devnull poles = index[k].poles;
185 28994509 2004-04-21 devnull for(i=0;i<index[k].npar;i++) {
186 28994509 2004-04-21 devnull if(i>=argc||option(argv[i])) {
187 28994509 2004-04-21 devnull fprintf(stderr,"%s needs %d params\n",index[k].name,index[k].npar);
188 28994509 2004-04-21 devnull exits("error");
190 28994509 2004-04-21 devnull params[i] = atof(argv[i]);
192 28994509 2004-04-21 devnull argv += i;
193 28994509 2004-04-21 devnull argc -= i;
194 28994509 2004-04-21 devnull while(argc>0&&option(argv[0])) {
197 28994509 2004-04-21 devnull switch(argv[-1][1]) {
198 28994509 2004-04-21 devnull case 'm':
199 28994509 2004-04-21 devnull if(file == &dfltfile) {
200 28994509 2004-04-21 devnull file = 0;
201 28994509 2004-04-21 devnull nfile = 0;
203 28994509 2004-04-21 devnull while(argc && !option(*argv)) {
204 28994509 2004-04-21 devnull file = realloc(file,(nfile+1)*sizeof(*file));
205 28994509 2004-04-21 devnull file[nfile].name = *argv;
206 28994509 2004-04-21 devnull file[nfile].color = currcolor;
207 28994509 2004-04-21 devnull file[nfile].style = SOLID;
208 28994509 2004-04-21 devnull nfile++;
213 28994509 2004-04-21 devnull case 'b':
214 28994509 2004-04-21 devnull bflag = 0;
215 28994509 2004-04-21 devnull for(nvert=0;nvert<NVERT&&argc>=2;nvert++) {
216 28994509 2004-04-21 devnull if(option(*argv))
218 28994509 2004-04-21 devnull v[nvert].x = atof(*argv++);
220 28994509 2004-04-21 devnull if(option(*argv))
222 28994509 2004-04-21 devnull v[nvert].y = atof(*argv++);
225 28994509 2004-04-21 devnull if(nvert>=NVERT)
226 28994509 2004-04-21 devnull error("too many clipping vertices");
228 28994509 2004-04-21 devnull case 'g':
229 28994509 2004-04-21 devnull gridcolor = currcolor;
230 28994509 2004-04-21 devnull for(i=0;i<3&&argc>i&&!option(argv[i]);i++)
231 28994509 2004-04-21 devnull grid[i] = atof(argv[i]);
232 28994509 2004-04-21 devnull switch(i) {
234 28994509 2004-04-21 devnull grid[0] = grid[1] = 0.;
237 28994509 2004-04-21 devnull grid[1] = grid[0];
239 28994509 2004-04-21 devnull argc -= i;
240 28994509 2004-04-21 devnull argv += i;
242 28994509 2004-04-21 devnull case 't':
243 28994509 2004-04-21 devnull style = SOLID;
244 28994509 2004-04-21 devnull goto casetu;
245 28994509 2004-04-21 devnull case 'u':
246 28994509 2004-04-21 devnull style = DOTDASH;
248 28994509 2004-04-21 devnull while(argc && !option(*argv)) {
249 28994509 2004-04-21 devnull track = realloc(track,(ntrack+1)*sizeof(*track));
250 28994509 2004-04-21 devnull track[ntrack].name = *argv;
251 28994509 2004-04-21 devnull track[ntrack].color = currcolor;
252 28994509 2004-04-21 devnull track[ntrack].style = style;
253 28994509 2004-04-21 devnull ntrack++;
258 28994509 2004-04-21 devnull case 'r':
259 28994509 2004-04-21 devnull rflag++;
261 28994509 2004-04-21 devnull case 's':
262 28994509 2004-04-21 devnull switch(argv[-1][2]) {
263 28994509 2004-04-21 devnull case '1':
264 28994509 2004-04-21 devnull s1flag++;
266 28994509 2004-04-21 devnull case 0: /* compatibility */
267 28994509 2004-04-21 devnull case '2':
268 28994509 2004-04-21 devnull s2flag++;
271 28994509 2004-04-21 devnull case 'o':
272 28994509 2004-04-21 devnull for(i=0;i<3&&i<argc&&!option(argv[i]);i++)
273 28994509 2004-04-21 devnull orientation[i] = atof(argv[i]);
274 28994509 2004-04-21 devnull oriented++;
275 28994509 2004-04-21 devnull argv += i;
276 28994509 2004-04-21 devnull argc -= i;
278 28994509 2004-04-21 devnull case 'l':
279 28994509 2004-04-21 devnull bordcolor = currcolor;
280 28994509 2004-04-21 devnull for(i=0;i<argc&&i<4&&!option(argv[i]);i++)
281 28994509 2004-04-21 devnull limits[i] = atof(argv[i]);
282 28994509 2004-04-21 devnull argv += i;
283 28994509 2004-04-21 devnull argc -= i;
285 28994509 2004-04-21 devnull case 'k':
286 28994509 2004-04-21 devnull kflag++;
287 28994509 2004-04-21 devnull for(i=0;i<argc&&i<4&&!option(argv[i]);i++)
288 28994509 2004-04-21 devnull klimits[i] = atof(argv[i]);
289 28994509 2004-04-21 devnull argv += i;
290 28994509 2004-04-21 devnull argc -= i;
292 28994509 2004-04-21 devnull case 'd':
293 28994509 2004-04-21 devnull if(argc>0&&!option(argv[0])) {
294 28994509 2004-04-21 devnull delta = atoi(argv[0]);
299 28994509 2004-04-21 devnull case 'w':
300 28994509 2004-04-21 devnull bordcolor = currcolor;
301 28994509 2004-04-21 devnull windowed++;
302 28994509 2004-04-21 devnull for(i=0;i<argc&&i<4&&!option(argv[i]);i++)
303 28994509 2004-04-21 devnull window[i] = atof(argv[i]);
304 28994509 2004-04-21 devnull argv += i;
305 28994509 2004-04-21 devnull argc -= i;
307 28994509 2004-04-21 devnull case 'c':
308 fa325e9b 2020-01-10 cross for(i=0;i<3&&argc>i&&!option(argv[i]);i++)
309 28994509 2004-04-21 devnull center[i] = atof(argv[i]);
310 28994509 2004-04-21 devnull argc -= i;
311 28994509 2004-04-21 devnull argv += i;
313 28994509 2004-04-21 devnull case 'p':
314 28994509 2004-04-21 devnull for(i=0;i<3&&argc>i&&!option(argv[i]);i++)
315 28994509 2004-04-21 devnull position[i] = atof(argv[i]);
316 28994509 2004-04-21 devnull argc -= i;
317 28994509 2004-04-21 devnull argv += i;
318 fa325e9b 2020-01-10 cross if(i!=3||position[2]<=0)
319 28994509 2004-04-21 devnull error("incomplete positioning");
321 28994509 2004-04-21 devnull case 'y':
322 28994509 2004-04-21 devnull if(argc>0&&!option(argv[0])) {
323 28994509 2004-04-21 devnull symbolfile = argv[0];
328 28994509 2004-04-21 devnull case 'v':
329 28994509 2004-04-21 devnull if(index[k].limb == 0)
330 28994509 2004-04-21 devnull error("-v does not apply here");
331 28994509 2004-04-21 devnull vflag = -1;
333 28994509 2004-04-21 devnull case 'x':
334 28994509 2004-04-21 devnull xflag = 1;
336 28994509 2004-04-21 devnull case 'C':
337 28994509 2004-04-21 devnull if(argc && !option(*argv)) {
338 28994509 2004-04-21 devnull currcolor = colorcode(*argv);
345 28994509 2004-04-21 devnull if(argc>0)
346 28994509 2004-04-21 devnull error("error in arguments");
347 28994509 2004-04-21 devnull pathnames();
348 28994509 2004-04-21 devnull clamp(&limits[0],-90.);
349 28994509 2004-04-21 devnull clamp(&limits[1],90.);
350 28994509 2004-04-21 devnull clamp(&klimits[0],-90.);
351 28994509 2004-04-21 devnull clamp(&klimits[1],90.);
352 28994509 2004-04-21 devnull clamp(&window[0],-90.);
353 28994509 2004-04-21 devnull clamp(&window[1],90.);
354 28994509 2004-04-21 devnull radbds(limits,rlimits);
355 28994509 2004-04-21 devnull limcase = limits[2]<-180.?0:
356 28994509 2004-04-21 devnull limits[3]>180.?2:
359 28994509 2004-04-21 devnull window[0]>=window[1]||
360 28994509 2004-04-21 devnull window[2]>=window[3]||
361 28994509 2004-04-21 devnull window[0]>90.||
362 28994509 2004-04-21 devnull window[1]<-90.||
363 28994509 2004-04-21 devnull window[2]>180.||
364 28994509 2004-04-21 devnull window[3]<-180.)
365 28994509 2004-04-21 devnull error("unreasonable window");
366 28994509 2004-04-21 devnull windlim();
367 28994509 2004-04-21 devnull radbds(window,rwindow);
368 28994509 2004-04-21 devnull upright = orientation[0]==90? 1: orientation[0]==-90? -1: 0;
369 28994509 2004-04-21 devnull if(index[k].spheroid && !upright)
370 28994509 2004-04-21 devnull error("can't tilt the spheroid");
371 28994509 2004-04-21 devnull if(limits[2]>limits[3])
372 28994509 2004-04-21 devnull limits[3] += 360;
373 28994509 2004-04-21 devnull if(!oriented)
374 28994509 2004-04-21 devnull orientation[2] = (limits[2]+limits[3])/2;
375 28994509 2004-04-21 devnull orient(orientation[0],orientation[1],orientation[2]);
376 28994509 2004-04-21 devnull projection = (*index[k].prog)(params[0],params[1]);
377 28994509 2004-04-21 devnull if(projection == 0)
378 28994509 2004-04-21 devnull error("unreasonable projection parameters");
379 28994509 2004-04-21 devnull clipinit();
380 28994509 2004-04-21 devnull grid[0] = fabs(grid[0]);
381 28994509 2004-04-21 devnull grid[1] = fabs(grid[1]);
382 28994509 2004-04-21 devnull if(!kflag)
383 28994509 2004-04-21 devnull for(i=0;i<4;i++)
384 28994509 2004-04-21 devnull klimits[i] = limits[i];
385 28994509 2004-04-21 devnull if(klimits[2]>klimits[3])
386 28994509 2004-04-21 devnull klimits[3] += 360;
387 28994509 2004-04-21 devnull lolat = limits[0];
388 28994509 2004-04-21 devnull hilat = limits[1];
389 28994509 2004-04-21 devnull lolon = limits[2];
390 28994509 2004-04-21 devnull hilon = limits[3];
391 28994509 2004-04-21 devnull if(lolon>=hilon||lolat>=hilat||lolat<-90.||hilat>90.)
392 28994509 2004-04-21 devnull error("unreasonable limits");
393 28994509 2004-04-21 devnull wlim = kflag? klimits: window;
394 28994509 2004-04-21 devnull dlat = fmin(hilat-lolat,wlim[1]-wlim[0])/16;
395 28994509 2004-04-21 devnull dlon = fmin(hilon-lolon,wlim[3]-wlim[2])/32;
396 28994509 2004-04-21 devnull dd = fmax(dlat,dlon);
397 28994509 2004-04-21 devnull while(grid[2]>fmin(dlat,dlon)/2)
398 28994509 2004-04-21 devnull grid[2] /= 2;
399 28994509 2004-04-21 devnull realcut();
400 28994509 2004-04-21 devnull if(nvert<=0) {
401 28994509 2004-04-21 devnull for(lat=klimits[0];lat<klimits[1]+dd-FUZZ;lat+=dd) {
402 28994509 2004-04-21 devnull if(lat>klimits[1])
403 28994509 2004-04-21 devnull lat = klimits[1];
404 28994509 2004-04-21 devnull for(lon=klimits[2];lon<klimits[3]+dd-FUZZ;lon+=dd) {
405 28994509 2004-04-21 devnull i = (kflag?posproj:normproj)
406 28994509 2004-04-21 devnull (lat,lon+(lon<klimits[3]?FUZZ:-FUZZ),
408 28994509 2004-04-21 devnull if(i*vflag <= 0)
409 28994509 2004-04-21 devnull continue;
410 28994509 2004-04-21 devnull if(x<xmin) xmin = x;
411 28994509 2004-04-21 devnull if(x>xmax) xmax = x;
412 28994509 2004-04-21 devnull if(y<ymin) ymin = y;
413 28994509 2004-04-21 devnull if(y>ymax) ymax = y;
416 28994509 2004-04-21 devnull } else {
417 28994509 2004-04-21 devnull for(i=0; i<nvert; i++) {
418 28994509 2004-04-21 devnull x = v[i].x;
419 28994509 2004-04-21 devnull y = v[i].y;
420 28994509 2004-04-21 devnull if(x<xmin) xmin = x;
421 28994509 2004-04-21 devnull if(x>xmax) xmax = x;
422 28994509 2004-04-21 devnull if(y<ymin) ymin = y;
423 28994509 2004-04-21 devnull if(y>ymax) ymax = y;
426 28994509 2004-04-21 devnull xrange = xmax - xmin;
427 28994509 2004-04-21 devnull yrange = ymax - ymin;
428 28994509 2004-04-21 devnull if(xrange<=0||yrange<=0)
429 28994509 2004-04-21 devnull error("map seems to be empty");
430 28994509 2004-04-21 devnull scaling = 2; /*plotting area from -1 to 1*/
431 28994509 2004-04-21 devnull if(position[2]!=0) {
432 28994509 2004-04-21 devnull if(posproj(position[0]-.5,position[1],&xcent,&ycent)==0||
433 28994509 2004-04-21 devnull posproj(position[0]+.5,position[1],&x,&y)==0)
434 28994509 2004-04-21 devnull error("unreasonable position");
435 28994509 2004-04-21 devnull scaling /= (position[2]*hypot(x-xcent,y-ycent));
436 28994509 2004-04-21 devnull if(posproj(position[0],position[1],&xcent,&ycent)==0)
437 28994509 2004-04-21 devnull error("unreasonable position");
438 28994509 2004-04-21 devnull } else {
439 28994509 2004-04-21 devnull scaling /= (xrange>yrange?xrange:yrange);
440 28994509 2004-04-21 devnull xcent = (xmin+xmax)/2;
441 28994509 2004-04-21 devnull ycent = (ymin+ymax)/2;
443 28994509 2004-04-21 devnull xoff = center[0]/scaling;
444 28994509 2004-04-21 devnull yoff = center[1]/scaling;
445 28994509 2004-04-21 devnull crot.l = center[2]*RAD;
446 28994509 2004-04-21 devnull sincos(&crot);
447 28994509 2004-04-21 devnull scaling *= HALFWIDTH*0.9;
448 fa325e9b 2020-01-10 cross if(symbolfile)
449 28994509 2004-04-21 devnull getsyms(symbolfile);
450 28994509 2004-04-21 devnull if(!s2flag) {
451 28994509 2004-04-21 devnull openpl();
452 28994509 2004-04-21 devnull erase();
454 28994509 2004-04-21 devnull range(left,bottom,right,top);
455 28994509 2004-04-21 devnull comment("grid","");
456 28994509 2004-04-21 devnull colorx(gridcolor);
457 28994509 2004-04-21 devnull pen(DOTTED);
458 28994509 2004-04-21 devnull if(grid[0]>0.)
459 28994509 2004-04-21 devnull for(lat=ceil(lolat/grid[0])*grid[0];
460 fa325e9b 2020-01-10 cross lat<=hilat;lat+=grid[0])
461 28994509 2004-04-21 devnull dogrid(lat,lat,lolon,hilon);
462 28994509 2004-04-21 devnull if(grid[1]>0.)
463 28994509 2004-04-21 devnull for(lon=ceil(lolon/grid[1])*grid[1];
464 fa325e9b 2020-01-10 cross lon<=hilon;lon+=grid[1])
465 28994509 2004-04-21 devnull dogrid(lolat,hilat,lon,lon);
466 28994509 2004-04-21 devnull comment("border","");
467 28994509 2004-04-21 devnull colorx(bordcolor);
468 28994509 2004-04-21 devnull pen(SOLID);
469 28994509 2004-04-21 devnull if(bflag) {
470 28994509 2004-04-21 devnull dolimb();
471 28994509 2004-04-21 devnull dobounds(lolat,hilat,lolon,hilon,0);
472 28994509 2004-04-21 devnull dobounds(window[0],window[1],window[2],window[3],1);
474 28994509 2004-04-21 devnull lolat = floor(limits[0]/10)*10;
475 28994509 2004-04-21 devnull hilat = ceil(limits[1]/10)*10;
476 28994509 2004-04-21 devnull lolon = floor(limits[2]/10)*10;
477 28994509 2004-04-21 devnull hilon = ceil(limits[3]/10)*10;
478 28994509 2004-04-21 devnull if(lolon>hilon)
479 28994509 2004-04-21 devnull hilon += 360.;
480 28994509 2004-04-21 devnull /*do tracks first so as not to lose the standard input*/
481 28994509 2004-04-21 devnull for(i=0;i<ntrack;i++) {
482 28994509 2004-04-21 devnull longlines = LONGLINES;
483 28994509 2004-04-21 devnull satellite(&track[i]);
484 28994509 2004-04-21 devnull longlines = shortlines;
486 28994509 2004-04-21 devnull for(i=0;i<nfile;i++) {
487 28994509 2004-04-21 devnull comment("mapfile",file[i].name);
488 28994509 2004-04-21 devnull colorx(file[i].color);
489 28994509 2004-04-21 devnull pen(file[i].style);
490 28994509 2004-04-21 devnull getdata(file[i].name);
492 28994509 2004-04-21 devnull move(right,bottom);
493 28994509 2004-04-21 devnull if(!s1flag)
494 28994509 2004-04-21 devnull closepl();
495 28994509 2004-04-21 devnull return 0;
498 28994509 2004-04-21 devnull /* Out of perverseness (really to recover from a dubious,
499 28994509 2004-04-21 devnull but documented, convention) the returns from projection
500 28994509 2004-04-21 devnull functions (-1 unplottable, 0 wrong sheet, 1 good) are
501 28994509 2004-04-21 devnull recoded into -1 wrong sheet, 0 unplottable, 1 good. */
504 28994509 2004-04-21 devnull fixproj(struct place *g, double *x, double *y)
506 28994509 2004-04-21 devnull int i = (*projection)(g,x,y);
507 28994509 2004-04-21 devnull return i<0? 0: i==0? -1: 1;
511 28994509 2004-04-21 devnull normproj(double lat, double lon, double *x, double *y)
514 28994509 2004-04-21 devnull struct place geog;
515 28994509 2004-04-21 devnull latlon(lat,lon,&geog);
517 28994509 2004-04-21 devnull printp(&geog);
519 28994509 2004-04-21 devnull normalize(&geog);
520 28994509 2004-04-21 devnull if(!inwindow(&geog))
521 28994509 2004-04-21 devnull return(-1);
522 28994509 2004-04-21 devnull i = fixproj(&geog,x,y);
524 28994509 2004-04-21 devnull *x = -*x;
526 28994509 2004-04-21 devnull printp(&geog);
527 28994509 2004-04-21 devnull fprintf(stderr,"%d %.3f %.3f\n",i,*x,*y);
529 28994509 2004-04-21 devnull return(i);
533 28994509 2004-04-21 devnull posproj(double lat, double lon, double *x, double *y)
536 28994509 2004-04-21 devnull struct place geog;
537 28994509 2004-04-21 devnull latlon(lat,lon,&geog);
538 28994509 2004-04-21 devnull normalize(&geog);
539 28994509 2004-04-21 devnull i = fixproj(&geog,x,y);
541 28994509 2004-04-21 devnull *x = -*x;
542 28994509 2004-04-21 devnull return(i);
546 28994509 2004-04-21 devnull inwindow(struct place *geog)
548 28994509 2004-04-21 devnull if(geog->nlat.l<rwindow[0]-FUZZ||
549 28994509 2004-04-21 devnull geog->nlat.l>rwindow[1]+FUZZ||
550 28994509 2004-04-21 devnull geog->wlon.l<rwindow[2]-FUZZ||
551 28994509 2004-04-21 devnull geog->wlon.l>rwindow[3]+FUZZ)
552 28994509 2004-04-21 devnull return(0);
553 28994509 2004-04-21 devnull else return(1);
557 28994509 2004-04-21 devnull inlimits(struct place *g)
559 28994509 2004-04-21 devnull if(rlimits[0]-FUZZ>g->nlat.l||
560 28994509 2004-04-21 devnull rlimits[1]+FUZZ<g->nlat.l)
561 28994509 2004-04-21 devnull return(0);
562 28994509 2004-04-21 devnull switch(limcase) {
564 28994509 2004-04-21 devnull if(rlimits[2]+TWOPI-FUZZ>g->wlon.l&&
565 28994509 2004-04-21 devnull rlimits[3]+FUZZ<g->wlon.l)
566 28994509 2004-04-21 devnull return(0);
569 28994509 2004-04-21 devnull if(rlimits[2]-FUZZ>g->wlon.l||
570 28994509 2004-04-21 devnull rlimits[3]+FUZZ<g->wlon.l)
571 28994509 2004-04-21 devnull return(0);
574 28994509 2004-04-21 devnull if(rlimits[2]>g->wlon.l&&
575 28994509 2004-04-21 devnull rlimits[3]-TWOPI+FUZZ<g->wlon.l)
576 28994509 2004-04-21 devnull return(0);
579 28994509 2004-04-21 devnull return(1);
583 28994509 2004-04-21 devnull long patch[18][36];
586 28994509 2004-04-21 devnull getdata(char *mapfile)
588 28994509 2004-04-21 devnull char *indexfile;
589 28994509 2004-04-21 devnull int kx,ky,c;
592 28994509 2004-04-21 devnull long *p;
593 28994509 2004-04-21 devnull int ip, jp;
595 28994509 2004-04-21 devnull struct place g;
596 28994509 2004-04-21 devnull int i, j;
597 28994509 2004-04-21 devnull double lat, lon;
598 28994509 2004-04-21 devnull int conn;
599 28994509 2004-04-21 devnull FILE *ifile, *xfile;
601 28994509 2004-04-21 devnull indexfile = mapindex(mapfile);
602 28994509 2004-04-21 devnull xfile = fopen(indexfile,"r");
603 28994509 2004-04-21 devnull if(xfile==NULL)
604 28994509 2004-04-21 devnull filerror("can't find map index", indexfile);
605 28994509 2004-04-21 devnull free(indexfile);
606 28994509 2004-04-21 devnull for(i=0,p=patch[0];i<18*36;i++,p++)
608 28994509 2004-04-21 devnull while(!feof(xfile) && fscanf(xfile,"%d%d%ld",&i,&j,&b)==3)
609 28994509 2004-04-21 devnull patch[i+9][j+18] = b;
610 28994509 2004-04-21 devnull fclose(xfile);
611 28994509 2004-04-21 devnull ifile = fopen(mapfile,"r");
612 28994509 2004-04-21 devnull if(ifile==NULL)
613 28994509 2004-04-21 devnull filerror("can't find map data", mapfile);
614 28994509 2004-04-21 devnull for(lat=lolat;lat<hilat;lat+=10.)
615 28994509 2004-04-21 devnull for(lon=lolon;lon<hilon;lon+=10.) {
616 28994509 2004-04-21 devnull if(!seeable(lat,lon))
617 28994509 2004-04-21 devnull continue;
618 28994509 2004-04-21 devnull i = pnorm(lat);
619 28994509 2004-04-21 devnull j = pnorm(lon);
620 28994509 2004-04-21 devnull if((b=patch[i+9][j+18])&1)
621 28994509 2004-04-21 devnull continue;
622 28994509 2004-04-21 devnull fseek(ifile,b,0);
623 28994509 2004-04-21 devnull while((ip=getc(ifile))>=0&&(jp=getc(ifile))>=0){
624 28994509 2004-04-21 devnull if(ip!=(i&0377)||jp!=(j&0377))
626 28994509 2004-04-21 devnull n = getshort(ifile);
627 28994509 2004-04-21 devnull conn = 0;
628 28994509 2004-04-21 devnull if(n > 0) { /* absolute coordinates */
629 28994509 2004-04-21 devnull kx = ky = 0; /* set */
630 28994509 2004-04-21 devnull for(k=0;k<n;k++){
631 28994509 2004-04-21 devnull kx = SCALERATIO*getshort(ifile);
632 28994509 2004-04-21 devnull ky = SCALERATIO*getshort(ifile);
633 28994509 2004-04-21 devnull if (((k%delta) != 0) && (k != (n-1)))
634 28994509 2004-04-21 devnull continue;
635 28994509 2004-04-21 devnull conv(kx,&g.nlat);
636 28994509 2004-04-21 devnull conv(ky,&g.wlon);
637 28994509 2004-04-21 devnull conn = plotpt(&g,conn);
639 28994509 2004-04-21 devnull } else { /* differential, scaled by SCALERATI0 */
641 28994509 2004-04-21 devnull kx = SCALERATIO*getshort(ifile);
642 28994509 2004-04-21 devnull ky = SCALERATIO*getshort(ifile);
643 28994509 2004-04-21 devnull for(k=0; k<n; k++) {
644 28994509 2004-04-21 devnull c = getc(ifile);
645 28994509 2004-04-21 devnull if(c&0200) c|= ~0177;
646 28994509 2004-04-21 devnull kx += c;
647 28994509 2004-04-21 devnull c = getc(ifile);
648 28994509 2004-04-21 devnull if(c&0200) c|= ~0177;
649 28994509 2004-04-21 devnull ky += c;
650 28994509 2004-04-21 devnull if(k%delta!=0&&k!=n-1)
651 28994509 2004-04-21 devnull continue;
652 28994509 2004-04-21 devnull conv(kx,&g.nlat);
653 28994509 2004-04-21 devnull conv(ky,&g.wlon);
654 28994509 2004-04-21 devnull conn = plotpt(&g,conn);
657 28994509 2004-04-21 devnull if(k==1) {
658 28994509 2004-04-21 devnull conv(kx,&g.nlat);
659 28994509 2004-04-21 devnull conv(ky,&g.wlon);
660 28994509 2004-04-21 devnull plotpt(&g,conn);
664 28994509 2004-04-21 devnull fclose(ifile);
668 28994509 2004-04-21 devnull seeable(double lat0, double lon0)
670 28994509 2004-04-21 devnull double x, y;
671 28994509 2004-04-21 devnull double lat, lon;
672 28994509 2004-04-21 devnull for(lat=lat0;lat<=lat0+10;lat+=2*grid[2])
673 28994509 2004-04-21 devnull for(lon=lon0;lon<=lon0+10;lon+=2*grid[2])
674 28994509 2004-04-21 devnull if(normproj(lat,lon,&x,&y)*vflag>0)
675 28994509 2004-04-21 devnull return(1);
676 28994509 2004-04-21 devnull return(0);
680 28994509 2004-04-21 devnull satellite(struct file *t)
682 28994509 2004-04-21 devnull char sym[50];
683 28994509 2004-04-21 devnull char lbl[50];
684 28994509 2004-04-21 devnull double scale;
685 28994509 2004-04-21 devnull int conn;
686 28994509 2004-04-21 devnull double lat,lon;
687 28994509 2004-04-21 devnull struct place place;
688 28994509 2004-04-21 devnull static FILE *ifile;
690 28994509 2004-04-21 devnull if(ifile == nil)
691 28994509 2004-04-21 devnull ifile = stdin;
692 28994509 2004-04-21 devnull if(t->name[0]!='-'||t->name[1]!=0) {
693 28994509 2004-04-21 devnull fclose(ifile);
694 28994509 2004-04-21 devnull if((ifile=fopen(t->name,"r"))==NULL)
695 28994509 2004-04-21 devnull filerror("can't find track", t->name);
697 28994509 2004-04-21 devnull comment("track",t->name);
698 28994509 2004-04-21 devnull colorx(t->color);
699 28994509 2004-04-21 devnull pen(t->style);
700 28994509 2004-04-21 devnull for(;;) {
701 28994509 2004-04-21 devnull conn = 0;
702 28994509 2004-04-21 devnull while(!feof(ifile) && fscanf(ifile,"%lf%lf",&lat,&lon)==2){
703 28994509 2004-04-21 devnull latlon(lat,lon,&place);
704 28994509 2004-04-21 devnull if(fscanf(ifile,"%1s",lbl) == 1) {
705 28994509 2004-04-21 devnull if(strchr("+-.0123456789",*lbl)==0)
707 28994509 2004-04-21 devnull ungetc(*lbl,ifile);
709 28994509 2004-04-21 devnull conn = plotpt(&place,conn);
711 28994509 2004-04-21 devnull if(feof(ifile))
713 28994509 2004-04-21 devnull fscanf(ifile,"%[^\n]",lbl+1);
714 28994509 2004-04-21 devnull switch(*lbl) {
715 28994509 2004-04-21 devnull case '"':
716 28994509 2004-04-21 devnull if(plotpt(&place,conn))
717 28994509 2004-04-21 devnull text(lbl+1);
719 28994509 2004-04-21 devnull case ':':
720 28994509 2004-04-21 devnull case '!':
721 28994509 2004-04-21 devnull if(sscanf(lbl+1,"%s %lf",sym,&scale) <= 1)
722 28994509 2004-04-21 devnull scale = 1;
723 28994509 2004-04-21 devnull if(plotpt(&place,conn?conn:-1)) {
724 28994509 2004-04-21 devnull int r = *lbl=='!'?0:rflag?-1:1;
725 28994509 2004-04-21 devnull pen(SOLID);
726 28994509 2004-04-21 devnull if(putsym(&place,sym,scale,r) == 0)
727 28994509 2004-04-21 devnull text(lbl);
728 28994509 2004-04-21 devnull pen(t->style);
731 28994509 2004-04-21 devnull default:
732 28994509 2004-04-21 devnull if(plotpt(&place,conn))
733 28994509 2004-04-21 devnull text(lbl);
740 28994509 2004-04-21 devnull pnorm(double x)
743 28994509 2004-04-21 devnull i = x/10.;
744 28994509 2004-04-21 devnull i %= 36;
745 28994509 2004-04-21 devnull if(i>=18) return(i-36);
746 28994509 2004-04-21 devnull if(i<-18) return(i+36);
747 28994509 2004-04-21 devnull return(i);
751 28994509 2004-04-21 devnull error(char *s)
753 28994509 2004-04-21 devnull fprintf(stderr,"map: \r\n%s\n",s);
754 28994509 2004-04-21 devnull exits("error");
758 28994509 2004-04-21 devnull filerror(char *s, char *f)
760 28994509 2004-04-21 devnull fprintf(stderr,"\r\n%s %s\n",s,f);
761 28994509 2004-04-21 devnull exits("error");
765 28994509 2004-04-21 devnull mapindex(char *s)
767 28994509 2004-04-21 devnull char *t = malloc(strlen(s)+3);
768 28994509 2004-04-21 devnull strcpy(t,s);
769 28994509 2004-04-21 devnull strcat(t,".x");
770 28994509 2004-04-21 devnull return t;
773 28994509 2004-04-21 devnull #define NOPT 32767
774 28994509 2004-04-21 devnull static int ox = NOPT;
775 28994509 2004-04-21 devnull static int oy = NOPT;
778 28994509 2004-04-21 devnull cpoint(int xi, int yi, int conn)
780 28994509 2004-04-21 devnull int dx = abs(ox-xi);
781 28994509 2004-04-21 devnull int dy = abs(oy-yi);
782 28994509 2004-04-21 devnull if(!xflag && (xi<left||xi>=right || yi<bottom||yi>=top)) {
783 28994509 2004-04-21 devnull ox = oy = NOPT;
784 28994509 2004-04-21 devnull return 0;
786 28994509 2004-04-21 devnull if(conn == -1) /* isolated plotting symbol */
788 28994509 2004-04-21 devnull else if(!conn)
789 28994509 2004-04-21 devnull point(xi,yi);
791 28994509 2004-04-21 devnull if(dx+dy>longlines) {
792 28994509 2004-04-21 devnull ox = oy = NOPT; /* don't leap across cuts */
793 28994509 2004-04-21 devnull return 0;
795 28994509 2004-04-21 devnull if(dx || dy)
796 28994509 2004-04-21 devnull vec(xi,yi);
798 28994509 2004-04-21 devnull ox = xi, oy = yi;
799 28994509 2004-04-21 devnull return dx+dy<=2? 2: 1; /* 2=very near; see dogrid */
803 28994509 2004-04-21 devnull struct place oldg;
806 28994509 2004-04-21 devnull plotpt(struct place *g, int conn)
808 28994509 2004-04-21 devnull int kx,ky;
809 28994509 2004-04-21 devnull int ret;
810 28994509 2004-04-21 devnull double cutlon;
811 28994509 2004-04-21 devnull if(!inlimits(g)) {
812 28994509 2004-04-21 devnull return(0);
814 28994509 2004-04-21 devnull normalize(g);
815 28994509 2004-04-21 devnull if(!inwindow(g)) {
816 28994509 2004-04-21 devnull return(0);
818 28994509 2004-04-21 devnull switch((*cut)(g,&oldg,&cutlon)) {
820 28994509 2004-04-21 devnull if(conn) {
821 28994509 2004-04-21 devnull ret = duple(g,cutlon)|duple(g,cutlon);
822 28994509 2004-04-21 devnull oldg = *g;
823 28994509 2004-04-21 devnull return(ret);
826 28994509 2004-04-21 devnull conn = 0;
827 28994509 2004-04-21 devnull default: /* prevent diags about bad return value */
829 28994509 2004-04-21 devnull oldg = *g;
830 28994509 2004-04-21 devnull ret = doproj(g,&kx,&ky);
831 28994509 2004-04-21 devnull if(ret==0 || !onlimb && ret*vflag<=0)
832 28994509 2004-04-21 devnull return(0);
833 28994509 2004-04-21 devnull ret = cpoint(kx,ky,conn);
834 28994509 2004-04-21 devnull return ret;
839 28994509 2004-04-21 devnull doproj(struct place *g, int *kx, int *ky)
842 28994509 2004-04-21 devnull double x,y,x1,y1;
843 28994509 2004-04-21 devnull /*fprintf(stderr,"dopr1 %f %f \n",g->nlat.l,g->wlon.l);*/
844 28994509 2004-04-21 devnull i = fixproj(g,&x,&y);
845 28994509 2004-04-21 devnull if(i == 0)
846 28994509 2004-04-21 devnull return(0);
847 28994509 2004-04-21 devnull if(rflag)
849 28994509 2004-04-21 devnull /*fprintf(stderr,"dopr2 %f %f\n",x,y);*/
850 28994509 2004-04-21 devnull if(!inpoly(x,y)) {
851 28994509 2004-04-21 devnull return 0;
853 28994509 2004-04-21 devnull x1 = x - xcent;
854 28994509 2004-04-21 devnull y1 = y - ycent;
855 28994509 2004-04-21 devnull x = (x1*crot.c - y1*crot.s + xoff)*scaling;
856 28994509 2004-04-21 devnull y = (x1*crot.s + y1*crot.c + yoff)*scaling;
857 28994509 2004-04-21 devnull *kx = x + (x>0?.5:-.5);
858 28994509 2004-04-21 devnull *ky = y + (y>0?.5:-.5);
859 28994509 2004-04-21 devnull return(i);
863 28994509 2004-04-21 devnull duple(struct place *g, double cutlon)
865 28994509 2004-04-21 devnull int kx,ky;
866 28994509 2004-04-21 devnull int okx,oky;
867 28994509 2004-04-21 devnull struct place ig;
868 28994509 2004-04-21 devnull revlon(g,cutlon);
869 28994509 2004-04-21 devnull revlon(&oldg,cutlon);
870 28994509 2004-04-21 devnull ig = *g;
871 28994509 2004-04-21 devnull invert(&ig);
872 28994509 2004-04-21 devnull if(!inlimits(&ig))
873 28994509 2004-04-21 devnull return(0);
874 28994509 2004-04-21 devnull if(doproj(g,&kx,&ky)*vflag<=0 ||
875 28994509 2004-04-21 devnull doproj(&oldg,&okx,&oky)*vflag<=0)
876 28994509 2004-04-21 devnull return(0);
877 28994509 2004-04-21 devnull cpoint(okx,oky,0);
878 28994509 2004-04-21 devnull cpoint(kx,ky,1);
879 28994509 2004-04-21 devnull return(1);
883 28994509 2004-04-21 devnull revlon(struct place *g, double cutlon)
885 28994509 2004-04-21 devnull g->wlon.l = reduce(cutlon-reduce(g->wlon.l-cutlon));
886 28994509 2004-04-21 devnull sincos(&g->wlon);
890 28994509 2004-04-21 devnull /* recognize problems of cuts
891 28994509 2004-04-21 devnull * move a point across cut to side of its predecessor
892 28994509 2004-04-21 devnull * if its very close to the cut
893 28994509 2004-04-21 devnull * return(0) if cut interrupts the line
894 28994509 2004-04-21 devnull * return(1) if line is to be drawn normally
895 28994509 2004-04-21 devnull * return(2) if line is so close to cut as to
896 28994509 2004-04-21 devnull * be properly drawn on both sheets
900 28994509 2004-04-21 devnull picut(struct place *g, struct place *og, double *cutlon)
902 28994509 2004-04-21 devnull *cutlon = PI;
903 28994509 2004-04-21 devnull return(ckcut(g,og,PI));
907 28994509 2004-04-21 devnull nocut(struct place *g, struct place *og, double *cutlon)
909 28994509 2004-04-21 devnull USED(g);
910 28994509 2004-04-21 devnull USED(og);
911 28994509 2004-04-21 devnull USED(cutlon);
913 28994509 2004-04-21 devnull #pragma ref g
914 28994509 2004-04-21 devnull #pragma ref og
915 28994509 2004-04-21 devnull #pragma ref cutlon
917 28994509 2004-04-21 devnull return(1);
921 28994509 2004-04-21 devnull ckcut(struct place *g1, struct place *g2, double lon)
923 28994509 2004-04-21 devnull double d1, d2;
924 28994509 2004-04-21 devnull double f1, f2;
925 28994509 2004-04-21 devnull int kx,ky;
926 28994509 2004-04-21 devnull d1 = reduce(g1->wlon.l -lon);
927 28994509 2004-04-21 devnull d2 = reduce(g2->wlon.l -lon);
928 28994509 2004-04-21 devnull if((f1=fabs(d1))<FUZZ)
929 28994509 2004-04-21 devnull d1 = diddle(g1,lon,d2);
930 28994509 2004-04-21 devnull if((f2=fabs(d2))<FUZZ) {
931 28994509 2004-04-21 devnull d2 = diddle(g2,lon,d1);
932 28994509 2004-04-21 devnull if(doproj(g2,&kx,&ky)*vflag>0)
933 28994509 2004-04-21 devnull cpoint(kx,ky,0);
935 28994509 2004-04-21 devnull if(f1<FUZZ&&f2<FUZZ)
936 28994509 2004-04-21 devnull return(2);
937 28994509 2004-04-21 devnull if(f1>PI*TWO_THRD||f2>PI*TWO_THRD)
938 28994509 2004-04-21 devnull return(1);
939 28994509 2004-04-21 devnull return(d1*d2>=0);
943 28994509 2004-04-21 devnull diddle(struct place *g, double lon, double d)
945 28994509 2004-04-21 devnull double d1;
946 28994509 2004-04-21 devnull d1 = FUZZ/2;
948 28994509 2004-04-21 devnull d1 = -d1;
949 28994509 2004-04-21 devnull g->wlon.l = reduce(lon+d1);
950 28994509 2004-04-21 devnull sincos(&g->wlon);
951 28994509 2004-04-21 devnull return(d1);
955 28994509 2004-04-21 devnull reduce(double lon)
957 28994509 2004-04-21 devnull if(lon>PI)
958 28994509 2004-04-21 devnull lon -= 2*PI;
959 28994509 2004-04-21 devnull else if(lon<-PI)
960 28994509 2004-04-21 devnull lon += 2*PI;
961 28994509 2004-04-21 devnull return(lon);
965 28994509 2004-04-21 devnull double tetrapt = 35.26438968; /* atan(1/sqrt(2)) */
968 28994509 2004-04-21 devnull dogrid(double lat0, double lat1, double lon0, double lon1)
970 28994509 2004-04-21 devnull double slat,slon,tlat,tlon;
971 28994509 2004-04-21 devnull register int conn, oconn;
972 28994509 2004-04-21 devnull slat = tlat = slon = tlon = 0;
973 28994509 2004-04-21 devnull if(lat1>lat0)
974 28994509 2004-04-21 devnull slat = tlat = fmin(grid[2],dlat);
976 28994509 2004-04-21 devnull slon = tlon = fmin(grid[2],dlon);;
977 28994509 2004-04-21 devnull conn = oconn = 0;
978 28994509 2004-04-21 devnull while(lat0<=lat1&&lon0<=lon1) {
979 28994509 2004-04-21 devnull conn = gridpt(lat0,lon0,conn);
980 28994509 2004-04-21 devnull if(projection==Xguyou&&slat>0) {
981 28994509 2004-04-21 devnull if(lat0<-45&&lat0+slat>-45)
982 28994509 2004-04-21 devnull conn = gridpt(-45.,lon0,conn);
983 28994509 2004-04-21 devnull else if(lat0<45&&lat0+slat>45)
984 28994509 2004-04-21 devnull conn = gridpt(45.,lon0,conn);
985 28994509 2004-04-21 devnull } else if(projection==Xtetra&&slat>0) {
986 28994509 2004-04-21 devnull if(lat0<-tetrapt&&lat0+slat>-tetrapt) {
987 28994509 2004-04-21 devnull gridpt(-tetrapt-.001,lon0,conn);
988 28994509 2004-04-21 devnull conn = gridpt(-tetrapt+.001,lon0,0);
990 28994509 2004-04-21 devnull else if(lat0<tetrapt&&lat0+slat>tetrapt) {
991 28994509 2004-04-21 devnull gridpt(tetrapt-.001,lon0,conn);
992 28994509 2004-04-21 devnull conn = gridpt(tetrapt+.001,lon0,0);
995 28994509 2004-04-21 devnull if(conn==0 && oconn!=0) {
996 28994509 2004-04-21 devnull if(slat+slon>.05) {
997 28994509 2004-04-21 devnull lat0 -= slat; /* steps too big */
998 28994509 2004-04-21 devnull lon0 -= slon; /* or near bdry */
999 28994509 2004-04-21 devnull slat /= 2;
1000 28994509 2004-04-21 devnull slon /= 2;
1001 28994509 2004-04-21 devnull conn = oconn = gridpt(lat0,lon0,conn);
1003 28994509 2004-04-21 devnull oconn = 0;
1004 28994509 2004-04-21 devnull } else {
1005 28994509 2004-04-21 devnull if(conn==2) {
1006 28994509 2004-04-21 devnull slat = tlat;
1007 28994509 2004-04-21 devnull slon = tlon;
1008 28994509 2004-04-21 devnull conn = 1;
1010 28994509 2004-04-21 devnull oconn = conn;
1012 28994509 2004-04-21 devnull lat0 += slat;
1013 28994509 2004-04-21 devnull lon0 += slon;
1015 28994509 2004-04-21 devnull gridpt(lat1,lon1,conn);
1018 28994509 2004-04-21 devnull static int gridinv; /* nonzero when doing window bounds */
1021 28994509 2004-04-21 devnull gridpt(double lat, double lon, int conn)
1023 28994509 2004-04-21 devnull struct place g;
1024 28994509 2004-04-21 devnull /*fprintf(stderr,"%f %f\n",lat,lon);*/
1025 28994509 2004-04-21 devnull latlon(lat,lon,&g);
1026 28994509 2004-04-21 devnull if(gridinv)
1027 28994509 2004-04-21 devnull invert(&g);
1028 28994509 2004-04-21 devnull return(plotpt(&g,conn));
1031 28994509 2004-04-21 devnull /* win=0 ordinary grid lines, win=1 window lines */
1034 28994509 2004-04-21 devnull dobounds(double lolat, double hilat, double lolon, double hilon, int win)
1036 28994509 2004-04-21 devnull gridinv = win;
1037 28994509 2004-04-21 devnull if(lolat>-90 || win && (poles&1)!=0)
1038 28994509 2004-04-21 devnull dogrid(lolat+FUZZ,lolat+FUZZ,lolon,hilon);
1039 28994509 2004-04-21 devnull if(hilat<90 || win && (poles&2)!=0)
1040 28994509 2004-04-21 devnull dogrid(hilat-FUZZ,hilat-FUZZ,lolon,hilon);
1041 28994509 2004-04-21 devnull if(hilon-lolon<360 || win && cut==picut) {
1042 28994509 2004-04-21 devnull dogrid(lolat,hilat,lolon+FUZZ,lolon+FUZZ);
1043 28994509 2004-04-21 devnull dogrid(lolat,hilat,hilon-FUZZ,hilon-FUZZ);
1045 28994509 2004-04-21 devnull gridinv = 0;
1048 28994509 2004-04-21 devnull static void
1049 28994509 2004-04-21 devnull dolimb(void)
1051 28994509 2004-04-21 devnull double lat, lon;
1052 28994509 2004-04-21 devnull double res = fmin(dlat, dlon)/4;
1053 28994509 2004-04-21 devnull int conn = 0;
1054 28994509 2004-04-21 devnull int newconn;
1055 28994509 2004-04-21 devnull if(limb == 0)
1056 28994509 2004-04-21 devnull return;
1057 28994509 2004-04-21 devnull onlimb = gridinv = 1;
1058 28994509 2004-04-21 devnull for(;;) {
1059 28994509 2004-04-21 devnull newconn = (*limb)(&lat, &lon, res);
1060 28994509 2004-04-21 devnull if(newconn == -1)
1062 28994509 2004-04-21 devnull conn = gridpt(lat, lon, conn*newconn);
1064 28994509 2004-04-21 devnull onlimb = gridinv = 0;
1069 28994509 2004-04-21 devnull radbds(double *w, double *rw)
1072 28994509 2004-04-21 devnull for(i=0;i<4;i++)
1073 28994509 2004-04-21 devnull rw[i] = w[i]*RAD;
1074 28994509 2004-04-21 devnull rw[0] -= FUZZ;
1075 28994509 2004-04-21 devnull rw[1] += FUZZ;
1076 28994509 2004-04-21 devnull rw[2] -= FUZZ;
1077 28994509 2004-04-21 devnull rw[3] += FUZZ;
1081 28994509 2004-04-21 devnull windlim(void)
1083 28994509 2004-04-21 devnull double center = orientation[0];
1084 28994509 2004-04-21 devnull double colat;
1085 28994509 2004-04-21 devnull if(center>90)
1086 28994509 2004-04-21 devnull center = 180 - center;
1087 28994509 2004-04-21 devnull if(center<-90)
1088 28994509 2004-04-21 devnull center = -180 - center;
1089 28994509 2004-04-21 devnull if(fabs(center)>90)
1090 28994509 2004-04-21 devnull error("unreasonable orientation");
1091 28994509 2004-04-21 devnull colat = 90 - window[0];
1092 28994509 2004-04-21 devnull if(center-colat>limits[0])
1093 28994509 2004-04-21 devnull limits[0] = center - colat;
1094 28994509 2004-04-21 devnull if(center+colat<limits[1])
1095 28994509 2004-04-21 devnull limits[1] = center + colat;
1100 28994509 2004-04-21 devnull getshort(FILE *f)
1102 28994509 2004-04-21 devnull int c, r;
1103 28994509 2004-04-21 devnull c = getc(f);
1104 28994509 2004-04-21 devnull r = (c | getc(f)<<8);
1105 28994509 2004-04-21 devnull if (r&0x8000)
1106 28994509 2004-04-21 devnull r |= ~0xFFFF; /* in case short > 16 bits */
1107 28994509 2004-04-21 devnull return r;
1111 28994509 2004-04-21 devnull fmin(double x, double y)
1113 28994509 2004-04-21 devnull return(x<y?x:y);
1117 28994509 2004-04-21 devnull fmax(double x, double y)
1119 28994509 2004-04-21 devnull return(x>y?x:y);
1123 28994509 2004-04-21 devnull clamp(double *px, double v)
1125 28994509 2004-04-21 devnull *px = (v<0?fmax:fmin)(*px,v);
1129 28994509 2004-04-21 devnull pathnames(void)
1132 28994509 2004-04-21 devnull char *t, *indexfile, *name;
1133 28994509 2004-04-21 devnull FILE *f, *fx;
1134 28994509 2004-04-21 devnull for(i=0; i<nfile; i++) {
1135 28994509 2004-04-21 devnull name = file[i].name;
1136 28994509 2004-04-21 devnull if(*name=='/')
1137 28994509 2004-04-21 devnull continue;
1138 28994509 2004-04-21 devnull indexfile = mapindex(name);
1139 28994509 2004-04-21 devnull /* ansi equiv of unix access() call */
1140 28994509 2004-04-21 devnull f = fopen(name, "r");
1141 28994509 2004-04-21 devnull fx = fopen(indexfile, "r");
1142 28994509 2004-04-21 devnull if(f) fclose(f);
1143 28994509 2004-04-21 devnull if(fx) fclose(fx);
1144 28994509 2004-04-21 devnull free(indexfile);
1145 28994509 2004-04-21 devnull if(f && fx)
1146 28994509 2004-04-21 devnull continue;
1147 28994509 2004-04-21 devnull t = malloc(strlen(name)+strlen(mapdir)+2);
1148 28994509 2004-04-21 devnull strcpy(t,mapdir);
1149 28994509 2004-04-21 devnull strcat(t,"/");
1150 28994509 2004-04-21 devnull strcat(t,name);
1151 28994509 2004-04-21 devnull file[i].name = t;
1156 28994509 2004-04-21 devnull clipinit(void)
1159 28994509 2004-04-21 devnull double s,t;
1160 28994509 2004-04-21 devnull if(nvert<=0)
1161 28994509 2004-04-21 devnull return;
1162 28994509 2004-04-21 devnull for(i=0; i<nvert; i++) { /*convert latlon to xy*/
1163 28994509 2004-04-21 devnull if(normproj(v[i].x,v[i].y,&v[i].x,&v[i].y)==0)
1164 28994509 2004-04-21 devnull error("invisible clipping vertex");
1166 28994509 2004-04-21 devnull if(nvert==2) { /*rectangle with diag specified*/
1167 28994509 2004-04-21 devnull nvert = 4;
1168 28994509 2004-04-21 devnull v[2] = v[1];
1169 28994509 2004-04-21 devnull v[1].x=v[0].x, v[1].y=v[2].y, v[3].x=v[2].x, v[3].y=v[0].y;
1171 28994509 2004-04-21 devnull v[nvert] = v[0];
1172 28994509 2004-04-21 devnull v[nvert+1] = v[1];
1174 28994509 2004-04-21 devnull for(i=1; i<=nvert; i++) { /*test for convexity*/
1175 28994509 2004-04-21 devnull t = (v[i-1].x-v[i].x)*(v[i+1].y-v[i].y) -
1176 28994509 2004-04-21 devnull (v[i-1].y-v[i].y)*(v[i+1].x-v[i].x);
1177 28994509 2004-04-21 devnull if(t<-FUZZ && s>=0) s = 1;
1178 28994509 2004-04-21 devnull if(t>FUZZ && s<=0) s = -1;
1179 28994509 2004-04-21 devnull if(-FUZZ<=t&&t<=FUZZ || t*s>0) {
1184 28994509 2004-04-21 devnull if(s==0)
1185 28994509 2004-04-21 devnull error("improper clipping polygon");
1186 28994509 2004-04-21 devnull for(i=0; i<nvert; i++) { /*edge equation ax+by=c*/
1187 28994509 2004-04-21 devnull e[i].a = s*(v[i+1].y - v[i].y);
1188 28994509 2004-04-21 devnull e[i].b = s*(v[i].x - v[i+1].x);
1189 28994509 2004-04-21 devnull e[i].c = s*(v[i].x*v[i+1].y - v[i].y*v[i+1].x);
1194 28994509 2004-04-21 devnull inpoly(double x, double y)
1197 28994509 2004-04-21 devnull for(i=0; i<nvert; i++) {
1198 28994509 2004-04-21 devnull register struct edge *ei = &e[i];
1199 28994509 2004-04-21 devnull double val = x*ei->a + y*ei->b - ei->c;
1200 28994509 2004-04-21 devnull if(val>10*FUZZ)
1201 28994509 2004-04-21 devnull return(0);
1203 28994509 2004-04-21 devnull return 1;
1207 2b604081 2005-05-07 devnull realcut(void)
1209 28994509 2004-04-21 devnull struct place g;
1210 28994509 2004-04-21 devnull double lat;
1212 28994509 2004-04-21 devnull if(cut != picut) /* punt on unusual cuts */
1213 28994509 2004-04-21 devnull return;
1214 28994509 2004-04-21 devnull for(lat=window[0]; lat<=window[1]; lat+=grid[2]) {
1215 28994509 2004-04-21 devnull g.wlon.l = PI;
1216 28994509 2004-04-21 devnull sincos(&g.wlon);
1217 28994509 2004-04-21 devnull g.nlat.l = lat*RAD;
1218 28994509 2004-04-21 devnull sincos(&g.nlat);
1219 28994509 2004-04-21 devnull if(!inwindow(&g)) {
1222 28994509 2004-04-21 devnull invert(&g);
1223 28994509 2004-04-21 devnull if(inlimits(&g)) {
1224 28994509 2004-04-21 devnull return;
1227 28994509 2004-04-21 devnull longlines = shortlines = LONGLINES;
1228 28994509 2004-04-21 devnull cut = nocut; /* not necessary; small eff. gain */
