ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng)/ProfitPress-MegaCDROM2.B6I/GRAPHICS/MISC/WORLD10.ZIP/WORLDMAP.TXT

Text File | 1987-12-06 | 30.5 KB | 988 lines

==============================

          READ.ME

==============================


To use Cartog.pas with an EGA or Hercules board, you will need to use
the following constants.


EGA Constants:                Hercules Constants:

CONST XCENTER = 320;          CONST XCENTER = 360;
      YCENTER = 174;                YCENTER = 174;
      ASPECT  = 1.37;               ASPECT  = 1.5;
      R       = 70;                 R       = 70;

Also you will need a software driver for these boards.  For example, if
you want to use Turbo Graphics Toolbox, you will need to make the
following changes: The following three lines go at the top of the
program, right after the line "Program Cartog;"

{$I Typedef.sys}
{$I graphix.sys}
{$I Kernel.sys}

replace "HiRes;  HiResColor(15);" with the following four lines:

Initgraphic;
defineWorld(1,0,348,719,0);
selectWorld(1);
selectWindow(1);

Finally replace  TEXTMODE(BW80); with  Leavegraphic;


---------
Setup for the plotter:

Our plots were made with a Hewlett-Packard 7475. A separate
version of CARTOG was made to output to the plotter. The
constants that we used were:

CONST XCENTER = 6.2; { Page center for B size plot }
      YCENTER = 5.1;
      ASPECT  = 1;  { Scale factors are the same in both axes    }
      R       = 1;  { Unit radius, scale in the plotting routine }

We also added a plotter initialization procedure:

PROCEDURE Initialize;
{ PLT is a text file containing the plotter commands.
  At the termination of this program, the file is copied to the
  plotter
}
BEGIN
   WRITE(' ':13, 'Enter map size (A/B): '); READLN(SIZE);
   SIZE:=UPCASE(SIZE);

   ASSIGN(PLT, 'WORLD.PLT');  REWRITE(PLT);

   WRITELN(PLT, 'IN;');  { Command to initialize plotter }

             { Set paper size and scaling points, so
               plotter effectively works in inches
               for B size plots and a little less for
               A size plots }
   IF SIZE = 'A' THEN
    WRITELN(PLT,'PS4; IP 100,100,868,868;')
   ELSE
    WRITELN(PLT, 'PS0; IP 100,100,1124,1124;');
   WRITELN(PLT, 'SC 0,1,0,1;');
            { Set pen speed (slow) and select pen 3 }
   WRITELN(PLT, 'VS 4; SP 3;');

END;                    { Initialize. }

Finally, we substituted the following for PROCEDURE PlotPt:

PROCEDURE PlotPt(VAR LastPtVis: BOOLEAN);
{  Draws a line from the last point to the current (XP,YP) if it
   is visible. }
BEGIN
    IF LastPtVis THEN           { Pen down, go to XP, YP }
     WRITELN(PLT, 'PD', XP:7:3, ',', YP:7:3,';')
    ELSE                        { Pen up,   go to XP, YP }
     WRITELN(PLT, 'PU' , XP:7:3, ',', YP:7:3,';');

    LastX:=; LastY:=YP;
        LastPtVis:=TRUE;
END; { PlotPt. }




==============================

          CARTOG.PAS

==============================


_______________________________________________________________________

CARTOG.PAS Accompanies "Mapping the World in Pascal" by Robert Miller and
Francis Reddy, BYTE, December 1987, page 329
_______________________________________________________________________

PROGRAM Cartog;
{ This program plots geographic data from the file
  WORLD.DAT and coordinate grids on the Mercator,
  Equidistant Cylindrical, Sinusoidal, Hammer, and
  Orthographic map projections.
}
CONST Sqrt2        = 1.4142135623731;
      PI           = 3.1415926535898;
      HalfPI       = 1.5707963267949;
      TwoPI        = 6.2831853071796;
      Radian       = 1.7453292519943E-2;
      RadianDiv100 = 1.7453292519943E-4; { PI/180/100, needed to convert }
                                         { data in WORLD.DAT to radians  }

CONST XCENTER    : INTEGER = 320;        { CGA Graphics constants.       }
      YCENTER    : INTEGER =  99;        { Screen center X and Y         }
      ASPECT     : REAL    = 2.4;        { 640x200 aspect ratio          }
      R          : REAL    =  40;        { Default map radius            }
      NotVisible : INTEGER = -32767;     { Flag for point visibility     }

TYPE  LLREC      = RECORD
      CODE       : ARRAY[0..1] OF CHAR;
      LONGI, LATI: INTEGER; END;

VAR   LL         : LLREC;
      LLF        : FILE OF LLREC;

VAR   LastX, LastY, XP, YP    : INTEGER; { Save variables for plotting }
      COLOR_GLB               : INTEGER;

VAR   I, J, K, MapType, M, X1,Y1,
      X2, Y2, SX, SY, CENTER  : INTEGER;

VAR   L, L1, LONGR, LSTEP,
      B, LATR, BSTEP, X, Y,
      PHI1, Lambda0           : REAL;

VAR   XX, YY, SA, SB          : REAL;

VAR   LastPtVis, GRID         : BOOLEAN;
VAR   CH                      : CHAR;

