mysql_function

BEGIN
	DECLARE a FLOAT;
   DECLARE b FLOAT;
	DECLARE pi FLOAT;
   DECLARE lon0 FLOAT;
   DECLARE k0 FLOAT;
   DECLARE dx FLOAT;
   DECLARE dy FLOAT;
   DECLARE e FLOAT;
   DECLARE M FLOAT;
   DECLARE mu FLOAT;
   DECLARE e1 FLOAT;
   DECLARE J1 FLOAT;
   DECLARE J2 FLOAT;
   DECLARE J3 FLOAT;
   DECLARE J4 FLOAT;
   DECLARE fp FLOAT;
	DECLARE e2 FLOAT;
	DECLARE C1 FLOAT;
	DECLARE T1 FLOAT;
	DECLARE R1 FLOAT;
	DECLARE N1 FLOAT;
	DECLARE D FLOAT;
	DECLARE Q1 FLOAT;
	DECLARE Q2 FLOAT;
	DECLARE Q3 FLOAT;
	DECLARE Q4 FLOAT;
	DECLARE lat FLOAT;
	DECLARE Q5 FLOAT;
	DECLARE Q6 FLOAT;
	DECLARE Q7 FLOAT;
	DECLARE lon FLOAT;
	DECLARE lonlat FLOAT;
	
	
	
	set a = 6378137.0;
   set b = 6356752.314245;
	set pi = 3.14159265358979323846;
   set lon0 = 121 * pi / 180;
   set k0 = 0.9999;
   set dx = 250000;
	set dy = 0;
   set e = pow((1- pow(b,2)/pow(a,2)), 0.5);

   set x = x-dx;
   set y = y-dy;

   set M = y/k0;


         set mu = M/(a*(1.0 - pow(e, 2)/4.0 - 3*pow(e, 4)/64.0 - 5*pow(e, 6)/256.0));
         set e1 = (1.0 - pow((1.0 - pow(e, 2)), 0.5)) / (1.0 + pow((1.0 - pow(e, 2)), 0.5));

         set J1 = (3*e1/2 - 27*pow(e1, 3)/32.0);
         set J2 = (21*pow(e1, 2)/16 - 55*pow(e1, 4)/32.0);
         set J3 = (151*pow(e1, 3)/96.0);
         set J4 = (1097*pow(e1, 4)/512.0);

         set fp = mu + J1*sin(2*mu) + J2*sin(4*mu) + J3*sin(6*mu) + J4*sin(8*mu);

         set e2 = pow((e*a/b), 2);
         set C1 = pow(e2*cos(fp), 2);
         set T1 = pow(tan(fp), 2);
         set R1 = a*(1-pow(e, 2))/pow((1-pow(e, 2)*pow(sin(fp), 2)), (3.0/2.0));
         set N1 = a/pow((1-pow(e, 2)*pow(sin(fp), 2)), 0.5);

         set D = x/(N1*k0);

#        // 計算緯度
         set Q1 = N1*tan(fp)/R1;
         set Q2 = (pow(D, 2)/2.0);
         set Q3 = (5 + 3*T1 + 10*C1 - 4*pow(C1, 2) - 9*e2)*pow(D, 4)/24.0;
         set Q4 = (61 + 90*T1 + 298*C1 + 45*pow(T1, 2) - 3*pow(C1, 2) - 252*e2)*pow(D, 6)/720.0;
         set lat = fp - Q1*(Q2 - Q3 + Q4);

#        // 計算經度
         set Q5 = D;
         set Q6 = (1 + 2*T1 + C1)*pow(D, 3)/6;
         set Q7 = (5 - 2*C1 + 28*T1 - 3*pow(C1, 2) + 8*e2 + 24*pow(T1, 2))*pow(D, 5)/120.0;
         set lon = lon0 + (Q5 - Q6 + Q7)/cos(fp);
# //緯
        set lat = (lat * 180) / pi; 
# //經        
        set lon = (lon * 180) / pi; 


        #set lonlat = lon + ",";
        RETURN CONCAT(lat,",",lon);
END