home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / msdos / turbo_c / mapper17.c < prev    next >
Encoding:
C/C++ Source or Header  |  1988-12-06  |  13.6 KB  |  663 lines
  1. /*
  2.  *    mapper.c
  3.  *
  4.  *    Version 1.7 by Steven R. Sampson, November 1988
  5.  *
  6.  *    Based on a program and article by William D. Johnston
  7.  *    Copyright (c) May-June 1979 BYTE, All Rights Reserved
  8.  *
  9.  *    This program draws three types of map projections:
  10.  *    Perspective, Modified Perspective, and Azimuthal Equidistant.
  11.  *
  12.  *    Compiled with Turbo-C V1.5
  13.  */
  14.  
  15. #include <dos.h>
  16. #include <math.h>
  17. #include <stdio.h>
  18. #include <conio.h>
  19. #include <string.h>
  20. #include <stdlib.h>
  21. #include <graphics.h>
  22.  
  23. typedef int    bool;
  24.  
  25. /* Program Constants */
  26.  
  27. #define    FALSE    (bool) 0
  28. #define    TRUE    (bool) ~FALSE
  29.  
  30. #define    PI    (3.141593F)
  31. #define    HALFPI    (1.570796F)
  32. #define    TWOPI    (2.0F * PI)        /* Two Pi alias 360 Degrees         */
  33.  
  34. #define    RADIAN    (180.0F / PI )        /* One radian                 */
  35. #define    TWO    (2.0F / RADIAN)        /* 2 degrees in radians             */
  36. #define    TEN    (10.0F / RADIAN)    /* 10 degrees in radians         */
  37. #define    EIGHTY    (80.0F / RADIAN)    /* 80 degrees in radians         */
  38. #define    EARTH    (6378.0F)        /* Mean radius of earth (Kilometers) */
  39.  
  40. /* Program Globals */
  41.  
  42. FILE    *fp;
  43.  
  44. float    angle, maxplot, center_lat, center_lon, lat, lon, distance,
  45.     sin_of_distance, cos_of_distance, sin_of_center_lat, cos_of_center_lat,
  46.     scale, g, h2, facing_azimuth, aspect;
  47.  
  48. int    option, center_x, center_y, grid_color, level = 5;
  49. int    GraphDriver = DETECT, GraphMode;
  50.  
  51. char    optstring[] = "bcd:gilm:rsx?";
  52. char    database[128] = "mwdbii";    /* default name 'MWDBII'         */
  53.                     /* leave room for pathname!         */
  54. bool    boundaries = TRUE,        /* defaults to Boundaries, Islands   */
  55.     countries  = FALSE,
  56.     grids      = FALSE,
  57.     islands    = TRUE,
  58.     lakes      = FALSE,
  59.     rivers     = FALSE,
  60.     states     = FALSE,
  61.     colors     = FALSE;
  62.  
  63. /* Forward Declarations, Prototypes */
  64.  
  65. extern    int    getopt(int, char **, char *);
  66. extern    int    optind, opterr;
  67. extern    char    *optarg;
  68.  
  69. float    parse(char *);
  70. void    grid(void), plotmap(void), prompts(void), quit(void);
  71. bool    compute(float *, float *, int *, int *);
  72.  
  73.  
  74. main(argc, argv)
  75. int    argc;
  76. char    *argv[];
  77. {
  78.     register int    i;
  79.     int        err, xasp, yasp;
  80.  
  81.     registerbgidriver(EGAVGA_driver);
  82.     registerbgidriver(CGA_driver);
  83.  
  84.     setcbrk(TRUE);        /* Allow Control-C checking             */
  85.     ctrlbrk(quit);        /* Execute quit() if Control-C detected         */
  86.  
  87.     while ((i = getopt(argc, argv, optstring)) != -1)  {
  88.         switch (i)  {
  89.         case 'b':
  90.             boundaries = FALSE;
  91.             break;
  92.         case 'c':
  93.             countries = TRUE;
  94.             break;
  95.         case 'd':
  96.             strcpy(database, optarg);
  97.             break;
  98.         case 'g':
  99.             grids = TRUE;
  100.             break;
  101.         case 'i':
  102.             islands = FALSE;
  103.             break;
  104.         case 'l':
  105.             lakes = TRUE;
  106.             break;
  107.         case 'm':
  108.             level = atoi(optarg);
  109.             break;
  110.         case 'r':
  111.             rivers = TRUE;
  112.             break;
  113.         case 's':
  114.             states = TRUE;
  115.             break;
  116.         case 'x':
  117.             colors = FALSE;
  118.             break;
  119.         case '?':
  120.         default:
  121.               printf("Usage: mapper [/bcdgilmrsx]\n\n");
  122.               printf("  /b   Boundaries Off\n");
  123.               printf("  /c   Countries On\n");
  124.               printf("  /dn  Database ('MWDBII' Default)\n");
  125.               printf("  /g   Grid lines On\n");
  126.               printf("  /i   Islands Off\n");
  127.               printf("  /l   Lakes On\n");
  128.               printf("  /mn  Map Resolution (5 Default)\n");
  129.               printf("  /r   Rivers On\n");
  130.               printf("  /s   States On\n");
  131.               printf("  /x   Colors On\n\n");
  132.               printf("Defaults to Boundaries and Islands On\n");
  133.               exit(0);
  134.         }
  135.     }
  136.  
  137.     if ((fp = fopen(database, "rb")) == (FILE *)NULL)  {
  138.         printf("\007Error: Can't locate Database '%s'\n", database);
  139.         exit(1);
  140.     }
  141.  
  142.     initgraph(&GraphDriver, &GraphMode, "");/* initialize graphics         */
  143.     err = graphresult();
  144.  
  145.     restorecrtmode();            /* get back to text mode     */
  146.  
  147.     if (err < 0)  {
  148.         printf("Graphics Error - %s\n", grapherrormsg(err));
  149.         exit(-1);
  150.     }
  151.  
  152.     center_x = getmaxx() / 2;        /* get screen size for x, y  */
  153.     center_y = getmaxy() / 2;
  154.     getaspectratio(&xasp, &yasp);        /* squish factor for y axis  */
  155.     aspect = (float)xasp / (float)yasp;
  156.  
  157.     prompts();                /* get the basic map info    */
  158.     setgraphmode(GraphMode);        /*  and go to graphics mode  */
  159.  
  160.     if (GraphMode != CGAHI)  {
  161.         setbkcolor(BLACK);        /* must be EGA or VGA then   */
  162.         if (colors)
  163.             grid_color = EGA_LIGHTRED;
  164.         else
  165.             grid_color = EGA_LIGHTGRAY;
  166.     } else
  167.         grid_color = LIGHTGRAY;        /* CGA only has two colors   */
  168.  
  169.     setcolor(LIGHTGRAY);
  170.  
  171.     /*
  172.      *    See if data plotting is even needed
  173.      */
  174.  
  175.     if (boundaries || countries || islands || lakes || rivers || states)
  176.         plotmap();            /* display map on screen     */
  177.  
  178.     if (grids)
  179.         grid();                /* draw lat & long ref lines */
  180.  
  181.     if (print)
  182.         printscreen();            /* relay screen to printer   */
  183.  
  184.     sound(800);                /* 800 Hz for 1/4 a second   */
  185.     delay(125);
  186.     nosound();
  187.  
  188.     getch();                /* pause until key pressed   */
  189.     closegraph();                /* graphics off             */
  190.     fclose(fp);                /* close database file         */
  191.  
  192.     exit(0);
  193. }
  194.  
  195. /*
  196.  *    Return to operator following Control-C
  197.  */
  198.  
  199. void quit()
  200. {
  201.     closegraph();
  202.     fclose(fp);
  203.  
  204.     exit(0);
  205. }
  206.  
  207. /*
  208.  *    Operator prompts and input.
  209.  */
  210.  
  211. void prompts()
  212. {
  213.     char    temp[16];
  214.     float    ret, altitude;
  215.  
  216.     printf("West Longitudes and South Lattitudes are negative\n");
  217.  
  218.     /* input the world Lat & Long that is to be centered on */
  219.     /*   then convert the human notation to radians         */
  220.  
  221.     do  {
  222.         printf("\nLatitude of the map center [+-]dd.mm : ");
  223.         scanf("%s", temp);
  224.         ret = parse(temp);
  225.     } while (ret > 90.0F || ret < -90.0F);
  226.  
  227.     /* the arcosine function has trouble at 90 degrees */
  228.  
  229.     if (ret == 90.0F)
  230.         ret = 89.9F;
  231.  
  232.     if (ret == -90.0F)
  233.         ret = -89.9F;
  234.  
  235.     center_lat = ret / RADIAN;
  236.     sin_of_center_lat = sin(center_lat);
  237.     cos_of_center_lat = cos(center_lat);
  238.  
  239.     do  {
  240.         printf("Longitude of the map center [+-]ddd.mm : ");
  241.         scanf("%s", temp);
  242.         ret = parse(temp);
  243.     } while (ret > 180.0F || ret < -180.0F);
  244.  
  245.     center_lon = ret / RADIAN;
  246.  
  247.     do  {
  248.         printf("\nSelect from the following options:\n\n");
  249.         printf("  1 - Perspective Projection\n");
  250.         printf("  2 - Modified Perspective Projection\n");
  251.         printf("  3 - Azimuthal Equidistant Projection\n\n");
  252.         printf("Choice : ");
  253.         scanf("%d", &option);
  254.     } while (option < 1 || option > 3);
  255.  
  256.     if (option == 3)  {
  257.         maxplot = PI;        /* use HALFPI for less area        */
  258.         scale = (float)center_y / maxplot;
  259.         return;
  260.     }
  261.  
  262.     /* input the height above the terrain */
  263.  
  264.     printf("\nObserver altitude (km) : ");
  265.     scanf("%f", &altitude);
  266.  
  267.     h2 = EARTH + altitude;
  268.     maxplot = acos(EARTH / h2);
  269.  
  270.     /* the operator can orient the world upside down if they want */
  271.  
  272.     do  {
  273.         printf("Observer facing azimuth (0 - 359 degrees) : ");
  274.         scanf("%f", &facing_azimuth);
  275.     } while (facing_azimuth < 0.0F || facing_azimuth > 360.0F);
  276.  
  277.     facing_azimuth = -facing_azimuth / RADIAN;
  278.  
  279.     /* calculate the scale for the polar coordinates */
  280.  
  281.     scale = (float)center_y / (EARTH * sin(maxplot));
  282.  
  283.     /* for the perspective projection */
  284.  
  285.     g = EARTH * (h2 - EARTH * cos(maxplot));
  286. }
  287.  
  288.  
  289. /*
  290.  *    Convert the database to the desired projection by computation.
  291.  *
  292.  *    This code is a hand translation from BASIC to C based on Mr. Johnstons
  293.  *    code.  It is a non-mathematicians idea of what he meant.
  294.  *
  295.  *    Return TRUE if offscale else FALSE.
  296.  */
  297.  
  298. bool compute(lat, lon, x, y)
  299. register float    *lat, *lon;
  300. register int    *x, *y;
  301. {
  302.     float    sin_of_lat,
  303.         cos_of_lat,
  304.         abs_delta_lon,            /* absolute value         */
  305.         delta_lon,            /* x distance from center    */
  306.         delta_lat,            /* y distance from center    */
  307.         temp;                /* temporary storage         */
  308.  
  309.     /* normalize the longitude to +/- PI */
  310.  
  311.     delta_lon = *lon - center_lon;
  312.  
  313.     if (delta_lon < -PI)
  314.         delta_lon = delta_lon + TWOPI;
  315.  
  316.     if (delta_lon > PI)
  317.         delta_lon = delta_lon - TWOPI;
  318.  
  319.     abs_delta_lon = fabs(delta_lon);
  320.  
  321.     /*
  322.      *    If the delta_lon is within .00015 radians of 0 then
  323.      *    the difference is considered exactly zero.
  324.      *
  325.      *    This also simplifys the great circle bearing calculation.
  326.      */
  327.  
  328.     if (abs_delta_lon <= 0.00015F)  {
  329.         delta_lat = fabs(center_lat - *lat);
  330.  
  331.         if (delta_lat > maxplot)
  332.             return TRUE;        /* offscale             */
  333.  
  334.         if (*lat < center_lat)
  335.             angle = PI;
  336.         else
  337.             angle = 0.0F;
  338.  
  339.         sin_of_distance = sin(delta_lat);
  340.         cos_of_distance = cos(delta_lat);
  341.     }
  342.  
  343.     /*
  344.      *    Check if delta_lon is within .00015 radians of PI.
  345.      */
  346.  
  347.     else if (fabs(PI - abs_delta_lon) <= 0.00015F)  {
  348.         delta_lat = PI - center_lat - *lat;
  349.  
  350.         if (delta_lat > PI)  {
  351.             delta_lat = TWOPI - delta_lat;
  352.             angle = PI;
  353.         } else
  354.             angle = 0.0F;
  355.  
  356.         if (delta_lat > maxplot)
  357.             return TRUE;        /* offscale             */
  358.  
  359.         sin_of_distance = sin(delta_lat);
  360.         cos_of_distance = cos(delta_lat);
  361.     }
  362.  
  363.     /*
  364.      *    Simple calculations are out, now get cosmic.
  365.      */
  366.  
  367.     else  {
  368.         sin_of_lat = sin(*lat);
  369.         cos_of_lat = cos(*lat);
  370.  
  371.         cos_of_distance = sin_of_center_lat * sin_of_lat +
  372.                     cos_of_center_lat * cos_of_lat *
  373.                       cos(delta_lon);
  374.  
  375.         distance = acos(cos_of_distance);
  376.  
  377.         if (distance > maxplot)
  378.             return TRUE;        /* offscale             */
  379.  
  380.         sin_of_distance = sin(distance);
  381.  
  382.         temp = (sin_of_lat - sin_of_center_lat * cos_of_distance) /
  383.             (cos_of_center_lat * sin_of_distance);
  384.  
  385.         if (temp < -1.0F || temp > 1.0F)
  386.             return TRUE;        /* offscale             */
  387.  
  388.         angle = acos(temp);
  389.  
  390.         if (delta_lon < 0.0F)
  391.             angle = TWOPI - angle;
  392.     }
  393.  
  394.     if (facing_azimuth != 0.0F)  {
  395.         angle = angle - facing_azimuth;
  396.         if (angle < 0.0F)
  397.             angle = TWOPI + angle;
  398.     }
  399.  
  400.     angle = HALFPI - angle;
  401.  
  402.     if (angle < -PI)
  403.         angle = angle + TWOPI;
  404.  
  405.     switch (option)  {
  406.     case 1:
  407.         temp  = (scale * (g * sin_of_distance)) /
  408.                 (h2 - EARTH * cos_of_distance);
  409.         break;
  410.     case 2:
  411.         temp = scale * EARTH * sin_of_distance;
  412.         break;
  413.     case 3:
  414.         temp = scale * distance;
  415.     }
  416.  
  417.     /* convert polar to rectangular, correct for screen aspect */
  418.  
  419.     *x = center_x + (int)(temp * cos(angle));
  420.     *y = center_y - (int)(temp * sin(angle) * aspect);
  421.  
  422.     return FALSE;
  423. }
  424.  
  425. /*
  426.  *    Read the database and plot points or lines.
  427.  *
  428.  *    The database is Micro World Data Bank II.  It's based on the
  429.  *    CIA WDB-II tape available from NTIS.  Micro WDB-II was created
  430.  *    by Micro Doc.  Placed in the public domain by Fred Pospeschil
  431.  *    and Antonio Riveria.  Check on availability at:
  432.  *    1-402-291-0795  (6-9 PM Central)
  433.  *
  434.  *    Austin Code Works has something called: The World Digitized
  435.  *    that sounds like the same thing ($30.00), 1-512-258-0785
  436.  *
  437.  *    Lone Star Software has something called: The World Digitized
  438.  *    that sounds like the same thing ($6.00), 1-800-445-6172.
  439.  *
  440.  *    Database is in Intel word order:
  441.  *    code_lsb, code_msb, lat_lsb, lat_msb, lon_lsb, lon_msb
  442.  *
  443.  *    Code:    Integer, two meanings:
  444.  *        1.  Detail Level (1 Highest - 5 Lowest)
  445.  *
  446.  *        2.  Header (1xxx - 7xxx)    Command Line Options
  447.  *
  448.  *            1xxx    Boundaries        /b
  449.  *            2xxx    Countries        /c
  450.  *    (decimal)    4xxx    States            /s
  451.  *            5xxx    Islands            /i
  452.  *            6xxx    Lakes            /l
  453.  *            7xxx    Rivers            /r
  454.  *
  455.  *    Lat & Long:  Integer
  456.  *        Representing Minutes of degree
  457.  */
  458.  
  459. void plotmap()
  460. {
  461.     struct    { short code, lat, lon; } coord;
  462.     float    lat, lon;
  463.     int    x, y;
  464.     bool    point;
  465.  
  466.     point = TRUE;
  467.     while (fread(&coord, sizeof coord, 1, fp) > 0)  {
  468.  
  469.         if (kbhit())  {
  470.             grids = print = FALSE;
  471.             getch();
  472.             return;
  473.         }
  474.             
  475.         /*
  476.          *    Skip data that has been optioned out.
  477.          */
  478.  
  479.         if (coord.code < level)
  480.             continue;
  481.  
  482.         if (coord.code > 5)  {        /* must be a header         */
  483.  
  484.             point = TRUE;
  485.  
  486.             switch (coord.code / 1000)  {
  487.             case 1:
  488.                 if (boundaries)  {
  489.                     if (colors)
  490.                         setcolor(EGA_LIGHTGRAY);
  491.                     break;
  492.                 }
  493.                 else
  494.                     continue;
  495.             case 2:
  496.                 if (countries)  {
  497.                     if (colors)
  498.                         setcolor(EGA_BROWN);
  499.                     break;
  500.                 }
  501.                 else
  502.                     continue;
  503.             case 4:
  504.                 if (states)  {
  505.                     if (colors)
  506.                         setcolor(EGA_BROWN);
  507.                     break;
  508.                 }
  509.                 else
  510.                     continue;
  511.             case 5:
  512.                 if (islands)  {
  513.                     if (colors)
  514.                         setcolor(EGA_LIGHTGRAY);
  515.                     break;
  516.                 }
  517.                 else
  518.                     continue;
  519.             case 6:
  520.                 if (lakes)  {
  521.                     if (colors)
  522.                         setcolor(EGA_BLUE);
  523.                     break;
  524.                 }
  525.                 else
  526.                     continue;
  527.             case 7:
  528.                 if (rivers)  {
  529.                     if (colors)
  530.                         setcolor(EGA_GREEN);
  531.                     break;
  532.                 }
  533.                 else
  534.                     continue;
  535.             }
  536.         }
  537.  
  538.         /*  Convert database minutes of a degree to radians */
  539.  
  540.         lat =  (float) coord.lat / 60.0F / RADIAN;
  541.         lon =  (float) coord.lon / 60.0F / RADIAN;
  542.  
  543.         if (compute(&lat, &lon, &x, &y))  {
  544.             point = TRUE;        /* offscale             */
  545.             continue;
  546.         }
  547.  
  548.         if (point)  {
  549.             putpixel(x, y, getcolor());/* put down a dot         */
  550.             moveto(x, y);
  551.             point = FALSE;
  552.         }
  553.         else
  554.             lineto(x, y);        /* connect the dots         */
  555.     }
  556. }
  557.  
  558. /*
  559.  *    parse +-ddd.mm
  560.  *
  561.  *    Change human degrees, and minutes to computer decimal.
  562.  *    Probably designed a monster for a simple solution here...
  563.  */
  564.  
  565. float parse(string)
  566. char    *string;
  567. {
  568.     char    *ptr, degrees[8], minutes[8];
  569.     float    num;
  570.  
  571.     strcpy(degrees, "       ");        /* pre-load with blanks      */
  572.     strcpy(minutes, "       ");
  573.  
  574.     /* if no decimal point we assume a whole number */
  575.  
  576.     if ( (ptr = strchr(string, '.')) == (char *)NULL )
  577.         return atof(string);
  578.  
  579.     /* else use the decimal point to offset */
  580.  
  581.     *ptr++ = '\0';
  582.  
  583.     strcpy(degrees, string);
  584.     num = atof(degrees);
  585.  
  586.     switch (strlen(ptr))  {
  587.     case 0:
  588.         return atof(string);
  589.     case 1:
  590.     case 2:
  591.         strcpy(minutes, ptr);
  592.         break;
  593.     default:
  594.         return 361.0F;    /* This will produce an error             */
  595.     }
  596.  
  597.     if (num >= 0.0F)
  598.         num += atof(minutes) / 60.0F;
  599.     else
  600.         num -= atof(minutes) / 60.0F;
  601.  
  602.     return num;
  603. }
  604.  
  605.  
  606. /*
  607.  *    Draw grid lines from -180 to +180 Degrees (Longitude Lines),
  608.  *    as well as +80 to -80 Degrees (Lattitude Lines).
  609.  */
  610.  
  611. void grid()
  612. {
  613.     float    lat, lon;
  614.     int    x, y, pass1;
  615.  
  616.     setcolor(grid_color);
  617.  
  618.     for (lon = -PI; lon <= PI; lon += TEN)  {
  619.         pass1 = TRUE;
  620.         for (lat = EIGHTY; lat > -EIGHTY; lat -= TEN)  {
  621.             if (!compute(&lat, &lon, &x, &y))  {
  622.                 if (pass1)  {
  623.                     putpixel(x, y, grid_color);
  624.                     moveto(x, y);
  625.                     pass1 = FALSE;
  626.                 } else
  627.                     lineto(x, y);
  628.             } else
  629.                 pass1 = TRUE;
  630.         }
  631.  
  632.         if (kbhit())  {
  633.             print = FALSE;
  634.             getch();
  635.             return;
  636.         }
  637.     }
  638.  
  639.     for (lat = EIGHTY; lat > -EIGHTY; lat -= TEN)  {
  640.         pass1 = TRUE;
  641.         for (lon = -PI; lon <= PI; lon += TEN)  {
  642.             if (!compute(&lat, &lon, &x, &y))  {
  643.                 if (pass1)  {
  644.                     putpixel(x, y, grid_color);
  645.                     moveto(x, y);
  646.                     pass1 = FALSE;
  647.                 } else
  648.                     lineto(x, y);
  649.             } else
  650.                 pass1 = TRUE;
  651.  
  652.         }
  653.  
  654.         if (kbhit())  {
  655.             print = FALSE;
  656.             getch();
  657.             return;
  658.         }
  659.     }
  660. }
  661.  
  662. /* EOF */
  663.