From 56c624f24bff638a4b1ef22168e84507b49ff8ec Mon Sep 17 00:00:00 2001
From: Gordon McCann <gmccann@fsu.edu>
Date: Wed, 8 Jun 2022 09:37:32 -0400
Subject: [PATCH] Add forward and reverse energy loss integrators using similar
 methods to those of the Yale SPANC integrators

---
 gwm_integrators.cpp   | 86 +++++++++++++++++++++++++++++++++++++++++++
 gwm_integrators.h     | 13 +++++++
 gwm_test/Makefile     | 18 +++++++++
 gwm_test/gwm_test.cpp | 17 +++++++++
 4 files changed, 134 insertions(+)
 create mode 100644 gwm_integrators.cpp
 create mode 100644 gwm_integrators.h
 create mode 100644 gwm_test/Makefile
 create mode 100644 gwm_test/gwm_test.cpp

diff --git a/gwm_integrators.cpp b/gwm_integrators.cpp
new file mode 100644
index 0000000..3bae853
--- /dev/null
+++ b/gwm_integrators.cpp
@@ -0,0 +1,86 @@
+#include "catima/gwm_integrators.h"
+
+namespace catima {
+
+    double integrate_energyloss(Projectile& proj, const Material& mat, const Config& c)
+    {
+        static double s_estep_max = 0.001;
+        static int s_depth_max = 100;
+        int depth = 0;
+        double e_in = proj.T; // MeV/u
+        double e_final = e_in;
+        double x_step = 0.25*mat.thickness(); //g/cm^2
+        double x_traversed = 0.0;
+        double e_step = dedx(proj, mat, c)*x_step;
+
+        while(true)
+        {
+            if(e_step/e_final > s_estep_max && depth < s_depth_max)
+            {
+                ++depth;
+                x_step *= 0.5;
+                e_step = dedx(proj, mat, c)*x_step;
+            }
+            else if(x_step + x_traversed >= mat.thickness())
+            {
+                x_step = mat.thickness() - x_traversed;
+                e_step = dedx(proj, mat, c)*x_step;
+                e_final -= e_step;
+                proj.T = e_final;
+                return (e_in - e_final)*proj.A;
+            }
+            else if(depth == s_depth_max)
+            {
+                return e_in*proj.A;
+            }
+            else
+            {
+                e_step = dedx(proj, mat, c)*x_step;
+                e_final -= e_step;
+                proj.T = e_final;
+                x_traversed += x_step;
+            }
+        }
+    }
+
+    double reverse_integrate_energyloss(Projectile& proj, const Material& mat, const Config& c)
+    {
+        static double s_estep_max = 0.001;
+        static int s_depth_max = 100;
+        int depth = 0;
+        double e_out = proj.T; //MeV/u
+        double e_initial = e_out;
+        double x_step = 0.25*mat.thickness(); //g/cm^2
+        double x_traversed = 0.0;
+        double e_step = dedx(proj, mat, c)*x_step;
+        
+        while(true)
+        {
+            if(e_step/e_initial > s_estep_max && depth < s_depth_max)
+            {
+                ++depth;
+                x_step *= 0.5;
+                e_step = dedx(proj, mat, c)*x_step;
+            }
+            else if(x_step + x_traversed >= mat.thickness())
+            {
+                x_step = mat.thickness() -  x_traversed;
+                e_step = dedx(proj, mat, c)*x_step;
+                e_initial += e_step;
+                proj.T = e_initial;
+                return (e_initial - e_out)*proj.A;
+            }
+            else if(depth == s_depth_max)
+            {
+                return e_out*proj.A;
+            }
+            else
+            {
+                e_step = dedx(proj, mat, c)*x_step;
+                e_initial += e_step;
+                proj.T = e_initial;
+                x_traversed += x_step;
+            }
+        }
+    }
+}
\ No newline at end of file
diff --git a/gwm_integrators.h b/gwm_integrators.h
new file mode 100644
index 0000000..6ccf939
--- /dev/null
+++ b/gwm_integrators.h
@@ -0,0 +1,13 @@
+#ifndef GWM_INTEGRATORS_H
+#define GWM_INTEGRATORS_H
+
+#include "catima/catima.h"
+
+namespace catima {
+
+    double integrate_energyloss(Projectile& proj, const Material& mat, const Config& c=default_config);
+
+    double reverse_integrate_energyloss(Projectile& proj, const Material& mat, const Config& c=default_config);
+}
+
+#endif
\ No newline at end of file
diff --git a/gwm_test/Makefile b/gwm_test/Makefile
new file mode 100644
index 0000000..84350f4
--- /dev/null
+++ b/gwm_test/Makefile
@@ -0,0 +1,18 @@
+CC=g++
+EXE=gwm_test
+CXX_FLAGS= -std=c++17 -g -Wall
+
+CATIMA_PATH=../build/
+LIB_PATH=$(CATIMA_PATH)lib/
+INCLUDE_PATH=$(CATIMA_PATH)include/
+CPP=gwm_test.cpp
+
+.PHONY: all clean
+
+all: $(EXE)
+
+$(EXE): $(CPP)
+	$(CC) $(CXX_FLAGS) -I$(INCLUDE_PATH) $^ $(LIB_PATH)libcatima.a -o $@
+
+clean:
+	$(RM) $(EXE)
diff --git a/gwm_test/gwm_test.cpp b/gwm_test/gwm_test.cpp
new file mode 100644
index 0000000..29e78eb
--- /dev/null
+++ b/gwm_test/gwm_test.cpp
@@ -0,0 +1,17 @@
+#include "catima/gwm_integrators.h"
+#include "catima/nucdata.h"
+#include <iostream>
+
+int main()
+{
+    std::cout<<"-------Testing GWM Energy Loss Integration-------"<<std::endl;
+    catima::Projectile p1(catima::element_atomic_weight(1), 1.0, 0, 3.0);
+    catima::Material mat1(catima::get_material(6));
+    mat1.density(2.23).thickness(500.0*1e-6);
+
+    double result = catima::integrate_energyloss(p1, mat1);
+    std::cout<<"Energy loss (MeV): "<<result<<" Final energy: "<<p1.T<<std::endl;
+    result = catima::reverse_integrate_energyloss(p1, mat1);
+    std::cout<<"Reverse Energy loss (MeV): "<<result<<" Initial energy: "<<p1.T<<std::endl;
+
+}
\ No newline at end of file