高斯投影正反算matlab程序

高斯-克吕格投影是由德国数学家、物理学家、天文学家高斯于19 世纪20 年代拟定，后经德国大地测量学家克吕格于1912 年对投影公式加以补充，故称为高斯-克吕格投影，又名"等角横切椭圆柱投影”，是地球椭球面和平面间正形投影的一种。 程序绝对可靠

function [BL] = Gaussion_inverse_calculation（xyL0）
%WGS-84椭球参数
a = 6378137;  %长半轴
b = 6356752.3142;  %短半轴
%f = 1/298.257223563;  %扁率
%b = a*（1-f）;
e = （sqrt（a^2-b^2））/a;  %椭球第一偏心率
e_ = （sqrt（a^2-b^2））/b;  %椭球第二偏心率

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m0 = a*（1-e^2）;
m2 = （3*e^2*m0）/2;
m4 = （5*e^2*m2）/4;
m6 = （7*e^2*m4）/6;
m8 = （9*e^2*m6）/8;
a0 = m0+m2/2+3*m4/8+5*m6/16+35*m8/128;
a2 = m2/2+m4/2+15*m6/32+7*m8/16;
a4 = m4/8+3*m6/16+7*m8/32;
a6 = m6/32+m8/16;
a8 = m8/128;

Bf = x/a0;%初始化的底点纬度
i = 0;
%迭代法求底点纬度
format long;
while 1
    New_Bf = （x-（-a2*sin（2*Bf）/2+a4*sin（4*Bf）/4-a6*sin（6*Bf）/6+a8*sin（8*Bf）/8））./a0;
        disp（Bf）;
        disp（New_Bf）;
        disp（i）;
    if （New_Bf-Bf）<=2.777778*10^-9
        break;
    end
    Bf= New_Bf;
    i = i+1;
end
t_f = tan（New_Bf）;
n_f = e_*cos（New_Bf）;
M_f = a*（1-e^2）/（sqrt（1-e^2*sin（New_Bf）^2））^3;
N_f = a/sqrt（1-e^2*sin（New_Bf）^2）;
%维度
B = New_Bf-t_f*y^2/（2*M_f*N_f）+...
    t_f^3*（5+3*t_f^2+n_f^2-9*n_f^2*t_f^2）*y^4/（24*M_f*N_f^3）;
%-t_f*（61+90*t_f^2+45*t_f^4）*y^6/（720*M_f*N_f^5）;
B = B*180/pi;
%经差
l = y/（N_f*cos（New_Bf））-...
    （1+2*t_f^2+n_f^2）*y^3/（6*N_f^3*cos（New_Bf））+...
    （5+28*t_f^2+24*t_f^4）*y^5/（120*N_f^5*cos（New_Bf））;
L = l*180/pi+L0;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        1329  2019-03-01 16:00  测绘程序\Gaussion_inverse_calculation.m
     文件        1065  2019-03-01 16:35  测绘程序\Gaussion_normal_calculation.m
     文件         319  2018-12-08 16:49  测绘程序\solve_decimal_angle2dms.m
     文件         245  2019-02-25 19:17  测绘程序\solve_dms2decimal_angle.m
     文件         519  2019-02-25 20:49  测绘程序\solve_meridian_length.m