FUNCTION ArcCos(X: REAL): REAL;
BEGIN
    IF  ABS(X) < 1  THEN ArcCos:= ARCTAN(SQRT(1-SQR(X))/X)
    ELSE  IF X = 1  THEN ArcCos:=  0
    ELSE  IF X =-1  THEN ArcCos:= PI;
END;   { ArcCos. }

FUNCTION ArcSin(X: REAL): REAL;
BEGIN
    IF ABS(X) < 1  THEN ArcSin:= ARCTAN(X/SQRT(1-SQR(X)))
    ELSE IF X = 1  THEN ArcSin:= HalfPI
    ELSE IF X =-1  THEN ArcSin:=-HalfPI;
END;    { ArcSin. }

FUNCTION ArcTanH(X : Real):  Real;
VAR A,T :  REAL;
BEGIN
    T:=ABS(X);
    IF T < 1 THEN
    BEGIN
       A := 0.5 * LN((1 + T)/(1 - T));
       IF X < 0 THEN ArcTanH := -A ELSE ArcTanH :=A;
    END;
END;  { ArcTanH. }

FUNCTION Meridian(Lambda, Lambda0: REAL):REAL;
{ Returns difference between current longitude and map center. }
VAR DelLam : REAL;
BEGIN
    DelLam := Lambda - Lambda0;
    IF DelLam < -PI THEN DelLam := DelLam + TwoPI
    ELSE
    IF DelLam >  PI THEN DelLam := DelLam - TwoPI;
    Meridian:=DelLam;
END;   { Meridian. }

PROCEDURE Mercator(Lambda, Lambda0, Phi, R : REAL; VAR X, Y : REAL);
{ For R = 1: -Pi <= X <= Pi, -Pi/2 <= Y <= Pi/2.  }
CONST MaxLat : REAL = 1.397;   {~80 degrees. }
     { REAL = 1.483;    ~85 degrees. }
BEGIN
    IF ABS(Phi) < MaxLat THEN
    BEGIN
       Lambda := Meridian(Lambda, Lambda0);
       X := R * Lambda;
       Y := R * ArcTanH(SIN(Phi));
    END
    ELSE X := NotVisible;
END;   { Mercator. }

PROCEDURE EquiCyl(Lambda, Lambda0, Phi, Phi1, R : REAL; VAR X, Y : REAL);
{ For R = 1: -Pi <= X <= Pi, -Pi/2 <= Y <= Pi/2. }
BEGIN
    Lambda := Meridian(Lambda, Lambda0);
    X := R * Lambda * COS(Phi1);
    Y := R * Phi;
END;   { EquiCyl. }

PROCEDURE Sinusoidal(Lambda, Lambda0, Phi, R : REAL; VAR X, Y : REAL);
{ For R = 1:  -Pi <= X <= Pi  and  -Pi/2 <= Y <= Pi/2. }
BEGIN
     Lambda := Meridian(Lambda, Lambda0);
     X := R * Cos(Phi) * Lambda ;
     Y := R * Phi;
END;  { Sinusoidal. }

PROCEDURE Hammer(Lambda, Lambda0, Phi, R : REAL; VAR X, Y : REAL);
{ For R = 1: -2√2 <= X <=2√2  and  - √2 <= Y <= √2. }
VAR K, CosPhi, HalfLambda        : REAL;
BEGIN
    HalfLambda := 0.5*Meridian(Lambda, Lambda0);
    CosPhi:=COS(Phi);
    K := R * SQRT2 / SQRT(1 +CosPhi * COS(HalfLambda));
    X := 2 * K * CosPhi * (SIN(HalfLambda));
    Y := K * SIN(Phi);
END;  { Hammer. }

PROCEDURE Orthographic(Lambda, Lambda0, Phi, Phi1, R: REAL; VAR X, Y : REAL);
{ For R = 1: -2 <= X,Y <= 2. }
VAR CosC, CosL, SinPhi1, CosPhi1, SinPhi, CosPhi, R2 :  Real;
BEGIN
    Lambda :=Meridian(Lambda, Lambda0);  R2:=R+R;
    CosPhi1:=COS(Phi1);   SinPhi1:=SIN(Phi1);
    CosPhi :=COS(Phi);    SinPhi:= SIN(Phi);
    CosL   :=COS(Lambda)*CosPhi;
    CosC   :=SinPhi1 * SinPhi + CosPhi1 * COSL;
    IF CosC >= 0 THEN
    BEGIN
       X :=R2 * CosPhi * SIN(Lambda);
       Y :=R2 * (CosPhi1 * SinPhi - SinPhi1 * COSL);
    END ELSE X:=NotVisible;
END;  { Orthographic. }

PROCEDURE Beep;
{ Sounds a tone when map is complete. }
BEGIN
    Sound(880);  Delay(250);  NoSound;
END;

PROCEDURE PlotPt(VAR LastPtVis: BOOLEAN);
{ Draws a line from the last point to the current (XP,YP) if it is visible. }
VAR IX,IY: INTEGER;
LABEL XIT;
BEGIN
    IX:=ROUND(XP); IY:=ROUND(YP);
    IF LastPtVis THEN DRAW(LastX,LastY,IX,IY,1);
    LastX:=IX; LastY:=IY;
    LastPtVis:=TRUE;
  XIT:
END; { PlotPt. }

