地籍调查

%白塞尔大地主题正算

%B2, L2, A1, B1, L1, A2 单位 rad

%S 单位 m

function [B2, L2, A2] = BesselGeodeticDirect(B1, L1, A1, S, eps, RefEllipsoid)

a = RefEllipsoid(1);

b = RefEllipsoid(2);

c = RefEllipsoid(3);

alpha = RefEllipsoid(4);

e_2 = RefEllipsoid(5);

e_2_ = RefEllipsoid(6); %e'_2

e = e_2^0.5;

%step1 project elements from ellipsoid to sphere

tanu1 = (1-e_2)^0.5*tan(B1);

u1 = atan(tanu1);

sinA0 = cos(u1)*sin(A1);

tanr1 = tanu1*sec(A1);

r1 = atan(tanr1);

pk_2 = e_2_ * ( 1 - sinA0*sinA0 );

pk_4 = pk_2*pk_2;

pk_6 = pk_4*pk_2;

pA = (1-pk_2/4+7*pk_4/64-15*pk_6/256)/b;

pB = pk_2/4-pk_4/8+37*pk_6/512;

pC = pk_4/128-pk_6/128;

r0 = pA*S;

while(1)

r = pA*S + pB*sin(r0)*cos(2*r1+r0)+pC*sin(2*r0)*cos(4*r1+2*r0);

if(abs(r-r0)

break;

end

r0 = r;

end

%step2 solve spherical trigonometry

tanA2 = cos(u1)*sin(A1)/( cos(u1)*cos(r)*cos(A1) - sin(u1)*sin(r) ); A2 = atan(tanA2);

sinu2 = sin(u1)*cos(r)+cos(u1)*cos(A1)*sin(r);

u2 = asin(sinu2);

tanfy = sin(A1)*sin(r)/( cos(u1)*cos(r) - sin(u1)*sin(r)*cos(A1) );

fy = atan(tanfy);

%step3 remap sphere element to ellipsoid

tanB2 = (1+e_2_)^0.5*tan(u2);

B2 = atan(tanB2);

sinA1 = sin(A1);

if(sinA1>0)

if(tanfy>0)

fy = fy;

else

fy = pi - fy;

end

else

if(tanfy<0)

fy = -fy;

else

fy = fy - pi;

end

end

pk_2_ = e_2 * ( 1 - sinA0*sinA0 );

pk_4_ = pk_2_*pk_2_;

e_4 = e_2*e_2;

e_6 = e_2*e_4;

pa = (e_2/2+e_4/8+e_6/16)-e_2*(1+e_2)*pk_2_/16+3*e_2*pk_4_/128;

pb = e_2*(1+e_2)*pk_2_/16 - e_2*pk_4_/32;

pc = e_2*pk_4_/256;

l = fy - sinA0*(pa*r+pb*sin(r)*cos(2*r1+r)+pc*sin(2*r)*cos(4*r1+2*r)); L2 = L1+l;

if(sinA1>0)

if(tanA2>0)

A2 = pi+A2;

else

A2 = 2*pi - A2;

end

else

if(tanA2>0)

A2 = A2;

else

A2 = pi - A2;

end

end

%-----------------------------------------------------BesselGeodeticI nverse.m-----------------------------------

%白塞尔大地主题反算

%B2, L2, A1, B1, L1, A2 单位 rad

%S 单位 m

function [S, A1, A2] = BesselGeodeticInverse(B1, L1, B2, L2, eps, RefEllipsoid)

a = RefEllipsoid(1);

b = RefEllipsoid(2);

c = RefEllipsoid(3);

alpha = RefEllipsoid(4);

e_2 = RefEllipsoid(5);

e_2_ = RefEllipsoid(6); %e'_2

e = e_2^0.5;

%step1 auxiliary computation

sinB1 = sin(B1);

sinB1_2 = sinB1*sinB1;

sinB2 = sin(B2);

sinB2_2 = sinB2*sinB2;

cosB1 = cos(B1);

cosB2 = cos(B2);

W1 = (1-e_2*sinB1_2)^0.5;

W2 = (1-e_2*sinB2_2)^0.5;

sinu1 = sinB1*(1-e_2)^0.5/W1;

sinu2 = sinB2*(1-e_2)^0.5/W2;

cosu1 = cosB1/W1;

cosu2 = cosB2/W2;

l = L2 - L1;

a1 = sinu1*sinu2;

a2 = cosu1*cosu2;