diff --git a/dart/src/wNewton.cpp b/dart/src/wNewton.cpp
index 1d801953e86634863fba0c643f8c2c3c633c975f..28efba8418e0a040d7100405ff4a5f9f630188e8 100644
--- a/dart/src/wNewton.cpp
+++ b/dart/src/wNewton.cpp
@@ -87,11 +87,14 @@ bool Newton::run()
     tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
 
     // Display current freestream conditions
-    std::cout << std::setprecision(2)
-              << "- Mach " << pbl->M_inf << ", "
-              << pbl->alpha * 180 / 3.14159 << " deg AoA, "
-              << pbl->beta * 180 / 3.14159 << " deg AoS"
-              << std::endl;
+    if (printOn)
+    {
+      std::cout << std::setprecision(2)
+                << "- Mach " << pbl->M_inf << ", "
+                << pbl->alpha * 180 / 3.14159 << " deg AoA, "
+                << pbl->beta * 180 / 3.14159 << " deg AoS"
+                << std::endl;
+    }
 
     // Initialize solver loop
     nIt = 0;           // iteration counter
@@ -117,21 +120,24 @@ bool Newton::run()
     ls.set(maxLsIt, lsTol, verbose);
 
     // Display residual
-    std::cout << std::setw(8) << "N_Iter"
-              << std::setw(8) << "L_Iter"
-              << std::setw(8) << "f_eval"
-              << std::setw(12) << "Cl"
-              << std::setw(12) << "Cd"
-              << std::setw(15) << "Rel_Res[phi]"
-              << std::setw(15) << "Abs_Res[phi]" << std::endl;
-    std::cout << std::fixed << std::setprecision(5);
-    std::cout << std::setw(8) << nIt
-              << std::setw(8) << 0
-              << std::setw(8) << 1
-              << std::setw(12) << "-"
-              << std::setw(12) << "-"
-              << std::setw(15) << log10(relRes)
-              << std::setw(15) << log10(absRes) << "\n";
+    if (printOn)
+    {
+      std::cout << std::setw(8) << "N_Iter"
+                << std::setw(8) << "L_Iter"
+                << std::setw(8) << "f_eval"
+                << std::setw(12) << "Cl"
+                << std::setw(12) << "Cd"
+                << std::setw(15) << "Rel_Res[phi]"
+                << std::setw(15) << "Abs_Res[phi]" << std::endl;
+      std::cout << std::fixed << std::setprecision(5);
+      std::cout << std::setw(8) << nIt
+                << std::setw(8) << 0
+                << std::setw(8) << 1
+                << std::setw(12) << "-"
+                << std::setw(12) << "-"
+                << std::setw(15) << log10(relRes)
+                << std::setw(15) << log10(absRes) << "\n";
+    }
 
     do
     {
@@ -197,7 +203,7 @@ bool Newton::run()
         Solver::computeLoad();
 
         // Display residual (at each iteration)
-        if (verbose > 0)
+        if (verbose > 0 && printOn)
         {
             std::cout << std::fixed << std::setprecision(5);
             std::cout << std::setw(8) << nIt
@@ -219,7 +225,7 @@ bool Newton::run()
     } while (nIt < maxIt);
 
     // Display residual (only last iteration)
-    if (verbose == 0)
+    if (verbose == 0 && printOn)
     {
         std::cout << std::fixed << std::setprecision(5);
         std::cout << std::setw(8) << nIt
@@ -237,11 +243,13 @@ bool Newton::run()
     // Check the solution
     if (relRes < relTol || absRes < absTol)
     {
-        std::cout << ANSI_COLOR_GREEN << "Newton solver converged!" << ANSI_COLOR_RESET << std::endl;
-        if (verbose > 1)
-            std::cout << "Newton solver CPU" << std::endl
-                      << tms;
-        std::cout << std::endl;
+      if (printOn){
+          std::cout << ANSI_COLOR_GREEN << "Newton solver converged!" << ANSI_COLOR_RESET << std::endl;
+          if (verbose > 1)
+              std::cout << "Newton solver CPU" << std::endl
+                        << tms;
+          std::cout << std::endl;
+      }
         return true;
     }
     else if (std::isnan(relRes))
diff --git a/dart/src/wNewton.h b/dart/src/wNewton.h
index 76cac5358ecb724d91b41afe7777c97382d2f832..f7b632f6090c1aef6216571229f76c5c7e4286b6 100644
--- a/dart/src/wNewton.h
+++ b/dart/src/wNewton.h
@@ -42,6 +42,8 @@ public:
     int maxLsIt;    ///< max number of line search
     double avThrsh; ///< residual threshold to deacrease artificial viscosity
 
+    bool printOn=true; 
+
 private:
     double mCOv; ///< variable cut-off Mach number
     double muCv; ///< variable artificial viscosity scaling factor
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index 5ffef5261bb69d22b908e74f76a21194480fc778..ebdc62d2bb103aa90dcb57f5bf1fe8c2bc59d272 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -109,13 +109,13 @@ void ViscSolver::InitMeshes(size_t region, std::vector<double> locM, std::vector
     std::vector<double> fGlobMesh = meshConverter->CreateFineMesh(globM);
     bl[region]->fineSize = fLocMesh.size();
     bl[region]->SetMesh(fLocMesh, fGlobMesh);
-    std::cout << "- " << bl[region]->name << ": Mesh mapped from " << locM.size() << " to " << bl[region]->GetnMarkers() << " elements." << std::endl;
+    //std::cout << "- " << bl[region]->name << ": Mesh mapped from " << locM.size() << " to " << bl[region]->GetnMarkers() << " elements." << std::endl;
   }
   else
   {
     bl[region]->SetMesh(locM, globM);
     bl[region]->fineSize = locM.size();
-    std::cout << "- " << bl[region]->name << ": Mesh elements " << bl[region]->GetnMarkers() << std::endl;
+    //std::cout << "- " << bl[region]->name << ": Mesh elements " << bl[region]->GetnMarkers() << std::endl;
   }
 }
 
@@ -134,7 +134,7 @@ void ViscSolver::SendInvBC(size_t region, std::vector<double> _Ue,
     std::vector<double> RhoeFine = meshConverter->InterpolateToFineMesh(_Rhoe, bl[region]->fineSize);
 
     bl[region]->SetInvBC(UeFine, MeFine, RhoeFine);
-    std::cout << "- " << bl[region]->name << ": Inviscid boundary interpolated." << std::endl;
+    //std::cout << "- " << bl[region]->name << ": Inviscid boundary interpolated." << std::endl;
   }
   else
   {
@@ -174,12 +174,16 @@ int ViscSolver::Run(unsigned int couplIter){
     tSolver->SetinitSln(true);
     if (couplIter != 0 && (nbElmUpper != bl[0]->GetnMarkers())){
       stagPtMvmt = true;
-      std::cout << "Stagnation point moved" << std::endl;
+      if (printOn){
+        std::cout << "Stagnation point moved" << std::endl;
+      }
     }
     else{
       stagPtMvmt = false;
     }
-    std::cout << "Restart solution" << std::endl;
+    if (printOn){
+      std::cout << "Restart solution" << std::endl;
+    }
   }
 
   else{
@@ -187,7 +191,9 @@ int ViscSolver::Run(unsigned int couplIter){
     tSolver->SetitFrozenJac(5);
     tSolver->SetinitSln(false);
     stagPtMvmt = false;
-    std::cout << "Updating solution" << std::endl;
+    if (printOn){
+      std::cout << "Updating solution" << std::endl;
+    }
   }
 
   nbElmUpper = bl[0]->GetnMarkers();
@@ -218,7 +224,9 @@ int ViscSolver::Run(unsigned int couplIter){
     }
     lockTrans = false;
 
-    std::cout << "- " << bl[iRegion]->name << " region";
+    if (printOn){
+      std::cout << "- " << bl[iRegion]->name << " region";
+    }
 
     if (iRegion == 2){
       /* Impose wake boundary condition */
@@ -227,7 +235,9 @@ int ViscSolver::Run(unsigned int couplIter){
       UpperCond = std::vector<double>(bl[0]->U.end() - bl[0]->GetnVar(), bl[0]->U.end()); /* Base std::vector constructor to slice */
       LowerCond = std::vector<double>(bl[1]->U.end() - bl[1]->GetnVar(), bl[1]->U.end()); /* Base std::vector constructor to slice */
       bl[iRegion]->SetWakeBC(Re, UpperCond, LowerCond);
-      std::cout << "\n" << std::endl;
+      if (printOn){
+        std::cout << "\n" << std::endl;
+      }
 
       lockTrans = true;
     }
@@ -245,7 +255,9 @@ int ViscSolver::Run(unsigned int couplIter){
       pointExitCode = tSolver->Integration(iPoint, bl[iRegion]);
       if (pointExitCode!=0)
       {
-        std::cout << "Point " << iPoint << ": Convergence terminated with exitcode " << pointExitCode << std::endl;
+        if (printOn){
+          std::cout << "Point " << iPoint << ": Convergence terminated with exitcode " << pointExitCode << std::endl;
+        }
         solverExitCode = -1; // Convergence failed
       }
 
@@ -258,7 +270,9 @@ int ViscSolver::Run(unsigned int couplIter){
         // Free transition.
         if (bl[iRegion]->U[iPoint*bl[iRegion]->GetnVar() + 2] >= bl[iRegion]->GetnCrit()){
           AverageTransition(iPoint, bl[iRegion], 0);
-          std::cout << ", free transition x/c = " << bl[iRegion]->xtr << ". Marker " << bl[iRegion]->transMarker << std::endl;
+          if (printOn){
+            std::cout << ", free transition x/c = " << bl[iRegion]->xtr << ". Marker " << bl[iRegion]->transMarker << std::endl;
+          }
           lockTrans = true;
         } // End if
       } // End if (!lockTrans).
@@ -393,7 +407,7 @@ std::vector<double> ViscSolver::ExtractDeltaStar(size_t iRegion)
   std::vector<double> cDStar;
   if (mapMesh)
   { 
-    std::cout << "DeltaStar size" << bl[iRegion]->DeltaStar.size() << std::endl;
+    // std::cout << "DeltaStar size" << bl[iRegion]->DeltaStar.size() << std::endl;
     cDStar = meshConverter->VarsFineToCoarse(bl[iRegion]->DeltaStar, bl[iRegion]->coarseSize);
   }
   else
diff --git a/dart/src/wViscSolver.h b/dart/src/wViscSolver.h
index a8286b46153e204b0b5658c40e950e55537fe3f1..e5c72884c7712fdd382eb35cb238aa2e5e4eeaf0 100644
--- a/dart/src/wViscSolver.h
+++ b/dart/src/wViscSolver.h
@@ -13,9 +13,9 @@ public:
   double Cdt;
   double Cdf;
   double Cdp;
+  double Re;
 
 private:
-  double Re;
   double Minf;
   double CFL0;
   size_t nbElmUpper;
@@ -28,6 +28,8 @@ private:
 public:
   std::vector<BLRegion *> bl;
 
+  bool printOn = true;
+
   ViscSolver(double _Re, double _Minf, double _CFL0, unsigned int meshFactor);
 
   ~ViscSolver();
diff --git a/dart/tests/blicpp.py b/dart/tests/blicpp.py
index 3b756977788a3a050226708880c302a04bb748dc..a2eefe0be72b588f15ecf2004dab1b75f08ed966 100644
--- a/dart/tests/blicpp.py
+++ b/dart/tests/blicpp.py
@@ -39,6 +39,7 @@ import dart.default as floD
 import dart.viscUtils as viscU
 
 import dart.vCoupler as floVC
+import numpy as np
 
 import tbox
 import tbox.utils as tboxU
@@ -56,10 +57,11 @@ def main():
 
     # define flow variables
     Re = 1e7
-    alpha = 15*math.pi/180
-    M_inf = 0.2
+    alpha = 2*math.pi/180
+    alphaSeq = np.linspace(0, 10, 10)
+    M_inf = 0.
 
-    meshFactor = 2
+    meshFactor = 1
     CFL0 = 1
 
     plotVar = 'Cf'
@@ -72,7 +74,7 @@ def main():
     # mesh the geometry
     print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
     tms['msh'].start()
-    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.0001}
+    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
     msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
     tms['msh'].stop()
 
@@ -93,7 +95,7 @@ def main():
 
     coupler.CreateProblem(blw[0], blw[1])
 
-    coupler.Run()
+    coupler.RunPolar(alphaSeq)
     tms['solver'].stop()
 
     # extract Cp
diff --git a/dart/tests/blicppPolar.py b/dart/tests/blicppPolar.py
new file mode 100644
index 0000000000000000000000000000000000000000..63a5ebbfaab8dac39a97c49692f33b8d7e2d64a6
--- /dev/null
+++ b/dart/tests/blicppPolar.py
@@ -0,0 +1,147 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2020 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+## Compute lifting (linear or nonlinear) viscous flow around a NACA 0012
+# Amaury Bilocq
+#
+# Test the viscous-inviscid interaction scheme
+# Reference to the master's thesis: http://hdl.handle.net/2268/252195
+# Reference test cases with Naca0012 (different from master's thesis):
+# 1) Incompressible: Re = 1e7, M_inf = 0, alpha = 5°, msTE = 0.01, msLE = 0.01
+#    -> nIt = 27, Cl = 0.55, Cd = 0.0058, xtrTop = 0.087, xtrBot = 0.741
+#    -> msLe = 0.001, xtrTop = 0.061
+# 2) Compressible: Re = 1e7, M_inf = 5, alpha = 5°, msTE = 0.01, msLE = 0.01
+#    -> nIt = 33, Cl = 0.65, Cd = 0.0063, xtrTop = 0.057, xtrBot = 0.740
+# 3) Separated: Re = 1e7, M_inf = 0, alpha = 12°, msTE = 0.01, msLE = 0.001
+#    -> nIt = 43, Cl = 1.30 , Cd = 0.0108, xtrTop = 0.010, xtrBot = 1
+#
+# CAUTION
+# This test is provided to ensure that the solver works properly.
+# Mesh refinement may have to be performed to obtain physical results.
+
+import dart.utils as floU
+import dart.default as floD
+import dart.viscUtils as viscU
+
+import dart.vCoupler as floVC
+
+import tbox
+import tbox.utils as tboxU
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+import numpy as np
+
+import math
+
+def main():
+    print('Start')
+    # timer
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    # define flow variables
+    Re = 1e7
+    alpha = np.linspace(-5, 5, 1)*math.pi/180
+    M_inf = 0.
+
+    meshFactor = 1
+    CFL0 = 1
+
+    plotVar = 'Cf'
+
+    # ========================================================================================
+
+    c_ref = 1
+    dim = 2
+
+    # mesh the geometry
+    print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+    tms['msh'].start()
+    pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.0001}
+    msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
+    tms['msh'].stop()
+
+    # set the problem and add medium, initial/boundary conditions
+
+    tms['pre'].start()
+    pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha[0], 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
+    tms['pre'].stop()
+
+    # solve viscous problem
+
+    print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
+    tms['solver'].start()
+    
+    isolver = floD.newton(pbl)
+    vsolver = floD.initBL(Re, M_inf, CFL0, meshFactor)
+
+    coupler = floVC.Coupler(isolver, vsolver, alpha)
+
+    coupler.CreateProblem(blw[0], blw[1])
+
+    coupler.Run()
+    tms['solver'].stop()
+
+    # extract Cp
+    Cp = floU.extract(bnd.groups[0].tag.elems, isolver.Cp)
+    tboxU.write(Cp, 'Cp_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
+    # display results
+    print(ccolors.ANSI_BLUE + 'PyRes...' + ccolors.ANSI_RESET)
+    print('      Re        M    alpha       Cl       Cd      Cdp      Cdf       Cm')
+    print('{0:6.1f}e6 {1:8.2f} {2:8.1f} {3:8.4f} {4:8.4f} {5:8.4f} {6:8.4f} {7:8.4f}'.format(Re/1e6, M_inf, alpha*180/math.pi, isolver.Cl, vsolver.Cdt, vsolver.Cdp, vsolver.Cdf, isolver.Cm))
+
+    # display timers
+    tms['total'].stop()
+    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+    print('CPU statistics')
+    print(tms)
+    
+    
+    
+    # visualize solution and plot results
+    #floD.initViewer(pbl)
+    tboxU.plot(Cp[:,0], Cp[:,3], 'x', 'Cp', 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(isolver.Cl, vsolver.Cdt, isolver.Cm, 4), True)
+
+    blScalars, blVectors = viscU.ExtractVars(vsolver)
+    xtr=vsolver.Getxtr()
+    wData = viscU.Dictionary(blScalars, blVectors, xtr)
+    viscU.PlotVars(wData, plotVar)
+
+    
+    # check results
+    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+    tests = CTests()
+    if Re == 1e7 and M_inf == 0.2 and alpha == 5*math.pi/180:
+        tests.add(CTest('Cl', isolver.Cl, 0.56, 5e-2)) # XFOIL 0.58
+        tests.add(CTest('Cd', vsolver.Cdt, 0.0063, 1e-3, forceabs=True)) # XFOIL 0.00617
+        tests.add(CTest('Cdp', vsolver.Cdp, 0.0019, 1e-3, forceabs=True)) # XFOIL 0.00176
+        tests.add(CTest('xtr_top', vsolver.xtr[0], 0.059, 0.03, forceabs=True)) # XFOIL 0.0510
+        tests.add(CTest('xtr_bot', vsolver.xtr[1], 0.738, 0.03, forceabs=True)) # XFOIL 0.7473
+    else:
+        raise Exception('Test not defined for this flow')
+
+    tests.run()
+
+
+    # eof
+    print('')
+
+if __name__ == "__main__":
+    main()
diff --git a/dart/vCoupler.py b/dart/vCoupler.py
index a93f9ece5f9cf530d8ad0f3c90820a3cdac1db59..6224e6e078b8df1094e5f38468a934af5e31c1de 100644
--- a/dart/vCoupler.py
+++ b/dart/vCoupler.py
@@ -19,11 +19,13 @@
 ## VII Coupler
 # Paul Dechamps
 
+from multiprocessing import set_forkserver_preload
 from dart.viscous.Master.DAirfoil import Airfoil
 from dart.viscous.Master.DWake import Wake
 from dart.viscous.Drivers.DGroupConstructor import GroupConstructor
 
 from fwk.coloring import ccolors
+import math
 
 import fwk
 import dart.viscUtils as viscU
@@ -56,8 +58,11 @@ class Coupler:
       nElmUpperPrev = 0
       couplIter = 0
       cdPrev = 0
+      print('- AoA: {:<2.2f} Mach: {:<1.1f} Re: {:<2.1f}e6'
+      .format(self.iSolver.pbl.alpha*180/math.pi, self.iSolver.pbl.M_inf, self.vSolver.Re/1e6))
+      print('{:>5s}| {:>10s} {:>10s} | {:>8s} {:>6s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
+      print(' --------------------------------------------------------')
       while couplIter < self.maxCouplIter:
-        
 
         # Run inviscid solver
         self.tms['inv'].start()
@@ -83,7 +88,7 @@ class Coupler:
 
         if ((len(fMeshDict['locCoord'][:fMeshDict['bNodes'][1]]) != nElmUpperPrev) or couplIter==0):
 
-          print('STAG POINT MOVED')
+          # print('STAG POINT MOVED')
           
           # Initialize mesh          
           self.vSolver.InitMeshes(0, fMeshDict['locCoord'][:fMeshDict['bNodes'][1]], fMeshDict['globCoord'][:fMeshDict['bNodes'][1]])
@@ -104,13 +109,13 @@ class Coupler:
         self.vSolver.SendInvBC(0, fMeshDict['vtExt'][:fMeshDict['bNodes'][1]], fMeshDict['Me'][:fMeshDict['bNodes'][1]], fMeshDict['rho'][:fMeshDict['bNodes'][1]])
         self.vSolver.SendInvBC(1, fMeshDict['vtExt'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['Me'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]], fMeshDict['rho'][fMeshDict['bNodes'][1]:fMeshDict['bNodes'][2]])
         self.vSolver.SendInvBC(2, fMeshDict['vtExt'][fMeshDict['bNodes'][2]:], fMeshDict['Me'][fMeshDict['bNodes'][2]:], fMeshDict['rho'][fMeshDict['bNodes'][2]:])
-        print('Inviscid boundary imposed')
+        # print('Inviscid boundary imposed')
         # Run viscous solver
 
         self.tms['visc'].start()
         exitCode = self.vSolver.Run(couplIter)
         self.tms['visc'].stop()
-        print('Viscous ops terminated with exit code ', exitCode)
+        # print('Viscous ops terminated with exit code ', exitCode)
 
         error = abs((self.vSolver.Cdt - cdPrev)/self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1
         
@@ -124,19 +129,27 @@ class Coupler:
 
         if error  <= self.couplTol and exitCode==0:
 
-          print(self.tms)
-          print('')
+          #print(self.tms)
+          print(ccolors.ANSI_GREEN,'{:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.2f} {:>6.2f} {:>8.2f}\n'.format(couplIter, self.iSolver.Cl, self.vSolver.Cdt, xtr[0], xtr[1], np.log10(error)),ccolors.ANSI_RESET)
+
+          """print('')
           print(ccolors.ANSI_GREEN, 'It: {:<3.0f} Cd = {:<6.5f}. Error = {:<2.3f} < {:<2.3f}'
           .format(couplIter, self.vSolver.Cdt, np.log10(error), np.log10(self.couplTol)), ccolors.ANSI_RESET)
-          print('')
+          print('')"""
           return 0
         
         cdPrev = self.vSolver.Cdt
+
+        xtr=self.vSolver.Getxtr()
+
+      #print('{:>5s} {:>12s} {:>12s} {:>12s} {:>6s} {:>12s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
+        if couplIter%5==0:
+          print(' {:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.2f} {:>6.2f} {:>8.2f}'.format(couplIter, self.iSolver.Cl, self.vSolver.Cdt, xtr[0], xtr[1], np.log10(error)))
         
-        print('')
+        """print('')
         print(ccolors.ANSI_BLUE, 'It: {:<3.0f} Cd = {:<6.5f}. Error = {:<2.3f} > {:<2.3f}'
         .format(couplIter, self.vSolver.Cdt, np.log10(error), np.log10(self.couplTol)), ccolors.ANSI_RESET)
-        print('')
+        print('')"""
 
         # Retreive and set blowing velocities.
 
@@ -198,3 +211,30 @@ class Coupler:
         print(ccolors.ANSI_RED, 'Maximum number of {:<.0f} coupling iterations reached'.format(self.maxCouplIter), ccolors.ANSI_RESET)
         return couplIter
 
+
+    def RunPolar(self, alphaSeq):
+      
+      import math
+
+      self.iSolver.printOn = False
+      self.vSolver.printOn = False
+
+      coeffTable = np.empty((2,len(alphaSeq)))
+
+      for ind, alpha in enumerate(alphaSeq):
+        self.iSolver.pbl.update(alpha*math.pi/180)
+
+        self.Run()
+
+        # Collect Coeff
+
+        coeffTable[0,ind] = self.iSolver.Cl
+        coeffTable[1,ind] = self.vSolver.Cdt
+
+      plt.plot(coeffTable[0,:], coeffTable[1,:], 'x-')
+      plt.xlabel('$c_l$ (-)')
+      plt.ylabel('$c_d$ (-)')
+      plt.show()
+
+
+