$lib graph on
data 0,1 , 200,1 , 440,1 , 680,2 , 880,2 , 1040,1 , 1200,1
data 1360,0 , 1520,0 , 1720,0 , 1880,3 , 1960,3 , 2040,4
dim rowdata%(12 , 1) , spotheights(800) , groundheight(800), spotngr$(800,1) , oldly(800 , 1)
pi = 3.141592653589
k = 1.33
distbreak = 120         'LOS lines or mouse cursor
def fnacs(x) = atn(sqr((1 - x * x) / x))
screen 12
for rowcount% = 0 to 12
	read rowdata%(rowcount% , 0)
	read rowdata%(rowcount% , 1)
next rowcount%
f$ = "heights2.dat"
if dir$(f$) = "" then
	f$ = "d:" + f$
	if dir$(f$) > "" then exit if
	sound 500 , 2
	color 12
	print "Height database missing"
	end
end if
if command$ > "" then
	insp% = instr(command$ , " ")
	if insp% = 0 then
		baselo$ = command$
		exit if
	else
		baselo$ = left$(command$ , insp% - 1)
		distlo$ = mid$(command$ , insp% + 1)
	end if
	print baselo$ , distlo$
	insp% = -1      'now a flag
end if

open "b" , 1 , f$
'==================================================================
cls
color 15
line (145 , 60) - (500 , 82) , 13 , b
locate 5 , 20
print "TERRAIN PLOT BETWEEN TWO SPECIFIED STATIONS"
color 7
line (50 , 105) - (595 , 225) , 10 , b
line (51 , 106) - (594 , 224) , 2 , b
line (52 , 107) - (593 , 223) , 2 , b
locate 8 , 10
if insp% = 0 then
	input; "BASE LOCATION   (NGR  ,  Locator   or   Sitename) " , baselo$
	if baselo$ = "" then end
	locate csrlin + 1 , 10
end if

call locationIn(baselo$ , he , hn , basemast, basename$)
if len(basename$) > 30 then
	for yy% = 30 to 1 step -1
		if mid$(basename$ , yy% , 1) = " " then exit for
	next yy%
	basename$ = left$(basename$ , yy%)
end if
baselabel$ = baselo$
if hn < 0 or hn > 1300 or he < rowdata%(hn \ 100 , 1) then
	print "No height data stored,  enter  as mast height"
	baseheight = 0
else
	call GetHeight(he , hn , BaseHeight)
end if
locate 9 , 10
if basemast <= 0 then input  "MAST HEIGHT       " , basemast else print "Mast height";basemast
locate 8 , 50
print using "        Base Height ####m";baseHeight
on timer(5) gosub cl30
do
	kcount% = 0
	locate 12,10
	if distlo$ = "" then
		input; "DISTANT LOCATION   [rtn] to end  ",distlo$
		locate csrlin + 1 , 10
	end if
	if distlo$ = "" then end
	distname$ = ""
	call LocationIn(distlo$ , de , dn , distmast , distname$)
	if len(distname$) > 30 then
		for yy% = 30 to 1 step -1
			if mid$(distname$ , yy% , 1) = " " then exit for
		next yy%
		distname$ = left$(distname$ , yy%)
	end if
	distlabel$ = distlo$
	locate 13 , 10
	if distmast <= 0 then input "DISTANT MAST HEIGHT " , distmast else print "Mast height";distmast
	if dn < 0 or dn > 1300 or de < rowdata%(dn \ 100 , 1) then
		Print "No height data stored,  enter  as mast height"
		distheight = 0
	else
		call GetHeight(de , dn , DistHeight)
	end if
	cls
	x = de - he
	y = dn - hn
	dist = sqr(x * x + y * y)
	if dist = -x then
		angle = pi / 2
	else
		angle = 2 * atn(y / (x + dist))
	end if
	bearing = 90 - angle * 180 / pi
	if bearing < 0 then bearing = bearing + 360
	if bearing > 360 then bearing = bearing - 360
	increment = max(dist / 400 , .5)
	alongcount% = 0
	maxheight = 0
	locate 25 , 30
	color 4
	print "Reading data  . . . . . . . ";
	color 7
	gosub calculations
	locate 25 , 30
	print space$(30);
	reg 1 , 0
	call interrupt &h33
	if dist < distbreak then reg 1 , 2 else reg 1 , 1
	call interrupt &h33
	reg 1 , 7
	reg 3 , pmap(0 , 0)
	reg 4 , pmap(alongcount% -1 , 0)
	call interrupt &h33
	reg 1 , 8
	reg 3 , pmap(0 , 1)
	reg 4 , pmap(maxheight * 1.5 , 1)
	call interrupt &h33

	fromy = baseheight + basemast
	toy  =  spotheights(alongcount% - 1) + groundheight(alongcount% - 1) + distmast
