You are here

lib-gputmcon.bas


Title: lib-gputmcon
Description: Package of GP/UTM programs

Category: Utilities/GP/UTM
Copyright: Bob Riess
Email: riessb@usa.net
Created: Nov 2002
Version: 1.0
Sbver: 0.8


'#sec:Main
SUB G2UINIT
LOCAL i,j
CONST nhz = "NPQRSTUVWXX?"
CONST shz = "MLKJHGFEDCC?"
CONST dats = 4
ndats = dats
DIM datums(1 TO dats + 1,1 TO 5)
DATA 6378206.4,294.978698000,"Clarke 1866 (NAD27)"
DATA 6378137.0,298.257222101,"WGS84/GRS80 (NAD83)"
DATA 6378388.0,297.000000000,"International 1910"
DATA 6378135.0,298.260000000,"WGS72"
FOR i = 1 TO dats
READ datums(i,1)
READ datums(i,2)
datums(i,3) = 1 / datums(i,2)
datums(i,4) = datums(i,3) * 2 - datums(i,3) ^ 2
READ datums(i,5)
NEXT
END


SUB
NDATUM(er,rf,desc)
ndats = dats + 1
datums(ndats,1) = er
datums(ndats,2) = rf
datums(ndats,3) = 1 / rf
datums(ndats,4) = datums(ndats,3) * 2 - datums(ndats,3) ^ 2
datums(ndats,5) = desc
END


SUB
GP2UTM(datnum,lat,ns,lon,ew,BYREF zone,BYREF east,BYREF north)
LOCAL phi,lam,iz,cm,fn,er,rf,f,esq,eps,pr,en,a,b,c,r,om,s,sinfi,cosfi,tn,ts,ets,l,ls,rn

IF ew = 0 THEN
ew = SGN(lon) + (lon = 0)
ENDIF
IF ew > 0 THEN
lam = 360 - ABS(lon)
iz = MAX(90 - lam \\ 6,31)
cm = RAD(543 - 6 * iz)
ELSE
lam = ABS(lon)
iz = MAX(30 - lam \\ 6,1)
cm = RAD(183 - 6 * iz)
ENDIF
lam = RAD(lam)

IF ns = 0 THEN
ns = SGN(lat) + (lat = 0)
ENDIF
phi = RAD(ABS(lat))
IF ns < 0 THEN
fn = 10000000
ELSE
fn = 0
ENDIF
zone = UTMZONE(iz,lat,ns)

er = datums(datnum,1)
rf = datums(datnum,2)
f = datums(datnum,3)
esq = datums(datnum,4)

eps = esq / (1 - esq)
pr = (1 - f) * er
en = (er - pr) / (er + pr)
r = er * (1 - en) * (1 - en ^ 2) * (1 + 2.25 * en ^ 2 + (225 / 64) * en ^ 4)
a = -1.5 * en + (9 / 16) * en ^ 3
b = 0.9375 * en ^ 2 - (15 / 32) * en ^ 4
c = -(35 / 48) * en ^ 3

om = phi + a * SIN(2 * phi) + b * SIN(4 * phi) + c * SIN(6 * phi)
s = r * om * 0.9996
sinfi = SIN(phi)
cosfi = COS(phi)
tn = sinfi / cosfi
ts = tn ^ 2
ets = eps * cosfi ^ 2
l = (lam - cm) * cosfi
ls = l * l
rn = 0.9996 * er / SQR(1 - esq * sinfi ^ 2)
a = rn * tn / 2
b = (5 - ts + ets * (9 + 4 * ets)) / 12
c = (61 + ts * (ts - 58) + ets * (270 - 330 * ts)) / 360
north = ABS(s + a * ls * (1 + ls * (b + c * ls)) - fn)
a = (1 - ts + ets) / 6
b = (5 + ts * (ts - 18) + ets * (14 - 58 * ts)) / 120
c = (61 - 479 * ts + 179 * ts ^ 2 - ts ^ 3) / 5040
east = 500000 - rn * l * (1 + ls * (a + ls * (b + c * ls)))

END


SUB
UTM2GP(datnum,BYREF lat,BYREF ns,BYREF lon,BYREF ew,BYREF zone,east,north)
LOCAL iz,cm,fn,er,rf,f,esq,eps,pr,en,c2,c4,c6,c8,v0,v2,v4,v6,r,om,cosom,foot,sinf,cosf,tn,ts,ets,rn,q,qs,l

iz = VAL(LEFT$(zone,LEN(zone) - 1))
IF iz < 30 THEN
cm = RAD(183 - 6 * iz)
ew = -1
ELSE
cm = RAD(543 - 6 * iz)
ew = 1
ENDIF

IF INSTR(shz,UCASE$(RIGHT$(zone,1))) THEN
north = 20000000 - north
ENDIF
IF north > 10000000 THEN
fn = 10000000
ns = -1
ELSE
fn = 0
ns = 1
ENDIF

er = datums(datnum,1)
rf = datums(datnum,2)
f = datums(datnum,3)
esq = datums(datnum,4)

eps = esq / (1 - esq)
pr = (1 - f) * er
en = (er - pr) / (er + pr)
r = er * (1 - en) * (1 - en ^ 2) * (1 + 2.25 * en ^ 2 + (225 / 64) * en ^ 4)
c2 = 3 * en / 2 - 27 * en ^ 3 / 32
c4 = 21 * en ^ 2 / 16 - 55 * en ^ 4 / 32
c6 = 151 * en ^ 3 / 96
c8 = 1097 * en ^ 4 / 512
v0 = 2 * (c2 - 2 * c4 + 3 * c6 - 4 * c8)
v2 = 8 * (c4 - 4 * c6 + 10 * c8)
v4 = 32 * (c6 - 6 * c8)
v6 = 128 * c8

om = (north - fn) / (r * 0.9996)
cosom = COS(om)
foot = om + SIN(om) * cosom * (v0 + v2 * cosom ^ 2 + v4 * cosom ^ 4 + v6 * cosom ^ 6)
sinf = SIN(foot)
cosf = COS(foot)
tn = sinf / cosf
ts = tn * tn
ets = eps * cosf ^ 2
rn = er * 0.9996 / SQR(1 - esq * sinf ^ 2)
q = (east - 500000) / rn
qs = q * q
c2 = -tn * (1 + ets) / 2
c4 = -(5 + 3 * ts + ets * (1 - 9 * ts) - 4 * ets ^ 2) / 12
c6 = (61 + 45 * ts * (2 + ts) + ets * (46 - 252 * ts - 60 * ts ^ 2)) / 360
lat = DEG(foot + c2 * qs * (1 + qs * (c4 + c6 * qs)))
zone = UTMZONE(iz,lat,ns)

c2 = -(1 + ts + ts + ets) / 6
c4 = (5 + ts * (28 + 24 * ts) + ets * (6 + 8 * ts)) / 120
c6 = -(61 + 662 * ts + 1320 * ts ^ 2 + 720 * ts ^ 3) / 5040
l = q * (1 + qs * (c2 + qs * (c4 + c6 * qs)))
lon = DEG(-l / cosf + cm)
IF lon > 180 THEN
lon = 360 - lon
ENDIF

END

FUNC
UTMZONE(iz,lat,ns)
IF ns < 0 THEN
UTMZONE = STR$(iz) + MID$(shz,MIN(ABS(lat) \\ 8 + 1,LEN(shz)),1)
ELSE
UTMZONE = STR$(iz) + MID$(nhz,MIN(ABS(lat) \\ 8 + 1,LEN(nhz)),1)
ENDIF
END
'