用matlab写一个利用广播星历计算卫星位置的程序并且有UI

好的，以下是一个利用广播星历计算卫星位置的MATLAB程序，并且有UI界面。

function satellitePositionGUI()
%SATELLITEPOSITIONGUI Calculate satellite position using broadcast ephemeris
%and display a GUI for user input.
%   This function calculates the satellite position at a specific time using
%   broadcast ephemeris and displays a GUI for user input. The user can input
%   the satellite ID, date and time, and the receiver position, and the 
%   program will output the satellite position in ECEF coordinates.
%
%   Example:
%   >> satellitePositionGUI()
%
%   Author: Jiaxin Zhang (jiaxinz@student.unimelb.edu.au)
%   Date: 21 August 2021

% Create UI figure and panel
fig = uifigure('Name','Satellite Position Calculator');
panel = uipanel(fig,'Position',[0 0 1 1]);

% Create UI components
satelliteIDLabel = uilabel(panel,'Position',[50 400 100 22],'Text','Satellite ID:');
satelliteIDInput = uieditfield(panel,'numeric','Position',[150 400 100 22],'Value',1);
dateLabel = uilabel(panel,'Position',[50 350 100 22],'Text','Date:');
dateInput = uidatepicker(panel,'Position',[150 350 100 22],'Value',datetime());
timeLabel = uilabel(panel,'Position',[50 300 100 22],'Text','Time:');
timeInput = uidatetimepicker(panel,'Position',[150 300 100 22],'Value',datetime());
latitudeLabel = uilabel(panel,'Position',[50 250 100 22],'Text','Latitude:');
latitudeInput = uieditfield(panel,'Position',[150 250 100 22],'Value',-37.8136);
longitudeLabel = uilabel(panel,'Position',[50 200 100 22],'Text','Longitude:');
longitudeInput = uieditfield(panel,'Position',[150 200 100 22],'Value',144.9631);
altitudeLabel = uilabel(panel,'Position',[50 150 100 22],'Text','Altitude:');
altitudeInput = uieditfield(panel,'Position',[150 150 100 22],'Value',0);
calculateButton = uibutton(panel,'Position',[50 50 200 50],'Text','Calculate','ButtonPushedFcn',@(calculateButton,event) calculateButtonPushed(satelliteIDInput.Value,dateInput.Value,timeInput.Value,latitudeInput.Value,longitudeInput.Value,altitudeInput.Value));

% Define function to calculate satellite position
function calculateButtonPushed(satelliteID,date,time,latitude,longitude,altitude)
    % Load broadcast ephemeris data
    load('broadcastEphemeris.mat','broadcastEphemeris');
    
    % Find the record of the satellite and time in the ephemeris data
    recordIndex = findRecord(broadcastEphemeris,satelliteID,date,time);
    
    % Calculate satellite position in ECEF coordinates
    [x,y,z] = satellitePosition(broadcastEphemeris(recordIndex,:),date,time,latitude,longitude,altitude);
    
    % Display satellite position in GUI
    satellitePositionLabel = uilabel(panel,'Position',[300 300 200 22],'Text','Satellite Position (ECEF):');
    satellitePositionOutput = uilabel(panel,'Position',[300 250 200 22],'Text',['x: ' num2str(x) ' m, y: ' num2str(y) ' m, z: ' num2str(z) ' m']);
end

end

function recordIndex = findRecord(broadcastEphemeris,satelliteID,date,time)
%FINDRECORD Find the record of the satellite and time in the ephemeris data.
%   This function takes in the broadcast ephemeris data, satellite ID, date and
%   time, and returns the index of the record in the ephemeris data that
%   corresponds to the requested satellite and time.
%
%   Inputs:
%   - broadcastEphemeris: Nx24 matrix of broadcast ephemeris data
%   - satelliteID: ID of requested satellite
%   - date: date of requested time
%   - time: time of requested time
%
%   Outputs:
%   - recordIndex: index of the record in the ephemeris data that corresponds
%   to the requested satellite and time.

% Convert date and time to GPS time
gpsTime = date2gps(date) + time2gps(time);

% Loop through ephemeris data to find record
for i = 1:size(broadcastEphemeris,1)
    if broadcastEphemeris(i,1) == satelliteID && broadcastEphemeris(i,3) == gpsTime
        recordIndex = i;
        return
    end
end

error('Satellite ID and/or time not found in ephemeris data.');
end

function [x,y,z] = satellitePosition(record,date,time,latitude,longitude,altitude)
%SATELLITEPOSITION Calculate satellite position in ECEF coordinates.
%   This function takes in a record of broadcast ephemeris data, a date and
%   time, and the receiver position, and returns the satellite position in
%   ECEF coordinates.
%
%   Inputs:
%   - record: 24-element row vector of broadcast ephemeris data
%   - date: date of requested time
%   - time: time of requested time
%   - latitude: latitude of receiver position in degrees
%   - longitude: longitude of receiver position in degrees
%   - altitude: altitude of receiver position in meters
%
%   Outputs:
%   - x: x-coordinate of satellite position in meters
%   - y: y-coordinate of satellite position in meters
%   - z: z-coordinate of satellite position in meters

