diff --git a/qucs/qucs-wideband-matching/GRABIM.cpp b/qucs/qucs-wideband-matching/GRABIM.cpp
index fc2fd7ee79..399484acf1 100644
--- a/qucs/qucs-wideband-matching/GRABIM.cpp
+++ b/qucs/qucs-wideband-matching/GRABIM.cpp
@@ -45,24 +45,13 @@ GRABIM::~GRABIM()
 {}
 
 
-// It sets the source impedance vs frequency
-int GRABIM::SetSourceImpedance(vector<complex<double>> zs)
+int GRABIM::SetData(MatchingData data)
 {
-    ZS = zs;
-    return 0;
-}
-
-// It loads a s2p file and calculates the input impedance
-int GRABIM::SetLoadImpedance(vector<complex<double>> zl)
-{
-    ZL = zl;
-    return 0;
-}
-
-// It sets the load impedance vs frequency
-int GRABIM::SetFrequency(vector<double> f)
-{
-    freq = f;
+    ZS = data.ZS;
+    ZL = data.ZL;
+    Zin_maxg = data.Zin_maxg;
+    Zout_maxg = data.Zout_maxg;
+    freq = data.freq;
     return 0;
 }
 
@@ -109,11 +98,15 @@ double GRABIM::GetThreshold()
     return MatchingThreshold;
 }
 
+
 // Starts the engine
