Why is my code getting stuck in interp1?

2 次查看（过去 30 天）
Megan Fiddleton
Megan Fiddleton 2017-12-6
评论： Walter Roberson 2017-12-6
I have a final project where I must calculate the power production of wind turbines in a city. Here is my code and its functions. When I run the code it just runs forever, and when i press pause, I see that it is always in the interp1 function.
MAIN FILE:
clc
clear
close all
% The data is saved as another .mat file. load the file
load('WindSpeed.mat')
% Convert this data to m/s and round to nearest integer
ws = (WindSpeed*0.27777);
B = 3;
% Calculate the wind speed at the hub height assigned
H = 10;
h = 95;
V = ws;
% To calculate windspeed v0 at any height given the windspeed at a certain height
v0 = V*(h/H)^(1/7);
% define the population
population = 316701;
% Load in the omega values
w = load('omega.dat');
% Convert rpm to rad/s
w = w*(2*pi/60);
% Load in the radial positions
r = load('radius_medium.dat');
% Load in the twist angle theta
theta = load('twist.dat')*(pi/180); %converting to radians
chord = load('chord.dat');
DU30 = load('DU30.dat');
DU30(:,1) = DU30(:,1)*(pi/180);
% Initialize a and a', phi, and ct are needed to find a and a' values
aprime = zeros(25,17);
a = zeros(25,17);
phi = zeros(25,17);
Ct = zeros(25,17);
for i = 1:25
[a(i,:),aprime(i,:),phi(i,:),Ct(i,:)]=A_func(i,w,r,theta,chord,DU30);
end
A_func:
function [ a,aprime,Ct,phi ] = A_func(i,w,r,theta,chord,DU30)
%This function is to calculate a and a' values
Avec=zeros(size(r)); %makes a 17x1 matrix
Aprimevec=zeros(size(r));
phivec=zeros(size(r));
% Initialize a Ct vec, as it must be used later on unlike Cn.
Ctvec=zeros(size(r));
for n = 1:length(r)
aold=1;aprimeold=1;a=2;aprime=2; %Initialize values to be different so it will enter the while loop
while(abs(aold-a)>=1e-6 || abs(aprimeold-aprime)>=1e-6)
aold=a;
aprimeold=aprime;
phi=atan(((1-aold)*i)/((1+aprimeold)*(w(i)*r(n))));
alpha = phi-theta(n);
% computing cl and cd
Cl=interp1(DU30(:,1),DU30(:,2),alpha);
Cd=interp1(DU30(:,1),DU30(:,3),alpha);
% now able to compute Cn and Ct
Cn = Cl*cos(phi)+Cd*sin(phi);
Ct = Cl*sin(phi)-Cd*cos(phi);
% Must define sigma for a and a' calculations
sig = chord(n)*3/(2*pi*r(n));
a = (1/((4*sin(phi)^2)/(sig*Cn)+1));
aprime = 1/((4*sin(phi)*cos(phi))/(sig*Ct)-1);
% Apply the Gaulert correction for when a >ac
if(a>=0.2)
% Define K
K = (4*sin(phi)^2)/(sig*Cn);
ac = 0.2;
a = 0.5*(2+K*(1-2*ac)-sqrt((K*(1-2*ac)+2)^2+4*(K*ac^2-1)));
end
end
Avec(n)=aold; Aprimevec(n)=aprimeold; phivec(n)= phi; Ctvec(n)=Ct;
end
a = Avec; aprime = Aprimevec; phi = phivec; Ct = Ctvec;
end
  2 个评论
ANKUR KUMAR
ANKUR KUMAR 2017-12-6
Mention the error which you are getting.
Megan Fiddleton
Megan Fiddleton 2017-12-6
It's not an error im getting. The code just runs forever. When I press pause, I see that it is always in the interp1 function. It never gets out of it.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Star Strider
Star Strider 2017-12-6
I cannot follow your code, and I don’t have the data. The interp1 function is likely not the problem, unless it’s actually throwing an error.
See the documentation on Debugging (link), since you’re the only one amongst us who has access to everything. The debugging steps will allow you to examine your variables to see what the problem is. The dbstop (link) function is particularly helpful.
Walter Roberson
Walter Roberson 2017-12-6
编辑：Walter Roberson 2017-12-6
It just is not efficient to keep calling interp1() to interpolate single points.
I suggest you do a first pass of interpolating every 0.05 to fill the gaps. After that, given any particular alpha, you can subtract the minimum and divide by 0.05 and floor() and add 1 to figure out which index number in the table to use. Then use that entry in the table and the next entry to do a simple linear interpolation to get your Cl and Cd values.
  1 个评论
Walter Roberson
Walter Roberson 2017-12-6
floor "round towards negative infinity". For example, floor(17.39) is 17, floor(3.99999999999999) is 3, floor(2) is 2 . floor() is the largest integer that does not exceed the input.
There is a corresponding function ceil() which is the smallest integer that equals or exceeds the input. For example ceil(17.39) is 18, ceil(3.99999999999999) is 4, ceil(2) is 2.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Logical 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by