2015-11-10 33 views
0

如何將一個for循環中的時間序列結果(As)用作第二個for循環的輸入?你能在一個.m文件中有2個for-loops嗎?或者我將它們分成兩個.m文件。如果我將它們分開,那麼如何在第二個迴路方程中調用第一個函數的輸出?Matlab:使用一個for-loop的時間序列結果作爲第二個for-loop的輸入。

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); 

回答

0

您可以在腳本或函數中儘可能多地使用for-loops。避免他們在這裏你可以和你的矢量化代碼,而不是

例如爲:

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循環,但將大大更快地運行。

看來,作爲一個載體,但在你的第二個循環,你把它當作一個標量,好像你缺少時間指數出現(即

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

應該

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

我想。

這個循環可能是不容易的向量化,因爲在它的計算依賴於以前的結果,所以可以離開它,因爲它是和剛修復錯誤。 一旦您修復了第二個循環中的錯誤,您應該在同一個腳本中同時運行兩個循環。

+0

謝謝Felix!我感謝幫助! – Wesser