voigt_with_derivs

def voigt_with_derivs (   self,   sigma,   alpha,   k_points = None,   xfullwidth = 10.0,   kexptail = 25   ) Compute a block of Voigt function values Parameters: sigma  The Gaussian width for the convolution, such that the Gaussian is y=exp[-x^2 / (2 sigma^2)] alpha  The half-width of the Lorentzian, i.e. y=1/(alpha^2 + x^2), suitably normalized k_points  The number of points in k-space to compute. If it is None, autorange based on kexptail (see below) xfullwidth  The maximum value of X for which the function is calculated. The domain is (-xfullwidth, xfullwidth) kexptail  The number of exponential decrements to compute in k-space. kexptail=25 gives about 10^-11 amplitude at the Nyquist frequency. Larger values give finer sampling of the peak in x-space. Return values: xvals  The grid of x values on which the function is computed. yft  The array of function values corresponding to the xvals given. The function is normalized so that sum(yft*dx)=1 dsig  The derivative of the Voigt function with respect to sigma dalpha  The derivative of the Voigt function with respect to alpha. return an x grid, voigt function with specified params, and d/dsigma and d/dalpha of this. Function is normalized to unit integral of points Definition at line 57 of file voigt_profile.py. 00057 : 00058 """return an x grid, voigt function with specified params, and d/dsigma and d/dalpha of this. 00059 Function is normalized to unit integral of points""" 00060 if k_points is None: 00061 #autorange k by solving for k such that 0.5*sigma**2*k**2+alpha*k==kexptail 00062 # so exp(-kexptail) is smallest point in spectrum 00063 q=-0.5*(alpha+math.sqrt(alpha*alpha+4.0*kexptail*(0.5*sigma**2))) #the Numerical Recipes quadratic trick 00064 k1=-kexptail/q 00065 kp=k1*xfullwidth/math.pi 00066 k_points=2**int(1+math.floor(math.log(kp)/math.log(2.0))) #next power of two up 00067 00068 if self.saved_params[:2] != (k_points, xfullwidth): 00069 kmax=math.pi*k_points/(xfullwidth) 00070 kvals=Numeric.array(range(k_points+1), _numeric_float)*(kmax/k_points) 00071 dx=xfullwidth/(2.0*k_points) 00072 xvals=Numeric.array(range(2*k_points), _numeric_float)*(xfullwidth/float(k_points))-xfullwidth 00073 k2=kvals*kvals 00074 x2=xvals*xvals 00075 cosarray=Numeric.cos((math.pi/xfullwidth)*xvals) 00076 self.saved_params=(k_points, xfullwidth, kvals, xvals, k2, cosarray, x2) 00077 else: 00078 k_points, xfullwidth, kvals, xvals, k2, cosarray, x2=self.saved_params 00079 kmax=math.pi*k_points/(xfullwidth) 00080 dx=xfullwidth/(2.0*k_points) 00081 00082 yvals=Numeric.exp(Numeric.clip((-0.5*sigma*sigma)*k2 - alpha*kvals, -100, 0)) 00083 yvals[1::2]*=-1.0 #modulate phase to put peak in center for convenience 00084 yft=inverse_real_fft(yvals) 00085 yft *= 0.5/dx #scale to unit area 00086 00087 #apply periodic boundary correction to function 00088 delta=2*xfullwidth 00089 k1=2*math.pi*alpha/delta 00090 shk1=math.sinh(k1)/delta 00091 chk1=math.cosh(k1) 00092 uu1=1.0/(chk1-cosarray ) 00093 uu2=Numeric.array(1.0/(x2+alpha**2)) 00094 deltav=shk1*uu1-(alpha/math.pi)*uu2 00095 00096 #this is an empirical correction on top of the analytical one, to get the next order correct 00097 #it is derived from a scaling argument (much hand waving) about 00098 #the effect of a gaussian convolution on a second derivative 00099 yft -= deltav*(1+(2*sigma**2/xfullwidth**4)*x2) 00100 00101 #transform of df/d(alpha) = d(transform of f)/d(alpha) = -k*transform 00102 yvals*=kvals 00103 dalpha=inverse_real_fft(yvals) 00104 dalpha*= -0.5/dx 00105 00106 #apply periodicity correction to d/dalpha 00107 deltad=(-2.0*math.pi*shk1*shk1)*(uu1*uu1) 00108 deltad += ( 2.0*math.pi*chk1/(delta*delta))*uu1 00109 deltad += (alpha*alpha - x2)*(uu2*uu2)/math.pi 00110 dalpha -= deltad 00111 00112 #transform of df/d(sigma) = d(transform of f)/d(sigma) = -k*k*sigma*transform 00113 yvals*=kvals 00114 dsig=inverse_real_fft(yvals) 00115 dsig*= -0.5*sigma/dx 00116 00117 #d/dsigma is very localized... no correction needed 00118 00119 return (xvals, yft, dsig, dalpha) 00120 00121 ##

---