diff --git a/Conductors3D/w3d.geo b/Conductors3D/w3d.geo
index d103979ff78346c5d89fb5c552ecd49531e9ce0a..58a8812340665ad09d1b5a6d447561ff828022ef 100644
--- a/Conductors3D/w3d.geo
+++ b/Conductors3D/w3d.geo
@@ -6,19 +6,19 @@ Mesh.VolumeEdges = 0; // hide volume edges
 Geometry.ExactExtrusion = 0; // to allow rotation of extruded shapes
 Solver.AutoMesh = 2; // always remesh if necessary (don't reuse mesh on disk)
 
-lc1 = (Flag_Thin) ? MeshSizeThinWire*mm : MeshSizeThinWire*mm ;
+lc1 = (Flag_Thin) ? MeshSizeThinWire*mm : MeshSizeWire*mm ;
 lc2 = box/10; // mesh size at outer surface
 
 llWires[] = {};
 centerWires[] = {};
 surfWires[] = {};
 
-If( !Flag_Thin || (Flag_Thin && Flag_Corr==2)) 
+If( !Flag_Thin || (Flag_Thin && Flag_RegSleeve)) 
   // then circles are in the geometry
   // their radius is R=rw if Flag_Thin != 0
   // or R=rs in the "regular sleeve case" (Flag_Corr == 2)
   
-  R = ( !Flag_Thin ) ? rw : rs;
+  R = ( Flag_Thin ) ? rs : rw;
       
   For i In {1:NumWires}
     pC = newp; Point(pC) = {WX~{i}  , WY~{i}       , 0, lc1};
@@ -96,7 +96,7 @@ If( !Flag_Thin ) // Wires with finite radius
   e[] = Extrude {0,0,Lz} { Surface{ s1 } ; /*Layers{1}; Recombine;*/ } ;
   Physical Volume ("AIR", 1) = { e[1] };
 
-ElseIf( Flag_Corr == 2 ) // regular sleeve
+ElseIf( Flag_RegSleeve ) // regular sleeve
   
   Wires[]=  {};
   For i In {1:NumWires}
@@ -122,19 +122,18 @@ Else
   e[] = Extrude {0,0,Lz} { Surface{ s1 } ; /*Layers{1}; Recombine;*/ } ;
   s2 = e[0];
   v1 = e[1];
+  Physical Volume ("AIR", 1) = { v1 };
   For i In {1:NumWires}
     e[] = Extrude {0,0,Lz} { Point{ centerWires[i-1] } ; /*Layers{1};*/ } ;
     //Printf("e=", e[]);
     Physical Line  (Sprintf("LWIRE_%g",i), 50+i) = { e[1] };
     Physical Point (Sprintf("PANODE_%g",i), 60+i) = { centerWires[i-1] };
     Physical Point (Sprintf("PCATHODE_%g",i), 70+i) = { e[0] };
-    Physical Volume ("AIR", 1) = { v1 };
     // embed line conductors in mesh
     Point { centerWires[i-1] } In Surface { s1 };
     Point { e[0] } In Surface { s2 };
     Curve { e[1] } In Volume { v1 };
   EndFor
-
 EndIf 
 
 Physical Surface("INF", 2) = { CombinedBoundary{ Volume{ : }; } };
diff --git a/Conductors3D/w3d.pro b/Conductors3D/w3d.pro
index 14417682c54d50d999024d9eb37660626b2eb4e7..c4c5a47ad3d66002d0f32c666ccc54c4c8007aef 100644
--- a/Conductors3D/w3d.pro
+++ b/Conductors3D/w3d.pro
@@ -22,7 +22,7 @@ DefineConstant[
 
 DefineConstant[
   R_ = {"Dynamic", Name "GetDP/1ResolutionChoices", Visible 0},
-  C_ = {"-sol -pos", Name "GetDP/9ComputeCommand", Visible 0},
+  C_ = {"-solve -v2", Name "GetDP/9ComputeCommand", Visible 0},
   P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}
   ResDir = "Res/" 
 ];
@@ -63,43 +63,33 @@ Group{
 
 
   // Abstract regions
-  // Cette formulation peut contenir des conducteurs volumiques Vol_C
-  // et lin�iques lVol_C.
 
   If( !Flag_Thin )
     Vol_C   = Region[ VWIRES ];
     Vol_CC  = Region[ AIR ];
-    //lVol_C  = Region[ {} ];
   Else
     Vol_C   = Region[ {LWIRES} ];
     Vol_CC  = Region[ {VWIRES, AIR} ];
-    //lVol_C  = Region[ LWIRES ];
   EndIf
 
   Vol_nu = Region[ {Vol_C, Vol_CC} ];
   Sur_Dirichlet_a = Region[ INF ];
   Electrodes = Region[ {ANODES, CATHODES} ];   
   
-  Vol_Tree = Region[ { Vol_nu } ];
-  Sur_Tree = Region[ { Sur_Dirichlet_a /*, SKIN*/ } ];
-  // If( Flag_Corr!= 0 )
-  //   Lin_Tree = Region[ { LWIRES /*, LINTREE*/ } ];
-  // Else
-  Lin_Tree = Region[ {} ];
-  
-
   Dom_Hcurl_a = Region[ Vol_nu ];
   Dom_Hgrad_u = Region[ Vol_C ];
-  //Dom_Hregion_i = Region[ Vol_nu ];
-  
-  // additional Groups for the local correction method
-  //Dom_Hthin_a = ElementsOf[ Vol_nu, OnOneSideOf LWIRES ];
-  
-  Dom_Hthin_a = Region[ Vol_nu ];
-  
-  // For integration on the volume elements (excluding line elements) 
-  // of the sleeve
-  Vol_Sleeve = ElementsOf[ {Vol_CC}, OnOneSideOf LWIRES ]; 
+  Dom_Hregion_i = Region[ Vol_C ];
+
+  // additional Groups for the semi_analytic approach
+  Dom_Hthin_a = ElementsOf[ Vol_nu, OnOneSideOf LWIRES ];
+  //Vol_Tree = ElementsOf[ Vol_nu, Not Dom_Hthin_a ];
+  Vol_Tree = Region[ { Vol_nu } ];
+  Sur_Tree = Region[ { Sur_Dirichlet_a /*, SKIN*/ } ];
+  If( !Flag_SemiAnalytic )
+    Lin_Tree = Region[ {} ];
+  Else
+    Lin_Tree = Region[ { LWIRES } ];
+  EndIf
 }
   
 
