# Changeset ae32bb8 in sasmodels

Ignore:
Timestamp:
Mar 2, 2017 4:32:31 AM (3 years ago)
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
37c7d5e
Parents:
1896dff (diff), e1aa129 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge branch 'master' of https://github.com/SasView/sasmodels

Files:
2 edited

### Legend:

Unmodified
 r7ae2b7f 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 scipy.special import j0 class SesansTransform(object): """ from sas.sascalc.data_util.nxsunit import Converter 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 q_min = dq = 0.1 * 2*pi / Rmax return np.arange(q_min, Converter(q_max)(q_max, units="1/A"), dq) def make_all_q(data): *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. """ 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] #: 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 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 result = call_HankelAccept(data, q_calc, Iq_calc, q, Iq_mono) else: Qx, Qy = [qmono, qmono] 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) #: q values to calculate when computing transform q_calc = None  # type: np.ndarray return result # transform arrays _H = None  # type: np.ndarray _H0 = None # type: np.ndarray 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 __init__(self, z, SElength, zaccept, Rmax): # type: (np.ndarray, float, float) -> None #import logging; logging.info("creating SESANS transform") self.q = z self._set_hankel(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 return P 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 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') 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-allq)*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 #Rmax = #value in text box somewhere in FitPage? q_max = 2*pi / (SElength - SElength) 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 sd[i] += Iq f[i] = 1-s[i]+sd[i] P[i] = (1-sd[i]/f[i])+1/f[i]*Gprime[i] 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 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-allq)*1e10 sigma[i]=wavelength[i]^2*thickness/2/pi*np.sum(allIq*allq*alldq) s[i]=1-exp(-sigma) dq = (q-q)*1e10 a = (x