Category : Science and Education
Archive   : WDBDISK5.ZIP
Filename : SATPLOT.PAS

PROGRAM SAT_PLOT(INPUT, OUTPUT);
{This program is unclassified and can be freely adapted and incorporated by
others. Appropriate credits should be given.

Original version created at the Applied Physics Lab by George Swartout.

Nov 81 - Adapted by Ed smith (JSCS) to the HP-9845.

Nov 83 - Reloaded, modified, and enhanced by Bob Cole (SAC/NRA).

Jun - Jul 86 - Converted from HP BASIC to MS Pascal and Flexi_graph by Fred
Pospeschil with math analysis support by Gerald Wendling.

COMMENTS ON PROGRAM STRUCTURE.

Symbolic analysis of the original formulas and program code indicated that
certain portions of the code always produced the same values once initial
setup was performed. Two formulas always produced a value of zero and did
not influence the overall operation of the program.

Based on this analysis several computations were eliminated and five of the
equasions in the SCALE_COORDS procedure were restructured to minimize the
total calculation load. all code which produced constant values was
collected into the COMPUTE_CONSTANTS procedure.

To further enhance performance all primary variables were declared global
to avoid loading and unloading up to 650,000 parameters on and off the
stack per execution of the program. Although these actions result in code
which is a little more difficult to read it produced a runtime performance
improvement of about 45%.

The essence of the program is contained in the procedures
COMPUTE_CONSTANTS, SCALE_COORDS, DRAW_GRID, AND DISPLAY_MAP.

This program is designed to use the Micro World Data Bank II data base.
}

CONST
R5 = 3432; {earth radius in NM}
PI = 3.14159265;

TYPE
FGREAL = REAL4; {Symbolic definition of the REAL precision}
{to make it easy to change from single to}
{double precsion throughout the program and}
{the graphics routines.}
TURTLE_TYPE = REAL4; {Symbolic definition as above. needed by}
{the external graphic routine definitions}
{which are in the 'included' file FGPASDEF.INC}

FGSTRING = LSTRING(80); {String definition needed by the text I/O}
{routine in the 'included' file FGREADLN.INC}

POINT = RECORD {Micro World Data Base II input record}
CODE : INTEGER;{Identifies the beginning of lines and line}
{Types. See the Micro World Data Base Doc}
LAT : INTEGER;{Latitude in signed minutes}
LON : INTEGER;{Longitude in signed minutes}
END;


