Annotation of researchv10no/cmd/map/export/map.cpio, revision 1.1
1.1 ! root 1: 0707070035051032341006640000030000040000011211070533573712000001100000001722Makefile CFLAGS = -O
! 2: CC=cc
! 3: OBJ = map.o index.o symbol.o libmap.a
! 4:
! 5: lplot: maplplot index.o symbol.o libmap.a
! 6: $(CC) $(OBJ) -lplot -lm -o map
! 7:
! 8: v10: mapv10 index.o symbol.o libmap.a
! 9: $(CC) $(OBJ) -lm -o map
! 10:
! 11: ps: maplplot plotPS.o index.o symbol.o libmap.a
! 12: $(CC) $(OBJ) plotPS.o -lm -o map
! 13:
! 14:
! 15: mapv10: map.c iplot.h
! 16: $(CC) -c $(CFLAGS) -DPLOT='"iplot.h"' map.c
! 17:
! 18: maplplot: map.c plot.h
! 19: $(CC) -c $(CFLAGS) -DPLOT='"plot.h"' map.c
! 20:
! 21:
! 22: route: route.o libmap.a
! 23: $(CC) $(CFLAGS) route.o libmap.a -lm -o route
! 24:
! 25: libmap.a:
! 26: cd libmap; make CC=$(CC) CFLAGS="-I.. $(CFLAGS)"
! 27:
! 28: install: map
! 29: strip map
! 30: test -d /usr/lib/map || mkdir /usr/lib/map
! 31: cp map.sh /usr/bin/map
! 32: cp map mapdata/* /usr/lib/map
! 33:
! 34: clean:
! 35: rm -f map route *.o libmap.a new.results
! 36: cd libmap; make clean
! 37:
! 38: quicktest: v10
! 39: MAPPROG=./map MAPDIR=./mapdata map.sh mercator -l 0 10 0 10 -g -b >new.results
! 40: :
! 41: : If any "diff" output follows, it should show only roundoff differences.
! 42: diff test.results new.results || true
! 43: rm -f map map.o
! 44: 0707070035051032331006640000030000040000010633140533562750600000700000007126README MATERIALS
! 45:
! 46: Manuals are included in two forms.
! 47: formatted ASCII file:
! 48: map.man
! 49: input for troff -man:
! 50: map.5, map.7, proj.3, proj.5, route.1
! 51:
! 52: The source directory has two subdirectories:
! 53: libmap source for all the projection subroutines
! 54: mapdata World Data Bank I, etc.
! 55:
! 56: The source is written in ANSI C. Check CFLAGS and CC in the
! 57: Makefile and make them reflect your system's conventions.
! 58:
! 59: QUICK TEST
! 60:
! 61: For a quick test, to see whether you can compile and compute, try this
! 62: test, which makes "map" and checks its output for a bit of Africa:
! 63:
! 64: make quicktest
! 65:
! 66: PLOTTING FILTERS
! 67:
! 68: Map produces output for a plotting filter (not included).
! 69: Unfortunately there is no Unix standard for plotting.
! 70: Here are ways to compile for various filters.
! 71:
! 72: make For System V and SunOS filters plot(1)
! 73: or tplot(1).
! 74:
! 75: make v10 For v10 research system plot(1), not
! 76: compatible with System V.
! 77:
! 78: make ps PostScript. Maps are drawn in a 6.5-inch
! 79: square centered 1 inch from the top of an
! 80: 8.5x11 page. To change, edit the
! 81: PostScript or plotPS.c.
! 82:
! 83: As map uses only simple plotting features, it is usually
! 84: easy to interface to other plotting packages. See
! 85: iplot.h (used with "make v10") for details.
! 86:
! 87: When you have a plotting filter, you can test map in the current
! 88: directory as follows. For sample arguments, see EXAMPLES on the
! 89: map(7) man page.
! 90:
! 91: MAPDIR=./mapdata MAPPROG=./map map.sh arguments | filter
! 92:
! 93: INSTALLATION
! 94:
! 95: For real installation, examine the recipe for
! 96:
! 97: make install
! 98:
! 99: It puts the map shell script in /usr/bin and everything else
! 100: in /usr/lib/map. If you want to put things elsewhere, adjust
! 101: variables MAPDIR and MAPPROG in map.sh. map.sh will become
! 102: the command "map". After installation it will be run thus:
! 103:
! 104: map arguments | filter
! 105:
! 106: POSSIBLE PROBLEMS
! 107:
! 108: Options
! 109: Option -f oes not work because World Data Bank II is
! 110: not in this distribution on account of its size (13Mb).
! 111:
! 112: -y files are like v10, not Sys V, plot files.
! 113: The plot(5) man page from v10 is included for reference,
! 114: although most of it is irrelevant for this application.
! 115:
! 116: Library
! 117: At least one version of tplot(1) has been seen to garble
! 118: the output. In this case the trouble went away by
! 119: compiling map with -l4014 instead of -lplot; see plot(3).
! 120:
! 121: Collisions with a nonstandard library function sincos()
! 122: have been observed. sincos will have to be renamed
! 123: in map.h and all .c files.
! 124:
! 125: Man pages
! 126: The man pages were produced with a -man macros slightly
! 127: different from Sys V. If you change font L to B throughout
! 128: you will probably get a passable result.
! 129:
! 130: REFERENCES
! 131:
! 132: Most standard texts on cartography discuss the major projections
! 133: and their uses. Some thorough works:
! 134:
! 135: J. P. Snyder, An Album of Map Projections, US Geological Survey
! 136: Professional Paper 1453, USGPO Washington (1989)
! 137:
! 138: J. P. Snyder, Map Projections - A Working Manual, US Geological
! 139: Survey Professional Paper 1395, USGPO Washington (1987)
! 140:
! 141: D. H. Maling, Coordinate Systems and Map Projections, George
! 142: Philip and Son, London (1973)
! 143:
! 144: J. A. Steers, An Introduction to the Study of Map Projections,
! 145: Univ. London Press (1970)
! 146:
! 147: A short introduction, with 24 sample maps drawn by the "map"
! 148: program itself and the commands to draw them:
! 149:
! 150: M. D. McIlroy, Projections: Mapmakers' Answers to the Riddle of
! 151: Presenting a Round Earth on Flat Paper, AT&T Bell Labs Computing
! 152: Science Tech. Report 140 (1987), order from Computing Science
! 153: Reports, Room 2C579, AT&T Bell Labs Murray Hill, NJ 07974
! 154:
! 155: [In this report, option -s is used in an obsolete way. It
! 156: should always be replaced by -s2, with -s1 added to the
! 157: preceding command.]
! 158:
! 159: Doug McIlroy
! 160: 201-582-6050
! 161: research!doug
! 162: [email protected]
! 163: 0707070035051032721006640000030000040000020633020533472714100001000000010405index.c #include "map.h"
! 164:
! 165: /* p0;p1 evaluated to shut off "not used" diags */
! 166: static proj Yaitoff(double p0,double p1){p0;p1;return aitoff();}
! 167: static proj Yalbers(double p0,double p1){return albers(p0,p1);}
! 168: static proj Yazequalarea(double p0,double p1){p0;p1;return azequalarea();}
! 169: static proj Yazequidistant(double p0,double p1){p0;p1;return azequidistant();}
! 170: static proj Ybicentric(double p0,double p1){p1;return bicentric(p0);}
! 171: static proj Ybonne(double p0,double p1){p1;return bonne(p0);}
! 172: static proj Yconic(double p0,double p1){p1;return conic(p0);}
! 173: static proj Ycylequalarea(double p0,double p1){p1;return cylequalarea(p0);}
! 174: static proj Ycylindrical(double p0,double p1){p0;p1;return cylindrical();}
! 175: static proj Yelliptic(double p0,double p1){p1;return elliptic(p0);}
! 176: static proj Yfisheye(double p0,double p1){p1;return fisheye(p0);}
! 177: static proj Ygall(double p0,double p1){p1;return gall(p0);}
! 178: static proj Yglobular(double p0,double p1){p0;p1;return globular();}
! 179: static proj Ygnomonic(double p0,double p1){p0;p1;return gnomonic();}
! 180: static proj Yguyou(double p0,double p1){p0;p1;return guyou();}
! 181: static proj Yharrison(double p0,double p1){return harrison(p0,p1);}
! 182: static proj Yhex(double p0,double p1){p0;p1;return hex();}
! 183: static proj Yhoming(double p0,double p1){p1;return homing(p0);}
! 184: static proj Ylagrange(double p0,double p1){p0;p1;return lagrange();}
! 185: static proj Ylambert(double p0,double p1){return lambert(p0,p1);}
! 186: static proj Ylaue(double p0,double p1){p0;p1;return laue();}
! 187: static proj Ymecca(double p0,double p1){p1;return mecca(p0);}
! 188: static proj Ymercator(double p0,double p1){p0;p1;return mercator();}
! 189: static proj Ymollweide(double p0,double p1){p0;p1;return mollweide();}
! 190: static proj Ynewyorker(double p0,double p1){p1;return newyorker(p0);}
! 191: static proj Yorthographic(double p0,double p1){p0;p1;return orthographic();}
! 192: static proj Yperspective(double p0,double p1){p1;return perspective(p0);}
! 193: static proj Ypolyconic(double p0,double p1){p0;p1;return polyconic();}
! 194: static proj Yrectangular(double p0,double p1){p1;return rectangular(p0);}
! 195: static proj Ysimpleconic(double p0,double p1){return simpleconic(p0,p1);}
! 196: static proj Ysinusoidal(double p0,double p1){p0;p1;return sinusoidal();}
! 197: static proj Ysp_albers(double p0,double p1){return sp_albers(p0,p1);}
! 198: static proj Ysp_mercator(double p0,double p1){p0;p1;return sp_mercator();}
! 199: static proj Ysquare(double p0,double p1){p0;p1;return square();}
! 200: static proj Ystereographic(double p0,double p1){p0;p1;return stereographic();}
! 201: static proj Ytetra(double p0,double p1){p0;p1;return tetra();}
! 202: static proj Ytrapezoidal(double p0,double p1){return trapezoidal(p0,p1);}
! 203: static proj Yvandergrinten(double p0,double p1){p0;p1;return vandergrinten();}
! 204:
! 205: struct index index[] = {
! 206: {"aitoff", Yaitoff, 0, picut, 0, 0},
! 207: {"albers", Yalbers, 2, picut, 3, 0},
! 208: {"azequalarea", Yazequalarea, 0, nocut, 1, 0},
! 209: {"azequidistant", Yazequidistant, 0, nocut, 1, 0},
! 210: {"bicentric", Ybicentric, 1, nocut, 0, 0},
! 211: {"bonne", Ybonne, 1, picut, 0, 0},
! 212: {"conic", Yconic, 1, picut, 0, 0},
! 213: {"cylequalarea", Ycylequalarea, 1, picut, 3, 0},
! 214: {"cylindrical", Ycylindrical, 0, picut, 0, 0},
! 215: {"elliptic", Yelliptic, 1, picut, 0, 0},
! 216: {"fisheye", Yfisheye, 1, nocut, 0, 0},
! 217: {"gall", Ygall, 1, picut, 3, 0},
! 218: {"globular", Yglobular, 0, picut, 0, 0},
! 219: {"gnomonic", Ygnomonic, 0, nocut, 0, 0},
! 220: {"guyou", Yguyou, 0, guycut, 0, 0},
! 221: {"harrison", Yharrison, 2, nocut, 0, 0},
! 222: {"hex", Yhex, 0, hexcut, 0, 0},
! 223: {"homing", Yhoming, 1, picut, 0, 0},
! 224: {"lagrange", Ylagrange,0,picut,0, 0},
! 225: {"lambert", Ylambert, 2, picut, 0, 0},
! 226: {"laue", Ylaue, 0, nocut, 0, 0},
! 227: {"mecca", Ymecca, 1, picut, 0, 0},
! 228: {"mercator", Ymercator, 0, picut, 0, 0},
! 229: {"mollweide", Ymollweide, 0, picut, 0, 0},
! 230: {"newyorker", Ynewyorker, 1, nocut, 0, 0},
! 231: {"orthographic", Yorthographic, 0, nocut, 0, 0},
! 232: {"perspective", Yperspective, 1, nocut, 0, 0},
! 233: {"polyconic", Ypolyconic, 0, picut, 0, 0},
! 234: {"rectangular", Yrectangular, 1, picut, 3, 0},
! 235: {"simpleconic", Ysimpleconic, 2, picut, 3, 0},
! 236: {"sinusoidal", Ysinusoidal, 0, picut, 0, 0},
! 237: {"sp_albers", Ysp_albers, 2, picut, 3, 1},
! 238: {"sp_mercator", Ysp_mercator, 0, picut, 0, 1},
! 239: {"square", Ysquare, 0, picut, 0, 0},
! 240: {"stereographic", Ystereographic, 0, nocut, 0, 0},
! 241: {"tetra", Ytetra, 0, tetracut, 0, 0},
! 242: {"trapezoidal", Ytrapezoidal, 2, picut, 3, 0},
! 243: {"vandergrinten", Yvandergrinten, 0, picut, 0, 0},
! 244: 0
! 245: };
! 246: 0707070035051032701006640000030000040000020035130417420602000001000000002016iplot.h /* Plotting functions for v8 and v9 systems */
! 247: /* This file is an alternative to plot.h */
! 248:
! 249: /* open the plotting output */
! 250: #define openpl() printf("o\n")
! 251:
! 252: /* close the plotting output */
! 253: #define closepl() printf("cl\n")
! 254:
! 255: /* make sure the page or screen is clear */
! 256: #define erase() printf("e\n")
! 257:
! 258: /* plot a point at _x,_y, which becomes current */
! 259: #define point(_x,_y) printf("poi %d %d\n", _x,_y)
! 260:
! 261: /* coordinates to be assigned to lower left and upper right
! 262: corners of (square) plotting area */
! 263: #define range(_x,_y,_X,_Y) printf("ra %d %d %d %d\n", _x,_y,_X,_Y)
! 264:
! 265: /* place text, first letter at current point, which does not change */
! 266: #define text(_s) {if(*(_s) == ' ')printf("t \"%s\"\n",_s); else printf("t %s\n", _s); }
! 267:
! 268: /* draw line from current point to _x,_y, which becomes current */
! 269: #define vec(_x,_y) printf("v %d %d\n", _x,_y)
! 270:
! 271: /* _x,_y becomes current point */
! 272: #define move(_x, _y) printf("m %d %d\n", _x, _y)
! 273:
! 274: /* specify style for drawing lines: "dotted", "solid", "dotdash" */
! 275: #define pen(_s) printf("pe %s\n", _s)
! 276: 0707070035051031101006640000030000040000011401020525002214000002000000001206libmap/Makefile
! 277: # .o's ordered for a nonrandom library
! 278: OBJ= aitoff.o albers.o azequalarea.o elliptic.o azequidist.o \
! 279: bicentric.o bonne.o conic.o cylequalarea.o cylindrical.o fisheye.o \
! 280: gall.o guyou.o harrison.o hex.o homing.o lagrange.o lambert.o \
! 281: laue.o mercator.o mollweide.o newyorker.o polyconic.o simpleconic.o \
! 282: sinusoidal.o tetra.o perspective.o orthographic.o trapezoidal.o \
! 283: rectangular.o twocirc.o cuts.o ccubrt.o cubrt.o elco2.o complex.o zcoord.o
! 284:
! 285: # ignore error on systems without ranlib
! 286:
! 287: ../libmap.a: ../map.h $(OBJ)
! 288: ar cr ../libmap.a $(OBJ)
! 289: ranlib ../libmap.a 2>/dev/null || true
! 290:
! 291: clean:
! 292: rm -f *.o
! 293:
! 294: %.o: %.c
! 295: $(CC) $(CFLAGS) -I.. -c $*.c
! 296: 0707070035051031071006640000030000040000011231100533472766300002000000000554libmap/aitoff.c #include "map.h"
! 297:
! 298: #define Xaitwist Xaitpole.nlat
! 299: static struct place Xaitpole;
! 300:
! 301: static int
! 302: Xaitoff(struct place *place, double *x, double *y)
! 303: {
! 304: struct place p;
! 305: copyplace(place,&p);
! 306: p.wlon.l /= 2.;
! 307: sincos(&p.wlon);
! 308: norm(&p,&Xaitpole,&Xaitwist);
! 309: Xazequalarea(&p,x,y);
! 310: *x *= 2.;
! 311: return(1);
! 312: }
! 313:
! 314: proj
! 315: aitoff(void)
! 316: {
! 317: latlon(0.,0.,&Xaitpole);
! 318: return(Xaitoff);
! 319: }
! 320: 0707070035051031061006640000030000040000011400540533472766300002000000004555libmap/albers.c #include "map.h"
! 321:
! 322: /* For Albers formulas see Deetz and Adams "Elements of Map Projection", */
! 323: /* USGS Special Publication No. 68, GPO 1921 */
! 324:
! 325: static double r0sq, r1sq, d2, n, den, sinb1, sinb2;
! 326: static struct coord plat1, plat2;
! 327: static southpole;
! 328:
! 329: static double num(double s)
! 330: {
! 331: if(d2==0)
! 332: return(1);
! 333: s = d2*s*s;
! 334: return(1+s*(2./3+s*(3./5+s*(4./7+s*5./9))));
! 335: }
! 336:
! 337: /* Albers projection for a spheroid, good only when N pole is fixed */
! 338:
! 339: static int
! 340: Xspalbers(struct place *place, double *x, double *y)
! 341: {
! 342: double r = sqrt(r0sq-2*(1-d2)*place->nlat.s*num(place->nlat.s)/n);
! 343: double t = n*place->wlon.l;
! 344: *y = r*cos(t);
! 345: *x = -r*sin(t);
! 346: if(!southpole)
! 347: *y = -*y;
! 348: else
! 349: *x = -*x;
! 350: return(1);
! 351: }
! 352:
! 353: /* lat1, lat2: std parallels; e2: squared eccentricity */
! 354:
! 355: static proj albinit(double lat1, double lat2, double e2)
! 356: {
! 357: double r1,r2;
! 358: double t;
! 359: for(;;) {
! 360: if(lat1 < -90)
! 361: lat1 = -180 - lat1;
! 362: if(lat2 > 90)
! 363: lat2 = 180 - lat2;
! 364: if(lat1 <= lat2)
! 365: break;
! 366: t = lat1; lat1 = lat2; lat2 = t;
! 367: }
! 368: if(lat2-lat1 < 1) {
! 369: if(lat1 > 89)
! 370: return(azequalarea());
! 371: return(0);
! 372: }
! 373: if(fabs(lat2+lat1) < 1)
! 374: return(cylequalarea(lat1));
! 375: d2 = e2;
! 376: den = num(1.);
! 377: deg2rad(lat1,&plat1);
! 378: deg2rad(lat2,&plat2);
! 379: sinb1 = plat1.s*num(plat1.s)/den;
! 380: sinb2 = plat2.s*num(plat2.s)/den;
! 381: n = (plat1.c*plat1.c/(1-e2*plat1.s*plat1.s) -
! 382: plat2.c*plat2.c/(1-e2*plat2.s*plat2.s)) /
! 383: (2*(1-e2)*den*(sinb2-sinb1));
! 384: r1 = plat1.c/(n*sqrt(1-e2*plat1.s*plat1.s));
! 385: r2 = plat2.c/(n*sqrt(2-e2*plat2.s*plat2.s));
! 386: r1sq = r1*r1;
! 387: r0sq = r1sq + 2*(1-e2)*den*sinb1/n;
! 388: southpole = lat1<0 && plat2.c>plat1.c;
! 389: return(Xspalbers);
! 390: }
! 391:
! 392: proj
! 393: sp_albers(double lat1, double lat2)
! 394: {
! 395: return(albinit(lat1,lat2,EC2));
! 396: }
! 397:
! 398: proj
! 399: albers(double lat1, double lat2)
! 400: {
! 401: return(albinit(lat1,lat2,0.));
! 402: }
! 403:
! 404: static double scale = 1;
! 405: static double twist = 0;
! 406:
! 407: void
! 408: albscale(double x, double y, double lat, double lon)
! 409: {
! 410: struct place place;
! 411: double alat, alon, x1,y1;
! 412: scale = 1;
! 413: twist = 0;
! 414: invalb(x,y,&alat,&alon);
! 415: twist = lon - alon;
! 416: deg2rad(lat,&place.nlat);
! 417: deg2rad(lon,&place.wlon);
! 418: Xspalbers(&place,&x1,&y1);
! 419: scale = sqrt((x1*x1+y1*y1)/(x*x+y*y));
! 420: }
! 421:
! 422: void
! 423: invalb(double x, double y, double *lat, double *lon)
! 424: {
! 425: int i;
! 426: double sinb_den, sinp;
! 427: x *= scale;
! 428: y *= scale;
! 429: *lon = atan2(-x,fabs(y))/(RAD*n) + twist;
! 430: sinb_den = (r0sq - x*x - y*y)*n/(2*(1-d2));
! 431: sinp = sinb_den;
! 432: for(i=0; i<5; i++)
! 433: sinp = sinb_den/num(sinp);
! 434: *lat = asin(sinp)/RAD;
! 435: }
! 436: 0707070035051031051006640000030000040000011400560533472766300002500000000361libmap/azequalarea.c #include "map.h"
! 437:
! 438: int
! 439: Xazequalarea(struct place *place, double *x, double *y)
! 440: {
! 441: double r;
! 442: r = sqrt(1. - place->nlat.s);
! 443: *x = - r * place->wlon.s;
! 444: *y = - r * place->wlon.c;
! 445: return(1);
! 446: }
! 447:
! 448: proj
! 449: azequalarea(void)
! 450: {
! 451: return(Xazequalarea);
! 452: }
! 453: 0707070035051031041006640000030000040000011400600533472766400002400000000401libmap/azequidist.c #include "map.h"
! 454:
! 455: int
! 456: Xazequidistant(struct place *place, double *x, double *y)
! 457: {
! 458: double colat;
! 459: colat = PI/2 - place->nlat.l;
! 460: *x = -colat * place->wlon.s;
! 461: *y = -colat * place->wlon.c;
! 462: return(1);
! 463: }
! 464:
! 465: proj
! 466: azequidistant(void)
! 467: {
! 468: return(Xazequidistant);
! 469: }
! 470: 0707070035051031031006640000030000040000011400620533472766400002300000000624libmap/bicentric.c #include "map.h"
! 471:
! 472: static struct coord center;
! 473:
! 474: static int
! 475: Xbicentric(struct place *place, double *x, double *y)
! 476: {
! 477: if(place->wlon.c<=.01||place->nlat.c<=.01)
! 478: return(-1);
! 479: *x = -center.c*place->wlon.s/place->wlon.c;
! 480: *y = place->nlat.s/(place->nlat.c*place->wlon.c);
! 481: return(*x**x+*y**y<=9);
! 482: }
! 483:
! 484: proj
! 485: bicentric(double l)
! 486: {
! 487: l = fabs(l);
! 488: if(l>89)
! 489: return(0);
! 490: deg2rad(l,¢er);
! 491: return(Xbicentric);
! 492: }
! 493: 0707070035051031021006640000030000040000011400640533472766400001700000001164libmap/bonne.c #include "map.h"
! 494:
! 495: static struct coord stdpar;
! 496: static double r0;
! 497:
! 498: static
! 499: Xbonne(struct place *place, double *x, double *y)
! 500: {
! 501: double r, alpha;
! 502: r = r0 - place->nlat.l;
! 503: if(r<.001)
! 504: if(fabs(stdpar.c)<1e-10)
! 505: alpha = place->wlon.l;
! 506: else if(fabs(place->nlat.c)==0)
! 507: alpha = 0;
! 508: else
! 509: alpha = place->wlon.l/(1+
! 510: stdpar.c*stdpar.c*stdpar.c/place->nlat.c/3);
! 511: else
! 512: alpha = place->wlon.l * place->nlat.c / r;
! 513: *x = - r*sin(alpha);
! 514: *y = - r*cos(alpha);
! 515: return(1);
! 516: }
! 517:
! 518: proj
! 519: bonne(double par)
! 520: {
! 521: if(fabs(par*RAD) < .01)
! 522: return(Xsinusoidal);
! 523: deg2rad(par, &stdpar);
! 524: r0 = stdpar.c/stdpar.s + stdpar.l;
! 525: return(Xbonne);
! 526: }
! 527: 0707070035051031011006640000030000040000011400660470474632500002000000000301libmap/ccubrt.c #include "map.h"
! 528:
! 529: void
! 530: ccubrt(double zr, double zi, double *wr, double *wi)
! 531: {
! 532: double r, theta;
! 533: theta = atan2(zi,zr);
! 534: r = cubrt(hypot(zr,zi));
! 535: *wr = r*cos(theta/3);
! 536: *wi = r*sin(theta/3);
! 537: }
! 538: 0707070035051031001006640000030000040000011400700470474632500002100000002535libmap/complex.c #include "map.h"
! 539:
! 540: /*complex divide, defensive against overflow from
! 541: * * and /, but not from + and -
! 542: * assumes underflow yields 0.0
! 543: * uses identities:
! 544: * (a + bi)/(c + di) = ((a + bd/c) + (b - ad/c)i)/(c + dd/c)
! 545: * (a + bi)/(c + di) = (b - ai)/(d - ci)
! 546: */
! 547: void
! 548: cdiv(double a, double b, double c, double d, double *u, double *v)
! 549: {
! 550: double r,t;
! 551: if(fabs(c)<fabs(d)) {
! 552: t = -c; c = d; d = t;
! 553: t = -a; a = b; b = t;
! 554: }
! 555: r = d/c;
! 556: t = c + r*d;
! 557: *u = (a + r*b)/t;
! 558: *v = (b - r*a)/t;
! 559: }
! 560:
! 561: void
! 562: cmul(double c1, double c2, double d1, double d2, double *e1, double *e2)
! 563: {
! 564: *e1 = c1*d1 - c2*d2;
! 565: *e2 = c1*d2 + c2*d1;
! 566: }
! 567:
! 568: void
! 569: csq(double c1, double c2, double *e1, double *e2)
! 570: {
! 571: *e1 = c1*c1 - c2*c2;
! 572: *e2 = c1*c2*2;
! 573: }
! 574:
! 575: /* complex square root
! 576: * assumes underflow yields 0.0
! 577: * uses these identities:
! 578: * sqrt(x+_iy) = sqrt(r(cos(t)+_isin(t))
! 579: * = sqrt(r)(cos(t/2)+_isin(t/2))
! 580: * cos(t/2) = sin(t)/2sin(t/2) = sqrt((1+cos(t)/2)
! 581: * sin(t/2) = sin(t)/2cos(t/2) = sqrt((1-cos(t)/2)
! 582: */
! 583: void
! 584: csqrt(double c1, double c2, double *e1, double *e2)
! 585: {
! 586: double r,s;
! 587: double x,y;
! 588: x = fabs(c1);
! 589: y = fabs(c2);
! 590: if(x>=y) {
! 591: if(x==0) {
! 592: *e1 = *e2 = 0;
! 593: return;
! 594: }
! 595: r = x;
! 596: s = y/x;
! 597: } else {
! 598: r = y;
! 599: s = x/y;
! 600: }
! 601: r *= sqrt(1+ s*s);
! 602: if(c1>0) {
! 603: *e1 = sqrt((r+c1)/2);
! 604: *e2 = c2/(2* *e1);
! 605: } else {
! 606: *e2 = sqrt((r-c1)/2);
! 607: if(c2<0)
! 608: *e2 = -*e2;
! 609: *e1 = c2/(2* *e2);
! 610: }
! 611: }
! 612: 0707070035051030771006640000030000040000011400720533472766500001700000000717libmap/conic.c #include "map.h"
! 613:
! 614: static struct coord stdpar;
! 615:
! 616: static int
! 617: Xconic(struct place *place, double *x, double *y)
! 618: {
! 619: double r;
! 620: if(fabs(place->nlat.l-stdpar.l) > 80.*RAD)
! 621: return(-1);
! 622: r = stdpar.c/stdpar.s - tan(place->nlat.l - stdpar.l);
! 623: *x = - r*sin(place->wlon.l * stdpar.s);
! 624: *y = - r*cos(place->wlon.l * stdpar.s);
! 625: if(r>3) return(0);
! 626: return(1);
! 627: }
! 628:
! 629: proj
! 630: conic(double par)
! 631: {
! 632: if(fabs(par) <.1)
! 633: return(Xcylindrical);
! 634: deg2rad(par, &stdpar);
! 635: return(Xconic);
! 636: }
! 637: 0707070035051030761006640000030000040000011400740470474632600001700000000450libmap/cubrt.c #include "map.h"
! 638:
! 639: double
! 640: cubrt(double a)
! 641: {
! 642: double x,y,x1;
! 643: if(a==0)
! 644: return(0.);
! 645: y = 1;
! 646: if(a<0) {
! 647: y = -y;
! 648: a = -a;
! 649: }
! 650: while(a<1) {
! 651: a *= 8;
! 652: y /= 2;
! 653: }
! 654: while(a>1) {
! 655: a /= 8;
! 656: y *= 2;
! 657: }
! 658: x = 1;
! 659: do {
! 660: x1 = x;
! 661: x = (2*x1+a/(x1*x1))/3;
! 662: } while(fabs(x-x1)>10.e-15);
! 663: return(x*y);
! 664: }
! 665: 0707070035051030451006640000030000040000011400760533472766500001600000001453libmap/cuts.c #include "map.h"
! 666: extern void abort(void);
! 667:
! 668: /* these routines duplicate names found in map.c. they are
! 669: called from routines in hex.c, guyou.c, and tetra.c, which
! 670: are in turn invoked directly from map.c. this bad organization
! 671: arises from data hiding; only these three files know stuff
! 672: that's necessary for the proper handling of the unusual cuts
! 673: involved in these projections.
! 674:
! 675: the calling routines are not advertised as part of the library,
! 676: and the library duplicates should never get loaded, however they
! 677: are included to make the libary self-standing.*/
! 678:
! 679: int
! 680: picut(struct place *g, struct place *og, double *cutlon)
! 681: {
! 682: g; og; cutlon;
! 683: abort();
! 684: return 0;
! 685: }
! 686:
! 687: int
! 688: ckcut(struct place *g1, struct place *g2, double lon)
! 689: {
! 690: g1; g2; lon;
! 691: abort();
! 692: return 0;
! 693: }
! 694:
! 695: double
! 696: reduce(double x)
! 697: {
! 698: x;
! 699: abort();
! 700: return 0;
! 701: }
! 702: 0707070035051030751006640000030000040000011401000533472766500002600000000502libmap/cylequalarea.c #include "map.h"
! 703:
! 704: static double a;
! 705:
! 706: static int
! 707: Xcylequalarea(struct place *place, double *x, double *y)
! 708: {
! 709: *x = - place->wlon.l * a;
! 710: *y = place->nlat.s;
! 711: return(1);
! 712: }
! 713:
! 714: proj
! 715: cylequalarea(double par)
! 716: {
! 717: struct coord stdp0;
! 718: if(par > 89.0)
! 719: return(0);
! 720: deg2rad(par, &stdp0);
! 721: a = stdp0.c*stdp0.c;
! 722: return(Xcylequalarea);
! 723: }
! 724: 0707070035051030741006640000030000040000011401100533472766600002500000000376libmap/cylindrical.c #include "map.h"
! 725:
! 726: int
! 727: Xcylindrical(struct place *place, double *x, double *y)
! 728: {
! 729: if(fabs(place->nlat.l) > 80.*RAD)
! 730: return(-1);
! 731: *x = - place->wlon.l;
! 732: *y = place->nlat.s / place->nlat.c;
! 733: return(1);
! 734: }
! 735:
! 736: proj
! 737: cylindrical(void)
! 738: {
! 739: return(Xcylindrical);
! 740: }
! 741: 0707070035051030731006640000030000040000011401120533472766600001700000004704libmap/elco2.c #include "map.h"
! 742:
! 743: /* elliptic integral routine, R.Bulirsch,
! 744: * Numerische Mathematik 7(1965) 78-90
! 745: * calculate integral from 0 to x+iy of
! 746: * (a+b*t^2)/((1+t^2)*sqrt((1+t^2)*(1+kc^2*t^2)))
! 747: * yields about D valid figures, where CC=10e-D
! 748: * for a*b>=0, except at branchpoints x=0,y=+-i,+-i/kc;
! 749: * there the accuracy may be reduced.
! 750: * fails for kc=0 or x<0
! 751: * return(1) for success, return(0) for fail
! 752: *
! 753: * special case a=b=1 is equivalent to
! 754: * standard elliptic integral of first kind
! 755: * from 0 to atan(x+iy) of
! 756: * 1/sqrt(1-k^2*(sin(t))^2) where k^2=1-kc^2
! 757: */
! 758:
! 759: #define ROOTINF 10.e18
! 760: #define CC 1.e-6
! 761:
! 762: int
! 763: elco2(double x, double y, double kc, double a, double b, double *u, double *v)
! 764: {
! 765: double c,d,dn1,dn2,e,e1,e2,f,f1,f2,h,k,m,m1,m2,sy;
! 766: double d1[13],d2[13];
! 767: int i,l;
! 768: if(kc==0||x<0)
! 769: return(0);
! 770: sy = y>0? 1: y==0? 0: -1;
! 771: y = fabs(y);
! 772: csq(x,y,&c,&e2);
! 773: d = kc*kc;
! 774: k = 1-d;
! 775: e1 = 1+c;
! 776: cdiv2(1+d*c,d*e2,e1,e2,&f1,&f2);
! 777: f2 = -k*x*y*2/f2;
! 778: csqr(f1,f2,&dn1,&dn2);
! 779: if(f1<0) {
! 780: f1 = dn1;
! 781: dn1 = -dn2;
! 782: dn2 = -f1;
! 783: }
! 784: if(k<0) {
! 785: dn1 = fabs(dn1);
! 786: dn2 = fabs(dn2);
! 787: }
! 788: c = 1+dn1;
! 789: cmul(e1,e2,c,dn2,&f1,&f2);
! 790: cdiv(x,y,f1,f2,&d1[0],&d2[0]);
! 791: h = a-b;
! 792: d = f = m = 1;
! 793: kc = fabs(kc);
! 794: e = a;
! 795: a += b;
! 796: l = 4;
! 797: for(i=1;;i++) {
! 798: m1 = (kc+m)/2;
! 799: m2 = m1*m1;
! 800: k *= f/(m2*4);
! 801: b += e*kc;
! 802: e = a;
! 803: cdiv2(kc+m*dn1,m*dn2,c,dn2,&f1,&f2);
! 804: csqr(f1/m1,k*dn2*2/f2,&dn1,&dn2);
! 805: cmul(dn1,dn2,x,y,&f1,&f2);
! 806: x = fabs(f1);
! 807: y = fabs(f2);
! 808: a += b/m1;
! 809: l *= 2;
! 810: c = 1 +dn1;
! 811: d *= k/2;
! 812: cmul(x,y,x,y,&e1,&e2);
! 813: k *= k;
! 814:
! 815: cmul(c,dn2,1+e1*m2,e2*m2,&f1,&f2);
! 816: cdiv(d*x,d*y,f1,f2,&d1[i],&d2[i]);
! 817: if(k<=CC)
! 818: break;
! 819: kc = sqrt(m*kc);
! 820: f = m2;
! 821: m = m1;
! 822: }
! 823: f1 = f2 = 0;
! 824: for(;i>=0;i--) {
! 825: f1 += d1[i];
! 826: f2 += d2[i];
! 827: }
! 828: x *= m1;
! 829: y *= m1;
! 830: cdiv2(1-y,x,1+y,-x,&e1,&e2);
! 831: e2 = x*2/e2;
! 832: d = a/(m1*l);
! 833: *u = atan2(e2,e1);
! 834: if(*u<0)
! 835: *u += PI;
! 836: a = d*sy/2;
! 837: *u = d*(*u) + f1*h;
! 838: *v = (-1-log(e1*e1+e2*e2))*a + f2*h*sy + a;
! 839: return(1);
! 840: }
! 841:
! 842: void
! 843: cdiv2(double c1, double c2, double d1, double d2, double *e1, double *e2)
! 844: {
! 845: double t;
! 846: if(fabs(d2)>fabs(d1)) {
! 847: t = d1, d1 = d2, d2 = t;
! 848: t = c1, c1 = c2, c2 = t;
! 849: }
! 850: if(fabs(d1)>ROOTINF)
! 851: *e2 = ROOTINF*ROOTINF;
! 852: else
! 853: *e2 = d1*d1 + d2*d2;
! 854: t = d2/d1;
! 855: *e1 = (c1+t*c2)/(d1+t*d2); /* (c1*d1+c2*d2)/(d1*d1+d2*d2) */
! 856: }
! 857:
! 858: /* complex square root of |x|+iy */
! 859: void
! 860: csqr(double c1, double c2, double *e1, double *e2)
! 861: {
! 862: double r2;
! 863: r2 = c1*c1 + c2*c2;
! 864: if(r2<=0) {
! 865: *e1 = *e2 = 0;
! 866: return;
! 867: }
! 868: *e1 = sqrt((sqrt(r2) + fabs(c1))/2);
! 869: *e2 = c2/(*e1*2);
! 870: }
! 871: 0707070035051030721006640000030000040000010612660533472766600002200000001122libmap/elliptic.c #include "map.h"
! 872:
! 873: struct coord center;
! 874:
! 875: static int
! 876: Xelliptic(struct place *place, double *x, double *y)
! 877: {
! 878: double r1,r2;
! 879: r1 = acos(place->nlat.c*(place->wlon.c*center.c
! 880: - place->wlon.s*center.s));
! 881: r2 = acos(place->nlat.c*(place->wlon.c*center.c
! 882: + place->wlon.s*center.s));
! 883: *x = -(r1*r1 - r2*r2)/(4*center.l);
! 884: *y = (r1*r1+r2*r2)/2 - (center.l*center.l+*x**x);
! 885: if(*y < 0)
! 886: *y = 0;
! 887: *y = sqrt(*y);
! 888: if(place->nlat.l<0)
! 889: *y = -*y;
! 890: return(1);
! 891: }
! 892:
! 893: proj
! 894: elliptic(double l)
! 895: {
! 896: l = fabs(l);
! 897: if(l>89)
! 898: return(0);
! 899: if(l<1)
! 900: return(Xazequidistant);
! 901: deg2rad(l,¢er);
! 902: return(Xelliptic);
! 903: }
! 904: 0707070035050603501006640000100000040000011210700533472766600002100000000566libmap/fisheye.c #include "map.h"
! 905: /* refractive fisheye, not logarithmic */
! 906:
! 907: static double n;
! 908:
! 909: static int
! 910: Xfisheye(struct place *place, double *x, double *y)
! 911: {
! 912: double r;
! 913: double u = sin(PI/4-place->nlat.l/2)/n;
! 914: if(fabs(u) > .97)
! 915: return -1;
! 916: r = tan(asin(u));
! 917: *x = -r*place->wlon.s;
! 918: *y = -r*place->wlon.c;
! 919: return 1;
! 920: }
! 921:
! 922: proj
! 923: fisheye(double par)
! 924: {
! 925: n = par;
! 926: return n<.1? 0: Xfisheye;
! 927: }
! 928: 0707070035051030401006640000030000040000011401140533472766700001600000000737libmap/gall.c #include "map.h"
! 929:
! 930: static double scale;
! 931:
! 932: static int
! 933: Xgall(struct place *place, double *x, double *y)
! 934: {
! 935: /* two ways to compute tan(place->nlat.l/2) */
! 936: if(fabs(place->nlat.s)<.1)
! 937: *y = sin(place->nlat.l/2)/cos(place->nlat.l/2);
! 938: else
! 939: *y = (1-place->nlat.c)/place->nlat.s;
! 940: *x = -scale*place->wlon.l;
! 941: return 1;
! 942: }
! 943:
! 944: proj
! 945: gall(double par)
! 946: {
! 947: double coshalf;
! 948: if(fabs(par)>80)
! 949: return 0;
! 950: par *= RAD;
! 951: coshalf = cos(par/2);
! 952: scale = cos(par)/(2*coshalf*coshalf);
! 953: return Xgall;
! 954: }
! 955: 0707070035051030701006640000030000040000011402000533472766700001700000003271libmap/guyou.c #include "map.h"
! 956:
! 957: static struct place gywhem, gyehem;
! 958: static struct coord gytwist;
! 959: static double gyconst, gykc, gyside;
! 960:
! 961:
! 962: static void
! 963: dosquare(double z1, double z2, double *x, double *y)
! 964: {
! 965: double w1,w2;
! 966: w1 = z1 -1;
! 967: if(fabs(w1*w1+z2*z2)>.000001) {
! 968: cdiv(z1+1,z2,w1,z2,&w1,&w2);
! 969: w1 *= gyconst;
! 970: w2 *= gyconst;
! 971: if(w1<0)
! 972: w1 = 0;
! 973: elco2(w1,w2,gykc,1.,1.,x,y);
! 974: } else {
! 975: *x = gyside;
! 976: *y = 0;
! 977: }
! 978: }
! 979:
! 980: int
! 981: Xguyou(struct place *place, double *x, double *y)
! 982: {
! 983: int ew; /*which hemisphere*/
! 984: double z1,z2;
! 985: struct place pl;
! 986: ew = place->wlon.l<0;
! 987: copyplace(place,&pl);
! 988: norm(&pl,ew?&gyehem:&gywhem,&gytwist);
! 989: Xstereographic(&pl,&z1,&z2);
! 990: dosquare(z1/2,z2/2,x,y);
! 991: if(!ew)
! 992: *x -= gyside;
! 993: return(1);
! 994: }
! 995:
! 996: proj
! 997: guyou(void)
! 998: {
! 999: double junk;
! 1000: gykc = 1/(3+2*sqrt(2.));
! 1001: gyconst = -(1+sqrt(2.));
! 1002: elco2(-gyconst,0.,gykc,1.,1.,&gyside,&junk);
! 1003: gyside *= 2;
! 1004: latlon(0.,90.,&gywhem);
! 1005: latlon(0.,-90.,&gyehem);
! 1006: deg2rad(0.,&gytwist);
! 1007: return(Xguyou);
! 1008: }
! 1009:
! 1010: int
! 1011: guycut(struct place *g, struct place *og, double *cutlon)
! 1012: {
! 1013: int c;
! 1014: c = picut(g,og,cutlon);
! 1015: if(c!=1)
! 1016: return(c);
! 1017: *cutlon = 0.;
! 1018: if(g->nlat.c<.7071||og->nlat.c<.7071)
! 1019: return(ckcut(g,og,0.));
! 1020: return(1);
! 1021: }
! 1022:
! 1023: static int
! 1024: Xsquare(struct place *place, double *x, double *y)
! 1025: {
! 1026: double z1,z2;
! 1027: double r, theta;
! 1028: struct place p;
! 1029: copyplace(place,&p);
! 1030: if(place->nlat.l<0) {
! 1031: p.nlat.l = -p.nlat.l;
! 1032: p.nlat.s = -p.nlat.s;
! 1033: }
! 1034: if(p.nlat.l<FUZZ && fabs(p.wlon.l)>PI-FUZZ){
! 1035: *y = -gyside/2;
! 1036: *x = p.wlon.l>0?0:gyside;
! 1037: return(1);
! 1038: }
! 1039: Xstereographic(&p,&z1,&z2);
! 1040: r = sqrt(sqrt(hypot(z1,z2)/2));
! 1041: theta = atan2(z1,-z2)/4;
! 1042: dosquare(r*sin(theta),-r*cos(theta),x,y);
! 1043: if(place->nlat.l<0)
! 1044: *y = -gyside - *y;
! 1045: return(1);
! 1046: }
! 1047:
! 1048: proj
! 1049: square(void)
! 1050: {
! 1051: guyou();
! 1052: return(Xsquare);
! 1053: }
! 1054: 0707070035051030501006640000030000040000010631010533472766700002200000001241libmap/harrison.c #include "map.h"
! 1055:
! 1056: static double v3,u2,u3,a,b; /*v=view,p=obj,u=unit.y*/
! 1057:
! 1058: static int
! 1059: Xharrison(struct place *place, double *x, double *y)
! 1060: {
! 1061: double p1 = -place->nlat.c*place->wlon.s;
! 1062: double p2 = -place->nlat.c*place->wlon.c;
! 1063: double p3 = place->nlat.s;
! 1064: double d = b + u3*p2 - u2*p3;
! 1065: double t;
! 1066: if(d < .01)
! 1067: return 0;
! 1068: t = a/d;
! 1069: if(t < 0)
! 1070: return 0;
! 1071: if(v3*place->nlat.s < 1.)
! 1072: return -1;
! 1073: *y = t*p2*u2 + (v3-t*(v3-p3))*u3;
! 1074: *x = t*p1;
! 1075: if(*x * *x + *y * *y > 16)
! 1076: return 0;
! 1077: return 1;
! 1078: }
! 1079:
! 1080: proj
! 1081: harrison(double r, double alpha)
! 1082: {
! 1083: u2 = cos(alpha*RAD);
! 1084: u3 = sin(alpha*RAD);
! 1085: v3 = r;
! 1086: b = r*u2;
! 1087: a = 1 + b;
! 1088: if(r<1.001 || a<sqrt(r*r-1))
! 1089: return 0;
! 1090: return Xharrison;
! 1091: }
! 1092: 0707070035051030671006640000030000040000011402020533472766700001500000004312libmap/hex.c #include "map.h"
! 1093:
! 1094: #define BIG 1.e15
! 1095: #define HFUZZ .0001
! 1096:
! 1097: static double hcut[3] ;
! 1098: static double kr[3] = { .5, -1., .5 };
! 1099: static double ki[3] = { -1., 0., 1. }; /*to multiply by sqrt(3)/2*/
! 1100: static double cr[3];
! 1101: static double ci[3];
! 1102: static struct place hem;
! 1103: static struct coord twist;
! 1104: static double rootroot3, hkc;
! 1105: static double w2;
! 1106: static double rootk;
! 1107:
! 1108: static void
! 1109: reflect(int i, double wr, double wi, double *x, double *y)
! 1110: {
! 1111: double pr,pi,l;
! 1112: pr = cr[i]-wr;
! 1113: pi = ci[i]-wi;
! 1114: l = 2*(kr[i]*pr + ki[i]*pi);
! 1115: *x = wr + l*kr[i];
! 1116: *y = wi + l*ki[i];
! 1117: }
! 1118:
! 1119: static int
! 1120: Xhex(struct place *place, double *x, double *y)
! 1121: {
! 1122: int ns;
! 1123: register i;
! 1124: double zr,zi;
! 1125: double sr,si,tr,ti,ur,ui,vr,vi,yr,yi;
! 1126: struct place p;
! 1127: copyplace(place,&p);
! 1128: ns = place->nlat.l >= 0;
! 1129: if(!ns) {
! 1130: p.nlat.l = -p.nlat.l;
! 1131: p.nlat.s = -p.nlat.s;
! 1132: }
! 1133: if(p.nlat.l<HFUZZ) {
! 1134: for(i=0;i<3;i++)
! 1135: if(fabs(reduce(p.wlon.l-hcut[i]))<HFUZZ) {
! 1136: if(i==2) {
! 1137: *x = 2*cr[0] - cr[1];
! 1138: *y = 0;
! 1139: } else {
! 1140: *x = cr[1];
! 1141: *y = 2*ci[2*i];
! 1142: }
! 1143: return(1);
! 1144: }
! 1145: p.nlat.l = HFUZZ;
! 1146: sincos(&p.nlat);
! 1147: }
! 1148: norm(&p,&hem,&twist);
! 1149: Xstereographic(&p,&zr,&zi);
! 1150: zr /= 2;
! 1151: zi /= 2;
! 1152: cdiv(1-zr,-zi,1+zr,zi,&sr,&si);
! 1153: csq(sr,si,&tr,&ti);
! 1154: ccubrt(1+3*tr,3*ti,&ur,&ui);
! 1155: csqrt(ur-1,ui,&vr,&vi);
! 1156: cdiv(rootroot3+vr,vi,rootroot3-vr,-vi,&yr,&yi);
! 1157: yr /= rootk;
! 1158: yi /= rootk;
! 1159: elco2(fabs(yr),yi,hkc,1.,1.,x,y);
! 1160: if(yr < 0)
! 1161: *x = w2 - *x;
! 1162: if(!ns) reflect(hcut[0]>place->wlon.l?0:
! 1163: hcut[1]>=place->wlon.l?1:
! 1164: 2,*x,*y,x,y);
! 1165: return(1);
! 1166: }
! 1167:
! 1168: proj
! 1169: hex(void)
! 1170: {
! 1171: register i;
! 1172: double t;
! 1173: double root3;
! 1174: double c,d;
! 1175: struct place p;
! 1176: hcut[2] = PI;
! 1177: hcut[1] = hcut[2]/3;
! 1178: hcut[0] = -hcut[1];
! 1179: root3 = sqrt(3.);
! 1180: rootroot3 = sqrt(root3);
! 1181: t = 15 -8*root3;
! 1182: hkc = t*(1-sqrt(1-1/(t*t)));
! 1183: elco2(BIG,0.,hkc,1.,1.,&w2,&t);
! 1184: w2 *= 2;
! 1185: rootk = sqrt(hkc);
! 1186: latlon(90.,90.,&hem);
! 1187: latlon(90.,0.,&p);
! 1188: Xhex(&p,&c,&t);
! 1189: latlon(0.,0.,&p);
! 1190: Xhex(&p,&d,&t);
! 1191: for(i=0;i<3;i++) {
! 1192: ki[i] *= root3/2;
! 1193: cr[i] = c + (c-d)*kr[i];
! 1194: ci[i] = (c-d)*ki[i];
! 1195: }
! 1196: deg2rad(0.,&twist);
! 1197: return(Xhex);
! 1198: }
! 1199:
! 1200: int
! 1201: hexcut(struct place *g, struct place *og, double *cutlon)
! 1202: {
! 1203: int t,i;
! 1204: if(g->nlat.l>=-HFUZZ&&og->nlat.l>=-HFUZZ)
! 1205: return(1);
! 1206: for(i=0;i<3;i++) {
! 1207: t = ckcut(g,og,*cutlon=hcut[i]);
! 1208: if(t!=1) return(t);
! 1209: }
! 1210: return(1);
! 1211: }
! 1212: 0707070035051030511006640000030000040000011402040533472767300002000000002430libmap/homing.c #include "map.h"
! 1213:
! 1214: static struct coord p0; /* standard parallel */
! 1215:
! 1216: static double
! 1217: trigclamp(double x)
! 1218: {
! 1219: return x>1? 1: x<-1? -1: x;
! 1220: }
! 1221:
! 1222: static struct coord az; /* azimuth of p0 seen from place */
! 1223: static struct coord rad; /* angular dist from place to p0 */
! 1224:
! 1225: static int
! 1226: azimuth(struct place *place)
! 1227: {
! 1228: if(place->nlat.c < .01)
! 1229: return 0;
! 1230: rad.c = trigclamp(p0.s*place->nlat.s + /* law of cosines */
! 1231: p0.c*place->nlat.c*place->wlon.c);
! 1232: rad.s = sqrt(1 - rad.c*rad.c);
! 1233: if(fabs(rad.s) < .001) {
! 1234: az.s = 0;
! 1235: az.c = 1;
! 1236: } else {
! 1237: az.s = trigclamp(p0.c*place->wlon.s/rad.s); /* sines */
! 1238: az.c = trigclamp((p0.s - rad.c*place->nlat.s)
! 1239: /(rad.s*place->nlat.c));
! 1240: }
! 1241: rad.l = atan2(rad.s, rad.c);
! 1242: return 1;
! 1243: }
! 1244:
! 1245: static int
! 1246: Xmecca(struct place *place, double *x, double *y)
! 1247: {
! 1248: if(!azimuth(place))
! 1249: return 0;
! 1250: *x = -place->wlon.l;
! 1251: *y = fabs(az.s)<.02? -az.c*rad.s/p0.c: *x*az.c/az.s;
! 1252: return fabs(*y)>2? 0:
! 1253: rad.c<0? -1:
! 1254: 1;
! 1255: }
! 1256:
! 1257: proj
! 1258: mecca(double par)
! 1259: {
! 1260: if(fabs(par)>80.)
! 1261: return(0);
! 1262: deg2rad(par,&p0);
! 1263: return(Xmecca);
! 1264: }
! 1265:
! 1266: static int
! 1267: Xhoming(struct place *place, double *x, double *y)
! 1268: {
! 1269: if(!azimuth(place))
! 1270: return 0;
! 1271: *x = -rad.l*az.s;
! 1272: *y = -rad.l*az.c;
! 1273: return place->wlon.c<0? -1: 1;
! 1274: }
! 1275:
! 1276: proj
! 1277: homing(double par)
! 1278: {
! 1279: if(fabs(par)>80.)
! 1280: return(0);
! 1281: deg2rad(par,&p0);
! 1282: return(Xhoming);
! 1283: }
! 1284:
! 1285: 0707070035051030351006640000030000040000011402060533472767300002200000000663libmap/lagrange.c #include "map.h"
! 1286:
! 1287: static int
! 1288: Xlagrange(struct place *place, double *x, double *y)
! 1289: {
! 1290: double z1,z2;
! 1291: double w1,w2,t1,t2;
! 1292: struct place p;
! 1293: copyplace(place,&p);
! 1294: if(place->nlat.l<0) {
! 1295: p.nlat.l = -p.nlat.l;
! 1296: p.nlat.s = -p.nlat.s;
! 1297: }
! 1298: Xstereographic(&p,&z1,&z2);
! 1299: csqrt(-z2/2,z1/2,&w1,&w2);
! 1300: cdiv(w1-1,w2,w1+1,w2,&t1,&t2);
! 1301: *y = -t1;
! 1302: *x = t2;
! 1303: if(place->nlat.l<0)
! 1304: *y = -*y;
! 1305: return(1);
! 1306: }
! 1307:
! 1308: proj
! 1309: lagrange(void)
! 1310: {
! 1311: return(Xlagrange);
! 1312: }
! 1313: 0707070035051030661006640000030000040000011402100533472767300002100000001576libmap/lambert.c #include "map.h"
! 1314:
! 1315: static struct coord stdp0, stdp1;
! 1316: static double k;
! 1317:
! 1318: static int
! 1319: Xlambert(struct place *place, double *x, double *y)
! 1320: {
! 1321: double r;
! 1322: if(place->nlat.l < -80.*RAD)
! 1323: return(-1);
! 1324: if(place->nlat.l > 89.*RAD)
! 1325: r = 0; /* slovenly */
! 1326: else
! 1327: r = stdp0.c*exp(0.5*k*log(
! 1328: (1+stdp0.s)*(1-place->nlat.s)/((1-stdp0.s)*(1+place->nlat.s))));
! 1329: if(stdp1.l<0.)
! 1330: r = -r;
! 1331: *x = - r*sin(k * place->wlon.l);
! 1332: *y = - r*cos(k * place->wlon.l);
! 1333: return(1);
! 1334: }
! 1335:
! 1336: proj
! 1337: lambert(double par0, double par1)
! 1338: {
! 1339: double temp;
! 1340: if(fabs(par0)>fabs(par1)){
! 1341: temp = par0;
! 1342: par0 = par1;
! 1343: par1 = temp;
! 1344: }
! 1345: deg2rad(par0, &stdp0);
! 1346: deg2rad(par1, &stdp1);
! 1347: if(fabs(par1+par0)<.1)
! 1348: return(mercator());
! 1349: if(fabs(par1-par0)<.1)
! 1350: return(perspective(-1.));
! 1351: if(fabs(par0)>89.5||fabs(par1)>89.5)
! 1352: return(0);
! 1353: k = 2*log(stdp1.c/stdp0.c)/log(
! 1354: (1+stdp0.s)*(1-stdp1.s)/((1-stdp0.s)*(1+stdp1.s)));
! 1355: return(Xlambert);
! 1356: }
! 1357: 0707070035051030651006640000030000040000011402120533472767500001600000000445libmap/laue.c #include "map.h"
! 1358:
! 1359:
! 1360: static int
! 1361: Xlaue(struct place *place, double *x, double *y)
! 1362: {
! 1363: double r;
! 1364: if(place->nlat.l<PI/4+FUZZ)
! 1365: return(-1);
! 1366: r = tan(PI-2*place->nlat.l);
! 1367: if(r>3)
! 1368: return(-1);
! 1369: *x = - r * place->wlon.s;
! 1370: *y = - r * place->wlon.c;
! 1371: return(1);
! 1372: }
! 1373:
! 1374: proj
! 1375: laue(void)
! 1376: {
! 1377: return(Xlaue);
! 1378: }
! 1379: 0707070035051030631006640000030000040000011402140533472767500002200000001026libmap/mercator.c #include "map.h"
! 1380:
! 1381: static int
! 1382: Xmercator(struct place *place, double *x, double *y)
! 1383: {
! 1384: if(fabs(place->nlat.l) > 80.*RAD)
! 1385: return(-1);
! 1386: *x = -place->wlon.l;
! 1387: *y = 0.5*log((1+place->nlat.s)/(1-place->nlat.s));
! 1388: return(1);
! 1389: }
! 1390:
! 1391: proj
! 1392: mercator(void)
! 1393: {
! 1394: return(Xmercator);
! 1395: }
! 1396:
! 1397: static double ecc = ECC;
! 1398:
! 1399: static int
! 1400: Xspmercator(struct place *place, double *x, double *y)
! 1401: {
! 1402: if(Xmercator(place,x,y) < 0)
! 1403: return(-1);
! 1404: *y += 0.5*ecc*log((1-ecc*place->nlat.s)/(1+ecc*place->nlat.s));
! 1405: return(1);
! 1406: }
! 1407:
! 1408: proj
! 1409: sp_mercator(void)
! 1410: {
! 1411: return(Xspmercator);
! 1412: }
! 1413: 0707070035051030621006640000030000040000011402160533472767500002300000000616libmap/mollweide.c #include "map.h"
! 1414:
! 1415: static int
! 1416: Xmollweide(struct place *place, double *x, double *y)
! 1417: {
! 1418: double z;
! 1419: double w;
! 1420: z = place->nlat.l;
! 1421: if(fabs(z)<89.9*RAD)
! 1422: do { /*newton for 2z+sin2z=pi*sin(lat)*/
! 1423: w = (2*z+sin(2*z)-PI*place->nlat.s)/(2+2*cos(2*z));
! 1424: z -= w;
! 1425: } while(fabs(w)>=.00001);
! 1426: *y = sin(z);
! 1427: *x = - (2/PI)*cos(z)*place->wlon.l;
! 1428: return(1);
! 1429: }
! 1430:
! 1431: proj
! 1432: mollweide(void)
! 1433: {
! 1434: return(Xmollweide);
! 1435: }
! 1436: 0707070035051030521006640000030000040000010633000533472767600002300000000570libmap/newyorker.c #include "map.h"
! 1437:
! 1438: static double a;
! 1439:
! 1440: static int
! 1441: Xnewyorker(struct place *place, double *x, double *y)
! 1442: {
! 1443: double r = PI/2 - place->nlat.l;
! 1444: double s;
! 1445: if(r<.001) /* cheat to plot center */
! 1446: s = 0;
! 1447: else if(r<a)
! 1448: return -1;
! 1449: else
! 1450: s = log(r/a);
! 1451: *x = -s * place->wlon.s;
! 1452: *y = -s * place->wlon.c;
! 1453: return(1);
! 1454: }
! 1455:
! 1456: proj
! 1457: newyorker(double a0)
! 1458: {
! 1459: a = a0*RAD;
! 1460: return(Xnewyorker);
! 1461: }
! 1462: 0707070035051030611006640000030000040000011402220533472767600002600000000370libmap/orthographic.c #include "map.h"
! 1463:
! 1464: int
! 1465: Xorthographic(struct place *place, double *x, double *y)
! 1466: {
! 1467: *x = - place->nlat.c * place->wlon.s;
! 1468: *y = - place->nlat.c * place->wlon.c;
! 1469: return(place->nlat.l<0.? 0 : 1);
! 1470: }
! 1471:
! 1472: proj
! 1473: orthographic(void)
! 1474: {
! 1475: return(Xorthographic);
! 1476: }
! 1477: 0707070035051030601006640000030000040000011402240533472767600002500000001435libmap/perspective.c #include "map.h"
! 1478:
! 1479: static double viewpt;
! 1480:
! 1481: static int
! 1482: Xperspective(struct place *place, double *x, double *y)
! 1483: {
! 1484: double r;
! 1485: if(viewpt<=1. && place->nlat.s<=viewpt+.01)
! 1486: return(-1);
! 1487: r = place->nlat.c*(viewpt - 1.)/(viewpt - place->nlat.s);
! 1488: *x = - r*place->wlon.s;
! 1489: *y = - r*place->wlon.c;
! 1490: if(r>4.)
! 1491: return(0);
! 1492: if(fabs(viewpt)>1. && place->nlat.s<=1./viewpt)
! 1493: return(0);
! 1494: return(1);
! 1495: }
! 1496:
! 1497: int
! 1498: Xstereographic(struct place *place, double *x, double *y)
! 1499: {
! 1500: viewpt = -1;
! 1501: return(Xperspective(place, x, y));
! 1502: }
! 1503:
! 1504: proj
! 1505: perspective(double radius)
! 1506: {
! 1507: viewpt = radius;
! 1508: if(radius>=1000.)
! 1509: return(Xorthographic);
! 1510: if(fabs(radius-1.)<.0001)
! 1511: return(0);
! 1512: return(Xperspective);
! 1513: }
! 1514:
! 1515: proj
! 1516: stereographic(void)
! 1517: {
! 1518: viewpt = -1.;
! 1519: return(Xperspective);
! 1520: }
! 1521:
! 1522: proj
! 1523: gnomonic(void)
! 1524: {
! 1525: viewpt = 0.;
! 1526: return(Xperspective);
! 1527: }
! 1528: 0707070035051030571006640000030000040000011402260533472767700002300000001047libmap/polyconic.c #include "map.h"
! 1529:
! 1530: int
! 1531: Xpolyconic(struct place *place, double *x, double *y)
! 1532: {
! 1533: double r, alpha;
! 1534: double lat2, lon2;
! 1535: if(fabs(place->nlat.l) > .01) {
! 1536: r = place->nlat.c / place->nlat.s;
! 1537: alpha = place->wlon.l * place->nlat.s;
! 1538: *y = place->nlat.l + r*(1 - cos(alpha));
! 1539: *x = - r*sin(alpha);
! 1540: } else {
! 1541: lon2 = place->wlon.l * place->wlon.l;
! 1542: lat2 = place->nlat.l * place->nlat.l;
! 1543: *y = place->nlat.l * (1+(lon2/2)*(1-(8+lon2)*lat2/12));
! 1544: *x = - place->wlon.l * (1-lat2*(3+lon2)/6);
! 1545: }
! 1546: return(1);
! 1547: }
! 1548:
! 1549: proj
! 1550: polyconic(void)
! 1551: {
! 1552: return(Xpolyconic);
! 1553: }
! 1554: 0707070035051030561006640000030000040000011402300533472770000002500000000426libmap/rectangular.c #include "map.h"
! 1555:
! 1556: static double scale;
! 1557:
! 1558: static int
! 1559: Xrectangular(struct place *place, double *x, double *y)
! 1560: {
! 1561: *x = -scale*place->wlon.l;
! 1562: *y = place->nlat.l;
! 1563: return(1);
! 1564: }
! 1565:
! 1566: proj
! 1567: rectangular(double par)
! 1568: {
! 1569: scale = cos(par*RAD);
! 1570: if(scale<.1)
! 1571: return 0;
! 1572: return(Xrectangular);
! 1573: }
! 1574: 0707070035051030361006640000030000040000011402320533472770000002500000001143libmap/simpleconic.c #include "map.h"
! 1575:
! 1576: static double r0, a;
! 1577:
! 1578: static int
! 1579: Xsimpleconic(struct place *place, double *x, double *y)
! 1580: {
! 1581: double r = r0 - place->nlat.l;
! 1582: double t = a*place->wlon.l;
! 1583: *x = -r*sin(t);
! 1584: *y = -r*cos(t);
! 1585: return 1;
! 1586: }
! 1587:
! 1588: proj
! 1589: simpleconic(double par0, double par1)
! 1590: {
! 1591: struct coord lat0;
! 1592: struct coord lat1;
! 1593: deg2rad(par0,&lat0);
! 1594: deg2rad(par1,&lat1);
! 1595: if(fabs(lat0.l+lat1.l)<.01)
! 1596: return rectangular(par0);
! 1597: if(fabs(lat0.l-lat1.l)<.01) {
! 1598: a = lat0.s/lat0.l;
! 1599: r0 = lat0.c/lat0.s + lat0.l;
! 1600: } else {
! 1601: a = (lat1.c-lat0.c)/(lat0.l-lat1.l);
! 1602: r0 = ((lat0.c+lat1.c)/a + lat1.l + lat0.l)/2;
! 1603: }
! 1604: return Xsimpleconic;
! 1605: }
! 1606: 0707070035051030551006640000030000040000011402360533472770000002400000000312libmap/sinusoidal.c #include "map.h"
! 1607:
! 1608: int
! 1609: Xsinusoidal(struct place *place, double *x, double *y)
! 1610: {
! 1611: *x = - place->wlon.l * place->nlat.c;
! 1612: *y = place->nlat.l;
! 1613: return(1);
! 1614: }
! 1615:
! 1616: proj
! 1617: sinusoidal(void)
! 1618: {
! 1619: return(Xsinusoidal);
! 1620: }
! 1621: 0707070035051030541006640000030000040000011402400533472770100001700000011457libmap/tetra.c #include "map.h"
! 1622:
! 1623: /*
! 1624: * conformal map of earth onto tetrahedron
! 1625: * the stages of mapping are
! 1626: * (a) stereo projection of tetrahedral face onto
! 1627: * isosceles curvilinear triangle with 3 120-degree
! 1628: * angles and one straight side
! 1629: * (b) map of this triangle onto half plane cut along
! 1630: * 3 rays from the roots of unity to infinity
! 1631: * formula (z^4+2*3^.5*z^2-1)/(z^4-2*3^.5*z^2-1)
! 1632: * (c) do 3 times for each sector of plane:
! 1633: * map of |arg z|<=pi/6, cut along z>1 into
! 1634: * triangle |arg z|<=pi/6, Re z<=const,
! 1635: * with upper side of cut going into upper half of
! 1636: * of vertical side of triangle and lowere into lower
! 1637: * formula int from 0 to z dz/sqrt(1-z^3)
! 1638: *
! 1639: * int from u to 1 3^.25*du/sqrt(1-u^3) =
! 1640: F(acos((rt3-1+u)/(rt3+1-u)),sqrt(1/2+rt3/4))
! 1641: * int from 1 to u 3^.25*du/sqrt(u^3-1) =
! 1642: * F(acos((rt3+1-u)/(rt3-1+u)),sqrt(1/2-rt3/4))
! 1643: * this latter formula extends analytically down to
! 1644: * u=0 and is the basis of this routine, with the
! 1645: * argument of complex elliptic integral elco2
! 1646: * being tan(acos...)
! 1647: * the formula F(pi-x,k) = 2*F(pi/2,k)-F(x,k) is
! 1648: * used to cross over into the region where Re(acos...)>pi/2
! 1649: * f0 and fpi are suitably scaled complete integrals
! 1650: */
! 1651:
! 1652: #define TFUZZ 0.00001
! 1653:
! 1654: static struct place tpole[4]; /* point of tangency of tetrahedron face*/
! 1655: static double tpoleinit[4][2] = {
! 1656: 1., 0.,
! 1657: 1., 180.,
! 1658: -1., 90.,
! 1659: -1., -90.
! 1660: };
! 1661: static struct tproj {
! 1662: double tlat,tlon; /* center of stereo projection*/
! 1663: double ttwist; /* rotatn before stereo*/
! 1664: double trot; /*rotate after projection*/
! 1665: struct place projpl; /*same as tlat,tlon*/
! 1666: struct coord projtw; /*same as ttwist*/
! 1667: struct coord postrot; /*same as trot*/
! 1668: } tproj[4][4] = {
! 1669: {/*00*/ {0.},
! 1670: /*01*/ {90., 0., 90., -90.},
! 1671: /*02*/ {0., 45., -45., 150.},
! 1672: /*03*/ {0., -45., -135., 30.}
! 1673: },
! 1674: {/*10*/ {90., 0., -90., 90.},
! 1675: /*11*/ {0.},
! 1676: /*12*/ {0., 135., -135., -150.},
! 1677: /*13*/ {0., -135., -45., -30.}
! 1678: },
! 1679: {/*20*/ {0., 45., 135., -30.},
! 1680: /*21*/ {0., 135., 45., -150.},
! 1681: /*22*/ {0.},
! 1682: /*23*/ {-90., 0., 180., 90.}
! 1683: },
! 1684: {/*30*/ {0., -45., 45., -150.},
! 1685: /*31*/ {0., -135., 135., -30.},
! 1686: /*32*/ {-90., 0., 0., 90.},
! 1687: /*33*/ {0.}
! 1688: }};
! 1689: static double tx[4] = { /*where to move facet after final rotation*/
! 1690: 0., 0., -1., 1. /*-1,1 to be sqrt(3)*/
! 1691: };
! 1692: static double ty[4] = {
! 1693: 0., 2., -1., -1.
! 1694: };
! 1695: static double root3;
! 1696: static double rt3inv;
! 1697: static double two_rt3;
! 1698: static double tkc,tk,tcon;
! 1699: static double f0r,f0i,fpir,fpii;
! 1700:
! 1701: static void
! 1702: twhichp(struct place *g, int *p, int *q)
! 1703: {
! 1704: int i,j,k;
! 1705: double cosdist[4];
! 1706: struct place *tp;
! 1707: for(i=0;i<4;i++) {
! 1708: tp = &tpole[i];
! 1709: cosdist[i] = g->nlat.s*tp->nlat.s +
! 1710: g->nlat.c*tp->nlat.c*(
! 1711: g->wlon.s*tp->wlon.s +
! 1712: g->wlon.c*tp->wlon.c);
! 1713: }
! 1714: j = 0;
! 1715: for(i=1;i<4;i++)
! 1716: if(cosdist[i] > cosdist[j])
! 1717: j = i;
! 1718: *p = j;
! 1719: k = j==0?1:0;
! 1720: for(i=0;i<4;i++)
! 1721: if(i!=j&&cosdist[i]>cosdist[k])
! 1722: k = i;
! 1723: *q = k;
! 1724: }
! 1725:
! 1726: int
! 1727: Xtetra(struct place *place, double *x, double *y)
! 1728: {
! 1729: int i,j;
! 1730: struct place pl;
! 1731: register struct tproj *tpp;
! 1732: double vr, vi;
! 1733: double br, bi;
! 1734: double zr,zi,z2r,z2i,z4r,z4i,sr,si,tr,ti;
! 1735: twhichp(place,&i,&j);
! 1736: copyplace(place,&pl);
! 1737: norm(&pl,&tproj[i][j].projpl,&tproj[i][j].projtw);
! 1738: Xstereographic(&pl,&vr,&vi);
! 1739: zr = vr/2;
! 1740: zi = vi/2;
! 1741: if(zr<=TFUZZ)
! 1742: zr = TFUZZ;
! 1743: csq(zr,zi,&z2r,&z2i);
! 1744: csq(z2r,z2i,&z4r,&z4i);
! 1745: z2r *= two_rt3;
! 1746: z2i *= two_rt3;
! 1747: cdiv(z4r+z2r-1,z4i+z2i,z4r-z2r-1,z4i-z2i,&sr,&si);
! 1748: csqrt(sr-1,si,&tr,&ti);
! 1749: cdiv(tcon*tr,tcon*ti,root3+1-sr,-si,&br,&bi);
! 1750: if(br<0) {
! 1751: br = -br;
! 1752: bi = -bi;
! 1753: if(!elco2(br,bi,tk,1.,1.,&vr,&vi))
! 1754: return 0;
! 1755: vr = fpir - vr;
! 1756: vi = fpii - vi;
! 1757: } else
! 1758: if(!elco2(br,bi,tk,1.,1.,&vr,&vi))
! 1759: return 0;
! 1760: if(si>=0) {
! 1761: tr = f0r - vi;
! 1762: ti = f0i + vr;
! 1763: } else {
! 1764: tr = f0r + vi;
! 1765: ti = f0i - vr;
! 1766: }
! 1767: tpp = &tproj[i][j];
! 1768: *x = tr*tpp->postrot.c +
! 1769: ti*tpp->postrot.s + tx[i];
! 1770: *y = ti*tpp->postrot.c -
! 1771: tr*tpp->postrot.s + ty[i];
! 1772: return(1);
! 1773: }
! 1774:
! 1775: int
! 1776: tetracut(struct place *g, struct place *og, double *cutlon)
! 1777: {
! 1778: int i,j,k;
! 1779: if((g->nlat.s<=-rt3inv&&og->nlat.s<=-rt3inv) &&
! 1780: (ckcut(g,og,*cutlon=0.)==2||ckcut(g,og,*cutlon=PI)==2))
! 1781: return(2);
! 1782: twhichp(g,&i,&k);
! 1783: twhichp(og,&j,&k);
! 1784: if(i==j||i==0||j==0)
! 1785: return(1);
! 1786: return(0);
! 1787: }
! 1788:
! 1789: proj
! 1790: tetra(void)
! 1791: {
! 1792: register i;
! 1793: int j;
! 1794: register struct place *tp;
! 1795: register struct tproj *tpp;
! 1796: double t;
! 1797: root3 = sqrt(3.);
! 1798: rt3inv = 1/root3;
! 1799: two_rt3 = 2*root3;
! 1800: tkc = sqrt(.5-.25*root3);
! 1801: tk = sqrt(.5+.25*root3);
! 1802: tcon = 2*sqrt(root3);
! 1803: elco2(tcon/(root3-1),0.,tkc,1.,1.,&f0r,&f0i);
! 1804: elco2(1.e15,0.,tk,1.,1.,&fpir,&fpii);
! 1805: fpir *= 2;
! 1806: fpii *= 2;
! 1807: for(i=0;i<4;i++) {
! 1808: tx[i] *= f0r*root3;
! 1809: ty[i] *= f0r;
! 1810: tp = &tpole[i];
! 1811: t = tp->nlat.s = tpoleinit[i][0]/root3;
! 1812: tp->nlat.c = sqrt(1 - t*t);
! 1813: tp->nlat.l = atan2(tp->nlat.s,tp->nlat.c);
! 1814: deg2rad(tpoleinit[i][1],&tp->wlon);
! 1815: for(j=0;j<4;j++) {
! 1816: tpp = &tproj[i][j];
! 1817: latlon(tpp->tlat,tpp->tlon,&tpp->projpl);
! 1818: deg2rad(tpp->ttwist,&tpp->projtw);
! 1819: deg2rad(tpp->trot,&tpp->postrot);
! 1820: }
! 1821: }
! 1822: return(Xtetra);
! 1823: }
! 1824:
! 1825: 0707070035051030711006640000030000040000011402440533472770100002500000001041libmap/trapezoidal.c #include "map.h"
! 1826:
! 1827: static struct coord stdpar0, stdpar1;
! 1828: static double k;
! 1829: static double yeq;
! 1830:
! 1831: static int
! 1832: Xtrapezoidal(struct place *place, double *x, double *y)
! 1833: {
! 1834: *y = yeq + place->nlat.l;
! 1835: *x = *y*k*place->wlon.l;
! 1836: return 1;
! 1837: }
! 1838:
! 1839: proj
! 1840: trapezoidal(double par0, double par1)
! 1841: {
! 1842: if(fabs(fabs(par0)-fabs(par1))<.1)
! 1843: return rectangular(par0);
! 1844: deg2rad(par0,&stdpar0);
! 1845: deg2rad(par1,&stdpar1);
! 1846: if(fabs(par1-par0) < .1)
! 1847: k = stdpar1.s;
! 1848: else
! 1849: k = (stdpar1.c-stdpar0.c)/(stdpar0.l-stdpar1.l);
! 1850: yeq = -stdpar1.l - stdpar1.c/k;
! 1851: return Xtrapezoidal;
! 1852: }
! 1853: 0707070035051030431006640000030000040000011400500533472770100002100000003153libmap/twocirc.c #include "map.h"
! 1854:
! 1855: static double
! 1856: quadratic(double a, double b, double c)
! 1857: {
! 1858: double disc = b*b - 4*a*c;
! 1859: return disc<0? 0: (-b-sqrt(disc))/(2*a);
! 1860: }
! 1861:
! 1862: /* for projections with meridians being circles centered
! 1863: on the x axis and parallels being circles centered on the y
! 1864: axis. Find the intersection of the meridian thru (m,0), (90,0),
! 1865: with the parallel thru (0,p), (p1,p2) */
! 1866:
! 1867: static int
! 1868: twocircles(double m, double p, double p1, double p2, double *x, double *y)
! 1869: {
! 1870: double a; /* center of meridian circle, a>0 */
! 1871: double b; /* center of parallel circle, b>0 */
! 1872: double t,bb;
! 1873: if(m > 0) {
! 1874: twocircles(-m,p,p1,p2,x,y);
! 1875: *x = -*x;
! 1876: } else if(p < 0) {
! 1877: twocircles(m,-p,p1,-p2,x,y);
! 1878: *y = -*y;
! 1879: } else if(p < .01) {
! 1880: *x = m;
! 1881: t = m/p1;
! 1882: *y = p + (p2-p)*t*t;
! 1883: } else if(m > -.01) {
! 1884: *y = p;
! 1885: *x = m - m*p*p;
! 1886: } else {
! 1887: b = p>=1? 1: p>.99? 0.5*(p+1 + p1*p1/(1-p)):
! 1888: 0.5*(p*p-p1*p1-p2*p2)/(p-p2);
! 1889: a = .5*(m - 1/m);
! 1890: t = m*m-p*p+2*(b*p-a*m);
! 1891: bb = b*b;
! 1892: *x = quadratic(1+a*a/bb, -2*a + a*t/bb,
! 1893: t*t/(4*bb) - m*m + 2*a*m);
! 1894: *y = (*x*a+t/2)/b;
! 1895: }
! 1896: return 1;
! 1897: }
! 1898:
! 1899: static int
! 1900: Xglobular(struct place *place, double *x, double *y)
! 1901: {
! 1902: twocircles(-2*place->wlon.l/PI,
! 1903: 2*place->nlat.l/PI, place->nlat.c, place->nlat.s, x, y);
! 1904: return 1;
! 1905: }
! 1906:
! 1907: proj
! 1908: globular(void)
! 1909: {
! 1910: return Xglobular;
! 1911: }
! 1912:
! 1913: static int
! 1914: Xvandergrinten(struct place *place, double *x, double *y)
! 1915: {
! 1916: double t = 2*place->nlat.l/PI;
! 1917: double abst = fabs(t);
! 1918: double pval = abst>=1? 1: abst/(1+sqrt(1-t*t));
! 1919: double p2 = 2*pval/(1+pval);
! 1920: twocircles(-place->wlon.l/PI, pval, sqrt(1-p2*p2), p2, x, y);
! 1921: if(t < 0)
! 1922: *y = -*y;
! 1923: return 1;
! 1924: }
! 1925:
! 1926: proj
! 1927: vandergrinten(void)
! 1928: {
! 1929: return Xvandergrinten;
! 1930: }
! 1931: 0707070035051030531006640000030000040000011402460533501406400002000000004626libmap/zcoord.c #include "map.h"
! 1932:
! 1933: static double cirmod(double);
! 1934:
! 1935: static struct place pole; /* map pole is tilted to here */
! 1936: static struct coord twist; /* then twisted this much */
! 1937: static struct place ipole; /* inverse transfrom */
! 1938: static struct coord itwist;
! 1939:
! 1940: void
! 1941: orient(double lat, double lon, double theta)
! 1942: {
! 1943: lat = cirmod(lat);
! 1944: if(lat>90.) {
! 1945: lat = 180. - lat;
! 1946: lon -= 180.;
! 1947: theta -= 180.;
! 1948: } else if(lat < -90.) {
! 1949: lat = -180. - lat;
! 1950: lon -= 180.;
! 1951: theta -= 180.;
! 1952: }
! 1953: latlon(lat,lon,&pole);
! 1954: deg2rad(theta, &twist);
! 1955: latlon(lat,180.-theta,&ipole);
! 1956: deg2rad(180.-lon, &itwist);
! 1957: }
! 1958:
! 1959: void
! 1960: latlon(double lat, double lon, struct place *p)
! 1961: {
! 1962: lat = cirmod(lat);
! 1963: if(lat>90.) {
! 1964: lat = 180. - lat;
! 1965: lon -= 180.;
! 1966: } else if(lat < -90.) {
! 1967: lat = -180. - lat;
! 1968: lon -= 180.;
! 1969: }
! 1970: deg2rad(lat,&p->nlat);
! 1971: deg2rad(lon,&p->wlon);
! 1972: }
! 1973:
! 1974: void
! 1975: deg2rad(double theta, struct coord *coord)
! 1976: {
! 1977: theta = cirmod(theta);
! 1978: coord->l = theta*RAD;
! 1979: if(theta==90) {
! 1980: coord->s = 1;
! 1981: coord->c = 0;
! 1982: } else if(theta== -90) {
! 1983: coord->s = -1;
! 1984: coord->c = 0;
! 1985: } else
! 1986: sincos(coord);
! 1987: }
! 1988:
! 1989: static double
! 1990: cirmod(double theta)
! 1991: {
! 1992: while(theta >= 180.)
! 1993: theta -= 360;
! 1994: while(theta<-180.)
! 1995: theta += 360.;
! 1996: return(theta);
! 1997: }
! 1998:
! 1999: void
! 2000: sincos(struct coord *coord)
! 2001: {
! 2002: coord->s = sin(coord->l);
! 2003: coord->c = cos(coord->l);
! 2004: }
! 2005:
! 2006: void
! 2007: normalize(struct place *gg)
! 2008: {
! 2009: norm(gg,&pole,&twist);
! 2010: }
! 2011:
! 2012: void
! 2013: invert(struct place *g)
! 2014: {
! 2015: norm(g,&ipole,&itwist);
! 2016: }
! 2017:
! 2018: void
! 2019: norm(struct place *gg, struct place *pp, struct coord *tw)
! 2020: {
! 2021: register struct place *g; /*geographic coords */
! 2022: register struct place *p; /* new pole in old coords*/
! 2023: struct place m; /* standard map coords*/
! 2024: g = gg;
! 2025: p = pp;
! 2026: if(p->nlat.s == 1.) {
! 2027: if(p->wlon.l+tw->l == 0.)
! 2028: return;
! 2029: g->wlon.l -= p->wlon.l+tw->l;
! 2030: } else {
! 2031: if(p->wlon.l != 0) {
! 2032: g->wlon.l -= p->wlon.l;
! 2033: sincos(&g->wlon);
! 2034: }
! 2035: m.nlat.s = p->nlat.s * g->nlat.s
! 2036: + p->nlat.c * g->nlat.c * g->wlon.c;
! 2037: m.nlat.c = sqrt(1. - m.nlat.s * m.nlat.s);
! 2038: m.nlat.l = atan2(m.nlat.s, m.nlat.c);
! 2039: m.wlon.s = g->nlat.c * g->wlon.s;
! 2040: m.wlon.c = p->nlat.c * g->nlat.s
! 2041: - p->nlat.s * g->nlat.c * g->wlon.c;
! 2042: m.wlon.l = atan2(m.wlon.s, - m.wlon.c)
! 2043: - tw->l;
! 2044: *g = m;
! 2045: }
! 2046: sincos(&g->wlon);
! 2047: if(g->wlon.l>PI)
! 2048: g->wlon.l -= 2*PI;
! 2049: else if(g->wlon.l<-PI)
! 2050: g->wlon.l += 2*PI;
! 2051: }
! 2052:
! 2053: void
! 2054: printp(struct place *g)
! 2055: {
! 2056: printf("%.3f %.3f %.3f %.3f %.3f %.3f\n",
! 2057: g->nlat.l,g->nlat.s,g->nlat.c,g->wlon.l,g->wlon.s,g->wlon.c);
! 2058: }
! 2059:
! 2060: void
! 2061: copyplace(struct place *g1, struct place *g2)
! 2062: {
! 2063: *g2 = *g1;
! 2064: }
! 2065: 0707070035050172611006640000030000040000010336560447134002200000600000003321map.5