diff --git a/Superconductors/helix.pro b/Superconductors/helix.pro
index fe60e7bd7527ff62c284206392639a8ba1e55b3c..426d1756794288a71935fb2ef776d3c83adee17f 100644
--- a/Superconductors/helix.pro
+++ b/Superconductors/helix.pro
@@ -5,6 +5,7 @@ Group {
   AirInf = Region[INF];
   Matrix = Region[MATRIX];
   BndMatrix = Region[BND_MATRIX];
+  BndInf = Region[BND_INF];
   Filaments = Region[{}];
   BndFilaments = Region[{}];
   For i In {1:NumLayers}
@@ -38,7 +39,7 @@ Function {
     sigmaMatrix = {6e7, Min 1e3, Max 6e7, Step 1e3, Visible ConductingMatrix,
       Name "Input/4Materials/Matrix conductivity [Sm⁻¹]"},
     Itot = {800/2, Min 100, Max 1000, Step 100,
-      Name "Input/3Source/Total current [A]"},
+      Name "Input/3Source/2Total current [A]"},
     Ec = {1e-4,
       Name "Input/4Materials/Critical electric field [Vm⁻¹]"},
     Jc = {5e8,
@@ -47,7 +48,7 @@ Function {
       ReadOnly Preset, Highlight "LightYellow",
       Name "Input/4Materials/Exponent (n) value"},
     Freq = {50, Min 1, Max 100, Step 1,
-      Name "Input/3Source/Frequency [Hz]"},
+      Name "Input/3Source/1Frequency [Hz]"},
     periods = {1., Min 0.1, Max 2.0, Step 0.05,
       Name "Input/Solver/0Periods to simulate"},
     time0 = 0, // initial time
@@ -65,7 +66,11 @@ Function {
     iter_max = {12,
       Name "Input/Solver/Maximum number of nonlinear iterations"},
     visu = {1, Choices{0, 1}, AutoCheck 0,
-      Name "Input/Solver/Visu", Label "Real-time visualization"}
+      Name "Input/Solver/Visu", Label "Real-time visualization"},
+    sourceField = {0, Choices{0, 1},
+      Name "Input/3Source/3Background transverse magnetic field?"},
+    hsValue = {1e5, Min 1e3, Max 1e7, Step 1e3, Visible sourceField,
+      Name "Input/3Source/4Source field amplitude [Am⁻¹]"}
   ];
 
   dt_max = adaptive ? dt_max : dt;
@@ -85,15 +90,25 @@ Function {
       Tensor[CompX[$1]^2, CompX[$1] * CompY[$1], CompX[$1] * CompZ[$1],
              CompY[$1] * CompX[$1], CompY[$1]^2, CompY[$1] * CompZ[$1],
              CompZ[$1] * CompX[$1], CompZ[$1] * CompY[$1], CompZ[$1]^2];
+
+  // test: source magnetic flux density
+  DtBs[] = hsValue * mu0 / Scaling * Sin_wt_p[]{2*Pi*Freq, 0.} * Vector[1,0,0];
 }
 
 Jacobian {
   { Name Vol ;
     Case {
-      { Region AirInf ; Jacobian VolCylShell{AirRadius*mm, InfRadius*mm} ; }
+      If(!sourceField)
+        { Region AirInf ; Jacobian VolCylShell{AirRadius*mm, InfRadius*mm} ; }
+      EndIf
       { Region All ; Jacobian Vol ; }
     }
   }
+  { Name Sur ;
+    Case {
+      { Region All ; Jacobian Sur ; }
+    }
+  }
 }
 
 Integration {
@@ -101,6 +116,7 @@ Integration {
     Case {
       { Type Gauss ;
 	Case {
+	  { GeoElement Line ; NumberOfPoints  4 ; }
 	  { GeoElement Triangle ; NumberOfPoints  4 ; }
           { GeoElement Quadrangle ; NumberOfPoints  4 ; }
 	  { GeoElement Tetrahedron ; NumberOfPoints  5 ; }
@@ -126,7 +142,7 @@ FunctionSpace {
   { Name HSpace; Type Form1;
     BasisFunction {
       { Name sn; NameOfCoef phin; Function BF_GradNode;
-        Support Omega; Entity NodesOf[OmegaCC]; }
+        Support Region[{Omega, BndInf}]; Entity NodesOf[OmegaCC]; }
       { Name se; NameOfCoef he; Function BF_Edge;
         Support OmegaC; Entity EdgesOf[All, Not BndOmegaC]; }
       { Name sc1; NameOfCoef I1; Function BF_GroupOfEdges;
@@ -174,8 +190,10 @@ Formulation {
       Galerkin { DtDof [ mu[] * Dof{h} , {h} ];
         In Omega; Integration Int; Jacobian Vol;  }
 
-      //Galerkin { [ mu[] * DtHs[] , {h} ];
-      //  In Omega; Integration Int; Jacobian Vol;  }
+      If(sourceField)
+        Galerkin { [ -DtBs[] * Normal[] , {dInv h} ];
+          In BndInf; Integration Int; Jacobian Sur;  }
+      EndIf
 
       Galerkin { [ rho[] * Dof{d h} , {d h} ];
         In Matrix; Integration Int; Jacobian Vol;  }