%Program to compute SLD99 values from accepted parameters of SLD99 report. %Coordinate Transformation From WGS84 To SLD99 system; %Entered Latitudes and Longitutes value in degrees, Minutes, Seconds formatin %a text file clc clear %Reading the text file format long [PI,PM,PS,LI,LM,LS,HI]=textread('wgsLP.txt','%f %f %f %f %f %f %f'); NPoints = min(length(PI)); for i=1:NPoints; phi=(PI(i)+(PM(i)/60)+(PS(i)/3600)); lambda=(LI(i)+(LM(i)/60)+(LS(i)/3600)); h=HI(i); % Constants for WGS84 ellipsoid a = 6378137.0; % semi-major axis in meters f = 1/298.257223563; % flattening e2 = (2*f-f^2); % square of first eccentricity b = a * (1 - f); % semi-minor axis

% Transformation parameters from WGS84 to Kandawala
deltaX =   0.2933; % Change in X (m)
deltaY = -766.9499; % Change in Y (m)
deltaZ =-87.713; % Change in Z (m)
delta_a = -860.655; % Change in semi-major axis (m)
delta_f = -0.2836136808*10^(-4); % Change in flattening


% Calculate radii of curvature
RM = a * (1 - e2) / (1 - e2 * sind(phi)^2)^(3/2);
RN = a / sqrt(1 - e2 * sind(phi)^2);

% Apply the transformation formulas
delta_phi = (1 / (RM + h)) * (-deltaX * sind(phi) * cosd(lambda) - deltaY * sind(phi) * sind(lambda) + deltaZ * cosd(phi) + ...
            (RN / a) * e2 * sind(phi) * cosd(phi) * delta_a  + ...
            sind(phi) * cosd(phi) * ((RN * b) / a + (RM * a) / b) * delta_f);
delta_lambda = (1 / ((RN + h) * cosd(phi))) * (-deltaX * sind(lambda) + deltaY * cosd(lambda));
delta_h = deltaX * cosd(phi) * cosd(lambda) + deltaY * cosd(phi) * sind(lambda) + deltaZ * sind(phi) - ...
          (a / RN) * delta_a + (RN * b) * (sind(phi)^2 / a) * delta_f;

% Calculate new coordinates in Kandawala
phi_kandawala = phi + rad2deg(delta_phi); % Convert back to degrees
lambda_kandawala = lambda + rad2deg(delta_lambda); % Convert back to degrees
h_kandawala = h + delta_h;




fprintf('\n %f %f %.3f  ',phi_kandawala,lambda_kandawala,h_kandawala);

end