Problem with precision in trigonometric functions

13 次查看（过去 30 天）
bradzik
bradzik 2016-10-26
编辑： Jan 2016-10-26
I am encountering problem with converting coordinates from geographic (GCS) to cartesian coordinate system (CCS) and vice versa. Below is my code, if you run it you will see that there is a loss of precision during two conversions. I'm starting with a point P in CCS, converting it to GCS and then back to CCS. As a result I am getting a bit different coordinates from initial point P coordinates. Moreover the result differs depending on the trigonometric functions used. The error seems small but it causes a significant difference in position calculated by trillateration. Is there any way to convert such coordinates with greater precision i Matlab?
% cartesian coordinates of point P
P = [3554.6 ; 1205.7 ; 5150.9]
% conversion to spherical geographic coordinate system method I
R = 6371;
lat1 = asin(P(3)/R); % lattitude of point P
lon1 = atan2(P(2),P(1)); % longitutde of point P
% conversion back to cartesian system for method I
P_converted1 = [R*cos(lon1) * cos(lat1); R*sin(lon1) * cos(lat1); R*sin(lat1)]
error1 = (abs(P-P_converted1)./P) * 100 % percentage error of two conversions
% conversion to spherical geographic coordinate system method II
R = 6371;
lat2 = acos(sqrt(P(1)^2 + P(2)^2)/R) % lattitude of point P
lon2 = atan2(P(2),P(1)); % longitutde of point P
% conversion back to cartesian system for method II
P_converted2 = [R*cos(lon2) * cos(lat2); R*sin(lon2) * cos(lat2); R*sin(lat2)]
error2 = (abs(P-P_converted2)./P) * 100 % percentage error of two conversions

请先登录，再回答此问题。

采纳的回答

Jan
Jan 2016-10-26
编辑：Jan 2016-10-26
The precision is reduced by the weak definition of the radius R:
P = [3554.6 ; 1205.7 ; 5150.9]
norm(P)
% >> 6373.43427517692
R = 6371; % !!!
With a correct radius, the error vanishes:
P = [3554.6 ; 1205.7 ; 5150.9]
R = norm(P);
lat1 = asin(P(3)/R); % lattitude of point P
lon1 = atan2(P(2),P(1)); % longitutde of point P
P_converted1 = [R*cos(lon1) * cos(lat1); R*sin(lon1) * cos(lat1); R*sin(lat1)]
error1 = (abs(P-P_converted1)./P) * 100
error1 = [1.27932074181754e-014, 0, 0] and the same for error2.

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Geographic Plots 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by