Commit c2250502 authored by cortes_c's avatar cortes_c
Browse files

I added the sorting of the Eigenvectors in SigmaGenerator.h. I also...

I added the sorting of the Eigenvectors in SigmaGenerator.h. I also implemented the renormalization method to ensure numerical consistence in Distribution.cpp. The code is working now, althought I dont know if the permutation plays a role in the matched distribution that is calculated
parent 848043d3
......@@ -1384,7 +1384,7 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
CyclotronElement->getFieldMapFN(),
Attributes::getReal(itsAttr[Attrib::Distribution::ORDERMAPS]),
CyclotronElement->getBScale(),
writeMap);$
writeMap);
if (siggen->match(accuracy,
Attributes::getReal(itsAttr[Attrib::Distribution::MAXSTEPSSI]),
......@@ -1406,7 +1406,13 @@ void Distribution::createMatchedGaussDistribution(size_t numberOfParticles, doub
for (unsigned int i = 0; i < siggen->getSigma().size1(); ++ i) {
*gmsg << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,0);
for (unsigned int j = 1; j < siggen->getSigma().size2(); ++ j) {
*gmsg << " & " << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,j);
if (siggen->getSigma()(i,j) < 10e-12){
*gmsg << " & " << std::setprecision(4) << std::setw(11) << 0.0;
}
else{
*gmsg << " & " << std::setprecision(4) << std::setw(11) << siggen->getSigma()(i,j);
}
}
*gmsg << " \\\\" << endl;
}
......@@ -2535,9 +2541,35 @@ void Distribution::generateGaussZ(size_t numberOfParticles) {
*gmsg << " \\\\" << endl;
}
#endif
//Sets the GSL error handler off, exception will be handled internaly with a renormalization method
gsl_set_error_handler_off();
int errcode = gsl_linalg_cholesky_decomp(corMat);
double rn = 1e-12;
while (errcode == GSL_EDOM) {
// Resets the correlation matrix
for (unsigned int i = 0; i < 6; ++ i) {
gsl_matrix_set(corMat, i, i, correlationMatrix_m(i, i));
for (unsigned int j = 0; j < i; ++ j) {
gsl_matrix_set(corMat, i, j, correlationMatrix_m(i, j));
gsl_matrix_set(corMat, j, i, correlationMatrix_m(i, j));
}
}
// Applying a renormalization method corMat = corMat + rn*Unitymatrix
// This is the renormalization
for(int i = 0; i < 6; i++){
double corMati = gsl_matrix_get(corMat, i , i);
gsl_matrix_set(corMat, i, i, corMati + rn);
}
//Just to be sure that the renormalization worked
errcode = gsl_linalg_cholesky_decomp(corMat);
if(errcode != GSL_EDOM) *gmsg << "* WARNING: Correlation matrix had to be renormalized"<<endl;
else rn *= 10;
}
//Sets again the standard GSL error handler on
gsl_set_error_handler(NULL);
//Just to be sure
if (errcode == GSL_EDOM) {
throw OpalException("Distribution::GenerateGaussZ",
"gsl_linalg_cholesky_decomp failed");
......@@ -3525,7 +3557,7 @@ void Distribution::setAttributes() {
itsAttr[Attrib::Distribution::MAXSTEPSCO]
= Attributes::makeReal("MAXSTEPSCO", "Maximum steps used to find closed orbit ", 100);
itsAttr[Attrib::Distribution::MAXSTEPSSI]
= Attributes::makeReal("MAXSTEPSSI", "Maximum steps used to find matched distribution ",2000);
= Attributes::makeReal("MAXSTEPSSI", "Maximum steps used to find matched distribution ",500);
itsAttr[Attrib::Distribution::ORDERMAPS]
= Attributes::makeReal("ORDERMAPS", "Order used in the field expansion ", 7);
itsAttr[Attrib::Distribution::SECTOR]
......
......@@ -499,12 +499,12 @@ template<typename Value_type, typename Size_type>
// write to terminal
*gmsg << "* ----------------------------" << endl
<< "* Closed orbit info (Gordon units):" << endl
<< "* Closed orbit info:" << endl
<< "*" << endl
<< "* average radius: " << ravg << " [m]" << endl
<< "* initial radius: " << r_turn[0] << " [m]" << endl
<< "* initial momentum: " << peo[0] << " [Beta Gamma]" << endl
<< "* frequency error: " << cof.getFrequencyError() << endl
<< "* frequency error: " << cof.getFrequencyError()*100 <<" [ % ] "<< endl
<< "* horizontal tune: " << tunes.first << endl
<< "* vertical tune: " << tunes.second << endl
<< "* ----------------------------" << endl << endl;
......@@ -552,7 +552,7 @@ template<typename Value_type, typename Size_type>
Mcycs[i] = mapgen.generateMap(H_m(h[i],
h[i]*h[i]+fidx[i],
-fidx[i]),
ds[i],truncOrder_m);
ds[i],truncOrder_m);
Mscs[i] = mapgen.generateMap(Hsc_m(sigmas_m[i](0,0),
sigmas_m[i](2,2),
......@@ -699,6 +699,7 @@ void SigmaGenerator<Value_type, Size_type>::eigsolve_m(const matrix_type& Mturn,
typedef gsl_vector_complex* gsl_cvector_t;
typedef gsl_matrix_complex* gsl_cmatrix_t;
typedef gsl_eigen_nonsymmv_workspace* gsl_wspace_t;
typedef boost::numeric::ublas::vector<complex_t> complex_vector_type;
gsl_cvector_t evals = gsl_vector_complex_alloc(6);
gsl_cmatrix_t evecs = gsl_matrix_complex_alloc(6, 6);
......@@ -718,7 +719,7 @@ void SigmaGenerator<Value_type, Size_type>::eigsolve_m(const matrix_type& Mturn,
// throw OpalException("SigmaGenerator::eigsolve_m()",
// "Couldn't perform eigendecomposition!");
/*status = */gsl_eigen_nonsymmv_sort(evals, evecs, GSL_EIGEN_SORT_ABS_ASC);
/*status = *///gsl_eigen_nonsymmv_sort(evals, evecs, GSL_EIGEN_SORT_ABS_ASC);
// if ( !status )
// throw OpalException("SigmaGenerator::eigsolve_m()",
......@@ -735,6 +736,61 @@ void SigmaGenerator<Value_type, Size_type>::eigsolve_m(const matrix_type& Mturn,
}
}
// Sorting the Eigenvectors
// This is an arbitrary threshold that has worked for me. (We should fix this)
value_type threshold = 10e-12;
bool isZdirection = false;
std::vector<complex_vector_type> zVectors{};
std::vector<complex_vector_type> xyVectors{};
for(size_type i = 0; i < 6; i++){
complex_t z = R(i,0);
if(std::abs(z) < threshold) z = 0.;
if(z == 0.) isZdirection = true;
complex_vector_type v(6);
if(isZdirection){
for(size_type j = 0;j < 6; j++){
complex_t z = R(i,j);
v(j) = z;
}
zVectors.push_back(v);
}
else{
for(size_type j = 0; j < 6; j++){
complex_t z = R(i,j);
v(j) = z;
}
xyVectors.push_back(v);
}
isZdirection = false;
}
//if z-direction not found, then the system does not behave as expected
if(zVectors.size() != 2)
throw OpalException("SigmaGenerator::eigsolve_m()",
"Couldn't find the vertical Eigenvectors.");
// Norms the Eigenvectors
for(size_type i = 0; i < 4; i++){
value_type norm{0};
for(size_type j = 0; j < 6; j++) norm += std::norm(xyVectors[i](j));
for(size_type j = 0; j < 6; j++) xyVectors[i](j) /= std::sqrt(norm);
}
for(size_type i = 0; i < 2; i++){
value_type norm{0};
for(size_type j = 0; j < 6; j++) norm += std::norm(zVectors[i](j));
for(size_type j = 0; j < 6; j++) zVectors[i](j) /= std::sqrt(norm);
}
for(value_type i = 0; i < 6; i++){
R(i,0) = xyVectors[0](i);
R(i,1) = xyVectors[1](i);
R(i,2) = zVectors[0](i);
R(i,3) = zVectors[1](i);
R(i,4) = xyVectors[2](i);
R(i,5) = xyVectors[3](i);
}
gsl_vector_complex_free(evals);
gsl_matrix_complex_free(evecs);
gsl_eigen_nonsymmv_free(wspace);
......
This diff is collapsed.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment