Merge pull request #191 from shirokobrod/master

Address a case where energy E is less than the lowest value in the energies array of Kissel's cross section Fixes #187
tschoonj · Feb 8, 2022 · 9e8ab95 · 9e8ab95
2 parents 1fd9e5c + ed31560
commit 9e8ab95
Show file tree

Hide file tree

Showing 4 changed files with 10 additions and 22 deletions.
diff --git a/java/Xraylib.java b/java/Xraylib.java
@@ -935,7 +935,10 @@ public static double CSb_Photo_Partial(int Z, int shell, double E) {
     }
 
     ln_E = Math.log(E);
-    if (EdgeEnergy_Kissel_arr[Z * SHELLNUM_K + shell] > EdgeEnergy_arr[Z * SHELLNUM + shell] && E < EdgeEnergy_Kissel_arr[Z * SHELLNUM_K + shell]) {
+    if (ln_E < E_Photo_Partial_Kissel_arr[Z][shell][0]) {
+	  /* Address a case where energy E is less than the lowest value in the energies array of Kissel's cross section
+         Fixes https://github.com/tschoonj/xraylib/issues/187 
+      */
       /*
        * use log-log extrapolation
        */
@@ -953,15 +956,6 @@ else if (m < -1.0)
         m=-1.0;
       ln_sigma = y0+m*(ln_E-x0);
     }
-    else if (ln_E < E_Photo_Partial_Kissel_arr[Z][shell][0]) {
-      /* Address edge case where Kissel edge energy is less than the lowest value in the energies array
-          Fixes https://github.com/tschoonj/xraylib/issues/187 
-          A better fix would involve setting the Kissel edge energies to the lowest value present in the energies array,
-          but this needs to be done in the script that extracts the data from the Kissel files.
-          This script is currently written in IDL, and I am in no mood to translate to Python right now
-      */
-      ln_sigma = Photo_Partial_Kissel_arr[Z][shell][0];
-    }
     else {
       ln_sigma = splint(E_Photo_Partial_Kissel_arr[Z][shell], Photo_Partial_Kissel_arr[Z][shell], Photo_Partial_Kissel_arr2[Z][shell],NE_Photo_Partial_Kissel_arr[Z][shell], ln_E);
     }

diff --git a/java/tests/TestKisselPE.java b/java/tests/TestKisselPE.java
@@ -108,7 +108,7 @@ public void test_cs_photo_partial_good() {
 		assertEquals(cs, 0.0024892741933504824, 1E-6);
 		/* see https://github.com/tschoonj/xraylib/issues/187 */
 		cs = Xraylib.CSb_Photo_Partial(47, Xraylib.L2_SHELL, 3.5282);
-		assertEquals(cs, 1.56958E+04, 1E-2);
+		assertEquals(cs, 1.569549E+04, 1E-2);
 	}
 
 	static Stream<Arguments> badCSPhotoPartialsValuesProvider() {

diff --git a/src/kissel_pe.c b/src/kissel_pe.c
@@ -132,7 +132,10 @@ double CSb_Photo_Partial(int Z, int shell, double E, xrl_error **error) {
   } 
 
   ln_E = log(E);
-  if (EdgeEnergy_Kissel[Z][shell] > EdgeEnergy_arr[Z][shell] && E < EdgeEnergy_Kissel[Z][shell]) {
+  if(ln_E < E_Photo_Partial_Kissel[Z][shell][0]) {
+  	/* Address a case where energy E is less than the lowest value in the energies array of Kissel's cross section
+       Fixes https://github.com/tschoonj/xraylib/issues/187 
+    */
     /*
      * use log-log extrapolation 
      */
@@ -150,15 +153,6 @@ double CSb_Photo_Partial(int Z, int shell, double E, xrl_error **error) {
       m = -1.0;
     ln_sigma = y0 + m * (ln_E - x0);
   }
-  else if (ln_E < E_Photo_Partial_Kissel[Z][shell][0]) {
-    /* Address edge case where Kissel edge energy is less than the lowest value in the energies array
-       Fixes https://github.com/tschoonj/xraylib/issues/187 
-       A better fix would involve setting the Kissel edge energies to the lowest value present in the energies array,
-       but this needs to be done in the script that extracts the data from the Kissel files.
-       This script is currently written in IDL, and I am in no mood to translate to Python right now
-    */
-    ln_sigma = Photo_Partial_Kissel[Z][shell][0];
-  }
   else {
     int splint_rv = splint(E_Photo_Partial_Kissel[Z][shell] - 1, Photo_Partial_Kissel[Z][shell] - 1, Photo_Partial_Kissel2[Z][shell] - 1, NE_Photo_Partial_Kissel[Z][shell], ln_E, &ln_sigma, error);
     if (!splint_rv)

diff --git a/tests/test-kissel_pe.c b/tests/test-kissel_pe.c
@@ -200,7 +200,7 @@ int main(int argc, char **argv) {
 	/* see https://github.com/tschoonj/xraylib/issues/187 */
 	cs = CSb_Photo_Partial(47, L2_SHELL, 3.5282, &error);
 	assert(error == NULL);
-	assert(fabs(cs - 1.56958E+04)/cs < 1E-6);
+	assert(fabs(cs - 1.569549E+04)/cs < 1E-6);
 
 	cs = CS_Photo_Partial(26, K_SHELL, 6.0, &error);
 	assert(error != NULL);