diff --git a/Plugin/NearToFarField.cpp b/Plugin/NearToFarField.cpp
index 2632a7d3897c8cd5ad20a9d4d561a002cb57da2e..99f1d9df6541fc24c78122cef2f825c2e8355c36 100644
--- a/Plugin/NearToFarField.cpp
+++ b/Plugin/NearToFarField.cpp
@@ -4,10 +4,6 @@
 // bugs and problems to <gmsh@geuz.org>.
 
 #include "NearToFarField.h"
-#include "Numeric.h"
-#include "PViewOptions.h"
-#include "MElement.h"
-#include "GModel.h"
 
 StringXNumber NearToFarFieldOptions_Number[] = {
   {GMSH_FULLRC, "Wavenumber",     NULL, 1.},
@@ -50,127 +46,70 @@ StringXNumber *GMSH_NearToFarFieldPlugin::getOption(int iopt)
   return &NearToFarFieldOptions_Number[iopt];
 }
 
-void GMSH_NearToFarFieldPlugin::CartesianToSpherical(int numSteps, double theta, double phi, double **Fc, double **Fsp)
-{
-  double sTheta = sin(theta) ;
-  double cTheta = cos(theta) ;
-  double sPhi   = sin(phi) ;
-  double cPhi   = cos(phi) ;
-
-   for(int step = 0; step < numSteps; step++){
-     Fsp[step][0] = Fc[step][0] * sTheta * cPhi + Fc[step][1] * sTheta * sPhi + Fc[step][2] *cTheta ; 
-     Fsp[step][1] = Fc[step][0] * cTheta * cPhi + Fc[step][1] * cTheta * sPhi - Fc[step][2]* sTheta ; 
-     Fsp[step][2] =-Fc[step][0] * sPhi          + Fc[step][1] * cPhi ;  
-   }
-}
 
