Why is my code getting stuck in interp1?

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

[ANKUR KUMAR]

Mention the error which you are getting.

[Megan Fiddleton]

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.

Answer 1

[Star Strider]

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.

Answer 2

[Walter Roberson]

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.

[Walter Roberson]

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.