/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/ /* [ Created with wxMaxima version 12.01.0 ] */ /* [wxMaxima: title start ] Finger-Cox-Jephcoat asymmetry correction [wxMaxima: title end ] */ /* [wxMaxima: section start ] Preliminary definitions [wxMaxima: section end ] */ /* [wxMaxima: input start ] */ hfunc(twopsi,twotheta) := sqrt(cos(twopsi)^2/cos(twotheta)^2 - 1); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ diff(hfunc(twopsi,twotheta),twotheta); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ hfunc_deriv(twopsi,twotheta) := ''%; /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ dfunc(twopsi,twotheta,h_l) := 1/(2*h_l*hfunc(twopsi,twotheta)*cos(twopsi)); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ Dfunc(twopsi,twotheta,h_l,s_l):=dfunc(twopsi,twotheta,h_l)*wfunc(twopsi,twotheta,h_l,s_l); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ wfunc(twopsi,twotheta,h_l,s_l) := block([], if (twopsi < twopsimin(h_l,s_l,twotheta)) then 0 elseif (twopsi < twopsiinfl(h_l,s_l,twotheta)) then (h_l + s_l - hfunc(twopsi,twotheta)) elseif (twopsi < twotheta) then 2 * min(h_l,s_l) else 0); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ twopsimin(h_l,s_l,twotheta) := block([bb], bb: cos(twotheta)*sqrt((h_l+s_l)^2 + 1),if bb>=1 then 0 else acos(bb)); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ twopsiinfl(h_l,s_l,twotheta) := acos(cos(twotheta)*sqrt((h_l-s_l)^2 + 1)); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ev(twopsiinfl(0.04,0.015,4.0*%pi/180),numer); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ plot2d (Dfunc(psi,0.3,0.125,0.02),[psi,0.2,0.2999],[y,0,10.0]); /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] Pseudo-Voight function is defined with constants used in powder diffraction to make the width parameter consistent between the two peaks (and unit area?) width(Gauss) = width/(2 sqrt(ln 2)) width(Lorenz) = width/2 The following equation matches Larry Finger's deposited output for PROFVAL, when multiplied by %pi/180 to account for the normalisation factor from radians to degrees [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ pseudov(twopsi,location,height,width,mixing) := height * ((1-mixing)*sqrt(4*log(2)/%pi)/(width)*exp(-4*log(2)*((location-twopsi)/width)^2) + 2*mixing/(%pi*width*(1 + 4*((location-twopsi)/width)^2))); /* [wxMaxima: input end ] */ /* [wxMaxima: section start ] Create and check analytic form for FCJ denominator [wxMaxima: section end ] */ /* [wxMaxima: comment start ] We believe that the integral over Dfunc has an analytic representation, which means that we can calculate the denominator exactly. Checking... [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ integrate(dfunc(twopsi,twotheta,h_l),twopsi); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ratsimp(%); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ trigsimp(%); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ dfunc_int(twopsi,twotheta,h_l) := ''%; /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ dfunc_int(0.1999999999,0.2,0.125)*4*0.125; /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ normval(centre,h_l,s_l) := quad_qag(Dfunc(delta,centre,h_l,s_l),delta,twopsimin(h_l,s_l,centre),centre-0.000001,6)[1]; /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ part_normval_upper(centre,h_l,s_l) := quad_qag(Dfunc(delta,centre,h_l,s_l),delta, twopsiinfl(h_l,s_l,centre),centre-0.0000001,6)[1]; /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ part_normval_lower(centre,h_l,s_l) := quad_qag(Dfunc(delta,centre,h_l,s_l),delta, twopsimin(h_l,s_l,centre),twopsiinfl(h_l,s_l,centre),6)[1]; /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] In the following, we use the identity that the limit of dfunc_int is pi/(4*h_l) as angle->2theta [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ normval_upper(twotheta,h_l,s_l) := 2*min(s_l,h_l)*(%pi/(4*h_l) - dfunc_int(twopsiinfl(h_l,s_l,twotheta),twotheta,h_l)); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ev(part_normval_upper(4.0*%pi/180,0.015,0.04),numer); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ev(part_normval_upper(4.0*%pi/180,0.015,0.04)-normval_upper(4.0*%pi/180,0.015,0.04),numer); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ extra_int(x) := 1/2 * (log(1/cos(x) + tan(x))); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ other_extra_int(x):= 1/4 * (log(abs(sin(x)+1)) - log(abs(sin(x)-1))); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ev(extra_int(5.0*%pi/180),numer); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ev(other_extra_int(5.0*%pi/180),numer); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ normval_lower(twotheta,h_l,s_l) := block([u,v], v:dfunc_int(twopsiinfl(h_l,s_l,twotheta),twotheta,h_l), u:dfunc_int(twopsimin(h_l,s_l,twotheta),twotheta,h_l), h_l*(v-u) + s_l*(v-u) - (1/h_l)*( other_extra_int(twopsiinfl(h_l,s_l,twotheta)) - other_extra_int(twopsimin(h_l,s_l,twotheta)))); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ev(normval_lower(4.0*%pi/180,0.015,0.04),numer); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ part_normval_lower(4.0*%pi/180,0.015,0.04); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ normval_analytic(twotheta,h_l,s_l) := normval_lower(twotheta,h_l,s_l) + normval_upper(twotheta,h_l,s_l); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ev(normval_analytic(4.0*%pi/180,0.015,0.04),numer)/normval(4.0*%pi/180,0.015,0.04); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ normval(4.0*%pi/180,0.04,0.015); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ev(twopsiinfl(0.015,0.04,4.0*%pi/180),numer); /* [wxMaxima: input end ] */ /* [wxMaxima: section start ] Define the FCJ function [wxMaxima: section end ] */ /* [wxMaxima: comment start ] Now we can define the Finger-Cox-Jephcoat asymmetry function. The multiplicative factor at the front arises from normalising the PV function on a degree-based x-axis: [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ fcj(psi,centre,height,fwhm,mixing,h_l,s_l) := %pi/180 * quad_qag (Dfunc(delta,centre,h_l,s_l)*pseudov(psi-delta,0.0,height,fwhm,mixing), delta,twopsimin(h_l,s_l,centre),centre-0.0000000001,6)[1]/normval_analytic(centre,h_l,s_l); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ev(fcj(3.5*3.14159/180,4*3.14159/180,1,0.06*3.14159/180,0.5,0.04,0.015),numer); /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] Define the FCJ function in alternative representation: [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ fcj_other(psi,centre,height,fwhm,mixing,h_l,s_l) := %pi/180 * quad_qag(Dfunc(psi-delta+centre,centre,h_l,s_l)* pseudov(delta,centre,height,fwhm,mixing),delta,psi+0.00000001,psi-twopsimin(h_l,s_l,centre)+centre,6)[1]/ normval_analytic(centre,h_l,s_l); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ev(fcj_other(3.5*3.14159/180,4*3.14159/180,1,0.06*3.14159/180,0.5,0.04,0.015),numer); /* [wxMaxima: input end ] */ /* [wxMaxima: section start ] Derive the derivatives [wxMaxima: section end ] */ /* [wxMaxima: comment start ] The FCJ derivatives. We need derivatives with respect to height, position, h_l, s_l and mixing All parameters (except peak centre) are separable, so we can use the Liebniz integral rule and the product rule to work through the derivatives in the convolution. d(fcj)/d(mixing) = convolution(D,dR/dmixing)/integral(D) Note that we introduce a factor of pi/180 to account for the 1/width normalisation factor in the Lorentzian definition; if we express our values in radians we need this when reporting against a degree-based x axis. [wxMaxima: comment end ] */ /* [wxMaxima: subsect start ] Derivative of peak shape parameters to check our approach [wxMaxima: subsect end ] */ /* [wxMaxima: input start ] */ diff(pseudov(twopsi,0.0,height,width,mixing),mixing); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ pvdiff_mixing(twopsi,height,width,mixing) := ''%; /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ df_dmixing(twopsi,centre,height,fwhm,mixing,h_l,s_l) := %pi/180 * quad_qag(Dfunc(delta,centre,h_l,s_l)* pvdiff_mixing(twopsi-delta,height,fwhm,mixing),delta,twopsimin(h_l,s_l,centre),centre-0.0000001,6)[1] /normval_analytic(centre,h_l,s_l); /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] Now check that we have the right value for the derivative: [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ ev(df_dmixing(0.15,0.2,100.0,0.01,0.5,0.125,0.02),numer); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ch1 : fcj(0.15,0.2,100.0,0.01,0.5,0.125,0.02); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ch2 : fcj(0.15,0.2,100.0,0.01,0.51,0.125,0.02); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ev((ch2-ch1)/0.01,numer); /* [wxMaxima: input end ] */ /* [wxMaxima: comment start ] Likewise for fwhm and height (height is just a normalised pseudo-voight, left out here): [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ diff(pseudov(delta,0.0,height,fwhm,mixing),fwhm); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ pvdiff_fwhm(delta,height,fwhm,mixing) := ''%; /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ df_dfwhm(twopsi,centre,height,fwhm,mixing,h_l,s_l) := %pi/180 * quad_qag(Dfunc(delta,centre,h_l,s_l)* pvdiff_fwhm(twopsi-delta,height,fwhm,mixing),delta,twopsimin(h_l,s_l,centre),centre-0.00001,6)[1] /normval(centre,h_l,s_l); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ev(df_dfwhm(0.15,0.2,100.0,0.01,0.5,0.125,0.02),numer); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ch1 : fcj(0.15,0.2,100.0,0.01,0.5,0.125,0.02); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ch2 : fcj(0.15,0.2,100.0,0.010001,0.5,0.125,0.02); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ ev((ch2-ch1)/0.000001,numer); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ diff(pseudov(delta,0.0,height,fwhm,mixing),delta); /* [wxMaxima: input end ] */ /* [wxMaxima: subsect start ] Derivatives with respect to h and s [wxMaxima: subsect end ] */ /* [wxMaxima: comment start ] The derivatives with respect to h_l,s_l and centre require evaluating the derivatives of the denominator [wxMaxima: comment end ] */ /* [wxMaxima: comment start ] The derivative of the denominator can be worked through using the Liebniz integral rule, with h(twopsimin)=0. We are ultimately left evaluating the integral of 1/(h(twopsi)*cos(twopsi)) from twopsimin to twopsiinfl [wxMaxima: comment end ] */ /* [wxMaxima: comment start ] Now we put together the full derivative for detector height: dF/dh_l = [convolve(d * R)/normval - 1/h_l FCJ] - df_dh_factor * FCJ where d is dfunc, not Dfunc - i.e. no weighting function. df_dh_factor, defined below, is the derivative of the denominator given by normval [wxMaxima: comment end ] */ /* [wxMaxima: input start ] */ wfunc_hl_diff(twopsi,twotheta,h_l,s_l) := block([],if (twopsi < twopsimin(h_l,s_l,twotheta)) then 0 elseif (twopsi < twopsiinfl(h_l,s_l,twotheta)) then 1 elseif (twopsi < twotheta and h_l < s_l) then 2 else 0); /* [wxMaxima: input end ] */ /* [wxMaxima: input start ] */ df_dh_factor(centre,h_l,s_l):= block([u],u:dfunc_int(twopsiinfl(h_l,s_l,centre),centre,h_l), (u - dfunc_int(twopsimin(h_l,s_l,centre),centre,h_l)) + if (h_l