osgb3source

'
' Program to convert between OS GB references, Latitude/Longitude
' and Bartholomew references.
'
' Lat/long are with respect to OSGB datum (49N 2W) NOT WGS84.
' Accuracy believed to be to about nearest metre.
'
'
' Change history:
'
' 27 July 98.    fixed error of about 2 metres.
' Accuracy probably about 6cm now.
' Value of TMscale refined.
' Value of minor axis adjusted.
' value of Northing&(origin) corrected
'      - now computed in setup rather than preset.
'
' 6 Jan 99
' allowed negative values of Bartholomew northings
' allow input of degrees, minutes, seconds
' corrected rounding error in printed values of degrees, minutes, seconds.
'    - dms$ routine rewritten.
'
'
'   Phil Brady, January 1999
'
DEFINT A-Z
DECLARE SUB bartref (index)
DECLARE FUNCTION dig7$ (r&)
DECLARE FUNCTION dms$ (a!, eastwest%)
DECLARE FUNCTION getdms! (LorL$, compass$)
DECLARE FUNCTION ReftoCOORDS (in$, index)
DECLARE FUNCTION OSReference$ (index)
DECLARE SUB setup ()
DECLARE SUB COORDStoLL (index)
DECLARE SUB LLtoCOORDS (index)

CONST current = 0
CONST work = -1
CONST origin = -2
CONST lastarray = 10
CONST TMscale! = .9996012717# 'was 0.9996

'the items work (-1) and origin (-2) in the following arrays
'are for internal use within the subroutines only.
'Items 0 to 10 are for user data though this example program
'only uses item 0.
DIM SHARED easting&(-2 TO lastarray), northing&(-2 TO lastarray) 'in metres
DIM SHARED longitude!(-2 TO lastarray), latitude!(-2 TO lastarray)

DIM SHARED gridletterx(26), gridlettery(26)
DIM SHARED majoraxis!, minoraxis!, rad#
DIM SHARED eccsq!, ecc2sq!, PA!, PB!, PC!

CALL setup
DO
    bad = 0
    CLS
    PRINT
    PRINT "Phil Brady's Long and Lat <-> OSGB reference demo program."
    PRINT
    PRINT " 1 to input Latitude and Longitude in decimal degrees"
    PRINT " 2 to input Latitude and Longitude in degrees, minutes, seconds"
    PRINT " 3 to input OS grid reference"
    PRINT " 4 to input a Bartholomew reference"
    PRINT " 5 TO QUIT"
    PRINT
    INPUT "Which: ", choice

    SELECT CASE choice
    CASE 1
        PRINT "Give latitude and longitude:"
        INPUT "Latitude (eg 51.6) : ", latitude!(current)
        INPUT "Longitude (eg -4.0 for 4ø west): ", longitude!(current)
        CALL LLtoCOORDS(current)
    CASE 2
        CLS
        latitude!(current) = getdms!("Latitude", "NS")
        longitude!(current) = getdms!("Longitude", "EW")
        CALL LLtoCOORDS(current)
    CASE 3
        INPUT "Give OS reference: ", os$
        IF ReftoCOORDS(os$, current) = 0 THEN
            PRINT "Invalid reference - "
            PRINT "needs 2 letters and an even number of digits"
            BEEP
            bad = 1
        ELSE
            CALL COORDStoLL(current)
        END IF
    CASE 4
        INPUT "Give Bartholomew eastings: ", easting&(current)
        INPUT "Give Bartholomew northings: ", northing&(current)
        CALL COORDStoLL(current)
    CASE 5: END
    CASE ELSE: bad = 1
    END SELECT
    PRINT
    IF bad = 0 THEN
        PRINT "OS ref= "; OSReference$(current)
        CALL bartref(current)
        PRINT USING "Lat=###.####øN   Long=###.####ø"; latitude!(current); ABS(longitude!(current));
        PRINT MID$("W E", SGN(longitude!(current)) + 2, 1)
        PRINT "Lat= "; dms$(latitude!(current), 0);
        PRINT " Long= "; dms$(longitude!(current), 1)
        PRINT
    END IF
    PRINT "Press any key to continue"
    WHILE INKEY$ = ""
    WEND
LOOP UNTIL 0
END

badosref:
PRINT
PRINT "Bad reference"
RESUME NEXT

