%The 7-parameter Molodensky transformation involves converting coordinates between two different geodetic systems with slightly different reference ellipsoids using parameters for translation, rotation, and scale change. Here is the code modified to use the 7-parameter Molodensky transformation: % Program to compute Transformation parameters from SL old to SLD99 clc; clear; format long;

% Computed Latitudes and longitudes values from earlier programs are used % with orthometric heights of the points

[PI, PM, PS, LI, LM, LS, HI] = textread('known-points.txt', '%f %f %f %f %f %f %f'); [PRR, LRR, Hv] = textread('Everest.txt', '%f %f %f'); % Determine the number of points NPoints = min(length(PI), length(PRR)); % Constants for WGS84 ag = 6378137.0; bg = 6356752.3142; eg = sqrt(((ag^2) - (bg^2)) / (ag^2));

for i = 1:NPoints % Convert degree, minute, second to decimal degrees PR(i) = PI(i) + (PM(i) / 60) + (PS(i) / 3600); LR(i) = LI(i) + (LM(i) / 60) + (LS(i) / 3600);

ng(i) = (eg * sind(PR(i)))^2;
vg(i) = ag / sqrt(1 - ng(i));
XI(i) = (vg(i) + HI(i)) * cosd(PR(i)) * cosd(LR(i));
YI(i) = (vg(i) + HI(i)) * cosd(PR(i)) * sind(LR(i));
ZI(i) = (vg(i) * (1 - (eg^2)) + HI(i)) * sind(PR(i));

end

% Conversion of (Everest) Lat, lon, and ellipsoidal heights to Geocentric (X, Y, Z) coordinates av = 6377276.345; bv = 6356075.413; ev = sqrt(((av^2) - (bv^2)) / (av^2));

for i = 1:NPoints Nv = av / sqrt(1 - (ev * sind(PRR(i)))^2); Xv(i) = (Nv + Hv(i)) * cosd(PRR(i)) * cosd(LRR(i)); Yv(i) = (Nv + Hv(i)) * cosd(PRR(i)) * sind(LRR(i)); Zv(i) = (Nv * (1 - (ev^2)) + Hv(i)) * sind(PRR(i)); end

C = zeros(NPoints * 3, 7); for i = 1:NPoints C(i * 3 - 2, 1) = XI(i); C(i * 3 - 2, 3) = -ZI(i); C(i * 3 - 2, 4) = YI(i); C(i * 3 - 2, 5) = 1; C(i * 3 - 1, 1) = YI(i); C(i * 3 - 1, 2) = ZI(i); C(i * 3 - 1, 4) = -XI(i); C(i * 3 - 1, 6) = 1; C(i * 3, 1) = ZI(i); C(i * 3, 2) = -YI(i); C(i * 3, 3) = XI(i); C(i * 3, 7) = 1; end

y = zeros(NPoints * 3, 1); for i = 1:NPoints y(i * 3 - 2, 1) = Xv(i) - XI(i); y(i * 3 - 1, 1) = Yv(i) - YI(i); y(i * 3, 1) = Zv(i) - ZI(i); end

% Bursa wolf 7 parameters are computed using least square principal. S = (inv(C' * C)) * (C' * y); Rx = ((180 / pi) * 3600) * S(2); Ry = ((180 / pi) * 3600) * S(3); Rz = ((180 / pi) * 3600) * S(4); Dx = S(5); Dy = S(6); Dz = S(7); Sf = 1000000 * S(1);

disp(S); M = [Dx; Dy; Dz; Rx; Ry; Rz; Sf]; fprintf('%.6f \n%.6f \n%.6f \n%.6f \n%.6f \n%.6f \n%.6f\n', Dx, Dy, Dz, Rx, Ry, Rz, Sf);