-double GMSH_NearToFarFieldPlugin::getFarField(PViewData *eData, PViewData *hData, double k0, double r_far, double theta, double phi)
+double GMSH_NearToFarFieldPlugin::getFarField(std::vector<element*> allElems,
+                                              std::vector<std::vector<double> > js, std::vector<std::vector<double> > ms,
+                                              double k0, double r_far, double theta, double phi)
 {
 
-  // theta in [0, pi] (elevation/polar angle)  
+  // theta in [0, pi] (elevation/polar angle)
   // phi in [0, 2*pi] (azimuthal angle)
- 
-  double r[3] = { sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta) }; // Unit vector position
-  double Z0 = 120 * M_PI ; // free-space impedance
-
-  int numSteps = eData->getNumTimeSteps() ;
-  int numEntities = eData->getNumEntities(0) ;
-
-  double **N  = new double * [numSteps] ;
-  double **Ns = new double * [numSteps] ;
-  double **L  = new double * [numSteps] ;
-  double **Ls = new double * [numSteps] ;
-  for(int step = 0; step < numSteps; step++){
-    N [step] = new double[3] ;
-    Ns[step] = new double[3] ;
-    L [step] = new double[3] ;
-    Ls[step] = new double[3] ;
-  }
-
-  for (int step = 0; step < numSteps; step++)
-    for(int comp = 0; comp < 3; comp++){
-      N[step][comp]= Ns[step][comp]= 0. ; 
-      L[step][comp]= Ls[step][comp]= 0. ;
-    }
- 
-  // tag all the nodes with "0" (the default tag)
-  for(int step = 0; step < numSteps; step++){
-    for(int ent = 0; ent < numEntities; ent++){
-      for(int ele = 0; ele < eData->getNumElements(step, ent); ele++){
-        if(eData->skipElement(step, ent, ele)) continue;
-        for(int nod = 0; nod < eData->getNumNodes(step, ent, ele); nod++)
-          eData->tagNode(step, ent, ele, nod, 0);
-      }
-    }
-  }
-
-  for(int ent = 0; ent < eData->getNumEntities(0); ent++){
-    for(int ele = 0; ele < eData->getNumElements(0, ent); ele++){
-      if(eData->skipElement(0, ent, ele)) continue;
-      int numComp  = eData->getNumComponents(0, ent, ele);
-      if(numComp != 3) continue ;
-      int numNodes = eData->getNumNodes(0, ent, ele);
-      std::vector<double> x(numNodes), y(numNodes), z(numNodes);
-      std::vector<int> tag(numNodes);
 
-      for(int nod = 0; nod < numNodes; nod++)
-        tag[nod] = eData->getNode(0, ent, ele, nod, x[nod], y[nod], z[nod]);
-      double n[3] = {0.,0.,0.};
-      normal3points(x[0], y[0], z[0], x[1], y[1], z[1], x[2], y[2], z[2], n); 
+  double sTheta = sin(theta) ; double cTheta = cos(theta) ;
+  double sPhi = sin(phi) ; double cPhi = cos(phi) ;
+  double r[3] = {sTheta*cPhi, sTheta*sPhi, cTheta}; // Unit vector position
 
-      Msg::Error("FIXME JS and MS computation need to be recoded");
-      // cannot allocate like this + fix numSteps
-      /* 
-      double Js[numSteps][numNodes*numComp], Ms[numSteps][numNodes*numComp] ;
+  double Z0 = 120 * M_PI ; // free-space impedance
 
-      for(int step = 0; step < numSteps; step++){
-        for(int nod = 0; nod < numNodes; nod++){
-          if(tag[nod]) continue ; // already condisered in integration
-          for(int comp = 0; comp < numComp; comp++){
-            eData->getValue(step, ent, ele, nod, comp, Ms[numSteps][numComp*nod + comp ]);
-            hData->getValue(step, ent, ele, nod, comp, Js[numSteps][numComp * nod + comp]);
-          }
+  int numComps = 3, numSteps = 2 ;
+  std::vector<std::vector<double> > N(numSteps,numComps), Ns(numSteps,numComps) ;
+  std::vector<std::vector<double> > L(numSteps,numComps), Ls(numSteps,numComps) ;
+
+  int i = 0 ;
+  for(unsigned int ele = 0; ele < allElems.size(); ele++){
+    element* e = allElems[ele] ;
+    int numNodes = e->getNumNodes() ;
+    int numEdges = e->getNumEdges() ;
+
+    std::vector<double > valN0(numNodes*numComps),valN1(numNodes*numComps), valL0(numNodes*numComps), valL1(numNodes*numComps) ;
+
+    for(int nod = 0; nod < numNodes; nod++){
+      double x, y, z;
+      e->getXYZ(nod, x, y, z) ;
+      double rr, r_nod[3] = {x, y, z};
+      prosca(r_nod, r, &rr) ;
+      double e_jk0rr[2] = {cos(k0*rr), sin(k0*rr)} ;
+
+      for(int comp = 0; comp < numComps; comp++){
+        if(i < js[0].size()){
+          valN0[numComps * nod + comp] = js[0][i] * e_jk0rr[0] - js[1][i] * e_jk0rr[1];
+          valN1[numComps * nod + comp] = js[0][i] * e_jk0rr[1] + js[1][i] * e_jk0rr[0];
+          valL0[numComps * nod + comp] = ms[0][i] * e_jk0rr[0] - ms[1][i] * e_jk0rr[1];
+          valL1[numComps * nod + comp] = ms[0][i] * e_jk0rr[1] + ms[1][i] * e_jk0rr[0];
+          i++;
         }
       }
+    }
 
-      // Integration 
-      double P0[3] = {x[0], y[0], z[0]} ;
-      double P1[3] = {x[1], y[1], z[1]} ;
-      double P2[3] = {x[2], y[2], z[2]} ;
-      double quad_area = triangle_area(P0,P1,P2);
-      
-      for(int nod = 0; nod < numNodes; nod++){
-        double rr, r_nod[3] = {x[nod], y[nod], z[nod]}; 
-        prosca(r_nod, r, &rr) ;
-        double cos_k0rr = quad_area*cos(k0*rr) ; 
-        double sin_k0rr = quad_area*sin(k0*rr) ;
-        
-        N[0][0] += Js[0][numComp*nod + 0] * cos_k0rr - Js[1][numComp*nod + 0] * sin_k0rr  ;   
-        N[0][1] += Js[0][numComp*nod + 1] * cos_k0rr - Js[1][numComp*nod + 1] * sin_k0rr  ;   
-        N[0][2] += Js[0][numComp*nod + 2] * cos_k0rr - Js[1][numComp*nod + 2] * sin_k0rr  ;   
-
-        N[1][0] += Js[0][numComp*nod + 0] * sin_k0rr + Js[1][numComp*nod + 0] * cos_k0rr  ;   
-        N[1][1] += Js[0][numComp*nod + 1] * sin_k0rr + Js[1][numComp*nod + 1] * cos_k0rr  ;   
-        N[1][2] += Js[0][numComp*nod + 2] * sin_k0rr + Js[1][numComp*nod + 2] * cos_k0rr  ;   
-
-        L[0][0] += Ms[0][numComp*nod + 0] * cos_k0rr - Ms[1][numComp*nod + 0] * sin_k0rr  ;   
-        L[0][1] += Ms[0][numComp*nod + 1] * cos_k0rr - Ms[1][numComp*nod + 1] * sin_k0rr  ;   
-        L[0][2] += Ms[0][numComp*nod + 2] * cos_k0rr - Ms[1][numComp*nod + 2] * sin_k0rr  ;   
-
-        L[1][0] += Ms[0][numComp*nod + 0] * sin_k0rr + Ms[1][numComp*nod + 0] * cos_k0rr  ;   
-        L[1][1] += Ms[0][numComp*nod + 1] * sin_k0rr + Ms[1][numComp*nod + 1] * cos_k0rr  ;   
-        L[1][2] += Ms[0][numComp*nod + 2] * sin_k0rr + Ms[1][numComp*nod + 2] * cos_k0rr  ;   
-
-        eData->tagNode(0, ent, ele, nod, 1);
-      }
-      */
+    N[0][0] += e->integrate(&valN0[0], 3) ; N[1][0] += e->integrate(&valN1[0], 3) ;
+    N[0][1] += e->integrate(&valN0[1], 3) ; N[1][1] += e->integrate(&valN1[1], 3) ;
+    N[0][2] += e->integrate(&valN0[2], 3) ; N[1][2] += e->integrate(&valN1[2], 3) ;
 
-    }
-  }
-  
+    L[0][0] += e->integrate(&valL0[0], 3) ; L[1][0] += e->integrate(&valL1[0], 3) ;
+    L[0][1] += e->integrate(&valL0[1], 3) ; L[1][1] += e->integrate(&valL1[1], 3) ;
+    L[0][2] += e->integrate(&valL0[2], 3) ; L[1][2] += e->integrate(&valL1[2], 3) ;
+ }
 
