rgps2agd.m

来自「此功能包用于各种GPS坐标和时间的转换」· M 代码 · 共 47 行

function rgps2agd (x1,y1,z1,x2,y2,z2)

% this function transforms the input Cartesian coordinates (x1,y1,z1)
% (in the WGS84 datum)
% to the Cartesian coordinates on the AGD66/84 datum (x2,y2,z2)
% Note the 7 parameter similarity transformation is hardcoded, based on values
% found at http://www.anzlic.org.au/icsm/gdatm/chapter7.htm (medium accuracy
% AGD84 to GDA94 similarity transformation.
%
% C. Rizos
%

format long;
in=[x1 y1 z1]';
out=[x2,y2,z2];

% parameters (translation of AGD ellipsoid centre from geocentre)
Dx=117.763;
Dy=51.51;
Dz=-139.061;

T=[Dx Dy Dz]';

% 3 small rotation angles about x, y, z-axes (WGS84 to AGD84)
Rx=0.292*pi/648000; %radians
Ry=0.443*pi/648000; %radians
Rz=0.277*pi/648000; %radians
% scale difference between two datums (as ppm)
Sc=1+0.191/1000000; %scale

% set up 3x3 rotation matrix (assume small angle approximation)
R1=[1 0 0
   0 cos(Rx) sin(Rx)
   0 -sin(Rx) cos(Rx)];

R2=[cos(Ry) 0 -sin(Ry)
   0 1 0
   sin(Ry) 0 cos(Ry)];

R3=[cos(Rz) sin(Rz) 0
   -sin(Rz) cos(Rz) 0
	0 0 1];

% similarity transformation
out=T+Sc*R3*R2*R1*in;
return