Oriented sesans transform #536

pkienzle · 2022-11-04T22:35:51Z

Allow modelling of SESANS for oriented samples.

This is an old set of changes that was never submitted for review.

I'm opening this PR as a draft so that SESANS developers can discuss it.

Note: there is a corresponding branch in the sasview repository. https://github.com/SasView/sasview/tree/costrafo411

smk78 · 2022-11-07T10:12:42Z

I will bring this to the attention of our SESANS team.

I note that the allied SasView branch is >4000 commits behand main...

rmdalgliesh · 2022-11-07T12:18:15Z

I suspect that this is the code that should have been worked on in the last code camp. The sesanscosine.py code contains a lot of modifications that make a lot more sense for TOF-SESANS. They answer quite a lot of the questions I had regarding how the sesans.py code was working for TOF.
I do have some concerns. The call to HankelAccept contains a reference to magfield in order to define the spin-echo length. This doesn't make sense in TOF. If this is required then I would suggest the calculation for SEL should be
SEL= SELC*wavelength^2 where SELC is a constant defined by the parameters of the beamline.

My other concern is the need for a reference for the cosine transform.
The hankel transform is well established and references but I'm not aware of anything where the cosine transform for anisotropic scattering is specifically discussed.

So... my comment is that I think this work and the work from the code camp need to be merged to remove the bugs that were ironed out but it is probably best if those patches are made and then we merge in these changes.

Rob

lucas-wilkins · 2022-11-07T13:35:21Z

Note, I've opened the corresponding PR in sasview as a draft.

lucas-wilkins

Definitely a few things here I'm unsure about, or that could be done better.

lucas-wilkins · 2022-11-07T14:06:53Z

sasmodels/sesanscos.py

+
+        repq = np.tile(q, (SElength.size, 1)).T
+        repSE = np.tile(SElength, (q.size, 1))
+        cosmat = np.float32(dq / (2 * pi)) * np.cos(repSE * repq)


I'm pretty sure an out of the box discrete cosine transform will have time complexity of $N \log N$ rather than the $N^2$ complexity of the transform here.

lucas-wilkins · 2022-11-07T14:12:38Z

sasmodels/sesanscosine.py

+
+    P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0]))
+
+def hankel(SElength, wavelength, thickness, q, Iq):


Also, the Hankel transform can be done in $N \log N$ time too.

lucas-wilkins · 2022-11-07T14:14:44Z

sasmodels/sesanscosine.py

+#     HankelTransform with fixed circular acceptance angle (circular aperture) for Time of Flight SESANS
+#==============================================================================
+#acceptq is the q-space needed to create limited acceptance effect
+    SElength= wavelength*magfield


@rmdalgliesh Mentioned that there are some difficulties with specifying things in terms of the magnetic field, as this is not as general as it can be.

lucas-wilkins · 2022-11-07T14:23:15Z

sasmodels/sesanscosine.py

+        # data.has_z_acceptance
+        return [q]
+
+def transform(data, q_calc, Iq_calc, qmono, Iq_mono):


Probably want to keep this under user control to some extent.

lucas-wilkins · 2022-11-07T14:23:55Z

sasmodels/sesanscosine.py

+Wim Bouwman ([email protected]), June 2013
+"""
+
+from __future__ import division


The future is here!

lucas-wilkins · 2022-11-07T14:24:39Z

sasmodels/sesanscosine.py

+from __future__ import division
+
+import numpy as np  # type: ignore
+from numpy import pi, exp  # type: ignore


Minor: should really use np.pi and np.exp

lucas-wilkins · 2022-11-07T14:25:20Z

sasmodels/sesanscosine.py

+
+    from sas.sascalc.data_util.nxsunit import Converter
+    wavelength = Converter(wavelength[1])(wavelength[0],"A")
+    thickness = Converter(thickness[1])(thickness[0],"A")


Are these up to date with the latest unit conversion stuff?

lucas-wilkins · 2022-11-07T14:27:56Z

sasmodels/sesanscosine.py

+    If the instrument has a circular acceptance, 1 all_q vector is needed
+
+    """
+    if not data.has_no_finite_acceptance:


is .has_no_finite_acceptance an attribute of all data classes - and should it be?

lucas-wilkins · 2022-11-07T14:30:57Z

sasmodels/sesanscosine.py

