# Remove this line and the comment symbol on the next line to get proper patch # diff -u -x '*.o' -x '*.lo' wojdyr-fityk-b558418/src/bfunc.cpp wojdyr-fityk-b558418-mod/src/bfunc.cpp --- wojdyr-fityk-b558418/src/bfunc.cpp 2011-09-29 08:41:42.000000000 +1000 +++ wojdyr-fityk-b558418-mod/src/bfunc.cpp 2012-04-26 12:08:41.000000000 +1000 @@ -491,6 +492,242 @@ } /////////////////////////////////////////////////////////////////////// +/* The FCJAsymm peakshape is that described in Finger, Cox and Jephcoat (1994) + J. Appl. Cryst. vol 27, pp 892-900. */ + +realt FuncFCJAsymm::dfunc_int(realt twopsi, realt twotheta) const { +/* The integral of the FCJ weight function: +0.5(asin((2 cos (twotheta)^2 + 2 sin(twopsi)-2)/|2 sin twopsi - 2|sin(twotheta)) - + asin((2 cos (twotheta)^2 - 2 sin(twopsi)-2)/|2 sin twopsi + 2|sin(twotheta)) ) +Twotheta and twopsi are in radians. We note that the limit as twopsi -> twotheta +is pi/2. + +Note that callers will need to apply the 1/(2_hl) factor found in the FCJ paper. +*/ +realt stwopsi = sin(twopsi); +realt stwoth = sin(twotheta); +realt ctwoth = cos(twotheta); +if(twopsi == 0) return 0.0; +if(twopsi == twotheta) return M_PI/2.0; +return 0.5 * (asin((2.0*ctwoth*ctwoth + 2*stwopsi -2)/(abs(2*stwopsi-2)*stwoth)) - + asin((2.0*ctwoth*ctwoth - 2*stwopsi -2)/(abs(2*stwopsi+2)*stwoth))); +} + +void FuncFCJAsymm::more_precomputations() +{ +denom=0.0; +realt hfunc_max = 0.0; +radians = M_PI/180.0; +cent_rad = av_[1]*radians; +realt hfunc_neg, hfunc_pos; +// If either of the below give a cosine greater than one, set to 0 +// Handle extrema by setting twopsimin to appropriate limit + twopsimin = 0.0; + if (cent_rad > M_PI/2) twopsimin = M_PI; +realt cospsimin = cos(cent_rad)*sqrt(pow(av_[4]+av_[5],2) + 1.0); + if(fabs(cospsimin)<1.0) twopsimin = acos(cospsimin); +twopsiinfl = 0.0; +realt cospsiinfl = cos(cent_rad)*sqrt(pow(av_[4]-av_[5],2) + 1.0); + if(fabs(cospsiinfl)<1.0) twopsiinfl = acos(cospsiinfl); + if(av_[4] == 0 && av_[5] == 0) denom = 1.0; + else { +/* The denominator for the FCJ expression can be calculated analytically. We define it in terms +of the integral of the weight function, dfunc_int: + denom = 2* min(h_l,s_l) * (pi/2 - dfunc_int (twopsiinfl,twotheta)) + + h_l * (v-u) + s_l * (v-u) - (extra_int (twopsiinfl,twotheta) - extra_int (twopsimin,twotheta)) + +where v = 1/(2h_l) * dfunc_int(twopsiinfl,twotheta) and u = 1/2h_l * dfunc_int(twopsimin,twotheta). +extra_int is the integral of 1/cos(psi). +*/ +realt u = 0.5*dfunc_int(twopsimin,cent_rad)/av_[4]; +realt v = 0.5*dfunc_int(twopsiinfl,cent_rad)/av_[4]; +denom_unscaled = 2.0 * fmin(av_[5],av_[4]) * (M_PI/(4.0*av_[4]) - v) + (av_[4] + av_[5])* (v - u) - + (1.0/(2*av_[4]))*0.5*(log(fabs(sin(twopsiinfl) + 1)) - log(fabs(sin(twopsiinfl)-1)) - + log(fabs(sin(twopsimin) + 1)) + log(fabs(sin(twopsimin)-1))); +denom = denom_unscaled * 2.0/fabs(cent_rad-twopsimin); //Scale to [-1,1] interval of G-L integration +// The following two factors are the analytic derivatives of the integral of D with respect to +// h_l and s_l. +df_dh_factor = (1.0/(2*av_[4]))*(dfunc_int(twopsiinfl,cent_rad) - dfunc_int(twopsimin,cent_rad)) - + (1.0/av_[4])*denom_unscaled; +df_ds_factor = (1.0/(2*av_[4]))*(dfunc_int(twopsiinfl,cent_rad) - dfunc_int(twopsimin,cent_rad)) + + (1.0/av_[4])*(M_PI/2 - dfunc_int(twopsiinfl,cent_rad)); +// Precalculate the x and hfunc coefficients +for(int pt=0;pt<512;pt++) { + delta_n_neg[pt] = (cent_rad + twopsimin)/2.0 - (cent_rad-twopsimin)*x1024[pt]/2; + delta_n_pos[pt] = (cent_rad + twopsimin)/2.0 + (cent_rad-twopsimin)*x1024[pt]/2; + hfunc_neg = sqrt(pow(cos(delta_n_neg[pt]),2)/pow(cos(cent_rad),2) -1); + hfunc_pos = sqrt(pow(cos(delta_n_pos[pt]),2)/pow(cos(cent_rad),2)-1); + //Weights for negative half of G-L interval + if(fabs(cos(delta_n_neg[pt])) > fabs(cos(twopsiinfl))) { + weight_neg[pt] = av_[4] + av_[5] - hfunc_neg; + } else { + weight_neg[pt] = 2 * min(av_[4],av_[5]); + } + //Weights for positive half of G-L interval + weight_neg[pt] = weight_neg[pt]/(2.0*av_[4]*hfunc_neg*abs(cos(delta_n_neg[pt]))); + if(fabs(cos(delta_n_pos[pt])) > fabs(cos(twopsiinfl))) { + weight_pos[pt] = av_[4] + av_[5] - hfunc_pos; + } else { + weight_pos[pt] = 2 * min(av_[4],av_[5]); + } + weight_pos[pt] = weight_pos[pt]/(2.0*av_[4]*hfunc_pos*abs(cos(delta_n_pos[pt]))); + // Apply Gauss-Legendre weights + weight_pos[pt]*=w1024[pt]; + weight_neg[pt]*=w1024[pt]; + } + } +} + +bool FuncFCJAsymm::get_nonzero_range(double level, + realt &left, realt &right) const +{ +if (level == 0) + return false; + else if (fabs(level) >= fabs(av_[0])) + left = right = 0; + else { + // As for Pseudo-Voight, allow 4 half-widths for Gaussian + realt pvw = av_[2]*(sqrt(fabs(av_[0]/(level*M_PI*av_[2])-1))+4.); //halfwidths for Lorenzian + // The first non-zero point occurs when the convoluting PV reaches twopsimin. The + // last value occurs when the convoluting PV moves completely past the centre + if(av_[1] < 90) { + left = twopsimin*180/M_PI - pvw; + right = av_[1] + pvw; + } else { + left = av_[1] - pvw; + right = twopsimin*180/M_PI + pvw; + } + } + return true; +} + +//Pseudo-Voight with scaling factors as used in crystallography. +realt FuncFCJAsymm::fcj_psv(realt x, realt location, realt fwhm, realt mixing) const { + realt xa1a2 = (location - x) / fwhm; + realt ex = exp(- 4.0 * M_LN2 * xa1a2 * xa1a2); + ex *= sqrt(4.0*M_LN2/M_PI)/fwhm; + realt lor = 2. / (M_PI * fwhm * (1 + 4* xa1a2 * xa1a2)); + return (1-mixing) * ex + mixing * lor; +} + +CALCULATE_VALUE_BEGIN(FuncFCJAsymm) + realt numer = 0.0; + realt fwhm_rad = av_[2]*2*M_PI/180.0; // Fityk uses hwhm, we use fwhm +if((av_[4]==0 && av_[5]==0) || cent_rad==M_PI/2) { // Plain PseudoVoight + numer = fcj_psv(x*radians,cent_rad,fwhm_rad, av_[3]); + } else { + //do the sum over 1024 Gauss-Legendre points + for(int pt=0; pt < 512; pt++) { + /* Note that the Pseudo-Voight equation for this calculation is that used + in powder diffraction, where the coefficients are chosen to give matching + width parameters, i.e. only one width parameter is necessary and the + width is expressed in degrees (which means a normalised height in + degrees */ + // negative and positive sides + realt psvval = 0.0; + psvval = fcj_psv(delta_n_neg[pt],x*radians,fwhm_rad,av_[3]); + numer += weight_neg[pt] * psvval; + psvval = fcj_psv(delta_n_pos[pt],x*radians,fwhm_rad,av_[3]); + numer += weight_pos[pt] * psvval; + } + } + //Radians scale below to make up for fwhm_rad in denominator of PV function +CALCULATE_VALUE_END(av_[0]*M_PI/180 * numer/denom) + +CALCULATE_DERIV_BEGIN(FuncFCJAsymm) + realt hfunc = 0.0; + realt fwhm_rad = av_[2]*2*M_PI/180.0; + realt numer = 0.0; +realt hfunc_neg,hfunc_pos; //FCJ hfunc for neg, pos regions of G-L integration + realt sumWdGdh = 0.0; // derivative with respect to H/L + realt sumWdGds = 0.0; // derivative with respect to S/L + realt sumWdRdG = 0.0; // sum w_i * dR/dgamma * W(delta,twotheta)/H + realt sumWdRde = 0.0 ;// as above with dR/deta + realt sumWdRdt = 0.0; // as above with dR/d2theta + + //do the sum over 1024 points + // parameters are height,centre,hwhm,eta (mixing),H/L,S/L + // 0 1 2 3 4 5 + for(int pt=0; pt < 512; pt++) { + for(int side=0;side < 2; side++) { + realt xa1a2 = 0.0; + if(side == 0) xa1a2 = (x*radians - delta_n_neg[pt]) / fwhm_rad; + else xa1a2 = (x*radians - delta_n_pos[pt]) / fwhm_rad ; + realt facta = -4.0 * M_LN2 * xa1a2 ; + realt ex = exp(facta * xa1a2); + ex *= sqrt(4.0*M_LN2/M_PI)/fwhm_rad; + realt lor = 2. / (M_PI * fwhm_rad *(1 + 4* xa1a2 * xa1a2)); + realt without_height = (1-av_[3]) * ex + av_[3] * lor; + realt psvval = av_[0] * without_height; + if(side == 0) { + numer += weight_neg[pt] * psvval; + hfunc_neg = 1/(2.0*av_[4]*sqrt(pow(cos(delta_n_neg[pt]),2)/pow(cos(cent_rad),2)-1)); + } else + if(side == 1) { //So hfunc nonzero + numer += weight_pos[pt] * psvval; + hfunc_pos = 1.0/(2.0*av_[4]*sqrt(pow(cos(delta_n_pos[pt]),2)/pow(cos(cent_rad),2)-1)); + } + // pseudo-voight derivatives: first fwhm. The parameter is expressed in + // degrees, but our calculations use radians,so we remember to scale at the end. + realt dRdg = ex/fwhm_rad * (-1.0 + -2.0*facta*xa1a2); //gaussian part:checked + dRdg = av_[0] * ((1-av_[3])* dRdg + av_[3]*(-1.0*lor/fwhm_rad + + 16*xa1a2*xa1a2/(M_PI*(fwhm_rad*fwhm_rad))*1.0/pow(1.0+4*xa1a2*xa1a2,2))); + realt dRde = av_[0] * (lor - ex); //with respect to mixing + realt dRdtt = -1.0* av_[0] * ((1.0-av_[3])*2.0*ex*facta/fwhm_rad - av_[3]*lor*8*xa1a2/(fwhm_rad*(1+4*xa1a2*xa1a2))); + if(side==0) { + /* We know that d(FCJ)/d(param) = sum w[i]*dR/d(param) * W/H(delta) */ + sumWdRdG += weight_neg[pt] * dRdg; + sumWdRde += weight_neg[pt] * dRde; + sumWdRdt += weight_neg[pt] * dRdtt; + } else { + sumWdRdG += weight_pos[pt] * dRdg; + sumWdRde += weight_pos[pt] * dRde; + sumWdRdt += weight_pos[pt] * dRdtt; + } + /* The derivative for h_l includes the convolution of dfunc with PV only up + to twopsiinfl when s_l < h_l as it is zero above this limit. To save + program bytes, we keep the same G-L interval and therefore quadrature + points. This is defensible as it is equivalent to including the zero + valued points in the integration. + The derivative for peak "centre" given here ignores the contribution + from changes in the weighting function, as it is not possible to + numerically integrate the derivative of dfunc that appears in the + full expression. It seems to work. + W is 2 min(s_l,h_l) above twopsiinfl. If W = 2h_l, the 2h_l factor + cancels at top and bottom, leaving a derivative of zero with respect + to both s_l and h_l. If W=2s_l, the derivative with respect to + s_l only is non-zero.*/ + if(fabs(cos(delta_n_pos[pt])) > fabs(cos(twopsiinfl)) && side == 1) { + realt dconvol = w1024[pt] * psvval * hfunc_pos / abs(cos(delta_n_pos[pt])); + sumWdGdh += dconvol; + sumWdGds += dconvol; + } + else if(side == 1 && av_[5] < av_[4]) { + sumWdGds += 2.0* w1024[pt] * psvval *hfunc_pos / abs(cos(delta_n_pos[pt])); + } + if(fabs(cos(delta_n_neg[pt])) > fabs(cos(twopsiinfl)) && side == 0) { + realt dconvol = w1024[pt] * psvval * hfunc_neg / abs(cos(delta_n_neg[pt])); + sumWdGdh += dconvol; + sumWdGds += dconvol; + } + else if (side == 0 && av_[5] < av_[4]) { + sumWdGds += 2.0* w1024[pt] * psvval * hfunc_neg / abs(cos(delta_n_neg[pt])); + } + } + } + // Note that we must scale any numerically integrated terms back to the correct interval + // This scale factor cancels for numer/denom, but not for e.g. numer/denom^2 +dy_dv[0] = M_PI/180 * numer/(av_[0]*denom); // height derivative(note numer contains height) +dy_dv[1] = M_PI*M_PI/(180*180) * sumWdRdt/denom; // peak position in degrees +dy_dv[2] = 2*M_PI*M_PI/(180*180) * sumWdRdG/denom; // fwhm/2 is hwhm, in degrees +dy_dv[3] = M_PI/180 * sumWdRde/denom; // mixing +dy_dv[4] = M_PI/180 * (sumWdGdh/denom - 1.0/av_[4] * numer/denom - + df_dh_factor*numer/(denom_unscaled*denom)); // h_l +dy_dv[5] = M_PI/180* (sumWdGds/denom - df_ds_factor * numer/(denom*denom_unscaled)); // s_l +dy_dx = -1.0*dy_dv[1]; +CALCULATE_DERIV_END(M_PI/180 * numer/denom) + +/////////////////////////////////////////////////////////////////////// void FuncVoigt::more_precomputations() { Only in wojdyr-fityk-b558418-mod/src: bfunc.cpp.orig diff -u -x '*.o' -x '*.lo' wojdyr-fityk-b558418/src/bfunc.h wojdyr-fityk-b558418-mod/src/bfunc.h --- wojdyr-fityk-b558418/src/bfunc.h 2011-09-29 08:41:42.000000000 +1000 +++ wojdyr-fityk-b558418-mod/src/bfunc.h 2012-03-28 16:38:11.000000000 +1100 @@ -134,6 +134,34 @@ bool get_area(realt* a) const; }; +class FuncFCJAsymm : public Function +{ + DECLARE_FUNC_OBLIGATORY_METHODS(FCJAsymm, Function) + void more_precomputations(); + bool get_nonzero_range(double level, realt &left, realt &right) const; + bool get_center(realt* a) const { *a = av_[1]; return true; } + bool get_height(realt* a) const { *a = av_[0]; return true; } + bool get_fwhm(realt* a) const { *a = 2 * fabs(av_[2]); return true; } + realt dfunc_int(realt angle1, realt angle2) const; + realt fcj_psv(realt x, realt location, realt fwhm, realt mixing) const; + static const double x100[]; + static const double w100[]; + static const double x1024[]; + static const double w1024[]; + realt twopsiinfl; + realt twopsimin; + realt cent_rad; + realt radians; + realt delta_n_neg[1024]; //same number of points as x1024 and w1024 + realt delta_n_pos[1024]; + realt weight_neg[1024]; + realt weight_pos[1024]; + realt denom; //denominator constant for given parameters + realt denom_unscaled; //denominator for x-axis in radians + realt df_ds_factor; //derivative with respect to denominator + realt df_dh_factor; //derivative with respect to denominator +}; + class FuncVoigt : public Function { DECLARE_FUNC_OBLIGATORY_METHODS(Voigt, Function) Common subdirectories: wojdyr-fityk-b558418/src/cli and wojdyr-fityk-b558418-mod/src/cli Only in wojdyr-fityk-b558418-mod/src: .deps Only in wojdyr-fityk-b558418-mod/src: gauss-legendre.cpp Only in wojdyr-fityk-b558418-mod/src: gauss-legendre.h Only in wojdyr-fityk-b558418-mod/src: .hg diff -u -x '*.o' -x '*.lo' wojdyr-fityk-b558418/src/info.cpp wojdyr-fityk-b558418-mod/src/info.cpp --- wojdyr-fityk-b558418/src/info.cpp 2011-09-29 08:41:42.000000000 +1000 +++ wojdyr-fityk-b558418-mod/src/info.cpp 2012-03-29 14:34:25.000000000 +1100 @@ -717,7 +717,7 @@ realt x = ep.calculate(); const Model* model = F->get_model(ds); vector symb = model->get_symbolic_derivatives(x); - vector num = model->get_numeric_derivatives(x, 1e-4); + vector num = model->get_numeric_derivatives(x, 1e-8); assert (symb.size() == num.size()); int n = symb.size() - 1; r += "F(" + S(x) + ")=" + S(model->value(x)); Only in wojdyr-fityk-b558418-mod/src: libfityk.la Only in wojdyr-fityk-b558418-mod/src: .libs Only in wojdyr-fityk-b558418-mod/src: Makefile Only in wojdyr-fityk-b558418-mod/src: Makefile.in Only in wojdyr-fityk-b558418-mod/src: swigluarun.h diff -u -x '*.o' -x '*.lo' wojdyr-fityk-b558418/src/tplate.cpp wojdyr-fityk-b558418-mod/src/tplate.cpp --- wojdyr-fityk-b558418/src/tplate.cpp 2011-09-29 08:41:42.000000000 +1000 +++ wojdyr-fityk-b558418-mod/src/tplate.cpp 2012-03-29 16:49:47.000000000 +1100 @@ -83,6 +83,7 @@ FACTORY_FUNC(FuncLogNormal) FACTORY_FUNC(FuncSpline) FACTORY_FUNC(FuncPolyline) +FACTORY_FUNC(FuncFCJAsymm) FACTORY_FUNC(CustomFunction) FACTORY_FUNC(CompoundFunction) @@ -178,6 +179,10 @@ "+shape/(1+((x-center)/hwhm)^2))", /*linear_d=*/false, /*peak_d=*/true, &create_FuncPseudoVoigt); + add("FCJAsymm","height,center,hwhm,shape,h_l,s_l",",,,0.5,,", + "Finger-Cox-Jephcoat asymmetry with PseudoVoight peakshape", + false,true,&create_FuncFCJAsymm); + add("Voigt", "height,center,gwidth,shape", ",,hwhm*0.8,0.1", "convolution of Gaussian and Lorentzian #", /*linear_d=*/false, /*peak_d=*/true, &create_FuncVoigt); Common subdirectories: wojdyr-fityk-b558418/src/wxgui and wojdyr-fityk-b558418-mod/src/wxgui