bl2gk.m

这是国外关于卫星导航方面一书的源代码

function [R,H]=bl2gk(bin,lin,l0,ell)
% [R,H]=bl2gk(B,L,L0)
% [R,H]=bl2gk(B,L,L0,ELLIPS)
%
% Umwandlung von Geodaetischen ellipsoidischen Koordinaten
%            in  Gauss-Krueger-Koordinaten (deutscher Konvention)
% Input:  Breite [grad.decimale], Laenge [grad.decimale], Nullmeridian [grad],
%         ELLIPS ist Textstring: 'wgs84', 'grs80', oder 'bessel' (default)
% Output: Rechtswert [m], Hochwert [m]
%
% Geschlossene Formeln, Reihenentwicklung nur fuer Meridianbogenlaenge aus
% Fusspunktsbreite nach K. Schnaedelbach, Skript Landesvermessung 
%
% (c) iapg - Zeb - 11/96

if nargin < 4, ell = 'bessel'; end	% Bessel default
ell = lower(deblank(ell));
if strcmp(ell,'bessel')			% Bessel-Ellipsoid 
   a  = 6377397.155;         		   % grosse Halbachse
   c  = 6398786.849;         		   % a*a/b
   es = 0.006719218798;     		   % (a*a-b*b)/(b*b)
   n  = 0.001674184801;      		   % (a-b)/(a+b)
elseif strcmp(ell,'wgs84') 		% WGS84-Ellipsoid
   a  = 6378137;
   c  = 6399593.626;
   es = 0.00673949674226;
   n  = 0.00167922038638;
elseif strcmp(ell,'grs80') 		% GRS80-Ellipsoid
   a  = 6378137;
   c  = 6399593.626;
   es = 0.00673949677548;
   n  = 0.00167922039463;
else
   error('Nicht-definiertes Ellipsoid')
end

   rho=180 / pi;

% GK in der Bundesrepublik Deutschland

   masstab = 1;
   y0 = 500000;


   b = bin / rho;
   l = (lin-l0) / rho;

% Berechnung Fusspunktsbreite

   eta2 = es .* (cos(b)).^2;
   v = sqrt(1+eta2);
   bf = atan(tan(b) ./ cos(v.*l) .* (1 + eta2/6 .* (1-3*sin(b).^2) .* l.^4));

% Berechnung der Meridianbogenlaenge aus Fusspunktsbreite

   b0 = (1 + n^2 / 4 + n^4 / 64) * bf;
   b2 = -1.5*n * (1-n^2/8) * sin(2.*bf);
   b4 = 15/16 * n^2 * (1 - n^2/4)*sin(4*bf);
   b6 = -35 / 48 * n^3 * sin(6*bf);
   b8 = 315 / 512 * n^4 * sin(8*bf);

   g = a / (1+n) * (b0 + b2 + b4 + b6 + b8);

% Berechnung Hochwert

   H = g * masstab;

% Berechnung Rechtswert

   eta2f = es .* cos(bf).^2;
   vf = sqrt(1+eta2f);
   ys0 = tan(l) .* cos(bf) ./ vf;
   ys0 = ys0 .* (1 + eta2f .* l.^2 .* cos(bf).^2 .* (eta2f/6 + l.^2/10) );
   ys = asinh(ys0);
   R = (ys .* c * masstab) + (l0/3)*1e6 + y0;

% -------------------Ende----------------------------------Zeb,11/96--------