Skip to content

Commit

Permalink
Added some unit tests. Refs #7653.
Browse files Browse the repository at this point in the history
  • Loading branch information
@wdzhou
wdzhou committed Aug 21, 2013
1 parent e29d710 commit 340b00a
Show file tree
Hide file tree
Showing 3 changed files with 85 additions and 45 deletions.
2 changes: 1 addition & 1 deletion Code/Mantid/Framework/API/src/IPowderDiffPeakFunction.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -174,7 +174,7 @@ namespace API
{
// Throw exception if tried to reset the miller index
stringstream errss;
errss << "ThermalNeutronBk2BkExpConvPVoigt Peak cannot have (HKL) reset.";
errss << "Profile function " << name() << "cannot have (HKL) reset.";
g_log.error(errss.str());
throw runtime_error(errss.str());
}
Expand Down
33 changes: 13 additions & 20 deletions Code/Mantid/Framework/CurveFitting/src/NeutronBk2BkExpConvPVoigt.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -138,6 +138,7 @@ namespace CurveFitting
*/
NeutronBk2BkExpConvPVoigt::NeutronBk2BkExpConvPVoigt()
{
mHKLSet = false;
}

//----------------------------------------------------------------------------------------------
Expand Down Expand Up @@ -249,27 +250,13 @@ namespace CurveFitting
beta = beta0 + beta1/pow(dh, 4.);
tof_h = zero + dtt1*dh + dtt2*dh*dh;

// - Start to calculate alpha, beta, sigma2, gamma,
#if 0
double n = 0.5*gsl_sf_erfc(wcross*(Tcross-1/dh));

double alpha_e = alph0 + alph1*dh;
double alpha_t = alph0t - alph1t/dh;
alpha = 1/(n*alpha_e + (1-n)*alpha_t);

double beta_e = beta0 + beta1*dh;
double beta_t = beta0t - beta1t/dh;
beta = 1/(n*beta_e + (1-n)*beta_t);

double Th_e = zero + dtt1*dh;
double Th_t = zerot + dtt1t*dh - dtt2t/dh;
tof_h = n*Th_e + (1-n)*Th_t;
#endif

sigma2 = sig0*sig0 + sig1*sig1*std::pow(dh, 2) + sig2*sig2*std::pow(dh, 4);
gamma = gam0 + gam1*dh + gam2*std::pow(dh, 2);

// - Calcualte H for the peak
g_log.notice() << "[DBx123] Gam-0 = " << gam0 << ", Gam-1 = " << gam1 << ", Gam-2 = " << gam2
<< ", Gamma = " << gamma << ".\n";

// Calcualte H for the peak
calHandEta(sigma2, gamma, H, eta);

N = alpha*beta*0.5/(alpha+beta);
Expand Down Expand Up @@ -494,6 +481,10 @@ namespace CurveFitting
{
g_log.warning() << "Calculated eta = " << eta << " is out of range [0, 1].\n";
}
else
{
g_log.notice() << "[DBx121] Eta = " << eta << "; Gamma = " << gamma << ".\n";
}

