Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- %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
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement