Fault calculation implemented

7 files changed, 620 insertions, 68 deletions

diff --git a/.gitignore b/.gitignore
index a46d72b..29383c9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,6 +4,7 @@
 .codelite/
 Project/Release/
 Project/Debug/
+doxygen/
 
 # Compiled source #
 ###################
diff --git a/Project/ElectricCalculation.cpp b/Project/ElectricCalculation.cpp
index c78864c..8f8fa35 100644
--- a/Project/ElectricCalculation.cpp
+++ b/Project/ElectricCalculation.cpp
@@ -37,7 +37,10 @@ void ElectricCalculation::GetElementsFromList(std::vector<Element*> elementList)
     }
 }
 
-bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > >& yBus, double systemPowerBase)
+bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> > >& yBus,
+    double systemPowerBase,
+    YBusSequence sequence,
+    bool includeSyncMachines)
 {
     if(m_busList.size() == 0) return false;
 
@@ -62,8 +65,8 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> >
     }
 
     // Load
-    for(auto itlo = m_loadList.begin(); itlo != m_loadList.end(); itlo++) {
-        Load* load = *itlo;
+    for(auto it = m_loadList.begin(), itEnd = m_loadList.end(); it != itEnd; ++it) {
+        Load* load = *it;
         if(load->IsOnline()) {
             int n = static_cast<Bus*>(load->GetParentList()[0])->GetEletricalData().number;
             LoadElectricalData data = load->GetPUElectricalData(systemPowerBase);
@@ -73,8 +76,8 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> >
     }
 
     // Capacitor
-    for(auto itca = m_capacitorList.begin(); itca != m_capacitorList.end(); itca++) {
-        Capacitor* capacitor = *itca;
+    for(auto it = m_capacitorList.begin(), itEnd = m_capacitorList.end(); it != itEnd; ++it) {
+        Capacitor* capacitor = *it;
         if(capacitor->IsOnline()) {
             int n = static_cast<Bus*>(capacitor->GetParentList()[0])->GetEletricalData().number;
             CapacitorElectricalData data = capacitor->GetPUElectricalData(systemPowerBase);
@@ -83,8 +86,8 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> >
     }
 
     // Inductor
-    for(auto itin = m_inductorList.begin(); itin != m_inductorList.end(); itin++) {
-        Inductor* inductor = *itin;
+    for(auto it = m_inductorList.begin(), itEnd = m_inductorList.end(); it != itEnd; ++it) {
+        Inductor* inductor = *it;
         if(inductor->IsOnline()) {
             int n = static_cast<Bus*>(inductor->GetParentList()[0])->GetEletricalData().number;
             InductorElectricalData data = inductor->GetPUElectricalData(systemPowerBase);
@@ -93,28 +96,43 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> >
     }
 
     // Power line
-    for(auto itl = m_lineList.begin(); itl != m_lineList.end(); itl++) {
-        Line* line = *itl;
+    for(auto it = m_lineList.begin(), itEnd = m_lineList.end(); it != itEnd; ++it) {
+        Line* line = *it;
         if(line->IsOnline()) {
             LineElectricalData data = line->GetElectricalData();
 
             int n1 = static_cast<Bus*>(line->GetParentList()[0])->GetEletricalData().number;
             int n2 = static_cast<Bus*>(line->GetParentList()[1])->GetEletricalData().number;
 
-            yBus[n1][n2] -= 1.0 / std::complex<double>(data.resistance, data.indReactance);
-            yBus[n2][n1] -= 1.0 / std::complex<double>(data.resistance, data.indReactance);
+            switch(sequence) {
+                case POSITIVE_SEQ:
+                case NEGATIVE_SEQ: {
+                    yBus[n1][n2] -= 1.0 / std::complex<double>(data.resistance, data.indReactance);
+                    yBus[n2][n1] -= 1.0 / std::complex<double>(data.resistance, data.indReactance);
 
-            yBus[n1][n1] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
-            yBus[n2][n2] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
+                    yBus[n1][n1] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
+                    yBus[n2][n2] += 1.0 / std::complex<double>(data.resistance, data.indReactance);
 
-            yBus[n1][n1] += std::complex<double>(0.0, data.capSusceptance / 2.0);
-            yBus[n2][n2] += std::complex<double>(0.0, data.capSusceptance / 2.0);
+                    yBus[n1][n1] += std::complex<double>(0.0, data.capSusceptance / 2.0);
+                    yBus[n2][n2] += std::complex<double>(0.0, data.capSusceptance / 2.0);
+                } break;
+                case ZERO_SEQ: {
+                    yBus[n1][n2] -= 1.0 / std::complex<double>(data.zeroResistance, data.zeroIndReactance);
+                    yBus[n2][n1] -= 1.0 / std::complex<double>(data.zeroResistance, data.zeroIndReactance);
+
+                    yBus[n1][n1] += 1.0 / std::complex<double>(data.zeroResistance, data.zeroIndReactance);
+                    yBus[n2][n2] += 1.0 / std::complex<double>(data.zeroResistance, data.zeroIndReactance);
+
+                    yBus[n1][n1] += std::complex<double>(0.0, data.zeroCapSusceptance / 2.0);
+                    yBus[n2][n2] += std::complex<double>(0.0, data.zeroCapSusceptance / 2.0);
+                }
+            }
         }
     }
 
     // Transformer
-    for(auto itt = m_transformerList.begin(); itt != m_transformerList.end(); ++itt) {
-        Transformer* transformer = *itt;
+    for(auto it = m_transformerList.begin(), itEnd = m_transformerList.end(); it != itEnd; ++it) {
+        Transformer* transformer = *it;
 
         if(transformer->IsOnline()) {
             TransformerElectricalData data = transformer->GetElectricalData();
@@ -124,7 +142,7 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> >
 
             // If the transformer have nominal turns ratio (1.0) and no phase shifting, it will be modelled as series
             // impedance.
-            if(data.turnsRatio == 1.0 && data.phaseShift == 0.0) {
+            if(data.turnsRatio == 1.0 && data.phaseShift == 0.0 && sequence != ZERO_SEQ) {
                 yBus[n1][n2] += -1.0 / std::complex<double>(data.resistance, data.indReactance);
                 yBus[n2][n1] += -1.0 / std::complex<double>(data.resistance, data.indReactance);
 
@@ -136,7 +154,7 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> >
             //[Ref. 1: Elementos de analise de sistemas de potencia - Stevenson - pg. 232]
             //[Ref. 2: http://www.ee.washington.edu/research/real/Library/Reports/Tap_Adjustments_in_AC_Load_Flows.pdf]
             //[Ref. 3: http://www.columbia.edu/~dano/courses/power/notes/power/andersson1.pdf]
-            else {
+            else if(sequence != ZERO_SEQ) {
                 // Complex turns ratio
                 double radPhaseShift = wxDegToRad(data.phaseShift);
                 std::complex<double> a = std::complex<double>(
@@ -145,10 +163,90 @@ bool ElectricCalculation::GetYBus(std::vector<std::vector<std::complex<double> >
                 // Transformer admitance
                 std::complex<double> y = 1.0 / std::complex<double>(data.resistance, data.indReactance);
 
-                yBus[n1][n1] += y / std::pow(std::abs(a), 2.0);
-                yBus[n1][n2] += -(y / std::conj(a));
-                yBus[n2][n1] += -(y / a);
-                yBus[n2][n2] += y;
+                if(sequence == POSITIVE_SEQ) {
+                    yBus[n1][n1] += y / std::pow(std::abs(a), 2.0);
+                    yBus[n1][n2] += -(y / std::conj(a));
+                    yBus[n2][n1] += -(y / a);
+                    yBus[n2][n2] += y;
+                } else if(sequence == NEGATIVE_SEQ) {
+                    yBus[n1][n1] += y / std::pow(std::abs(a), 2.0);
+                    yBus[n1][n2] += -(y / a);
+                    yBus[n2][n1] += -(y / std::conj(a));
+                    yBus[n2][n2] += y;
+                }
+            } else if(sequence == ZERO_SEQ) {
+                switch(data.connection) {
+                    case GWYE_GWYE: {
+                        std::complex<double> y =
+                            1.0 / std::complex<double>(data.zeroResistance +
+                                          3.0 * (data.primaryGrndResistance + data.secondaryGrndResistance),
+                                      data.zeroIndReactance +
+                                          3.0 * (data.primaryGrndReactance + data.secondaryGrndReactance));
+                        std::complex<double> a = std::complex<double>(data.turnsRatio, 0.0);
+
+                        yBus[n1][n1] += y / (a * a);
+                        yBus[n1][n2] += -(y / a);
+                        yBus[n2][n1] += -(y / a);
+                        yBus[n2][n2] += y;
+                    } break;
+                    case DELTA_GWYE: {
+                        std::complex<double> y =
+                            1.0 / std::complex<double>(data.zeroResistance + 3.0 * (data.secondaryGrndResistance),
+                                      data.zeroIndReactance + 3.0 * (data.secondaryGrndReactance));
+                        yBus[n2][n2] += y;
+                        break;
+                    }
+                    case GWYE_DELTA: {
+                        std::complex<double> y =
+                            1.0 / std::complex<double>(data.zeroResistance + 3.0 * (data.primaryGrndResistance),
+                                      data.zeroIndReactance + 3.0 * (data.primaryGrndReactance));
+                        yBus[n1][n1] += y;
+                        break;
+                    }
+                    default:
+                        break;
+                }
+            }
+        }
+    }
+
+    if(includeSyncMachines) {
+        // Synchronous generator
+        for(auto it = m_syncGeneratorList.begin(), itEnd = m_syncGeneratorList.end(); it != itEnd; ++it) {
+            SyncGenerator* syncGenerator = *it;
+            if(syncGenerator->IsOnline()) {
+                int n = static_cast<Bus*>(syncGenerator->GetParentList()[0])->GetEletricalData().number;
+                SyncGeneratorElectricalData data = syncGenerator->GetPUElectricalData(systemPowerBase);
+                switch(sequence) {
+                    case POSITIVE_SEQ: {
+                        yBus[n][n] += 1.0 / std::complex<double>(data.positiveResistance, data.positiveReactance);
+                    } break;
+                    case NEGATIVE_SEQ: {
+                        yBus[n][n] += 1.0 / std::complex<double>(data.negativeResistance, data.negativeReactance);
+                    } break;
+                    case ZERO_SEQ: {
+                        yBus[n][n] += 1.0 / std::complex<double>(data.zeroResistance, data.zeroReactance);
+                    } break;
+                }
+            }
+        }
+        // Synchronous motor
+        for(auto it = m_syncMotorList.begin(), itEnd = m_syncMotorList.end(); it != itEnd; ++it) {
+            SyncMotor* syncMotor = *it;
+            if(syncMotor->IsOnline()) {
+                int n = static_cast<Bus*>(syncMotor->GetParentList()[0])->GetEletricalData().number;
+                SyncMotorElectricalData data = syncMotor->GetPUElectricalData(systemPowerBase);
+                switch(sequence) {
+                    case POSITIVE_SEQ: {
+                        yBus[n][n] += 1.0 / std::complex<double>(data.positiveResistance, data.positiveReactance);
+                    } break;
+                    case NEGATIVE_SEQ: {
+                        yBus[n][n] += 1.0 / std::complex<double>(data.negativeResistance, data.negativeReactance);
+                    } break;
+                    case ZERO_SEQ: {
+                        yBus[n][n] += 1.0 / std::complex<double>(data.zeroResistance, data.zeroReactance);
+                    } break;
+                }
             }
         }
     }
@@ -429,3 +527,72 @@ void ElectricCalculation::UpdateElementsPowerFlow(std::vector<std::complex<doubl
         }
     }
 }
+
+bool ElectricCalculation::InvertMatrix(std::vector<std::vector<std::complex<double> > > matrix,
+    std::vector<std::vector<std::complex<double> > >& inverse)
+{
+    int order = static_cast<int>(matrix.size());
+
+    inverse.clear();
+    // Fill the inverse matrix with identity.
+    for(int i = 0; i < order; ++i) {
+        std::vector<std::complex<double> > line;
+        for(int j = 0; j < order; ++j) {
+            line.push_back(i == j ? std::complex<double>(1.0, 0.0) : std::complex<double>(0.0, 0.0));
+        }
+        inverse.push_back(line);
+    }
+
+    // Check if a main diagonal value of the matrix is zero, if one is zero, try a linear combination to remove it.
+    for(int i = 0; i < order; ++i) {
+        for(int j = 0; j < order; ++j) {
+            if(i == j && matrix[i][j] == std::complex<double>(0.0, 0.0)) {
+                int row = 0;
+                while(row < order) {
+                    if(matrix[row][j] != std::complex<double>(0.0, 0.0)) {
+                        for(int k = 0; k < order; ++k) {
+                            matrix[i][k] += matrix[row][k];
+                            inverse[i][k] += inverse[row][k];
+                        }
+                        break;
+                    }
+                    row++;
+                }
+                // If all line values are zero, the matrix is singular and the solution is impossible.
+                if(row == order) return false;
+            }
+        }
+    }
+
+    // Linear combinations are made in both matrices, the goal is the input matrix become the identity. The final result
+    // have two matrices: the identity and the inverse of the input.
+    for(int i = 0; i < order; ++i) {
+        for(int j = 0; j < order; ++j) {
+            if(i != j) {
+                if(matrix[i][i] == std::complex<double>(0.0, 0.0)) return false;
+
+                std::complex<double> factor = matrix[j][i] / matrix[i][i];
+                for(int k = 0; k < order; ++k) {
+                    matrix[j][k] -= factor * matrix[i][k];
+                    inverse[j][k] -= factor * inverse[i][k];
+                }
+            }
+        }
+    }
+    // Main diagonal calculation.
+    for(int i = 0; i < order; ++i) {
+        for(int j = 0; j < order; ++j) {
+            if(i == j) {
+                if(matrix[i][j] == std::complex<double>(0.0, 0.0)) return false;
+
+                std::complex<double> factor = (matrix[i][j] - std::complex<double>(1.0, 0.0)) / matrix[i][j];
+                for(int k = 0; k < order; ++k) {
+                    matrix[j][k] -= factor * matrix[i][k];
+                    inverse[j][k] -= factor * inverse[i][k];
+                }
+            }
+        }
+    }
+
+    return true;
+}
diff --git a/Project/ElectricCalculation.h b/Project/ElectricCalculation.h
index e4e4ef7..3eea5dd 100644
--- a/Project/ElectricCalculation.h
+++ b/Project/ElectricCalculation.h
@@ -18,14 +18,16 @@
 enum BusType { BUS_SLACK = 0, BUS_PV, BUS_PQ };
 
 enum ReactiveLimitsType {
-    RL_UNLIMITED = 0,     // The bus can generate any ammount of reactive power.
-    RL_LIMITED,           // The bus reactive power generation is limited.
-    RL_UNLIMITED_SOURCE,  // The bus have at least one source of infinite reative power.
-    RL_MAX_REACHED,       // Max limit reached
-    RL_MIN_REACHED,       // Min limit reached
-    RL_NONE_REACHED       // No limits reached
+    RL_UNLIMITED = 0,    // The bus can generate any ammount of reactive power.
+    RL_LIMITED,          // The bus reactive power generation is limited.
+    RL_UNLIMITED_SOURCE, // The bus have at least one source of infinite reative power.
+    RL_MAX_REACHED,      // Max limit reached
+    RL_MIN_REACHED,      // Min limit reached
+    RL_NONE_REACHED      // No limits reached
 };
 
+enum YBusSequence { POSITIVE_SEQ = 0, NEGATIVE_SEQ, ZERO_SEQ };
+
 struct ReactiveLimits {
     double maxLimit = 0.0;
     double minLimit = 0.0;
@@ -34,30 +36,123 @@ struct ReactiveLimits {
     ReactiveLimitsType limitReached = RL_NONE_REACHED;
 };
 
+/**
+ * @class ElectricCalculation
+ * @author Thales Lima Oliveira
+ * @date 09/01/2017
+ * @file ElectricCalculation.h
+ * @brief Base class of electric calculations, with general methods.
+ */
 class ElectricCalculation
 {
-   public:
+public:
+    /**
+     * @brief Constructor.
+     */
     ElectricCalculation();
+
+    /**
+     * @brief Destructor.
+     */
     ~ElectricCalculation();
+
+    /**
+     * @brief Separate the elements from a generic list.
+     * @param elementList List of generic elements.
+     */
     virtual void GetElementsFromList(std::vector<Element*> elementList);
-    virtual bool GetYBus(std::vector<std::vector<std::complex<double> > >& yBus, double systemPowerBase);
+
+    /**
+     * @brief Get the admittance matrix from the list of elements (use GetElementsFromList first).
+     * @param yBus Admittance matrix. The previous content will be erased.
+     * @param systemPowerBase Base power of the system.
+     * @param sequence Sequence of admittance matrix (positive, negative and zero).
+     * @param includeSyncMachines Include the synchronous machines on calculation.
+     * @return Return true if was possible to build the admittance matrix.
+     */
+    virtual bool GetYBus(std::vector<std::vector<std::complex<double> > >& yBus,
+        double systemPowerBase,
+        YBusSequence sequence = POSITIVE_SEQ,
+        bool includeSyncMachines = false);
+
+    /**
+     * @brief Invert a matrix.
+     * @param matrix Matrix to invert.
+     * @param inverse Inverted matrix. The previous content will be erased.
+     * @return Return true if was possible to invert the matrix.
+     */
+    virtual bool InvertMatrix(std::vector<std::vector<std::complex<double> > > matrix,
+        std::vector<std::vector<std::complex<double> > >& inverse);
+
+    /**
+     * @brief Update the elements after the power flow calculation.
+     * @param voltage Array with the buses voltages.
+     * @param power Array with the buses injected power.
+     * @param busType Array with the buses type.
+     * @param reactiveLimit Array with the reactive limit data.
+     * @param systemPowerBase Base power of the system.
+     */
     virtual void UpdateElementsPowerFlow(std::vector<std::complex<double> > voltage,
-                                         std::vector<std::complex<double> > power,
-                                         std::vector<BusType> busType,
-                                         std::vector<ReactiveLimits> reactiveLimit,
-                                         double systemPowerBase);
+        std::vector<std::complex<double> > power,
+        std::vector<BusType> busType,
+        std::vector<ReactiveLimits> reactiveLimit,
+        double systemPowerBase);
 
+    /**
+     * @brief Get the buses of the system (use GetElementsFromList first).
+     * @return A list of bus elements.
+     */
     const std::vector<Bus*> GetBusList() const { return m_busList; }
+
+    /**
+     * @brief Get the capacitors of the system (use GetElementsFromList first).
+     * @return A list of capacitor elements.
+     */
     const std::vector<Capacitor*> GetCapacitorList() const { return m_capacitorList; }
+    
+    /**
+     * @brief Get the induction motors of the system (use GetElementsFromList first).
+     * @return A list of induction motor elements.
+     */
     const std::vector<IndMotor*> GetIndMotorList() const { return m_indMotorList; }
+    
+    /**
+     * @brief Get the inductors of the system (use GetElementsFromList first).
+     * @return A list of inductor elements.
+     */
     const std::vector<Inductor*> GetInductorList() const { return m_inductorList; }
+    
+    /**
+     * @brief Get the lines of the system (use GetElementsFromList first).
+     * @return A list of line elements.
+     */
     const std::vector<Line*> GetLineList() const { return m_lineList; }
+    
+    /**
+     * @brief Get the loads of the system (use GetElementsFromList first).
+     * @return A list of load elements.
+     */
     const std::vector<Load*> GetLoadList() const { return m_loadList; }
+    
+    /**
+     * @brief Get the synchronous generators of the system (use GetElementsFromList first).
+     * @return A list of synchronous generator elements.
+     */
     const std::vector<SyncGenerator*> GetSyncGeneratorList() const { return m_syncGeneratorList; }
+    
+    /**
+     * @brief Get the synchronous motors of the system (use GetElementsFromList first).
+     * @return A list of synchronous motor elements.
+     */
     const std::vector<SyncMotor*> GetSyncMotorList() const { return m_syncMotorList; }
-    const std::vector<Transformer*> GetTransformerList() const { return m_transformerList; }
     
-   protected:
+    /**
+     * @brief Get the transformers of the system (use GetElementsFromList first).
+     * @return A list of transformer elements.
+     */
+    const std::vector<Transformer*> GetTransformerList() const { return m_transformerList; }
+
+protected:
     std::vector<Bus*> m_busList;
     std::vector<Capacitor*> m_capacitorList;
     std::vector<IndMotor*> m_indMotorList;
@@ -69,4 +164,4 @@ class ElectricCalculation
     std::vector<Transformer*> m_transformerList;
 };
 
-#endif  // ELECTRICCALCULATION_H
+#endif // ELECTRICCALCULATION_H
diff --git a/Project/Fault.cpp b/Project/Fault.cpp
new file mode 100644
index 0000000..058b4ea
--- /dev/null
+++ b/Project/Fault.cpp
@@ -0,0 +1,238 @@
+#include "Fault.h"
+
+Fault::Fault()
+    : ElectricCalculation()
+{
+}
+
+Fault::Fault(std::vector<Element*> elementList) { GetElementsFromList(elementList); }
+
+Fault::~Fault() {}
+
+bool Fault::RunFaultCalculation(double systemPowerBase)
+{
+    int numberOfBuses = static_cast<int>(m_busList.size());
+    if(numberOfBuses == 0) {
+        m_errorMsg = _("There is no buses in the system.");
+        return false;
+    }
+
+    // Get adimittance matrices.
+    std::vector<std::vector<std::complex<double> > > yBusPos;
+    GetYBus(yBusPos, systemPowerBase, POSITIVE_SEQ, true);
+    std::vector<std::vector<std::complex<double> > > yBusNeg;
+    GetYBus(yBusNeg, systemPowerBase, NEGATIVE_SEQ, true);
+    std::vector<std::vector<std::complex<double> > > yBusZero;
+    GetYBus(yBusZero, systemPowerBase, ZERO_SEQ, true);
+
+    // Calculate the impedance matrices.
+    if(!InvertMatrix(yBusPos, m_zBusPos)) {
+        m_errorMsg = _("Fail to invert the positive sequence admittance matrix.");
+        return false;
+    }
+    if(!InvertMatrix(yBusNeg, m_zBusNeg)) {
+        m_errorMsg = _("Fail to invert the negative sequence admittance matrix.");
+        return false;
+    }
+    if(!InvertMatrix(yBusZero, m_zBusZero)) {
+        m_errorMsg = _("Fail to invert the zero sequence admittance matrix.");
+        return false;
+    }
+
+    // Pre-fault voltages (power flow solution).
+    std::vector<std::complex<double> > preFaultVoltages;
+
+    // Get fault parameters.
+    int fNumber = -1;
+    FaultData fType = FAULT_THREEPHASE;
+    FaultData fLocation = FAULT_LINE_A;
+    std::complex<double> fImpedance = std::complex<double>(0.0, 0.0);
+    for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) {
+        Bus* bus = *it;
+        BusElectricalData data = bus->GetEletricalData();
+        preFaultVoltages.push_back(data.voltage);
+        if(data.hasFault) {
+            fNumber = data.number;
+            fType = data.faultType;
+            fLocation = data.faultLocation;
+            fImpedance = std::complex<double>(data.faultResistance, data.faultReactance);
+        }
+    }
+
+    if(fNumber == -1) {
+        m_errorMsg = _("There is no fault in the system.");
+        return false;
+    }
+
+    // Fault calculation.
+    std::complex<double> fCurrentPos = std::complex<double>(0.0, 0.0);
+    std::complex<double> fCurrentNeg = std::complex<double>(0.0, 0.0);
+    std::complex<double> fCurrentZero = std::complex<double>(0.0, 0.0);
+
+    std::complex<double> preFaultVoltage = preFaultVoltages[fNumber];
+    std::complex<double> a = std::complex<double>(-0.5, 0.866025403784);
+    std::complex<double> a2 = std::complex<double>(-0.5, -0.866025403784);
+
+    switch(fType) {
+        case FAULT_THREEPHASE: {
+            fCurrentPos = preFaultVoltage / (m_zBusPos[fNumber][fNumber] + fImpedance);
+        } break;
+        case FAULT_2LINE: {
+            fCurrentPos = preFaultVoltage / (m_zBusPos[fNumber][fNumber] + m_zBusNeg[fNumber][fNumber] + fImpedance);
+
+            switch(fLocation) {
+                case FAULT_LINE_A: {
+                    fCurrentNeg = -a2 * fCurrentPos;
+                } break;
+                case FAULT_LINE_B: {
+                    fCurrentNeg = -fCurrentPos;
+                } break;
+                case FAULT_LINE_C: {
+                    fCurrentNeg = -a * fCurrentPos;
+                } break;
+                default:
+                    break;
+            }
+        } break;
+        case FAULT_2LINE_GROUND: {
+            std::complex<double> z1 = m_zBusPos[fNumber][fNumber];
+            std::complex<double> z2 = m_zBusNeg[fNumber][fNumber];
+            std::complex<double> z0 = m_zBusZero[fNumber][fNumber];
+            std::complex<double> zf_3 = std::complex<double>(3.0, 0.0) * fImpedance;
+
+            fCurrentPos = (preFaultVoltage * (z2 + z0 + zf_3)) / (z1 * z2 + z2 * z0 + z2 * zf_3 + z1 * z0 + z1 * zf_3);
+
+            switch(fLocation) {
+                case FAULT_LINE_A: {
+                    fCurrentNeg = -a2 * ((preFaultVoltage - z1 * fCurrentPos) / z2);
+                    fCurrentZero = -a * ((preFaultVoltage - z1 * fCurrentPos) / (z0 + zf_3));
+                } break;
+                case FAULT_LINE_B: {
+                    fCurrentNeg = -((preFaultVoltage - z1 * fCurrentPos) / z2);
+                    fCurrentZero = -((preFaultVoltage - z1 * fCurrentPos) / (z0 + zf_3));
+                } break;
+                case FAULT_LINE_C: {
+                    fCurrentNeg = -a * ((preFaultVoltage - z1 * fCurrentPos) / z2);
+                    fCurrentZero = -a2 * ((preFaultVoltage - z1 * fCurrentPos) / (z0 + zf_3));
+                } break;
+                default:
+                    break;
+            }
+        } break;
+        case FAULT_LINE_GROUND: {
+            fCurrentPos =
+                preFaultVoltage / (m_zBusPos[fNumber][fNumber] + m_zBusNeg[fNumber][fNumber] +
+                                      m_zBusZero[fNumber][fNumber] + std::complex<double>(3.0, 0.0) * fImpedance);
+            switch(fLocation) {
+                case FAULT_LINE_A: {
+                    fCurrentNeg = fCurrentPos;
+                    fCurrentZero = fCurrentPos;
+                } break;
+                case FAULT_LINE_B: {
+                    fCurrentNeg = a * fCurrentPos;
+                    fCurrentZero = a2 * fCurrentPos;
+                } break;
+                case FAULT_LINE_C: {
+                    fCurrentNeg = a2 * fCurrentPos;
+                    fCurrentZero = a * fCurrentPos;
+                } break;
+                default:
+                    break;
+            }
+        } break;
+        default:
+            break;
+    }
+
+    // Convert sequence currents to ABC. [Iabc] = [A]*[I012]
+    m_fCurrentA = fCurrentPos + fCurrentNeg + fCurrentZero;
+    m_fCurrentB = fCurrentPos + a2 * fCurrentNeg + a * fCurrentZero;
+    m_fCurrentC = fCurrentPos + a * fCurrentNeg + a2 * fCurrentZero;
+
+    // Pos-fault voltages calculation
+    m_posFaultVoltagePos.clear();
+    m_posFaultVoltageNeg.clear();
+    m_posFaultVoltageZero.clear();
+    m_posFaultVoltageA.clear();
+    m_posFaultVoltageB.clear();
+    m_posFaultVoltageC.clear();
+
+    for(int i = 0; i < numberOfBuses; ++i) {
+        m_posFaultVoltagePos.push_back(preFaultVoltages[i] - m_zBusPos[i][fNumber] * fCurrentPos);
+        m_posFaultVoltageNeg.push_back(-m_zBusNeg[i][fNumber] * fCurrentNeg);
+        m_posFaultVoltageZero.push_back(-m_zBusZero[i][fNumber] * fCurrentZero);
+
+        // V012 -> Vabc
+        m_posFaultVoltageA.push_back(m_posFaultVoltagePos[i] + m_posFaultVoltageNeg[i] + m_posFaultVoltageZero[i]);
+        m_posFaultVoltageB.push_back(
+            m_posFaultVoltagePos[i] + a2 * m_posFaultVoltageNeg[i] + a * m_posFaultVoltageZero[i]);
+        m_posFaultVoltageC.push_back(
+            m_posFaultVoltagePos[i] + a * m_posFaultVoltageNeg[i] + a2 * m_posFaultVoltageZero[i]);
+    }
+
+    UpdateElementsFault(systemPowerBase);
+    return true;
+}
+
+void Fault::UpdateElementsFault(double systemPowerBase)
+{
+    std::complex<double> a = std::complex<double>(-0.5, 0.866025403784);
+    std::complex<double> a2 = std::complex<double>(-0.5, -0.866025403784);
+
+    for(auto it = m_busList.begin(), itEnd = m_busList.end(); it != itEnd; ++it) {
+        Bus* bus = *it;
+        auto data = bus->GetEletricalData();
+        if(data.hasFault) {
+            data.faultCurrent[0] = m_fCurrentA;
+            data.faultCurrent[1] = m_fCurrentB;
+            data.faultCurrent[2] = m_fCurrentC;
+        } else {
+            data.faultCurrent[0] = data.faultCurrent[1] = data.faultCurrent[2] = std::complex<double>(0.0, 0.0);
+        }
+        data.faultVoltage[0] = m_posFaultVoltageA[data.number];
+        data.faultVoltage[1] = m_posFaultVoltageB[data.number];
+        data.faultVoltage[2] = m_posFaultVoltageC[data.number];
+        bus->SetElectricalData(data);
+    }
+
+    for(auto it = m_lineList.begin(), itEnd = m_lineList.end(); it != itEnd; ++it) {
+        Line* line = *it;
+        if(line->IsOnline()) {
+            int n1 = static_cast<Bus*>(line->GetParentList()[0])->GetEletricalData().number;
+            int n2 = static_cast<Bus*>(line->GetParentList()[1])->GetEletricalData().number;
+            auto data = line->GetElectricalData();
+            std::complex<double> vPos[2] = { m_posFaultVoltagePos[n1], m_posFaultVoltagePos[n2] };
+            std::complex<double> vNeg[2] = { m_posFaultVoltageNeg[n1], m_posFaultVoltageNeg[n2] };
+            std::complex<double> vZero[2] = { m_posFaultVoltageZero[n1], m_posFaultVoltageZero[n2] };
+            std::complex<double> zPos(data.resistance, data.indReactance);
+            std::complex<double> bPos(0.0, data.capSusceptance / 2.0);
+            std::complex<double> zZero(data.zeroResistance, data.zeroIndReactance);
+            std::complex<double> bZero(0.0, data.zeroCapSusceptance / 2.0);
+
+            std::complex<double> lineCurrentPos[2];
+            std::complex<double> lineCurrentNeg[2];
+            std::complex<double> lineCurrentZero[2];
+
+            lineCurrentPos[0] = ((vPos[0] - vPos[1]) / zPos) + (vPos[0] * bPos);
+            lineCurrentNeg[0] = ((vNeg[0] - vNeg[1]) / zPos) + (vNeg[0] * bPos);
+            lineCurrentZero[0] = ((vZero[0] - vZero[1]) / zZero) + (vZero[0] * bZero);
+            lineCurrentPos[1] = ((vPos[1] - vPos[0]) / zPos) + (vPos[1] * bPos);
+            lineCurrentNeg[1] = ((vNeg[1] - vNeg[0]) / zPos) + (vNeg[1] * bPos);
+            lineCurrentZero[1] = ((vZero[1] - vZero[0]) / zZero) + (vZero[1] * bZero);
+
+            data.faultCurrent[0][0] = lineCurrentPos[0] + lineCurrentNeg[0] + lineCurrentZero[0];
+            data.faultCurrent[1][0] = lineCurrentPos[0] + a2 * lineCurrentNeg[0] + a * lineCurrentZero[0];
+            data.faultCurrent[2][0] = lineCurrentPos[0] + a * lineCurrentNeg[0] + a2 * lineCurrentZero[0];
+            data.faultCurrent[0][1] = lineCurrentPos[1] + lineCurrentNeg[1] + lineCurrentZero[1];
+            data.faultCurrent[1][1] = lineCurrentPos[1] + a2 * lineCurrentNeg[1] + a * lineCurrentZero[1];
+            data.faultCurrent[2][1] = lineCurrentPos[1] + a * lineCurrentNeg[1] + a2 * lineCurrentZero[1];
+
+            line->SetElectricalData(data);
+        }
+    }
+    
+    for(auto it = m_transformerList.begin(), itEnd = m_transformerList.end(); it != itEnd; ++it) {
+        Transformer* transformer = *it;
+        
+    }
+}
diff --git a/Project/Fault.h b/Project/Fault.h
new file mode 100644
index 0000000..dfe0e75
--- /dev/null
+++ b/Project/Fault.h
@@ -0,0 +1,37 @@
+#ifndef FAULT_H
+#define FAULT_H
+
+#include "ElectricCalculation.h"
+
+class Fault : public ElectricCalculation
+{
+public:
+    Fault(std::vector<Element*> elementList);
+    Fault();
+    ~Fault();
+
+    virtual bool RunFaultCalculation(double systemPowerBase);
+    virtual void UpdateElementsFault(double systemPowerBase);
+    virtual wxString GetErrorMessage() { return m_errorMsg; }
+
+protected:
+    wxString m_errorMsg = "";
+
+    std::vector<std::vector<std::complex<double> > > m_zBusPos;
+    std::vector<std::vector<std::complex<double> > > m_zBusNeg;
+    std::vector<std::vector<std::complex<double> > > m_zBusZero;
+
+    std::vector<std::complex<double> > m_posFaultVoltagePos;
+    std::vector<std::complex<double> > m_posFaultVoltageNeg;
+    std::vector<std::complex<double> > m_posFaultVoltageZero;
+
+    std::complex<double> m_fCurrentA;
+    std::complex<double> m_fCurrentB;
+    std::complex<double> m_fCurrentC;
+
+    std::vector<std::complex<double> > m_posFaultVoltageA;
+    std::vector<std::complex<double> > m_posFaultVoltageB;
+    std::vector<std::complex<double> > m_posFaultVoltageC;
+};
+
+#endif // FAULT_H
diff --git a/Project/PowerFlow.cpp b/Project/PowerFlow.cpp
index c94d8d8..bfee2af 100644
--- a/Project/PowerFlow.cpp
+++ b/Project/PowerFlow.cpp
@@ -1,12 +1,23 @@
 #include "PowerFlow.h"
 
-PowerFlow::PowerFlow(std::vector<Element*> elementList) : ElectricCalculation() { GetElementsFromList(elementList); }
+PowerFlow::PowerFlow()
+    : ElectricCalculation()
+{
+}
+
+PowerFlow::PowerFlow(std::vector<Element*> elementList)
+    : ElectricCalculation()
+{
+    GetElementsFromList(elementList);
+}
+
 PowerFlow::~PowerFlow() {}
+
 bool PowerFlow::RunGaussSeidel(double systemPowerBase,
-                               int maxIteration,
-                               double error,
-                               double initAngle,
-                               double accFactor)
+    int maxIteration,
+    double error,
+    double initAngle,
+    double accFactor)
 {
     // Calculate the Ybus.
     if(!GetYBus(m_yBus, systemPowerBase)) {
@@ -15,13 +26,13 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase,
     }
 
     // Number of buses on the system.
-    int numberOfBuses = (int)m_busList.size();
+    int numberOfBuses = static_cast<int>(m_busList.size());
 
-    std::vector<BusType> busType;                  // Bus type
-    std::vector<std::complex<double> > voltage;    // Voltage of buses
-    std::vector<std::complex<double> > power;      // Injected power
-    std::vector<std::complex<double> > loadPower;  // Only the load power
-    std::vector<ReactiveLimits> reactiveLimit;     // Limit of reactive power on PV buses
+    std::vector<BusType> busType;                 // Bus type
+    std::vector<std::complex<double> > voltage;   // Voltage of buses
+    std::vector<std::complex<double> > power;     // Injected power
+    std::vector<std::complex<double> > loadPower; // Only the load power
+    std::vector<ReactiveLimits> reactiveLimit;    // Limit of reactive power on PV buses
 
     reactiveLimit.resize(numberOfBuses);
 
@@ -61,7 +72,7 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase,
         }
 
         // Fill the power array
-        power.push_back(std::complex<double>(0.0, 0.0));  // Initial value
+        power.push_back(std::complex<double>(0.0, 0.0)); // Initial value
         loadPower.push_back(std::complex<double>(0.0, 0.0));
 
         // Synchronous generator
@@ -168,12 +179,12 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase,
     }
 
     // Gauss-Seidel method
-    std::vector<std::complex<double> > oldVoltage;  // Old voltage array.
+    std::vector<std::complex<double> > oldVoltage; // Old voltage array.
     oldVoltage.resize(voltage.size());
 
     auto oldBusType = busType;
 
-    int iteration = 0;  // Current itaration number.
+    int iteration = 0; // Current itaration number.
 
     while(true) {
         // Reach the max number of iterations.
@@ -203,7 +214,7 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase,
 
                 // Apply the acceleration factor.
                 newVolt = std::complex<double>(accFactor * (newVolt.real() - voltage[i].real()) + voltage[i].real(),
-                                               accFactor * (newVolt.imag() - voltage[i].imag()) + voltage[i].imag());
+                    accFactor * (newVolt.imag() - voltage[i].imag()) + voltage[i].imag());
 
                 voltage[i] = newVolt;
             }
@@ -227,15 +238,15 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase,
 
                 // Apply the acceleration factor.
                 newVolt = std::complex<double>(accFactor * (newVolt.real() - voltage[i].real()) + voltage[i].real(),
-                                               accFactor * (newVolt.imag() - voltage[i].imag()) + voltage[i].imag());
+                    accFactor * (newVolt.imag() - voltage[i].imag()) + voltage[i].imag());
 
                 // Keep the same voltage magnitude.
                 voltage[i] = std::complex<double>(std::abs(voltage[i]) * std::cos(std::arg(newVolt)),
-                                                  std::abs(voltage[i]) * std::sin(std::arg(newVolt)));
+                    std::abs(voltage[i]) * std::sin(std::arg(newVolt)));
             }
 
-            double busError = std::max(std::abs(voltage[i].real() - oldVoltage[i].real()),
-                                       std::abs(voltage[i].imag() - oldVoltage[i].imag()));
+            double busError = std::max(
+                std::abs(voltage[i].real() - oldVoltage[i].real()), std::abs(voltage[i].imag() - oldVoltage[i].imag()));
 
             if(busError > iterationError) iterationError = busError;
         }
@@ -246,7 +257,8 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase,
                 if(busType[i] == BUS_PV) {
                     if(reactiveLimit[i].maxLimitType == RL_LIMITED) {
                         if(power[i].imag() - loadPower[i].imag() > reactiveLimit[i].maxLimit) {
-                            power[i] = std::complex<double>(power[i].real(), reactiveLimit[i].maxLimit + loadPower[i].imag());
+                            power[i] =
+                                std::complex<double>(power[i].real(), reactiveLimit[i].maxLimit + loadPower[i].imag());
                             busType[i] = BUS_PQ;
                             reactiveLimit[i].limitReached = RL_MAX_REACHED;
                             limitReach = true;
@@ -254,7 +266,8 @@ bool PowerFlow::RunGaussSeidel(double systemPowerBase,
                     }
                     if(reactiveLimit[i].minLimitType == RL_LIMITED) {
                         if(power[i].imag() - loadPower[i].imag() < reactiveLimit[i].minLimit) {
-                            power[i] = std::complex<double>(power[i].real(), reactiveLimit[i].minLimit + loadPower[i].imag());
+                            power[i] =
+                                std::complex<double>(power[i].real(), reactiveLimit[i].minLimit + loadPower[i].imag());
                             busType[i] = BUS_PQ;
                             reactiveLimit[i].limitReached = RL_MIN_REACHED;
                             limitReach = true;
diff --git a/Project/PowerFlow.h b/Project/PowerFlow.h
index 1e5c621..29e205a 100644
--- a/Project/PowerFlow.h
+++ b/Project/PowerFlow.h
@@ -4,24 +4,25 @@
 #include "ElectricCalculation.h"
 
 #include <wx/string.h>
-#include <wx/intl.h>  //_()
-#include <wx/log.h>   //temp
+#include <wx/intl.h> //_()
 
 class PowerFlow : public ElectricCalculation
 {
-   public:
+public:
+    PowerFlow();
     PowerFlow(std::vector<Element*> elementList);
     ~PowerFlow();
     virtual bool RunGaussSeidel(double systemPowerBase = 100e6,
-                                int maxIteration = 5000,
-                                double error = 1e-6,
-                                double initAngle = 0.0,
-                                double accFactor = 1.0);
+        int maxIteration = 5000,
+        double error = 1e-6,
+        double initAngle = 0.0,
+        double accFactor = 1.0);
 
     virtual wxString GetErrorMessage() { return m_errorMsg; }
-   protected:
+
+protected:
     std::vector<std::vector<std::complex<double> > > m_yBus;
     wxString m_errorMsg = "";
 };
 
-#endif  // POWERFLOW_H
+#endif // POWERFLOW_H