From bc5a0faa47025547a179b61860089b8875b3e2c5 Mon Sep 17 00:00:00 2001
From: Ali Aghaeifar <ali.aghaeifar@tuebingen.mpg.de>
Date: Thu, 15 Feb 2024 12:53:53 +0100
Subject: [PATCH] reading gradient

---
 config/config_default.ini |  9 +++++
 src/file_utils.h          | 81 ++++++++++++++++++++++++---------------
 src/miscellaneous.h       | 26 +++++++++++--
 3 files changed, 82 insertions(+), 34 deletions(-)

diff --git a/config/config_default.ini b/config/config_default.ini
index af6c537..da8d0e1 100644
--- a/config/config_default.ini
+++ b/config/config_default.ini
@@ -61,6 +61,15 @@ DEPHASING[2] =
 DEPHASING_T[0] = 
 DEPHASING_T[1] =
 DEPHASING_T[2] =
+; Apply gradient in T/m for each axis. Gradients are active for one DWELL_TIME
+GRADIENT_XYZ[0] = 0.0 0.0 0.0
+GRADIENT_XYZ[1] = 0.0 0.0 0.0
+GRADIENT_XYZ[2] = 0.0 0.0 0.0
+; Time (sec) to apply gradient.
+GRADIENT_T[0] = 0.0
+GRADIENT_T[1] = 0.0
+GRADIENT_T[2] = 0.0
+
 
 DWELL_TIME  = 50e-6
 
diff --git a/src/file_utils.h b/src/file_utils.h
index 8ef55e6..e80ba33 100644
--- a/src/file_utils.h
+++ b/src/file_utils.h
@@ -162,16 +162,7 @@ bool read_config(std::string config_filename, simulation_parameters& param, std:
         for(i=0; i<MAX_RF && ini.get("SCAN_PARAMETERS").has("RF_ST[" + std::to_string(i) + "]"); i++)           
             param.RF_ST[i] = std::stof(ini.get("SCAN_PARAMETERS").get("RF_ST[" + std::to_string(i) + "]")) / param.dt;
         
-        // check RF start time conditions
-        if (std::is_sorted(param.RF_ST, param.RF_ST + i) == false || 
-            std::adjacent_find(param.RF_ST, param.RF_ST + i) != param.RF_ST + i || 
-            param.RF_ST[0] != 0 || 
-            (param.n_RF = i) == 0)
-        {
-            std::copy(param.RF_ST, param.RF_ST + param.n_RF, std::ostream_iterator<int>(std::cout, " ")); std::cout << std::endl;
-            std::cout << ERR_MSG << "RF Times must be in ascending order, starts with 0 and must not have duplicates values" << std::endl;
-            return false;
-        }
+        param.n_RF = i;
         // RF flip angles
         for(i=0; i<param.n_RF && ini.get("SCAN_PARAMETERS").has("RF_FA[" + std::to_string(i) + "]"); i++)
             param.RF_FA[i] = std::stof(ini.get("SCAN_PARAMETERS").get("RF_FA[" + std::to_string(i) + "]")) ;
@@ -191,31 +182,24 @@ bool read_config(std::string config_filename, simulation_parameters& param, std:
             return false;
         }
 
-        // ---------------- dephasing (start times, Flip angles ) ----------------
-        // Dephase start times
-        for(i=0; i<MAX_RF && ini.get("SCAN_PARAMETERS").has("DEPHASING_T[" + std::to_string(i) + "]"); i++)     
-        {
-            try    
-            {        
-                param.dephasing_T[i] = std::stof(ini.get("SCAN_PARAMETERS").get("DEPHASING_T[" + std::to_string(i) + "]")) / param.dt;    
-            }        
-            catch(const std::exception& e)
-            {
-                break; 
-            }                      
-        }      
-            
-        param.n_dephasing = i;
-
         // check RF start time conditions
-        if (std::is_sorted(param.dephasing_T, param.dephasing_T + i) == false || 
-            std::adjacent_find(param.dephasing_T, param.dephasing_T + i) != param.dephasing_T + i )
+        if (std::is_sorted(param.RF_ST, param.RF_ST + i) == false || 
+            std::adjacent_find(param.RF_ST, param.RF_ST + i) != param.RF_ST + i || 
+            param.RF_ST[0] != 0 || 
+            param.n_RF == 0)
         {
             std::copy(param.RF_ST, param.RF_ST + param.n_RF, std::ostream_iterator<int>(std::cout, " ")); std::cout << std::endl;
-            std::cout << ERR_MSG << "RF Times must be in ascending order and must not have duplicates values" << std::endl;
+            std::cout << ERR_MSG << "RF Times must be in ascending order, starts with 0 and must not have duplicates values" << std::endl;
             return false;
-        }        
+        }
 
