太阳ECEF坐标计算

在电离层研究、导航精密单点定位过程中，经常需要计算太阳的地心地固ECEF坐标下的坐标，下面给出matlab计算代码，供大家参考：

function [sunECEF] = sunPositionECEF(y, m, D, UT)
% calculate approximate position of sun in Earth-Centered-Fixed-Frame [km]
% 
% INPUT:
%   y               year, 4-digit
%   m               month
%   D               day of month
%   UT              hour of day [UTC]
% OUTPUT:
%	sunECEF         3x1, sun position in ECEF [km]


deg2rad = pi/180;                   % conversion from degree to radiant
AU = 1.49597870e8;                  % [km], 1 Astronomical Unit - mean distance earth-sun
JD = cal2jd_GT(y, m, D + UT/24);    % julian date
MJD = JD - 2400000.5;               % modified julian date

% Formulas adapted from glab
fday = MJD - floor(MJD);
JDN  = MJD - 15019.5;

vl   = mod(279.696678 + 0.9856473354*JDN,360);
gstr = mod(279.690983 + 0.9856473354*JDN + 360*fday + 180,360);
g    = mod(358.475845 + 0.985600267*JDN,360)*deg2rad;

slong = vl + (1.91946-0.004789*JDN/36525)*sin(g) + 0.020094*sin(2*g);
obliq = (23.45229-0.0130125*JDN/36525)*deg2rad;

slp  = (slong-0.005686)*deg2rad;
sind = sin(obliq)*sin(slp);
cosd = sqrt(1-sind*sind);
sdec = atan2(sind,cosd)/deg2rad;
	
sra = 180 - atan2(sind/cosd/tan(obliq),-cos(slp)/cosd)/deg2rad;
	
sunECI = [  cos(sdec*deg2rad) * cos((sra)*deg2rad) * AU;
            cos(sdec*deg2rad) * sin((sra)*deg2rad) * AU;
            sin(sdec*deg2rad) * AU];
        
sunECEF = rot(3,gstr*deg2rad)*sunECI;           % transform to ECEF


function jd = cal2jd_GT(yr, mn, dy)
% cal2jd_GT  Converts calendar date to Julian date using algorithm
%   from "Practical Ephemeris Calculations" by Oliver Montenbruck
%   (Springer-Verlag, 1989). Uses astronomical year for B.C. dates
%   (2 BC = -1 yr). 
% Usage:   jd=cal2jd_GT(yr,mn,dy)
% Input:   yr - calendar year (4-digit including century)
%          mn - calendar month
%          dy - calendar day (including fractional day)
% Output:  jd - jJulian date

jd = 0;     % MFG: added initialization of output argument

if nargin ~= 3
  warning('Incorrect number of input arguments');
  return;
end
if mn < 1 | mn > 12
  warning('Invalid input month');
  return
end
if dy < 1
  if (mn == 2 & dy > 29) | (any(mn == [3 5 9 11]) & dy > 30) | (dy > 31)
    warning('Invalid input day');
    return
  end
end

if mn > 2
  y = yr;
  m = mn;
else
  y = yr - 1;
  m = mn + 12;
end

date1=4.5+31*(10+12*1582);   % Last day of Julian calendar (1582.10.04 Noon)
date2=15.5+31*(10+12*1582);  % First day of Gregorian calendar (1582.10.15 Noon)
date=dy+31*(mn+12*yr);

if date <= date1
  b = -2;
elseif date >= date2
  b = fix(y/400) - fix(y/100);
else
  warning('Dates between October 5 & 15, 1582 do not exist');
  return;
end

if y > 0
  jd = fix(365.25*y) + fix(30.6001*(m+1)) + b + 1720996.5 + dy;
else
  jd = fix(365.25*y-0.75) + fix(30.6001*(m+1)) + b + 1720996.5 + dy;
end

function R = rot(ax, alpha)
% Creates rotation matrix
% 
% INPUT:
%   ax        number of axis
%   alpha     angle to rotate
% OUTPUT:
%   rotation matrix

switch ax
    case 3      % z-axis
        R = [ cos(alpha)  sin(alpha)  0
             -sin(alpha)  cos(alpha)  0
                  0            0      1];
    case 2      % y-axis
        R = [ cos(alpha)  0  sin(alpha)
                  0       1      0
             -sin(alpha)  0  cos(alpha)];
    case 1      % x-axis
        R = [1       0           0
             0   cos(alpha)  sin(alpha)
             0  -sin(alpha)  cos(alpha)];
end

测试结果如下：

 [sunECEF] = sunPositionECEF(2024, 6, 24, 5.4)

  sunECEF =

   1.0e+08 *

   -0.2296
    1.3536
    0.5941