写一个完整可运行的Matlab代码，计算某经纬度一天的太阳高度角和方位角

以下是一个计算某经纬度一天的太阳高度角和方位角的Matlab代码：

clc;
clear all;
close all;

% 地球半径
R = 6371;

% 经纬度
lat = 31.23; % 纬度
lon = 121.47; % 经度

% 时间
date = datetime('now','TimeZone','local','Format','yyyy-MM-dd');
yr = year(date);
mon = month(date);
day = day(date);

% 一天的时间范围
t = linspace(0,24,24*60);

% 计算太阳位置
for i = 1:length(t)
    [az, el] = sun_position(yr, mon, day, t(i), lat, lon, R);
    azimuth(i) = az;
    elevation(i) = el;
end

% 计算太阳高度角和方位角
h = 90 - elevation; % 太阳高度角
A = azimuth; % 太阳方位角

% 绘制太阳高度角和方位角图像
figure
subplot(2,1,1)
plot(t,h)
title('太阳高度角')
xlabel('时间（小时）')
ylabel('高度角（度）')
ylim([0,90])

subplot(2,1,2)
plot(t,A)
title('太阳方位角')
xlabel('时间（小时）')
ylabel('方位角（度）')
ylim([0,360])

function [az, el] = sun_position(yr, mon, day, hr, lat, lon, R)
    % 计算儒略日
    J2000 = 2451545;
    J = J2000 + (datenum(yr,mon,day)-datenum(2000,1,1)) + (hr-12)/24;

    % 计算黄道坐标
    T = (J - 2451545)/36525;
    L = 280.460 + 36000.772*T;
    M = 357.528 + 35999.050*T;
    lambda = L + 1.915*sin(deg2rad(M)) + 0.020*sin(2*deg2rad(M));
    epsilon = 23.439 - 0.013*T;
    alpha = atan2(cos(deg2rad(epsilon))*sin(deg2rad(lambda)), cos(deg2rad(lambda)));
    delta = asin(sin(deg2rad(epsilon))*sin(deg2rad(lambda)));

    % 计算赤道坐标
    GMST = mod(280.46061837 + 360.98564736629*(J-2451545), 360);
    H = mod(GMST + lon, 360);
    alpha = rad2deg(alpha);
    delta = rad2deg(delta);
    H = H*15;
    [x, y, z] = sph2cart(deg2rad(H), deg2rad(delta), R);
    [x1, y1, z1] = sph2cart(deg2rad(lat), deg2rad(90), R);
    [x2, y2, z2] = sph2cart(deg2rad(alpha), deg2rad(90-delta), R);
    x3 = x2-x1;
    y3 = y2-y1;
    z3 = z2-z1;
    [az, el, r] = cart2sph(x3, y3, z3);
    az = rad2deg(az);
    el = rad2deg(el);
end

这个代码使用了一个名为 `sun_position` 的函数，用于计算太阳位置的赤道坐标，然后再计算出太阳在某个地点的高度角和方位角。这个函数的输入参数包括年、月、日、小时、纬度、经度和地球半径。输出参数是太阳的方位角和高度角。主函数则使用 `linspace` 函数生成一天的时间范围，然后遍历这个范围，计算每个时刻的太阳高度角和方位角，并绘制出图像。