Matlab：使用一个for循环的时间序列结果作为第二个for循环的输入。？

Question

您如何使用来自一个 for 循环的时间序列结果 (As) 作为第二个 for 循环的输入？你可以在一个 .m 文件中有 2 个 for 循环吗？还是我将它们分成两个 .m 文件。如果我将它们分开，您如何在第二个循环的方程中调用第一个函数的输出？

x=load('C:....dat'); %Load the file with data of precipitation, temperature and snowpack
td=x(:,4);              %hour of the day
GlobalRad=x(:,5);          %Global Radiation
T=x(:,6);                %Air Temperature
P=x(:,7);                %Precipitation
Snow_obs=x(:,8);           %SWE Observed
SnowDepth=x(:,9);        %Snow Depth Observed
ts=x(:,10);                %start time of daylight on day d (t0)
te=x(:,11);               %end time of daylight on day d (t1)
dTd=x(:,12);              %difference between the max and min daily temperatures on day d

snow_sim(1)=0;
runoff=zeros(length(P));

%Parameters
Scf=1.3;                        %Snowfall correction factor
TT=1;                          %Threshold temperature
C=3.5;                           %Degree-day factor (mm day-1 C-1)(Ac)
Cfr=0.05;                        %Refreezing coefficient (use default value of 0.05)
Cwh=0.1;                         %Water holding capacity (use default value of 0.1)
B=0.05;         % Factor to convert the temp amplitude into a degree day factor


snow_sim(1)=0;                   %Simulated snowpack is 0 mm for day 1
snow_sim_water(1)=0;             %Liquid water in snowpack for day 1 the water is 0 mm


for t=1 : length(td);       %**FIRST time series loop**

    ln(t)=24-te(t)+ts(t);                    %length of the night
    Z(t)=2*((te(t)-ts(t))/(3.14*ln(t)));        %factor ensuring that daily mean vaules of As equals Ac

    if ts(t)<=td(t)<te(t);

      As(t)=C+(B*dTd(t)*(sin(3.14*((td(t)-ts(t))/(te(t)-ts(t))))));       %equation for time variant degree day factor secnario1

    else As(t)=C-(B*dTd(t)*Z(t));            %equation for time variant degree day factor senario 2

    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



for t=2 : length(P);             %**SECOND time series loop**

    if T(t)< TT;                 %If the temperature for day t is lower than the threshold value for melting (=below TT degrees) then refreezing of melt water will occur
        refreez(t)=Cfr*As*(TT-T(t));      %Equation for refreezing of meltwater
        if refreez(t) > snow_sim_water(t-1);     %If the refreezing for day t is larger than the water in the snowpack the day before, then the refreezing is limited
            refreez(t)= snow_sim_water(t-1);      %Can't freez more water than is accumulated in the snowpack
        end

        snow_sim(t)=P(t)*Scf+snow_sim(t-1)+refreez(t);     %The total simulated snowpack for any given day is the precipitation that day together with snow pack from day before and refreeze of that day.
        snow_sim_water(t)=snow_sim_water(t-1)-refreez(t);        %The total simulated amount of water in the snowpack is the water in the snowpack the day before minus the water refrozen the same day


    else %T(t) > TT                   %temperature above threshold temperature, snowmelt will occur
        Melt(t)=As*(T(t)-TT);           %Equations for melting rate of existing snowpack
        if Melt(t) > snow_sim(t-1);    %If the melting rate for day t is larger than the snowpack the day before, then the melting is limited
            Melt(t)= snow_sim(t-1);     %Because it can't melt more snow than is available
        end
        snow_sim(t)=snow_sim(t-1)-Melt(t);                  %Total simulated snow is the simulated snowpack for the day before minus the melted snow
        snow_sim_water(t)=snow_sim_water(t-1)+P(t)+Melt(t);       %Total water amount in snow is the water amount in snow for the day before plus the precipitation and the melted snow
        if Cwh*snow_sim(t) < snow_sim_water(t);                 %The snowpack can retain as much as 10% of its water equivalent, but not more
            runoff(t)=snow_sim_water(t)-0.1*snow_sim(t);        %if there is more liquid water, this goes to runoff (note:if there is no snowpack all water will go to runoff
            snow_sim_water(t)=0.1*snow_sim(t);
        end
    end
end



snow_sim_total=snow_sim+snow_sim_water;                   %The total simulated snowpack is the water in snow and the simulated snowpack

daynr=1:length(P);

Answer 1

您可以在脚本或函数中使用任意数量的 for 循环。尽可能避免使用它们，而是将代码矢量化

例如：

As=zeros(size(te));
diff1=te-ts; % this difference gets used multiple times
ln=24-diff1;
Z=2*diff1./(pi*ln);
ix=ts<=td<te;
As(ix)=C+(B*dTd(ix).*sin(pi*((td(ix)-ts(ix))./diff1(ix))));
As(~ix)=C-(B*dTd(~ix).*Z(~ix));

与第一个 for 循环的作用相同，但运行速度要快得多。

看起来 As 是一个向量，但在你的第二个循环中，你把它当作一个标量，似乎你错过了那里的时间索引（即

refreez(t)=Cfr*As*(TT-T(t)); 

应该是

refreez(t)=Cfr*As(t)*(TT-T(t));

我认为。

这个循环可能不容易向量化，因为它内部的计算取决于以前的结果，所以可以保持原样并修复错误。 修复第二个循环中的错误后，在同一个脚本中运行两个循环应该没有问题。