$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
'defdbl a - z
def fnacs(x) = atn(sqr((1 - x * x) / x))
dim rowdata%(12 , 1) , sitenames$(1000 , 2)

pi = 3.141592653589
k = 1.33
re = 6371.4 * k
ang = 6
ReadInDBase
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 + 16
	print "Height database missing"
	end
end if
if command$ > "" then
	baselo$ = command$
	print baselo$
end if
screen 12
open "b" , 1 , f$
'==================================================================
cls
color 15
locate 7 , 20
print "VISIBLE GROUND FROM SPECIFIED LOCATION"
color 7
line (60 , 200) - (580 , 390) , 10 , b
line (61 , 201) - (579 , 389) , 2 , b

locate 14 , 10
if baselo$ = "" then input; "BASE LOCATION   (NGR ,  Locator   or  Sitename)   " , baselo$
if baselo$ = "" then end
call locationIn(baselo$ , he , hn , mastheight , basename$)
locate csrlin + 1 , 10
print"MAST HEIGHT";
locate csrlin , 60
if mastheight = 0 then
	input , mastheight
else
	print using "### m";mastheight
end if
call GetHeight(he , hn , BaseHeight)
zhorizdist = sqr(2 * re * (baseheight + mastheight) / 1000)
defaultlimit = 10 * int((zhorizdist * 1.4) / 10 + .5)
color 2
locate csrlin + 1 , 10
print using "Base Height #### m   Zero horizon at #### km";baseHeight , zhorizdist
color 7
locate csrlin + 1, 10
print "MAXIMUM RADIAL DISTANCE OUT   [rtn] for ";using$("#### km   " , defaultlimit);
input ,limit
if limit = 0 then	limit = defaultlimit

bearingincr = round(180 / pi * atn(1 / limit) ,1)   '1km
locate csrlin + 1, 10
input "LIMIT SECTOR           y / N                      ",sect$
locate csrlin , 10
if ucase$(sect$) = "Y" then
	input "BEARING RANGE,   clockwise       From  ,  To      ",bearingfrom , bearingto
	if bearingto < bearingfrom then bearingto = bearingto + 360
else
	bearingfrom = 0
	bearingto = 360
end if
cls

SetUpPalette
ChangePalette
line (0 , 58) - (100 , 100) , 15 , b
locate 5 , 2 : print " ' VIEW '"
locate 6 , 2 : print "G4JNT 1997"
window (-limit * 1.33  , -limit) - (1.33 * limit , limit)
axes
locate 4 , 68
print "Horizon angle";

for co% = 13 to 2 step -1
	th = (co% - 8) / 7
	locate co% + 5 , 72
	color 15
	print using "+#.##ø "; th;
	color co%
	print "Û";
next co%

ts = timer
for bearing = bearingfrom to bearingto + bearingincr step bearingincr
	increment = .5
	alongcount% = 0
	maxheight = 0
	maxtheta1 = -99
	angle = (90 - bearing) * pi / 180
	xmult = cos(angle)
	ymult = sin(angle)
	for along = increment to limit step increment
		ae = along * xmult
		an = along * ymult
		call GetHeight(ae + he , an + hn , aheight)
		theta1 = ( (aheight - baseheight - mastheight) / 1000 / along - along / 2 / re ) * 180 / pi
		if theta1 > maxtheta1 then
			maxtheta1 = theta1
			horizondist = along
			co% = theta1 * 7 + 8
			if co% < 2 then co% = 2
			if co% > 13 then co% = 13
			if limit > 40 then
				pset (ae , an) , co%
			else
				line (ae , an) - (ae + increment / 2 , an + increment / 2) , co% , bf
			end if
		end if
	next along
	color 15
	locate 30 , 1
	print using "Bearing ###ø   Horizon ##.##ø   at ###.# km";bearing , maxtheta1 , horizondist;
	if inkey$ = chr$(27) then exit for