@@ -109,13 +99,13 @@ Function{
   mur = 1.;
   sigma = 5.96e7; 
 
-  i[] = Complex[0,1];
-
   mu[] = mu0;
   nu[] = 1/mu[];
   sigma[] = sigma;
   skin_depth = Sqrt[ 2/(omega*sigma*mu0*mur) ];
-  
+
+  // Functions for the semi-analytic approach
+  i[] = Complex[0,1];
   tau[] = Complex[-1,1]/skin_depth*rw;
   J0[] = JnComplex[0,$1];
   J1[] = JnComplex[1,$1];
@@ -131,7 +121,7 @@ Function{
   EndFor
 
   // shape function of the current, wI[r] 
-  wI[] = tau[] * J0[tau[]*$1/rw] / J1[tau[]] /(2*A_c); 
+  wI[] = tau[] * J0[tau[]*$1/rw]/J1[tau[]]/(2*A_c); 
 
   //radial analytical solutions with a zero flux condition imposed at $1(=r)=rs
   Analytic_A[] = // per Amp
@@ -140,8 +130,8 @@ Function{
   AnalyticStatic_A[] = mu0/(2*Pi) * // per Amp
 	(($1>rw) ? Log[rs/$1] : Log[rs/rw] + mur*(1-($1/rw)^2)/2 );
   Analytic_B[] = // per Amp
-    (($1>rw) ?  mu0/(2*Pi*$1) :
-      mu0*mur/(2*Pi) * J1[tau[]*$1/rw] / J1[tau[]] );
+    (($1>=rw) ?  mu0/(2*Pi*$1) :
+      mu0*mur/(2*Pi*rw) * J1[tau[]*$1/rw] / J1[tau[]] );
  
   // Impedance of thin wire p.u. length
   R_DC = 1./(sigma*A_c);
@@ -155,11 +145,11 @@ Function{
   If ( NumWires == 2 )
     Exact_B[] = Analytic_B[R_1[]] + Analytic_B[R_2[]] ;
     Correction_B[] = ((R_1[]<rs) ? Analytic_B[R_1[]] :
-      ((R_2[]<rs) ? Analytic_B[R_2[]] : 0 ));
+                     ((R_2[]<rs) ? Analytic_B[R_2[]] : 0 ));
   EndIf
   If ( NumWires == 3 )
     Exact_B[] = Analytic_B[R_1[]] + Analytic_B[R_2[]] + Analytic_B[R_3[]] ;
-    Correction_A[] = ((R_1[]<rs) ? Analytic_B[R_1[]] : 
+    Correction_B[] = ((R_1[]<rs) ? Analytic_B[R_1[]] : 
                      ((R_2[]<rs) ? Analytic_B[R_2[]] : 
                      ((R_3[]<rs) ? Analytic_B[R_3[]] : 0 )));
   EndIf
@@ -180,14 +170,14 @@ Integration {
     Case {
       { Type Gauss ;
         Case {
-	  { GeoElement Point       ; NumberOfPoints  1 ; }
-	  { GeoElement Line        ; NumberOfPoints  3 ; }
-	  { GeoElement Triangle    ; NumberOfPoints  3 ; }
-	  { GeoElement Quadrangle  ; NumberOfPoints  4 ; }
-	  { GeoElement Prism       ; NumberOfPoints 21 ; }
-	  { GeoElement Tetrahedron ; NumberOfPoints  4 ; }
-	  { GeoElement Hexahedron  ; NumberOfPoints  6 ; }
-	  { GeoElement Pyramid     ; NumberOfPoints  8 ; }
+          { GeoElement Point       ; NumberOfPoints  1 ; }
+          { GeoElement Line        ; NumberOfPoints  3 ; }
+          { GeoElement Triangle    ; NumberOfPoints  3 ; }
+          { GeoElement Quadrangle  ; NumberOfPoints  4 ; }
+          { GeoElement Prism       ; NumberOfPoints 21 ; }
+          { GeoElement Tetrahedron ; NumberOfPoints  4 ; }
+          { GeoElement Hexahedron  ; NumberOfPoints  6 ; }
+          { GeoElement Pyramid     ; NumberOfPoints  8 ; }
         }
       }
     }
@@ -197,19 +187,19 @@ Integration {
 Constraint{
   { Name Impose_U ;
     Case {
-      For i In {1:NumWires}
-        If( Flag_Form == 1) // massive
+      For i In {1:NumWires}  
+        If( !Flag_Stranded ) // massive
           { Region ANODE~{i} ; Value  0 ; }
         EndIf
-        If( Flag_U )
-          
-          If( Flag_Form == 1) // massive
+        
+        If( Flag_U )  
+          If( !Flag_Stranded ) // massive
             { Region CATHODE~{i} ; Value -R_DC*Lz*WI~{i} ; }
           Else // stranded
             { Region VWIRE~{i} ; Value -R_DC*Lz*WI~{i} ; }
           EndIf
-          
         EndIf
+        
       EndFor
     }
   }
@@ -217,30 +207,48 @@ Constraint{
     Case {
       For i In {1:NumWires}
         If( !Flag_U )
-          
-          If( Flag_Form == 1) // massive
+          If( !Flag_Stranded ) // massive
             { Region ANODE~{i}   ; Value  WI~{i} ; }
             { Region CATHODE~{i} ; Value  WI~{i} ; }
           Else // stranded
             { Region VWIRE~{i} ; Value  WI~{i} ; }
             { Region LWIRE~{i} ; Value  WI~{i} ; }
-          EndIf
-            
+          EndIf            
         EndIf
       EndFor
     }
   }
 
-  { Name Hcurl_a_3D_ac; Type Assign;
+  { Name Hcurl_a_3D; Type Assign;
     Case {
       { Region Region[ Sur_Dirichlet_a ]; Value 0.; }
     }
   }
-
+    
   { Name GaugeCondition_a; Type Assign;
     Case {
       { Region Vol_Tree; SubRegion Region[ { Sur_Tree, Lin_Tree } ]; 
-        Value 0.; }
+        Value 0; }
+    }
+  }
+  
+  // Semi-analytic approach
+  { Name Impose_corr ;
+    Case {
+      For i In {1:NumWires}
+        { Region LWIRE~{i} ; Value  1. ; }
+      EndFor
+    }
+  }
+
+  { Name Hcurl_acorr_3D; Type Assign;
+    Case {
+      { Region Region[ Sur_Dirichlet_a ]; Value 0.; }
+      If( Flag_Thin )
+        For i In {1:NumWires}
+          { Region Region[ LWIRE~{i} ]; Value 1 ; }
+        EndFor
+      EndIf
     }
   }
 }
@@ -256,7 +264,7 @@ FunctionSpace {
     }
     Constraint {
       { NameOfCoef ac;
-		EntityType Auto; NameOfConstraint Hcurl_a_3D_ac; }
+		EntityType Auto; NameOfConstraint Hcurl_a_3D; }
 	  { NameOfCoef ac; EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
 		NameOfConstraint GaugeCondition_a ; }
     }
@@ -280,10 +288,11 @@ FunctionSpace {
     }
   }
 
+  // Electric current per regionl for stranded conductors
   { Name Hregion_i_3D; Type Vector;
     BasisFunction {
       { Name sr; NameOfCoef ir; Function BF_RegionZ;
-        Support Dom_Hthin_a; Entity Vol_C; }
+        Support Dom_Hregion_i; Entity Vol_C; }
     }
     GlobalQuantity {
       { Name Is; Type AliasOf       ; NameOfCoef ir; }
@@ -295,52 +304,93 @@ FunctionSpace {
     }
   }
 
-  // For the local correction method
-  { Name Hcurl_acorr_3D; Type Form1;
+  // Function spaces for the semi-analytical approach
+  
+  // { Name Hcurl_acorr_3D; Type Form1;
+  //   BasisFunction {
+  //     { Name sw; NameOfCoef aw; Function BF_Edge;
+  //       Support Dom_Hthin_a; Entity EdgesOf[ Vol_nu, ConnectedTo LWIRES ]; }
+  //   }
+  //   Constraint {
+  //     { NameOfCoef aw;
+  //   	EntityType Auto; NameOfConstraint Hcurl_a_3D; }
+  //     { NameOfCoef aw; EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
+  //     	NameOfConstraint GaugeCondition_a ; }
+  //   }
+  // }
+  
+  { Name Hcurl_athin_3D; Type Form1;
     BasisFunction {
-      { Name sw; NameOfCoef aw; Function BF_Edge;
-        Support Vol_nu; Entity EdgesOf[ Vol_nu, ConnectedTo LWIRES ]; }
+      { Name si; NameOfCoef ai; Function BF_GroupOfEdges;
+        Support Dom_Hcurl_a; Entity GroupsOfEdgesOf[ LWIRES ]; }
+      { Name sj; NameOfCoef aj; Function BF_Edge;
+        Support Dom_Hcurl_a; Entity EdgesOf[ Vol_nu, Not LWIRES ]; }
+    }
+    GlobalQuantity {
+      { Name F; Type AliasOf       ; NameOfCoef ai; }
+      { Name I; Type AssociatedWith; NameOfCoef ai; }
     }
     Constraint {
-	  // { NameOfCoef aw; EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
-	  // 	NameOfConstraint GaugeCondition_a ; }
+      { NameOfCoef I; EntityType Region; NameOfConstraint Impose_corr; }
+      { NameOfCoef aj; EntityType Auto; NameOfConstraint Hcurl_a_3D; }
+      { NameOfCoef aj; EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
+       	NameOfConstraint GaugeCondition_a ; }
+    }
+  }
+  { Name Hcurl_asleeve_3D; Type Form1;
+    BasisFunction {
+      { Name si; NameOfCoef ai; Function BF_GroupOfEdges;
+        Support Dom_Hthin_a; Entity GroupsOfEdgesOf[ LWIRES ]; }
+      { Name sj; NameOfCoef aj; Function BF_Edge;
+        Support Dom_Hthin_a; Entity EdgesOf[ Vol_nu, ConnectedTo LWIRES ]; }
+    }
+    GlobalQuantity {
+      { Name F; Type AliasOf       ; NameOfCoef ai; }
+      { Name I; Type AssociatedWith; NameOfCoef ai; }
+    }
+    Constraint {
+      { NameOfCoef I; EntityType Region; NameOfConstraint Impose_corr; }
+      { NameOfCoef aj; EntityType Auto; NameOfConstraint Hcurl_a_3D; }
     }
   }
 }
 
 Formulation {
   { Name MagnetoDynamics; Type FemEquation;
+    
+    If( !Flag_SemiAnalytic )
+      // Conventional massive and stranded conductor formulation
+      // If Flag_Thin is true, naive thine wire formulations
+      
+      If( !Flag_Stranded ) // Massive conductor
+        Quantity {
+          { Name a;  Type Local;  NameOfSpace Hcurl_a_3D; }
 
-    If( Flag_Form == 1) // Massive conductor
-      Quantity {
-        { Name a;  Type Local;  NameOfSpace Hcurl_a_3D; }
-
-        { Name v;  Type Local;  NameOfSpace Hgrad_u_3D; }
-        { Name U;  Type Global; NameOfSpace Hgrad_u_3D [U]; }
-        { Name I;  Type Global; NameOfSpace Hgrad_u_3D [I]; }
-      }
+          { Name v;  Type Local;  NameOfSpace Hgrad_u_3D; }
+          { Name U;  Type Global; NameOfSpace Hgrad_u_3D [U]; }
+          { Name I;  Type Global; NameOfSpace Hgrad_u_3D [I]; }
+        }
 
-      Equation {   
-        Integral { [ nu[] * Dof{d a} , {d a} ];
-          In Vol_nu; Jacobian Vol; Integration I1; }
+        Equation {   
+          Integral { [ nu[] * Dof{d a} , {d a} ];
+            In Vol_nu; Jacobian Vol; Integration I1; }
 
-        Integral { DtDof[ sigma[] * Dof{a} , {a} ];
-          In Vol_C; Jacobian Vol; Integration I1; }
-        Integral { [ sigma[] * Dof{d v} , {a} ];
-          In Vol_C; Jacobian Vol; Integration I1; }
-        Integral { DtDof[ sigma[] * Dof{a} , {d v} ];
-          In Vol_C; Jacobian Vol; Integration I1; }      
-        Integral { [ sigma[] * Dof{d v} , {d v} ];
-          In Vol_C; Jacobian Vol; Integration I1; } 
+          Integral { DtDof[ sigma[] * Dof{a} , {a} ];
+            In Vol_C; Jacobian Vol; Integration I1; }
+          Integral { [ sigma[] * Dof{d v} , {a} ];
+            In Vol_C; Jacobian Vol; Integration I1; }
+          Integral { DtDof[ sigma[] * Dof{a} , {d v} ];
+            In Vol_C; Jacobian Vol; Integration I1; }      
+          Integral { [ sigma[] * Dof{d v} , {d v} ];
+            In Vol_C; Jacobian Vol; Integration I1; } 
         
-        GlobalTerm { [Dof{I} , {U} ]; In CATHODES; }  
-      }
+          GlobalTerm { [Dof{I} , {U} ]; In CATHODES; }  
+        }
       
-    Else // stranded conductor   
-      If( Flag_Corr==0 ) // naive thin wire model
+        Else // stranded conductor
         Quantity {
           { Name a;  Type Local;  NameOfSpace Hcurl_a_3D; }
-
+            
           { Name  i; Type Local;  NameOfSpace Hregion_i_3D; }
           { Name Is; Type Global; NameOfSpace Hregion_i_3D [Is]; }
           { Name Us; Type Global; NameOfSpace Hregion_i_3D [Us]; }
@@ -351,7 +401,7 @@ Formulation {
             In Vol_nu; Jacobian Vol; Integration I1; }
           Integral { [ -Dof{i}/A_c , {a} ];
             In Vol_C; Jacobian Vol; Integration I1; }
-        
+            
           // Integral { [ -J[] , {a} ];
           //   In Vol_C; Jacobian Vol; Integration I1; }
           
@@ -361,44 +411,53 @@ Formulation {
             In Vol_C; Jacobian Vol; Integration I1;}
           Integral { [ Dof{i}/(sigma[]*A_c) , {i} ];
             In Vol_C; Jacobian Vol; Integration I1;}
-            GlobalTerm { [ Dof{Us}*A_c , {Is} ]; In Vol_C; }
-        }
-      Else // use local correction method
-        Quantity {
-          { Name a;  Type Local;  NameOfSpace Hcurl_a_3D; }
-          { Name as; Type Local;  NameOfSpace Hcurl_acorr_3D; }
-
-          { Name  i; Type Local;  NameOfSpace Hregion_i_3D; }
-          { Name Is; Type Global; NameOfSpace Hregion_i_3D [Is]; }
-          { Name Us; Type Global; NameOfSpace Hregion_i_3D [Us]; }
+          GlobalTerm { [ Dof{Us}*A_c , {Is} ]; In Vol_C; }
         }
+      EndIf
+        
+    EndIf
       
-        Equation {
-          /*
-          Integral { [ nu[] * Dof{d a} , {d a} ];
-            In Vol_nu; Jacobian Vol; Integration I1; }
-          Integral { [ -Dof{i}/A_c , {a} ];
-            In Vol_C; Jacobian Vol; Integration I1; }
-            */
-          Integral { [ nu[] * Dof{d as} , {d as} ];
-            In Vol_nu; Jacobian Vol; Integration I1; }
-          Integral { [ -Dof{i}/A_c , {as} ];
-            In Vol_C; Jacobian Vol; Integration I1; }
+    If( Flag_SemiAnalytic )
+        // Semi-analytic correction methods 
+        
+      Quantity {
+        { Name a;  Type Local; NameOfSpace Hcurl_athin_3D; }
+        { Name F; Type Global; NameOfSpace Hcurl_athin_3D [F]; }
+        { Name I; Type Global; NameOfSpace Hcurl_athin_3D [I]; }
           
-          /*    
-          GlobalTerm { [ Analytic_R[] * Dof{Is} , {Is} ]; 
-            In Vol_C;}
-          Integral { DtDof [ Dof{a}/A_c , {i} ];
-            In Vol_C; Jacobian Vol; Integration I1; }
-          Integral { DtDof [ -Dof{as}/A_c , {i} ];
-            In Vol_C; Jacobian Vol; Integration I1;}
-          GlobalTerm { DtDof [ Analytic_L[] * Dof{Is} , {Is} ]; 
-            In Vol_C;}
-          GlobalTerm { [ Dof{Us} , {Is} ]; 
-            In Vol_C; }
-            */
-        }
-      EndIf
+        { Name as; Type Local;  NameOfSpace Hcurl_asleeve_3D; }
+        { Name Fs; Type Global; NameOfSpace Hcurl_asleeve_3D [F]; }
+        { Name Is; Type Global; NameOfSpace Hcurl_asleeve_3D [I]; }
+      }
+      
+      Equation {
+          
+        Integral { [ nu[] * Dof{d a} , {d a} ];
+          In Vol_nu; Jacobian Vol; Integration I1; }
+        GlobalTerm { [ -Dof{I}*1e1 , {F} ]; 
+          In Vol_nu; }
+          
+        Integral { [ nu[] * Dof{d as} , {d as} ];
+          In Vol_nu; Jacobian Vol; Integration I1; }
+        GlobalTerm { [ -Dof{Is}*1e1 , {Fs} ]; 
+          In Vol_nu; }
+            
+        // Integral { [ -Dof{i}/A_c , {as} ];
+        //   In Vol_C; Jacobian Vol; Integration I1; }
+          
+        /*    
+              GlobalTerm { [ Analytic_R[] * Dof{Is} , {Is} ]; 
+              In Vol_C;}
+              Integral { DtDof [ Dof{a}/A_c , {i} ];
+              In Vol_C; Jacobian Vol; Integration I1; }
+              Integral { DtDof [ -Dof{as}/A_c , {i} ];
+              In Vol_C; Jacobian Vol; Integration I1;}
+              GlobalTerm { DtDof [ Analytic_L[] * Dof{Is} , {Is} ]; 
+              In Vol_C;}
+              GlobalTerm { [ Dof{Us} , {Is} ]; 
+              In Vol_C; }
+        */
+      }
     EndIf
   }
 }
@@ -443,41 +502,65 @@ PostProcessing {
       { Name b;
         Value { Local { [ {d a}]; In  Dom_Hcurl_a; Jacobian Vol; }}}
       { Name by;
-        Value { Local { [ Re[Cart2Pol[CompY [ {d a}]]]]; 
+        Value { Local { [ Cart2Pol[ Norm[ {d a}]] ]; 
             In  Dom_Hcurl_a; Jacobian Vol; }}}
-
-      If( Flag_Form == 1 ) // massive
-        { Name J;
-          Value { Local { [ -sigma[] * ( Dt[{a}] + {d v} ) ]; 
-			In Vol_C; Jacobian Vol; }}}
-        { Name U;
-          Value { Term { [ {U} ]; In Electrodes; }}}
+      
+      If( Flag_SemiAnalytic )
+        { Name bs;
+          Value { Local { [ {d as} ]; In Vol_nu; Jacobian Vol; }}}
+        { Name bsy;
+          Value { Local { [ Norm[ {d as} ] ]; In Vol_nu; Jacobian Vol; }}}
+        { Name bbsy;
+          Value { Local { [ Norm[ {d a} - {d as} ] ]; In Vol_nu; Jacobian Vol; }}}
+        { Name bcorry;
+          Value { Local { [ Norm[ {d a} - {d as} ] + Correction_B[] ]; 
+              In Vol_nu; Jacobian Vol; }}}
+      EndIf
+      
+      
+      
+      If( !Flag_SemiAnalytic )
+        If( !Flag_Stranded ) // massive
+          { Name J;
+            Value { Local { [ -sigma[] * ( Dt[{a}] + {d v} ) ]; 
+                In Vol_C; Jacobian Vol; }}}
+          { Name U;
+            Value { Term { [ {U} ]; In Electrodes; }}}
+          { Name I;
+            Value { Term { [ {I} ]; In Electrodes; }}}
+          { Name R;
+            Value { Term { [ Re[ -{U}/{I}/Lz ] ]; In Electrodes; }}}
+          { Name L; // actually total flux/I
+            Value { Term { [ Im[ -{U}/{I}/omega/Lz ] ]; In Electrodes; }}}
+        Else // stranded
+          { Name J;
+            Value { Local { [ {Is}/A_c ]; 
+                In Vol_C; Jacobian Vol; }}}
+          { Name U;
+            Value { Term { [ {Us} ]; In Vol_C; }}}
+          { Name I;
+            Value { Term { [ {Is} ]; In Vol_C; }}}
+          { Name R;
+            Value { Term { [ Re[ -{Us}/{Is}/Lz ] ]; In Vol_C; }}}
+          { Name L; // actually total flux/I
+            Value { Term { [ Im[ -{Us}/{Is}/omega/Lz ] ]; In Vol_C; }}}
+        EndIf
+      EndIf
+      
+      If( Flag_SemiAnalytic )        
+        { Name F;
+          Value { Term { [ {F} ]; In Vol_C; }}}
         { Name I;
-          Value { Term { [ {I} ]; In Electrodes; }}}
-        { Name R;
-          Value { Term { [ Re[ -{U}/{I}/Lz ] ]; In Electrodes; }}}
-        { Name L; // actually total flux/I
-          Value { Term { [ Im[ -{U}/{I}/omega/Lz ] ]; In Electrodes; }}}
-      Else // stranded
+          Value { Term { [ {I} ]; In Vol_C; }}}
         { Name J;
-          Value { Local { [ {i}/A_c ]; 
+          Value { Local { [ {I}/A_c ]; 
               In Vol_C; Jacobian Vol; }}}
         { Name U;
-          Value { Term { [ {Us} ]; In Vol_C; }}}
-        { Name I;
-          Value { Term { [ {Is} ]; In Vol_C; }}}
+          Value { Term { [ i[]*omega*{F}*10 + Analytic_R[]*{I} ]; In Vol_C; }}}
         { Name R;
-          Value { Term { [ Re[ -{Us}/{Is}/Lz ] ]; In Vol_C; }}}
+          Value { Term { [ Analytic_R[] + {I}*0 ]; In Vol_C; }}}
         { Name L; // actually total flux/I
-          Value { Term { [ Im[ -{Us}/{Is}/omega/Lz ] ]; In Vol_C; }}}
-      EndIf
-      
-      If( (Flag_Form == 2) && Flag_Corr )
-        { Name bs;
-          Value { Local { [ {d as}]; In Dom_Hthin_a; Jacobian Vol; }}}
-        { Name bcorry;
-          Value { Local { [ Re[Cart2Pol[CompY [{d a}-{d as}]]] + Correction_B[] ]; 
-              In  Dom_Hthin_a; Jacobian Vol; }}}
+          Value { Term { [ {F}/{I}/Lz ]; In Vol_C; }}}
       EndIf
         
     }
@@ -485,28 +568,48 @@ PostProcessing {
 }
 
 PostOperation map UsingPost MagnetoDynamics {
-  Print[ J, OnElementsOf Vol_C, File "j.pos"];
+     
   Print[ b, OnElementsOf Vol_nu, File "b.pos"];
-  If( (Flag_Form == 2) && Flag_Corr )
+    
+  If( !Flag_SemiAnalytic )
+    Print[ J, OnElementsOf Vol_C, File "j.pos"];
+  Else
     Print[ bs, OnElementsOf Vol_nu, File "bs.pos"];
     PrintGroup[ EdgesOf[ Vol_nu, ConnectedTo LWIRES ],
       In Vol_nu, File "group.pos"];
   EndIf
+  PrintGroup[ EdgesOfTreeIn[ { Vol_Tree }, StartingOn { Sur_Tree, Lin_Tree } ],
+    In Vol_Tree, File "Tree.pos"];
+  //PrintGroup[ _CO_Entity_44, In Vol_nu, File "Tree.pos"];
+
 }
 
 PostOperation integaz UsingPost MagnetoDynamics {
   For i In {1:NumWires}
-    If( Flag_Form == 1 ) // massive
-        
-      Print[ U, OnRegion CATHODES, File > "U.dat", Format Table,
-        SendToServer Sprintf["Results/1Voltage/Wire %g",i]];
-      Print[ I, OnRegion CATHODES, File > "I.dat", Format Table,
-        SendToServer Sprintf["Results/2Current/Wire %g",i]];
-      Print[ R, OnRegion CATHODES, File > "Z.dat", Format Table,
-        SendToServer Sprintf["Results/3Resistance p.u.l./Wire %g",i]];
-      Print[ L, OnRegion CATHODES, File > "Z.dat", Format Table,
-        SendToServer Sprintf["Results/4Inductance p.u.l./Wire %g",i]];
-    Else
+    If( !Flag_SemiAnalytic )
+      If( !Flag_Stranded ) // massive
+        Print[ U, OnRegion CATHODES, File > "U.dat", Format Table,
+          SendToServer Sprintf["Results/1Voltage/Wire %g",i]];
+        Print[ I, OnRegion CATHODES, File > "I.dat", Format Table,
+          SendToServer Sprintf["Results/2Current/Wire %g",i]];
+        Print[ R, OnRegion CATHODES, File > "Z.dat", Format Table,
+          SendToServer Sprintf["Results/3Resistance p.u.l./Wire %g",i]];
+        Print[ L, OnRegion CATHODES, File > "Z.dat", Format Table,
+          SendToServer Sprintf["Results/4Inductance p.u.l./Wire %g",i]];
+      Else
+        Print[ U, OnRegion Vol_C, File > "U.dat", Format Table,
+          SendToServer Sprintf["Results/1Voltage/Wire %g",i]];
+        Print[ I, OnRegion Vol_C, File > "I.dat", Format Table,
+          SendToServer Sprintf["Results/2Current/Wire %g",i]];
+        Print[ R, OnRegion Vol_C, File > "Z.dat", Format Table,
+          SendToServer Sprintf["Results/3Resistance p.u.l./Wire %g",i]];
+        Print[ L, OnRegion Vol_C, File > "Z.dat", Format Table,
+          SendToServer Sprintf["Results/4Inductance p.u.l./Wire %g",i]];
+      EndIf
+    EndIf
+    If( Flag_SemiAnalytic )
+      Print[ F, OnRegion Vol_C, File > "Flux.dat", Format Table,
+        SendToServer Sprintf["Results/0Flux/Wire %g",i]];
       Print[ U, OnRegion Vol_C, File > "U.dat", Format Table,
         SendToServer Sprintf["Results/1Voltage/Wire %g",i]];
       Print[ I, OnRegion Vol_C, File > "I.dat", Format Table,
@@ -519,27 +622,66 @@ PostOperation integaz UsingPost MagnetoDynamics {
   EndFor
 }
 
-NbPoints = 400;
-Xmax = 0.15*box;
+NbPoints = 1000;
+Xcut = 0.1*box;
+Zcut = Lz/2;
 
 PostOperation cut UsingPost MagnetoDynamics {
-    Print [ by, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
-			Format Gmsh, File "by.pos" ];
+  Print [ by, OnLine { {-Xcut,0,0} {Xcut,0,Zcut} }{NbPoints},
+          Format Gmsh, File "by.pos" ];
+  Echo[ StrCat["l=PostProcessing.NbViews-1;",
+      "View[l].Name = 'cut by';",
+      "View[l].Axes = 3;",
+      "View[l].LineWidth = 3;",
+      "View[l].Type = 2;"],
+    File "tmp.geo", LastTimeStepOnly];
+      
+  If( Flag_SemiAnalytic )
+    Print [ bsy, OnLine { {-Xcut,0,0} {Xcut,0,Zcut} }{NbPoints},
+      Format Gmsh, File "bcorry.pos" ];
     Echo[ StrCat["l=PostProcessing.NbViews-1;",
-				 "View[l].Name = 'cut by';",
-				 "View[l].Axes = 3;",
-				 "View[l].LineWidth = 3;",
-				 "View[l].Type = 2;"],
-	File "tmp.geo", LastTimeStepOnly];
-  If( (Flag_Form == 2) && Flag_Corr )
-    Print [ bcorry, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
-      Format Gmsh, File "by.pos" ];
+        "View[l].Name = 'cut bsy ';",
+        "View[l].Axes = 3;",
+        "View[l].LineWidth = 3;",
+        "View[l].Type = 2;"],
+      File "tmp.geo", LastTimeStepOnly];
+
+    Print [ bbsy, OnLine { {-Xcut,0,0} {Xcut,0,Zcut} }{NbPoints},
+      Format Gmsh, File "bcorry.pos" ];
     Echo[ StrCat["l=PostProcessing.NbViews-1;",
-        "View[l].Name = 'cut by';",
+        "View[l].Name = 'cut by-bsy ';",
         "View[l].Axes = 3;",
         "View[l].LineWidth = 3;",
         "View[l].Type = 2;"],
       File "tmp.geo", LastTimeStepOnly];
+
+    Print [ bcorry, OnLine { {-Xcut,0,0} {Xcut,0,Zcut} }{NbPoints},
+      Format Gmsh, File "bcorry.pos" ];
+    Echo[ StrCat["l=PostProcessing.NbViews-1;",
+        "View[l].Name = 'cut by corrected';",
+        "View[l].Axes = 3;",
+        "View[l].LineWidth = 3;",
+        "View[l].Type = 2;"],
+      File "tmp.geo", LastTimeStepOnly];
+  EndIf
+  
+  // cuts in txt format for Gnuplot
+  // Print [ exact, OnLine { {0,0,0} {Xcut,0,0} } {NbPoints},
+  //         Format SimpleTable, File "Cut_analytic.txt" ];
+
+  
+  If( !Flag_SemiAnalytic )
+    Print [ by, OnLine { {0,0,0} {Xcut,0,0} } {NbPoints},
+      Format SimpleTable, File "Cut_byref.txt" ];
+  Else
+    Print [ by, OnLine { {0,0,0} {Xcut,0,0} } {NbPoints},
+      Format SimpleTable, File "Cut_by.txt" ];
+    Print [ bsy, OnLine { {0,0,0} {Xcut,0,0} } {NbPoints},
+      Format SimpleTable, File "Cut_bsy.txt" ];
+    Print [ bbsy, OnLine { {0,0,0} {Xcut,0,0} } {NbPoints},
+      Format SimpleTable, File "Cut_bbsy.txt" ];
+    Print [ bcorry, OnLine { {0,0,0} {Xcut,0,0} } {NbPoints},
+      Format SimpleTable, File "Cut_bcorry.txt" ];
   EndIf
 }
 
diff --git a/Conductors3D/w3d_common.pro b/Conductors3D/w3d_common.pro
index e8dc467f0dfc6edd0cfbabd4403ae550b31be86b..3f34eed820047969ea14c92178c5a3f4c2ca3992 100644
--- a/Conductors3D/w3d_common.pro
+++ b/Conductors3D/w3d_common.pro
@@ -4,32 +4,43 @@ deg = 180/Pi;
 
 /*
  Expected inductance of a rectilinear 1mm radius wire: 4 nH/cm 
+
+ git/documentation/models/Inductor/magstadyn_av_js0_3d.pro 
  */
 
 DefineConstant[
-  NumWires = {1, Name "Parameters/0Number of wires", 
+  NumWires = {1, Name "Parameters/00Number of wires", 
 			  Choices {1="1", 2="2", 3="3"}, Visible 1}
-  Flag_Form = {2, Name "Parameters/1Conductor type", Visible 1, 
-               Choices {1="Massive", 2="Stranded"}}
-  Flag_Thin = {1, Name "Parameters/2Use thin wires", Choices {0,1}, Visible 1}
-  Flag_Corr = {1, Name "Parameters/3Correction method", 
-    Visible (Flag_Form==2) && Flag_Thin, 
-    Choices {0="None", 1="Raw sleeve", 2="regular sleeve"}}
-  Flag_U = {0, Name "Parameters/4Impose U", Choices {0,1}, Visible 1}
+  
+  Flag_Thin = {1, Name "Parameters/01Thin wires", 
+    Choices {0,1}, Visible 1}
+  Flag_RegSleeve = {0, Name "Parameters/02Regular sleeves", 
+    Choices {0,1}, Visible Flag_Thin}
+  
+  Flag_SemiAnalytic = {1, Name "Parameters/03Semi-analytic approach", 
+    Choices {0,1}, Visible Flag_Thin}
+  
+  Flag_Stranded = {0, Name "Parameters/04Stranded conductors", 
+                   Choices {0,1}, Visible !Flag_SemiAnalytic || !Flag_Thin}
+  Flag_Dual = {0, Name "Parameters/05Dual approach", 
+    Choices {0,1}, Visible Flag_SemiAnalytic && Flag_Thin}
+  Flag_U = {0, Name "Parameters/06Impose U", Choices {0,1}, Visible 1}
+  
+  WireRadius   = {1, Name "Parameters/10Wire radius [mm]"}
+  MeshSizeWire = {0.25, Name "Parameters/11Mesh size in wire [mm]", 
+    Visible !Flag_Thin}
+  MeshSizeThinWire = {2, Name "Parameters/12Mesh size on thin wire [mm]", 
+    Visible Flag_Thin}
+  //SleeveRadius = {2, Name "Parameters/13Sleeve radius [mm]", Visible Flag_RegSleeve}
+  box = {100*mm, Name "Parameters/5box half-width [m]"}
+  
   Flag_Maps = {1, Name "Input/2Display maps", Choices {0,1}, Visible 1}
   
   LogFreq = {4,  Min 1, Max 8, Step 0.25, Name "Input/4Log of frequency"}
   Freq = 10^LogFreq
-  WireRadius   = {1, Name "Parameters/5Wire radius [mm]"}
-  MeshSizeWire = {1, Name "Parameters/6Mesh size in wire [mm]", 
-    Visible !Flag_Thin}
-  MeshSizeThinWire = {1, Name "Parameters/7Mesh size on thin wire [mm]", 
-    Visible Flag_Thin}
-  SleeveRadius = {2, Name "Parameters/8Sleeve radius [mm]"}
-  box = {100*mm, Name "Parameters/9box half-width [m]"}
-
+  
   rw = WireRadius*mm 
-  rs = SleeveRadius*mm
+  rs = MeshSizeThinWire*mm
   A_c = Pi*rw^2 // wire cross section
   rout = box
   Lz = 10*mm
@@ -41,6 +52,10 @@ DefineConstant[
    */
 ];
 
+If( Flag_SemiAnalytic ) 
+  SetNumber("Parameters/01Thin wires", 1);
+  Flag_Thin=1;
+EndIf
 
 
 Xlocs() = { 0, 5*mm, -5*mm};
@@ -50,7 +65,7 @@ For i In {1:NumWires}
 	WX~{i} = { Xlocs(i-1), // place wires symmetrically around x=0
 			   Name Sprintf("Parameters/Wire %g/0X position [m]", i), Closed }
     WY~{i} = { Ylocs(i-1), Name Sprintf("Parameters/Wire %g/0Y position [m]", i) }
-    WI~{i} = { 1.23456, Name Sprintf("Parameters/Wire %g/2Current [A]", i) }
+    WI~{i} = { 1, Name Sprintf("Parameters/Wire %g/2Current [A]", i) }
     WP~{i} = { 0*NumWires*(i-1),
       Name Sprintf("Parameters/Wire %g/3Phase [deg]", i) }
     ];