'	if dist < distbreak then line (0 , fromy) - (alongcount% , toy) , 8 ,, &hF0F0
	locate 30 , 64
	print "F1 - .BMP file";
	do
		reg 1 , 3
		call interrupt &h33
		button% = reg(2)
		mx% = reg(3)
		my% = reg(4)

		if mx% <> oldmx% or my% <> oldmy% then moved% = -1 else moved% = 0
		oldmx% = mx%
		oldmy% = my%
		a% = max(pmap(mx% , 2) , 0)
		mh = pmap (my% , 3)

      if button% = 1 then
			locngr$ = spotngr$(a% , 0)
			call LocationIn(locngr$ , me , mn , dummy, sitename$)
			me1 = me : mn1 = mn
			call NearestSite(me , mn , nearestname$)
			searchdist = sqr((me1 - me) ^ 2 + (mn1 - mn) ^ 2)
			timer on
			if searchdist < 2 then
				ns$ = "Nearest site " + left$(nearestname$ , 40)_
						 + using$(" #.#km" , searchdist)
			else
				ns$ = "No site near" + space$(68)
			end if
			locate 30 , 1
			print ns$;
		end if

		along = max(a% / alongcount% * dist , 1)
		theta1 = ((mh - groundheight(a%) - baseheight - basemast) / 1000 / along - along / 2 / re ) * 180 / pi
		theta2 = ((mh - groundheight(a%) - distheight - distmast) / 1000 / (dist - along) - (dist - along) / 2 / re) * 180 / pi
		if abs(theta1) > 9.99 then theta1 = 9.99
		if abs(theta2) > 9.99 then theta2 = 9.99

		if moved% and dist < distbreak then
			for lx% = 0 to alongcount% - 1
				pset(lx% , oldly(lx% , 0)) , oldly(lx% , 1)  'delete last line
				if lx% < a% then
					ly = lx% / a% * (mh - fromy) + fromy
				else
					ly = (alongcount% - lx%) / (alongcount% - a%) * (mh - toy) + toy
				end if
				oldly(lx% , 0) = ly
				oldly(lx% , 1) = point(lx% , ly)
				if lx% = a% then col% = 14 else col% = 6
				pset (lx% , ly) , col%
			next lx%
		end if

		locate 27 , 20
		print using "####.# km    \      \ \      \   ####m ASL";a% * increment / 2 , spotngr$(a%,0) , spotngr$(a% , 1) , mh - groundheight(a%);
		locate 23 , 5 : print using "+#.##ø";theta1
		locate 23 , 70: print using "+#.##ø";theta2

		ik$ = ucase$(inkey$)
		if ik$ = "K" then
			kcount% = (kcount% + 1) mod 4
			k = val (mid$("1.331.501.151.00" , kcount% * 4 + 1 , 4))
			cls
			locate 30 , 30
			color 10
			print "Reading data  . . . . . . . ";
			color 7
			gosub calculations
			locate 30 , 30
			print space$(30);
         locate 30 , 64
			print "F1 - .BMP file";
		end if
		if ik$ = chr$(0 , 59) then
			window
			DumpVGAToBmp
			window (-alongcount% / 5 , -maxheight / 2 ) - ( alongcount% * 1.2 , maxheight * 2)
		end if
		if ik$ = chr$(27) then exit
		locate 1 , 20 : color 14
		print "' TERRAIN '      G4JNT  1997"
		color 7
	loop
	cls
	distlo$ = ""          'to stop continuous looping
