Merge pull request #189 from tschoonj/fix-issue-187

Fix issue 187
tschoonj · Feb 1, 2022 · 1fd9e5c · 1fd9e5c
2 parents 3fe41fb + d1b6403
commit 1fd9e5c
Show file tree

Hide file tree

Showing 4 changed files with 30 additions and 2 deletions.
diff --git a/java/Xraylib.java b/java/Xraylib.java
@@ -953,6 +953,15 @@ 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
@@ -106,6 +106,9 @@ public void test_cs_photo_partial_good() {
 		assertEquals(cs, 22.40452459077649, 1E-6);
 		cs = Xraylib.CS_Photo_Partial(26, Xraylib.K_SHELL, 300.0);
 		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);
 	}
 
 	static Stream<Arguments> badCSPhotoPartialsValuesProvider() {

diff --git a/src/kissel_pe.c b/src/kissel_pe.c
@@ -150,6 +150,15 @@ 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
@@ -92,6 +92,10 @@ int main(int argc, char **argv) {
 	assert(error == NULL);
 	assert(fabs(0.102639909656483 - cs) < 1E-6);
 
+	/* see https://github.com/tschoonj/xraylib/issues/187 */
+	cs = CS_FluorLine_Kissel(47, L2M4_LINE, 3.5282, &error);
+	assert(error == NULL);
+
 	/* lets try some bad input */
 	cs = CS_FluorLine_Kissel(0, KL3_LINE, 10.0, &error);
 	assert(error != NULL);
@@ -186,15 +190,18 @@ int main(int argc, char **argv) {
 	cs += CS_FluorLine_Kissel(92, L3N7_LINE, 30.0, &error);
 	assert(error == NULL);
 	cs2 = CS_FluorLine_Kissel(92, LB_LINE, 30.0, NULL);
-	fprintf(stderr, "cs: %f\n", cs);
-	fprintf(stderr, "cs2: %f\n", cs2);
 	assert(fabs(cs2 - cs) < 1E-6);
 
 	/* CS_Photo_Partial tests */
 	cs = CS_Photo_Partial(26, K_SHELL, 20.0, &error);
 	assert(error == NULL);
 	assert(fabs(cs - 22.40452459077649) < 1E-6);
 
+	/* 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);
+
 	cs = CS_Photo_Partial(26, K_SHELL, 6.0, &error);
 	assert(error != NULL);
 	assert(cs == 0.0);