function [Xpuwg, Ypuwg]=WGS84_do_PUWG92(B_stopnie, L_stopnie);
% Konwersja współrzędnych układu WGS84 do układu PUWG92
% Opis:
%    konwersja wspolrzednych z ukladu WGS 84 do ukladu PUWG 1992
% Parametry:
%    B_stopnie - szerokosc geograficzna wyrazona w stopniach
%    L_stopnie - dlugosc geograficzna wyrazona w stopniach
% Wartosc zwracana
%    Xpuwg - wspolrzedna X ukladu PUWG 1992 (UWAGA - wspolrzedna pionowa) 
%            lub NaN, jesli parametry wejsciowe sa poza zakresem
%    Ypuwg - wspolrzedna Y ukladu PUWG 1992 (UWAGA - wspolrzedna pozioma)
%            lub NaN, jesli parametry wejsciowe sa poza zakresem

%
% Autor: Zbigniew Szymanski
% E-mail: z.szymanski@szymanski-net.eu
% Wersja: 1.0
% Data modyfikacji: 2012-11-17
% Uwagi: Oprogramowanie darmowe. Dozwolone jest wykorzystanie i modyfikacja 
%       niniejszego oprogramowania do wlasnych celow pod warunkiem 
%       pozostawienia wszystkich informacji z naglowka. W przypadku 
%       wykorzystania niniejszego oprogramowania we wszelkich projektach
%       naukowo-badawczych, rozwojowych, wdrozeniowych i dydaktycznych prosze
%       o zacytowanie nastepujacego artykulu:
       
%       Zbigniew Szymanski, Stanislaw Jankowski, Jan Szczyrek, 
%       "Reconstruction of environment model by using radar vector field histograms.",
%       Photonics Applications in Astronomy, Communications, Industry, and 
%       High-Energy Physics Experiments 2012, Proc. of SPIE Vol. 8454, pp. 845422 - 1-8,
%       doi:10.1117/12.2001354
%       
%Literatura:
%       Uriasz, J., “Wybrane odwzorowania kartograficzne”, Akademia Morska w Szczecinie,
%       http://uriasz.am.szczecin.pl/naw_bezp/odwzorowania.html


Xpuwg=NaN;
Ypuwg=NaN;

%% Parametry elipsoidy GRS-80
e=0.0818191910428;  % pierwszymimośród elipsoidy
R0=6367449.14577; %promień sfery Lagrange’a
Snorm=2.0E-6;   %parametr normujący
xo=5760000; %parametr centrujący
a0=5765181.11148097;
a1=499800.81713800;
a2=-63.81145283;
a3=0.83537915;
a4=0.13046891;
a5=-0.00111138;
a6=-0.00010504;

%% Parametry odwzorowania Gaussa-Kruegera dla układu PUWG92
L0_stopnie=19; %Początek układu wsp. PUWG92 (długość)
m0=0.9993;
x0=-5300000;
y0= 500000;


%% Zakres stosowalnosci metody
Bmin=48*pi/180;
Bmax=56*pi/180;
dLmin=-6*pi/180;
dLmax=6*pi/180;

%% Weryfikacja danych wejsciowych
B=B_stopnie*pi/180;
dL_stopnie=L_stopnie-L0_stopnie;
dL=dL_stopnie*pi/180;

if ((B<Bmin) || (B>Bmax))
   disp('Blad: niepoprawny zakres parametru B');
   return;
end

if ((dL<dLmin) || (dL>dLmax))
   disp('Blad: niepoprawny zakres parametru L');
   return;
end

%% etap I - elipsoida na kulę
% Lagrange’a

U=1-e*sin(B);
V=1+e*sin(B);
K=(U/V)^(e/2);
C=K*tan(B/2+pi/4);
fi=2*atan(C)-pi/2;
d_lambda=dL;

%% etap II - kula na walec
% Merkatora

p=sin(fi);
q=cos(fi)*cos(d_lambda);
r=1+cos(fi)*sin(d_lambda);
s=1-cos(fi)*sin(d_lambda);
XMERC=R0*atan(p/q);
YMERC=0.5*R0*log(r/s);

%% etap III - walec na płaszczyznę
% (XMERC,YMERC)->(XGK,YGK)
Z=complex((XMERC-xo)*Snorm, YMERC*Snorm);
Zgk=a0+Z*(a1+Z*(a2+Z*(a3+Z*(a4+Z*(a5+Z*a6)))));
Xgk=real(Zgk);
Ygk=imag(Zgk);

%% Przejście do układu aplikacyjnego
Xpuwg=m0*Xgk+x0;
Ypuwg=m0*Ygk+y0;