Skip to content

Commit

Permalink
second batch refs #7808
Browse files Browse the repository at this point in the history
  • Loading branch information
@NickDraper
NickDraper committed Aug 22, 2013
1 parent 9587df8 commit b62c13f
Show file tree
Hide file tree
Showing 9 changed files with 65 additions and 84 deletions.
3 changes: 1 addition & 2 deletions Code/Mantid/Framework/Algorithms/src/StripVanadiumPeaks.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,7 +132,6 @@ void StripVanadiumPeaks::exec()
//Middle of each X bin
MantidVec midX;
convertToBinCentre(X, midX);
double totY;

//This'll be the output
MantidVec outY = Y;
Expand All @@ -150,7 +149,7 @@ void StripVanadiumPeaks::exec()
int L1 = getBinIndex(X, center - width * 0.75);
int L2 = getBinIndex(X, center - width * 0.25);
double leftX = (midX[L1] + midX[L2])/2;
totY = 0;
double totY = 0;
for (int i=L1; i<=L2; i++)
totY += Y[i];
double leftY = totY / (L2-L1+1);
Expand Down
8 changes: 3 additions & 5 deletions Code/Mantid/Framework/Crystal/src/CalculateUMatrix.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -182,8 +182,6 @@ namespace Crystal
OrientedLattice o(a,b,c,alpha,beta,gamma);
Matrix<double> B=o.getB();

double H,K,L;

PeaksWorkspace_sptr ws;
ws = AnalysisDataService::Instance().retrieveWS<PeaksWorkspace>(this->getProperty("PeaksWorkspace") );
if (!ws) throw std::runtime_error("Problems reading the peaks workspace");
Expand All @@ -195,9 +193,9 @@ namespace Crystal
for (int i=0;i<ws->getNumberPeaks();i++)
{
Peak p=ws->getPeaks()[i];
H=p.getH();
K=p.getK();
L=p.getL();
double H=p.getH();
double K=p.getK();
double L=p.getL();
if(H*H+K*K+L*L>0)
{
nIndexedpeaks++;
Expand Down
9 changes: 6 additions & 3 deletions Code/Mantid/Framework/Crystal/src/LoadIsawPeaks.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -172,10 +172,11 @@ namespace Crystal
g_log.error()<<"Missing L1 or Time offset"<<std::endl;
throw std::invalid_argument("Missing L1 or Time offset");
}
double L1;

try
{
std::istringstream iss( L1s+" "+T0s, std::istringstream::in);
double L1;
iss>>L1;
iss>>T0;
V3D sampPos=instr->getSample()->getPos();
Expand Down Expand Up @@ -222,11 +223,13 @@ namespace Crystal
startChar = getWord( in, false);// blank lines ?? and # lines ignore

std::istringstream iss( line, std::istringstream::in);
int bankNum,nrows,ncols;
double width,height,depth,detD,Centx,Centy,Centz,Basex,Basey,Basez,
int bankNum;
double width,height,Centx,Centy,Centz,Basex,Basey,Basez,
Upx,Upy,Upz;
try
{
int nrows,ncols;
double depth,detD;
iss>>bankNum>>nrows>>ncols>>width>>height>>depth>>detD
>>Centx>>Centy>>Centz>>Basex>>Basey>>Basez
>>Upx>>Upy>>Upz;
Expand Down
2 changes: 1 addition & 1 deletion Code/Mantid/Framework/Crystal/src/MaskPeaksWorkspace.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,8 +113,8 @@ namespace Mantid
Geometry::IComponent_const_sptr comp = inst->getComponentByName(bankName);
if (!comp)
{
continue;
g_log.debug() << "Component "+bankName+" does not exist in instrument\n";
continue;
}

// determine the range in time-of-flight
Expand Down
3 changes: 1 addition & 2 deletions Code/Mantid/Framework/Crystal/src/OptimizeExtinctionParameters.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,7 +143,6 @@ namespace Mantid
if(type.compare(7,2,"II")==0)nopt = 2;
size_t iter = 0;
int status = 0;
double size;

/* Starting point */
x = gsl_vector_alloc (nopt);
Expand All @@ -170,7 +169,7 @@ namespace Mantid
if (status)
break;

size = gsl_multimin_fminimizer_size (s);
double size = gsl_multimin_fminimizer_size (s);
status = gsl_multimin_test_size (size, 1e-4);

}
Expand Down
2 changes: 1 addition & 1 deletion Code/Mantid/Framework/Crystal/src/SaveIsawUB.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -112,7 +112,7 @@ namespace Crystal
dV+= U*U;
}

Volume = Volume/2.0/(1-xA*xA-xB*xB-xC*xC+2*xA*xB*xC);

double U=(lattice_errors[3])*(sin(2*latticeParams[3]/180.*M_PI)-
sin(latticeParams[3]/180.*M_PI)*cos(latticeParams[4]/180*M_PI)*
cos(latticeParams[5]/180*M_PI));
Expand Down
2 changes: 1 addition & 1 deletion Code/Mantid/Framework/Crystal/src/SortPeaksWorkspace.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -145,7 +145,7 @@ namespace Mantid
} catch (std::invalid_argument& ex)
{
this->g_log.error("Specified ColumnToSortBy does not exist");
throw ex;
throw;
}

}
Expand Down
3 changes: 1 addition & 2 deletions Code/Mantid/Framework/CurveFitting/src/Fit1D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -496,7 +496,6 @@ void Fit1D::exec()

