summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--gui.py89
-rw-r--r--wing_scripts/eye_beam_example.m70
-rw-r--r--wing_scripts/get_dp.m4
-rw-r--r--wing_scripts/get_ds.m20
-rw-r--r--wing_scripts/get_int.m35
-rw-r--r--wing_scripts/get_z.m34
-rw-r--r--wing_scripts/my_progress.m459
-rw-r--r--wing_scripts/stringersBeamExample.m47
-rw-r--r--wing_scripts/wingAnalysis_190422.m579
9 files changed, 1298 insertions, 39 deletions
diff --git a/gui.py b/gui.py
index cd036b1..a1b8d8e 100644
--- a/gui.py
+++ b/gui.py
@@ -13,9 +13,10 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
-import creator
-import evaluator
-import generator
+from tools import creator, evaluator, generator
+# import creator
+# import evaluator
+# import generator
import tkinter as tk
import tkinter.ttk as ttk
@@ -23,6 +24,51 @@ from matplotlib.backends.backend_tkagg import (
FigureCanvasTkAgg, NavigationToolbar2Tk)
+class MainWindow(tk.Frame):
+ """Main editor window."""
+
+ def __init__(self, *args, **kwargs):
+ tk.Frame.__init__(self, *args, **kwargs)
+ root = tk.Tk()
+ root.wm_title('MAE 154B - Airfoil Design, Evaluation, Optimization')
+
+ # self.button = tk.Button(self, text="Create new window",
+ # command=self.create_window)
+ # self.button.pack(side="top")
+ frame_1 = ttk.Frame(root)
+ l_naca, e_naca = new_field(frame_1, 'naca')
+ l_chord, e_chord = new_field(frame_1, 'chord')
+ l_semi_span, e_semi_span = new_field(frame_1, 'semi_span')
+ af = generator.default_airfoil()
+ # Graph window
+ frame_2 = ttk.Frame(root)
+ fig, ax = creator.plot_geom(af, False)
+ plot = FigureCanvasTkAgg(fig, frame_2)
+ # plot.draw()
+ toolbar = NavigationToolbar2Tk(plot, frame_2)
+ # toolbar.update()
+
+ l_naca.grid(row=0, column=0)
+ e_naca.grid(row=0, column=1, padx=4)
+ # b_naca.grid(row=0, column=2)
+ l_chord.grid(row=1, column=0)
+ e_chord.grid(row=1, column=1, padx=4)
+ l_semi_span.grid(row=2, column=0, padx=4)
+ e_semi_span.grid(row=2, column=1, padx=4)
+ frame_1.pack(side=tk.LEFT)
+ # Graph window
+ plot.get_tk_widget().pack(expand=1, fill=tk.BOTH)
+ toolbar.pack()
+ frame_2.pack(side=tk.LEFT)
+
+ def create_window(self):
+ self.counter += 1
+ window = tk.Toplevel(self)
+ window.wm_title("Window #%s" % self.counter)
+ label = tk.Label(window, text="This is window #%s" % self.counter)
+ label.pack(side="top", fill="both", expand=True, padx=100, pady=100)
+
+
def new_field(parent, name):
"""Add a new user input field."""
@@ -46,40 +92,5 @@ def set_semi_span(name):
print(semi_span)
-root = tk.Tk()
-root.wm_title('MAE 154B - Airfoil Design, Evaluation, Optimization')
-# root.geometry('1000x400')
-
-# Object definition
-# User inputs
-frame_1 = ttk.Frame(root)
-l_naca, e_naca = new_field(frame_1, 'naca')
-l_chord, e_chord = new_field(frame_1, 'chord')
-l_semi_span, e_semi_span = new_field(frame_1, 'semi_span')
-af = generator.default_airfoil()
-# Graph window
-frame_2 = ttk.Frame(root)
-fig, ax = creator.plot_geom(af, False)
-plot = FigureCanvasTkAgg(fig, frame_2)
-# plot.draw()
-toolbar = NavigationToolbar2Tk(plot, frame_2)
-# toolbar.update()
-
-# Layout
-# User input
-l_naca.grid(row=0, column=0)
-e_naca.grid(row=0, column=1, padx=4)
-# b_naca.grid(row=0, column=2)
-l_chord.grid(row=1, column=0)
-e_chord.grid(row=1, column=1, padx=4)
-l_semi_span.grid(row=2, column=0, padx=4)
-e_semi_span.grid(row=2, column=1, padx=4)
-frame_1.pack(side=tk.LEFT)
-# Graph window
-plot.get_tk_widget().pack(expand=1, fill=tk.BOTH)
-toolbar.pack()
-frame_2.pack(side=tk.LEFT)
-
# plot.get_tk_widget().pack()
-
-root.mainloop()
+MainWindow().mainloop()
diff --git a/wing_scripts/eye_beam_example.m b/wing_scripts/eye_beam_example.m
new file mode 100644
index 0000000..70b4d92
--- /dev/null
+++ b/wing_scripts/eye_beam_example.m
@@ -0,0 +1,70 @@
+% Bending/Shear stress example
+close all;
+
+length = 20; % in
+force = 10000; %lbs
+
+%eye-beam dimensions
+
+max_width = 4; % in
+min_width = 1; % in
+y_max = 4; % in
+center_y = 2; % in
+
+
+%max bending moment at the root...
+
+M = force*length;
+
+I = min_width*(2*center_y)^3/12 + 2*( max_width*(y_max-center_y)^3/12 + ...
+ max_width*(y_max-center_y)*((y_max+center_y)/2)^2);
+
+sigma_max = M * y_max / I;
+
+
+% solve for shear stress distribution
+% V / (I * t) * int(y*da)
+
+% Point 1: evaluated at location just before thickness changes from 4 to 1 in
+tempCoeff = force / (I * max_width);
+int_y_da = ((y_max+center_y)/2) * max_width*(y_max-center_y);
+shear_1 = tempCoeff*int_y_da;
+
+% Point 2: evaluated at location just after thickness changes from 4 to 1 in
+tempCoeff = force / (I * min_width);
+shear_2 = tempCoeff*int_y_da;
+
+
+% Point 3: evaluated at center of beam
+tempCoeff = force / (I * min_width);
+int_y_da = (center_y/2) * min_width*center_y;
+shear_3 = shear_2+tempCoeff*int_y_da;
+
+%evaluating continous integral for width of 4..
+int_y_da_4 = force / (I * max_width)*4*(y_max^2/2 - (center_y:.1:y_max).^2/2);
+
+%evaluating continous integral for width of 1..
+int_y_da_1 = shear_2 + force / (I * min_width)*1*(center_y^2/2 - (0:.1:center_y).^2/2);
+
+figure; grid on; hold on;set(gcf,'color',[1 1 1]);
+plot(int_y_da_4,center_y:.1:y_max,'linewidth',2)
+plot(int_y_da_1,0:.1:center_y,'linewidth',2)
+plot(int_y_da_1,0:-.1:-center_y,'linewidth',2)
+plot(int_y_da_4,-center_y:-.1:-y_max,'linewidth',2)
+plot([shear_1 shear_2],[center_y center_y],'linewidth',2)
+plot([shear_1 shear_2],[-center_y -center_y],'linewidth',2)
+
+plot(shear_1,center_y,'o')
+plot(shear_2,center_y,'o')
+plot(shear_3,0,'o')
+plot(shear_2,-center_y,'o')
+plot(shear_1,-center_y,'o')
+
+xlabel('shear stress (lb/in^2)','fontsize',16,'fontweight','bold');ylabel('Distance from Center (in)','fontsize',16,'fontweight','bold')
+set(gca,'FontSize',16,'fontweight','bold');
+
+
+figure; grid on; hold on;set(gcf,'color',[1 1 1]);
+plot([0 4 4 2.5 2.5 4 4 0 0 1.5 1.5 0 0],[4 4 2 2 -2 -2 -4 -4 -2 -2 2 2 4],'linewidth',2)
+
+
diff --git a/wing_scripts/get_dp.m b/wing_scripts/get_dp.m
new file mode 100644
index 0000000..2a3281d
--- /dev/null
+++ b/wing_scripts/get_dp.m
@@ -0,0 +1,4 @@
+function z = get_dp(xDist,zDist,Vx,Vz,Ix,Iz,Ixz,A)
+
+denom = (Ix*Iz-Ixz^2);
+z = -A*xDist*(Ix*Vx-Ixz*Vz)/denom - A*zDist*(Iz*Vz-Ixz*Vx)/denom;
diff --git a/wing_scripts/get_ds.m b/wing_scripts/get_ds.m
new file mode 100644
index 0000000..2f0eb9d
--- /dev/null
+++ b/wing_scripts/get_ds.m
@@ -0,0 +1,20 @@
+function ds = get_ds(xi,xf,u)
+
+dist = 0;
+numSteps = 10;
+dx = (xf-xi)/numSteps;
+z0 = get_z(xi,u);
+x0 = xi;
+for i=1:10
+ tempX = x0+dx;
+ if tempX > 0
+ tempZ = get_z(tempX,u);
+ else
+ tempZ = 0;
+ end
+ dist = dist + (dx^2+(tempZ-z0)^2)^.5;
+ z0 = tempZ;
+ x0 = tempX;
+end
+
+ds =dist; \ No newline at end of file
diff --git a/wing_scripts/get_int.m b/wing_scripts/get_int.m
new file mode 100644
index 0000000..edbfda3
--- /dev/null
+++ b/wing_scripts/get_int.m
@@ -0,0 +1,35 @@
+function z = get_int(xi,xf,u)
+
+M = 0.02;
+P = 0.4;
+T = 0.12;
+a0 = 0.2969;
+a1 = -0.126;
+a2 = -0.3516;
+a3 = 0.2843;
+a4 = -0.1015;
+
+
+%evaluate the integral of camber line, depending on xi and xf related to P
+
+if xf <P
+ intCamb = M/P^2*(2*P*xf^2/2 - xf^3/3) - M/P^2*(2*P*xi^2/2 - xi^3/3);
+elseif xi<P
+ intCamb = (M/(1-P)^2)*((1 - 2*P)*xf +2*P*xf^2/2 - xf^3/3) - (M/(1-P)^2)*((1 - 2*P)*P +2*P*P^2/2 - P^3/3);
+ intCamb = intCamb + M/P^2*(2*P*P^2/2 - P^3/3) - M/P^2*(2*P*xi^2/2 - xi^3/3);
+else
+ intCamb = (M/(1-P)^2)*((1 - 2*P)*xf +2*P*xf^2/2 - xf^3/3) - (M/(1-P)^2)*((1 - 2*P)*xi +2*P*xi^2/2 - xi^3/3);
+end
+
+% do integral on thickness line
+%z_thickness = (T/0.2)*(a0*x^.5+a1*x+a2*x^2+a3*x^3+a4*x^4);
+
+intThickness = (T/0.2)*(a0*xf^1.5/1.5 + a1*xf^2/2 + a2*xf^3/3 + a3*xf^4/4 +a4*xf^5/5);
+intThickness = intThickness - (T/0.2)*(a0*xi^1.5/1.5 + a1*xi^2/2 + a2*xi^3/3 + a3*xi^4/4 +a4*xi^5/5);
+
+% combine both integral results to get total integral
+if u == 1
+ z = intCamb + intThickness;
+else
+ z = abs(intCamb - intThickness);
+end \ No newline at end of file
diff --git a/wing_scripts/get_z.m b/wing_scripts/get_z.m
new file mode 100644
index 0000000..5387b52
--- /dev/null
+++ b/wing_scripts/get_z.m
@@ -0,0 +1,34 @@
+function z = get_z(x,u)
+
+
+
+if (x < 0 )
+ disp('invalid X')
+end
+
+M = 0.02;
+P = 0.4;
+T = 0.12;
+a0 = 0.2969;
+a1 = -0.126;
+a2 = -0.3516;
+a3 = 0.2843;
+a4 = -0.1015;
+
+if x <P
+ z_camber = M/P^2*(2*P*x - x^2);
+else
+ z_camber = (M/(1-P)^2)*(1 - 2*P +2*P*x - x^2);
+end
+
+%z_camber = M/P^2*(2*P*x - x^2);
+z_thickness = (T/0.2)*(a0*x^.5+a1*x+a2*x^2+a3*x^3+a4*x^4);
+
+if u==1
+ z = z_camber + z_thickness;
+else
+ z = z_camber - z_thickness;
+end
+
+
+
diff --git a/wing_scripts/my_progress.m b/wing_scripts/my_progress.m
new file mode 100644
index 0000000..e87ef23
--- /dev/null
+++ b/wing_scripts/my_progress.m
@@ -0,0 +1,459 @@
+%wing shear flow
+clear all;
+close all;
+
+Vx = 1; Vz = 1; My = 1; %test loads will be applied individually
+
+
+%Ixz = -Ixz;
+
+%define webs
+
+%% web cell 1
+
+%upper webs
+numStringers = numTopStringers;
+stringerGap = upperStringerGap;
+webThickness = t_upper;
+tempStringers = topStringers;
+
+for i=1:(numStringers+1)
+ web(i).xStart = sparCaps(1).posX + stringerGap*(i-1);
+ web(i).xEnd = sparCaps(1).posX + stringerGap*(i);
+ web(i).thickness = webThickness;
+ web(i).zStart = get_z(web(i).xStart/chord,1)*chord;
+ web(i).zEnd = get_z(web(i).xEnd/chord,1)*chord;
+ if i==1
+ web(i).dp_area = sparCaps(1).area;
+ web(i).dP_X = 0;
+ web(i).dP_Z = 0;
+ web(i).qPrime_X = 0;
+ web(i).qPrime_Z = 0;
+ else
+ web(i).dp_area = tempStringers(i-1).area;
+ dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
+ web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area); %just Vx
+ web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area); %just Vz
+ web(i).qPrime_X = web(i-1).qPrime_X - web(i).dP_X;
+ web(i).qPrime_Z = web(i-1).qPrime_Z - web(i).dP_Z;
+ end
+ tempInt = get_int(web(i).xStart/chord,web(i).xEnd/chord,1)*chord^2; %integral of airfoil function
+ triangle1 = abs( (web(i).xStart - sparCaps(1).posX)*web(i).zStart/2);
+ triangle2 = abs((web(i).xEnd - sparCaps(1).posX)*web(i).zEnd/2);
+ web(i).Area = tempInt + triangle1 - triangle2;
+ web(i).ds = get_ds(web(i).xStart/chord,web(i).xEnd/chord,1)*chord;
+ web(i).dS_over_t = web(i).ds / web(i).thickness;
+
+ web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
+ web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
+ web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
+ web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
+ web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
+ web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
+ web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
+ web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
+end
+webTop = web;
+web = [];
+
+%rear spar
+i=1;
+web(i).xStart = sparCaps(3).posX;
+web(i).xEnd = sparCaps(4).posX;
+web(i).thickness = t_rearSpar;
+web(i).zStart = sparCaps(3).posZ;
+web(i).zEnd = sparCaps(4).posZ;
+web(i).dp_area = sparCaps(3).area;
+dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
+web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
+web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
+web(i).qPrime_X = webTop(numTopStringers+1).qPrime_X - web(i).dP_X;
+web(i).qPrime_Z = webTop(numTopStringers+1).qPrime_Z - web(i).dP_Z;
+
+web(i).Area = (sparCaps(3).posX-sparCaps(1).posX)*sparCaps(3).posZ/2 + ...
+ abs((sparCaps(3).posX-sparCaps(1).posX)*sparCaps(4).posZ/2);
+web(i).ds = abs(sparCaps(3).posZ - sparCaps(4).posZ);
+web(i).dS_over_t = web(i).ds / web(i).thickness;
+
+web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
+web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
+web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
+web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
+web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
+web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
+web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
+web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
+
+webRearSpar = web;
+web = [];
+
+
+%lower webs
+numStringers = numBottomStringers;
+stringerGap = lowerStringerGap;
+webThickness = t_lower;
+tempStringers = bottomStringers;
+
+for i=1:(numStringers+1)
+ web(i).xStart = sparCaps(4).posX - stringerGap*(i-1);
+ web(i).xEnd = sparCaps(4).posX - stringerGap*(i);
+ web(i).thickness = webThickness;
+ web(i).zStart = get_z(web(i).xStart/chord,0)*chord;
+ web(i).zEnd = get_z(web(i).xEnd/chord,0)*chord;
+ dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
+ if i==1
+ web(i).dp_area = sparCaps(4).area;
+ web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
+ web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
+ web(i).qPrime_X = webRearSpar.qPrime_X - web(i).dP_X;
+ web(i).qPrime_Z = webRearSpar.qPrime_Z - web(i).dP_Z;
+ else
+ web(i).dp_area = tempStringers(i-1).area;
+ web(i).dP_X = get_dp(dx,dz, Vx,0,Ix,Iz,Ixz,web(i).dp_area);
+ web(i).dP_Z = get_dp(dx,dz, 0,Vz,Ix,Iz,Ixz,web(i).dp_area);
+ web(i).qPrime_X = web(i-1).qPrime_X - web(i).dP_X;
+ web(i).qPrime_Z = web(i-1).qPrime_Z - web(i).dP_Z;
+ end
+
+ tempInt = get_int(web(i).xEnd/chord,web(i).xStart/chord,0)*chord^2; %integral of airfoil function
+ triangle2 = abs((web(i).xStart - sparCaps(1).posX)*web(i).zStart/2);
+ triangle1 = abs((web(i).xEnd - sparCaps(1).posX)*web(i).zEnd/2);
+ web(i).Area = tempInt + triangle1 - triangle2;
+ web(i).ds = get_ds(web(i).xStart/chord,web(i).xEnd/chord,0)*chord;
+ web(i).dS_over_t = web(i).ds / web(i).thickness;
+
+ web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
+ web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
+ web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
+ web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
+ web(i).qp_dx_X = web(i).qPrime_X*(web(i).xEnd-web(i).xStart);
+ web(i).qp_dx_Z = web(i).qPrime_Z*(web(i).xEnd-web(i).xStart);
+ web(i).qp_dz_X = web(i).qPrime_X*(web(i).zEnd-web(i).zStart);
+ web(i).qp_dz_Z = web(i).qPrime_Z*(web(i).zEnd-web(i).zStart);
+
+ %web(i).radCurv = ... Example: get_curve(web(i).xStart,web(i).xEnd,1)
+end
+webBottom = web;
+web = [];
+
+%front Spar
+i=1;
+web(i).xStart = sparCaps(2).posX;
+web(i).xEnd = sparCaps(1).posX;
+web(i).thickness = t_frontSpar;
+web(i).zStart = sparCaps(2).posZ;
+web(i).zEnd = sparCaps(1).posZ;
+web(i).dp_area = sparCaps(2).area;
+dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
+web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
+web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
+web(i).qPrime_X = webBottom(numBottomStringers+1).qPrime_X - web(i).dP_X;
+web(i).qPrime_Z = webBottom(numBottomStringers+1).qPrime_Z - web(i).dP_Z;
+web(i).Area = 0;
+web(i).ds = abs(sparCaps(2).posZ - sparCaps(1).posZ);
+web(i).dS_over_t = web(i).ds / web(i).thickness;
+
+web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
+web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
+web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
+web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
+web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
+web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
+web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
+web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
+
+webFrontSpar = web;
+web = [];
+
+
+
+
+%% web cell 2
+
+%lower nose webs
+numStringers = numNoseBottomStringers;
+stringerGap = lowerNoseStringerGap;
+webThickness = t_lower_front;
+tempStringers = noseBottomStringers;
+
+for i=1:(numStringers+1)
+ web(i).xStart = sparCaps(2).posX - stringerGap*(i-1);
+ web(i).xEnd = sparCaps(2).posX - stringerGap*(i);
+ web(i).thickness = webThickness;
+ web(i).zStart = get_z(web(i).xStart/chord,0)*chord;
+ web(i).zEnd = get_z(web(i).xEnd/chord,0)*chord;
+ dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
+
+ if i==1
+ web(i).dp_area = sparCaps(2).area;
+ web(i).dP_X = 0;
+ web(i).dP_Z = 0;
+ web(i).qPrime_X = 0;
+ web(i).qPrime_Z = 0;
+ else
+ web(i).dp_area = tempStringers(i-1).area;
+ web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
+ web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
+ web(i).qPrime_X = web(i-1).qPrime_X - web(i).dP_X;
+ web(i).qPrime_Z = web(i-1).qPrime_Z - web(i).dP_Z;
+ end
+ tempInt = get_int(web(i).xEnd/chord,web(i).xStart/chord,0)*chord^2; %integral of airfoil function
+ triangle1 = abs((web(i).xStart - sparCaps(2).posX)*web(i).zStart/2);
+ triangle2 = abs((web(i).xEnd - sparCaps(2).posX)*web(i).zEnd/2);
+ web(i).Area = tempInt + triangle1 - triangle2;
+ web(i).ds = get_ds(web(i).xStart/chord,web(i).xEnd/chord,0)*chord;
+ web(i).dS_over_t = web(i).ds / web(i).thickness;
+
+ web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
+ web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
+ web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
+ web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
+ web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
+ web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
+ web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
+ web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
+
+ %web(i).radCurv = ... Example: get_curve(web(i).xStart,web(i).xEnd,1)
+end
+webLowerNose = web;
+web = [];
+
+%upper nose webs
+numStringers = numNoseTopStringers;
+stringerGap = upperNoseStringerGap;
+webThickness = t_upper_front;
+tempStringers = noseTopStringers;
+
+for i=1:(numStringers+1)
+ web(i).xStart = stringerGap*(i-1);
+ web(i).xEnd = stringerGap*(i);
+ web(i).thickness = webThickness;
+ web(i).zStart = get_z(web(i).xStart/chord,1)*chord;
+ web(i).zEnd = get_z(web(i).xEnd/chord,1)*chord;
+ dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
+ if i==1
+ web(i).dp_area = 0;
+ web(i).dP_X = 0;
+ web(i).dP_Z = 0;
+ web(i).qPrime_X = webLowerNose(numNoseBottomStringers+1).qPrime_X - web(i).dP_X;
+ web(i).qPrime_Z = webLowerNose(numNoseBottomStringers+1).qPrime_Z - web(i).dP_Z;
+ else
+ web(i).dp_area = tempStringers(i-1).area;
+ web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
+ web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
+ web(i).qPrime_X = web(i-1).qPrime_X - web(i).dP_X;
+ web(i).qPrime_Z = web(i-1).qPrime_Z - web(i).dP_Z;
+ end
+ tempInt = get_int(web(i).xStart/chord,web(i).xEnd/chord,1)*chord^2; %integral of airfoil function
+ triangle2 = abs((web(i).xStart - sparCaps(2).posX)*web(i).zStart/2);
+ triangle1 = abs((web(i).xEnd - sparCaps(2).posX)*web(i).zEnd/2);
+ web(i).Area = tempInt + triangle1 - triangle2;
+ web(i).ds = get_ds(web(i).xStart/chord,web(i).xEnd/chord,1)*chord;
+ web(i).dS_over_t = web(i).ds / web(i).thickness;
+
+ web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
+ web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
+ web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
+ web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
+ web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
+ web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
+ web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
+ web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
+
+end
+webUpperNose = web;
+web = [];
+
+
+%front Spar
+i=1;
+web(i).xStart = sparCaps(1).posX;
+web(i).xEnd = sparCaps(2).posX;
+web(i).thickness = t_frontSpar;
+web(i).zStart = sparCaps(1).posZ;
+web(i).zEnd = sparCaps(2).posZ;
+web(i).dp_area = sparCaps(1).area;
+dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
+
+web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
+web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
+web(i).qPrime_X = webUpperNose(numNoseTopStringers+1).qPrime_X - web(i).dP_X;
+web(i).qPrime_Z = webUpperNose(numNoseTopStringers+1).qPrime_Z - web(i).dP_Z;
+web(i).Area = 0;
+web(i).ds = abs(sparCaps(1).posZ - sparCaps(2).posZ);
+web(i).dS_over_t = web(i).ds / web(i).thickness;
+web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
+web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
+web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
+web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
+web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
+web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
+web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
+web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
+
+webFrontSparCell2 = web;
+web = [];
+
+
+%check that q'*dx sums up to Vx
+
+Fx = sum([webTop.qp_dx_X])+webRearSpar.qp_dx_X+ sum([webBottom.qp_dx_X])+webFrontSpar.qp_dx_X; %cell 1
+Fx = Fx + sum([webLowerNose.qp_dx_X])+ sum([webUpperNose.qp_dx_X]); %cell 2
+Fx
+Fz = sum([webTop.qp_dz_X])+webRearSpar.qp_dz_X+ sum([webBottom.qp_dz_X])+webFrontSpar.qp_dz_X; %cell 1
+Fz = Fz + sum([webLowerNose.qp_dz_X])+ sum([webUpperNose.qp_dz_X]); %cell 2
+Fz
+
+%check that q'*dz sums up to Vz
+
+
+Fx = sum([webTop.qp_dx_Z])+webRearSpar.qp_dx_Z+ sum([webBottom.qp_dx_Z])+webFrontSpar.qp_dx_Z; %cell 1
+Fx = Fx + sum([webLowerNose.qp_dx_Z])+ sum([webUpperNose.qp_dx_Z]); %cell 2
+Fx
+Fz = sum([webTop.qp_dz_Z])+webRearSpar.qp_dz_Z+ sum([webBottom.qp_dz_Z])+webFrontSpar.qp_dz_Z; %cell 1
+Fz = Fz + sum([webLowerNose.qp_dz_Z])+ sum([webUpperNose.qp_dz_Z]); %cell 2
+Fz
+
+%%
+
+% sum up the ds/t and q*ds/t to solve 2 equations, 2 unknowns
+
+% [A]*[q1s q2s] = B
+
+A11 = sum([webTop.dS_over_t])+webRearSpar.dS_over_t+ sum([webBottom.dS_over_t])+webFrontSpar.dS_over_t;
+A22 = sum([webLowerNose.dS_over_t])+ sum([webUpperNose.dS_over_t])+webFrontSparCell2.dS_over_t;
+A12 = -webFrontSpar.dS_over_t;
+A21 = -webFrontSparCell2.dS_over_t;
+
+B1_X = sum([webTop.q_dS_over_t_X])+webRearSpar.q_dS_over_t_X+ sum([webBottom.q_dS_over_t_X])+webFrontSpar.q_dS_over_t_X;
+B2_X = sum([webLowerNose.q_dS_over_t_X])+ sum([webUpperNose.q_dS_over_t_X])+webFrontSparCell2.q_dS_over_t_X;
+B1_Z = sum([webTop.q_dS_over_t_Z])+webRearSpar.q_dS_over_t_Z+ sum([webBottom.q_dS_over_t_Z])+webFrontSpar.q_dS_over_t_Z;
+B2_Z = sum([webLowerNose.q_dS_over_t_Z])+ sum([webUpperNose.q_dS_over_t_Z])+webFrontSparCell2.q_dS_over_t_Z;
+
+Amat = [A11 A12; A21 A22];
+Bmat_X = -[B1_X;B2_X];
+Bmat_Z = -[B1_Z;B2_Z];
+
+qs_X = inv(Amat)*Bmat_X;
+qs_Z = inv(Amat)*Bmat_Z;
+
+
+
+sum_2_a_q_X = sum([webTop.two_A_qprime_X])+webRearSpar.two_A_qprime_X+ sum([webBottom.two_A_qprime_X]); %cell 1 qprimes
+sum_2_a_q_X = sum_2_a_q_X + sum([webLowerNose.two_A_qprime_X])+ sum([webUpperNose.two_A_qprime_X]); %cell 2 qprimes
+sum_2_a_q_X = sum_2_a_q_X + 2*qs_X(1)*(sum([webTop.Area])+webRearSpar.Area+ sum([webBottom.Area]));
+sum_2_a_q_X = sum_2_a_q_X + 2*qs_X(2)*(sum([webLowerNose.Area])+ sum([webUpperNose.Area]));
+
+sum_2_a_q_Z = sum([webTop.two_A_qprime_Z])+webRearSpar.two_A_qprime_Z+ sum([webBottom.two_A_qprime_Z]); %cell 1 qprimes
+sum_2_a_q_Z = sum_2_a_q_Z + sum([webLowerNose.two_A_qprime_Z])+ sum([webUpperNose.two_A_qprime_Z]); %cell 2 qprimes
+sum_2_a_q_Z = sum_2_a_q_Z + 2*qs_Z(1)*(sum([webTop.Area])+webRearSpar.Area+ sum([webBottom.Area]));
+sum_2_a_q_Z = sum_2_a_q_Z + 2*qs_Z(2)*(sum([webLowerNose.Area])+ sum([webUpperNose.Area]));
+
+%shear center
+sc.posX = sum_2_a_q_Z / Vz + frontSpar*chord;
+sc.posZ = - sum_2_a_q_X / Vx;
+
+
+% now consider the torque representing shifting the load from the quarter
+% chord to the SC (need to check signs on these moments)
+
+torque_Z = Vz*(sc.posX - 0.25*chord);
+torque_X = -Vx*sc.posZ;
+
+
+Area1 = sum([webTop.Area]) + webRearSpar.Area + sum([webBottom.Area]);
+%check area
+Area1_check = get_int(frontSpar,backSpar,1)*chord^2 + get_int(frontSpar,backSpar,0)*chord^2;
+
+Area2 = sum([webLowerNose.Area]) + sum([webUpperNose.Area]);
+Area2_check = get_int(0,frontSpar,1)*chord^2 + get_int(0,frontSpar,0)*chord^2;
+
+
+%for twist equation (see excel spreadsheet example)
+
+q1t_over_q2t = (A22/Area2 + webFrontSpar.dS_over_t/Area1)/(A11/Area1 + webFrontSpar.dS_over_t/Area2);
+
+q2t = torque_X/(2*Area1*q1t_over_q2t + 2*Area2);
+q1t = q2t*q1t_over_q2t;
+qt_X = [q1t;q2t];
+
+q2t = torque_Z/(2*Area1*q1t_over_q2t + 2*Area2);
+q1t = q2t*q1t_over_q2t;
+qt_Z = [q1t;q2t];
+
+
+
+% --- - add up all shear flows: qtot = (qPrime + qs) + qt
+
+
+
+
+%--- insert force balance to check total shear flows ---
+
+% --- --
+
+
+%end
+
+sc
+
+
+%plotting airfoil cross-section
+
+xChord = 0:.01:1;
+xChord = xChord*chord;
+upperSurface = zeros(1,length(xChord));
+lowerSurface = zeros(1,length(xChord));
+
+for i=1:length(xChord)
+ upperSurface(i) = get_z(xChord(i)/chord,1)*chord;
+ lowerSurface(i) = get_z(xChord(i)/chord,0)*chord;
+end
+
+figure; hold on; axis equal; grid on;
+%plot(xChord,z_camber,'-')
+plot(xChord,upperSurface,'-k','linewidth',2)
+plot(xChord,lowerSurface,'-k','linewidth',2)
+plot([0 1],[0 0],'--k','linewidth',1)
+
+
+for i = 1:length(webTop)
+ vecX = [frontSpar*chord webTop(i).xStart webTop(i).xEnd];
+ vecZ = [0 webTop(i).zStart webTop(i).zEnd];
+ fill(vecX,vecZ,[0.9 0.9 0.9])
+end
+
+for i = 1:length(webBottom)
+ vecX = [frontSpar*chord webBottom(i).xStart webBottom(i).xEnd];
+ vecZ = [0 webBottom(i).zStart webBottom(i).zEnd];
+ fill(vecX,vecZ,[0.9 0.9 0.9])
+end
+
+for i = 1:length(webUpperNose)
+ vecX = [frontSpar*chord webUpperNose(i).xStart webUpperNose(i).xEnd];
+ vecZ = [0 webUpperNose(i).zStart webUpperNose(i).zEnd];
+ fill(vecX,vecZ,[0.7 0.9 1.0])
+end
+
+for i = 1:length(webLowerNose)
+ vecX = [frontSpar*chord webLowerNose(i).xStart webLowerNose(i).xEnd];
+ vecZ = [0 webLowerNose(i).zStart webLowerNose(i).zEnd];
+ fill(vecX,vecZ,[0.7 0.9 1.0])
+end
+
+ vecX = [frontSpar*chord sparCaps(3).posX sparCaps(4).posX];
+ vecZ = [0 sparCaps(3).posZ sparCaps(4).posZ];
+ fill(vecX,vecZ,[0.9 0.9 0.9])
+
+
+sparCapSize = 18;
+stringerSize = 18;
+plot([sparCaps(1).posX sparCaps(2).posX],[sparCaps(1).posZ sparCaps(2).posZ],'-k','linewidth',2)
+plot([sparCaps(3).posX sparCaps(4).posX],[sparCaps(3).posZ sparCaps(4).posZ],'-k','linewidth',2)
+plot([sparCaps.posX],[sparCaps.posZ],'.b','markersize',sparCapSize)
+plot([topStringers.posX],[topStringers.posZ],'.r','markersize',stringerSize)
+plot([bottomStringers.posX],[bottomStringers.posZ],'.r','markersize',stringerSize)
+plot([noseTopStringers.posX],[noseTopStringers.posZ],'.r','markersize',stringerSize)
+plot([noseBottomStringers.posX],[noseBottomStringers.posZ],'.r','markersize',stringerSize)
+plot(centroid.posX,centroid.posZ,'.k','markerSize',18)
+plot(sc.posX,sc.posZ,'.g','markersize',18)
diff --git a/wing_scripts/stringersBeamExample.m b/wing_scripts/stringersBeamExample.m
new file mode 100644
index 0000000..cd3bcb6
--- /dev/null
+++ b/wing_scripts/stringersBeamExample.m
@@ -0,0 +1,47 @@
+close all;
+force = 8000; % lbs
+stringer_A = 0.5; % in^2
+thickness = 0.04; % in
+
+top_stringers_y = 6; % in
+middle_stringers_y = 2; % in
+
+I = 2*stringer_A*top_stringers_y^2 + 2*stringer_A*middle_stringers_y^2;
+
+% solve for shear stress distribution. this calc ignores the thickness of
+% the web between teh stringers (assumes bending taken by stringers)
+% V / (I * t) * int(y*da)
+
+shear_top_web = force / (I*thickness) * top_stringers_y * stringer_A;
+shear_middle_web = shear_top_web + (force / (I*thickness) * middle_stringers_y * stringer_A);
+
+figure; grid on; hold on;set(gcf,'color',[1 1 1]);
+
+
+plot([shear_top_web shear_top_web],[middle_stringers_y top_stringers_y],'linewidth',2);
+plot([shear_middle_web shear_middle_web],[-middle_stringers_y middle_stringers_y],'linewidth',2);
+plot([shear_top_web shear_top_web],[-middle_stringers_y -top_stringers_y],'linewidth',2);
+
+plot([0 shear_top_web],[top_stringers_y top_stringers_y],'linewidth',2);
+plot([0 shear_top_web],[-top_stringers_y -top_stringers_y],'linewidth',2);
+plot([shear_middle_web shear_top_web],[middle_stringers_y middle_stringers_y],'linewidth',2);
+plot([shear_middle_web shear_top_web],[-middle_stringers_y -middle_stringers_y],'linewidth',2);
+xlabel('shear stress (lb/in^2)','fontsize',16,'fontweight','bold');ylabel('Distance from Center (in)','fontsize',16,'fontweight','bold')
+set(gca,'FontSize',16,'fontweight','bold');
+
+%Alternate approach.. compute change in bending stress at each stringer to
+%find the change in shear load
+
+%at top stringer
+d_sigma = force * top_stringers_y / I; %(lbs/in^2)
+d_force_top = d_sigma * stringer_A;
+
+%at middle stringer..
+d_sigma = force * middle_stringers_y / I; %(lbs/in^2)
+d_force_middle = d_force_top + d_sigma*stringer_A;
+
+%check if load balances
+check_load = 2*d_force_top*4 + d_force_middle*4;
+
+
+
diff --git a/wing_scripts/wingAnalysis_190422.m b/wing_scripts/wingAnalysis_190422.m
new file mode 100644
index 0000000..a7d65e2
--- /dev/null
+++ b/wing_scripts/wingAnalysis_190422.m
@@ -0,0 +1,579 @@
+%wing shear flow
+clear all;
+close all;
+
+
+
+
+Vx = 1; Vz = 1; My = 1; %test loads will be applied individually
+
+%define a few
+numTopStringers = 6;
+numBottomStringers = 8;
+numNoseTopStringers = 4;
+numNoseBottomStringers = 4;
+
+t_upper = 0.02/12;
+t_lower = 0.02/12;
+t_upper_front = 0.02/12;
+t_lower_front = 0.02/12;
+t_frontSpar = 0.04/12;
+t_rearSpar = 0.04/12;
+
+frontSpar = 0.2;
+backSpar = 0.7;
+chord = 5;
+
+sparCaps(1).posX = frontSpar*chord;
+sparCaps(2).posX = frontSpar*chord;
+sparCaps(3).posX = backSpar*chord;
+sparCaps(4).posX = backSpar*chord;
+
+sparCaps(1).posZ = get_z(frontSpar,1)*chord;
+sparCaps(2).posZ = get_z(frontSpar,0)*chord;
+sparCaps(3).posZ = get_z(backSpar,1)*chord;
+sparCaps(4).posZ = get_z(backSpar,0)*chord;
+
+sparCaps(1).area = .1;
+sparCaps(2).area = .1;
+sparCaps(3).area = .1;
+sparCaps(4).area = .1;
+
+upperStringerGap = (sparCaps(3).posX - sparCaps(1).posX)/(numTopStringers + 1);
+lowerStringerGap = (sparCaps(3).posX - sparCaps(1).posX)/(numBottomStringers + 1);
+upperNoseStringerGap = (sparCaps(1).posX - 0)/(numNoseTopStringers + 1);
+lowerNoseStringerGap = (sparCaps(1).posX - 0)/(numNoseBottomStringers + 1);
+
+
+%set stringers spaced evenly along X axis betwen Spars
+%top Stringers
+for i=1:numTopStringers
+ topStringers(i).posX = sparCaps(1).posX + upperStringerGap*i;
+ topStringers(i).posZ = get_z(topStringers(i).posX/chord,1)*chord;
+ topStringers(i).area = .1;
+end
+
+%bottom Stringers
+for i=1:numBottomStringers
+ bottomStringers(i).posX = sparCaps(4).posX - lowerStringerGap*i;
+ bottomStringers(i).posZ = get_z(bottomStringers(i).posX/chord,0)*chord;
+ bottomStringers(i).area = .1;
+
+end
+
+%nose bottom Stringers
+for i=1:numNoseBottomStringers
+ noseBottomStringers(i).posX = sparCaps(2).posX - lowerNoseStringerGap*i;
+ noseBottomStringers(i).posZ = get_z(noseBottomStringers(i).posX/chord,0)*chord;
+ noseBottomStringers(i).area = .1;
+end
+
+%nose top Stringers
+for i=1:numNoseTopStringers
+ noseTopStringers(i).posX = upperNoseStringerGap*i;
+ noseTopStringers(i).posZ = get_z(noseTopStringers(i).posX/chord,1)*chord;
+ noseTopStringers(i).area = .1;
+end
+
+
+centroid.posX = sum([sparCaps.posX].*[sparCaps.area]) + ...
+ sum([topStringers.posX].*[topStringers.area]) + ...
+ sum([bottomStringers.posX].*[bottomStringers.area]) + ...
+ sum([noseTopStringers.posX].*[noseTopStringers.area]) + ...
+ sum([noseBottomStringers.posX].*[noseBottomStringers.area]);
+
+centroid.posX = centroid.posX / ( sum([sparCaps.area]) + sum([topStringers.area]) + ...
+ sum([bottomStringers.area]) + sum([noseTopStringers.area]) + sum([noseBottomStringers.area]));
+
+centroid.posZ = sum([sparCaps.posZ].*[sparCaps.area]) + ...
+ sum([topStringers.posZ].*[topStringers.area]) + ...
+ sum([bottomStringers.posZ].*[bottomStringers.area]) + ...
+ sum([noseTopStringers.posZ].*[noseTopStringers.area]) + ...
+ sum([noseBottomStringers.posZ].*[noseBottomStringers.area]);
+
+centroid.posZ = centroid.posZ / ( sum([sparCaps.area]) + sum([topStringers.area]) + ...
+ sum([bottomStringers.area]) + sum([noseTopStringers.area]) + sum([noseBottomStringers.area]));
+
+%summing contributions for inertia terms
+Ix = 0; Iz = 0; Ixz = 0;
+
+for i=1:4 %spar caps
+ Ix = Ix + sparCaps(i).area*(sparCaps(i).posZ-centroid.posZ)^2;
+ Iz = Iz + sparCaps(i).area*(sparCaps(i).posX-centroid.posX)^2;
+ Ixz = Ixz + sparCaps(i).area*(sparCaps(i).posX-centroid.posX)*(sparCaps(i).posZ-centroid.posZ);
+end
+
+
+for i=1:numTopStringers %top stringers
+ Ix = Ix + topStringers(i).area*(topStringers(i).posZ-centroid.posZ)^2;
+ Iz = Iz + topStringers(i).area*(topStringers(i).posX-centroid.posX)^2;
+ Ixz = Ixz + topStringers(i).area*(topStringers(i).posX-centroid.posX)*(topStringers(i).posZ-centroid.posZ);
+end
+for i=1:numBottomStringers %bottom stringers
+ Ix = Ix + bottomStringers(i).area*(bottomStringers(i).posZ-centroid.posZ)^2;
+ Iz = Iz + bottomStringers(i).area*(bottomStringers(i).posX-centroid.posX)^2;
+ Ixz = Ixz + bottomStringers(i).area*(bottomStringers(i).posX-centroid.posX)*(bottomStringers(i).posZ-centroid.posZ);
+end
+for i=1:numNoseTopStringers %nose top stringers
+ Ix = Ix + noseTopStringers(i).area*(noseTopStringers(i).posZ-centroid.posZ)^2;
+ Iz = Iz + noseTopStringers(i).area*(noseTopStringers(i).posX-centroid.posX)^2;
+ Ixz = Ixz + noseTopStringers(i).area*(noseTopStringers(i).posX-centroid.posX)*(noseTopStringers(i).posZ-centroid.posZ);
+end
+for i=1:numNoseBottomStringers %nose bottom stringers
+ Ix = Ix + noseBottomStringers(i).area*(noseBottomStringers(i).posZ-centroid.posZ)^2;
+ Iz = Iz + noseBottomStringers(i).area*(noseBottomStringers(i).posX-centroid.posX)^2;
+ Ixz = Ixz + noseBottomStringers(i).area*(noseBottomStringers(i).posX-centroid.posX)*(noseBottomStringers(i).posZ-centroid.posZ);
+end
+
+%Ixz = -Ixz;
+
+%define webs
+
+%% web cell 1
+
+%upper webs
+numStringers = numTopStringers;
+stringerGap = upperStringerGap;
+webThickness = t_upper;
+tempStringers = topStringers;
+
+for i=1:(numStringers+1)
+ web(i).xStart = sparCaps(1).posX + stringerGap*(i-1);
+ web(i).xEnd = sparCaps(1).posX + stringerGap*(i);
+ web(i).thickness = webThickness;
+ web(i).zStart = get_z(web(i).xStart/chord,1)*chord;
+ web(i).zEnd = get_z(web(i).xEnd/chord,1)*chord;
+ if i==1
+ web(i).dp_area = sparCaps(1).area;
+ web(i).dP_X = 0;
+ web(i).dP_Z = 0;
+ web(i).qPrime_X = 0;
+ web(i).qPrime_Z = 0;
+ else
+ web(i).dp_area = tempStringers(i-1).area;
+ dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
+ web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area); %just Vx
+ web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area); %just Vz
+ web(i).qPrime_X = web(i-1).qPrime_X - web(i).dP_X;
+ web(i).qPrime_Z = web(i-1).qPrime_Z - web(i).dP_Z;
+ end
+ tempInt = get_int(web(i).xStart/chord,web(i).xEnd/chord,1)*chord^2; %integral of airfoil function
+ triangle1 = abs( (web(i).xStart - sparCaps(1).posX)*web(i).zStart/2);
+ triangle2 = abs((web(i).xEnd - sparCaps(1).posX)*web(i).zEnd/2);
+ web(i).Area = tempInt + triangle1 - triangle2;
+ web(i).ds = get_ds(web(i).xStart/chord,web(i).xEnd/chord,1)*chord;
+ web(i).dS_over_t = web(i).ds / web(i).thickness;
+
+ web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
+ web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
+ web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
+ web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
+ web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
+ web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
+ web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
+ web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
+end
+webTop = web;
+web = [];
+
+%rear spar
+i=1;
+web(i).xStart = sparCaps(3).posX;
+web(i).xEnd = sparCaps(4).posX;
+web(i).thickness = t_rearSpar;
+web(i).zStart = sparCaps(3).posZ;
+web(i).zEnd = sparCaps(4).posZ;
+web(i).dp_area = sparCaps(3).area;
+dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
+web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
+web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
+web(i).qPrime_X = webTop(numTopStringers+1).qPrime_X - web(i).dP_X;
+web(i).qPrime_Z = webTop(numTopStringers+1).qPrime_Z - web(i).dP_Z;
+
+web(i).Area = (sparCaps(3).posX-sparCaps(1).posX)*sparCaps(3).posZ/2 + ...
+ abs((sparCaps(3).posX-sparCaps(1).posX)*sparCaps(4).posZ/2);
+web(i).ds = abs(sparCaps(3).posZ - sparCaps(4).posZ);
+web(i).dS_over_t = web(i).ds / web(i).thickness;
+
+web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
+web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
+web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
+web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
+web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
+web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
+web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
+web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
+
+webRearSpar = web;
+web = [];
+
+
+%lower webs
+numStringers = numBottomStringers;
+stringerGap = lowerStringerGap;
+webThickness = t_lower;
+tempStringers = bottomStringers;
+
+for i=1:(numStringers+1)
+ web(i).xStart = sparCaps(4).posX - stringerGap*(i-1);
+ web(i).xEnd = sparCaps(4).posX - stringerGap*(i);
+ web(i).thickness = webThickness;
+ web(i).zStart = get_z(web(i).xStart/chord,0)*chord;
+ web(i).zEnd = get_z(web(i).xEnd/chord,0)*chord;
+ dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
+ if i==1
+ web(i).dp_area = sparCaps(4).area;
+ web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
+ web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
+ web(i).qPrime_X = webRearSpar.qPrime_X - web(i).dP_X;
+ web(i).qPrime_Z = webRearSpar.qPrime_Z - web(i).dP_Z;
+ else
+ web(i).dp_area = tempStringers(i-1).area;
+ web(i).dP_X = get_dp(dx,dz, Vx,0,Ix,Iz,Ixz,web(i).dp_area);
+ web(i).dP_Z = get_dp(dx,dz, 0,Vz,Ix,Iz,Ixz,web(i).dp_area);
+ web(i).qPrime_X = web(i-1).qPrime_X - web(i).dP_X;
+ web(i).qPrime_Z = web(i-1).qPrime_Z - web(i).dP_Z;
+ end
+
+ tempInt = get_int(web(i).xEnd/chord,web(i).xStart/chord,0)*chord^2; %integral of airfoil function
+ triangle2 = abs((web(i).xStart - sparCaps(1).posX)*web(i).zStart/2);
+ triangle1 = abs((web(i).xEnd - sparCaps(1).posX)*web(i).zEnd/2);
+ web(i).Area = tempInt + triangle1 - triangle2;
+ web(i).ds = get_ds(web(i).xStart/chord,web(i).xEnd/chord,0)*chord;
+ web(i).dS_over_t = web(i).ds / web(i).thickness;
+
+ web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
+ web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
+ web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
+ web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
+ web(i).qp_dx_X = web(i).qPrime_X*(web(i).xEnd-web(i).xStart);
+ web(i).qp_dx_Z = web(i).qPrime_Z*(web(i).xEnd-web(i).xStart);
+ web(i).qp_dz_X = web(i).qPrime_X*(web(i).zEnd-web(i).zStart);
+ web(i).qp_dz_Z = web(i).qPrime_Z*(web(i).zEnd-web(i).zStart);
+
+ %web(i).radCurv = ... Example: get_curve(web(i).xStart,web(i).xEnd,1)
+end
+webBottom = web;
+web = [];
+
+%front Spar
+i=1;
+web(i).xStart = sparCaps(2).posX;
+web(i).xEnd = sparCaps(1).posX;
+web(i).thickness = t_frontSpar;
+web(i).zStart = sparCaps(2).posZ;
+web(i).zEnd = sparCaps(1).posZ;
+web(i).dp_area = sparCaps(2).area;
+dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
+web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
+web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
+web(i).qPrime_X = webBottom(numBottomStringers+1).qPrime_X - web(i).dP_X;
+web(i).qPrime_Z = webBottom(numBottomStringers+1).qPrime_Z - web(i).dP_Z;
+web(i).Area = 0;
+web(i).ds = abs(sparCaps(2).posZ - sparCaps(1).posZ);
+web(i).dS_over_t = web(i).ds / web(i).thickness;
+
+web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
+web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
+web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
+web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
+web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
+web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
+web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
+web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
+
+webFrontSpar = web;
+web = [];
+
+
+
+
+%% web cell 2
+
+%lower nose webs
+numStringers = numNoseBottomStringers;
+stringerGap = lowerNoseStringerGap;
+webThickness = t_lower_front;
+tempStringers = noseBottomStringers;
+
+for i=1:(numStringers+1)
+ web(i).xStart = sparCaps(2).posX - stringerGap*(i-1);
+ web(i).xEnd = sparCaps(2).posX - stringerGap*(i);
+ web(i).thickness = webThickness;
+ web(i).zStart = get_z(web(i).xStart/chord,0)*chord;
+ web(i).zEnd = get_z(web(i).xEnd/chord,0)*chord;
+ dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
+
+ if i==1
+ web(i).dp_area = sparCaps(2).area;
+ web(i).dP_X = 0;
+ web(i).dP_Z = 0;
+ web(i).qPrime_X = 0;
+ web(i).qPrime_Z = 0;
+ else
+ web(i).dp_area = tempStringers(i-1).area;
+ web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
+ web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
+ web(i).qPrime_X = web(i-1).qPrime_X - web(i).dP_X;
+ web(i).qPrime_Z = web(i-1).qPrime_Z - web(i).dP_Z;
+ end
+ tempInt = get_int(web(i).xEnd/chord,web(i).xStart/chord,0)*chord^2; %integral of airfoil function
+ triangle1 = abs((web(i).xStart - sparCaps(2).posX)*web(i).zStart/2);
+ triangle2 = abs((web(i).xEnd - sparCaps(2).posX)*web(i).zEnd/2);
+ web(i).Area = tempInt + triangle1 - triangle2;
+ web(i).ds = get_ds(web(i).xStart/chord,web(i).xEnd/chord,0)*chord;
+ web(i).dS_over_t = web(i).ds / web(i).thickness;
+
+ web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
+ web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
+ web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
+ web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
+ web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
+ web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
+ web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
+ web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
+
+ %web(i).radCurv = ... Example: get_curve(web(i).xStart,web(i).xEnd,1)
+end
+webLowerNose = web;
+web = [];
+
+%upper nose webs
+numStringers = numNoseTopStringers;
+stringerGap = upperNoseStringerGap;
+webThickness = t_upper_front;
+tempStringers = noseTopStringers;
+
+for i=1:(numStringers+1)
+ web(i).xStart = stringerGap*(i-1);
+ web(i).xEnd = stringerGap*(i);
+ web(i).thickness = webThickness;
+ web(i).zStart = get_z(web(i).xStart/chord,1)*chord;
+ web(i).zEnd = get_z(web(i).xEnd/chord,1)*chord;
+ dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
+ if i==1
+ web(i).dp_area = 0;
+ web(i).dP_X = 0;
+ web(i).dP_Z = 0;
+ web(i).qPrime_X = webLowerNose(numNoseBottomStringers+1).qPrime_X - web(i).dP_X;
+ web(i).qPrime_Z = webLowerNose(numNoseBottomStringers+1).qPrime_Z - web(i).dP_Z;
+ else
+ web(i).dp_area = tempStringers(i-1).area;
+ web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
+ web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
+ web(i).qPrime_X = web(i-1).qPrime_X - web(i).dP_X;
+ web(i).qPrime_Z = web(i-1).qPrime_Z - web(i).dP_Z;
+ end
+ tempInt = get_int(web(i).xStart/chord,web(i).xEnd/chord,1)*chord^2; %integral of airfoil function
+ triangle2 = abs((web(i).xStart - sparCaps(2).posX)*web(i).zStart/2);
+ triangle1 = abs((web(i).xEnd - sparCaps(2).posX)*web(i).zEnd/2);
+ web(i).Area = tempInt + triangle1 - triangle2;
+ web(i).ds = get_ds(web(i).xStart/chord,web(i).xEnd/chord,1)*chord;
+ web(i).dS_over_t = web(i).ds / web(i).thickness;
+
+ web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
+ web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
+ web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
+ web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
+ web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
+ web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
+ web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
+ web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
+
+end
+webUpperNose = web;
+web = [];
+
+
+%front Spar
+i=1;
+web(i).xStart = sparCaps(1).posX;
+web(i).xEnd = sparCaps(2).posX;
+web(i).thickness = t_frontSpar;
+web(i).zStart = sparCaps(1).posZ;
+web(i).zEnd = sparCaps(2).posZ;
+web(i).dp_area = sparCaps(1).area;
+dx = web(i).xStart-centroid.posX; dz = web(i).zStart-centroid.posZ;
+
+web(i).dP_X = get_dp(dx,dz,Vx,0,Ix,Iz,Ixz,web(i).dp_area);
+web(i).dP_Z = get_dp(dx,dz,0,Vz,Ix,Iz,Ixz,web(i).dp_area);
+web(i).qPrime_X = webUpperNose(numNoseTopStringers+1).qPrime_X - web(i).dP_X;
+web(i).qPrime_Z = webUpperNose(numNoseTopStringers+1).qPrime_Z - web(i).dP_Z;
+web(i).Area = 0;
+web(i).ds = abs(sparCaps(1).posZ - sparCaps(2).posZ);
+web(i).dS_over_t = web(i).ds / web(i).thickness;
+web(i).q_dS_over_t_X = web(i).qPrime_X * web(i).dS_over_t;
+web(i).q_dS_over_t_Z = web(i).qPrime_Z * web(i).dS_over_t;
+web(i).two_A_qprime_X = 2*web(i).Area*web(i).qPrime_X;
+web(i).two_A_qprime_Z = 2*web(i).Area*web(i).qPrime_Z;
+web(i).qp_dx_X = web(i).qPrime_X *(web(i).xEnd-web(i).xStart);
+web(i).qp_dx_Z = web(i).qPrime_Z *(web(i).xEnd-web(i).xStart);
+web(i).qp_dz_X = web(i).qPrime_X *(web(i).zEnd-web(i).zStart);
+web(i).qp_dz_Z = web(i).qPrime_Z *(web(i).zEnd-web(i).zStart);
+
+webFrontSparCell2 = web;
+web = [];
+
+
+%check that q'*dx sums up to Vx
+
+Fx = sum([webTop.qp_dx_X])+webRearSpar.qp_dx_X+ sum([webBottom.qp_dx_X])+webFrontSpar.qp_dx_X; %cell 1
+Fx = Fx + sum([webLowerNose.qp_dx_X])+ sum([webUpperNose.qp_dx_X]); %cell 2
+Fx
+Fz = sum([webTop.qp_dz_X])+webRearSpar.qp_dz_X+ sum([webBottom.qp_dz_X])+webFrontSpar.qp_dz_X; %cell 1
+Fz = Fz + sum([webLowerNose.qp_dz_X])+ sum([webUpperNose.qp_dz_X]); %cell 2
+Fz
+
+%check that q'*dz sums up to Vz
+
+
+Fx = sum([webTop.qp_dx_Z])+webRearSpar.qp_dx_Z+ sum([webBottom.qp_dx_Z])+webFrontSpar.qp_dx_Z; %cell 1
+Fx = Fx + sum([webLowerNose.qp_dx_Z])+ sum([webUpperNose.qp_dx_Z]); %cell 2
+Fx
+Fz = sum([webTop.qp_dz_Z])+webRearSpar.qp_dz_Z+ sum([webBottom.qp_dz_Z])+webFrontSpar.qp_dz_Z; %cell 1
+Fz = Fz + sum([webLowerNose.qp_dz_Z])+ sum([webUpperNose.qp_dz_Z]); %cell 2
+Fz
+
+%%
+
+% sum up the ds/t and q*ds/t to solve 2 equations, 2 unknowns
+
+% [A]*[q1s q2s] = B
+
+A11 = sum([webTop.dS_over_t])+webRearSpar.dS_over_t+ sum([webBottom.dS_over_t])+webFrontSpar.dS_over_t;
+A22 = sum([webLowerNose.dS_over_t])+ sum([webUpperNose.dS_over_t])+webFrontSparCell2.dS_over_t;
+A12 = -webFrontSpar.dS_over_t;
+A21 = -webFrontSparCell2.dS_over_t;
+
+B1_X = sum([webTop.q_dS_over_t_X])+webRearSpar.q_dS_over_t_X+ sum([webBottom.q_dS_over_t_X])+webFrontSpar.q_dS_over_t_X;
+B2_X = sum([webLowerNose.q_dS_over_t_X])+ sum([webUpperNose.q_dS_over_t_X])+webFrontSparCell2.q_dS_over_t_X;
+B1_Z = sum([webTop.q_dS_over_t_Z])+webRearSpar.q_dS_over_t_Z+ sum([webBottom.q_dS_over_t_Z])+webFrontSpar.q_dS_over_t_Z;
+B2_Z = sum([webLowerNose.q_dS_over_t_Z])+ sum([webUpperNose.q_dS_over_t_Z])+webFrontSparCell2.q_dS_over_t_Z;
+
+Amat = [A11 A12; A21 A22];
+Bmat_X = -[B1_X;B2_X];
+Bmat_Z = -[B1_Z;B2_Z];
+
+qs_X = inv(Amat)*Bmat_X;
+qs_Z = inv(Amat)*Bmat_Z;
+
+
+
+sum_2_a_q_X = sum([webTop.two_A_qprime_X])+webRearSpar.two_A_qprime_X+ sum([webBottom.two_A_qprime_X]); %cell 1 qprimes
+sum_2_a_q_X = sum_2_a_q_X + sum([webLowerNose.two_A_qprime_X])+ sum([webUpperNose.two_A_qprime_X]); %cell 2 qprimes
+sum_2_a_q_X = sum_2_a_q_X + 2*qs_X(1)*(sum([webTop.Area])+webRearSpar.Area+ sum([webBottom.Area]));
+sum_2_a_q_X = sum_2_a_q_X + 2*qs_X(2)*(sum([webLowerNose.Area])+ sum([webUpperNose.Area]));
+
+sum_2_a_q_Z = sum([webTop.two_A_qprime_Z])+webRearSpar.two_A_qprime_Z+ sum([webBottom.two_A_qprime_Z]); %cell 1 qprimes
+sum_2_a_q_Z = sum_2_a_q_Z + sum([webLowerNose.two_A_qprime_Z])+ sum([webUpperNose.two_A_qprime_Z]); %cell 2 qprimes
+sum_2_a_q_Z = sum_2_a_q_Z + 2*qs_Z(1)*(sum([webTop.Area])+webRearSpar.Area+ sum([webBottom.Area]));
+sum_2_a_q_Z = sum_2_a_q_Z + 2*qs_Z(2)*(sum([webLowerNose.Area])+ sum([webUpperNose.Area]));
+
+%shear center
+sc.posX = sum_2_a_q_Z / Vz + frontSpar*chord;
+sc.posZ = - sum_2_a_q_X / Vx;
+
+
+% now consider the torque representing shifting the load from the quarter
+% chord to the SC (need to check signs on these moments)
+
+torque_Z = Vz*(sc.posX - 0.25*chord);
+torque_X = -Vx*sc.posZ;
+
+
+Area1 = sum([webTop.Area]) + webRearSpar.Area + sum([webBottom.Area]);
+%check area
+Area1_check = get_int(frontSpar,backSpar,1)*chord^2 + get_int(frontSpar,backSpar,0)*chord^2;
+
+Area2 = sum([webLowerNose.Area]) + sum([webUpperNose.Area]);
+Area2_check = get_int(0,frontSpar,1)*chord^2 + get_int(0,frontSpar,0)*chord^2;
+
+
+%for twist equation (see excel spreadsheet example)
+
+q1t_over_q2t = (A22/Area2 + webFrontSpar.dS_over_t/Area1)/(A11/Area1 + webFrontSpar.dS_over_t/Area2);
+
+q2t = torque_X/(2*Area1*q1t_over_q2t + 2*Area2);
+q1t = q2t*q1t_over_q2t;
+qt_X = [q1t;q2t];
+
+q2t = torque_Z/(2*Area1*q1t_over_q2t + 2*Area2);
+q1t = q2t*q1t_over_q2t;
+qt_Z = [q1t;q2t];
+
+
+
+% --- - add up all shear flows: qtot = (qPrime + qs) + qt
+
+
+
+
+%--- insert force balance to check total shear flows ---
+
+% --- --
+
+
+%end
+
+sc
+
+
+%plotting airfoil cross-section
+
+xChord = 0:.01:1;
+xChord = xChord*chord;
+upperSurface = zeros(1,length(xChord));
+lowerSurface = zeros(1,length(xChord));
+
+for i=1:length(xChord)
+ upperSurface(i) = get_z(xChord(i)/chord,1)*chord;
+ lowerSurface(i) = get_z(xChord(i)/chord,0)*chord;
+end
+
+figure; hold on; axis equal; grid on;
+%plot(xChord,z_camber,'-')
+plot(xChord,upperSurface,'-k','linewidth',2)
+plot(xChord,lowerSurface,'-k','linewidth',2)
+plot([0 1],[0 0],'--k','linewidth',1)
+
+
+for i = 1:length(webTop)
+ vecX = [frontSpar*chord webTop(i).xStart webTop(i).xEnd];
+ vecZ = [0 webTop(i).zStart webTop(i).zEnd];
+ fill(vecX,vecZ,[0.9 0.9 0.9])
+end
+
+for i = 1:length(webBottom)
+ vecX = [frontSpar*chord webBottom(i).xStart webBottom(i).xEnd];
+ vecZ = [0 webBottom(i).zStart webBottom(i).zEnd];
+ fill(vecX,vecZ,[0.9 0.9 0.9])
+end
+
+for i = 1:length(webUpperNose)
+ vecX = [frontSpar*chord webUpperNose(i).xStart webUpperNose(i).xEnd];
+ vecZ = [0 webUpperNose(i).zStart webUpperNose(i).zEnd];
+ fill(vecX,vecZ,[0.7 0.9 1.0])
+end
+
+for i = 1:length(webLowerNose)
+ vecX = [frontSpar*chord webLowerNose(i).xStart webLowerNose(i).xEnd];
+ vecZ = [0 webLowerNose(i).zStart webLowerNose(i).zEnd];
+ fill(vecX,vecZ,[0.7 0.9 1.0])
+end
+
+ vecX = [frontSpar*chord sparCaps(3).posX sparCaps(4).posX];
+ vecZ = [0 sparCaps(3).posZ sparCaps(4).posZ];
+ fill(vecX,vecZ,[0.9 0.9 0.9])
+
+
+sparCapSize = 18;
+stringerSize = 18;
+plot([sparCaps(1).posX sparCaps(2).posX],[sparCaps(1).posZ sparCaps(2).posZ],'-k','linewidth',2)
+plot([sparCaps(3).posX sparCaps(4).posX],[sparCaps(3).posZ sparCaps(4).posZ],'-k','linewidth',2)
+plot([sparCaps.posX],[sparCaps.posZ],'.b','markersize',sparCapSize)
+plot([topStringers.posX],[topStringers.posZ],'.r','markersize',stringerSize)
+plot([bottomStringers.posX],[bottomStringers.posZ],'.r','markersize',stringerSize)
+plot([noseTopStringers.posX],[noseTopStringers.posZ],'.r','markersize',stringerSize)
+plot([noseBottomStringers.posX],[noseBottomStringers.posZ],'.r','markersize',stringerSize)
+plot(centroid.posX,centroid.posZ,'.k','markerSize',18)
+plot(sc.posX,sc.posZ,'.g','markersize',18)
Copyright 2019--2024 Marius PETER