next bearing
axes
ts = timer - ts
locate 30 , 1
print space$(70);
mins = ts \ 60
secs = ts mod 60
locate 30 , 1
print using "Plot time #####m ##s";mins, secs;
on timer (5) gosub cl30
timer on
locate 1 , 1
print using " Visibility from   \                            \";ucase$(basename$);
locate 2 , 1
print using "####m ASL  Mast Height ###m";baseheight , mastheight;
locate 29 , 66
print "F1  - .BMP File";
reg 1 , 0 : mouse
reg 1 , 1 : mouse
color 15
do
	ik$ = ucase$(inkey$)
	reg 1 , 3
	mouse
	button% = reg(2)
	mx = reg(3) - 320
	my = reg(4) - 240
   select case ik$
	case chr$(0 , 75)
		decr mx
	case chr$(0 , 77)
		incr mx
	case chr$(0 , 72)
		decr my
	case chr$(0 , 80)
		incr my
	end select
	reg 3 , mx + 320
	reg 4 , my + 240
	reg 1 , 4
	mouse
	me = mx / 320 * limit * 1.333333 + he
	mn = hn - my / 240 * limit
	call ENtoNGR(me , mn , ngr$)
	call getheight(me , mn , height)
	call ENtoLOC(me , mn , locator$)
	call DistBear(me , mn , he , hn , distance , bearing)
	if button% = 1 or ik$ = "N" then
		timer off
		call NearestSite(me , mn , sitename$)
		locate 30 , 1
		print space$(41);
		locate 30 , 1
		print sitename$;
		mx = (me - he) * 320 / limit / 1.33333 + 320
		my = (hn - mn) * 240 / limit + 240
		reg 1 , 4
		reg 3 , mx
		reg 4 , my
		mouse
	end if
	locate 30 , 42
	print using "\      \  \      \  ####m ###.#km ###ø";ngr$ , locator$ , height , distance , bearing;
'on timer(5) gosub cl30
	timer on
	if ik$ = "P" then changepalette
	if ik$ = chr$(0 , 59) then
		window
		DumpVGAToBmp
		window (-limit * 1.33  , -limit) - (1.33 * limit , limit)
      locate 29 , 66
		print "F1  - .BMP File";
	end if
loop until ik$ = chr$(27)
palette
end
'==========================================================================
Sub LocationIn(lo$ , e , n , mastheight , basename$)
	shared sitenames$() , maxcounter%
	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%
	basename$ = lo$
	if (notloc% and notngr%) or len(lo$) < 6 then
		found%=0
		OSITE$=lo$
		searchloop%=0
		for n% = 0 to maxcounter%
			incr searchloop%
			sitename$=sitenames$(n% , 0)
			stloc$ = ucase$(sitenames$(n% , 1))
			readline$ = ucase$(sitenames$(n% , 2))
			posnhgt% = instr(readline$ , "H")
			if posnhgt% > 0 then mastheight = val(mid$(readline$ , posnhgt% + 1))_
						else mastheight = 0
			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 for
			end if
		next n%
		lo$ = stloc$

		if found%=0 then
			sound 300,1
			print "Not stored"
			end
		else
			l% = len(sitename$)
			possp% = instr(30 , sitename$ , " ")
			if l% < 30 then possp% = 30
			basename$ = left$(sitename$ , possp%)
		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 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&
	if pointer& > 0 then
		seek 1 , pointer&
		get$ 1 , 1 , ch1$
		height = asc(ch1$) ^ 1.32
	else
		height = 0
	end if
end sub

sub axes
	shared limit , pi
	circle (0 , 0) , .1 , 1
	locate 1 , 65
	color 15
	print "10 km circles"
	for dgrid = 0 to limit step 10
		circle (0 , 0) , dgrid , 1
	next dgrid
	dlimit = dgrid - 10
	for grid = 0 to 360 step 30
		line (0 , 0) - (dlimit * cos(grid * pi / 180) , dlimit * sin(grid * pi / 180)) , 1 ,, &hAAAA
	next grid
end sub

sub mouse
	call interrupt &h33
end sub

sub ENtoNGR(e,n,ngr$)
	 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 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 DistBear (de , dn , he , hn , dist , bear)
	dx = de - he
	dy = dn - hn
	dist = sqr(dx * dx + dy * dy)
	if (dx + dist) > 0 then bear = 360 / 3.1415927 * atn(dy / (dx + dist)) else bear = 0
	bear = 90 - bear
	if bear < 0 then bear = bear + 360
end sub

sub NearestSite(me , mn , foundname$)
	shared sitenames$() , maxcounter%
	found%=0
	nearest=9999
	foundname$=""
	for n% = 0 to maxcounter%
		sitename$= sitenames$(n% , 0)
		location$ = sitenames$(n% , 1)
		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
	next n%
	close 3
	me = founde
	mn = foundn               'move mouse to nearest site on return
	sound 2000,.1
end sub