return;
}
Expand All @@ -508,14 +499,14 @@ namespace CurveFitting
const double sigma2, const double invert_sqrt2sigma,
const bool explicitoutput) const
{
// Transform to variable u, v, y, z
const double u = 0.5*alpha*(alpha*sigma2+2.*x);
const double y = (alpha*sigma2 + x)*invert_sqrt2sigma;

const double v = 0.5*beta*(beta*sigma2 - 2.*x);
const double z = (beta*sigma2 - x)*invert_sqrt2sigma;

// 2. Calculate

// Calculate Gaussian part
const double erfcy = gsl_sf_erfc(y);
double part1(0.);
if (fabs(erfcy) > DBL_MIN)
Expand All @@ -527,6 +518,8 @@ namespace CurveFitting
part2 = exp(v)*erfcz;

const double omega1 = (1.-eta)*N*(part1 + part2);

// Calculate Lorenzian part
double omega2(0.);
if (eta >= 1.0E-8)
{
Expand Down
95 changes: 71 additions & 24 deletions Code/Mantid/Framework/CurveFitting/test/NeutronBk2BkExpConvPVoigtTest.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,11 @@ class NeutronBk2BkExpConvPVoigtTest : public CxxTest::TestSuite
//----------------------------------------------------------------------------------------------
/** Set and get parameter
*/
void test_accessParameter()
void Ptest_accessParameter()
{
NeutronBk2BkExpConvPVoigt func;
func.initialize();

func.setParameter("Dtt1", 1.0);
double dtt1 = func.getParameter("Dtt1");
TS_ASSERT_EQUALS(dtt1, 1.0);
Expand All @@ -38,53 +40,73 @@ class NeutronBk2BkExpConvPVoigtTest : public CxxTest::TestSuite
//----------------------------------------------------------------------------------------------
/** Calculate peak positions: data is from Fullprof's sample: arg_si
*/
void test_calculatePeakPositions()
void Ptest_calculatePeakPositions()
{
// (1,1,1)
NeutronBk2BkExpConvPVoigt func;
func.initialize();

func.setParameter("Dtt1", 7476.910);
func.setParameter("Dtt2", -1.540);
func.setParameter("Zero", -9.227);
func.setParameter("LatticeConstant", 5.431363); // Silicon

// (1,1,1)
func.setMillerIndex(1, 1, 1);
func.calculateParameters(false);
double tofh1 = func.centre();
TS_ASSERT_DELTA(tofh1, 23421.7207, 0.01);

// (2, 2, 0)
func.setMillerIndex(2, 2, 0);
func.calculateParameters(false);
double tofh2 = func.centre();
NeutronBk2BkExpConvPVoigt func220;
func220.initialize();

func220.setParameter("Dtt1", 7476.910);
func220.setParameter("Dtt2", -1.540);
func220.setParameter("Zero", -9.227);
func220.setParameter("LatticeConstant", 5.431363); // Silicon
func220.setMillerIndex(2, 2, 0);
func220.calculateParameters(false);
double tofh2 = func220.centre();
TS_ASSERT_DELTA(tofh2, 14342.8350, 0.01);

// (3,1,1)
func.setMillerIndex(3, 1, 1);
func.calculateParameters(false);
double tofh3 = func.centre();
NeutronBk2BkExpConvPVoigt func311;
func311.initialize();

func311.setParameter("Dtt1", 7476.910);
func311.setParameter("Dtt2", -1.540);
func311.setParameter("Zero", -9.227);
func311.setParameter("LatticeConstant", 5.431363); // Silicon

func311.setMillerIndex(3, 1, 1);
func311.calculateParameters(false);
double tofh3 = func311.centre();
TS_ASSERT_DELTA(tofh3, 12230.9648, 0.01);

// (2, 2, 2)
func.setMillerIndex(2, 2, 2);
func.calculateParameters(false);
double tofh4 = func.centre();
NeutronBk2BkExpConvPVoigt func222;
func222.initialize();

func222.setParameter("Dtt1", 7476.910);
func222.setParameter("Dtt2", -1.540);
func222.setParameter("Zero", -9.227);
func222.setParameter("LatticeConstant", 5.431363); // Silicon
func222.setMillerIndex(2, 2, 2);

func222.calculateParameters(false);
double tofh4 = func222.centre();
TS_ASSERT_DELTA(tofh4, 11710.0332, 0.01);


}
/* From arg_si.prf
11710.0332 0 ( 2 2 2) 0 1
12230.9648 0 ( 3 1 1) 0 1
14342.8350 0 ( 2 2 0) 0 1
23421.7207 0 ( 1 1 1) 0 1
*/

//----------------------------------------------------------------------------------------------
/** Calculate peak positions: data is from Fullprof's sample: arg_si
*/
void test_calculatePeakShape()
{
NeutronBk2BkExpConvPVoigt func;
func.initialize();

func.setParameter("Dtt1", 7476.910);
func.setParameter("Dtt2", -1.540);
Expand All @@ -95,18 +117,43 @@ class NeutronBk2BkExpConvPVoigtTest : public CxxTest::TestSuite
func.setParameter("Alph1", 0.597100);
func.setParameter("Beta0", 0.042210);
func.setParameter("Beta1", 0.009460);
func.setParameter("Sig0", 3.032);
func.setParameter("Sig1", 33.027);
func.setParameter("Sig0", sqrt(3.032));
func.setParameter("Sig1", sqrt(33.027));
func.setParameter("Sig2", 0.000);
func.setParameter("Gam0", 0.000);
func.setParameter("Gam0", 2.604);
func.setParameter("Gam0", 0.000);
func.setParameter("Gam1", 2.604);
func.setParameter("Gam2", 0.000);

func.setMillerIndex(1, 1, 0);
func.setMillerIndex(1, 1, 1);
func.calculateParameters(false);

// Peak centre
double tofh1 = func.centre();
TS_ASSERT_DELTA(tofh1, 0.0, 0.01);
TS_ASSERT_DELTA(tofh1, 23421.7207, 0.01);

// Peak shape
func.setParameter("Height", 1.0);

double fwhm = func.fwhm();
TS_ASSERT_DELTA(fwhm, 47.049, 0.001);

/*
cout << "FWHM = " << fwhm << ".\n";
vector<double> vecX;
double tof = tofh1 - 10*fwhm;
while (tof < tofh1 + 10*fwhm)
{
vecX.push_back(tof);
tof += fwhm*0.1;
}
vector<double> vecY(vecX.size(), 0.0);
func.function(vecY, vecX);
for (size_t i = 0; i < vecX.size(); ++i)
cout << vecX[i] << "\t\t" << vecY[i] << "\n";
TS_ASSERT_EQUALS(1, 432);
*/
}

};
Expand Down

0 comments on commit 340b00a

Please sign in to comment.