BASIC programs for surveying and mapping
* * OR * *
BASIC programs for geomatics
By Ian W. Douglas, December 2011
Revised late December 2011 (the revisions are to the program CORDINV.BAS ; see the notes below for that
program)
About 14 years ago, I was a somewhat unsuccessful geomatics student at a local tech school, where I wrote the
following programs. My lack of academic success was due to a combination of a lack of paid experience in surveying
or related endeavours, a lack of experience in a tech school environment (for starters, some of my classmates were
veterans of expulsion from other technical programs and knew more than I did how to navigate through tech school
courses), excessive dependence on calculators and other electronic devices, unwillingness to devote sufficient time
to studies, lack of diplomacy in dealing with other students and instructors, and a combination of these and other
factors. (Regarding this second last point: I have still not learned diplomacy, and my career in the foreign service
would last about a week if I were to pass the entrance exam.) However, I did write programs that worked in
exam/assignment/field work/other test condition situations, whether I was a team member or on my own. The programs
here are sort of "boiler plate" problem solving tools, but I recommend their use by trained geomatics technicians only.
I assume no responsibility if one or more of my programs cause financial loss, territorial disputes, misalignment of
a road, etc.
Acknowledgements: The general guidance of these instructors at the Southern Alberta Institute of Technology (SAIT) made
these programs possible, although the instructors were not involved directly or specifically in the writing of these
programs: Tom Erdman, Ed Scovill, Neil Ormerud, Bob Pearce, Yee Leung, and M. Willman.
The references to the "PC1403H" mean a hand held calculator manufactured by Sharp, which accepts and executes BASIC
programs. I found it useful to have programs for this calculator in the field, for example when repeating angles on a
repeating theodolite.
I have written the programs for BASIC interpreters which do trigonometric calculations in radians, such as Microsoft's
GW-BASIC interpreter. This is not an issue for programs such as my DIFFLEV.BAS, since in a program like this, the only
arithmetic done is addition and subtraction. I have tested these programs successfully using the Microsoft GW-BASIC
interpreter.
My email address is idouglas1304@yahoo.ca .
Table of contents:
[1] DIFFLEV.BAS : differential levelling
[2] CORDINV.BAS : coordinate inversion
[3] MEANREP.BAS : mean of a repeated horizontal angle on a repeating theodolite
[4] TRAV.BAS : traverse computation
[5] CIRCCV. BAS : circular highway curves
[6] AREATRV.BAS : area of a traverse by coordinates
[7] LINELIN.BAS : line-line intersection
[8] CONTR.BAS : contour line interpolation
[1] [Levelling program] This program is for differential levelling, by means of a dumpy level or a similar instrument
shooting at a vertical levelling rod. The program can be modified easily (probably even used as-is) for laser
levelling, three wire levelling, and precise levelling. The program assumes the standard method of entering data in
field notes for levelling, that is, column 1 is for station numbers (e.g. benchmarks and turning points), column 2 is
for backsights, column 3 is for heights of instrument (H.I.s), column 4 is for foresights, and column 5 is for
elevations. (This pattern of assigning columns is not necessary for successful use of the program; I am only citing
a field note setup which is generally accepted in the industry. On the other hand, the program does need the above
types of basic information in order to run.) "Foresight #3" as a prompt by the program is for the user to enter the
foresight in the row in the field note table corresponding to station #3; this row is the third row in the table.
The same principle applies when entering backsights from field notes into the program. The program calculates H.I.s,
elevations, the sum of the foresights, and the sum of backsights. When reducing a differential levelling field note
table which contains intermediate foresights, treat the last station as if it had a foresight instead of an intermediate
foresight, using the intermediate foresight figure as the foresight figure. Program name: DIFFLEV.BAS .
10 CLS:CLEAR:PRINT "Reduce field data, differential levelling."
20 PRINT
30 INPUT "How many stations ";S
40 PRINT
50 DIM BS(S),FS(S),HI(S),EL(S)
60 REM above are arrays for backsights, foresights, heights of instrument, elevations respectively
70 PRINT
80 FOR N=1 TO S
90 IF N=S THEN PRINT "Backsight #";S;"left blank":N=N+1:GOTO 130
100 PRINT "Backsight #";N;
110 INPUT BS(N)
120 NEXT N
130 PRINT
140 FOR O=1 TO S
150 IF O=1 THEN PRINT "Foresight # 1 left blank":O=O+1:GOTO 160
160 PRINT "Foresight #";O;
170 INPUT FS(O)
180 NEXT O
190 PRINT
200 INPUT "Starting elevation ";SE
210 PRINT
220 EL(1)=SE
230 HI(1)=SE+BS(1)
240 FOR E=1 TO S-1
250 EL(E+1)=HI(E)-FS(E+1)
260 HI(E+1)=EL(E+1)+BS(E+1)
270 NEXT E
280 PRINT "H.I.s:":PRINT
290 PRINT
300 FOR X=1 TO S
310 IF X=S THEN PRINT " H.I.";S;"left blank":GOTO 340
320 PRINT HI(X)
330 NEXT X
340 PRINT
350 PRINT "Elevations:"
360 PRINT
370 FOR Y=2 TO S
380 IF Y=2 THEN PRINT " First elevation was given"
390 PRINT EL(Y)
400 NEXT Y
410 PRINT
420 FOR P=1 TO S
430 SB=SB+BS(P)
440 NEXT P
450 FOR Q=1 TO S
460 SF=SF+FS(Q)
470 NEXT Q
480 PRINT "Sum of backsights: ";SB
490 PRINT "Sum of foresights: ";SF
500 PRINT:PRINT"(Ignore this last figure if"
510 PRINT "table contained intermediate foresights)"
520 END
[2] This program calculates (on the one hand) the magnitude and direction of a vector connecting 2 Cartesian points,
and (on the other hand) calculates the end point of a vector of known starting point, magnitude, and direction. All
this is sometimes called "coordinate inversion." The angles are not, strictly speaking, measured from the "standard position," as the remarks in the program indicate, but from lines which could be the x axis or a line parallel to the
x axis. The program also deals with angles expressed as azimuths. Angles are entered as input and generated as output
in decimal degrees, and not in degrees-minutes-seconds or radians. This program assumes that the BASIC interpreter
used is Microsoft GW-BASIC or a similar interpreter; this type of interpreter only works in angles (e.g. for trig
functions) expressed as radians. I have written my program accordingly. Program name: CORDINV.BAS .
I recently revised this program to better handle situations where there was an exact 90 degree, 180 degree, 270 degree,
and 360 degree course between two points, when converting points to a vector. The angles I mentioned above are
both for azimuth and angle measured from the standard position (the actual courses are different, depending on the
system used; I am referring in general to dead on compass rose situations and how my program handles them).
10 CLS:CLEAR:PRINT "Conversion of (on the one hand) two 2 dimensional Cartesian points"
20 PRINT "to vector, or (on the other hand) conversion of a vector of known starting point, magnitude,"
25 PRINT "to an end point."
30 PRINT:PRINT "Convert points to vector (1)"
40 INPUT"or vector and starting point to end point (2) ";CH
50 ON CH GOTO 70, 390
60 PRINT:END
70 CLS:PRINT "Convert points to vector."
80 PRINT "(angle given as standard"
90 PRINT "angle and as azimuth.)"
100 PRINT:PRINT"Enter the 1st point"
110 INPUT "xa (Easting) = ";XA
120 INPUT "ya (Northing) = ";YA
130 PRINT:PRINT "Enter the 2nd point"
140 INPUT "xb (Easting) = ";XB
150 INPUT "yb (Northing) = ";YB
160 DX=XB-XA:DY=YB-YA
170 D=SQR((DX^2)+(DY^2)):REM Distance
180 IF DX=0 THEN GOTO 250
190 REM Above line prevents division by zero in next line
200 S=DY/DX:REM Slope
210 AN=ATN(S):AN=AN*57.29577951#:REM Angle referred to standard position, calculated from slope
220 REM Leave out conversion to degrees on PC1403H
230 IF DX<0 THEN AN=AN+180
240 IF DX>0 AND DY<0 THEN AN=AN+360
250 IF DX=0 AND DY=0 THEN AN=0
260 IF DX=0 AND DY>0 THEN AN=90
270 IF DX<0 AND DY=0 THEN AN=180
280 IF DX=0 AND DY<0 THEN AN=270
290 IF DX>0 AND DY=0 THEN AN=360
300 PRINT "Distance = ";D
310 PRINT "Angle, counterclockwise"
320 PRINT "from standard position:"
330 PRINT AN:REM Note again that this angle is in decimal degrees, not in degrees-minutes-seconds.
340 AA=360-AN+90
350 IF AA>360 THEN AA=AA-360
360 PRINT:PRINT "Angle expressed as"
370 PRINT "azimuth: ";AA
380 GOTO 60
390 CLS:PRINT "Calculate end coordinates"
400 PRINT "of vector, from angle, distance"
410 PRINT "and starting coordinates."
420 PRINT:INPUT "Distance ";D
430 PRINT:PRINT "Enter the starting coordinates"
440 INPUT "xa (Easting) = ";XA
450 INPUT "ya (Northing) = ";YA
460 PRINT:INPUT "Angle = ";AN:AN=AN/57.29577951#
470 PRINT:PRINT "Is angle standard angle (s)"
480 INPUT "or azimuth (a) ";C$
490 IF C$="s" OR C$="S" THEN GOTO 520
500 IF C$="a" OR C$="A" THEN GOTO 540
510 GOTO 470
520 XB=XA+(D*COS(AN)):YB=YA+(D*SIN(AN))
530 GOTO 550
540 XB=XA+(D*SIN(AN)):YB=YA+(D*COS(AN))
550 PRINT:PRINT "xb (Easting) = ";XB
560 PRINT "yb (Northing) = ";YB
570 GOTO 60
[3] [Mean of repeated horizontal angles on repeating theodolite] I never did master the proper field procedure for
using a repeating theodolite, but the basic idea is this: you measure an angle several times, using an instrument
that permits averaging of the readings without the scale suggesting to the instrument operator what the reading
"ought" to be. The procedure (explained at greater length in surveying textbooks and in formal courses) involves
managing circular and measurement errors on the horizontal circle by means of plunging the scope and measuring an
angle repeatedly; only the first and last scale readings need be recorded, together with the number of repetitions.
For example, if you are at a station shooting two prisms exactly 116 degrees 19 min 00 sec apart from your station,
your first reading is 116 degrees 20 min 00 sec, your last reading is 232 degrees 38 min 00 sec, and you repeat the
angle twice, the mean angle as determined by you and your instrument will be 116 degrees 19 min 00 sec. Program
name: MEANREP.BAS .
10 PRINT "Mean of repeated angle on repeating theodolite."
20 INPUT "1st approximation (decimal degrees) ";FA
30 INPUT "Number of repetitions ";NR
40 INPUT "Final horizontal circle reading ";FH
50 A=FA*NR:A=A/360:A=INT(A)
60 B=A*360
70 MA=(FH+B)/NR
80 PRINT "Mean angle (decimal degrees): ";MA
90 END
[4] Program to compute traverses. Program name: TRAV.BAS
Input data: 1 azimuth (decimal degrees)
Northing and easting of first station
Interior angles of traverse (decimal degrees)
Output: Northing and easting of stations
Latitudes and departures
Forward and back azimuths (decimal degrees)
Traverse length
Note: in some traverses which can be assumed to be closed, starting and ending on the same station, the first and last
station may have to counted as two stations, depending on the kind pf problem encountered. For example, in a closed
traverse having four stations, one might have to enter "5" after the "How many stations?" prompt in the program.
10 PRINT:PRINT "Traverse Computation."
20 PRINT "With Distances, Northings,"
30 PRINT "and Eastings."
40 PRINT:PRINT"Clockwise (1) or"
50 INPUT "Anticlockwise (2) ";CH
60 ON CH GOSUB 80, 610
70 END
80 CLEAR:PRINT:PRINT "Clockwise Traverse."
90 PRINT:INPUT "How many stations ";NS
100 DIM FA(NS), BA(NS), IA(NS), D(NS), N(NS), E(NS), LA(NS), DE(NS):REM Arrays of Forward Azimuths,
110 REM Back Azimuths, Interior Angles, Distances, Northings, Eastings, Latitudes,
120 REM Departures respectively
130 DIM RF(NS), RB(NS): REM Arrays for Radian Values of Forward Azimuths and Back Azimuths respectively; leave out
on PC1403H
140 PRINT:PRINT "Forward azimuth of"
150 INPUT "first station ";FA(1)
160 PRINT:PRINT "Northing of"
170 INPUT "station #1 ";N(1)
180 PRINT:PRINT "Easting of"
190 INPUT "station #1 ";E(1)
200 FOR I=2 TO NS
210 PRINT:PRINT "Interior angle of"
220 PRINT "station # ";I;" ";:INPUT IA(I)
230 NEXT I
240 FOR L=1 TO (NS-1)
250 PRINT:PRINT "Distance from station # ";L
260 PRINT "to station # ";L+1;" ";:INPUT D(L)
270 NEXT L
280 FOR J=2 TO NS
290 BA(J)=FA(J-1)-180:IF BA(J)<180 THEN BA(J)=360+BA(J)
300 IF BA(J)>360 THEN BA(J)=BA(J)-360
310 FA(J)=BA(J)-IA(J)
320 IF FA(J)<0 THEN FA(J)=FA(J)+360
330 IF FA(J)>360 THEN FA(J)=FA(J)-360
340 NEXT J
350 FOR M=2 TO NS
360 REM Calculate Latitudes and Departures
370 F(M-1)=FA(M-1)/57.29577951#:LA(M)=D(M-1)*COS(F(M-1)):DE(M)=D(M-1)*SIN(F(M-1)):REM Modify on PC1403H; radians not used
380 NEXT M
390 FOR O=2 TO NS
400 REM Calculate Northings and Eastings
410 N(O)=N(O-1)+LA(O):E(O)=E(O-1)+DE(O)
420 P=P+D(O-1)
430 NEXT O
440 FOR K=2 TO NS
450 PRINT:PRINT "Back Azimuth of station # ";K;" = ";BA(K):REM Arrange for output in Deg/min/sec on PC1403H
460 PRINT "Forward Azimuth of station # ";K;" = ";FA(K):REM Arrange for output in Deg/min/sec on PC1403H
470 PRINT "Northing of station # ";K;" = ";N(K)
480 PRINT "Easting of station # ";K;" = ";E(K)
490 NEXT K
500 PRINT:PRINT "Print latitudes, departures,"
510 INPUT "and total distance ";Y$
520 IF Y$="Y" OR Y$="y" THEN GOTO 530 ELSE GOTO 70
530 PRINT
540 FOR A=2 TO NS
550 PRINT:PRINT "Latitude to station ";A;" : ";LA(A)
560 PRINT "Departure to station ";A;" : ";DE(A)
570 NEXT A
580 PRINT:PRINT "Total length of traverse: ";P
590 GOTO 70
600 REM
610 CLEAR:PRINT:PRINT "Anticlockwise traverse."
620 PRINT:INPUT "How many stations ";NS
630 DIM FA(NS), BA(NS), IA(NS), D(NS), N(NS), E(NS), LA(NS), DE(NS):REM Arrays for Forward Azimuths,
640 REM Back Azimuths, Interior Angles, Distances, Northings, Eastings, Latitudes,
650 REM Departures respectively
660 DIM RF(NS), RB(NS):REM Arrays of Radian Values of Forward Azimuths and Back Azimuths respectively; leave out on PC1403H
670 PRINT:PRINT"Forward azimuth of"
680 INPUT "first station ";FA(1):REM FA(1)=DEG FA(1); arrange for deg/min/sec entry on PC1403H
690 PRINT:PRINT "Northing of"
700 INPUT "station #1 ";N(1)
710 PRINT:PRINT "Easting of"
720 INPUT "station #1";E(1)
730 FOR I=2 TO NS
740 PRINT:PRINT "Interior angle of"
750 PRINT "station # ";I;" " ;:INPUT IA(I):REM IA(I) = DEG IA(I); arrange for deg/min/sec entry on PC1403H
760 NEXT I
770 FOR L=1 TO (NS-1)
780 PRINT:PRINT "Distance from station # ";L
790 PRINT "to station # ";L+1;" ";:INPUT D(L)
800 NEXT L
810 FOR J=2 TO NS
820 BA(J)=FA(J-1)-180:IF BA(J)<180 THEN BA(J)=360+BA(J)
830 IF BA(J)>360 THEN BA(J)=BA(J)-360
840 FA(J)=BA(J)+IA(J)
850 IF FA(J)<0 THEN FA(J)=FA(J)+360
860 IF FA(J)>360 THEN FA(J)=FA(J)-360
870 NEXT J
880 FOR M=2 TO NS
890 REM Calculate Latitudes and Departures
900 F(M-1)=FA(M-1)/57.29577951#:LA(M)=D(M-1)*COS(F(M-1)):DE(M)=D(M-1)*SIN(F(M-1)):REM Modify on PC1403H; radians
not used on that calculator
910 NEXT M
920 FOR O=2 TO NS
930 REM Calculate Northings and Eastings
940 N(O)=N(O-1)+LA(O):E(O)=E(O-1)+DE(O)
950 P=P+D(O-1)
960 NEXT O
970 FOR K=2 TO NS
980 PRINT:PRINT "Back Azimuth of station # ";K;" = ";BA(K):REM arrange for deg/min/sec output on PC1403H
990 PRINT "Forward Azimuth of station # ";K;" = "FA(K):REM arrange for deg/min/sec output on PC1403H
1000 PRINT "Northing of station # ";K;" = ";N(K)
1010 PRINT "Easting of station # ";K;" = ";E(K)
1020 NEXT K
1030 PRINT:PRINT "Print latitudes, departures"
1040 INPUT "and total distance ";Y$
1050 IF Y$="Y" OR Y$="y" THEN GOTO 1060 ELSE GOTO 70
1060 PRINT
1070 FOR A=2 TO NS
1080 PRINT:PRINT "Latitude to station ";A;" : ";LA(A)
1090 PRINT "Departure to station ";A;" : ";DE(A)
1100 NEXT A
1110 PRINT:PRINT "Total length of traverse: ";P
1120 GOTO 70
[5] Program to lay out circular curves, by deflection azimuths and chords, from the Beginning of Curve (BC). This
program will give a "next without for in 590" error message if interval of chainages and/or chords is longer than
the curve arc length. Program name: CIRCCV.BAS .
10 PRINT "Circular highway curve layout"
20 PRINT
30 PRINT "Chainage of Beginning of"
40 PRINT "Curve (assume 0+ and input"
50 PRINT "numbers to right of '+' sign)"
60 INPUT CB
70 INPUT "Radius ";R
80 PRINT "Deflection angle"
90 INPUT "(decimal degrees) ";DA:DA=DA/57.29577951#
100 INPUT "Interval of stations on curve ";IN
110 PRINT "Does curve go left (l) or right (r)"
120 INPUT "from BC ";LR$
130 DT=DA/2
140 L=R*DA:REM Arc Length
150 REM Calculate # of stations on curve
160 NS=(INT(L/IN))+3
170 DIM CH(NS), A(NS), D(NS), C(NS), N(NS), E(NS), AZ(NS)
180 REM Compute chainages of 1st, last stations
190 CH(0)=CB:CH(1)=((INT(CB/IN))*IN)+IN
200 CH(NS-1)=CB+L
210 REM Compute chainages, arcs of remaining stations on curve
220 FOR I=1 TO NS-1
230 IF I=1 THEN GOTO 270
240 IF I=NS-1 THEN GOTO 280
250 CH(I)=CH(I-1)+IN
260 IF CH(I)>CH(NS-1) THEN CH(I)=0:GOTO 290
270 A(I)=CH(I)-CH(0):GOTO 290:REM Arcs
280 A(NS-1)=CH(NS-1)
290 NEXT I
300 REM Compute deflection azimuths
310 FOR J=1 TO NS-1
320 IF J=NS-1 THEN GOTO 340
330 D(J)=(DT/L)*A(J):GOTO 350
340 D(NS-1)=DT
350 IF LR$="L" OR LR$="l" THEN D(J)=6.283185307#-D(J)
360 NEXT J
370 REM Compute chords
380 FOR K=1 TO NS-1
390 C(K)=2*R*SIN(D(K))
400 NEXT K
410 REM Print results
420 PRINT:PRINT "Radius: ";R
430 PRINT:PRINT "Chainage of BC: ";CB
440 PRINT:PRINT "Deflection angle: ";DA*57.29577951#;" decimal degrees"
450 PRINT: PRINT "Chainage" TAB(23) "Chord" TAB(43) "Defl. azimuth":PRINT
460 FOR M=1 TO NS-1
470 IF CH(M)=0 THEN GOTO 540
480 D(M)=D(M)*57.29577951#
490 REM Convert decimal degrees to deg, min, sec
500 DD=INT(D(M)):FD=ABS((DD-D(M)))
510 RM=FD*60:DM=INT(RM):FM=RM-DM:DS=INT(FM*60)
520 REM
530 PRINT CH(M) TAB(21) ABS(C(M)) TAB(41) DD; TAB(45) " d "; TAB(51) DM; TAB(55) " m "; TAB(61) DS; TAB(65) " s "
540 NEXT M
550 PRINT:PRINT "Deflection azimuths are in relation to 'straight-ahead' line from"
560 PRINT "the Beginning of Curve to the Point of Intersection, with this line assumed"
570 PRINT "to have an azimuth of 0 degrees."
580 PRINT:PRINT "Last station is End of Curve."
590 END
[6] Compute, by means of coordinates, the area of a traverse having any number of sides. Coordinates can be in
any quadrant in the Cartesian plane. See Kavanagh and Bird, "Surveying," pp. 197-199. Program name: AREATRAV.BAS .
10 PRINT "Area of a traverse":REM "A"
20 PRINT "by coordinates."
30 PRINT:INPUT "How many stations";NS
40 DIM N(NS),E(NS):REM Arrays for northings and eastings
50 PRINT:FOR I=1 TO NS
60 PRINT "Northing of station ";I;:INPUT N(I)
70 PRINT "Easting of station ";I;:INPUT E(I)
80 PRINT
90 NEXT I
100 A=E(1)*(N(2)-N(NS))
110 FOR J=2 TO NS-1
120 A=A+(E(J)*(N(J+1)-N(J-1)))
130 NEXT J
140 A=A+(E(NS)*(N(1)-N(NS-1)))
150 A=ABS(A/2)
160 PRINT:PRINT "Area = ";A
170 END
[7] Line-line intersection. Program name: LINELIN.BAS .
10 CLS:CLEAR:PRINT "Intersection coordinates"
20 PRINT "of two lines of unknown"
30 PRINT "distance, known starting coordinates,"
40 PRINT "and known angles."
50 PRINT:PRINT "Enter data for the 1st line."
60 PRINT:INPUT "angle = ";AA
70 PRINT:PRINT "Starting Coordinates:"
80 INPUT "xa (easting) = ";XA
90 INPUT "ya (northing) = ";YA
100 PRINT:PRINT "Enter data for the 2nd line."
110 PRINT:INPUT "angle = ";AB
120 PRINT:PRINT "Starting Coordinates:"
130 PRINT:INPUT "xb (easting) = ";XB
140 INPUT "yb (northing) = ";YB
150 PRINT:PRINT "Are angles standard"
160 INPUT "angles (s) or azimuths (a)";A$
170 IF A$="S" OR A$="s" THEN AA=360-AA+90:IF AA>360 THEN AA=AA-360
180 IF A$="S" OR A$="s" THEN AB=360-AB+90:IF AB>360 THEN AB=AB-360
190 AA=AA/57.29577951#:AB=AB/57.29577951#:REM leave out conversion to radians on PC1403H
200 T=AB-AA:DX=XB-XA:DY=YB-YA
210 S1=((DX*COS(AB))-(DY*SIN(AB)))/-SIN(T)
220 XI=XA+(S1*SIN(AA)):YI=YA+(S1*COS(AA))
230 PRINT:PRINT:PRINT "Coordinates of intersection point:"
240 PRINT:PRINT "x (easting)= ";XI:PRINT "y (northing)= ";YI
250 PRINT:END
[8] [Contour line interpolation program] I figured out how to write this program long after first encountering the
problem I intended to solve with it, and too late to enable me to pass this part of the course. This program deals
with the radial line topography survey problem, for example, making a topographic map of a small hill with the transit
set up on top of the hill. You shoot along several horizontal rays, for example every 20 degrees in azimuth (e.g.
North, 20 degrees azimuth, 40 degrees, 60 degrees, 80 degrees--10 more degrees would mean shooting due East--etc.)
Along each ray, you use a tape or Electronic Distance Meter (EDM) to get the distance from the transit and the angle
of depression to each station shot thus along the ray. The raw data, using the following or a similar program, can
form the basis of a topographic map with a contour interval of a round number, e.g. every 2 metres, 10 metres, etc.
If the height of instrument above ground level is the same as the height of the pogo stick prism (again above ground
level), then you are virtually measuring ground level to ground level characteristics of the land, namely azimuth
versus amount of depression/elevation, and distance. By the way, this program forms the nucleus of a program I wrote
later to do normalization of data on 2 dimensional line/scatter graphs for work in a drug testing lab. Program
name: CONTR.BAS .
10 CLS:CLEAR:PRINT "Contour line interpolation"
20 PRINT:INPUT "Horizontal distance between transit and target ";HD
30 INPUT "Desired contour interval ";CI
40 INPUT "Elevation of ground at transit ";ER
50 INPUT "Elevation of ground at target ";EA
60 DE=EA-ER:IF DE=0 THEN EA=EA+.0000001: REM Last statement prevents division by zero later on in program
70 REM by making program "believe" that there is still an extremely small difference in elevation
80 Q=ER/CI:QI=INT(Q):QM=QI*CI
90 IF DE<0 THEN GOTO 120
100 A=INT(QM)+CI:IF A>=EA THEN GOTO 410
110 GOTO 130
120 A=INT(QM):IF A<=EA THEN GOTO 410
130 NC=1:CP=A
140 IF DE<0 THEN GOTO 170
150 CP=CP+CI:IF CP>EA THEN GOTO 200
160 NC=NC+1:GOTO 150
170 CP=CP-CI:IF CP0 AND DY<0 THEN XYAN=XYAN+360
340 IF DX=0 AND DY=0 THEN XYAN=0
350 IF DX=0 AND DY>0 THEN XYAN=90
360 IF DX<0 AND DY=0 THEN XYAN=180
370 IF DX=0 AND DY<0 THEN XYAN=270
380 IF DX>0 AND DY=0 THEN XYAN=360
390 REM End getting XY angle into correct quadrant, referred to standard position
400 IF DX=0 THEN GOTO 480
410 REM Above line prevents division by zero in next line
420 XZSLOPE=DZ/DX:REM Slope in, or reduced to, XZ plane
430 XZAN=ATN(XZSLOPE):XZAN=XZAN*57.29577951#:REM XZ Angle referred to standard position, calculated from slope
440 REM Leave out conversion to degrees on PC1403H
450 REM Begin getting XZ angle into correct quadrant, referred to standard position
460 IF DX<0 THEN XZAN=XZAN+180
470 IF DX>0 AND DZ<0 THEN XZAN=XZAN+360
480 IF DX=0 AND DZ=0 THEN XZAN=0
490 IF DX=0 AND DZ>0 THEN XZAN=90
500 IF DX<0 AND DZ=0 THEN XZAN=180
510 IF DX=0 AND DZ<0 THEN XZAN=270
520 IF DX>0 AND DZ=0 THEN XZAN=360
530 REM End getting XZ angle into correct quadrant, referred to standard position
540 REM Begin calculating horizontal and vertical distances of 2nd point from 1st point.
550 REM Horizontal separation is absolute horizontal distance as seen when looking onto
560 REM XY plane, but otherwise not reduced to any other plane.
570 VESEP=DZ:HOSEP=SQR((DX^2)+(DY^2))
580 REM End calculating horizontal and vertical distances of 2nd point from 1st point
590 REM Begin calculating absolute angle of elevation of 2nd point from 1st point, that is,
600 REM angle above XY plane between the two points reduced to no plane but its own (see below).
610 IF HOSEP=0 THEN GOTO 680
620 REM Above line prevents division by zero in next line
630 ELANG=ATN(VESEP/HOSEP):ELANG=ELANG*57.29577951#:REM vector's absolute angle of elevation above
640 REM a plane 'parallel' (in 3 dimensions) to XY plane and touching the vector's start point. This angle is
650 REM reduced to no plane but its own, and is between +90 degrees and -90 degrees; negative angles of
660 REM elevation are angles of depression.
670 REM Begin traps of situations where HOSEP is zero.
680 IF HOSEP=0 AND VESEP=0 THEN ELANG=0
690 IF HOSEP=0 AND VESEP>0 THEN ELANG=90
700 IF HOSEP=0 AND VESEP<0 THEN ELANG=-90
710 REM End traps of situations where HOSEP is zero.
720 REM End calculation of vector's absolute angle of elevation above a plane 'parallel' (in 3 dimensions) above
730 REM XY plane (see above)
740 PRINT:PRINT "Distance = ";DIST
750 PRINT:PRINT "XY Angle, counterclockwise from standard position:"
760 PRINT XYAN;" decimal degrees"
770 XYAAAN=360-XYAN+90
780 IF XYAAAN>360 THEN XYAAAN=XYAAAN-360
790 PRINT "XY Angle expressed as azimuth:"
800 PRINT XYAAAN;" decimal degrees"
810 PRINT:PRINT "XZ Angle, counterclockwise from standard position:"
820 PRINT XZAN;" decimal degrees"
830 PRINT:PRINT "Vector's absolute angle of elevation above a plane touching the vector's start point, the plane"
840 PRINT "being 'parallel' (in 3 dimensions) to the XY plane; this angle of elevation is reduced to its own"
850 PRINT "plane but no other except by coincidence; this angle must be between +90 degrees and -90 degrees;"
860 PRINT "negative angles of elevation are angles of depression = ";ELANG;" decimal degrees":PRINT
870 GOTO 120
880 PRINT "Calculate end coordinates of 3D vector, given distance, XY angle, angle of elevation/depression,"
890 PRINT "and starting coordinates. Program accepts XY (horizontal plane) angles as either azimuths or angles"
900 PRINT "referred to standard position. Program only accepts angles of elevation between +90 and -90 degrees;"
910 PRINT "negative angles of elevation are angles of depression."
920 PRINT:INPUT "Distance ";DIST
930 PRINT:PRINT "Enter the starting coordinates"
940 INPUT "xa (Easting) = ";XA
950 INPUT "ya (Northing) = ";YA
960 INPUT "za (vertical distance from z = 0) = ";ZA
970 PRINT:INPUT "XY Angle (decimal degrees) = ";XYAN:XYAN=XYAN/57.29577951#
980 PRINT:PRINT "Vector's absolute angle of elevation above a plane touching the vector's start point, the"
990 PRINT "just-previously-mentioned plane being 'parallel' (in 3 dimensions) to XY plane (negative angles are"
1000 PRINT "angles of depression) (this angle reduced to its own plane but no other except by coincidence) (decimal"
1010 INPUT "degrees) = ";ELANG:PRINT
1020 IF 90