PROCEDURE CoordinateGrid(OUTLINE: BOOLEAN; MapType: INTEGER);
CONST LatitudeSpacing  = 30;
      LongitudeSpacing = 30;

VAR Longitude, Latitude, LatLimit,
    MaxLat, LongIncr, LatIncr      : INTEGER;
VAR LL, PP, A, R2, RA, XN, YN,
    SINDT, COSDT                   : REAL;
BEGIN
    CASE MapType OF
       1: BEGIN   MaxLat:=80;  LongIncr:=360;  LatIncr:=160;    END;
       2: BEGIN   MaxLat:=90;  LongIncr:=360;  LatIncr:=180;    END;
       3: BEGIN   MaxLat:=90;  LongIncr:=360;  LatIncr:=5;      END;
    4..5: BEGIN   Maxlat:=90;  LongIncr:=5;    LatIncr:=5;      END;
    END;  { CASE...}

    LL:=0;  PP:=Phi1;
    IF OUTLINE THEN
    BEGIN
        IF MapType = 5 THEN PP:=0;
        LatLimit:=MaxLat;              { Draw only extreme latitudes }
                                       { to make map outline         }
    END
    ELSE LatLimit:= MaxLat DIV LatitudeSpacing*LatitudeSpacing;

    Latitude:=LatLimit;

    WHILE Latitude >= -LatLimit DO      { Draw parallels }
    BEGIN
       LATR:=Latitude*Radian;
       LastPtVis:=FALSE;

       Longitude:=-180;
       WHILE Longitude <= 180 DO
       BEGIN
          LONGR:=Longitude*Radian;

          CASE MapType OF
          1: BEGIN MERCATOR(LONGR, LL, LATR, R, X, Y);      END;
          2: BEGIN EQUICYL(LONGR, LL, LATR, PP, R, X, Y);   END;
          3: BEGIN SINUSOIDAL(LONGR, LL, LATR, R, X, Y);    END;
          4: BEGIN HAMMER(LONGR, LL, LATR, R, X, Y);        END;
          5: BEGIN ORTHOGRAPHIC (LONGR, LL, LATR, PP, R, X, Y); END;
          END;  { CASE...}

          IF X > -300 THEN
          BEGIN
             XP:=ROUND(X*ASPECT)+XCENTER;
             YP:=YCENTER-ROUND(Y);
             PlotPt(LastPtVis);
          END ELSE LastPtVis:=FALSE;

          Longitude:=Longitude+LongIncr;
       END;

       IF OUTLINE THEN
            Latitude:=Latitude-2*MaxLat
       ELSE
            Latitude:=Latitude-LatitudeSpacing;
       END;

       IF OUTLINE THEN LL:=0 ELSE LL:=Lambda0;

       Longitude:=-180;                   { Draw meridians }

       IF MapType >= 4 THEN MaxLat:=90;
       WHILE Longitude <= 180 DO
       BEGIN
          LONGR:=Longitude*Radian;
          LastPtVis:=FALSE;
          Latitude:=MaxLat;
          WHILE Latitude >= -MaxLat DO
          BEGIN
             LATR:=Latitude*Radian;

             CASE MapType OF
             1: BEGIN MERCATOR(LONGR, LL, LATR, R, X, Y);      END;
             2: BEGIN EQUICYL(LONGR, LL, LATR, PP, R, X, Y);   END;
             3: BEGIN SINUSOIDAL(LONGR, LL, LATR, R, X, Y);    END;
             4: BEGIN HAMMER(LONGR, LL, LATR, R, X, Y);        END;
             5: BEGIN ORTHOGRAPHIC( LONGR, LL, LATR, PP, R, X, Y); END;
             END;  { CASE...}

          IF X > -300 THEN
          BEGIN
             XP:=ROUND(X*ASPECT)+XCENTER;
             YP:=YCENTER-ROUND(Y);
             PlotPt(LastPtVis);
          END ELSE LastPtVis:=FALSE;

          Latitude:=Latitude-LatIncr;
     END;

     IF OUTLINE THEN
          Longitude:=Longitude+360
     ELSE
          Longitude:=Longitude+LongitudeSpacing;
 END;

 IF OUTLINE AND (MapType=5) THEN
 BEGIN
    A:=0;                           { Draw circular outline }
    LastPtVis:=False;
    R2:=R + R;
    RA:= R2 * Aspect;
    SINDT:= 0.05996400648;
    COSDT:= 0.99820053993;
    X:=1;   Y:=0;
    XP:= ROUND(XCENTER + RA);
    YP:= ROUND(YCENTER);
    PlotPt(LastPtVis);
    WHILE A <= TwoPI DO
    BEGIN                           { Compute points on the circle }
       XN:= X * COSDT - Y * SINDT;
       YN:= X * SINDT + Y * COSDT;
       X:= XN;  Y:= YN;
       XP:= XCENTER + ROUND(X*RA);
       YP:= YCENTER + ROUND(Y*R2);
       PlotPt(LastPtVis);
       A:= A+0.06;
    END; { While. }
  END;
END;  { CoordinateGrid. }