+        // ---------------- dephasing (start times, Flip angles ) ----------------
+        // Dephase start times
+        for(i=0; i<MAX_RF && ini.get("SCAN_PARAMETERS").has("DEPHASING_T[" + std::to_string(i) + "]"); i++)       
+            param.dephasing_T[i] = std::stof(ini.get("SCAN_PARAMETERS").get("DEPHASING_T[" + std::to_string(i) + "]")) / param.dt;    
+ 
+        param.n_dephasing = i;   
+        // Dephase flip angles
         for(i=0; i<param.n_dephasing && ini.get("SCAN_PARAMETERS").has("DEPHASING[" + std::to_string(i) + "]"); i++)
             param.dephasing[i] = std::stof(ini.get("SCAN_PARAMETERS").get("DEPHASING[" + std::to_string(i) + "]")) ;
         
@@ -223,7 +207,42 @@ bool read_config(std::string config_filename, simulation_parameters& param, std:
         {
             std::cout << ERR_MSG << "DEPHASING and DEPHASING_T must have the same number of elements" << std::endl;
             return false;
-        }
+        } 
+
+        // check Dephase start time conditions
+        if (std::is_sorted(param.dephasing_T, param.dephasing_T + i) == false || 
+            std::adjacent_find(param.dephasing_T, param.dephasing_T + i) != param.dephasing_T + i )
+        {
+            std::copy(param.dephasing_T, param.dephasing_T + param.n_dephasing, std::ostream_iterator<int>(std::cout, " ")); std::cout << std::endl;
+            std::cout << ERR_MSG << "dephasing Times must be in ascending order and must not have duplicates values" << std::endl;
+            return false;
+        } 
+            
+        // ---------------- Gradients (start times, strength (T/m) ) ----------------
+        // Gradient start times
+        for(i=0; i<MAX_GRADIENT && ini.get("SCAN_PARAMETERS").has("GRADIENT_T[" + std::to_string(i) + "]"); i++)        
+            param.gradient_T[i] = std::stof(ini.get("SCAN_PARAMETERS").get("GRADIENT_T[" + std::to_string(i) + "]")) / param.dt;    
+
+        param.n_gradient = i;
+        // Gradient strength
+        for(i=0; i<param.n_gradient && ini.get("SCAN_PARAMETERS").has("GRADIENT_XYZ[" + std::to_string(i) + "]"); i++)
+            param.gradient_xyz[i] = std::stof(ini.get("SCAN_PARAMETERS").get("GRADIENT_XYZ[" + std::to_string(i) + "]")) ;
+               
+        if(i != param.n_gradient)
+        {
+            std::cout << ERR_MSG << "GRADIENT_XYZ and GRADIENT_T must have the same number of elements" << std::endl;
+            return false;
+        } 
+
+        // check Dephase start time conditions
+        if (std::is_sorted(param.gradient_T, param.gradient_T + i) == false || 
+            std::adjacent_find(param.gradient_T, param.gradient_T + i) != param.gradient_T + i )
+        {
+            std::copy(param.gradient_T, param.gradient_T + param.n_gradient, std::ostream_iterator<int>(std::cout, " ")); std::cout << std::endl;
+            std::cout << ERR_MSG << "Gradient Times must be in ascending order and must not have duplicates values" << std::endl;
+            return false;
+        } 
+        
     }
 
     // ============== reading section SCAN_PARAMETERS ==============
diff --git a/src/miscellaneous.h b/src/miscellaneous.h
index fdf7fd0..c824a36 100644
--- a/src/miscellaneous.h
+++ b/src/miscellaneous.h
@@ -25,20 +25,36 @@
 #define MAX_RF 128          // maximum number of RF
 #define MAX_TE 256          // maximum number of echo times
 #define MAX_DEPHASE 128     // maximum number of dephasing