sub SetUpPalette
	for a% = 0 to 15
		palette a% , a%
	next a%
   reg 1 , &h1010
	reg 2 , 1: reg 3 , 256 * 32 + 32: reg 4 , 32 * 256: call interrupt &h10
   reg 1 , &h1010
	reg 2 , 15: reg 3 , 256 * 63 + 63: reg 4 , 63 * 256: call interrupt &h10

end sub

sub ChangePalette
	static pal%
	incr pal%
	pal% = pal% mod 5
	for co% = 2 to 13
		select case pal%
		case 0
			r% = co% * 4 + 11
			g% = co% * 4 + 11
			b% = co% * 4 + 11
		case 4
			if co% >= 9 then r% = (co% - 9) * 9
			g% = min((co%) * 7 , 63)
			if co% <= 5 then b% = co% * 12 else b% = max((10 - co%) * 10 , 0)
			if co% >= 12 then b% = (co% - 12) * 18
		case 2
			if co% > 7 then r% = min ((co% - 4) * 9 , 63) else r% = 0
			if co% > 2 then g% = min(co% * 8  , 63)
			if co% > 5 then g% = max((12 - co%) * 6 , 0)
			if co% > 13 then g% = min((co% - 13) * 32 , 63)
			b% = min(co% * 12  , 63)
			if co% > 3 then b% = max((7 - co%) * 5 , 0)
			if co% > 11 then b% = min((co% - 11) * 15 , 63)
		case 3
			if co% > 7 then g% = min ((co% - 6) * 8 , 63) else g% = 0
			if co% > 2 then r% = min(co% * 5  , 63)
			b% = min(co% * 12  , 63)
			if co% > 3 then b% = max((7 - co%) * 5 , 0)
			if co% > 11 then b% = min((co% - 11) * 15 , 63)
		case 1  'startup
			b% = 63 - co% * 3
			g% = 63
			r% = 24 + co% * 3
		end select
		reg 1 , &h1010
		reg 2 , co%
		reg 3 , 256 * g% + b%
		reg 4 , r% * 256
		call interrupt &h10
	next co%
	restore
end sub

sub DumpVGAToBmp       'Palette table assumes EGA compatibility
	basename$ = "VIEW"
	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 29 , 66
	print using "\             \";f$;
	bwidth& = 640           'previously set by   PALETTE CO%,CO%
	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
	for rgb% = 0 to 15
		reg 1 , &h1015
		reg 2 , rgb%
		call interrupt &h10
		b% = (reg(3) mod 256) * 4
		g% = (reg(3) \ 256) * 4
		r% = (reg(4) \ 256) * 4
		pa$ = chr$(b%) + chr$(g%) + chr$(r%) + chr$(0)
		put$ fi% , pa$
	next rgb%

	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
end sub

sub ReadInDBase
	shared sitenames$(), maxcounter%
	open "i" , 1 , "sites.dat"
	while not eof(1)
		line input #1 , l$
		if l$ = "" then exit loop
		if instr(l$ , "ENDEND") > 0 then exit loop
		sp1% = instr(l$ , "~")
		sitenames$(maxcounter% , 0) = left$(l$ , sp1% - 1)
		sp2% = instr(sp1% + 1 , l$ , "~")
		if sp2% > 0 then
			sitenames$(maxcounter% , 1) = ucase$(mid$(l$ , sp1% + 1 , sp2% - sp1% - 1))
			sitenames$(maxcounter% , 2) = mid$(l$ , sp2% + 1)
		else
			sitenames$(maxcounter% , 1) = ucase$(mid$(l$ , sp1% + 1))
		end if
		incr maxcounter%
	wend
	close 1
	f$ = dir$ ("Auxsites.dat")
	if f$ > "" then
		open "i" , 1 , f$
		while not eof(1)
			line input #1 , l$
			if l$ = "" then exit loop
			sp1% = instr(l$ , "~")
			sitenames$(maxcounter% , 0) = left$(l$ , sp1% - 1)
			sp2% = instr(sp1% + 1 , l$ , "~")
			if sp2% > 0 then
				sitenames$(maxcounter% , 1) = ucase$(mid$(l$ , sp1% + 1 , sp2% - sp1% - 1))
				sitenames$(maxcounter% , 2) = mid$(l$ , sp2% + 1)
			else
				sitenames$(maxcounter% , 1) = ucase$(mid$(l$ , sp1% + 1))
			end if
			incr maxcounter%
		wend
	end if
	close 1
	decr maxcounter%
end sub

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