sampson@killer.DALLAS.TX.US (Steve Sampson) (12/07/88)
Posting-number: Volume 5, Issue 72
Submitted-by: "Steve Sampson" <sampson@killer.DALLAS.TX.US>
Archive-name: mapper
This is for Turbo-C users, no binaries included. Easily modified for
other compilers.
#! /bin/sh
# Contents: mapper.doc mapper.c
# Wrapped by sampson on Mon Dec 05 15:07:18 1988
PATH=/bin:/usr/bin:/usr/ucb ; export PATH
if test -f mapper.doc -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"mapper.doc\"
else
echo shar: Extracting \"mapper.doc\" \(4532 characters\)
sed "s/^X//" >mapper.doc <<'END_OF_mapper.doc'
X
X The Mapper Program
X
X Version 1.7
X
X 22 November 1988
X
X
X I was interested in techniques for producing maps, and found
X the article by William D. Johnston in the May and June 1979 Byte
X Magazine. This two part article provided an excellent
X introduction and source code in Basic Language. His code was
X restricted to the algorithms and did not get involved with user
X interface. To evaluate his algorithms and try out the displays I
X coded the program and a simple interface in Turbo-C Version 1.5.
X The program in its current form is highly based on Mr. Johnston's
X algorithms and provides no significant additional capabilities.
X
X I also found a high resolution database called the Micro World
X Data Bank II (MWDBII). This database was 1 megabyte in length
X and good down to minutes of a degree. See the C source code
X comments for availability.
X
X To run the program and receive help you use the DOS common
X method of the question option "/?". Just type "mapper/?" and the
X following usage help will be displayed:
X
X Usage: mapper [/bcdgilmrsx]
X
X /b Boundaries Off
X /c Countries On
X /dn Database ('MWDBII' Default)
X /g Grid lines On
X /i Islands Off
X /l Lakes Off
X /mn Map Resolution (5 Default)
X /r Rivers On
X /s States On
X /x Colors On
X
X Defaults to Boundaries and Islands On
X
X The defaults are what I thought should be fairly common. The map
X database has 5 resolutions, and can be selected with the 'm'
X option. 5 is the lowest resolution and 1 is the greatest. If
X you have several different databases you can use the 'd' option
X and provide the path and filename (128 Characters max). The 'm'
X and 'd' options should be placed at the end. They can be put
X anywhere but it's a little easier at the end. Example:
X mapper/glrsm1. If you use the option in the middle you will need
X to put a space between it and the remaining options. Example:
X mapper/glddata /rs. These are the most foolproof methods. Note:
X The level 5 database included doesn't really use the options yet.
X The program works as advertised on level 1. There are some errors
X with the database as you'll see. I've converted the database to
X ASCII, and am working on cleaning up the errors and redundancies.
X
X A little about the speed of the result. The program is quite
X slow on an 8088 without a math coprocessor, and speed is getting
X acceptable on an 80286. The C language standard uses double
X precision math. This is a waste with the current database
X resolution. An integer version of the math routines would sure
X speed up the program quite a bit. The mapper program uses
X Turbo-C auto detect of a math coprocessor and graphics device
X type (CGA, EGA, and VGA).
X
X If you want to quit the plotting on the screen, just hit any
X key and the bell will sound and exit you back to DOS. You can
X also use Control-C to get out.
X
X The C program lists three sources for the Micro World Data
X Bank II database. The database is 1 Megabyte (800K Compressed)
X and is just too much data for reasonable downloading. To see if
X the program would be useful for you I included a level 5
X resolution map (the lowest resolution). This particular database
X has all the data thrown together so the command line options
X aren't fully functional. They become functional at about level
X 3 I believe.
X
X This program was tested on a PC XT Clone, 640K, and NSI EGA
X board. Also a Zenith Z-248 with 512k and CGA was tested. Other
X configurations will need to be tested.
X
X Due to the grid method used, it shouldn't be used with an
X Azimuthal Equidistant map. You can try it once to see what I
X mean. There's lots of room for improvements, "Handle It!".
X
X Thanks to Mr. Johnston for his article and algorithms.
X
X USMail: Steve R. Sampson, Box 45668, Tinker AFB, Oklahoma, 73145
X Compuserve: 75136,626 Unix: sampson@killer.dallas.tx
END_OF_mapper.doc
if test 4532 -ne `wc -c <mapper.doc`; then
echo shar: \"mapper.doc\" unpacked with wrong size!
fi
# end of overwriting check
fi
if test -f mapper.c -a "${1}" != "-c" ; then
echo shar: Will not over-write existing file \"mapper.c\"
else
echo shar: Extracting \"mapper.c\" \(13949 characters\)
sed "s/^X//" >mapper.c <<'END_OF_mapper.c'
X/*
X * mapper.c
X *
X * Version 1.7 by Steven R. Sampson, November 1988
X *
X * Based on a program and article by William D. Johnston
X * Copyright (c) May-June 1979 BYTE, All Rights Reserved
X *
X * This program draws three types of map projections:
X * Perspective, Modified Perspective, and Azimuthal Equidistant.
X *
X * Compiled with Turbo-C V1.5
X */
X
X#include <dos.h>
X#include <math.h>
X#include <stdio.h>
X#include <conio.h>
X#include <string.h>
X#include <stdlib.h>
X#include <graphics.h>
X
Xtypedef int bool;
X
X/* Program Constants */
X
X#define FALSE (bool) 0
X#define TRUE (bool) ~FALSE
X
X#define PI (3.141593F)
X#define HALFPI (1.570796F)
X#define TWOPI (2.0F * PI) /* Two Pi alias 360 Degrees */
X
X#define RADIAN (180.0F / PI ) /* One radian */
X#define TWO (2.0F / RADIAN) /* 2 degrees in radians */
X#define TEN (10.0F / RADIAN) /* 10 degrees in radians */
X#define EIGHTY (80.0F / RADIAN) /* 80 degrees in radians */
X#define EARTH (6378.0F) /* Mean radius of earth (Kilometers) */
X
X/* Program Globals */
X
XFILE *fp;
X
Xfloat angle, maxplot, center_lat, center_lon, lat, lon, distance,
X sin_of_distance, cos_of_distance, sin_of_center_lat, cos_of_center_lat,
X scale, g, h2, facing_azimuth, aspect;
X
Xint option, center_x, center_y, grid_color, level = 5;
Xint GraphDriver = DETECT, GraphMode;
X
Xchar optstring[] = "bcd:gilm:rsx?";
Xchar database[128] = "mwdbii"; /* default name 'MWDBII' */
X /* leave room for pathname! */
Xbool boundaries = TRUE, /* defaults to Boundaries, Islands */
X countries = FALSE,
X grids = FALSE,
X islands = TRUE,
X lakes = FALSE,
X rivers = FALSE,
X states = FALSE,
X colors = FALSE;
X
X/* Forward Declarations, Prototypes */
X
Xextern int getopt(int, char **, char *);
Xextern int optind, opterr;
Xextern char *optarg;
X
Xfloat parse(char *);
Xvoid grid(void), plotmap(void), prompts(void), quit(void);
Xbool compute(float *, float *, int *, int *);
X
X
Xmain(argc, argv)
Xint argc;
Xchar *argv[];
X{
X register int i;
X int err, xasp, yasp;
X
X registerbgidriver(EGAVGA_driver);
X registerbgidriver(CGA_driver);
X
X setcbrk(TRUE); /* Allow Control-C checking */
X ctrlbrk(quit); /* Execute quit() if Control-C detected */
X
X while ((i = getopt(argc, argv, optstring)) != -1) {
X switch (i) {
X case 'b':
X boundaries = FALSE;
X break;
X case 'c':
X countries = TRUE;
X break;
X case 'd':
X strcpy(database, optarg);
X break;
X case 'g':
X grids = TRUE;
X break;
X case 'i':
X islands = FALSE;
X break;
X case 'l':
X lakes = TRUE;
X break;
X case 'm':
X level = atoi(optarg);
X break;
X case 'r':
X rivers = TRUE;
X break;
X case 's':
X states = TRUE;
X break;
X case 'x':
X colors = FALSE;
X break;
X case '?':
X default:
X printf("Usage: mapper [/bcdgilmrsx]\n\n");
X printf(" /b Boundaries Off\n");
X printf(" /c Countries On\n");
X printf(" /dn Database ('MWDBII' Default)\n");
X printf(" /g Grid lines On\n");
X printf(" /i Islands Off\n");
X printf(" /l Lakes On\n");
X printf(" /mn Map Resolution (5 Default)\n");
X printf(" /r Rivers On\n");
X printf(" /s States On\n");
X printf(" /x Colors On\n\n");
X printf("Defaults to Boundaries and Islands On\n");
X exit(0);
X }
X }
X
X if ((fp = fopen(database, "rb")) == (FILE *)NULL) {
X printf("\007Error: Can't locate Database '%s'\n", database);
X exit(1);
X }
X
X initgraph(&GraphDriver, &GraphMode, "");/* initialize graphics */
X err = graphresult();
X
X restorecrtmode(); /* get back to text mode */
X
X if (err < 0) {
X printf("Graphics Error - %s\n", grapherrormsg(err));
X exit(-1);
X }
X
X center_x = getmaxx() / 2; /* get screen size for x, y */
X center_y = getmaxy() / 2;
X getaspectratio(&xasp, &yasp); /* squish factor for y axis */
X aspect = (float)xasp / (float)yasp;
X
X prompts(); /* get the basic map info */
X setgraphmode(GraphMode); /* and go to graphics mode */
X
X if (GraphMode != CGAHI) {
X setbkcolor(BLACK); /* must be EGA or VGA then */
X if (colors)
X grid_color = EGA_LIGHTRED;
X else
X grid_color = EGA_LIGHTGRAY;
X } else
X grid_color = LIGHTGRAY; /* CGA only has two colors */
X
X setcolor(LIGHTGRAY);
X
X /*
X * See if data plotting is even needed
X */
X
X if (boundaries || countries || islands || lakes || rivers || states)
X plotmap(); /* display map on screen */
X
X if (grids)
X grid(); /* draw lat & long ref lines */
X
X if (print)
X printscreen(); /* relay screen to printer */
X
X sound(800); /* 800 Hz for 1/4 a second */
X delay(125);
X nosound();
X
X getch(); /* pause until key pressed */
X closegraph(); /* graphics off */
X fclose(fp); /* close database file */
X
X exit(0);
X}
X
X/*
X * Return to operator following Control-C
X */
X
Xvoid quit()
X{
X closegraph();
X fclose(fp);
X
X exit(0);
X}
X
X/*
X * Operator prompts and input.
X */
X
Xvoid prompts()
X{
X char temp[16];
X float ret, altitude;
X
X printf("West Longitudes and South Lattitudes are negative\n");
X
X /* input the world Lat & Long that is to be centered on */
X /* then convert the human notation to radians */
X
X do {
X printf("\nLatitude of the map center [+-]dd.mm : ");
X scanf("%s", temp);
X ret = parse(temp);
X } while (ret > 90.0F || ret < -90.0F);
X
X /* the arcosine function has trouble at 90 degrees */
X
X if (ret == 90.0F)
X ret = 89.9F;
X
X if (ret == -90.0F)
X ret = -89.9F;
X
X center_lat = ret / RADIAN;
X sin_of_center_lat = sin(center_lat);
X cos_of_center_lat = cos(center_lat);
X
X do {
X printf("Longitude of the map center [+-]ddd.mm : ");
X scanf("%s", temp);
X ret = parse(temp);
X } while (ret > 180.0F || ret < -180.0F);
X
X center_lon = ret / RADIAN;
X
X do {
X printf("\nSelect from the following options:\n\n");
X printf(" 1 - Perspective Projection\n");
X printf(" 2 - Modified Perspective Projection\n");
X printf(" 3 - Azimuthal Equidistant Projection\n\n");
X printf("Choice : ");
X scanf("%d", &option);
X } while (option < 1 || option > 3);
X
X if (option == 3) {
X maxplot = PI; /* use HALFPI for less area */
X scale = (float)center_y / maxplot;
X return;
X }
X
X /* input the height above the terrain */
X
X printf("\nObserver altitude (km) : ");
X scanf("%f", &altitude);
X
X h2 = EARTH + altitude;
X maxplot = acos(EARTH / h2);
X
X /* the operator can orient the world upside down if they want */
X
X do {
X printf("Observer facing azimuth (0 - 359 degrees) : ");
X scanf("%f", &facing_azimuth);
X } while (facing_azimuth < 0.0F || facing_azimuth > 360.0F);
X
X facing_azimuth = -facing_azimuth / RADIAN;
X
X /* calculate the scale for the polar coordinates */
X
X scale = (float)center_y / (EARTH * sin(maxplot));
X
X /* for the perspective projection */
X
X g = EARTH * (h2 - EARTH * cos(maxplot));
X}
X
X
X/*
X * Convert the database to the desired projection by computation.
X *
X * This code is a hand translation from BASIC to C based on Mr. Johnstons
X * code. It is a non-mathematicians idea of what he meant.
X *
X * Return TRUE if offscale else FALSE.
X */
X
Xbool compute(lat, lon, x, y)
Xregister float *lat, *lon;
Xregister int *x, *y;
X{
X float sin_of_lat,
X cos_of_lat,
X abs_delta_lon, /* absolute value */
X delta_lon, /* x distance from center */
X delta_lat, /* y distance from center */
X temp; /* temporary storage */
X
X /* normalize the longitude to +/- PI */
X
X delta_lon = *lon - center_lon;
X
X if (delta_lon < -PI)
X delta_lon = delta_lon + TWOPI;
X
X if (delta_lon > PI)
X delta_lon = delta_lon - TWOPI;
X
X abs_delta_lon = fabs(delta_lon);
X
X /*
X * If the delta_lon is within .00015 radians of 0 then
X * the difference is considered exactly zero.
X *
X * This also simplifys the great circle bearing calculation.
X */
X
X if (abs_delta_lon <= 0.00015F) {
X delta_lat = fabs(center_lat - *lat);
X
X if (delta_lat > maxplot)
X return TRUE; /* offscale */
X
X if (*lat < center_lat)
X angle = PI;
X else
X angle = 0.0F;
X
X sin_of_distance = sin(delta_lat);
X cos_of_distance = cos(delta_lat);
X }
X
X /*
X * Check if delta_lon is within .00015 radians of PI.
X */
X
X else if (fabs(PI - abs_delta_lon) <= 0.00015F) {
X delta_lat = PI - center_lat - *lat;
X
X if (delta_lat > PI) {
X delta_lat = TWOPI - delta_lat;
X angle = PI;
X } else
X angle = 0.0F;
X
X if (delta_lat > maxplot)
X return TRUE; /* offscale */
X
X sin_of_distance = sin(delta_lat);
X cos_of_distance = cos(delta_lat);
X }
X
X /*
X * Simple calculations are out, now get cosmic.
X */
X
X else {
X sin_of_lat = sin(*lat);
X cos_of_lat = cos(*lat);
X
X cos_of_distance = sin_of_center_lat * sin_of_lat +
X cos_of_center_lat * cos_of_lat *
X cos(delta_lon);
X
X distance = acos(cos_of_distance);
X
X if (distance > maxplot)
X return TRUE; /* offscale */
X
X sin_of_distance = sin(distance);
X
X temp = (sin_of_lat - sin_of_center_lat * cos_of_distance) /
X (cos_of_center_lat * sin_of_distance);
X
X if (temp < -1.0F || temp > 1.0F)
X return TRUE; /* offscale */
X
X angle = acos(temp);
X
X if (delta_lon < 0.0F)
X angle = TWOPI - angle;
X }
X
X if (facing_azimuth != 0.0F) {
X angle = angle - facing_azimuth;
X if (angle < 0.0F)
X angle = TWOPI + angle;
X }
X
X angle = HALFPI - angle;
X
X if (angle < -PI)
X angle = angle + TWOPI;
X
X switch (option) {
X case 1:
X temp = (scale * (g * sin_of_distance)) /
X (h2 - EARTH * cos_of_distance);
X break;
X case 2:
X temp = scale * EARTH * sin_of_distance;
X break;
X case 3:
X temp = scale * distance;
X }
X
X /* convert polar to rectangular, correct for screen aspect */
X
X *x = center_x + (int)(temp * cos(angle));
X *y = center_y - (int)(temp * sin(angle) * aspect);
X
X return FALSE;
X}
X
X/*
X * Read the database and plot points or lines.
X *
X * The database is Micro World Data Bank II. It's based on the
X * CIA WDB-II tape available from NTIS. Micro WDB-II was created
X * by Micro Doc. Placed in the public domain by Fred Pospeschil
X * and Antonio Riveria. Check on availability at:
X * 1-402-291-0795 (6-9 PM Central)
X *
X * Austin Code Works has something called: The World Digitized
X * that sounds like the same thing ($30.00), 1-512-258-0785
X *
X * Lone Star Software has something called: The World Digitized
X * that sounds like the same thing ($6.00), 1-800-445-6172.
X *
X * Database is in Intel word order:
X * code_lsb, code_msb, lat_lsb, lat_msb, lon_lsb, lon_msb
X *
X * Code: Integer, two meanings:
X * 1. Detail Level (1 Highest - 5 Lowest)
X *
X * 2. Header (1xxx - 7xxx) Command Line Options
X *
X * 1xxx Boundaries /b
X * 2xxx Countries /c
X * (decimal) 4xxx States /s
X * 5xxx Islands /i
X * 6xxx Lakes /l
X * 7xxx Rivers /r
X *
X * Lat & Long: Integer
X * Representing Minutes of degree
X */
X
Xvoid plotmap()
X{
X struct { short code, lat, lon; } coord;
X float lat, lon;
X int x, y;
X bool point;
X
X point = TRUE;
X while (fread(&coord, sizeof coord, 1, fp) > 0) {
X
X if (kbhit()) {
X grids = print = FALSE;
X getch();
X return;
X }
X
X /*
X * Skip data that has been optioned out.
X */
X
X if (coord.code < level)
X continue;
X
X if (coord.code > 5) { /* must be a header */
X
X point = TRUE;
X
X switch (coord.code / 1000) {
X case 1:
X if (boundaries) {
X if (colors)
X setcolor(EGA_LIGHTGRAY);
X break;
X }
X else
X continue;
X case 2:
X if (countries) {
X if (colors)
X setcolor(EGA_BROWN);
X break;
X }
X else
X continue;
X case 4:
X if (states) {
X if (colors)
X setcolor(EGA_BROWN);
X break;
X }
X else
X continue;
X case 5:
X if (islands) {
X if (colors)
X setcolor(EGA_LIGHTGRAY);
X break;
X }
X else
X continue;
X case 6:
X if (lakes) {
X if (colors)
X setcolor(EGA_BLUE);
X break;
X }
X else
X continue;
X case 7:
X if (rivers) {
X if (colors)
X setcolor(EGA_GREEN);
X break;
X }
X else
X continue;
X }
X }
X
X /* Convert database minutes of a degree to radians */
X
X lat = (float) coord.lat / 60.0F / RADIAN;
X lon = (float) coord.lon / 60.0F / RADIAN;
X
X if (compute(&lat, &lon, &x, &y)) {
X point = TRUE; /* offscale */
X continue;
X }
X
X if (point) {
X putpixel(x, y, getcolor());/* put down a dot */
X moveto(x, y);
X point = FALSE;
X }
X else
X lineto(x, y); /* connect the dots */
X }
X}
X
X/*
X * parse +-ddd.mm
X *
X * Change human degrees, and minutes to computer decimal.
X * Probably designed a monster for a simple solution here...
X */
X
Xfloat parse(string)
Xchar *string;
X{
X char *ptr, degrees[8], minutes[8];
X float num;
X
X strcpy(degrees, " "); /* pre-load with blanks */
X strcpy(minutes, " ");
X
X /* if no decimal point we assume a whole number */
X
X if ( (ptr = strchr(string, '.')) == (char *)NULL )
X return atof(string);
X
X /* else use the decimal point to offset */
X
X *ptr++ = '\0';
X
X strcpy(degrees, string);
X num = atof(degrees);
X
X switch (strlen(ptr)) {
X case 0:
X return atof(string);
X case 1:
X case 2:
X strcpy(minutes, ptr);
X break;
X default:
X return 361.0F; /* This will produce an error */
X }
X
X if (num >= 0.0F)
X num += atof(minutes) / 60.0F;
X else
X num -= atof(minutes) / 60.0F;
X
X return num;
X}
X
X
X/*
X * Draw grid lines from -180 to +180 Degrees (Longitude Lines),
X * as well as +80 to -80 Degrees (Lattitude Lines).
X */
X
Xvoid grid()
X{
X float lat, lon;
X int x, y, pass1;
X
X setcolor(grid_color);
X
X for (lon = -PI; lon <= PI; lon += TEN) {
X pass1 = TRUE;
X for (lat = EIGHTY; lat > -EIGHTY; lat -= TEN) {
X if (!compute(&lat, &lon, &x, &y)) {
X if (pass1) {
X putpixel(x, y, grid_color);
X moveto(x, y);
X pass1 = FALSE;
X } else
X lineto(x, y);
X } else
X pass1 = TRUE;
X }
X
X if (kbhit()) {
X print = FALSE;
X getch();
X return;
X }
X }
X
X for (lat = EIGHTY; lat > -EIGHTY; lat -= TEN) {
X pass1 = TRUE;
X for (lon = -PI; lon <= PI; lon += TEN) {
X if (!compute(&lat, &lon, &x, &y)) {
X if (pass1) {
X putpixel(x, y, grid_color);
X moveto(x, y);
X pass1 = FALSE;
X } else
X lineto(x, y);
X } else
X pass1 = TRUE;
X
X }
X
X if (kbhit()) {
X print = FALSE;
X getch();
X return;
X }
X }
X}
X
X/* EOF */
END_OF_mapper.c
if test 13949 -ne `wc -c <mapper.c`; then
echo shar: \"mapper.c\" unpacked with wrong size!
fi
# end of overwriting check
fi
echo shar: End of shell archive.
exit 0