Source code for refl1d.staj

# This program is in the public domain
# Author: Paul Kienzle
r"""
Read and write staj files

Staj files are the model files for the mlayer and gj2 programs, which are
used as the calculation engine for the reflpak suite. Mlayer supports
unpolarized beam with multilayer models,  and has files ending in
**.staj**. GJ2 supports polarized beam without multilayer models, and
has files ending in **.sta**.
"""

from math import pi
import numpy
from bumps.wsolve import wsolve

ERF_FWHM = 2.35482004503095 # 2 * sqrt(2*log(2))
TANH_FWHM = 0.47320111770856327 # 1/2 atanh(erf(1/sqrt(2))) / acosh(sqrt(2))
# Derivation
# ==========
# Converting from error function FWHM to 1-sigma requires solving
# exp(-0.5*(z/1)**2) = 0.5 for z and doubling it, giving 2*sqrt(2*log(2))
#
# Converting from tanh FWHM to 1-sigma error function is more complicated.
# First find C where w is defined as 1-sigma equivalent of tanh, using the
# identity Erf.CDF(z=sigma;w=sigma) = tanh.CDF(z=sigma;w=sigma).
# This simplifies to::
#
#    erf.CDF  = (1+erf(z/(w*sqrt(2)))/2 = (1+erf(1/sqrt(2)))/2
#    tanh.CDF = (1+tanh(C/w*z))/2       = (1+tanh(C))/2
#
#    erf.CDF = tanh.CDF => C = atanh(erf(1/sqrt(2)))
#
# Next find C where w is defined as FWHM, using the equivalent probability
# density function::
#
#    PDF(z) = C/2w * sech(C/w*z)**2
#
# Solving PDF(w/2) = PDF(0)/2 yields::
#
#    Pw = PDF(w/2) = C/2w * sech(C/2)**2
#    Po = PDF(0) = C/2w * sech(0)**2/2 = C/2w
#
#    Pw = Po/2 => sech(C/2)**2 = 1/2
#              => C = 2 acosh(sqrt(2))
#
# To find 1-sigma width given tanh FWHM of w, use the scale factor
# s = C_1_sigma/C_fwhm = 1/2 atanh(erf(1/sqrt(2)))/acosh(sqrt(2))

#----------------------------------------------------------------------------
#The format of a Non-Magnetic Staj file is shown below.  It is an ASCII text
#file composed of lines as follows:
#
#Note that nLayers = nTLayers + nMLayers + nBLayers + 3
#This is the number of top, middle, and bottom layers plus extra profile lines.
#
#Note that nRepeat is the number of repeats of the mLayers profile lines.
#This parameter is ignored by KsRefl.
#
#Line 1:         nTLayers  nMLayers  nBLayers  nRepeat  nFitParams  nRoughSteps
#Line 2:         wavelength  wavelengthDiv  angularDiv  [theta_offset]
#Line 3:         intensity  background  Qmin  Qmax  nQ (data points)
#Line 4:         profileType ('E' for error function, 'H' for tanh)
#Line 5:         datafileName
#Line 6:         <optional outputfileName>
##Sections have an ignored layer followed by the layers for the section
##The ignored layer of the top section contains vacuum SLD
##Layers have rho  mrho  depth  rough  mu
#Line 7 to 7+nL: sections
#Line 7+nL+1:    p1 p2 p3 ... (fitted parameter numbers)
#Line 7+nL+2 to end: constraints
## Constraints often end with a line of garbage characters.
#
#On reading a Staj file the following conversions are made:
#- scattering_length_density = 1e6 rho / (16 pi)
#- magnetic_scat_len_density = 1e6 mrho / (16 pi)
#- roughness_between_layers  = rough
#- absorption                = mu / (2 wavelength)
#
#-------------------------------------------------------------------------------
#The following shows the contents of an sample Non-Magnetic Staj file:
#            1            1            1            1            0           21
#        6.00000       0.160000    0.000300000       0.000000
#        1.00000   1.00000E-010     0.00741408       0.109617            267
#E
#SNS_Ni_Si_3col.txt
#
#  0.000000E+000  0.000000E+000  0.000000E+000  1.000000E-010  0.000000E+000
#  0.000000E+000  0.000000E+000  9.000000E+002  0.000000E+000  0.000000E+000
#  0.000000E+000  0.000000E+000  0.000000E+000  1.000000E-010  0.000000E+000
#  4.452130E-004  0.000000E+000  7.439860E+002  5.308710E+001  0.000000E+000
#  4.452130E-004  0.000000E+000  0.000000E+000  5.308710E+001  0.000000E+000
#  1.279860E-004  0.000000E+000  6.250000E+002  1.814250E+001  0.000000E+000
#
#MQC1=MQC2
#]
#
#-------------------------------------------------------------------------------
#The sample Non-Magnetic Staj File shown above (after conversion of data
#elements) consists of the following layers:
#
#- T0 (Air)  rho =  0.0, depth =   0.0, mu = 0.0, mrho = 0.0
#                                          roughness between layers =  0.0
#- T1 (Air)  rho =  0.0, depth = 900.0, mu = 0.0, mrho = 0.0
#                                          roughness between layers = 22.5
#- M1 (Ni)   rho = 8.86, depth = 744.0, mu = 0.0, mrho = 0.0
#                                          roughness between layers = 7.70
#- B1 (Si)   rho = 2.55, depth = 625.0, mu = 0.0, mrho = 0.0
#
#Note that comments or additional blank lines are not permitted except after
#all required lines (i.e. line 1 through line 6+nLayers).
#
#Note that the number of values on the first line determines if the format is
#for a non-magnetic (6 values) or a magnetic (4 values) Staj file.

