xy2latlon(x,y,lat0,​lon0,azimuth)

This is a simple transformation between local planar coordinates and longitude latitude coordinates, such as WGS84.
I did not use complicate transformation of coordinate systems of different earth model. Instead I only using coordinate system translation, rotation, and scaling. I tested using google earth and it seems the accuracy is pretty reliable. But notice this way will be inaccurate for large scale transformation.
.
The reverse transformation latlon2xy is posted here:
% From latitude longitude geographical coordinate system to local planar
% coordinate system involves translation, rotation, and scaling of
% coordinate system.
%Input lat,lon of the points;
% (lat0,lon0): the latitude longitude for the origin point of the local planar coordinate system
% azimuth: the azimuth angle of x-axis of the local coordinate system,
% y-axis is clockwise from x-axis.
%
% output x=xy(2); y=xy(1);
% (x,y): the points location of local planar coordinate system. unit in meter
%
% Notice: Not suitable for calculing large scale points. For points within 1 km,
% the accuracy depends on the azimuth angle.

% example: Set the local origin point (0,0) at (55.709264 N, 13.201449 E), a building corner
% the azimuth angle of x-axis is 104.5 deg
% the the point 55.709224,13.202661) of the geographic coordinate
% system has local coordinates at
% xy = latlon2xy(55.70922447,13.20266068,55.709264,13.201449,104.5);
%
% See also xy2latlon
% Hongxiao Jin Nov 29 2014
%
function xy = latlon2xy(lat,lon,lat0,lon0,azimuth)
xy = [];
if size(lat,1) >size(lat,2), lat=lat'; end
if size(lon,1) >size(lon,2), lon=lon'; end
if size(lat,1) ~= 1
msgbox('latitude should be a vector!');return;
end
if size(lon,1) ~= 1
msgbox('longitude should be a vector!');return;
end
if size(lat,2) ~= size(lon,2)
msgbox('Latitude Longitude should be vectors of equal length !');return;
end

R = 6371007.181; % earth radius Spheric model, same as Sinosoidal projection system
PI=3.141592653589793238;
theta=azimuth-90;
xy = [cos(PI*theta/180) sin(PI*theta/180); -sin(PI*theta/180) cos(PI*theta/180)]*[1 0;0 cos(PI*lat0/180)]*([lat;lon] - [lat0 lon0]'*ones(size(lat))) *PI*R/180;