park.util.peaks

1 # This program is public domain 2 3 """ 4 Peaks and backgrounds. 5 6 Available peaks and backgrounds:: 7 8 gauss, lorentz, voigt, quadratic 9 constant, linear 10 11 The functions are standardized to unit height at zero, but can be adjusted 12 with center, scale and width parameters. 13 14 Models corresponding to these use uppercase names. E.g.:: 15 16 gauss(x, sigma=3, center=2, scale=4) 17 18 corresponds to:: 19 20 m = Gauss(sigma=3, center=2, scale=4) 21 m(x) 22 23 The models define derivatives with respect to the parameters for the 24 individual evaluation points x. E.g.:: 25 26 m.dsigma(x) 27 """ 28 29 import numpy 30 from math import sqrt, pi 31 32 # Peak definitions

33 -def gauss(x, sigma=1, center=0, scale=1):

34 """ 35 Gauss peak. 36 37 The formula used to calculate this is:: 38 39 G(x) = exp( -0.5 (x-center)**2 / sigma**2 ) 40 41 :Parameters: 42 sigma : real 43 The 1-sigma width of the Gaussian 44 center : real 45 Location of the center of the peak 46 scale : real 47 Value at the highest point of the peak 48 49 """ 50 return scale*numpy.exp(-((numpy.asarray(x)-center)/sigma)**2)

51

52 -def lorentz(x, gamma=1, center=0, scale=1):

53 """ 54 Lorentz peak. 55 56 The formula used to calculate this is:: 57 58 L(x) = scale/pi gamma/((x-center)**2 + gamma**2) 59 60 :Parameters: 61 gamma : real 62 Half-width at half-maximum peak width 63 center : real 64 Location of the center of the peak 65 scale : real 66 Value at the highest point of the peak 67 """ 68 return scale/numpy.pi * gamma/((x-center)**2 + gamma**2)

69

70 -def voigt(x, sigma=1, gamma=1, center=0, scale=1):

71 """ 72 Voigt peak. 73 74 The voigt peak is the convolution of a Lorentz peak with a Gaussian peak. 75 76 The formula used to calculate this is:: 77 78 z(x) = (x + 1j gamma) / (sqrt(2) sigma) 79 w(z) = exp(-z**2) erfc(-1j z) / (sqrt(2 pi) sigma) 80 81 V(x) = scale Re(w(z(x-center))) 82 83 :Parameters: 84 gamma : real 85 The half-width half-maximum of the Lorentzian 86 sigma : real 87 The 1-sigma width of the Gaussian, which is one standard deviation. 88 center : real 89 Location of the center of the peak 90 scale : real 91 Value at the highest point of the peak 92 93 Ref: W.I.F. David, J. Appl. Cryst. (1986). 19, 63-64 94 95 Note: adjusted to use stddev and HWHM rather than FWHM parameters 96 """ 97 # wofz function = w(z) = Fad[d][e][y]eva function = exp(-z**2)erfc(-iz) 98 from scipy.special import wofz 99 z = (numpy.asarray(x)-center+1j*gamma)/(sigma*sqrt(2)) 100 V = wofz(z)/(sqrt(2*pi)*sigma) 101 return scale*V.real

102

103 -def quadratic(x, width=1, center=0, scale=0):

104 """ 105 Quadratic peak. 106 107 The formula used to calculate this is:: 108 109 Q(x) = A ( 1 - ( (x - c)/w )**2 / 2 ) if |x-c| < w sqrt(2) 110 Q(x) = 0 otherwise 111 112 :Parameters: 113 width : real 114 The half-width half-maximum of the quadratic 115 center : real 116 Location of the center of the peak 117 scale : real 118 Value at the highest point of the peak 119 """ 120 v = scale*(1 - ((numpy.asarray(x)-center)/width)**2/2) 121 v[v<0] = 0 122 return v

123

124 -def linear(x, center=0, scale=1):

125 """ 126 Sloping background. 127 128 The background has unit height at zero. 129 130 The formula used to calculate this is:: 131 132 S(x) = scale ( x - center + 1) 133 134 :Parameters: 135 center : real 136 Location of the zero crossing 137 scale : real 138 Slope of the line 139 """ 140 return scale*(numpy.asarray(x) - center + 1)

141

142 -def constant(x, scale=0):