[docs]class MlayerModel(object): r""" Model definition used by MLayer program. **Attributes:** Q values and reflectivity come from a data file with Q, R, dR or from simulation with linear spacing from Qmin to Qmax in equal steps: *data_file* name of the data file, or None if this is simulation only *Qmin*, *Qmax*, *num_Q* for simulation, Q sample points Resolution is defined by wavelength and by incident angle: *wavelength*, *wavelength_dispersion*, *angular_divergence* resolution is calculated as $\Delta Q/Q = \Delta\lambda/\lambda + \Delta\theta/\theta$ Additional beam parameters correct for intensity, background and possibly sample alignment: *intensity*, *background* incident beam intensity and sample background *theta_offset* alignment angle correction The model is defined in terms of layers, with three sections. The top and bottom section correspond to the fixed layers at the surface and the substrate. The middle section layers can be repeated an arbitrary number of times, as defined by the number of repeats attribute. The attributes defining the sections are: *num_top* *num_middle* *num_bottom* section sizes *num_repeats* number of times middle section repeats Interfaces are split into discrete steps according to a profile, either error function or hyperbolic tangent. For sharp interfaces which do not overlap within a layer, the interface is broken into a fixed number of slabs with slabs having different widths, but equal changes in height. For broad interfaces, the whole layer is split into the same fixed number of slabs, but with each slab having the same width. The following attributes are used: *roughness_steps* number of roughness steps (13 is coarse; 51 is fine) *roughness_profile* roughness profile is either 'E' for error function or 'H' for tanh Layers have thickness, interface roughness and real and imaginary scattering length density (SLD). Roughness is stored in the file using full width at half maximum (FWHM) for the given profile type. For convenience, roughness can also be set or queried using a 1-\ $\sigma$ equivalent roughness on an error function profile. Regardless, layer parameters are represented as vectors with one entry for each top, middle and bottom layer using the following attributes: *thickness*, *roughness* : float | |Ang| layer thickness and FWHM roughness *rho*, *irho*, *incoh* : float | |1e-6/Ang^2| complex coherent $\rho + j \rho_i$ and incoherent SLD Computed attributes are provided for convenience: *sigma_roughness* : float | |Ang| 1-\ $\sigma$ equivalent roughness for erf profile *mu* absorption cross section (2*wavelength*irho + incoh) .. Note:: The staj files store SLD as $16\pi\rho$, $2\lambda\rho_i$ with an additional column of 0 for magnetic SLD. This conversion happens automatically on read/write. The incoherent cross section is assumed to be zero. The layers are ordered from surface to substrate. Additional attributes are as follows: *fitpars* individual fit parameter numbers *constraints* constraints between layers *output_file* name of the output file These can be safely ignored, except perhaps if you want to try to compile the constraints into something that can be used by your system. **Methods:** model = MlayerModel(attribute=value, ...) Construct a new MLayer model with the given attributes set. model = MlayerModel.load(filename) Construct a new MLayer model from a staj file. model.set(attribute=value, ...) Replace a set of attribute values. model.fit_resolution(Q,dQ) Choose the best resolution parameters to match the given Q,dQ resolution. Returns the object so that calls can be chained. model.resolution(Q) Return the resolution at Q for the current resolution parameters. model.split_sections() Assign top, middle, bottom and repeats to distribute the layers across sections. Returns the object so that calls can be chained. model.save(filename) Write the model to the given named file. Raises ValueError if the model is invalid. **Constructing new files:** Staj files can be constructed directly. The MlayerModel constructor can accept all data attributes as key word arguments. Models require at least *data_file*, *wavelength*, *thickness*, *roughness* and *rho*. Resolution parameters can be set using model.fit_resolution(Q,dQ). Section sizes can be set using model.split_sections(). Everything else has reasonable defaults. """ data_file = "" Qmin = 0 Qmax = 0.5 num_Q = 200 wavelength = 1 wavelength_dispersion = 0.01 angular_divergence = 0.001 intensity = 1 background = 0 theta_offset = 0 num_top, num_middle, num_bottom = 0, 0, 0 num_repeats = 1 roughness_steps = 13 roughness_profile = 'E' thickness = roughness = rho = None irho = incoh = 0 fitpars = [] constraints = "" output_file = "" def __init__(self, **kw): self.set(**kw)
[docs] def set(self, **kw): valid = ('data_file','Qmin','Qmax','num_Q','wavelength', 'wavelength_dispersion','angular_divergence','intensity', 'background','theta_offset', 'num_top','num_middle','num_bottom','num_repeats', 'roughness_steps','roughness_profile','sigma_roughness', 'thickness','roughness','rho','irho','incoh', 'fitpars','constraints','output_file') for k,v in kw.items(): if k not in valid: raise TypeError("Unexpected attribute '%s' in Mlayer Model"%k) setattr(self,k,v)
@classmethod
[docs] def load(cls, filename): """ Load a staj file, returning an MlayerModel object """ fin = open(filename,'r') lines = fin.readlines() fin.close() self = cls() self._parse(lines) return self
[docs] def save(self, filename): """ Save the staj file """ self._check() fid = open(filename,'w') try: self._write(fid) finally: fid.close()
[docs] def FWHMresolution(self, Q): r""" Return the resolution at Q for mlayer with the current settings for wavelength, wavelength divergence and angular divergence. Resolution is full-width at half maximum (FWHM), not 1-\ $\sigma$. """ return (abs(Q) * self.wavelength_dispersion + 4 * pi * self.angular_divergence) / self.wavelength
[docs] def fit_FWHMresolution(self, Q, dQ, weight=1): r""" Choose the best dL and dT to match the resolution dQ. Given that mlayer uses the following resolution function: .. math:: \Delta Q_k = (|Q_k| \Delta\lambda + 4 \pi \Delta\theta)/\lambda_k we can use a linear system solver to find the optimal $\Delta \lambda$ and $\Delta \theta$ across our dataset from the over-determined system: .. math:: [|Q_k|/\lambda_k, 4\pi/\lambda_k][\Delta\lambda, \Delta\theta]^T = \Delta Q_k If weights are provided (e.g., $\Delta R_k/R_k$), then weigh each point during the fit. Given that the experiment is often run with fixed slits at the start and end, you may choose to match the resolution across the entire $Q$ range, or instead restrict it to just the region where the slits are opening. You will generally want to get the resolution correct at the critical edge since that's where it will have the largest effect on the fit. Returns the object so that operations can be chained. """ A = numpy.array([abs(Q)/self.wavelength, numpy.ones_like(Q)*(4*pi/self.wavelength)]).T s = wsolve(A,y=dQ,dy=weight) self.wavelength_dispersion = s.x self.angular_divergence = s.x return self
[docs] def split_sections(self): """ Split the given set of layers into sections, putting as many layers as possible into the middle section, then the bottom and finally the top. Returns the object so that operations can be chained. """ self.num_repeats = 1 n = len(self.thickness) if n > 28: raise ValueError("A maximum of 28 layers is allowed") if n < 2: raise ValueError("Must have at least two layers") n -= 1 # Incident medium layer if n >= 11: self.num_middle = 9 elif n > 2: self.num_middle = n-2 elif n > 1: self.num_middle = 1 else: self.num_middle = 0 n -= self.num_middle if n >= 10: self.num_bottom = 9 else: self.num_bottom = n-1 n -= self.num_bottom self.num_top = n # This may be zero if there are less than 3 layers. return self
def __str__(self): line = [] if self.data_file is not "": line.append("Data: %s"%self.data_file) else: line.append("Q: %g to %g in %d steps" %(self.Qmin,self.Qmax,self.num_Q)) line.append("Wavelength L: %g dL/L: %g, dTheta: %g" %(self.wavelength, self.wavelength_dispersion, self.angular_divergence)) line.append("Beam intensity: %g Background: %g, Theta offset: %g" %(self.intensity, self.background, self.theta_offset)) profile = 'error function' if self.roughness_profile=='E' else 'tanh' line.append("Interface: %s in %d steps" %(profile, self.roughness_steps)) w,s = self.thickness, self.roughness rho,mu = self.rho, self.mu line.append("Layers:") line.append((" " + "%15s "*4) %("Width (A)","Interface (FWHM)", "Rho (1e-6/A)","Mu (1e-6/A)")) for i in range(len(self.rho)): if i == 0: name = 'V' elif i < self.num_top + 1: name = 'T%d'%(i) elif i < self.num_top + self.num_middle + 1: name = 'M%d'%(i-self.num_top) else: name = 'B%d'%(i-self.num_top-self.num_middle) line.append(("%s:"+("%15g "*4))%(name,w[i],s[i],rho[i],mu[i])) if self.constraints != "": line.append("Constraints:") line.append(self.constraints) return "\n".join(line) def _check(self): """ Verify that the staj file is correct and ready for writing, filling in the details that are missing. """ if ((self.irho is not None and len(self.rho) != len(self.irho)) or len(self.rho) != len(self.thickness) or len(self.rho) != len(self.roughness)): raise ValueError("layer parameters have different lengths") # Could check if Qmin/Qmax/num_Q matches data, but I don't think # it matters, so skip it if self.num_top+self.num_middle+self.num_bottom+1 != len(self.rho): raise ValueError("section sizes do not match number of layers") def _get_mu(self): mu = 2*self.wavelength*self.irho + self.incoh return mu*numpy.ones_like(self.rho) def _set_mu(self, mu): self.irho = mu/(2*self.wavelength) mu = property(fget=_get_mu, fset=_set_mu) def _get_sigma(self): if self.roughness_profile == 'H': scale = TANH_FWHM else: scale = ERF_FWHM return self.roughness/scale def _set_sigma(self, v): if self.roughness_profile == 'H': scale = TANH_FWHM else: scale = ERF_FWHM self.roughness = v * scale sigma_roughness = property(fget=_get_sigma, fset=_set_sigma) def _parse(self, lines): #1: num_top num_middle num_bottom num_repeats num_fit rough_steps nums = [int(s) for s in lines.split()] self.num_top, self.num_middle, self.num_bottom, self.num_repeats, \ _, self.roughness_steps = nums #2: wavelength wavelength_dispersion angular_divergence [theta_offset] nums = [float(s) for s in lines.split()] self.wavelength, self.wavelength_dispersion, self.angular_divergence \ = nums[:3] self.theta_offset = nums if len(nums) > 3 else 0 #3: intensity background Qmin Qmax num_Q nums = [float(s) for s in lines.split()] self.intensity, self.background, self.Qmin, self.Qmax = nums[:4] self.num_Q = int(nums) #4: profile_type ('E' for error function, 'H' for tanh) self.roughness_profile = lines.strip().upper()[:1] if self.roughness_profile not in ('E','H'): raise ValueError("Expected roughness profile type E or H") #5: data_file #6: output_file self.data_file = lines.strip() self.output_file = lines.strip() #nL = num_top+num_middle+num_bottom+3 #7 to 7+nL: rho mrho depth rough mu #ignore the layer before each section nL = self.num_top+self.num_middle+self.num_bottom+3 layers = [[float(v) for v in line.split()] for line in lines[6:6+nL]] del layers[self.num_top+1] del layers[self.num_top+self.num_middle+1] A = numpy.array(layers) self.rho = A[:,0] * (1e6/16/pi) self.irho = A[:,4] * (1e6/2/self.wavelength) self.incoh = A[:,0] * 0 self.thickness = A[:,2] self.roughness = A[:,3] #7+nL+1: P1 P2 P3 ... (fit parameters) self.fitpars = [int(s) for s in lines[6+nL].split()] #7+nL+2 to end: constraints self.constraints = "".join(lines[6+nL+1:]) def _write(self, fid): #1: num_top num_middle num_bottom num_repeats num_fit rough_steps fid.write("%d %d %d %d %d %d\n"%(self.num_top, self.num_middle, self.num_bottom, self.num_repeats, len(self.fitpars), self.roughness_steps)) #2: wavelength wavelength_dispersion angular_divergence [theta_offset] fid.write("%g %g %g %g\n"%(self.wavelength, self.wavelength_dispersion, self.angular_divergence, self.theta_offset)) #3: intensity background Qmin Qmax num_Q fid.write("%g %g %g %g %d\n"%(self.intensity, self.background, self.Qmin, self.Qmax, self.num_Q)) #4: profile_type ('E' for error function, 'H' for tanh) fid.write("%s\n"%self.roughness_profile) #5: data_file #6: output_file fid.write("%s\n%s\n"%(self.data_file, self.output_file)) #nL = num_top+num_middle+num_bottom+3 #7 to 7+nL: rho mrho depth rough mu #ignore the layer before each section rho = self.rho*(16*pi*1e-6) mu = self.mu*1e-6 w,s = self.thickness, self.roughness def _write_layer(idx): fid.write("%g %g %g %g %g\n"%(rho[idx], 0., w[idx], s[idx], mu[idx])) offset = 0 for n in [self.num_top, self.num_middle, self.num_bottom]: _write_layer(offset) if n == 0: # In the case of only two or three layers, some sections # may have no layers, and need to be filled with a repeated # value from the next section. _write_layer(offset) for i in range(n): _write_layer(i+offset+1) offset += n #7+nL+1: P1 P2 P3 ... (fit parameters) fid.write(" ".join(str(p) for p in self.fitpars)+"\n") #7+nL+2 to end: constraints fid.write(self.constraints)
#============================================================================== # Translate the Staj file into model representation (magnetic case). #============================================================================== #------------------------------------------------------------------------------- #The format of a Magnetic Staj file is shown below. It is an ASCII text file #composed of lines as follows: # #Note that nLayers does not count the top layer, just the middle and bottom #layers (though data for the top (incident) layer is provided) in the file. # #Line 1: wavelength wavelengthDiv angularDiv [aguide] #Line 2: intensity background #Line 3: nLayers nRoughSteps nFitParam #Line 4: Qmin Qmax nQ (data points in 'a' datafle) #Line 5: Qmin Qmax nQ (data points in 'b' datafle) #Line 6: Qmin Qmax nQ (data points in 'c' datafle) #Line 7: Qmin Qmax nQ (data points in 'd' datafle) #Line 8: profileType ('E' for error function, 'H' for tanh) #Line 9: active cross sections (usually 'abcd' or 'ABCD') #Line 10: data file (base name without suffix char such as test.refl) #Line 11: output file ## Layer information follows in sets of 3 lines x (nL+1). #Line next (1): rho depth rough mu #Line next (2): mrho mdepth mrough (mrho is also known as phi) #Line next (3): mtheta ## Fit information fills the remainder of the file #Line 11+3*(nL+1)+1: Fit parameters (integers) #Line 11+3*(nL+1)+2 to end-of-file: Constraint program # #On reading a Staj file the following conversions are made: #- scattering length density = rho * (1e6 / 16 / pi ) #- complex SLD = irho * (1e6 / 2 / wavelength) #- magnetic SLD = mrho * (1e6 / 16 / pi ) #- structural roughness = rough #- magnetic roughness = mrough # #------------------------------------------------------------------------------- #The following shows the contents of an sample Magnetic Staj file: # # 5.00000 0.0500000 1.00000E-005 -90.0000 # 1.00000 1.00000E-010 # 7 7 0 # 0.00613985 0.174845 215 # 0.00613985 0.174845 215 # 0.00613985 0.174845 215 # 0.00613985 0.174845 215 #E # abcd #Test.refl # # 0.000000 0.000000 0.000000 0.000000 # 0.000000 0.000000 0.000000 # 270.000 # 0.000318683 50.0000 2.00000 0.000000 # 0.000000 50.0000 2.00000 # 270.000 # 0.000329239 80.0000 2.00000 0.000000 # 0.000000 80.0000 2.00000 # 270.000 # 0.000352864 216.000 2.00000 0.000000 # 0.000117621 216.000 2.00000 # 270.000 # 0.000329239 80.0000 2.00000 0.000000 # 0.000000 80.0000 2.00000 # 270.000 # 0.000318683 50.0000 2.00000 0.000000 # 0.000000 50.0000 2.00000 # 270.000 # 9.70124E-005 10.0000 2.00000 0.000000 # 0.000000 10.0000 2.00000 # 270.000 # 0.000104050 100.000 2.00000 1.00000E-009 # 0.000000 100.000 2.00000 # 270.000 # #------------------------------------------------------------------------------- #The sample Magnetic Staj File shown above (after conversion of data elements) #consists of the following layers (with mdepth and mrough not shown): # #- MV (Air) rho = 0.0, depth = 0.0, mu = 0.0, mrho = 0.0, mtheta = 270.0 # roughness between layers = 0.0 #- M1 (Pt) rho = 6.34, depth = 50.0, mu = 0.0, mrho = 0.0, mtheta = 270.0 # roughness between layers = 0.849 #- M2 (Cu) rho = 6.55, depth = 80.0, mu = 0.0, mrho = 0.0, mtheta = 270.0 # roughness between layers = 0.849 #- M3 (Co+Ni) rho = 7.02, depth = 216.0, mu = 0.0, mrho = 2.34, mtheta = 270.0 # roughness between layers = 0.849 #- M4 (Cu) rho = 6.55, depth = 80.0, mu = 0.0, mrho = 0.0, mtheta = 270.0 # roughness between layers = 0.849 #- M5 (Pt) rho = 6.34, depth = 50.0, mu = 0.0, mrho = 0.0, mtheta = 270.0 # roughness between layers = 0.849 #- M6 (SiO2) rho = 1.93, depth = 10.0, mu = 0.0, mrho = 0.0, mtheta = 270.0 # roughness between layers = 0.849 #- M7 (Si) rho = 2.07, depth = 100.0, mu = 0.0, mrho = 0.0, mtheta = 270.0 # #Note that comments or additional blank lines are not permitted except after #all required lines (i.e. line 1 through line 11+3*(nLayers+1)). # #Note that the number of values on the first line determines if the format is #for a non-magnetic (6 values) or a magnetic (4 values) Staj file. #
[docs]class MlayerMagnetic(object): r""" Model definition used by GJ2 program. **Attributes:** Q values and reflectivity come from a data file with Q, R, dR or from simulation with linear spacing from Qmin to Qmax in equal steps: *data_file* base name of the data file, or None if this is simulation only *active_xsec* active cross sections (usually 'abcd' for all cross sections) *Qmin*, *Qmax*, *num_Q* for simulation, Q sample points Resolution is defined by wavelength and by incident angle: *wavelength*, *wavelength_dispersion*, *angular_divergence* resolution is calculated as $\Delta Q/Q = \Delta\lambda/\lambda + \Delta\theta/\theta$ Additional beam parameters correct for intensity, background and possibly guide field angle: *intensity*, *background* incident beam intensity and sample background *guide_angle* angle of the guide field Unlike pure structural models, magnetic models are in one large section with no repeats. The single parameter is the number of layers, which is implicit in the length of the layer data and does not need to be an explicit attribute. Interfaces are split into discrete steps according to a profile, either error function or hyperbolic tangent. For sharp interfaces which do not overlap within a layer, the interface is broken into a fixed number of slabs with slabs having different widths, but equal changes in height. For broad interfaces, the whole layer is split into the same fixed number of slabs, but with each slab having the same width. The following attributes are used: *roughness_steps* number of roughness steps (13 is coarse; 51 is fine) *roughness_profile* roughness profile is either 'E' for error function or 'H' for tanh Layers have thickness, interface roughness and real and imaginary scattering length density (SLD). Roughness is stored in the file using full width at half maximum (FWHM) for the given profile type. For convenience, roughness can also be set or queried using a 1-\ $\sigma$ equivalent roughness on an error function profile. Regardless, layer parameters are represented as vectors with one entry for each top, middle and bottom layer using the following attributes: *thickness*, *roughness* : float | |Ang| layer thickness and FWHM roughness *rho*, *irho* : float, float | $16 \pi \rho$, $2\lambda\rho_i$ complex scattering length density *mthickness*, *mroughness* : float | |Ang| magnetic thickness and roughness *mrho* : float | $16 \pi \rho_M$ magnetic scattering length density *mtheta* : float | |deg| magnetic angle *sigma_roughness*, *sigma_mroughness* : float | |Ang| computed 1-\ $\sigma$ equivalent roughness for erf profile The conversion from stored $16 \pi \rho$, $2\lambda\rho_i$ to in memory $10^6 \rho$, $10^6 \rho_i$ happens automatically on read/write. The layers are ordered from surface to substrate. Additional attributes are as follows: *fitpars* individual fit parameter numbers *constraints* constraints between layers *output_file* name of the output file These can be safely ignored, except perhaps if you want to try to compile the constraints into something that can be used by your system. **Methods:** model = MlayerMagnetic(attribute=value, ...) Construct a new MLayer model with the given attributes set. model = MlayerMagnetic.load(filename) Construct a new MLayer model from a sta file. model.set(attribute=value, ...) Replace a set of attribute values. model.fit_resolution(Q,dQ) Choose the best resolution parameters to match the given Q,dQ resolution. Returns the object so that calls can be chained. model.resolution(Q) Return the resolution at Q for the current resolution parameters. model.save(filename) Write the model to the given named file. Raises ValueError if the model is invalid. **Constructing new files:** Staj files can be constructed directly. The MlayerModel constructor can accept all data attributes as key word arguments. Models require at least *data_file*, *wavelength*, *thickness*, *roughness* and *rho*. Resolution parameters can be set using model.fit_resolution(Q,dQ). Everything else has reasonable defaults. """ data_file = "" active_xsec = "abcd" Qmin = 0 Qmax = 0.5 num_Q = 200 wavelength = 1 wavelength_dispersion = 0.01 angular_divergence = 0.001 intensity = 1 background = 0 guide_angle = 270 roughness_steps = 13 roughness_profile = 'E' thickness = roughness = rho = irho = None mthickness = mroughness = mrho = mtheta = None fitpars = [] constraints = "" output_file = "" def __init__(self, **kw): self.set(**kw)
[docs] def set(self, **kw): valid = ('data_file','active_xsec','Qmin','Qmax','num_Q','wavelength', 'wavelength_dispersion','angular_divergence','intensity', 'background','guide_angle', 'roughness_steps','roughness_profile', 'thickness','roughness','rho','irho','sigma_roughness', 'mthickness','mroughness','mrho','mtheta','sigma_mroughness', 'fitpars','constraints','output_file') for k,v in kw.items(): if k not in valid: raise TypeError("Unexpected attribute '%s' in Mlayer Model"%k) setattr(self,k,v)
@classmethod
[docs] def load(cls, filename): """ Load a staj file, returning an MlayerModel object """ fin = open(filename,'r') lines = fin.readlines() fin.close() self = cls() try: self._parse(lines) except: raise ValueError("Improper staj file") return self
[docs] def save(self, filename): """ Save the staj file """ self._check() fid = open(filename,'w') try: self._write(fid) finally: fid.close()
[docs] def FWHMresolution(self, Q): return (abs(Q) * self.wavelength_dispersion + 4 * pi * self.angular_divergence) / self.wavelength
FWHMresolution.__doc__ = MlayerModel.FWHMresolution.__doc__
[docs] def fit_FWHMresolution(self, Q, dQ, weight=1): A = numpy.array([abs(Q)/self.wavelength, numpy.ones_like(Q)*(4*pi/self.wavelength)]) s = wsolve(A,y=dQ,dy=weight) self.wavelength_dispersion = s.x self.angular_divergence = s.x return self
fit_FWHMresolution.__doc__ = MlayerModel.fit_FWHMresolution.__doc__ def __str__(self): line = [] if self.data_file is not "": line.append("Data: %s[%s]" %(self.data_file,self.active_xsec.upper())) else: line.append("Q: %g to %g in %d steps" %(self.Qmin,self.Qmax,self.num_Q)) line.append("Wavelength L: %g dL/L: %g, dTheta: %g" %(self.wavelength, self.wavelength_dispersion, self.angular_divergence)) line.append("Beam intensity: %g Background: %g, Guide angle: %g" %(self.intensity, self.background, self.guide_angle)) profile = 'error function' if self.roughness_profile=='E' else 'tanh' line.append("Interface: %s in %d steps" %(profile, self.roughness_steps)) w,s = self.thickness, self.roughness wm = self.mthickness if self.mthickness is not None else w sm = self.mroughness if self.mroughness is not None else s rho = self.rho irho = self.irho if self.irho is not None else numpy.zeros_like(w) mrho = self.mrho if self.mrho is not None else numpy.zeros_like(w) mtheta = self.mtheta if self.mtheta is not None else 270*numpy.ones_like(w) line.append("Layers:") line.append((" " + ("%15s "*4) + "\n " + ("%15s "*4)) %("Width (A)","Interface (FWHM)", "Rho (1e-6/A)","iRho (1e-6/A)", "Mag width","Mag interface","Mag rho","Mag angle (deg)")) for i in range(len(self.rho)): line.append(("%3d:"+("%15g "*4)+"\n " + ("%15g "*4)) %(i, w[i],s[i],rho[i],irho[i], wm[i],sm[i],mrho[i],mtheta[i])) if self.constraints != "": line.append("Constraints:") line.append(self.constraints) return "\n".join(line) def _check(self): """ Verify that the staj file is correct and ready for writing, filling in the details that are missing. """ ns = [(len(v) if v is not None else 0) for v in (self.rho, self.thickness, self.roughness, self.irho, self.mrho, self.mthickness, self.mroughness, self.mtheta)] if any((n != ns and n != 0) for n in ns[1:]): raise ValueError("layer parameters have different lengths") if any((n==0) for n in ns[:3]): raise ValueError("rho, thickness and roughness are required") # Could check if Qmin/Qmax/num_Q matches data, but I don't think # it matters, so skip it def _get_sigma(self): if self.roughness_profile == 'H': scale = TANH_FWHM else: scale = ERF_FWHM return self.roughness/scale def _set_sigma(self, v): if self.roughness_profile == 'H': scale = TANH_FWHM else: scale = ERF_FWHM self.roughness = v * scale sigma_roughness = property(_get_sigma, _set_sigma) def _get_msigma(self): if self.roughness_profile == 'H': scale = TANH_FWHM else: scale = ERF_FWHM return self.mroughness/scale def _set_msigma(self, v): if self.mroughness_profile == 'H': scale = TANH_FWHM else: scale = ERF_FWHM self.mroughness = v * scale sigma_mroughness = property(_get_msigma, _set_msigma) def _parse(self, lines): #1: wavelength wavelength_dispersion angular_divergence [aguide] nums = [float(s) for s in lines.split()] self.wavelength, self.wavelength_dispersion, self.angular_divergence \ = nums[:3] self.guide_angle = nums if len(nums) > 3 else 0 #2: intensity background nums = [float(s) for s in lines.split()] self.intensity, self.background = nums #3: maxLayer nRoughSteps nFitParam nums = [int(s) for s in lines.split()] maxLayer, self.roughness_steps, _ = nums #4-7: Qmin Qmax nQ (data points in a, b, c and d) # Note that we are only keeping the first one; for simulation all the # others should be the same, and for loading, the datafile will tell # us what Q to use. nums = [float(s) for s in lines.split()] self.Qmin, self.Qmax, self.num_Q = nums, nums, int(nums) #8: profile_type ('E' for error function, 'H' for tanh) self.roughness_profile = lines.strip() if self.roughness_profile not in ('E','H'): raise ValueError("Expected roughness profile type E or H") #9: active cross sections (usually 'abcd' or 'ABCD') self.active_xsec = lines.strip().lower() #10: data file (base name without suffix char such as test.refl) #11: output_file self.data_file = lines.strip() self.output_file = lines.strip() ## Layer information follows in sets of 3 lines x (nL+1). #Line next (1): rho depth rough mu #Line next (2): mrho mdepth mrough (mrho is also known as phi) #Line next (3): mtheta nL = maxLayer + 1 layers = [[float(v) for v in " ".join(lines[11+3*i:11+3*(i+1)]).split()] for i in range(nL)] A = numpy.array(layers) self.rho = A[:,0]* (1e6/16/pi) self.irho = A[:,3] * (1e6/2/self.wavelength) self.thickness = A[:,1] self.roughness = A[:,2] self.mrho = A[:,4] * (1e6/16/pi) self.mtheta = A[:,7] self.mthickness = A[:,5] self.mroughness = A[:,6] ## Fit information fills the remainder of the file #footer+1: P1 P2 P3 ... (fit parameters) self.fitpars = [int(s) for s in lines[10+3*nL+1].split()] #footer+2 to end: constraints self.constraints = "".join(lines[10+3*nL+2:]) def _write(self, fid): #1: wavelength wavelength_dispersion angular_divergence [aguide] fid.write("%g %g %g %g\n"%(self.wavelength, self.wavelength_dispersion, self.angular_divergence, self.guide_angle)) #2: intensity background fid.write("%g %g\n"%(self.intensity, self.background)) #3: maxLayer nRoughSteps nFitParam fid.write("%d %d %d\n"%(len(self.rho)-1, self.roughness_steps, len(self.fitpars))) #4-7: Qmin Qmax nQ (data points in a, b, c and d) fid.write("%g %g %d\n"%(self.Qmin, self.Qmax, self.num_Q)) fid.write("%g %g %d\n"%(self.Qmin, self.Qmax, self.num_Q)) fid.write("%g %g %d\n"%(self.Qmin, self.Qmax, self.num_Q)) fid.write("%g %g %d\n"%(self.Qmin, self.Qmax, self.num_Q)) #8: profile_type ('E' for error function, 'H' for tanh) fid.write("%s\n"%self.roughness_profile) #9: active cross sections (usually 'abcd' or 'ABCD') fid.write(" %s\n"%self.active_xsec) #10: data file (base name without suffix char such as test.refl) #11: output_file fid.write("%s\n%s\n"%(self.data_file, self.output_file)) ## Layer information follows in sets of 3 lines x (nL+1). #Line next (1): rho depth rough mu #Line next (2): mrho mdepth mrough (mrho is also known as phi) #Line next (3): mtheta rho = self.rho*(16*pi/1e6) w,s = self.thickness, self.roughness wm = self.mthickness if self.mthickness is not None else self.thickness sm = self.mroughness if self.mroughness is not None else self.roughness if self.irho is not None: mu = self.irho*(2*self.wavelength/1e6) else: mu = numpy.zeros_like(rho) if self.mrho is not None: mrho = self.mrho * (16*pi/1e6) else: mrho = numpy.zeros_like(rho) if self.mtheta is not None: mtheta = self.mtheta else: mtheta = 270*numpy.ones_like(rho) for i in range(len(rho)): fid.write("%g %g %g %g\n%g %g %g\n%g\n" %(rho[i],w[i],s[i],mu[i],mrho[i],wm[i],sm[i],mtheta[i])) ## Fit information fills the remainder of the file #footer+1: P1 P2 P3 ... (fit parameters) fid.write(" ".join(str(p) for p in self.fitpars)+"\n") #footer+2 to end: constraints fid.write(self.constraints)