Repository: tschoonj/xraylib (language: python)

Issue:
CS_FluorLine_Kissel fails on some energies
Hi!

I've been using xraylib a lot and I really appreciate it very much! Thank you so much for continuing the work on it, Tom.

Only recently it happened that I came across this problem. So out of millions of calls only one failed so far! So it definitely does not happen often, but nevertheless it causes of a crash, so I thought better to report it in case there is a nice fix for it. So if you call the CS_FluorLine_Kissel with certain energy near the edge Ag L2, the function call fails. Just a little bit different energy will work. I did not test systematically all elements and lines and energies down to fifth decimal, but that might be a fun exercise to run (and would take some calculation time!).

xrl.CS_FluorLine_Kissel(47, -63, 3.5282)
Traceback (most recent call last):
File ..., in <module>xrl.CS_FluorLine_Kissel(47, -63, 3.5282)
File "C:\ProgramData\Anaconda3\lib\site-packages\xraylib.py", line 4424, in CS_FluorLine_Kissel
return _xraylib.CS_FluorLine_Kissel(Z, line, E)
ValueError: Spline extrapolation is not allowed

Kind regards,

Ulla
CS_FluorLine_Kissel fails on some energies
Hi!

I've been using xraylib a lot and I really appreciate it very much! Thank you so much for continuing the work on it, Tom.

Only recently it happened that I came across this problem. So out of millions of calls only one failed so far! So it definitely does not happen often, but nevertheless it causes of a crash, so I thought better to report it in case there is a nice fix for it. So if you call the CS_FluorLine_Kissel with certain energy near the edge Ag L2, the function call fails. Just a little bit different energy will work. I did not test systematically all elements and lines and energies down to fifth decimal, but that might be a fun exercise to run (and would take some calculation time!).

xrl.CS_FluorLine_Kissel(47, -63, 3.5282)
Traceback (most recent call last):
File ..., in <module>xrl.CS_FluorLine_Kissel(47, -63, 3.5282)
File "C:\ProgramData\Anaconda3\lib\site-packages\xraylib.py", line 4424, in CS_FluorLine_Kissel
return _xraylib.CS_FluorLine_Kissel(Z, line, E)
ValueError: Spline extrapolation is not allowed

Kind regards,

Ulla


Discussion / hints:
Hi Ulla, this looks like a bug indeed. I will try and come up with a fix, but it may take several weeks before I have the time to do so.

Best, Tom
Hi Ulla, this looks like a bug indeed. I will try and come up with a fix, but it may take several weeks before I have the time to do so.

Best, Tom

Fix patch:
diff --git a/java/Xraylib.java b/java/Xraylib.java
index a36f3a29..1c21a5f4 100644
--- 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/src/kissel_pe.c b/src/kissel_pe.c
index 1e8932a2..844bb1bc 100644
--- 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)


Test patch:
diff --git a/java/tests/TestKisselPE.java b/java/tests/TestKisselPE.java
index 7ac683d9..66371e54 100644
--- 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/tests/test-kissel_pe.c b/tests/test-kissel_pe.c
index 77a8c85f..775a9f6f 100644
--- 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,8 +190,6 @@ 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 */
@@ -195,6 +197,11 @@ int main(int argc, char **argv) {
 	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);