BESSEL_2.M
上传用户:sfyaiting
上传日期:2009-10-25
资源大小:320k
文件大小:4k
源码类别:

GPS编程

开发平台:

Matlab

  1. function [s12,A1,A2] = bessel_2(phi1d,phi1m,phi1s,l1d,l1m,...
  2.                         l1s,phi2d,phi2m,phi2s,l2d,l2m,l2s,a,finv)
  3. %BESSEL_2 Solution of the inverse geodetic problem according to
  4. %         the Bessel-Helmert method as described in Zhan Xue-Lian.
  5. %         Given two points with coordinates (phi1, l1) and 
  6. %         (phi2,l2). Per definition we always have l2 > l1.
  7. %         The given reference ellipsoid has semi-major
  8. %         axis a and inverse flattening finv. The coordinates
  9. %         are in the format degree, minute, and second with decimals. 
  10. %
  11. %         Zhan Xue-Lian (1985) The nested coefficient method for
  12. %            accurate solutions of direct and inverse geodetic 
  13. %            problems with any length. Proceedings of the 7th
  14. %            International Symposium on Geodetic Computations.
  15. %            Cracow, June 18--21, 1985. Institute of Geodesy and
  16. %            Cartography. Wasaw, Poland, ul. Jasna 2/4
  17. %
  18. %         For a good background of the problems, see 
  19. %
  20. %         Bodem"uller, H.(1954): Die geod"atischen Linien des 
  21. %          Rotationsellipsoides und die L"osung der geod"atischen
  22. %          Hauptaufgaben f"ur gross{}e Strecken unter 
  23. %            besonderer Ber"ucksichtigung der Bessel-Helmertschen
  24. %            L"osungsmethode.
  25. %          Deutsche Geod"atische Kommission, Reihe B, Nr. 13.
  26. %Kai Borre, January 25, 1999
  27. %Copyright (c) by Kai Borre
  28. %$Revision 1.0 $  $Date:1999/01/25  $
  29. if nargin == 0
  30.    phi1  = dms2rad(50,0,0);
  31.    l1 = dms2rad(10,0,0);
  32.    phi2 = dms2rad(-62,57,3.203824);
  33.    l2 = dms2rad(105,5,38.299430);
  34.    a = 6378388;
  35.    finv = 297;
  36. else
  37.    phi1 = dms2rad(phi1d,phi1m,phi1s);
  38.    l1 = dms2rad(l1d,l1m,l1s);
  39.    phi2 = dms2rad(phi2d,phi2m,phi2s);
  40.    l2 = dms2rad(l2d,l2m,l2s);
  41. end
  42. dtr = pi/180;
  43. f = 1/finv; 
  44. ex2 = (2-f)*f/(1-f)^2;
  45. %Reduced latitudes
  46. u1 = atan((1-f)*tan(phi1)); 
  47. u2 = atan((1-f)*tan(phi2));
  48. deltaomega_old = 0;
  49. deltaomega_new = 1;
  50. while abs(deltaomega_old-deltaomega_new) > 1.e-12
  51.    deltaomega_old = deltaomega_new;
  52.    omega = l2-l1+deltaomega_old;
  53.    sigma = atan2(sqrt((cos(u2)*sin(omega))^2+...
  54.       (cos(u1)*sin(u2)-sin(u1)*cos(u2)*cos(omega))^2),...
  55.       sin(u1)*sin(u2)+cos(u1)*cos(u2)*cos(omega)); 
  56.    cosun = cos(u1)*cos(u2)*sin(omega)/sin(sigma);   
  57.    sinun2 = 1-cosun^2;   
  58.    if sinun2 == 0
  59.       sigmam = acos(cos(sigma)-2);
  60.    else
  61.       sigmam = acos(cos(sigma)-2*sin(u1)*sin(u2)/sinun2);   
  62.    end
  63.    v = f*sinun2/4;
  64.    K3 = v*(1+f+f^2-v*(3+7*f-13*v));
  65.    deltaomega_new = (1-K3)*f*cosun*(sigma+K3*sin(sigma)*(cos(sigmam)+...
  66.       K3*cos(sigma)*cos(2*sigmam)));
  67. end
  68. t = ex2*sinun2/4;
  69. K1 = 1+t*(1-t*(3-t*(5-11*t))/4);
  70. K2 = t*(1-t*(2-t*(37-94*t)/8));
  71. deltasigma = K2*sin(sigma)*(cos(sigmam)+...
  72.    K2*(cos(sigma)*cos(2*sigmam)+...
  73.    K2*(1+2*cos(2*sigma))*cos(3*sigmam)/6)/4);
  74. s12 = K1*(1-f)*a*(sigma-deltasigma);
  75. A1 = atan2(cos(u2)*sin(omega),...
  76.    cos(u1)*sin(u2)-sin(u1)*cos(u2)*cos(omega));
  77. A2 = atan2(cos(u1)*sin(omega),...
  78.    cos(u1)*sin(u2)*cos(omega)-sin(u1)*cos(u2));
  79.        
  80. %----------------------------------------------
  81. function result = dms2rad(deg,min,sec);
  82. % Conversion of degrees, minutes, and seconds to radians
  83. neg_arg = 'FALSE';
  84. if deg < 0
  85.    neg_arg = 'TRUE ';
  86.    deg = -deg;
  87. end
  88. arg = deg+min/60+sec/3600;
  89. result = arg*pi/180;
  90. if neg_arg == 'TRUE ';
  91.    result = -result;
  92. end
  93. %------------------------------------------
  94. function result = rad2dms(arg)
  95. %RAD2DMS Conversion of radians to degrees, minutes, and seconds
  96. neg_arg = 'FALSE';
  97. if arg < 0
  98.    neg_arg = 'TRUE ';
  99.    arg = -arg;
  100. end
  101. arg = arg*180/pi;
  102. result = zeros(1,3);
  103. result(1) = fix(arg);
  104. if result(1) == 0
  105.    result(2) = fix(arg*60);
  106. else
  107.    result(2) = fix(rem(arg,result(1))*60);
  108. end
  109. result(3) = (arg-result(1)-result(2)/60)*3600;
  110. if neg_arg == 'TRUE '
  111.    result(1) = -result(1);
  112. end
  113. fprintf('   %3.0f %2.0f %8.6fn',result(1),result(2),result(3))
  114. %%%%%%%%%%%%%%%%% end bessel_2.m  %%%%%%%%%%%%%%%%%%%%%%%%%%%