Consider following matlab code,% 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; eg = 0.081819190842622; bg = 6356752.3142; % Convert degrees, minutes, seconds to decimal degrees and compute geocentric coordinates for i = 1:NPoints 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 % Constants for Everest Ellipsoid av = 6377276.345; bv = 6356075.413; ev = sqrt(((av^2) - (bv^2)) / (av^2)); % Conversion of (Everest) Lat, lon, and ellipsoidal heights to Geocentric (X, Y, Z) coordinates for i = 1:NPoints Nv(i) = av / sqrt(1 - (ev * sind(PRR(i)))^2); Xv(i) = (Nv(i) + Hv(i)) * cosd(PRR(i)) * cosd(LRR(i)); Yv(i) = (Nv(i) + Hv(i)) * cosd(PRR(i)) * sind(LRR(i)); Zv(i) = (Nv(i) * (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); % Conversion to the correct units for the rotation parameters 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); % Display the results fprintf('%.13f\n ',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); % Calculate residuals r = y - C * S; % Residual sum of squares RSS = sum(r.^2); % Number of observations n = length(y); % Number of parameters p = length(S); % Estimate of the variance of the residuals sigma2 = RSS / (n - p); % Compute the covariance matrix of the parameter estimates CovS = sigma2 * inv(C' * C); % Standard deviations of the parameter estimates are the square roots of the diagonal elements of the covariance matrix stdS = sqrt(diag(CovS)); % Display the standard deviations disp('Standard deviations of the parameters:'); disp(stdS); % Print each one with labels fprintf('Standard deviation of Scale Factor: %.6f\n', stdS(1)); fprintf('Standard deviation of Rx: %.6f\n', stdS(2) * ((180 / pi) * 3600)); fprintf('Standard deviation of Ry: %.6f\n', stdS(3) * ((180 / pi) * 3600)); fprintf('Standard deviation of Rz: %.6f\n', stdS(4) * ((180 / pi) * 3600)); fprintf('Standard deviation of Dx: %.6f\n', stdS(5)); fprintf('Standard deviation of Dy: %.6f\n', stdS(6)); fprintf('Standard deviation of Dz: %.6f\n', stdS(7));