diff --git a/aerodynamic_analysis/CMakeLists.txt b/aerodynamic_analysis/CMakeLists.txt
index 0c1ea812d0fc097b979adbd377364beadf19964c..f5c099c506ed947993058ad0259cb4c7971c6faa 100644
--- a/aerodynamic_analysis/CMakeLists.txt
+++ b/aerodynamic_analysis/CMakeLists.txt
@@ -30,6 +30,7 @@ set(MODULE_SOURCES
 	src/bwb/bwbcalculatePolarConfig.cpp
 	src/bwb/bwbDefaultStrategy.cpp
 	src/methods/viscDragRaymer.cpp
+	src/methods/viscDragForTandem.cpp
 	src/methods/liftingLinePolar.cpp
 	src/methods/liftingLineForTAW.cpp
 	src/methods/semiEmpiricalHighLiftAdaption.cpp
diff --git a/aerodynamic_analysis/src/methods/viscDragForTandem.cpp b/aerodynamic_analysis/src/methods/viscDragForTandem.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8624f5b1867eab0cad71b796a53398485aeaf1e4
--- /dev/null
+++ b/aerodynamic_analysis/src/methods/viscDragForTandem.cpp
@@ -0,0 +1,412 @@
+/*  Copyright (C) 2010-2024 Institute of Aerospace Systems, RWTH Aachen University - All rights reserved.
+    This file is part of UNICADO.
+
+    UNICADO 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.
+
+    UNICADO 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 UNICADO.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+/*
+ * Raymer, D. P. - Aircraft Design: A Conceptual Approach, American Institute of Aeronautics and Astronautics, 2018, ISBN: 978-1-62410-574-6
+ * p. 280 ff. (pdf p. 294 ff.) 
+ */
+
+#include "viscDragForTandem.h"
+
+#include <aixml/node.h>
+#include <aircraftGeometry2/fuselage.h>
+#include <aircraftGeometry2/airfoil_surface.h>
+#include <aircraftGeometry2/processing/measure.h>
+#include <moduleBasics/module.h>
+#include <moduleBasics/runtimeIO.h>
+#include <runtimeInfo/runtimeInfo.h>
+#include <moduleBasics/strategySelector.h>
+
+#include <filesystem>
+#include <vector>
+#include <memory>
+#include <algorithm>
+#include <exception>
+#include <iomanip>
+#include <string>
+#include <utility>  // make_pair
+
+#include "../polar.h"
+#include "generalMethods/auxFunctions.h"
+
+using std::endl;
+using std::exception;
+using std::fixed;
+using std::ofstream;
+using std::setprecision;
+using std::string;
+using std::stringstream;
+using std::isinf;
+using std::isnan;
+using std::vector;
+/* TODO:    -find way to substract fuselage from wetted surface wing
+            -find way to calculate fuselage upsweep angle
+            -check if wing drag is for entire wing or just right side
+*/
+viscDragForTandemInitialization::viscDragForTandemInitialization()
+    :
+    calculateWing(false),
+    calculateFuselage(false),
+    calculateVTP(false),
+    calculateBackWing(false),
+    calculateNacelles(false),
+    calculateMisc(false),
+    manualTransitionLocation(false),
+    transitionLocationWing(NAN),
+    transitionLocationBackWing(NAN),
+    transitionLocationVTP(NAN),
+    doCalibration(false),
+    CDvisc_CDsum(NAN),
+    CDvisc_CLfact(NAN),
+    CDvisc_CLexp(NAN),
+    CDvisc_CDsum_lowMa(NAN),
+    CDvisc_CLfact_lowMa(NAN),
+    CDvisc_CLexp_lowMa(NAN),
+    jetEngine(false) {
+}
+
+viscDragForTandem::viscDragForTandem(const std::shared_ptr<RuntimeIO>& rtIO, const viscDragForTandemInitialization& viscDragInit)
+    :
+    // the following variables need to be initialized externally
+    calculateWing(viscDragInit.calculateWing),
+    calculateFuselage(viscDragInit.calculateFuselage),
+    calculateVTP(viscDragInit.calculateVTP),
+    calculateBackWing(viscDragInit.calculateBackWing),
+    calculateNacelles(viscDragInit.calculateNacelles),
+    calculateMisc(viscDragInit.calculateMisc),
+    manualTransitionLocation(viscDragInit.manualTransitionLocation),
+    transitionLocationWing(viscDragInit.transitionLocationWing),
+    transitionLocationBackWing(viscDragInit.transitionLocationBackWing),
+    transitionLocationVTP(viscDragInit.transitionLocationVTP),
+    doCalibration(viscDragInit.doCalibration),
+    CDvisc_CDsum(viscDragInit.CDvisc_CDsum),
+    CDvisc_CLfact(viscDragInit.CDvisc_CLfact),
+    CDvisc_CLexp(viscDragInit.CDvisc_CLexp),
+    CDvisc_CDsum_lowMa(viscDragInit.CDvisc_CDsum_lowMa),
+    CDvisc_CLfact_lowMa(viscDragInit.CDvisc_CLfact_lowMa),
+    CDvisc_CLexp_lowMa(viscDragInit.CDvisc_CLexp_lowMa),
+    jetEngine(viscDragInit.jetEngine),
+    Sref_Wing(NAN),
+    Swet_Wing(NAN),
+    MAC_Wing(NAN),
+    avgTtoC_Wing(NAN),
+    mXtoC_Wing(NAN),
+    eta_ref_Wing(NAN),
+    phiAtmXtoC_Wing(NAN),
+    QC_Wing(NAN),
+    Swet_BackWing(NAN),
+    MAC_BackWing(NAN),
+    avgTtoC_BackWing(NAN),
+    mXtoC_BackWing(NAN),
+    eta_ref_BackWing(NAN),
+    phiAtmXtoC_BackWing(NAN),
+    QC_BackWing(NAN),
+    Swet_VTP(NAN),
+    MAC_VTP(NAN),
+    avgTtoC_VTP(NAN),
+    mXtoC_VTP(NAN),
+    eta_ref_VTP(NAN),
+    phiAtmXtoC_VTP(NAN),
+    QC_VTP(NAN),
+    Swet_Fuselage(NAN),
+    l_Fuselage(NAN),
+    maxWidthFuselage(NAN),
+    maxHeightFuselage(NAN),
+    QC_Fuselage(NAN),
+    fuselageUpsweepAngle(NAN),
+    Swet_Nacelles(NAN),
+    l_Nacelles(NAN),
+    maxWidthNacelles(NAN),
+    maxHeightNacelles(NAN),
+    QC_Nacelles(NAN),
+    factorProtuberances(NAN) {
+    if (calculateWing) {
+        theAircraftComponents.myWing = viscDragInit.theWing;
+    }
+    if (calculateFuselage) {
+        theAircraftComponents.myFuselage = viscDragInit.theFuselage;
+    }
+    if (calculateNacelles) {
+        theAircraftComponents.myNacelles = viscDragInit.theNacelles;
+    }
+    if (calculateBackWing) {
+        theAircraftComponents.myBackWing = viscDragInit.theBackWing;
+    }
+    if (calculateVTP) {
+        theAircraftComponents.myVTP = viscDragInit.theVTP;
+    }
+    if (!manualTransitionLocation) {
+        transitionLocationWing = 0.0;
+        transitionLocationBackWing = 0.0;
+        transitionLocationVTP = 0.0;
+    }
+}
+
+void viscDragForTandem::preprocessGeometry() {
+    myRuntimeInfo->debug << "viscDragForTandem: preprocess the aircraft geometry" << endl;
+    if (calculateWing) {
+        Sref_Wing = geom2::measure::reference_area(theAircraftComponents.myWing);
+        Swet_Wing = geom2::measure::area(theAircraftComponents.myWing);
+        MAC_Wing = geom2::measure::mean_aerodynamic_chord(theAircraftComponents.myWing);
+        avgTtoC_Wing = getAverageTtoC(theAircraftComponents.myWing);
+        mXtoC_Wing = 0.5;  // Assumption by Raymer for chordwise location of airfoil maximum thickness point for high-speed airfoils
+        eta_ref_Wing = 0.65;  // Reference position used for local sweep <- cant  find this value in the source
+        phiAtmXtoC_Wing = getLocalSweep(theAircraftComponents.myWing, eta_ref_Wing * getHalfSpan(theAircraftComponents.myWing), mXtoC_Wing, "RADIAN");
+        QC_Wing = 1.;  // Assumtion for high, mid, and well fitted low wing
+        /*
+        // Wetted Surface 
+        if (!myAcftPt->theFuselage.empty()) {
+            wFuselage = myAcftPt->theFuselage.at(0).getLocalWidth(myAcftPt->theWing.at(i).referencePtx);
+            hFuselage = myAcftPt->theFuselage.at(0).getLocalHeight(myAcftPt->theWing.at(i).referencePtx);
+            zFuselage = myAcftPt->theFuselage.at(0).getLocalOffset(myAcftPt->theWing.at(i).referencePtx);
+        }
+        // Check if wing object is located in fuselage
+        if (myAcftPt->theWing.at(i).referencePtz < (zFuselage + hFuselage / 2) &&
+                myAcftPt->theWing.at(i).referencePtz > (zFuselage - hFuselage / 2)) {  // If tailplane Z-coordinate inside fuselage
+            Swet.push_back(myAcftPt->theWing.at(i).getWingSurface() - 2.*(wFuselage * myAcftPt->theWing.at(i).l_i_right));
+        } else {
+            Swet.push_back(myAcftPt->theWing.at(i).getWingSurface());
+        }
+        */
+    }
+    if (calculateBackWing) { //TODO kayra: her türlü tandem icin yap bir tane,  damit nutzer nicht verwirrt wird.
+        Swet_BackWing = geom2::measure::area(theAircraftComponents.myBackWing);
+        MAC_BackWing = geom2::measure::mean_aerodynamic_chord(theAircraftComponents.myBackWing);
+        avgTtoC_BackWing = getAverageTtoC(theAircraftComponents.myBackWing);
+        mXtoC_BackWing = 0.4;  // Assumption by Raymer for chordwise location of airfoil maximum thickness point for BackWing
+        eta_ref_BackWing = 0.65;  // Reference position used for local sweep
+        phiAtmXtoC_BackWing = getLocalSweep(theAircraftComponents.myBackWing, eta_ref_BackWing * getHalfSpan(theAircraftComponents.myBackWing), mXtoC_BackWing, "RADIAN");
+        QC_BackWing = 1.04;  // Form factor BackWing
+    }
+    if (calculateVTP) {
+        Swet_VTP = geom2::measure::area(theAircraftComponents.myVTP);
+        MAC_VTP = geom2::measure::mean_aerodynamic_chord(theAircraftComponents.myVTP);
+        avgTtoC_VTP = getAverageTtoC(theAircraftComponents.myVTP);
+        mXtoC_VTP = 0.4;  // Assumption by Raymer for chordwise location of airfoil maximum thickness point for high-speed airfoils
+        eta_ref_VTP = 0.65;
+        phiAtmXtoC_VTP = getLocalSweep(theAircraftComponents.myVTP, eta_ref_VTP * getHalfSpan(theAircraftComponents.myVTP), mXtoC_VTP, "RADIAN");
+        QC_VTP = 1.04;  // Form factor VTP
+    }
+    if (calculateFuselage) {
+        Swet_Fuselage = geom2::measure::area(theAircraftComponents.myFuselage);
+        l_Fuselage = geom2::measure::length(theAircraftComponents.myFuselage);
+        maxWidthFuselage = geom2::measure::width_max(theAircraftComponents.myFuselage);
+        maxHeightFuselage = geom2::measure::height_max(theAircraftComponents.myFuselage);
+        QC_Fuselage = 1.0;
+        fuselageUpsweepAngle = 3 * PI / 180;  // be creative!!!
+    }
+    if (calculateNacelles) {
+        Swet_Nacelles = getNacelleWettedArea(theAircraftComponents.myNacelles.front());
+        l_Nacelles = geom2::measure::length(theAircraftComponents.myNacelles.front());
+        maxWidthNacelles = geom2::measure::width_max(theAircraftComponents.myNacelles.front());
+        maxHeightNacelles = geom2::measure::height_max(theAircraftComponents.myNacelles.front());
+        QC_Nacelles = 1.4;
+    }
+    // according to [Ray18] (chap. Aerodynamics - Leakage and Protuberances): 2-5% for jet, 5-10% for prop
+    if (jetEngine) {
+        factorProtuberances = 0.03;
+    } else {
+        factorProtuberances = 0.07;
+    }
+    myRuntimeInfo->out << "Viscous drag calculation by Raymer successfully initialized for:" << endl;
+    if (calculateWing)
+        myRuntimeInfo->out << " -wing" << endl;
+    if (calculateBackWing)
+        myRuntimeInfo->out << " -BackWing" << endl;
+    if (calculateVTP)
+        myRuntimeInfo->out << " -VTP" << endl;
+    if (calculateFuselage)
+        myRuntimeInfo->out << " -fuselage" << endl;
+    if (calculateNacelles)
+        myRuntimeInfo->out << " -nacelles" << endl;
+    if (calculateMisc)
+        myRuntimeInfo->out << " -miscellaneous" << endl;
+}
+
+std::vector<double> viscDragForTandem::calculateViscDrag(const Polar& aPolar, const double transitionLocationWing,
+                                                    const double transitionLocationBackWing, const double transitionLocationVTP) {
+    // [Ray18] p. 280 ff. (pdf p. 294 ff.)
+    if (calculateWing == false) {
+        std::stringstream errorMsg;
+        errorMsg << "No wing in viscous drag raymer. Unable to calculate drag";
+        throwError(__FILE__, __func__, __LINE__, errorMsg.str());
+    }
+    // vector with component order: 0 = Wing, 1 = BackWing, 2 = VTP, 3 = Fuselage, 4 = Nacelles, 5 = Miscellaneous
+    vector <double> theComponentDrags;
+    theComponentDrags.push_back(calculateViscDragWing(aPolar, transitionLocationWing));
+    // calculate BackWing drag if switch is on
+    if (calculateBackWing) {
+        theComponentDrags.push_back(calculateViscDragBackWing(aPolar, transitionLocationBackWing));
+    } else {
+        theComponentDrags.push_back(0.);
+    }
+    // calculate VTP drag if switch is on
+    if (calculateVTP) {
+        theComponentDrags.push_back(calculateViscDragVTP(aPolar, transitionLocationVTP));
+    } else {
+        theComponentDrags.push_back(0.);
+    }
+    // calculate fuselage drag if switch is on
+    if (calculateFuselage) {
+        theComponentDrags.push_back(calculateViscDragFuselage(aPolar));
+    } else {
+        theComponentDrags.push_back(0.);
+    }
+    // calculate nacelle drag if switch is on
+    if (calculateNacelles) {
+        theComponentDrags.push_back(calculateViscDragNacelles(aPolar));
+    } else {
+        theComponentDrags.push_back(0.);
+    }
+    // calculate miscellaneous and protuberance drag if switch is on
+    if (calculateMisc) {
+        theComponentDrags.push_back(calculateViscDragMisc(aPolar));
+        // according to [Ray18] (chap. Aerodynamics - Leakage and Protuberances): 2-5% for jet, 5-10% for prop (Page 289)
+        double CD_protub = std::accumulate(theComponentDrags.begin(), theComponentDrags.end(), 0.0) * factorProtuberances;
+        theComponentDrags.back() += CD_protub;
+    } else {
+        theComponentDrags.push_back(0.);
+    }
+    // maybe do a plausibility check here..
+    return theComponentDrags;
+}
+
+std::vector<double> viscDragForTandem::calculateCorrectedViscDrag(const std::vector<double>& componentDrags, const double aMachNumber, const double CL) {
+    /* Drag calibration. First introduced in commit [4eab519c2] based on experimental data in [Lam2010], 23.06.2010. After several changes, 
+       the current definition and values in config are from Risse in commit [2ec91580] (02.08.2015) */
+    std::vector <double> correctedComponentDrags = componentDrags;
+    if (doCalibration) {
+        if (aMachNumber > 0.5) {
+            correctedComponentDrags.front() += CDvisc_CDsum + CDvisc_CLfact * pow(CL, CDvisc_CLexp);
+        } else if (aMachNumber <= 0.5) {
+            correctedComponentDrags.front() += CDvisc_CDsum_lowMa + CDvisc_CLfact_lowMa * pow(CL, CDvisc_CLexp_lowMa);
+        }
+    }
+    return correctedComponentDrags;
+}
+
+double viscDragForTandem::calculateViscDragWing(const Polar& aPolar, const double transitionLocationWing) {
+    // [Ray18] p. 280 ff. (pdf p. 294 ff.)
+    myRuntimeInfo->debug << "      - Get viscous drag wing ..." << endl;
+    double CDvisc_Wing(0.);
+    double reynoldNumberWing = getReynoldsNumber(MAC_Wing, aPolar.altitude, aPolar.MachNumber);
+    double Cf;
+    double Cflaminar;
+    double FF;
+    // Friction Coefficient (turbulent) for the wing
+    Cf = (getFrictionCoefficient(aPolar.MachNumber, reynoldNumberWing));
+    // Friction Coefficient (laminar) for the wing
+    Cflaminar = (1.328 / sqrt(reynoldNumberWing));
+    // Form factor for the components (eq. 12.30)
+    FF = (1. + 0.6 / mXtoC_Wing * avgTtoC_Wing + 100.*pow(avgTtoC_Wing, 4.))
+                 * (1.34 * pow(aPolar.MachNumber, 0.18) * pow(cos(phiAtmXtoC_Wing), 0.28));
+    // fully turbulent for transitionLocation = 0.0
+    CDvisc_Wing = (Cflaminar * transitionLocationWing + Cf * (1. - transitionLocationWing))
+                                                            * FF * QC_Wing * Swet_Wing / Sref_Wing;
+    return CDvisc_Wing;
+}
+
+double viscDragForTandem::calculateViscDragBackWing(const Polar& aPolar, const double transitionLocationBackWing) {
+    // [Ray18] p. 280 ff. (pdf p. 294 ff.)
+    myRuntimeInfo->debug << "      - Get viscous drag BackWing ..." << endl;
+    double CDvisc_BackWing(0.0);
+    double reynoldNumberBackWing = getReynoldsNumber(MAC_BackWing, aPolar.altitude, aPolar.MachNumber);
+    double Cf;
+    double Cflaminar;
+    double FF;
+    // Friction Coefficient (turbulent) for the BackWing
+    Cf = getFrictionCoefficient(aPolar.MachNumber, reynoldNumberBackWing);
+    // Friction Coefficient (laminar) for the BackWing
+    Cflaminar = 1.328 / sqrt(reynoldNumberBackWing);
+    // Form factor for the BackWing (eq. 12.30)
+    FF = ((1. + 0.6 / mXtoC_BackWing * avgTtoC_BackWing + 100.*pow(avgTtoC_BackWing, 4.))
+                 * (1.34 * pow(aPolar.MachNumber, 0.18) * pow(cos(phiAtmXtoC_BackWing), 0.28)));
+    CDvisc_BackWing = (Cflaminar * transitionLocationBackWing + Cf * (1. - transitionLocationBackWing))
+                                                                * FF * QC_BackWing * Swet_BackWing / Sref_Wing;
+    return CDvisc_BackWing;
+}
+
+double viscDragForTandem::calculateViscDragVTP(const Polar& aPolar, const double transitionLocationVTP) {
+    // [Ray18] p. 280 ff. (pdf p. 294 ff.)
+    myRuntimeInfo->debug << "      - Get viscous drag VTP ..." << endl;
+    double CDvisc_VTP(0.0);
+    double reynoldNumberVTP = getReynoldsNumber(MAC_VTP, aPolar.altitude, aPolar.MachNumber);
+    double Cf;
+    double Cflaminar;
+    double FF;
+    // Friction Coefficient (turbulent) for component
+    Cf = getFrictionCoefficient(aPolar.MachNumber, reynoldNumberVTP);
+    Cflaminar = 1.328 / sqrt(reynoldNumberVTP);
+    // Form factor for components
+    FF = (1. + 0.6 / mXtoC_VTP * avgTtoC_VTP + 100.*pow(avgTtoC_VTP, 4.))
+                 * (1.34 * pow(aPolar.MachNumber, 0.18) * pow(cos(phiAtmXtoC_VTP), 0.28));
+    // Store components
+    CDvisc_VTP = (Cflaminar * transitionLocationVTP + Cf * (1. - transitionLocationVTP))
+                                                            * FF * QC_VTP * Swet_VTP / Sref_Wing;
+    return CDvisc_VTP;
+}
+
+double viscDragForTandem::calculateViscDragFuselage(const Polar& aPolar) {
+    // [Ray18] p. 280 ff. (pdf p. 294 ff.)
+    myRuntimeInfo->debug << "      - Get viscous drag Fuselage ..." << endl;
+    double CDvisc_Fuselage(0.0);
+    double reynoldsNumberFuselage = getReynoldsNumber(l_Fuselage, aPolar.altitude, aPolar.MachNumber);
+    double Cf;
+    double FF;
+    // Friction Coefficient (turbulent) for component
+    Cf = getFrictionCoefficient(aPolar.MachNumber, reynoldsNumberFuselage);
+    // Form factor for fuselage
+    double f_Fuselage = l_Fuselage / sqrt(maxHeightFuselage * maxWidthFuselage);
+    FF = 1. + (60. / pow(f_Fuselage, 3.)) + (f_Fuselage / 400.);
+    CDvisc_Fuselage = Cf * FF * QC_Fuselage * Swet_Fuselage / Sref_Wing;
+    return CDvisc_Fuselage;
+}
+
+double viscDragForTandem::calculateViscDragNacelles(const Polar& aPolar) {
+    // [Ray18] p. 280 ff. (pdf p. 294 ff.)
+    myRuntimeInfo->debug << "      - Get viscous drag Nacelles ..." << endl;
+    double CDvisc_Nacelles(0.0);
+    double reynoldsNumberNacelles = getReynoldsNumber(l_Nacelles, aPolar.altitude, aPolar.MachNumber);
+    double Cf;
+    double FF;
+    // Reynolds number of component
+    reynoldsNumberNacelles = getReynoldsNumber(l_Nacelles, aPolar.altitude, aPolar.MachNumber);
+    // Friction Coefficient (turbulent) for nacelles
+    Cf = getFrictionCoefficient(aPolar.MachNumber, reynoldsNumberNacelles);
+    // Form factor for components
+    double f_Nacelles = l_Nacelles / sqrt(maxHeightNacelles * maxWidthNacelles);
+    FF = 1 + (0.35 / f_Nacelles);
+    // Store components
+    CDvisc_Nacelles = theAircraftComponents.myNacelles.size() * Cf * FF * QC_Nacelles * Swet_Nacelles / Sref_Wing;
+    return CDvisc_Nacelles;
+}
+
+double viscDragForTandem::calculateViscDragMisc(const Polar& aPolar) {
+    // [Ray18] p. 280 ff. (pdf p. 294 ff.)
+    myRuntimeInfo->debug << "      - Get viscous drag Miscellaneous ..." << endl;
+    double CDvisc_misc = 0.0;
+    if (calculateFuselage) {
+        CDvisc_misc = 3.83 * pow(fuselageUpsweepAngle, 2.5) * (PI * 0.5 * maxWidthFuselage * 0.5 * maxHeightFuselage) / Sref_Wing;
+    }
+    return CDvisc_misc;
+}
+
+double viscDragForTandem::getFrictionCoefficient(double flightMach, double ReynoldsNumber) {
+    // [Ray18] p. 280 ff. (pdf p. 294 ff.)
+    // Friction Coefficient (turbulent) for dcomponent
+    return 0.455 / (pow(log10(ReynoldsNumber), 2.58) * pow(1 + 0.144 * pow(flightMach, 2.), 0.65));
+}
diff --git a/aerodynamic_analysis/src/methods/viscDragForTandem.h b/aerodynamic_analysis/src/methods/viscDragForTandem.h
new file mode 100644
index 0000000000000000000000000000000000000000..10a5b905b6e7d037a8ba647393d9e682712e205b
--- /dev/null
+++ b/aerodynamic_analysis/src/methods/viscDragForTandem.h
@@ -0,0 +1,193 @@
+/*  Copyright (C) 2010-2024 Institute of Aerospace Systems, RWTH Aachen University - All rights reserved.
+    This file is part of UNICADO.
+
+    UNICADO 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.
+
+    UNICADO 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 UNICADO.  If not, see <http://www.gnu.org/licenses/>.
+*/
+#ifndef CALCULATEPOLAR_SRC_METHODS_VISCDRAGFORTANDEM_H_
+#define CALCULATEPOLAR_SRC_METHODS_VISCDRAGFORTANDEM_H_
+
+#include <moduleBasics/runtimeIO.h>
+#include <aircraftGeometry2/airfoil_surface.h>
+#include <aircraftGeometry2/hull_surface.h>
+#include <aircraftGeometry2/fuselage.h>
+
+#include <memory>
+#include <vector>
+#include <string>
+
+#include "../polar.h"
+
+class viscDragForTandemInitialization {
+ public:
+    // geometry switches
+    bool calculateWing;                 /**< calculate viscous drag of the wing switch */
+    bool calculateFuselage;             /**< calculate viscous drag of the fuselage switch */
+    bool calculateVTP;                  /**< calculate viscous drag of the VTP switch */
+    bool calculateBackWing;                  /**< calcualte viscous drag of the  BackWing switch */
+    bool calculateNacelles;             /**< calculate viscous drag of the nacelles switch */
+    bool calculateMisc;                 /**< calculate miscellaneous drag switch */
+    // transition parameters
+    bool manualTransitionLocation;          /**< switch wether transition location is given manually or flow is considered turbulent over all surfaces */
+    double transitionLocationWing;          /**< relative chord position on wing, where the flow goes from laminar to turbulent [-] */
+    double transitionLocationBackWing;           /**< relative chord position on BackWing, where the flow goes from laminar to turbulent [-] */
+    double transitionLocationVTP;           /**< relative chord position on VTP, where the flow goes from laminar to turbulent [-] */
+    // calibration parameters
+    bool doCalibration;
+    double CDvisc_CDsum;
+    double CDvisc_CLfact;
+    double CDvisc_CLexp;
+    double CDvisc_CDsum_lowMa;
+    double CDvisc_CLfact_lowMa;
+    double CDvisc_CLexp_lowMa;
+    // other
+    bool jetEngine;                         /**< switch if aircraft flies with jet engines or propellers */
+    // geometry
+    geom2::MultisectionSurface<geom2::AirfoilSection> theWing;                      /**< geometry of the main wing */
+    geom2::MultisectionSurface<geom2::AirfoilSection> theBackWing;                       /**< geometry of the horzontal tail plane */
+    geom2::MultisectionSurface<geom2::AirfoilSection> theVTP;                       /**< geometry of the vertical tail plaine */
+    geom2::MultisectionSurface<geom2::PolygonSection> theFuselage;                  /**< geometry of the fuselage */
+    std::vector<geom2::MultisectionSurface<geom2::PolygonSection>> theNacelles;     /**< vector containing all nacelle geometries */
+    viscDragForTandemInitialization();
+    virtual ~viscDragForTandemInitialization() {}
+};
+class viscDragForTandem {
+ public:
+    // general variables
+    std::string config;
+    std::filesystem::path aircraftPath;
+    std::string airfoilPath;
+    // geometry switches
+    bool calculateWing;                 /**< calculate viscous drag of the wing switch */
+    bool calculateFuselage;             /**< calculate viscous drag of the fuselage switch */
+    bool calculateVTP;                  /**< calculate viscous drag of the VTP switch */
+    bool calculateBackWing;                  /**< calcualte viscous drag of the  BackWing switch */
+    bool calculateNacelles;             /**< calculate viscous drag of the nacelles switch */
+    bool calculateMisc;                 /**< calculate miscellaneous drag switch */
+    // transition parameters
+    bool manualTransitionLocation;
+    double transitionLocationWing;          /**< relative chord position on wing, where the flow goes from laminar to turbulent [-] */
+    double transitionLocationBackWing;    /**< relative chord position on BackWing, where the flow goes from laminar to turbulent [-] */
+    double transitionLocationVTP;           /**< relative chord position on VTP, where the flow goes from laminar to turbulent [-] */
+    // calibration parameters
+    bool doCalibration;
+    double CDvisc_CDsum;
+    double CDvisc_CLfact;
+    double CDvisc_CLexp;
+    double CDvisc_CDsum_lowMa;
+    double CDvisc_CLfact_lowMa;
+    double CDvisc_CLexp_lowMa;
+    // sonstiges
+    bool jetEngine;                         /**< switch if aircraft flies with jet engines or propellers */
+    // geometry
+    struct aircraftComponents {                                                         /**< struct which stores all possible aircraft geometry elements */
+        geom2::MultisectionSurface<geom2::AirfoilSection> myWing;                      /**< geometry of the main wing */
+        geom2::MultisectionSurface<geom2::AirfoilSection> myBackWing;                       /**< geometry of the horzontal tail plane */
+        geom2::MultisectionSurface<geom2::AirfoilSection> myVTP;                       /**< geometry of the vertical tail plaine */
+        geom2::MultisectionSurface<geom2::PolygonSection> myFuselage;                  /**< geometry of the fuselage */
+        std::vector<geom2::MultisectionSurface<geom2::PolygonSection>> myNacelles;     /**< vector containing all nacelle geometries */
+    };
+    aircraftComponents theAircraftComponents;
+    // geometry dependent wing parameters
+    double Sref_Wing;           /**< Reference area of the wing [m^2] */
+    double Swet_Wing;           /**< Wetted area of the wing [m^2] */
+    double MAC_Wing;            /**< Mean aerodynamic chord of the wing [m] */
+    double avgTtoC_Wing;        /**< Average thickness to chord ratio of the wing [-] */
+    double eta_ref_Wing;        /**< Relative span position used for local sweep [-] */
+    double phiAtmXtoC_Wing;     /**< Local sweep at reference position [rad] */
+    double mXtoC_Wing;          /**< Chordwise location of airfoil maximum thickness point [-] */
+    double QC_Wing;             /**< Form factor wing [-] */
+    // geometry dependent BackWing parameters
+    double Swet_BackWing;            /**< Wetted area of the BackWing [m^2] */       //TODO kayra: buradaki parametreleri de degistirmem lazim mi?
+    double MAC_BackWing;             /**< Mean aerodynamic chord of the BackWing [m] */
+    double avgTtoC_BackWing;         /**< Average thickness to chord ratio of the BackWing [-] */
+    double mXtoC_BackWing;           /**< Chordwise location of airfoil maximum thickness point [-] */
+    double eta_ref_BackWing;         /**< Relative span position used for local sweep [-] */
+    double phiAtmXtoC_BackWing;      /**< Local sweep at reference position [rad] */
+    double QC_BackWing;              /**< Form factor BackWing [-] */
+    // geometry dependent VTP parameters
+    double Swet_VTP;            /**< Wetted area of the VTP [m^2] */
+    double MAC_VTP;             /**< Mean aerodynamic chord of the VTP [m] */
+    double avgTtoC_VTP;         /**< Average thickness to chord ratio of the VTP [-] */
+    double mXtoC_VTP;           /**< Chordwise location of airfoil maximum thickness point [-] */
+    double eta_ref_VTP;         /**< Relative span position used for local sweep [-] */
+    double phiAtmXtoC_VTP;      /**< Local sweep at reference position [rad] */
+    double QC_VTP;              /**< Form factor VTP [-] */
+    // geometry dependent fuselage parameters
+    double Swet_Fuselage;       /**< Wetted area of the fuselage [m^2] */
+    double l_Fuselage;          /**< lenght of the fuselage [m] */
+    double maxWidthFuselage;    /**< maximum width of the fuselage [m] */
+    double maxHeightFuselage;   /**< maximum height of the fuselage [m] */
+    double QC_Fuselage;         /**< Form factor fuselage [-] */
+    double fuselageUpsweepAngle;  /**< Up sweep angle of the fuselage tail [rad] */
+    // geometry dependent nacelle parameters
+    double Swet_Nacelles;       /**< Wetted area of one nacelle [m^2] */
+    double l_Nacelles;          /**< Length of one nacelle [m] */
+    double maxWidthNacelles;    /**< Maximum width of one nacelle [m] */
+    double maxHeightNacelles;   /**< Maximum height of one nacelle [m] */
+    double QC_Nacelles;         /**< Form factor nacelles [-] */
+    double factorProtuberances; /**< protuberance factor of the aircraft [-] */
+    // calculate parameters that are only dependent from the geometry
+    void preprocessGeometry();
+    // calculations
+    /** \brief Overall function to calculate viscous drag
+     *  \param aPolar polar/flight condition for which the drag gets calculated
+     *  \param transitionLocationWing relative transition location where the flow turns turbulent on the wing
+     *  \param transitionLocationBackWing relative transition location where the flow turns turbulent on the BackWing
+     *  \param transitionLocationVTP relative transition location where the flow turns turbulent on the VTP
+     */   
+    std::vector <double> calculateViscDrag(const Polar& aPolar, const double transitionLocationWing,
+                                    const double transitionLocationBackWing, const double transitionLocationVTP);
+    /** \brief Overall function to calculate viscous drag corrections
+     *  \param componentDrags vector of the component drags which need to be corrected
+     *  \param aMachNumber mach number for the calibration
+     *  \param CL CL value for the calibration
+     */                                       
+    std::vector <double> calculateCorrectedViscDrag(const std::vector <double>& componentDrags, const double aMachNumber, const double CL);
+    /** \brief Function to calculate the wing component drag
+     *  \param aPolar polar/flight condition for which the drag gets calculated
+     *  \param transitionLocationWing relative transition location where the flow turns turbulent on the wing
+     */          
+    double calculateViscDragWing(const Polar& aPolar, const double transitionLocationWing);
+    /** \brief Function to calculate the BackWing component drag
+     *  \param aPolar polar/flight condition for which the drag gets calculated
+     *  \param transitionLocationBackWing relative transition location where the flow turns turbulent on the BackWing
+     */  
+    double calculateViscDragBackWing(const Polar& aPolar, const double transitionLocationBackWing);
+    /** \brief Function to calculate the VTP component drag
+     *  \param aPolar polar/flight condition for which the drag gets calculated
+     *  \param transitionLocationVTP relative transition location where the flow turns turbulent on the VTP
+     */  
+    double calculateViscDragVTP(const Polar& aPolar, const double transitionLocationVTP);
+    /** \brief Function to calculate the fuselage component drag
+     *  \param aPolar polar/flight condition for which the drag gets calculated
+     */  
+    double calculateViscDragFuselage(const Polar& aPolar);
+    /** \brief Function to calculate the nacelle component drag
+     *  \param aPolar polar/flight condition for which the drag gets calculated
+     */  
+    double calculateViscDragNacelles(const Polar& aPolar);
+    /** \brief Function to calculate the miscellaneous component drag
+     *  \param aPolar polar/flight condition for which the drag gets calculated
+     */     
+    double calculateViscDragMisc(const Polar& aPolar);
+    /** \brief Function to calculate the friction coefficient
+     *  \param flightMach mach number
+     *  \param ReynoldsNumber reynolds number
+     */ 
+    double getFrictionCoefficient(double flightMach, double ReynoldsNumber);
+    explicit viscDragForTandem(const std::shared_ptr<RuntimeIO>& rtIO, const viscDragForTandemInitialization& viscDragInit);
+    virtual ~viscDragForTandem() {}
+};
+
+#endif  // CALCULATEPOLAR_SRC_METHODS_VISCDRAGFORTANDEM_H_
diff --git a/aerodynamic_analysis/src/polar.h b/aerodynamic_analysis/src/polar.h
index c3389134478760354c394ce72fb9b88c894abc2c..145823abafc5ec5fba02cc3cc00761b901723dad 100644
--- a/aerodynamic_analysis/src/polar.h
+++ b/aerodynamic_analysis/src/polar.h
@@ -139,8 +139,9 @@ class tandemDefaultPolar : public Polar {
     std::vector<defaultPolarPoint> PointsMisc;
 
     //TODO kayra:  Two main wings instead of one wing + horizontal stab
-    std::vector<defaultPolarPointLiftSurface> PointsWing; //PointsFrontWing
-    std::vector<defaultPolarPointLiftSurface> PointsBackWing;
+    std::vector<defaultPolarPointLiftSurface> PointsFrontWing; //PointsFrontWing
+    std::vector<defaultPolarPointLiftSurface> PointsBackWing; //PointsBackWing
+    std::vector<defaultPolarPointLiftSurface> PointsStab;  //TODO kayra: gecici olarak diger strategydeki kodlari silmeden silme
 
     tandemDefaultPolar();
     virtual ~tandemDefaultPolar() {}
diff --git a/aerodynamic_analysis/src/tandem/tandemCalculatePolarConfig.cpp b/aerodynamic_analysis/src/tandem/tandemCalculatePolarConfig.cpp
index fcd1158bbf78e7042b3893324384f124e4310bd8..8b0f77980080e9d70c5d839b8a1294deb897286a 100644
--- a/aerodynamic_analysis/src/tandem/tandemCalculatePolarConfig.cpp
+++ b/aerodynamic_analysis/src/tandem/tandemCalculatePolarConfig.cpp
@@ -39,7 +39,7 @@ tandemCalculatePolarConfig::tandemCalculatePolarConfig(const node& config)
     // execution mode
     execution_mode(EndnodeReadOnly<std::string>("program_settings/execution_mode")),
     // lifting line settings
-    stepWidthCL(EndnodeReadOnly<double>("LiftingLineForTAW/stepWidthCL")),      //TODO kayra: LiftingLineForTandem gibi bir copy paste yapayim mi ayni yeri mi kullanayim
+    stepWidthCL(EndnodeReadOnly<double>("LiftingLineForTAW/stepWidthCL")),      
     // other settings
     CLModeAirfoilSelection(EndnodeReadOnly<std::string>("program_settings/CLModeAirfoilSelection")),
     setFlightConditionsMode(EndnodeReadOnly<std::string>("program_settings/FlightConditions/AdaptMachNumbersToCruiseReq")),
@@ -53,12 +53,12 @@ tandemCalculatePolarConfig::tandemCalculatePolarConfig(const node& config)
     delta_CM0(EndnodeReadOnly<double>("LiftingLineForTAW/PitchingMoment/delta_CM0")),
     delta_dCMdCL(EndnodeReadOnly<double>("LiftingLineForTAW/PitchingMoment/delta_dCMdCL")),
     // lifiting line for TAW induced drag corrections
-    indDragCtCorrForCalibration(EndnodeReadOnly<double>("LiftingLineForTAW/InducedDragCorrection/indDragCtCorrForCalibration")), //TODO kayra: bu correctionlardan iki tane tanimlama gerek var mi mesela yoksa ortak kullanilir mi  
+    indDragCtCorrForCalibration(EndnodeReadOnly<double>("LiftingLineForTAW/InducedDragCorrection/indDragCtCorrForCalibration")),  
     factorIndDragCleanPolar(EndnodeReadOnly<double>("LiftingLineForTAW/InducedDragCorrection/factorIndDragCleanPolar")),
     // visc drag raymer transition parameters
     manualTransition(EndnodeReadOnly<bool>("program_settings/ViscDragRaymer/ManualTransition")),
     manualTransitionWing(EndnodeReadOnly<double>("program_settings/ViscDragRaymer/TransitionLocationWing")),
-    manualTransitionStabilizer(EndnodeReadOnly<double>("program_settings/ViscDragRaymer/TransitionLocationStabilizer")), //TODO kayra:bunun degismesi lazim mesela ama araya back wing komple copy paste mi olsun
+    manualTransitionBackWing(EndnodeReadOnly<double>("program_settings/ViscDragRaymer/TransitionLocationStabilizer")), //TODO kayra:bunu stab mi wing mi yapayim sor
     manualTransitionFin(EndnodeReadOnly<double>("program_settings/ViscDragRaymer/TransitionLocationFin")),
     // visc drag raymer calibration parameters
     doViscDragCalibration(EndnodeReadOnly<bool>("program_settings/ViscDragRaymer/UseCalibration")),
@@ -83,8 +83,8 @@ tandemCalculatePolarConfig::tandemCalculatePolarConfig(const node& config)
     doTrim(EndnodeReadOnly<bool>("program_settings/Trim/DoTrim")),
     trimAdditionalCGPositions(EndnodeReadOnly<bool>("program_settings/Trim/TrimWithAdditionalCoGPositions")),
     trimHighLiftPolars(EndnodeReadOnly<bool>("program_settings/Trim/TrimHighLift")),
-    stabAngleGridMode(EndnodeReadOnly<std::string>("program_settings/Trim/StabAngleGrid")),
-    customStabAngleGrid(EndnodeReadOnly<std::string>("program_settings/Trim/customStabAngleGrid")),       //TODO kayra: HLW/nur Höhenrudder?(aeleron ile böyle birsey yapilabilir) zaten set degil mi empenage design icin
+    backWingAngleGridMode(EndnodeReadOnly<std::string>("program_settings/Trim/StabAngleGrid")),
+    customBackWingAngleGrid(EndnodeReadOnly<std::string>("program_settings/Trim/customStabAngleGrid")),       
     deltaTotalDragHL(EndnodeReadOnly<double>("program_settings/SemiEmpiricalHighLiftAdaptions/deltaTotalDrag")),
     factorDragHL(EndnodeReadOnly<double>("program_settings/SemiEmpiricalHighLiftAdaptions/factorDrag")),
     // calibration settings
@@ -102,7 +102,7 @@ tandemCalculatePolarConfig::tandemCalculatePolarConfig(const node& config)
     pathToLiftingLine.read(config);
     manualTransition.read(config);
     manualTransitionWing.read(config);
-    manualTransitionStabilizer.read(config);
+    manualTransitionBackWing.read(config);
     manualTransitionFin.read(config);
     doViscDragCalibration.read(config);
     CDvisc_CDsum.read(config);
@@ -124,8 +124,8 @@ tandemCalculatePolarConfig::tandemCalculatePolarConfig(const node& config)
     doTrim.read(config);
     trimAdditionalCGPositions.read(config);
     trimHighLiftPolars.read(config);
-    stabAngleGridMode.read(config);
-    customStabAngleGrid.read(config);
+    backWingAngleGridMode.read(config);
+    customBackWingAngleGrid.read(config);
     deltaTotalDragHL.read(config);
     factorDragHL.read(config);
     CM_corr_fuselage_mode.read(config);
@@ -170,19 +170,20 @@ void tandemCalculatePolarConfig::readflightConditions(const node& config, size_t
 
 void tandemCalculatePolarConfig::initializeTrimSettings() {
     if (doTrim.value()) {
-        if (stabAngleGridMode.value() == "mode_0") {
-            stabAngles = {-6.0, 0.0, 6.0};
-            neutralStabAngleID = 1;
-        } else if (stabAngleGridMode.value() == "mode_1") {
-            stabAngles = {-6.0, -4.0, -2.0, 0.0, 0.2, 0.4, 0.6};
-            neutralStabAngleID = 3;
-        } else if (stabAngleGridMode.value() == "mode_2") {
+        if (backWingAngleGridMode.value() == "mode_0") {
+            backWingAngles = {-6.0, 0.0, 6.0}; //ekside angle of attack daha fazla oluyor. TODO kayra
+            //backWingAngles = {90, 90, 90}; //TODO kayra:sadece fuselage ici kisma twist koyuyordu kodu degistirdim
+            neutralBackWingAngleID = 1;
+        } else if (backWingAngleGridMode.value() == "mode_1") {
+            backWingAngles = {-6.0, -4.0, -2.0, 0.0, 0.2, 0.4, 0.6};
+            neutralBackWingAngleID = 3;
+        } else if (backWingAngleGridMode.value() == "mode_2") {
             myRuntimeInfo->out << "read out custom angles not implemented yet, tschüss!" << endl;
             exit(1);
         }
     } else {
-        stabAngles = {0.0};
-        neutralStabAngleID = 0;
+        backWingAngles = {0.0};
+        neutralBackWingAngleID = 0;
     }
     if (doTrim.value() && trimAdditionalCGPositions.value()) {
         CGPositions = {"design", "forward", "aft"};
diff --git a/aerodynamic_analysis/src/tandem/tandemCalculatePolarConfig.h b/aerodynamic_analysis/src/tandem/tandemCalculatePolarConfig.h
index e2624b5a579875a7525c18dea0ee3bfe5c8e286d..ed64e8ed0c2fb3fffb75ea4eb2353180f2129be7 100644
--- a/aerodynamic_analysis/src/tandem/tandemCalculatePolarConfig.h
+++ b/aerodynamic_analysis/src/tandem/tandemCalculatePolarConfig.h
@@ -56,7 +56,7 @@ class tandemCalculatePolarConfig{
     // visc drag raymer
     EndnodeReadOnly<bool> manualTransition;
     EndnodeReadOnly<double> manualTransitionWing;
-    EndnodeReadOnly<double> manualTransitionStabilizer;
+    EndnodeReadOnly<double> manualTransitionBackWing;
     EndnodeReadOnly<double> manualTransitionFin;
     EndnodeReadOnly<bool> doViscDragCalibration;
     EndnodeReadOnly<double> CDvisc_CDsum;
@@ -80,8 +80,8 @@ class tandemCalculatePolarConfig{
     EndnodeReadOnly<bool> doTrim;
     EndnodeReadOnly<bool> trimAdditionalCGPositions;
     EndnodeReadOnly<bool> trimHighLiftPolars;
-    EndnodeReadOnly<std::string> stabAngleGridMode;
-    EndnodeReadOnly<std::string> customStabAngleGrid;
+    EndnodeReadOnly<std::string> backWingAngleGridMode;
+    EndnodeReadOnly<std::string> customBackWingAngleGrid;
     // semi empirical high lift adaptions
     EndnodeReadOnly<double> deltaTotalDragHL;
     EndnodeReadOnly<double> factorDragHL;
@@ -106,9 +106,9 @@ class tandemCalculatePolarConfig{
     std::string stabtoolFolderPath; /**< Path to Stabtool folder */
 
     // trim settings
-    std::vector<double> stabAngles; /**< vector of stabilizer angles in [deg] */
+    std::vector<double> backWingAngles; /**< vector of stabilizer angles in [deg] */
     std::vector<std::string> CGPositions;  /**< vector of CG positions, order: design, fwd, aft*/
-    size_t neutralStabAngleID;
+    size_t neutralBackWingAngleID;
 
     /* Polar Configuration */
     /** \brief  Class which contains configurations for the polar calculations
diff --git a/aerodynamic_analysis/src/tandem/tandemDefaultStrategy.cpp b/aerodynamic_analysis/src/tandem/tandemDefaultStrategy.cpp
index dc1861771f3eeb67e1875963cbd993f7676c5b8f..73eeb4e7edac0bace78885a8b9244a3bb4ccf793 100644
--- a/aerodynamic_analysis/src/tandem/tandemDefaultStrategy.cpp
+++ b/aerodynamic_analysis/src/tandem/tandemDefaultStrategy.cpp
@@ -41,7 +41,7 @@
 #include "../methods/generalMethods/auxFunctions.h"
 #include "../methods/liftingLinePolar.h"
 #include "../methods/liftingLineForTAW.h"
-#include "../methods/viscDragRaymer.h"
+#include "../methods/viscDragForTandem.h"
 #include "../methods/waveDragMason.h"
 #include "../methods/semiEmpiricalHighLiftAdaption.h"
 #include "../methods/trimInterpolation.h"
@@ -69,18 +69,18 @@ tandemDefaultStrategy::tandemDefaultStrategy(const std::shared_ptr<RuntimeIO>& r
     config_(config) {
 }
 
-void tandemDefaultStrategy::run() {
+void tandemDefaultStrategy::run() {         //TODO kayra: stab angle silinmesi grekebilir oan göre fkt übergabeleri de iptallenebilir,händisch angle geben und gucken ob ganze flügel dreht oder nicht
     // loop over all cg positions and stabilizer angles
     for (size_t CGID(0); CGID < config_->CGPositions.size(); ++CGID) {
-        for (size_t stabID(0); stabID < config_->stabAngles.size(); ++stabID) {
+        for (size_t backWingID(0); backWingID < config_->backWingAngles.size(); ++backWingID) {
             this->clearDataForNextCalculationRun();
             myRuntimeInfo-> out << endl;
             myRuntimeInfo->out << "Run aero calculations for " << config_->CGPositions.at(CGID) <<
-                            " CoG and stabilizer angle of " << config_->stabAngles.at(stabID) << " deg..." << endl;
+                            " CoG and stabilizer angle of " << config_->backWingAngles.at(backWingID) << " deg..." << endl;
             myRuntimeInfo-> out << endl;
-            this->createAircraftGeometry(config_->stabAngles.at(stabID));
-            this->runAeroCalculations(data_->CoGPositions.at(config_->CGPositions.at(CGID)), config_->stabAngles.at(stabID));
-            this->copyPolarsToCGMap(CGID, stabID);
+            this->createAircraftGeometry(config_->backWingAngles.at(backWingID));
+            this->runAeroCalculations(data_->CoGPositions.at(config_->CGPositions.at(CGID)), config_->backWingAngles.at(backWingID));
+            this->copyPolarsToCGMap(CGID, backWingID);
         }
     }
     // interpolate trimed polars, or just copy polars if trim is disabled
@@ -98,7 +98,7 @@ void tandemDefaultStrategy::run() {
     this->calibratePolars();
 }
 
-void tandemDefaultStrategy::createAircraftGeometry(const double stabAngle) {
+void tandemDefaultStrategy::createAircraftGeometry(const double backWingAngle) {
     auto get_position = [](const node &component) -> geom2::Point_3 {
         return geom2::Point_3{double{component.at("position/x/value")},
                               double{component.at("position/y/value")},
@@ -106,23 +106,34 @@ void tandemDefaultStrategy::createAircraftGeometry(const double stabAngle) {
     };
     std::filesystem::path file{rtIO_->acxml.getName()};
     std::shared_ptr<node> AcXML = aixml::openDocument(file);
-    // create wing
+    // create wings
     this->aircraft_xml = aixml::openDocument(rtIO_->acxmlAccess);
-    const auto wing_position = get_position(this->aircraft_xml->at("aircraft_exchange_file/component_design/wing"));
+    const auto tandem_wing_set_position = get_position(this->aircraft_xml->at("aircraft_exchange_file/component_design/wing"));
     geom2::WingFactory wingFactory{AcXML, rtIO_->getAirfoilDataDir()};
-    geom2::MultisectionSurface<geom2::AirfoilSection> mainWing = wingFactory.create("wing/specific/geometry/aerodynamic_surface@0");
-    auto offset = wing_position - CGAL::ORIGIN;
-    mainWing.origin += offset;
-    liftingSurfaces.push_back(mainWing);
-    // create HTP
+    // front wing
+    geom2::MultisectionSurface<geom2::AirfoilSection> frontWing = wingFactory.create("wing/specific/geometry/aerodynamic_surface@0");
+    auto offset = tandem_wing_set_position - CGAL::ORIGIN;
+    frontWing.origin += offset;
+    liftingSurfaces.push_back(frontWing);
+    // back wing
+    geom2::MultisectionSurface<geom2::AirfoilSection> backWing = wingFactory.create("wing/specific/geometry/aerodynamic_surface@1");
+    backWing.origin += offset; 
+    for (auto& section : backWing.sections) {
+        double designBackWingAngle = section.get_twist_angle();
+        section.set_twist_angle(designBackWingAngle + convertUnit(DEGREE, RADIAN, backWingAngle));
+    }
+    /* double designBackWingAngle = backWing.sections.front().get_twist_angle(); //TODO KAyra: ausprobierern mit stl für backwing //TODO kayra:denemek icin bu eskisi
+    backWing.sections.front().set_twist_angle(designBackWingAngle + convertUnit(DEGREE, RADIAN, backWingAngle)); */
+    liftingSurfaces.push_back(backWing);
+    /*/ create HTP  //TODO kayra: back wing icin olustur acilari
     const auto htp_position = get_position(this->aircraft_xml->at("aircraft_exchange_file/component_design/empennage"));
     geom2::WingFactory HTPFactory{AcXML, rtIO_->getAirfoilDataDir()};
     geom2::MultisectionSurface<geom2::AirfoilSection> HTP = HTPFactory.create("empennage/specific/geometry/aerodynamic_surface@1");
     offset = htp_position - CGAL::ORIGIN;
     HTP.origin += offset;
-    double designStabAngle = HTP.sections.front().get_twist_angle();
+    double designStabAngle = HTP.sections.front().get_twist_angle(); //TODO KAyra: ausprobierern mit stl für backwing
     HTP.sections.front().set_twist_angle(designStabAngle + convertUnit(DEGREE, RADIAN, stabAngle));
-    liftingSurfaces.push_back(HTP);
+    liftingSurfaces.push_back(HTP);*/
     // create fuselage
     geom2::FuselageFactory fuselageFactory{AcXML, rtIO_->getGeometryDir()};
     fuselage = fuselageFactory.create("fuselage/specific/geometry/fuselage@0");
@@ -132,13 +143,15 @@ void tandemDefaultStrategy::createAircraftGeometry(const double stabAngle) {
     VTP = VTPFactory.create("empennage/specific/geometry/aerodynamic_surface@0");
     offset = vtp_position - CGAL::ORIGIN;
     VTP.origin += offset;
-    geom2::Mesh wing_mesh = geom2::transform::to_mesh(mainWing);
+    geom2::Mesh frontWing_mesh = geom2::transform::to_mesh(frontWing);
+    geom2::Mesh backWing_mesh  = geom2::transform::to_mesh(backWing);
+    backWing_mesh += frontWing_mesh;
     geom2::Mesh vtp_mesh = geom2::transform::to_mesh(VTP);
-    vtp_mesh += wing_mesh;
-    geom2::Mesh htp_mesh = geom2::transform::to_mesh(HTP);
-    htp_mesh += vtp_mesh;
+    vtp_mesh += backWing_mesh;
+    //geom2::Mesh htp_mesh = geom2::transform::to_mesh(HTP);  //TODO kayra:silmem gerekicek
+    //htp_mesh += vtp_mesh;
     geom2::Mesh fuselage_mesh = geom2::transform::to_mesh(fuselage);
-    fuselage_mesh += htp_mesh;
+    fuselage_mesh += vtp_mesh;
     // create nacelles
     const auto nacelle_position = get_position(this->aircraft_xml->at("aircraft_exchange_file/component_design/propulsion"));
     offset = nacelle_position - CGAL::ORIGIN;
@@ -151,21 +164,21 @@ void tandemDefaultStrategy::createAircraftGeometry(const double stabAngle) {
         fuselage_mesh += aNacelle_mesh;
         nacelles.push_back(aNacelle);
     }
-    CGAL::IO::write_STL("TAW_mesh.stl", fuselage_mesh); //TODO kayra: cok yanlis suan:update geldi mi?
+    CGAL::IO::write_STL("tandemWing_mesh.stl", fuselage_mesh); //TODO kayra: cok yanlis suan:update geldi mi?
 }
 
 void tandemDefaultStrategy::runAeroCalculations(const Point CoGPosition, const double currentIStabAngle) {
-    // lifting line for TAW
-    LiLiforTAWInitialization myLiLiInit;
-    this->initializeLiLiforTAW(&myLiLiInit);
-    liftingLineForTAW* myLiftingLineForTAW = new liftingLineForTAW(rtIO_, config_->liftingLineFolderPath, myLiLiInit);
-    myLiftingLineForTAW->run(CoGPosition);
-    this->getPolarDataFromLiftingLine(*myLiftingLineForTAW);
+    // lifting line polar
+    LiLiPolarInitialization myLiLiInit;
+    this->initializeLiLiPolar(&myLiLiInit);
+    liftingLinePolar* myLiftingLinePolar = new liftingLinePolar(rtIO_, config_->liftingLineFolderPath, myLiLiInit);
+    myLiftingLinePolar->run(CoGPosition);
+    this->getPolarDataFromLiftingLine(*myLiftingLinePolar);
     // viscous drag by raymer
-    viscDragRaymerInitialization myViscDragInit;
-    this->initializeViscDragRaymer(&myViscDragInit);
-    viscDragRaymer* myViscDragRamyer = new viscDragRaymer(rtIO_, myViscDragInit);
-    this->calculateViscDragRaymer(myViscDragRamyer);
+    viscDragForTandemInitialization myViscDragInit;
+    this->initializeViscDragForTandem(&myViscDragInit);
+    viscDragForTandem* myViscDragForTandem = new viscDragForTandem(rtIO_, myViscDragInit);
+    this->calculateViscDragForTandem(myViscDragForTandem);
     // wave drag by mason
     waveDragInitialization myWaveDragInit;
     this->initializeWaveDragMason(&myWaveDragInit);
@@ -176,23 +189,23 @@ void tandemDefaultStrategy::runAeroCalculations(const Point CoGPosition, const d
     // high lift calculation
     semiEmpiricalHighLiftAdaption* mySemiEmpHLAdaption = new semiEmpiricalHighLiftAdaption(rtIO_, liftingSurfaces.front());
     this->copyWingDevicesAndStabilizer(mySemiEmpHLAdaption);
-    this->copyCleanPolarForHLCalculation(mySemiEmpHLAdaption, myLiftingLineForTAW);
+    this->copyCleanPolarForHLCalculation(mySemiEmpHLAdaption, myLiftingLinePolar);
     mySemiEmpHLAdaption->runHighLiftAdaption();
     this->calculateHighLiftPolars(mySemiEmpHLAdaption, CoGPosition.xCoordinate, currentIStabAngle);
 }
 
-void tandemDefaultStrategy::initializeLiLiforTAW(LiLiforTAWInitialization *myLiLiInit) {
+void tandemDefaultStrategy::initializeLiLiPolar(LiLiPolarInitialization *myLiLiInit) {
     myLiLiInit->CLModeAirfoilSelection = config_->CLModeAirfoilSelection.value();
     myLiLiInit->stepWidthCL = config_->stepWidthCL.value();
-    myLiLiInit->CM_corr_fuselage_mode = config_->CM_corr_fuselage_mode.value();
-    myLiLiInit->CM_corr_nacelle_mode = config_->CM_corr_nacelle_mode.value();
+    //myLiLiInit->CM_corr_fuselage_mode = config_->CM_corr_fuselage_mode.value(); TODO kayra: diese waren für lilifortaw.
+    //myLiLiInit->CM_corr_nacelle_mode = config_->CM_corr_nacelle_mode.value();
     myLiLiInit->delta_CM0 = config_->delta_CM0.value();
     myLiLiInit->delta_dCMdCL = config_->delta_dCMdCL.value();
     myLiLiInit->indDragCtCorrForCalibration = config_->indDragCtCorrForCalibration.value();
     myLiLiInit->factorDragCleanPolar = config_->factorIndDragCleanPolar.value();
     myLiLiInit->bestCruiseCL = data_->bestCruiseCL;
     myLiLiInit->initialMachCruise = data_->initialCruiseMachNumber;
-    myLiLiInit->engineType = data_->engineType;
+    //myLiLiInit->engineType = data_->engineType;
     for (size_t flightConditionID(0); flightConditionID < config_->flightConditions.size(); ++flightConditionID) {
         flightCondition* aLiLiFlightCondition = new flightCondition;
         aLiLiFlightCondition->Mach = config_->flightConditions.at(flightConditionID).Mach;
@@ -206,16 +219,16 @@ void tandemDefaultStrategy::initializeLiLiforTAW(LiLiforTAWInitialization *myLiL
         delete aLiLiFlightCondition;
     }
     myLiLiInit->liftingSurfaces = liftingSurfaces;
-    myLiLiInit->theFuselage = fuselage;
-    myLiLiInit->theNacelles = nacelles;
+    //myLiLiInit->theFuselage = fuselage;       //TODO kayra:sil
+    //myLiLiInit->theNacelles = nacelles;
 }
 
-void tandemDefaultStrategy::initializeViscDragRaymer(viscDragRaymerInitialization *myViscDragInit) {
+void tandemDefaultStrategy::initializeViscDragForTandem(viscDragForTandemInitialization *myViscDragInit) {
     // copy geometry
     myViscDragInit->theWing = liftingSurfaces.at(0);
     myViscDragInit->calculateWing = true;
-    myViscDragInit->theHTP = liftingSurfaces.at(1);
-    myViscDragInit->calculateHTP = true;
+    myViscDragInit->theBackWing = liftingSurfaces.at(1);
+    myViscDragInit->calculateBackWing = true;
     myViscDragInit->theVTP = VTP;
     myViscDragInit->calculateVTP = true;
     myViscDragInit->theFuselage = fuselage;
@@ -226,7 +239,7 @@ void tandemDefaultStrategy::initializeViscDragRaymer(viscDragRaymerInitializatio
     // initialize transition parameters
     myViscDragInit->manualTransitionLocation = config_->manualTransition.value();
     myViscDragInit->transitionLocationWing = config_->manualTransitionWing.value();
-    myViscDragInit->transitionLocationHTP = config_->manualTransitionStabilizer.value();
+    myViscDragInit->transitionLocationBackWing = config_->manualTransitionBackWing.value();
     myViscDragInit->transitionLocationVTP = config_->manualTransitionFin.value();
     // initialize calibration parameters
     myViscDragInit->doCalibration = config_->doViscDragCalibration.value();
@@ -252,59 +265,66 @@ void tandemDefaultStrategy::initializeWaveDragMason(waveDragInitialization *myWa
     myWaveDragInit->liftingSurfaces = liftingSurfaces;
 }
 
-void tandemDefaultStrategy::getPolarDataFromLiftingLine(const liftingLineForTAW& liftingLineForTAW) {
+void tandemDefaultStrategy::getPolarDataFromLiftingLine(const liftingLinePolar& liftingLinePolar) {
     for (size_t machID(0); machID < config_->flightConditions.size(); ++machID) {
-        tawDefaultPolar aCleanPolar;
-        aCleanPolar.extrapolationMargin = liftingLineForTAW.theCleanPolars.at(machID).extrapolationMargin;
-        aCleanPolar.allowGridChange = liftingLineForTAW.theCleanPolars.at(machID).allowGridChange;
-        aCleanPolar.altitude = liftingLineForTAW.theCleanPolars.at(machID).altitude;
-        aCleanPolar.MachNumber = liftingLineForTAW.theCleanPolars.at(machID).MachNumber;
-        aCleanPolar.ReynoldsNumber = liftingLineForTAW.theCleanPolars.at(machID).ReynoldsNumber;
-        for (size_t pointID(0); pointID < liftingLineForTAW.theCleanPolars.at(machID).Points.size(); ++pointID) {
+        tandemDefaultPolar aCleanPolar;
+        aCleanPolar.extrapolationMargin = liftingLinePolar.theCleanPolars.at(machID).extrapolationMargin;
+        aCleanPolar.allowGridChange = liftingLinePolar.theCleanPolars.at(machID).allowGridChange;
+        aCleanPolar.altitude = liftingLinePolar.theCleanPolars.at(machID).altitude;
+        aCleanPolar.MachNumber = liftingLinePolar.theCleanPolars.at(machID).MachNumber;
+        aCleanPolar.ReynoldsNumber = liftingLinePolar.theCleanPolars.at(machID).ReynoldsNumber;
+        for (size_t pointID(0); pointID < liftingLinePolar.theCleanPolars.at(machID).Points.size(); ++pointID) {
             aCleanPolar.Points.push_back(defaultPolarPoint());
-            aCleanPolar.Points.back().CL = liftingLineForTAW.theCleanPolars.at(machID).Points.at(pointID).CL;
-            aCleanPolar.Points.back().AoA = liftingLineForTAW.theCleanPolars.at(machID).Points.at(pointID).AoA;
-            aCleanPolar.Points.back().CD_ind = liftingLineForTAW.theCleanPolars.at(machID).Points.at(pointID).CDind;
-            aCleanPolar.Points.back().CM = liftingLineForTAW.theCleanPolars.at(machID).Points.at(pointID).CM;
-            aCleanPolar.Points.back().iStabPolar = liftingLineForTAW.theCleanPolars.at(machID).Points.at(pointID).iStabPolar;
-            aCleanPolar.PointsWing.push_back(defaultPolarPointLiftSurface());
-            aCleanPolar.PointsWing.back().yStations = liftingLineForTAW.theCleanPolars.at(machID).liftingSurfacePolars.at(0).at(pointID).yStations;
-            aCleanPolar.PointsWing.back().Cl_distr = liftingLineForTAW.theCleanPolars.at(machID).liftingSurfacePolars.at(0).at(pointID).Cl_distr;
-            aCleanPolar.PointsWing.back().Lift_distr = liftingLineForTAW.theCleanPolars.at(machID).liftingSurfacePolars.at(0).at(pointID).Lift_distr;
-            aCleanPolar.PointsWing.back().CL = liftingLineForTAW.theCleanPolars.at(machID).liftingSurfacePolars.at(0).at(pointID).CL;
-            aCleanPolar.PointsWing.back().CD_ind = liftingLineForTAW.theCleanPolars.at(machID).liftingSurfacePolars.at(0).at(pointID).CDind;
-            aCleanPolar.PointsStab.push_back(defaultPolarPointLiftSurface());
-            aCleanPolar.PointsStab.back().yStations = liftingLineForTAW.theCleanPolars.at(machID).liftingSurfacePolars.at(1).at(pointID).yStations;
-            aCleanPolar.PointsStab.back().Cl_distr = liftingLineForTAW.theCleanPolars.at(machID).liftingSurfacePolars.at(1).at(pointID).Cl_distr;
-            aCleanPolar.PointsStab.back().Lift_distr = liftingLineForTAW.theCleanPolars.at(machID).liftingSurfacePolars.at(1).at(pointID).Lift_distr;
-            aCleanPolar.PointsStab.back().CL = liftingLineForTAW.theCleanPolars.at(machID).liftingSurfacePolars.at(1).at(pointID).CL;
-            aCleanPolar.PointsStab.back().CD_ind = liftingLineForTAW.theCleanPolars.at(machID).liftingSurfacePolars.at(1).at(pointID).CDind;
-            aCleanPolar.PointsNacelles.push_back(defaultPolarPoint());
-            aCleanPolar.PointsNacelles.back().CD_ind = liftingLineForTAW.theCleanPolars.at(machID).PointsNacelles.at(pointID).CDind;
+            aCleanPolar.Points.back().CL = liftingLinePolar.theCleanPolars.at(machID).Points.at(pointID).CL;
+            aCleanPolar.Points.back().AoA = liftingLinePolar.theCleanPolars.at(machID).Points.at(pointID).AoA;
+            aCleanPolar.Points.back().CD_ind = liftingLinePolar.theCleanPolars.at(machID).Points.at(pointID).CDind;
+            aCleanPolar.Points.back().CM = liftingLinePolar.theCleanPolars.at(machID).Points.at(pointID).CM;
+            aCleanPolar.Points.back().iStabPolar = liftingLinePolar.theCleanPolars.at(machID).Points.at(pointID).iStabPolar;
+            aCleanPolar.PointsFrontWing.push_back(defaultPolarPointLiftSurface());
+            aCleanPolar.PointsFrontWing.back().yStations = liftingLinePolar.theCleanPolars.at(machID).liftingSurfacePolars.at(0).at(pointID).yStations;
+            aCleanPolar.PointsFrontWing.back().Cl_distr = liftingLinePolar.theCleanPolars.at(machID).liftingSurfacePolars.at(0).at(pointID).Cl_distr;
+            aCleanPolar.PointsFrontWing.back().Lift_distr = liftingLinePolar.theCleanPolars.at(machID).liftingSurfacePolars.at(0).at(pointID).Lift_distr;
+            aCleanPolar.PointsFrontWing.back().CL = liftingLinePolar.theCleanPolars.at(machID).liftingSurfacePolars.at(0).at(pointID).CL;
+            aCleanPolar.PointsFrontWing.back().CD_ind = liftingLinePolar.theCleanPolars.at(machID).liftingSurfacePolars.at(0).at(pointID).CDind;
+            aCleanPolar.PointsBackWing.push_back(defaultPolarPointLiftSurface());
+            aCleanPolar.PointsBackWing.back().yStations = liftingLinePolar.theCleanPolars.at(machID).liftingSurfacePolars.at(0).at(pointID).yStations;
+            aCleanPolar.PointsBackWing.back().Cl_distr = liftingLinePolar.theCleanPolars.at(machID).liftingSurfacePolars.at(0).at(pointID).Cl_distr;
+            aCleanPolar.PointsBackWing.back().Lift_distr = liftingLinePolar.theCleanPolars.at(machID).liftingSurfacePolars.at(0).at(pointID).Lift_distr;
+            aCleanPolar.PointsBackWing.back().CL = liftingLinePolar.theCleanPolars.at(machID).liftingSurfacePolars.at(0).at(pointID).CL;
+            aCleanPolar.PointsBackWing.back().CD_ind = liftingLinePolar.theCleanPolars.at(machID).liftingSurfacePolars.at(0).at(pointID).CDind;
+            /*aCleanPolar.PointsStab.push_back(defaultPolarPointLiftSurface());
+            aCleanPolar.PointsStab.back().yStations = liftingLinePolar.theCleanPolars.at(machID).liftingSurfacePolars.at(1).at(pointID).yStations;
+            aCleanPolar.PointsStab.back().Cl_distr = liftingLinePolar.theCleanPolars.at(machID).liftingSurfacePolars.at(1).at(pointID).Cl_distr;
+            aCleanPolar.PointsStab.back().Lift_distr = liftingLinePolar.theCleanPolars.at(machID).liftingSurfacePolars.at(1).at(pointID).Lift_distr;
+            aCleanPolar.PointsStab.back().CL = liftingLinePolar.theCleanPolars.at(machID).liftingSurfacePolars.at(1).at(pointID).CL;
+            aCleanPolar.PointsStab.back().CD_ind = liftingLinePolar.theCleanPolars.at(machID).liftingSurfacePolars.at(1).at(pointID).CDind;*/
+            aCleanPolar.PointsNacelles.push_back(defaultPolarPoint());  
+            //TODO kayra:Even though liftingLine has no fuselage/nacelle data, i create an entry so that we can store our Raymer results later
+            //aCleanPolar.PointsNacelles.back().CD_ind = liftingLinePolar.theCleanPolars.at(machID).PointsNacelles.at(pointID).CDind;   
             aCleanPolar.PointsFuselage.push_back(defaultPolarPoint());
-            aCleanPolar.PointsFuselage.back().CD_ind = liftingLineForTAW.theCleanPolars.at(machID).PointsFuselage.at(pointID).CDind;
+            //aCleanPolar.PointsFuselage.back().CD_ind = liftingLinePolar.theCleanPolars.at(machID).PointsFuselage.at(pointID).CDind;
             aCleanPolar.PointsFin.push_back(defaultPolarPoint());
             aCleanPolar.PointsMisc.push_back(defaultPolarPoint());
         }
         currentPolarSet.cleanPolars.insert(std::make_pair(machID, aCleanPolar));
     }
-    LILI_Version = liftingLineForTAW.LILI_version;
+    LILI_Version = liftingLinePolar.LILI_version;
 }
 
-void tandemDefaultStrategy::calculateViscDragRaymer(viscDragRaymer *myViscDragRaymer) {
-    myViscDragRaymer->preprocessGeometry();
+void tandemDefaultStrategy::calculateViscDragForTandem(viscDragForTandem *myViscDragForTandem) {
+    myViscDragForTandem->preprocessGeometry();
     // run viscous drag calculation by Raymer over the current polar set
-    for (std::pair<size_t, tawDefaultPolar> polarIterator : currentPolarSet.cleanPolars) {
+    for (std::pair<size_t, tandemDefaultPolar> polarIterator : currentPolarSet.cleanPolars) {
         size_t polarID = polarIterator.first;
-        vector<double> theComponentDrags = myViscDragRaymer->calculateViscDrag(currentPolarSet.cleanPolars.at(polarID), 0.0, 0.0, 0.0);
+        vector<double> theComponentDrags = myViscDragForTandem->calculateViscDrag(currentPolarSet.cleanPolars.at(polarID), 0.0, 0.0, 0.0);
         // loop over all points to calculate CL dependent visc drag correction
         for (size_t pointID(0); pointID < currentPolarSet.cleanPolars.at(polarID).Points.size(); ++pointID) {
-            vector <double> theCorrectedComponentDrags = myViscDragRaymer->calculateCorrectedViscDrag(theComponentDrags, currentPolarSet.cleanPolars.at(polarID).MachNumber,
+            vector <double> theCorrectedComponentDrags = myViscDragForTandem->calculateCorrectedViscDrag(theComponentDrags, currentPolarSet.cleanPolars.at(polarID).MachNumber,
                                                                                                         currentPolarSet.cleanPolars.at(polarID).Points.at(pointID).CL);
             // assign total and componentwise viscous drag to the polar map
             currentPolarSet.cleanPolars.at(polarID).Points.at(pointID).CD_visc = std::accumulate(theCorrectedComponentDrags.begin(), theCorrectedComponentDrags.end(), 0.0);
-            currentPolarSet.cleanPolars.at(polarID).PointsWing.at(pointID).CD_visc = theCorrectedComponentDrags.at(0);
-            currentPolarSet.cleanPolars.at(polarID).PointsStab.at(pointID).CD_visc = theCorrectedComponentDrags.at(1);
+            currentPolarSet.cleanPolars.at(polarID).PointsFrontWing.at(pointID).CD_visc = theCorrectedComponentDrags.at(0);
+            currentPolarSet.cleanPolars.at(polarID).PointsBackWing.at(pointID).CD_visc = theCorrectedComponentDrags.at(1);
             currentPolarSet.cleanPolars.at(polarID).PointsFin.at(pointID).CD_visc = theCorrectedComponentDrags.at(2);
             currentPolarSet.cleanPolars.at(polarID).PointsFuselage.at(pointID).CD_visc = theCorrectedComponentDrags.at(3);
             currentPolarSet.cleanPolars.at(polarID).PointsNacelles.at(pointID).CD_visc = theCorrectedComponentDrags.at(4);
@@ -324,18 +344,18 @@ void tandemDefaultStrategy::calculateWaveDragMason(waveDragMason *myWaveDragMaso
             sweepPosition = 0.5;
         }
         for (size_t pointID(0); pointID < currentPolarSet.cleanPolars.at(polarID).Points.size(); ++pointID) {
-            currentPolarSet.cleanPolars.at(polarID).PointsWing.at(pointID).CD_wave
+            currentPolarSet.cleanPolars.at(polarID).PointsFrontWing.at(pointID).CD_wave
                                             = myWaveDragMason->calculateWaveDrag(currentPolarSet.cleanPolars.at(polarID).Points.at(pointID),
                                                                                 currentPolarSet.cleanPolars.at(polarID).MachNumber,
-                                                                                0, sweepPosition, currentPolarSet.cleanPolars.at(polarID).PointsWing.at(pointID).yStations,
-                                                                                currentPolarSet.cleanPolars.at(polarID).PointsWing.at(pointID).Cl_distr, true);
-            currentPolarSet.cleanPolars.at(polarID).PointsStab.at(pointID).CD_wave
+                                                                                0, sweepPosition, currentPolarSet.cleanPolars.at(polarID).PointsFrontWing.at(pointID).yStations,
+                                                                                currentPolarSet.cleanPolars.at(polarID).PointsFrontWing.at(pointID).Cl_distr, true);
+            currentPolarSet.cleanPolars.at(polarID).PointsBackWing.at(pointID).CD_wave
                                             = myWaveDragMason->calculateWaveDrag(currentPolarSet.cleanPolars.at(polarID).Points.at(pointID),
                                                                                 currentPolarSet.cleanPolars.at(polarID).MachNumber,
-                                                                                1, sweepPosition, currentPolarSet.cleanPolars.at(polarID).PointsStab.at(pointID).yStations,
-                                                                                currentPolarSet.cleanPolars.at(polarID).PointsStab.at(pointID).Cl_distr, false);
-            currentPolarSet.cleanPolars.at(polarID).Points.at(pointID).CD_wave = currentPolarSet.cleanPolars.at(polarID).PointsWing.at(pointID).CD_wave
-                                                                               + currentPolarSet.cleanPolars.at(polarID).PointsStab.at(pointID).CD_wave;
+                                                                                1, sweepPosition, currentPolarSet.cleanPolars.at(polarID).PointsBackWing.at(pointID).yStations,
+                                                                                currentPolarSet.cleanPolars.at(polarID).PointsBackWing.at(pointID).Cl_distr, false);
+            currentPolarSet.cleanPolars.at(polarID).Points.at(pointID).CD_wave = currentPolarSet.cleanPolars.at(polarID).PointsFrontWing.at(pointID).CD_wave
+                                                                               + currentPolarSet.cleanPolars.at(polarID).PointsBackWing.at(pointID).CD_wave;
         }
     }
     myRuntimeInfo->out << "... Wave drag has been calculated sucessfully!" << endl;
@@ -373,11 +393,11 @@ void tandemDefaultStrategy::copyWingDevicesAndStabilizer(semiEmpiricalHighLiftAd
     }
 }
 
-void tandemDefaultStrategy::copyCleanPolarForHLCalculation(semiEmpiricalHighLiftAdaption *mySemiEmpHLAdaption, liftingLineForTAW *myLiftingLineForTAW) {
+void tandemDefaultStrategy::copyCleanPolarForHLCalculation(semiEmpiricalHighLiftAdaption *mySemiEmpHLAdaption, liftingLinePolar *myLiftingLinePolar) {
     // find polar with minimum mach number
     double minMachNumber(1.0);
     size_t minMachID;
-    for (std::pair<size_t, tawDefaultPolar> polarIterator : currentPolarSet.cleanPolars) {
+    for (std::pair<size_t, tandemDefaultPolar> polarIterator : currentPolarSet.cleanPolars) {
         if (polarIterator.second.MachNumber < minMachNumber) {
             minMachID = polarIterator.first;
             minMachNumber = polarIterator.second.MachNumber;
@@ -404,22 +424,22 @@ void tandemDefaultStrategy::copyCleanPolarForHLCalculation(semiEmpiricalHighLift
     mySemiEmpHLAdaption->stepWidthCL = config_->stepWidthCL.value();
     mySemiEmpHLAdaption->reductionDragCountsHighLift = config_->deltaTotalDragHL.value();
     mySemiEmpHLAdaption->factorDragHighLift = config_->factorDragHL.value();
-    mySemiEmpHLAdaption->CDmin_inc = myLiftingLineForTAW->theAnalyticalPolarData.CDmin_inc.at(0);
-    mySemiEmpHLAdaption->kFactor_inc = myLiftingLineForTAW->theAnalyticalPolarData.kFactor_inc.at(0);
-    mySemiEmpHLAdaption->CLatCDmin_inc = myLiftingLineForTAW->theAnalyticalPolarData.CLatCDmin_inc.at(0);
-    mySemiEmpHLAdaption->dCMtodCL_inc = myLiftingLineForTAW->theAnalyticalPolarData.dCMtodCL_inc.at(0);
-    mySemiEmpHLAdaption->CMatCL0_inc = myLiftingLineForTAW->theAnalyticalPolarData.CMatCL0_inc.at(0);
-    mySemiEmpHLAdaption->dCMtodAoA = myLiftingLineForTAW->theAnalyticalPolarData.dCMtodAoA.at(0);
-    mySemiEmpHLAdaption->dCLtodAoA = myLiftingLineForTAW->theAnalyticalPolarData.dCLtodAoA.at(0);
-    mySemiEmpHLAdaption->CLatAoA0 = myLiftingLineForTAW->theAnalyticalPolarData.CLatAoA0.at(0);
-    mySemiEmpHLAdaption->CMatAoA0 = myLiftingLineForTAW->theAnalyticalPolarData.CMatAoA0.at(0);
+    mySemiEmpHLAdaption->CDmin_inc = myLiftingLinePolar->theAnalyticalPolarData.CDmin_inc.at(0);
+    mySemiEmpHLAdaption->kFactor_inc = myLiftingLinePolar->theAnalyticalPolarData.kFactor_inc.at(0);
+    mySemiEmpHLAdaption->CLatCDmin_inc = myLiftingLinePolar->theAnalyticalPolarData.CLatCDmin_inc.at(0);
+    mySemiEmpHLAdaption->dCMtodCL_inc = myLiftingLinePolar->theAnalyticalPolarData.dCMtodCL_inc.at(0);
+    mySemiEmpHLAdaption->CMatCL0_inc = myLiftingLinePolar->theAnalyticalPolarData.CMatCL0_inc.at(0);
+    mySemiEmpHLAdaption->dCMtodAoA = myLiftingLinePolar->theAnalyticalPolarData.dCMtodAoA.at(0);
+    mySemiEmpHLAdaption->dCLtodAoA = myLiftingLinePolar->theAnalyticalPolarData.dCLtodAoA.at(0);
+    mySemiEmpHLAdaption->CLatAoA0 = myLiftingLinePolar->theAnalyticalPolarData.CLatAoA0.at(0);
+    mySemiEmpHLAdaption->CMatAoA0 = myLiftingLinePolar->theAnalyticalPolarData.CMatAoA0.at(0);
 }
 
 void tandemDefaultStrategy::calculateHighLiftPolars(semiEmpiricalHighLiftAdaption *mySemiEmpHLAdaption, const double CoG_X_Position, const double currentIStabAngle) {
     for (size_t configID(0); configID < HLconfigs.size(); ++configID) {
         myRuntimeInfo->out << "Calculate high lift polar for " << HLconfigs.at(configID) << endl;
         Polar aPolar = mySemiEmpHLAdaption->calculateHighLiftPolar(HLconfigs.at(configID), CoG_X_Position, currentIStabAngle);
-        currentPolarSet.highLiftPolars[HLconfigs.at(configID)] = tawDefaultPolar();
+        currentPolarSet.highLiftPolars[HLconfigs.at(configID)] = tandemDefaultPolar();
         // copy polar data
         currentPolarSet.highLiftPolars.at(HLconfigs.at(configID)).Config = HLconfigs.at(configID);
         currentPolarSet.highLiftPolars.at(HLconfigs.at(configID)).extrapolationMargin = aPolar.extrapolationMargin;
@@ -454,7 +474,7 @@ void tandemDefaultStrategy::doTrimInterpolation() {
         // copy polar data into the generalized interpolation container
         for (size_t polarID(0); polarID < polarDataCGMap.at(CGID).at(0).cleanPolars.size(); ++polarID) {
             myRuntimeInfo->out << "Interpolate trim for clean polars at " << config_->CGPositions.at(CGID) << " CoG at Ma " << config_->machNumbers.at(polarID) << "..." << endl;
-            for (size_t stabID(0); stabID < config_->stabAngles.size(); ++stabID) {
+            for (size_t stabID(0); stabID < config_->backWingAngles.size(); ++stabID) {
                 vector<trimInterpolation::interpolationPoint> aStabPolar;
                 for (size_t pointID(0); pointID < polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(polarID).Points.size(); ++pointID) {
                     trimInterpolation::interpolationPoint anInterpolationPoint;
@@ -466,14 +486,14 @@ void tandemDefaultStrategy::doTrimInterpolation() {
                     anInterpolationPoint.values.push_back(polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(polarID).Points.at(pointID).CD_ind);
                     anInterpolationPoint.values.push_back(polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(polarID).Points.at(pointID).CD_visc);
                     anInterpolationPoint.values.push_back(polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(polarID).Points.at(pointID).CD_wave);
-                    anInterpolationPoint.values.push_back(polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(polarID).PointsWing.at(pointID).CD_wave);
+                    anInterpolationPoint.values.push_back(polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(polarID).PointsFrontWing.at(pointID).CD_wave);
                     anInterpolationPoint.values.push_back(polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(polarID).PointsStab.at(pointID).CD_ind);
                     anInterpolationPoint.values.push_back(polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(polarID).PointsStab.at(pointID).CD_wave);
                     // interpolate additional values for cruise mach polar
                     if (accuracyCheck(data_->initialCruiseMachNumber, polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(polarID).MachNumber, ACCURACY_HIGH)) {
-                        anInterpolationPoint.vectorValues.push_back(polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(polarID).PointsWing.at(pointID).Cl_distr);
-                        anInterpolationPoint.vectorValues.push_back(polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(polarID).PointsWing.at(pointID).Lift_distr);
-                        anInterpolationPoint.vectorValues.push_back(polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(polarID).PointsWing.at(pointID).yStations);
+                        anInterpolationPoint.vectorValues.push_back(polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(polarID).PointsFrontWing.at(pointID).Cl_distr);
+                        anInterpolationPoint.vectorValues.push_back(polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(polarID).PointsFrontWing.at(pointID).Lift_distr);
+                        anInterpolationPoint.vectorValues.push_back(polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(polarID).PointsFrontWing.at(pointID).yStations);
                     }
                     aStabPolar.push_back(anInterpolationPoint);
                 }
@@ -491,16 +511,16 @@ void tandemDefaultStrategy::doTrimInterpolation() {
             polarDataTrimMap.at(CGID).cleanPolars.at(polarID).allowGridChange = polarDataCGMap.at(CGID).at(0).cleanPolars.at(polarID).allowGridChange;
             for (size_t pointID(0); pointID < interpolatedTrimPolar.size(); ++pointID) {
                 polarDataTrimMap.at(CGID).cleanPolars.at(polarID).Points.push_back(defaultPolarPoint());
-                polarDataTrimMap.at(CGID).cleanPolars.at(polarID).PointsWing.push_back(defaultPolarPointLiftSurface());
+                polarDataTrimMap.at(CGID).cleanPolars.at(polarID).PointsFrontWing.push_back(defaultPolarPointLiftSurface());
                 polarDataTrimMap.at(CGID).cleanPolars.at(polarID).Points.back().CL = interpolatedTrimPolar.at(pointID).CL;
                 polarDataTrimMap.at(CGID).cleanPolars.at(polarID).Points.back().CM = interpolatedTrimPolar.at(pointID).CM;
                 polarDataTrimMap.at(CGID).cleanPolars.at(polarID).Points.back().CD = interpolatedTrimPolar.at(pointID).values.at(0);
                 polarDataTrimMap.at(CGID).cleanPolars.at(polarID).Points.back().AoA = interpolatedTrimPolar.at(pointID).values.at(1);
                 polarDataTrimMap.at(CGID).cleanPolars.at(polarID).Points.back().iStabPolar = interpolatedTrimPolar.at(pointID).values.at(2);
                 if (accuracyCheck(data_->initialCruiseMachNumber, polarDataTrimMap.at(CGID).cleanPolars.at(polarID).MachNumber, ACCURACY_HIGH)) {
-                    polarDataTrimMap.at(CGID).cleanPolars.at(polarID).PointsWing.back().Cl_distr = interpolatedTrimPolar.at(pointID).vectorValues.at(0);
-                    polarDataTrimMap.at(CGID).cleanPolars.at(polarID).PointsWing.back().Lift_distr = interpolatedTrimPolar.at(pointID).vectorValues.at(1);
-                    polarDataTrimMap.at(CGID).cleanPolars.at(polarID).PointsWing.back().yStations = interpolatedTrimPolar.at(pointID).vectorValues.at(2);
+                    polarDataTrimMap.at(CGID).cleanPolars.at(polarID).PointsFrontWing.back().Cl_distr = interpolatedTrimPolar.at(pointID).vectorValues.at(0);
+                    polarDataTrimMap.at(CGID).cleanPolars.at(polarID).PointsFrontWing.back().Lift_distr = interpolatedTrimPolar.at(pointID).vectorValues.at(1);
+                    polarDataTrimMap.at(CGID).cleanPolars.at(polarID).PointsFrontWing.back().yStations = interpolatedTrimPolar.at(pointID).vectorValues.at(2);
                 }
             }
         }
@@ -523,7 +543,7 @@ void tandemDefaultStrategy::copyCleanPolars() {
 void tandemDefaultStrategy::copyHighLiftPolars() {
     for (size_t CGID(0); CGID < config_->CGPositions.size(); ++CGID) {
         for (size_t configID(0); configID < HLconfigs.size(); ++configID) {
-            polarDataTrimMap[CGID].highLiftPolars[HLconfigs.at(configID)] = polarDataCGMap.at(CGID).at(config_->neutralStabAngleID).highLiftPolars.at(HLconfigs.at(configID));
+            polarDataTrimMap[CGID].highLiftPolars[HLconfigs.at(configID)] = polarDataCGMap.at(CGID).at(config_->neutralBackWingAngleID).highLiftPolars.at(HLconfigs.at(configID));
         }
     }
 }
@@ -644,7 +664,7 @@ void tandemDefaultStrategy::writePolarXML() {
             }
             XML << "            </polars>" << endl;
             XML << "        </configuration>" << endl;
-            for (const pair<const string, tawDefaultPolar>& thePolarIterator : polarDataTrimMap.at(CGID).highLiftPolars) {
+            for (const pair<const string, tandemDefaultPolar>& thePolarIterator : polarDataTrimMap.at(CGID).highLiftPolars) {
                 XML << "        <configuration ID=\"" << thePolarIterator.first << "\">" << endl;
                 XML << "            <name>" << thePolarIterator.first << "</name>" << endl;
                 XML << "            <comment>" << thePolarIterator.first << " Configuration</comment>" << endl;
@@ -696,11 +716,11 @@ void tandemDefaultStrategy::writePolarData() {
     for (size_t CGID(0); CGID < polarDataCGMap.size(); ++CGID) {
         for (const pair<const double, PolarData>& thePolarIterator : polarDataCGMap.at(CGID)) {
             /* plot all clean polars */
-            for (const pair<const int, tawDefaultPolar>& theMachIDIterator : thePolarIterator.second.cleanPolars) {
+            for (const pair<const int, tandemDefaultPolar>& theMachIDIterator : thePolarIterator.second.cleanPolars) {
                 stringstream machNumber;
                 stringstream iStabStream;
                 machNumber << setprecision(2) << fixed << theMachIDIterator.second.MachNumber;
-                iStabStream << setprecision(2) << fixed << Rounding(config_->stabAngles.at(thePolarIterator.first), 2);
+                iStabStream << setprecision(2) << fixed << Rounding(config_->backWingAngles.at(thePolarIterator.first), 2);
                 std::string csvFiles = "tawCalculatePolar_M" + replaceAll(machNumber.str(), ".", "") + "_"
                                        + CoGnames.at(CGID) + "_" + replaceAll(iStabStream.str(), ".", "") +"_plot.csv";
                 plot.open(csvFiles.c_str());
@@ -723,8 +743,8 @@ void tandemDefaultStrategy::writePolarData() {
                         plot << setprecision(6) << fixed << theMachIDIterator.second.Points.at(i).CM << ";";
                         plot << setprecision(6) << fixed << theMachIDIterator.second.Points.at(i).CL / theMachIDIterator.second.Points.at(i).CD << ";";
                         plot << setprecision(6) << fixed << theMachIDIterator.second.Points.at(i).iStabPolar << ";";
-                        plot << setprecision(6) << fixed << theMachIDIterator.second.PointsWing.at(i).CD_visc << ";";
-                        plot << setprecision(6) << fixed << theMachIDIterator.second.PointsWing.at(i).CD_wave << ";";
+                        plot << setprecision(6) << fixed << theMachIDIterator.second.PointsFrontWing.at(i).CD_visc << ";";
+                        plot << setprecision(6) << fixed << theMachIDIterator.second.PointsFrontWing.at(i).CD_wave << ";";
                         plot << setprecision(6) << fixed << theMachIDIterator.second.PointsStab.at(i).CD_visc << ";";
                         plot << setprecision(6) << fixed << theMachIDIterator.second.PointsStab.at(i).CD_wave << ";";
                         plot << setprecision(6) << fixed << theMachIDIterator.second.PointsStab.at(i).CD_ind << ";";
@@ -739,11 +759,11 @@ void tandemDefaultStrategy::writePolarData() {
     for (size_t CGID(0); CGID < polarDataCGMap.size(); ++CGID) {
         for (const pair<const double, PolarData>& thePolarIterator : polarDataCGMap.at(CGID)) {
             /* plot all clean polars */
-            for (const pair<const std::string, tawDefaultPolar>& theConfigIterator : thePolarIterator.second.highLiftPolars) {
+            for (const pair<const std::string, tandemDefaultPolar>& theConfigIterator : thePolarIterator.second.highLiftPolars) {
                 stringstream machNumber;
                 stringstream iStabStream;
                 machNumber << setprecision(2) << fixed << theConfigIterator.second.MachNumber;
-                iStabStream << setprecision(2) << fixed << Rounding(config_->stabAngles.at(thePolarIterator.first), 2);
+                iStabStream << setprecision(2) << fixed << Rounding(config_->backWingAngles.at(thePolarIterator.first), 2);
                 std::string csvFiles = "tawCalculatePolar_" + theConfigIterator.first + "_"
                                        + CoGnames.at(CGID) + "_" + replaceAll(iStabStream.str(), ".", "") +"_plot.csv";
                 plot.open(csvFiles.c_str());
@@ -786,9 +806,9 @@ void tandemDefaultStrategy::writePolarData() {
     } else {
         plot << "# " + csvFiles << endl;
         plot << "# (1)C_L; (2)C_Lw; (3)y, m; (4)C_l; (5)Cl*l, m" << endl;
-        for (size_t pointID(0); pointID < polarDataTrimMap.at(0).cleanPolars.at(cruiseMachID).PointsWing.size(); ++pointID) {
-            vector <double> clDistribution_tmp = polarDataTrimMap.at(0).cleanPolars.at(cruiseMachID).PointsWing.at(pointID).Cl_distr;
-            vector <double> yStations_tmp      = polarDataTrimMap.at(0).cleanPolars.at(cruiseMachID).PointsWing.at(pointID).yStations;
+        for (size_t pointID(0); pointID < polarDataTrimMap.at(0).cleanPolars.at(cruiseMachID).PointsFrontWing.size(); ++pointID) {
+            vector <double> clDistribution_tmp = polarDataTrimMap.at(0).cleanPolars.at(cruiseMachID).PointsFrontWing.at(pointID).Cl_distr;
+            vector <double> yStations_tmp      = polarDataTrimMap.at(0).cleanPolars.at(cruiseMachID).PointsFrontWing.at(pointID).yStations;
             reverse(clDistribution_tmp.begin(), clDistribution_tmp.end());
             reverse(yStations_tmp.begin(), yStations_tmp.end());
             // for Cl distr analysis at single CL
@@ -797,7 +817,7 @@ void tandemDefaultStrategy::writePolarData() {
                 string delimiter(";");
                 for (size_t j(0); j < clDistribution_tmp.size(); ++j) {
                     plot << polarDataTrimMap.at(0).cleanPolars.at(cruiseMachID).Points.at(pointID).CL << delimiter;
-                    plot << polarDataTrimMap.at(0).cleanPolars.at(cruiseMachID).PointsWing.at(pointID).CL << delimiter;
+                    plot << polarDataTrimMap.at(0).cleanPolars.at(cruiseMachID).PointsFrontWing.at(pointID).CL << delimiter;
                     plot << yStations_tmp.at(j) << delimiter;
                     plot << clDistribution_tmp.at(j) << delimiter;
                     plot << getLocalChordLength(liftingSurfaces.front(), yStations_tmp.at(j)) << endl;
@@ -840,7 +860,7 @@ void tandemDefaultStrategy::writeViscDragData() {
                 stringstream machNumber;
                 stringstream iStabStream;
                 machNumber << setprecision(2) << fixed << polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(machID).MachNumber;
-                iStabStream << setprecision(2) << fixed << Rounding(config_->stabAngles.at(stabID), 2);
+                iStabStream << setprecision(2) << fixed << Rounding(config_->backWingAngles.at(stabID), 2);
                 std::string csvFile = "viscDrag_M" + replaceAll(machNumber.str(), ".", "") + "_"
                                        + CoGnames.at(CGID) + "_" + replaceAll(iStabStream.str(), ".", "") +"_plot.csv";
                 plot.open(csvFile.c_str());
@@ -857,7 +877,7 @@ void tandemDefaultStrategy::writeViscDragData() {
                         plot << setprecision(3) << fixed << polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(machID).Points.at(i).AoA << ";";
                         plot << setprecision(3) << fixed << polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(machID).Points.at(i).CL << ";";
                         plot << setprecision(6) << fixed << polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(machID).Points.at(i).CD_visc << ";";
-                        plot << setprecision(6) << fixed << polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(machID).PointsWing.at(i).CD_visc << ";";
+                        plot << setprecision(6) << fixed << polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(machID).PointsFrontWing.at(i).CD_visc << ";";
                         plot << setprecision(6) << fixed << polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(machID).PointsStab.at(i).CD_visc << ";";
                         plot << setprecision(6) << fixed << polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(machID).PointsFin.at(i).CD_visc << ";";
                         plot << setprecision(6) << fixed << polarDataCGMap.at(CGID).at(stabID).cleanPolars.at(machID).PointsFuselage.at(i).CD_visc << ";";
@@ -876,7 +896,7 @@ void tandemDefaultStrategy::setReferenceValues() {
     data_->span.set_value(getHalfSpan(liftingSurfaces.at(0)) * 2);
     data_->MAC.set_value(geom2::measure::mean_aerodynamic_chord(liftingSurfaces.at(0)));
     data_->S_ref.set_value(geom2::measure::reference_area(liftingSurfaces.at(0)));
-    data_->neutral_point_xCoordinate.set_value(getNeutralPointXposition(getTrimmedCruisePolar(polarDataCGMap.at(0).at(0).cleanPolars, data_->initialCruiseMachNumber),
+    data_->neutral_point_xCoordinate.set_value(getNeutralPointXposition(getTrimmedTandemCruisePolar(polarDataCGMap.at(0).at(0).cleanPolars, data_->initialCruiseMachNumber),
                                                                         data_->MAC.value(), data_->CoGPositions["design"].xCoordinate));
 }
 
@@ -917,7 +937,7 @@ void tandemDefaultStrategy::update_xml_entries(std::string mode) {
 double tandemDefaultStrategy::getCLoptCruise() {
     double LoverDmax(0.0);
     double CLopt(0.0);
-    std::vector<defaultPolarPoint> theTrimmedCruisePolar(getTrimmedCruisePolar(polarDataTrimMap.at(0).cleanPolars, data_->initialCruiseMachNumber));
+    std::vector<defaultPolarPoint> theTrimmedCruisePolar(getTrimmedTandemCruisePolar(polarDataTrimMap.at(0).cleanPolars, data_->initialCruiseMachNumber));
     for (size_t pointID(0); pointID < theTrimmedCruisePolar.size(); ++pointID) {
         double LoverD = theTrimmedCruisePolar.at(pointID).CL / theTrimmedCruisePolar.at(pointID).CD;
             if (LoverD > LoverDmax) {
@@ -928,7 +948,7 @@ double tandemDefaultStrategy::getCLoptCruise() {
     return CLopt;
 }
 
-double tandemDefaultStrategy::getMaxCL(const tawDefaultPolar& aPolar) {
+double tandemDefaultStrategy::getMaxCL(const tandemDefaultPolar& aPolar) {
     double maxCL = -10000;
     for (size_t ID = 0; ID < aPolar.Points.size(); ++ID) {
         if (aPolar.Points.at(ID).CL > maxCL) {
diff --git a/aerodynamic_analysis/src/tandem/tandemDefaultStrategy.h b/aerodynamic_analysis/src/tandem/tandemDefaultStrategy.h
index efb51500ae619648fd5eddbcb43530fe2b38caf2..568fa0a001496937845826f37c4d151068edeb59 100644
--- a/aerodynamic_analysis/src/tandem/tandemDefaultStrategy.h
+++ b/aerodynamic_analysis/src/tandem/tandemDefaultStrategy.h
@@ -32,7 +32,7 @@
 #include "../methods/liftingLineForTAW.h"
 #include "../methods/liftingLinePolar.h"
 
-#include "../methods/viscDragRaymer.h"
+#include "../methods/viscDragForTandem.h"
 #include "../methods/waveDragMason.h"
 #include "../methods/semiEmpiricalHighLiftAdaption.h"
 
@@ -40,7 +40,7 @@ class tandemDefaultStrategy {
  public:
     // geometry bools
     // containers for the aircraft geometry
-    std::vector <geom2::MultisectionSurface<geom2::AirfoilSection>> liftingSurfaces;    /**< vector which stores all lifting surface geometries (wing and HTP)*/
+    std::vector <geom2::MultisectionSurface<geom2::AirfoilSection>> liftingSurfaces;    /**< vector which stores all lifting surface geometries (wings and HTP)*/
     geom2::MultisectionSurface<geom2::AirfoilSection> VTP;                              /**< the VTP geometry */
     geom2::MultisectionSurface<geom2::PolygonSection> fuselage;                         /**< the fuselage geometry */
     std::vector <geom2::MultisectionSurface<geom2::PolygonSection>> nacelles;           /**< vector which stores the nacelle geometries*/
@@ -75,14 +75,14 @@ class tandemDefaultStrategy {
     double C_LoptimumCruise;
 
     void run();
-    /** \brief Function initializes liftingLineForTAW with the required variables
-     *  \param myLiLiInit the lifting line initializer for TAW method with gets initialized
+    /** \brief Function initializes liftingLinePolar with the required variables
+     *  \param myLiLiInit the lifting line polar initializer method with gets initialized
      */
-    void initializeLiLiforTAW(LiLiforTAWInitialization *myLiLiInit);
+    void initializeLiLiPolar(LiLiPolarInitialization *myLiLiInit);
     /** \brief Function initializes viscDragRaymer with the required variables
      *  \param myViscDragRaymer the viscous drag raymer method with gets initialized
      */
-    void initializeViscDragRaymer(viscDragRaymerInitialization *myViscDragInit);
+    void initializeViscDragForTandem(viscDragForTandemInitialization *myViscDragInit);
     /** \brief Function initializes waveDragMason with the required variables
      *  \param myWaveDragMason the wave drag mason method with gets initialized
      */
@@ -90,26 +90,26 @@ class tandemDefaultStrategy {
     /** \brief Function with builds and saves the entire aircraft geometry
      *  \param stabAngle the angle of the stabilizer relative to the design angle, which gets altered for the trimm setting [deg]
      */
-    void createAircraftGeometry(const double stabAngle);
+    void createAircraftGeometry(const double backWingAngle);
     /** \brief Function witch executes all aero calculation methods
      *  \param CoGPosition the 3D center of gravity coordinate as reference for the aero moments [m]
      */    
     void runAeroCalculations(const Point CoGPosition, const double currentIStabAngle);
-    /** \brief Function to copy the polar set from liftingLineForTAW to the strategyClass
-     *  \param liftingLineForTAW the lifting line for TAW method, which generated the clean polar set
+    /** \brief Function to copy the polar set from liftingLinePolar to the strategyClass
+     *  \param liftingLinePolar the lifting line for TAW method, which generated the clean polar set
      */
-    void getPolarDataFromLiftingLine(const liftingLineForTAW& liftingLineForTAW);
+    void getPolarDataFromLiftingLine(const liftingLinePolar& liftingLinePolar);
     /** \brief Function to calculate the viscous drag for the current polar set
-     *  \param myViscDragRaymer the viscous drag raymer method to calculate the drag
+     *  \param myViscDragForTandem the viscous drag raymer method to calculate the drag
      */
-    void calculateViscDragRaymer(viscDragRaymer *myViscDragRaymer);
+    void calculateViscDragForTandem(viscDragForTandem *myViscDragForTandem);
     /** \brief Function to calculate the wave drag for the current polar set
      *  \param myWaveDragMason the wave drag mason method to calculate the drag
      */
     void calculateWaveDragMason(waveDragMason *myWaveDragMason);
     void sumUpDragComponents();
     void copyWingDevicesAndStabilizer(semiEmpiricalHighLiftAdaption *mySemiEmpHLAdaption);
-    void copyCleanPolarForHLCalculation(semiEmpiricalHighLiftAdaption *mySemiEmpHLAdaption, liftingLineForTAW *myLiftingLineForTAW);
+    void copyCleanPolarForHLCalculation(semiEmpiricalHighLiftAdaption *mySemiEmpHLAdaption, liftingLinePolar *myLiftingLinePolar);
     void calculateHighLiftPolars(semiEmpiricalHighLiftAdaption *mySemiEmpHLAdaption, const double CoG_X_Position, const double currentIStabAngle);
     /** \brief Copy the temporary current polar set to the permanent polarDataCGMap
      *  \param CGID index of the current center of gravity position
@@ -146,7 +146,7 @@ class tandemDefaultStrategy {
 
     double getCLoptCruise();
 
-    double getMaxCL(const tawDefaultPolar& aPolar);
+    double getMaxCL(const tandemDefaultPolar& aPolar);
 
     void setReferenceValues();
 
diff --git a/wing_design/src/tandem/cantilever/cantileverTandemWingDesign.cpp b/wing_design/src/tandem/cantilever/cantileverTandemWingDesign.cpp
index 65de7d509ef03d95c0de0b917d489c24cbc900b4..9729e2dd80761e720776f95f1a9bbb0d3868df48 100644
--- a/wing_design/src/tandem/cantilever/cantileverTandemWingDesign.cpp
+++ b/wing_design/src/tandem/cantilever/cantileverTandemWingDesign.cpp
@@ -53,8 +53,7 @@ Wing::Wing(const std::shared_ptr<RuntimeIO>& rtIO)
   /* Setup mass modes */
   mass_mode_runner["mode_0"] = [this]() { flops_mass(); };    //TODO kayra: buradan ayari
   mass_mode_runner["mode_1"] = [this]() { chiozzotto_wer_mass(); };
-  back_mass_mode_runner["mode_0"] = [this]() { back_flops_mass(); };    //TODO kayra: buradan ayari
-  back_mass_mode_runner["mode_1"] = [this]() { back_chiozzotto_wer_mass(); };
+
 }
 
 void Wing::initialize() {
@@ -75,7 +74,6 @@ void Wing::initialize() {
   /* Store method to variable */
   selected_design_mode = config->wing_design_mode.value();
   selected_mass_mode = config->wing_mass_mode.value();
-  selected_back_mass_mode = config->back_wing_mass_mode.value(); //TODO kayra :altta sikintisi var:bunlardan biri gerekmeyebilir
 
 
   airfoils_library = std::make_shared<Airfoils>(static_cast<std::filesystem::path>(config->common_airfoil_data_path.value()));
@@ -91,8 +89,7 @@ void Wing::run() {
   myRuntimeInfo->out << "Selected design method ... [FINISHED]" << std::endl;
 
   myRuntimeInfo->out << "Selected mass method   ... [START]" << std::endl;
-  mass_mode_runner[selected_mass_mode]();             //TODO kayra: burada bir eksik var back mass mode runner,
-  back_mass_mode_runner[selected_back_mass_mode]();            //belki ikisi icin ortak mass mode daha mantikli
+  mass_mode_runner[selected_mass_mode]();    
   myRuntimeInfo->out << "Selected mass method   ... [FINISHED]" << std::endl;
 
   myRuntimeInfo->out << "Run wing design ... [FINISHED]" << std::endl;
@@ -186,7 +183,7 @@ void Wing::standard_design() {
   myRuntimeInfo->out << "1/4 sweep of rear Wing: " << std::endl;    
   // Compute back wing quarter chord based on selected mode          
   const double installed_back_quarter_chord_sweep = get_wing_quarter_chord_sweep(
-      config->back_sweep_calculation_mode.value(), config->user_back_sweep_angle_at_quarter_chord.value(),
+      config->sweep_calculation_mode.value(), config->user_back_sweep_angle_at_quarter_chord.value(),
       data->specification_design_mach.value(), data->sizing_point_wing_loading.value(),
       data->specification_design_altitude.value(), back_max_thickness_to_chord_ratio,
       config->back_delta_drag_divergence_to_mach_design.value(), config->back_korn_technology_factor.value());
@@ -196,7 +193,7 @@ void Wing::standard_design() {
   std::map<std::string, std::tuple<double, double>> limits_field_length_and_wing_span =
       icao_easa_aerodrome_reference_code(data->specification_icao_aerodrome_reference_code.value());
 
-  /* Current limits of wing span according to reference code data*/
+  /* Current limits of wing span according to reference code data*/   //TODO kayra: lower limit neden var
   const auto [lower_wing_span, upper_wing_span] = limits_field_length_and_wing_span["wing_span"];
 
   /* Compute aspect ratios based on selected modes*/
@@ -206,19 +203,19 @@ void Wing::standard_design() {
   //myRuntimeInfo->out <<"WSR: "<< wing_span_ratio << std::endl;//TODO kayra:suan yanlislik yaratiyor
   //TODO kayra: bu design icin gereksiz bir parametre olarak kalabilir eger farkli yapmazsak.
 
-  myRuntimeInfo->out <<"AR_1 of forward wing: "<< std::endl;          
+  myRuntimeInfo->out <<"initial AR_1 of front wing: "<< std::endl;          
   double aspect_ratio = get_design_aspect_ratio(config->aspect_ratio_calculation_mode.value(),
                                                 config->user_aspect_ratio.value(), installed_quarter_chord_sweep);
                                                 
-  myRuntimeInfo->out <<"AR_2 of rear wing: "<< std::endl; 
-  double back_aspect_ratio = get_design_aspect_ratio(config->back_aspect_ratio_calculation_mode.value(),
+  myRuntimeInfo->out <<"initial AR_2 of back wing: "<< std::endl; 
+  double back_aspect_ratio = get_design_aspect_ratio(config->aspect_ratio_calculation_mode.value(),
                                                 config->user_back_aspect_ratio.value(), installed_back_quarter_chord_sweep);
   
   /* Determine spans and possible aspect ratio limits */
   double span = sqrt(aspect_ratio * installed_front_wing_area);
-  myRuntimeInfo->out <<"b_1: "<< span << std::endl;
+  myRuntimeInfo->out <<"initial calculated b_1 for front wing: "<< span << std::endl;
   double back_span = sqrt(back_aspect_ratio * installed_back_wing_area);
-  myRuntimeInfo->out <<"b_2: "<< back_span << std::endl;
+  myRuntimeInfo->out <<"initial calculated b_2 for back wing: "<< back_span << std::endl;
 
 
   double maximum_possible_aspect_ratio = pow(upper_wing_span, 2.) / installed_front_wing_area;
@@ -229,13 +226,13 @@ void Wing::standard_design() {
   /* if current aspect ratios exceed maximum possible aspect ratio -> set aspect ratios to maximum possible aspect ratio
    * and set span to according wing span */
   if (aspect_ratio > maximum_possible_aspect_ratio) {
-    myRuntimeInfo->warn << "Calculated aspect ratio for forward wing > maximum_possible_aspect_ratio" << std::endl;
+    myRuntimeInfo->warn << "Initial aspect ratio for front wing > maximum_possible_aspect_ratio for front wing" << std::endl;
     myRuntimeInfo->warn << "Switch to maximum possible aspect_ratio" << std::endl;
     aspect_ratio = maximum_possible_aspect_ratio;
     span = upper_wing_span;
   }
   if (back_aspect_ratio > maximum_possible_back_aspect_ratio) {
-    myRuntimeInfo->warn << "Calculated aspect ratio for rear wing > maximum_possible_back_aspect_ratio" << std::endl;
+    myRuntimeInfo->warn << "Initial aspect ratio for back wing > maximum_possible_back_aspect_ratio for back wing" << std::endl;
     myRuntimeInfo->warn << "Switch to maximum possible back_aspect_ratio" << std::endl;
     back_aspect_ratio = maximum_possible_back_aspect_ratio;
     back_span = upper_wing_span;
@@ -243,13 +240,13 @@ void Wing::standard_design() {
   // Check if current aspect ratios are below lower span limits -> set aspect ratios to minimum possible aspect ratios and
   // set spans to according wing spans */
   if (aspect_ratio < minimum_possible_aspect_ratio) {
-    myRuntimeInfo->warn << "Current aspect ratio for forward wing breach lower span limit ..." << std::endl;
+    myRuntimeInfo->warn << "Current aspect ratio for front wing breach lower span limit ..." << std::endl;
     myRuntimeInfo->warn << "Switch to minimum possible aspect_ratio" << std::endl;
     aspect_ratio = minimum_possible_aspect_ratio;
     span = lower_wing_span;
   }
   if (back_aspect_ratio < minimum_possible_back_aspect_ratio) {
-    myRuntimeInfo->warn << "Current aspect ratio for rear wing breach lower span limit ..." << std::endl;
+    myRuntimeInfo->warn << "Current aspect ratio for back wing breach lower span limit ..." << std::endl;
     myRuntimeInfo->warn << "Switch to minimum possible back_aspect_ratio" << std::endl;
     back_aspect_ratio = minimum_possible_back_aspect_ratio;
     back_span = lower_wing_span;
@@ -261,17 +258,22 @@ void Wing::standard_design() {
   const double installed_span = span;
   const double installed_back_span = back_span;
   
-  myRuntimeInfo->out <<"jetzt korrigierte Werte: "<< std::endl;
+  myRuntimeInfo->out <<"Corrected final values for aspect ratios and wing spans according to the airport limits: "<< std::endl;
   myRuntimeInfo->out <<"AR_1: "<< installed_aspect_ratio << std::endl;
   myRuntimeInfo->out <<"AR_2: "<< installed_back_aspect_ratio << std::endl;
   myRuntimeInfo->out <<"b_1: "<< installed_span << std::endl;
   myRuntimeInfo->out <<"b_2: "<< installed_back_span << std::endl;
 
-  double wing_span_ratio = installed_back_span / installed_span; //TODO kayra: just to show
-  myRuntimeInfo->out <<"Wing Span Ratio: "<< wing_span_ratio << std::endl;
+  const double installed_wing_span_ratio = installed_back_span / installed_span; //TODO kayra: just to show
+  myRuntimeInfo->out <<"Wing Span Ratio: "<< installed_wing_span_ratio << std::endl;
+
+  const double installed_aspect_ratio_ratio = installed_back_aspect_ratio / installed_aspect_ratio; //TODO kayra: just to show
+  myRuntimeInfo->out <<"Aspect Ratio Ratio: "<< installed_aspect_ratio_ratio << std::endl;
 
-  double aspect_ratio_ratio = installed_back_aspect_ratio / installed_aspect_ratio; //TODO kayra: just to show
-  myRuntimeInfo->out <<"Aspect Ratio Ratio: "<< aspect_ratio_ratio << std::endl;
+  if(installed_wing_span_ratio > 1.2 || installed_wing_span_ratio < 0.83 || 
+      installed_aspect_ratio_ratio  > 2|| installed_aspect_ratio_ratio < 0.83) {
+        myRuntimeInfo->warn << "Please consider that generally for a properly designed tandem wing configuration Wing Area Ratio, Wing Span Ratio and Aspect Ratio Ratio should be between 0.83 and 1.2" << std::endl;
+  }
 
   // Compute taper ratios based on selected modes
   myRuntimeInfo->out <<"Taper Ratio of forward wing: "<< std::endl;
@@ -280,18 +282,18 @@ void Wing::standard_design() {
                              installed_aspect_ratio, installed_quarter_chord_sweep);
   myRuntimeInfo->out <<"Taper Ratio of rear wing: "<< std::endl;
   const double installed_back_taper_ratio =
-      get_design_taper_ratio(config->back_taper_ratio_calculation_mode.value(), config->user_back_taper_ratio.value(),
+      get_design_taper_ratio(config->taper_ratio_calculation_mode.value(), config->user_back_taper_ratio.value(),
                              installed_back_aspect_ratio, installed_back_quarter_chord_sweep);
 
   // Comptue dihedrals based on selected modes 
-  //TODO kayra: iste burada rolling stability acisindan ciddi sikintilar cikabilir. dihedral diff gerekebilir. Paper var
+  //TODO kayra: paper var buna göre yap
   myRuntimeInfo->out <<"Dihedral Angle of forward wing: "<< std::endl;
   const double installed_dihedral = get_design_dihedral(
       config->dihedral_calculation_mode.value(), data->specification_wing_mounting.value(),
       config->user_dihedral.value(), installed_quarter_chord_sweep, data->specification_design_mach.value());
   myRuntimeInfo->out <<"Dihedral Angle of rear wing: "<< std::endl;
   const double installed_back_dihedral = get_design_dihedral(
-      config->back_dihedral_calculation_mode.value(), data->specification_back_wing_mounting.value(),
+      config->dihedral_calculation_mode.value(), data->specification_back_wing_mounting.value(),
       config->user_back_dihedral.value(), installed_back_quarter_chord_sweep, data->specification_design_mach.value());
 
   // Generate wings
@@ -472,7 +474,7 @@ geom2::Point_3 Wing::calculate_wing_position() {
 
   double y_mac_position = geom2::measure::mean_aerodynamic_chord_position(data->wing);
   double x_mac_position_quarter_chord_referred_to_centerline =
-      geom2::measure::offset_LE(data->wing, y_mac_position).x() +
+      geom2::measure::offset_LE(data->wing, y_mac_position).x() +   //TODO kayra: if twist this should be adjusted
       0.25 * geom2::measure::chord(data->wing, y_mac_position);
 
   const double fuselage_length = geom2::measure::length(data->fuselage);
@@ -515,7 +517,7 @@ geom2::Point_3 Wing::calculate_back_wing_position() {
 
   double y_mac_position = geom2::measure::mean_aerodynamic_chord_position(data->back_wing);
   double x_mac_position_quarter_chord_referred_to_centerline =
-      geom2::measure::offset_LE(data->back_wing, y_mac_position).x() +
+      geom2::measure::offset_LE(data->back_wing, y_mac_position).x() + //TODO kayra: if twist this should be adjusted
       0.25 * geom2::measure::chord(data->back_wing, y_mac_position);
 
   const double fuselage_length = geom2::measure::length(data->fuselage);
@@ -1005,8 +1007,7 @@ geom2::MultisectionSurface<geom2::AirfoilSection> Wing::calculate_unkinked_back_
 
 
 void Wing::flops_mass() {
-  myRuntimeInfo->out << "Front Wing mass estimation using ... FLOPS" << std::endl;
-  myRuntimeInfo->out << "Calculating ... " << std::endl;
+
   auto get_max_thickness_to_chord_ratio = [](const geom2::AirfoilSection section) {
     std::vector<double> t_to_c;
     const double chord = section.get_chord_length();
@@ -1018,10 +1019,7 @@ void Wing::flops_mass() {
     return *std::max_element(t_to_c.begin(), t_to_c.end());
   };
 
-  myRuntimeInfo->warn << "Following values are Hard coded!" << std::endl;
-  myRuntimeInfo->warn << "================================" << std::endl;
-  myRuntimeInfo->warn << "\twing_mounted_engines" << std::endl;
-  myRuntimeInfo->warn << "================================" << std::endl;
+  // Unit conversions and reference values for front and back wing
   const double to_foot = convertUnit(METER, FOOT, 1.0);
   const double to_meter = 1. / to_foot;
   const double taper_ratio = geom2::measure::taper_ratio(data->wing);
@@ -1036,32 +1034,65 @@ void Wing::flops_mass() {
   const double to_lb = convertUnit(KILO, GRAM, NOPREFIX, POUND, 1.0);
   const double to_kg = 1 / to_lb;
 
-  const double dg = data->mtom_mass_properties.data["mass"].value() * to_lb;
+  const double back_taper_ratio = geom2::measure::taper_ratio(data->back_wing);
+  const double back_span_ft = geom2::measure::span(data->back_wing) * to_foot;
+  const double back_wing_area_sqrft = geom2::measure::reference_area(data->back_wing) * to_foot * to_foot;
+  const double back_aspect_ratio = geom2::measure::aspect_ratio(data->back_wing);
+  const double back_t_to_c_root = get_max_thickness_to_chord_ratio(data->back_wing.sections.front());
+  const double back_t_to_c_tip = get_max_thickness_to_chord_ratio(data->back_wing.sections.back());
+  const double back_quarter_chord_sweep =
+      geom2::measure::sweep(data->back_wing, 0.25 * back_span_ft * to_meter, 0.25);  // measure at 50% half span -> 25% of span
+
+  
+  // It is assumed for preliminary design stage that each wing carries load proportinal to its wing area  //TODO kayra
+      
+  // PCTL for front wing
+  const double PCTL = wing_area_sqrft / (wing_area_sqrft + back_wing_area_sqrft); //TODO kayra:degismeli
+
+  // PCTL for back wing  
+  const double back_PCTL = back_wing_area_sqrft / (wing_area_sqrft + back_wing_area_sqrft);   //TODO kayra: degismeli
+
+  const double dg = data->mtom_mass_properties.data["mass"].value() * to_lb;  //* PCTL; //TODO kayra: bu isi de sor
 
-  // Tangent of 3/4 chord wing sweep
+  const double back_dg = data->mtom_mass_properties.data["mass"].value() * to_lb; //* back_PCTL; //TODO kayra: bu isi de sor
+
+
+
+  myRuntimeInfo->out << "Front Wing mass estimation using ... FLOPS" << std::endl;
+  myRuntimeInfo->out << "Calculating ... " << std::endl;
+
+  myRuntimeInfo->warn << "Following values are Hard coded!" << std::endl;
+  myRuntimeInfo->warn << "================================" << std::endl;
+  myRuntimeInfo->warn << "\twing_mounted_engines" << std::endl;
+  myRuntimeInfo->warn << "================================" << std::endl;
+
+
+
+
+  // Tangent of 3/4 chord wing sweep of front wing
   double TLAM = tan(quarter_chord_sweep) - 2.0 * (1.0 - taper_ratio) / (aspect_ratio * (1.0 + taper_ratio));
-  // Sine of 3/4 chord wing sweep angle
+  // Sine of 3/4 chord wing sweep angle of front wing
   double SLAM = TLAM / (sqrt(1 + TLAM * TLAM));
 
-  // Factor Eq. 15
+  // Factor Eq. 15 of front wing
   double C4 = 1.0 - 0.5 * config->flops_faert.value();
 
-  // Factor Eq. 16
+  // Factor Eq. 16 of front wing
   double C6 = 0.5 * config->flops_faert.value() - 0.16 * config->flops_fstrt.value();
 
-  // Aspect ratio defined factor Eq. 17
+  // Aspect ratio defined factor Eq. 17 of front wing
   double CAYA = aspect_ratio > 5. ? aspect_ratio - 5.0 : 0.;
-  // Wing sweep and aeroelastic factor
+  // Wing sweep and aeroelastic factor of front wing
   double CAYL = (1.0 - SLAM * SLAM) * (1.0 + C6 * SLAM * SLAM + 0.03 * CAYA * C4 * SLAM);
 
-  // Wing strut bracing factor (for cantilever w/o wing strut)
+  // Wing strut bracing factor of front wing (for cantilever w/o wing strut)
   double EMS = 1.0 - 0.25 * config->flops_fstrt.value();
-  // Weighted average of the wing thickness to chord ratio (Jenkinson 1999 formulation)
+  // Weighted average of the wing thickness to chord ratio of front wing (Jenkinson 1999 formulation)
   double TCA = 0.25 * (3. * t_to_c_tip + t_to_c_root);
-  // wing equivalent bending material factor Dimension [-]
+  // wing equivalent bending material factor of front wing Dimension [-]
   double BT = 0.215 * (0.37 + 0.7 * taper_ratio) * pow(span_ft * span_ft / wing_area_sqrft, EMS) / (CAYL * TCA);
 
-  // Total wing bending material weight
+  // Total wing bending material weight of front wing
 
   // A1 - Table 1 Flops wing weight constants for transport and hwb
   const double A1 = 8.80;
@@ -1079,33 +1110,32 @@ void Wing::flops_mass() {
   const double VARSWP = 0.0;
   // VFACT
   const double VFACT = 1.0;
-  // PCTL
-  const double PCTL = 1.0;
 
-  // Wing bending material weight, not including the effects of inertia [lb]
+
+  // Wing bending material weight of front wing, not including the effects of inertia [lb]
   const double W1NIR = A1 * BT * (1.0 + sqrt(A2 / span_ft)) * ULF * span_ft *
                        (1.0 - 0.4 * config->flops_fcomp.value()) * (1.0 - 0.1 * config->flops_faert.value()) * CAYF *
                        VFACT * PCTL / 1000000.0;
 
-  // Total wing shear material and control surface weight
+  // Total wing shear material and control surface weight of front wing
   const double A3 = 0.68;
   const double A4 = 0.34;
   const double A5 = 0.6;
 
-  // total movable wing surface area including flaps, elevators, spoilers etc [ft^2]
+  // total movable wing surface area of front wing including flaps, elevators, spoilers etc [ft^2]
 
   const double SFLAP = get_control_device_area_in_fraction_of_wing_area() *
                        wing_area_sqrft;  // set to 20% of wing_area - testing only TODO
 
   const double W2 = A3 * (1.0 - 0.17 * config->flops_fcomp.value()) * pow(SFLAP, A4) * pow(dg, A5);
 
-  // Total wing miscellaneous items weight
+  // Total wing miscellaneous items weight of front wing
   const double A6 = 0.035;
   const double A7 = 1.50;
 
   const double W3 = A6 * (1.0 - 0.3 * config->flops_fcomp.value()) * pow(wing_area_sqrft, A7);
 
-  // Wing bending material weight intertia relief adjustment
+  // Wing bending material weight intertia relief adjustment of front wing
   const int NEW = data->specification_number_of_wing_mounted_engines;
   const double CAYE = 1.0 - 0.03 * NEW;
 
@@ -1114,179 +1144,211 @@ void Wing::flops_mass() {
   // Total aft body weight for hybrid wing body aircraft - here no HWB
   const double W4 = 0.0;
 
-  // Total wing weight
+  // Total wing weight of front wing
   const double wing_mass = config->flops_technology_factor.value() * (W1 + W2 + W3 + W4) * to_kg;
 
-  myRuntimeInfo->out << "Calculated wing mass ... " << wing_mass << "[kg]" << std::endl;
-  myRuntimeInfo->out << "Wing bending material mass ... " << config->flops_technology_factor.value() * W1 * to_kg << "[kg]" << std::endl;
-  myRuntimeInfo->out << "Wing shear material and control surface mass ... " << config->flops_technology_factor.value() * W2 * to_kg << "[kg]" << std::endl;
-  myRuntimeInfo->out << "Wing miscellaneous items mass ... " << config->flops_technology_factor.value() * W3 * to_kg << "[kg]" << std::endl;
-  data->wing_mass_properties.data["mass"] = wing_mass;
-
-  // Calculate wing center of gravity according to Howe 2005 page 252
-  myRuntimeInfo->out << "Calculate Center of Gravity ..." << std::endl;
-  const double half_span_mac_position = geom2::measure::mean_aerodynamic_chord_position(data->wing);
-  const double mean_chord = geom2::measure::reference_area(data->wing) / geom2::measure::span(data->wing);
-  const double mac = geom2::measure::chord(data->wing, half_span_mac_position);
-  const geom2::Point_3 leading_edge_point = geom2::measure::offset_LE(data->wing, half_span_mac_position);
-  data->wing_mass_properties.data["x"].set_value(data->wing_position.x.value() + leading_edge_point.x() + 0.25 * mac +
-                                                 0.1 * mean_chord);
-  // Transform values into global coordinate system z-> -y z
-  data->wing_mass_properties.data["y"].set_value(data->wing_position.y.value() + fabs(leading_edge_point.z()) -
-                                                 fabs(leading_edge_point.z()));
-  data->wing_mass_properties.data["z"].set_value(data->wing_position.z.value() + leading_edge_point.y());
+  myRuntimeInfo->out << "Calculated front wing mass ... " << wing_mass << "[kg]" << std::endl;
+  myRuntimeInfo->out << "Wing bending material mass of front wing... " << config->flops_technology_factor.value() * W1 * to_kg << "[kg]" << std::endl;
+  myRuntimeInfo->out << "Wing shear material and control surface mass of front wing ... " << config->flops_technology_factor.value() * W2 * to_kg << "[kg]" << std::endl;
+  myRuntimeInfo->out << "Wing miscellaneous items mass of front wing ... " << config->flops_technology_factor.value() * W3 * to_kg << "[kg]" << std::endl;
 
-  myRuntimeInfo->out << "CG Position..." << std::endl;
-  myRuntimeInfo->out << "x ... " << data->wing_mass_properties.data["x"].value() << std::endl;
-  myRuntimeInfo->out << "y ... " << data->wing_mass_properties.data["y"].value() << std::endl;
-  myRuntimeInfo->out << "z ... " << data->wing_mass_properties.data["z"].value() << std::endl;
 
-}
 
-void Wing::back_flops_mass() {
   myRuntimeInfo->out << "Back Wing mass estimation using ... FLOPS" << std::endl;
   myRuntimeInfo->out << "Calculating ... " << std::endl;
-  auto get_max_thickness_to_chord_ratio = [](const geom2::AirfoilSection section) {
-    std::vector<double> t_to_c;
-    const double chord = section.get_chord_length();
-    const int steps = 50;
-    const double step_width = chord / steps;
-    for (size_t i = 0; i < steps; ++i) {
-      t_to_c.push_back(geom2::measure::thickness(section, i * step_width) / chord);
-    }
-    return *std::max_element(t_to_c.begin(), t_to_c.end());
-  };
 
   myRuntimeInfo->warn << "Following values are Hard coded!" << std::endl;
   myRuntimeInfo->warn << "================================" << std::endl;
   myRuntimeInfo->warn << "\twing_mounted_engines" << std::endl;
   myRuntimeInfo->warn << "================================" << std::endl;
-  const double to_foot = convertUnit(METER, FOOT, 1.0);
-  const double to_meter = 1. / to_foot;
-  const double taper_ratio = geom2::measure::taper_ratio(data->back_wing);
-  const double span_ft = geom2::measure::span(data->back_wing) * to_foot;
-  const double wing_area_sqrft = geom2::measure::reference_area(data->back_wing) * to_foot * to_foot;
-  const double aspect_ratio = geom2::measure::aspect_ratio(data->back_wing);
-  const double t_to_c_root = get_max_thickness_to_chord_ratio(data->back_wing.sections.front());
-  const double t_to_c_tip = get_max_thickness_to_chord_ratio(data->back_wing.sections.back());
-  const double quarter_chord_sweep =
-      geom2::measure::sweep(data->back_wing, 0.25 * span_ft * to_meter, 0.25);  // measure at 50% half span -> 25% of span
 
-  const double to_lb = convertUnit(KILO, GRAM, NOPREFIX, POUND, 1.0);
-  const double to_kg = 1 / to_lb;
 
-  const double dg = data->mtom_mass_properties.data["mass"].value() * to_lb;
+  // Tangent of 3/4 chord wing sweep of back wing
+  double back_TLAM = tan(back_quarter_chord_sweep) - 2.0 * (1.0 - back_taper_ratio) / (back_aspect_ratio * (1.0 + back_taper_ratio));
+  // Sine of 3/4 chord wing sweep angle of back wing
+  double back_SLAM = back_TLAM / (sqrt(1 + back_TLAM * back_TLAM));
 
-  // Tangent of 3/4 chord wing sweep
-  double TLAM = tan(quarter_chord_sweep) - 2.0 * (1.0 - taper_ratio) / (aspect_ratio * (1.0 + taper_ratio));
-  // Sine of 3/4 chord wing sweep angle
-  double SLAM = TLAM / (sqrt(1 + TLAM * TLAM));
+  // Factor Eq. 15 of back wing
+  double back_C4 = 1.0 - 0.5 * config->back_flops_faert.value();
 
-  // Factor Eq. 15
-  double C4 = 1.0 - 0.5 * config->back_flops_faert.value();
+  // Factor Eq. 16 of back wing
+  double back_C6 = 0.5 * config->back_flops_faert.value() - 0.16 * config->back_flops_fstrt.value();
 
-  // Factor Eq. 16
-  double C6 = 0.5 * config->back_flops_faert.value() - 0.16 * config->back_flops_fstrt.value();
+  // Aspect ratio defined factor Eq. 17 of back wing
+  double back_CAYA = back_aspect_ratio > 5. ? back_aspect_ratio - 5.0 : 0.;
+  // Wing sweep and aeroelastic factor of back wing
+  double back_CAYL = (1.0 - back_SLAM * back_SLAM) * (1.0 + back_C6 * back_SLAM * back_SLAM + 0.03 * back_CAYA * back_C4 * back_SLAM);
 
-  // Aspect ratio defined factor Eq. 17
-  double CAYA = aspect_ratio > 5. ? aspect_ratio - 5.0 : 0.;
-  // Wing sweep and aeroelastic factor
-  double CAYL = (1.0 - SLAM * SLAM) * (1.0 + C6 * SLAM * SLAM + 0.03 * CAYA * C4 * SLAM);
-
-  // Wing strut bracing factor (for cantilever w/o wing strut)
-  double EMS = 1.0 - 0.25 * config->back_flops_fstrt.value();
-  // Weighted average of the wing thickness to chord ratio (Jenkinson 1999 formulation)
-  double TCA = 0.25 * (3. * t_to_c_tip + t_to_c_root);
-  // wing equivalent bending material factor Dimension [-]
-  double BT = 0.215 * (0.37 + 0.7 * taper_ratio) * pow(span_ft * span_ft / wing_area_sqrft, EMS) / (CAYL * TCA);
-
-  // Total wing bending material weight
+  // Wing strut bracing factor of back wing(for cantilever w/o wing strut)
+  double back_EMS = 1.0 - 0.25 * config->back_flops_fstrt.value();
+  // Weighted average of the wing thickness to chord ratio of back wing(Jenkinson 1999 formulation)
+  double back_TCA = 0.25 * (3. * back_t_to_c_tip + back_t_to_c_root);
+  // wing equivalent bending material factor of back wing Dimension [-]
+  double back_BT = 0.215 * (0.37 + 0.7 * back_taper_ratio) * pow(back_span_ft * back_span_ft / back_wing_area_sqrft, back_EMS) / (back_CAYL * back_TCA);
 
+  // Total wing bending material weight of back wing
+  //TODO kayra: I keep these variables with different names altough they are same with front wing values, in order to keep maintenance easier if some of them should change in the future.
   // A1 - Table 1 Flops wing weight constants for transport and hwb
-  const double A1 = 8.80;
+  const double back_A1 = 8.80;
 
   // A2
-  const double A2 = 6.25;
+  const double back_A2 = 6.25;
 
   // ULF - Structural Ultimate load factor - default 3.75
-  const double ULF = 3.75;
+  const double back_ULF = 3.75;
 
   // CAYF - only single fuselage usage
-  const double CAYF = 1.0;
+  const double back_CAYF = 1.0;
 
   // VARSWP - no variable sweep allowed
-  const double VARSWP = 0.0;
+  const double back_VARSWP = 0.0;
   // VFACT
-  const double VFACT = 1.0;
-  // PCTL
-  const double PCTL = 1.0;
+  const double back_VFACT = 1.0;
 
-  // Wing bending material weight, not including the effects of inertia [lb]
-  const double W1NIR = A1 * BT * (1.0 + sqrt(A2 / span_ft)) * ULF * span_ft *
-                       (1.0 - 0.4 * config->back_flops_fcomp.value()) * (1.0 - 0.1 * config->back_flops_faert.value()) * CAYF *
-                       VFACT * PCTL / 1000000.0;
 
-  // Total wing shear material and control surface weight
-  const double A3 = 0.68;
-  const double A4 = 0.34;
-  const double A5 = 0.6;
+  // Wing bending material weight of back wing, not including the effects of inertia [lb]
+  const double back_W1NIR = back_A1 * back_BT * (1.0 + sqrt(back_A2 / back_span_ft)) * back_ULF * back_span_ft *
+                       (1.0 - 0.4 * config->back_flops_fcomp.value()) * (1.0 - 0.1 * config->back_flops_faert.value()) * back_CAYF *
+                       back_VFACT * back_PCTL / 1000000.0;
 
-  // total movable wing surface area including flaps, elevators, spoilers etc [ft^2]
+  // Total wing shear material and control surface weight of back wing
+  const double back_A3 = 0.68;
+  const double back_A4 = 0.34;
+  const double back_A5 = 0.6;
 
-  const double SFLAP = get_back_control_device_area_in_fraction_of_back_wing_area() *
-                       wing_area_sqrft;  // set to 20% of wing_area - testing only TODO
+  // total movable wing surface area of back wing including flaps, elevators, spoilers etc [ft^2]
 
-  const double W2 = A3 * (1.0 - 0.17 * config->back_flops_fcomp.value()) * pow(SFLAP, A4) * pow(dg, A5);
+  const double back_SFLAP = get_back_control_device_area_in_fraction_of_back_wing_area() *
+                       back_wing_area_sqrft;  // set to 20% of wing_area - testing only TODO
 
-  // Total wing miscellaneous items weight
-  const double A6 = 0.035;
-  const double A7 = 1.50;
+  const double back_W2 = back_A3 * (1.0 - 0.17 * config->back_flops_fcomp.value()) * pow(back_SFLAP, back_A4) * pow(back_dg, back_A5);
 
-  const double W3 = A6 * (1.0 - 0.3 * config->back_flops_fcomp.value()) * pow(wing_area_sqrft, A7);
+  // Total wing miscellaneous items weight of back wing
+  const double back_A6 = 0.035;
+  const double back_A7 = 1.50;
 
-  // Wing bending material weight intertia relief adjustment
-  const int NEW = data->specification_number_of_wing_mounted_engines;
-  const double CAYE = 1.0 - 0.03 * NEW;
+  const double back_W3 = back_A6 * (1.0 - 0.3 * config->back_flops_fcomp.value()) * pow(back_wing_area_sqrft, back_A7);
 
-  const double W1 = (dg * CAYE * W1NIR + W2 + W3) / (1.0 + W1NIR) - W2 - W3;
+  // Wing bending material weight intertia relief adjustment of back wing
+  const int back_NEW = data->specification_number_of_wing_mounted_engines;
+  const double back_CAYE = 1.0 - 0.03 * back_NEW;
+
+  const double back_W1 = (back_dg * back_CAYE * back_W1NIR + back_W2 + back_W3) / (1.0 + back_W1NIR) - back_W2 - back_W3;
 
   // Total aft body weight for hybrid wing body aircraft - here no HWB
-  const double W4 = 0.0;
+  const double back_W4 = 0.0;
+
+  // Total wing weight of back wing
+  const double back_wing_mass = config->back_flops_technology_factor.value() * (back_W1 + back_W2 + back_W3 + back_W4) * to_kg;
+
+  myRuntimeInfo->out << "Calculated back wing mass ... " << back_wing_mass << "[kg]" << std::endl;
+  myRuntimeInfo->out << "Wing bending material mass of back wing ... " << config->back_flops_technology_factor.value() * back_W1 * to_kg << "[kg]" << std::endl;
+  myRuntimeInfo->out << "Wing shear material and control surface mass of back wing ... " << config->back_flops_technology_factor.value() * back_W2 * to_kg << "[kg]" << std::endl;
+  myRuntimeInfo->out << "Wing miscellaneous items mass of back wing ... " << config->back_flops_technology_factor.value() * back_W3 * to_kg << "[kg]" << std::endl;
+
+  myRuntimeInfo->out << "Overall Tandem Wing Set mass estimation by adding both wings masses calculated with FLOPS method ..." << std::endl;
+  myRuntimeInfo->out << "Calculating ... " << std::endl;
+
+  const double tandem_wing_set_mass = wing_mass + back_wing_mass;
+
+  myRuntimeInfo->out << "Calculated overall tandem wing set mass ... " << tandem_wing_set_mass << "[kg]" << std::endl;
+
+  data->wing_mass_properties.data["mass"] = tandem_wing_set_mass;
+  //data->back_wing_mass_properties.data["mass"] = back_wing_mass;  //TODO kayra: olmamasi lazim
+
+
+
+  // Calculate wing center of gravity of front wing according to Howe 2005 page 252
+  myRuntimeInfo->out << "Calculate Center of Gravity of front wing ..." << std::endl;
+  const double half_span_mac_position = geom2::measure::mean_aerodynamic_chord_position(data->wing);
+  const double mean_chord = geom2::measure::reference_area(data->wing) / geom2::measure::span(data->wing);
+  const double mac = geom2::measure::chord(data->wing, half_span_mac_position);
+  const geom2::Point_3 leading_edge_point = geom2::measure::offset_LE(data->wing, half_span_mac_position); //TODO kayra: if twist this should be adjusted
+  const double CoG_x = data->wing_position.x.value() + leading_edge_point.x() + 0.25 * mac +
+                                                 0.1 * mean_chord;
+  // Transform values into global coordinate system z-> -y z
+  const double CoG_y = data->wing_position.y.value() + fabs(leading_edge_point.z()) -
+                                                 fabs(leading_edge_point.z());
+  const double CoG_z = data->wing_position.z.value() + leading_edge_point.y();
+
+  myRuntimeInfo->out << "CG Position of front wing..." << std::endl;
+  myRuntimeInfo->out << "x ... " << CoG_x << std::endl;
+  myRuntimeInfo->out << "y ... " << CoG_y << std::endl;
+  myRuntimeInfo->out << "z ... " << CoG_z << std::endl;
+
+
+
+
+  // Calculate wing center of gravity of back wing according to Howe 2005 page 252
+  myRuntimeInfo->out << "Calculate Center of Gravity of back wing..." << std::endl;
+  const double back_half_span_mac_position = geom2::measure::mean_aerodynamic_chord_position(data->back_wing);
+  const double back_mean_chord = geom2::measure::reference_area(data->back_wing) / geom2::measure::span(data->back_wing);
+  const double back_mac = geom2::measure::chord(data->back_wing, back_half_span_mac_position);
+  const geom2::Point_3 back_leading_edge_point = geom2::measure::offset_LE(data->back_wing, back_half_span_mac_position); //TODO kayra: if twist this should be adjusted
+  const double back_CoG_x = data->wing_position.x.value() + data->back_wing.origin.x() + back_leading_edge_point.x() + 0.25 * back_mac +
+                                                 0.1 * back_mean_chord;
+  // Transform values into global coordinate system z-> -y z
+  const double back_CoG_y = data->wing_position.y.value() + data->back_wing.origin.y() + fabs(back_leading_edge_point.z()) -
+                                                 fabs(back_leading_edge_point.z());
+  const double back_CoG_z = data->wing_position.z.value() + data->back_wing.origin.z() + back_leading_edge_point.y();
+
+  myRuntimeInfo->out << "CG Position of back wing..." << std::endl;
+  myRuntimeInfo->out << "x ... " << back_CoG_x << std::endl;
+  myRuntimeInfo->out << "y ... " << back_CoG_y << std::endl;
+  myRuntimeInfo->out << "z ... " << back_CoG_z << std::endl;
 
-  // Total wing weight
-  const double wing_mass = config->back_flops_technology_factor.value() * (W1 + W2 + W3 + W4) * to_kg;
-
-  myRuntimeInfo->out << "Calculated wing mass ... " << wing_mass << "[kg]" << std::endl;
-  myRuntimeInfo->out << "Wing bending material mass ... " << config->back_flops_technology_factor.value() * W1 * to_kg << "[kg]" << std::endl;
-  myRuntimeInfo->out << "Wing shear material and control surface mass ... " << config->back_flops_technology_factor.value() * W2 * to_kg << "[kg]" << std::endl;
-  myRuntimeInfo->out << "Wing miscellaneous items mass ... " << config->back_flops_technology_factor.value() * W3 * to_kg << "[kg]" << std::endl;
-  data->back_wing_mass_properties.data["mass"] = wing_mass;
-
-  // Calculate wing center of gravity according to Howe 2005 page 252
-  myRuntimeInfo->out << "Calculate Center of Gravity ..." << std::endl;
-  const double half_span_mac_position = geom2::measure::mean_aerodynamic_chord_position(data->back_wing);
-  const double mean_chord = geom2::measure::reference_area(data->back_wing) / geom2::measure::span(data->back_wing);
-  const double mac = geom2::measure::chord(data->back_wing, half_span_mac_position);
-  const geom2::Point_3 leading_edge_point = geom2::measure::offset_LE(data->back_wing, half_span_mac_position);
-  data->back_wing_mass_properties.data["x"].set_value(data->back_wing.origin.x() + leading_edge_point.x() + 0.25 * mac +
+ 
+  // Calculate center of gravity of tandem wing set using mass weighted CoG formula //TODO kayra
+  myRuntimeInfo->out << "Calculate Center of Gravity of overall tandem wing set ..." << std::endl;
+  const double wingset_CoG_x = (CoG_x * wing_mass + back_CoG_x * back_wing_mass) / tandem_wing_set_mass;
+  const double wingset_CoG_y = (CoG_y * wing_mass + back_CoG_y * back_wing_mass) / tandem_wing_set_mass;
+  const double wingset_CoG_z = (CoG_z * wing_mass + back_CoG_z * back_wing_mass) / tandem_wing_set_mass;
+
+  
+  data->wing_mass_properties.data["x"].set_value(wingset_CoG_x);
+  data->wing_mass_properties.data["y"].set_value(wingset_CoG_y);
+  data->wing_mass_properties.data["z"].set_value(wingset_CoG_z);
+
+  myRuntimeInfo->out << "CG Position of Tandem Wing Set ..." << std::endl;
+  myRuntimeInfo->out << "x ... " << data->wing_mass_properties.data["x"].value() << std::endl;
+  myRuntimeInfo->out << "y ... " << data->wing_mass_properties.data["y"].value() << std::endl;
+  myRuntimeInfo->out << "z ... " << data->wing_mass_properties.data["z"].value() << std::endl;
+
+
+  //TODO kayra: yedek olarak biraktim
+  /*data->wing_mass_properties.data["x"].set_value(data->wing_position.x.value() + leading_edge_point.x() + 0.25 * mac +
                                                  0.1 * mean_chord);
   // Transform values into global coordinate system z-> -y z
+  data->wing_mass_properties.data["y"].set_value(data->wing_position.y.value() + fabs(leading_edge_point.z()) -
+                                                 fabs(leading_edge_point.z()));
+  data->wing_mass_properties.data["z"].set_value(data->wing_position.z.value() + leading_edge_point.y());
+
+  myRuntimeInfo->out << "CG Position of front wing..." << std::endl;
+  myRuntimeInfo->out << "x ... " << data->wing_mass_properties.data["x"].value() << std::endl;
+  myRuntimeInfo->out << "y ... " << data->wing_mass_properties.data["y"].value() << std::endl;
+  myRuntimeInfo->out << "z ... " << data->wing_mass_properties.data["z"].value() << std::endl;
+
+
+  data->back_wing_mass_properties.data["x"].set_value(data->back_wing.origin.x() + back_leading_edge_point.x() + 0.25 * back_mac +
+                                                 0.1 * back_mean_chord);
+  // Transform values into global coordinate system z-> -y z
   data->back_wing_mass_properties.data["y"].set_value(data->back_wing.origin.y() + fabs(leading_edge_point.z()) -
                                                  fabs(leading_edge_point.z()));
   data->back_wing_mass_properties.data["z"].set_value(data->back_wing.origin.z() + leading_edge_point.y());
 
-  myRuntimeInfo->out << "CG Position..." << std::endl;
+  myRuntimeInfo->out << "CG Position of back wing..." << std::endl;
   myRuntimeInfo->out << "x ... " << data->back_wing_mass_properties.data["x"].value() << std::endl;
   myRuntimeInfo->out << "y ... " << data->back_wing_mass_properties.data["y"].value() << std::endl;
-  myRuntimeInfo->out << "z ... " << data->back_wing_mass_properties.data["z"].value() << std::endl;
+  myRuntimeInfo->out << "z ... " << data->back_wing_mass_properties.data["z"].value() << std::endl;*/
 
 }
 
+
 /**
  * \brief Chiozzotto weight estimation relationship mass estimation - assumes symmetric engine distribution on wings
  *
  */
-// TODO kayra: bu wing loadingli falan bir garip olabilir tandem icin, mtom da gariplik yaratir.
+// TODO kayra: bunda asla kink olmazsa sonuc cikmiyor sikintili bir yöntem
 void Wing::chiozzotto_wer_mass() {
   myRuntimeInfo->out << "Front Wing mass estimation using ... CHIOZZOTTO" << std::endl;
   myRuntimeInfo->out << "Calculating ... " << std::endl;
@@ -1337,19 +1399,26 @@ void Wing::chiozzotto_wer_mass() {
       {"AL ASW", {8.41, 3.017, 0.817}},   {"AL USW", {12.5, 2.177, 1.011}},   {"AL FSW", {25.4, 1.941, 1.204}},
       {"CFRP ASW", {15.3, 3.104, 0.909}}, {"CFRP USW", {13.6, 2.569, 0.910}}, {"CFRP FSW", {24.7, 2.117, 1.148}}};
 
+
+  const double ref_area = geom2::measure::reference_area(data->wing);             //TODO kayra: düzenleme
+  const double back_ref_area = geom2::measure::reference_area(data->back_wing);
+
+  const double area_portion = ref_area / (ref_area + back_ref_area);
+  const double back_area_portion = back_ref_area / (ref_area + back_ref_area);
+
   // Chiozzotto Engine relief factor k_e (5.5 - Page 116)
   // IF NUMBER OF ENGINES ON WING > 2 ->
   // double engine_relative_spanwise_position = 0.3;  // BOOKMARK
   // Vector with data for one half ->
   // More than one engine -> k_e = k_e_1 * k_e_2
 
-  double wing_loading = data->sizing_point_wing_loading.value();
+  double wing_loading = data->sizing_point_wing_loading.value();  //TODO kayra: düzenleme yok zaten kendini buluyor
   if (data->sizing_point_wing_loading.unit() == "kg/m^2") {
     wing_loading *= G_FORCE;
   }
 
   atmosphere isa;
-  const double mtom = data->mtom_mass_properties.data["mass"].value();
+  const double mtom = data->mtom_mass_properties.data["mass"].value() * area_portion;  //TODO kayra: düzenleme  : lift generation prop. to its area portion: first assumption
   const double aspect_ratio = geom2::measure::aspect_ratio(data->wing);
   const double taper_ratio = geom2::measure::taper_ratio(data->wing);
   const double half_span = geom2::measure::span(data->wing) * 0.5;
@@ -1375,10 +1444,10 @@ void Wing::chiozzotto_wer_mass() {
   } else {
     material += "USW";
   }
-  myRuntimeInfo->out << "Wing material used ... " << material << std::endl;
+  myRuntimeInfo->out << "Wing material used for front wing ... " << material << std::endl;
 
   // Get design relative kink -
-  double t_to_c_kink = 0.0;
+  double t_to_c_kink = 0.0;   //TODO kayra : düzeltilecek
   double eta_kink = 0.35;
   if (is_wing_kinked_check()) {
     eta_kink = get_design_relative_kink_position(
@@ -1413,27 +1482,27 @@ void Wing::chiozzotto_wer_mass() {
 
   // Check inputs
   if (!(chiozzotto_wer_check_limits(min_max_mtom.at(material), mtom))) {
-    myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... mtom ... " << mtom << std::endl;
+    myRuntimeInfo->warn << "Chiozzotto limit check for front wing out of bounds ... mtom front wing portion ... " << mtom << std::endl;
   }
   if (!(chiozzotto_wer_check_limits(min_max_w_to_s.at(material), wing_loading))) {
-    myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... wing_loading ... " << wing_loading << std::endl;
+    myRuntimeInfo->warn << "Chiozzotto limit check for front wing out of bounds ... front wing_loading ... " << wing_loading << std::endl;
   }
   if (!(chiozzotto_wer_check_limits(min_max_aspect_ratio.at(material), aspect_ratio))) {
-    myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... aspect ratio ... " << aspect_ratio << std::endl;
+    myRuntimeInfo->warn << "Chiozzotto limit check for front wing out of bounds ... front wings aspect ratio ... " << aspect_ratio << std::endl;
   }
   if (!(chiozzotto_wer_check_limits(min_max_taper_ratio.at(material), taper_ratio))) {
-    myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... taper ratio ... " << taper_ratio << std::endl;
+    myRuntimeInfo->warn << "Chiozzotto limit check for front wing out of bounds ... front wings taper ratio ... " << taper_ratio << std::endl;
   }
   if (!(chiozzotto_wer_check_limits(min_max_sweep.at(material), quarter_chord_sweep * geom2::detail::to_degrees))) {
-    myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... sweep ... "
+    myRuntimeInfo->warn << "Chiozzotto limit check for front wing out of bounds ... front wings sweep ... "
                         << quarter_chord_sweep * geom2::detail::to_degrees << std::endl;
   }
   if (!(chiozzotto_wer_check_limits(min_max_t_to_c.at(material), t_to_c_kink))) {
-    myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... thickness to chord ... " << t_to_c_kink
+    myRuntimeInfo->warn << "Chiozzotto limit check for front wing out of bounds ... front wings thickness to chord ... " << t_to_c_kink
                         << std::endl;
   }
   if (!(chiozzotto_wer_check_limits(min_max_v_eas.at(material), v_eas_design))) {
-    myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... equivalent airspeed at design point ... "
+    myRuntimeInfo->warn << "Chiozzotto limit check for front wing out of bounds ... equivalent airspeed at design point ... "
                         << v_eas_design << std::endl;
   }
 
@@ -1468,238 +1537,220 @@ void Wing::chiozzotto_wer_mass() {
 
   double wing_mass = mass_covers + mass_spars_and_ribs + mass_secondary_structure;
 
-  myRuntimeInfo->out << "Calculated wing mass ... " << wing_mass << "[kg]" << std::endl;
-  myRuntimeInfo->out << "Wing mass spars and ribs ..." << mass_spars_and_ribs << "[kg]" << std::endl;
-  myRuntimeInfo->out << "Wing mass covers  ..." << mass_covers << "[kg]" << std::endl;
-  myRuntimeInfo->out << "Wing mass secondary structure ..." << mass_secondary_structure << "[kg]" << std::endl;
-  data->wing_mass_properties.data["mass"] = wing_mass;
+  myRuntimeInfo->out << "Calculated front wing mass ... " << wing_mass << "[kg]" << std::endl;
+  myRuntimeInfo->out << "Front Wing mass spars and ribs ..." << mass_spars_and_ribs << "[kg]" << std::endl;
+  myRuntimeInfo->out << "Front Wing mass covers  ..." << mass_covers << "[kg]" << std::endl;
+  myRuntimeInfo->out << "Front Wing mass secondary structure ..." << mass_secondary_structure << "[kg]" << std::endl;
 
-  // Calculate wing center of gravity according to Howe 2005 page 252
-  myRuntimeInfo->out << "Calculate Center of Gravity ..." << std::endl;
-  const double half_span_mac_position = geom2::measure::mean_aerodynamic_chord_position(data->wing);
-  const double mean_chord = geom2::measure::reference_area(data->wing) / geom2::measure::span(data->wing);
-  const double mac = geom2::measure::chord(data->wing, half_span_mac_position);
-  const geom2::Point_3 leading_edge_point = geom2::measure::offset_LE(data->wing, half_span_mac_position);
-  data->wing_mass_properties.data["x"].set_value(data->wing_position.x.value() + leading_edge_point.x() + 0.25 * mac +
-                                                 0.1 * mean_chord);
-  // Transform values into global coordinate system z-> -y z
-  data->wing_mass_properties.data["y"].set_value(data->wing_position.y.value() + fabs(leading_edge_point.z()) -
-                                                 fabs(leading_edge_point.z()));
-  data->wing_mass_properties.data["z"].set_value(data->wing_position.z.value() + leading_edge_point.y());
 
-  myRuntimeInfo->out << "CG Position..." << std::endl;
-  myRuntimeInfo->out << "x ... " << data->wing_mass_properties.data["x"].value() << std::endl;
-  myRuntimeInfo->out << "y ... " << data->wing_mass_properties.data["y"].value() << std::endl;
-  myRuntimeInfo->out << "z ... " << data->wing_mass_properties.data["z"].value() << std::endl;
-}
-
-void Wing::back_chiozzotto_wer_mass() {
   myRuntimeInfo->out << "Back Wing mass estimation using ... CHIOZZOTTO" << std::endl;
   myRuntimeInfo->out << "Calculating ... " << std::endl;
-
-
-  std::map<std::string, std::tuple<double, double>> min_max_mtom = {
-      {"AL ASW", {20000, 250000}},   {"AL USW", {20000, 250000}},   {"AL FSW", {20000, 250000}},
-      {"CFRP ASW", {20000, 250000}}, {"CFRP USW", {20000, 250000}}, {"CFRP FSW", {20000, 250000}}};
-  std::map<std::string, std::tuple<double, double>> min_max_w_to_s = {
-      {"AL ASW", {3000, 8000}},   {"AL USW", {3000, 8000}},   {"AL FSW", {3000, 8000}},
-      {"CFRP ASW", {3000, 8000}}, {"CFRP USW", {3000, 8000}}, {"CFRP FSW", {3000, 8000}}};
-  std::map<std::string, std::tuple<double, double>> min_max_aspect_ratio = {
-      {"AL ASW", {10, 20}},   {"AL USW", {10, 20}},   {"AL FSW", {8, 16}},
-      {"CFRP ASW", {10, 20}}, {"CFRP USW", {10, 20}}, {"CFRP FSW", {8, 16}}};
-  std::map<std::string, std::tuple<double, double>> min_max_sweep = {{"AL ASW", {10, 40}},  {"AL USW", {0, 10}},
-                                                                     {"AL FSW", {-5, -25}}, {"CFRP ASW", {10, 40}},
-                                                                     {"CFRP USW", {0, 10}}, {"CFRP FSW", {-5, -25}}};
-  std::map<std::string, std::tuple<double, double>> min_max_t_to_c = {
-      {"AL ASW", {0.08, 0.14}},   {"AL USW", {0.08, 0.18}},   {"AL FSW", {0.08, 0.16}},
-      {"CFRP ASW", {0.08, 0.14}}, {"CFRP USW", {0.08, 0.18}}, {"CFRP FSW", {0.08, 0.16}}};
-  std::map<std::string, std::tuple<double, double>> min_max_v_eas = {
-      {"AL ASW", {130, 200}},   {"AL USW", {130, 200}},   {"AL FSW", {130, 200}},
-      {"CFRP ASW", {130, 200}}, {"CFRP USW", {130, 200}}, {"CFRP FSW", {130, 200}}};
-  std::map<std::string, std::tuple<double, double>> min_max_taper_ratio = {
-      {"AL ASW", {0.1, 0.5}},   {"AL USW", {0.1, 0.5}},   {"AL FSW", {0.1, 0.5}},
-      {"CFRP ASW", {0.1, 0.5}}, {"CFRP USW", {0.1, 0.5}}, {"CFRP FSW", {0.1, 0.5}}};
-
-  /**< Parameter containing a table of the constants and exponents of conventional Aircraft Covers*/
-  std::map<std::string, std::tuple<double, double, double, double, double, double, double, double>>
-      lookup_table_covers = {{"AL ASW", {2.14e-3, 1.345, -0.623, 1.286, -1.601, -0.830, -0.010, 0.104}},
-                             {"AL USW", {2.66e-3, 1.320, -0.731, 1.370, 6.733, -0.945, 0.110, 0.305}},
-                             {"AL FSW", {2.82e-5, 1.373, -1.152, 2.006, -5.507, -1.549, 1.036, 0.698}},
-                             {"CFRP ASW", {9.47e-4, 1.440, -0.699, 1.233, -1.044, -0.808, 0.017, 0.518}},
-                             {"CFRP USW", {1.29e-3, 1.396, -0.929, 1.4481, 0.667, -1.051, 0.266, 0.633}},
-                             {"CFRP FSW", {3.92e-5, 1.401, -1.111, 1.874, -4.084, -1.436, 0.889, 0.751}}};
-
-  /**< Parameter containing a table of constants and exponents of conventional Aircraft spars and ribs*/
-  std::map<std::string, std::tuple<double, double, double, double, double, double, double, double>>
-      lookup_table_spars_and_ribs = {{"AL ASW", {1.86, 1.384, -0.972, -0.072, -0.072, 0.486, 0.043, -0.126}},
-                                     {"AL USW", {2.99, 1.376, -1.090, 0.008, 5.376, 0.423, 0.112, -0.070}},
-                                     {"AL FSW", {9.82e-2, 1.362, -1.319, 0.598, -3.711, -0.126, 0.673, 0.391}},
-                                     {"CFRP ASW", {1.50e-1, 1.411, -0.832, 0.105, -0.235, 0.301, 0.057, 0.101}},
-                                     {"CFRP USW", {2.85e-1, 1.391, -0.994, 0.210, 5.854, 0.194, 0.164, 0.189}},
-                                     {"CFRP FSW", {2.31e-2, 1.370, -1.174, 0.650, -2.750, -0.204, 0.620, 0.454}}};
-
-  /**< Parameter containing a table of constants  and exponents of conventional Aircraft engine relief*/
-  std::map<std::string, std::tuple<double, double, double>> lookup_table_engine_relief = {
-      {"AL ASW", {8.41, 3.017, 0.817}},   {"AL USW", {12.5, 2.177, 1.011}},   {"AL FSW", {25.4, 1.941, 1.204}},
-      {"CFRP ASW", {15.3, 3.104, 0.909}}, {"CFRP USW", {13.6, 2.569, 0.910}}, {"CFRP FSW", {24.7, 2.117, 1.148}}};
-
+  
+  
   // Chiozzotto Engine relief factor k_e (5.5 - Page 116)
   // IF NUMBER OF ENGINES ON WING > 2 ->
   // double engine_relative_spanwise_position = 0.3;  // BOOKMARK
   // Vector with data for one half ->
   // More than one engine -> k_e = k_e_1 * k_e_2
 
-  double wing_loading = data->sizing_point_wing_loading.value();
+  double back_wing_loading = data->sizing_point_wing_loading.value();    //TODO kayra: düzenleme
   if (data->sizing_point_wing_loading.unit() == "kg/m^2") {
-    wing_loading *= G_FORCE;
+    back_wing_loading *= G_FORCE;
   }
 
-  atmosphere isa;
-  const double mtom = data->mtom_mass_properties.data["mass"].value();
-  const double aspect_ratio = geom2::measure::aspect_ratio(data->back_wing);
-  const double taper_ratio = geom2::measure::taper_ratio(data->back_wing);
-  const double half_span = geom2::measure::span(data->back_wing) * 0.5;
-  const double speed_of_sound_at_zero = isa.getSpeedOfSoundISA(0.0);
-  const double pressure_at_zero = isa.getPressureISA(0.0);
-  const double pressure_at_design_altitude = isa.getPressureISA(data->specification_design_altitude.value());
+  const double back_mtom = data->mtom_mass_properties.data["mass"].value() * back_area_portion;    //TODO kayra: düzenleme
+  const double back_aspect_ratio = geom2::measure::aspect_ratio(data->back_wing);
+  const double back_taper_ratio = geom2::measure::taper_ratio(data->back_wing);
+  const double back_half_span = geom2::measure::span(data->back_wing) * 0.5;
 
-  std::vector<std::tuple<double, double>> engines_mass_and_eta = data->wing_mounted_engine_data;
+  std::vector<std::tuple<double, double>> back_engines_mass_and_eta = data->wing_mounted_engine_data;
 
-  // Calculation of EAS Design
-  const double v_eas_design = data->specification_design_mach.value() * speed_of_sound_at_zero *
-                              sqrt(pressure_at_design_altitude / pressure_at_zero);
+  // Calculation of EAS Design same with front wing
 
   // No LRA sweep existent - use 1/4 chord sweep (deviation from chiozzotto P.112 -> referred to figure 5.3)
-  const double quarter_chord_sweep = -geom2::measure::sweep(data->back_wing, 0.5 * half_span, 0.25);
+  const double back_quarter_chord_sweep = -geom2::measure::sweep(data->back_wing, 0.5 * back_half_span, 0.25);
 
   // Get material selection
-  std::string material = config->back_chiozzotto_material.value() + " ";
-  if (quarter_chord_sweep > ACCURACY_LOW) {
-    material += "ASW";
-  } else if (quarter_chord_sweep < -ACCURACY_LOW) {
-    material += "FSW";
+  std::string back_material = config->back_chiozzotto_material.value() + " ";
+  if (back_quarter_chord_sweep > ACCURACY_LOW) {
+    back_material += "ASW";
+  } else if (back_quarter_chord_sweep < -ACCURACY_LOW) {
+    back_material += "FSW";
   } else {
-    material += "USW";
+    back_material += "USW";
   }
-  myRuntimeInfo->out << "Wing material used ... " << material << std::endl;
+  myRuntimeInfo->out << "Wing material used for back wing ... " << back_material << std::endl;
 
   // Get design relative kink -
-  double t_to_c_kink = 0.0;
-  double eta_kink = 0.35;
+  double back_t_to_c_kink = 0.0;
+  double back_eta_kink = 0.35;
   if (is_back_wing_kinked_check()) {
-    eta_kink = get_design_relative_kink_position(
+    back_eta_kink = get_design_relative_kink_position(
         config->back_relative_kink_position_calculation_mode.value(), config->user_back_relative_kink.value(),
         config->track_based_back_initial_relative_kink.value(), data->track_based_back_relative_kink.value());
   }
 
   // Find thickness to chord ratio
-  const double chord_at_kink = geom2::measure::chord(data->back_wing, eta_kink * half_span);
+  const double back_chord_at_kink = geom2::measure::chord(data->back_wing, back_eta_kink * back_half_span);
 
-  double eta_position_last = 0.;
-  double max_thickness_last = 0.0;
-  double eta_position_next = 0.;
-  double max_thickness_next = 0.0;
+  double back_eta_position_last = 0.;
+  double back_max_thickness_last = 0.0;
+  double back_eta_position_next = 0.;
+  double back_max_thickness_next = 0.0;
   for (size_t idx = 0; idx < data->back_wing.sections.size() - 1; ++idx) {
-    eta_position_last = fabs(data->back_wing.sections.at(idx).origin.z()) / half_span;
-    eta_position_next = fabs(data->back_wing.sections.at(idx + 1).origin.z()) / half_span;
+    back_eta_position_last = fabs(data->back_wing.sections.at(idx).origin.z()) / back_half_span;
+    back_eta_position_next = fabs(data->back_wing.sections.at(idx + 1).origin.z()) / back_half_span;
     // If kink at a specified section
-    if (fabs(eta_position_last - eta_kink) < ACCURACY_MEDIUM) {
-      t_to_c_kink = geom2::measure::thickness_max(data->back_wing.sections.at(idx)) / chord_at_kink;
-    } else if (fabs(eta_position_next - eta_kink) < ACCURACY_MEDIUM) {
-      t_to_c_kink = geom2::measure::thickness_max(data->back_wing.sections.at(idx)) / chord_at_kink;
-    } else if (eta_position_last > eta_kink && eta_position_next < eta_kink) {
-      double max_thickness_last = geom2::measure::thickness_max(data->back_wing.sections.at(idx));
-      double max_thickness_next = geom2::measure::thickness_max(data->back_wing.sections.at(idx));
-      double dimensionless_position = (eta_kink - eta_position_last) / (eta_position_last - eta_position_next);
-      t_to_c_kink = std::lerp(max_thickness_last, max_thickness_next, dimensionless_position) / chord_at_kink;
+    if (fabs(back_eta_position_last - back_eta_kink) < ACCURACY_MEDIUM) {
+      back_t_to_c_kink = geom2::measure::thickness_max(data->back_wing.sections.at(idx)) / back_chord_at_kink;
+    } else if (fabs(back_eta_position_next - back_eta_kink) < ACCURACY_MEDIUM) {
+      back_t_to_c_kink = geom2::measure::thickness_max(data->back_wing.sections.at(idx)) / back_chord_at_kink;
+    } else if (back_eta_position_last > back_eta_kink && back_eta_position_next < back_eta_kink) {
+      double back_max_thickness_last = geom2::measure::thickness_max(data->back_wing.sections.at(idx));
+      double back_max_thickness_next = geom2::measure::thickness_max(data->back_wing.sections.at(idx));
+      double back_dimensionless_position = (back_eta_kink - back_eta_position_last) / (back_eta_position_last - back_eta_position_next);
+      back_t_to_c_kink = std::lerp(back_max_thickness_last, back_max_thickness_next, back_dimensionless_position) / back_chord_at_kink;
     } else {
       // Do nothing
     }
   }
 
   // Check inputs       //TODO kayra: bundan assagida bir seyi ikiyle carpip sonuc alabilirsin sanirim iki wing massi ile
-  if (!(chiozzotto_wer_check_limits(min_max_mtom.at(material), mtom))) {
-    myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... mtom ... " << mtom << std::endl;
+  if (!(chiozzotto_wer_check_limits(min_max_mtom.at(back_material), back_mtom))) {
+    myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... mtom ... " << back_mtom << std::endl;
   }
-  if (!(chiozzotto_wer_check_limits(min_max_w_to_s.at(material), wing_loading))) {
-    myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... wing_loading ... " << wing_loading << std::endl;
+  if (!(chiozzotto_wer_check_limits(min_max_w_to_s.at(back_material), back_wing_loading))) {
+    myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... wing_loading ... " << back_wing_loading << std::endl;
   }
-  if (!(chiozzotto_wer_check_limits(min_max_aspect_ratio.at(material), aspect_ratio))) {
-    myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... aspect ratio ... " << aspect_ratio << std::endl;
+  if (!(chiozzotto_wer_check_limits(min_max_aspect_ratio.at(back_material), back_aspect_ratio))) {
+    myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... aspect ratio ... " << back_aspect_ratio << std::endl;
   }
-  if (!(chiozzotto_wer_check_limits(min_max_taper_ratio.at(material), taper_ratio))) {
-    myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... taper ratio ... " << taper_ratio << std::endl;
+  if (!(chiozzotto_wer_check_limits(min_max_taper_ratio.at(back_material), back_taper_ratio))) {
+    myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... taper ratio ... " << back_taper_ratio << std::endl;
   }
-  if (!(chiozzotto_wer_check_limits(min_max_sweep.at(material), quarter_chord_sweep * geom2::detail::to_degrees))) {
+  if (!(chiozzotto_wer_check_limits(min_max_sweep.at(back_material), back_quarter_chord_sweep * geom2::detail::to_degrees))) {
     myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... sweep ... "
-                        << quarter_chord_sweep * geom2::detail::to_degrees << std::endl;
+                        << back_quarter_chord_sweep * geom2::detail::to_degrees << std::endl;
   }
-  if (!(chiozzotto_wer_check_limits(min_max_t_to_c.at(material), t_to_c_kink))) {
-    myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... thickness to chord ... " << t_to_c_kink
+  if (!(chiozzotto_wer_check_limits(min_max_t_to_c.at(back_material), back_t_to_c_kink))) {
+    myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... thickness to chord ... " << back_t_to_c_kink
                         << std::endl;
   }
-  if (!(chiozzotto_wer_check_limits(min_max_v_eas.at(material), v_eas_design))) {
+  if (!(chiozzotto_wer_check_limits(min_max_v_eas.at(back_material), v_eas_design))) {
     myRuntimeInfo->warn << "Chiozzotto limit check out of bounds ... equivalent airspeed at design point ... "
                         << v_eas_design << std::endl;
   }
 
-  const auto [C_er, E_eta_er, E_K_er] = lookup_table_engine_relief[material];
-
+  const auto [back_C_er, back_E_eta_er, back_E_K_er] = lookup_table_engine_relief[back_material];
+ 
+  
   // Engine relief factor
-  double k_e = 1.;
-  for (auto engine : engines_mass_and_eta) {
-    const auto [engine_mass, engine_spanwise_position] = engine;
+  double back_k_e = 1.;
+  for (auto engine : back_engines_mass_and_eta) {
+    const auto [back_engine_mass, back_engine_spanwise_position] = engine;
     /* If engines exist, the spanwise position is absolute, if engines not exist, the default values are in rel. half span */
-    double engine_relative_spanwise_position = data->engines_exist ? std::fabs(engine_spanwise_position)/half_span : engine_spanwise_position;
+    double back_engine_relative_spanwise_position = data->engines_exist ? std::fabs(back_engine_spanwise_position)/back_half_span : back_engine_spanwise_position;
     //TODO kayra: engine existin arka kanat icin farkli olmasi gerekebilir
-    k_e *= (1 - C_er * pow(engine_relative_spanwise_position, E_eta_er) * pow(engine_mass / mtom, E_K_er));
+    back_k_e *= (1 - back_C_er * pow(back_engine_relative_spanwise_position, back_E_eta_er) * pow(back_engine_mass / back_mtom, back_E_K_er));
   }
   // Secondary structure mass
-  double mass_secondary_structure = 0.0443 * mtom * config->back_chiozzotto_technology_factor.value();
+  double back_mass_secondary_structure = 0.0443 * back_mtom * config->back_chiozzotto_technology_factor.value();
 
   // Mass of covers
-  const auto [C_cov, mtom_cov, w_to_s_cov, ar_cov, cos_lambda_cov, t_to_c_kink_cov, v_eas_design_cov,
-              one_plus_taper_cov] = lookup_table_covers[material];
-  double mass_covers = (k_e * C_cov * pow(mtom, mtom_cov) * pow(wing_loading, w_to_s_cov) * pow(aspect_ratio, ar_cov) *
-                        pow(cos(quarter_chord_sweep), cos_lambda_cov) * pow(t_to_c_kink, t_to_c_kink_cov) *
-                        pow(v_eas_design, v_eas_design_cov) * pow(1 + taper_ratio, one_plus_taper_cov)) *
+  const auto [back_C_cov, back_mtom_cov, back_w_to_s_cov, back_ar_cov, back_cos_lambda_cov, back_t_to_c_kink_cov, back_v_eas_design_cov,
+              back_one_plus_taper_cov] = lookup_table_covers[back_material];
+  double back_mass_covers = (back_k_e * back_C_cov * pow(back_mtom, back_mtom_cov) * pow(back_wing_loading, back_w_to_s_cov) * pow(back_aspect_ratio, back_ar_cov) *
+                        pow(cos(back_quarter_chord_sweep), back_cos_lambda_cov) * pow(back_t_to_c_kink, back_t_to_c_kink_cov) *
+                        pow(v_eas_design, back_v_eas_design_cov) * pow(1 + back_taper_ratio, back_one_plus_taper_cov)) *
                        config->back_chiozzotto_technology_factor.value();
 
   // Mass of spars and ribs
-  const auto [C_sr, mtom_sr, w_to_s_sr, ar_sr, cos_lambda_sr, t_to_c_kink_sr, v_eas_design_sr, one_plus_taper_sr] =
-      lookup_table_spars_and_ribs[material];
-  double mass_spars_and_ribs = (C_sr * pow(mtom, mtom_sr) * pow(wing_loading, w_to_s_sr) * pow(aspect_ratio, ar_sr) *
-                                pow(cos(quarter_chord_sweep), cos_lambda_sr) * pow(t_to_c_kink, t_to_c_kink_sr) *
-                                pow(v_eas_design, v_eas_design_sr) * pow(1 + taper_ratio, one_plus_taper_sr)) *
+  const auto [back_C_sr, back_mtom_sr, back_w_to_s_sr, back_ar_sr, back_cos_lambda_sr, back_t_to_c_kink_sr, back_v_eas_design_sr, back_one_plus_taper_sr] =
+      lookup_table_spars_and_ribs[back_material];
+  double back_mass_spars_and_ribs = (back_C_sr * pow(back_mtom, back_mtom_sr) * pow(back_wing_loading, back_w_to_s_sr) * pow(back_aspect_ratio, back_ar_sr) *
+                                pow(cos(back_quarter_chord_sweep), back_cos_lambda_sr) * pow(back_t_to_c_kink, back_t_to_c_kink_sr) *
+                                pow(v_eas_design, back_v_eas_design_sr) * pow(1 + back_taper_ratio, back_one_plus_taper_sr)) *
                                config->back_chiozzotto_technology_factor.value();
 
-  double wing_mass = mass_covers + mass_spars_and_ribs + mass_secondary_structure;
+  double back_wing_mass = back_mass_covers + back_mass_spars_and_ribs + back_mass_secondary_structure;
 
-  myRuntimeInfo->out << "Calculated wing mass ... " << wing_mass << "[kg]" << std::endl;
-  myRuntimeInfo->out << "Wing mass spars and ribs ..." << mass_spars_and_ribs << "[kg]" << std::endl;
-  myRuntimeInfo->out << "Wing mass covers  ..." << mass_covers << "[kg]" << std::endl;
-  myRuntimeInfo->out << "Wing mass secondary structure ..." << mass_secondary_structure << "[kg]" << std::endl;
-  data->back_wing_mass_properties.data["mass"] = wing_mass;
-
-  // Calculate wing center of gravity according to Howe 2005 page 252
-  myRuntimeInfo->out << "Calculate Center of Gravity ..." << std::endl;
-  const double half_span_mac_position = geom2::measure::mean_aerodynamic_chord_position(data->back_wing);
-  const double mean_chord = geom2::measure::reference_area(data->back_wing) / geom2::measure::span(data->back_wing);
-  const double mac = geom2::measure::chord(data->back_wing, half_span_mac_position);
-  const geom2::Point_3 leading_edge_point = geom2::measure::offset_LE(data->back_wing, half_span_mac_position);
-  data->back_wing_mass_properties.data["x"].set_value(data->back_wing.origin.x() + leading_edge_point.x() + 0.25 * mac +
-                                                 0.1 * mean_chord);
+  myRuntimeInfo->out << "Calculated back wing mass ... " << back_wing_mass << "[kg]" << std::endl;
+  myRuntimeInfo->out << "Back Wing mass spars and ribs ..." << back_mass_spars_and_ribs << "[kg]" << std::endl;
+  myRuntimeInfo->out << "Back Wing mass covers  ..." << back_mass_covers << "[kg]" << std::endl;
+  myRuntimeInfo->out << "Back Wing mass secondary structure ..." << back_mass_secondary_structure << "[kg]" << std::endl;
+
+  
+  myRuntimeInfo->out << "Overall Tandem Wing Set mass estimation by adding both wings masses calculated with Chiozotto method ..." << std::endl;
+  myRuntimeInfo->out << "Calculating ... " << std::endl;
+
+  const double tandem_wing_set_mass = wing_mass + back_wing_mass;
+
+  myRuntimeInfo->out << "Calculated overall tandem wing set mass ... " << tandem_wing_set_mass << "[kg]" << std::endl;
+
+  data->wing_mass_properties.data["mass"] = tandem_wing_set_mass;
+  //data->back_wing_mass_properties.data["mass"] = back_wing_mass;  //TODO kayra: olmamasi lazim
+  
+  
+  // Calculate wing center of gravity of front wing according to Howe 2005 page 252
+  myRuntimeInfo->out << "Calculate Center of Gravity of front wing ..." << std::endl;
+  const double half_span_mac_position = geom2::measure::mean_aerodynamic_chord_position(data->wing);
+  const double mean_chord = geom2::measure::reference_area(data->wing) / geom2::measure::span(data->wing);
+  const double mac = geom2::measure::chord(data->wing, half_span_mac_position);
+  const geom2::Point_3 leading_edge_point = geom2::measure::offset_LE(data->wing, half_span_mac_position); //TODO kayra: if twist this should be adjusted
+  const double CoG_x = data->wing_position.x.value() + leading_edge_point.x() + 0.25 * mac +
+                                                 0.1 * mean_chord;
   // Transform values into global coordinate system z-> -y z
-  data->back_wing_mass_properties.data["y"].set_value(data->back_wing.origin.y() + fabs(leading_edge_point.z()) -
-                                                 fabs(leading_edge_point.z()));
-  data->back_wing_mass_properties.data["z"].set_value(data->back_wing.origin.z() + leading_edge_point.y());
+  const double CoG_y = data->wing_position.y.value() + fabs(leading_edge_point.z()) -
+                                                 fabs(leading_edge_point.z());
+  const double CoG_z = data->wing_position.z.value() + leading_edge_point.y();
 
-  myRuntimeInfo->out << "CG Position..." << std::endl;
-  myRuntimeInfo->out << "x ... " << data->back_wing_mass_properties.data["x"].value() << std::endl;
-  myRuntimeInfo->out << "y ... " << data->back_wing_mass_properties.data["y"].value() << std::endl;
-  myRuntimeInfo->out << "z ... " << data->back_wing_mass_properties.data["z"].value() << std::endl;
+  myRuntimeInfo->out << "CG Position of front wing..." << std::endl;
+  myRuntimeInfo->out << "x ... " << CoG_x << std::endl;
+  myRuntimeInfo->out << "y ... " << CoG_y << std::endl;
+  myRuntimeInfo->out << "z ... " << CoG_z << std::endl;
+
+
+
+
+  // Calculate wing center of gravity of back wing according to Howe 2005 page 252
+  myRuntimeInfo->out << "Calculate Center of Gravity of back wing..." << std::endl;
+  const double back_half_span_mac_position = geom2::measure::mean_aerodynamic_chord_position(data->back_wing);
+  const double back_mean_chord = geom2::measure::reference_area(data->back_wing) / geom2::measure::span(data->back_wing);
+  const double back_mac = geom2::measure::chord(data->back_wing, back_half_span_mac_position);
+  const geom2::Point_3 back_leading_edge_point = geom2::measure::offset_LE(data->back_wing, back_half_span_mac_position); //TODO kayra: if twist this should be adjusted
+  const double back_CoG_x = data->wing_position.x.value() + data->back_wing.origin.x() + back_leading_edge_point.x() + 0.25 * back_mac +
+                                                 0.1 * back_mean_chord;
+  // Transform values into global coordinate system z-> -y z
+  const double back_CoG_y = data->wing_position.y.value() + data->back_wing.origin.y() + fabs(back_leading_edge_point.z()) -
+                                                 fabs(back_leading_edge_point.z());
+  const double back_CoG_z = data->wing_position.z.value() + data->back_wing.origin.z() + back_leading_edge_point.y();
+
+  myRuntimeInfo->out << "CG Position of back wing..." << std::endl;
+  myRuntimeInfo->out << "x ... " << back_CoG_x << std::endl;
+  myRuntimeInfo->out << "y ... " << back_CoG_y << std::endl;
+  myRuntimeInfo->out << "z ... " << back_CoG_z << std::endl;
+
+ 
+  // Calculate center of gravity of tandem wing set using mass weighted CoG formula //TODO kayra
+  myRuntimeInfo->out << "Calculate Center of Gravity of overall tandem wing set ..." << std::endl;
+  const double wingset_CoG_x = (CoG_x * wing_mass + back_CoG_x * back_wing_mass) / tandem_wing_set_mass;
+  const double wingset_CoG_y = (CoG_y * wing_mass + back_CoG_y * back_wing_mass) / tandem_wing_set_mass;
+  const double wingset_CoG_z = (CoG_z * wing_mass + back_CoG_z * back_wing_mass) / tandem_wing_set_mass;
+
+  
+  data->wing_mass_properties.data["x"].set_value(wingset_CoG_x);
+  data->wing_mass_properties.data["y"].set_value(wingset_CoG_y);
+  data->wing_mass_properties.data["z"].set_value(wingset_CoG_z);
+
+  myRuntimeInfo->out << "CG Position of Tandem Wing Set ..." << std::endl;
+  myRuntimeInfo->out << "x ... " << data->wing_mass_properties.data["x"].value() << std::endl;
+  myRuntimeInfo->out << "y ... " << data->wing_mass_properties.data["y"].value() << std::endl;
+  myRuntimeInfo->out << "z ... " << data->wing_mass_properties.data["z"].value() << std::endl;
+
+    
 }
 
+
 bool Wing::chiozzotto_wer_check_limits(std::tuple<double, double>& min_max, double value_to_check) {
   const auto [min, max] = min_max;
   if ((value_to_check > min && value_to_check < max) ||
diff --git a/wing_design/src/tandem/cantilever/cantileverTandemWingDesign.h b/wing_design/src/tandem/cantilever/cantileverTandemWingDesign.h
index c725f6756e431ecbe0c058fd793974018db0647c..4c02725c76b9a3cb4ee4158e63c8922a92b6f3a5 100644
--- a/wing_design/src/tandem/cantilever/cantileverTandemWingDesign.h
+++ b/wing_design/src/tandem/cantilever/cantileverTandemWingDesign.h
@@ -65,30 +65,18 @@ class Wing : public Strategy {
   void standard_design();
 
   /**
-   * \brief mass method of front wing -> flops
+   * \brief mass method for front and back wing -> flops
    * 
    */
   void flops_mass();
 
-  /**
-   * \brief mass method of back wing -> flops
-   * 
-   */
-  void back_flops_mass();
-
 
   /**
-   * \brief mass method of front wing -> chiozzotto weight estimation relationship
+   * \brief mass method for front and back wing -> chiozzotto weight estimation relationship
    * 
    */
   void chiozzotto_wer_mass();
 
-  /**
-   * \brief mass method of back wing -> chiozzotto weight estimation relationship
-   * 
-   */
-  void back_chiozzotto_wer_mass();
-
 
   /**
    * \brief Chiozzotto utility check function
@@ -242,26 +230,47 @@ class Wing : public Strategy {
   void set_html_body(void);
 
   /**
-   * \brief Generate plot for wing planform
+   * \brief Generate plot for front wing planform
    * 
    * \return fs::path relative path between reportpage and wing planform plot
    */
   fs::path plot_wing_planform();
 
+    /**
+   * \brief Generate plot for back wing planform
+   * 
+   * \return fs::path relative path between reportpage and wing planform plot
+   */
+  fs::path plot_back_wing_planform();
+
   /**
-   * \brief Generate plots for thickness and twist 
+   * \brief Generate plots for thickness and twist of front wing
    * 
    * \return fs::path relative path between reportpage and thickness and twist distribution plot
    */
   fs::path plot_thickness_and_twist_distribution();
 
   /**
-   * \brief Generate airfoil plots
+   * \brief Generate plots for thickness and twist of back wing
+   * 
+   * \return fs::path relative path between reportpage and thickness and twist distribution plot
+   */
+  fs::path plot_back_thickness_and_twist_distribution();
+
+  /**
+   * \brief Generate airfoil plots for front wing
    * 
    * \return std::vector<fs::path> vector of relative paths between reportpage and thickness and twist distribution plot
    */
   std::vector<fs::path> plot_airfoils();
 
+  /**
+   * \brief Generate airfoil plots for back wing
+   * 
+   * \return std::vector<fs::path> vector of relative paths between reportpage and thickness and twist distribution plot
+   */
+  std::vector<fs::path> plot_back_airfoils();
+
   /**
    * \brief Method to combine two meshes into one
    * 
diff --git a/wing_design/src/tandem/cantilever/cantileverTandemWingDesignPlot.cpp b/wing_design/src/tandem/cantilever/cantileverTandemWingDesignPlot.cpp
index 63df5fba08b4e6446f5fdbd666589bc0baa02c3d..c37043e0990cb66a942e6a9bb488a36b2f5ffc1c 100644
--- a/wing_design/src/tandem/cantilever/cantileverTandemWingDesignPlot.cpp
+++ b/wing_design/src/tandem/cantilever/cantileverTandemWingDesignPlot.cpp
@@ -45,7 +45,7 @@ std::array<float,4UL> rgba_tandem(float r, float b, float g, float a) {
   return {1-a,r/255,b/255,g/255};
 }
 
-
+//TODO kayra: error : altough there is no kink. Calculated with kink
 /* */
 std::array<float,4UL> spar_color_tandem = rgba_tandem(99, 110, 114,1.0);
 
@@ -121,7 +121,7 @@ fs::path Wing::plot_wing_planform() {
   std::deque<double> chordwise;
   auto half_span = 0.5 * geom2::measure::span(data->wing);
   auto sections = data->wing.sections;
-  auto fuselage_max_half = geom2::measure::width_max(data->fuselage)*0.5;
+  auto fuselage_half = geom2::measure::width(data->fuselage , data->wing_position.x.value())*0.5;
 
   /* Loop coordinates from out to inside for LE and TE */
   for (auto it = sections.rbegin(); it != sections.rend(); ++it) {
@@ -162,8 +162,8 @@ fs::path Wing::plot_wing_planform() {
     legend_strings.push_back(device.name);
   }
   auto limits = ax->ylim();
-  ax->plot(std::vector<double>({fuselage_max_half,fuselage_max_half}),std::vector<double>({limits[0]-2,limits[1]+2}))->color({0, 1, 0,0}).line_style("--").line_width(1.0);
-  legend_strings.push_back("fuselage max width");
+  ax->plot(std::vector<double>({fuselage_half,fuselage_half}),std::vector<double>({limits[0]-2,limits[1]+2}))->color({0, 1, 0,0}).line_style("--").line_width(1.0);
+  legend_strings.push_back("fuselage width");
   matplot::legend_handle lh;
   ax->xlabel("y [m]");
   ax->ylabel("x [m]");
@@ -184,6 +184,80 @@ fs::path Wing::plot_wing_planform() {
   return fs::relative(plot_path, fs::path(rtIO->getHtmlDir()));
 }
 
+fs::path Wing::plot_back_wing_planform() {
+  /* Plot active */
+  if (!rtIO->plotOn) {
+    return "";
+  }
+
+  auto coordinates_xy = geom2::transform::outline_2d(data->back_wing, geom2::Direction_3(0, 0, 1));
+  std::deque<double> spanwise;
+  std::deque<double> chordwise;
+  auto half_span = 0.5 * geom2::measure::span(data->back_wing);
+  auto sections = data->back_wing.sections;
+  auto fuselage_half = geom2::measure::width(data->fuselage , (data->wing_position.x.value()+data->back_wing.origin.x()))*0.5;
+
+  /* Loop coordinates from out to inside for LE and TE */
+  for (auto it = sections.rbegin(); it != sections.rend(); ++it) {
+    spanwise.push_front(-it->origin.z());
+    spanwise.push_back(-it->origin.z());
+    chordwise.push_front(it->origin.x());
+    chordwise.push_back(it->origin.x()+it->get_chord_length());
+  }
+
+  /* Get mac coordinates */
+  auto mac_position_y = geom2::measure::mean_aerodynamic_chord_position(data->back_wing);
+  auto mac_position = geom2::measure::offset_LE(data->back_wing,mac_position_y);
+  auto mac = geom2::measure::chord(data->back_wing,mac_position_y);
+
+
+
+  /* Generate plot */
+  std::vector<std::string> legend_strings;
+  auto f = matplot::figure(true);
+  f->title("Wing planform");
+  auto ax = f->add_axes(false);
+  ax->plot(std::vector<double>(spanwise.begin(),spanwise.end()), std::vector<double>(chordwise.begin(),chordwise.end()))->color("black").line_width(1.5);
+  legend_strings.push_back("wing contour");
+  matplot::hold(true);
+  ax->plot({-mac_position.z(),-mac_position.z()},{mac_position.x(),mac_position.x()+mac})->color({0,128,128,128}).line_style(":").line_width(1.5);
+  legend_strings.push_back("MAC");
+
+  /* Plot spars */
+  for(auto& spar : data->back_spars) {
+    const auto [x,y,_] = gen_polygon_tandem(data->back_wing, spar);
+    (void)_;
+    ax->plot(y,x)->color(spar_color_tandem).line_width(1).line_style("--");
+    legend_strings.push_back("spar");
+  }
+  for(auto& device : data->back_control_devices) {
+    const auto [x,y,name] = gen_polygon_tandem(data->back_wing, device);
+    ax->plot(y,x)->color(control_device_colormap_tandem[name]).line_width(1.5);
+    legend_strings.push_back(device.name);
+  }
+  auto limits = ax->ylim();
+  ax->plot(std::vector<double>({fuselage_half,fuselage_half}),std::vector<double>({limits[0]-2,limits[1]+2}))->color({0, 1, 0,0}).line_style("--").line_width(1.0);
+  legend_strings.push_back("fuselage width");
+  matplot::legend_handle lh;
+  ax->xlabel("y [m]");
+  ax->ylabel("x [m]");
+  ax->limits_mode(matplot::manual);
+  ax->axis(matplot::equal);
+  ax->xlim({0,*std::max_element(spanwise.begin(), spanwise.end())*1.1});
+  ax->legend(legend_strings);
+  ax->legend()->location(matplot::legend::general_alignment::bottomright);
+  ax->legend()->inside(true);
+  ax->legend()->box(true);
+  ax->legend()->font_size(8);
+
+
+  fs::path plot_path = rtIO->getPlotDir() + "/back_wing_planfrom.svg";
+  matplot::save(f, plot_path.string(), "svg");
+  std::this_thread::sleep_for(std::chrono::milliseconds(500));
+  system("pkill gnuplot");
+  return fs::relative(plot_path, fs::path(rtIO->getHtmlDir()));
+}
+
 fs::path Wing::plot_thickness_and_twist_distribution() {
   /* Plot active */
   if (!rtIO->plotOn) {
@@ -240,15 +314,77 @@ fs::path Wing::plot_thickness_and_twist_distribution() {
   ax2->ylim({0, 0.2});
   ax2->y2lim({-1, 1});
   */
-  fs::path plot_path = rtIO->getPlotDir() + "/thickness_and_twist_distribution.svg";
+  fs::path plot_path = rtIO->getPlotDir() + "/front_wing_thickness_and_twist_distribution.svg";
   matplot::save(f, plot_path.string(), "svg");
   std::this_thread::sleep_for(std::chrono::milliseconds(500));
   system("pkill gnuplot");
   return fs::relative(plot_path, fs::path(rtIO->getHtmlDir()));
 }
 
+fs::path Wing::plot_back_thickness_and_twist_distribution() {
+  /* Plot active */
+  if (!rtIO->plotOn) {
+    return "";
+  }
+  std::vector<double> eta;
+  std::vector<double> thickness_to_chord_ratio;
+  std::vector<double> twist;
+  double half_span = geom2::measure::span(data->back_wing) * 0.5;
+  for (auto section : data->back_wing.sections) {
+    eta.push_back(-section.origin.z() / half_span);
+    thickness_to_chord_ratio.push_back(section.get_thickness_scale() * geom2::measure::thickness_max(section) /
+                                       section.get_chord_length());
+    twist.push_back(section.get_twist_angle());
+  }
+
+  auto f = matplot::figure(true);
+
+  auto ax1 = f->add_subplot(2, 1, 0);
+  ax1->plot(eta, thickness_to_chord_ratio, "g-")->color("black").line_width(1.5);
+  // ax1->xlabel("\u03B7");
+  ax1->x_axis().visible(true);
+  ax1->grid(true);
+  ax1->ylabel("t/c");
+
+  auto ax2 = f->add_subplot(2, 1, 1);
+  ax2->plot(eta, twist, "b-")->color("black").line_style("--").line_width(1.5);
+  ax2->ylabel("\u03F5 [\u00b0]");
+  ax2->xlabel("\u03B7");
+  ax2->grid(true);
+  /*
+  // Plot the first set of data
+  ax1->plot(eta, thickness_to_chord_ratio, "g-")->use_y2(false).color("black").line_width(1.5);
+  ax1->xlabel("\u03B7");
+  ax1->ylabel("t/c");
+  ax1->limits_mode_automatic();
+  ax1->legend({"t/c"});
+
+  matplot::hold(true);
+  // Create a secondary y-axis
+  auto ax2 = f->add_axes(ax1, false, false);
+  ax2->plot(eta, twist, "b-")->use_y2(true).color("green").line_style("--").line_width(1.5);
+  ax2->y2_axis().limits({-1, 1});
+  // ax2->visible(false);
+  ax2->grid(true);
+  ax2->y_axis().color("black");
+  ax2->y2_axis().label("\u03F5 [\u00b0]");
+  ax2->legend({"t/c", "\u03F5"});
+
+  // Set the secondary y-axis invisible except for the y-axis itself
+  ax2->x_axis().visible(true);
+  ax2->y2_axis().color("black");
+  // Adjust limits and ticks if necessary
+  ax2->ylim({0, 0.2});
+  ax2->y2lim({-1, 1});
+  */
+  fs::path plot_path = rtIO->getPlotDir() + "/back_wing_thickness_and_twist_distribution.svg";
+  matplot::save(f, plot_path.string(), "svg");
+  std::this_thread::sleep_for(std::chrono::milliseconds(500));
+  system("pkill gnuplot");
+  return fs::relative(plot_path, fs::path(rtIO->getHtmlDir()));
+}
 /**
- * \brief Plot airfoils
+ * \brief Plot airfoils of front wing
  *
  * \return fs::path
  */
@@ -301,5 +437,60 @@ std::vector<fs::path> Wing::plot_airfoils() {
   return airfoil_plots_path;
 }
 
+/**
+ * \brief Plot airfoils of back wing
+ *
+ * \return fs::path
+ */
+std::vector<fs::path> Wing::plot_back_airfoils() {
+  /* Plot active */
+  if (!rtIO->plotOn) {
+    return {""};
+  }
+  std::map<std::string, std::vector<geom2::Point_2>> airfoils_data{};
+  std::vector<fs::path> airfoil_plots_path;
+
+  for (auto section : data->back_wing.sections) {
+    if (airfoils_data.find(section.name) != airfoils_data.end()) {
+      continue;
+    } else {
+      airfoils_data.insert({section.name, section.get_contour().vertices()});
+    }
+  }
+  /* loop over airfoils */
+
+  for (const auto& airfoil_data : airfoils_data) {
+    std::vector<double> x{}, y{};
+    for (auto coord : airfoil_data.second) {
+      x.push_back(coord.x());
+      y.push_back(coord.y());
+    }
+
+    {
+      auto f = matplot::figure(true);
+      auto h = f->current_axes()->plot(x, y);
+      h->line_width(1.5).color("black");
+
+      auto ax = f->current_axes();
+      ax->title(std::format("Airfoil - {}", airfoil_data.first));
+      matplot::xlabel("\u03B7");
+      ax->ylabel("z/c");
+      ax->xlabel("\u03B7");
+      ax->grid(true);
+      ax->ylim({-0.2, 0.2});
+      ax->xlim({0, 1});
+      ax->axis(matplot::equal);
+
+      fs::path plot_path = rtIO->getPlotDir() + "/" + airfoil_data.first + "_airfoil.svg";
+      matplot::save(f, plot_path.string(), "svg");
+      airfoil_plots_path.push_back(fs::relative(plot_path, fs::path(rtIO->getHtmlDir())));
+      std::this_thread::sleep_for(std::chrono::milliseconds(500));
+      system("pkill gnuplot");
+    }
+  }
+  return airfoil_plots_path;
+}
+
+
 }  // namespace cantilever
 }  // namespace tandem
diff --git a/wing_design/src/tandem/cantilever/cantileverTandemWingDesignReport.cpp b/wing_design/src/tandem/cantilever/cantileverTandemWingDesignReport.cpp
index 9a88f7a35f484ce5e480e4296f2018cb0285766c..67c68b6e16146ca094995eb2bb577b3f5eacaa1c 100644
--- a/wing_design/src/tandem/cantilever/cantileverTandemWingDesignReport.cpp
+++ b/wing_design/src/tandem/cantilever/cantileverTandemWingDesignReport.cpp
@@ -21,10 +21,123 @@ void Wing::set_html_body() {
   reporter.htmlReportStream() << "<div class=\"box data\">\n";
   /* Add headline Data*/
   reporter.htmlReportStream() << "<h2>Data</h2>\n";
-  /* Table 1 - General parameters */
+
+/* Table 1 - Positioning paramters */
+
+  reporter.htmlReportStream()
+  << "<table class=\"content-table\">\n"
+  << "<caption>Positioning parameters of tandem wing set</caption>\n"
+  << "<thead>\n"
+  << "<tr>\n"
+  << "<th>parameter</th><th>symbol</th><th>unit</th><th>value</th>\n"
+  << "</tr>\n"
+  << "</thead>\n"
+  << "<tbody>\n"
+  << "<tr>\n"
+  << std::format("<td>Horizontal LE position of front wing</td><td>x<sub>front,LE</sub></td><td>m</td><td>{:.4f}</td>\n",
+                 (data->wing_position.x.value() + data->wing.origin.x()))
+  << "</tr>\n"
+  << std::format("<td>Horizontal LE position of back wing</td><td>x<sub>back,LE</sub></td><td>m</td><td>{:.4f}</td>\n",
+                 (data->wing_position.x.value() + data->back_wing.origin.x()))
+  << "</tr>\n"
+  << "<tr>\n"
+  << std::format("<td>Verical LE position of front wing</td><td>z<sub>front,LE</sub></td><td>m</td><td>{:.4f}</td>\n",
+                 (data->wing_position.z.value() + data->wing.origin.z()))
+  << "</tr>\n"
+  << "<tr>\n"
+  << std::format("<td>Verical LE position of back wing</td><td>z<sub>back,LE</sub></td><td>m</td><td>{:.4f}</td>\n",
+                 (data->wing_position.z.value() + data->back_wing.origin.z()))
+  << "</tr>\n"
+  << "<tr>\n"
+  << std::format("<td>Leading edge stagger</td><td>St<sub>LE</sub></td><td>m</td><td>{:.4f}</td>\n",
+                 (data->back_wing.origin.x() - data->wing.origin.x()))
+  << "</tr>\n"
+  << "<tr>\n"
+  << std::format("<td>Relative leading edge stagger</td><td>St<sub>rel,LE</sub></td><td>1</td><td>{:.4f}</td>\n",
+                 (data->back_wing.origin.x() - data->wing.origin.x()) / geom2::measure::length(data->fuselage))
+  << "</tr>\n"
+  << "<tr>\n"
+  << std::format("<td>1/4 MAC stagger</td><td>St<sub>25_MAC</sub></td><td>m</td><td>{:.4f}</td>\n",
+                ((data->back_wing.origin.x()
+                + geom2::measure::offset_LE(data->back_wing, geom2::measure::mean_aerodynamic_chord_position(data->back_wing)).x() 
+                + 0.25 * geom2::measure::chord(data->back_wing, geom2::measure::mean_aerodynamic_chord_position(data->back_wing))
+                * cos(geom2::measure::twist(data->back_wing, geom2::measure::mean_aerodynamic_chord_position(data->back_wing))))
+                - (data->wing.origin.x()
+                + geom2::measure::offset_LE(data->wing, geom2::measure::mean_aerodynamic_chord_position(data->wing)).x() 
+                + 0.25 * geom2::measure::chord(data->wing, geom2::measure::mean_aerodynamic_chord_position(data->wing))
+                * cos(geom2::measure::twist(data->wing, geom2::measure::mean_aerodynamic_chord_position(data->wing))))) 
+                )  //TODO kayra: twist diger hesaplamalarda yoktu. Ben ekledim
+                            
+  << "</tr>\n"
+  << "<tr>\n"
+  << std::format("<td>Relative 1/4 MAC stagger</td><td>St<sub>rel,25_MAC</sub></td><td>1</td><td>{:.4f}</td>\n",
+                ((data->back_wing.origin.x()
+                + geom2::measure::offset_LE(data->back_wing, geom2::measure::mean_aerodynamic_chord_position(data->back_wing)).x() 
+                + 0.25 * geom2::measure::chord(data->back_wing, geom2::measure::mean_aerodynamic_chord_position(data->back_wing))
+                * cos(geom2::measure::twist(data->back_wing, geom2::measure::mean_aerodynamic_chord_position(data->back_wing))))
+                - (data->wing.origin.x()
+                + geom2::measure::offset_LE(data->wing, geom2::measure::mean_aerodynamic_chord_position(data->wing)).x() 
+                + 0.25 * geom2::measure::chord(data->wing, geom2::measure::mean_aerodynamic_chord_position(data->wing))
+                * cos(geom2::measure::twist(data->wing, geom2::measure::mean_aerodynamic_chord_position(data->wing))))) 
+                / geom2::measure::length(data->fuselage))
+                            
+  << "</tr>\n"
+  << "<tr>\n"
+  << std::format("<td>Vertical leading edge gap</td><td>G<sub>LE</sub></td><td>m</td><td>{:.4f}</td>\n", 
+                data->back_wing.origin.z() - data->wing.origin.z())
+  << "</tr>\n"
+  << "<tr>\n"
+  << std::format("<td>Relative vertical leading edge gap</td><td>G<sub>rel,LE</sub></td><td>1</td><td>{:.4f}</td>\n", 
+                (data->back_wing.origin.z() - data->wing.origin.z()) / geom2::measure::height_max(data->fuselage))
+  << "</tr>\n"
+  << "<tr>\n"
+  << std::format("<td>Vertical 1/4 MAC gap</td><td>G<sub>25_MAC</sub></td><td>m</td><td>{:.4f}</td>\n", 
+                ((data->back_wing.origin.z()
+                + geom2::measure::offset_LE(data->back_wing, geom2::measure::mean_aerodynamic_chord_position(data->back_wing)).y() 
+                + 0.25 * geom2::measure::chord(data->back_wing, geom2::measure::mean_aerodynamic_chord_position(data->back_wing))
+                * sin(geom2::measure::twist(data->back_wing, geom2::measure::mean_aerodynamic_chord_position(data->back_wing))))
+                - (data->wing.origin.z()
+                + geom2::measure::offset_LE(data->wing, geom2::measure::mean_aerodynamic_chord_position(data->wing)).y() 
+                + 0.25 * geom2::measure::chord(data->wing, geom2::measure::mean_aerodynamic_chord_position(data->wing))
+                * sin(geom2::measure::twist(data->wing, geom2::measure::mean_aerodynamic_chord_position(data->wing))))) 
+                )  
+  << "</tr>\n"
+
+  << "<tr>\n"
+  << std::format("<td>Relative vertical 1/4 MAC gap</td><td>G<sub>rel,25_MAC</sub></td><td>1</td><td>{:.4f}</td>\n", 
+                ((data->back_wing.origin.z()
+                + geom2::measure::offset_LE(data->back_wing, geom2::measure::mean_aerodynamic_chord_position(data->back_wing)).y() 
+                + 0.25 * geom2::measure::chord(data->back_wing, geom2::measure::mean_aerodynamic_chord_position(data->back_wing))
+                * sin(geom2::measure::twist(data->back_wing, geom2::measure::mean_aerodynamic_chord_position(data->back_wing))))
+                - (data->wing.origin.z()
+                + geom2::measure::offset_LE(data->wing, geom2::measure::mean_aerodynamic_chord_position(data->wing)).y() 
+                + 0.25 * geom2::measure::chord(data->wing, geom2::measure::mean_aerodynamic_chord_position(data->wing))
+                * sin(geom2::measure::twist(data->wing, geom2::measure::mean_aerodynamic_chord_position(data->wing))))) 
+                / geom2::measure::height_max(data->fuselage))
+  << "</tr>\n"
+
+  << "<tr>\n"
+  << std::format("<td>Fuselage length</td><td>l<sub>fus</sub></td><td>m</td><td>{:.4f}</td>\n",
+                 geom2::measure::length(data->fuselage))
+  << "</tr>\n"
+  << "<tr>\n"
+  << std::format("<td>Fuselage max. height</td><td>h<sub>fus,max</sub></td><td>m</td><td>{:.4f}</td>\n",
+                 geom2::measure::height_max(data->fuselage))
+  << "</tr>\n"
+  << "</tbody>\n"
+  << "</table>\n";
+  reporter.htmlReportStream()
+  << "<p style=\"font-size: 14px; font-style: italic;\">"
+  << "Note: Horizontal relative values are calculated as fractions of the fuselage length. <br>"
+  << "Vertical relative values are calculated as fractions of the maximum fuselage height."
+  << "</p>\n";
+
+
+
+  /* Table 2 - General parameters for Front Wing*/
   reporter.htmlReportStream()
       << "<table class=\"content-table\">\n"
-      << "<caption>General parameters</caption>\n"
+      << "<caption>General parameters of front wing</caption>\n"
       << "<thead>\n"
       << "<tr>\n"
       << "<th>parameter</th><th>symbol</th><th>unit</th><th>value</th>\n"
@@ -53,9 +166,41 @@ void Wing::set_html_body() {
       << "</tbody>\n"
       << "</table>\n";
 
-  /* Table 2 - Section parameters*/
+/* Table 3 - General parameters for Back Wing */
+  reporter.htmlReportStream()
+    << "<table class=\"content-table\">\n"
+    << "<caption>General parameters of back wing</caption>\n"
+    << "<thead>\n"
+    << "<tr>\n"
+    << "<th>parameter</th><th>symbol</th><th>unit</th><th>value</th>\n"
+    << "</tr>\n"
+    << "</thead>\n"
+    << "<tbody>\n"
+    << "<tr>\n"
+    << std::format("<td>Reference area</td><td>S<sub>ref</sub></td><td>m<sup>2</sup></td><td>{:.2f}</td>\n",
+                   geom2::measure::reference_area(data->back_wing))
+    << "</tr>\n"
+    << "<tr>\n"
+    << std::format("<td>Mean aero. chord</td><td>MAC</td><td>m</td><td>{:.2f}</td>\n",
+                   geom2::measure::mean_aerodynamic_chord(data->back_wing))
+    << "</tr>\n"
+    << "<tr>\n"
+    << std::format("<td>Span</td><td>b</td><td>m</td><td>{:.2f}</td>\n", geom2::measure::span(data->back_wing))
+    << "</tr>\n"
+    << "<tr>\n"
+    << std::format("<td>Aspect ratio</td><td>AR</td><td>1</td><td>{:.2f}</td>\n",
+                   geom2::measure::aspect_ratio(data->back_wing))
+    << "</tr>\n"
+    << "<tr>\n"
+    << std::format("<td>Taper ratio</td><td>&lambda;</td><td>1</td><td>{:.2f}</td>\n",
+                   geom2::measure::taper_ratio(data->back_wing))
+    << "</tr>\n"
+    << "</tbody>\n"
+    << "</table>\n";
+
+  /* Table 4 - Section parameters of Front Wing*/
   reporter.htmlReportStream() << "<table class=\"content-table\">\n"
-                              << "<caption>Section parameters</caption>\n"
+                              << "<caption>Section parameters of front wing</caption>\n"
                               << "<thead>\n"
                               << "<tr>\n"
                               << "<th>parameter</th><th>symbol</th><th>unit</th>";
@@ -143,11 +288,102 @@ void Wing::set_html_body() {
                               << "</tbody>\n"
                               << "</table>\n";
 
-  /* Table 3 - Spar parameters*/
+  /* Table 2 - Section parameters of Back Wing*/
+  reporter.htmlReportStream() << "<table class=\"content-table\">\n"
+                              << "<caption>Section parameters of back wing</caption>\n"
+                              << "<thead>\n"
+                              << "<tr>\n"
+                              << "<th>parameter</th><th>symbol</th><th>unit</th>";
+  for (size_t i = 0; i < data->back_wing.sections.size(); ++i) {
+    reporter.htmlReportStream() << std::format("<th>S<sub>{:02d}</sub></th>", i);
+  }
+  reporter.htmlReportStream() << "\n</tr>\n"
+                              << "</thead>\n"
+                              << "<tbody>\n"
+                              /* Dimensionless half span - eta*/
+                              << "<tr>\n"
+                              << "<td>Dimensionless half span</td><td>&eta;</td><td>1</td>";
+  for (auto section : data->back_wing.sections) {
+    reporter.htmlReportStream() << std::format("<td>{:.2f}</td>",
+                                               -2 * section.origin.z() / geom2::measure::span(data->back_wing));
+  }
+  /* Spanwise coordinate - y*/
+  reporter.htmlReportStream() << "\n</tr>\n"
+                              << "<td>Spanwise coordinate</td><td>y</td><td>m</td>";
+  for (auto section : data->back_wing.sections) {
+    reporter.htmlReportStream() << std::format("<td>{:.2f}</td>", -section.origin.z());
+  }
+  reporter.htmlReportStream() << "\n</tr>\n"
+                              /* Chord - c*/
+                              << "<tr>\n"
+                              << "<td>Chord</td><td>c</td><td>m</td>";
+  for (auto section : data->back_wing.sections) {
+    reporter.htmlReportStream() << std::format("<td>{:.2f}</td>", section.get_chord_length());
+  }
+  reporter.htmlReportStream() << "\n</tr>\n"
+                              /* Sweep - Leading edge - phi_LE*/
+                              << "<tr>\n"
+                              << "<td>Sweep leading edge</td><td>&phi;<sub>LE</sub></td><td>&deg;</td>";
+  for (auto section : data->back_wing.sections) {
+    reporter.htmlReportStream() << std::format(
+        "<td>{:.2f}</td>", -geom2::detail::to_degrees * geom2::measure::sweep(data->back_wing, section.origin.z(), 0.0));
+  }
+  reporter.htmlReportStream() << "\n</tr>\n"
+                              /* Sweep - 25 percent - phi_25%*/
+                              << "<tr>\n"
+                              << "<td>Sweep quarter chord</td><td>&phi;<sub>25</sub></td><td>&deg;</td>";
+  for (auto section : data->back_wing.sections) {
+    reporter.htmlReportStream() << std::format(
+        "<td>{:.2f}</td>", -geom2::detail::to_degrees * geom2::measure::sweep(data->back_wing, section.origin.z(), 0.25));
+  }
+  reporter.htmlReportStream() << "\n</tr>\n"
+                              /* Sweep - Trailing edge percent - phi_TE*/
+                              << "<tr>\n"
+                              << "<td>Sweep Trailing edge</td><td>&phi;<sub>TE</sub></td><td>&deg;</td>";
+  for (auto section : data->back_wing.sections) {
+    reporter.htmlReportStream() << std::format(
+        "<td>{:.2f}</td>", -geom2::detail::to_degrees * geom2::measure::sweep(data->back_wing, section.origin.z(), 1.0));
+  }
+  reporter.htmlReportStream() << "\n</tr>\n"
+                              /* Dihedral - nu*/
+                              << "<tr>\n"
+                              << "<td>Dihedral</td><td>&nu;</td><td>&deg;</td>";
+  for (auto section : data->back_wing.sections) {
+    reporter.htmlReportStream() << std::format(
+        "<td>{:.2f}</td>", geom2::detail::to_degrees * geom2::measure::dihedral(data->back_wing, section.origin.z()));
+  }
+  reporter.htmlReportStream() << "\n</tr>\n"
+                              /* Twist - epsilon*/
+                              << "<tr>\n"
+                              << "<td>Twist angle</td><td>&epsilon;</td><td>&deg;</td>";
+  for (auto section : data->back_wing.sections) {
+    reporter.htmlReportStream() << std::format("<td>{:.2f}</td>", geom2::detail::to_degrees * section.rotation_z);
+  }
+  reporter.htmlReportStream() << "\n</tr>\n"
+                              /* Thickness to chord - epsilon*/
+                              << "<tr>\n"
+                              << "<td>Thickness ratio</td><td>t/c</td><td>%</td>";
+  for (auto section : data->back_wing.sections) {
+    reporter.htmlReportStream() << std::format(
+        "<td>{:.2f}</td>", 100 * geom2::measure::thickness_max(section) / section.get_chord_length());
+  }
+  reporter.htmlReportStream() << "\n</tr>\n"
+                              /* Airfoil - none */
+                              << "<tr>\n"
+                              << "<td>Airfoil</td><td>-</td><td>-</td>";
+  for (auto section : data->back_wing.sections) {
+    reporter.htmlReportStream() << std::format("<td>{}</td>", section.name);
+  }
+  reporter.htmlReportStream() << "\n</tr>\n"
+                              << "</tbody>\n"
+                              << "</table>\n";
+
+
+  /* Table 5 - Spar parameters of Front Wing*/
   for (size_t spar_id = 0; spar_id < data->spars.size(); ++spar_id) {
     reporter.htmlReportStream() << "<table class=\"content-table\">\n";
     if (spar_id == 0) {
-      reporter.htmlReportStream() << "<caption>Spar parameters</caption>\n";
+      reporter.htmlReportStream() << "<caption>Spar parameters of front wing</caption>\n";
     }
     reporter.htmlReportStream() << "<thead>\n"
                                 << "<tr>\n"
@@ -175,7 +411,39 @@ void Wing::set_html_body() {
                                 << "</table>\n";
   }
 
-  /* Sort devices to leading, trailing and other devices*/
+  /* Table 6 - Spar parameters of Back Wing*/
+  for (size_t spar_id = 0; spar_id < data->back_spars.size(); ++spar_id) {
+    reporter.htmlReportStream() << "<table class=\"content-table\">\n";
+    if (spar_id == 0) {
+      reporter.htmlReportStream() << "<caption>Spar parameters of back wing</caption>\n";
+    }
+    reporter.htmlReportStream() << "<thead>\n"
+                                << "<tr>\n"
+                                << std::format("<th>{:<10}</th><th>&eta; [1]<th>x/c [%]</th>",
+                                               data->back_spars.at(spar_id).name)
+                                << "\n</tr>\n"
+                                << "</thead>\n"
+                                << "<tbody>\n";
+    /* Spar name - eta, x/c*/
+    for (size_t i = 0; i < data->back_spars.at(spar_id).sections.size(); ++i) {
+      reporter.htmlReportStream() << "<tr>\n";
+      std::string tag = "";
+      if (i == 0) {
+        tag = "from";
+      } else if (i == data->back_spars.at(spar_id).sections.size() - 1) {
+        tag = "to";
+      }
+      reporter.htmlReportStream() << std::format("<td>{}</td><td>{:.2f}</td><td>{:.2f}</td>\n", tag,
+                                                 data->back_spars.at(spar_id).sections.at(i).origin.z(),
+                                                 data->back_spars.at(spar_id).sections.at(i).get_contour().vertex(0).x());
+      reporter.htmlReportStream() << "\n</tr>\n";
+    }
+    /* Spanwise coordinate - y*/
+    reporter.htmlReportStream() << "</tbody>\n"
+                                << "</table>\n";
+  }
+
+  /* Sort devices to leading, trailing and other devices for front wing*/
   /* device data - name, eta_from, rel_chord_from_start, rel_chord_from_end, eta_to, rel_chord_to_start,
    * rel_chord_to_end */
   using device_data = std::tuple<std::string, double, double, double, double, double, double>;
@@ -209,7 +477,7 @@ void Wing::set_html_body() {
   if (!leading_edge_devices.empty()) {
     reporter.htmlReportStream()
         << "<table class=\"content-table\">\n"
-        << "<caption>Leading edge control devices</caption>\n"
+        << "<caption>Leading edge control devices of front wing</caption>\n"
         << "<thead>\n"
         << "<tr>\n"
         << "<th>Device</th><th>&eta;<sub>from</sub></th><th>x/c<sub>from,start</sub></th><th>x/c<sub>from,end</sub></"
@@ -235,7 +503,7 @@ void Wing::set_html_body() {
   if (!trailing_edge_devices.empty()) {
     reporter.htmlReportStream()
         << "<table class=\"content-table\">\n"
-        << "<caption>Trailing edge control devices</caption>\n"
+        << "<caption>Trailing edge control devices of front wing</caption>\n"
         << "<thead>\n"
         << "<tr>\n"
         << "<th>Device</th><th>&eta;<sub>from</sub></th><th>x/c<sub>from,start</sub></th><th>x/c<sub>from,end</sub></"
@@ -260,7 +528,7 @@ void Wing::set_html_body() {
   if (!other_devices.empty()) {
     reporter.htmlReportStream()
         << "<table class=\"content-table\">\n"
-        << "<caption>Other control devices</caption>\n"
+        << "<caption>Other control devices of front wing</caption>\n"
         << "<thead>\n"
         << "<tr>\n"
         << "<th>Device</th><th>&eta;<sub>from</sub></th><th>x/c<sub>from,start</sub></th><th>x/c<sub>from,end</sub></"
@@ -283,8 +551,116 @@ void Wing::set_html_body() {
                                 << "</table>\n";
   }
 
+    /* Sort devices to leading, trailing and other devices for back wing*/
+  /* device data - name, eta_from, rel_chord_from_start, rel_chord_from_end, eta_to, rel_chord_to_start,
+   * rel_chord_to_end */
+  using back_device_data = std::tuple<std::string, double, double, double, double, double, double>;
+  std::vector<back_device_data> back_leading_edge_devices;
+  std::vector<back_device_data> back_trailing_edge_devices;
+  std::vector<back_device_data> back_other_devices;
+  for (auto device : data->back_control_devices) {
+    /* Device has only two sections (from + to) spanwise */
+    double eta_from = device.sections.at(0).origin.z();
+    double rel_chord_from_start = device.sections.at(0).get_contour().vertex(0).x(),
+           rel_chord_from_end = device.sections.at(0).get_contour().vertex(1).x();
+    double eta_to = device.sections.at(1).origin.z();
+    double rel_chord_to_start = device.sections.at(1).get_contour().vertex(0).x(),
+           rel_chord_to_end = device.sections.at(1).get_contour().vertex(1).x();
+
+    if (std::fabs(rel_chord_from_start - rel_chord_to_start) < ACCURACY_MEDIUM &&
+        std::fabs(rel_chord_from_start) < ACCURACY_MEDIUM) {
+      back_leading_edge_devices.push_back({device.name, eta_from, rel_chord_from_start, rel_chord_from_end, eta_to,
+                                      rel_chord_to_start, rel_chord_to_end});
+    } else if (std::fabs(rel_chord_from_end - rel_chord_to_end) < ACCURACY_MEDIUM &&
+               std::fabs(rel_chord_from_end - 1) < ACCURACY_MEDIUM) {
+      back_trailing_edge_devices.push_back({device.name, eta_from, rel_chord_from_start, rel_chord_from_end, eta_to,
+                                       rel_chord_to_start, rel_chord_to_end});
+    } else {
+      back_other_devices.push_back({device.name, eta_from, rel_chord_from_start, rel_chord_from_end, eta_to,
+                               rel_chord_to_start, rel_chord_to_end});
+    }
+  }
+
+  /* Leading edge table*/
+  if (!back_leading_edge_devices.empty()) {
+    reporter.htmlReportStream()
+        << "<table class=\"content-table\">\n"
+        << "<caption>Leading edge control devices of back wing</caption>\n"
+        << "<thead>\n"
+        << "<tr>\n"
+        << "<th>Device</th><th>&eta;<sub>from</sub></th><th>x/c<sub>from,start</sub></th><th>x/c<sub>from,end</sub></"
+           "th><th>&eta;<sub>to</sub></th><th>x/c<sub>to,start</sub></th><th>x/c<sub>to,end</sub></th>\n"
+        << "</tr>\n"
+        << "</thead>\n"
+        << "<tbody>\n";
+    for (auto device : back_leading_edge_devices) {
+      const auto [name, eta_from, rel_chord_from_start, rel_chord_from_end, eta_to, rel_chord_to_start,
+                  rel_chord_to_end] = device;
+      reporter.htmlReportStream() << "<tr>\n";
+      reporter.htmlReportStream() << std::format(
+          "<td>{}</td><td>{:.2f}</td><td>{:.2f}</td><td>{:.2f}</td><td>{:.2f}</td><td>{:.2f}</td><td>{:.2f}</td>\n",
+          name, eta_from, rel_chord_from_start, rel_chord_from_end, eta_to, rel_chord_to_start,
+          rel_chord_to_end);
+      reporter.htmlReportStream() << "\n</tr>\n";
+    }
+    /* Spanwise coordinate - y*/
+    reporter.htmlReportStream() << "</tbody>\n"
+                                << "</table>\n";
+  }
+
+  if (!back_trailing_edge_devices.empty()) {
+    reporter.htmlReportStream()
+        << "<table class=\"content-table\">\n"
+        << "<caption>Trailing edge control devices of back wing</caption>\n"
+        << "<thead>\n"
+        << "<tr>\n"
+        << "<th>Device</th><th>&eta;<sub>from</sub></th><th>x/c<sub>from,start</sub></th><th>x/c<sub>from,end</sub></"
+           "th><th>&eta;<sub>to</sub></th><th>x/c<sub>to,start</sub></th><th>x/c<sub>to,end</sub></th>\n"
+        << "</tr>\n"
+        << "</thead>\n"
+        << "<tbody>\n";
+    for (auto device : back_trailing_edge_devices) {
+      const auto [name, eta_from, rel_chord_from_start, rel_chord_from_end, eta_to, rel_chord_to_start,
+                  rel_chord_to_end] = device;
+      reporter.htmlReportStream() << "<tr>\n";
+      reporter.htmlReportStream() << std::format(
+          "<td>{}</td><td>{:.2f}</td><td>{:.2f}</td><td>{:.2f}</td><td>{:.2f}</td><td>{:.2f}</td><td>{:.2f}</td>\n",
+          name, eta_from, rel_chord_from_start, rel_chord_from_end, eta_to, rel_chord_to_start,
+          rel_chord_to_end);
+      reporter.htmlReportStream() << "\n</tr>\n";
+    }
+    reporter.htmlReportStream() << "</tbody>\n"
+                                << "</table>\n";
+  }
+
+  if (!back_other_devices.empty()) {
+    reporter.htmlReportStream()
+        << "<table class=\"content-table\">\n"
+        << "<caption>Other control devices of back wing</caption>\n"
+        << "<thead>\n"
+        << "<tr>\n"
+        << "<th>Device</th><th>&eta;<sub>from</sub></th><th>x/c<sub>from,start</sub></th><th>x/c<sub>from,end</sub></"
+           "th><th>&eta;<sub>to</sub></th><th>x/c<sub>to,start</sub></th><th>x/c<sub>to,end</sub></th>\n"
+        << "</tr>\n"
+        << "</thead>\n"
+        << "<tbody>\n";
+    for (auto device : back_other_devices) {
+      const auto [name, eta_from, rel_chord_from_start, rel_chord_from_end, eta_to, rel_chord_to_start,
+                  rel_chord_to_end] = device;
+      reporter.htmlReportStream() << "<tr>\n";
+      reporter.htmlReportStream() << std::format(
+          "<td>{}</td><td>{:.2f}</td><td>{:.2f}</td><td>{:.2f}</td><td>{:.2f}</td><td>{:.2f}</td><td>{:.2f}</td>\n",
+          name, eta_from, rel_chord_from_start, rel_chord_from_end, eta_to, rel_chord_to_start,
+          rel_chord_to_end);
+      reporter.htmlReportStream() << "\n</tr>\n";
+    }
+    /* Spanwise coordinate - y*/
+    reporter.htmlReportStream() << "</tbody>\n"
+                                << "</table>\n";
+  }
+
   reporter.htmlReportStream() << "<table class=\"content-table\">\n"
-                              << "<caption>Mass and CoG</caption>\n"
+                              << "<caption>Mass and CoG of tandem wing set</caption>\n"
                               << "<thead>\n"
                               << "<tr>\n"
                               << "<th>Mass</th><th>CoG<sub>x</sub></th><th>CoG<sub>y</sub></th><th>CoG<sub>z</sub></th>\n"
@@ -308,19 +684,34 @@ void Wing::set_html_body() {
   /* Add headline Data*/
   reporter.htmlReportStream() << "<h2>Plots</h2>\n";
 
-  reporter.htmlReportStream() << "<h3>Geometry</h3>\n";
+  reporter.htmlReportStream() << "<h3>Geometry of front wing</h3>\n";
   auto planform_plot_path = plot_wing_planform();
   reporter.htmlReportStream() << std::format("<img class=\"image-plot\" src=\"{}\"/>\n",planform_plot_path.string());
 
+  reporter.htmlReportStream() << "<h3>Geometry of back wing</h3>\n";
+  auto back_planform_plot_path = plot_back_wing_planform();
+  reporter.htmlReportStream() << std::format("<img class=\"image-plot\" src=\"{}\"/>\n",back_planform_plot_path.string());
+
   auto airfoil_plot_paths = plot_airfoils();
-  reporter.htmlReportStream() << "<h3>Airfoils</h3>\n";
+  reporter.htmlReportStream() << "<h3>Airfoils of front wing</h3>\n";
   for( auto airfoil_path : airfoil_plot_paths) {
     reporter.htmlReportStream() << std::format("<img class=\"image-plot\" src=\"{}\"/>\n", airfoil_path.string());
   }
-  reporter.htmlReportStream() << "<h3>Thickness & Geometric Twist</h3>\n";
+
+  auto back_airfoil_plot_paths = plot_back_airfoils();
+  reporter.htmlReportStream() << "<h3>Airfoils of back wing</h3>\n";
+  for( auto airfoil_path : back_airfoil_plot_paths) {
+    reporter.htmlReportStream() << std::format("<img class=\"image-plot\" src=\"{}\"/>\n", airfoil_path.string());
+  }
+
+  reporter.htmlReportStream() << "<h3>Thickness & Geometric Twist of front wing</h3>\n";
   auto thickness_and_twist_plot_path = plot_thickness_and_twist_distribution();
   reporter.htmlReportStream() << std::format("<img class=\"image-plot\" src=\"{}\"/>\n", thickness_and_twist_plot_path.string());
 
+  reporter.htmlReportStream() << "<h3>Thickness & Geometric Twist of back wing</h3>\n";
+  auto back_thickness_and_twist_plot_path = plot_thickness_and_twist_distribution();
+  reporter.htmlReportStream() << std::format("<img class=\"image-plot\" src=\"{}\"/>\n", back_thickness_and_twist_plot_path.string());
+
   /* close box plot */
   reporter.htmlReportStream() << "</div>\n";
 }
diff --git a/wing_design/src/tandem/tandemWingDesignConfig.cpp b/wing_design/src/tandem/tandemWingDesignConfig.cpp
index c83c2141abfcbab1af91e8e81dedbc0ccda7b3b3..9694f1ed0e125db845b3fb001f49941d5d012ef7 100644
--- a/wing_design/src/tandem/tandemWingDesignConfig.cpp
+++ b/wing_design/src/tandem/tandemWingDesignConfig.cpp
@@ -5,9 +5,9 @@
 tandem::cantilever::WingDesignConfig::WingDesignConfig()
     : wing_design_mode(EndnodeReadOnly<std::string>("program_settings/modes/design_mode")),
       wing_mass_mode(
-          EndnodeReadOnly<std::string>("program_settings/tandem/cantilever/calculation_methods/front_wing/mass/method")),
-      back_wing_mass_mode(
-          EndnodeReadOnly<std::string>("program_settings/tandem/cantilever/calculation_methods/back_wing/mass/method")),         
+          EndnodeReadOnly<std::string>("program_settings/tandem/cantilever/calculation_methods/mass_method/method")),
+      //back_wing_mass_mode(
+          //EndnodeReadOnly<std::string>("program_settings/tandem/cantilever/calculation_methods/back_wing/mass/method")),         
       wing_configuration_mode(EndnodeReadOnly<std::string>("program_settings/tandem/wing_configuration")),
       total_wing_area_calculation_mode(EndnodeReadOnly<std::string>(
           "program_settings/tandem/cantilever/calculation_methods/total_wing_area/method")),
@@ -21,7 +21,7 @@ tandem::cantilever::WingDesignConfig::WingDesignConfig()
           "program_settings/tandem/cantilever/calculation_methods/relative_LE_stagger/parameters/mode_0/relative_LE_stagger")),
 
       sweep_calculation_mode(
-          EndnodeReadOnly<std::string>("program_settings/tandem/cantilever/calculation_methods/front_wing/sweep/method")),
+          EndnodeReadOnly<std::string>("program_settings/tandem/cantilever/calculation_methods/sweep_method/method")),
       user_sweep_angle_at_quarter_chord(EndnodeReadOnly<double>(
           "program_settings/tandem/cantilever/calculation_methods/front_wing/sweep/parameters/mode_0/sweep_angle")),
       korn_technology_factor(EndnodeReadOnly<double>("program_settings/tandem/cantilever/calculation_methods/front_wing/"
@@ -29,8 +29,8 @@ tandem::cantilever::WingDesignConfig::WingDesignConfig()
       delta_drag_divergence_to_mach_design(
           EndnodeReadOnly<double>("program_settings/tandem/cantilever/calculation_methods/front_wing/sweep/parameters/"
                                   "mode_1/delta_drag_divergence_to_mach_design")),
-      back_sweep_calculation_mode(
-          EndnodeReadOnly<std::string>("program_settings/tandem/cantilever/calculation_methods/back_wing/sweep/method")),
+      //back_sweep_calculation_mode(
+          //EndnodeReadOnly<std::string>("program_settings/tandem/cantilever/calculation_methods/back_wing/sweep/method")), //TODO kayra:delete
       user_back_sweep_angle_at_quarter_chord(EndnodeReadOnly<double>(
           "program_settings/tandem/cantilever/calculation_methods/back_wing/sweep/parameters/mode_0/sweep_angle")),
       back_korn_technology_factor(EndnodeReadOnly<double>("program_settings/tandem/cantilever/calculation_methods/back_wing/"
@@ -40,32 +40,32 @@ tandem::cantilever::WingDesignConfig::WingDesignConfig()
                                   "mode_1/delta_drag_divergence_to_mach_design")),
 
       taper_ratio_calculation_mode(EndnodeReadOnly<std::string>(
-          "program_settings/tandem/cantilever/calculation_methods/front_wing/taper_ratio/method")),
+          "program_settings/tandem/cantilever/calculation_methods/taper_ratio_method/method")),
       user_taper_ratio(EndnodeReadOnly<double>(
           "program_settings/tandem/cantilever/calculation_methods/front_wing/taper_ratio/parameters/mode_0/taper_ratio")),
-      back_taper_ratio_calculation_mode(EndnodeReadOnly<std::string>(
-          "program_settings/tandem/cantilever/calculation_methods/back_wing/taper_ratio/method")),
+      //back_taper_ratio_calculation_mode(EndnodeReadOnly<std::string>(
+          //"program_settings/tandem/cantilever/calculation_methods/back_wing/taper_ratio/method")), //TODO kayra:delete
       user_back_taper_ratio(EndnodeReadOnly<double>(
           "program_settings/tandem/cantilever/calculation_methods/back_wing/taper_ratio/parameters/mode_0/taper_ratio")),          
       
       aspect_ratio_calculation_mode(EndnodeReadOnly<std::string>(
-          "program_settings/tandem/cantilever/calculation_methods/front_wing/aspect_ratio/method")),
+          "program_settings/tandem/cantilever/calculation_methods/aspect_ratio_method/method")),
       user_aspect_ratio(EndnodeReadOnly<double>(
           "program_settings/tandem/cantilever/calculation_methods/front_wing/aspect_ratio/parameters/mode_0/aspect_ratio")),
-      back_aspect_ratio_calculation_mode(EndnodeReadOnly<std::string>(
-          "program_settings/tandem/cantilever/calculation_methods/back_wing/aspect_ratio/method")),
+      //back_aspect_ratio_calculation_mode(EndnodeReadOnly<std::string>(
+          //"program_settings/tandem/cantilever/calculation_methods/back_wing/aspect_ratio/method")),  //TODO kayra:delete
       user_back_aspect_ratio(EndnodeReadOnly<double>(
           "program_settings/tandem/cantilever/calculation_methods/back_wing/aspect_ratio/parameters/mode_0/aspect_ratio")),
 
       dihedral_calculation_mode(EndnodeReadOnly<std::string>(
-          "program_settings/tandem/cantilever/calculation_methods/front_wing/dihedral/method")),
+          "program_settings/tandem/cantilever/calculation_methods/dihedral_method/method")),
       user_dihedral(EndnodeReadOnly<double>(
           "program_settings/tandem/cantilever/calculation_methods/front_wing/dihedral/parameters/mode_0/dihedral_angle")),
       dihedral_limitation_mode(
           EndnodeReadOnly<std::string>("program_settings/tandem/cantilever/calculation_methods/front_wing/dihedral/"
                                        "parameters/mode_1/dihedral_limitation")),
-      back_dihedral_calculation_mode(EndnodeReadOnly<std::string>(
-          "program_settings/tandem/cantilever/calculation_methods/back_wing/dihedral/method")),
+      //back_dihedral_calculation_mode(EndnodeReadOnly<std::string>(
+          //"program_settings/tandem/cantilever/calculation_methods/back_wing/dihedral/method")),  //TODO kayra:delete
       user_back_dihedral(EndnodeReadOnly<double>(
           "program_settings/tandem/cantilever/calculation_methods/back_wing/dihedral/parameters/mode_0/dihedral_angle")),
       back_dihedral_limitation_mode(
@@ -180,7 +180,7 @@ tandem::cantilever::WingDesignConfig::WingDesignConfig()
 void tandem::cantilever::WingDesignConfig::read(const node& xml) {
   wing_design_mode.read(xml);
   wing_mass_mode.read(xml);
-  back_wing_mass_mode.read(xml);
+  //back_wing_mass_mode.read(xml);
 
   wing_configuration_mode.read(xml);
   total_wing_area_calculation_mode.read(xml);
@@ -197,25 +197,25 @@ void tandem::cantilever::WingDesignConfig::read(const node& xml) {
   korn_technology_factor.read(xml);
   delta_drag_divergence_to_mach_design.read(xml);
   
-  back_sweep_calculation_mode.read(xml);
+  //back_sweep_calculation_mode.read(xml); //TODO kayra:delete
   user_back_sweep_angle_at_quarter_chord.read(xml);
   back_korn_technology_factor.read(xml);
   back_delta_drag_divergence_to_mach_design.read(xml);
 
   taper_ratio_calculation_mode.read(xml);
   user_taper_ratio.read(xml);
-  back_taper_ratio_calculation_mode.read(xml);
+  //back_taper_ratio_calculation_mode.read(xml); //TODO kayra:delete
   user_back_taper_ratio.read(xml);
 
   aspect_ratio_calculation_mode.read(xml);
   user_aspect_ratio.read(xml);
-  back_aspect_ratio_calculation_mode.read(xml);
+  //back_aspect_ratio_calculation_mode.read(xml); //TODO kayra:delete
   user_back_aspect_ratio.read(xml);
 
   dihedral_calculation_mode.read(xml);
   user_dihedral.read(xml);
   dihedral_limitation_mode.read(xml);
-  back_dihedral_calculation_mode.read(xml);
+  //back_dihedral_calculation_mode.read(xml);  //TODO kayra:delete
   user_back_dihedral.read(xml);
   back_dihedral_limitation_mode.read(xml);
 
diff --git a/wing_design/src/tandem/tandemWingDesignConfig.h b/wing_design/src/tandem/tandemWingDesignConfig.h
index 9a7910106ab8e1a660ef32dfb1f3a91c7c08bc99..909e9d9b7b4196084981fa452ed7744b20dacbfc 100644
--- a/wing_design/src/tandem/tandemWingDesignConfig.h
+++ b/wing_design/src/tandem/tandemWingDesignConfig.h
@@ -34,7 +34,7 @@ class WingDesignConfig {
   EndnodeReadOnly<std::string> wing_design_mode;
   EndnodeReadOnly<std::string> wing_configuration_mode;
   EndnodeReadOnly<std::string> wing_mass_mode;
-  EndnodeReadOnly<std::string> back_wing_mass_mode;
+  //EndnodeReadOnly<std::string> back_wing_mass_mode;
 
   /* Design */
   /* Data for method - wing area calculation */
@@ -54,7 +54,7 @@ class WingDesignConfig {
   EndnodeReadOnly<double> korn_technology_factor;
   EndnodeReadOnly<double> delta_drag_divergence_to_mach_design;
 
-  EndnodeReadOnly<std::string> back_sweep_calculation_mode;
+  //EndnodeReadOnly<std::string> back_sweep_calculation_mode; //TODO kayra:delete
   EndnodeReadOnly<double> user_back_sweep_angle_at_quarter_chord;
   EndnodeReadOnly<double> back_korn_technology_factor;
   EndnodeReadOnly<double> back_delta_drag_divergence_to_mach_design;
@@ -64,22 +64,22 @@ class WingDesignConfig {
   EndnodeReadOnly<std::string> taper_ratio_calculation_mode;
   EndnodeReadOnly<double> user_taper_ratio;
 
-  EndnodeReadOnly<std::string> back_taper_ratio_calculation_mode;
+  //EndnodeReadOnly<std::string> back_taper_ratio_calculation_mode; //TODO kayra:delete
   EndnodeReadOnly<double> user_back_taper_ratio;
 
   /* Data for method - aspect ratio calculation */
   EndnodeReadOnly<std::string> aspect_ratio_calculation_mode;
   EndnodeReadOnly<double> user_aspect_ratio;
 
-  EndnodeReadOnly<std::string> back_aspect_ratio_calculation_mode;
+  //EndnodeReadOnly<std::string> back_aspect_ratio_calculation_mode;  //TODO kayra:delete
   EndnodeReadOnly<double> user_back_aspect_ratio;
 
   /* Data for method - dihedral calculation*/
   EndnodeReadOnly<std::string> dihedral_calculation_mode;
-  EndnodeReadOnly<std::string> dihedral_limitation_mode;
+  EndnodeReadOnly<std::string> dihedral_limitation_mode;      //TODO kayra: bunu d birlestirmen lazim bir ara
   EndnodeReadOnly<double> user_dihedral;
 
-  EndnodeReadOnly<std::string> back_dihedral_calculation_mode;
+  //EndnodeReadOnly<std::string> back_dihedral_calculation_mode;   //TODO kayra:delete
   EndnodeReadOnly<std::string> back_dihedral_limitation_mode;
   EndnodeReadOnly<double> user_back_dihedral;
 
diff --git a/wing_design/src/tandem/tandemWingDesignData.cpp b/wing_design/src/tandem/tandemWingDesignData.cpp
index 115234e70c7e1ac02d23509cf2e60a535f1560e7..dce7bd117a359254ffe54fc03c19c4aa5233464a 100644
--- a/wing_design/src/tandem/tandemWingDesignData.cpp
+++ b/wing_design/src/tandem/tandemWingDesignData.cpp
@@ -37,9 +37,9 @@ tandem::cantilever::WingDesignData::WingDesignData()
       specification_number_of_wing_mounted_engines(0),
       sizing_point_wing_loading(EndnodeReadOnly<double>("sizing_point/wing_loading")),
       mtom_mass_properties(MassPropertiesIO("analysis/masses_cg_inertia/maximum_takeoff_mass","maximum takeoff mass")),
-      wing_mass_properties(MassPropertiesIO("component_design/wing", "wing")),
-      back_wing_mass_properties(MassPropertiesIO("component_design/back_wing", "Back Wing")),
-      wing_position(PositionIO("component_design/wing/position", "tandem wing set")) {}
+      wing_mass_properties(MassPropertiesIO("component_design/wing", "Tandem Wing Set")),
+      //back_wing_mass_properties(MassPropertiesIO("component_design/back_wing", "Back Wing")), //TODO kayra:silinecek
+      wing_position(PositionIO("component_design/wing/position", "Tandem Wing Set")) {}
 
 
 void tandem::cantilever::WingDesignData::read(const node& xml) {
@@ -202,7 +202,7 @@ void tandem::cantilever::WingDesignData::update(node& xml) {
   wing_mass_properties.update(xml); 
   
 
-  back_wing_mass_properties.update(xml); 
+  //back_wing_mass_properties.update(xml);  //TODO kayra: silinecek
   
 
   wing_position.update(xml);
diff --git a/wing_design/src/tandem/tandemWingDesignData.h b/wing_design/src/tandem/tandemWingDesignData.h
index 57baccd4f234f1a2d1b81e6d31b0b861fdf3b4bf..a30a354fc958482e13c7afca9988f391623f6830 100644
--- a/wing_design/src/tandem/tandemWingDesignData.h
+++ b/wing_design/src/tandem/tandemWingDesignData.h
@@ -56,13 +56,13 @@ class WingDesignData {
   EndnodeReadOnly<double> track_based_back_relative_kink;
 
   Endnode<double> spanwise_kink;
-  size_t specification_number_of_wing_mounted_engines;     //TODO kayra: bununla ilgili seyleri back winge aynen aldim ayarlama lazim olabilir
+  size_t specification_number_of_wing_mounted_engines;     //TODO kayra: I got everything for engines exactly same to back wing
   std::vector<std::tuple<double,double>> wing_mounted_engine_data;
   bool engines_exist = true;
 
   MassPropertiesIO mtom_mass_properties;
   MassPropertiesIO wing_mass_properties;
-  MassPropertiesIO back_wing_mass_properties;
+  //MassPropertiesIO back_wing_mass_properties; //TODO kayra: silinecek
   
   PositionIO wing_position;
   //PositionIO back_wing_position; //TODO kayra: silinecek
diff --git a/wing_design/wing_design_conf.xml b/wing_design/wing_design_conf.xml
index 10c83f691362941a97076c95fc04e58d50884dc2..91fbda97ef0b1c15f80374b5853fc5d8f6e20db3 100644
--- a/wing_design/wing_design_conf.xml
+++ b/wing_design/wing_design_conf.xml
@@ -2,10 +2,10 @@
 <module_configuration_file name="wing_design_conf.xml">
 	<control_settings description="General control settings for this tool">
 		<aircraft_exchange_file_name description="Specify the name of the exchange file">
-            <value>csmr-2020.xml</value>
+            <value>csmr-2020-voll.xml</value>
         </aircraft_exchange_file_name>
 		<aircraft_exchange_file_directory description="Specify the direction in which the aircraft exchange file can be found">
-			<value>../projects/</value>
+			<value>../projects/CSMR_voll/</value>
 		</aircraft_exchange_file_directory>
 		<own_tool_level description="Specify the tool level of this tool">
 			<value>1</value>
@@ -756,7 +756,7 @@
 								<wing_area_ratio description="wing area ratio">
 									<value>1</value>
 									<unit>1</unit>
-									<lower_boundary>0.9</lower_boundary>
+									<lower_boundary>0.83</lower_boundary>
 									<upper_boundary>1.2</upper_boundary>
 								</wing_area_ratio>
 							</mode_0>
@@ -770,7 +770,7 @@
 						<parameters>
 							<mode_0 description="user_defined">
 								<relative_LE_stagger description="relative leading edge stagger">
-									<value>0.5</value>
+									<value>0.3</value>
 									<unit>1</unit>
 									<lower_boundary>0</lower_boundary>
 									<upper_boundary>1</upper_boundary>
@@ -778,12 +778,39 @@
 							</mode_0>
 						</parameters>
 					</relative_LE_stagger>
+					
+					<mass_method description="switch for mass calculation method, look below for parameters">
+						<method description="selector mode_0: flops, mode_1: chiozzotto_wer">
+							<value>mode_0</value>
+						</method>
+					</mass_method>
+					
+					<sweep_method description="switch for sweep calculation method, look below for parameters">
+						<method description="selector mode_0: user_defined, mode_1: drag_divergence">
+							<value>mode_0</value>
+						</method>
+					</sweep_method>
+
+					<taper_ratio_method description="switch for taper ratio calculation method, look below for parameters">
+						<method description="selector mode_0: user_defined, mode_1: howe">
+							<value>mode_0</value>
+						</method>
+					</taper_ratio_method>
+
+					<dihedral_method description="switch for dihedral calculation method, look below for parameters">
+						<method description="selector mode_0: user_defined, mode_1: by_wing_position_and_quarter_chord_sweep">
+							<value>mode_0</value>
+						</method>
+					</dihedral_method>
+
+					<aspect_ratio_method description="switch for aspect ratio calculation method, look below for parameters">
+						<method description="selector mode_0: user_defined, mode_1: by_pitch_up_limit_function">
+							<value>mode_0</value>
+						</method>
+					</aspect_ratio_method>
 
 					<front_wing>
-						<sweep description="sweep calculation method">
-							<method description="selector mode_0: user_defined, mode_1: drag_divergence">
-								<value>mode_0</value>
-							</method>
+						<sweep description="parameters for sweep calculation methods">
 							<parameters description="sweep method parameters">
 								<mode_0 description="user_defined">
 									<sweep_angle description="sweep angle at quarter chord">
@@ -809,10 +836,7 @@
 								</mode_1>
 							</parameters>
 						</sweep>
-						<taper_ratio description="taper ratio calculation method">
-							<method description="selector mode_0: user_defined, mode_1: howe">
-								<value>mode_0</value>
-							</method>
+						<taper_ratio description="parameters for taper ratio calculation methods">
 							<parameters>
 								<mode_0 description="user_defined">
 									<taper_ratio description="taper ratio of overall wing">
@@ -824,10 +848,7 @@
 								</mode_0>
 							</parameters>
 						</taper_ratio>
-						<dihedral description="dihedral calculation method">
-							<method description="selector mode_0: user_defined, mode_1: by_wing_position_and_quarter_chord_sweep">
-								<value>mode_0</value>
-							</method>
+						<dihedral description="parameters for dihedral calculation methods">
 							<parameters description="parameters for methods">
 								<mode_0 description="user_defined">
 									<dihedral_angle description="dihedral angle">
@@ -844,10 +865,7 @@
 								</mode_1>
 							</parameters>
 						</dihedral>
-						<aspect_ratio description="aspect ratio calculation method">
-							<method description="selector mode_0: user_defined, mode_1: by_pitch_up_limit_function">
-								<value>mode_0</value>
-							</method>
+						<aspect_ratio description="parameters for aspect ratio calculation method">
 							<parameters description="Mode parameters">
 								<mode_0 description="user_defined">
 									<aspect_ratio description="user aspect ratio">
@@ -996,10 +1014,7 @@
 								</mode_1>
 							</parameters>
 						</wing_profile_and_thickness_distribution>
-						<mass description="mass calculation methods">
-							<method description="selector mode_0: flops, mode_1: chiozzotto_wer">
-								<value>mode_0</value>
-							</method>
+						<mass description="parameters for mass calculation methods">
 							<parameters>
 								<mode_0 description="flops">
 									<fstrt description="wing strut bracing factor from 0.0 (no strut) to 1.0 (full benefit from strut bracing)">
@@ -1430,10 +1445,7 @@
 						</spars>
 					</front_wing>
 					<back_wing>
-						<sweep description="sweep calculation method">
-							<method description="selector mode_0: user_defined, mode_1: drag_divergence">
-								<value>mode_0</value>
-							</method>
+						<sweep description="parameters for sweep calculation methods">
 							<parameters description="sweep method parameters">
 								<mode_0 description="user_defined">
 									<sweep_angle description="sweep angle at quarter chord">
@@ -1459,10 +1471,7 @@
 								</mode_1>
 							</parameters>
 						</sweep>
-						<taper_ratio description="taper ratio calculation method">
-							<method description="selector mode_0: user_defined, mode_1: howe">
-								<value>mode_0</value>
-							</method>
+						<taper_ratio description="parameters for taper ratio calculation methods">
 							<parameters>
 								<mode_0 description="user_defined">
 									<taper_ratio description="taper ratio of overall wing">
@@ -1474,10 +1483,7 @@
 								</mode_0>
 							</parameters>
 						</taper_ratio>
-						<dihedral description="dihedral calculation method">
-							<method description="selector mode_0: user_defined, mode_1: by_wing_position_and_quarter_chord_sweep">
-								<value>mode_0</value>
-							</method>
+						<dihedral description="parameters for dihedral calculation methods">
 							<parameters description="parameters for methods">
 								<mode_0 description="user_defined">
 									<dihedral_angle description="dihedral angle">
@@ -1494,10 +1500,7 @@
 								</mode_1>
 							</parameters>
 						</dihedral>
-						<aspect_ratio description="aspect ratio calculation method">
-							<method description="selector mode_0: user_defined, mode_1: by_pitch_up_limit_function">
-								<value>mode_0</value>
-							</method>
+						<aspect_ratio description="parameters for aspect ratio calculation methods">
 							<parameters description="Mode parameters">
 								<mode_0 description="user_defined">
 									<aspect_ratio description="user aspect ratio">
@@ -1646,10 +1649,7 @@
 								</mode_1>
 							</parameters>
 						</wing_profile_and_thickness_distribution>
-						<mass description="mass calculation methods">
-							<method description="selector mode_0: flops, mode_1: chiozzotto_wer">
-								<value>mode_0</value>
-							</method>
+						<mass description="parameters for mass calculation methods">
 							<parameters>
 								<mode_0 description="flops">
 									<fstrt description="wing strut bracing factor from 0.0 (no strut) to 1.0 (full benefit from strut bracing)">