diff --git a/include/maxwell/maxwell_common.h b/include/maxwell/maxwell_common.h
index d4cfc6be3889b81c37d22f5f671dcbb67ae94c0c..463481727547b89a96f73ba32061c372878f857d 100644
--- a/include/maxwell/maxwell_common.h
+++ b/include/maxwell/maxwell_common.h
@@ -62,6 +62,9 @@ eval_boundary_sources(const model& mod, const parameter_loader& mpl,
                 case boundary_condition::H_FIELD:
                     break;
 
+                case boundary_condition::EH_FIELD:
+                    break;
+
                 case boundary_condition::IMPEDANCE:
                     break;
 
@@ -164,6 +167,24 @@ eval_boundary_sources_new(const model& mod, const parameter_loader& mpl,
             case boundary_condition::SURFACE_CURRENT:
                     break;
 
+            case boundary_condition::EH_FIELD: {
+                for (auto& [local, global] : sd.face_numbers)
+                {
+                    auto fEH = [&](const point_3d& pt) -> std::pair<vec3d, vec3d> {
+                        return mpl.eval_boundary_source_eh(mtag, pt, state.curr_time);
+                    };
+
+                    matxd prj = project_EH_to_face(e, local, fEH);
+                    auto num_bf = prj.rows();
+                    srcfield.Ex.segment(global*num_bf, num_bf) = prj.col(0);
+                    srcfield.Ey.segment(global*num_bf, num_bf) = prj.col(1);
+                    srcfield.Ez.segment(global*num_bf, num_bf) = prj.col(2);
+                    srcfield.Hx.segment(global*num_bf, num_bf) = prj.col(3);
+                    srcfield.Hy.segment(global*num_bf, num_bf) = prj.col(4);
+                    srcfield.Hz.segment(global*num_bf, num_bf) = prj.col(5);
+                }
+            } break; /* case boundary_condition::EH_FIELD */
+
             case boundary_condition::PLANE_WAVE_E: {
                 for (auto& [local, global] : sd.face_numbers)
                 {
diff --git a/include/maxwell/maxwell_interface.h b/include/maxwell/maxwell_interface.h
index af0379ec1113d414241eae07122ac44c3c8ab7f2..e9bcfde5fa535f6a89fe738fffb6cf471d0c518a 100644
--- a/include/maxwell/maxwell_interface.h
+++ b/include/maxwell/maxwell_interface.h
@@ -43,6 +43,7 @@ enum class boundary_condition
     PLANE_WAVE_H,
     E_FIELD,
     H_FIELD,
+    EH_FIELD,
     SURFACE_CURRENT
 };
 
@@ -80,6 +81,7 @@ public:
     bool    volume_has_source(int) const;
     vec3d   eval_volume_source(int, const point_3d&, double) const;
     vec3d   eval_boundary_source(int, const point_3d&, double) const;
+    std::pair<vec3d,vec3d> eval_boundary_source_eh(int, const point_3d&, double) const;
     vec3d   eval_interface_source(int, const point_3d&, double) const;
 
     bool    has_analytical_solution(void) const;
diff --git a/src/maxwell/maxwell_gpu.cpp b/src/maxwell/maxwell_gpu.cpp
index e13aff21bc11f7e49aee6360227941af541aba79..25b7a288ddc92d05855d630838243a9ad31d424c 100644
--- a/src/maxwell/maxwell_gpu.cpp
+++ b/src/maxwell/maxwell_gpu.cpp
@@ -105,6 +105,7 @@ build_source_compression_tables(const model& mod, solver_state_gpu& state, const
                 switch ( mpl.boundary_type(bd.gmsh_tag()) )
                 {
                     case boundary_condition::PLANE_WAVE_E:
+                    case boundary_condition::EH_FIELD:
                         for (size_t i = 0; i < num_fluxes; i++)
                             state.bndsrcs_decomp_table_cpu.push_back(gfnum*num_fluxes+i);
                         break;
diff --git a/src/maxwell/maxwell_interface.cpp b/src/maxwell/maxwell_interface.cpp
index 32899c4ea9f4cb4a96eff973df1865a2bfda578c..5f3921dc6c12f0e76c41ff46348c437bb9768955 100644
--- a/src/maxwell/maxwell_interface.cpp
+++ b/src/maxwell/maxwell_interface.cpp
@@ -489,6 +489,7 @@ material_params_gpu::copyin(const material_params& mp)
 #define BCOND_TYPE_HPLANEW      "plane_wave_H"
 #define BCOND_TYPE_EFIELD       "E_field"
 #define BCOND_TYPE_HFIELD       "H_field"
+#define BCOND_TYPE_EHFIELD      "EH_field"
 #define BCOND_TYPE_SURFCURR     "surface_current"
 
 #define SEC_IFCONDS             "ifaceconds"
@@ -581,7 +582,7 @@ parameter_loader::validate_boundary_condition(const model&, int tag) const
     if (bc_kind == BCOND_TYPE_IMPEDANCE)
         return true;
 
-    if (bc_kind == BCOND_TYPE_EPLANEW)
+    if (bc_kind == BCOND_TYPE_EPLANEW || bc_kind == BCOND_TYPE_EHFIELD)
     {
         if (bc[BC_SOURCE].valid())
             return true;
@@ -791,6 +792,9 @@ parameter_loader::boundary_type(int tag) const
     
     if (kind_str == BCOND_TYPE_HFIELD)
         return boundary_condition::H_FIELD;
+    
+    if (kind_str == BCOND_TYPE_EHFIELD)
+        return boundary_condition::EH_FIELD;
 
     if (kind_str == BCOND_TYPE_SURFCURR)
         return boundary_condition::SURFACE_CURRENT;
@@ -846,6 +850,15 @@ parameter_loader::eval_boundary_source(int tag, const point_3d& pt, double t) co
     return ret;
 }
 
+std::pair<vec3d, vec3d>
+parameter_loader::eval_boundary_source_eh(int tag, const point_3d& pt, double t) const
+{
+    vec3d retE, retH;
+    sol::tie(retE(0), retE(1), retE(2), retH(0), retH(1), retH(2)) =
+        lua[SEC_BCONDS][tag][BC_SOURCE](tag, pt.x(), pt.y(), pt.z(), t);
+    return std::make_pair(retE, retH);
+}
+
 vec3d
 parameter_loader::eval_interface_source(int tag, const point_3d& pt, double t) const
 {