VAR
SYSTEM : INTEGER; {Stores system identifier 1-5}
PRINTER : TEXT; {Printer is TEXT output device}
PLOTTER : TEXT; {Plotter is written to like an ASCII file}
PLOTTER_ON : BOOLEAN; {TRUE = output to plotter, FALSE = CRT}
KEYBOARD : FGSTRING; {Buffer for keyboard input}
STRBUFF : FGSTRING; {Buffer to hold output text strings}
MAP_FILE : FILE OF POINT; {Symbolic name for input map data file}
FILE_NAME : LSTRING(30); {Stores map file name input by user}
NEXT_POINT : POINT; {Holds next point to be plotted}
LAST_POINT : POINT; {Holds last point plotted}
INIT_VEC : ARRAY [1..7] OF INTEGER; {graphics initialization parms}
DEG_TO_RAD : FGREAL; {Factor to convert from degrees to radians}
MIN_TO_RAD : FGREAL; {Factor to convert from minutes to radians}
X_MIN : FGREAL; {Min world X or LONGITUDE value}
X_MAX : FGREAL; {Max " " " " " }
Y_MIN : FGREAL; {Min " Y " LATITUDE " }
Y_MAX : FGREAL; {Max " " " " " }
X_Y_ASPECT : FGREAL; {X / Y aspect ratio of the display screen}

WORLDX_TO_PLOTTER : FGREAL; {World to plotter scale factor for X / LON}
WORLDY_TO_PLOTTER : FGREAL; {World to plotter scale factor for Y / LAT}
PLOTTER_X_CENTER : INTEGER; {Center of the plotter X axis}
PLOTTER_Y_CENTER : INTEGER; {Center of the plotter Y axis}

L1 : FGREAL; {Satellite latitude in radians}
L2 : FGREAL; {Satellite longitude in radians}
H0 : FGREAL; {Satellite altitude in nautical miles}

SAT_LAT : FGREAL; {Satellite latitude in degrees}
SAT_LON : FGREAL; {Satellite longitude in degrees}
SAT_ALT : FGREAL; {Satellite altitude in nautical miles}
S : FGREAL; {Distance from point on earth to satellite}
S0 : FGREAL; {Distance from the satellite to the horizon}
S1 : FGREAL; {Projected X component}
S2 : FGREAL; {Projected Y component}
S3 : FGREAL; {Projected Z component}
LAT : FGREAL; {Unprojected latitude in radians}
LON : FGREAL; {Unprojected longitude in radians}
PROJ_LAT : FGREAL; {Projected latitude in radians}
PROJ_LON : FGREAL; {Projected longitude in radians}
G0 : FGREAL; {Largest distance from satellite coords at}
{which points will be plotted. i.e. the }
{maximum projected LAT/LON distance in NM}

{Constants}
SINL1 : FGREAL; {which are computed}
COSL1 : FGREAL; {and used by the coordinate}
SINL2 : FGREAL; {projection and scaling routine}
COSL2 : FGREAL;
SINL1SINL2 : FGREAL;
SINL1COSL2 : FGREAL;
S1ADJUST : FGREAL;
S2ADJUST : FGREAL;
S3ADJUST : FGREAL;
VISABLE : BOOLEAN;


{$INCLUDE:'FGPASDEF.INC'} {Flexi-Graph graphics routine definitions}
{$INCLUDE:'IDSHDUMP.INC'} {Screen dump routine for IDS Microprism printer}
{$INCLUDE:'FGREADLN.INC'} {Routine to do non-distructive text I/O in the}
{graphics mode}

{NOTE : All of the above included routines are provided as part of the}
{Flexi-Graph graphics support system and are provided along with this}
{program to improve the utility of this public domain program. These}
{files are provided with the permission of NOGDS and Micro Doc}

{*************************************}
{This procedure computes a number of values which remain constant
throughout the execution of the program. By pre-computing these values
considerable performance improvements were realized.}

PROCEDURE COMPUTE_CONSTANTS;
VAR
X, Y : INTEGER;

BEGIN
DEG_TO_RAD := PI / 180.0;
MIN_TO_RAD := DEG_TO_RAD / 60.0;
L1 := L1 * DEG_TO_RAD;
L2 := L2 * DEG_TO_RAD;

S0 := SQRT(H0 * (2*R5 + H0)); {Compute distance from satellite}
{to horizon}
{Compute distance from horizon to line from the satellite to the}
{the center of the earth}
G0 := (R5*SQRT(SQR(R5+H0) - SQR(R5))) / (R5+H0);

FGMAXS(X, Y); {Get max X and Y for the selected CRT display}
{Compute X/Y CRT ratio so that the world coordinate window can be
adjusted to the same ratio}
X_Y_ASPECT := FLOAT(X) / FLOAT(Y);
X_MIN := -G0 * X_Y_ASPECT; {Here we compute the min and max X/Y world}
X_MAX := G0 * X_Y_ASPECT; {coordinates}
Y_MIN := -G0;
Y_MAX := G0;

SINL1 := SIN(L1); {Compute static values used by the}
COSL1 := COS(L1); {SCALE_COORDS procedure}
SINL2 := SIN(L2);
COSL2 := COS(L2);
SINL1SINL2 := SINL1 * SINL2;
SINL1COSL2 := SINL1 * COSL2;
S1ADJUST := (R5+H0) * COSL1 * SINL2;
S2ADJUST := (R5+H0) * COSL1 * COSL2;
S3ADJUST := (R5+H0) * SINL1;
END;

{*************************************}
PROCEDURE SCALE_COORDS;
{Here we scale the coordinates to reflect the curvature of the earth}
{and determine if the points are visable.}

BEGIN
S1 := R5 * COS(LAT) * SIN(LON) - S1ADJUST;
S2 := R5 * COS(LAT) * COS(LON) - S2ADJUST;
S3 := R5 * SIN(LAT) - S3ADJUST;
S := SQRT(SQR(S1) + SQR(S2) + SQR(S3));

IF S <= S0 THEN {If the distance from the satellite to the point}
VISABLE := TRUE {is less than or equal to the distance from the}

ELSE {satellite to the computed horizon then the point}
VISABLE := FALSE; {is visable.}

PROJ_LAT := (-S1 * SINL1SINL2) - (S2 * SINL1COSL2) + (S3 * COSL1);
PROJ_LON := (S1 * COSL2) - (S2 * SINL2);

END;

{*************************************}
PROCEDURE SETUP_PLOTTER(PLOTTER_X_EXTENT, PLOTTER_Y_EXTENT : INTEGER);
{This procedure sets up values needed to scale the points for display}
{on a plotter. It assumes that the plotter also uses the standard}
{convention of X=0 and Y=0 is at the lower left corner of the plotting}
{surface}

VAR
WORLD_X_EXTENT : FGREAL;
WORLD_Y_EXTENT : FGREAL;

BEGIN
WORLD_X_EXTENT := X_MAX - X_MIN;
WORLD_Y_EXTENT := Y_MAX - Y_MIN;
WORLDX_TO_PLOTTER := PLOTTER_X_EXTENT / WORLD_X_EXTENT;
WORLDY_TO_PLOTTER := PLOTTER_Y_EXTENT / WORLD_Y_EXTENT;
PLOTTER_X_CENTER := PLOTTER_X_EXTENT DIV 2;
PLOTTER_Y_CENTER := PLOTTER_Y_EXTENT DIV 2;
END;


{The following procedures PLOTTER_SCALE, PLOTTER_MOVE, PLOTTER_PLOT,
PLOTTER_DRAW and PLOTTER_CIRC perform the operations needed by the plotter
which are automatically handled by the Flexi-Graph CRT graphics routines
The WRITELN statements write the necessary text strings for a Mannesmann
Tally PIXY3 plotter and will have to be changed for other plotters Only the
most basic plotter commands were used so as to make this code as portable
to other plotters as possible.}

{*************************************}
PROCEDURE PLOT_SCALE(WORLD_X, WORLD_Y : FGREAL;
VAR PLOTTER_X, PLOTTER_Y : INTEGER);
BEGIN
PLOTTER_X := TRUNC(PLOTTER_X_CENTER + (WORLD_X * WORLDX_TO_PLOTTER));
PLOTTER_Y := TRUNC(PLOTTER_Y_CENTER + (WORLD_Y * WORLDY_TO_PLOTTER));
END;

{*************************************}
PROCEDURE PLOTTER_MOVE(WORLD_X, WORLD_Y : FGREAL);
VAR
PLOTTER_X : INTEGER;
PLOTTER_Y : INTEGER;

BEGIN
PLOT_SCALE(WORLD_X, WORLD_Y, PLOTTER_X, PLOTTER_Y);
WRITELN(PLOTTER,'M ',PLOTTER_X:5, PLOTTER_Y:5);
END;

{*************************************}
PROCEDURE PLOTTER_PLOT(WORLD_X, WORLD_Y : FGREAL);
VAR
PLOTTER_X : INTEGER;
PLOTTER_Y : INTEGER;

BEGIN
PLOT_SCALE(WORLD_X, WORLD_Y, PLOTTER_X, PLOTTER_Y);
WRITELN(PLOTTER,'M ',PLOTTER_X:5, PLOTTER_Y:5);
WRITELN(PLOTTER,'I ',0:5,0:5);
END;

{*************************************}
PROCEDURE PLOTTER_DRAW(WORLD_X, WORLD_Y : FGREAL);
VAR
PLOTTER_X : INTEGER;
PLOTTER_Y : INTEGER;

BEGIN
PLOT_SCALE(WORLD_X, WORLD_Y, PLOTTER_X, PLOTTER_Y);
WRITELN(PLOTTER,'D ',PLOTTER_X:5, PLOTTER_Y:5);
END;

{*************************************}
PROCEDURE PLOTTER_CIRC(WORLD_X, WORLD_Y, WORLD_RADIUS : FGREAL);
VAR
PLOTTER_X : INTEGER;
PLOTTER_Y : INTEGER;
PLOTTER_R : INTEGER;

BEGIN
PLOT_SCALE(WORLD_X, WORLD_Y, PLOTTER_X, PLOTTER_Y);
PLOTTER_R := TRUNC(WORLD_RADIUS * WORLDX_TO_PLOTTER);
WRITELN(PLOTTER,'W ',PLOTTER_X:5, PLOTTER_Y:5,
PLOTTER_R:5, PLOTTER_R:5,
0:5, 3600:5);
END;


{*************************************}
PROCEDURE DRAW_GRID;
{This procedure draws the horizon as a solid circle and the latitude}
{and longitude lines as dots spaced at three degree intervals.}

VAR
I : INTEGER;
J : INTEGER;

BEGIN

IF PLOTTER_ON THEN
BEGIN
WRITELN(PLOTTER,'J3'); {Select plotter pen 3}
PLOTTER_CIRC(0.0, 0.0, G0);
END
ELSE
BEGIN
LINCLR(3); {Select CRT color 3}
CIRC(0.0, 0.0, G0);
END;

FOR I := - 5 TO 5 DO
BEGIN
LAT := FLOAT(I * 15) * DEG_TO_RAD;
FOR J := 0 TO 120 DO
BEGIN
LON := FLOAT(J * 3) * DEG_TO_RAD;
SCALE_COORDS;
IF VISABLE THEN
IF PLOTTER_ON THEN
PLOTTER_PLOT(PROJ_LON, PROJ_LAT)
ELSE
PLOT(PROJ_LON, PROJ_LAT);
END;
END;

FOR I := 0 TO 23 DO
BEGIN
LON := FLOAT(I * 15) * DEG_TO_RAD;
FOR J := -29 TO 29 DO
BEGIN
LAT := FLOAT(J * 3) * DEG_TO_RAD;
SCALE_COORDS;
IF VISABLE THEN
IF PLOTTER_ON THEN
PLOTTER_PLOT(PROJ_LON, PROJ_LAT)
ELSE
PLOT(PROJ_LON, PROJ_LAT);
END;
END;
END;

{*************************************}
PROCEDURE DISPLAY_MAP;


PROCEDURE OUTPUT_TO_CRT;
VAR
LINE_TYPE : INTEGER;

BEGIN
{now we build the display}
WHILE NOT EOF(MAP_FILE) DO {process all records in file}
BEGIN
READ(MAP_FILE, NEXT_POINT); {get a record}
LAT := FLOAT(NEXT_POINT.LAT) * MIN_TO_RAD;
LON := FLOAT(NEXT_POINT.LON) * MIN_TO_RAD;
SCALE_COORDS;

IF VISABLE THEN
BEGIN
IF (NEXT_POINT.CODE > 5) THEN {check for first point in each line}
BEGIN
LINE_TYPE := NEXT_POINT.CODE DIV 1000;

IF LINE_TYPE = 2 THEN {country boundries}
LINCLR(2)
ELSE IF LINE_TYPE IN [4,7] THEN {state boundries, and rivers}
LINCLR(3)
ELSE
LINCLR(1); {coast lines, islands, lakes}

PLOT(PROJ_LON, PROJ_LAT); {plot the starting point}
END
ELSE
DRAW(PROJ_LON, PROJ_LAT); {draw line}
END

ELSE {the point / line is not visable}
MOVE(PROJ_LON, PROJ_LAT);

END; {display completed}
END; {end procedure OUTPUT_TO_CRT}



PROCEDURE OUTPUT_TO_PLOTTER;
VAR
LINE_TYPE : INTEGER;

BEGIN
{now we build the display on the plotter}
WHILE NOT EOF(MAP_FILE) DO {process all records in file}
BEGIN
READ(MAP_FILE, NEXT_POINT); {get a record}
LAT := FLOAT(NEXT_POINT.LAT) * MIN_TO_RAD;
LON := FLOAT(NEXT_POINT.LON) * MIN_TO_RAD;
SCALE_COORDS;

IF VISABLE THEN
BEGIN
IF (NEXT_POINT.CODE > 5) THEN {check for first point in each line}
BEGIN
LINE_TYPE := NEXT_POINT.CODE DIV 1000;

IF LINE_TYPE = 2 THEN {country boundries}
WRITELN(PLOTTER,'J2')
ELSE IF LINE_TYPE IN [4,7] THEN {state boundries, and rivers}
WRITELN(PLOTTER,'J3')
ELSE
WRITELN(PLOTTER,'J1'); {coast lines, islands, lakes}

PLOTTER_PLOT(PROJ_LON, PROJ_LAT); {plot the starting point}
END
ELSE
PLOTTER_DRAW(PROJ_LON, PROJ_LAT) {draw to the point}
END
ELSE {the point / line is not visable}
PLOTTER_MOVE(PROJ_LON, PROJ_LAT); {move to the point}
END; {display completed}
END; {end procedure OUTPUT_TO_PLOTTER}

BEGIN {main code}
IF PLOTTER_ON THEN {If plotter output was selected}
OUTPUT_TO_PLOTTER {output to plotter}
ELSE {else}
OUTPUT_TO_CRT; {output to the CRT}

WRITE(CHR(7)); {signal plot is done}

END;

{*************************************}
BEGIN
{have the user identify the system configuration}
REPEAT
WRITELN(' Satellite Projections');
WRITELN(' using');
WRITELN(' Micro World Data Bank II ');
WRITELN;
WRITELN(' What type of microcomputer are you running ?');
WRITELN;
WRITELN(' 1 = IBM PC or compatable with color graphics adapter');
WRITELN(' 2 = IBM AT with enhanced color graphics adapter');
WRITELN(' 3 = Z110 or Z120 with 32K video memory chips');
WRITELN(' 4 = Z110 or Z120 with 64K video memory chips');
WRITELN(' 5 = IBM PC with IDS 440 video display card');

READLN(SYSTEM); {get type of computer and display system}
UNTIL SYSTEM IN [1..5];

WRITELN;
WRITELN;
WRITELN;

REPEAT
WRITE('Enter satellite latitude <-90.0 thru 90.0> ');
READLN(L1);
UNTIL (L1 >= -90.0) AND (L1 <= 90.0);
SAT_LAT := L1;

REPEAT
WRITE('Enter satellite longitude <-180.0 thru 180.0> ');
READLN(L2);
UNTIL (L2 >= -180.0) AND (L2 <= 180.0);
SAT_LON := L2;

REPEAT
WRITE('Enter satellite altitude in NM <1.0 thru 30000.0> ');
READLN(H0);
UNTIL (H0 >= 1.0) AND (H0 <= 30000.0);

WRITELN;
WRITE('Enter the name of the map data file > ');
READLN(FILE_NAME);
ASSIGN(MAP_FILE, FILE_NAME); {identify map file to operating system}
RESET(MAP_FILE); {position to beginning of file}


{initialize the graphics display and Flexi-Graph routines}
CASE SYSTEM OF
1 : GOPEN(0, 1); {IBM PC 320x200}
2 : GOPEN(1, 0); {EGA 640 x 350}
3 : GOPEN(2, 1); {Z100 640x240}
4 : GRSTRT(1,1,1); {Z100 640x480}
5 : GRSTRT(2,0,0); {PC / IDS 640x400}
END;

COMPUTE_CONSTANTS;

NOCLIP; {Turn off Flexi-Graph clipping as}
{this program keeps all points}
{within the world coordinate space}
VECABS; {Plot and Draw with absolute coords}
FGSFNT; {Select standard Flexi-Graph font}
WINDOW(X_MIN, X_MAX, Y_MIN, Y_MAX); {Set the bounds of the world coords}

PLOTTER_ON := FALSE; {On initial pass, outuput to the CRT}
DRAW_GRID;
DISPLAY_MAP;
CLOSE(MAP_FILE); {close map data input file}

SETTOG; {Enter toggle pixels mode}
MOVE(X_MIN + 200, Y_MIN + 200); {Move to lower left corner}
STRBUFF := 'Dump to printer '; {Set up prompt message}
TYPET(1,0,ORD(STRBUFF.LEN),STRBUFF[1]); {Display it on the CRT}
FGREADLN(KEYBOARD); {See if dump of screen to printer}
{is desired}
MOVE(X_MIN + 200, Y_MIN + 200); {Move back to beginning of prompt}
TYPET(1,0,ORD(STRBUFF.LEN),STRBUFF[1]); {Toggle/erase prompt message}

IF KEYBOARD[1] = 'Y' THEN {is desired}
BEGIN
ASSIGN(PRINTER, 'PRN');
REWRITE(PRINTER);
WRITELN(PRINTER, ' SATELLITE PLOT OF EARTH');
WRITELN(PRINTER);
WRITELN(PRINTER, ' SATELLITE LATITUDE = ',SAT_LAT:8:2, ' DEGREES');
WRITELN(PRINTER, ' SATELLITE LONGITUDE = ', SAT_LON:8:2, ' DEGREES');
WRITELN(PRINTER, ' SATELLITE ALTITUDE = ', H0:8:2, ' NAUTICAL MILES');
WRITELN(PRINTER, ' MAP DATA BASE USED WAS ', FILE_NAME);
WRITELN(PRINTER);
CLOSE(PRINTER);
IDSHDUMP(0); {Dump the image to IDS Microprism printer. The}
END; {"0" tells the dump routine that the}
{background conor of the display is 0.}

MOVE(X_MIN + 200, Y_MIN + 200);
STRBUFF := 'Dump to plotter ';
TYPET(1,0,ORD(STRBUFF.LEN),STRBUFF[1]);
FGREADLN(KEYBOARD); {See if a plotter copy of the}
MOVE(X_MIN + 200, Y_MIN + 200); {image is wanted}
TYPET(1,0,ORD(STRBUFF.LEN),STRBUFF[1]);
CLRTOG; {Exit the toggle mode}

IF KEYBOARD[1] = 'Y' THEN
BEGIN
ASSIGN(MAP_FILE, FILE_NAME); {Identify map file to operating system}
RESET(MAP_FILE); {position to beginning of file}
PLOTTER_ON := TRUE;
ASSIGN(PLOTTER, 'AUX:');
REWRITE(PLOTTER);
SETUP_PLOTTER(2400, 1800); {PIXY3 plotter has a physical address}
{space of X = 0 thru 2400 and}
{Y = 0 thru 1800}
WRITELN(PLOTTER); {These two CRs tell the plotter that}
WRITELN(PLOTTER); {command strings will be delimited}
{with a single CR}
WRITELN(PLOTTER,'T0'); {Set plotter pen speed to maximum}
WRITELN(PLOTTER, 'J1'); {Select plotter pen #1}
WRITELN(PLOTTER,'H'); {Home the plotter pen}
DRAW_GRID;
DISPLAY_MAP;
WRITELN(PLOTTER,'J0'); {Park all pens at the end of plotting}
CLOSE(PLOTTER); {Flush output buffer to plotter}
CLOSE(MAP_FILE); {close map data input file}
END;

{exit the graphics mode and leave the system in normal text mode}
CASE SYSTEM OF
1 : GCLOSE;
2 : GCLOSE;
3 : GCLOSE;
4 : GRSTOP;
5 : GRSTOP;
END;
END.