请问怎么把这段代码修改成循环呢，希望实现的是每个时刻的图像按时间顺序依次出图。而我现在是最粗糙的办法，每个时刻去修改，一张一张的出图
filename = 'molandi1.nc';
% 读取经纬度数据
lon = ncread(filename, 'longitude');
lat = ncread(filename, 'latitude');
u_wind = ncread(filename, 'u10');
v_wind = ncread(filename, 'v10');
t = ncread(filename,'time');
lona = double(lon);lata = double(lat);
k1 = size(u_wind,1),k2 = size(u_wind,2),k3 = size(u_wind,3)
ws = [];
for i1 = 1:k3
for i2 = 1:k1
for i3 = 1:k2
ws(i2,i3,i1) = sqrt(u_wind(i2,i3,i1).^2+v_wind(i2,i3,i1).^2);
end
end
end
[lat1,lon1] = meshgrid(lata,lona)
tf = load('CH2016BSTmolandi.txt')
% 读取CMA数据集中的台风路径数据
a = importdata('CH2016BSTmolandi.txt',' ',1);
time = a.data(:, 3); % 时间数据
lonc = a.data(:, 1); % 经度数据
latc = a.data(:, 2); % 纬度数据
Vmax = a.data(:, 5); %CMA最大风速
Jm = [];
for i4 = 2:19
m = acosd((tf(i4,1)-tf(i4+1,1))/sqrt((tf(i4+1,1)-tf(i4,1))^2+(tf(i4+1,2)-tf(i4,2))^2));
n = acosd((tf(i4-1,1)-tf(i4,1))/sqrt((tf(i4,1)-tf(i4-1,1))^2+(tf(i4,2)-tf(i4-1,2))^2));
x1 = lona;
y1 = tand(90-(m+n)/2)*(x1-tf(i4,1))+tf(i4,2);
y2(:,i4-1) = y1;
for i5 = 1:k1
Jm(i5,i4-1) = interp2(lat1,lon1,ws(:,:,79),y1(i5,1),x1(i5,1),"spline");
end
end
s1 = [];
s2 = [];
s3 = [];
%(限定了角平分线的范围)
for i6 = 1:k1
if y2(i6,13)>5&&y2(i6,13)<45
s1 = [s1;lona(i6,1)];
s2 = [s2;y2(i6,13)];
s3 = [s3;Jm(i6,13)];
end
end
time_diff = 6;
%（要改）
distance = distance_on_unit_sphere(tf(14,1), tf(14,2), tf(13,1), tf(13,2));
Vmc = distance / time_diff/3.6;
Vs = max(s3)
Vm = Vs-Vmc
Rmax1 = 51.6*exp(-0.0223*Vm+0.0281*latc(13,:));
% 圆心坐标
center_lat = latc(13,:);
center_lon = lonc(13,:);
Rmax11 = Rmax1/(111*cosd(center_lat));
% 圆半径
radius = 9.5*Rmax11;
p = polyfit(s1, s2, 1);
slope = p(1);
p_slope = -1/slope;
angle_rad = atan(p_slope);
angle = 180+rad2deg(angle_rad);
% 网格间隔
grid_interval = 0.25; % 0.25°
% 计算网格的经度和纬度范围
lon_range = center_lon - radius : grid_interval : center_lon + radius;
lat_range = center_lat - radius : grid_interval : center_lat + radius;
% 计算网格的经度和纬度
[lon12, lat12] = meshgrid(lon_range, lat_range);
% 计算每个网格点到圆心的距离
r = sqrt((lon12 - center_lon).^2 + (lat12 - center_lat).^2);
Vmov = Vmc*((r*Rmax11)./(r.^2+Rmax11^2));
vx = Vmov*cosd(angle);
vy = Vmov*sind(angle);
theta = linspace(0, 2*pi, 100);
xm = center_lon + radius * cos(theta);
ym = center_lat + radius * sin(theta);
% 计算圆内部的索引
dist = sqrt((lon12 - center_lon).^2 + (lat12 - center_lat).^2);
inside_circle = dist <= radius;
outside = dist >= radius;
U1 = u_wind(:,:,79);
V1 = v_wind(:,:,79);
U2 = vx(inside_circle);
V2 = vy(inside_circle);
% 调整矢量场的尺寸
U2 = imresize(U2, size(U1));
V2 = imresize(V2, size(V1));
% 创建网格
[YY,XX] = meshgrid(lata,lona);
% 判断网格点是否属于重叠部分，并进行相减操作
idx = (XX >= (lonc(13,:)-radius)) & (XX <= (lonc(13,:)+radius)) & (YY >= (latc(13,:)-radius)) & (YY <= (latc(13,:)+radius));
U = U1;
V = V1;
U(idx) = U1(idx) - U2(idx);
V(idx) = V1(idx) - V2(idx);
figure
plot(lonc, latc, 'r-', 'LineWidth', 1.5);
hold on;
scatter(lonc, latc, 'b');
hold on;
quiver(lon1,lat1,u_wind(:,:,79), v_wind(:,:,79));
hold on;
scatter(x,y,'g');
hold on;
plot(xm, ym,'--');
quiver(lon12(inside_circle), lat12(inside_circle), vx(inside_circle), vy(inside_circle), 'r');
hold on;
plot(s1,s2);
hold off;
axis equal;
xlabel('经度');
ylabel('纬度');
title('风场图');
fname = '16莫兰蒂风场.png';
saveas(gcf, fname, 'png');
figure
plot(lonc, latc, 'r-', 'LineWidth', 1.5);
hold on;
scatter(lonc, latc, 'b');
hold on;
quiver(XX, YY, U, V);
hold on;
scatter(x,y,'g');
hold on;
plot(xm, ym,'--');
hold on;
plot(s1,s2);
hold off;
axis equal;
xlabel('经度');
ylabel('纬度');
title('风场图');
fname = '16莫兰蒂风场（后）.png';
saveas(gcf, fname, 'png');
figure;
hold on;
axis equal;
plot(xm, ym,'--');
quiver(lon12(inside_circle), lat12(inside_circle), vx(inside_circle), vy(inside_circle), 'r');
title('圆内各点速度大小');
xlabel('经度');
ylabel('纬度');
Fs = [];
vrot = sqrt(U.^2+V.^2);
Fs(:,:) = interp2(lat1,lon1,vrot,s2,s1,"spline");
figure
plot(s1,s3);
hold on;
plot(s1,Fs);
hold off;
title('风速剖面图');
xlabel('经度');
ylabel('风速');
figure
plot(s1,s3,'blue');
hold on;
plot(s1,Fs,'yellow');
hold on;
yc = 0.93*Vmax(13,:) - max(Fs);
ycc = 0.93*Vmax(13,:) - 14.9231;
index1 = find(Fs == max(Fs));
index2 = 34;
index3 = 69;
index4 = 87;
xx = s1(index3:index4);
yy = (ycc/(xx(1,:)-xx(19,:)))*(xx-xx(1,:))+ycc ;
xx1 = s1(index2:index1);
yy1 = (yc/(xx1(25,:)-xx1(1,:)))*(xx1-xx1(1,:));
Fs(index3:index4) = Fs(index3:index4) + yy;
Fs(index2:index1) = Fs(index2:index1) + yy1;
plot(s1,Fs,'r');
hold off;
title('风速校正图');
xlabel('经度');
ylabel('风速');
zFC = scatteredInterpolant(s2, s1, Fs, 'natural','nearest');
AA = zFC(lat1,lon1);
function distance = distance_on_unit_sphere(lon1, lat1, lon2, lat2)
lon1_rad = deg2rad(lon1);
lat1_rad = deg2rad(lat1);
lon2_rad = deg2rad(lon2);
lat2_rad = deg2rad(lat2);
delta_lon = lon2_rad - lon1_rad;
delta_lat = lat2_rad - lat1_rad;
a = sin(delta_lat/2)^2 + cos(lat1_rad) * cos(lat2_rad) * sin(delta_lon/2)^2;
c = 2 * atan2(sqrt(a), sqrt(1-a));
distance = c * 6372.8;
end
当我把ws里面的第三列数字79换成循环里的i4，出的图则是很多线连成一坨，应该是所有时刻的线都画出来了，但是都在同一张图上，而且还首尾相连。怎么把这些线按时刻区分并绘制在不同的图上。
for i5 = 1:k1
Jm(i5,i4-1) = interp2(lat1,lon1,ws(:,:,79),y1(i5,1),x1(i5,1),"spline");
end
