Code chương trình:
'********************************************************************************************* _
Nguyen Bach Thao _
Phong Ky thuat-Cong ty Dien luc Quang Nam _
Tel: 0983130081 _
Email: thaoeqn1983@gmail.com _
*********************************************************************************************
Function Helmert_X(X, Y, Z, DX, Y_rot, Z_rot, S)
'Computed Helmert transformed X coordinate.
'Input: - _
cartesian XYZ coords (X,Y,Z), X translation (DX) all in meters ; _
Y and Z rotations in seconds of arc (Y_Rot, Z_Rot) and scale in ppm (s).
'Convert rotations to radians and ppm scale to a factor
Pi = 3.14159265358979
sfactor = S * 0.000001
RadY_Rot = (Y_rot / 3600) * (Pi / 180)
RadZ_Rot = (Z_rot / 3600) * (Pi / 180)
'Compute transformed X coord
Helmert_X = X + (X * sfactor) - (Y * RadZ_Rot) + (Z * RadY_Rot) + DX
End Function
Function Helmert_Y(X, Y, Z, DY, X_rot, Z_rot, S)
'Computed Helmert transformed Y coordinate.
'Input: - _
cartesian XYZ coords (X,Y,Z), Y translation (DY) all in meters ; _
X and Z rotations in seconds of arc (X_Rot, Z_Rot) and scale in ppm (s).
'Convert rotations to radians and ppm scale to a factor
Pi = 3.14159265358979
sfactor = S * 0.000001
RadX_Rot = (X_rot / 3600) * (Pi / 180)
RadZ_Rot = (Z_rot / 3600) * (Pi / 180)
'Compute transformed Y coord
Helmert_Y = (X * RadZ_Rot) + Y + (Y * sfactor) - (Z * RadX_Rot) + DY
End Function
Function Helmert_Z(X, Y, Z, DZ, X_rot, Y_rot, S)
'Computed Helmert transformed Z coordinate.
'Input: - _
cartesian XYZ coords (X,Y,Z), Z translation (DZ) all in meters ; _
X and Y rotations in seconds of arc (X_Rot, Y_Rot) and scale in ppm (s).
'Convert rotations to radians and ppm scale to a factor
Pi = 3.14159265358979
sfactor = S * 0.000001
RadX_Rot = (X_rot / 3600) * (Pi / 180)
RadY_Rot = (Y_rot / 3600) * (Pi / 180)
'Compute transformed Z coord
Helmert_Z = (-1 * X * RadY_Rot) + (Y * RadX_Rot) + Z + (Z * sfactor) + DZ
End Function
Function XYZ_to_Lat(X, Y, Z, a, b)
'Convert XYZ to Latitude (PHI) in Dec Degrees.
'Input: - _
XYZ cartesian coords (X,Y,Z) and ellipsoid axis dimensions (a & b), all in meters.
'THIS FUNCTION REQUIRES THE "Iterate_XYZ_to_Lat" FUNCTION
'THIS FUNCTION IS CALLED BY THE "XYZ_to_H" FUNCTION
RootXYSqr = Sqr((X ^ 2) + (Y ^ 2))
e2 = ((a ^ 2) - (b ^ 2)) / (a ^ 2)
PHI1 = Atn(Z / (RootXYSqr * (1 - e2)))
PHI = Iterate_XYZ_to_Lat(a, e2, PHI1, Z, RootXYSqr)
Pi = 3.14159265358979
XYZ_to_Lat = PHI * (180 / Pi)
End Function
Function Iterate_XYZ_to_Lat(a, e2, PHI1, Z, RootXYSqr)
'Iteratively computes Latitude (PHI).
'Input: - _
ellipsoid semi major axis (a) in meters; _
eta squared (e2); _
estimated value for latitude (PHI1) in radians; _
cartesian Z coordinate (Z) in meters; _
RootXYSqr computed from X & Y in meters.
'THIS FUNCTION IS CALLED BY THE "XYZ_to_PHI" FUNCTION
'THIS FUNCTION IS ALSO USED ON IT'S OWN IN THE _
"Projection and Transformation Calculations.xls" SPREADSHEET
V = a / (Sqr(1 - (e2 * ((Sin(PHI1)) ^ 2))))
PHI2 = Atn((Z + (e2 * V * (Sin(PHI1)))) / RootXYSqr)
Do While Abs(PHI1 - PHI2) > 0.000000001
PHI1 = PHI2
V = a / (Sqr(1 - (e2 * ((Sin(PHI1)) ^ 2))))
PHI2 = Atn((Z + (e2 * V * (Sin(PHI1)))) / RootXYSqr)
Loop
Iterate_XYZ_to_Lat = PHI2
End Function
Function XYZ_to_Long(X, Y)
'Convert XYZ to Longitude (LAM) in Dec Degrees.
'Input: - _
X and Y cartesian coords in meters.
Pi = 3.14159265358979
'tailor the output to fit the equatorial quadrant as determined by the signs of X and Y
If X >= 0 Then 'longitude is in the W90 thru 0 to E90 hemisphere
XYZ_to_Long = (Atn(Y / X)) * (180 / Pi)
End If
If X < 0 And Y >= 0 Then 'longitude is in the E90 to E180 quadrant
XYZ_to_Long = ((Atn(Y / X)) * (180 / Pi)) + 180
End If
If X < 0 And Y < 0 Then 'longitude is in the E180 to W90 quadrant
XYZ_to_Long = ((Atn(Y / X)) * (180 / Pi)) - 180
End If
End Function
Function XYZ_to_H(X, Y, Z, a, b)
'Convert XYZ to Ellipsoidal Height.
'Input: - _
XYZ cartesian coords (X,Y,Z) and ellipsoid axis dimensions (a & b), all in meters.
'REQUIRES THE "XYZ_to_Lat" FUNCTION
'Compute PHI (Dec Degrees) first
PHI = XYZ_to_Lat(X, Y, Z, a, b)
'Convert PHI radians
Pi = 3.14159265358979
RadPHI = PHI * (Pi / 180)
'Compute H
RootXYSqr = Sqr((X ^ 2) + (Y ^ 2))
e2 = ((a ^ 2) - (b ^ 2)) / (a ^ 2)
V = a / (Sqr(1 - (e2 * ((Sin(RadPHI)) ^ 2))))
H = (RootXYSqr / Cos(RadPHI)) - V
XYZ_to_H = H
End Function
Function Lat_Long_H_to_X(PHI, LAM, H, a, b)
'Convert geodetic coords lat (PHI), long (LAM) and height (H) to cartesian X coordinate.
'Input: - _
Latitude (PHI)& Longitude (LAM) both in decimal degrees; _
Ellipsoidal height (H) and ellipsoid axis dimensions (a & b) all in meters.
'Convert angle measures to radians
Pi = 3.14159265358979
RadPHI = PHI * (Pi / 180)
RadLAM = LAM * (Pi / 180)
'Compute eccentricity squared and nu
e2 = ((a ^ 2) - (b ^ 2)) / (a ^ 2)
V = a / (Sqr(1 - (e2 * ((Sin(RadPHI)) ^ 2))))
'Compute X
Lat_Long_H_to_X = (V + H) * (Cos(RadPHI)) * (Cos(RadLAM))
End Function
Function Lat_Long_H_to_Y(PHI, LAM, H, a, b)
'Convert geodetic coords lat (PHI), long (LAM) and height (H) to cartesian Y coordinate.
'Input: - _
Latitude (PHI)& Longitude (LAM) both in decimal degrees; _
Ellipsoidal height (H) and ellipsoid axis dimensions (a & b) all in meters.
'Convert angle measures to radians
Pi = 3.14159265358979
RadPHI = PHI * (Pi / 180)
RadLAM = LAM * (Pi / 180)
'Compute eccentricity squared and nu
e2 = ((a ^ 2) - (b ^ 2)) / (a ^ 2)
V = a / (Sqr(1 - (e2 * ((Sin(RadPHI)) ^ 2))))
'Compute Y
Lat_Long_H_to_Y = (V + H) * (Cos(RadPHI)) * (Sin(RadLAM))
End Function
Function Lat_H_to_Z(PHI, H, a, b)
'Convert geodetic coord components latitude (PHI) and height (H) to cartesian Z coordinate.
'Input: - _
Latitude (PHI) decimal degrees; _
Ellipsoidal height (H) and ellipsoid axis dimensions (a & b) all in meters.
'Convert angle measures to radians
Pi = 3.14159265358979
RadPHI = PHI * (Pi / 180)
'Compute eccentricity squared and nu
e2 = ((a ^ 2) - (b ^ 2)) / (a ^ 2)
V = a / (Sqr(1 - (e2 * ((Sin(RadPHI)) ^ 2))))
'Compute X
Lat_H_to_Z = ((V * (1 - e2)) + H) * (Sin(RadPHI))
End Function
Function Lat_Long_to_East(PHI, LAM, a, b, E0, f0, PHI0, LAM0)
'Project Latitude and longitude to Transverse Mercator eastings.
'Input: - _
Latitude (PHI) and Longitude (LAM) in decimal degrees; _
ellipsoid axis dimensions (a & b) in meters; _
eastings of false origin (e0) in meters; _
central meridian scale factor (f0); _
latitude (PHI0) and longitude (LAM0) of false origin in decimal degrees.
'Convert angle measures to radians
Pi = 3.14159265358979
RadPHI = PHI * (Pi / 180)
RadLAM = LAM * (Pi / 180)
RadPHI0 = PHI0 * (Pi / 180)
RadLAM0 = LAM0 * (Pi / 180)
af0 = a * f0
bf0 = b * f0
e2 = ((af0 ^ 2) - (bf0 ^ 2)) / (af0 ^ 2)
N = (af0 - bf0) / (af0 + bf0)
nu = af0 / (Sqr(1 - (e2 * ((Sin(RadPHI)) ^ 2))))
rho = (nu * (1 - e2)) / (1 - (e2 * (Sin(RadPHI)) ^ 2))
eta2 = (nu / rho) - 1
p = RadLAM - RadLAM0
IV = nu * (Cos(RadPHI))
V = (nu / 6) * ((Cos(RadPHI)) ^ 3) * ((nu / rho) - ((Tan(RadPHI) ^ 2)))
VI = (nu / 120) * ((Cos(RadPHI)) ^ 5) * (5 - (18 * ((Tan(RadPHI)) ^ 2)) + ((Tan(RadPHI)) ^ 4) + (14 * eta2) - (58 * ((Tan(RadPHI)) ^ 2) * eta2))
Lat_Long_to_East = E0 + (p * IV) + ((p ^ 3) * V) + ((p ^ 5) * VI)
End Function
Function Lat_Long_to_North(PHI, LAM, a, b, E0, N0, f0, PHI0, LAM0)
'Project Latitude and longitude to Transverse Mercator northings
'Input: - _
Latitude (PHI) and Longitude (LAM) in decimal degrees; _
ellipsoid axis dimensions (a & b) in meters; _
eastings (e0) and northings (n0) of false origin in meters; _
central meridian scale factor (f0); _
latitude (PHI0) and longitude (LAM0) of false origin in decimal degrees.
'REQUIRES THE "Marc" FUNCTION
'Convert angle measures to radians
Pi = 3.14159265358979
RadPHI = PHI * (Pi / 180)
RadLAM = LAM * (Pi / 180)
RadPHI0 = PHI0 * (Pi / 180)
RadLAM0 = LAM0 * (Pi / 180)
af0 = a * f0
bf0 = b * f0
e2 = ((af0 ^ 2) - (bf0 ^ 2)) / (af0 ^ 2)
N = (af0 - bf0) / (af0 + bf0)
nu = af0 / (Sqr(1 - (e2 * ((Sin(RadPHI)) ^ 2))))
rho = (nu * (1 - e2)) / (1 - (e2 * (Sin(RadPHI)) ^ 2))
eta2 = (nu / rho) - 1
p = RadLAM - RadLAM0
m = Marc(bf0, N, RadPHI0, RadPHI)
i = m + N0
II = (nu / 2) * (Sin(RadPHI)) * (Cos(RadPHI))
III = ((nu / 24) * (Sin(RadPHI)) * ((Cos(RadPHI)) ^ 3)) * (5 - ((Tan(RadPHI)) ^ 2) + (9 * eta2))
IIIA = ((nu / 720) * (Sin(RadPHI)) * ((Cos(RadPHI)) ^ 5)) * (61 - (58 * ((Tan(RadPHI)) ^ 2)) + ((Tan(RadPHI)) ^ 4))
Lat_Long_to_North = i + ((p ^ 2) * II) + ((p ^ 4) * III) + ((p ^ 6) * IIIA)
End Function
Function E_N_to_Lat(East, North, a, b, E0, N0, f0, PHI0, LAM0)
'Un-project Transverse Mercator eastings and northings back to latitude.
'Input: - _
eastings (East) and northings (North) in meters; _
ellipsoid axis dimensions (a & b) in meters; _
eastings (e0) and northings (n0) of false origin in meters; _
central meridian scale factor (f0) and _
latitude (PHI0) and longitude (LAM0) of false origin in decimal degrees.
'REQUIRES THE "Marc" AND "InitialLat" FUNCTIONS
'Convert angle measures to radians
Pi = 3.14159265358979
RadPHI0 = PHI0 * (Pi / 180)
RadLAM0 = LAM0 * (Pi / 180)
'Compute af0, bf0, e squared (e2), n and Et
af0 = a * f0
bf0 = b * f0
e2 = ((af0 ^ 2) - (bf0 ^ 2)) / (af0 ^ 2)
N = (af0 - bf0) / (af0 + bf0)
Et = East - E0
'Compute initial value for latitude (PHI) in radians
PHId = InitialLat(North, N0, af0, RadPHI0, N, bf0)
'Compute nu, rho and eta2 using value for PHId
nu = af0 / (Sqr(1 - (e2 * ((Sin(PHId)) ^ 2))))
rho = (nu * (1 - e2)) / (1 - (e2 * (Sin(PHId)) ^ 2))
eta2 = (nu / rho) - 1
'Compute Latitude
VII = (Tan(PHId)) / (2 * rho * nu)
VIII = ((Tan(PHId)) / (24 * rho * (nu ^ 3))) * (5 + (3 * ((Tan(PHId)) ^ 2)) + eta2 - (9 * eta2 * ((Tan(PHId)) ^ 2)))
IX = ((Tan(PHId)) / (720 * rho * (nu ^ 5))) * (61 + (90 * ((Tan(PHId)) ^ 2)) + (45 * ((Tan(PHId)) ^ 4)))
E_N_to_Lat = (180 / Pi) * (PHId - ((Et ^ 2) * VII) + ((Et ^ 4) * VIII) - ((Et ^ 6) * IX))
End Function
Function E_N_to_Long(East, North, a, b, E0, N0, f0, PHI0, LAM0)
'Un-project Transverse Mercator eastings and northings back to longitude.
'Input: - _
eastings (East) and northings (North) in meters; _
ellipsoid axis dimensions (a & b) in meters; _
eastings (e0) and northings (n0) of false origin in meters; _
central meridian scale factor (f0) and _
latitude (PHI0) and longitude (LAM0) of false origin in decimal degrees.
'REQUIRES THE "Marc" AND "InitialLat" FUNCTIONS
'Convert angle measures to radians
Pi = 3.14159265358979
RadPHI0 = PHI0 * (Pi / 180)
RadLAM0 = LAM0 * (Pi / 180)
'Compute af0, bf0, e squared (e2), n and Et
af0 = a * f0
bf0 = b * f0
e2 = ((af0 ^ 2) - (bf0 ^ 2)) / (af0 ^ 2)
N = (af0 - bf0) / (af0 + bf0)
Et = East - E0
'Compute initial value for latitude (PHI) in radians
PHId = InitialLat(North, N0, af0, RadPHI0, N, bf0)
'Compute nu, rho and eta2 using value for PHId
nu = af0 / (Sqr(1 - (e2 * ((Sin(PHId)) ^ 2))))
rho = (nu * (1 - e2)) / (1 - (e2 * (Sin(PHId)) ^ 2))
eta2 = (nu / rho) - 1
'Compute Longitude
X = ((Cos(PHId)) ^ -1) / nu
XI = (((Cos(PHId)) ^ -1) / (6 * (nu ^ 3))) * ((nu / rho) + (2 * ((Tan(PHId)) ^ 2)))
XII = (((Cos(PHId)) ^ -1) / (120 * (nu ^ 5))) * (5 + (28 * ((Tan(PHId)) ^ 2)) + (24 * ((Tan(PHId)) ^ 4)))
XIIA = (((Cos(PHId)) ^ -1) / (5040 * (nu ^ 7))) * (61 + (662 * ((Tan(PHId)) ^ 2)) + (1320 * ((Tan(PHId)) ^ 4)) + (720 * ((Tan(PHId)) ^ 6)))
E_N_to_Long = (180 / Pi) * (RadLAM0 + (Et * X) - ((Et ^ 3) * XI) + ((Et ^ 5) * XII) - ((Et ^ 7) * XIIA))
End Function
Function InitialLat(North, N0, afo, PHI0, N, bfo)
'Compute initial value for Latitude (PHI) IN RADIANS.
'Input: - _
northing of point (North) and northing of false origin (n0) in meters; _
semi major axis multiplied by central meridian scale factor (af0) in meters; _
latitude of false origin (PHI0) IN RADIANS; _
n (computed from a, b and f0) and _
ellipsoid semi major axis multiplied by central meridian scale factor (bf0) in meters.
'REQUIRES THE "Marc" FUNCTION
'THIS FUNCTION IS CALLED BY THE "E_N_to_Lat", "E_N_to_Long" and "E_N_to_C" FUNCTIONS
'THIS FUNCTION IS ALSO USED ON IT'S OWN IN THE "Projection and Transformation Calculations.xls" SPREADSHEET
'First PHI value (PHI1)
PHI1 = ((North - N0) / afo) + PHI0
'Calculate M
m = Marc(bfo, N, PHI0, PHI1)
'Calculate new PHI value (PHI2)
PHI2 = ((North - N0 - m) / afo) + PHI1
'Iterate to get final value for InitialLat
Do While Abs(North - N0 - m) > 0.00001
PHI2 = ((North - N0 - m) / afo) + PHI1
m = Marc(bfo, N, PHI0, PHI2)
PHI1 = PHI2
Loop
InitialLat = PHI2
End Function
Function Marc(bf0, N, PHI0, PHI)
'Compute meridional arc.
'Input: - _
ellipsoid semi major axis multiplied by central meridian scale factor (bf0) in meters; _
n (computed from a, b and f0); _
lat of false origin (PHI0) and initial or final latitude of point (PHI) IN RADIANS.
'THIS FUNCTION IS CALLED BY THE - _
"Lat_Long_to_North" and "InitialLat" FUNCTIONS
'THIS FUNCTION IS ALSO USED ON IT'S OWN IN THE "Projection and Transformation Calculations.xls" SPREADSHEET
Marc = bf0 * (((1 + N + ((5 / 4) * (N ^ 2)) + ((5 / 4) * (N ^ 3))) * (PHI - PHI0)) _
- (((3 * N) + (3 * (N ^ 2)) + ((21 / 8) * (N ^ 3))) * (Sin(PHI - PHI0)) * (Cos(PHI + PHI0))) _
+ ((((15 / 8) * (N ^ 2)) + ((15 / 8) * (N ^ 3))) * (Sin(2 * (PHI - PHI0))) * (Cos(2 * (PHI + PHI0)))) _
- (((35 / 24) * (N ^ 3)) * (Sin(3 * (PHI - PHI0))) * (Cos(3 * (PHI + PHI0)))))
End Function
Function Lat_Long_to_C(PHI, LAM, LAM0, a, b, f0)
'Compute convergence (in decimal degrees) from latitude and longitude
'Input: - _
latitude (PHI), longitude (LAM) and longitude (LAM0) of false origin in decimal degrees; _
ellipsoid axis dimensions (a & b) in meters; _
central meridian scale factor (f0).
'Convert angle measures to radians
Pi = 3.14159265358979
RadPHI = PHI * (Pi / 180)
RadLAM = LAM * (Pi / 180)
RadLAM0 = LAM0 * (Pi / 180)
'Compute af0, bf0 and e squared (e2)
af0 = a * f0
bf0 = b * f0
e2 = ((af0 ^ 2) - (bf0 ^ 2)) / (af0 ^ 2)
'Compute nu, rho, eta2 and p
nu = af0 / (Sqr(1 - (e2 * ((Sin(RadPHI)) ^ 2))))
rho = (nu * (1 - e2)) / (1 - (e2 * (Sin(RadPHI)) ^ 2))
eta2 = (nu / rho) - 1
p = RadLAM - RadLAM0
'Compute Convergence
XIII = Sin(RadPHI)
XIV = ((Sin(RadPHI) * ((Cos(RadPHI)) ^ 2)) / 3) * (1 + (3 * eta2) + (2 * (eta2 ^ 2)))
XV = ((Sin(RadPHI) * ((Cos(RadPHI)) ^ 4)) / 15) * (2 - ((Tan(RadPHI)) ^ 2))
Lat_Long_to_C = (180 / Pi) * ((p * XIII) + ((p ^ 3) * XIV) + ((p ^ 5) * XV))
End Function
Function E_N_to_C(East, North, a, b, E0, N0, f0, PHI0)
'Compute convergence (in decimal degrees) from easting and northing
'Input: - _
Eastings (East) and Northings (North) in meters; _
ellipsoid axis dimensions (a & b) in meters; _
easting (e0) and northing (n0) of true origin in meters; _
central meridian scale factor (f0); _
latitude of central meridian (PHI0) in decimal degrees.
'REQUIRES THE "Marc" AND "InitialLat" FUNCTIONS
'Convert angle measures to radians
Pi = 3.14159265358979
RadPHI0 = PHI0 * (Pi / 180)
'Compute af0, bf0, e squared (e2), n and Et
af0 = a * f0
bf0 = b * f0
e2 = ((af0 ^ 2) - (bf0 ^ 2)) / (af0 ^ 2)
N = (af0 - bf0) / (af0 + bf0)
Et = East - E0
'Compute initial value for latitude (PHI) in radians
PHId = InitialLat(North, N0, af0, RadPHI0, N, bf0)
'Compute nu, rho and eta2 using value for PHId
nu = af0 / (Sqr(1 - (e2 * ((Sin(PHId)) ^ 2))))
rho = (nu * (1 - e2)) / (1 - (e2 * (Sin(PHId)) ^ 2))
eta2 = (nu / rho) - 1
'Compute Convergence
XVI = (Tan(PHId)) / nu
XVII = ((Tan(PHId)) / (3 * (nu ^ 3))) * (1 + ((Tan(PHId)) ^ 2) - eta2 - (2 * (eta2 ^ 2)))
XVIII = ((Tan(PHId)) / (15 * (nu ^ 5))) * (2 + (5 * ((Tan(PHId)) ^ 2)) + (3 * ((Tan(PHId)) ^ 4)))
E_N_to_C = (180 / Pi) * ((Et * XVI) - ((Et ^ 3) * XVII) + ((Et ^ 5) * XVIII))
End Function
Function Lat_Long_to_LSF(PHI, LAM, LAM0, a, b, f0)
'Compute local scale factor from latitude and longitude
'Input: - _
latitude (PHI), longitude (LAM) and longitude (LAM0) of false origin in decimal degrees; _
ellipsoid axis dimensions (a & b) in meters; _
central meridian scale factor (f0).
'Convert angle measures to radians
Pi = 3.14159265358979
RadPHI = PHI * (Pi / 180)
RadLAM = LAM * (Pi / 180)
RadLAM0 = LAM0 * (Pi / 180)
'Compute af0, bf0 and e squared (e2)
af0 = a * f0
bf0 = b * f0
e2 = ((af0 ^ 2) - (bf0 ^ 2)) / (af0 ^ 2)
'Compute nu, rho, eta2 and p
nu = af0 / (Sqr(1 - (e2 * ((Sin(RadPHI)) ^ 2))))
rho = (nu * (1 - e2)) / (1 - (e2 * (Sin(RadPHI)) ^ 2))
eta2 = (nu / rho) - 1
p = RadLAM - RadLAM0
'Compute local scale factor
XIX = ((Cos(RadPHI) ^ 2) / 2) * (1 + eta2)
XX = ((Cos(RadPHI) ^ 4) / 24) * (5 - (4 * ((Tan(RadPHI)) ^ 2)) + (14 * eta2) - (28 * ((Tan(RadPHI * eta2)) ^ 2)))
Lat_Long_to_LSF = f0 * (1 + ((p ^ 2) * XIX) + ((p ^ 4) * XX))
End Function
Function E_N_to_LSF(East, North, a, b, E0, N0, f0, PHI0)
'Compute local scale factor from from easting and northing
'Input: - _
Eastings (East) and Northings (North) in meters; _
ellipsoid axis dimensions (a & b) in meters; _
easting (e0) and northing (n0) of true origin in meters; _
central meridian scale factor (f0); _
latitude of central meridian (PHI0) in decimal degrees.
'REQUIRES THE "Marc" AND "InitialLat" FUNCTIONS
'Convert angle measures to radians
Pi = 3.14159265358979
RadPHI0 = PHI0 * (Pi / 180)
'Compute af0, bf0, e squared (e2), n and Et
af0 = a * f0
bf0 = b * f0
e2 = ((af0 ^ 2) - (bf0 ^ 2)) / (af0 ^ 2)
N = (af0 - bf0) / (af0 + bf0)
Et = East - E0
'Compute initial value for latitude (PHI) in radians
PHId = InitialLat(North, N0, af0, RadPHI0, N, bf0)
'Compute nu, rho and eta2 using value for PHId
nu = af0 / (Sqr(1 - (e2 * ((Sin(PHId)) ^ 2))))
rho = (nu * (1 - e2)) / (1 - (e2 * (Sin(PHId)) ^ 2))
eta2 = (nu / rho) - 1
'Compute local scale factor
XXI = 1 / (2 * rho * nu)
XXII = (1 + (4 * eta2)) / (24 * (rho ^ 2) * (nu ^ 2))
E_N_to_LSF = f0 * (1 + ((Et ^ 2) * XXI) + ((Et ^ 4) * XXII))
End Function
Function E_N_to_t_minus_T(AtEast, AtNorth, ToEast, ToNorth, a, b, E0, N0, f0, PHI0)
'Compute (t-T) correction in decimal degrees at point (AtEast, AtNorth) to point (ToEast,ToNorth)
'Input: - _
Eastings (AtEast) and Northings (AtNorth) in meters, of point where (t-T) is being computed; _
Eastings (ToEast) and Northings (ToNorth) in meters, of point at other end of line to which (t-T) is being computed; _
ellipsoid axis dimensions (a & b) and easting & northing (e0 & n0) of true origin in meters; _
central meridian scale factor (f0); _
latitude of central meridian (PHI0) in decimal degrees.
'REQUIRES THE "Marc" AND "InitialLat" FUNCTIONS
'Convert angle measures to radians
Pi = 3.14159265358979
RadPHI0 = PHI0 * (Pi / 180)
'Compute af0, bf0, e squared (e2), n and Nm (Northing of mid point)
af0 = a * f0
bf0 = b * f0
e2 = ((af0 ^ 2) - (bf0 ^ 2)) / (af0 ^ 2)
N = (af0 - bf0) / (af0 + bf0)
Nm = (AtNorth + ToNorth) / 2
'Compute initial value for latitude (PHI) in radians
PHId = InitialLat(Nm, N0, af0, RadPHI0, N, bf0)
'Compute nu, rho and eta2 using value for PHId
nu = af0 / (Sqr(1 - (e2 * ((Sin(PHId)) ^ 2))))
rho = (nu * (1 - e2)) / (1 - (e2 * (Sin(PHId)) ^ 2))
'Compute (t-T)
XXIII = 1 / (6 * nu * rho)
E_N_to_t_minus_T = (180 / Pi) * ((2 * (AtEast - E0)) + (ToEast - E0)) * (AtNorth - ToNorth) * XXIII
End Function
Private Function TrueAzimuth(AtEast, AtNorth, ToEast, ToNorth, a, b, E0, N0, f0, PHI0)
'Compute true azimuth in decimal degrees at point (AtEast, AtNorth) to point (ToEast,ToNorth)
'Input: - _
Eastings (AtEast) and Northings (AtNorth) in meters, of point where true azimuth is being computed; _
Eastings (ToEast) and Northings (ToNorth) in meters, of point at other end of line to which true azimuth is being computed; _
ellipsoid axis dimensions (a & b) and easting & northing (e0 & n0) of true origin in meters; _
central meridian scale factor (f0); _
latitude of central meridian (PHI0) in decimal degrees.
'REQUIRES THE "Marc", "InitialLat", "E_N_to_t_minus_T" and "E_N_to_C" FUNCTIONS
'Compute eastings and northings differences
Diffe = ToEast - AtEast
Diffn = ToNorth - AtNorth
'Compute grid bearing
If Diffe = 0 Then
If Diffn < 0 Then
GridBearing = 180
Else
GridBearing = 0
End If
GoTo EndOfComputeBearing
End If
Ratio = Diffn / Diffe
Pi = 3.14159265358979
GridAngle = (180 / Pi) * Atn(Ratio)
If Diffe > 0 Then
GridBearing = 90 - GridAngle
End If
If Diffe < 0 Then
GridBearing = 270 - GridAngle
End If
EndOfComputeBearing:
'Compute convergence
Convergence = E_N_to_C(AtEast, AtNorth, a, b, E0, N0, f0, PHI0)
'Compute (t-T) correction
t_minus_T = E_N_to_t_minus_T(AtEast, AtNorth, ToEast, ToNorth, a, b, E0, N0, f0, PHI0)
'Compute initial azimuth
InitAzimuth = GridBearing + Convergence - t_minus_T
'Set TrueAzimuth >=0 and <=360
If InitAzimuth < 0 Then
TrueAzimuth = InitAzimuth + 360
ElseIf InitAzimuth > 360 Then
TrueAzimuth = InitAzimuth - 360
Else
TrueAzimuth = InitAzimuth
End If
End Function
Public Function NBT_to_VN2000_X(Latitude As Double, longitude As Double, Height As Double)
Dim a, b, F As Double
Dim KTtruc, VTruc As Double
Dim muichieu As Double
Dim E0, N0 As Double
Dim DX, DY, DZ As Double
Dim X_rot, Y_rot, Z_rot As Double
Dim S As Double
Dim X, Y, Z As Double
Dim X1, Y1, Z1 As Double
Dim B1, L1, H1 As Double
Dim VN2000X, VN2000Y, VN2000H As Double
' Lay cac thong so co ban cua VN2000
Dim slieu As Worksheet
Set slieu = ActiveWorkbook.Worksheets("So lieu")
a = slieu.Range("a")
b = slieu.Range("b")
F = slieu.Range("f")
If b = 0 Then b = a * (1 - F)
KTtruc = slieu.Range("LAM0")
vttruc = slieu.Range("PHI0")
E0 = slieu.Range("E0")
N0 = slieu.Range("N0")
muichieu = slieu.Range("F0")
DX = slieu.Range("DX")
DY = slieu.Range("DY")
DZ = slieu.Range("DZ")
X_rot = slieu.Range("X_Rotation")
Y_rot = slieu.Range("Y_Rotation")
Z_rot = slieu.Range("Z_Rotation")
S = slieu.Range("Scale")
' Lat long H to XYZ WGS84
X = Lat_Long_H_to_X(Latitude, longitude, Height, a, b)
Y = Lat_Long_H_to_Y(Latitude, longitude, Height, a, b)
Z = Lat_H_to_Z(Latitude, Height, a, b)
'XYZ to X1Y1Z1 VN2000
X1 = Helmert_X(X, Y, Z, DX, Y_rot, Z_rot, S)
Y1 = Helmert_Y(X, Y, Z, DY, X_rot, Z_rot, S)
Z1 = Helmert_Z(X, Y, Z, DZ, X_rot, Y_rot, S)
'X1Y1Z1 to Lat Long H VN2000
B1 = XYZ_to_Lat(X1, Y1, Z1, a, b)
L1 = XYZ_to_Long(X1, Y1)
H1 = XYZ_to_H(X1, Y1, Z1, a, b)
'Lat Long H to E N
VN2000Y = Lat_Long_to_East(B1, L1, a, b, E0, muichieu, vttruc, KTtruc)
VN2000X = Lat_Long_to_North(B1, L1, a, b, E0, N0, muichieu, vttruc, KTtruc)
VN2000H = H1
NBT_to_VN2000_X = VN2000X
End Function
Public Function NBT_to_VN2000_Y(Latitude As Double, longitude As Double, Height As Double)
Dim a, b, F As Double
Dim KTtruc, VTruc As Double
Dim muichieu As Double
Dim E0, N0 As Double
Dim DX, DY, DZ As Double
Dim X_rot, Y_rot, Z_rot As Double
Dim S As Double
Dim X, Y, Z As Double
Dim X1, Y1, Z1 As Double
Dim B1, L1, H1 As Double
Dim VN2000X, VN2000Y, VN2000H As Double
' Lay cac thong so co ban cua VN2000
Dim slieu As Worksheet
Set slieu = ActiveWorkbook.Worksheets("So lieu")
a = slieu.Range("a")
b = slieu.Range("b")
F = slieu.Range("f")
If b = 0 Then b = a * (1 - F)
KTtruc = slieu.Range("LAM0")
vttruc = slieu.Range("PHI0")
E0 = slieu.Range("E0")
N0 = slieu.Range("N0")
muichieu = slieu.Range("F0")
DX = slieu.Range("DX")
DY = slieu.Range("DY")
DZ = slieu.Range("DZ")
X_rot = slieu.Range("X_Rotation")
Y_rot = slieu.Range("Y_Rotation")
Z_rot = slieu.Range("Z_Rotation")
S = slieu.Range("Scale")
' Lat long H to XYZ WGS84
X = Lat_Long_H_to_X(Latitude, longitude, Height, a, b)
Y = Lat_Long_H_to_Y(Latitude, longitude, Height, a, b)
Z = Lat_H_to_Z(Latitude, Height, a, b)
'XYZ to X1Y1Z1 VN2000
X1 = Helmert_X(X, Y, Z, DX, Y_rot, Z_rot, S)
Y1 = Helmert_Y(X, Y, Z, DY, X_rot, Z_rot, S)
Z1 = Helmert_Z(X, Y, Z, DZ, X_rot, Y_rot, S)
'X1Y1Z1 to Lat Long H VN2000
B1 = XYZ_to_Lat(X1, Y1, Z1, a, b)
L1 = XYZ_to_Long(X1, Y1)
H1 = XYZ_to_H(X1, Y1, Z1, a, b)
'Lat Long H to E N
VN2000Y = Lat_Long_to_East(B1, L1, a, b, E0, muichieu, vttruc, KTtruc)
VN2000X = Lat_Long_to_North(B1, L1, a, b, E0, N0, muichieu, vttruc, KTtruc)
VN2000H = H1
NBT_to_VN2000_Y = VN2000Y
End Function
Public Function NBT_to_VN2000_Z(Latitude As Double, longitude As Double, Height As Double)
Dim a, b, F As Double
Dim KTtruc, VTruc As Double
Dim muichieu As Double
Dim E0, N0 As Double
Dim DX, DY, DZ As Double
Dim X_rot, Y_rot, Z_rot As Double
Dim S As Double
Dim X, Y, Z As Double
Dim X1, Y1, Z1 As Double
Dim B1, L1, H1 As Double
Dim VN2000X, VN2000Y, VN2000H As Double
' Lay cac thong so co ban cua VN2000
Dim slieu As Worksheet
Set slieu = ActiveWorkbook.Worksheets("So lieu")
a = slieu.Range("a")
b = slieu.Range("b")
F = slieu.Range("f")
If b = 0 Then b = a * (1 - F)
KTtruc = slieu.Range("LAM0")
vttruc = slieu.Range("PHI0")
E0 = slieu.Range("E0")
N0 = slieu.Range("N0")
muichieu = slieu.Range("F0")
DX = slieu.Range("DX")
DY = slieu.Range("DY")
DZ = slieu.Range("DZ")
X_rot = slieu.Range("X_Rotation")
Y_rot = slieu.Range("Y_Rotation")
Z_rot = slieu.Range("Z_Rotation")
S = slieu.Range("Scale")
' Lat long H to XYZ WGS84
X = Lat_Long_H_to_X(Latitude, longitude, Height, a, b)
Y = Lat_Long_H_to_Y(Latitude, longitude, Height, a, b)
Z = Lat_H_to_Z(Latitude, Height, a, b)
'XYZ to X1Y1Z1 VN2000
X1 = Helmert_X(X, Y, Z, DX, Y_rot, Z_rot, S)
Y1 = Helmert_Y(X, Y, Z, DY, X_rot, Z_rot, S)
Z1 = Helmert_Z(X, Y, Z, DZ, X_rot, Y_rot, S)
'X1Y1Z1 to Lat Long H VN2000
B1 = XYZ_to_Lat(X1, Y1, Z1, a, b)
L1 = XYZ_to_Long(X1, Y1)
H1 = XYZ_to_H(X1, Y1, Z1, a, b)
'Lat Long H to E N
VN2000Y = Lat_Long_to_East(B1, L1, a, b, E0, muichieu, vttruc, KTtruc)
VN2000X = Lat_Long_to_North(B1, L1, a, b, E0, N0, muichieu, vttruc, KTtruc)
VN2000H = H1
NBT_to_VN2000_Z = VN2000H
End Function
Public Function NBT_to_WGS84_Lat(North As Double, East As Double, Height As Double)
Dim a, b, F As Double
Dim KTtruc, VTruc As Double
Dim muichieu As Double
Dim E0, N0 As Double
Dim DX, DY, DZ As Double
Dim X_rot, Y_rot, Z_rot As Double
Dim S As Double
Dim X, Y, Z As Double
Dim X1, Y1, Z1 As Double
Dim B1, L1, H1 As Double
Dim WGS84Lat, WGS84Long, WGS84H As Double
' Lay cac thong so co ban cua VN2000
Dim slieu As Worksheet
Set slieu = ActiveWorkbook.Worksheets("So lieu")
a = slieu.Range("a")
b = slieu.Range("b")
F = slieu.Range("f")
If b = 0 Then b = a * (1 - F)
KTtruc = slieu.Range("LAM0")
vttruc = slieu.Range("PHI0")
E0 = slieu.Range("E0")
N0 = slieu.Range("N0")
muichieu = slieu.Range("F0")
DX = slieu.Range("DX")
DY = slieu.Range("DY")
DZ = slieu.Range("DZ")
X_rot = slieu.Range("X_Rotation")
Y_rot = slieu.Range("Y_Rotation")
Z_rot = slieu.Range("Z_Rotation")
S = slieu.Range("Scale")
' E N to Lat Long VN2000
B1 = E_N_to_Lat(East, North, a, b, E0, N0, muichieu, vttruc, KTtruc)
L1 = E_N_to_Long(East, North, a, b, E0, N0, muichieu, vttruc, KTtruc)
H1 = Height
' Lat long H to XYZ VN2000
X1 = Lat_Long_H_to_X(B1, L1, Height, a, b)
Y1 = Lat_Long_H_to_Y(B1, L1, Height, a, b)
Z1 = Lat_H_to_Z(B1, L1, a, b)
'X1Y1Z1 to XYZ WGS84
X = Helmert_X(X1, Y1, Z1, DX, Y_rot, Z_rot, S)
Y = Helmert_Y(X1, Y1, Z1, DY, X_rot, Z_rot, S)
Z = Helmert_Z(X1, Y1, Z1, DX, X_rot, Y_rot, S)
'XYZ to Lat Long H WGS84
WGS84Lat = XYZ_to_Lat(X, Y, Z, a, b)
WGS84Long = XYZ_to_Long(X, Y)
WGS84H = XYZ_to_H(X, Y, Z, a, b)
NBT_to_WGS84_Lat = WGS84Lat - 2.90417781181418E-03 - 0.00006
End Function
Public Function NBT_to_WGS84_Long(North As Double, East As Double, Height As Double)
Dim a, b, F As Double
Dim KTtruc, VTruc As Double
Dim muichieu As Double
Dim E0, N0 As Double
Dim DX, DY, DZ As Double
Dim X_rot, Y_rot, Z_rot As Double
Dim S As Double
Dim X, Y, Z As Double
Dim X1, Y1, Z1 As Double
Dim B1, L1, H1 As Double
Dim WGS84Lat, WGS84Long, WGS84H As Double
' Lay cac thong so co ban cua VN2000
Dim slieu As Worksheet
Set slieu = ActiveWorkbook.Worksheets("So lieu")
a = slieu.Range("a")
b = slieu.Range("b")
F = slieu.Range("f")
If b = 0 Then b = a * (1 - F)
KTtruc = slieu.Range("LAM0")
vttruc = slieu.Range("PHI0")
E0 = slieu.Range("E0")
N0 = slieu.Range("N0")
muichieu = slieu.Range("F0")
DX = slieu.Range("DX")
DY = slieu.Range("DY")
DZ = slieu.Range("DZ")
X_rot = slieu.Range("X_Rotation")
Y_rot = slieu.Range("Y_Rotation")
Z_rot = slieu.Range("Z_Rotation")
S = slieu.Range("Scale")
' E N to Lat Long VN2000
B1 = E_N_to_Lat(East, North, a, b, E0, N0, muichieu, vttruc, KTtruc)
L1 = E_N_to_Long(East, North, a, b, E0, N0, muichieu, vttruc, KTtruc)
H1 = Height
' Lat long H to XYZ VN2000
X1 = Lat_Long_H_to_X(B1, L1, Height, a, b)
Y1 = Lat_Long_H_to_Y(B1, L1, Height, a, b)
Z1 = Lat_H_to_Z(B1, L1, a, b)
'X1Y1Z1 to XYZ WGS84
X = Helmert_X(X1, Y1, Z1, DX, Y_rot, Z_rot, S)
Y = Helmert_Y(X1, Y1, Z1, DY, X_rot, Z_rot, S)
Z = Helmert_Z(X1, Y1, Z1, DX, X_rot, Y_rot, S)
'XYZ to Lat Long H WGS84
WGS84Lat = XYZ_to_Lat(X, Y, Z, a, b)
WGS84Long = XYZ_to_Long(X, Y)
WGS84H = XYZ_to_H(X, Y, Z, a, b)
NBT_to_WGS84_Long = WGS84Long + 3.76088254250817E-03 - 0.00013
End Function
Public Function NBT_to_WGS84_H(North As Double, East As Double, Height As Double)
Dim a, b, F As Double
Dim KTtruc, VTruc As Double
Dim muichieu As Double
Dim E0, N0 As Double
Dim DX, DY, DZ As Double
Dim X_rot, Y_rot, Z_rot As Double
Dim S As Double
Dim X, Y, Z As Double
Dim X1, Y1, Z1 As Double
Dim B1, L1, H1 As Double
Dim WGS84Lat, WGS84Long, WGS84H As Double
' Lay cac thong so co ban cua VN2000
Dim slieu As Worksheet
Set slieu = ActiveWorkbook.Worksheets("So lieu")
a = slieu.Range("a")
b = slieu.Range("b")
F = slieu.Range("f")
If b = 0 Then b = a * (1 - F)
KTtruc = slieu.Range("LAM0")
vttruc = slieu.Range("PHI0")
E0 = slieu.Range("E0")
N0 = slieu.Range("N0")
muichieu = slieu.Range("F0")
DX = slieu.Range("DX")
DY = slieu.Range("DY")
DZ = slieu.Range("DZ")
X_rot = slieu.Range("X_Rotation")
Y_rot = slieu.Range("Y_Rotation")
Z_rot = slieu.Range("Z_Rotation")
S = slieu.Range("Scale")
' E N to Lat Long VN2000
B1 = E_N_to_Lat(East, North, a, b, E0, N0, muichieu, vttruc, KTtruc)
L1 = E_N_to_Long(East, North, a, b, E0, N0, muichieu, vttruc, KTtruc)
H1 = Height
' Lat long H to XYZ VN2000
X1 = Lat_Long_H_to_X(B1, L1, Height, a, b)
Y1 = Lat_Long_H_to_Y(B1, L1, Height, a, b)
Z1 = Lat_H_to_Z(B1, L1, a, b)
'X1Y1Z1 to XYZ WGS84
X = Helmert_X(X1, Y1, Z1, DX, Y_rot, Z_rot, S)
Y = Helmert_Y(X1, Y1, Z1, DY, X_rot, Z_rot, S)
Z = Helmert_Z(X1, Y1, Z1, DX, X_rot, Y_rot, S)
'XYZ to Lat Long H WGS84
WGS84Lat = XYZ_to_Lat(X, Y, Z, a, b)
WGS84Long = XYZ_to_Long(X, Y)
WGS84H = XYZ_to_H(X, Y, Z, a, b)
NBT_to_WGS84_H = WGS84H - 91.4520812714472 - 4.2632564145606E-14 + 47
End Function
