用matlab写一个GPS时间与坐标系统转换程序

好的，下面是用MATLAB编写的GPS时间与坐标系统转换程序：

function [x,y,z] = gps2xyz(t, lat, lon, h)

% t: GPS时间，单位为秒
% lat: 纬度，单位为度
% lon: 经度，单位为度
% h: 海拔高度，单位为米

% 常数定义
a = 6378137; % WGS84椭球体长半轴
f = 1/298.257223563; % WGS84椭球体扁率
e2 = f*(2-f); % WGS84椭球体第一偏心率的平方

% 计算地球半径
N = a / sqrt(1 - e2 * sin(lat*pi/180)^2);
R = [N+h, (N+h)*cos(lat*pi/180)*sin(lon*pi/180), (N+h)*cos(lat*pi/180)*cos(lon*pi/180)];

% 计算UTC时间
t = t - 18; % 时区为北京时间东八区，需要减去18秒
y = str2double(datestr(datenum([1980, 1, 6, 0, 0, t]), 'yyyy'));
m = str2double(datestr(datenum([1980, 1, 6, 0, 0, t]), 'mm'));
d = str2double(datestr(datenum([1980, 1, 6, 0, 0, t]), 'dd'));
h = str2double(datestr(datenum([1980, 1, 6, 0, 0, t]), 'HH'));
min = str2double(datestr(datenum([1980, 1, 6, 0, 0, t]), 'MM'));
s = str2double(datestr(datenum([1980, 1, 6, 0, 0, t]), 'SS'));

% 计算儒略日
JD = 367*y - fix(7*(y+fix((m+9)/12))/4) + fix(275*m/9) + d + 1721013.5 + ((s/60+min)/60+h)/24;

% 计算GPS周数
GPS_Week = fix((JD-2444244.5)/7);

% 计算GPS周内秒数
GPS_Sec = (JD-2444244.5-GPS_Week*7)*24*3600;

% 计算GPS周内秒数对应的UTC时间
UTC_Year = y;
UTC_Month = m;
UTC_Day = d;
UTC_Hour = h;
UTC_Minute = min;
UTC_Second = s-GPS_Sec;

if UTC_Second < 0
    UTC_Second = UTC_Second + 60;
    UTC_Minute = UTC_Minute - 1;
end

if UTC_Minute < 0
    UTC_Minute = UTC_Minute + 60;
    UTC_Hour = UTC_Hour - 1;
end

if UTC_Hour < 0
    UTC_Hour = UTC_Hour + 24;
    UTC_Day = UTC_Day - 1;
end

if UTC_Day < 1
    if UTC_Month == 3
        if mod(UTC_Year, 4) == 0
            if mod(UTC_Year, 100) == 0
                if mod(UTC_Year, 400) == 0
                    days = 29;
                else
                    days = 28;
                end
            else
                days = 29;
            end
        else
            days = 28;
        end
    else
        if UTC_Month == 1 || UTC_Month == 2 || UTC_Month == 4 || UTC_Month == 6 || UTC_Month == 8 || UTC_Month == 9 || UTC_Month == 11
            days = 31;
        else
            days = 30;
        end
    end
    
    UTC_Day = UTC_Day + days;
    UTC_Month = UTC_Month - 1;
end

if UTC_Month < 1
    UTC_Month = UTC_Month + 12;
    UTC_Year = UTC_Year - 1;
end

% 计算UTC时间对应的GPS周数和秒数
GPS_Week = fix((JD-2444244.5)/7);
GPS_Sec = (JD-2444244.5-GPS_Week*7)*24*3600;

% 计算GPS周内秒数对应的GPS时间
TOW = mod(GPS_Sec, 604800);

% 计算GPS周内秒数对应的GPS日期
GPS_Year = UTC_Year;
GPS_Month = UTC_Month;
GPS_Day = UTC_Day - 1;

if TOW < 0
    TOW = TOW + 604800;
    GPS_Day = GPS_Day - 1;
end

if GPS_Day < 1
    if GPS_Month == 3
        if mod(GPS_Year, 4) == 0
            if mod(GPS_Year, 100) == 0
                if mod(GPS_Year, 400) == 0
                    days = 29;
                else
                    days = 28;
                end
            else
                days = 29;
            end
        else
            days = 28;
        end
    else
        if GPS_Month == 1 || GPS_Month == 2 || GPS_Month == 4 || GPS_Month == 6 || GPS_Month == 8 || GPS_Month == 9 || GPS_Month == 11
            days = 31;
        else
            days = 30;
        end
    end
    
    GPS_Day = GPS_Day + days;
    GPS_Month = GPS_Month - 1;
end

if GPS_Month < 1
    GPS_Month = GPS_Month + 12;
    GPS_Year = GPS_Year - 1;
end

% 计算GPS周内秒数对应的GPS时间
GPS_Time = GPS_Week*604800 + TOW;

% 计算WGS84坐标
WGS84 = [0, 0, 0];
WGS84(1) = R(1)*cos(GPS_Time*7.2921151467e-5);
WGS84(2) = R(1)*sin(GPS_Time*7.2921151467e-5);
WGS84(3) = R(2);

% 计算ENU坐标
phi = lat;
lambda = lon;
R_N = a / sqrt(1 - e2 * sin(phi*pi/180)^2);
R_M = a * (1 - e2) / ((1 - e2 * sin(phi*pi/180)^2)^1.5);
N = [ -sin(phi*pi/180) * cos(lambda*pi/180); -sin(phi*pi/180) * sin(lambda*pi/180); cos(phi*pi/180)];
E = [ -sin(lambda*pi/180); cos(lambda*pi/180); 0];
U = [ cos(phi*pi/180) * cos(lambda*pi/180); cos(phi*pi/180) * sin(lambda*pi/180); sin(phi*pi/180)];
ENU = [E'; N'; U'] * (WGS84' - R');

% 输出结果
x = ENU(1);
y = ENU(2);
z = ENU(3);
end

使用方法：

[x,y,z] = gps2xyz(t, lat, lon, h);

其中，t为GPS时间，单位为秒；lat为纬度，单位为度；lon为经度，单位为度；h为海拔高度，单位为米。函数返回值x、y、z为ENU坐标系下的坐标，单位为米。