PROCEDURE DrawMap(MapType: INTEGER);
VAR Latitude, Longitude : REAL;
VAR LastX : INTEGER;
LABEL XIT;
BEGIN
    LastPtVis:=FALSE; LastX:=0;

    ASSIGN(LLF, 'WORLD.DAT');  RESET(LLF);
    WHILE NOT EOF(LLF) DO
    BEGIN
       READ(LLF, LL);
       IF KeyPressed THEN GOTO XIT;
       LONGR:=LL.LONGI * RadianDiv100;
       LATR :=LL.LATI  * RadianDiv100;

       IF LL.CODE = 'LS' THEN LastPtVis:=FALSE;
       IF (LL.CODE = 'S ') OR (LL.CODE = 'LS') THEN
       BEGIN
          CASE MapType OF
          1: BEGIN MERCATOR(LONGR, Lambda0, LATR, R, X, Y);      END;
          2: BEGIN EQUICYL(LONGR, Lambda0, LATR, Phi1, R, X, Y); END;
          3: BEGIN SINUSOIDAL(LONGR, Lambda0, LATR, R, X, Y);    END;
          4: BEGIN HAMMER(LONGR, Lambda0, LATR, R, X, Y);        END;
          5: BEGIN ORTHOGRAPHIC(LONGR, Lambda0, LATR, Phi1, R, X, Y);  END;
          END;  { CASE...}

          IF X > -300 THEN
          BEGIN
             XP:=ROUND(X*ASPECT)+XCENTER;
             IF ABS(LastX-XP) > 100 THEN LastPtVis:=FALSE;
             YP:= YCENTER-ROUND(Y);
             PlotPt(LastPtVis);  LastX:=XP;
          END ELSE LastPtVis:=FALSE;
       END;
    END;
 XIT:
END;  { DrawMap. }

(* ---------------------  MAIN  PROGRAM  ------------------ *)

VAR RESP : CHAR;
LABEL XIT;
BEGIN
    MapType:=1;
    WHILE MapType > 0 DO      (*       MENU       *)
    BEGIN
        ClrScr;

        GOTOXY(24,1);  WRITE('C A R T O G');
        LowVideo;
        GOTOXY(1,24);
        WRITE('':4,'Copyright 1987 by Robert Miller and Francis Reddy');
        GOTOXY(1,3);
        WRITELN(' ':4,'To PLOT:  Choose a projection.  Enter the Central ');
        WRITELN(' ':4,'Meridian of the map (180 to -180 degrees, longitudes');
        WRITELN(' ':4,'west of Greenwich negative). If applicable, enter');
        WRITELN(' ':4,'the Standard Parallel (90 to -90 degrees, southern');
        WRITELN(' ':4,'latitudes negative). The file WORLD.DAT must also be');
        WRITELN(' ':4,'on the logged drive. A tone means the map is done.');
        WRITELN;
        WRITELN(' ':4,'Any key ABORTS plot.  Hit return to restore MENU.');
        NormVideo;
        WRITELN;
        WRITELN;
        WRITE(' ':6,'1. Mercator');
        WRITELN(' ':21,'4. Hammer');
        WRITE(' ':6,'2. Equidistant Cylindrical');
        WRITELN(' ':6,'5. Orthographic');
        WRITELN(' ':6,'3. Sinusoidal');
        WRITELN;
        WRITE(' ':8,'Projection number (1-5) or 0 to quit:   ');
        READLN(MapType);
        If MapType = 0 THEN GOTO XIT;

        WRITELN;
        WRITE(' ':8,'Central Longitude of Map (default = 0): ');
        Lambda0:=0;
        READLN(Lambda0);  Lambda0:=Lambda0*Radian;

        IF (MapType = 2) OR (MapType = 5) THEN
        BEGIN
            WRITE(' ':8,'Central Latitude  of Map (default = 0): ');
            Phi1:=0;  READLN(Phi1);
            IF Phi1 = 90 THEN Phi1 := HalfPI
            ELSE
            Phi1:=Phi1*Radian;
        END;

        IF MapType >= 4 THEN R:=83 ELSE R:=70;
        WRITE(' ':8,'Plot grid, continents or both (G/C/B)?  ');
        READLN(RESP);  RESP:=UPCASE(RESP);
        GRID:=(RESP ='G') OR (RESP = 'B');

        HiRes;  HiResColor(15);            { Set CGA Graphics Mode }

        IF GRID THEN CoordinateGrid(FALSE,  MapType);
        CoordinateGrid(TRUE, MapType);

        IF (RESP = 'B') OR (RESP = 'C') THEN DrawMap(MapType);
      Beep;
    XIT:

    IF MapType > 0 THEN
       While NOT KeyPressed DO ;           { Wait for key strike }

    TEXTMODE(BW80);                        { Return to Text Mode }
    ClrScr;
    END;  { WHILE MapType > 0...}
END.



==============================

          CPLOT.PAS

==============================



_______________________________________________________________________

CPLOT.PAS Accompanies "Mapping the World in Pascal" by Robert Miller and Francis Reddy, BYTE, December 1987, page 329
_______________________________________________________________________

PROGRAM MapProjections;
{  Map projection program for Hewlett-Packard 7475 plotter. }
{$V-}
CONST Sqrt2   = 1.4142135623710;
      PI      = 3.141592653589793238;
      HalfPI  = 1.570796326794897;
      ThreePI2= 4.71238898038469;
      TwoPI   = 6.283185307179587;
      Degree  = 57.29577951308232088;
      Radian  = 0.01745329251994329577;