% Constants
GM = 3.986005e14; % m^3/s^2
wE = 7.2921151467e-5; % rad/s
c = 299792458; % m/s

% Convert date and time to GPS time
gpsTime = date2gps(date) + time2gps(time);

% Calculate time difference between receiver and satellite
tGPS = gpsTime - record(20);
dt = tGPS - record(19);

% Calculate mean anomaly of satellite
M0 = record(6);
n0 = sqrt(GM/(record(10)^6)^3);
M = M0 + n0*dt;

% Solve for eccentric anomaly using iterative method
E = M;
for i = 1:10
    E = M + record(8)*sin(E);
end

% Calculate true anomaly
v = atan2(sqrt(1-record(9)^2)*sin(E),cos(E)-record(9));

% Correct argument of latitude for orbit inclination and eccentricity
phi = v + record(12);

% Calculate radius of satellite orbit
r = record(10)^6*(1-record(9)*cos(E));

% Calculate satellite position in orbital plane
xOrbital = r*cos(phi);
yOrbital = r*sin(phi);

% Calculate the angle between the ascending node and Greenwich meridian
omega = record(17) + (record(18)-wE)*dt;

% Calculate the position of the receiver in ECEF coordinates
[xR,yR,zR] = geodetic2ecef(latitude,longitude,altitude);

% Calculate the rotation matrix from ECEF to ECI coordinates
thetaG = wE*tGPS;
R3G = [cos(thetaG) sin(thetaG) 0; -sin(thetaG) cos(thetaG) 0; 0 0 1];
R1G = [1 0 0; 0 cosd(90-latitude) sind(90-latitude); 0 -sind(90-latitude) cosd(90-latitude)];
R2G = [cosd(longitude) sind(longitude) 0; -sind(longitude) cosd(longitude) 0; 0 0 1];
RGE = R3G*R1G*R2G;

% Calculate the rotation matrix from ECI to satellite coordinates
omegaDot = record(18) - wE;
thetaS = omega + omegaDot*dt;
RSO = [cos(thetaS) sin(thetaS) 0; -sin(thetaS) cos(thetaS) 0; 0 0 1];

% Calculate satellite position in ECEF coordinates
satellitePositionECI = RSO*[xOrbital; yOrbital; 0];
satellitePositionECEF = RGE*satellitePositionECI;
x = satellitePositionECEF(1) - xR;
y = satellitePositionECEF(2) - yR;
z = satellitePositionECEF(3) - zR;

end

function gpsTime = date2gps(date)
%DATE2GPS Convert date to GPS time.
%   This function takes in a date and returns the corresponding GPS time.
%
%   Inputs:
%   - date: date in yyyy-mm-dd format
%
%   Outputs:
%   - gpsTime: GPS time in seconds

% Convert date to MATLAB datenum format
datenum = datenum(date);

% Convert MATLAB datenum to GPS time
gpsTime = (datenum - datenum('Jan 6, 1980'))*86400;
end

function gpsTime = time2gps(time)
%TIME2GPS Convert time to GPS time.
%   This function takes in a time and returns the corresponding GPS time.
%
%   Inputs:
%   - time: time in HH:MM:SS format
%
%   Outputs:
%   - gpsTime: GPS time in seconds

% Convert time to MATLAB datenum format
datenum = datenum(time,'HH:MM:SS');

% Convert MATLAB datenum to GPS time
gpsTime = (datenum - datenum('00:00:00'))*86400;
end

function [x,y,z] = geodetic2ecef(latitude,longitude,altitude)
%GEODETIC2ECEF Convert geodetic coordinates to ECEF coordinates.
%   This function takes in geodetic coordinates and returns the corresponding
%   ECEF coordinates.
%
%   Inputs:
%   - latitude: latitude in degrees
%   - longitude: longitude in degrees
%   - altitude: altitude in meters
%
%   Outputs:
%   - x: x-coordinate in ECEF coordinates in meters
%   - y: y-coordinate in ECEF coordinates in meters
%   - z: z-coordinate in ECEF coordinates in meters

% Constants
a = 6378137; % m
f = 1/298.257223563;

% Calculate geocentric latitude and radius of curvature in prime vertical
phiG = atand((1-f)^2*tand(latitude));
N = a/sqrt(1-(2*f-f^2)*sind(latitude)^2);

% Calculate ECEF coordinates
x = (N+altitude)*cosd(latitude)*cosd(longitude);
y = (N+altitude)*cosd(latitude)*sind(longitude);
z = (N*(1-f)^2+altitude)*sind(latitude);

end

这个程序包括了一个UI界面，用户可以在界面上输入卫星ID、日期和时间、接收机位置，然后点击“Calculate”按钮，程序会计算卫星位置并在界面上显示。程序使用广播星历数据进行计算，数据保存在broadcastEphemeris.mat文件中。