First commit!

Changed coordinate lists into numpy arrays.

18 files changed, 2153 insertions, 0 deletions

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..c796902
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,4 @@
+# .gitignore
+**/__pycache__/
+**/log.txt
+save/
\ No newline at end of file
diff --git a/README.org b/README.org
new file mode 100644
index 0000000..caa04a9
--- /dev/null
+++ b/README.org
@@ -0,0 +1,7 @@
+#+TITLE: UCLA MAE 154B
+#+SUBTITLE: Spring 2019 Final Project
+
+This program enables the creation of NACA airfoils;
+the analysis of the airfoil's structural properties;
+the optimization via genetic algorithm of a population of airfoils;
+With the final objective of designing a lightweight FAR 23 compliant airfoil.
diff --git a/creator/fuselage.py b/creator/fuselage.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/creator/fuselage.py
diff --git a/creator/propulsion.py b/creator/propulsion.py
new file mode 100644
index 0000000..e69de29
--- /dev/null
+++ b/creator/propulsion.py
diff --git a/creator/wing.py b/creator/wing.py
new file mode 100644
index 0000000..4988cb5
--- /dev/null
+++ b/creator/wing.py
@@ -0,0 +1,375 @@
+"""
+The wing.py module contains class definitions for and various components
+we add to an airfoil (spars, stringers, and ribs).
+
+Classes:
+    Airfoil: instantiated with class method to provide coordinates to heirs.
+    Spar: inherits from Airfoil.
+    Stringer: also inherits from Airfoil.
+
+Functions:
+    plot_geom(airfoil): generates a 2D plot of the airfoil & any components.
+"""
+
+import sys
+import os.path
+import logging
+import numpy as np
+from math import sin, cos, atan
+import bisect as bi
+import matplotlib.pyplot as plt
+
+logging.basicConfig(filename='log.txt',
+                    level=logging.DEBUG,
+                    format='%(asctime)s - %(levelname)s - %(message)s')
+
+
+class Component:
+    """Basic component providing coordinates and tools."""
+
+    # TODO: define defaults in separate module
+    def __init__(self):
+        self.x = np.array([])
+        self.z = np.array([])
+        self.material = str()
+        self.mass = float()
+
+    def set_material(self, material):
+        """Set the component bulk material."""
+        self.material = material
+
+    def info_print(self, round):
+        """Print all the component's coordinates to the terminal."""
+        name = f'    CREATOR DATA FOR {str(self).upper()}    '
+        num_of_dashes = len(name)
+        print(num_of_dashes * '-')
+        print(name)
+        for k, v in self.__dict__.items():
+            if type(v) != list:
+                print('{}:\n'.format(k), v)
+                print(num_of_dashes * '-')
+        for k, v in self.__dict__.items():
+            if type(v) == list:
+                print('{}:\n'.format(k), np.around(v, round))
+        return None
+
+    def info_save(self, save_path, number):
+        """Save all the object's coordinates (must be full path)."""
+        file_name = f'{str(self).lower()}_{number}.txt'
+        full_path = os.path.join(save_path, file_name)
+        try:
+            with open(full_path, 'w') as sys.stdout:
+                self.info_print(6)
+                # This line required to reset behavior of sys.stdout
+                sys.stdout = sys.__stdout__
+                logging.debug(f'Successfully wrote to file {full_path}')
+        except IOError:
+            print(f'Unable to write {file_name} to specified directory.\n',
+                  'Was the full path passed to the function?')
+        return None
+
+
+class Airfoil(Component):
+    """This class represents a single NACA airfoil.
+
+    The coordinates are saved as two lists
+    for the x- and z-coordinates. The coordinates start at
+    the leading edge, travel over the airfoil's upper edge,
+    then loop back to the leading edge via the lower edge.
+
+    This method was chosen for easier future exports
+    to 3D CAD packages like SolidWorks, which can import such
+    geometry as coordinates written in a CSV file.
+    """
+
+    # TODO: default values in separate module
+    def __init__(self, chord, semi_span, material):
+        super().__init__()
+        # self.x = np.array([])
+        # self.z = np.array([])
+        # self.chord = chord
+        """Create airfoil from its chord and semi-span."""
+        self.chord = chord if chord > 20 else 20
+        if chord <= 20:
+            logging.debug('Chord too small, using minimum value of 20.')
+        self.semi_span = semi_span
+        self.material = material
+        self.naca_num = int()
+
+    def __str__(self):
+        return type(self).__name__
+
+    def add_naca(self, naca_num):
+        """Generate surface geometry for a NACA airfoil.
+
+        The nested functions perform the required steps to generate geometry,
+        and can be called to solve the geometry y-coordinate for any 'x' input.
+        Equation coefficients were retrieved from Wikipedia.org.
+
+        Parameters:
+        naca_num: 4-digit NACA wing
+
+        Return:
+        None
+        """
+        self.naca_num = naca_num
+        # Variables extracted from naca_num argument passed to the function
+        m = int(str(naca_num)[0]) / 100
+        p = int(str(naca_num)[1]) / 10
+        t = int(str(naca_num)[2:]) / 100
+        # x-coordinate of maximum camber
+        p_c = p * self.chord
+
+        def get_camber(x):
+            """
+            Returns camber z-coordinate from 1 'x' along the airfoil chord.
+            """
+            z_c = float()
+            if 0 <= x < p_c:
+                z_c = (m / (p**2)) * (2 * p * (x / self.chord) -
+                                      (x / self.chord)**2)
+            elif p_c <= x <= self.chord:
+                z_c = (m /
+                       ((1 - p)**2)) * ((1 - 2 * p) + 2 * p *
+                                        (x / self.chord) - (x / self.chord)**2)
+            return (z_c * self.chord)
+
+        def get_thickness(x):
+            """Return thickness from 1 'x' along the airfoil chord."""
+            x = 0 if x < 0 else x
+            z_t = 5 * t * self.chord * (+0.2969 *
+                                        (x / self.chord)**0.5 - 0.1260 *
+                                        (x / self.chord)**1 - 0.3516 *
+                                        (x / self.chord)**2 + 0.2843 *
+                                        (x / self.chord)**3 - 0.1015 *
+                                        (x / self.chord)**4)
+            return z_t
+
+        def get_theta(x):
+            dz_c = float()
+            if 0 <= x < p_c:
+                dz_c = ((2 * m) / p**2) * (p - x / self.chord)
+            elif p_c <= x <= self.chord:
+                dz_c = (2 * m) / ((1 - p)**2) * (p - x / self.chord)
+
+            theta = atan(dz_c)
+            return theta
+
+        def get_coord_u(x):
+            x = x - get_thickness(x) * sin(get_theta(x))
+            z = get_camber(x) + get_thickness(x) * cos(get_theta(x))
+            return (x, z)
+
+        def get_coord_l(x):
+            x = x + get_thickness(x) * sin(get_theta(x))
+            z = get_camber(x) - get_thickness(x) * cos(get_theta(x))
+            return (x, z)
+
+        # Densify x-coordinates 10 times for first 1/4 chord length
+        x_chord_25_percent = round(self.chord / 4)
+        x_chord = [i / 10 for i in range(x_chord_25_percent * 10)]
+        x_chord.extend(i for i in range(x_chord_25_percent, self.chord + 1))
+        # Generate our airfoil skin geometry from previous sub-functions
+        self.x_c = np.array([])
+        self.z_c = np.array([])
+        # Upper surface and camber line
+        for x in x_chord:
+            self.x_c = np.append(self.x_c, x)
+            self.z_c = np.append(self.z_c, get_camber(x))
+            self.x = np.append(self.x, get_coord_u(x)[0])
+            self.z = np.append(self.z, get_coord_u(x)[1])
+        # Lower surface
+        for x in x_chord[::-1]:
+            self.x = np.append(self.x, get_coord_l(x)[0])
+            self.z = np.append(self.z, get_coord_l(x)[1])
+        return None
+
+
+class Spar(Component):
+    """Contains a single spar's data."""
+    def __init__(self, airfoil, loc_percent, material):
+        """Set spar location as percent of chord length."""
+        super().__init__()
+        super().set_material(material)
+        self.cap_area = float()
+        loc = loc_percent * airfoil.chord
+        # bi.bisect_left: returns index of first value in airfoil.x > loc
+        # This ensures that spar geom intersects with airfoil geom.
+        # Spar upper coordinates
+        spar_u = bi.bisect_left(airfoil.x, loc) - 1
+        self.x = np.append(self.x, airfoil.x[spar_u])
+        self.z = np.append(self.z, airfoil.z[spar_u])
+        # Spar lower coordinates
+        spar_l = bi.bisect_left(airfoil.x[::-1], loc)
+        self.x = np.append(self.x, airfoil.x[-spar_l])
+        self.z = np.append(self.z, airfoil.z[-spar_l])
+        return None
+
+    def set_cap_area(self, cap_area):
+        self.cap_area = cap_area
+        return None
+
+    def set_mass(self, mass):
+        self.mass = mass
+        return None
+
+
+class Stringer(Component):
+    """Contains the coordinates of all stringers."""
+    def __init__(self):
+        super().__init__()
+        self.x_start = []
+        self.x_end = []
+        self.z_start = []
+        self.z_end = []
+        self.diameter = float()
+        self.area = float()
+
+    def add_coord(self, airfoil, spars, stringer_u_1, stringer_u_2,
+                  stringer_l_1, stringer_l_2):
+        """Add equally distributed stringers to four airfoil locations
+        (upper nose, lower nose, upper surface, lower surface).
+
+        Parameters:
+        airfoil_coord: packed airfoil coordinates
+        spar_coord: packed spar coordinates
+        stringer_u_1: upper nose number of stringers
+        stringer_u_2: upper surface number of stringers
+        stringer_l_1: lower nose number of stringers
+        stringer_l_2: lower surface number of stringers
+
+        Returns:
+        None
+        """
+
+        # Find distance between leading edge and first upper stringer
+        interval = spars.x[0][0] / (stringer_u_1 + 1)
+        # initialise first self.stringer_x at first interval
+        x = interval
+        # Add upper stringers from leading edge until first spar.
+        for _ in range(0, stringer_u_1):
+            # Index of the first value of airfoil.x > x
+            i = bi.bisect_left(airfoil.x, x)
+            self.x.append(airfoil.x[i])
+            self.z.append(airfoil.z[i])
+            x += interval
+            # Add upper stringers from first spar until last spar
+            # TODO: stringer placement if only one spar is created
+        interval = (airfoil.spar.x[-1][0] -
+                    airfoil.spar.x[0][0]) / (stringer_u_2 + 1)
+        x = interval + airfoil.spar.x[0][0]
+        for _ in range(0, stringer_u_2):
+            i = bi.bisect_left(airfoil.x, x)
+            self.x.append(airfoil.x[i])
+            self.z.append(airfoil.z[i])
+            x += interval
+
+        # Find distance between leading edge and first lower stringer
+        interval = airfoil.spar.x[0][1] / (stringer_l_1 + 1)
+        x = interval
+        # Add lower stringers from leading edge until first spar.
+        for _ in range(0, stringer_l_1):
+            i = bi.bisect_left(airfoil.x[::-1], x)
+            self.x.append(airfoil.x[-i])
+            self.z.append(airfoil.z[-i])
+            x += interval
+            # Add lower stringers from first spar until last spar
+        interval = (airfoil.spar.x[-1][1] -
+                    airfoil.spar.x[0][1]) / (stringer_l_2 + 1)
+        x = interval + airfoil.spar.x[0][1]
+        for _ in range(0, stringer_l_2):
+            i = bi.bisect_left(airfoil.x[::-1], x)
+            self.x.append(airfoil.x[-i])
+            self.z.append(airfoil.z[-i])
+            x += interval
+        return None
+
+    def add_area(self, area):
+        self.area = area
+        return None
+
+    def add_mass(self, mass):
+        self.mass = len(self.x) * mass + len(self.x) * mass
+        return None
+
+    def add_webs(self, thickness):
+        """Add webs to stringers."""
+        for _ in range(len(self.x) // 2):
+            self.x_start.append(self.x[_])
+            self.x_end.append(self.x[_ + 1])
+            self.z_start.append(self.z[_])
+            self.z_end.append(self.z[_ + 1])
+            self.thickness = thickness
+        return None
+
+    def info_print(self, round):
+        super().info_print(round)
+        print('Stringer Area:\n', np.around(self.area, round))
+        return None
+
+
+def plot_geom(airfoil, spars, stringers):
+    """This function plots the airfoil's + sub-components' geometry."""
+    fig, ax = plt.subplots()
+
+    # Plot chord
+    x = [0, airfoil.chord]
+    y = [0, 0]
+    ax.plot(x, y, linewidth='1')
+    # Plot quarter chord
+    ax.plot(airfoil.chord / 4,
+            0,
+            '.',
+            color='g',
+            markersize=24,
+            label='Quarter-chord')
+    # Plot mean camber line
+    ax.plot(airfoil.x_c,
+            airfoil.z_c,
+            '-.',
+            color='r',
+            linewidth='2',
+            label='Mean camber line')
+    # Plot airfoil surfaces
+    ax.plot(airfoil.x, airfoil.z, color='b', linewidth='1')
+
+    # Plot spars
+    try:
+        for spar in spars:
+            x = (spar.x)
+            y = (spar.z)
+            ax.plot(x, y, '-', color='y', linewidth='4')
+    except AttributeError:
+        print('No spars to plot.')
+        # Plot stringers
+    try:
+        for _ in range(0, len(airfoil.stringer.x)):
+            x = airfoil.stringer.x[_]
+            y = airfoil.stringer.z[_]
+            ax.plot(x, y, '.', color='y', markersize=12)
+    except AttributeError:
+        print('No stringers to plot.')
+
+    # Graph formatting
+    # plot_bound = np.amax(airfoil.x)
+    ax.set(
+        title='NACA ' + str(airfoil.naca_num) + ' airfoil',
+        xlabel='X axis',
+        # xlim=[-0.10 * plot_bound, 1.10 * plot_bound],
+        ylabel='Z axis')
+    # ylim=[-(1.10 * plot_bound / 2), (1.10 * plot_bound / 2)])
+
+    plt.grid(axis='both', linestyle=':', linewidth=1)
+    plt.gca().set_aspect('equal', adjustable='box')
+    plt.gca().legend(bbox_to_anchor=(1, 1),
+                     bbox_transform=plt.gcf().transFigure)
+    plt.show()
+    return fig, ax
+
+
+def main():
+    return None
+
+
+if __name__ == '__main__':
+    main()
diff --git a/evaluator.py b/evaluator.py
new file mode 100644
index 0000000..85325a9
--- /dev/null
+++ b/evaluator.py
@@ -0,0 +1,284 @@
+"""
+The evaluator.py module contains a single Evaluator class,
+which knows all the attributes of a specified Airfoil instance,
+and contains functions to analyse the airfoil's geometrical
+& structural properties.
+"""
+
+import sys
+import os.path
+import numpy as np
+from math import sqrt
+import matplotlib.pyplot as plt
+
+
+class Evaluator:
+    """Performs structural evaluations for the airfoil passed as argument."""
+    def __init__(self, airfoil):
+        # Evaluator knows all geometrical info from evaluated airfoil
+        self.airfoil = airfoil
+        self.spar = airfoil.spar
+        self.stringer = airfoil.stringer
+        # Global dimensions
+        self.chord = airfoil.chord
+        self.semi_span = airfoil.semi_span
+        # Mass & spanwise distribution
+        self.mass_total = float(airfoil.mass + airfoil.spar.mass +
+                                airfoil.stringer.mass)
+        self.mass_dist = []
+        # Lift
+        self.lift_rectangular = []
+        self.lift_elliptical = []
+        self.lift_total = []
+        # Drag
+        self.drag = []
+        # centroid
+        self.centroid = []
+        # Inertia terms:
+        self.I_ = {'x': 0, 'z': 0, 'xz': 0}
+
+    def __str__(self):
+        return type(self).__name__
+
+    def info_print(self, round):
+        """Print all the component's evaluated data to the terminal."""
+        name = '    EVALUATOR DATA FOR {}    '.format(str(self).upper())
+        num_of_dashes = len(name)
+        print(num_of_dashes * '-')
+        print(name)
+        for k, v in self.__dict__.items():
+            if type(v) != list:
+                print('{}:\n'.format(k), v)
+        print(num_of_dashes * '-')
+        for k, v in self.__dict__.items():
+            if type(v) == list:
+                print('{}:\n'.format(k), np.around(v, round))
+        return None
+
+    def info_save(self, save_path, number):
+        """Save all the object's coordinates (must be full path)."""
+        file_name = 'airfoil_{}_eval.txt'.format(number)
+        full_path = os.path.join(save_path, file_name)
+        try:
+            with open(full_path, 'w') as sys.stdout:
+                self.info_print(6)
+                # This line required to reset behavior of sys.stdout
+                sys.stdout = sys.__stdout__
+                print('Successfully wrote to file {}'.format(full_path))
+        except IOError:
+            print(
+                'Unable to write {} to specified directory.\n'.format(
+                    file_name), 'Was the full path passed to the function?')
+        return None
+
+    # All these functions take integer arguments and return lists.
+
+    def get_lift_rectangular(self, lift):
+        L_prime = [lift / (self.semi_span * 2) for x in range(self.semi_span)]
+        return L_prime
+
+    def get_lift_elliptical(self, L_0):
+        L_prime = [
+            L_0 / (self.semi_span * 2) * sqrt(1 - (y / self.semi_span)**2)
+            for y in range(self.semi_span)
+        ]
+        return L_prime
+
+    def get_lift_total(self):
+        F_z = [(self.lift_rectangular[_] + self.lift_elliptical[_]) / 2
+               for _ in range(len(self.lift_rectangular))]
+        return F_z
+
+    def get_mass_distribution(self, total_mass):
+        F_z = [total_mass / self.semi_span for x in range(0, self.semi_span)]
+        return F_z
+
+    def get_drag(self, drag):
+        # Transform semi-span integer into list
+        semi_span = [x for x in range(0, self.semi_span)]
+
+        # Drag increases after 80% of the semi_span
+        cutoff = round(0.8 * self.semi_span)
+
+        # Drag increases by 25% after 80% of the semi_span
+        F_x = [drag for x in semi_span[0:cutoff]]
+        F_x.extend([1.25 * drag for x in semi_span[cutoff:]])
+        return F_x
+
+    def get_centroid(self):
+        """Return the coordinates of the centroid."""
+        stringer_area = self.stringer.area
+        cap_area = self.spar.cap_area
+
+        caps_x = [value for spar in self.spar.x for value in spar]
+        caps_z = [value for spar in self.spar.z for value in spar]
+        stringers_x = self.stringer.x
+        stringers_z = self.stringer.z
+
+        denominator = float(
+            len(caps_x) * cap_area + len(stringers_x) * stringer_area)
+
+        centroid_x = float(
+            sum([x * cap_area for x in caps_x]) +
+            sum([x * stringer_area for x in stringers_x]))
+        centroid_x = centroid_x / denominator
+
+        centroid_z = float(
+            sum([z * cap_area for z in caps_z]) +
+            sum([z * stringer_area for z in stringers_z]))
+        centroid_z = centroid_z / denominator
+
+        return (centroid_x, centroid_z)
+
+    def get_inertia_terms(self):
+        """Obtain all inertia terms."""
+        stringer_area = self.stringer.area
+        cap_area = self.spar.cap_area
+
+        # Adds upper and lower components' coordinates to list
+        x_stringers = self.stringer.x
+        z_stringers = self.stringer.z
+        x_spars = self.spar.x[:][0] + self.spar.x[:][1]
+        z_spars = self.spar.z[:][0] + self.spar.z[:][1]
+        stringer_count = range(len(x_stringers))
+        spar_count = range(len(self.spar.x))
+
+        # I_x is the sum of the contributions of the spar caps and stringers
+        # TODO: replace list indices with dictionary value
+        I_x = sum([
+            cap_area * (z_spars[i] - self.centroid[1])**2 for i in spar_count
+        ])
+        I_x += sum([
+            stringer_area * (z_stringers[i] - self.centroid[1])**2
+            for i in stringer_count
+        ])
+
+        I_z = sum([
+            cap_area * (x_spars[i] - self.centroid[0])**2 for i in spar_count
+        ])
+        I_z += sum([
+            stringer_area * (x_stringers[i] - self.centroid[0])**2
+            for i in stringer_count
+        ])
+
+        I_xz = sum([
+            cap_area * (x_spars[i] - self.centroid[0]) *
+            (z_spars[i] - self.centroid[1]) for i in spar_count
+        ])
+        I_xz += sum([
+            stringer_area * (x_stringers[i] - self.centroid[0]) *
+            (z_stringers[i] - self.centroid[1]) for i in stringer_count
+        ])
+        return (I_x, I_z, I_xz)
+
+    def get_dx(self, component):
+        return [x - self.centroid[0] for x in component.x_start]
+
+    def get_dz(self, component):
+        return [x - self.centroid[1] for x in component.x_start]
+
+    def get_dP(self, xDist, zDist, V_x, V_z, area):
+        I_x = self.I_['x']
+        I_z = self.I_['z']
+        I_xz = self.I_['xz']
+        denom = float(I_x * I_z - I_xz**2)
+        z = float()
+        for _ in range(len(xDist)):
+            z += float(-area * xDist[_] * (I_x * V_x - I_xz * V_z) / denom -
+                       area * zDist[_] * (I_z * V_z - I_xz * V_x) / denom)
+        return z
+
+    def analysis(self, V_x, V_z):
+        """Perform all analysis calculations and store in class instance."""
+        self.drag = self.get_drag(10)
+        self.lift_rectangular = self.get_lift_rectangular(13.7)
+        self.lift_elliptical = self.get_lift_elliptical(15)
+        self.lift_total = self.get_lift_total()
+        self.mass_dist = self.get_mass_distribution(self.mass_total)
+        self.centroid = self.get_centroid()
+        self.I_['x'] = self.get_inertia_terms()[0]
+        self.I_['z'] = self.get_inertia_terms()[1]
+        self.I_['xz'] = self.get_inertia_terms()[2]
+        spar_dx = self.get_dx(self.spar)
+        spar_dz = self.get_dz(self.spar)
+        self.spar.dP_x = self.get_dP(spar_dx, spar_dz, V_x, 0,
+                                     self.spar.cap_area)
+        self.spar.dP_z = self.get_dP(spar_dx, spar_dz, 0, V_z,
+                                     self.spar.cap_area)
+        return None
+
+
+def plot_geom(evaluator):
+    """This function plots analysis results over the airfoil's geometry."""
+    # Plot chord
+    x_chord = [0, evaluator.chord]
+    y_chord = [0, 0]
+    plt.plot(x_chord, y_chord, linewidth='1')
+    # Plot quarter chord
+    plt.plot(evaluator.chord / 4,
+             0,
+             '.',
+             color='g',
+             markersize=24,
+             label='Quarter-chord')
+    # Plot airfoil surfaces
+    x = [0.98 * x for x in evaluator.airfoil.x]
+    y = [0.98 * z for z in evaluator.airfoil.z]
+    plt.fill(x, y, color='w', linewidth='1', fill=False)
+    x = [1.02 * x for x in evaluator.airfoil.x]
+    y = [1.02 * z for z in evaluator.airfoil.z]
+    plt.fill(x, y, color='b', linewidth='1', fill=False)
+
+    # Plot spars
+    try:
+        for _ in range(len(evaluator.spar.x)):
+            x = (evaluator.spar.x[_])
+            y = (evaluator.spar.z[_])
+            plt.plot(x, y, '-', color='b')
+    except AttributeError:
+        print('No spars to plot.')
+    # Plot stringers
+    try:
+        for _ in range(0, len(evaluator.stringer.x)):
+            x = evaluator.stringer.x[_]
+            y = evaluator.stringer.z[_]
+            plt.plot(x, y, '.', color='y', markersize=12)
+    except AttributeError:
+        print('No stringers to plot.')
+
+    # Plot centroid
+    x = evaluator.centroid[0]
+    y = evaluator.centroid[1]
+    plt.plot(x, y, '.', color='r', markersize=24, label='centroid')
+
+    # Graph formatting
+    plt.xlabel('X axis')
+    plt.ylabel('Z axis')
+
+    plot_bound = max(evaluator.airfoil.x)
+    plt.xlim(-0.10 * plot_bound, 1.10 * plot_bound)
+    plt.ylim(-(1.10 * plot_bound / 2), (1.10 * plot_bound / 2))
+    plt.gca().set_aspect('equal', adjustable='box')
+    plt.gca().legend()
+    plt.grid(axis='both', linestyle=':', linewidth=1)
+    plt.show()
+    return None
+
+
+def plot_lift(evaluator):
+    x = range(evaluator.semi_span)
+    y_1 = evaluator.lift_rectangular
+    y_2 = evaluator.lift_elliptical
+    y_3 = evaluator.lift_total
+    plt.plot(x, y_1, '.', color='b', markersize=4, label='Rectangular lift')
+    plt.plot(x, y_2, '.', color='g', markersize=4, label='Elliptical lift')
+    plt.plot(x, y_3, '.', color='r', markersize=4, label='Total lift')
+
+    # Graph formatting
+    plt.xlabel('Semi-span location')
+    plt.ylabel('Lift')
+
+    plt.gca().legend()
+    plt.grid(axis='both', linestyle=':', linewidth=1)
+    plt.show()
+    return None
diff --git a/example_airfoil.py b/example_airfoil.py
new file mode 100644
index 0000000..292c1e9
--- /dev/null
+++ b/example_airfoil.py
@@ -0,0 +1,82 @@
+"""This example illustrates the usage of creator, evaluator and generator.
+
+All the steps of airfoil creation & evaluation are detailed here;
+furthermore, the generator.py module contains certain presets
+(default airfoils).
+
+Create an airfoil;
+Evaluate an airfoil;
+Generate a population of airfoils & optimize.
+"""
+
+from resources import materials as mt
+from creator import wing, fuselage, propulsion
+# from evaluator import
+# from generator import
+
+import time
+start_time = time.time()
+
+# Airfoil dimensions (in)
+NACA_NUM = 2412
+
+# Thicknesses
+SPAR_THICKNESS = 0.4
+SKIN_THICKNESS = 0.1
+
+# Component masses (lbs)
+AIRFOIL_MASS = 10
+SPAR_MASS = 10
+STRINGER_MASS = 5
+
+# Area (sqin)
+SPAR_CAP_AREA = 0.3
+STRINGER_AREA = 0.1
+
+# Amount of stringers
+TOP_STRINGERS = 6
+BOTTOM_STRINGERS = 4
+NOSE_TOP_STRINGERS = 3
+NOSE_BOTTOM_STRINGERS = 5
+
+SAVE_PATH = '/home/blendux/Projects/Aircraft_Studio/save'
+
+# Create airfoil instance
+af = wing.Airfoil(68, 150, mt.aluminium)
+af.add_naca(NACA_NUM)
+# af.info_print(2)
+af.info_save(SAVE_PATH, 'foo_name')
+
+# Create spar instances
+af.spar1 = wing.Spar(af, 0.23, mt.aluminium)
+af.spar2 = wing.Spar(af, 0.57, mt.aluminium)
+# af.spar1.info_print(2)
+# af.spar2.info_print(2)
+af.spar1.info_save(SAVE_PATH, 'spar1')
+af.spar2.info_save(SAVE_PATH, 'spar2')
+
+# # Create stringer instance
+# af.stringer = wing.Stringer()
+# # Compute the stringer coordinates from their quantity in each zone
+# af.stringer.add_coord(af, [af.spar1, af.spar2], NOSE_TOP_STRINGERS, TOP_STRINGERS,
+# NOSE_BOTTOM_STRINGERS, BOTTOM_STRINGERS)
+# af.stringer.add_area(STRINGER_AREA)
+# af.stringer.add_webs(SKIN_THICKNESS)
+# af.stringer.info_print(2)
+# af.stringer.info_save(SAVE_PATH, 'foo_name')
+
+# Plot components with matplotlib
+wing.plot_geom(af, [af.spar1, af.spar2], None)
+
+# Evaluator object contains airfoil analysis results.
+# eval = evaluator.Evaluator(af)
+# The analysis is performed in the evaluator.py module.
+# eval.analysis(1, 1)
+# eval.info_print(2)
+# eval.info_save(SAVE_PATH, 'foo_name')
+# evaluator.plot_geom(eval)
+# evaluator.plot_lift(eval)
+
+# Final execution time
+final_time = time.time() - start_time
+print(f"--- {round(final_time, 4)}s seconds ---")
diff --git a/generator.py b/generator.py
new file mode 100644
index 0000000..0213828
--- /dev/null
+++ b/generator.py
@@ -0,0 +1,64 @@
+# This file is part of Marius Peter's airfoil analysis package (this program).
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+"""
+The generator.py module contains a single Population class,
+which represents a collection of randomized airfoils.
+"""
+
+from tools import creator
+
+
+def default_airfoil():
+    """Generate the default airfoil."""
+    airfoil = creator.Airfoil.from_dimensions(100, 200)
+    airfoil.add_naca(2412)
+    airfoil.add_mass(10)
+
+    airfoil.spar = creator.Spar()
+    airfoil.spar.add_coord(airfoil, 0.23)
+    airfoil.spar.add_coord(airfoil, 0.57)
+    airfoil.spar.add_spar_caps(0.3)
+    airfoil.spar.add_mass(10)
+    airfoil.spar.add_webs(0.4)
+
+    airfoil.stringer = creator.Stringer()
+    airfoil.stringer.add_coord(airfoil, 3, 6, 5, 4)
+    airfoil.stringer.add_area(0.1)
+    airfoil.stringer.add_mass(5)
+    airfoil.stringer.add_webs(0.1)
+
+    return airfoil
+
+
+class Population(creator.Airfoil):
+    """Collection of random airfoils."""
+
+    def __init__(self, size):
+        af = creator.Airfoil
+        # print(af)
+        self.size = size
+        self.gen_number = 0  # incremented for every generation
+
+    def mutate(self, prob_mt):
+        """Randomly mutate the genes of prob_mt % of the population."""
+
+    def crossover(self, prob_cx):
+        """Combine the genes of prob_cx % of the population."""
+
+    def reproduce(self, prob_rp):
+        """Pass on the genes of the fittest prob_rp % of the population."""
+
+    def fitness():
+        """Rate the fitness of an individual on a relative scale (0-100)"""
diff --git a/gui.py b/gui.py
new file mode 100644
index 0000000..2eea281
--- /dev/null
+++ b/gui.py
@@ -0,0 +1,81 @@
+from tools import creator, evaluator, generator
+# import creator
+# import evaluator
+# import generator
+import tkinter as tk
+import tkinter.ttk as ttk
+
+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."""
+
+    label = ttk.Label(parent, text=name)
+    entry = ttk.Entry(parent)
+    return label, entry
+
+
+def set_naca(name):
+    naca_num = name.get()
+    print(naca_num)
+
+
+def set_chord(name):
+    chord = name.get()
+    print(chord)
+
+
+def set_semi_span(name):
+    semi_span = name.get()
+    print(semi_span)
+
+
+# plot.get_tk_widget().pack()
+MainWindow().mainloop()
diff --git a/resources/materials.py b/resources/materials.py
new file mode 100644
index 0000000..480b518
--- /dev/null
+++ b/resources/materials.py
@@ -0,0 +1,8 @@
+aluminium = {
+    "name": "aluminium",
+    "category": "metal",
+    "density": 2.70,
+    "mod_young": 70E9,
+    "mod_shear": 26E9,
+    "mod_bulk": 76E9
+}
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)