I have to use the source panel method to calculate the Cp about a given airfoil. I have done my best to follow my professors instructions, but my Matlab code is giving me incorrect results. Can someone help me troubleshoot this code?

oad Original Airfoil Points

data $=$ load$(\text{'}$Airfoil$\_$Data.mat'$)$;

x $=$ data.unnamed$($:$,1)$;

y $=$ data.unnamed$($:$,2)$;

V$\_0=1$;

alpha $=$ deg$2$rad$([-5])$;

Interpolate new points

ds $=$ sqrt$($diff$($x$).^2+$ diff$($y$).^2)$;

s $=[0$; cumsum$($ds$)]$;

N $=128$;

s$\_$new $=$ linspace$(0,$ s$($end$),$ N $+1)$;

x$\_$new $=$ interp$1($s$,$ x$,$ s$\_$new, 'linear'$)$;

y$\_$new $=$ interp$1($s$,$ y$,$ s$\_$new, 'linear'$)$;

points$=[$x$\_$new$($:$),$ y$\_$new$($:$)]$;

figure;

plot$($x$,$ y$,\text{'}-$o$\text{'},$ 'DisplayName', 'Original Airfoil'$)$;

hold on;

plot$($x$\_$new, y$\_$new, $\text{'}-$x$\text{'},$ 'DisplayName', 'Interpolated Airfoil'$)$;

legend;

xlabel$(\text{'}$x$\text{'})$;

ylabel$(\text{'}$y$\text{'})$;

title$(\text{'}$Original and Interpolated Airfoil'$)$;

grid on;

axis equal;

Control Points

x$\_$control $=$ zeros$($N$,1)$;

y$\_$control $=$ zeros$($N$,1)$;

for i $=1$:N

next$\_$point $=$ mod$($i$,$ N$)+1$;

x$\_$control$($i$)=($points$($i$,1)+$ points$($next$\_$point, $1))/2$;

y$\_$control$($i$)=($points$($i$,2)+$ points$($next$\_$point, $2))/2$;

end

controlPoints $=[$x$\_$control, y$\_$control$]$;

$\%$ Plot the airfoil and control points

figure;

plot$($points$($:$,1),$ points$($:$,2),\text{'}-$o$\text{'},$ 'DisplayName', 'Airfoil'$)$; $\%$ Airfoil points

hold on;

plot$($controlPoints$($:$,1),$ controlPoints$($:$,2),\text{'}$xr$\text{'},$ 'DisplayName', 'Control Points'$)$; $\%$ Control points

legend;

xlabel$(\text{'}$x$\text{'})$;

ylabel$(\text{'}$y$\text{'})$;

title$(\text{'}$Airfoil and Control Points'$)$;

grid on;

axis equal;

Panel Lengths

deltaS $=$ zeros$($N$,1)$;

for i $=1$:N

next$\_$point $=$ mod$($i$,$ N$)+1$;

dx $=$ points$($next$\_$point, $1)-$ points$($i$,1)$;

dy $=$ points$($next$\_$point, $2)-$ points$($i$,2)$;

deltaS$($i$)=$ sqrt$($dx$^2+$dy$^2)$;

end

Panel Angle w$/$ respect to x axis

theta $=$ zeros$($N$,1)$;

for i $=1$:N

next$\_$point $=$ mod$($i$,$ N$)+1$;

dx $=$ points$($next$\_$point, $1)-$ points$($i$,1)$;

dy $=$ points$($next$\_$point, $2)-$ points$($i$,2)$;

theta$($i$)=$ atan$2($dy$,$ dx$)$;

if theta$($i$)<0$

theta$($i$)=$ theta$($i$)+2*$pi;

end

end

Length from control point to control point

r $=$ zeros$($N$,$ N$)$;

for i $=1$:N

for j $=1$:N

dzeta $=$ controlPoints$($j$,1)-$ controlPoints$($i$,1)$;

deta $=$ controlPoints$($j$,2)-$ controlPoints$($i$,2)$;

r$($i$,$ j$)=$ sqrt$($dzeta$^2+$ deta$^2)$;

end

end

Angle between control point i and j w$/$ respect to the x axis

phi $=$ zeros$($N$,$ N$)$;

for i $=1$:N

for j $=1$:N

dzeta $=$ controlPoints$($j$,1)-$ controlPoints$($i$,1)$;

deta $=$ controlPoints$($j$,2)-$ controlPoints$($i$,2)$;

phi$($i$,$ j$)=$ atan$2($deta$,$ dzeta$)$;

if phi$($i$,$ j$)<0$

phi$($i$,$ j$)=$ phi$($i$,$ j$)+2*$pi;

end

end

end

Circulation and Normal influence coefficients

C $=$ zeros$($N$,$ N$)$;

C$\_$bar $=$ zeros$($N$,$ N$)$;

for i $=1$:N

for j $=1$:N

if i ~$=$ j

C$($i$,$ j$)=($sin$($theta$($i$)-$ phi$($i$,$ j$))/(2*$ pi $*$ r$($i$,$ j$)))*$ deltaS$($j$)$;

C$\_$bar$($i$,$ j$)=($cos$($theta$($i$)-$ phi$($i$,$ j$))/(2*$ pi $*$ r$($i$,$ j$)))*$ deltaS$($j$)$;

else

C$($i$,$ j$)=0$;

C$\_$bar$($i$,$ j$)=0$;

end

end

end

Source Strength

RHS $=$ zeros$($N$,1)$;

for i $=1$:N

RHS$($i$)=$ V$\_0*$ sin$($theta$($i$)-$ alpha$)$;

end

A $=$ zeros$($N$,$ N$)$;

for i $=1$:N

for j $=1$:N

if i $==$ j

A$($i$,$ j$)=1/2$;

else

A$($i$,$ j$)=$ C$($i$,$ j$)$;

end

end

end

q $=$ A $\backslash$ RHS;

Tangential Velocity

Vt $=$ zeros$($N$,1)$;

for i $=1$:N

Vt$($i$)=$ V$\_0*$ cos$($theta$($i$)-$ alpha$)$;

for j $=1$:N

if i ~$=$ j

Vt$($i$)=$ Vt$($i$)-$ q$($j$)*$ C$\_$bar$($i$,$ j$)$;

end

end

end

Pressure Coefficient

$\%$ Ensure Cp is computed

Cp $=1-($Vt $/$ V$\_0).^2$; $\%$ Vectorized computation of Cp

$\%$ Plot Cp distribution

figure;

plot$($controlPoints$($:$,1),$ Cp$,\text{'}-$o$\text{'})$; $\%$ Cp vs$.$ x$-$coordinates of control points

xlabel$(\text{'}$x $($Chordwise Position$)\text{'})$;

ylabel$(\text{'}$C$\_$p$\text{'})$;

title$(\text{'}$Pressure Coefficient Distribution'$)$;

grid on;