function [panel, Vol] = DeterminePanelGeometry(inputgeo, figs)
% find coordinates for horseshoe vortices and control points and plot
% Aircraft-Fixed coordinate system:
% x: forward along fuselage axis
% y: out starboard wing
% z: down
% Local Profile coordinate system:
% x: along chord from leading edge toward trailing edge
% z: up
plot_on = 1;
%% constants
halfSpan=inputgeo.b/2; % half wingspan
sweep=inputgeo.sweep; % wing sweep angle in radians
dihedral=inputgeo.dih; %dihedral angle in radians
nc=inputgeo.nc; %number of panels on camber line, upper and lower surfaces
ns = inputgeo.ns; % number of span segments per side
M = inputgeo.M; %Freestream mach number
beta = sqrt(1-M^2); %Prandtl-Glauert correction
%% Root Airfoil data
y=0; % span station
chord=inputgeo.c_r; % chord length
alphaRoot=inputgeo.i_r; % geometric angle of incidence at root in radians
% call function
[Croot,Troot,Aroot]=NacaCoord(inputgeo.root,y,nc); % Croot = nondimensional camber line, Troot = nondimensional surface;
%Croot coordinates = [x location in fraction chord, y location in fraction
% chord, z location in fraction chord]
%Troot coordinates = [x location in fraction chord, y location in fraction
% chord, z location in fraction chord]
%% Scale to chord length
Croot(:,1)=chord*Croot(:,1);Croot(:,3)=chord*Croot(:,3);
Troot(:,1)=chord*Troot(:,1);Troot(:,3)=chord*Troot(:,3);
Aroot=Aroot*chord^2;
%Calculate slope of camber line for each panel at the root
%Calculate the chord along each elemental panel
dzdxRoot = zeros(1,nc);
ccRoot = zeros(1,nc);
if nc == 1 %If there is only one chordwise element
i=1;
dzdxRoot(i)=-(Croot(i+1,3)-Croot(i,3))/(Croot(i+1,1)-Croot(i,1));
ccRoot(i)=0.75*(Croot(i+1,1)-Croot(i,1));
else %if there are more than one chordwise elements
for i=1:nc
dzdxRoot(i)=(Croot(i+1,3)-Croot(i,3))/(Croot(i+1,1)-Croot(i,1));
if i==nc
ccRoot(i)=0.75*(Croot(i+1,1)-Croot(i,1));
else
ccRoot(i)=0.75*(Croot(i+1,1)-Croot(i,1))+0.25*(Croot(i+2,1)-Croot(i+1,1));
end
end
end
if plot_on
axes(figs.root); cla;
plot(Croot(:,1)*-3.2808,Croot(:,3)*-3.2808,'r') %Convert to feet; plot with z-axis positive up and x-axis to the left
hold on
plot(Troot(:,1)*-3.2808,Troot(:,3)*-3.2808,'b') %Convert to feet; plot with z-axis positive up and x-axis to the left
axis equal
title('Root airfoil')
end
%% Tip Airfoil data
y=-halfSpan; % span station; left wing
chord=inputgeo.taper*inputgeo.c_r; % chord length
alphaTip=inputgeo.i_r+inputgeo.twist; % geometric angle of incidence at tip in radians
%call function
[Ctip,Ttip,Atip]=NacaCoord(inputgeo.tip,y,nc);
%% scale to chord length
Ctip(:,1)=chord*Ctip(:,1);Ctip(:,3)=chord*Ctip(:,3);
Ttip(:,1)=chord*Ttip(:,1);Ttip(:,3)=chord*Ttip(:,3);
Atip=Atip*chord^2;
%Calculate slope of camber line for each panel at the tip
%Calculate the chord along each elemental panel
dzdxTip = zeros(1,nc);
ccTip = zeros(1,nc);
if nc == 1 %If there is only one chordwise element
dzdxTip(i)=-(Ctip(i+1,3)-Ctip(i,3))/(Ctip(i+1,1)-Ctip(i,1));
ccTip(i)=0.75*(Ctip(i+1,1)-Ctip(i,1));
else %if there are more than one chordwise elements
for i=1:nc
dzdxTip(i)=(Ctip(i+1,3)-Ctip(i,3))/(Ctip(i+1,1)-Ctip(i,1));
if i==nc
ccTip(i)=0.75*(Ctip(i+1,1)-Ctip(i,1));
else
ccTip(i)=0.75*(Ctip(i+1,1)-Ctip(i,1))+0.25*(Ctip(i+2,1)-Ctip(i+1,1));
end
end
end
if plot_on
axes(figs.tip); cla;
plot(Ctip(:,1)*-3.2808,Ctip(:,3)*-3.2808,'r') %Convert to feet; plot with z-axis positive up and x-axis to the left
hold on
plot(Ttip(:,1)*-3.2808,Ttip(:,3)*-3.2808,'b') %Convert to feet; plot with z-axis positive up and x-axis to the left
axis equal
title('Tip airfoil')
end
%% Root: Transform to Aircraft-Fixed Coordinate System:
% skin:
%Skin(number of spanwise sections, number of chordwise sections, coordinates in 3 dimensions (x,y,z))
Skin(1,1:2*nc+1,1:3) = Troot*[cos(alphaRoot) 0 sin(alphaRoot); 0 1 0; -sin(alphaRoot) 0 cos(alphaRoot)];
% Mean camber line:
Camber(1,1:nc+1,1:3) = Croot*[cos(alphaRoot) 0 sin(alphaRoot); 0 1 0; -sin(alphaRoot) 0 cos(alphaRoot)];
%% Tip: Transform to Aircraft-Fixed Coordinate System:
% sweep and dihedral
% Negative sign occurs because local axes are opposite to global
xLETip = -halfSpan*(tan(sweep));% x coordinate of Leading Edge at Tip
zLETip = -halfSpan*(tan(dihedral));% z coordinate of Leading Edge at Tip
% Skin:
LETip=ones(2*nc+1,1)*[xLETip 0 zLETip];
Skin(ns+1,1:2*nc+1,1:3) = Ttip*[cos(alphaTip) 0 sin(alphaTip); 0 1 0; -sin(alphaTip) 0 cos(alphaTip)]+LETip;
% Mean camber line
LETip=ones(nc+1,1)*[xLETip 0 zLETip];
Camber(ns+1,1:nc+1,1:3) = Ctip*[cos(alphaTip) 0 sin(alphaTip); 0 1 0; -sin(alphaTip) 0 cos(alphaTip)]+LETip;
%% Intermediate stations
for k=1:ns-1 % k = 0 corresponds to root, k = ns corresponds to tip
eta=k/ns;
for i=1:nc+1 % along camber
Camber(k+1,i,:) = Camber(1,i,:)+eta*(Camber(ns+1,i,:)-Camber(1,i,:));
end
for i=1:2*nc+1 %along surface
Skin(k+1,i,:) = Skin(1,i,:)+eta*(Skin(ns+1,i,:)-Skin(1,i,:));
end
end
if plot_on
%Plot wing surface
axes(figs.fig); cla; hold on; axis equal;
for k=1:ns+1
plot3(Skin(k,:,1)*-3.2808,Skin(k,:,2)*3.2808,Skin(k,:,3)*-3.2808); %Convert to feet; plot with z-axis positive up and x-axis to the left
plot3(Skin(k,:,1)*-3.2808,-Skin(k,:,2)*3.2808,Skin(k,:,3)*-3.2808); %Convert to feet; plot with z-axis positive up and x-axis to the left
plot3(Camber(k,:,1)*-3.2808,Camber(k,:,2)*3.2808,Camber(k,:,3)*-3.2808,'r') %Convert to feet; plot with z-axis positive up and x-axis to the left
end
for i=1:size(Skin,2)
plot3(Skin(:,i,1)*-3.2808,Skin(:,i,2)*3.2808,Skin(:,i,3)*-3.2808);
plot3(Skin(:,i,1)*-3.2808,-Skin(:,i,2)*3.2808,Skin(:,i,3)*-3.2808);
end
end
%Determine panel properties
twist = zeros(1,ns);
dzdx = zeros(ns,nc);
cc = zeros(ns,nc);
CP = zeros(ns,nc,3);
BV = zeros(ns,nc,3);
BV1 = zeros(ns,nc,3);
BV2 = zeros(ns,nc,3);
BVhalfspan = zeros(ns,nc);
sweep_c4 = zeros(ns,nc);
for k = 1:ns
%Determine local twist at each panel, dz/dx at control points, and
%chord along left trailing leg of elemental panel
eta=k/(ns+1);
twist(k)=alphaRoot+eta*(alphaTip-alphaRoot);
dzdx(k,:)=dzdxRoot+eta*(dzdxTip-dzdxRoot);
cc(k,:)=ccRoot+eta*(ccTip-ccRoot); %Validated
%Determine control point coordinates, coordinates to center of bound
%vortex, half span of bound vortex, and quarter-chord sweep angle
for i = 1:nc
CP(k,i,:)=0.5*(Camber(k,i,:)+Camber(k+1,i,:))+0.75*(0.5*(Camber(k,i+1,:)+Camber(k+1,i+1,:))-0.5*(Camber(k,i,:)+Camber(k+1,i,:)));
BV(k,i,:)=0.5*(Camber(k,i,:)+Camber(k+1,i,:))+0.25*(0.5*(Camber(k,i+1,:)+Camber(k+1,i+1,:))-0.5*(Camber(k,i,:)+Camber(k+1,i,:)));
BV1(k,i,:)=0.75*Camber(k+1,i,:)+0.25*Camber(k+1,i+1,:); %Left bound vortex coordinate
BV2(k,i,:)=0.75*Camber(k,i,:)+0.25*Camber(k,i+1,:); %Right bound vortex coordinate
BVhalfspan(k,i)= 0.5*(0.75*norm(shiftdim(Camber(k,i,2:3)-Camber(k+1,i,2:3)))+0.25*norm(shiftdim(Camber(k,i+1,2:3)-Camber(k+1,i+1,2:3))));
%Only consider [Y,Z] distances due to definition of halfspan in
%NASA paper; considering X distance as well makes halfspan increase
%greatly as you increase sweep and invalidates calculations
%Validated BVhalfspan calculation through (1:ns-1,1:nc)
sweep_c4(k,i)=atan((Camber(k,i,1)-Camber(k+1,i,1))/(Camber(k,i,2)-Camber(k+1,i,2))); %Validated
end
end
%Prandtl-Glauert corrections
CPprime = CP(:,:,1)/beta;
BVprime = BV(:,:,1)/beta;
BV1prime = BV1(:,:,1)/beta;
BV2prime = BV2(:,:,1)/beta;
sweep_c4_prime = atan(tan(sweep_c4)/beta);
for k = 1:ns
for i = 1:nc
panel(k,i).CP=shiftdim(CP(k,i,:)); %Removing singleton dimensions
panel(k,i).BV=shiftdim(BV(k,i,:)); %Removing singleton dimensions
panel(k,i).BV1=shiftdim(BV1(k,i,:)); %Removing singleton dimensions
panel(k,i).BV2=shiftdim(BV2(k,i,:)); %Removing singleton dimensions
panel(k,i).dzdx=dzdx(k,i);
panel(k,i).twist=twist(k);
panel(k,i).s=BVhalfspan(k,i);
panel(k,i).sweep=sweep_c4(k,i);
panel(k,i).sweepprime=sweep_c4_prime(k,i);
panel(k,i).CPprime=CPprime(k,i);
panel(k,i).BVprime=BVprime(k,i);
panel(k,i).BV1prime=BV1prime(k,i);
panel(k,i).BV2prime=BV2prime(k,i);
panel(k,i).cc=cc(k,i);
end
end
%Determine wing volume
Vol = 1/3*(Aroot +sqrt(Aroot*Atip) + Atip)*inputgeo.b;