CONST XCENTER : REAL    = 6.2;      { CGA Graphics constants. }
      YCENTER : REAL    = 5.1;
      ASPECT  : REAL    =   1;
      R       : REAL    =   1;

TYPE  S80       = STRING[80];
TYPE  LLREC     = RECORD
      CODE        : ARRAY[0..1] OF CHAR;
      LONGI, LATI : INTEGER; END;

VAR   LL         : LLREC;
      LLF        : FILE OF LLREC;
      PLT        : TEXT[$1000];
VAR   FNAME      : STRING[64];

VAR   LastX, LastY, XP, YP  : REAL;  { Save variables for plotting. }
      COLOR_GLB, IPEN       : INTEGER;
      SIZE                  : CHAR;

VAR   I, J, K, MapType, M, X1,Y1, X2, Y2,
      SX, SY, CENTER        : INTEGER;

VAR   L, L1, LONGR, LSTEP,
      B, LATR, BSTEP, X, Y,
      PHI1, Lambda0         : REAL;

VAR XX, YY, SA, SB          : REAL;

VAR   LastPtVis, GRID       : BOOLEAN;
VAR   CH                    : CHAR;

PROCEDURE Initialize;
BEGIN
    WRITE(' ':13, 'Enter map size (A/B): '); READLN(SIZE);
    SIZE:=UPCASE(SIZE);;

    ASSIGN(PLT, 'WORLD.PLT');  REWRITE(PLT);
    WRITELN(PLT, 'IN;');
    IF SIZE = 'A' THEN WRITELN(PLT,'PS4; IP 100,100,868,868;')
    ELSE  WRITELN(PLT, 'PS0; IP 100,100,1124,1124;');
    WRITELN(PLT, 'SC 0,1,0,1;');
    WRITELN(PLT, 'VS 4; SP 3;');
END;  { Initialize. }

PROCEDURE LabelMap(MapType : INTEGER);
VAR TITLE  : STRING[80];
VAR XT, YT : REAL;
BEGIN
     CASE MapType OF
     1: BEGIN TITLE:='Mercator';                  END;
     2: BEGIN TITLE:='Equidistant Cylindrical';   END;
     3: BEGIN TITLE:='Miller';                    END;
     4: BEGIN TITLE:='Sinusoidal';                END;
     5: BEGIN TITLE:='Hammer';                    END;
     6: BEGIN TITLE:='Orthographic' ;             END;
     7: BEGIN TITLE:='Stereographic';             END;
     END;

     TITLE:=TITLE + ' Map Projection';
     IF SIZE = 'A' THEN XT:=12.38
     ELSE XT:= 14;
     YT:= YCENTER-4.5;
     WRITELN(PLT, 'DI 0,1; SR 18,18; PU',XT:7:2, ',', YT:7:2);
     WRITELN(PLT, 'LB' + TITLE + CHR(3)+';');

     XT:=XT +0.24;
     WRITELN(PLT, 'PU ',XT:7:2, ',', YT:7:2, '; SR 14,14;');
     WRITELN(PLT, 'LBCentral Meridian ',ROUND(Lambda0*Degree):4,
                  ' degrees'+CHR(3)+';');

     IF (MapType = 2) OR (MapType = 6) OR (MapType = 7) THEN
     BEGIN
         XT:=XT +0.24;
         WRITELN(PLT, 'PU ',XT:7:2, ',', YT:7:2, ';');
         WRITELN(PLT, 'LBCentral Latitude ',ROUND(Phi1*Degree):4,
                      ' degrees'+CHR(3)+';');
     END;

     XT:=XT +0.24;
     WRITELN(PLT, 'SR 10,10; PU', XT:7:2 ,',', YT:7:2,';');
     WRITELN(PLT, 'LBRobert D. Miller, 1986.' +CHR(3)+';');
END;  { LabelMap. }

FUNCTION ArcCos(X: REAL): REAL;
BEGIN
    IF  ABS(X) < 1  THEN ArcCos:= ARCTAN(SQRT(1-SQR(X))/X)
    ELSE  IF X = 1  THEN ArcCos:=  0
    ELSE  IF X =-1  THEN ArcCos:= PI;
END;   { ArcCos. }

FUNCTION ArcSin(X: REAL): REAL;
BEGIN
    IF ABS(X) < 1  THEN ArcSin:= ARCTAN(X/SQRT(1-SQR(X)))
    ELSE IF X = 1  THEN ArcSin:= HalfPI
    ELSE IF X =-1  THEN ArcSin:=-HalfPI;
END;    { ArcSin. }

FUNCTION ArcTanH(TERM : Real):  Real;   (* Inverse hyperbolic tangent *)
VAR A,T :  real;
BEGIN
    T:=ABS(TERM);
    IF T < 1 THEN
    BEGIN
        A := 0.5 * LN((1 +T)/(1 -T));
        IF TERM < 0 THEN ArcTanH := -A ELSE ArcTanH :=A;
    END;
END;  { ArcTanH. }

(*
        Map Projection Library
        by Robert Miller and Francis Reddy, 20 October 1986.
        The following routines are based on equations found
        in the U.S. Geological Survey's Bulletin 1532,
        "Map Projections Used by the U.S. Geological Survey"
        by John P. Snyder and "Introduction to Map Projections"
        by Porter McDonnell, Jr.
*)
FUNCTION Meridian(Lambda, Lambda0: REAL):REAL;
{ Returns difference between current longitude and map center. }
VAR DelLam : REAL;
BEGIN
    DelLam := Lambda - Lambda0;
    IF DelLam < -PI THEN DelLam := DelLam + TwoPi
    ELSE
    IF DelLam >  PI THEN DelLam := DelLam - TwoPi;
    Meridian:=DelLam;
