function mut = getmut(mjd2000,theta_NGP,phi_NGP)

% mut = getmut(mjd2000,theta_NGP,phi_NGP)
%
% Calculate the Magnetic Universal Time for a given MJD2000 and a given 
% position of the North Geomagnetic Pole.
%
% Input:    mjd2000         modified julian day 2000 (days)
%           theta_NGP       geographic co-latitude of NGP (deg)
%           phi_NGP         geographic longitude of NGP (deg)
% Output:   mut             magnetic universal time (hr)
%
% A. Chulliat, October 2008

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rad = pi/180;

theta_0 = theta_NGP;
phi_0   = phi_NGP;


[rasc, decl] = sun_md2000(mjd2000);
theta_s = 90 - decl/rad;

% old: 
% ut = mod(mjd2000, 1)*24;
% phi_s = 15*(12-ut);

% NIO: 
gst = gast_md2000(mjd2000);
phi_s = (rasc - gst)/rad;

eopp = sind(theta_s).*sind(phi_s-phi_0);
eadj = cosd(theta_0).*sind(theta_s).*cosd(phi_s-phi_0) ...
    - sind(theta_0).* cosd(theta_s);

mut = 12 - atan2(eopp,eadj)/rad/15;





            
      