143 """ 144 Constant background. 145 146 The formula used to calculate this is:: 147 148 C(x) = scale 149 150 :Parameters: 151 scale : real 152 Value at the highest point of the peak 153 """ 154 return scale*numpy.ones(numpy.asarray(x).shape)

155 156 157 # Models with derivatives

158 -class Gauss(object):

159 """Gauss peak: scale exp ( -0.5 (x-center)**2 / sigma**2 )"""

160 - def __init__(self, scale=1., center=0, sigma=1):

161 self.scale,self.center = scale,center 162 self.sigma = sigma

163 - def __call__(self, x):

164 return gauss(x, scale=self.scale, center=self.center, sigma=self.sigma)

165 - def dscale(self, x):

166 return self(x)/self.scale

167 - def dcenter(self, x):

168 return self(x) * (x - self.center) / self.sigma**2

169 - def dsigma(self, x):

170 return self(x) * (x - self.center)**2 / self.sigma**3

171

172 -class Lorentz(object):

173 """Lorentz peak (HWHM): scale/pi gamma/((x-center)**2 + gamma**2)"""

174 - def __init__(self, scale=1., center=0, gamma=1):

175 self.scale,self.center = scale,center 176 self.gamma = gamma

177 - def __call__(self, x):

178 return lorentz(x, scale=self.scale, center=self.center, 179 gamma=self.gamma)

180 - def dscale(self, x):

181 return self(x)/self.scale

182 - def dcenter(self, x):

183 return self(x) * 2*(x-self.center)/((x-self.center)**2+self.gamma**2)

184 - def dgamma(self, x):

185 return self(x) * (1/self.gamma - self.gamma/((x-self.center)**2 186 + self.gamma**2))

187

188 -class Voigt:

189 """Voigt peak (HWHM,sigma): A [G(sigma) * L(gamma)](x-center)"""

190 - def __init__(self, scale=1, center=0, sigma=1, gamma=1):

191 self.scale,self.center = scale,center 192 self.sigma = sigma 193 self.gamma = gamma

194 - def __call__(self, x):

195 return voigt(x, scale=self.scale, center=self.center, 196 sigma=self.sigma, gamma=self.gamma)

197 - def dscale(self, x):

198 return self(x)/self.scale

199 - def dcenter(self, x):

200 # Not hard given that d/dz erfc = 2 e(-z**2)/sqrt(pi) 201 raise NotImplementedError

202 - def dsigma(self, x):

203 # Not hard given that d/dz erfc = 2 e(-z**2)/sqrt(pi) 204 raise NotImplementedError

205 - def dgamma(self, x):

206 # Not hard given that d/dz erfc = 2 e(-z**2)/sqrt(pi) 207 raise NotImplementedError

208

209 -class Quadratic:

210 """Quadratic peak (HWHM): 1 - 1/2 ((x-c)/w)**2"""

211 - def __init__(self, scale=1, center=0, width=1):

212 self.scale, self.center = scale, center 213 self.width = width

214 - def __call__(self, x):

215 return quadratic(x, scale=self.scale, center=self.center, 216 width=self.width)

217 - def dscale(self, x):

218 return self(x)/self.scale

219 - def dcenter(self, x):

220 v = self.scale/self.width**2 * (x-self.center) 221 v[abs(x-self.center) >= self.width*sqrt(2)] = 0 222 return v

223 - def dwidth(self, x):

224 v = self.scale/self.width**3 * (x-self.center)**2 225 v[abs(x-self.center) >= self.width*sqrt(2)] = 0 226 return v

227

228 -class Linear:

229 """Linear background: scale (x - center + 1)"""

230 - def __init__(self, scale=1, center=0):

231 self.scale, self.center = scale, center

232 - def __call__(self, x):

233 return linear(x, scale=self.scale, center=self.center)

234 - def dscale(self, x):

235 return self(x)/self.scale

236 - def dcenter(self, x):

237 return constant(x, scale=-self.scale)

238

239 -class Constant:

240 """Constant background: scale"""

241 - def __init__(self, scale=1):

242 self.scale = scale

243 - def __call__(self, x):

244 return constant(x, scale=self.scale)

245 - def dscale(self, x):

246 return constant(x, scale=1)

247

Source Code for Module park.util.peaks