END;   { Meridian. }

PROCEDURE Mercator(Lambda, Lambda0, Phi, R : REAL; VAR X, Y : REAL);
{  For R = 1: -Pi <= X <= Pi, -Pi/2 <= Y <= Pi/2.  }
CONST MaxLat : REAL = 1.397;   {~80 degrees. }
             { REAL = 1.483;    ~85 degrees. }
BEGIN
    IF ABS(Phi) < MaxLat THEN
    BEGIN
        Lambda := Meridian(Lambda, Lambda0);
        X := R * Lambda;
        Y := R * ArcTanH(SIN(Phi));
    END
    ELSE X := -32767;
END;   { Mercator. }

PROCEDURE EquiCyl(Lambda, Lambda0, Phi, Phi1, R : REAL; VAR X, Y : REAL);
{  For R = 1: -Pi <= X <= Pi, -Pi/2 <= Y <= Pi/2. }
BEGIN
    Lambda := Meridian(Lambda, Lambda0);
    X := R * Lambda * COS(Phi1);
    Y := R * Phi;
END;   { EquiCyl. }

PROCEDURE Miller(Lambda, Lambda0, Phi, R : REAL; VAR X, Y : REAL);
{ For R = 1: -Pi <= X <= Pi, -Pi/2 <= Y <= Pi/2. }
BEGIN
     Lambda := Meridian(Lambda, Lambda0);
     X := R * Lambda;
     Y := R * ArcTanH(SIN(0.8 * Phi)) * 1.25;
END;   { Miller. }

PROCEDURE Sinusoidal(Lambda, Lambda0, Phi, R : REAL; VAR X, Y : REAL);
{ For R = 1:  -Pi <= X <= Pi  and  -Pi/2 <= Y <= Pi/2. }
BEGIN
     Lambda := Meridian(Lambda, Lambda0);
     X := R * Cos(Phi) * Lambda ;
     Y := R * Phi;
END;  { Sinusoidal. }

PROCEDURE Hammer(Lambda, Lambda0, Phi, R : REAL; VAR X, Y : REAL);
{  For R = 1: -2√2 <= X <=2√2  and  -√2 <= Y <= √2. }
VAR K, CosPhi, HalfLambda        : REAL;
BEGIN
    HalfLambda := 0.5*Meridian(Lambda, Lambda0);
    CosPhi:=COS(Phi);
    K := R * SQRT2 / SQRT(1 +CosPhi * COS(HalfLambda));
    X := 2 * K * CosPhi * (SIN(HalfLambda));
    Y := K * SIN(Phi);
END;  { Hammer. }

PROCEDURE Orthographic(Lambda, Lambda0, Phi, Phi1, R: REAL; VAR X, Y : REAL);
{  For unit radius R of generating globe, R = 1,  -2 <= X,Y <= 2. }
VAR CosC, CosL, SinPhi1, CosPhi1, SinPhi, CosPhi, R2 :  Real;
BEGIN
    Lambda :=Meridian(Lambda, Lambda0);  R2:=R+R;
    CosPhi1:=COS(Phi1);   SinPhi1:=SIN(Phi1);
    CosPhi :=COS(Phi);    SinPhi:= SIN(Phi);
    CosL   :=COS(Lambda)*CosPhi;
    CosC   :=SinPhi1 * SinPhi + CosPhi1 * COSL;
    IF CosC >= 0 THEN
    BEGIN
        X :=R2 * CosPhi * SIN(Lambda);
        Y :=R2 * (CosPhi1 * SinPhi - SinPhi1 * COSL);
    END ELSE X:=-32767;
END;  { Orthographic. }

PROCEDURE Stereographic(Lambda, Lambda0, Phi, Phi1, R : REAL; VAR X, Y :REAL);
{  For R = 1, -2 <= X,Y <= 2. }
VAR K, CosC, CosL, CosPhi, SinPhi, SinPhi1, CosPhi1  : Real;
BEGIN
    Lambda :=Meridian(Lambda, Lambda0);
    IF (Lambda <> Lambda0 + TwoPi) OR (Lambda <> Lambda0-TwoPi) THEN
    BEGIN
        SinPhi:= SIN(Phi);    CosPhi:= COS(Phi);
        SinPhi1:=SIN(Phi1);   CosPhi1:=COS(Phi1);
        CosL:=COS(Lambda) * CosPhi;
        CosC :=SinPhi1 * SinPhi + CosPhi1 * COSL;
        IF CosC >= 0 THEN
        BEGIN
            K := R * 2 / (1 + CosC);
            X := K * CosPhi * SIN(Lambda);
            Y := K * (CosPhi1 * SinPhi - SinPhi1 * COSL) ;
        END
        ELSE X:=-32767;
    END;
END;  { Stereographic. }

