[記事] WGS84 轉 UMT (TWD97TM2)

 public double[] WGS84toTWD97TM2(double lat, double lon){
 // Ref1 : http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.htm
 // Ref2 : http://blog.ez2learn.com/2009/08/15/lat-lon-to-twd97/

 // convert from degrees to radians
 lat = Math.toRadians(lat);
 lon = Math.toRadians(lon);

 double a =  6378137.0;
 double b = 6356752.314245;
 double long0 = Math.toRadians(121);
 double k0 = 0.9999;
 double dx = 250000;

 double e = Math.sqrt(1- Math.pow(b, 2)/Math.pow(a,2));

 double e2 = Math.pow(e,2)/(1- Math.pow(e,2));

 double n = (a-b) / (a+b);

 double n2 = n*n;
 double n3 = n*n*n;
 double n4 = n*n*n*n;
 double n5 = n*n*n*n*n;

 double nu = a / Math.sqrt(1- Math.pow(e,2) * Math.pow(Math.sin(lat),2));

 double p = lon - long0;

 double A = a * (1 - n + (5/4.0)*(n2-n3) + (81/64.0)*(n4 - n5));
 double B = (3*a*n/2.0)*(1 - n + (7/8.0)*(n2 -n3) + (55/64.0)*(n4 - n5));
 double C = (15*a*n2/16.0)*(1 - n + (3/4.0)*(n2 - n3));
 double D = (35*a*n3/48.0)*(1 - n + (11/16.0)*(n2 - n3));
 double E = (315*a*n4/51.0)*(1-n);
 
 double S = A*lat - B*Math.sin(2*lat) + C*Math.sin(4*lat) - D*Math.sin(6*lat) + E*Math.sin(8*lat);

 double K1 = S*k0;
 double K2 = k0*nu*Math.sin(2*lat)/4.0;
 double K3 = (k0*nu*Math.sin(lat)*(Math.pow(Math.cos(lat),3)/24.0)) * (5 - Math.pow(Math.tan(lat), 2) + 9*e2*Math.pow(Math.cos(lat),2)) + 4*Math.pow(e2,2)*Math.pow(Math.cos(lat), 4);
 double y = K1+K2*Math.pow(p,2) + K3*Math.pow(p, 4);

 double K4 = k0*nu*Math.cos(lat);
 double K5 = (k0*nu*Math.pow(Math.cos(lat), 3)/6.0) * (1-Math.pow(Math.tan(lat),2) + e2*Math.pow(Math.cos(lat), 2));

 double x = K4*p + K5*Math.pow(p, 3) + dx;

 return new double[] {x,y};
 }