$lib graph off
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)
pi = 3.141592653589
k = 1.33
re = 6371.4 * k

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

open "b" , 1 , f$
'==================================================================
cls
locate 1 , 10
color 11
print "TROPOSCATTER PATH LOSS PREDICTIONS FROM HEIGHT DATABASE"
color 3
locate 4 , 15
input; "Base Location  NGR  or  Locator    " , baselo$
locate 4 , 50
color 11
print ucase$(baselo$);
call LocationIn(baselo$ , he , hn , basemast)
print " ";baselo$;
call GetHeight(he , hn , BaseHeight)
color 3
print using " #### m ASL";baseHeight
locate 5 , 15
if basemast = 0 then
	input; "Base mast height                   ",basemast
else
	print  "Base mast height                   ";basemast;
end if
print "m"
if basemast = 0 then basemast = 1
locate 7 , 15
input; "Base station antenna beamwidth     " , bw1
print "ø"
locate 9 , 15
input "Frequency   (MHz)                  " , freq
color 7
locate 11 , 1
print string$(80 , 205);
do
	for y% = 13 to 24
		locate y% , 1
		print space$(79)
	next y%
	color 7
	locate 13 , 15
	input; "Distant Location                   ",distlo$
	if distlo$ = "" then distlo$ = olddistlo$
	olddistlo$ = distlo$
	locate 13 , 50
	print ucase$(distlo$);
	color 7
	call LocationIn(distlo$ , de , dn , farmast)
	print " ";distlo$;
	call GetHeight(de , dn , farheight)
	print; using " ####m ASL";farheight

	locate 14 , 15
	if farmast = 0 then
		input;  "Mast height                        ",farmast
	else
		print   "Mast Height                        ";farmast;
	end if
	print "m"
	oldfarmast = farmast
	locate 16 , 15
	input; "Distant station antenna beamwidth  " , bw2
	print "ø"
	if bw2 = 0 then bw2 = oldbw2
	oldbw2 = bw2
	locate 16 , 49
	print bw2
	if farmast = 0 then farmast = 1
	x = de - he
	y = dn - hn
	distance = sqr(x * x + y * y)
	if distance = -x then
		angle = pi / 2
	else
		angle = 2 * atn(y / (x + distance))
	end if
	bearing = 90 - angle * 180 / pi
	if bearing < 0 then bearing = bearing + 360
	if bearing > 360 then bearing = bearing - 360
	increment = .7
	alongcount% = 0
	maxheight = 0

'find base horizon-----------------------------
	maxneartheta = -99
	zhorizdist = sqr(2 * re * (baseheight + basemast) / 1000)
	for along = increment to distance step increment
		ae = he + along * cos(angle)
		an = hn + along * sin(angle)
		call GetHeight(ae , an , aheight)
		arc = along / re          'great circle angle in radians

		call Thetas(arc , baseheight + basemast , aheight , basetheta1 , basetheta2)

		if basetheta1 > maxneartheta then
			maxneartheta = basetheta1
			basehorizondist = along
			call ENtoNGR(ae , an , wherenearhoriz$)
			basehorhgt = aheight
		end if
	next along

'find far horizon-----------------------------
	maxfartheta = -99
	zhorizdist = sqr(2 * re * (farheight + farmast) / 1000)
	for along =  distance - increment to increment step -increment
		ae = he + along * cos(angle)
		an = hn + along * sin(angle)
      call GetHeight(ae , an , aheight)
		arc = (distance - along) / re
		call Thetas(arc , farheight + farmast , aheight , fartheta1 , fartheta2)
		if fartheta1 > maxfartheta then
			maxfartheta = fartheta1
			farhorizondist = distance - along
			call ENtoNGR(ae , an , wherefarhoriz$)
			farhorhgt = aheight
		end if
	next along

'now find arc of direct path to check for LOS -----------------------------
	arc = distance / re
	call Thetas(arc , baseheight + basemast , farheight + farmast , alltheta1 , alltheta2)
	thetadirect = alltheta1

	color 14
	locate 18 , 20
	print using "Distance  ####.#km    Bearing = ####ø          ";distance , bearing
	color 7
	locate 20 , 1
	print using "Base Horiz ####m  at ###.# km  \      \    Distant ####m at ###.# km  \      \";basehorhgt , basehorizondist , wherenearhoriz$ , farhorhgt , farhorizondist , wherefarhoriz$

	scatterangle = maxneartheta + maxfartheta + arc * 180 / pi
	fsl = 32.4 + 20 * log10(freq * distance)
	meanbeamwidth = sqr(bw1 * bw2)
	lam = 2 + 2 * scatterangle / meanbeamwidth
	color 8
	print using "é1 +#.##ø   é2 +#.##ø    és #.###ø";maxneartheta , maxfartheta , scatterangle
	color 7
	if thetadirect > maxneartheta then
		color 14
		locate 24 , 17
		print using "Path is Line of sight      Path loss  ### dB"; fsl
	else
		ls = 21 + 10 * scatterangle + 10 * log10(freq)
		ltot=fsl+ls+lam
		print using "Scattering loss ###.# db   Coupling Loss ###.# dB"; ls , lam
		color 14
		locate 24 , 17
		print using "Total Troposcatter path loss    ### dB"; ltot
	end if
	while not instat:wend
	if inkey$ = chr$(27) then end
loop
end



'==========================================================================
Sub LocationIn(lo$ , e , n , mastheight)
	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
		lo$ = stloc$
		if found%=0 then	sound 300,1 : end
	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$
	height = asc(ch1$) ^ 1.32
end sub


Sub Thetas(arc , baseheight , aheight , theta1 , theta2)
	shared pi , re
	aa = sqr( (baseheight  / 1000 + re) ^ 2 + (aheight / 1000 + re) ^ 2 - 2 * ((baseheight + basemast) / 1000 + re) * (aheight / 1000 + re) * cos(arc))
	theta1 = 180 / pi * fnacs ((aheight / 1000 + re) / aa * sin(arc))
	theta2 = 180 / pi * fnacs ((baseheight / 1000 + re) / aa * sin(arc))
	if theta1 < arc * 180 / pi and theta2 < arc * 180 / pi then
		theta1 = -theta1
		theta2 = -theta2
	else
		if aheight > baseheight then theta2 = -theta2
		if baseheight > aheight then theta1 = -theta1
	end if
end sub