PROCEDURE PlotPt(VAR LastPtVis: BOOLEAN);
{   Draws a line from the last point to the current (XP,YP) if it is visible. }
LABEL XIT;
BEGIN
    IF LastPtVis THEN  WRITELN(PLT, 'PD',XP:7:3, ',', YP:7:3,';')
    ELSE WRITELN(PLT, 'PU',XP:7:3, ',', YP:7:3,';');
    LastX:=XP; LastY:=YP; LastPtVis:=TRUE;
  XIT:
END; { PlotPt. }

PROCEDURE CoordinateGrid(OUTLINE: BOOLEAN; MapType: INTEGER);
VAR Longitude, Latitude, MaxLat, LongIncr, LatIncr : INTEGER;
VAR LL, PP, A, RA, R2   : REAL;
VAR X, Y, XN, YN, SINDT, COSDT  : REAL;
BEGIN
    CASE MapType OF
       1: BEGIN   MaxLat:=75;  LongIncr:=360;  LatIncr:=160;    END;
       2: BEGIN   MaxLat:=90;  LongIncr:=360;  LatIncr:=180;    END;
       3: BEGIN   MaxLat:=90;  LongIncr:=360;  LatIncr:=5;      END;
    4..5: BEGIN   MaxLat:=75;  LongIncr:=5;    LatIncr:=5;      END;
    END;  { CASE...}

    LL:=0;  PP:=Phi1;
    IF OUTLINE THEN
    BEGIN IF MapType >  1 THEN MaxLat:=90
          ELSE MaxLat:=80;
          IF MapType >= 5 THEN PP:=0;
    END;

    Latitude:=MaxLat;

    WHILE Latitude >= -MaxLat DO       { Draw parallels }
    BEGIN
        LATR:=Latitude*Radian;
        LastPtVis:=FALSE;

        Longitude:=-180;
        WHILE Longitude <= 180 DO
        BEGIN
            LONGR:=Longitude*Radian;

            CASE MapType OF
            1: BEGIN MERCATOR(LONGR, LL, LATR, R, X, Y);      END;
            2: BEGIN EQUICYL(LONGR, LL, LATR, PP, R, X, Y);   END;
            3: BEGIN MILLER(LONGR, LL, LATR, R, X, Y);        END;
            4: BEGIN SINUSOIDAL(LONGR, LL, LATR, R, X, Y);    END;
            5: BEGIN HAMMER(LONGR, LL, LATR, R, X, Y);        END;
            6: BEGIN ORTHOGRAPHIC( LONGR, LL, LATR, PP, R, X, Y); END;
            7: BEGIN STEREOGRAPHIC(LONGR, LL, LATR, PP, R, X, Y); END;
            END;  { CASE...}

            IF X > -300 THEN
            BEGIN
                XP:=X+XCENTER;
                YP:=Y+YCENTER;
                PlotPt(LastPtVis);
            END ELSE LastPtVis:=FALSE;

            Longitude:=Longitude+LongIncr;
        END;

        IF OUTLINE THEN Latitude:=Latitude-2*MaxLat
        ELSE Latitude:=Latitude-15;    { Latitude:=Latitude-10; }
    END;

    IF OUTLINE THEN LL:=0 ELSE LL:=Lambda0;

    Longitude:=-180;                   { Draw meridians }
    IF OUTLINE AND (MapType > 5) THEN Longitude:=-90;

    IF MapType >= 4 THEN MaxLat:=90;
    WHILE Longitude <= 180 DO
    BEGIN
        LONGR:=Longitude*Radian;
        LastPtVis:=FALSE;
        Latitude:=MaxLat;
        WHILE Latitude >= -MaxLat DO
        BEGIN
            LATR:=Latitude*Radian;

            CASE MapType OF
            1: BEGIN MERCATOR(LONGR, LL, LATR, R, X, Y);      END;
            2: BEGIN EQUICYL(LONGR, LL, LATR, PP, R, X, Y);   END;
            3: BEGIN MILLER(LONGR, LL, LATR, R, X, Y);        END;
            4: BEGIN SINUSOIDAL(LONGR, LL, LATR, R, X, Y);    END;
            5: BEGIN HAMMER(LONGR, LL, LATR, R, X, Y);        END;
            6: BEGIN ORTHOGRAPHIC( LONGR, LL, LATR, PP, R, X, Y); END;
            7: BEGIN STEREOGRAPHIC(LONGR, LL, LATR, PP, R, X, Y); END;
            END;  { CASE...}

            IF X > -300 THEN
            BEGIN
                XP:=X+XCENTER;
                YP:=Y+YCENTER;
                PlotPt(LastPtVis);
            END ELSE LastPtVis:=FALSE;

            Latitude:=Latitude-LatIncr;
        END;

        IF OUTLINE THEN
            IF MapType <= 5 THEN Longitude:=Longitude+360
            ELSE
            Longitude:=Longitude+180
        ELSE Longitude:=Longitude+15;
    END;

    IF NOT OUTLINE AND (MapType = 6) THEN { Draw circular outline. }
    BEGIN
        A:=0;     LastPtVis:=FALSE;
        R2:=R+R;  RA:=R2*Aspect;
        SINDT:=0.05996400648;
        COSDT:=0.99820053993;

        X:=1;  Y:=0;
        XP:=ROUND(XCENTER+RA);  YP:=ROUND(YCENTER);
        PlotPt(LastPtVis);

        WHILE A <= 6.28318 DO
        BEGIN
            XN:= X*COSDT - Y*SINDT;
            YN:= X*SINDT + Y*COSDT;
            X:=XN;  Y:=YN;
            XP:=XCENTER + ROUND(X*RA);  YP:=YCENTER + ROUND(Y*R2);
            PlotPt(LastPtVis);
            A:=A+0.06;
        END;
    END;
