diff --git a/CMakeLists.txt b/CMakeLists.txt
index 43483efaf..fceda4550 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -402,6 +402,7 @@ add_library(crpropa SHARED
   src/magneticField/turbulentField/SimpleGridTurbulence.cpp
   src/magneticField/TF17Field.cpp
   src/magneticField/UF23Field.cpp
+  src/magneticField/KST24Field.cpp
   src/magneticField/CMZField.cpp
   src/advectionField/AdvectionField.cpp
   src/massDistribution/ConstantDensity.cpp
diff --git a/include/CRPropa.h b/include/CRPropa.h
index bac540bb9..e8491cdb8 100644
--- a/include/CRPropa.h
+++ b/include/CRPropa.h
@@ -65,6 +65,7 @@
 #include "crpropa/magneticField/QuimbyMagneticField.h"
 #include "crpropa/magneticField/TF17Field.h"
 #include "crpropa/magneticField/UF23Field.h"
+#include "crpropa/magneticField/KST24Field.h"
 #include "crpropa/magneticField/CMZField.h"
 #include "crpropa/magneticField/turbulentField/GridTurbulence.h"
 #include "crpropa/magneticField/turbulentField/HelicalGridTurbulence.h"
diff --git a/include/crpropa/magneticField/KST24Field.h b/include/crpropa/magneticField/KST24Field.h
new file mode 100644
index 000000000..46edb889b
--- /dev/null
+++ b/include/crpropa/magneticField/KST24Field.h
@@ -0,0 +1,154 @@
+#ifndef _KST24_GMF_H_
+#define _KST24_GMF_H_
+
+#include <vector>
+#include "crpropa/magneticField/MagneticField.h"
+
+namespace crpropa {
+
+/* 
+The C++ implementation of the GMF model KST24 (A.Korochkin, D.Semikoz, P.Tinyakov 2024)
+The model was presented in arXiv:2407.02148 and published in A&A
+If you use the model, please cite A&A, 693, A284 (2025)
+
+In KST24 GMF model the position of the Solar System is at {-8.2 kpc, 0, 0}
+The Galactic Center is at {0, 0, 0}
+z-axis points to the North pole
+*/
+
+
+class KST24Field : public MagneticField
+{
+private:
+	double north_tor_B_gauss;
+	double north_tor_zmin_kpc;
+	double north_tor_zmax_kpc;
+	double north_tor_rmin_kpc;
+	double north_tor_rmax_kpc;
+
+	double south_tor_B_gauss;
+	double south_tor_zmin_kpc;
+	double south_tor_zmax_kpc;
+	double south_tor_rmin_kpc;
+	double south_tor_rmax_kpc;
+
+	double Xfield_B_gauss;
+	double Xfield_rmin_kpc;
+	double Xfield_rmax_kpc;
+	double Xfield_theta_deg;
+	double Xfield_theta_rad;
+
+	double LB_B_gauss;
+	double LB_lB_deg;
+	double LB_bB_deg;
+	double LB_rmin_kpc; 
+	double LB_dr_kpc;
+	double LB_x0_kpc;
+	double LB_y0_kpc;
+	double LB_z0_kpc;
+	double LB_Bdir_x;
+	double LB_Bdir_y;
+	double LB_Bdir_z;
+
+	double ScutumArm_B_gauss;
+	double ScutumArm_pitch_deg;
+	double ScutumArm_phi0_deg;
+	double ScutumArm_x_shift_kpc;
+	double ScutumArm_y_shift_kpc;
+	double ScutumArm_arc_radius1_kpc;
+	double ScutumArm_arc_radius2_kpc;
+	double ScutumArm_arc_eps;
+	double ScutumArm_arc_div_deg; 
+	double ScutumArm_rmin_kpc;
+	double ScutumArm_rmax_kpc;
+	double ScutumArm_zmin_kpc;
+	double ScutumArm_zmax_kpc;
+
+	double CarinaSagittariusArm_B_gauss;
+	double CarinaSagittariusArm_pitch_deg;
+	double CarinaSagittariusArm_phi0_deg;
+	double CarinaSagittariusArm_x_shift_kpc;
+	double CarinaSagittariusArm_y_shift_kpc;
+	double CarinaSagittariusArm_arc_radius1_kpc;
+	double CarinaSagittariusArm_arc_radius2_kpc;
+	double CarinaSagittariusArm_arc_eps;
+	double CarinaSagittariusArm_arc_div_deg; 
+	double CarinaSagittariusArm_rmin_kpc;
+	double CarinaSagittariusArm_rmax_kpc;
+	double CarinaSagittariusArm_zmin_kpc;
+	double CarinaSagittariusArm_zmax_kpc;
+
+	double LocalArm_B_gauss;
+	double LocalArm_pitch_deg;
+	double LocalArm_phi0_deg;
+	double LocalArm_x_shift_kpc;
+	double LocalArm_y_shift_kpc;
+	double LocalArm_arc_radius1_kpc;
+	double LocalArm_arc_radius2_kpc;
+	double LocalArm_arc_eps;
+	double LocalArm_arc_div_deg; 
+	double LocalArm_rmin_kpc;
+	double LocalArm_rmax_kpc;
+	double LocalArm_zmin_kpc;
+	double LocalArm_zmax_kpc;
+
+	double PerseusArm1_B_gauss;
+	double PerseusArm1_pitch_deg;
+	double PerseusArm1_phi0_deg;
+	double PerseusArm1_x_shift_kpc;
+	double PerseusArm1_y_shift_kpc;
+	double PerseusArm1_arc_radius1_kpc;
+	double PerseusArm1_arc_radius2_kpc;
+	double PerseusArm1_arc_eps;
+	double PerseusArm1_arc_div_deg; 
+	double PerseusArm1_rmin_kpc;
+	double PerseusArm1_rmax_kpc;
+	double PerseusArm1_zmin_kpc;
+	double PerseusArm1_zmax_kpc;
+
+	double PerseusArm2_B_gauss;
+	double PerseusArm2_pitch_deg;
+	double PerseusArm2_phi0_deg;
+	double PerseusArm2_x_shift_kpc;
+	double PerseusArm2_y_shift_kpc;
+	double PerseusArm2_arc_radius1_kpc;
+	double PerseusArm2_arc_radius2_kpc;
+	double PerseusArm2_arc_eps;
+	double PerseusArm2_arc_div_deg; 
+	double PerseusArm2_rmin_kpc;
+	double PerseusArm2_rmax_kpc;
+	double PerseusArm2_zmin_kpc;
+	double PerseusArm2_zmax_kpc;
+
+public:
+	Vector3d getField(const Vector3d& pos) const;
+	KST24Field();
+
+	Vector3d get_toroidal(const Vector3d pos_kpc, const double tor_B_gauss, 
+						  const double tor_zmin_kpc, const double tor_zmax_kpc, 
+						  const double tor_rmin_kpc, const double tor_rmax_kpc) const;
+
+	Vector3d get_Xfield(const Vector3d pos_kpc, const double Xfield_B_gauss, 
+						const double Xfield_rmin_kpc, const double Xfield_rmax_kpc, 
+						const double Xfield_theta_rad) const;
+
+	bool is_LB(const Vector3d pos_kpc, const double LB_rmin_kpc, const double LB_dr_kpc,
+			   const double LB_x0_kpc, const double LB_y0_kpc, const double LB_z0_kpc) const;
+
+	Vector3d get_LB(const Vector3d pos_kpc, const double LB_B_gauss, 
+					const double LB_lB_deg, const double LB_bB_deg, 
+					const double LB_rmin_kpc, const double LB_dr_kpc,
+					const double LB_x0_kpc, const double LB_y0_kpc, const double LB_z0_kpc) const;
+
+	Vector3d get_logspiral(const Vector3d pos_kpc, const double B_gauss,  
+						   const double pitch_deg, const double phi0_deg,
+						   const double x_shift_kpc, const double y_shift_kpc,
+						   const double arc_radius1_kpc, const double arc_radius2_kpc, 
+						   const double arc_eps , const double arc_div_deg, 
+						   const double rmin_kpc, const double rmax_kpc, 
+						   const double zmin_kpc, const double zmax_kpc) const;
+};
+
+}// namespace crpropa 
+
+#endif /* _KST24_GMF_H_ */
\ No newline at end of file
diff --git a/python/2_headers.i b/python/2_headers.i
index d34665838..073cd23b1 100644
--- a/python/2_headers.i
+++ b/python/2_headers.i
@@ -341,6 +341,7 @@
 %include "crpropa/magneticField/PT11Field.h"
 %include "crpropa/magneticField/TF17Field.h"
 %include "crpropa/magneticField/UF23Field.h"
+%include "crpropa/magneticField/KST24Field.h"
 %include "crpropa/magneticField/ArchimedeanSpiralField.h"
 %include "crpropa/magneticField/CMZField.h"
 %include "crpropa/magneticField/turbulentField/TurbulentField.h"
diff --git a/src/magneticField/KST24Field.cpp b/src/magneticField/KST24Field.cpp
new file mode 100644
index 000000000..16f688b6b
--- /dev/null
+++ b/src/magneticField/KST24Field.cpp
@@ -0,0 +1,375 @@
+#include "crpropa/magneticField/KST24Field.h"
+
+/* 
+The C++ implementation of the GMF model KST24 (A.Korochkin, D.Semikoz, P.Tinyakov 2024)
+The model was presented in arXiv:2407.02148 and published in A&A
+If you use the model, please cite A&A, 693, A284 (2025)
+
+In KST24 GMF model the position of the Solar System is at {-8.2 kpc, 0, 0}
+The Galactic Center is at {0, 0, 0}
+z-axis points to the North pole
+*/
+
+
+namespace crpropa {
+KST24Field::KST24Field()
+{
+	north_tor_B_gauss  = -3.2e-6;
+	north_tor_zmin_kpc = 1.185;
+	north_tor_zmax_kpc = 2.1;
+	north_tor_rmin_kpc = 1.0;
+	north_tor_rmax_kpc = 9.1;
+
+	south_tor_B_gauss  = 3.2e-6;
+	south_tor_zmin_kpc = -2.5;
+	south_tor_zmax_kpc = -1.22;
+	south_tor_rmin_kpc = 1.0;
+	south_tor_rmax_kpc = 14.0;
+
+	Xfield_B_gauss   = 1.8e-6;
+	Xfield_rmin_kpc  = 1.0;
+	Xfield_rmax_kpc  = 6.2;
+	Xfield_theta_deg = 28.0;
+	Xfield_theta_rad = Xfield_theta_deg*M_PI/180.;
+
+	LB_B_gauss  = -3.5e-6;
+	LB_lB_deg   = 50.0;
+	LB_bB_deg   = 2.0;
+	LB_rmin_kpc = 0.2; 
+	LB_dr_kpc   = 0.03;
+	LB_x0_kpc   = -8.2;
+	LB_y0_kpc   = 0.095;
+	LB_z0_kpc   = -0.05;
+	LB_Bdir_x   = cos(LB_bB_deg*M_PI/180)*cos(LB_lB_deg*M_PI/180);
+	LB_Bdir_y   = cos(LB_bB_deg*M_PI/180)*sin(LB_lB_deg*M_PI/180);
+	LB_Bdir_z   = sin(LB_bB_deg*M_PI/180);
+
+
+	ScutumArm_B_gauss         = 4.9e-6;
+	ScutumArm_pitch_deg       = 20;
+	ScutumArm_phi0_deg        = -134.0;
+	ScutumArm_x_shift_kpc     = 1.0;
+	ScutumArm_y_shift_kpc     = 0.0;
+	ScutumArm_arc_radius1_kpc = 0.8;
+	ScutumArm_arc_radius2_kpc = 1.0;
+	ScutumArm_arc_eps         = 2.0;
+	ScutumArm_arc_div_deg     = 0.0;
+	ScutumArm_rmin_kpc        = 3.0;
+	ScutumArm_rmax_kpc        = 16.0;
+	ScutumArm_zmin_kpc        = -1.0;
+	ScutumArm_zmax_kpc        = 1.0;	
+
+	CarinaSagittariusArm_B_gauss         = 1.3e-6;
+	CarinaSagittariusArm_pitch_deg       = 20;
+	CarinaSagittariusArm_phi0_deg        = -80.0;
+	CarinaSagittariusArm_x_shift_kpc     = 1.37;
+	CarinaSagittariusArm_y_shift_kpc     = 0.0;
+	CarinaSagittariusArm_arc_radius1_kpc = 1.0;
+	CarinaSagittariusArm_arc_radius2_kpc = 0.79;
+	CarinaSagittariusArm_arc_eps         = 2.3;
+	CarinaSagittariusArm_arc_div_deg     = 3.0;
+	CarinaSagittariusArm_rmin_kpc        = 3.0;
+	CarinaSagittariusArm_rmax_kpc        = 15.0;
+	CarinaSagittariusArm_zmin_kpc        = -1.2;
+	CarinaSagittariusArm_zmax_kpc        = 1.2;
+
+	LocalArm_B_gauss         = -3.5e-6;
+	LocalArm_pitch_deg       = 20;
+	LocalArm_phi0_deg        = -2.2;
+	LocalArm_x_shift_kpc     = -0.15;
+	LocalArm_y_shift_kpc     = 0.0;
+	LocalArm_arc_radius1_kpc = 0.73;
+	LocalArm_arc_radius2_kpc = 1.0;
+	LocalArm_arc_eps         = 1.45;
+	LocalArm_arc_div_deg     = 0.0;
+	LocalArm_rmin_kpc        = 3.3;
+	LocalArm_rmax_kpc        = 14.0;
+	LocalArm_zmin_kpc        = -1.0;
+	LocalArm_zmax_kpc        = 1.0;
+
+	PerseusArm1_B_gauss         = -2.0e-6;
+	PerseusArm1_pitch_deg       = 20;
+	PerseusArm1_phi0_deg        = 46.0;
+	PerseusArm1_x_shift_kpc     = -1.0;
+	PerseusArm1_y_shift_kpc     = 0.0;
+	PerseusArm1_arc_radius1_kpc = 0.4;
+	PerseusArm1_arc_radius2_kpc = 1.1;
+	PerseusArm1_arc_eps         = 2.0;
+	PerseusArm1_arc_div_deg     = 0.0;
+	PerseusArm1_rmin_kpc        = 3.0;
+	PerseusArm1_rmax_kpc        = 17.0;
+	PerseusArm1_zmin_kpc        = -1.0;
+	PerseusArm1_zmax_kpc        = 1.0;
+
+	PerseusArm2_B_gauss         = -3.5e-6;
+	PerseusArm2_pitch_deg       = 20;
+	PerseusArm2_phi0_deg        = 46.0;
+	PerseusArm2_x_shift_kpc     = -1.0;
+	PerseusArm2_y_shift_kpc     = 0.0;
+	PerseusArm2_arc_radius1_kpc = 1.2;
+	PerseusArm2_arc_radius2_kpc = 1.1;
+	PerseusArm2_arc_eps         = 2.0;
+	PerseusArm2_arc_div_deg     = 0.0;
+	PerseusArm2_rmin_kpc        = 3.0;
+	PerseusArm2_rmax_kpc        = 17.0;
+	PerseusArm2_zmin_kpc        = -2.0;
+	PerseusArm2_zmax_kpc        = 2.0;
+}
+
+
+
+Vector3d KST24Field::getField(const Vector3d& pos) const
+{
+	Vector3d vals(0, 0, 0);
+	Vector3d pos_kpc = pos/kpc;
+
+	if (pos_kpc.z > north_tor_zmin_kpc) // northern halo height lower boundary
+	{
+		vals += get_toroidal(pos_kpc, north_tor_B_gauss, north_tor_zmin_kpc, north_tor_zmax_kpc, 
+							 north_tor_rmin_kpc, north_tor_rmax_kpc);
+		vals += get_Xfield(pos_kpc, Xfield_B_gauss, Xfield_rmin_kpc, Xfield_rmax_kpc, Xfield_theta_rad);
+		return vals*gauss;
+	}
+
+	if (pos_kpc.z < south_tor_zmax_kpc) // southern halo height upper boundary
+	{
+		vals += get_toroidal(pos_kpc, south_tor_B_gauss, south_tor_zmin_kpc, south_tor_zmax_kpc, 
+							 south_tor_rmin_kpc, south_tor_rmax_kpc);
+		vals += get_Xfield(pos_kpc, Xfield_B_gauss, Xfield_rmin_kpc, Xfield_rmax_kpc, Xfield_theta_rad);
+		return vals*gauss;
+	}
+
+	if (is_LB(pos_kpc, LB_rmin_kpc, LB_dr_kpc, LB_x0_kpc, LB_y0_kpc, LB_z0_kpc))
+	{
+		vals += get_LB(pos_kpc, LB_B_gauss, LB_lB_deg, LB_bB_deg, 
+					   LB_rmin_kpc, LB_dr_kpc, LB_x0_kpc, LB_y0_kpc, LB_z0_kpc);
+		return vals*gauss;
+	}
+
+	vals += get_logspiral(pos_kpc, ScutumArm_B_gauss, ScutumArm_pitch_deg, ScutumArm_phi0_deg,
+						  ScutumArm_x_shift_kpc, ScutumArm_y_shift_kpc, ScutumArm_arc_radius1_kpc,
+						  ScutumArm_arc_radius2_kpc, ScutumArm_arc_eps, ScutumArm_arc_div_deg,
+						  ScutumArm_rmin_kpc, ScutumArm_rmax_kpc, ScutumArm_zmin_kpc, ScutumArm_zmax_kpc);
+
+	vals += get_logspiral(pos_kpc, CarinaSagittariusArm_B_gauss, CarinaSagittariusArm_pitch_deg, CarinaSagittariusArm_phi0_deg,
+						  CarinaSagittariusArm_x_shift_kpc, CarinaSagittariusArm_y_shift_kpc, CarinaSagittariusArm_arc_radius1_kpc,
+						  CarinaSagittariusArm_arc_radius2_kpc, CarinaSagittariusArm_arc_eps, CarinaSagittariusArm_arc_div_deg,
+						  CarinaSagittariusArm_rmin_kpc, CarinaSagittariusArm_rmax_kpc, CarinaSagittariusArm_zmin_kpc, CarinaSagittariusArm_zmax_kpc);
+
+	vals += get_logspiral(pos_kpc, LocalArm_B_gauss, LocalArm_pitch_deg, LocalArm_phi0_deg,
+						  LocalArm_x_shift_kpc, LocalArm_y_shift_kpc, LocalArm_arc_radius1_kpc,
+						  LocalArm_arc_radius2_kpc, LocalArm_arc_eps, LocalArm_arc_div_deg,
+						  LocalArm_rmin_kpc, LocalArm_rmax_kpc, LocalArm_zmin_kpc, LocalArm_zmax_kpc);
+
+	vals += get_logspiral(pos_kpc, PerseusArm1_B_gauss, PerseusArm1_pitch_deg, PerseusArm1_phi0_deg,
+						  PerseusArm1_x_shift_kpc, PerseusArm1_y_shift_kpc, PerseusArm1_arc_radius1_kpc,
+						  PerseusArm1_arc_radius2_kpc, PerseusArm1_arc_eps, PerseusArm1_arc_div_deg,
+						  PerseusArm1_rmin_kpc, PerseusArm1_rmax_kpc, PerseusArm1_zmin_kpc, PerseusArm1_zmax_kpc);
+
+	vals += get_logspiral(pos_kpc, PerseusArm2_B_gauss, PerseusArm2_pitch_deg, PerseusArm2_phi0_deg,
+						  PerseusArm2_x_shift_kpc, PerseusArm2_y_shift_kpc, PerseusArm2_arc_radius1_kpc,
+						  PerseusArm2_arc_radius2_kpc, PerseusArm2_arc_eps, PerseusArm2_arc_div_deg,
+						  PerseusArm2_rmin_kpc, PerseusArm2_rmax_kpc, PerseusArm2_zmin_kpc, PerseusArm2_zmax_kpc);
+
+	vals += get_Xfield(pos_kpc, Xfield_B_gauss, Xfield_rmin_kpc, Xfield_rmax_kpc, Xfield_theta_rad);
+
+	return vals*gauss;
+}
+
+bool KST24Field::is_LB(const Vector3d pos_kpc, const double LB_rmin_kpc, const double LB_dr_kpc,
+					   const double LB_x0_kpc, const double LB_y0_kpc, const double LB_z0_kpc) const
+{
+	double eRx, eRy, eRz;
+	double cR;
+
+	eRx = pos_kpc.x - LB_x0_kpc;
+	eRy = pos_kpc.y - LB_y0_kpc;
+	eRz = pos_kpc.z - LB_z0_kpc;
+	cR = sqrt(eRx*eRx + eRy*eRy + eRz*eRz);
+	if (cR < LB_rmin_kpc + LB_dr_kpc)
+		return true;
+	return false;
+}
+
+Vector3d KST24Field::get_toroidal(const Vector3d pos_kpc, const double tor_B_gauss, 
+								  const double tor_zmin_kpc, const double tor_zmax_kpc, 
+								  const double tor_rmin_kpc, const double tor_rmax_kpc) const 
+{
+	Vector3d vals_gauss(0, 0, 0);
+	if ((pos_kpc.z <= tor_zmin_kpc) or (tor_zmax_kpc <= pos_kpc.z))
+		return vals_gauss;
+	
+	double cR, theta;
+	cR = sqrt(pow(pos_kpc.x, 2) + pow(pos_kpc.y, 2));
+	if ((tor_rmin_kpc <= cR) and (cR <= tor_rmax_kpc))
+	{
+		theta = atan2(pos_kpc.y, pos_kpc.x);
+		vals_gauss.x =  tor_B_gauss*sin(theta);
+		vals_gauss.y = -tor_B_gauss*cos(theta);
+		vals_gauss.z =  0;
+	}
+	return vals_gauss;
+}
+
+Vector3d KST24Field::get_Xfield(const Vector3d pos_kpc, const double Xfield_B_gauss, 
+								const double Xfield_rmin_kpc, const double Xfield_rmax_kpc, 
+								const double Xfield_theta_rad) const
+{
+	Vector3d vals_gauss(0, 0, 0);
+
+	int sign = 1;
+	if (pos_kpc.z < 0)
+		sign = -1;
+	if (fabs(pos_kpc.z) > 10)
+		return vals_gauss;
+
+	double cR, cR_0;
+	cR = sqrt(pos_kpc.x*pos_kpc.x + pos_kpc.y*pos_kpc.y);
+	cR_0 = cR - pos_kpc.z*sign*tan(Xfield_theta_rad);
+
+	if ((cR_0 < Xfield_rmin_kpc) or (Xfield_rmax_kpc < cR_0))
+		return vals_gauss;
+
+	double mgn_frc = cR_0/cR;
+	double xy_theta = atan2(pos_kpc.y, pos_kpc.x);
+	vals_gauss.x =  Xfield_B_gauss*mgn_frc*cos(xy_theta)*sign*sin(Xfield_theta_rad);
+	vals_gauss.y =  Xfield_B_gauss*mgn_frc*sin(xy_theta)*sign*sin(Xfield_theta_rad);
+	vals_gauss.z =  Xfield_B_gauss*mgn_frc*cos(Xfield_theta_rad);
+
+	return vals_gauss;
+}
+
+Vector3d KST24Field::get_LB(const Vector3d pos_kpc, const double LB_B_gauss, 
+							const double LB_lB_deg, const double LB_bB_deg, 
+							const double LB_rmin_kpc, const double LB_dr_kpc,
+							const double LB_x0_kpc, const double LB_y0_kpc, const double LB_z0_kpc) const
+{
+	Vector3d vals_gauss(0, 0, 0);
+
+	double eRx, eRy, eRz, cR;
+	eRx = pos_kpc.x - LB_x0_kpc;
+	eRy = pos_kpc.y - LB_y0_kpc;
+	eRz = pos_kpc.z - LB_z0_kpc;
+	cR = sqrt(eRx*eRx + eRy*eRy + eRz*eRz);
+ 	if ((cR < LB_rmin_kpc) or (LB_rmin_kpc + LB_dr_kpc < cR))
+		return vals_gauss;
+	eRx /= cR;
+	eRy /= cR;
+	eRz /= cR;
+
+	double fdir_x, fdir_y, fdir_z, cosTheta;
+	cosTheta = eRx*LB_Bdir_x + eRy*LB_Bdir_y + eRz*LB_Bdir_z; // expanding cross product: b(ac) - c(ab)
+	fdir_x = LB_Bdir_x - eRx*cosTheta;
+	fdir_y = LB_Bdir_y - eRy*cosTheta;
+	fdir_z = LB_Bdir_z - eRz*cosTheta;
+
+	// magnetic field amplification factor
+	double ampl, ampl_elec;
+	ampl = (1 + pow(LB_rmin_kpc, 2)/(2*LB_rmin_kpc*LB_dr_kpc + LB_dr_kpc*LB_dr_kpc));
+
+	vals_gauss.x = LB_B_gauss*ampl*fdir_x;
+	vals_gauss.y = LB_B_gauss*ampl*fdir_y;
+	vals_gauss.z = LB_B_gauss*ampl*fdir_z;
+
+	return vals_gauss;
+}
+
+Vector3d KST24Field::get_logspiral(const Vector3d pos_kpc, const double B_gauss,  
+								   const double pitch_deg, const double phi0_deg,
+								   const double x_shift_kpc, const double y_shift_kpc,
+								   const double arc_radius1_kpc, const double arc_radius2_kpc, 
+								   const double arc_eps , const double arc_div_deg, 
+								   const double rmin_kpc, const double rmax_kpc, 
+								   const double zmin_kpc, const double zmax_kpc) const
+{
+	Vector3d vals_gauss(0, 0, 0);
+
+	if ((pos_kpc.z < zmin_kpc) or (zmax_kpc < pos_kpc.z))
+		return vals_gauss;
+
+	std::vector<double> pos_v;
+	pos_v.push_back(pos_kpc.x + x_shift_kpc);
+	pos_v.push_back(pos_kpc.y + y_shift_kpc);
+	pos_v.push_back(pos_kpc.z);
+
+	double r = sqrt(pos_v[0]*pos_v[0] + pos_v[1]*pos_v[1]);
+	if ((r < rmin_kpc) or (rmax_kpc < r))
+		return vals_gauss;
+
+	double a_kpc = 3;
+	double k = tan(pitch_deg*M_PI/180.);
+	double cos_pitch = cos(pitch_deg*M_PI/180.);
+	double sin_pitch = sin(pitch_deg*M_PI/180.);
+	double phi0 = phi0_deg*M_PI/180.;
+	double arc_div_rad = arc_div_deg*M_PI/180.;
+
+	int nn;
+	double phi = atan2(pos_v[1], pos_v[0]);	
+	nn = floor((log(r/a_kpc) - k*(phi + phi0))/(2*M_PI*k));
+
+	double r1, r2;
+	r1 = a_kpc*exp(k*(phi + phi0 + 2*M_PI*nn));
+	r2 = r1*exp(k*2*M_PI);
+
+	if (fabs(r - r1) > fabs(r - r2))
+	{
+		r1 = r2;
+		nn++;  // we need outer arm
+	}
+
+
+	// perpendicular to the spiral arm axis
+	double xi, yi, dir1, dir2, dirtmp;
+	xi = r1*cos(phi);
+	yi = r1*sin(phi);
+	dir1 = k*xi - yi;
+	dir2 = k*yi + xi;
+	dirtmp = sqrt(dir1*dir1 + dir2*dir2);
+	dir1 /= dirtmp;
+	dir2 /= dirtmp;
+	// std::cout << ((xi - pos[0])*dir1 + (yi - pos[1])*dir2)/sqrt(pow((xi - pos[0]), 2) + pow((yi - pos[1]), 2)) << '\n';
+
+
+	// first order correction 
+	double delta_phi;
+	delta_phi = -((xi - pos_v[0])*dir1 + (yi - pos_v[1])*dir2)/(r1/cos_pitch);
+	r1 = a_kpc*exp(k*(phi + phi0 + delta_phi + 2*M_PI*nn));
+	xi = r1*cos(phi + delta_phi);
+	yi = r1*sin(phi + delta_phi);
+	dir1 = k*xi - yi;
+	dir2 = k*yi + xi;
+	dirtmp = sqrt(dir1*dir1 + dir2*dir2);
+	dir1 /= dirtmp;
+	dir2 /= dirtmp;
+
+
+	double d1, L1, d_at_L1;
+	d1 = sqrt(pow((xi - pos_v[0]), 2) + pow((yi - pos_v[1]), 2));
+	L1 = r1/sin_pitch;
+
+	double r_scale = 5;
+	double L1_at_r_scale, arc_r_plane, arc_r_z;
+	L1_at_r_scale = r_scale/sin_pitch;
+	arc_r_plane = std::max(arc_radius2_kpc, arc_radius2_kpc*(1 + (L1 - L1_at_r_scale)*arc_div_rad));
+	// arc_r_plane = arc_radius2_kpc;
+	arc_r_z     = std::max(arc_radius1_kpc, arc_radius1_kpc*(1 + (L1 - L1_at_r_scale)*arc_div_rad));
+	if (arc_r_plane > 1.2)
+		arc_r_plane = 1.2;
+	if (arc_r_z > 1.2)
+		arc_r_z = 1.2;
+
+
+	d_at_L1 = sqrt(pow(fabs(d1/arc_r_plane), arc_eps) + pow(fabs(pos_v[2]/arc_r_z), arc_eps));
+	if ((d_at_L1 > 1) or (d_at_L1 < 0.))
+		return vals_gauss;
+	else
+	{
+		double scaling_factor = (arc_radius1_kpc*arc_radius2_kpc)/(arc_r_plane*arc_r_z);
+		vals_gauss.x = B_gauss*dir1*scaling_factor;
+		vals_gauss.y = B_gauss*dir2*scaling_factor;
+		vals_gauss.z = 0;
+	}
+
+	return vals_gauss;
+}
+}
\ No newline at end of file