function [l,b,mul,mub]=equ2gal(ra,de,mura,mude);
%%
%
%  Calculating galactic coordinates and
%  proper motion components:
%  Transformation: Equatorial ----> Galactic
%  Epoch and equinix J2000
%
%  Function call: [l,b,mul,mub]=equ2gal(ra,de,mura,mude)
%
%
%  Input:   ra, de - eequatorial coordinates (in degrees)
%           mura, mude - proper motion components
%           in the equatorial coordinates (arcsec/year, msec/year ,,,)
%  Output:  l, b - galactic coordinates (in degrees)
%           mul, mub - proper motion components
%           in the galactic coordinates, in the same units as mura, mude
%                                         (arcsec/year, msec/year ,,,)
%
% Programmer: A.S.Rastorguev, Lomonosov Moscow State University,
%                             Sternberg Astronomical Institute, 2013
% -------------------------------------------------------------------------

%% Calculating galactic coordinates:
% From: HIPPARCOS: Section 1.5
% "Transformation of Astrometric Data and Associated Error Propagation"
% Expression (1.5.11):
Ag = [ -0.0548755604  -0.8734370902   -0.4838350155,
       +0.4941094279  -0.4448296300   +0.7469822445,
       -0.8676661990  -0.1980763734   +0.4559837762 ];

rad=180/pi;
ra=ra/rad; de=de/rad;
% Unit radius-vector (column vector) in the equatorial coordinate frame:
r=[ cos(de)*cos(ra),
    cos(de)*sin(ra),
            sin(de) ];
% Unit radius-vector (column vector) in the galactic coordinate frame:
rg=Ag*r;
% Calculating galactic longitude and latitude:
b=asin(rg(3,1));
l=atan2(rg(2,1),rg(1,1));
if l < 0,
      l=l+2*pi;
end
l=l*rad;
b=b*rad;

%% Calculating proper motions in the galactic coordinates:
% From: HIPPARCOS: Section 1.5
% "Transformation of Astrometric Data and Associated Error Propagation"
% Expression (1.5.15)
Z=[0 0 1]';
Zg=Ag(3,:)';
% Zg=[-0.8676661990 -0.1980763734 +0.4559837762]';

P=cross(Z,r)/norm(cross(Z,r));
Pg=cross(Zg,r)/norm(cross(Zg,r));

Q=cross(r,P);
Qg=cross(r,Pg);

% Формула (1.5.19)
M=[ Pg'*P  Pg'*Q,
    Qg'*P  Qg'*Q ];

MUe=[ mura,
      mude ];

MUg=M*MUe;

mul=MUg(1);
mub=MUg(2);
% End