SUB bartref (index)
'print Bartholomew reference.
'PRINT "Easting "; easting&(index)
'PRINT "Northing "; northing&(index)
PRINT "Bartholomew E:"; dig7$(easting&(index)); " N:"; dig7$(northing&(index))
END SUB

SUB COORDStoLL (index)

'compute lat & long from easting/northings
'correction constant is in degrees/metre

IF index < 0 THEN
PRINT "Bad index for COORDStoLL"
EXIT SUB
END IF

k! = 1 / (rad# * majoraxis!)

'initial wild guess
latitude!(work) = 52
longitude!(work) = -3
count = 0
DO
    CALL LLtoCOORDS(work)
    dx& = easting&(index) - easting&(work)
    dy& = northing&(index) - northing&(work)
    latitude!(work) = latitude!(work) + dy& * k!
    longitude!(work) = longitude!(work) + dx& * k! / COS(rad# * latitude!(work))
    count = count + 1
LOOP UNTIL (ABS(dx&) + ABS(dy&) < 1) OR (count > 20)
latitude!(index) = latitude!(work)
longitude!(index) = longitude!(work)
END SUB

FUNCTION dig7$ (r&)
'return integer in 7 spaces
dig7$ = RIGHT$(" " + STR$(r&), 7)
END FUNCTION

FUNCTION dms$ (a!, eastwest)

' returns a! in degrees, seconds, minutes.
' eastwest=1 gives E/W rather than N/S

b& = INT(ABS(a!) * 3600 + .5)
d = b& \ 3600
b& = b& - 3600& * d
m = b& \ 60
b& = b& - 60 * m

out$ = RIGHT$(" " + STR$(d), 3) + "ø"
out$ = out$ + RIGHT$(STR$(m), 2) + "'"
out$ = out$ + RIGHT$(STR$(b&), 2) + CHR$(34)
s = SGN(a!)
out$ = out$ + MID$("SNNWEE", s + 2 + 3 * eastwest, 1)
dms$ = out$

END FUNCTION

FUNCTION getdms! (LorL$, compass$)
PRINT
PRINT "Please give "; LorL$
DO
    INPUT "Degrees: ", d
LOOP UNTIL ABS(d) < 90
DO
    INPUT "Minutes: ", m
LOOP UNTIL (m > -1) AND (m < 60)
DO
    INPUT "Seconds: ", s
LOOP UNTIL (s > -1) AND (s < 60)
IF compass$ = "EW" THEN
    DO
        INPUT "East or West?   ", ew$
        ew$ = UCASE$(LEFT$(ew$ + " ", 1))
    LOOP UNTIL INSTR("EW", ew$)
END IF
deg! = d + m / 60 + s / 3600
IF ew$ = "W" THEN deg! = -deg!
getdms! = deg!
END FUNCTION

SUB LLtoCOORDS (index)

longi! = rad# * (longitude!(index) - longitude!(origin))
lati! = rad# * latitude!(index)
t! = TAN(lati!)
tsq! = t! * t!

L! = longi! * COS(lati!)
Lsq! = L! * L!

etasq! = COS(lati!)
etasq! = etasq! * etasq! * ecc2sq!
nu! = SIN(lati!)
nu! = SQR(1 - eccsq! * nu! * nu!)
nu! = majoraxis! / nu!

'compute easting

A7! = (61 - tsq! * (479 - 179! * tsq! + tsq! * tsq!)) / 5040!
A5! = (5 - tsq! * (18 - tsq!) + etasq! * (14 - 58 * tsq!)) / 120
A3! = (1 - tsq! + etasq!) / 6
A1! = TMscale! * nu!
easting&(index) = A1! * L! * (1 + Lsq! * (A3! + Lsq! * (A5! + A7! * Lsq!)))
easting&(index) = easting&(index)

'now do northings
m! = majoraxis! * (1 - eccsq!) * (PA! * lati! - .5 * PB! * SIN(2 * lati!) + .25 * PC! * SIN(4 * lati!))
a2! = .5 * A1! * t!
A4! = (5 - tsq! + etasq! * (9 + 4 * etasq!)) / 12
A6! = (61 - tsq! * (58 - tsq!) + etasq! * (270 - 330 * tsq!)) / 360
northing&(index) = TMscale! * m! + a2! * Lsq! * (1 + Lsq! * (A4! + A6! * Lsq!))
IF index <> origin THEN
northing&(index) = northing&(index) - northing&(origin)
easting&(index) = easting&(index) - easting&(origin)
END IF
END SUB

FUNCTION OSReference$ (index)

'returns 10 digit reference for given easting and northing in km.
' eg SO 12345 67890

ON ERROR GOTO badosref

e& = 1000000 + easting&(index)
N& = northing&(index) + 500000

i = e& \ 500000
j = N& \ 500000
L$ = MID$("VWXYZQRSTULMNOPFGHJKABCDE", i + 1 + j * 5, 1)

e& = e& MOD 500000
N& = N& MOD 500000
i = e& \ 100000
j = N& \ 100000
L$ = L$ + MID$("VWXYZQRSTULMNOPFGHJKABCDE", j * 5 + i + 1, 1)

e& = e& MOD 100000
N& = N& MOD 100000

L$ = L$ + " " + RIGHT$(STR$(INT(100000# + e&)), 5) + " " + RIGHT$(STR$(INT(100000# + N&)), 5)

OSReference$ = L$
END FUNCTION

FUNCTION ReftoCOORDS (in$, index)

'Calculate eastings and northings for a grid reference.
' eg SO 123 456
' Note our OS origin is at reference VV 0 0 ie 1000 km W and
' 500 km S of OS false origin.
' results in km!
' Return -1 if successful, 0 if faulty input.

t$ = UCASE$(in$)
FOR i = 1 TO LEN(in$)
    c$ = MID$(t$, i, 1)
    IF INSTR("1234567890", c$) THEN
        numbers$ = numbers$ + c$
    ELSEIF INSTR("ABCDEFGHJKLMNOPQRSTUVWXYZ", c$) THEN
        refLetters$ = refLetters$ + c$
    END IF
NEXT
IF LEN(refLetters$) <> 2 THEN EXIT FUNCTION
L = LEN(numbers$)
IF (L = 0) OR (L AND 1) THEN EXIT FUNCTION
L = L / 2
refEasting$ = LEFT$(numbers$, L)
refEasting$ = LEFT$(refEasting$ + "00000", 5)

refNorthing$ = MID$(numbers$, L + 1)
refNorthing$ = LEFT$(refNorthing$ + "00000", 5)

ind = ASC(refLetters$) - 65

easting&(index) = 5000! * gridletterx(ind)
northing&(index) = 5000! * gridlettery(ind)

ind = ASC(MID$(refLetters$, 2, 1)) - 65
easting&(index) = easting&(index) + gridletterx(ind) * 1000!
northing&(index) = northing&(index) + gridlettery(ind) * 1000!
easting&(index) = easting&(index) + VAL(refEasting$) - 1000000
northing&(index) = northing&(index) + VAL(refNorthing$) - 500000

ReftoCOORDS = -1
END FUNCTION

SUB setup

CLS
'set up for 1st two letters of grid ref
char = ASC("A")
FOR v = 400 TO 0 STEP -100
    FOR h = 0 TO 400 STEP 100
        gridletterx(char - 65) = h
        gridlettery(char - 65) = v
        char = char + 1
        IF char = ASC("I") THEN char = char + 1
    NEXT
NEXT

'Airey 1849 global constants
majoraxis! = 6377563.396#
'minoraxis! = 6356256.9#
minoraxis! = 6356256.91# 'more correct?

'define origin
longitude!(origin) = -2
latitude!(origin) = 49

rad# = ATN(1) / 45
majsq! = majoraxis! * majoraxis!
minsq! = minoraxis! * minoraxis!
eccsq! = (majsq! - minsq!) / majsq! 'eccentricity ^ 2
ecc4! = eccsq! * eccsq!
ecc6! = eccsq! * ecc4!
ecc2sq! = (majsq! - minsq!) / minsq! 'second eccentricity ^2
'
'meridinal distance parameters
PA! = 1 + .75 * eccsq! + .703125 * ecc4! + .068359375# * ecc6!
PB! = .75 * eccsq! + .9375 * ecc4! + 1.025390625# * ecc6!
PC! = .234375 * ecc4! + .41015625# * ecc6!

CALL LLtoCOORDS(origin)
easting&(origin) = -400000!
northing&(origin) = northing&(origin) + 100000

' replaced code 27 July 98
'easting&(origin) = -400000!
' was northing&(origin) = 5527056 now computed as 5527063

END SUB