From e4e5e2926535e9bc7a2f21be024b517e10a388e5 Mon Sep 17 00:00:00 2001 From: jhbakker Date: Sat, 8 Apr 2017 16:13:19 +0200 Subject: [PATCH 1/4] include 2D cosine transform in sesans.py --- sasmodels/sesans.py | 40 ++++++++++++++++++++++++++++++++++++---- 1 file changed, 36 insertions(+), 4 deletions(-) diff --git a/sasmodels/sesans.py b/sasmodels/sesans.py index b7aabefab..9b8b482d0 100644 --- a/sasmodels/sesans.py +++ b/sasmodels/sesans.py @@ -44,13 +44,23 @@ def __init__(self, z, SElength, zaccept, Rmax): # type: (np.ndarray, float, float) -> None #import logging; logging.info("creating SESANS transform") self.q = z + # isoriented flag determines whether data is from an oriented sample or not, should be in the data or selectable in GUI. + #if self.isoriented==False self._set_hankel(SElength, zaccept, Rmax) + #if self.isoriented==True + self._set_cosmat(SElength, zaccept, Rmax) def apply(self, Iq): - # tye: (np.ndarray) -> np.ndarray - G0 = np.dot(self._H0, Iq) - G = np.dot(self._H.T, Iq) - P = G - G0 + if len(Iq.size) == 1: # if isotropic, do Hankel transform + # type: (np.ndarray) -> np.ndarray + G0 = np.dot(self._H0, Iq) + G = np.dot(self._H.T, Iq) + P = G - G0 + elif len(Iq.size) == 2: + dq=self.q_calc[0] + G0 = sum(np.dot(self._cos0, Iq)*dq) + G = sum(np.dot(self._cosmat.T, Iq)*dq) + P = G - G0 return P def _set_hankel(self, SElength, zaccept, Rmax): @@ -72,3 +82,25 @@ def _set_hankel(self, SElength, zaccept, Rmax): self.q_calc = q self._H, self._H0 = H, H0 + + def _set_cosmat(self, SElength, zaccept, Rmax): + # type: (np.ndarray, float, float) -> None + # Force float32 arrays, otherwise run into memory problems on some machines + SElength = np.asarray(SElength, dtype='float32') + #qymax and qzmax depend on detector shape + qzmax= + # Rmax = #value in text box somewhere in FitPage? + q_max = 2 * pi / (SElength[1] - SElength[0]) + q_min = 0.1 * 2 * pi / (np.size(SElength) * SElength[-1]) + + q = np.arange(q_min, q_max, q_min, dtype='float32') + dq = q_min + + cos0 = np.float32(dq / (2 * pi)) + + 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) + + self.q_calc = q + self._cosmat, self._cos0 = cosmat, cos0 \ No newline at end of file From a86e249e1edc9902d8497f521e4085aa171ee02e Mon Sep 17 00:00:00 2001 From: jhbakker Date: Sat, 8 Apr 2017 17:11:29 +0200 Subject: [PATCH 2/4] isgood --- example/6pc SD.txt | 65 +++++++++++ example/Beise009798_01.ses | 91 ++++++++++++++++ example/Beise009798_02.ses | 91 ++++++++++++++++ example/Beise009798_03.ses | 91 ++++++++++++++++ example/Gab_5pc_deb2p0.ses | 36 +++++++ example/Gab_5pc_deb2p0.txt | 36 +++++++ sasmodels/sesans.py | 29 ++--- sasmodels/sesanscos.py | 106 ++++++++++++++++++ sasmodels/sesanscosine.py | 214 +++++++++++++++++++++++++++++++++++++ 9 files changed, 745 insertions(+), 14 deletions(-) create mode 100644 example/6pc SD.txt create mode 100644 example/Beise009798_01.ses create mode 100644 example/Beise009798_02.ses create mode 100644 example/Beise009798_03.ses create mode 100644 example/Gab_5pc_deb2p0.ses create mode 100644 example/Gab_5pc_deb2p0.txt create mode 100644 sasmodels/sesanscos.py create mode 100644 sasmodels/sesanscosine.py diff --git a/example/6pc SD.txt b/example/6pc SD.txt new file mode 100644 index 000000000..956d7691a --- /dev/null +++ b/example/6pc SD.txt @@ -0,0 +1,65 @@ +# X Y E +1 +0.004056 1231.66 21.0814 +0.00438048 903.487 12.6614 +0.00473092 657.371 7.23208 +0.00510939 491.896 4.65659 +0.00551814 383.556 3.20568 +0.00595959 300.487 2.226 +0.00643636 236.132 1.61404 +0.00695127 191.016 1.19379 +0.00750737 154.596 0.896938 +0.00810796 126.281 0.679435 +0.0087566 103.86 0.527471 +0.00945713 85.2982 0.409737 +0.0102137 70.6566 0.323172 +0.0110308 56.9606 0.253332 +0.0119133 45.5839 0.199817 +0.0128663 36.0221 0.156881 +0.0138956 28.1599 0.123091 +0.0150073 21.7195 0.096631 +0.0162079 16.61 0.0757612 +0.0175045 12.7726 0.060008 +0.0189048 9.65055 0.0475836 +0.0204172 7.53418 0.038362 +0.0220506 5.90617 0.031078 +0.0238147 4.66578 0.0255239 +0.0257198 3.85815 0.0213397 +0.0277774 3.22471 0.0181047 +0.0299996 2.69734 0.015346 +0.0323996 2.29566 0.0131455 +0.0349915 1.92759 0.0111743 +0.0377909 1.56761 0.00943598 +0.0408141 1.27379 0.00796468 +0.0440793 0.988188 0.00664852 +0.0476056 0.789308 0.00564 +0.0514141 0.620302 0.00477223 +0.0555272 0.506472 0.00413974 +0.0599694 0.397452 0.00355314 +0.0647669 0.335179 0.00315824 +0.0699483 0.28453 0.00282602 +0.0755441 0.234323 0.00254771 +0.0815876 0.197999 0.00232224 +0.0881147 0.172155 0.00213735 +0.0951638 0.149657 0.00200684 +0.102777 0.129361 0.00188584 +0.110999 0.111647 0.00180023 +0.119879 0.0964134 0.00172185 +0.129469 0.0881419 0.00167425 +0.139827 0.0782058 0.00164565 +0.151013 0.0704272 0.00164037 +0.163094 0.063286 0.0016955 +0.176142 0.0563907 0.00176072 +0.190233 0.0579218 0.00190181 +0.205452 0.0530909 0.00206054 +0.221888 0.0474212 0.0022456 +0.239639 0.044957 0.00244592 +0.25881 0.0410137 0.00265731 +0.279515 0.0310335 0.00287899 +0.301876 0.0284024 0.00319526 +0.326026 0.0363374 0.00356354 +0.352108 0.0331389 0.00393627 +0.380277 0.0340982 0.00438344 +0.410699 0.0264766 0.00490149 +0.443555 0.035785 0.00570573 +0.480307 0.0220708 0.00660036 diff --git a/example/Beise009798_01.ses b/example/Beise009798_01.ses new file mode 100644 index 000000000..1b43c95bc --- /dev/null +++ b/example/Beise009798_01.ses @@ -0,0 +1,91 @@ +DataFileTitle SBA 5%, D2O etc, @ Short position, Sine ++ only, NOT shaken +Sample 6% CaCas +Settings D1=D2=20x8 mm,Ds' = 14x8 mm (WxH), GF1 =scanning, GF2 = 2.5 A. +Operator BT +Date wo 3 feb 2016 11:53:17 +ScanType sine one element scan +Thickness [cm]  +Q_zmax [\A^-1] 0.05 +Q_ymax [\A^-1] 0.05 + +spin echo length [A] Depolarisation [A-2 cm-1] error depol [A-2 cm-1] error SEL [A] wavelength [A] error wavelength [A] polarisation error pol +85.566 -0.00018001 0.00071881 4.2783 2.11 0.1055 0.9992 0.0031976 +134.65 -0.00078142 0.00074052 6.7326 2.11 0.1055 0.99653 0.0032854 +184.58 -0.0015094 0.00078507 9.2289 2.11 0.1055 0.9933 0.0034718 +234.52 0.00063399 0.00074764 11.726 2.11 0.1055 1.0028 0.003338 +285.02 -0.00023297 0.00078333 14.251 2.11 0.1055 0.99896 0.0034838 +336.04 -0.0016295 0.00079926 16.802 2.11 0.1055 0.99277 0.0035327 +387.53 -0.00184 0.0007789 19.376 2.11 0.1055 0.99184 0.0034395 +439.53 -0.001342 0.00074893 21.976 2.11 0.1055 0.99404 0.0033144 +492.41 -0.002132 0.00072869 24.621 2.11 0.1055 0.99055 0.0032136 +545.58 -0.00041763 0.00072165 27.279 2.11 0.1055 0.99814 0.0032069 +599.67 -0.0025597 0.00071984 29.983 2.11 0.1055 0.98867 0.0031685 +654.37 -0.0019647 0.00070481 32.719 2.11 0.1055 0.99129 0.0031105 +709.85 -0.0072549 0.00070816 35.493 2.11 0.1055 0.96822 0.0030526 +766.32 -0.0045379 0.00070593 38.316 2.11 0.1055 0.98 0.00308 +823.47 -0.0040602 0.00071259 41.173 2.11 0.1055 0.98209 0.0031157 +881.71 -0.0050426 0.00072632 44.085 2.11 0.1055 0.9778 0.0031619 +941.01 -0.0063506 0.00074944 47.051 2.11 0.1055 0.97212 0.0032436 +1001.4 -0.0060148 0.00076607 50.07 2.11 0.1055 0.97358 0.0033205 +1063.1 -0.0061236 0.00076025 53.157 2.11 0.1055 0.97311 0.0032937 +1126 -0.0079583 0.00074065 56.299 2.11 0.1055 0.96519 0.0031827 +1190.3 -0.0083306 0.00073057 59.516 2.11 0.1055 0.96359 0.0031342 +1256 -0.0093108 0.00071948 62.801 2.11 0.1055 0.95939 0.0030731 +1323.5 -0.0077527 0.00069781 66.174 2.11 0.1055 0.96607 0.0030013 +1392.7 -0.0093389 0.00069307 69.633 2.11 0.1055 0.95927 0.00296 +1463.5 -0.010127 0.00071306 73.176 2.11 0.1055 0.95591 0.0030347 +1536.4 -0.0098888 0.00073211 76.82 2.11 0.1055 0.95693 0.003119 +1611.6 -0.010577 0.00074926 80.58 2.11 0.1055 0.954 0.0031824 +1688.8 -0.010665 0.00077448 84.44 2.11 0.1055 0.95363 0.0032882 +1768.6 -0.011789 0.00076507 88.431 2.11 0.1055 0.94887 0.003232 +1851.1 -0.011337 0.00077711 92.557 2.11 0.1055 0.95078 0.0032895 +1936.2 -0.0116 0.00080642 96.811 2.11 0.1055 0.94966 0.0034095 +2024.4 -0.012613 0.00083683 101.22 2.11 0.1055 0.94539 0.0035222 +2116 -0.012266 0.00087268 105.8 2.11 0.1055 0.94686 0.0036788 +2210.8 -0.013226 0.00084958 110.54 2.11 0.1055 0.94282 0.0035661 +2309.5 -0.013611 0.00081592 115.47 2.11 0.1055 0.9412 0.003419 +2412.2 -0.014943 0.00078266 120.61 2.11 0.1055 0.93564 0.0032602 +2519.3 -0.016327 0.00079435 125.97 2.11 0.1055 0.92989 0.0032886 +2630.8 -0.01548 0.00082678 131.54 2.11 0.1055 0.9334 0.0034357 +2747.4 -0.014239 0.00087384 137.37 2.11 0.1055 0.93857 0.0036514 +2869.4 -0.01676 0.00083787 143.47 2.11 0.1055 0.9281 0.0034621 +2996.9 -0.016441 0.00079794 149.85 2.11 0.1055 0.92942 0.0033018 +3130.5 -0.016089 0.00080484 156.53 2.11 0.1055 0.93088 0.0033355 +3270.9 -0.017261 0.00084897 163.55 2.11 0.1055 0.92603 0.0035001 +3418.2 -0.016636 0.00085359 170.91 2.11 0.1055 0.92861 0.003529 +3573.1 -0.01796 0.00079916 178.66 2.11 0.1055 0.92315 0.0032845 +3736.1 -0.017599 0.00078142 186.81 2.11 0.1055 0.92464 0.0032168 +3908 -0.018449 0.0008443 195.4 2.11 0.1055 0.92115 0.0034625 +4089.2 -0.018261 0.00082252 204.46 2.11 0.1055 0.92192 0.003376 +4280.5 -0.018156 0.00077108 214.02 2.11 0.1055 0.92235 0.0031664 +4482.3 -0.017473 0.0008328 224.11 2.11 0.1055 0.92516 0.0034302 +4695.9 -0.017402 0.00084046 234.8 2.11 0.1055 0.92545 0.0034629 +4921.9 -0.018548 0.00077627 246.1 2.11 0.1055 0.92074 0.0031821 +5161.3 -0.020108 0.00074633 258.06 2.11 0.1055 0.91437 0.0030382 +5414.8 -0.019294 0.00075007 270.74 2.11 0.1055 0.91769 0.0030645 +5683.8 -0.01921 0.00079006 284.19 2.11 0.1055 0.91803 0.0032291 +5969.3 -0.017449 0.00074079 298.47 2.11 0.1055 0.92525 0.0030516 +6272.4 -0.01757 0.00077985 313.62 2.11 0.1055 0.92476 0.0032107 +6594.8 -0.018636 0.00075007 329.74 2.11 0.1055 0.92038 0.0030735 +6937.2 -0.017352 0.00077649 346.86 2.11 0.1055 0.92566 0.0032 +7301.8 -0.018904 0.00077855 365.09 2.11 0.1055 0.91928 0.0031864 +7689.6 -0.018071 0.0007507 384.48 2.11 0.1055 0.9227 0.0030838 +8102.7 -0.018453 0.00081005 405.14 2.11 0.1055 0.92113 0.003322 +8542.9 -0.017641 0.00086596 427.15 2.11 0.1055 0.92447 0.0035641 +9012 -0.019005 0.0009162 450.6 2.11 0.1055 0.91887 0.0037481 +9512.5 -0.017995 0.00089011 475.63 2.11 0.1055 0.92301 0.0036578 +10046 -0.017917 0.00088139 502.31 2.11 0.1055 0.92333 0.0036232 +10616 -0.019032 0.00087904 530.79 2.11 0.1055 0.91876 0.0035956 +11224 -0.018727 0.0008816 561.2 2.11 0.1055 0.92001 0.003611 +11874 -0.017395 0.00087074 593.69 2.11 0.1055 0.92548 0.0035877 +12568 -0.020032 0.00086947 628.38 2.11 0.1055 0.91468 0.0035407 +13309 -0.020078 0.00087639 665.46 2.11 0.1055 0.91449 0.0035681 +14102 -0.020734 0.0008901 705.09 2.11 0.1055 0.91182 0.0036134 +14949 -0.01854 0.000971 747.47 2.11 0.1055 0.92077 0.0039805 +15856 -0.018621 0.00105 792.78 2.11 0.1055 0.92044 0.0043029 +16825 -0.017531 0.0010941 841.25 2.11 0.1055 0.92492 0.0045054 +17863 -0.019811 0.0010575 893.13 2.11 0.1055 0.91558 0.0043108 +18973 -0.017963 0.0011122 948.63 2.11 0.1055 0.92314 0.0045711 +20161 -0.020882 0.0011029 1008 2.11 0.1055 0.91122 0.0044744 +21433 -0.01926 0.0012863 1071.7 2.11 0.1055 0.91782 0.0052561 +22795 -0.020119 0.0013739 1139.8 2.11 0.1055 0.91432 0.0055927 diff --git a/example/Beise009798_02.ses b/example/Beise009798_02.ses new file mode 100644 index 000000000..b170e1249 --- /dev/null +++ b/example/Beise009798_02.ses @@ -0,0 +1,91 @@ +DataFileTitle SBA 5%, D2O etc, @ Short position, Sine ++ only, NOT shaken +Sample 12% CaCas +Settings D1=D2=20x8 mm,Ds' = 14x8 mm (WxH), GF1 =scanning, GF2 = 2.5 A. +Operator BT +Date wo 3 feb 2016 11:53:17 +ScanType sine one element scan +Thickness [cm]  +Q_zmax [\A^-1] 0.05 +Q_ymax [\A^-1] 0.05 + +spin echo length [A] Depolarisation [A-2 cm-1] error depol [A-2 cm-1] error SEL [A] wavelength [A] error wavelength [A] polarisation error pol +85.566 -0.0013573 0.00076786 4.2783 2.11 0.1055 0.99398 0.003398 +134.65 -0.00042677 0.00078731 6.7326 2.11 0.1055 0.9981 0.0034985 +184.58 -0.0010147 0.00082955 9.2289 2.11 0.1055 0.99549 0.0036766 +234.52 -0.00171 0.0008052 11.726 2.11 0.1055 0.99242 0.0035577 +285.02 -0.0011597 0.00083535 14.251 2.11 0.1055 0.99485 0.0036999 +336.04 -0.0035924 0.00085393 16.802 2.11 0.1055 0.98413 0.0037415 +387.53 -0.0029474 0.00083597 19.376 2.11 0.1055 0.98696 0.0036733 +439.53 -0.004379 0.00080319 21.976 2.11 0.1055 0.98069 0.0035068 +492.41 -0.0032451 0.00078158 24.621 2.11 0.1055 0.98566 0.0034298 +545.58 -0.0054541 0.00078838 27.279 2.11 0.1055 0.97601 0.0034257 +599.67 -0.0069394 0.00078029 29.983 2.11 0.1055 0.96958 0.0033682 +654.37 -0.0067755 0.00076938 32.719 2.11 0.1055 0.97029 0.0033236 +709.85 -0.009938 0.00076426 35.493 2.11 0.1055 0.95672 0.0032553 +766.32 -0.0096871 0.00077257 38.316 2.11 0.1055 0.95779 0.0032944 +823.47 -0.0084367 0.00077256 41.173 2.11 0.1055 0.96314 0.0033127 +881.71 -0.0088836 0.00078965 44.085 2.11 0.1055 0.96122 0.0033793 +941.01 -0.010687 0.00081352 47.051 2.11 0.1055 0.95354 0.0034536 +1001.4 -0.01057 0.00082919 50.07 2.11 0.1055 0.95403 0.0035219 +1063.1 -0.011139 0.00082387 53.157 2.11 0.1055 0.95162 0.0034905 +1126 -0.014209 0.00081277 56.299 2.11 0.1055 0.9387 0.0033967 +1190.3 -0.015686 0.000804 59.516 2.11 0.1055 0.93255 0.003338 +1256 -0.015399 0.0007846 62.801 2.11 0.1055 0.93374 0.0032616 +1323.5 -0.015742 0.00077507 66.174 2.11 0.1055 0.93232 0.0032171 +1392.7 -0.016249 0.00076562 69.633 2.11 0.1055 0.93021 0.0031707 +1463.5 -0.018137 0.00079379 73.176 2.11 0.1055 0.92243 0.0032599 +1536.4 -0.01749 0.00080971 76.82 2.11 0.1055 0.92509 0.0033349 +1611.6 -0.01865 0.00083145 80.58 2.11 0.1055 0.92032 0.0034068 +1688.8 -0.019166 0.00085411 84.44 2.11 0.1055 0.91821 0.0034916 +1768.6 -0.018526 0.00084525 88.431 2.11 0.1055 0.92083 0.0034652 +1851.1 -0.020626 0.0008628 92.557 2.11 0.1055 0.91226 0.0035043 +1936.2 -0.020195 0.00088415 96.811 2.11 0.1055 0.91401 0.0035979 +2024.4 -0.02309 0.00093632 101.22 2.11 0.1055 0.90231 0.0037614 +2116 -0.021222 0.00096216 105.8 2.11 0.1055 0.90984 0.0038974 +2210.8 -0.021873 0.00093881 110.54 2.11 0.1055 0.90721 0.0037918 +2309.5 -0.02239 0.00090297 115.47 2.11 0.1055 0.90512 0.0036387 +2412.2 -0.023905 0.00087241 120.61 2.11 0.1055 0.89904 0.0034919 +2519.3 -0.02531 0.00088202 125.97 2.11 0.1055 0.89344 0.0035084 +2630.8 -0.024152 0.00091392 131.54 2.11 0.1055 0.89805 0.003654 +2747.4 -0.025189 0.0009793 137.37 2.11 0.1055 0.89392 0.0038974 +2869.4 -0.027767 0.00093671 143.47 2.11 0.1055 0.88371 0.0036854 +2996.9 -0.025304 0.00088801 149.85 2.11 0.1055 0.89346 0.0035323 +3130.5 -0.027575 0.00090803 156.53 2.11 0.1055 0.88447 0.0035756 +3270.9 -0.026555 0.00094156 163.55 2.11 0.1055 0.8885 0.0037245 +3418.2 -0.026484 0.00094427 170.91 2.11 0.1055 0.88878 0.0037364 +3573.1 -0.026045 0.00088943 178.66 2.11 0.1055 0.89051 0.0035263 +3736.1 -0.02714 0.000871 186.81 2.11 0.1055 0.88618 0.0034364 +3908 -0.027349 0.00093219 195.4 2.11 0.1055 0.88536 0.0036744 +4089.2 -0.028458 0.00091508 204.46 2.11 0.1055 0.881 0.0035892 +4280.5 -0.027419 0.00086098 214.02 2.11 0.1055 0.88508 0.0033927 +4482.3 -0.027101 0.00091969 224.11 2.11 0.1055 0.88634 0.0036292 +4695.9 -0.02763 0.00094596 234.8 2.11 0.1055 0.88426 0.003724 +4921.9 -0.027258 0.00086491 246.1 2.11 0.1055 0.88572 0.0034106 +5161.3 -0.027401 0.00082853 258.06 2.11 0.1055 0.88516 0.0032651 +5414.8 -0.027698 0.00083429 270.74 2.11 0.1055 0.88398 0.0032834 +5683.8 -0.027046 0.00087081 284.19 2.11 0.1055 0.88655 0.0034371 +5969.3 -0.027659 0.00082332 298.47 2.11 0.1055 0.88414 0.0032408 +6272.4 -0.02659 0.00086779 313.62 2.11 0.1055 0.88836 0.0034322 +6594.8 -0.027022 0.00083144 329.74 2.11 0.1055 0.88665 0.0032821 +6937.2 -0.026556 0.00086642 346.86 2.11 0.1055 0.88849 0.0034272 +7301.8 -0.02814 0.00086337 365.09 2.11 0.1055 0.88225 0.0033912 +7689.6 -0.027075 0.00083927 384.48 2.11 0.1055 0.88644 0.0033122 +8102.7 -0.02749 0.00090678 405.14 2.11 0.1055 0.88481 0.003572 +8542.9 -0.02817 0.00097409 427.15 2.11 0.1055 0.88213 0.0038256 +9012 -0.026146 0.00099712 450.6 2.11 0.1055 0.89012 0.0039515 +9512.5 -0.026982 0.00099049 475.63 2.11 0.1055 0.88681 0.0039106 +10046 -0.027148 0.00097204 502.31 2.11 0.1055 0.88615 0.0038349 +10616 -0.025813 0.00096382 530.79 2.11 0.1055 0.89144 0.0038252 +11224 -0.027438 0.00097686 561.2 2.11 0.1055 0.88501 0.003849 +11874 -0.028203 0.00098246 593.69 2.11 0.1055 0.882 0.0038579 +12568 -0.029499 0.00097006 628.38 2.11 0.1055 0.87693 0.0037873 +13309 -0.027996 0.00096756 665.46 2.11 0.1055 0.88281 0.0038029 +14102 -0.029864 0.00098435 705.09 2.11 0.1055 0.8755 0.0038368 +14949 -0.028606 0.0010866 747.47 2.11 0.1055 0.88042 0.0042592 +15856 -0.02789 0.0011605 792.78 2.11 0.1055 0.88323 0.0045634 +16825 -0.027305 0.0012079 841.25 2.11 0.1055 0.88553 0.0047621 +17863 -0.029264 0.0011722 893.13 2.11 0.1055 0.87784 0.0045813 +18973 -0.028019 0.0012295 948.63 2.11 0.1055 0.88272 0.004832 +20161 -0.031972 0.0012303 1008 2.11 0.1055 0.86732 0.0047508 +21433 -0.031008 0.0014282 1071.7 2.11 0.1055 0.87105 0.0055388 +22795 -0.030828 0.0015173 1139.8 2.11 0.1055 0.87175 0.005889 diff --git a/example/Beise009798_03.ses b/example/Beise009798_03.ses new file mode 100644 index 000000000..edabab6a1 --- /dev/null +++ b/example/Beise009798_03.ses @@ -0,0 +1,91 @@ +DataFileTitle SBA 5%, D2O etc, @ Short position, Sine ++ only, NOT shaken +Sample 18% CaCas +Settings D1=D2=20x8 mm,Ds' = 14x8 mm (WxH), GF1 =scanning, GF2 = 2.5 A. +Operator BT +Date wo 3 feb 2016 11:53:17 +ScanType sine one element scan +Thickness [cm]  +Q_zmax [\A^-1] 0.05 +Q_ymax [\A^-1] 0.05 + +spin echo length [A] Depolarisation [A-2 cm-1] error depol [A-2 cm-1] error SEL [A] wavelength [A] error wavelength [A] polarisation error pol +85.566 -0.00093161 0.00081941 4.2783 2.11 0.1055 0.99586 0.003633 +134.65 -0.0018955 0.00083937 6.7326 2.11 0.1055 0.9916 0.0037056 +184.58 -0.001239 0.00088616 9.2289 2.11 0.1055 0.9945 0.0039236 +234.52 -0.0016163 0.00085367 11.726 2.11 0.1055 0.99283 0.0037734 +285.02 -0.0043577 0.00089744 14.251 2.11 0.1055 0.98079 0.0039187 +336.04 -0.0046675 0.0009078 16.802 2.11 0.1055 0.97943 0.0039585 +387.53 -0.004882 0.00089861 19.376 2.11 0.1055 0.9785 0.0039147 +439.53 -0.0059567 0.00086489 21.976 2.11 0.1055 0.97383 0.0037498 +492.41 -0.0049538 0.0008391 24.621 2.11 0.1055 0.97819 0.0036543 +545.58 -0.0058798 0.00084218 27.279 2.11 0.1055 0.97416 0.0036526 +599.67 -0.0094648 0.00084814 29.983 2.11 0.1055 0.95874 0.0036202 +654.37 -0.0085715 0.00082442 32.719 2.11 0.1055 0.96256 0.003533 +709.85 -0.010518 0.00082071 35.493 2.11 0.1055 0.95425 0.0034867 +766.32 -0.0086817 0.00081638 38.316 2.11 0.1055 0.96209 0.0034968 +823.47 -0.010962 0.00083145 41.173 2.11 0.1055 0.95237 0.0035254 +881.71 -0.011815 0.00085069 44.085 2.11 0.1055 0.94876 0.0035933 +941.01 -0.011222 0.00087706 47.051 2.11 0.1055 0.95127 0.0037145 +1001.4 -0.013706 0.00090257 50.07 2.11 0.1055 0.94081 0.0037805 +1063.1 -0.013354 0.00089079 53.157 2.11 0.1055 0.94228 0.003737 +1126 -0.017477 0.00088102 56.299 2.11 0.1055 0.92514 0.0036288 +1190.3 -0.016154 0.00085934 59.516 2.11 0.1055 0.9306 0.0035604 +1256 -0.017295 0.0008411 62.801 2.11 0.1055 0.92589 0.0034672 +1323.5 -0.018177 0.00083528 66.174 2.11 0.1055 0.92226 0.0034297 +1392.7 -0.019546 0.00083231 69.633 2.11 0.1055 0.91666 0.0033967 +1463.5 -0.018236 0.00084938 73.176 2.11 0.1055 0.92202 0.0034867 +1536.4 -0.018868 0.00086985 76.82 2.11 0.1055 0.91943 0.0035606 +1611.6 -0.02081 0.00089671 80.58 2.11 0.1055 0.91151 0.003639 +1688.8 -0.021007 0.00092412 84.44 2.11 0.1055 0.91072 0.0037469 +1768.6 -0.02258 0.00091166 88.431 2.11 0.1055 0.90436 0.0036706 +1851.1 -0.023449 0.00093925 92.557 2.11 0.1055 0.90087 0.0037671 +1936.2 -0.021908 0.00094932 96.811 2.11 0.1055 0.90707 0.0038337 +2024.4 -0.02322 0.00099073 101.22 2.11 0.1055 0.90179 0.0039777 +2116 -0.023607 0.0010431 105.8 2.11 0.1055 0.90023 0.0041808 +2210.8 -0.022224 0.001001 110.54 2.11 0.1055 0.90579 0.0040368 +2309.5 -0.02295 0.00096699 115.47 2.11 0.1055 0.90287 0.003887 +2412.2 -0.024982 0.00093832 120.61 2.11 0.1055 0.89474 0.0037378 +2519.3 -0.0267 0.00094607 125.97 2.11 0.1055 0.88792 0.0037399 +2630.8 -0.025881 0.00099024 131.54 2.11 0.1055 0.89117 0.0039289 +2747.4 -0.02587 0.0010509 137.37 2.11 0.1055 0.89121 0.0041699 +2869.4 -0.026249 0.001001 143.47 2.11 0.1055 0.88971 0.003965 +2996.9 -0.025131 0.00095182 149.85 2.11 0.1055 0.89415 0.003789 +3130.5 -0.025426 0.00095098 156.53 2.11 0.1055 0.89297 0.0037807 +3270.9 -0.026334 0.0010089 163.55 2.11 0.1055 0.88937 0.0039949 +3418.2 -0.025836 0.0010176 170.91 2.11 0.1055 0.89134 0.0040384 +3573.1 -0.025482 0.00094325 178.66 2.11 0.1055 0.89275 0.0037491 +3736.1 -0.025594 0.00093094 186.81 2.11 0.1055 0.89231 0.0036983 +3908 -0.025207 0.00099237 195.4 2.11 0.1055 0.89384 0.0039491 +4089.2 -0.027043 0.00098036 204.46 2.11 0.1055 0.88657 0.0038696 +4280.5 -0.027299 0.00091603 214.02 2.11 0.1055 0.88556 0.0036116 +4482.3 -0.026166 0.0009785 224.11 2.11 0.1055 0.89003 0.0038773 +4695.9 -0.025901 0.001006 234.8 2.11 0.1055 0.89109 0.0039909 +4921.9 -0.026386 0.00092292 246.1 2.11 0.1055 0.88916 0.0036535 +5161.3 -0.028562 0.00089093 258.06 2.11 0.1055 0.88059 0.0034929 +5414.8 -0.026413 0.00088955 270.74 2.11 0.1055 0.88906 0.003521 +5683.8 -0.026904 0.00094025 284.19 2.11 0.1055 0.88711 0.0037135 +5969.3 -0.027137 0.00087956 298.47 2.11 0.1055 0.8862 0.0034703 +6272.4 -0.027597 0.0009387 313.62 2.11 0.1055 0.88439 0.003696 +6594.8 -0.02686 0.00088706 329.74 2.11 0.1055 0.88729 0.0035041 +6937.2 -0.025797 0.00093116 346.86 2.11 0.1055 0.8915 0.0036958 +7301.8 -0.026455 0.00091303 365.09 2.11 0.1055 0.88889 0.0036133 +7689.6 -0.025799 0.00088226 384.48 2.11 0.1055 0.89149 0.0035017 +8102.7 -0.028113 0.00097679 405.14 2.11 0.1055 0.88236 0.0038372 +8542.9 -0.026597 0.0010345 427.15 2.11 0.1055 0.88833 0.0040916 +9012 -0.026901 0.001078 450.6 2.11 0.1055 0.88713 0.0042576 +9512.5 -0.028878 0.0010695 475.63 2.11 0.1055 0.87936 0.0041869 +10046 -0.027146 0.0010342 502.31 2.11 0.1055 0.88616 0.0040802 +10616 -0.027076 0.0010354 530.79 2.11 0.1055 0.88644 0.0040864 +11224 -0.026394 0.0010397 561.2 2.11 0.1055 0.88913 0.0041157 +11874 -0.02899 0.0010494 593.69 2.11 0.1055 0.87891 0.0041065 +12568 -0.027433 0.0010354 628.38 2.11 0.1055 0.88503 0.0040799 +13309 -0.02857 0.0010465 665.46 2.11 0.1055 0.88056 0.0041028 +14102 -0.030266 0.0010601 705.09 2.11 0.1055 0.87394 0.0041247 +14949 -0.026457 0.0011507 747.47 2.11 0.1055 0.88888 0.0045536 +15856 -0.027132 0.001235 792.78 2.11 0.1055 0.88622 0.0048727 +16825 -0.028503 0.0013094 841.25 2.11 0.1055 0.88082 0.005135 +17863 -0.028169 0.0012591 893.13 2.11 0.1055 0.88213 0.0049449 +18973 -0.027315 0.0013031 948.63 2.11 0.1055 0.88549 0.0051371 +20161 -0.029182 0.0013068 1008 2.11 0.1055 0.87816 0.0051093 +21433 -0.030792 0.0015285 1071.7 2.11 0.1055 0.87189 0.0059334 +22795 -0.02962 0.0016099 1139.8 2.11 0.1055 0.87645 0.0062818 diff --git a/example/Gab_5pc_deb2p0.ses b/example/Gab_5pc_deb2p0.ses new file mode 100644 index 000000000..8852201ca --- /dev/null +++ b/example/Gab_5pc_deb2p0.ses @@ -0,0 +1,36 @@ +DataFileTitle 5pc_deb2p0 +Sample 5pc_deb2p0 +Settings blabla +Operator CPd +Date wo 28 jun 2016 14:30 +ScanType sine one element scan +Thickness[cm] 0.2 +Q_zmax[\A^-1] 0.05 +Q_ymax[\A^-1] 0.05 + +spin echo length [A] Depolarisation [A-2 cm-1] error depol [A-2 cm-1] error SEL [A] wavelength [A] error wavelength [A] polarization error pol +8.91E+02 -0.005000634 3.68E-03 4.46E+01 2.11 0.1055 0.995557234 3.68E-03 +2.57E+03 -0.030114103 3.52E-03 1.28E+02 2.11 0.1055 0.973542109 3.52E-03 +4.26E+03 -0.067322739 3.43E-03 2.13E+02 2.11 0.1055 0.941815849 3.43E-03 +5.99E+03 -0.100439625 3.23E-03 2.99E+02 2.11 0.1055 0.914449119 3.23E-03 +7.74E+03 -0.138541245 3.22E-03 3.87E+02 2.11 0.1055 0.883945477 3.22E-03 +9.54E+03 -0.163854409 3.33E-03 4.77E+02 2.11 0.1055 0.864244778 3.33E-03 +1.14E+04 -0.196263409 3.05E-03 5.70E+02 2.11 0.1055 0.839661147 3.05E-03 +1.33E+04 -0.224571679 3.06E-03 6.66E+02 2.11 0.1055 0.818760953 3.06E-03 +1.53E+04 -0.251404993 3.08E-03 7.67E+02 2.11 0.1055 0.799430217 3.08E-03 +1.75E+04 -0.274978996 3.11E-03 8.75E+02 2.11 0.1055 0.78282446 3.11E-03 +1.98E+04 -0.306175187 3.09E-03 9.90E+02 2.11 0.1055 0.761378624 3.09E-03 +2.23E+04 -0.333022802 3.12E-03 1.12E+03 2.11 0.1055 0.743393204 3.12E-03 +2.51E+04 -0.365092007 3.34E-03 1.26E+03 2.11 0.1055 0.722465779 3.34E-03 +2.82E+04 -0.394694216 3.40E-03 1.41E+03 2.11 0.1055 0.703671519 3.40E-03 +3.18E+04 -0.435265203 3.34E-03 1.59E+03 2.11 0.1055 0.678704909 3.34E-03 +3.60E+04 -0.466568709 3.28E-03 1.80E+03 2.11 0.1055 0.660048402 3.28E-03 +4.10E+04 -0.506015307 3.27E-03 2.05E+03 2.11 0.1055 0.637267261 3.27E-03 +4.69E+04 -0.538798205 3.49E-03 2.35E+03 2.11 0.1055 0.618933962 3.49E-03 +5.42E+04 -0.581251553 3.72E-03 2.71E+03 2.11 0.1055 0.595974137 3.72E-03 +6.31E+04 -0.633703036 3.83E-03 3.15E+03 2.11 0.1055 0.568779833 3.83E-03 +7.42E+04 -0.671815932 4.22E-03 3.71E+03 2.11 0.1055 0.549801298 4.22E-03 +8.81E+04 -0.70514095 4.86E-03 4.40E+03 2.11 0.1055 0.533726574 4.86E-03 +1.06E+05 -0.766341597 5.68E-03 5.28E+03 2.11 0.1055 0.505419812 5.68E-03 +1.28E+05 -0.811142639 6.29E-03 6.40E+03 2.11 0.1055 0.48565459 6.29E-03 +1.56E+05 -0.835418913 9.16E-03 7.82E+03 2.11 0.1055 0.47526929 9.16E-03 diff --git a/example/Gab_5pc_deb2p0.txt b/example/Gab_5pc_deb2p0.txt new file mode 100644 index 000000000..8852201ca --- /dev/null +++ b/example/Gab_5pc_deb2p0.txt @@ -0,0 +1,36 @@ +DataFileTitle 5pc_deb2p0 +Sample 5pc_deb2p0 +Settings blabla +Operator CPd +Date wo 28 jun 2016 14:30 +ScanType sine one element scan +Thickness[cm] 0.2 +Q_zmax[\A^-1] 0.05 +Q_ymax[\A^-1] 0.05 + +spin echo length [A] Depolarisation [A-2 cm-1] error depol [A-2 cm-1] error SEL [A] wavelength [A] error wavelength [A] polarization error pol +8.91E+02 -0.005000634 3.68E-03 4.46E+01 2.11 0.1055 0.995557234 3.68E-03 +2.57E+03 -0.030114103 3.52E-03 1.28E+02 2.11 0.1055 0.973542109 3.52E-03 +4.26E+03 -0.067322739 3.43E-03 2.13E+02 2.11 0.1055 0.941815849 3.43E-03 +5.99E+03 -0.100439625 3.23E-03 2.99E+02 2.11 0.1055 0.914449119 3.23E-03 +7.74E+03 -0.138541245 3.22E-03 3.87E+02 2.11 0.1055 0.883945477 3.22E-03 +9.54E+03 -0.163854409 3.33E-03 4.77E+02 2.11 0.1055 0.864244778 3.33E-03 +1.14E+04 -0.196263409 3.05E-03 5.70E+02 2.11 0.1055 0.839661147 3.05E-03 +1.33E+04 -0.224571679 3.06E-03 6.66E+02 2.11 0.1055 0.818760953 3.06E-03 +1.53E+04 -0.251404993 3.08E-03 7.67E+02 2.11 0.1055 0.799430217 3.08E-03 +1.75E+04 -0.274978996 3.11E-03 8.75E+02 2.11 0.1055 0.78282446 3.11E-03 +1.98E+04 -0.306175187 3.09E-03 9.90E+02 2.11 0.1055 0.761378624 3.09E-03 +2.23E+04 -0.333022802 3.12E-03 1.12E+03 2.11 0.1055 0.743393204 3.12E-03 +2.51E+04 -0.365092007 3.34E-03 1.26E+03 2.11 0.1055 0.722465779 3.34E-03 +2.82E+04 -0.394694216 3.40E-03 1.41E+03 2.11 0.1055 0.703671519 3.40E-03 +3.18E+04 -0.435265203 3.34E-03 1.59E+03 2.11 0.1055 0.678704909 3.34E-03 +3.60E+04 -0.466568709 3.28E-03 1.80E+03 2.11 0.1055 0.660048402 3.28E-03 +4.10E+04 -0.506015307 3.27E-03 2.05E+03 2.11 0.1055 0.637267261 3.27E-03 +4.69E+04 -0.538798205 3.49E-03 2.35E+03 2.11 0.1055 0.618933962 3.49E-03 +5.42E+04 -0.581251553 3.72E-03 2.71E+03 2.11 0.1055 0.595974137 3.72E-03 +6.31E+04 -0.633703036 3.83E-03 3.15E+03 2.11 0.1055 0.568779833 3.83E-03 +7.42E+04 -0.671815932 4.22E-03 3.71E+03 2.11 0.1055 0.549801298 4.22E-03 +8.81E+04 -0.70514095 4.86E-03 4.40E+03 2.11 0.1055 0.533726574 4.86E-03 +1.06E+05 -0.766341597 5.68E-03 5.28E+03 2.11 0.1055 0.505419812 5.68E-03 +1.28E+05 -0.811142639 6.29E-03 6.40E+03 2.11 0.1055 0.48565459 6.29E-03 +1.56E+05 -0.835418913 9.16E-03 7.82E+03 2.11 0.1055 0.47526929 9.16E-03 diff --git a/sasmodels/sesans.py b/sasmodels/sesans.py index 9b8b482d0..31bfa713a 100644 --- a/sasmodels/sesans.py +++ b/sasmodels/sesans.py @@ -44,23 +44,25 @@ def __init__(self, z, SElength, zaccept, Rmax): # type: (np.ndarray, float, float) -> None #import logging; logging.info("creating SESANS transform") self.q = z - # isoriented flag determines whether data is from an oriented sample or not, should be in the data or selectable in GUI. - #if self.isoriented==False - self._set_hankel(SElength, zaccept, Rmax) - #if self.isoriented==True - self._set_cosmat(SElength, zaccept, Rmax) + # isoriented flag determines whether data is from an oriented sample or not, should be a selection variable upon entering SESANS data. + if self.isoriented==False: + self._set_hankel(SElength, zaccept, Rmax) + if self.isoriented==True: + self._set_cosmat(SElength, zaccept, Rmax) def apply(self, Iq): - if len(Iq.size) == 1: # if isotropic, do Hankel transform - # type: (np.ndarray) -> np.ndarray + try: G0 = np.dot(self._H0, Iq) G = np.dot(self._H.T, Iq) P = G - G0 - elif len(Iq.size) == 2: - dq=self.q_calc[0] - G0 = sum(np.dot(self._cos0, Iq)*dq) - G = sum(np.dot(self._cosmat.T, Iq)*dq) - P = G - G0 + except: + try: + dq = self.q_calc[0] + G0 = sum(np.dot(self._cos0, Iq) * dq) + G = sum(np.dot(self._cosmat.T, Iq) * dq) + P = G - G0 + except: + raise ValueError('Sesanstransform.apply cannot generate either a Hankel transform or a cosine transform of you SESANS data') return P def _set_hankel(self, SElength, zaccept, Rmax): @@ -87,8 +89,7 @@ def _set_cosmat(self, SElength, zaccept, Rmax): # type: (np.ndarray, float, float) -> None # Force float32 arrays, otherwise run into memory problems on some machines SElength = np.asarray(SElength, dtype='float32') - #qymax and qzmax depend on detector shape - qzmax= + # Rmax = #value in text box somewhere in FitPage? q_max = 2 * pi / (SElength[1] - SElength[0]) q_min = 0.1 * 2 * pi / (np.size(SElength) * SElength[-1]) diff --git a/sasmodels/sesanscos.py b/sasmodels/sesanscos.py new file mode 100644 index 000000000..9b8b482d0 --- /dev/null +++ b/sasmodels/sesanscos.py @@ -0,0 +1,106 @@ +""" +Conversion of scattering cross section from SANS (I(q), or rather, ds/dO) in absolute +units (cm-1)into SESANS correlation function G using a Hankel transformation, then converting +the SESANS correlation function into polarisation from the SESANS experiment + +Everything is in units of metres except specified otherwise (NOT TRUE!!!) +Everything is in conventional units (nm for spin echo length) + +Wim Bouwman (w.g.bouwman@tudelft.nl), June 2013 +""" + +from __future__ import division + +import numpy as np # type: ignore +from numpy import pi, exp # type: ignore +from scipy.special import j0 + +class SesansTransform(object): + """ + Spin-Echo SANS transform calculator. Similar to a resolution function, + the SesansTransform object takes I(q) for the set of *q_calc* values and + produces a transformed dataset + + *SElength* (A) is the set of spin-echo lengths in the measured data. + + *zaccept* (1/A) is the maximum acceptance of scattering vector in the spin + echo encoding dimension (for ToF: Q of min(R) and max(lam)). + + *Rmax* (A) is the maximum size sensitivity; larger radius requires more + computation time. + """ + #: SElength from the data in the original data units; not used by transform + #: but the GUI uses it, so make sure that it is present. + q = None # type: np.ndarray + + #: q values to calculate when computing transform + q_calc = None # type: np.ndarray + + # transform arrays + _H = None # type: np.ndarray + _H0 = None # type: np.ndarray + + def __init__(self, z, SElength, zaccept, Rmax): + # type: (np.ndarray, float, float) -> None + #import logging; logging.info("creating SESANS transform") + self.q = z + # isoriented flag determines whether data is from an oriented sample or not, should be in the data or selectable in GUI. + #if self.isoriented==False + self._set_hankel(SElength, zaccept, Rmax) + #if self.isoriented==True + self._set_cosmat(SElength, zaccept, Rmax) + + def apply(self, Iq): + if len(Iq.size) == 1: # if isotropic, do Hankel transform + # type: (np.ndarray) -> np.ndarray + G0 = np.dot(self._H0, Iq) + G = np.dot(self._H.T, Iq) + P = G - G0 + elif len(Iq.size) == 2: + dq=self.q_calc[0] + G0 = sum(np.dot(self._cos0, Iq)*dq) + G = sum(np.dot(self._cosmat.T, Iq)*dq) + P = G - G0 + return P + + def _set_hankel(self, SElength, zaccept, Rmax): + # type: (np.ndarray, float, float) -> None + # Force float32 arrays, otherwise run into memory problems on some machines + SElength = np.asarray(SElength, dtype='float32') + + #Rmax = #value in text box somewhere in FitPage? + q_max = 2*pi / (SElength[1] - SElength[0]) + q_min = 0.1 * 2*pi / (np.size(SElength) * SElength[-1]) + q = np.arange(q_min, q_max, q_min, dtype='float32') + dq = q_min + + H0 = np.float32(dq/(2*pi)) * q + + repq = np.tile(q, (SElength.size, 1)).T + repSE = np.tile(SElength, (q.size, 1)) + H = np.float32(dq/(2*pi)) * j0(repSE*repq) * repq + + self.q_calc = q + self._H, self._H0 = H, H0 + + def _set_cosmat(self, SElength, zaccept, Rmax): + # type: (np.ndarray, float, float) -> None + # Force float32 arrays, otherwise run into memory problems on some machines + SElength = np.asarray(SElength, dtype='float32') + #qymax and qzmax depend on detector shape + qzmax= + # Rmax = #value in text box somewhere in FitPage? + q_max = 2 * pi / (SElength[1] - SElength[0]) + q_min = 0.1 * 2 * pi / (np.size(SElength) * SElength[-1]) + + q = np.arange(q_min, q_max, q_min, dtype='float32') + dq = q_min + + cos0 = np.float32(dq / (2 * pi)) + + 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) + + self.q_calc = q + self._cosmat, self._cos0 = cosmat, cos0 \ No newline at end of file diff --git a/sasmodels/sesanscosine.py b/sasmodels/sesanscosine.py new file mode 100644 index 000000000..62c4d0b0f --- /dev/null +++ b/sasmodels/sesanscosine.py @@ -0,0 +1,214 @@ +""" +Conversion of scattering cross section from SANS (I(q), or rather, ds/dO) in absolute +units (cm-1)into SESANS correlation function G using a Hankel transformation, then converting +the SESANS correlation function into polarisation from the SESANS experiment + +Everything is in units of metres except specified otherwise (NOT TRUE!!!) +Everything is in conventional units (nm for spin echo length) + +Wim Bouwman (w.g.bouwman@tudelft.nl), June 2013 +""" + +from __future__ import division + +import numpy as np # type: ignore +from numpy import pi, exp # type: ignore +from scipy.special import jv as besselj +#import direct_model.DataMixin as model + +def make_q(q_max, Rmax): + r""" + Return a $q$ vector suitable for SESANS covering from $2\pi/ (10 R_{\max})$ + to $q_max$. This is the integration range of the Hankel transform; bigger range and + more points makes a better numerical integration. + Smaller q_min will increase reliable spin echo length range. + Rmax is the "radius" of the largest expected object and can be set elsewhere. + q_max is determined by the acceptance angle of the SESANS instrument. + """ + from sas.sascalc.data_util.nxsunit import Converter + + q_min = dq = 0.1 * 2*pi / Rmax + return np.arange(q_min, + Converter(q_max[1])(q_max[0], + units="1/A"), + dq) + +def make_all_q(data): + """ + Return a $q$ vector suitable for calculating the total scattering cross section for + calculating the effect of finite acceptance angles on Time of Flight SESANS instruments. + If no acceptance is given, or unwanted (set "unwanted" flag in paramfile), no all_q vector is needed. + If the instrument has a rectangular acceptance, 2 all_q vectors are needed. + If the instrument has a circular acceptance, 1 all_q vector is needed + + """ + if not data.has_no_finite_acceptance: + return [] + elif data.has_yz_acceptance(data): + # compute qx, qy + Qx, Qy = np.meshgrid(qx, qy) + return [Qx, Qy] + else: + # else only need q + # data.has_z_acceptance + return [q] + +def transform(data, q_calc, Iq_calc, qmono, Iq_mono): + """ + Decides which transform type is to be used, based on the experiment data file contents (header) + (2016-03-19: currently controlled from parameters script) + nqmono is the number of q vectors to be used for the detector integration + """ + nqmono = len(qmono) + if nqmono == 0: + result = call_hankel(data, q_calc, Iq_calc) + elif nqmono == 1: + q = qmono[0] + result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono) + else: + Qx, Qy = [qmono[0], qmono[1]] + Qx = np.reshape(Qx, nqx, nqy) + Qy = np.reshape(Qy, nqx, nqy) + Iq_mono = np.reshape(Iq_mono, nqx, nqy) + qx = Qx[0, :] + qy = Qy[:, 0] + result = call_Cosine2D(data, q_calc, Iq_calc, qx, qy, Iq_mono) + + return result + +def call_hankel(data, q_calc, Iq_calc): + return hankel((data.x, data.x_unit), + (data.lam, data.lam_unit), + (data.sample.thickness, + data.sample.thickness_unit), + q_calc, Iq_calc) + +def call_HankelAccept(data, q_calc, Iq_calc, q_mono, Iq_mono): + return hankel(data.x, data.lam * 1e-9, + data.sample.thickness / 10, + 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, + data.sample.thickness / 10, + q_calc, Iq_calc) + +def TotalScatter(model, parameters): #Work in progress!! +# Calls a model with existing model parameters already in place, then integrate the product of q and I(q) from 0 to (4*pi/lambda) + allq = np.linspace(0,4*pi/wavelength,1000) + allIq = 1 + integral = allq*allIq + + + +def Cosine2D(wavelength, magfield, thickness, qy, qz, Iqy, Iqz, modelname): #Work in progress!! Needs to call model still +#============================================================================== +# 2D Cosine Transform if "wavelength" is a vector +#============================================================================== +#allq is the q-space needed to create the total scattering cross-section + + Gprime = np.zeros_like(wavelength, 'd') + s = np.zeros_like(wavelength, 'd') + sd = np.zeros_like(wavelength, 'd') + Gprime = np.zeros_like(wavelength, 'd') + f = np.zeros_like(wavelength, 'd') + for i, wavelength_i in enumerate(wavelength): + z = magfield*wavelength_i + allq=np.linspace() #for calculating the Q-range of the scattering power integral + 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) + for j, Iqy_j, qy_j in enumerate(qy): + for k, Iqz_k, qz_k in enumerate(qz): + Iq = np.sqrt(Iqy_j^2+Iqz_k^2) + q = np.sqrt(qy_j^2 + qz_k^2) + Gintegral = Iq*cos(z*Qz_k) + Gprime[i] += Gintegral +# sigma = wavelength^2*thickness/2/pi* allq[i]*allIq[i] +# s[i] += 1-exp(Totalscatter(modelname)*thickness) +# For now, work with standard 2-phase scatter + + + sd[i] += Iq + f[i] = 1-s[i]+sd[i] + P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i] + + + + +def HankelAccept(wavelength, magfield, thickness, q, Iq, theta, modelname): +#============================================================================== +# 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 + G = np.zeros_like(SElength, 'd') + threshold=2*pi*theta/wavelength + for i, SElength_i in enumerate(SElength): + allq=np.linspace() #for calculating the Q-range of the scattering power integral + 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) + + dq = (q[1]-q[0])*1e10 + a = (x Date: Mon, 10 Apr 2017 09:39:56 -0400 Subject: [PATCH 3/4] oriented sesans: split transform into oriented and unoriented versions --- sasmodels/sesans.py | 78 +++++++++++++++++++++++++++++++-------------- 1 file changed, 54 insertions(+), 24 deletions(-) diff --git a/sasmodels/sesans.py b/sasmodels/sesans.py index 31bfa713a..127079919 100644 --- a/sasmodels/sesans.py +++ b/sasmodels/sesans.py @@ -19,7 +19,7 @@ class SesansTransform(object): """ Spin-Echo SANS transform calculator. Similar to a resolution function, the SesansTransform object takes I(q) for the set of *q_calc* values and - produces a transformed dataset + produces a transformed dataset. *SElength* (A) is the set of spin-echo lengths in the measured data. @@ -44,25 +44,14 @@ def __init__(self, z, SElength, zaccept, Rmax): # type: (np.ndarray, float, float) -> None #import logging; logging.info("creating SESANS transform") self.q = z - # isoriented flag determines whether data is from an oriented sample or not, should be a selection variable upon entering SESANS data. - if self.isoriented==False: - self._set_hankel(SElength, zaccept, Rmax) - if self.isoriented==True: - self._set_cosmat(SElength, zaccept, Rmax) + # isoriented flag determines whether data is from an oriented sample or + # not, should be a selection variable upon entering SESANS data. + self._set_hankel(SElength, zaccept, Rmax) def apply(self, Iq): - try: - G0 = np.dot(self._H0, Iq) - G = np.dot(self._H.T, Iq) - P = G - G0 - except: - try: - dq = self.q_calc[0] - G0 = sum(np.dot(self._cos0, Iq) * dq) - G = sum(np.dot(self._cosmat.T, Iq) * dq) - P = G - G0 - except: - raise ValueError('Sesanstransform.apply cannot generate either a Hankel transform or a cosine transform of you SESANS data') + G0 = np.dot(self._H0, Iq) + G = np.dot(self._H.T, Iq) + P = G - G0 return P def _set_hankel(self, SElength, zaccept, Rmax): @@ -85,6 +74,48 @@ def _set_hankel(self, SElength, zaccept, Rmax): self.q_calc = q self._H, self._H0 = H, H0 +class OrientedSesansTransform(object): + """ + Oriented Spin-Echo SANS transform calculator. Similar to a resolution + function, the OrientedSesansTransform object takes I(q) for the set + of *q_calc* values and produces a transformed dataset. + + *SElength* (A) is the set of spin-echo lengths in the measured data. + + *zaccept* (1/A) is the maximum acceptance of scattering vector in the spin + echo encoding dimension (for ToF: Q of min(R) and max(lam)). + + *Rmax* (A) is the maximum size sensitivity; larger radius requires more + computation time. + """ + #: SElength from the data in the original data units; not used by transform + #: but the GUI uses it, so make sure that it is present. + q = None # type: np.ndarray + + #: q values to calculate when computing transform + q_calc = None # type: np.ndarray + + # transform arrays + _cosmat = None # type: np.ndarray + _cos0 = None # type: np.ndarray + _Iq_shape = None # type: Tuple[int, int] + + def __init__(self, z, SElength, zaccept, Rmax): + # type: (np.ndarray, float, float) -> None + #import logging; logging.info("creating SESANS transform") + self.q = z + # isoriented flag determines whether data is from an oriented sample or + # not, should be a selection variable upon entering SESANS data. + self._set_cosmat(SElength, zaccept, Rmax) + + def apply(self, Iq): + dq = self.q_calc[0][0] + Iq = np.reshape(Iq, self._Iq_shape) + G0 = self._cos0 * np.sum(Iq) * dq + G = np.sum(np.dot(Iq, self._cosmat.T), axis=1) * dq + P = G - G0 + return P + def _set_cosmat(self, SElength, zaccept, Rmax): # type: (np.ndarray, float, float) -> None # Force float32 arrays, otherwise run into memory problems on some machines @@ -98,10 +129,9 @@ def _set_cosmat(self, SElength, zaccept, Rmax): dq = q_min cos0 = np.float32(dq / (2 * pi)) + cosmat = np.float32(dq / (2 * pi)) * np.cos(q[:, None] * SE[None, :]) - 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) - - self.q_calc = q - self._cosmat, self._cos0 = cosmat, cos0 \ No newline at end of file + qx, qy = np.meshgrid(q, q) + self._Iq_shape = qx.shape + self.q_calc = qx.flatten(), qy.flatten() + self._cosmat, self._cos0 = cosmat, cos0 From 2cdc35b9c3dcac265c15de2e96a382526a4f45b8 Mon Sep 17 00:00:00 2001 From: Paul Kienzle Date: Tue, 11 Apr 2017 11:44:18 -0400 Subject: [PATCH 4/4] provide working oriented/unoriented sesans examples --- example/oriented_sesans.py | 32 ++++++++++++++++++++++ example/unoriented_sesans.py | 26 ++++++++++++++++++ sasmodels/direct_model.py | 51 ++++++++++++++++-------------------- sasmodels/sesans.py | 9 +++---- 4 files changed, 83 insertions(+), 35 deletions(-) create mode 100644 example/oriented_sesans.py create mode 100644 example/unoriented_sesans.py diff --git a/example/oriented_sesans.py b/example/oriented_sesans.py new file mode 100644 index 000000000..e152a9347 --- /dev/null +++ b/example/oriented_sesans.py @@ -0,0 +1,32 @@ +from bumps.names import * + +from sasmodels.data import load_data +from sasmodels.core import load_model +from sasmodels.bumps_model import Model, Experiment + +# Spherical particle data, not ellipsoids +data = load_data('../../sasview/sasview/test/sesans_data/sphere2micron.ses') +data.oriented = True + +kernel = load_model("ellipsoid") + +model = Model( + kernel, + scale=0.08, background=0, + sld=.291, sld_solvent=7.105, + radius_polar=10000, radius_polar_pd=0.222296, radius_polar_pd_n=0, + radius_equatorial=2600, radius_equatorial_pd=0.28, radius_equatorial_pd_n=0, + theta=60, theta_pd=0, theta_pd_n=0, + phi=60, phi_pd=0, phi_pd_n=0, + ) + +# SET THE FITTING PARAMETERS +model.radius_polar.range(1000, 100000) +model.radius_equatorial.range(1000, 100000) +model.theta.range(0, 90) +model.phi.range(0, 180) +model.background.range(0,1000) +model.scale.range(0, 10) + +M = Experiment(data=data, model=model) +problem = FitProblem(M) diff --git a/example/unoriented_sesans.py b/example/unoriented_sesans.py new file mode 100644 index 000000000..a3b154f6b --- /dev/null +++ b/example/unoriented_sesans.py @@ -0,0 +1,26 @@ +from bumps.names import * + +from sasmodels.data import load_data +from sasmodels.core import load_model +from sasmodels.bumps_model import Model, Experiment + +# Spherical particle data, not ellipsoids +data = load_data('../../sasview/sasview/test/sesans_data/sphere2micron.ses') +kernel = load_model("sphere*hardsphere") + +model = Model( + kernel, + scale=0.08, background=0, + sld=.291, sld_solvent=7.105, + radius=10000, radius_pd=0.222296, radius_pd_n=0, + volfraction=0.2, + ) + +# SET THE FITTING PARAMETERS +model.radius.range(1000, 100000) +model.volfraction.range(0, 1) +model.background.range(0, 1000) +model.scale.range(0, 10) + +M = Experiment(data=data, model=model) +problem = FitProblem(M) diff --git a/sasmodels/direct_model.py b/sasmodels/direct_model.py index 386e852f9..289b90fc2 100644 --- a/sasmodels/direct_model.py +++ b/sasmodels/direct_model.py @@ -201,16 +201,23 @@ def _interpret_data(self, data, model): self.data_type = 'Iq' if self.data_type == 'sesans': - q = sesans.make_q(data.sample.zacceptance, data.Rmax) - index = slice(None, None) - res = None - if data.y is not None: - Iq, dIq = data.y, data.dy + from sas.sascalc.data_util.nxsunit import Converter + qmax, qunits = data.sample.zacceptance + SElength = Converter(data._xunit)(data.x, "A") + zaccept = Converter(qunits)(qmax, "1/A"), + Rmax = 10000000 + index = slice(None, None) # equivalent to index [:] + Iq = data.y[index] + dIq = data.dy[index] + oriented = getattr(data, 'oriented', False) + if oriented: + res = sesans.OrientedSesansTransform(data.x[index], SElength, zaccept, Rmax) + # Oriented sesans transform produces q_calc = [qx, qy] + q_vectors = res.q_calc else: - Iq, dIq = None, None - #self._theory = np.zeros_like(q) - q_vectors = [q] - q_mono = sesans.make_all_q(data) + res = sesans.SesansTransform(data.x[index], SElength, zaccept, Rmax) + # Unoriented sesans transform produces q_calc = q + q_vectors = [res.q_calc] elif self.data_type == 'Iqxy': #if not model.info.parameters.has_2d: # raise ValueError("not 2D without orientation or magnetic parameters") @@ -229,7 +236,6 @@ def _interpret_data(self, data, model): nsigma=3.0, accuracy=accuracy) #self._theory = np.zeros_like(self.Iq) q_vectors = res.q_calc - q_mono = [] elif self.data_type == 'Iq': index = (data.x >= data.qmin) & (data.x <= data.qmax) if data.y is not None: @@ -254,7 +260,6 @@ def _interpret_data(self, data, model): #self._theory = np.zeros_like(self.Iq) q_vectors = [res.q_calc] - q_mono = [] elif self.data_type == 'Iq-oriented': index = (data.x >= data.qmin) & (data.x <= data.qmax) if data.y is not None: @@ -271,14 +276,12 @@ def _interpret_data(self, data, model): qx_width=data.dxw[index], qy_width=data.dxl[index]) q_vectors = res.q_calc - q_mono = [] else: raise ValueError("Unknown data type") # never gets here # Remember function inputs so we can delay loading the function and # so we can save/restore state self._kernel_inputs = q_vectors - self._kernel_mono_inputs = q_mono self._kernel = None self.Iq, self.dIq, self.index = Iq, dIq, index self.resolution = res @@ -305,9 +308,6 @@ def _calc_theory(self, pars, cutoff=0.0): # type: (ParameterSet, float) -> np.ndarray if self._kernel is None: self._kernel = self._model.make_kernel(self._kernel_inputs) - self._kernel_mono = ( - self._model.make_kernel(self._kernel_mono_inputs) - if self._kernel_mono_inputs else None) Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff) # Storing the calculated Iq values so that they can be plotted. @@ -315,19 +315,12 @@ def _calc_theory(self, pars, cutoff=0.0): # TODO: extend plotting of calculate Iq to other measurement types # TODO: refactor so we don't store the result in the model self.Iq_calc = None - if self.data_type == 'sesans': - Iq_mono = (call_kernel(self._kernel_mono, pars, mono=True) - if self._kernel_mono_inputs else None) - result = sesans.transform(self._data, - self._kernel_inputs[0], Iq_calc, - self._kernel_mono_inputs, Iq_mono) - else: - result = self.resolution.apply(Iq_calc) - if hasattr(self.resolution, 'nx'): - self.Iq_calc = ( - self.resolution.qx_calc, self.resolution.qy_calc, - np.reshape(Iq_calc, (self.resolution.ny, self.resolution.nx)) - ) + result = self.resolution.apply(Iq_calc) + if hasattr(self.resolution, 'nx'): + self.Iq_calc = ( + self.resolution.qx_calc, self.resolution.qy_calc, + np.reshape(Iq_calc, (self.resolution.ny, self.resolution.nx)) + ) return result diff --git a/sasmodels/sesans.py b/sasmodels/sesans.py index 127079919..fee7f08ca 100644 --- a/sasmodels/sesans.py +++ b/sasmodels/sesans.py @@ -44,8 +44,6 @@ def __init__(self, z, SElength, zaccept, Rmax): # type: (np.ndarray, float, float) -> None #import logging; logging.info("creating SESANS transform") self.q = z - # isoriented flag determines whether data is from an oriented sample or - # not, should be a selection variable upon entering SESANS data. self._set_hankel(SElength, zaccept, Rmax) def apply(self, Iq): @@ -104,15 +102,13 @@ def __init__(self, z, SElength, zaccept, Rmax): # type: (np.ndarray, float, float) -> None #import logging; logging.info("creating SESANS transform") self.q = z - # isoriented flag determines whether data is from an oriented sample or - # not, should be a selection variable upon entering SESANS data. self._set_cosmat(SElength, zaccept, Rmax) def apply(self, Iq): dq = self.q_calc[0][0] Iq = np.reshape(Iq, self._Iq_shape) G0 = self._cos0 * np.sum(Iq) * dq - G = np.sum(np.dot(Iq, self._cosmat.T), axis=1) * dq + G = np.sum(np.dot(Iq, self._cosmat), axis=0) * dq P = G - G0 return P @@ -124,12 +120,13 @@ def _set_cosmat(self, SElength, zaccept, Rmax): # Rmax = #value in text box somewhere in FitPage? q_max = 2 * pi / (SElength[1] - SElength[0]) q_min = 0.1 * 2 * pi / (np.size(SElength) * SElength[-1]) + q_min *= 100 q = np.arange(q_min, q_max, q_min, dtype='float32') dq = q_min cos0 = np.float32(dq / (2 * pi)) - cosmat = np.float32(dq / (2 * pi)) * np.cos(q[:, None] * SE[None, :]) + cosmat = np.float32(dq / (2 * pi)) * np.cos(q[:, None] * SElength[None, :]) qx, qy = np.meshgrid(q, q) self._Iq_shape = qx.shape