+                  q_calc, Iq_calc)
+
+def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono):
+    return hankel(data.x, data.lam * 1e-9,


Is calling hankel here correct?

lucas-wilkins · 2022-11-07T14:31:14Z

sasmodels/sesanscosine.py

+                  q_calc, Iq_calc)
+
+def call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono):
+    return hankel(data.x, data.y, data.lam * 1e-9,


Is calling hankel here correct, either?

smk78 · 2022-11-07T16:29:35Z

@rmdalgliesh commented:

My other concern is the need for a reference for the cosine transform. The hankel transform is well established and references but I'm not aware of anything where the cosine transform for anisotropic scattering is specifically discussed.

@wimbouwman replies by email:

The cosine transform is one of the first SESANS articles (1996) by Theo Rekveldt (https://doi.org/10.1016/0168-583X(96)00213-3). The numerical drawback of the cosine transform is that it is 2D. Only later we found out that it can be written as the Hankeltransform (2008, https://dx.doi.org/10.1107/S0021889808026770), which is only 1D.

butlerpd · 2023-11-21T23:56:44Z

@caitwolf , @dehoni @rmdalgliesh @gnsmith @wimbouwman - Should this be a priority for 6.0?

dehoni · 2023-11-25T15:46:04Z

sasmodels/sesanscosine.py

+        allIq=np.linspace()  # This is the model applied to the allq q-space. Needs to refference the model somehow
+        alldq = (allq[1]-allq[0])*1e10
+        sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq)
+        s[i]=1-exp(-sigma)


sigma[i] and s[i] not used further, would bring the code in line with the other calls to calculate the correlation function G.

dehoni · 2023-11-25T15:50:56Z

sasmodels/sesanscosine.py

+
+    G *= dq*1e10*2*pi
+
+    P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0]))


One would expect to get corrected data in reduced data (spin-echolength, correlation function G). Separate the conversion to polarisation into separate function to avoid code repetition. I would assume that the data format will develop with time and instrument specific information (wavelength, sample thickness, magnetic field) is retained in the data reduction.

wimbouwman · 2023-11-27T08:01:37Z