int iter = 0;
int status;
double size; // for simplex algorithm
double finalCostFuncVal;
double dof = static_cast<double>(l_data.n - l_data.p); // dof stands for degrees of freedom

Expand Down Expand Up @@ -536,7 +535,7 @@ void Fit1D::exec()
if (status) // break if error
break;

size = gsl_multimin_fminimizer_size(simplexMinimizer);
double size = gsl_multimin_fminimizer_size(simplexMinimizer);
status = gsl_multimin_test_size(size, 1e-2);
prog.report();
}
Expand Down
117 changes: 50 additions & 67 deletions Code/Mantid/Framework/CurveFitting/src/IkedaCarpenterPV.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -254,9 +254,6 @@ void IkedaCarpenterPV::constFunction(double* out, const double* xValues, const i

const double k = 0.05;

double u,v,s,r;
double yu, yv, ys, yr;

// Not entirely sure what to do if sigmaSquared ever negative
// for now just post a warning
double someConst = std::numeric_limits<double>::max() / 100.0;
Expand All @@ -267,48 +264,42 @@ void IkedaCarpenterPV::constFunction(double* out, const double* xValues, const i
g_log.warning() << "sigmaSquared negative in functionLocal.\n";
}

double R, Nu, Nv, Ns, Nr, N;

std::complex<double> zs, zu, zv, zr;

double alpha, a_minus, a_plus, x, y, z;

// update wavelength vector
calWavelengthAtEachDataPoint(xValues, nData);

for (int i = 0; i < nData; i++) {
double diff=xValues[i]-X0;

R = exp(-81.799/(m_waveLength[i]*m_waveLength[i]*kappa));
alpha = 1.0 / (alpha0+m_waveLength[i]*alpha1);
double R = exp(-81.799/(m_waveLength[i]*m_waveLength[i]*kappa));
double alpha = 1.0 / (alpha0+m_waveLength[i]*alpha1);

a_minus = alpha*(1-k);
a_plus = alpha*(1+k);
x=a_minus-beta;
y=alpha-beta;
z=a_plus-beta;
double a_minus = alpha*(1-k);
double a_plus = alpha*(1+k);
double x=a_minus-beta;
double y=alpha-beta;
double z=a_plus-beta;

Nu=1-R*a_minus/x;
Nv=1-R*a_plus/z;
Ns=-2*(1-R*alpha/y);
Nr=2*R*alpha*alpha*beta*k*k/(x*y*z);
double Nu=1-R*a_minus/x;
double Nv=1-R*a_plus/z;
double Ns=-2*(1-R*alpha/y);
double Nr=2*R*alpha*alpha*beta*k*k/(x*y*z);

u=a_minus*(a_minus*sigmaSquared-2*diff)/2.0;
v=a_plus*(a_plus*sigmaSquared-2*diff)/2.0;
s=alpha*(alpha*sigmaSquared-2*diff)/2.0;
r=beta*(beta*sigmaSquared-2*diff)/2.0;
double u=a_minus*(a_minus*sigmaSquared-2*diff)/2.0;
double v=a_plus*(a_plus*sigmaSquared-2*diff)/2.0;
double s=alpha*(alpha*sigmaSquared-2*diff)/2.0;
double r=beta*(beta*sigmaSquared-2*diff)/2.0;

yu = (a_minus*sigmaSquared-diff)*someConst;
yv = (a_plus*sigmaSquared-diff)*someConst;
ys = (alpha*sigmaSquared-diff)*someConst;
yr = (beta*sigmaSquared-diff)*someConst;
double yu = (a_minus*sigmaSquared-diff)*someConst;
double yv = (a_plus*sigmaSquared-diff)*someConst;
double ys = (alpha*sigmaSquared-diff)*someConst;
double yr = (beta*sigmaSquared-diff)*someConst;

zs = std::complex<double>(-alpha*diff, 0.5*alpha*gamma);
zu = (1-k)*zs;
zv = (1-k)*zs;
zr = std::complex<double>(-beta*diff, 0.5*beta*gamma);
std::complex<double> zs = std::complex<double>(-alpha*diff, 0.5*alpha*gamma);
std::complex<double> zu = (1-k)*zs;
std::complex<double> zv = (1-k)*zs;
std::complex<double> zr = std::complex<double>(-beta*diff, 0.5*beta*gamma);

N = 0.25*alpha*(1-k*k)/(k*k);
double N = 0.25*alpha*(1-k*k)/(k*k);

out[i] = I*N*( (1-eta)*(Nu*exp(u+gsl_sf_log_erfc(yu))+Nv*exp(v+gsl_sf_log_erfc(yv)) +
Ns*exp(s+gsl_sf_log_erfc(ys))+Nr*exp(r+gsl_sf_log_erfc(yr))) -
Expand Down Expand Up @@ -340,10 +331,7 @@ void IkedaCarpenterPV::functionLocal(double* out, const double* xValues, const s
// equations taken from Fullprof manual

const double k = 0.05;

double u,v,s,r;
double yu, yv, ys, yr;


// Not entirely sure what to do if sigmaSquared ever negative
// for now just post a warning
double someConst = std::numeric_limits<double>::max() / 100.0;
Expand All @@ -354,50 +342,45 @@ void IkedaCarpenterPV::functionLocal(double* out, const double* xValues, const s
g_log.warning() << "sigmaSquared negative in functionLocal.\n";
}

double R, Nu, Nv, Ns, Nr, N;

std::complex<double> zs, zu, zv, zr;

double alpha, a_minus, a_plus, x, y, z;


// update wavelength vector
calWavelengthAtEachDataPoint(xValues, nData);

for (size_t i = 0; i < nData; i++) {
double diff=xValues[i]-X0;


R = exp(-81.799/(m_waveLength[i]*m_waveLength[i]*kappa));
alpha = 1.0 / (alpha0+m_waveLength[i]*alpha1);
double R = exp(-81.799/(m_waveLength[i]*m_waveLength[i]*kappa));
double alpha = 1.0 / (alpha0+m_waveLength[i]*alpha1);


a_minus = alpha*(1-k);
a_plus = alpha*(1+k);
x=a_minus-beta;
y=alpha-beta;
z=a_plus-beta;
double a_minus = alpha*(1-k);
double a_plus = alpha*(1+k);
double x=a_minus-beta;
double y=alpha-beta;
double z=a_plus-beta;

Nu=1-R*a_minus/x;
Nv=1-R*a_plus/z;
Ns=-2*(1-R*alpha/y);
Nr=2*R*alpha*alpha*beta*k*k/(x*y*z);
double Nu=1-R*a_minus/x;
double Nv=1-R*a_plus/z;
double Ns=-2*(1-R*alpha/y);
double Nr=2*R*alpha*alpha*beta*k*k/(x*y*z);

u=a_minus*(a_minus*sigmaSquared-2*diff)/2.0;
v=a_plus*(a_plus*sigmaSquared-2*diff)/2.0;
s=alpha*(alpha*sigmaSquared-2*diff)/2.0;
r=beta*(beta*sigmaSquared-2*diff)/2.0;
double u=a_minus*(a_minus*sigmaSquared-2*diff)/2.0;
double v=a_plus*(a_plus*sigmaSquared-2*diff)/2.0;
double s=alpha*(alpha*sigmaSquared-2*diff)/2.0;
double r=beta*(beta*sigmaSquared-2*diff)/2.0;

yu = (a_minus*sigmaSquared-diff)*someConst;
yv = (a_plus*sigmaSquared-diff)*someConst;
ys = (alpha*sigmaSquared-diff)*someConst;
yr = (beta*sigmaSquared-diff)*someConst;
double yu = (a_minus*sigmaSquared-diff)*someConst;
double yv = (a_plus*sigmaSquared-diff)*someConst;
double ys = (alpha*sigmaSquared-diff)*someConst;
double yr = (beta*sigmaSquared-diff)*someConst;

zs = std::complex<double>(-alpha*diff, 0.5*alpha*gamma);
zu = (1-k)*zs;
zv = (1-k)*zs;
zr = std::complex<double>(-beta*diff, 0.5*beta*gamma);
std::complex<double> zs = std::complex<double>(-alpha*diff, 0.5*alpha*gamma);
std::complex<double> zu = (1-k)*zs;
std::complex<double> zv = (1-k)*zs;
std::complex<double> zr = std::complex<double>(-beta*diff, 0.5*beta*gamma);

N = 0.25*alpha*(1-k*k)/(k*k);
double N = 0.25*alpha*(1-k*k)/(k*k);

out[i] = I*N*( (1-eta)*(Nu*exp(u+gsl_sf_log_erfc(yu))+Nv*exp(v+gsl_sf_log_erfc(yv)) +
Ns*exp(s+gsl_sf_log_erfc(ys))+Nr*exp(r+gsl_sf_log_erfc(yr))) -
Expand Down

0 comments on commit b62c13f

Please sign in to comment.