-  CartesianToSpherical(numSteps, theta, phi, N, Ns) ;
-  CartesianToSpherical(numSteps, theta, phi, L, Ls) ;
+  // From Cartesian to spherical coordinates
+  for(int step = 0; step < 2; step++){
+    Ns[step][0] = N[step][0] * sTheta * cPhi + N[step][1] * sTheta * sPhi + N[step][2] *cTheta ;
+    Ns[step][1] = N[step][0] * cTheta * cPhi + N[step][1] * cTheta * sPhi - N[step][2]* sTheta ;
+    Ns[step][2] =-N[step][0] * sPhi          + N[step][1] * cPhi ;
+
+    Ls[step][0] = L[step][0] * sTheta * cPhi + L[step][1] * sTheta * sPhi + L[step][2] *cTheta ;
+    Ls[step][1] = L[step][0] * cTheta * cPhi + L[step][1] * cTheta * sPhi - L[step][2]* sTheta ;
+    Ls[step][2] =-L[step][0] * sPhi          + L[step][1] * cPhi ;
+   }
 
   // E_r radial component is negligible in far field
   double E_theta[2] ;
@@ -180,32 +119,17 @@ double GMSH_NearToFarFieldPlugin::getFarField(PViewData *eData, PViewData *hData
   double sin_k0r = sin(k0*r_far) ;
 
   // Elevation component
-  E_theta[0] = -k0_over_4pir * ( (Ls[0][2] + Z0 * Ns[0][1]) * sin_k0r -  
+  E_theta[0] = -k0_over_4pir * ( (Ls[0][2] + Z0 * Ns[0][1]) * sin_k0r -
                                  (Ls[1][2] + Z0 * Ns[1][1]) * cos_k0r ) ;
-  E_theta[1] = -k0_over_4pir * ( (Ls[0][2] + Z0 * Ns[0][1]) * cos_k0r +  
+  E_theta[1] = -k0_over_4pir * ( (Ls[0][2] + Z0 * Ns[0][1]) * cos_k0r +
                                  (Ls[1][2] + Z0 * Ns[1][1]) * sin_k0r ) ;
-  // Azimuthal component 
-  E_phi[0]   =  k0_over_4pir * ( (Ls[0][1] - Z0 * Ns[0][2]) * sin_k0r -  
+  // Azimuthal component
+  E_phi[0]   =  k0_over_4pir * ( (Ls[0][1] - Z0 * Ns[0][2]) * sin_k0r -
                                  (Ls[1][1] - Z0 * Ns[1][2]) * cos_k0r ) ;
-  E_phi[1]   =  k0_over_4pir * ( (Ls[0][1] - Z0 * Ns[0][2]) * cos_k0r +  
+  E_phi[1]   =  k0_over_4pir * ( (Ls[0][1] - Z0 * Ns[0][2]) * cos_k0r +
                                  (Ls[1][1] - Z0 * Ns[1][2]) * sin_k0r ) ;
 
-
-  //printf("Ephi %g %g \n ", E_phi[0], E_phi[1]) ;
-  //printf("Etheta %g %g\n ", E_theta[0], E_theta[1]) ;
-
- double farF =  1./2./Z0 * ( (E_theta[0]*E_theta[0] + E_theta[1]*E_theta[1]) + (E_phi[0]*E_phi[0]+E_phi[1] *E_phi[1]) ) ;
-
-  for (int step = 0; step < numSteps; step++){
-      delete [] N [step]  ; 
-      delete [] Ns[step] ; 
-      delete [] L [step]  ; 
-      delete [] Ls[step] ; 
-  }
-  delete [] N  ; 
-  delete [] Ns ; 
-  delete [] L  ; 
-  delete [] Ls ;   
+  double farF =  1./2./Z0 * ( (E_theta[0]*E_theta[0] + E_theta[1]*E_theta[1]) + (E_phi[0]*E_phi[0]+E_phi[1] *E_phi[1]) ) ;
 
   return farF ;
 }
@@ -225,8 +149,8 @@ PView *GMSH_NearToFarFieldPlugin::execute(PView * v)
   PView *ve = getView(_eView, v);
   if(!ve){
      Msg::Error("NearToFarField plugin could not find EView %i", _eView);
-    return v; 
-  } 
+    return v;
+  }
   PView *vh = getView(_hView, v);
   if(!vh){
     Msg::Error("NearToFarField plugin could not find HView %i", _hView);
@@ -242,8 +166,8 @@ PView *GMSH_NearToFarFieldPlugin::execute(PView * v)
     return v;
   }
 
-  if(eData->getNumTimeSteps()!= 2){
-    Msg::Error("NearToFarField Plugin only implemented for frequency domain, fields must be complex");
+  if(eData->getNumTimeSteps()!= 2 || hData->getNumTimeSteps() != 2){
+    Msg::Error("Invalid number of steps for EView or HView, fields must be complex");
     return v;
   }
 
@@ -257,95 +181,83 @@ PView *GMSH_NearToFarFieldPlugin::execute(PView * v)
     Msg::Error("EView %i and HView %i must be given on the same grid", _eView, _hView);
     return v;
   }
-  
 
   // View for far field: represented on a sphere of radious determined by the size of the BoundingBox
   PView *vf = new PView();
   PViewDataList *dataFar = getDataList(vf);
-   
+
+  std::vector<element*> allElems  ;
+  std::vector<std::vector<double> > js ;
+  std::vector<std::vector<double> > ms ;
+
+  int numSteps = eData->getNumTimeSteps() ;
+  js.resize(numSteps);
+  ms.resize(numSteps);
+
+  for(int ent = 0; ent < eData->getNumEntities(0); ent++){
+    for(int ele = 0; ele < eData->getNumElements(0, ent); ele++){
+      if(eData->skipElement(0, ent, ele)) continue;
+      if(hData->skipElement(0, ent, ele)) continue;
+      int numComp  = eData->getNumComponents(0, ent, ele);
+      if(numComp != 3) continue ;
+
+      int dim = eData->getDimension(0, ent, ele);
+      int numNodes = eData->getNumNodes(0, ent, ele);
+      std::vector<double> x(numNodes), y(numNodes), z(numNodes);
+      for(int nod = 0; nod < numNodes; nod++)
+        eData->getNode(0, ent, ele, nod, x[nod], y[nod], z[nod]);
+
+      elementFactory factory;
+      allElems.push_back(factory.create(numNodes, dim, &x[0], &y[0], &z[0], true));
+
+      double n[3] = {0.,0.,0.};
+      if(numNodes>2)
+        normal3points(x[0], y[0], z[0], x[1], y[1], z[1], x[2], y[2], z[2], n);
+      else
+        normal2points(x[0], y[0], z[0], x[1], y[1], z[1], n);
+
+      for(int step = 0; step < numSteps; step++){
+        for(int nod = 0; nod < numNodes; nod++){
+          double h[3], e[3];
+          for(int comp = 0; comp < numComp; comp++){
+            eData->getValue(step, ent, ele, nod, comp, e[comp]);
+            hData->getValue(step, ent, ele, nod, comp, h[comp]);
+          }
+          double j[3], m[3] ;
+          prodve(n, h, j) ; // Js =   n x H ; Surface electric current
+          prodve(e, n, m) ; // Ms = - n x E ; Surface magnetic current
+          js[step].push_back(j[0]) ; js[step].push_back(j[1]) ; js[step].push_back(j[2]) ;
+          ms[step].push_back(m[0]) ; ms[step].push_back(m[1]) ; ms[step].push_back(m[2]) ;
+        }
+      }
+
+    }
+  }
+
+  // -------------------------------------------------------------
+  // Generating radiation pattern
+  // -------------------------------------------------------------
   double phi,   dPhi   = 2*M_PI/_NbPhi ;
   double theta, dTheta =   M_PI/_NbThe ;
   double ffmax = 0.0 ;
 
-  double **allPhi = new double *[_NbPhi+1] ;
-  double **allThe = new double *[_NbPhi+1] ;
-  double **farF   = new double *[_NbPhi+1] ;
-  for (int i = 0; i <= _NbPhi; i++){
-    allPhi[i] = new double [_NbThe+1] ;
-    allThe[i] = new double [_NbThe+1] ;
-    farF[i]   = new double [_NbThe+1] ;
-  }
-
- ve->setChanged(true);
- vh->setChanged(true);
-
- for(int step = 0; step < eData->getNumTimeSteps(); step++){
-   // tag all the nodes with "0" (the default tag)
-   for(int ent = 0; ent < eData->getNumEntities(step); ent++){
-     for(int ele = 0; ele < eData->getNumElements(step, ent); ele++){
-       if(eData->skipElement(step, ent, ele)) continue;
-       if(hData->skipElement(step, ent, ele)) continue;
-       for(int nod = 0; nod < eData->getNumNodes(step, ent, ele); nod++){
-         eData->tagNode(step, ent, ele, nod, 0);
-         hData->tagNode(step, ent, ele, nod, 0);
-       }
-     }
-   }
-   
-   for(int ent = 0; ent < eData->getNumEntities(step); ent++){
-     for(int ele = 0; ele < eData->getNumElements(step, ent); ele++){
-       if(eData->skipElement(0, ent, ele)) continue;
-       if(hData->skipElement(0, ent, ele)) continue;
-       int numComp  = eData->getNumComponents(0, ent, ele);
-       if(numComp != 3) continue ;
-       int numNodes = eData->getNumNodes(0, ent, ele);
-       std::vector<double> x(numNodes), y(numNodes), z(numNodes);
-       std::vector<int> tag(numNodes);
-       
-       for(int nod = 0; nod < numNodes; nod++)
-         tag[nod] = eData->getNode(step, ent, ele, nod, x[nod], y[nod], z[nod]);
-       double n[3] = {0.,0.,0.};
-       normal3points(x[0], y[0], z[0], x[1], y[1], z[1], x[2], y[2], z[2], n); 
-       std::vector<double> valE(numNodes*numComp), valH(numNodes*numComp);
-       
-       for(int nod = 0; nod < numNodes; nod++){
-         if(tag[nod]) continue ; // already considered 
-         for(int comp = 0; comp < numComp; comp++){
-           eData->getValue(step, ent, ele, nod, comp, valE[numComp * nod + comp]);
-           hData->getValue(step, ent, ele, nod, comp, valH[numComp * nod + comp]);
-         }
-         double H[3] = { valH[numComp * nod + 0], valH[numComp * nod + 1], valH[numComp * nod + 2] } ;
-         double E[3] = { valE[numComp * nod + 0], valE[numComp * nod + 1], valE[numComp * nod + 2] } ;
-         double J[3], M[3] ;
-         prodve(n, H, J) ; // Js =  n x H ;  Surface electric current
-         prodve(E, n, M) ; // Ms = - n x E ; Surface magnetic current
-
-         for(int comp = 0; comp < numComp; comp++){
-           eData->setValue(step, ent, ele, nod, comp, M[comp]);
-           hData->setValue(step, ent, ele, nod, comp, J[comp]);
-           eData->tagNode(step, ent, ele, nod, 1);
-           hData->tagNode(step, ent, ele, nod, 1);
-         }
-       }
-     }
-   }
- }
+  std::vector<std::vector<double> > allPhi(_NbPhi+1,_NbThe+1) ;
+  std::vector<std::vector<double> > allThe(_NbPhi+1,_NbThe+1) ;
+  std::vector<std::vector<double> > farF(_NbPhi+1,_NbThe+1) ;
 
- eData->finalize(); 
- hData->finalize(); 
- 
-//****************************************
   for (int i = 0; i <= _NbPhi; i++){
     phi = i*dPhi ;
     for (int j = 0; j <= _NbThe; j++){
       theta = j * dTheta ;
       allPhi[i][j] = phi ;
       allThe[i][j] = theta ;
-      farF[i][j] = getFarField(eData, hData, _k0, _r_far, theta, phi) ; 
+
+      farF[i][j] = getFarField(allElems, js, ms, _k0, _r_far, theta, phi) ;
+
       ffmax = (ffmax < farF[i][j]) ? farF[i][j] : ffmax ;
     }
   }
-    
+
   if(_normalize){
       for (int i = 0; i <= _NbPhi; i++)
         for (int j = 0; j <= _NbThe; j++)
@@ -357,58 +269,49 @@ PView *GMSH_NearToFarFieldPlugin::execute(PView * v)
 
   // Constructing sphere for visualization
   // centered at center of bb and with radious relative to the bb size
-  double r_bb[3] = { eData->getBoundingBox().max().x()-eData->getBoundingBox().min().x(), 
-                     eData->getBoundingBox().max().y()-eData->getBoundingBox().min().y(), 
+  double r_bb[3] = { eData->getBoundingBox().max().x()-eData->getBoundingBox().min().x(),
+                     eData->getBoundingBox().max().y()-eData->getBoundingBox().min().y(),
                      eData->getBoundingBox().max().z()-eData->getBoundingBox().min().z() };
   double r_sph = norm3(r_bb) ;
   r_sph = (r_sph) ? r_sph/2 : 1./2. ; // radious of sphere for visu
 
   for (int i = 0; i < _NbPhi; i++){
     for (int j = 0; j < _NbThe; j++){
-      double P1[3] = { x0 + r_sph * farF[i  ][j  ] * sin(allThe[i  ][j  ]) * cos(allPhi[i  ][j  ]),
-                       y0 + r_sph * farF[i  ][j  ] * sin(allThe[i  ][j  ]) * sin(allPhi[i  ][j  ]),
-                       z0 + r_sph * farF[i  ][j  ] * cos(allThe[i  ][j  ]) } ;
-      double P2[3] = { x0 + r_sph * farF[i+1][j  ] * sin(allThe[i+1][j  ]) * cos(allPhi[i+1][j  ]),
-                       y0 + r_sph * farF[i+1][j  ] * sin(allThe[i+1][j  ]) * sin(allPhi[i+1][j  ]),
-                       z0 + r_sph * farF[i+1][j  ] * cos(allThe[i+1][j  ]) } ;
+      double P1[3] = { x0 + r_sph * farF[i ][j ] * sin(allThe[i ][j ]) * cos(allPhi[i ][j ]),
+                       y0 + r_sph * farF[i ][j ] * sin(allThe[i ][j ]) * sin(allPhi[i ][j ]),
+                       z0 + r_sph * farF[i ][j ] * cos(allThe[i ][j ]) } ;
+      double P2[3] = { x0 + r_sph * farF[i+1][j ] * sin(allThe[i+1][j ]) * cos(allPhi[i+1][j  ]),
+                       y0 + r_sph * farF[i+1][j ] * sin(allThe[i+1][j ]) * sin(allPhi[i+1][j  ]),
+                       z0 + r_sph * farF[i+1][j ] * cos(allThe[i+1][j ]) } ;
       double P3[3] = { x0 + r_sph * farF[i+1][j+1] * sin(allThe[i+1][j+1]) * cos(allPhi[i+1][j+1]),
                        y0 + r_sph * farF[i+1][j+1] * sin(allThe[i+1][j+1]) * sin(allPhi[i+1][j+1]),
                        z0 + r_sph * farF[i+1][j+1] * cos(allThe[i+1][j+1]) } ;
-      double P4[3] = { x0 + r_sph * farF[i  ][j+1] * sin(allThe[i  ][j+1]) * cos(allPhi[i  ][j+1]),
-                       y0 + r_sph * farF[i  ][j+1] * sin(allThe[i  ][j+1]) * sin(allPhi[i  ][j+1]),
-                       z0 + r_sph * farF[i  ][j+1] * cos(allThe[i  ][j+1]) } ;     
-
-      dataFar->SQ.push_back(P1[0]); dataFar->SQ.push_back(P2[0]);  dataFar->SQ.push_back(P3[0]); dataFar->SQ.push_back(P4[0]); 
-      dataFar->SQ.push_back(P1[1]); dataFar->SQ.push_back(P2[1]);  dataFar->SQ.push_back(P3[1]); dataFar->SQ.push_back(P4[1]); 
-      dataFar->SQ.push_back(P1[2]); dataFar->SQ.push_back(P2[2]);  dataFar->SQ.push_back(P3[2]); dataFar->SQ.push_back(P4[2]); 
-      (dataFar->NbSQ)++;      
+      double P4[3] = { x0 + r_sph * farF[i ][j+1] * sin(allThe[i ][j+1]) * cos(allPhi[i ][j+1]),
+                       y0 + r_sph * farF[i ][j+1] * sin(allThe[i ][j+1]) * sin(allPhi[i ][j+1]),
+                       z0 + r_sph * farF[i ][j+1] * cos(allThe[i ][j+1]) } ;
+
+      dataFar->SQ.push_back(P1[0]); dataFar->SQ.push_back(P2[0]);  dataFar->SQ.push_back(P3[0]); dataFar->SQ.push_back(P4[0]);
+      dataFar->SQ.push_back(P1[1]); dataFar->SQ.push_back(P2[1]);  dataFar->SQ.push_back(P3[1]); dataFar->SQ.push_back(P4[1]);
+      dataFar->SQ.push_back(P1[2]); dataFar->SQ.push_back(P2[2]);  dataFar->SQ.push_back(P3[2]); dataFar->SQ.push_back(P4[2]);
+      (dataFar->NbSQ)++;
       if(!_dB){
-        dataFar->SQ.push_back(farF[i  ][j  ]); 
-        dataFar->SQ.push_back(farF[i+1][j  ]); 
-        dataFar->SQ.push_back(farF[i+1][j+1]); 
-        dataFar->SQ.push_back(farF[i  ][j+1]);
+        dataFar->SQ.push_back(farF[i ][j  ]);
+        dataFar->SQ.push_back(farF[i+1][j  ]);
+        dataFar->SQ.push_back(farF[i+1][j+1]);
+        dataFar->SQ.push_back(farF[i ][j+1]);
       }
       else{
-        dataFar->SQ.push_back(log10(farF[i  ][j  ])); 
-        dataFar->SQ.push_back(log10(farF[i+1][j  ])); 
-        dataFar->SQ.push_back(log10(farF[i+1][j+1])); 
-        dataFar->SQ.push_back(log10(farF[i  ][j+1]));
+        dataFar->SQ.push_back(10*log10(farF[i ][j]));
+        dataFar->SQ.push_back(10*log10(farF[i+1][j ]));
+        dataFar->SQ.push_back(10*log10(farF[i+1][j+1]));
+        dataFar->SQ.push_back(10*log10(farF[i ][j+1]));
       }
     }
   }
 
-  for (int i = 0; i <= _NbPhi; i++){
-    delete[] allPhi[i] ;
-    delete[] allThe[i] ;
-    delete[] farF[i];
-  }
-  delete [] allPhi ;
-  delete [] allThe ;
-  delete [] farF ;
-
   dataFar->setName("_NearToFarField");
   dataFar->setFileName("_NearToFarField.pos");
   dataFar->finalize();
-  
+
   return vf;
 }
diff --git a/Plugin/NearToFarField.h b/Plugin/NearToFarField.h
index eed6474a3214753feeda853736c573dc6089f300..618982623dcd2f0032a449f877c936ad3a5e54c7 100644
--- a/Plugin/NearToFarField.h
+++ b/Plugin/NearToFarField.h
@@ -7,6 +7,11 @@
 #define _NEARTOFARFIELD_H_
 
 #include "Plugin.h"
+#include "shapeFunctions.h"
+#include "Numeric.h"
+#include "PViewOptions.h"
+#include "MElement.h"
+#include "GModel.h"
 
 extern "C"
 {
@@ -24,12 +29,12 @@ class GMSH_NearToFarFieldPlugin : public GMSH_PostPlugin
   }
   std::string getHelp() const;
   int getNbOptions() const;
-  StringXNumber* getOption(int iopt);  
+  StringXNumber* getOption(int iopt);
   PView *execute(PView *);
- 
-  static double getFarField(PViewData *eData, PViewData *hData, double k0, double r_far, double theta, double phi) ;
-  static void CartesianToSpherical(int numSteps, double theta, double phi, double **Fc, double **Fsp) ;
-    
-};
+
+  static double getFarField(std::vector<element*> allElems,
+                            std::vector<std::vector<double> >js, std::vector<std::vector<double> >ms,
+                            double k0, double r_far, double theta, double phi);
+ };
 
 #endif