loop until insp%         'stops if ngrs in command line
end
'==========================================================================
Sub LocationIn(lo$ , e , n , mastheight , sitename$)
	lo$ = ucase$(lo$)
	q1$="AA00AA  "
	q2$="RR99XX99"
	q3$="HA000000"
	q4$="YZ999999"
	whichtype%=0:notloc%=0:notngr%=0:found%=0
	for chk%=1 to len(lo$)
		a$=mid$(lo$,chk%,1)
		if a$ < mid$(q1$,chk%,1) or a$ > mid$(q2$,chk%,1) then notloc%=-1
		if a$ < mid$(q3$,chk%,1) or a$ > mid$(q4$,chk%,1) then notngr%=-1
	next chk%
	if (notloc% and notngr%) or len(lo$) < 6 then
		found%=0
		OSITE$=lo$
		searchloop%=0
		close 3             ' In case
		open "i",3,"sites.dat"
		do
			incr searchloop%
			line input #3,readline$
			posn1%=instr(readline$,"~")
			posn2%=instr(posn1%+7,readline$,"~")
			sitename$=left$(readline$,posn1%-1)
			if posn2%=0 then
				stloc$=mid$(readline$,posn1%+1,8)
			else
				stloc$=mid$(readline$,posn1%+1,posn2%-posn1%-1)
				posnhgt% = instr(posn2% , ucase$(readline$) , "H")
				if posnhgt% > 0 then mastheight = val(mid$(readline$ , posnhgt% + 1))_
							else mastheight = 0
			end if

			stloc$=ucase$(stloc$)
			if instr(stloc$,"~")>0 then stloc$=left$(stloc$,6)  ' 6 Digit loc stored
			b%=instr(ucase$(sitename$),osite$)
			IF B%>0 THEN
				found%=-1
				exit loop
			end if
		loop until eof(3)
		close 3
		if not found%  and dir$("auxsites.dat") > ""then
			open "i",3,"auxsites.dat"
			searchloop% = 0
			do
				incr searchloop%
				line input #3,readline$
				posn1%=instr(readline$,"~")
				posn2%=instr(posn1%+7,readline$,"~")
					sitename$=left$(readline$,posn1%-1)
				if posn2%=0 then
					stloc$=mid$(readline$,posn1%+1,8)
				else
					stloc$=mid$(readline$,posn1%+1,posn2%-posn1%-1)
					posnhgt% = instr(posn2% , ucase$(readline$) , "H")
					if posnhgt% > 0 then mastheight = val(mid$(readline$ , posnhgt% + 1))_
								else mastheight = 0
				end if

				stloc$=ucase$(stloc$)
				if instr(stloc$,"~")>0 then stloc$=left$(stloc$,6)  ' 6 Digit loc stored
				b%=instr(ucase$(sitename$),osite$)
				IF B%>0 THEN
					found%=-1
					exit loop
				end if
			loop until eof(3)
		end if
		close 3

		lo$ = stloc$
		if found%=0 then
			sound 300,1
			end
		else
      	print " ";lo$
		end if
	end if

	c6% = asc(mid$(lo$ , 6 , 1))
	c5% = asc(mid$(lo$ , 5 , 1))
	if c5% >= 48 and c5% <= 57 then
		HN%=VAL(MID$(lo$,6,3))
		HE%=VAL(MID$(lo$,3,3))
		T1%=ASC(MID$(lo$,1,1))-65:IF T1%>8 THEN T1%=T1%-1
		T2%=ASC(MID$(lo$,2,1))-65:IF T2%>8 THEN T2%=T2%-1
		E = 500 * (T1% MOD 5) + 100 * (T2% MOD 5) - 1000 + HE% / 10
		N = 1900 - 500 * (T1% \ 5) - 100 * (T2% \ 5) + HN% / 10
	else
		if len(lo$) <= 6 then lo$ = lo$ + "55"
		longd=20*(asc(mid$(lo$,1,1))-74)+2*(asc(mid$(lo$,3,1))-48)_
				+(asc(mid$(lo$,5,1))-65)/12_
				+(asc(mid$(lo$,7,1))-48)/120

		latd= 10*(asc(mid$(lo$,2,1))-74)+asc(mid$(lo$,4,1))-48_
				+(asc(mid$(lo$,6,1))-65)/24_
				+(asc(mid$(lo$,8,1))-48)/240

		DR = 180/3.141592654
		LATR = LATD / DR
		LO = LONGD / DR
		Y = SQR(1 - .00667054# * SIN(LATR) * SIN(LATR))
		N = -100 + 111.0793 * (LATD - 49)
		N = N - 31.947 * SIN((LATD - 49) / DR) * COS(LATR + 49 / DR) + .971 * (2 + LONGD) ^ 2 * SIN(LATR) * COS(LATR) / Y
		E = 400 + 111.2651 * (2 + LONGD) * COS(LATR) / Y
		E = E - .0056 * (2 + LONGD) ^ 3 * (COS(LATR)) ^ 3 * ((TAN(LATR)) ^ 2 - 1)
	end if
end sub

'=========================================================================

sub ENtoNGR(e,n,ngr$)             'e,n  in kmetres
	 if e<0 or n<-100 then
		ngr$="   --   "
		exit sub
	 end if
	 e1% = e\100
	 n1% = n\100
	 c1%=83-5*(n1%\5)+(e1%\5)
	 if c1%=73 then c1%=72
	 eh%=(e1% mod 5)
	 nh%=(n1% mod 5)
	 c2%=(4-nh%)*5+eh%+65
	 IF c2%>72 then incr c2%
	 e2&=(e-e1%*100)*10
	 n2&=(n-n1%*100)*10
	 ce$=mid$(STR$(e2&),2)
	 cn$=mid$(STR$(n2&),2)
	 while len(ce$)<3
		ce$="0"+ce$
	 wend
	 while len(cn$)<3
		cn$="0"+cn$
	 wend
	 NGR$=CHR$(C1%)+CHR$(C2%)+CE$+CN$
end sub


sub GetHeight (ein , nin , height)
	shared rowdata%()
	e = ein * 1000
	n = nin * 1000
	n100% = n \ 100000       '100 km row
	e100% = e \ 100000
	n10% = n \ 10000 mod 10
	e10% = e \ 10000 mod 10

	npoint% = n \ 500 mod 20
	epoint% = e \ 500 mod 20

	estart% = rowdata%(n100% , 1)
	filebreak& = rowdata%(n100% , 0) * 1000

	km100sq% = e100% - estart%
	km10sq% = n10% + e10% * 10

	offset10& = km100sq% * 40000  +  km10sq% * 400

	fineoffset% = epoint% + 380 - 20 * npoint%

	offset& = offset10& + fineoffset%
	pointer& = offset& + filebreak&

	seek 1 , pointer&
	get$ 1 , 1 , ch1$
	h% = asc(ch1$) ^ 1.32
	height = h%
end sub


sub ENtoLoc (e , n , locator$)
	dr = 180 / 3.141593
	T1=(N+5548.79)/6371.28
	T2=2*ATN(EXP((E-400)/6389.7))
	LONG=ATN(-COS(T2)/(COS(T1)*SIN(T2)))-2/DR
	LA=SIN(T1)*SIN(T2)
	LAT=ATN(LA/SQR(1-LA*LA))
	latd = lat*dr
	longd=long*dr
	locator$=chr$(65+int((longd+180)/20))+chr$(65+int((latd+90)/10))+_
	chr$(48+int((longd+180)/2 mod 10))+chr$(48+int((latd+90) mod 10))+_
	chr$(65+int((longd+180)*12 mod 24))+chr$(65+int((latd+90)*24 mod 24))+_
	chr$(48+int((longd+180)*120 mod 10))+chr$(48+int((latd+90)*240 mod 10))
end sub

sub NearestSite(me , mn , foundname$)
	found%=0
	nearest=9999
	foundname$=""
	open "i",3,"sites.dat"
	while not eof(3)
		line input #3,l$
		posn1%=instr(l$,"~")
		posn2%=instr(posn1%+7,l$,"~")
		sitename$=left$(l$,posn1%-1)
		if posn2%=0 then
			location$=mid$(l$,posn1%+1,8)
		else
			location$=mid$(l$,posn1%+1,posn2%-posn1%-1)
		end if
		lo$=ucase$(location$)
		c6% = asc(mid$(lo$ , 6 , 1))
		c5% = asc(mid$(lo$ , 5 , 1))
		if c5% >= 48 and c5% <= 57 then
			HN%=VAL(MID$(lo$,6,3))
			HE%=VAL(MID$(lo$,3,3))
			T1%=ASC(MID$(lo$,1,1))-65:IF T1%>8 THEN T1%=T1%-1
			T2%=ASC(MID$(lo$,2,1))-65:IF T2%>8 THEN T2%=T2%-1
			E = 500 * (T1% MOD 5) + 100 * (T2% MOD 5) - 1000 + HE% / 10
			N = 1900 - 500 * (T1% \ 5) - 100 * (T2% \ 5) + HN% / 10
		else
			if len(lo$) <= 6 then lo$ = lo$ + "55"
			longd=20*(asc(mid$(lo$,1,1))-74)+2*(asc(mid$(lo$,3,1))-48)_
					+(asc(mid$(lo$,5,1))-65)/12_
					+(asc(mid$(lo$,7,1))-48)/120

			latd= 10*(asc(mid$(lo$,2,1))-74)+asc(mid$(lo$,4,1))-48_
					+(asc(mid$(lo$,6,1))-65)/24_
					+(asc(mid$(lo$,8,1))-48)/240

			DR = 180/3.141592654
			LATR = LATD / DR
			LO = LONGD / DR
			Y = SQR(1 - .00667054# * SIN(LATR) * SIN(LATR))
			N = -100 + 111.0793 * (LATD - 49)
			N = N - 31.947 * SIN((LATD - 49) / DR) * COS(LATR + 49 / DR) + .971 * (2 + LONGD) ^ 2 * SIN(LATR) * COS(LATR) / Y
			E = 400 + 111.2651 * (2 + LONGD) * COS(LATR) / Y
			E = E - .0056 * (2 + LONGD) ^ 3 * (COS(LATR)) ^ 3 * ((TAN(LATR)) ^ 2 - 1)
		end if
		dist = sqr((e - me) ^ 2 + (n - mn) ^ 2)
		if dist <= nearest then
			nearest = dist
			foundname$=sitename$
			founde = e
			foundn = n
		end if
	wend
	close 3
	me = founde
	mn = foundn               'move mouse to nearest site on return
	sound 2000,1
end sub

sub DumpVGAToBmp
	basename$ = "TERRN"
	filename$ = dir$(basename$ + "???.bmp")
	n% = val(mid$(filename$ , 6, 3))
	do
		f$ = dir$
		n% = max(val(mid$(f$ , 6, 3)) , n%)
	loop until f$ = ""
	n$ = mid$(str$(min(n% + 1, 999)) , 2)
	do while len(n$) < 3
		n$ = "0" + n$
	loop
	f$ = basename$ + n$ + ".BMP"
	locate 30 , 64 : print using"\               \";f$;
	bwidth& = 640
	bheight& = 480
	ColoursInTable& = 16
	BitsPerPixel% = 4        '1 / 4 / 8 / 24  only
	NPlanes% = 1             '1 only
	InfoSize& = 40

	ImageSize& = bwidth& * bheight& * BitsperPixel% / 8
	dataoffset& = 54 + ColoursInTable& * 4
	filesize& = dataoffset& + ImageSize&

	t = timer
	if instr(ucase$(f$) , ".BMP") = 0 then f$ = f$ + ".BMP"
	fi% = freefile
	open "b" , fi% , f$
	i1% = 0 : i2% = 0 : compression& = 0
	bm$ = "BM"
	put$ fi% , bm$
	put fi% , 2 , filesize&
	put fi% ,, i1%
	put fi% ,, i2%
	put fi% ,, DataOffset&
	put fi% ,, InfoSize&
	put fi% ,, bwidth&
	put fi% ,, bheight&
	put fi% ,, NPlanes%
	put fi% ,, BitsPerPixel%
	put fi% ,, compression&
	put fi% ,, ImageSize&
	put fi% ,, XPelsPerMetre&
	put fi% ,, YPelsPerMetre&
	put fi% ,, ColoursInTable&
	put fi% ,, ImportantColours&

	seek fi% , 54
	pa$ =       chr$(0,0,0,0,128,0,0,0)
	pa$ = pa$ + chr$(0,128,0,0,128,128,0,0)
	pa$ = pa$ + chr$(0,0,0,128,128,0,128,0)
	pa$ = pa$ + chr$(0,63,128,0,128,128,128,0)
	pa$ = pa$ + chr$(63,63,63,0,191,63,63,0)
	pa$ = pa$ + chr$(63,191,63,0,191,191,63,0)
	pa$ = pa$ + chr$(63,63,191,0,0,140,210,0)
	pa$ = pa$ + chr$(63,191,191,0,191,191,191,0)
	put$ fi% , pa$

	seek fi% , dataoffset&         'data starts at bottom left of picture
	nibble% = 0
	for y% = bheight& - 1 to 0 step -1
		dat$ = ""
		for x% = 0 to bwidth& - 2 step 2
				d% = 16 * point (x% , y%) + point (x% + 1 , y%)
				dat$ = dat$ + chr$(d%)
		next x%
      put$ fi% , dat$
	next y%
	close fi%
	sound 1750 , 2
	locate 30 , 64
	print "F1 - .BMP file";
end sub

cl30:
	locate 30 , 1
	print space$(80);
	timer off
return

Calculations:
	maxheight = 0
	maxtheta1 = -9999
	maxtheta2 = -9999
	re = 6371.4 * k
	arc = dist / re
   y0 =  1000 * sqr(re * re - dist * dist / 4)
	alongcount% = 0
	for along = 0 to dist step increment / 2
		ae = he + along * cos(angle)
		an = hn + along * sin(angle)
		call ENtoNGR(ae , an , angr$)
		call ENtoLoc(ae , an , aloc$)
		spotngr$(alongcount% , 0) = angr$
		spotngr$(alongcount% , 1) = aloc$
		groundheight(alongcount%) = 1000 * sqr(re * re - (along - dist / 2) ^ 2) - y0
		if an < 0 or an > 1300 or ae < rowdata%(an \ 100 , 1) then
			aheight = 0
		else
			call GetHeight(ae , an , aheight)
		end if
		spotheights(alongcount%) = aheight
		maxheight = max(maxheight , aheight)
		if along > increment * 2 and along < dist - increment * 2 then
			theta1 = ((aheight - baseheight - basemast) / 1000 / along - along / 2 / re ) * 180 / pi
			theta2 = ((aheight - distheight - distmast) / 1000 / (dist - along) - (dist - along) / 2 / re) * 180 / pi
			if theta1 > maxtheta1 then
				maxtheta1 = theta1
				max1% = alongcount%
			end if
			if theta2 > maxtheta2 then
				maxtheta2 = theta2
				max2% = alongcount%
			end if
		end if
		incr alongcount%
	next along
	spotheights(alongcount% - 1) = distheight
	maxheight = maxheight + 125 * dist ^ 2 / re
   palette 13 , 38
	window (-alongcount% / 5 , -maxheight / 2 ) - ( alongcount% * 1.2 , maxheight * 2)
	line (0 , 0) - (alongcount% - 1 , 0) , 8
	line - (alongcount% - 1 , maxheight) , 8
	line (0 , 0) - (0 , maxheight) , 8
	locate 29 , 25 : print using "Distance ####.# km   Bearing ####ø";dist,bearing;
	locate 25 , 9 : print baselabel$
	locate 26 , 2 : print basename$
	locate 25 , 65: print left$(distlabel$ , 15);
	locate 26 , 80 - len(distname$) : print distname$;
	locate 20 , 72: color 13 :	print "K";
	color 7 : print using " (#.##)";k;
	lastground = 0
	lastheight = baseheight
	for a% = 1 to alongcount% - 1
		line (a% - 1 , 0) - (a% , groundheight(a%)) , 8 , bf
		line (a% - 1 , lastheight) - (a% , spotheights(a%) + groundheight(a%)) , 10
		lastheight = spotheights(a%) + groundheight(a%)
	next a%
	line - (a% - 1 , distmast + spotheights(a% - 1) + groundheight(a% - 1)) , 14
	line (0 , baseheight) - (0 , baseheight + basemast) , 14 'draw masts
	line (0 , baseheight + basemast) - (max1% , groundheight(max1%) + spotheights(max1%)) , 8 ,, &haaaa
	line (alongcount% , distheight + distmast) - (max2% , groundheight(max2%) + spotheights(max2%)) , 8 ,,&haaaa
return
