diff --git a/direfl/api/invert.py b/direfl/api/invert.py index 479044a..cfaad53 100755 --- a/direfl/api/invert.py +++ b/direfl/api/invert.py @@ -876,7 +876,7 @@ def plotamp(Q, r, dr=None, scaled=True, ylabel="Real R", **kw): pylab.fill_between(Q, scale*(r-dr), scale*(r+dr), color=h.get_color(), alpha=0.2) pylab.ylabel(ylabel) - pylab.xlabel("Q (inv A)") + pylab.xlabel("Q $[\AA^{-1}]$") class Interpolator(): @@ -908,8 +908,10 @@ def remesh(data, xmin, xmax, npts, left=None, right=None): x, y = data x, y = x[isfinite(x)], y[isfinite(y)] + if npts > len(x): npts = len(x) + newx = np.linspace(xmin, xmax, npts) newy = np.interp(newx, x, y, left=left, right=right) return np.array((newx, newy)) @@ -1533,6 +1535,9 @@ def _phase_reconstruction(Q, R1sq, R2sq, rho_u, rho_v1, rho_v2): Compute phase reconstruction from back reflectivity on paired samples with varying surface materials. + "Fixed Nonvacuum Fronting, Variable Backing" + Uses eq. (31), (32) from [Majkrzak2003]. + Inputs:: *Q* is the measurement positions @@ -1543,6 +1548,15 @@ def _phase_reconstruction(Q, R1sq, R2sq, rho_u, rho_v1, rho_v2): Returns RealR, ImagR """ + # The used notation here is different from the paper [Majkrzak2003]. + # To more easily understand the code, take a look at the following translation table + # + # Paper | Code + # f^2 = usq + # f^2(a^2 + f^2b^2) = alpha + # f^2(d^2 + c^2) = beta + # \Sigma^{fh_i} = sigmai with i = 1, 2 + # h_1^2, h_1^2 = v1sq, v2sq Qsq = Q**2 + 16.*pi*rho_u*1e-6 usq, v1sq, v2sq = [(1-16*pi*rho*1e-6/Qsq) for rho in (rho_u, rho_v1, rho_v2)] diff --git a/direfl/api/reference_layer.py b/direfl/api/reference_layer.py index 61119c0..8f8c17f 100644 --- a/direfl/api/reference_layer.py +++ b/direfl/api/reference_layer.py @@ -1,43 +1,69 @@ -from math import pi import cmath - import numpy as np -from .util import isstr -from .invert import SurroundVariation, remesh +from math import pi + +from .invert import SurroundVariation, reduce from .sld_profile import SLDProfile, refr_idx +try: # CRUFT: basestring isn't used in python3 + basestring +except: + basestring = str + +ZERO_TOL = 1e-10 +# The tolerance to decide, when the reflectivity is 1, i.e. |r - 1| < tol. +REFLECTIVITY_UNITY_TOL = 1e-10 """ References: + + [Majkrzak2003] C. F. Majkrzak, N. F. Berk and U. A. Perez-Salas. Langmuir (2003), 19, 7796-7810. + Phase-Sensitive Neutron Reflectometry - [Majkrzak2003] C. F. Majkrzak, N. F. Berk: Physica B 336 (2003) 27-38 - Phase sensitive reflectometry and the unambiguous determination - of scattering length density profiles """ class AbstractReferenceVariation(SurroundVariation): def __init__(self, fronting_sld, backing_sld): + """ + fronting and backing sld are in units of [AA^-2]. + + Example: + For a backing SLD of Si, use backing_sld = 2.1e-6 + + :param fronting_sld: + :param backing_sld: + """ self._f = float(fronting_sld) self._b = float(backing_sld) self._measurements = [] - # The tolerance to decide, when the reflectivity is 1, i.e. |r - 1| < tol. - self.REFLECTIVITY_UNITY_TOL = 1e-10 - self.ZERO_TOL = 1e-10 - self.MATRIX_ILL_CONDITIONED = 1e10 - self.dImagR = None self.dRealR = None self.plot_imaginary_offset = -5 def run(self): + """ + Runs the phase reconstruction algorithm on the loaded data set. + + The result (Real and Imaginary part of the reflection) can be access by the attributes + self.Q, self.RealR, self.ImagR + + :return: None + """ self._check_measurements() self._calc() self.Qin = self.Q def _check_measurements(self): + """ + Checks that the given measurements coincide on the q-grid. + Checks that the given SLDs are not equal (breaks the reconstruction algorithm) + + :raise: RuntimeWarning if the checks fail + :return: None + """ if len(self._measurements) <= 1: return @@ -45,54 +71,128 @@ def _check_measurements(self): for ms in self._measurements: if not ms['Qin'].shape == q1.shape or not all(q1 == ms['Qin']): - raise ValueError("Q points do not match in data files") + raise RuntimeWarning("Q points do not match in data files") slds = [ms['sld'] for ms in self._measurements] - # Checkts that there are no duplicate slds added + # Checks that there are no duplicate SLDs added # duplicate sld profile yield a singular matrix in the constraint system (obviously) if any([x == y for i, x in enumerate(slds) for j, y in enumerate(slds) if i != j]): raise RuntimeWarning("Two equal sld profiles found. The profiles have to be " "different.") - def remesh(self): + def remesh(self, interpolation=1, interpolation_kind=None): + """ + Re-meshes the loaded data onto a common grid by interpolation. + + Usually, the reflectivity data is not measured at the very same q points. Also the min/max range of + the q values might vary. But, to reconstruct the phase, the reflectivity needs to be measured at + the same q values. This method achieves this goal. - qmin, qmax, npts = [], [], [] - for measurement in self._measurements: - qmin.append(measurement['Qin'][0]) - qmax.append(measurement['Qin'][-1]) - npts.append(len(measurement['Qin'])) + The new grid is the coarsest possible grid. The min/max values are chosen such that every reflectivity + measurement contains the min/max values. - qmin = max(qmin) - qmax = min(qmax) - npts = min(npts) + :param interpolation: integer number. Defines the number of additional interpolations between + two q-grid points + :param interpolation_kind: interpolation of the function between the point (linear/quadratic/etc..) + See scipy.interp1d for possible values + :return: None + """ + + # coarsest possible grid + qmin = max([ms['Qin'][0] for ms in self._measurements]) + qmax = min([ms['Qin'][-1] for ms in self._measurements]) + npts = min([len(ms['Qin']) for ms in self._measurements]) + + new_mesh = np.linspace(qmin, qmax, npts + 1) for measurement in self._measurements: q, R, dR = measurement['Qin'], measurement['Rin'], measurement['dRin'] + try: + from skipi.function import Function + f = Function.to_function(q, R, interpolation=interpolation_kind).remesh(new_mesh).oversample( + interpolation) - q, R = remesh([q, R], qmin, qmax, npts, left=0, right=0) - if dR is not None: - q, dR = remesh([q, dR], qmin, qmax, npts, left=0, right=0) + if dR is not None: + df = Function.to_function(q, dR, interpolation=interpolation_kind).remesh( + new_mesh).oversample(interpolation) + dR = df.eval() - measurement['Qin'], measurement['Rin'], measurement['dRin'] = q, R, dR + measurement['Qin'], measurement['Rin'], measurement['dRin'] = f.get_domain(), f.eval(), dR + except: + # Fallback if skipi is not available + q, R = remesh([q, R], qmin, qmax, npts, left=0, right=0) + if dR is not None: + q, dR = remesh([q, dR], qmin, qmax, npts, left=0, right=0) + + measurement['Qin'], measurement['Rin'], measurement['dRin'] = q, R, dR def load_data(self, q, R, dq, dR, sld_profile, name='unknown'): + """ + Load data directly + + :param q: array of q + :param R: array of R(q) + :param dq: array of error in q or None + :param dR: array of error in R(q) or None + :param sld_profile: reference layer as SLDProfile + :param name: name of the measurement + :return: None + """ assert isinstance(sld_profile, SLDProfile) self._measurements.append( {'name': name, 'Qin': q, 'dQin': dq, 'Rin': R, 'dRin': dR, 'sld': sld_profile}) - self.number_measurements = len(self._measurements) + def load_function(self, function, sld_profile, name='unknown'): + """ + Load the reflectivity data by using the skipi function package. + + :param function: reflectivity measurement, including errors in dx, dy + :param sld_profile: reference layer as SLDProfile + :param name: name of the measurement + :return: None + """ + from skipi.function import Function - def load(self, file, sld_profile, use_columns=None): + assert isinstance(sld_profile, SLDProfile) + assert isinstance(function, Function) + + self._measurements.append({ + 'name': name, + 'Qin': function.get_domain(), + 'dQin': function.dx.eval(), + 'Rin': function.eval(), + 'dRin': function.dy.eval(), + 'sld': sld_profile + }) + + def load(self, file, sld_profile, use_columns=None, name=None, q0=0): + """ + Load the reflectivity data by reading from a file. + + File structure: + q, R(q), [dR(q)] for 2-3 column files + q, dq, R(q), dR(q), [wavelength] for 4-5 column files + + :param file: file name + :param sld_profile: reference layer as SLDProfile + :param use_columns: columns to read the data from + :param name: optional name of the measurement, if None then filename is used + :param q0: minimum q value to consider, all measurements below q0 are discarded + :return: None + """ assert isinstance(sld_profile, SLDProfile) - if isstr(file): + if isinstance(file, basestring): d = np.loadtxt(file, usecols=use_columns).T - name = file + _name = file else: d = file - name = "SimData{}".format(len(self._measurements) + 1) + _name = "Measurement {}".format(len(self._measurements) + 1) + + if name is None: + name = _name q, dq, r, dr = None, None, None, None @@ -110,14 +210,20 @@ def load(self, file, sld_profile, use_columns=None): elif ncols >= 5: q, dq, r, dr, lamb = d[0:5] + dq = dq[q > q0] if dq is not None else None + dr = dr[q > q0] if dr is not None else np.zeros(len(r)) + r = r[q > q0] + q = q[q > q0] + self._measurements.append( {'name': name, 'Qin': q, 'dQin': dq, 'Rin': r, 'dRin': dr, 'sld': sld_profile}) - self.number_measurements = len(self._measurements) - def _calc(self): - self.Q, self.Rall = self._phase_reconstruction() + self.Q, self.Rall, self.dR = self._phase_reconstruction() l = len(self.Rall) + + # If only two measurements are supplied, Rall contains then + # two branches of the reflection. This code splits the tuple up into R+ and R- self.Rp, self.Rm = np.zeros(l, dtype=complex), np.zeros(l, dtype=complex) for idx, el in enumerate(self.Rall): if type(el) is np.ndarray or type(el) is list: @@ -131,22 +237,72 @@ def _calc(self): self.R = self.Rp self.RealR, self.ImagR = self.R.real, self.R.imag - def _refl(self, alpha_u, beta_u, gamma_u): + @classmethod + def _refl(cls, alpha_u, beta_u, gamma_u): # Compute the reflection coefficient, based on the knowledge of alpha_u, beta_u, # gamma_u where these parameters are the solution of the matrix equation # # See eq (38)-(40) in [Majkrzak2003] return - (alpha_u - beta_u + 2 * 1j * gamma_u) / (alpha_u + beta_u + 2) - def _calc_refl_constraint(self, q, reflectivity, sld_reference, fronting, - backing): + @classmethod + def _drefl(cls, alpha, beta, gamma, cov): + r""" + Estimates the error within the reconstructed phase information (real and imaginary part are + calculated separately) by the propagation of error method of the formula: + + ..math:: + R = (\beta - \alpha - 2i\gamma) / (\alpha + \beta + 2) + + :param alpha: alpha_u: solution of the linear equation + :param beta: beta_u: solution of the linear equation + :param gamma: gamma_u: solution of the linear equation + :param cov: Covariance matrix of the linear least squares method, with the order + [alpha, beta, gamma] = [0, 1, 2] + :return: Estimated error as complex number + """ + a, b, g = alpha, beta, gamma + + # dR/d alpha, real part + dAlpha = - 2 * (b + 1) / (a + b + 2) ** 2 + # dR/d beta, real part + dBeta = 2 * (a + 1) / (a + b + 2) ** 2 + # dR/d gamma = 0 (real part is independent of gamma) + + # dR/d alpha, imag part + dAlphaIm = -2 * g / (a + b + 2) ** 2 + # dR/d beta, imag part + dBetaIm = -2 * g / (a + b + 2) ** 2 + # dR/d gamma, imag part + dGammaIm = -2 / (a + b + 2) + + r""" + Error propagation + + ..math:: + Re(dR)^2 = (\frac{dR}{d\alpha} \sigma_\alpha)^2 + (\frac{dR}{d\beta} \sigma_\beta)^2 + + 2 * \frac{dR}{d\alpha}{\frac{dR}{d\beta} \sigma_{\alpha, \beta} + + where :math:`\sigma_\alpha` denotes the error within :math:`\alpha` and :math:`\sigma_{\alpha, \beta}` + denotes the covariance of :math:`\alpha' and :math:`\beta` + + The error for the imaginary part is computed analogously. + """ + sResqr = dAlpha ** 2 * cov[0][0] + dBeta ** 2 * cov[1][1] + 2 * dAlpha * dBeta * cov[0][1] + sImsqr = dAlphaIm ** 2 * cov[0][0] + dBetaIm ** 2 * cov[1][1] + dGammaIm ** 2 * cov[2][2] + \ + 2 * dAlphaIm * dBetaIm * cov[0][1] + \ + 2 * dAlphaIm * dGammaIm * cov[0][2] + \ + 2 * dBetaIm * dGammaIm * cov[1][2] + + return np.sqrt(sResqr) + 1j * np.sqrt(sImsqr) + + @classmethod + def _calc_refl_constraint(cls, q, reflectivity, sld_reference, fronting, backing): # See the child classes for implementations raise NotImplementedError() - def _phase_reconstruction(self): + def _do_phase_reconstruction(self, q, Rs, dRs, SLDs): """ - Here, we reconstruct the reflection coefficients for every q. - The calculation is split up in multiple parts (to keep the code repetition low). First, we calculate the constraining linear equations for the reflection coefficient. Only this depends on the location of the reference layer (front or back). Next, @@ -155,53 +311,86 @@ def _phase_reconstruction(self): reflection coefficients). Using the solution of the linear system, we finally calculate the reflection coefficient. - :return: q, r(q) for each q - """ - qs = [] - r = [] - - for idx, q in enumerate(self._measurements[0]['Qin']): - - # TODO: check this - q = cmath.sqrt(q ** 2 + 16.0 * pi * self._f).real - - # Skip those q values which are too close to zero, this would break the - # refractive index calculation otherwise - if abs(q) < self.ZERO_TOL: - continue + Note that this reconstructs the reflection and also returns a new q value since this might have + changed do to a non-zero fronting medium. - f = refr_idx(q, self._f) - b = refr_idx(q, self._b) + :param q: The q value + :param Rs: The reflections measured at q + :param dRs: The uncertainties in the measured reflections at q + :param SLDs: The SLDs corresponding to the reflections + :return: q_new, Reflection, Error in Reflection + """ + # Shift the q vector if the incidence medium is not vacuum + # See eq (49) in [Majkrzak2003] + # + # Note that this also prohibits to measure the phase information + # below the critical edge by simply flipping the sample. + # The first q value you effectively measure is the first + # one direct _after_ the critical edge .. + q = cmath.sqrt(q ** 2 + 16.0 * pi * self._f).real + + # Skip those q values which are too close to zero, this would break the + # refractive index calculation otherwise + if abs(q) < ZERO_TOL: + return None + + f = refr_idx(q, self._f) + b = refr_idx(q, self._b) + + A = [] + c = [] + # Calculate for each measurement a linear constraint. Putting all of the + # constraints together enables us to solve for the reflection itself. How to + # calculate the linear constraint using a reference layer can be + # found in [Majkrzak2003] + + for R, dR, SLD in zip(Rs, dRs, SLDs): + # Don't use values close to the total refection regime. + # You can't reconstruct the reflection below there with this method. + if abs(R - 1) < REFLECTIVITY_UNITY_TOL: + return None + + lhs, rhs, drhs = self._calc_refl_constraint(q, R, SLD, f, b) + + sigma = 1e-10 + + # Note: the right hand side is a function of R and thus, the std deviation of the rhs is + # simply the derivative times the std deviation of R + if abs(dR) > ZERO_TOL: + sigma = drhs * dR + + # divide by sigma, so that we do a chi squared minimization. + A.append(np.array(lhs) / sigma) + c.append(rhs / sigma) + + try: + R, dR = self._solve_reference_layer(A, c) + return q, R, dR + except RuntimeWarning as e: + print("Could not reconstruct the phase for q = {}. Reason: {}".format(q, e)) - A = [] - c = [] - # Calculate for each measurement a linear constraint. Putting all of the - # constraints together enables us to solve for the reflection itself. How to - # calculate the linear constraint using a reference layer can be - # found in [Majkrzak2003] - for ms in self._measurements: + def _phase_reconstruction(self): + """ + Here, we reconstruct the reflection coefficients for every q. - # Don't use values close to the total refection regime. - # You can't reconstruct the reflection below there with this method. - if abs(ms['Rin'][idx] - 1) < self.REFLECTIVITY_UNITY_TOL: - continue + :return: q, r(q), dr(q) for each q + """ + qr = np.empty(len(self._measurements[0]['Qin']), dtype=tuple) - lhs, rhs = self._calc_refl_constraint(q, ms['Rin'][idx], ms['sld'], f, b) + SLDs = [ms['sld'] for ms in self._measurements] - A.append(lhs) - c.append(rhs) + for idx, q in enumerate(self._measurements[0]['Qin']): + Rs = [ms['Rin'][idx] for ms in self._measurements] + dRs = [ms['dRin'][idx] for ms in self._measurements] - try: - reflection = self._solve_reference_layer(A, c) + qr[idx] = self._do_phase_reconstruction(q, Rs, dRs, SLDs) - qs.append(q) - r.append(reflection) - except RuntimeWarning as e: - print("Could not reconstruct the phase for q = {}. Reason: {}".format(q, e)) + qs, rs, dRs = zip(*qr[qr != None]) - return np.array(qs), np.array(r) + return np.array(qs), np.array(rs), np.array(dRs) - def _solve_reference_layer(self, A, c): + @classmethod + def _solve_reference_layer(cls, A, c): """ Solving the linear system A x = c with x = [alpha_u, beta_u, gamma_u], being the unknown coefficients for the @@ -221,10 +410,9 @@ def _solve_reference_layer(self, A, c): The condition gamma^2 = alpha * beta - 1 will be used to construct two reflection coefficients which solve the equation. A list of two reflection coefficients is returned then. - N == 3: - A usual linear inversion is performed - N >= 4: + N >= 3: A least squares fit is performed (A^T A x = c) + (which is exact for N=3) If any of the operations is not possible, (bad matrix condition number, quadratic eq has no real solution) a RuntimeWarning exception is raised @@ -264,12 +452,12 @@ def _solve_reference_layer(self, A, c): # alpha = u - v * gamma, see above, the linear relationship alpha_beta = lambda u, v, g: u - v * g - if abs(det) < self.ZERO_TOL: + if abs(det) < ZERO_TOL: # Luckily, we get just one solution for gamma :D gamma_u = -b / (2 * a) alpha_u = alpha_beta(u[0], v[0], gamma_u) beta_u = alpha_beta(u[1], v[1], gamma_u) - return self._refl(alpha_u, beta_u, gamma_u) + return cls._refl(alpha_u, beta_u, gamma_u), 0 elif det > 0: reflection = [] # Compute first gamma using both branches of the quadratic solution @@ -279,31 +467,25 @@ def _solve_reference_layer(self, A, c): gamma_u = (-b + sign * cmath.sqrt(det).real) / (2 * a) alpha_u = alpha_beta(u[0], v[0], gamma_u) beta_u = alpha_beta(u[1], v[1], gamma_u) - reflection.append(self._refl(alpha_u, beta_u, gamma_u)) + reflection.append(cls._refl(alpha_u, beta_u, gamma_u)) # Returns the reflection branches, R+ and R- - return reflection + return reflection, 0 else: - # This usually happens is the reference sld's are not correct. + # This usually happens if the reference sld's are not correct. raise RuntimeWarning("The quadratic equation has no real solution.") - if len(A) == 3: - # Highly ill-conditioned, better throw away the solution than pretending it's - # good ... - # TODO: maybe least squares? - condition_number = np.linalg.cond(A) - if condition_number > self.MATRIX_ILL_CONDITIONED: - raise RuntimeWarning("Given linear constraints are ill conditioned. " - "Condition number {}".format(condition_number)) - alpha_u, beta_u, gamma_u = np.linalg.solve(A, c) - return self._refl(alpha_u, beta_u, gamma_u) - - if len(A) > 3: + if len(A) >= 3: + # least squares solves exact for 3x3 matrices # Silence the FutureWarning with rcond=None solution, residuals, rank, singular_values = np.linalg.lstsq(A, c, rcond=None) alpha_u, beta_u, gamma_u = solution - return self._refl(alpha_u, beta_u, gamma_u) + + # covariance matrix + C = np.linalg.inv(np.array(A).T.dot(A)) + + return cls._refl(alpha_u, beta_u, gamma_u), cls._drefl(alpha_u, beta_u, gamma_u, C) def choose(self, plus_or_minus): """ @@ -410,8 +592,8 @@ def plot_r_choose(self, branch_selection=1, plot_jumps_continuity=True, class BottomReferenceVariation(AbstractReferenceVariation): - - def _calc_refl_constraint(self, q, reflectivity, sld_reference, fronting, + @classmethod + def _calc_refl_constraint(cls, q, reflectivity, sld_reference, fronting, backing): """ Solving the linear system A x = c @@ -422,6 +604,7 @@ def _calc_refl_constraint(self, q, reflectivity, sld_reference, fronting, of the equation (38) in [Majkrzak2003] This method returns one row in this linear system as: lhs, rhs with lhs being one row in the matrix A and rhs being one scalar in the vector b + drhs is the "variance" of the rhs, i.e. just the first derivative of rhs """ w, x, y, z = sld_reference.as_matrix(q) f, b = fronting, backing @@ -432,12 +615,14 @@ def _calc_refl_constraint(self, q, reflectivity, sld_reference, fronting, lhs = [f ** 2 * beta, b ** 2 * alpha, 2 * f * b * gamma] rhs = 2 * f * b * (1 + reflectivity) / (1 - reflectivity) - return lhs, rhs + drhs = 4 * f * b * 1 / (1 - reflectivity) ** 2 + return lhs, rhs, drhs -class TopReferenceVariation(AbstractReferenceVariation): - def _calc_refl_constraint(self, q, reflectivity, sld_reference, fronting, +class TopReferenceVariation(AbstractReferenceVariation): + @classmethod + def _calc_refl_constraint(cls, q, reflectivity, sld_reference, fronting, backing): """ Solving the linear system A x = c @@ -448,6 +633,7 @@ def _calc_refl_constraint(self, q, reflectivity, sld_reference, fronting, of the equation (33) in [Majkrzak2003] This method returns one row in this linear system as: lhs, rhs with lhs being one row in the matrix A and rhs being one scalar in the vector b + drhs is the "variance" of the rhs, i.e. just the first derivative of rhs """ w, x, y, z = sld_reference.as_matrix(q) f, b = fronting, backing @@ -458,5 +644,6 @@ def _calc_refl_constraint(self, q, reflectivity, sld_reference, fronting, lhs = [b ** 2 * beta, f ** 2 * alpha, 2 * f * b * gamma] rhs = 2 * f * b * (1 + reflectivity) / (1 - reflectivity) + drhs = 4 * f * b * 1 / (1 - reflectivity) ** 2 - return lhs, rhs + return lhs, rhs, drhs diff --git a/direfl/api/sld_profile.py b/direfl/api/sld_profile.py index 56a1eae..a8f8a36 100644 --- a/direfl/api/sld_profile.py +++ b/direfl/api/sld_profile.py @@ -1,7 +1,9 @@ +import numpy as np import cmath + +from functools import reduce from math import pi, ceil -import numpy as np from numpy import sin, cos from scipy.interpolate import interp1d @@ -66,6 +68,9 @@ def as_matrix(self, q): class ConstantSLDProfile(SLDProfile): def __init__(self, sld, thickness, sigma=0): + if sld > 15: + raise RuntimeError("SLD seems to be unreasonable high") + self._sld = float(sld) self._d = float(thickness) self._r = float(sigma) @@ -141,13 +146,14 @@ def from_sample(cls, sample, dz=0.1, dA=1e-4, probe=None): @classmethod def from_slabs(cls, thickness, sld, roughness, precision=1): - # You should rather use the from_sample method, since this easier to + # You should rather use the from_sample method, since its easier to # understand. This method here is just a kind of 'fallback' - # if you dont wanna have the overhead of building the Stacks in refl1d + # if you don't wanna have the overhead of building the Stacks in refl1d # just to put the data in here.. - + # # WARNING: from_slabs and from_sample do not create the same slab profile - # they are shifted profiles (by I'd guess 3*roughness[0]) + # they are shifted profiles (by I'd guess 3*roughness[0]?) + from refl1d.profile import build_profile w = thickness @@ -186,7 +192,7 @@ def plot_profile(self, offset=0, reverse=False): pylab.plot(self._z + offset, rho) def as_matrix(self, q): - from functools import reduce + # len(dz) = len(self._z) - 1 dz = np.diff(self._z) m = len(dz) * [None]