agd2utm.m

此功能包用于各种GPS坐标和时间的转换

function agd2utm (lat,long,E,N)

% this function converts ellipsoidal coordinates
% to conformal transverse Mercator coordinates
% using formulas given in GDA manual
% input latitude and longitude in radians
% output East, North in metres


% deg conversion
longitude=(long/(2*pi))*360;

% find zone (6 degree zones, 60 total)
% zone 1 has 177 degrees west longitude
% as central meridian
% ie -177 degrees
zone=fix((longitude+180)/6+1);
% long0 is the zone central meridian
% convert to radians
long0=(((zone-1)*6)-177)*2*pi/360;

% set up constants
% ANS reference ellipsoid
invflat=298.25;
a=6378160;
% f is the flatening
f=1/invflat;
esq=2*f-f^2;
e=sqrt(esq);
% V is radius of curvature in meridian
% normal section on ellipsoid
V=a/(sqrt(1-esq*(sin(lat))^2));
% P is radius of curvature in prime
% vertical normal section on ellipsoid
P=a*(1-esq)/((1-esq*(sin(lat))^2)^(3/2));
% W (upper case)
W=V/P;

T=tan(lat);
% K0 is the central scale factor
K0=0.9996;
% w (lower case) is the change in longitude
% long0 (long-zero) is the origin longitude
w=long-long0;
% false easting and northing
fEast=500000;
fNorth=10000000;
% M meridian distance
A0=1-esq/4-(3/64)*(e^4)-(5/256)*(e^6);
A2=(3/8)*(esq+(e^4)/4+(15/128)*(e^6));
A4=(15/256)*(e^4+(3/4)*(e^6));
A6=(35/3072)*(e^6);
M=a*(A0*lat-A2*sin(2*lat)+A4*sin(4*lat)-A6*sin(6*lat));

%*********************
%* main calculations *
%*********************

% term calcs first
t1=((w^2)/6)*((cos(lat))^2)*(W-(T^2));
t2=((w^4)/120)*((cos(lat))^4)*((4*(W^3))*(1-6*(T^2))+(W^2)*(1+8*(T^2))-W*2*(T^2)+T^4);
t3=((w^6)/5040)*((cos(lat))^6)*(61-479*(T^2)+179*(T^4)-T^6);

East=(K0*V*w*cos(lat))*(1+t1+t2+t3)+fEast;
%*****************************************

% next bunch of terms
t4=((w^2)/2)*V*sin(lat)*cos(lat);
t5=((w^4)/24)*V*sin(lat)*((cos(lat))^3)*(4*(W^2)+W-T^2);
t6=((w^6)/720)*V*sin(lat)*(cos(lat)^5)*(8*(W^4)*(11-24*(T^2))-28*(W^3)*(1-6*(T^2))+(W^2)*(1-32*(T^2))-W*(2*(T^2))+T^4);
t7=((w^8)/40320)*V*sin(lat)*((cos(lat))^7)*(1385-3111*(T^2)+543*(T^4)-T^6);

North=K0*(M+t4+t5+t6+t7)+fNorth;
%*******************************************
EN=[East North];
return