g高斯反算及求hm

function[myll,mybb,mygm]=gsfs(myzyzwx,mycza,mybl,myxx,myyy)
jdn=floor(myyy/1000000);
mydzb=(1-mybl)*mycza;mydypxl=1-mydzb*mydzb/mycza/mycza;mydepxl=mycza*mycza/mydzb/mydzb-1;%%计算短轴b和第一第二偏心率e和e'
m0=mycza*(1-mydypxl);%%%??????
m0=mycza*(1-mydypxl);m2=1.5*mydypxl*m0;m4=1.25*mydypxl*m2;m6=7.0/6.0*mydypxl*m4;
m8=1.125*mydypxl*m6;
a0=m0+0.5*m2+3/8*m4+5/16*m6+35/128*m8;
a2=0.5*m2+0.5*m4+15/32*m6+7/16*m8;
a4=0.125*m4+3/16*m6+7/32*m8;
a6=m6/32+m8/16;a8=m8/128;
xxs=[a0;-0.5*a2;0.25*a4;-a6/6.0;0.125*a8];
yzb=myyy-jdn*1000000-500000;xzb=myxx;
bf0=xzb/xxs(1);
bff=xxs(2)*sin(2*bf0)+xxs(3)*sin(4*bf0)+xxs(4)*sin(6*bf0)+xxs(5)*sin(8*bf0);
bfj1=(xzb-bff)/xxs(1);
while abs(bfj1-bf0)>0.000000000001%%%%%%%%%%迭代法求大地纬度
    bf0=bfj1;
    bff=xxs(2)*sin(2*bf0)+xxs(3)*sin(4*bf0)+xxs(4)*sin(6*bf0)+xxs(5)*sin(8*bf0);
    bfj1=(xzb-bff)/xxs(1);
end
bf=bfj1;
w=sqrt(1-mydypxl*sin(bf)*sin(bf));
zwqlbj=mycza*(1-mydypxl)/w^3;myqlbj=mycza/w;%%%%求子午曲率半径和卯酉曲率半径
m1=zwqlbj;
n1=myqlbj;n3=n1*n1*n1;n5=n1*n1*n3;
t1=tan(bf);t2=t1*t1;t3=t2*t1;t4=t2*t2;
et2=mydepxl*cos(bf)*cos(bf);yyzb=yzb;
ycn=yzb/n1;
wdb=bf-t1/2*yzb/m1*yzb/n1+t1/24*yzb/m1*(yzb/n1)^3*(5+3*t2+et2-9*et2*t2)-t1/720*yzb/m1*(yzb/n1)^5*(61+90*t2+45*t4);
jdc=ycn/cos(bf)-ycn^3/(6*cos(bf))*(1+2*t2+et2)+ycn^5/120/cos(bf)*(5+28*t2+24*t4+6*et2+8*et2*t2);
zwxslj=1/n1*yzb*t1-yzb^2/(3*n3)*t1*(1+t2-et2)+yzb^5/(15*n5)*t1*(2+5*t2+3*t4);
jdl=jdc+deg2rad(myzyzwx);
mybb=rad2deg(wdb);myll=rad2deg(jdl);mygm=rad2deg(zwxslj)

%%%%运行结果
% [myll,mybb,mygm]=gsfs(117,6378245.00000000000000,0.00335232986925914,4064723.15217828,77034.839408148)
%
% mygm =
%
%           0.51542891950434
%
%
% myll =
%
%           117.862157027779
%
%
% mybb =
%
%           36.7097876111137
%
%
% mygm =
%
%           0.51542891950434

 

 

 

gszb=[4076088.839,597710.960;4063398.870,584606.600;4083934.536,606973.739;4086499.419,591151.579;4079756.943,599351.575;4077148.241,581837.429;4075127.860,590810.987;4073156.809,595793.135;4073850.585,598833.173;4071886.536,584064.934;4071351.614,588689.345;4069166.256,587651.599;4068055.790,596323.583;4068099.648,596378.712;4069251.225,598595.545];

[mm,nn]=size(gszb);

for i=1:mm

[myll,mybb,mygm]=gsfs(117,6378245.00000000000000,0.00335232986925914,gszb(i,1),gszb(i,2))

mylll(i,1)=myll;

mylll(i,2)=mybb;

end

aaa=mylll

lmax=max(aaa(:,1))

lmin=min(aaa(:,1))

bmax=max(aaa(:,2))

bmin=min(aaa(:,2))