I want to develop a mathematical equations for twisted curved surface. The cross section is nacaa 4 series aerofoil along xz plane .I want to twist the surface along Y axis but I am unable to get the right plot.

syms s t t1=0.12; p=0.4; m=.02; a1=(m/p^2)*(2*p*s-s^2); b1=(m/(1-p)^2)*((1-2*p)+2*p*s-s^2); z_t=5*t1*(0.2969*sqrt(s)-0.1260*s-0.351) ; z_c=piecewise(0<=s<=p,a1,p<=s<=1,b1); diif=diff(z_c); theta=atan(diif) x_u=s-z_t*sin(theta); z_u=z_c+z_t*cos(theta); z_tx_L=s+z_t*sin(theta); z_L=z_c-z_t*cos(theta);

a=[s^3 s^2 s 1]; chord=20;

span=37.1; x1=s*chord; y1=t*span;

b=[1 1 1 1]; b1=transpose(b)*b;

n=[-1 3 -3 1;3 -6 3 0;-3 3 0 0;1 0 0 0]; c=[t^3 t^2 y 1]; d=transpose(c);

z1=20*z_u*a*n*b1*n*d; z2=20*z_L*a*n*b1*n*d ; y=y1 x=x1*cos(t)-z1*sin(t); z11=x1*sin(t)+z1*cos(t);

figure(1) fsurf(x, y, z11,[0,1,0,1])

figure(2) y=y1 x=x1*cos(t)-z2*sin(t); z22=x1*sin(t)+z2*cos(t); fsurf(x, y, -z22,[0,1,0,1]) figure (3) fsurf(x, y, z11, [0,1,0,1]) hold on fsurf(x, y,- z22,[0,1,0,1])

zlim([-5 5])

xlabel('chord'); ylabel('span'); zlabel('z'); title('top surface') box on