-GRABIM_Result GRABIM::RunGRABIM()
+GRABIM_Result GRABIM::RunGRABIM(vector<complex<double>> Z1_, vector<complex<double>> Z2_)
 {
     GRABIM_Result Res;
-
+    //Set input and output impedances
+    Z1 = Z1_;
+    Z2 = Z2_;
+    CreateTopologiesList();
     vector<vector<double>> Vopt;//This data structure which contains the optimum vector for each circuit topology.
     vector<string> candidates;//List of topologies. The grid search will be done for each one of them
     GridSearch_DifferentTopologies(Vopt, candidates);//Runs the grid search for each topology
@@ -164,8 +157,8 @@ GRABIM_Result GRABIM::RunGRABIM()
     //At this point, it was found the best candidate so the code below arranges data so as to produce the output data
     topology = candidates[best_index];
     Res.topology = topology.c_str();//Save network topology to create a Qucs schematic
-    Res.ZS = ZS;
-    Res.ZL = ZL;
+    Res.ZS = Z1;
+    Res.ZL = Z2;
     Res.freq = freq;
 
     //Initialize S param vectors
@@ -192,7 +185,7 @@ GRABIM_Result GRABIM::RunGRABIM()
     for (unsigned int i = 0; i < freq.size(); i++)
     {
         // Grid search
-        S_gridsearch = S2PEngine.getSparams(Res.x_grid_search, ZS[i], ZL[i], freq[i], topology);
+        S_gridsearch = S2PEngine.getSparams(Res.x_grid_search, Z1[i], Z2[i], freq[i], topology);
         Res.S11_gridsearch[i] = S_gridsearch(0,0);
         Res.S21_gridsearch[i] = S_gridsearch(1,0);
         Res.S12_gridsearch[i] = S_gridsearch(0,1);
@@ -201,7 +194,7 @@ GRABIM_Result GRABIM::RunGRABIM()
         // Local optimizer
         if (refine)
         {
-        S_nlopt = S2PEngine.getSparams(Res.x_local_opt, ZS[i], ZL[i], freq[i], topology);
+        S_nlopt = S2PEngine.getSparams(Res.x_local_opt, Z1[i], Z2[i], freq[i], topology);
         Res.S11_nlopt[i] = S_nlopt(0,0);
         Res.S21_nlopt[i] = S_nlopt(1,0);
         Res.S12_nlopt[i] = S_nlopt(0,1);
@@ -338,7 +331,6 @@ int GRABIM::GridSearch_DifferentTopologies(vector<vector<double>> & Vopt, vector
     double gridtest, opttopo = 1e12;
     vector<double> Vaux;
     string tag, best;
-
     unsigned int n_topo = TopoList.size();
     for (unsigned int i = 0; i < n_topo; i++)
     {
@@ -347,8 +339,8 @@ int GRABIM::GridSearch_DifferentTopologies(vector<vector<double>> & Vopt, vector
         TopoList.pop_front();
         TagList.pop_front();
         AutoSetInitialPivot();;
-        Vaux = GridSearch(); 
-        gridtest = CandidateEval(Vaux);  
+        Vaux = GridSearch();
+        gridtest = CandidateEval(Vaux);
         cout << 100.*(i+1)/n_topo<< "% completed. " << tag << ": S11_min = " <<  gridtest << " dB" << endl;
 
         if (gridtest < -8.)
@@ -530,7 +522,7 @@ vector<double> GRABIM::InspectCandidate(vector<double> xk)
         default://It is a transmission line or a stub
             if (xk[index] < 0)//Something is wrong...
             {
-                xk[index] = .5*(mean(abs(ZS))+mean(abs(ZL)));
+                xk[index] = .5*(mean(abs(Z1))+mean(abs(Z2)));
             }
             if (xk[index+1] < 0)
             {
@@ -566,19 +558,18 @@ double GRABIM::CandidateEval(vector<double> x)
     double fobj = -1e3;
     Mat S, ABCD;
     SparEngine S2PEngine;
-    vector<complex<double>> s11 (freq.size());
     for (unsigned int i = 0; i < freq.size(); i++)
     {
         if (ObjFun == ObjectiveFunction::NINF_S11dB)
         {
-            S = S2PEngine.getSparams(x, ZS[i], ZL[i], freq[i], topology);
+            S = S2PEngine.getSparams(x, Z1[i], Z2[i], freq[i], topology);
             if (abs(S(0,0)) > fobj) fobj = abs(S(0,0));//N inf
   
         }
         if (ObjFun == ObjectiveFunction::NINF_POWERTRANS)
         {
             ABCD = S2PEngine.getABCDmatrix(x, freq.at(i), topology);
-            fobj = CalcInvPowerTransfer(ABCD, ZS.at(i), ZL.at(i));
+            fobj = CalcInvPowerTransfer(ABCD, Z1.at(i), Z2.at(i));
         }
     }
     if (ObjFun == ObjectiveFunction::NINF_S11dB)fobj = 20*log10(fobj);//|grad{log(x)}| > |grad{x}| when x < 1;
@@ -620,7 +611,7 @@ void GRABIM::AutoSetInitialPivot()
     double meanf = .5*(min(freq)+max(freq));
     double meanw = 2*pi*meanf;
     double lambda4 = c0/(4.*meanf);
-    double meanZ = 0.75*(mean(abs(ZL))+mean(abs(ZS)));
+    double meanZ = 0.75*(mean(abs(Z1))+mean(abs(Z2)));
     queue <double> XINI;
     for (unsigned int i = 0; i< topology.size();i++)
     {
@@ -628,16 +619,16 @@ void GRABIM::AutoSetInitialPivot()
         if (!topology.substr(i,1).compare("2")) XINI.push(meanZ/meanw);
         if (!topology.substr(i,1).compare("1")) XINI.push(1/(meanZ*meanw));
         if (!topology.substr(i,1).compare("3")) XINI.push(1/(meanZ*meanw));
-        if((!topology.substr(i,1).compare("5"))||(!topology.substr(i,1).compare("6")))XINI.push(mean(real(ZS)+real(ZL))),XINI.push(lambda4);
+        if((!topology.substr(i,1).compare("5"))||(!topology.substr(i,1).compare("6")))XINI.push(mean(real(Z1)+real(Z2))),XINI.push(lambda4);
 
         if (!topology.substr(i,1).compare("4"))//Transmission line
         {
             // Here it is defined a linear impedance profile. In case the user selected 444, this tries
             // works real-to-real impedance transformer. Luckily, GRABIM will find a better result
             double Zi;
-            double m = (mean(real(ZL))-mean(real(ZS)))/topology.length();
+            double m = (mean(real(Z2))-mean(real(Z1)))/topology.length();
 
-            Zi = m*i+mean(real(ZS));
+            Zi = m*i+mean(real(Z1));
             XINI.push(Zi);
             XINI.push(lambda4);
         }
@@ -684,7 +675,15 @@ void GRABIM::CheckDim(vector<vector<double>> & C, vector<vector<double>> & xkq,
 // 2: LC networks up to 4 elements
 // 3: LC + TL networks up to 6 elements
 // 4: LC + TL + stubs networks up to 6 elements
+
 void GRABIM::setSearchMode(int mode)
+{
+    search_mode = mode;
+}
+
+
+
+void GRABIM::CreateTopologiesList()
 {
     string str;
     TopoList.clear();
@@ -692,7 +691,7 @@ void GRABIM::setSearchMode(int mode)
     list<string>::iterator it_topo, it_tag;
     it_topo = TopoList.begin();
     it_tag = TagList.begin();
-    switch (mode)
+    switch (search_mode)
     {
     case -1://Unique network (Only quick command line checks)
     {
@@ -1055,8 +1054,8 @@ vector<double> GRABIM::NelderMead(vector<double> x)
  return getRow(X, 0);
 }
 
-void GRABIM::SetAmplifierS2P(vector<vector<complex<double>>> amps2p)
+void GRABIM::SetAmplifierS2P(AmpData amps2p)
 {
-
+  Amplifier_SPAR = amps2p;
 }
 
diff --git a/qucs/qucs-wideband-matching/GRABIM.h b/qucs/qucs-wideband-matching/GRABIM.h
index 395fe3a5b6..7cf345c635 100644
--- a/qucs/qucs-wideband-matching/GRABIM.h
+++ b/qucs/qucs-wideband-matching/GRABIM.h
@@ -33,6 +33,7 @@ using namespace std;
 
 enum ObjectiveFunction {NINF_S11dB, NINF_POWERTRANS};
 
+
 typedef struct GRABIM_Result {
     vector<double> x_grid_search;
     double grid_val;
@@ -45,15 +46,29 @@ typedef struct GRABIM_Result {
     std::string topology;
     std::string source_path;
     std::string load_path;
-    std::string QucsVersion;
+
 
 } GRABIM_Result;
 
+typedef struct TwoPortCircuit{
+    GRABIM_Result InputMatching;
+    GRABIM_Result OutputMatching;
+    std::string amplifier_path;
+    std::string QucsVersion;
+
+}TwoPortCircuit;
+
 typedef struct AmpData{
-       deque<complex<double>> S11, S21, S12, S22;
+       vector<complex<double>> S11, S21, S12, S22;
+       double Z0;//Reference impedance at which the SPAR data was extracted
 } AmpData;
 
 
+typedef struct MatchingData{
+       vector<complex<double>> ZS, ZL, Zin_maxg, Zout_maxg;
+       vector<double> freq;
+} MatchingData;
+
 
 //Reference:
 // [1] Broadband direct-coupled and RF matching networks. Thomas R. Cuthbert, 1999
@@ -66,15 +81,13 @@ class GRABIM
     GRABIM();
     ~GRABIM();
 
-    int SetSourceImpedance(vector<complex<double>>);
-    int SetLoadImpedance(vector<complex<double>>);
-    int SetFrequency(vector<double>);
+    int SetData(MatchingData);
 
     int SetInitialPivot(vector<double>);
     vector<double> GetInitialPivot();
     int SetMatchingBand(double, double);
     int SetTopology(std::string);
-    GRABIM_Result RunGRABIM();
+    GRABIM_Result RunGRABIM(vector<complex<double>>, vector<complex<double>>);
     int SetMaxIterGridSearch(int);
     int SetThreshold(double);
     double GetThreshold();
@@ -86,9 +99,10 @@ class GRABIM
 
     void setTopoScript(std::string);
     void setSearchMode(int);
+    void CreateTopologiesList();
     void SimplifyNetwork(bool);
     bool refine;//Runs a Nelder-Mead local optimizer to refine results
-    void SetAmplifierS2P(vector<vector<complex<double>>>);
+    void SetAmplifierS2P(AmpData);
 
 private:
     vector<double> GridSearch();
@@ -98,8 +112,12 @@ class GRABIM
     vector<double> InspectCandidate(vector<double>);
     vector<double> x_ini;
 
+    vector<complex<double>> ZS, ZL, Zin_maxg, Zout_maxg;
     vector<double> freq;
-    vector<complex<double>> ZS, ZL;
+    vector<complex<double>> Z1, Z2;//One-port: Z1 = ZS, Z2 = ZL
+                                   //Two-port: Z1 = ZS, Z2 = Zin_maxg
+                                   //          Z1 = Zout_maxg, Z2 = ZL
+    AmpData Amplifier_SPAR;
 
     std::string topology;
     unsigned int Grid_MaxIter;
diff --git a/qucs/qucs-wideband-matching/MathOperations.cpp b/qucs/qucs-wideband-matching/MathOperations.cpp
index a8e7392ca7..9ece108829 100644
--- a/qucs/qucs-wideband-matching/MathOperations.cpp
+++ b/qucs/qucs-wideband-matching/MathOperations.cpp
@@ -72,9 +72,8 @@ vector<complex<double>>  interp(vector<double> f1, vector<complex<double>> Z, ve
    vector<double> realpart = getRealPart(Z);
    vector<double> imagpart = getImagPart(Z);
    vector<complex<double>> Z_(f2.size());
-  
    //Check f2 \in f1
-   if ( (min(f2) < min(f1)) || (max(f2) > max(f1)) )
+   if ( (round(min(f2)) < round(min(f1))) || (round(max(f2)) > round(max(f1))) )
    {
     cout << "ERROR: Cannot interporlate a sequence defined outside the original one" << endl;
    }
@@ -87,10 +86,10 @@ vector<complex<double>>  interp(vector<double> f1, vector<complex<double>> Z, ve
 	{
           i1 = closestIndex(f2, f1[ind]);
 	  i2 = closestIndex(f2, f1[ind+1]);
-	  nslots = i2-i1;
-
-          step_r = (realpart[ind + 1] - realpart[ind])/nslots;
-          step_i = (imagpart[ind + 1] - imagpart[ind])/nslots;
+	  nslots = i2-i1;//Number of samples between two consecutive indexes of f1
+        
+          step_r = (realpart[ind + 1] - realpart[ind])/nslots;//Slope of the real part
+          step_i = (imagpart[ind + 1] - imagpart[ind])/nslots;//Slope of the imaginary part
 	
 	  for (unsigned int i = a; i <= nslots; i++, aux++)
 	  {
@@ -98,7 +97,7 @@ vector<complex<double>>  interp(vector<double> f1, vector<complex<double>> Z, ve
 	    imagpart_ = imagpart[ind] + i*step_i;
             Z_[aux] = complex<double>(realpart_, imagpart_);
           }
-          if (a==0)a=1;
+          if (a==0)a=1;//This trick is done in order to use the first sample of the first block
 
 	}
    }
@@ -135,12 +134,26 @@ vector<double> getImagPart(vector<complex<double> > Z)
   return imagpart;
 }
 
-// Conjugates a complex vector
+// Conjugates a complex number
 complex<double> conj(complex<double> Z)
 {
   return complex<double>(real(Z), -imag(Z));
 }
 
+// Conjugates a complex vector
+vector<complex<double>> conj(vector<complex<double>> Z)
+{
+  for (unsigned int i = 0; i < Z.size(); i++) Z[i] = complex<double>(real(Z[i]), -imag(Z[i]));
+  return Z;
+}
+
+//Calculates the square root of the elements of the vector
+vector<double> sqrt(vector<double> U)
+{
+  for (unsigned int i = 0; i < U.size(); i++) U[i] = sqrt(U[i]);
+  return U;
+}
+
 
 // Creates a N-dim ones vector
 vector<double> ones(unsigned int N)
@@ -150,14 +163,6 @@ vector<double> ones(unsigned int N)
   return V;
 }
 
-// double * vector operation
-vector<double> operator * (double d, vector<double> Z)
-{
-  vector<double> V(Z.size());
-  for (unsigned int i=0; i< Z.size(); i++) V[i] = d*Z[i];
-  return V;
-}
-
 // double + vector operation
 vector<double> operator + (double u, vector<double> V)
 {
@@ -182,14 +187,46 @@ vector<double> operator + (vector<double> U, vector<double> V)
   return V_;
 }
 
+
+vector<complex<double>> operator + (vector<double> U, vector<complex<double>> V)
+{
+  vector<complex<double>> Z(U.size());
+  for (unsigned int i=0; i< V.size(); i++) Z[i] = complex<double>(U[i], 0) + V[i];
+  return Z;
+}
+
+vector<complex<double>> operator - (vector<double> U, vector<complex<double>> V)
+{
+  vector<complex<double>> Z(U.size());
+  for (unsigned int i=0; i< V.size(); i++) Z[i] = complex<double>(U[i], 0) - V[i];
+  return Z;
+}
+
 // vector - vector operation
 vector<double> operator - (vector<double> U, vector<double> V)
 {
   vector<double> V_(U.size());
-  for (auto a = U.begin(), b=V.begin(), c = V_.begin(); a != U.begin()+(U.size()-1); ++a, ++b, ++c) *c = *a - *b;
+  for (unsigned int i=0; i< V.size(); i++) V_[i] = U[i] - V[i];
   return V_;
 }
 
+// vector<double> - vector<double> operation
+vector<complex<double>> operator - (vector<complex<double>> U, vector<complex<double>> V)
+{
+  vector<complex<double>> V_(U.size());
+  for (unsigned int i=0; i< V.size(); i++) V_[i] = U[i] - V[i];
+  return V_;
+}
+
+// double * vector operation
+vector<double> operator * (double d, vector<double> Z)
+{
+  vector<double> V(Z.size());
+  for (unsigned int i=0; i< Z.size(); i++) V[i] = d*Z[i];
+  return V;
+}
+
+
 // vector * vector dot operation
 vector<double> operator * (vector<double> U, vector<double> V)
 {
@@ -198,6 +235,42 @@ vector<double> operator * (vector<double> U, vector<double> V)
   return V_;
 }
 
+vector<complex<double>> operator * (double d, vector<complex<double>> Z)
+{
+    vector<complex<double>> V(Z.size());
+    for (unsigned int i=0; i< Z.size(); i++) V[i] = d*Z[i];
+    return V;
+}
+
+// vector<complex> * vector<complex> dot product operation
+vector<complex<double>> operator * (vector<complex<double>> U, vector<complex<double>> V)
+{
+  vector<complex<double>> V_(U.size());
+  for (unsigned int i=0; i< U.size(); i++) V_[i] = U[i]*V[i];
+  return V_;
+}
+
+complex<double> operator / (double d, complex<double> z)
+{
+  return d*conj(z)/(abs(z)*abs(z));
+}
+
+vector<complex<double>> operator / (vector<double> U, vector<complex<double>> V)
+{
+  vector<complex<double>> Z(U.size());
+  for (unsigned int i=0; i< Z.size(); i++) Z[i] = U[i]/V[i];
+  return Z;
+}
+
+
+vector<complex<double>> operator / (vector<complex<double>> U, vector<complex<double>> V)
+{
+    vector<complex<double>> Z(U.size());
+    for (unsigned int i=0; i< Z.size(); i++) Z[i] = U[i]/V[i];
+    return Z;
+}
+
+
 // Prints the content of a real vector in the standard output
 void print(vector<double> V)
 {
diff --git a/qucs/qucs-wideband-matching/MathOperations.h b/qucs/qucs-wideband-matching/MathOperations.h
index 74e349927c..62fd106056 100644
--- a/qucs/qucs-wideband-matching/MathOperations.h
+++ b/qucs/qucs-wideband-matching/MathOperations.h
@@ -47,13 +47,23 @@ using namespace std;
  vector<double> getRealPart(vector<complex<double>>);
  vector<double> getImagPart(vector<complex<double>>);
  complex<double> conj(complex<double>);
+ vector<complex<double>> conj(vector<complex<double>>);
+ vector<double> sqrt(vector<double>);
  vector<double> ones(unsigned int);
  vector<double> operator * (double, vector<double>);
  vector<double> operator * (vector<double>, vector<double>);
+ vector<complex<double>> operator * (double, vector<complex<double>>);
+ vector<complex<double>> operator * (vector<complex<double>>, vector<complex<double>>);
+ vector<complex<double>> operator / (vector<double>, vector<complex<double>>);
+ vector<complex<double>> operator / (vector<complex<double>>, vector<complex<double>>);
+ complex<double> operator / (double, complex<double>);
  vector<double> operator + (double, vector<double>);
  vector<double> operator + (vector<double>, vector<double>);
+ vector<complex<double>> operator + (vector<double>, vector<complex<double>>);
  vector<double> operator - (double, vector<double>);
  vector<double> operator - (vector<double>, vector<double>);
+ vector<complex<double>> operator - (vector<double>, vector<complex<double>>);
+ vector<complex<double>> operator - (vector<complex<double>>, vector<complex<double>>);
  vector<complex<double>> operator + (vector<complex<double>>, vector<complex<double>>);
  void print(vector<double>);
  void print(vector<complex<double>>);
diff --git a/qucs/qucs-wideband-matching/io.cpp b/qucs/qucs-wideband-matching/io.cpp
index b261276804..a8a0163cd4 100644
--- a/qucs/qucs-wideband-matching/io.cpp
+++ b/qucs/qucs-wideband-matching/io.cpp
@@ -312,14 +312,12 @@ int IO::loadS1Pdata(std::string filepath, terminal Port)
     {
         fS = freq;
         ZS = Z;
-        if (!ZL.empty())ResampleImpedances();
     }
 
     if (Port == LOAD)//Set load port properties
     {
         fL = freq;
         ZL = Z;
-        if (!ZS.empty())ResampleImpedances();
     }
     return 0;
 }
@@ -332,7 +330,6 @@ int IO::loadS2Pdata(std::string filepath)
     {
         return -1;
     }
-
     std::string line;
     double freq_scale = 1;
     double Zref = 50;
@@ -351,13 +348,10 @@ int IO::loadS2Pdata(std::string filepath)
 
     int Rindex = line.find_last_of("r");
     Rindex = line.find_first_not_of(" ", Rindex);
-    Zref = atof(line.substr(Rindex+1).c_str());
+    AmpS2P.Z0 = atof(line.substr(Rindex+1).c_str());
 
     bool is_indB = (line.find("db") != string::npos);
-    bool RI = (line.find("ma") == string::npos);
-    bool isS_par = (line.find(" s ") != string::npos);
-    bool isZ_par = (line.find(" z ") != string::npos);
-
+    bool RI = (line.find("ri") != string::npos);
 
     while( getline(s2pfile, line) )
     {//Looking for the start of the raw data
@@ -442,8 +436,7 @@ int IO::loadS2Pdata(std::string filepath)
         qsize++;
     }while (std::getline(s2pfile, line));
     s2pfile.close();
-    deque<double> freq(qsize);
-    AmpData AmpS2P;
+    vector<double> freq(qsize);
 
     double S11m, S11a, S21m, S21a, S12m, S12a, S22m, S22a;
     for (unsigned int i = 0; i < qsize; i++)
@@ -484,9 +477,7 @@ int IO::loadS2Pdata(std::string filepath)
     	    AmpS2P.S22.push_back(complex<double>(S22m, S22a)); 
         }
         else
-        {   
-
-           
+        {            
             S11a = (pi/180)*S11a;
             S21a = (pi/180)*S21a;
             S12a = (pi/180)*S12a;
@@ -496,16 +487,12 @@ int IO::loadS2Pdata(std::string filepath)
     	    AmpS2P.S21.push_back(complex<double>(S21m,0)*complex<double>(cos(S21a), sin(S21a))); 
     	    AmpS2P.S12.push_back(complex<double>(S12m,0)*complex<double>(cos(S12a), sin(S12a))); 
     	    AmpS2P.S22.push_back(complex<double>(S22m,0)*complex<double>(cos(S22a), sin(S22a)));
-/*
-           cout << freq.at(i) << " S11:" <<S11m << ", " << S11a << "==> " << AmpS2P.S11.at(i) << endl;
-           cout << freq.at(i) << " S21:" <<S21m << ", " << S21a << "==> " << AmpS2P.S21.at(i) << endl;
-           cout << freq.at(i) << " S12:" <<S12m << ", " << S12a << "==> " << AmpS2P.S12.at(i) << endl;
-           cout << freq.at(i) << " S22:" <<S22m << ", " << S22a << "==> " << AmpS2P.S22.at(i) << endl << endl;*/
         }
  
 
 
     }
+    fAMP = freq;
     return 0;
 }
 
@@ -514,123 +501,181 @@ int IO::loadS2Pdata(std::string filepath)
 // spline or cubic interpolation, but it seems that they are not implemented in Armadillo
 int IO::ResampleImpedances()
 {
-    bool ready = false;
+    if (ZS.size() == 0) return 0;//Not set
+    if (ZL.size() == 0) return 0;//Not set
+    if ((Two_Port_Matching)&&(fAMP.size() == 0)) return 0;//Not set
+
     if (fmatching_min == -1) return 0;
 
-    //Check whether the vectors are already resampled or not
-    if ((ZS.size() == ZL.size()) && (ZS.size() == freq.size()))
+    freq = linspace(fmatching_min, fmatching_max, Nsamples);//Frequency vector employed for matching
+
+    //Source impedance
+    if (ZS.size() == 1)
     {
-        return 0;
+       complex<double> zs_temp = ZS[0];
+       ZS.resize(freq.size());
+       ZS = ones(ZS);
+       ZS = Product(ZS,zs_temp);
     }
-
-    double fmin, fmax;
-    /*We have to discern the following cases
-    1: ZS constant, ZL constant versus frequency
-    2: ZS s2p/s1p, ZL constant vs freq
-    3: ZS constant vs freq, ZS s2p/s1p
-    4: ZS s2p/s1p vs ZL s2p/s1p
-    */
-
-    fmin = fmatching_min;
-    fmax = fmatching_max;
-
-    if ((ZS.size() == 1) && (ZL.size() == 1)) //ZS constant, ZL constant vs frequency
+    else
     {
-        freq = linspace(fmatching_min, fmatching_max, Nsamples);//Available freqs
-        complex<double> zs_temp = ZS[0];
-        complex<double> zl_temp = ZL[0];
-        //Create ZS and ZL vectors with the same length than freq
-        ZS.resize(freq.size());
-        ZS = ones(ZS);
-        ZS = Product(ZS,zs_temp);
-
-        ZL.resize(freq.size());
-        ZL = ones(ZL);
-        ZL = Product(ZL,zl_temp);
-
-        return 0;
+        //Extract impedance lying on the specified frequency band defined by the user in the UI
+        unsigned int i1s = closestIndex(fS, fmatching_min);
+        unsigned int i2s = closestIndex(fS, fmatching_max);
+        vector<complex<double>> zs = SubVector(ZS, i1s, i2s);//Useful impedance data
+        vector<double> fs = SubVector(fS, i1s, i2s);//Frequency points where the impedance zs is sampled
+        vector<complex<double>> ZS_ = interp(fs, zs, freq);//Interpolation using input data
+        ZS.resize(ZS_.size());//Resize ZS array
+        ZS = ZS_;//Copy data
     }
 
-    if ((ZS.size() == 1) && (ZL.size() != 1)) //ZS constant vs ZL s2p/s1p
+    //Load impedance
+    if (ZL.size() == 1)
     {
-        if (min(fL) > fmin ) fmin = min(fL);
-        if (max(fL) < fmax ) fmax = max(fL);
-        fS = fL;
-        complex<double> zs_temp = ZS[0];
-        ZS.resize(ZL.size());
-        ZS = ones(ZL);
-        ZS = Product(ZS,zs_temp);
-        ready = true;
+       complex<double> zl_temp = ZL[0];
+       ZL.resize(freq.size());
+       ZL = ones(ZL);
+       ZL = Product(ZL,zl_temp);
     }
-    if ((ZS.size() != 1) && (ZL.size() == 1) && (!ready)) //ZS s2p/s1p vs ZL constant
+    else
     {
-        if (min(fS) > fmin ) fmin = min(fS);
-        if (max(fS) < fmax ) fmax = max(fS);
-        fL = fS;
-        complex<double> zl_temp = ZL[0];
-        ZL.resize(ZS.size());
-        ZL = ones(ZL);
-        ZL = Product(ZL,zl_temp);
-        ready = true;
+        //Extract impedance lying on the specified frequency band defined by the user in the UI
+        unsigned int i1l = closestIndex(fL, fmatching_min);
+        unsigned int i2l = closestIndex(fL, fmatching_max);
+        vector<complex<double>> zl = SubVector(ZL, i1l, i2l);//Useful impedance data
+        vector<double> fl = SubVector(fL, i1l, i2l);//Frequency points where the impedance zs is sampled
+        vector<complex<double>> ZL_ = interp(fl, zl, freq);//Interpolation using input data
+        ZL.resize(ZL_.size());//Resize ZS array
+        ZL = ZL_;//Copy data
     }
-    if ((ZS.size() != 1) && (ZL.size() != 1) && (!ready)) //ZS s2p/s1p vs ZL s2p/s1p
-    {
-          double fmin_aux, fmax_aux;
-          //Define vector of available freqs
-          (min(fS) > min(fL)) ? fmin_aux = min(fS) : fmin_aux = min(fL);
-          (max(fS) > max(fL)) ? fmax_aux = max(fL) : fmax_aux = max(fS);
 
-          if (fmin_aux < fmin ) fmin = fmin_aux;
-          if (fmax_aux > fmax ) fmax = fmax_aux;
+    if (AmpS2P.S11.empty()) return 0;//One port matching
 
-          //int N;//Number of points for linear interpolation
-          //(ZS.n_elem > ZL.n_elem) ? N = 2*ZS.n_elem : N = 2*ZL.n_elem;
+    //Amplifier SPAR data
+    if (AmpS2P.S11.size() == 1)//Constant SPAR vs frequency. Corner case, but some user may input data like this
+    {
+       complex<double> S11_temp = AmpS2P.S11[0];
+       complex<double> S21_temp = AmpS2P.S21[0];
+       complex<double> S12_temp = AmpS2P.S12[0];
+       complex<double> S22_temp = AmpS2P.S22[0];
+
+       AmpS2P.S11.resize(freq.size());
+       AmpS2P.S21.resize(freq.size());
+       AmpS2P.S12.resize(freq.size());
+       AmpS2P.S22.resize(freq.size());
+
+       AmpS2P.S11 = ones(AmpS2P.S11);
+       AmpS2P.S21 = ones(AmpS2P.S21);
+       AmpS2P.S12 = ones(AmpS2P.S12);
+       AmpS2P.S22 = ones(AmpS2P.S22);
+
+       AmpS2P.S11 = Product(AmpS2P.S11,S11_temp);
+       AmpS2P.S21 = Product(AmpS2P.S21,S21_temp);
+       AmpS2P.S12 = Product(AmpS2P.S12,S12_temp);
+       AmpS2P.S22 = Product(AmpS2P.S22,S22_temp);
     }
+   else
+   {
+        //Extract impedance lying on the specified frequency band defined by the user in the UI
+        unsigned int i1amp = closestIndex(fAMP, fmatching_min);
+        unsigned int i2amp = closestIndex(fAMP, fmatching_max);
+        
+        vector<complex<double>> s11 = SubVector(AmpS2P.S11, i1amp, i2amp);
+        vector<complex<double>> s21 = SubVector(AmpS2P.S21, i1amp, i2amp);
+        vector<complex<double>> s12 = SubVector(AmpS2P.S12, i1amp, i2amp);
+        vector<complex<double>> s22 = SubVector(AmpS2P.S22, i1amp, i2amp);
+        vector<double> fa = SubVector(fAMP, i1amp, i2amp);//Frequency points where the impedance zs is sampled
+
+        vector<complex<double>> S11_ = interp(fa, s11, freq);//Interpolation using input data
+        vector<complex<double>> S21_ = interp(fa, s21, freq);
+        vector<complex<double>> S12_ = interp(fa, s12, freq);
+        vector<complex<double>> S22_ = interp(fa, s22, freq);
+
+        AmpS2P.S11.resize(S11_.size());
+        AmpS2P.S21.resize(S21_.size());
+        AmpS2P.S12.resize(S12_.size());
+        AmpS2P.S22.resize(S22_.size());
+
+        AmpS2P.S11 = S11_;//Copy data
+        AmpS2P.S21 = S21_;
+        AmpS2P.S12 = S12_;
+        AmpS2P.S22 = S22_;
+   }
 
-        freq = linspace(fmin, fmax, Nsamples);//Available freqs
-        unsigned int i1s = closestIndex(fS, fmin)-1;
-        unsigned int i2s = closestIndex(fS, fmax)+1;
-        unsigned int i1l = closestIndex(fL, fmin)-1;
-        unsigned int i2l = closestIndex(fL, fmax)+1;  
-
-        vector<complex<double>> zs = SubVector(ZS, i1s, i2s);
-        vector<complex<double>> zl = SubVector(ZL, i1l, i2l);
-        vector<double> fs = SubVector(fS, i1s, i2s);
-        vector<double> fl = SubVector(fL, i1l, i2l);
-
-        vector<complex<double>> ZS_ = interp(fs, zs, freq);
-        vector<complex<double>> ZL_ = interp(fl, zl, freq);  
-          
-        ZS.resize(ZS_.size());
-        ZL.resize(ZL_.size());
-
-        ZS = ZS_;
-        ZL = ZL_;
-        return 0;
+   //Finally, we calculate conj(Zin) and conj(Zout) for achieving maximum gain on a two-port device
+
+    /***********************************************************************
+      (gamma_S)                (gamma_in)    (gamma_out)                (gamma_L)
+              _________________       _______        __________________
+             |                 | Zin |       | Zout |                  |
+        ZS---| INPUT MATCHING  |-----|  TRT  |------| OUTPUT MATCHING  |---ZL
+             |_________________|     |_______|      |__________________|
+
+    ************************************************************************/
+
+   vector<complex<double>> delta = AmpS2P.S11*AmpS2P.S22 - AmpS2P.S21*AmpS2P.S12;
+   vector<double> B1 = 1 + abs(AmpS2P.S11)*abs(AmpS2P.S11) - abs(AmpS2P.S22)*abs(AmpS2P.S22) - abs(delta)*abs(delta);
+   vector<double> B2 = 1 + abs(AmpS2P.S22)*abs(AmpS2P.S22) - abs(AmpS2P.S11)*abs(AmpS2P.S11) - abs(delta)*abs(delta);
+   vector<complex<double>> C1 = AmpS2P.S11 - delta*conj(AmpS2P.S22);
+   vector<complex<double>> C2 = AmpS2P.S22 - delta*conj(AmpS2P.S11);
+   complex<double> gamma_S, gamma_L;
+   int ext_code = 0;
+   Zin_maxg.resize(C1.size());
+   Zout_maxg.resize(C1.size());
+   //The sign of the square root must be chosen for each frequency so here we cannot apply vector operations
+   for (int i = 0; i < C1.size(); i++)
+   {
+      gamma_S = (B1[i]-sqrt(B1[i]*B1[i] - 4.*abs(C1[i])*abs(C1[i])))/(2.*C1[i]);
+      gamma_L = (B2[i]-sqrt(B2[i]*B2[i] - 4.*abs(C2[i])*abs(C2[i])))/(2.*C2[i]);
+
+      if ((gamma_S != gamma_S) || (gamma_L != gamma_L))//Check if Nan
+      {
+          ext_code = -2;
+          Zin_maxg[i]  = AmpS2P.Z0*(complex<double>(1,0)+conj(AmpS2P.S11.at(i)))/(complex<double>(1,0)-conj(AmpS2P.S11.at(i)));
+          Zout_maxg[i] = AmpS2P.Z0*(complex<double>(1,0)+conj(AmpS2P.S22.at(i)))/(complex<double>(1,0)-conj(AmpS2P.S22.at(i)));
+      }
+      else
+      {
+         Zin_maxg[i]  = AmpS2P.Z0*(complex<double>(1,0)+conj(gamma_S))/(complex<double>(1,0)-conj(gamma_S));
+         Zout_maxg[i] = AmpS2P.Z0*(complex<double>(1,0)+conj(gamma_L))/(complex<double>(1,0)-conj(gamma_L));
+         continue;
+      }
+      if (real(Zin_maxg[i]) <= 0)
+      {
+        gamma_S = (B1[i]+sqrt(B1[i]*B1[i] - 4.*abs(C1[i])*abs(C1[i])))/(2.*C1[i]);
+        Zin_maxg[i]  = ZS[i]*(complex<double>(1,0)+conj(gamma_S))/(complex<double>(1,0)-conj(gamma_S));
+      }
 
-}
+      if (real(Zout_maxg[i]) <= 0)
+      {
+        gamma_L = (B2[i]+sqrt(B2[i]*B2[i] - 4.*abs(C2[i])*abs(C2[i])))/(2.*C2[i]);
+        Zout_maxg[i] = ZL[i]*(complex<double>(1,0)+conj(gamma_L))/(complex<double>(1,0)-conj(gamma_L));
+      }
+   }
 
-vector<complex<double>> IO::getSourceImpedance()
-{
-    ResampleImpedances();
-    return ZS;
+   //Notice that the equation below difer from the typical equations for two-port matching in the sense that
+   //here ZS and ZL are used instead of Z0. Although in general the amplifier is matched to Z0=50 or 75\Ohm
+   //this tool aims to treat the impedance matching problem in a more general fashion
+
+
+   return ext_code;
 }
-vector<complex<double>> IO::getLoadImpedance()
+
+MatchingData IO::getMatchingData()
 {
-    ResampleImpedances();
-    return ZL;
+   MatchingData data;
+   data.freq = freq;
+   data.ZS = ZS;
+   data.ZL = ZL;
+   data.Zin_maxg = Zin_maxg;
+   data.Zout_maxg = Zout_maxg;
+   return data;
 }
 
 
-vector<vector<complex<double>>> IO::getAmplifierS2P()
-{
-    ResampleImpedances();
-    return AmpS2P;
-}
 
 vector<double> IO::getFrequency()
 {
-    ResampleImpedances();
     return freq;
 }
 
@@ -656,26 +701,6 @@ void IO::set_matching_band(double fmin, double fmax)
     freq = linspace(fmatching_min, fmatching_max, Nsamples);//Available freqs
 }
 
-// Set load/source matching impedace
-int IO::setMatchingImpedances()
-{
-    //Select impedance values at the matching band
-    int index1 = getFreqIndex(fmatching_min);
-    int index2 = getFreqIndex(fmatching_max);
-
-    vector<complex<double>> ZS_(&ZS[index1], &ZS[index2+1]);
-    vector<complex<double>> ZL_(&ZL[index1], &ZL[index2+1]);
-
-    ZS_matching = ZS_;
-    ZL_matching = ZL_;
-
-    //Now, it's needed to update the matching band according to f_analysis values
-    vector<double>freq_(&freq[index1], &freq[index2]);
-    f_matching = freq_;
-
-    return 0;
-}
-
 
 // Gets the index of a given frequency
 int IO::getFreqIndex(double f1)
@@ -683,18 +708,6 @@ int IO::getFreqIndex(double f1)
     return closestIndex(freq, f1) ;
 }
 
-/*
-//Sets local optimisation algorithm
-void IO::setLocalOptimiser(nlopt::algorithm NL)
-{
-    LocalOptAlgo = NL;
-}
-
-
-nlopt::algorithm IO::getLocalOptimiser()
-{
-    return LocalOptAlgo;
-}*/
 
 
 //Get freq scale from a string line
@@ -726,15 +739,15 @@ double IO::getS2PfreqScale(string line)
 
 
 // This function exports the best topology found to a Qucs schematic
-int IO::ExportQucsSchematic(GRABIM_Result R, string QucsFile)
+int IO::ExportQucsSchematic(TwoPortCircuit TPC, string QucsFile)
 {
     std::string wirestr = "";
     std::string componentstr = "";
     std::string paintingstr = "";
     int x_pos = 0;
-    SchematicParser(R, x_pos, componentstr, wirestr, paintingstr);
-    if (R.QucsVersion.empty()) R.QucsVersion = "0.0.19";
-    CreateSchematic(componentstr, wirestr, paintingstr, R.QucsVersion, QucsFile);
+    SchematicParser(TPC, x_pos, componentstr, wirestr, paintingstr);
+    if (TPC.QucsVersion.empty()) TPC.QucsVersion = "0.0.19";
+    CreateSchematic(componentstr, wirestr, paintingstr, TPC.QucsVersion, QucsFile);
     return 0;
 }
 
@@ -748,10 +761,8 @@ void IO::generateConstant_s1p(string datapath, complex<double> Z)
 
 // Given a string code of inductors, capacitors and transmission lines, it generates the Qucs network. Notice that the schematic is split into
 // three part: components, wires and paintings, all of them are passed by reference.
-int IO::SchematicParser(GRABIM_Result R, int & x_pos, string & componentstr, string & wirestr, string & paintingstr)
+int IO::SchematicParser(TwoPortCircuit TPC, int & x_pos, string & componentstr, string & wirestr, string & paintingstr)
 {
-    string component;
-    int x_series = 150, x_shunt = 120;//x-axis spacing depending on whether the component is placed in a series or shunt configuration
     //Clear input strings (just in case)
     componentstr = "";
     wirestr = "";
@@ -764,6 +775,8 @@ int IO::SchematicParser(GRABIM_Result R, int & x_pos, string & componentstr, str
     // 2: Port 2
     // 3: Port 1, Port 2 and S parameter simulation
 
+    GRABIM_Result R = TPC.InputMatching;
+
     if ((R.source_path.empty()) && (abs(ZS.at(0).imag()) < 1e-3) && (ZS.at(0).real() > 1e-3))
     {//Conventional term
         componentstr += "<Pac P1 1 " + Num2String(x_pos) + " -30 18 -26 0 1 \"1\" 1 \"" + Num2String(ZS.at(0).real()) + " Ohm\" 1 \"0 dBm\" 0 \"1 GHz\" 0>\n";
@@ -789,6 +802,102 @@ int IO::SchematicParser(GRABIM_Result R, int & x_pos, string & componentstr, str
         x_pos +=100;
     }
 
+    //Draw ladder network
+    CreateLadder(R, componentstr, wirestr, x_pos);
+    if (!TPC.OutputMatching.freq.empty())
+    {//Two-port network
+        //Place S2P device
+        x_pos += 30;
+        componentstr += "<SPfile X1 1 " + Num2String(x_pos) + " -120 -26 -67 0 0 \"" + TPC.amplifier_path + "\" 1 \"rectangular\" 0 \"linear\" 0 \"open\" 0 \"2\" 0>\n";
+        componentstr += "<GND * 1 " + Num2String(x_pos) +" -90 0 0 0 0>\n";
+        x_pos += 30;
+        wirestr += "<" + Num2String(x_pos+30) + " -120 " + Num2String(x_pos) + " -120>\n";
+        R = TPC.OutputMatching;
+        CreateLadder(R, componentstr, wirestr, x_pos);
+    }
+
+    int spacing = 30;
+    x_pos += spacing;
+
+    if ((R.load_path.empty()) && (abs(ZL.at(0).imag()) < 1e-3) && (ZL.at(0).real() > 1e-3))
+    {//Conventional term
+        componentstr += "<Pac P1 1 " + Num2String(x_pos) + " -30 18 -26 0 1 \"1\" 1 \"" + Num2String(ZL.at(0).real()) + " Ohm\" 1 \"0 dBm\" 0 \"1 GHz\" 0>\n";
+        componentstr += "<GND * 1 "+Num2String(x_pos)+" 0 0 0 0 0>\n";
+
+        wirestr += "<" + Num2String(x_pos) + " -60 " + Num2String(x_pos) + " -120 "" 0 0 0 "">\n";
+        wirestr += "<" + Num2String(x_pos-spacing) + " -120 " + Num2String(x_pos) + " -120 "" 0 0 0 "">\n";
+        term2 = true;//Add S22_dB = dB(S[2,2]) equation
+    }
+    else
+    {//Place a S-param file
+        if (R.load_path.empty())
+        {
+          string s1p_path = tmp_path + "/ZL_data.s1p";
+          generateConstant_s1p(s1p_path, ZL.at(0));
+          R.load_path =  s1p_path;
+        }
+        componentstr += "<SPfile X1 1 " + Num2String(x_pos) + " -120 -26 -67 0 0 \"" + R.load_path + "\" 1 \"rectangular\" 0 \"linear\" 0 \"open\" 0 \"1\" 0>\n";
+        componentstr += "<GND * 1 " + Num2String(x_pos) + " -90 0 0 0 0>\n";
+    }
+
+    //Finally, add the S-par simulation block
+    double fstart = fmatching_min*0.3;
+    double fstop = fmatching_max*1.7;
+    componentstr += "<.SP SP1 1 200 200 0 67 0 0 \"lin\" 1 \"" + Num2String(fstart) + "\" 1 \"" + Num2String(fstop) + "\" 1 \"300\" 1 \"no\" 0 \"1\" 0 \"2\" 0>\n";
+   
+    //Add equations
+    string eqn = ""; 
+    if (term1) eqn = "\"S11_dB=dB(S[1,1])\" 1 ";
+    if (term2) eqn += "\"S22_dB=dB(S[2,2])\" 1 ";
+    if (term1 && term2) eqn += "\"S21_dB=dB(S[2,1])\" 1 ";
+    if (term1 || term2) componentstr += "<Eqn Eqn1 1 50 200 -28 15 0 0 " + eqn + "\"yes\" 0>\n";
+
+    return 0;
+}
+
+
+//-----------------------------------------------------------------------------
+// Given the components, wires and paintings, it creates the schematic and copies on the clipboard
+bool IO::CreateSchematic(string components, string wires, string paintings, string QucsVersion, string QucsFilePath)
+{
+    //Header
+    std::string Schematic = "<Qucs Schematic " + QucsVersion + ">\n";
+
+    //Add components
+    Schematic += "<Components>\n";
+    Schematic += components;
+    Schematic += "</Components>\n";
+
+    //Add wires
+    Schematic+= "<Wires>\n";
+    Schematic += wires;
+    Schematic+= "</Wires>\n";
+
+    //Add paintings
+    Schematic += "<Paintings>\n";
+    Schematic += paintings;
+    Schematic += "</Paintings>\n";
+
+    if (QucsFilePath.empty())
+    {
+       QApplication::clipboard()->setText(QString(Schematic.c_str()), QClipboard::Clipboard);//Copy into clipboard
+    }
+    else//Dump content into a file
+    {
+       std::ofstream QucsFile(QucsFilePath.c_str(), ios_base::out);
+       QucsFile << Schematic;
+       QucsFile.close();
+    }
+
+    return true;
+}
+
+
+
+void IO::CreateLadder(GRABIM_Result R, string & componentstr, string & wirestr, int & x_pos)
+{
+    string component;
+    int x_series = 150, x_shunt = 120;//x-axis spacing depending on whether the component is placed in a series or shunt configuration
 
     // The string format is as follows: "XX<value>;XX<value2>;...XX<valueN>;"
     // where XX, YY, ZZ define the type of component and its configuration.
@@ -866,82 +975,9 @@ int IO::SchematicParser(GRABIM_Result R, int & x_pos, string & componentstr, str
         }
 
     }
-    double spacing = 30;
-    x_pos += spacing;
-
-    if ((R.load_path.empty()) && (abs(ZL.at(0).imag()) < 1e-3) && (ZL.at(0).real() > 1e-3))
-    {//Conventional term
-        componentstr += "<Pac P1 1 " + Num2String(x_pos) + " -30 18 -26 0 1 \"1\" 1 \"" + Num2String(ZL.at(0).real()) + " Ohm\" 1 \"0 dBm\" 0 \"1 GHz\" 0>\n";
-        componentstr += "<GND * 1 "+Num2String(x_pos)+" 0 0 0 0 0>\n";
-
-        wirestr += "<" + Num2String(x_pos) + " -60 " + Num2String(x_pos) + " -120 "" 0 0 0 "">\n";
-        wirestr += "<" + Num2String(x_pos-spacing) + " -120 " + Num2String(x_pos) + " -120 "" 0 0 0 "">\n";
-        term2 = true;//Add S22_dB = dB(S[2,2]) equation
-    }
-    else
-    {//Place a S-param file
-        if (R.load_path.empty())
-        {
-          string s1p_path = tmp_path + "/ZL_data.s1p";
-          generateConstant_s1p(s1p_path, ZL.at(0));
-          R.load_path =  s1p_path;
-        }
-        componentstr += "<SPfile X1 1 " + Num2String(x_pos) + " -120 -26 -67 0 0 \"" + R.load_path + "\" 1 \"rectangular\" 0 \"linear\" 0 \"open\" 0 \"1\" 0>\n";
-        componentstr += "<GND * 1 " + Num2String(x_pos) + " -90 0 0 0 0>\n";
-    }
-
-    //Finally, add the S-par simulation block
-    double fstart = fmatching_min*0.3;
-    double fstop = fmatching_max*1.7;
-    componentstr += "<.SP SP1 1 200 200 0 67 0 0 \"lin\" 1 \"" + Num2String(fstart) + "\" 1 \"" + Num2String(fstop) + "\" 1 \"300\" 1 \"no\" 0 \"1\" 0 \"2\" 0>\n";
-   
-    //Add equations
-    string eqn = ""; 
-    if (term1) eqn = "\"S11_dB=dB(S[1,1])\" 1 ";
-    if (term2) eqn += "\"S22_dB=dB(S[2,2])\" 1 ";
-    if (term1 && term2) eqn += "\"S21_dB=dB(S[2,1])\" 1 ";
-    if (term1 || term2) componentstr += "<Eqn Eqn1 1 50 200 -28 15 0 0 " + eqn + "\"yes\" 0>\n";
-
-    return 0;
 }
 
 
-//-----------------------------------------------------------------------------
-// Given the components, wires and paintings, it creates the schematic and copies on the clipboard
-bool IO::CreateSchematic(string components, string wires, string paintings, string QucsVersion, string QucsFilePath)
-{
-    //Header
-    std::string Schematic = "<Qucs Schematic " + QucsVersion + ">\n";
-
-    //Add components
-    Schematic += "<Components>\n";
-    Schematic += components;
-    Schematic += "</Components>\n";
-
-    //Add wires
-    Schematic+= "<Wires>\n";
-    Schematic += wires;
-    Schematic+= "</Wires>\n";
-
-    //Add paintings
-    Schematic += "<Paintings>\n";
-    Schematic += paintings;
-    Schematic += "</Paintings>\n";
-
-    if (QucsFilePath.empty())
-    {
-       QApplication::clipboard()->setText(QString(Schematic.c_str()), QClipboard::Clipboard);//Copy into clipboard
-    }
-    else//Dump content into a file
-    {
-       std::ofstream QucsFile(QucsFilePath.c_str(), ios_base::out);
-       QucsFile << Schematic;
-       QucsFile.close();
-    }
-
-    return true;
-}
-
 string IO::Num2String(double Num)
 {
   char c = 0;
diff --git a/qucs/qucs-wideband-matching/io.h b/qucs/qucs-wideband-matching/io.h
index c397a824f2..c41ab78728 100644
--- a/qucs/qucs-wideband-matching/io.h
+++ b/qucs/qucs-wideband-matching/io.h
@@ -49,42 +49,38 @@ class IO
     int loadS1Pdata(std::string, terminal);
     int loadS2Pdata(std::string);
     int ResampleImpedances();
-    vector<complex<double>> getSourceImpedance();
-    vector<complex<double>> getLoadImpedance();
-    vector<vector<complex<double>>> getAmplifierS2P();
     vector<double> getFrequency();
     void set_constant_ZS_vs_freq(complex<double>);
     void set_constant_ZL_vs_freq(complex<double>);
     void set_matching_band(double, double);
-    int ExportQucsSchematic(GRABIM_Result, string);
+    int ExportQucsSchematic(TwoPortCircuit, string);
     void PrintNetwork_StandardOutput(GRABIM_Result);
     void UseClipboard(bool);
     bool UseGNUplot;
     string tmp_path;//Path to a temporary directory for data dumping
+    bool Two_Port_Matching;
 
-
-private:
     // ZS and ZL are the source and load impedances, respectively whereas fS and fL indicates the frequencies where
     // ZS and ZL were sampled
     vector<complex<double>> ZS, ZL;
-    vector<vector<complex<double>>> AmpS2P;
-    vector<double>  fS, fL;
-
+    vector<complex<double>> Zin_maxg, Zout_maxg;//Optimal Zin and Zout for achieving maximum gain on a two-port device
     vector<double> freq;//More often than not, ZS and ZL are sampled at different frecuencies, so it is necessary to have
     // common frequency vector for pairing ZS and ZL.
+    MatchingData getMatchingData();
+    AmpData AmpS2P;
 
+private:
+    vector<double>  fS, fL, fAMP;//Frequencies at which the input source, load and amplifier data are sampled
 
     double fmatching_min, fmatching_max;
     int getFreqIndex(double);
 
-    vector<complex<double>> ZS_matching, ZL_matching;
     vector<double> f_matching;
-    int setMatchingImpedances();
     double getS2PfreqScale(string line);
-   // nlopt::algorithm LocalOptAlgo;
 
-    int SchematicParser(GRABIM_Result, int &, string &, string &, string &);
+    int SchematicParser(TwoPortCircuit, int &, string &, string &, string &);
     bool CreateSchematic(string, string, string, string, string);
+    void CreateLadder(GRABIM_Result, string &, string &, int &);
 
     int Nsamples;//Impedance samples within matching band
 
diff --git a/qucs/qucs-wideband-matching/main.cpp b/qucs/qucs-wideband-matching/main.cpp
index a319f3503d..8b82b75497 100644
--- a/qucs/qucs-wideband-matching/main.cpp
+++ b/qucs/qucs-wideband-matching/main.cpp
@@ -262,16 +262,15 @@ int main(int argc, char *argv[])
 
         GRABIM * MatchingObject = new GRABIM();
         // Impedance and frequency settings
-        MatchingObject->SetSourceImpedance(inout_operations->getSourceImpedance());
-        MatchingObject->SetLoadImpedance(inout_operations->getLoadImpedance());
-        MatchingObject->SetFrequency(inout_operations->getFrequency());
+        MatchingData data = inout_operations->getMatchingData();
+        MatchingObject->SetData(data);
         MatchingObject->setTopoScript(TopoScript_path);
         MatchingObject->SetTopology(topo);
         MatchingObject->setSearchMode(search_mode);
         MatchingObject->SimplifyNetwork(!no_simplify);
 
 
-        GRABIM_Result R = MatchingObject->RunGRABIM();//Runs GRABIM. Please bear in mind that this is not a rigorous implementation of [1] 
+        GRABIM_Result R = MatchingObject->RunGRABIM(inout_operations->ZS, inout_operations->ZL);//Runs GRABIM. Please bear in mind that this is not a rigorous implementation of [1]
         // Biggest differences between this code and [1]:
         // 1) The candidate vector is always defined in natural units rather than logarithmic units.
         // 2) Frequency is not normalized.
@@ -296,8 +295,10 @@ int main(int argc, char *argv[])
         cout << "The GNUplot data was written at " << GNUplot_path << endl;
         cout << "You can find a sample script to plot S11 in the same folder as the source code" << endl;
 
-        inout_operations->exportGNUplot(R, GNUplot_path, false);
-        inout_operations->ExportQucsSchematic(R, QucsSchPath);
+        TwoPortCircuit TPC;
+        TPC.InputMatching = R;
+        inout_operations->exportGNUplot(TPC.InputMatching, GNUplot_path, false);
+        inout_operations->ExportQucsSchematic(TPC, QucsSchPath);
         delete inout_operations;
         delete MatchingObject;
     }
diff --git a/qucs/qucs-wideband-matching/ui.cpp b/qucs/qucs-wideband-matching/ui.cpp
index 5c3bc96d9c..2655af7283 100644
--- a/qucs/qucs-wideband-matching/ui.cpp
+++ b/qucs/qucs-wideband-matching/ui.cpp
@@ -51,6 +51,7 @@ ui::ui()
       ZSOhmLabel = new QLabel("Ohm");
       ZLOhmLabel = new QLabel("Ohm");
 
+
       ZSOhmLabel->setVisible(false);
       ZLOhmLabel->setVisible(false);
 
@@ -66,14 +67,14 @@ ui::ui()
       ConstantZLLayout->addWidget(FixedZLLineedit);
       ConstantZLLayout->addWidget(ZLOhmLabel);
 
-      QLabel *SourceLabel = new QLabel("Source");
+      SourceLabel = new QLabel("Source");
       SourceLabel->setStyleSheet("QLabel {font-weight: bold; border: 1px solid black;}");
       SourceLabel->setAlignment(Qt::AlignHCenter);
       SourceLayout->addWidget(SourceLabel);
       SourceLayout->addWidget(SourceFileButton);
       SourceLayout->addLayout(ConstantZSLayout);
 
-      QLabel *LoadLabel = new QLabel("Load");
+      LoadLabel = new QLabel("Load");
       LoadLabel->setStyleSheet("QLabel {font-weight: bold; border: 1px solid black;}");
       LoadLabel->setAlignment(Qt::AlignHCenter);
       LoadLayout->addWidget(LoadLabel);
@@ -90,9 +91,11 @@ ui::ui()
       LoadLayout->addWidget(FixedZLCheckbox);
 
       SourceImpGroupBox = new QGroupBox();
+      SourceImpGroupBox->setMinimumSize(130,20);
       SourceImpGroupBox->setStyleSheet("QGroupBox{  border: 1px solid black;}");
       LoadImpGroupBox = new QGroupBox();
       LoadImpGroupBox->setStyleSheet("QGroupBox{  border: 1px solid black;}");
+      LoadImpGroupBox->setMinimumSize(130,20);
 
       SourceImpGroupBox->setLayout(SourceLayout);
       LoadImpGroupBox->setLayout(LoadLayout);
@@ -178,6 +181,7 @@ ui::ui()
        OptionsLayout->addWidget(GNUplotButton, 2, 1);
        GNUplot_path = "/tmp";
 
+
       //Create go/cancel buttons
       RunButton = new QPushButton("Calculate and put into clipboard");
       ButtonsLayout->addWidget(RunButton);
@@ -267,34 +271,52 @@ void ui::go_clicked()
     if (LoadFile.contains(".s1p")) formatLoad = 0;
 
 
-    if ((formatSource != 0) && (!FixedZSCheckbox->isChecked()))
+    if ((formatSource != 0) && (!FixedZSCheckbox->isChecked()) && (!TwoPortMatchingCheckbox->isChecked()))
     {
         QMessageBox::warning(0, QObject::tr("Error"),
                              QObject::tr("Unknown source impedace"));
         return;
     }
 
-    if ((formatLoad != 0)  && (!FixedZLCheckbox->isChecked()))
+    if ((formatLoad != 0)  && (!FixedZLCheckbox->isChecked()) && (!TwoPortMatchingCheckbox->isChecked()))
     {
         QMessageBox::warning(0, QObject::tr("Error"),
                              QObject::tr("Unknown load impedance"));
         return;
     }
 
-    //Impedance data paths were already specified, let's proceed to bring the S-par data into memory
     IO * inout_operations = new IO();
 
-    if (!FixedZSCheckbox->isChecked())//Read source impedance from file
+    if (TwoPortMatchingCheckbox->isChecked())
+    {
+       inout_operations->Two_Port_Matching=true;
+       inout_operations->loadS2Pdata(AmpFile.toStdString());
+    }
+    else
+    {
+       inout_operations->Two_Port_Matching=false;
+    }
+
+
+    //Impedance data paths were already specified, let's proceed to bring the S-par data into memory
+
+    if (!FixedZSCheckbox->isChecked() && (!TwoPortMatchingCheckbox->isChecked()))//Read source impedance from file
     {
         inout_operations->loadS1Pdata(SourceFile.toStdString(), SOURCE);//s1p
     }
     else//Set constant source impedance
     {
        complex<double> zs_temp;
+       if (!TwoPortMatchingCheckbox->isChecked())
+       {
        QByteArray ba = FixedZSLineedit->text().toLatin1();
        char * text = ba.data();
        zs_temp = getComplexImpedanceFromText(text);
-
+       }
+       else
+       {
+         zs_temp = inout_operations->AmpS2P.Z0;
+       }
        if (real(zs_temp) == -1)//Check if the input value is correct
        {
            QMessageBox::warning(0, QObject::tr("Error"),
@@ -303,25 +325,24 @@ void ui::go_clicked()
        }
        inout_operations->set_constant_ZS_vs_freq(zs_temp);
     }
-    if (TwoPortMatchingCheckbox->isChecked())
-    {
-       inout_operations->loadS2Pdata(AmpFile.toStdString());
-    }
-    else
-    {
-       inout_operations->loadS2Pdata("");
-    } 
 
-    if (!FixedZLCheckbox->isChecked())
+    if (!FixedZLCheckbox->isChecked() && (!TwoPortMatchingCheckbox->isChecked()))
     {
        inout_operations->loadS1Pdata(LoadFile.toStdString(), LOAD);//s1p
     }
     else
     {
         complex<double> zl_temp;
+        if (!TwoPortMatchingCheckbox->isChecked())
+        {
         QByteArray ba = FixedZLLineedit->text().toLatin1();
         char * text = ba.data();
         zl_temp = getComplexImpedanceFromText(text);
+        }
+        else
+        {
+          zl_temp = inout_operations->AmpS2P.Z0;
+        }
 
         if (zl_temp.real() == -1)//Check if the input value is correct
         {
@@ -330,7 +351,6 @@ void ui::go_clicked()
             return;
         }
         inout_operations->set_constant_ZL_vs_freq(zl_temp);
-
     }
 
     //Check frequency specifications
@@ -343,56 +363,69 @@ void ui::go_clicked()
                              QObject::tr("Wrong frequency settings"));
         return;
     }
-    else//Everything correct... lets set frequency
+
+    //Scale frequency according to the combobox units
+    fmatching_min *= getFreqScale(minFUnitsCombo->currentIndex());
+    fmatching_max *= getFreqScale(maxFUnitsCombo->currentIndex());
+    inout_operations->set_matching_band(fmatching_min, fmatching_max);
+    if (fmatching_min < min(inout_operations->getFrequency())-1e6)//The lower freq is not present at s1p/s2p
     {
-        //Scale frequency according to the combobox units
-        fmatching_min *= getFreqScale(minFUnitsCombo->currentIndex());
-        fmatching_max *= getFreqScale(maxFUnitsCombo->currentIndex());
-        inout_operations->set_matching_band(fmatching_min, fmatching_max);
-        if (fmatching_min < min(inout_operations->getFrequency())-1e6)//The lower freq is not present at s1p/s2p
-        {
-            QMessageBox::warning(0, QObject::tr("Error"),
-                                 QObject::tr("One of the impedance data files does not contain the specified lower frequency"));
-            return;
-        }
-        if (fmatching_max > max(inout_operations->getFrequency())+1e6)//The maximum freq is not present at s1p/s2p
-        {
-            QMessageBox::warning(0, QObject::tr("Error"),
-                                 QObject::tr("One of the impedance data files does not contain the specified upper frequency"));
-            return;
-        }
-        //Check if the specified frequencies lie with the s1p/s2p data
-        inout_operations->ResampleImpedances();//Force data update
+        QMessageBox::warning(0, QObject::tr("Error"),
+                             QObject::tr("One of the impedance data files does not contain the specified lower frequency"));
+        return;
     }
+    if (fmatching_max > max(inout_operations->getFrequency())+1e6)//The maximum freq is not present at s1p/s2p
+    {
+        QMessageBox::warning(0, QObject::tr("Error"),
+                             QObject::tr("One of the impedance data files does not contain the specified upper frequency"));
+        return;
+    }
+
+    //Check if the specified frequencies lie with the s1p/s2p data
+    if (inout_operations->ResampleImpedances() == -2)
+    {
+        QMessageBox::warning(0, QObject::tr("Warning"),
+                             QObject::tr("It is not possible to match the source and load ports simultaneously.\nThe unilateral approach will be used..."));
+    }
+
     inout_operations->UseGNUplot = UseGNUplotCheckbox->isChecked();
 
     GRABIM * MatchingObject = new GRABIM();
     // Impedance and frequency settings
-    MatchingObject->SetSourceImpedance(inout_operations->getSourceImpedance());
-    MatchingObject->SetLoadImpedance(inout_operations->getLoadImpedance());
-    MatchingObject->SetAmplifierS2P(inout_operations->getAmplifierS2P());
-    MatchingObject->SetFrequency(inout_operations->getFrequency());
+    MatchingObject->SetData(inout_operations->getMatchingData());
     MatchingObject->refine = RefineCheckbox->isChecked();
     MatchingObject->setSearchMode(0);//Default mode
 
-    GRABIM_Result R = MatchingObject->RunGRABIM();//Runs GRABIM. Well, this is not exactly the algorithm
+    GRABIM_Result R, RS, RL;
+
+    if (TwoPortMatchingCheckbox->isChecked())
+    {
+         RS = MatchingObject->RunGRABIM(inout_operations->ZS, inout_operations->Zin_maxg);//Input matching network
+         RL = MatchingObject->RunGRABIM(inout_operations->Zout_maxg, inout_operations->ZL);//Output matching network
+         (FixedZSCheckbox->isChecked()) ? RS.source_path = "" : RL.source_path = SourceFile.toStdString();
+         (FixedZLCheckbox->isChecked()) ? RL.load_path = "": RL.load_path = LoadFile.toStdString();
+         TwoPortCircuit TPC;
+         TPC.InputMatching = RS;
+         TPC.OutputMatching = RL;
+         TPC.amplifier_path = AmpFile.toStdString();
+         inout_operations->ExportQucsSchematic(TPC, "");
+
+    }
+    else
+    {
+        R = MatchingObject->RunGRABIM(inout_operations->ZS, inout_operations->ZL);
+        (FixedZSCheckbox->isChecked()) ? R.source_path = "" : R.source_path = SourceFile.toStdString();
+        (FixedZLCheckbox->isChecked()) ? R.load_path = "": R.load_path = LoadFile.toStdString();
 
-    // detailed at [1] but from my point of view is functional and easier to code...
-    //Notes:
-    // 1) The candidate vector is not in log units. I do not see a good reason for doing so. Maybe I am missing something important
-    // 2) Frequency is not in rad/s.
-    // 3) The objective function is the magnitude of S11 expressed in dB. log(x) functions usually have strong
-    // gradients so it seem to suggest that this is good for derivative free opt algorithms
-    // 4) This code takes advantage from NLopt derivative-free local optimisers. NLopt is easy to link and it works
-    // fine. Despite the fact that the Nelder-Mead algorithm does not guarantee convergence (among other problems), it leads to achieve a good local (probably, global) optimum. NM is known to fail when the dimension of the problem is above 14-15. However, as far as matching networks are concerned, the lower order, the simpler the tuning. In practice it is not usual to find matching networks with more than 6 elements
+        TwoPortCircuit TPC;
+        TPC.InputMatching = R;
+        TPC.QucsVersion = PACKAGE_VERSION;
 
-    (FixedZSCheckbox->isChecked()) ? R.source_path = "" : R.source_path = SourceFile.toStdString();
-    (FixedZLCheckbox->isChecked()) ? R.load_path = "": R.load_path = LoadFile.toStdString();
-    R.QucsVersion = PACKAGE_VERSION;
+        inout_operations->tmp_path = "/tmp";
+        if (UseGNUplotCheckbox->isChecked())inout_operations->exportGNUplot(R, GNUplot_path.toStdString(), MatchingObject->refine);
+        inout_operations->ExportQucsSchematic(TPC, "");
+    }
 
-    inout_operations->tmp_path = "/tmp";
-    if (UseGNUplotCheckbox->isChecked())inout_operations->exportGNUplot(R, GNUplot_path.toStdString(), MatchingObject->refine);
-    inout_operations->ExportQucsSchematic(R, "");
     delete MatchingObject;
     delete inout_operations;
    
@@ -551,12 +584,41 @@ void ui::TwoPortMatchingCheckbox_state(int state)
 {
    if (state == Qt::Unchecked)
    { 
+      //Deactivate the S2P button
       imgWidget->setVisible(true);
       AmplifierS2Pbutton->setVisible(false);
+
+      //Show the source and load options
+      (FixedZSCheckbox->isChecked())? FixedZSLineedit->setVisible(true): SourceFileButton->setVisible(true);
+      (FixedZLCheckbox->isChecked())? FixedZLLineedit->setVisible(true): LoadFileButton->setVisible(true);
+
+      FixedZSCheckbox->setVisible(true);
+      FixedZLCheckbox->setVisible(true);
+
+      SourceLabel->setText("Source");
+      LoadLabel->setText("Load");
+
    }
    else
    {
+      //Activate S2P button
       imgWidget->setVisible(false);
       AmplifierS2Pbutton->setVisible(true);
+
+      //Hide terms. Two-port matching can only be done at the impedance
+      //at which the S2P device was measured:
+
+      SourceFileButton->setVisible(false);
+      LoadFileButton->setVisible(false);
+
+      FixedZSCheckbox->setVisible(false);
+      FixedZLCheckbox->setVisible(false);
+
+      FixedZSLineedit->setVisible(false);
+      FixedZLLineedit->setVisible(false);
+
+      SourceLabel->setText("Source\nZ0");
+      LoadLabel->setText("Load\nZ0");
+
    }
 }
diff --git a/qucs/qucs-wideband-matching/ui.h b/qucs/qucs-wideband-matching/ui.h
index 383583ffc1..fa0a992318 100644
--- a/qucs/qucs-wideband-matching/ui.h
+++ b/qucs/qucs-wideband-matching/ui.h
@@ -75,6 +75,7 @@ class ui :public QMainWindow
     QLineEdit * FixedZSLineedit, *FixedZLLineedit;
 
     QLabel *ZLOhmLabel, *ZSOhmLabel;
+    QLabel *SourceLabel, *LoadLabel;
 
     QVBoxLayout * mainLayout, * SourceLayout, * LoadLayout;
     QHBoxLayout * ConstantZSLayout, * ConstantZLLayout, *ButtonsLayout;
@@ -83,7 +84,6 @@ class ui :public QMainWindow
     QGroupBox *FreqgroupBox;
 
 
-
  private slots:
 
     void go_clicked();