+#define MAX_GRADIENT 128    // maximum number of gradient
 
 typedef struct simulation_parameters
 {
     float T1, T2, TR, dt, B0, e1, e12, e2, e22, c, s, c2, s2;
     float RF_FA[MAX_RF], RF_PH[MAX_RF]; // refocusing FA
     float dephasing[MAX_DEPHASE]; // dephasing in degree
-    int32_t RF_ST[MAX_RF], TE[MAX_TE], dephasing_T[MAX_DEPHASE]; // refocusing time in dt, echo times in dt, dephasing time in dt
+    float gradient_xyz[3*MAX_GRADIENT]; // gradient in T/m
+    int32_t RF_ST[MAX_RF], TE[MAX_TE], dephasing_T[MAX_DEPHASE], gradient_T[MAX_GRADIENT]; // refocusing time in dt, echo times in dt, dephasing time in dt
     float sample_length[3], scale2grid[3], diffusion_const, phase_cycling;
-    int32_t n_timepoints, n_sample_length_scales, n_fieldmaps, n_TE, n_RF, n_dephasing;
+    int32_t n_timepoints, n_sample_length_scales, n_fieldmaps, n_TE, n_RF, n_dephasing, n_gradient;
     int32_t n_dummy_scan;
     uint32_t n_spins, fieldmap_size[3], seed;
     uint64_t matrix_length;
     bool enDebug, enApplyFA2, enCrossBoundry;
-    simulation_parameters():T1(2.2),T2(0.04),TR(0.04),dt(5e-5),B0(9.4),n_TE(0),n_RF(0),n_dummy_scan(0),phase_cycling(0.),enApplyFA2(false),enDebug(false),enCrossBoundry(true)
+    simulation_parameters():
+        T1(2.2),
+        T2(0.04),
+        TR(0.04),
+        dt(5e-5),
+        B0(9.4),
+        n_TE(0),
+        n_RF(0),
+        n_dephasing(0),
+        n_gradient(0),
+        n_dummy_scan(0),
+        phase_cycling(0.),
+        enApplyFA2(false),
+        enDebug(false),
+        enCrossBoundry(true)
     {
         memset(fieldmap_size, 0, 3*sizeof(fieldmap_size[0])); 
         memset(sample_length, 0, 3*sizeof(sample_length[0]));
@@ -48,6 +64,8 @@ typedef struct simulation_parameters
         memset(RF_PH, 0, MAX_RF*sizeof(RF_PH[0]));
         memset(dephasing, 0, MAX_DEPHASE*sizeof(dephasing[0]));
         memset(dephasing_T, 0, MAX_DEPHASE*sizeof(dephasing_T[0]));
+        memset(gradient_xyz, 0, 3*MAX_GRADIENT*sizeof(gradient_xyz[0]));
+        memset(gradient_T, 0, MAX_GRADIENT*sizeof(gradient_T[0]));
     }
 
     void dump()
@@ -59,6 +77,8 @@ typedef struct simulation_parameters
         std::cout<<"RF time = "; for(int i=0; i<n_RF; i++) std::cout<<RF_ST[i]*dt<<' '; std::cout<<'\n';
         std::cout<<"dephasing degree = "; for(int i=0; i<n_dephasing; i++) std::cout<<dephasing[i]<<' '; std::cout<<'\n';
         std::cout<<"dephasing time = "; for(int i=0; i<n_dephasing; i++) std::cout<<dephasing_T[i]*dt<<' '; std::cout<<'\n';
+        std::cout<<"gradient T/m = "; for(int i=0; i<3*n_gradient; i++) std::cout<<gradient_xyz[i]<<' '; std::cout<<'\n';
+        std::cout<<"gradient time = "; for(int i=0; i<n_gradient; i++) std::cout<<gradient_T[i]<<' '; std::cout<<'\n';
         std::cout<<"sample length = "<< sample_length[0] << " x " << sample_length[1] << " x " << sample_length[2] << " m" << '\n';
         std::cout<<"scale2grid = "<< scale2grid[0] << " x " << scale2grid[1] << " x " << scale2grid[2] << '\n';
         std::cout<<"fieldmap size = "<< fieldmap_size[0] << " x " << fieldmap_size[1] << " x " << fieldmap_size[2] << '\n';