diff --git a/include/trackcpp/auxiliary.h b/include/trackcpp/auxiliary.h
index 35e2956..55d7c91 100644
--- a/include/trackcpp/auxiliary.h
+++ b/include/trackcpp/auxiliary.h
@@ -38,7 +38,8 @@ class PassMethodsClass {
     static const int pm_kickmap_pass                   = 8;
     static const int pm_matrix_pass                    = 9;
     static const int pm_drift_g2l_pass                 = 10;
-    static const int pm_nr_pms                         = 11;  // counter for number of passmethods
+    static const int pm_field3d_pass                   = 11;
+    static const int pm_nr_pms                         = 12;  // counter for number of passmethods
     static const std::vector<int> time_aware_passmethods;
     PassMethodsClass() {
         passmethods.push_back("identity_pass");
@@ -52,6 +53,7 @@ class PassMethodsClass {
         passmethods.push_back("kicktable_pass");
         passmethods.push_back("matrix_pass");
         passmethods.push_back("drift_g2l_pass");
+        passmethods.push_back("field3d_pass");
     }
     int size() const { return passmethods.size(); }
     std::string operator[](const int i) const { return passmethods[i]; }
@@ -76,7 +78,8 @@ struct PassMethod {
         pm_kickmap_pass                   = 8,
         pm_matrix_pass                    = 9,
         pm_drift_g2l_pass                 = 10,
-        pm_nr_pms                         = 11,
+        pm_field3d_pass                   = 11,
+        pm_nr_pms                         = 12,
     };
 };
 
@@ -94,6 +97,7 @@ const std::vector<std::string> pm_dict = {
     "kicktable_pass",
     "matrix_pass",
     "drift_g2l_pass",
+    "field3d_pass",
 };
 
 struct RadiationState {
diff --git a/include/trackcpp/coefmatrix.h b/include/trackcpp/coefmatrix.h
new file mode 100644
index 0000000..6886da2
--- /dev/null
+++ b/include/trackcpp/coefmatrix.h
@@ -0,0 +1,124 @@
+// TRACKCPP - Particle tracking code
+// Copyright (C) 2015  LNLS Accelerator Physics Group
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef COEFMATRIX_H
+#define COEFMATRIX_H
+
+#include <cstddef>
+#include <vector>
+
+#ifndef SWIG
+template<typename T>
+class MatrixRow {
+public:
+    explicit MatrixRow(T* ptr)
+        : ptr_(ptr)
+    {}
+
+    T& operator[](size_t j)
+    {
+        return ptr_[j];
+    }
+
+private:
+    T* ptr_;
+};
+#endif
+
+
+class CoefMatrix {
+public:
+
+#ifndef SWIG
+    using Row = MatrixRow<double>;
+    using ConstRow = MatrixRow<const double>;
+#endif
+
+    CoefMatrix() = default;
+
+    CoefMatrix(size_t rows, size_t cols)
+        : rows_(rows),
+          cols_(cols),
+          data_(rows * cols, 0.0)
+    {}
+
+    void resize(size_t rows, size_t cols)
+    {
+        rows_ = rows;
+        cols_ = cols;
+        data_.resize(rows * cols);
+    }
+
+    double& operator()(size_t i, size_t j)
+    {
+        return data_[i * cols_ + j];
+    }
+
+    const double& operator()(size_t i, size_t j) const
+    {
+        return data_[i * cols_ + j];
+    }
+
+#ifndef SWIG
+    Row operator[](size_t i)
+    {
+        return Row(data_.data() + i * cols_);
+    }
+
+    ConstRow operator[](size_t i) const
+    {
+        return ConstRow(data_.data() + i * cols_);
+    }
+#endif
+
+    size_t rows() const
+    {
+        return rows_;
+    }
+
+    size_t cols() const
+    {
+        return cols_;
+    }
+
+    size_t size() const
+    {
+        return data_.size();
+    }
+
+    bool empty() const
+    {
+        return data_.empty();
+    }
+
+    double* data()
+    {
+        return data_.data();
+    }
+
+    const double* data() const
+    {
+        return data_.data();
+    }
+
+private:
+
+    size_t rows_ = 0;
+    size_t cols_ = 0;
+    std::vector<double> data_;
+};
+
+#endif
\ No newline at end of file
diff --git a/include/trackcpp/elements.h b/include/trackcpp/elements.h
index 94942d9..4ab1a05 100644
--- a/include/trackcpp/elements.h
+++ b/include/trackcpp/elements.h
@@ -24,6 +24,7 @@
 #include <string>
 #include <fstream>
 #include <cfloat>
+#include "coefmatrix.h"
 
 
 class Element {
@@ -63,7 +64,12 @@ class Element {
   double        phase_lag   = 0;  // [rad]
   int           kicktable_idx = -1;   // index of kickmap object in kicktable_list
   double        rescale_kicks = 1.0;  // for kickmaps
-
+  double        field3d_hori_ks = 0;       // [1/m]
+  double        field3d_hori_kx = 0;       // [1/m]
+  double        field3d_hori_s0 = 0;       // [m]
+  CoefMatrix field3d_hori_coefs_cos;
+  CoefMatrix field3d_hori_coefs_sin;
+  
   std::vector<double> polynom_a = default_polynom;
   std::vector<double> polynom_b = default_polynom;
   Matrix              matrix66 = Matrix(6);
@@ -111,6 +117,10 @@ class Element {
   static Element sextupole  (const std::string& fam_name_, const double& length_, const double& S_, const int nr_steps_ = 5);
   static Element rfcavity   (const std::string& fam_name_, const double& length_, const double& frequency_, const double& voltage_, const double& phase_lag_);
   static Element kickmap    (const std::string& fam_name_, const std::string& kicktable_fname_, const int nr_steps_ = 20, const double& rescale_length_ = 1.0, const double& rescale_kicks_ = 1.0);
+  static Element field3d    (const std::string& fam_name_, const double& length_, const double& field3d_hori_s0_, const double& field3d_hori_kx_, const double& field3d_hori_ks_,
+                             const CoefMatrix& field3d_hori_coefs_cos_,
+                             const CoefMatrix& field3d_hori_coefs_sin_,
+                             const int nr_steps_ = 40);
 
   bool operator==(const Element& o) const;
   bool operator!=(const Element& o) const { return !(*this == o); };
@@ -132,5 +142,9 @@ void initialize_quadrupole(Element& element, const double& K, const int& nr_step
 void initialize_sextupole(Element& element, const double& S, const int& nr_steps);
 void initialize_rfcavity(Element& element, const double& frequency, const double& voltage, const double& phase_lag);
 void initialize_kickmap(Element& element, const int& kicktable_idx, const int& nr_steps, const double &rescale_kicks);
+void initialize_field3d(Element& element, const double& field3d_hori_s0_, const double& field3d_hori_kx_, const double& field3d_hori_ks_,
+                        const CoefMatrix& field3d_hori_coefs_cos_,
+                        const CoefMatrix& field3d_hori_coefs_sin_,
+                        const int nr_steps_);
 
 #endif
diff --git a/include/trackcpp/field3d.hpp b/include/trackcpp/field3d.hpp
new file mode 100644
index 0000000..b63ce5c
--- /dev/null
+++ b/include/trackcpp/field3d.hpp
@@ -0,0 +1,214 @@
+// TRACKCPP - Particle tracking code
+// Copyright (C) 2015  LNLS Accelerator Physics Group
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#ifndef _FIELD3D_H
+#define _FIELD3D_H
+
+#include <vector>
+#include <string>
+#include <fstream>
+#include <cmath>
+#include "auxiliary.h"
+#include "pos.h"
+#include "coefmatrix.h"
+
+
+inline void calc_sdependence(
+    const CoefMatrix& hori_coefs_cos,
+    const CoefMatrix& hori_coefs_sin,
+    const double& hori_ks,
+    const double& s,
+    CoefMatrix& sdependence
+)
+{
+    size_t M = hori_coefs_cos.rows();
+    size_t N = hori_coefs_cos.cols();
+    for (int n = 1; n <= N; ++n) {
+        const double nks = n * hori_ks;
+        const double coss = cos(nks * s);
+        const double sins = sin(nks * s);
+        for (int m = 1; m <= M; ++m) {
+            const double fac1 = hori_coefs_cos[m - 1][n - 1];
+            const double fac2 = -hori_coefs_sin[m - 1][n - 1];
+            sdependence[m-1][n-1] = fac1 * sins + fac2 * coss;
+        }
+    }
+}
+
+template <typename T>
+void ay_inty_day_dx(
+    const double& brho, const double& hori_kx, const double& hori_ks,
+    const CoefMatrix& sdependence,
+    const T& x, const T& y, const double& s,
+    T& ayn, T& idayn
+)
+{
+    size_t M = sdependence.rows();
+    size_t N = sdependence.cols();
+
+    std::vector<double> nks(N), nks2(N);
+    for (int n = 1; n <= N; ++n) {
+        nks[n-1] = n * hori_ks;
+        nks2[n-1] = nks[n-1] * nks[n-1];
+    }
+
+    for (int m = 1; m <= M; ++m) {
+        const double mkx = m * hori_kx;
+        const double mkx2 = mkx * mkx;
+        const T sinx = sin(mkx * x);
+        const T cosx = cos(mkx * x);
+        for (int n = 1; n <= N; ++n) {
+            const double hori_ky = std::sqrt(mkx2 + nks2[n-1]);
+            double ratio = mkx / (nks[n-1] * hori_ky);
+            ayn += sinx * sinh(hori_ky * y) * ratio * sdependence[m-1][n-1];
+
+            ratio = mkx2 / (nks[n-1] * hori_ky * hori_ky);
+            idayn += cosx * (cosh(hori_ky * y) - 1.0) * ratio * sdependence[m-1][n-1];
+        }
+    }
+    ayn /= brho;
+    idayn /= brho;
+}
+
+template <typename T>
+void ax_intx_dax_dy(
+    const double& brho, const double& hori_kx, const double& hori_ks,
+    const CoefMatrix& sdependence,
+    const T& x, const T& y, const double& s,
+    T& axn, T& idaxn
+)
+{
+    size_t M = sdependence.rows();
+    size_t N = sdependence.cols();
+
+    std::vector<double> nks(N), nks2(N);
+    for (int n = 1; n <= N; ++n) {
+        nks[n-1] = n * hori_ks;
+        nks2[n-1] = nks[n-1] * nks[n-1];
+    }
+
+    for (int m = 1; m <= M; ++m) {
+        const double mkx = m * hori_kx;
+        const double mkx2 = mkx * mkx;
+        const T sinx = sin(mkx * x);
+        const T cosx = cos(mkx * x);
+        for (int n = 1; n <= N; ++n) {
+            const double hori_ky = std::sqrt(mkx2 + nks2[n-1]);
+            double ratio = 1 / nks[n-1];
+            axn += cosx * cosh(hori_ky * y) * ratio * sdependence[m-1][n-1];
+
+            ratio = hori_ky / (nks[n-1] * mkx);
+            idaxn += sinx * sinh(hori_ky * y) * ratio * sdependence[m-1][n-1];
+        }
+    }
+    axn /= brho;
+    idaxn /= brho;
+}
+
+template <typename T>
+void exp_h1_s(T& s, T step) {
+    s += step / 2.0;
+}
+
+template <typename T>
+void exp_h2_y(Pos<T>& map, const T& pnorm, double step) {
+    map.ry += 0.5*step*pnorm*map.py;
+}
+
+
+template <typename T>
+void exp_h2_z(Pos<T>& map, const T& pnorm, double step) {
+    map.dl += 0.25*step*pnorm*pnorm*(map.py * map.py);
+}
+
+
+template <typename T>
+void exp_h3_x(Pos<T>& map, const T& pnorm, double step) {
+    map.rx += step*pnorm*map.px;
+}
+
+template <typename T>
+void exp_h3_z(Pos<T>& map, const T& pnorm, double step) {
+    map.dl += 0.5*step*pnorm*pnorm*(map.px * map.px);
+}
+
+template <typename T>
+void prop_h1(Pos<T>& map, double& s, double step) {
+    exp_h1_s(s, step);
+}
+
+template <typename T>
+void prop_h2(Pos<T>& map, const T& pnorm, double step) {
+    exp_h2_y(map, pnorm, step);
+    exp_h2_z(map, pnorm, step);
+}
+
+template <typename T>
+void prop_h3(Pos<T>& map, const T& pnorm, double step) {
+    exp_h3_x(map, pnorm, step);
+    exp_h3_z(map, pnorm, step);
+}
+
+template <typename T>
+void prop_ix(const double& brho, const double& hori_kx, const double& hori_ks,
+             const CoefMatrix& sdependence,
+             Pos<T>& map, double s, int sign, double step) {
+    T axn = 0.0;
+    T idaxn = 0.0;
+    ax_intx_dax_dy(
+        brho, hori_kx, hori_ks, sdependence,
+        map.rx, map.ry, s, axn, idaxn
+    );
+    map.px += axn * sign * -1.0;
+    map.py += idaxn * sign * -1.0;
+}
+
+template <typename T>
+void prop_iy(const double& brho, const double& hori_kx, const double& hori_ks,
+             const CoefMatrix& sdependence,
+             Pos<T>& map, double s, int sign, double step) {
+    T ayn = 0.0;
+    T idayn = 0.0;
+    ay_inty_day_dx(
+        brho, hori_kx, hori_ks, sdependence,
+        map.rx, map.ry, s, ayn, idayn
+    );
+    map.px += idayn * sign * -1.0;
+    map.py += ayn * sign * -1.0;
+}
+
+template <typename T>
+void prop_step(const double& brho, const double& hori_kx, const double& hori_ks,
+               const CoefMatrix& hori_coefs_cos,
+               const CoefMatrix& hori_coefs_sin,
+               Pos<T>& map, const T& pnorm, double& s, double step) {
+    prop_h1(map, s, step);
+
+    CoefMatrix sdependence(hori_coefs_cos.rows(), hori_coefs_cos.cols());
+    calc_sdependence(hori_coefs_cos, hori_coefs_sin, hori_ks, s, sdependence);
+
+    prop_iy(brho, hori_kx, hori_ks, sdependence, map, s, +1, step);
+    prop_h2(map, pnorm, step);
+    prop_iy(brho, hori_kx, hori_ks, sdependence, map, s, -1, step);
+    prop_ix(brho, hori_kx, hori_ks, sdependence, map, s, +1, step);
+    prop_h3(map, pnorm, step);
+    prop_ix(brho, hori_kx, hori_ks, sdependence, map, s, -1, step);
+    prop_iy(brho, hori_kx, hori_ks, sdependence, map, s, +1, step);
+    prop_h2(map, pnorm, step);
+    prop_iy(brho, hori_kx, hori_ks, sdependence, map, s, -1, step);
+    prop_h1(map, s, step);
+}
+
+#endif
diff --git a/include/trackcpp/passmethods.h b/include/trackcpp/passmethods.h
index 51a38aa..35f9e67 100644
--- a/include/trackcpp/passmethods.h
+++ b/include/trackcpp/passmethods.h
@@ -66,6 +66,7 @@ template <typename T> Status::type pm_thinsext_pass              (Pos<T> &pos, c
 template <typename T> Status::type pm_kickmap_pass               (Pos<T> &pos, const Element &elem, const Accelerator& accelerator);
 template <typename T> Status::type pm_matrix_pass                (Pos<T> &pos, const Element &elem, const Accelerator& accelerator);
 template <typename T> Status::type pm_drift_g2l_pass             (Pos<T> &pos, const Element &elem, const Accelerator& accelerator);
+template <typename T> Status::type pm_field3d_pass               (Pos<T> &pos, const Element &elem, const Accelerator& accelerator);
 
 #include "passmethods.hpp"
 
diff --git a/include/trackcpp/passmethods.hpp b/include/trackcpp/passmethods.hpp
index 0931064..9baccdb 100644
--- a/include/trackcpp/passmethods.hpp
+++ b/include/trackcpp/passmethods.hpp
@@ -31,6 +31,7 @@
 #include "auxiliary.h"
 #include "tpsa.h"
 #include "linalg.h"
+#include "field3d.hpp"
 #include <cmath>
 #include <algorithm>
 
@@ -114,6 +115,7 @@ Status::type kicktablethinkick(Pos<T>& pos, const int& kicktable_idx,
   return status;
 }
 
+
 template <typename T>
 void matthinkick(Pos<T> &pos, const Matrix &m) {
 
@@ -469,6 +471,23 @@ Status::type pm_kickmap_pass(Pos<T> &pos, const Element &elem,
   return status;
 }
 
+template <typename T>
+Status::type pm_field3d_pass(Pos<T> &pos, const Element &elem,
+                            const Accelerator& accelerator) {
+                  
+  global_2_local(pos, elem);
+  const double brho = get_magnetic_rigidity(accelerator.energy);
+  const double gamma = accelerator.energy / electron_rest_energy_eV;
+  double step   = elem.length / float(elem.nr_steps);
+  double s0 = elem.field3d_hori_s0;
+  T pnorm = 1 / (1 + pos.de);
+  for (int i=0; i<elem.nr_steps; ++i){
+    prop_step(-1*brho, elem.field3d_hori_kx, elem.field3d_hori_ks, elem.field3d_hori_coefs_cos, elem.field3d_hori_coefs_sin, pos, pnorm, s0, step);
+  }
+  local_2_global(pos, elem);
+  return Status::success;
+}
+
 template <typename T>
 Status::type pm_matrix_pass(Pos<T> &pos, const Element &elem,
                             const Accelerator& accelerator) {
diff --git a/include/trackcpp/tracking.h b/include/trackcpp/tracking.h
index 0542992..1c585ef 100644
--- a/include/trackcpp/tracking.h
+++ b/include/trackcpp/tracking.h
@@ -94,6 +94,9 @@ Status::type track_elementpass (
     case PassMethod::pm_drift_g2l_pass:
         if ((status = pm_drift_g2l_pass<T>(orig_pos, el, accelerator)) != Status::success) return status;
         break;
+    case PassMethod::pm_field3d_pass:
+        if ((status = pm_field3d_pass<T>(orig_pos, el, accelerator)) != Status::success) return status;
+        break;
     default:
         return Status::passmethod_not_defined;
     }
diff --git a/python_package/interface.cpp b/python_package/interface.cpp
index 7e99942..faf8481 100644
--- a/python_package/interface.cpp
+++ b/python_package/interface.cpp
@@ -241,6 +241,14 @@ Element kickmap_wrapper(const std::string& fam_name_,  const std::string& kickta
     return Element::kickmap(fam_name_, kicktable_fname_, nr_steps_, rescale_length_, rescale_kicks_);
 }
 
+Element field3d_wrapper(const std::string& fam_name_, const double& length_,
+                        const double& hori_s0_, const double& hori_kx_, const double& hori_ks_, 
+                        const CoefMatrix& hori_coefs_cos_,
+                        const CoefMatrix& hori_coefs_sin_,
+                        const int nr_steps_) {
+    return Element::field3d(fam_name_, length_, hori_s0_, hori_kx_, hori_ks_, hori_coefs_cos_, hori_coefs_sin_, nr_steps_);
+}
+
 Status::type read_flat_file_wrapper(String& fname, Accelerator& accelerator, bool file_flag) {
   return read_flat_file(fname.data, accelerator, file_flag);
 }
diff --git a/python_package/interface.h b/python_package/interface.h
index e35e092..40f6595 100644
--- a/python_package/interface.h
+++ b/python_package/interface.h
@@ -26,6 +26,7 @@
 #include <trackcpp/diffusion_matrix.h>
 #include <trackcpp/pos.h>
 #include <trackcpp/naff.h>
+#include <trackcpp/coefmatrix.h>
 
 struct LinePassArgs {
     unsigned int element_offset;
@@ -91,6 +92,10 @@ Element quadrupole_wrapper(const std::string& fam_name_, const double& length_,
 Element sextupole_wrapper(const std::string& fam_name_, const double& length_, const double& S_, const int nr_steps_ = 5);
 Element rfcavity_wrapper(const std::string& fam_name_, const double& length_, const double& frequency_, const double& voltage_, const double& phase_lag);
 Element kickmap_wrapper(const std::string& fam_name_,  const std::string& kicktable_fname_, const int nr_steps_ = 20, const double& rescale_length = 1.0, const double& rescale_kicks = 1.0);
+Element field3d_wrapper(const std::string& fam_name_, const double& length_, const double& field3d_hori_s0_, const double& field3d_hori_kx_,
+                        const double& field3d_hori_ks_, const CoefMatrix& field3d_hori_coefs_cos_,
+                        const CoefMatrix& field3d_hori_coefs_sin_,
+                        const int nr_steps_);
 Element rbend_wrapper(const std::string& fam_name_, const double& length_,
                       const double& angle_, const double& angle_in_, const double& angle_out_,
                       const double& gap_, const double& fint_in_, const double& fint_out_,
diff --git a/python_package/trackcpp.i b/python_package/trackcpp.i
index 766a189..e3b6a1e 100644
--- a/python_package/trackcpp.i
+++ b/python_package/trackcpp.i
@@ -29,6 +29,7 @@
 #include <trackcpp/diffusion_matrix.h>
 #include <trackcpp/naff.h>
 #include <trackcpp/linalg.h>
+#include <trackcpp/coefmatrix.h>
 #include "interface.h"
 %}
 
@@ -124,6 +125,7 @@ double get_double_max() {
 %include "../include/trackcpp/diffusion_matrix.h"
 %include "../include/trackcpp/naff.h"
 %include "../include/trackcpp/linalg.h"
+%include "../include/trackcpp/coefmatrix.h"
 %include "interface.h"
 
 %template(CppDoublePos) Pos<double>;
diff --git a/src/elements.cpp b/src/elements.cpp
index 5d96fb1..c97b403 100644
--- a/src/elements.cpp
+++ b/src/elements.cpp
@@ -19,6 +19,7 @@
 #include <trackcpp/kicktable.h>
 #include <cfloat>
 #include <cstring> // necessary for memcmp
+#include <trackcpp/coefmatrix.h>
 
 const std::vector<double> Element::default_polynom = std::vector<double>(3,0);
 
@@ -156,6 +157,15 @@ Element Element::kickmap (const std::string& fam_name_, const std::string& kickt
   return e;
 }
 
+Element Element::field3d (const std::string& fam_name_, const double& length_, const double& hori_s0_, const double& hori_kx_, const double& hori_ks_,
+                          const CoefMatrix& hori_coefs_cos_,
+                          const CoefMatrix& hori_coefs_sin_,
+                          const int nr_steps_) {
+  Element e = Element(fam_name_, length_);
+    initialize_field3d(e, hori_s0_, hori_kx_, hori_ks_, hori_coefs_cos_, hori_coefs_sin_, nr_steps_);
+  return e;
+}
+
 
 void print_polynom(std::ostream& out, const std::string& label, const std::vector<double>& polynom) {
   int order = 0;
@@ -308,3 +318,16 @@ void initialize_kickmap(Element& element, const int& kicktable_idx, const int& n
     element.kicktable_idx = kicktable_idx;
     element.rescale_kicks = rescale_kicks;
 }
+
+void initialize_field3d(Element& element, const double& hori_s0, const double& hori_kx, const double& hori_ks, 
+                        const CoefMatrix& hori_coefs_cos, 
+                        const CoefMatrix& hori_coefs_sin, 
+                        const int nr_steps) {
+    element.pass_method = PassMethod::pm_field3d_pass;
+    element.nr_steps = nr_steps;
+    element.field3d_hori_coefs_cos = hori_coefs_cos;
+    element.field3d_hori_coefs_sin = hori_coefs_sin;
+    element.field3d_hori_kx = hori_kx;
+    element.field3d_hori_ks = hori_ks;
+    element.field3d_hori_s0 = hori_s0;
+}
diff --git a/src/flat_file.cpp b/src/flat_file.cpp
index a6ecc4b..21be2cf 100644
--- a/src/flat_file.cpp
+++ b/src/flat_file.cpp
@@ -30,9 +30,11 @@ static int process_rad_property(std::istringstream& ss);
 static std::string get_boolean_string(bool value);
 static bool has_matrix66(const Matrix& r);
 static bool has_polynom(const std::vector<double>& p);
+static bool has_coeffs(const CoefMatrix& matrix);
 static void write_6d_vector(std::ostream& fp, const std::string& label, const double* t);
 static void write_6d_vector(std::ostream& fp, const std::string& label, const std::vector<double>& t);
 static void write_polynom(std::ostream& fp, const std::string& label, const std::vector<double>& p);
+static void write_coeffs(std::ostream& fp, const std::string& label, const CoefMatrix& matrix);
 static void synchronize_polynomials(Element& e);
 static void read_polynomials(std::ifstream& fp, Element& e);
 static void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator);
@@ -127,6 +129,13 @@ void write_flat_file_trackcpp(std::ostream& fp, const Accelerator& accelerator)
     if (e.angle_in != 0) { fp << std::setw(pw) << "angle_in" << e.angle_in << '\n'; }
     if (e.angle_out != 0) { fp << std::setw(pw) << "angle_out" << e.angle_out << '\n'; }
     if (e.rescale_kicks != 1.0) { fp << std::setw(pw) << "rescale_kicks" << e.rescale_kicks << '\n'; }
+    if (e.field3d_hori_kx != 0) { fp << std::setw(pw) << "kx" << e.field3d_hori_kx << '\n'; }
+    if (e.field3d_hori_ks != 0) { fp << std::setw(pw) << "ks" << e.field3d_hori_ks << '\n'; }
+    if (e.field3d_hori_s0 != 0) { fp << std::setw(pw) << "s0" << e.field3d_hori_s0 << '\n'; }
+    if (has_coeffs(e.field3d_hori_coefs_cos)) write_coeffs(fp, "coefs_cos", e.field3d_hori_coefs_cos);
+    if (has_coeffs(e.field3d_hori_coefs_sin)) write_coeffs(fp, "coefs_sin", e.field3d_hori_coefs_sin);
+    
+    
     if (e.has_t_in) write_6d_vector(fp, "t_in", e.t_in);
     if (e.has_t_out) write_6d_vector(fp, "t_out", e.t_out);
     if (e.has_r_in) {
@@ -226,6 +235,9 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator)
     if (cmd.compare("angle_in")    == 0) { ss >> e.angle_in;  continue; }
     if (cmd.compare("angle_out")   == 0) { ss >> e.angle_out; continue; }
     if (cmd.compare("rescale_kicks")   == 0) { ss >> e.rescale_kicks; continue; }
+    if (cmd.compare("kx")       == 0) { ss >> e.field3d_hori_kx;     continue; }
+    if (cmd.compare("ks")       == 0) { ss >> e.field3d_hori_ks;     continue; }
+    if (cmd.compare("s0")       == 0) { ss >> e.field3d_hori_s0;     continue; }
     if (cmd.compare("t_in")      == 0) { for(auto i=0; i<6; ++i) ss >> e.t_in[i]; e.reflag_t_in(); continue; }
     if (cmd.compare("t_out")     == 0) { for(auto i=0; i<6; ++i) ss >> e.t_out[i]; e.reflag_t_out(); continue; }
     if (cmd.compare("rx|r_in")   == 0) { for(auto i=0; i<6; ++i) ss >> e.r_in[0*6+i]; e.reflag_r_in(); continue; }
@@ -303,6 +315,48 @@ Status::type read_flat_file_trackcpp(std::istream& fp, Accelerator& accelerator)
       synchronize_polynomials(e);
       continue;
     }
+
+    if (cmd.compare("coefs_cos") == 0 || 
+        cmd.compare("coefs_sin") == 0) {
+        CoefMatrix Element::* M = nullptr;
+      
+        if (cmd == "coefs_cos")
+        M = &Element::field3d_hori_coefs_cos;
+
+        if (cmd == "coefs_sin")
+        M = &Element::field3d_hori_coefs_sin;
+
+        std::vector<unsigned int> row;
+        std::vector<unsigned int> col;
+        std::vector<double>       val;
+
+        unsigned int max_row = 0;
+        unsigned int max_col = 0;
+
+        while (!ss.eof()) {
+            unsigned int r, c;
+            double v;
+            ss >> r >> c >> v;
+
+            if (ss.eof()) break;
+
+            row.push_back(r);
+            col.push_back(c);
+            val.push_back(v);
+
+            if (r + 1 > max_row) max_row = r + 1;
+            if (c + 1 > max_col) max_col = c + 1;
+        }
+
+        if (max_row > 0 && max_col > 0) {
+            (e.*M).resize(max_row, max_col);
+
+          for (unsigned int k = 0; k < val.size(); ++k)
+              (e.*M)[row[k]][col[k]] = val[k];
+        }
+
+    continue;
+  }
     if (line.size()<2) continue;
     return Status::flat_file_error;
   }
@@ -490,6 +544,11 @@ static bool has_matrix66(const Matrix& m) {
   return false;
 }
 
+static bool has_coeffs(const CoefMatrix& matrix)
+{
+    return !matrix.empty();
+}
+
 static bool has_polynom(const std::vector<double>& p) {
   for (int i=0; i<p.size(); ++i)
     if (p[i] != 0)
@@ -512,6 +571,25 @@ static void write_6d_vector(std::ostream& fp, const std::string& label, const st
   fp << '\n';
 }
 
+static void write_coeffs(std::ostream& fp,
+                         const std::string& label,
+                         const CoefMatrix& M)
+{
+    fp << std::setw(pw) << label;
+
+    for (size_t i = 0; i < M.rows(); ++i) {
+        for (size_t j = 0; j < M.cols(); ++j) {
+
+            fp.unsetf(std::ios_base::showpos);
+            fp << i << ' ' << j << ' ';
+            fp.setf(std::ios_base::showpos);
+            fp << M[i][j] << ' ';
+        }
+    }
+
+    fp << '\n';
+}
+
 static void write_polynom(std::ostream& fp, const std::string& label, const std::vector<double>& p) {
   fp << std::setw(pw) << label;
   for (int i=0; i<p.size(); ++i)