%WGS84toTWD97
clear
WGS=input('緯度/精度：');
[i,j]=size(WGS);
latd=WGS(:,1);
lond=WGS(:,2);
a = 6378137.0;
b = 6356752.314245;
long0 = pi/180*121;
k0 = 0.9999;
dx = 250000;
lat=latd.*(pi/180);
lon=lond.*(pi/180);

e = (1-b^2/a^2)^(0.5);
e2 = e^2/(1-e^2);
n = (a-b)/(a+b);
for z=1:i
    nu(z) = a/(1-(e^2)*(sin(lat(z))^2))^0.5;
end
p = (lon-long0)';

A = a*(1 - n + (5/4.0)*(n^2 - n^3) + (81/64.0)*(n^4  - n^5));
B = (3*a*n/2.0)*(1 - n + (7/8.0)*(n^2 - n^3) + (55/64.0)*(n^4 - n^5));
C = (15*a*(n^2)/16.0)*(1 - n + (3/4.0)*(n^2 - n^3));
D = (35*a*(n^3)/48.0)*(1 - n + (11/16.0)*(n^2 - n^3));
E = (315*a*(n^4)/51.0)*(1 - n);
for z=1:i
    S(z) = A*lat(z) - B*sin(2*lat(z)) + C*sin(4*lat(z)) - D*sin(6*lat(z)) + E*sin(8*lat(z));
end
K1 = S.*k0;
for z=1:i
    K2(z) = k0*nu(z)*sin(2*lat(z))/4.0;
    K3(z) = (k0*nu(z)*sin(lat(z))*(cos(lat(z))^3)/24.0)*(5 - tan(lat(z))^2 + 9*e2*(cos(lat(z))^2) + 4*(e2^2)*(cos(lat(z))^4));
    K4(z) = k0*nu(z)*cos(lat(z));
    K5(z) = (k0*nu(z)*(cos(lat(z))^3)/6.0)*(1 - tan(lat(z))^2 + e2*(cos(lat(z))^2));
end
for z=1:i
    y(z) = K1(z) + K2(z)*(p(z)^2) + K3(z)*(p(z)^4);
    x(z) = K4(z)*p(z) + K5(z)*(p(z)^3) + dx;
end
for z=1:i
    TWD97(z,1)=x(z);
    TWD97(z,2)=y(z);
end