END;  { CoordinateGrid. }

PROCEDURE DrawMap(MapType: INTEGER);
VAR Latitude, Longitude : REAL;
VAR LastX : REAL;
BEGIN
    LastPtVis:=FALSE;  LastX:=0;  LastY:=0;

    IF Fname = '' THEN
    ASSIGN(LLF, 'EARTH.LAT')
    ELSE ASSIGN(LLF, Fname);

    RESET(LLF);

    WHILE NOT EOF(LLF) DO
    BEGIN
        READ(LLF, LL);
        LONGR:=LL.LONGI* 1.745329251994329577E-4;
        LATR :=LL.LATI * 1.745329251994329577E-4;

        IF LL.CODE = 'LS' THEN LastPtVis:=FALSE;
        IF (LL.CODE = 'S ') OR (LL.CODE = 'LS') THEN
        BEGIN
            CASE MapType OF
            1: BEGIN MERCATOR(LONGR, Lambda0, LATR, R, X, Y);      END;
            2: BEGIN EQUICYL(LONGR, Lambda0, LATR, Phi1, R, X, Y); END;
            3: BEGIN MILLER(LONGR, Lambda0, LATR, R, X, Y);        END;
            4: BEGIN SINUSOIDAL(LONGR, Lambda0, LATR, R, X, Y);    END;
            5: BEGIN HAMMER(LONGR, Lambda0, LATR, R, X, Y);        END;
            6: BEGIN ORTHOGRAPHIC(LONGR, Lambda0, LATR, Phi1, R, X, Y);  END;
            7: BEGIN STEREOGRAPHIC(LONGR, Lambda0, LATR, Phi1, R, X, Y); END;
            END;  { CASE...}

             IF X > -300 THEN
             BEGIN
                 XP:=X+XCENTER;
                 IF ABS(LastX-XP) > 0.4 THEN LastPtVis:=FALSE;
                 YP:=Y+YCENTER;
                 PlotPt(LastPtVis);  LastX:=XP;
              END ELSE LastPtVis:=FALSE;
        END;
    END;
END;  { DrawMap. }

(* ---------------------  MAIN  PROGRAM  ------------------ *)
VAR RESP : CHAR;
LABEL XIT;
BEGIN
    LastPtVis:=FALSE;  LastX:=0;  LastY:=0;  IPEN:=-1;
    MapType:=1;
    Fname:='';
    Fname:='WORLD4.DAT';

(*                             MENU                                  *)
    WHILE MapType > 0 DO
    BEGIN
        ClrScr;

        GOTOXY(30,1);  WRITE('GLOBE MAP PROJECTIONS');
        GOTOXY(32,3);  WRITELN('SELECT PARAMETERS');

        WRITELN;
        WRITELN(' ':10,'1. Mercator');
        WRITELN(' ':10,'2. Equidistant Cylindrical');
        WRITELN(' ':10,'3. Miller Cylindrical');
        WRITELN(' ':10,'4. Sinusoidal');
        WRITELN(' ':10,'5. Hammer-Aitoff');
        WRITELN(' ':10,'6. Orthographic');
        WRITELN(' ':10,'7. Stereographic');

        WRITELN;
        WRITE(' ':13,'Projection number (1-7) or 0 to quit: ');
        READLN(MapType);
        If MapType = 0 THEN GOTO XIT;

        CASE MapType OF
        1: BEGIN R:=1.82;   END;
        2: BEGIN R:=1.82;   END;
        3: BEGIN R:=1.82;   END;
        4: BEGIN R:=1.82;   END;
        5: BEGIN R:=2.04;   END;
        6: BEGIN R:=2.18;   END;
        7: BEGIN R:=2.18;   END;
        END;

        WRITELN;
        WRITE(' ':13,'Central Longitude of Map (degrees, default = 0): ');
        Lambda0:=0;
        READLN(Lambda0);  Lambda0:=Lambda0*Radian;

        IF (MapType = 2) OR (MapType = 6) OR (MapType = 7) THEN
        BEGIN
            WRITE(' ':13,'Central Latitude  of Map (degrees, -90 - 90): ');
            Phi1:=0;  READLN(Phi1);  Phi1:=Phi1*Radian;
        END;

        WRITE(' ':13, 'Plot grid, continents or both (G/C/B)? ');
        READLN(RESP); RESP:=UPCASE(RESP);
        GRID:=(RESP ='G') OR (RESP = 'B');

        Initialize;

        IF GRID THEN
        BEGIN
            WRITELN(PLT,'SP 2;');
            CoordinateGrid(FALSE,  MapType);
        END;

        WRITELN(PLT,'SP 3;');
        CoordinateGrid(TRUE, MapType);

        WRITELN(PLT, 'SP 1');
        IF (RESP = 'B') OR (RESP = 'C') THEN DrawMap(MapType);

     {  LabelMap(MapType);  }

     WRITELN(PLT, 'PU 0,0; SP 0;');
     CLOSE(PLT);
  XIT:
    END;  { WHILE MapType > 0...}
END.