经纬度坐标转为上海地方坐标代码（js代码）

转   http://blog.sina.com.cn/s/blog_67fa5ecd0101aefi.html

// 经纬度坐标转换为上海地方坐标

function wgs84Tosh(lat,lon) {

    var xx, yy, r2d, tolat, tolon, rearth, PI;

    PI = 3.141592653589793;

    r2d = 57.2957795131;

    tolat = (31 + (14.0 + 7.55996 / 60.0) / 60.0) / r2d;

    tolon = (121.0 + (28.0 + 1.80651 / 60.0) / 60) / r2d;

    rearth = 6371006.84;

    var hor, frlat, frlon, gcdist,clatf, clatt, slatf, slatt,gcbrg;

    var dlon, cdlon, sdlon, sdist, cdist, sbrg, cbrg, temp; 

    var intlat, intlon;

    intlat = lat;

    intlon = lon;

    frlat = lat / r2d;

    frlon = lon / r2d;

    clatt = Math.cos(frlat);

    clatf = Math.cos(tolat);

    slatt = Math.sin(frlat);

    slatf = Math.sin(tolat);

    dlon = frlon - tolon;

    cdlon = Math.cos(dlon);

    sdlon = Math.sin(dlon);

    cdist = slatf * slatt + clatf * clatt * cdlon;

    temp = (clatt * sdlon) * (clatt * sdlon) + (clatf * slatt - slatf * clatt * cdlon) * (clatf * slatt - slatf * clatt *cdlon);

    sdist = Math.sqrt(Math.abs(temp));

    if ((Math.abs(sdist) > 1e-7) || (Math.abs(cdist) > 1e-7))

        gcdist = Math.atan2(sdist, cdist);

    else

        gcdist = 0;

    sbrg = sdlon * clatt;

    cbrg = (clatf * slatt - slatf * clatt * cdlon);

    if ((Math.abs(sbrg) > 1e-7) || (Math.abs(cbrg) > 1e-7)) {

        temp = Math.atan2(sbrg, cbrg);

        while (temp < 0) {

            temp = temp + 2 * PI;

            gcbrg=temp;

        }

    } else

        gcbrg = 0;

    hor = gcdist * rearth;

    xx = hor * Math.sin(temp);

    yy = hor * Math.cos(temp); 

    return [xx,yy];

}

c#

class CWgs84Tosh
    {
        // 经纬度坐标转换为上海地方坐标
        public double[] wgs84Tosh(double lat, double lon)
        {
            double xx, yy, r2d, tolat, tolon, rearth, PI;
            PI = 3.141592653589793;
            r2d = 57.2957795131;
            tolat = (31 + (14.0 + 7.55996 / 60.0) / 60.0) / r2d;
            tolon = (121.0 + (28.0 + 1.80651 / 60.0) / 60) / r2d;
            rearth = 6371006.84;
            double hor, frlat, frlon, gcdist, clatf, clatt, slatf, slatt, gcbrg;
            double dlon, cdlon, sdlon, sdist, cdist, sbrg, cbrg, temp;
            double intlat, intlon;
            intlat = lat;
            intlon = lon;
            frlat = lat / r2d;
            frlon = lon / r2d;
            clatt = Math.Cos(frlat);
            clatf = Math.Cos(tolat);
            slatt = Math.Sin(frlat);
            slatf = Math.Sin(tolat);
            dlon = frlon - tolon;
            cdlon = Math.Cos(dlon);
            sdlon = Math.Sin(dlon);
            cdist = slatf * slatt + clatf * clatt * cdlon;
            temp = (clatt * sdlon) * (clatt * sdlon) + (clatf * slatt - slatf * clatt * cdlon) * (clatf * slatt - slatf * clatt * cdlon);
            sdist = Math.Sqrt(Math.Abs(temp));
            if ((Math.Abs(sdist) > 1e-7) || (Math.Abs(cdist) > 1e-7))
                gcdist = Math.Atan2(sdist, cdist);
            else
                gcdist = 0;
            sbrg = sdlon * clatt;
            cbrg = (clatf * slatt - slatf * clatt * cdlon);
            if ((Math.Abs(sbrg) > 1e-7) || (Math.Abs(cbrg) > 1e-7))
            {
                temp = Math.Atan2(sbrg, cbrg);
                while (temp < 0)
                {
                    temp = temp + 2 * PI;
                    gcbrg = temp;
                }
            }
            else
                gcbrg = 0;


            hor = gcdist * rearth;
            xx = hor * Math.Sin(temp);
            yy = hor * Math.Cos(temp);
            double[] res = new double[2];
            res[0] = xx;
            res[1] = yy;
            return res;
        }
    }