""" This module implements invariant and its related computations. @author: Gervaise B. Alina/UTK TODO: - intro / documentation - add unit tests for sufrace/volume computation with and without extrapolation. - replace the get_extra_data_* methods """ import math import numpy from DataLoader.data_info import Data1D as LoaderData1D # The minimum q-value to be used when extrapolating Q_MINIMUM = 1e-5 # The maximum q-value to be used when extrapolating Q_MAXIMUM = 10 # Number of steps in the extrapolation INTEGRATION_NSTEPS = 1000 class Transform(object): """ Define interface that need to compute a function or an inverse function given some x, y """ def linearize_data(self, data): """ Linearize data so that a linear fit can be performed. Filter out the data that can't be transformed. @param data : LoadData1D instance """ # Check that the vector lengths are equal assert(len(data.x)==len(data.y)) if data.dy is not None: assert(len(data.x)==len(data.dy)) dy = data.dy else: dy = numpy.ones(len(data.y)) # Transform the data data_points = zip(data.x, data.y, dy) output_points = [(self.linearize_q_value(p[0]), math.log(p[1]), p[2]/p[1]) for p in data_points if p[0]>0 and p[1]>0 and p[2]>0] x_out, y_out, dy_out = zip(*output_points) # Create Data1D object x_out = numpy.asarray(x_out) y_out = numpy.asarray(y_out) dy_out = numpy.asarray(dy_out) linear_data = LoaderData1D(x=x_out, y=y_out, dy=dy_out) return linear_data def get_allowed_bins(self, data): """ Goes through the data points and returns a list of boolean values to indicate whether each points is allowed by the model or not. @param data: Data1D object """ return [p[0]>0 and p[1]>0 and p[2]>0 for p in zip(data.x, data.y, data.dy)] def linearize_q_value(self, value): """ Transform the input q-value for linearization """ return NotImplemented def extract_model_parameters(self, constant, slope, dconstant=0, dslope=0): """ set private member """ return NotImplemented def evaluate_model(self, x): """ Returns an array f(x) values where f is the Transform function. """ return NotImplemented def evaluate_model_errors(self, x): """ Returns an array of I(q) errors """ return NotImplemented class Guinier(Transform): """ class of type Transform that performs operations related to guinier function """ def __init__(self, scale=1, radius=60): Transform.__init__(self) self.scale = scale self.radius = radius ## Uncertainty of scale parameter self.dscale = 0 ## Unvertainty of radius parameter self.dradius = 0 def linearize_q_value(self, value): """ Transform the input q-value for linearization @param value: q-value @return: q*q """ return value * value def extract_model_parameters(self, constant, slope, dconstant=0, dslope=0): """ assign new value to the scale and the radius """ self.scale = math.exp(constant) self.radius = math.sqrt(-3 * slope) # Errors self.dscale = math.exp(constant)*dconstant self.dradius = -3.0/2.0/math.sqrt(-3 * slope)*dslope def evaluate_model(self, x): """ return F(x)= scale* e-((radius*x)**2/3) """ return self._guinier(x) def evaluate_model_errors(self, x): """ Returns the error on I(q) for the given array of q-values @param x: array of q-values """ p1 = numpy.array([self.dscale * math.exp(-((self.radius * q)**2/3)) for q in x]) p2 = numpy.array([self.scale * math.exp(-((self.radius * q)**2/3)) * (-(q**2/3)) * 2 * self.radius * self.dradius for q in x]) diq2 = p1*p1 + p2*p2 return numpy.array( [math.sqrt(err) for err in diq2] ) def _guinier(self, x): """ Retrive the guinier function after apply an inverse guinier function to x Compute a F(x) = scale* e-((radius*x)**2/3). @param x: a vector of q values @param scale: the scale value @param radius: the guinier radius value @return F(x) """ # transform the radius of coming from the inverse guinier function to a # a radius of a guinier function if self.radius <= 0: raise ValueError("Rg expected positive value, but got %s"%self.radius) value = numpy.array([math.exp(-((self.radius * i)**2/3)) for i in x ]) return self.scale * value class PowerLaw(Transform): """ class of type transform that perform operation related to power_law function """ def __init__(self, scale=1, power=4): Transform.__init__(self) self.scale = scale self.power = power def linearize_q_value(self, value): """ Transform the input q-value for linearization @param value: q-value @return: log(q) """ return math.log(value) def extract_model_parameters(self, constant, slope, dconstant=0, dslope=0): """ Assign new value to the scale and the power """ self.power = -slope self.scale = math.exp(constant) # Errors self.dscale = math.exp(constant)*dconstant self.dradius = -dslope def evaluate_model(self, x): """ given a scale and a radius transform x, y using a power_law function """ return self._power_law(x) def evaluate_model_errors(self, x): """ Returns the error on I(q) for the given array of q-values @param x: array of q-values """ p1 = numpy.array([self.dscale * math.pow(q, -self.power) for q in x]) p2 = numpy.array([self.scale * self.power * math.pow(q, -self.power-1) * self.dradius for q in x]) diq2 = p1*p1 + p2*p2 return numpy.array( [math.sqrt(err) for err in diq2] ) def _power_law(self, x): """ F(x) = scale* (x)^(-power) when power= 4. the model is porod else power_law The model has three parameters: @param x: a vector of q values @param power: power of the function @param scale : scale factor value @param F(x) """ if self.power <= 0: raise ValueError("Power_law function expected positive power, but got %s"%self.power) if self.scale <= 0: raise ValueError("scale expected positive value, but got %s"%self.scale) value = numpy.array([ math.pow(i, -self.power) for i in x ]) return self.scale * value class Extrapolator: """ Extrapolate I(q) distribution using a given model """ def __init__(self, data, model=None): """ Determine a and b given a linear equation y = ax + b If a model is given, it will be used to linearize the data before the extrapolation is performed. If None, a simple linear fit will be done. @param data: data containing x and y such as y = ax + b @param model: optional Transform object """ self.data = data self.model = model # Set qmin as the lowest non-zero value self.qmin = Q_MINIMUM for q_value in self.data.x: if q_value > 0: self.qmin = q_value break self.qmax = max(self.data.x) def fit(self, power=None, qmin=None, qmax=None): """ Fit data for y = ax + b return a and b @param power: a fixed, otherwise None @param qmin: Minimum Q-value @param qmax: Maximum Q-value """ if qmin is None: qmin = self.qmin if qmax is None: qmax = self.qmax # Identify the bin range for the fit idx = (self.data.x >= qmin) & (self.data.x <= qmax) fx = numpy.zeros(len(self.data.x)) # Uncertainty if type(self.data.dy)==numpy.ndarray and len(self.data.dy)==len(self.data.x): sigma = self.data.dy else: sigma = numpy.ones(len(self.data.x)) # Compute theory data f(x) fx[idx] = self.data.y[idx] # Linearize the data if self.model is not None: linearized_data = self.model.linearize_data(LoaderData1D(self.data.x[idx], fx[idx], dy = sigma[idx])) else: linearized_data = LoaderData1D(self.data.x[idx], fx[idx], dy = sigma[idx]) ##power is given only for function = power_law if power != None: sigma2 = linearized_data.dy * linearized_data.dy a = -(power) b = (numpy.sum(linearized_data.y/sigma2) \ - a*numpy.sum(linearized_data.x/sigma2))/numpy.sum(1.0/sigma2) deltas = linearized_data.x*a+numpy.ones(len(linearized_data.x))*b-linearized_data.y residuals = numpy.sum(deltas*deltas/sigma2) err = math.fabs(residuals) / numpy.sum(1.0/sigma2) return [a, b], [0, math.sqrt(err)] else: A = numpy.vstack([ linearized_data.x/linearized_data.dy, 1.0/linearized_data.dy]).T (p, residuals, rank, s) = numpy.linalg.lstsq(A, linearized_data.y/linearized_data.dy) # Get the covariance matrix, defined as inv_cov = a_transposed * a err = numpy.zeros(2) try: inv_cov = numpy.dot(A.transpose(), A) cov = numpy.linalg.pinv(inv_cov) err_matrix = math.fabs(residuals) * cov err = [math.sqrt(err_matrix[0][0]), math.sqrt(err_matrix[1][1])] except: err = [-1.0, -1.0] return p, err class InvariantCalculator(object): """ Compute invariant if data is given. Can provide volume fraction and surface area if the user provides Porod constant and contrast values. @precondition: the user must send a data of type DataLoader.Data1D the user provide background and scale values. @note: Some computations depends on each others. """ def __init__(self, data, background=0, scale=1 ): """ Initialize variables @param data: data must be of type DataLoader.Data1D @param background: Background value. The data will be corrected before processing @param scale: Scaling factor for I(q). The data will be corrected before processing """ # Background and scale should be private data member if the only way to # change them are by instantiating a new object. self._background = background self._scale = scale # The data should be private self._data = self._get_data(data) # Since there are multiple variants of Q*, you should force the # user to use the get method and keep Q* a private data member self._qstar = None # You should keep the error on Q* so you can reuse it without # recomputing the whole thing. self._qstar_err = 0 # Extrapolation parameters self._low_extrapolation_npts = 4 self._low_extrapolation_function = Guinier() self._low_extrapolation_power = None self._high_extrapolation_npts = 4 self._high_extrapolation_function = PowerLaw() self._high_extrapolation_power = None def _get_data(self, data): """ @note this function must be call before computing any type of invariant @return data= self._scale *data - self._background """ if not issubclass(data.__class__, LoaderData1D): #Process only data that inherited from DataLoader.Data_info.Data1D raise ValueError,"Data must be of type DataLoader.Data1D" new_data = (self._scale * data) - self._background # Check that the vector lengths are equal assert(len(new_data.x)==len(new_data.y)) # Verify that the errors are set correctly if new_data.dy is None or len(new_data.x) != len(new_data.dy) or \ (min(new_data.dy)==0 and max(new_data.dy)==0): new_data.dy = numpy.ones(len(new_data.x)) return new_data def _fit(self, model, qmin=Q_MINIMUM, qmax=Q_MAXIMUM, power=None): """ fit data with function using data= self._get_data() fx= Functor(data , function) y = data.y slope, constant = linalg.lstsq(y,fx) @param qmin: data first q value to consider during the fit @param qmax: data last q value to consider during the fit @param power : power value to consider for power-law @param function: the function to use during the fit @return a: the scale of the function @return b: the other parameter of the function for guinier will be radius for power_law will be the power value """ extrapolator = Extrapolator(data=self._data, model=model) p, dp = extrapolator.fit(power=power, qmin=qmin, qmax=qmax) return model.extract_model_parameters(constant=p[1], slope=p[0], dconstant=dp[1], dslope=dp[0]) def _get_qstar(self, data): """ Compute invariant for pinhole data. This invariant is given by: q_star = x0**2 *y0 *dx0 +x1**2 *y1 *dx1 + ..+ xn**2 *yn *dxn where n >= len(data.x)-1 dxi = 1/2*(xi+1 - xi) + (xi - xi-1) dx0 = (x1 - x0)/2 dxn = (xn - xn-1)/2 @param data: the data to use to compute invariant. @return q_star: invariant value for pinhole data. q_star > 0 """ if len(data.x) <= 1 or len(data.y) <= 1 or len(data.x)!= len(data.y): msg = "Length x and y must be equal" msg += " and greater than 1; got x=%s, y=%s"%(len(data.x), len(data.y)) raise ValueError, msg else: n = len(data.x)- 1 #compute the first delta q dx0 = (data.x[1] - data.x[0])/2 #compute the last delta q dxn = (data.x[n] - data.x[n-1])/2 sum = 0 sum += data.x[0] * data.x[0] * data.y[0] * dx0 sum += data.x[n] * data.x[n] * data.y[n] * dxn if len(data.x) == 2: return sum else: #iterate between for element different from the first and the last for i in xrange(1, n-1): dxi = (data.x[i+1] - data.x[i-1])/2 sum += data.x[i] * data.x[i] * data.y[i] * dxi return sum def _get_qstar_uncertainty(self, data): """ Compute invariant uncertainty with with pinhole data. This uncertainty is given as follow: dq_star = math.sqrt[(x0**2*(dy0)*dx0)**2 + (x1**2 *(dy1)*dx1)**2 + ..+ (xn**2 *(dyn)*dxn)**2 ] where n >= len(data.x)-1 dxi = 1/2*(xi+1 - xi) + (xi - xi-1) dx0 = (x1 - x0)/2 dxn = (xn - xn-1)/2 dyn: error on dy @param data: note: if data doesn't contain dy assume dy= math.sqrt(data.y) """ if len(data.x) <= 1 or len(data.y) <= 1 or \ len(data.x) != len(data.y) or \ (data.dy is not None and (len(data.dy) != len(data.y))): msg = "Length of data.x and data.y must be equal" msg += " and greater than 1; got x=%s, y=%s"%(len(data.x), len(data.y)) raise ValueError, msg else: #Create error for data without dy error if data.dy is None: dy = math.sqrt(y) else: dy = data.dy n = len(data.x) - 1 #compute the first delta dx0 = (data.x[1] - data.x[0])/2 #compute the last delta dxn= (data.x[n] - data.x[n-1])/2 sum = 0 sum += (data.x[0] * data.x[0] * dy[0] * dx0)**2 sum += (data.x[n] * data.x[n] * dy[n] * dxn)**2 if len(data.x) == 2: return math.sqrt(sum) else: #iterate between for element different from the first and the last for i in xrange(1, n-1): dxi = (data.x[i+1] - data.x[i-1])/2 sum += (data.x[i] * data.x[i] * dy[i] * dxi)**2 return math.sqrt(sum) def _get_extrapolated_data(self, model, npts=INTEGRATION_NSTEPS, q_start=Q_MINIMUM, q_end=Q_MAXIMUM): """ @return extrapolate data create from data """ #create new Data1D to compute the invariant q = numpy.linspace(start=q_start, stop=q_end, num=npts, endpoint=True) iq = model.evaluate_model(q) diq = model.evaluate_model_errors(q) result_data = LoaderData1D(x=q, y=iq, dy=diq) return result_data def get_qstar_low(self): """ Compute the invariant for extrapolated data at low q range. Implementation: data = self._get_extra_data_low() return self._get_qstar() @return q_star: the invariant for data extrapolated at low q. """ # Data boundaries for fitting qmin = self._data.x[0] qmax = self._data.x[self._low_extrapolation_npts - 1] # Extrapolate the low-Q data self._fit(model=self._low_extrapolation_function, qmin=qmin, qmax=qmax, power=self._low_extrapolation_power) # Distribution starting point q_start = Q_MINIMUM if Q_MINIMUM >= qmin: q_start = qmin/10 data = self._get_extrapolated_data(model=self._low_extrapolation_function, npts=INTEGRATION_NSTEPS, q_start=q_start, q_end=qmin) # Systematic error # If we have smearing, the shape of the I(q) distribution at low Q will # may not be a Guinier or simple power law. The following is a conservative # estimation for the systematic error. err = qmin*qmin*math.fabs((qmin-q_start)*(data.y[0] - data.y[INTEGRATION_NSTEPS-1])) return self._get_qstar(data), self._get_qstar_uncertainty(data)+err def get_qstar_high(self): """ Compute the invariant for extrapolated data at high q range. Implementation: data = self._get_extra_data_high() return self._get_qstar() @return q_star: the invariant for data extrapolated at high q. """ # Data boundaries for fitting x_len = len(self._data.x) - 1 qmin = self._data.x[x_len - (self._high_extrapolation_npts - 1)] qmax = self._data.x[x_len] q_end = Q_MAXIMUM # fit the data with a model to get the appropriate parameters self._fit(model=self._high_extrapolation_function, qmin=qmin, qmax=qmax, power=self._high_extrapolation_power) #create new Data1D to compute the invariant data = self._get_extrapolated_data(model=self._high_extrapolation_function, npts=INTEGRATION_NSTEPS, q_start=qmax, q_end=q_end) return self._get_qstar(data), self._get_qstar_uncertainty(data) def get_extra_data_low(self, npts_in=None, q_start=Q_MINIMUM, nsteps=INTEGRATION_NSTEPS): """ This method generates 2 data sets , the first is a data created during low extrapolation . its y is generated from x in [ Q_MINIMUM - the minimum of data.x] and the outputs of the extrapolator . (data is the data used to compute invariant) the second is also data produced during the fit but the x range considered is within the reel range of data x. x uses is in [minimum of data.x up to npts_in points] @param npts_in: the number of first points of data to consider for computing y's coming out of the fit. @param q_start: is the minimum value to uses for extrapolated data @param npts: the number of point used to create extrapolated data """ # Create a data from result of the fit for a range outside of the data # at low q range q_start = max(Q_MINIMUM, q_start) qmin = min(self._data.x) if q_start < qmin: data_out_range = self._get_extrapolated_data(model=self._low_extrapolation_function, npts=nsteps, q_start=q_start, q_end=qmin) else: data_out_range = LoaderData1D(x=numpy.zeros(0), y=numpy.zeros(0)) # Create data from the result of the fit for a range inside data q range for # low q if npts_in is None : npts_in = self._low_extrapolation_npts x = self._data.x[:npts_in] y = self._low_extrapolation_function.evaluate_model(x=x) data_in_range = LoaderData1D(x=x, y=y) return data_out_range, data_in_range def get_extra_data_high(self, npts_in=None, q_end=Q_MAXIMUM, nsteps=INTEGRATION_NSTEPS ): """ This method generates 2 data sets , the first is a data created during low extrapolation . its y is generated from x in [ the maximum of data.x to Q_MAXIMUM] and the outputs of the extrapolator . (data is the data used to compute invariant) the second is also data produced during the fit but the x range considered is within the reel range of data x. x uses is from maximum of data.x up to npts_in points before data.x maximum. @param npts_in: the number of first points of data to consider for computing y's coming out of the fit. @param q_end: is the maximum value to uses for extrapolated data @param npts: the number of point used to create extrapolated data """ #Create a data from result of the fit for a range outside of the data # at low q range qmax = max(self._data.x) if q_end != Q_MAXIMUM or nsteps != INTEGRATION_NSTEPS: if q_end > Q_MAXIMUM: q_end = Q_MAXIMUM elif q_end <= qmax: q_end = qmax * 10 #compute the new data with the proper result of the fit for different #boundary and step, outside of data data_out_range = self._get_extrapolated_data(model=self._high_extrapolation_function, npts=nsteps, q_start=qmax, q_end=q_end) else: data_out_range = LoaderData1D(x=numpy.zeros(0), y=numpy.zeros(0)) #Create data from the result of the fit for a range inside data q range for #high q if npts_in is None : npts_in = self._high_extrapolation_npts x_len = len(self._data.x) x = self._data.x[(x_len-npts_in):] y = self._high_extrapolation_function.evaluate_model(x=x) data_in_range = LoaderData1D(x=x, y=y) return data_out_range, data_in_range def set_extrapolation(self, range, npts=4, function=None, power=None): """ Set the extrapolation parameters for the high or low Q-range. Note that this does not turn extrapolation on or off. @param range: a keyword set the type of extrapolation . type string @param npts: the numbers of q points of data to consider for extrapolation @param function: a keyword to select the function to use for extrapolation. of type string. @param power: an power to apply power_low function """ range = range.lower() if range not in ['high', 'low']: raise ValueError, "Extrapolation range should be 'high' or 'low'" function = function.lower() if function not in ['power_law', 'guinier']: raise ValueError, "Extrapolation function should be 'guinier' or 'power_law'" if range == 'high': if function != 'power_law': raise ValueError, "Extrapolation only allows a power law at high Q" self._high_extrapolation_npts = npts self._high_extrapolation_power = power else: if function == 'power_law': self._low_extrapolation_function = PowerLaw() else: self._low_extrapolation_function = Guinier() self._low_extrapolation_npts = npts self._low_extrapolation_power = power def get_qstar(self, extrapolation=None): """ Compute the invariant of the local copy of data. @param extrapolation: string to apply optional extrapolation @return q_star: invariant of the data within data's q range @warning: When using setting data to Data1D , the user is responsible of checking that the scale and the background are properly apply to the data """ self._qstar = self._get_qstar(self._data) self._qstar_err = self._get_qstar_uncertainty(self._data) if extrapolation is None: return self._qstar # Compute invariant plus invariant of extrapolated data extrapolation = extrapolation.lower() if extrapolation == "low": qs_low, dqs_low = self.get_qstar_low() qs_hi, dqs_hi = 0, 0 elif extrapolation == "high": qs_low, dqs_low = 0, 0 qs_hi, dqs_hi = self.get_qstar_high() elif extrapolation == "both": qs_low, dqs_low = self.get_qstar_low() qs_hi, dqs_hi = self.get_qstar_high() self._qstar += qs_low + qs_hi self._qstar_err = math.sqrt(self._qstar_err*self._qstar_err \ + dqs_low*dqs_low + dqs_hi*dqs_hi) return self._qstar def get_surface(self, contrast, porod_const, extrapolation=None): """ Compute the surface of the data. Implementation: V= self.get_volume_fraction(contrast, extrapolation) Compute the surface given by: surface = (2*pi *V(1- V)*porod_const)/ q_star @param contrast: contrast value to compute the volume @param porod_const: Porod constant to compute the surface @param extrapolation: string to apply optional extrapolation @return: specific surface """ # Compute the volume volume = self.get_volume_fraction(contrast, extrapolation) return 2 * math.pi * volume *(1 - volume) * float(porod_const)/self._qstar def get_volume_fraction(self, contrast, extrapolation=None): """ Compute volume fraction is deduced as follow: q_star = 2*(pi*contrast)**2* volume( 1- volume) for k = 10^(-8)*q_star/(2*(pi*|contrast|)**2) we get 2 values of volume: with 1 - 4 * k >= 0 volume1 = (1- sqrt(1- 4*k))/2 volume2 = (1+ sqrt(1- 4*k))/2 q_star: the invariant value included extrapolation is applied unit 1/A^(3)*1/cm q_star = self.get_qstar() the result returned will be 0 <= volume <= 1 @param contrast: contrast value provides by the user of type float. contrast unit is 1/A^(2)= 10^(16)cm^(2) @param extrapolation: string to apply optional extrapolation @return: volume fraction @note: volume fraction must have no unit """ if contrast <= 0: raise ValueError, "The contrast parameter must be greater than zero" # Make sure Q star is up to date self.get_qstar(extrapolation) if self._qstar <= 0: raise RuntimeError, "Invalid invariant: Invariant Q* must be greater than zero" # Compute intermediate constant k = 1.e-8 * self._qstar/(2 * (math.pi * math.fabs(float(contrast)))**2) # Check discriminant value discrim = 1 - 4 * k # Compute volume fraction if discrim < 0: raise RuntimeError, "Could not compute the volume fraction: negative discriminant" elif discrim == 0: return 1/2 else: volume1 = 0.5 * (1 - math.sqrt(discrim)) volume2 = 0.5 * (1 + math.sqrt(discrim)) if 0 <= volume1 and volume1 <= 1: return volume1 elif 0 <= volume2 and volume2 <= 1: return volume2 raise RuntimeError, "Could not compute the volume fraction: inconsistent results" def get_qstar_with_error(self, extrapolation=None): """ Compute the invariant uncertainty. This uncertainty computation depends on whether or not the data is smeared. @param extrapolation: string to apply optional extrapolation @return: invariant, the invariant uncertainty """ self.get_qstar(extrapolation) return self._qstar, self._qstar_err def get_volume_fraction_with_error(self, contrast, extrapolation=None): """ Compute uncertainty on volume value as well as the volume fraction This uncertainty is given by the following equation: dV = 0.5 * (4*k* dq_star) /(2* math.sqrt(1-k* q_star)) for k = 10^(-8)*q_star/(2*(pi*|contrast|)**2) q_star: the invariant value including extrapolated value if existing dq_star: the invariant uncertainty dV: the volume uncertainty The uncertainty will be set to -1 if it can't be computed. @param contrast: contrast value @param extrapolation: string to apply optional extrapolation @return: V, dV = volume fraction, error on volume fraction """ volume = self.get_volume_fraction(contrast, extrapolation) # Compute error k = 1.e-8 * self._qstar /(2 * (math.pi* math.fabs(float(contrast)))**2) # Check value inside the sqrt function value = 1 - k * self._qstar if (value) <= 0: uncertainty = -1 # Compute uncertainty uncertainty = math.fabs((0.5 * 4 * k * self._qstar_err)/(2 * math.sqrt(1 - k * self._qstar))) return volume, uncertainty def get_surface_with_error(self, contrast, porod_const, extrapolation=None): """ Compute uncertainty of the surface value as well as the surface value. The uncertainty is given as follow: dS = porod_const *2*pi[( dV -2*V*dV)/q_star + dq_star(v-v**2) q_star: the invariant value dq_star: the invariant uncertainty V: the volume fraction value dV: the volume uncertainty @param contrast: contrast value @param porod_const: porod constant value @param extrapolation: string to apply optional extrapolation @return S, dS: the surface, with its uncertainty """ # We get the volume fraction, with error # get_volume_fraction_with_error calls get_volume_fraction # get_volume_fraction calls get_qstar # which computes Qstar and dQstar v, dv = self.get_volume_fraction_with_error(contrast, extrapolation) s = self.get_surface(contrast=contrast, porod_const=porod_const) ds = porod_const * 2 * math.pi * (( dv - 2 * v * dv)/ self._qstar\ + self._qstar_err * ( v - v**2)) return s, ds