Problem with precision in trigonometric functions

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

Answer 1

[Jan]

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.