Hi Dirk, For your reference, it started even with the cosine transform: Rekveldt, M. T. (1996). Novel SANS instrument using neutron spin echo. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms, 114(3-4), 366-370. Only later with Timofei Kruglov we saw that the integral could be made 1D with the Hankel transform. Wim From: Dirk Honecker ***@***.***> Sent: zaterdag 25 november 2023 16:06 To: SasView/sasmodels ***@***.***> Cc: Wim Bouwman ***@***.***>; Mention ***@***.***> Subject: Re: [SasView/sasmodels] Oriented sesans transform (PR #536) I suspect that this is the code that should have been worked on in the last code camp. The sesanscosine.py code contains a lot of modifications that make a lot more sense for TOF-SESANS. They answer quite a lot of the questions I had regarding how the sesans.py code was working for TOF. I do have some concerns. The call to HankelAccept contains a reference to magfield in order to define the spin-echo length. This doesn't make sense in TOF. If this is required then I would suggest the calculation for SEL should be SEL= SELC*wavelength^2 where SELC is a constant defined by the parameters of the beamline. My other concern is the need for a reference for the cosine transform. The hankel transform is well established and references but I'm not aware of anything where the cosine transform for anisotropic scattering is specifically discussed. So... my comment is that I think this work and the work from the code camp need to be merged to remove the bugs that were ironed out but it is probably best if those patches are made and then we merge in these changes. Rob Only for spherically symmetric (i.e. isotropic) density distribution you will have the Fourier–Abel–Hankel. For an oriented system, one has to consider the projection of the scattering length density correlation function along the encoding direction J. Krouglov et al., Appl. Cryst. (2003). 36, 1417. — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https:/github.com/SasView/sasmodels/pull/536*issuecomment-1826353544__;Iw!!PAKc-5URQlI!4BvOZm5ggtDppH70C-Gr6pgdyUYJ-WCKtMZNcsNeKTgJA_KXLM6JS0V0hmBMfvX4fl52y4IL7zeAx5LdRAAt8Hd_rgJ1QQ$>, or unsubscribe<https://urldefense.com/v3/__https:/github.com/notifications/unsubscribe-auth/ACUIQEJ5KTGFEIKZKZJBPILYGICOPAVCNFSM6AAAAAARXTNHMOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQMRWGM2TGNJUGQ__;!!PAKc-5URQlI!4BvOZm5ggtDppH70C-Gr6pgdyUYJ-WCKtMZNcsNeKTgJA_KXLM6JS0V0hmBMfvX4fl52y4IL7zeAx5LdRAAt8HcajV2EIg$>. You are receiving this because you were mentioned.Message ID: ***@***.******@***.***>>

wimbouwman · 2023-11-27T08:05:59Z

Hi Dirk, I am indeed surprised by this polarisation calculation. I didn’t think we originally had a polarisation calculation to avoid instrument specific parameters in SasView. If anybody would really feel a strong need for it, then it should not be in the cosine transform to avoid code duplication. Wim From: Dirk Honecker ***@***.***> Sent: zaterdag 25 november 2023 16:51 To: SasView/sasmodels ***@***.***> Cc: Wim Bouwman ***@***.***>; Mention ***@***.***> Subject: Re: [SasView/sasmodels] Oriented sesans transform (PR #536) @dehoni commented on this pull request.

________________________________ In sasmodels/sesanscosine.py<https://urldefense.com/v3/__https:/github.com/SasView/sasmodels/pull/536*discussion_r1405146543__;Iw!!PAKc-5URQlI!7ijUgtX8YouLiMuTclrcWzts8b_GlkltLETyqGifuIHqV8lH1ufhYehSCJ6-P9CyGt01ilyt5YZ583udKVtYQCngJtcuaQ$>:

+ alldq = (allq[1]-allq[0])*1e10

+ sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) + s[i]=1-exp(-sigma) + + dq = (q[1]-q[0])*1e10 + a = (x<threshold) + acceptq = a*q + acceptIq = a*Iq + + G[i] = np.sum(besselj(0, acceptq*SElength_i)*acceptIq*acceptq*dq) + +# G[i]=np.sum(integral) + + G *= dq*1e10*2*pi + + P = exp(thickness*wavelength**2/(4*pi**2)*(G-G[0])) One would expect to get corrected data in reduced data (spin-echolength, correlation function G). Separate the conversion to polarisation into separate function to avoid code repetition. I would assume that the data format will develop with time and instrument specific information (wavelength, sample thickness, magnetic field) is retained in the data reduction. — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https:/github.com/SasView/sasmodels/pull/536*pullrequestreview-1749205098__;Iw!!PAKc-5URQlI!7ijUgtX8YouLiMuTclrcWzts8b_GlkltLETyqGifuIHqV8lH1ufhYehSCJ6-P9CyGt01ilyt5YZ583udKVtYQClasUAUOQ$>, or unsubscribe<https://urldefense.com/v3/__https:/github.com/notifications/unsubscribe-auth/ACUIQEKMXWZRFJKTAY3K6RLYGIHWZAVCNFSM6AAAAAARXTNHMOVHI2DSMVQWIX3LMV43YUDVNRWFEZLROVSXG5CSMV3GSZLXHMYTONBZGIYDKMBZHA__;!!PAKc-5URQlI!7ijUgtX8YouLiMuTclrcWzts8b_GlkltLETyqGifuIHqV8lH1ufhYehSCJ6-P9CyGt01ilyt5YZ583udKVtYQCmlrDpJig$>. You are receiving this because you were mentioned.Message ID: ***@***.******@***.***>>

jhbakker and others added 7 commits April 8, 2017 16:13

include 2D cosine transform in sesans.py

e4e5e29

isgood

a86e249

oriented sesans: split transform into oriented and unoriented versions

7c1cce3

provide working oriented/unoriented sesans examples

2cdc35b

Merge branch 'master' into costrafo411

2f9f1ec

Merge branch 'master' into costrafo411

65ceb7d

Merge branch 'master' into costrafo411

3a45c2c

pkienzle marked this pull request as draft November 4, 2022 22:37

pkienzle requested review from rprospero and wimbouwman November 4, 2022 22:38

smk78 requested review from dehoni and rmdalgliesh November 7, 2022 10:18

caitwolf self-requested a review November 7, 2022 14:06

lucas-wilkins requested changes Nov 7, 2022

View reviewed changes

caitwolf mentioned this pull request Nov 21, 2023

Fixed inverse q-spacing in the geometric extrapolation function #571

Merged

dehoni reviewed Nov 25, 2023

View reviewed changes

wpotrzebowski added DiscussAtTheCall and removed DiscussAtTheCall labels Dec 4, 2023

wpotrzebowski requested a review from gnsmith December 5, 2023 14:24

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Oriented sesans transform #536

Oriented sesans transform #536

pkienzle commented Nov 4, 2022 •

edited by smk78

Loading

smk78 commented Nov 7, 2022 •

edited

Loading

rmdalgliesh commented Nov 7, 2022

lucas-wilkins commented Nov 7, 2022

lucas-wilkins left a comment

lucas-wilkins Nov 7, 2022

lucas-wilkins Nov 7, 2022

lucas-wilkins Nov 7, 2022

lucas-wilkins Nov 7, 2022

lucas-wilkins Nov 7, 2022

lucas-wilkins Nov 7, 2022

lucas-wilkins Nov 7, 2022

lucas-wilkins Nov 7, 2022

lucas-wilkins Nov 7, 2022

lucas-wilkins Nov 7, 2022

smk78 commented Nov 7, 2022 •

edited

Loading

butlerpd commented Nov 21, 2023

dehoni Nov 25, 2023

dehoni Nov 25, 2023

wimbouwman commented Nov 27, 2023 via email

wimbouwman commented Nov 27, 2023 via email


		P = exp(thicknesswavelength2/(4pi*2)(G-G[0]))

		def hankel(SElength, wavelength, thickness, q, Iq):


		G = dq1e102pi

		P = exp(thicknesswavelength2/(4pi*2)(G-G[0]))

Oriented sesans transform #536

Are you sure you want to change the base?

Oriented sesans transform #536

Conversation

pkienzle commented Nov 4, 2022 • edited by smk78 Loading

smk78 commented Nov 7, 2022 • edited Loading

rmdalgliesh commented Nov 7, 2022

lucas-wilkins commented Nov 7, 2022

lucas-wilkins left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

smk78 commented Nov 7, 2022 • edited Loading

butlerpd commented Nov 21, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

wimbouwman commented Nov 27, 2023 via email

wimbouwman commented Nov 27, 2023 via email

pkienzle commented Nov 4, 2022 •

edited by smk78

Loading

smk78 commented Nov 7, 2022 •

edited

Loading

smk78 commented Nov 7, 2022 •

edited

Loading