Changeset 042f82 for src/vector.cpp


Ignore:
Timestamp:
Jul 23, 2009, 2:21:57 PM (16 years ago)
Author:
Frederik Heber <heber@…>
Branches:
Action_Thermostats, Add_AtomRandomPerturbation, Add_FitFragmentPartialChargesAction, Add_RotateAroundBondAction, Add_SelectAtomByNameAction, Added_ParseSaveFragmentResults, AddingActions_SaveParseParticleParameters, Adding_Graph_to_ChangeBondActions, Adding_MD_integration_tests, Adding_ParticleName_to_Atom, Adding_StructOpt_integration_tests, AtomFragments, Automaking_mpqc_open, AutomationFragmentation_failures, Candidate_v1.5.4, Candidate_v1.6.0, Candidate_v1.6.1, Candidate_v1.7.0, ChangeBugEmailaddress, ChangingTestPorts, ChemicalSpaceEvaluator, CombiningParticlePotentialParsing, Combining_Subpackages, Debian_Package_split, Debian_package_split_molecuildergui_only, Disabling_MemDebug, Docu_Python_wait, EmpiricalPotential_contain_HomologyGraph, EmpiricalPotential_contain_HomologyGraph_documentation, Enable_parallel_make_install, Enhance_userguide, Enhanced_StructuralOptimization, Enhanced_StructuralOptimization_continued, Example_ManyWaysToTranslateAtom, Exclude_Hydrogens_annealWithBondGraph, FitPartialCharges_GlobalError, Fix_BoundInBox_CenterInBox_MoleculeActions, Fix_ChargeSampling_PBC, Fix_ChronosMutex, Fix_FitPartialCharges, Fix_FitPotential_needs_atomicnumbers, Fix_ForceAnnealing, Fix_IndependentFragmentGrids, Fix_ParseParticles, Fix_ParseParticles_split_forward_backward_Actions, Fix_PopActions, Fix_QtFragmentList_sorted_selection, Fix_Restrictedkeyset_FragmentMolecule, Fix_StatusMsg, Fix_StepWorldTime_single_argument, Fix_Verbose_Codepatterns, Fix_fitting_potentials, Fixes, ForceAnnealing_goodresults, ForceAnnealing_oldresults, ForceAnnealing_tocheck, ForceAnnealing_with_BondGraph, ForceAnnealing_with_BondGraph_continued, ForceAnnealing_with_BondGraph_continued_betteresults, ForceAnnealing_with_BondGraph_contraction-expansion, FragmentAction_writes_AtomFragments, FragmentMolecule_checks_bonddegrees, GeometryObjects, Gui_Fixes, Gui_displays_atomic_force_velocity, ImplicitCharges, IndependentFragmentGrids, IndependentFragmentGrids_IndividualZeroInstances, IndependentFragmentGrids_IntegrationTest, IndependentFragmentGrids_Sole_NN_Calculation, JobMarket_RobustOnKillsSegFaults, JobMarket_StableWorkerPool, JobMarket_unresolvable_hostname_fix, MoreRobust_FragmentAutomation, ODR_violation_mpqc_open, PartialCharges_OrthogonalSummation, PdbParser_setsAtomName, PythonUI_with_named_parameters, QtGui_reactivate_TimeChanged_changes, Recreated_GuiChecks, Rewrite_FitPartialCharges, RotateToPrincipalAxisSystem_UndoRedo, SaturateAtoms_findBestMatching, SaturateAtoms_singleDegree, StoppableMakroAction, Subpackage_CodePatterns, Subpackage_JobMarket, Subpackage_LinearAlgebra, Subpackage_levmar, Subpackage_mpqc_open, Subpackage_vmg, Switchable_LogView, ThirdParty_MPQC_rebuilt_buildsystem, TrajectoryDependenant_MaxOrder, TremoloParser_IncreasedPrecision, TremoloParser_MultipleTimesteps, TremoloParser_setsAtomName, Ubuntu_1604_changes, stable
Children:
36ec71
Parents:
205ccd
Message:

fixed indentation from tabs to two spaces.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/vector.cpp

    r205ccd r042f82  
    2828double Vector::DistanceSquared(const Vector *y) const
    2929{
    30         double res = 0.;
    31         for (int i=NDIM;i--;)
    32                 res += (x[i]-y->x[i])*(x[i]-y->x[i]);
    33         return (res);
     30  double res = 0.;
     31  for (int i=NDIM;i--;)
     32    res += (x[i]-y->x[i])*(x[i]-y->x[i]);
     33  return (res);
    3434};
    3535
    … …  
    4040double Vector::Distance(const Vector *y) const
    4141{
    42         double res = 0.;
    43         for (int i=NDIM;i--;)
    44                 res += (x[i]-y->x[i])*(x[i]-y->x[i]);
    45         return (sqrt(res));
     42  double res = 0.;
     43  for (int i=NDIM;i--;)
     44    res += (x[i]-y->x[i])*(x[i]-y->x[i]);
     45  return (sqrt(res));
    4646};
    4747
    … …  
    5353double Vector::PeriodicDistance(const Vector *y, const double *cell_size) const
    5454{
    55         double res = Distance(y), tmp, matrix[NDIM*NDIM];
    56         Vector Shiftedy, TranslationVector;
    57         int N[NDIM];
    58         matrix[0] = cell_size[0];
    59         matrix[1] = cell_size[1];
    60         matrix[2] = cell_size[3];
    61         matrix[3] = cell_size[1];
    62         matrix[4] = cell_size[2];
    63         matrix[5] = cell_size[4];
    64         matrix[6] = cell_size[3];
    65         matrix[7] = cell_size[4];
    66         matrix[8] = cell_size[5];
    67         // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
    68         for (N[0]=-1;N[0]<=1;N[0]++)
    69                 for (N[1]=-1;N[1]<=1;N[1]++)
    70                         for (N[2]=-1;N[2]<=1;N[2]++) {
    71                                 // create the translation vector
    72                                 TranslationVector.Zero();
    73                                 for (int i=NDIM;i--;)
    74                                         TranslationVector.x[i] = (double)N[i];
    75                                 TranslationVector.MatrixMultiplication(matrix);
    76                                 // add onto the original vector to compare with
    77                                 Shiftedy.CopyVector(y);
    78                                 Shiftedy.AddVector(&TranslationVector);
    79                                 // get distance and compare with minimum so far
    80                                 tmp = Distance(&Shiftedy);
    81                                 if (tmp < res) res = tmp;
    82                         }
    83         return (res);
     55  double res = Distance(y), tmp, matrix[NDIM*NDIM];
     56  Vector Shiftedy, TranslationVector;
     57  int N[NDIM];
     58  matrix[0] = cell_size[0];
     59  matrix[1] = cell_size[1];
     60  matrix[2] = cell_size[3];
     61  matrix[3] = cell_size[1];
     62  matrix[4] = cell_size[2];
     63  matrix[5] = cell_size[4];
     64  matrix[6] = cell_size[3];
     65  matrix[7] = cell_size[4];
     66  matrix[8] = cell_size[5];
     67  // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
     68  for (N[0]=-1;N[0]<=1;N[0]++)
     69    for (N[1]=-1;N[1]<=1;N[1]++)
     70      for (N[2]=-1;N[2]<=1;N[2]++) {
     71        // create the translation vector
     72        TranslationVector.Zero();
     73        for (int i=NDIM;i--;)
     74          TranslationVector.x[i] = (double)N[i];
     75        TranslationVector.MatrixMultiplication(matrix);
     76        // add onto the original vector to compare with
     77        Shiftedy.CopyVector(y);
     78        Shiftedy.AddVector(&TranslationVector);
     79        // get distance and compare with minimum so far
     80        tmp = Distance(&Shiftedy);
     81        if (tmp < res) res = tmp;
     82      }
     83  return (res);
    8484};
    8585
    … …  
    9191double Vector::PeriodicDistanceSquared(const Vector *y, const double *cell_size) const
    9292{
    93         double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
    94         Vector Shiftedy, TranslationVector;
    95         int N[NDIM];
    96         matrix[0] = cell_size[0];
    97         matrix[1] = cell_size[1];
    98         matrix[2] = cell_size[3];
    99         matrix[3] = cell_size[1];
    100         matrix[4] = cell_size[2];
    101         matrix[5] = cell_size[4];
    102         matrix[6] = cell_size[3];
    103         matrix[7] = cell_size[4];
    104         matrix[8] = cell_size[5];
    105         // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
    106         for (N[0]=-1;N[0]<=1;N[0]++)
    107                 for (N[1]=-1;N[1]<=1;N[1]++)
    108                         for (N[2]=-1;N[2]<=1;N[2]++) {
    109                                 // create the translation vector
    110                                 TranslationVector.Zero();
    111                                 for (int i=NDIM;i--;)
    112                                         TranslationVector.x[i] = (double)N[i];
    113                                 TranslationVector.MatrixMultiplication(matrix);
    114                                 // add onto the original vector to compare with
    115                                 Shiftedy.CopyVector(y);
    116                                 Shiftedy.AddVector(&TranslationVector);
    117                                 // get distance and compare with minimum so far
    118                                 tmp = DistanceSquared(&Shiftedy);
    119                                 if (tmp < res) res = tmp;
    120                         }
    121         return (res);
     93  double res = DistanceSquared(y), tmp, matrix[NDIM*NDIM];
     94  Vector Shiftedy, TranslationVector;
     95  int N[NDIM];
     96  matrix[0] = cell_size[0];
     97  matrix[1] = cell_size[1];
     98  matrix[2] = cell_size[3];
     99  matrix[3] = cell_size[1];
     100  matrix[4] = cell_size[2];
     101  matrix[5] = cell_size[4];
     102  matrix[6] = cell_size[3];
     103  matrix[7] = cell_size[4];
     104  matrix[8] = cell_size[5];
     105  // in order to check the periodic distance, translate one of the vectors into each of the 27 neighbouring cells
     106  for (N[0]=-1;N[0]<=1;N[0]++)
     107    for (N[1]=-1;N[1]<=1;N[1]++)
     108      for (N[2]=-1;N[2]<=1;N[2]++) {
     109        // create the translation vector
     110        TranslationVector.Zero();
     111        for (int i=NDIM;i--;)
     112          TranslationVector.x[i] = (double)N[i];
     113        TranslationVector.MatrixMultiplication(matrix);
     114        // add onto the original vector to compare with
     115        Shiftedy.CopyVector(y);
     116        Shiftedy.AddVector(&TranslationVector);
     117        // get distance and compare with minimum so far
     118        tmp = DistanceSquared(&Shiftedy);
     119        if (tmp < res) res = tmp;
     120      }
     121  return (res);
    122122};
    123123
    … …  
    128128void Vector::KeepPeriodic(ofstream *out, double *matrix)
    129129{
    130 //      int N[NDIM];
    131 //      bool flag = false;
    132         //vector Shifted, TranslationVector;
    133         Vector TestVector;
    134 //      *out << Verbose(1) << "Begin of KeepPeriodic." << endl;
    135 //      *out << Verbose(2) << "Vector is: ";
    136 //      Output(out);
    137 //      *out << endl;
    138         TestVector.CopyVector(this);
    139         TestVector.InverseMatrixMultiplication(matrix);
    140         for(int i=NDIM;i--;) { // correct periodically
    141                 if (TestVector.x[i] < 0) {      // get every coefficient into the interval [0,1)
    142                         TestVector.x[i] += ceil(TestVector.x[i]);
    143                 } else {
    144                         TestVector.x[i] -= floor(TestVector.x[i]);
    145                 }
    146         }
    147         TestVector.MatrixMultiplication(matrix);
    148         CopyVector(&TestVector);
    149 //      *out << Verbose(2) << "New corrected vector is: ";
    150 //      Output(out);
    151 //      *out << endl;
    152 //      *out << Verbose(1) << "End of KeepPeriodic." << endl;
     130//  int N[NDIM];
     131//  bool flag = false;
     132  //vector Shifted, TranslationVector;
     133  Vector TestVector;
     134//  *out << Verbose(1) << "Begin of KeepPeriodic." << endl;
     135//  *out << Verbose(2) << "Vector is: ";
     136//  Output(out);
     137//  *out << endl;
     138  TestVector.CopyVector(this);
     139  TestVector.InverseMatrixMultiplication(matrix);
     140  for(int i=NDIM;i--;) { // correct periodically
     141    if (TestVector.x[i] < 0) {  // get every coefficient into the interval [0,1)
     142      TestVector.x[i] += ceil(TestVector.x[i]);
     143    } else {
     144      TestVector.x[i] -= floor(TestVector.x[i]);
     145    }
     146  }
     147  TestVector.MatrixMultiplication(matrix);
     148  CopyVector(&TestVector);
     149//  *out << Verbose(2) << "New corrected vector is: ";
     150//  Output(out);
     151//  *out << endl;
     152//  *out << Verbose(1) << "End of KeepPeriodic." << endl;
    153153};
    154154
    … …  
    159159double Vector::ScalarProduct(const Vector *y) const
    160160{
    161         double res = 0.;
    162         for (int i=NDIM;i--;)
    163                 res += x[i]*y->x[i];
    164         return (res);
     161  double res = 0.;
     162  for (int i=NDIM;i--;)
     163    res += x[i]*y->x[i];
     164  return (res);
    165165};
    166166
    167167
    168168/** Calculates VectorProduct between this and another vector.
    169  *      -# returns the Product in place of vector from which it was initiated
    170  *      -# ATTENTION: Only three dim.
    171  *      \param *y array to vector with which to calculate crossproduct
    172  *      \return \f$ x \times y \f&
     169 *  -# returns the Product in place of vector from which it was initiated
     170 *  -# ATTENTION: Only three dim.
     171 *  \param *y array to vector with which to calculate crossproduct
     172 *  \return \f$ x \times y \f&
    173173 */
    174174void Vector::VectorProduct(const Vector *y)
    175175{
    176         Vector tmp;
    177         tmp.x[0] = x[1]* (y->x[2]) - x[2]* (y->x[1]);
    178         tmp.x[1] = x[2]* (y->x[0]) - x[0]* (y->x[2]);
    179         tmp.x[2] = x[0]* (y->x[1]) - x[1]* (y->x[0]);
    180         this->CopyVector(&tmp);
     176  Vector tmp;
     177  tmp.x[0] = x[1]* (y->x[2]) - x[2]* (y->x[1]);
     178  tmp.x[1] = x[2]* (y->x[0]) - x[0]* (y->x[2]);
     179  tmp.x[2] = x[0]* (y->x[1]) - x[1]* (y->x[0]);
     180  this->CopyVector(&tmp);
    181181
    182182};
    … …  
    189189void Vector::ProjectOntoPlane(const Vector *y)
    190190{
    191         Vector tmp;
    192         tmp.CopyVector(y);
    193         tmp.Normalize();
    194         tmp.Scale(ScalarProduct(&tmp));
    195         this->SubtractVector(&tmp);
     191  Vector tmp;
     192  tmp.CopyVector(y);
     193  tmp.Normalize();
     194  tmp.Scale(ScalarProduct(&tmp));
     195  this->SubtractVector(&tmp);
    196196};
    197197
    … …  
    202202double Vector::Projection(const Vector *y) const
    203203{
    204         return (ScalarProduct(y));
     204  return (ScalarProduct(y));
    205205};
    206206
    … …  
    210210double Vector::Norm() const
    211211{
    212         double res = 0.;
    213         for (int i=NDIM;i--;)
    214                 res += this->x[i]*this->x[i];
    215         return (sqrt(res));
     212  double res = 0.;
     213  for (int i=NDIM;i--;)
     214    res += this->x[i]*this->x[i];
     215  return (sqrt(res));
    216216};
    217217
    … …  
    220220void Vector::Normalize()
    221221{
    222         double res = 0.;
    223         for (int i=NDIM;i--;)
    224                 res += this->x[i]*this->x[i];
    225         if (fabs(res) > MYEPSILON)
    226                 res = 1./sqrt(res);
    227         Scale(&res);
     222  double res = 0.;
     223  for (int i=NDIM;i--;)
     224    res += this->x[i]*this->x[i];
     225  if (fabs(res) > MYEPSILON)
     226    res = 1./sqrt(res);
     227  Scale(&res);
    228228};
    229229
    … …  
    232232void Vector::Zero()
    233233{
    234         for (int i=NDIM;i--;)
    235                 this->x[i] = 0.;
     234  for (int i=NDIM;i--;)
     235    this->x[i] = 0.;
    236236};
    237237
    … …  
    240240void Vector::One(double one)
    241241{
    242         for (int i=NDIM;i--;)
    243                 this->x[i] = one;
     242  for (int i=NDIM;i--;)
     243    this->x[i] = one;
    244244};
    245245
    … …  
    248248void Vector::Init(double x1, double x2, double x3)
    249249{
    250         x[0] = x1;
    251         x[1] = x2;
    252         x[2] = x3;
     250  x[0] = x1;
     251  x[1] = x2;
     252  x[2] = x3;
    253253};
    254254
    … …  
    266266  if (angle > 1)
    267267    angle = 1;
    268         return acos(angle);
     268  return acos(angle);
    269269};
    270270
    … …  
    275275void Vector::RotateVector(const Vector *axis, const double alpha)
    276276{
    277         Vector a,y;
    278         // normalise this vector with respect to axis
    279         a.CopyVector(this);
    280         a.Scale(Projection(axis));
    281         SubtractVector(&a);
    282         // construct normal vector
    283         y.MakeNormalVector(axis,this);
    284         y.Scale(Norm());
    285         // scale normal vector by sine and this vector by cosine
    286         y.Scale(sin(alpha));
    287         Scale(cos(alpha));
    288         // add scaled normal vector onto this vector
    289         AddVector(&y);
    290         // add part in axis direction
    291         AddVector(&a);
     277  Vector a,y;
     278  // normalise this vector with respect to axis
     279  a.CopyVector(this);
     280  a.Scale(Projection(axis));
     281  SubtractVector(&a);
     282  // construct normal vector
     283  y.MakeNormalVector(axis,this);
     284  y.Scale(Norm());
     285  // scale normal vector by sine and this vector by cosine
     286  y.Scale(sin(alpha));
     287  Scale(cos(alpha));
     288  // add scaled normal vector onto this vector
     289  AddVector(&y);
     290  // add part in axis direction
     291  AddVector(&a);
    292292};
    293293
    … …  
    299299Vector& operator+=(Vector& a, const Vector& b)
    300300{
    301         a.AddVector(&b);
    302         return a;
     301  a.AddVector(&b);
     302  return a;
    303303};
    304304/** factor each component of \a a times a double \a m.
    … …  
    309309Vector& operator*=(Vector& a, const double m)
    310310{
    311         a.Scale(m);
    312         return a;
    313 };
    314 
    315 /** Sums two vectors \a and \b component-wise.
     311  a.Scale(m);
     312  return a;
     313};
     314
     315/** Sums two vectors \a  and \b component-wise.
    316316 * \param a first vector
    317317 * \param b second vector
    … …  
    320320Vector& operator+(const Vector& a, const Vector& b)
    321321{
    322         Vector *x = new Vector;
    323         x->CopyVector(&a);
    324         x->AddVector(&b);
    325         return *x;
     322  Vector *x = new Vector;
     323  x->CopyVector(&a);
     324  x->AddVector(&b);
     325  return *x;
    326326};
    327327
    … …  
    333333Vector& operator*(const Vector& a, const double m)
    334334{
    335         Vector *x = new Vector;
    336         x->CopyVector(&a);
    337         x->Scale(m);
    338         return *x;
     335  Vector *x = new Vector;
     336  x->CopyVector(&a);
     337  x->Scale(m);
     338  return *x;
    339339};
    340340
    … …  
    345345bool Vector::Output(ofstream *out) const
    346346{
    347         if (out != NULL) {
    348                 *out << "(";
    349                 for (int i=0;i<NDIM;i++) {
    350                         *out << x[i];
    351                         if (i != 2)
    352                                 *out << ",";
    353                 }
    354                 *out << ")";
    355                 return true;
    356         } else
    357                 return false;
     347  if (out != NULL) {
     348    *out << "(";
     349    for (int i=0;i<NDIM;i++) {
     350      *out << x[i];
     351      if (i != 2)
     352        *out << ",";
     353    }
     354    *out << ")";
     355    return true;
     356  } else
     357    return false;
    358358};
    359359
    360360ostream& operator<<(ostream& ost,Vector& m)
    361361{
    362         ost << "(";
    363         for (int i=0;i<NDIM;i++) {
    364                 ost << m.x[i];
    365                 if (i != 2)
    366                         ost << ",";
    367         }
    368         ost << ")";
    369         return ost;
     362  ost << "(";
     363  for (int i=0;i<NDIM;i++) {
     364    ost << m.x[i];
     365    if (i != 2)
     366      ost << ",";
     367  }
     368  ost << ")";
     369  return ost;
    370370};
    371371
    … …  
    375375void Vector::Scale(double **factor)
    376376{
    377         for (int i=NDIM;i--;)
    378                 x[i] *= (*factor)[i];
     377  for (int i=NDIM;i--;)
     378    x[i] *= (*factor)[i];
    379379};
    380380
    381381void Vector::Scale(double *factor)
    382382{
    383         for (int i=NDIM;i--;)
    384                 x[i] *= *factor;
     383  for (int i=NDIM;i--;)
     384    x[i] *= *factor;
    385385};
    386386
    387387void Vector::Scale(double factor)
    388388{
    389         for (int i=NDIM;i--;)
    390                 x[i] *= factor;
     389  for (int i=NDIM;i--;)
     390    x[i] *= factor;
    391391};
    392392
    … …  
    396396void Vector::Translate(const Vector *trans)
    397397{
    398         for (int i=NDIM;i--;)
    399                 x[i] += trans->x[i];
     398  for (int i=NDIM;i--;)
     399    x[i] += trans->x[i];
    400400};
    401401
    … …  
    405405void Vector::MatrixMultiplication(double *M)
    406406{
    407         Vector C;
    408         // do the matrix multiplication
    409         C.x[0] = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
    410         C.x[1] = M[1]*x[0]+M[4]*x[1]+M[7]*x[2];
    411         C.x[2] = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
    412         // transfer the result into this
    413         for (int i=NDIM;i--;)
    414                 x[i] = C.x[i];
     407  Vector C;
     408  // do the matrix multiplication
     409  C.x[0] = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
     410  C.x[1] = M[1]*x[0]+M[4]*x[1]+M[7]*x[2];
     411  C.x[2] = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
     412  // transfer the result into this
     413  for (int i=NDIM;i--;)
     414    x[i] = C.x[i];
    415415};
    416416
    … …  
    447447void Vector::InverseMatrixMultiplication(double *A)
    448448{
    449         Vector C;
    450         double B[NDIM*NDIM];
    451         double detA = RDET3(A);
    452         double detAReci;
    453 
    454         // calculate the inverse B
    455         if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
    456                 detAReci = 1./detA;
    457                 B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);            // A_11
    458                 B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);            // A_12
    459                 B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);            // A_13
    460                 B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);            // A_21
    461                 B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);            // A_22
    462                 B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);            // A_23
    463                 B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);            // A_31
    464                 B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);            // A_32
    465                 B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);            // A_33
    466 
    467                 // do the matrix multiplication
    468                 C.x[0] = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
    469                 C.x[1] = B[1]*x[0]+B[4]*x[1]+B[7]*x[2];
    470                 C.x[2] = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
    471                 // transfer the result into this
    472                 for (int i=NDIM;i--;)
    473                         x[i] = C.x[i];
    474         } else {
    475                 cerr << "ERROR: inverse of matrix does not exists!" << endl;
    476         }
     449  Vector C;
     450  double B[NDIM*NDIM];
     451  double detA = RDET3(A);
     452  double detAReci;
     453
     454  // calculate the inverse B
     455  if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
     456    detAReci = 1./detA;
     457    B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
     458    B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
     459    B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
     460    B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
     461    B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
     462    B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
     463    B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
     464    B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
     465    B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
     466
     467    // do the matrix multiplication
     468    C.x[0] = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
     469    C.x[1] = B[1]*x[0]+B[4]*x[1]+B[7]*x[2];
     470    C.x[2] = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
     471    // transfer the result into this
     472    for (int i=NDIM;i--;)
     473      x[i] = C.x[i];
     474  } else {
     475    cerr << "ERROR: inverse of matrix does not exists!" << endl;
     476  }
    477477};
    478478
    … …  
    487487void Vector::LinearCombinationOfVectors(const Vector *x1, const Vector *x2, const Vector *x3, double *factors)
    488488{
    489         for(int i=NDIM;i--;)
    490                 x[i] = factors[0]*x1->x[i] + factors[1]*x2->x[i] + factors[2]*x3->x[i];
     489  for(int i=NDIM;i--;)
     490    x[i] = factors[0]*x1->x[i] + factors[1]*x2->x[i] + factors[2]*x3->x[i];
    491491};
    492492
    … …  
    496496void Vector::Mirror(const Vector *n)
    497497{
    498         double projection;
    499         projection = ScalarProduct(n)/n->ScalarProduct(n);              // remove constancy from n (keep as logical one)
    500         // withdraw projected vector twice from original one
    501         cout << Verbose(1) << "Vector: ";
    502         Output((ofstream *)&cout);
    503         cout << "\t";
    504         for (int i=NDIM;i--;)
    505                 x[i] -= 2.*projection*n->x[i];
    506         cout << "Projected vector: ";
    507         Output((ofstream *)&cout);
    508         cout << endl;
     498  double projection;
     499  projection = ScalarProduct(n)/n->ScalarProduct(n);    // remove constancy from n (keep as logical one)
     500  // withdraw projected vector twice from original one
     501  cout << Verbose(1) << "Vector: ";
     502  Output((ofstream *)&cout);
     503  cout << "\t";
     504  for (int i=NDIM;i--;)
     505    x[i] -= 2.*projection*n->x[i];
     506  cout << "Projected vector: ";
     507  Output((ofstream *)&cout);
     508  cout << endl;
    509509};
    510510
    … …  
    518518bool Vector::MakeNormalVector(const Vector *y1, const Vector *y2, const Vector *y3)
    519519{
    520         Vector x1, x2;
    521 
    522         x1.CopyVector(y1);
    523         x1.SubtractVector(y2);
    524         x2.CopyVector(y3);
    525         x2.SubtractVector(y2);
    526         if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
    527                 cout << Verbose(4) << "Given vectors are linear dependent." << endl;
    528                 return false;
    529         }
    530 //      cout << Verbose(4) << "relative, first plane coordinates:";
    531 //      x1.Output((ofstream *)&cout);
    532 //      cout << endl;
    533 //      cout << Verbose(4) << "second plane coordinates:";
    534 //      x2.Output((ofstream *)&cout);
    535 //      cout << endl;
    536 
    537         this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
    538         this->x[1] = (x1.x[2]*x2.x[0] - x1.x[0]*x2.x[2]);
    539         this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
    540         Normalize();
    541 
    542         return true;
     520  Vector x1, x2;
     521
     522  x1.CopyVector(y1);
     523  x1.SubtractVector(y2);
     524  x2.CopyVector(y3);
     525  x2.SubtractVector(y2);
     526  if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
     527    cout << Verbose(4) << "Given vectors are linear dependent." << endl;
     528    return false;
     529  }
     530//  cout << Verbose(4) << "relative, first plane coordinates:";
     531//  x1.Output((ofstream *)&cout);
     532//  cout << endl;
     533//  cout << Verbose(4) << "second plane coordinates:";
     534//  x2.Output((ofstream *)&cout);
     535//  cout << endl;
     536
     537  this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
     538  this->x[1] = (x1.x[2]*x2.x[0] - x1.x[0]*x2.x[2]);
     539  this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
     540  Normalize();
     541
     542  return true;
    543543};
    544544
    … …  
    554554bool Vector::MakeNormalVector(const Vector *y1, const Vector *y2)
    555555{
    556         Vector x1,x2;
    557         x1.CopyVector(y1);
    558         x2.CopyVector(y2);
    559         Zero();
    560         if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
    561                 cout << Verbose(4) << "Given vectors are linear dependent." << endl;
    562                 return false;
    563         }
    564 //      cout << Verbose(4) << "relative, first plane coordinates:";
    565 //      x1.Output((ofstream *)&cout);
    566 //      cout << endl;
    567 //      cout << Verbose(4) << "second plane coordinates:";
    568 //      x2.Output((ofstream *)&cout);
    569 //      cout << endl;
    570 
    571         this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
    572         this->x[1] = (x1.x[2]*x2.x[0] - x1.x[0]*x2.x[2]);
    573         this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
    574         Normalize();
    575 
    576         return true;
     556  Vector x1,x2;
     557  x1.CopyVector(y1);
     558  x2.CopyVector(y2);
     559  Zero();
     560  if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
     561    cout << Verbose(4) << "Given vectors are linear dependent." << endl;
     562    return false;
     563  }
     564//  cout << Verbose(4) << "relative, first plane coordinates:";
     565//  x1.Output((ofstream *)&cout);
     566//  cout << endl;
     567//  cout << Verbose(4) << "second plane coordinates:";
     568//  x2.Output((ofstream *)&cout);
     569//  cout << endl;
     570
     571  this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
     572  this->x[1] = (x1.x[2]*x2.x[0] - x1.x[0]*x2.x[2]);
     573  this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
     574  Normalize();
     575
     576  return true;
    577577};
    578578
    … …  
    584584bool Vector::MakeNormalVector(const Vector *y1)
    585585{
    586         bool result = false;
    587         Vector x1;
    588         x1.CopyVector(y1);
    589         x1.Scale(x1.Projection(this));
    590         SubtractVector(&x1);
    591         for (int i=NDIM;i--;)
    592                 result = result || (fabs(x[i]) > MYEPSILON);
    593 
    594         return result;
     586  bool result = false;
     587  Vector x1;
     588  x1.CopyVector(y1);
     589  x1.Scale(x1.Projection(this));
     590  SubtractVector(&x1);
     591  for (int i=NDIM;i--;)
     592    result = result || (fabs(x[i]) > MYEPSILON);
     593
     594  return result;
    595595};
    596596
    … …  
    603603bool Vector::GetOneNormalVector(const Vector *GivenVector)
    604604{
    605         int Components[NDIM]; // contains indices of non-zero components
    606         int Last = 0;    // count the number of non-zero entries in vector
    607         int j;  // loop variables
    608         double norm;
    609 
    610         cout << Verbose(4);
    611         GivenVector->Output((ofstream *)&cout);
    612         cout << endl;
    613         for (j=NDIM;j--;)
    614                 Components[j] = -1;
    615         // find two components != 0
    616         for (j=0;j<NDIM;j++)
    617                 if (fabs(GivenVector->x[j]) > MYEPSILON)
    618                         Components[Last++] = j;
    619         cout << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl;
    620 
    621         switch(Last) {
    622                 case 3: // threecomponent system
    623                 case 2: // two component system
    624                         norm = sqrt(1./(GivenVector->x[Components[1]]*GivenVector->x[Components[1]]) + 1./(GivenVector->x[Components[0]]*GivenVector->x[Components[0]]));
    625                         x[Components[2]] = 0.;
    626                         // in skp both remaining parts shall become zero but with opposite sign and third is zero
    627                         x[Components[1]] = -1./GivenVector->x[Components[1]] / norm;
    628                         x[Components[0]] = 1./GivenVector->x[Components[0]] / norm;
    629                         return true;
    630                         break;
    631                 case 1: // one component system
    632                         // set sole non-zero component to 0, and one of the other zero component pendants to 1
    633                         x[(Components[0]+2)%NDIM] = 0.;
    634                         x[(Components[0]+1)%NDIM] = 1.;
    635                         x[Components[0]] = 0.;
    636                         return true;
    637                         break;
    638                 default:
    639                         return false;
    640         }
     605  int Components[NDIM]; // contains indices of non-zero components
     606  int Last = 0;  // count the number of non-zero entries in vector
     607  int j;  // loop variables
     608  double norm;
     609
     610  cout << Verbose(4);
     611  GivenVector->Output((ofstream *)&cout);
     612  cout << endl;
     613  for (j=NDIM;j--;)
     614    Components[j] = -1;
     615  // find two components != 0
     616  for (j=0;j<NDIM;j++)
     617    if (fabs(GivenVector->x[j]) > MYEPSILON)
     618      Components[Last++] = j;
     619  cout << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl;
     620
     621  switch(Last) {
     622    case 3:  // threecomponent system
     623    case 2:  // two component system
     624      norm = sqrt(1./(GivenVector->x[Components[1]]*GivenVector->x[Components[1]]) + 1./(GivenVector->x[Components[0]]*GivenVector->x[Components[0]]));
     625      x[Components[2]] = 0.;
     626      // in skp both remaining parts shall become zero but with opposite sign and third is zero
     627      x[Components[1]] = -1./GivenVector->x[Components[1]] / norm;
     628      x[Components[0]] = 1./GivenVector->x[Components[0]] / norm;
     629      return true;
     630      break;
     631    case 1: // one component system
     632      // set sole non-zero component to 0, and one of the other zero component pendants to 1
     633      x[(Components[0]+2)%NDIM] = 0.;
     634      x[(Components[0]+1)%NDIM] = 1.;
     635      x[Components[0]] = 0.;
     636      return true;
     637      break;
     638    default:
     639      return false;
     640  }
    641641};
    642642
    … …  
    649649double Vector::CutsPlaneAt(Vector *A, Vector *B, Vector *C)
    650650{
    651 //      cout << Verbose(3) << "For comparison: ";
    652 //      cout << "A " << A->Projection(this) << "\t";
    653 //      cout << "B " << B->Projection(this) << "\t";
    654 //      cout << "C " << C->Projection(this) << "\t";
    655 //      cout << endl;
    656         return A->Projection(this);
     651//  cout << Verbose(3) << "For comparison: ";
     652//  cout << "A " << A->Projection(this) << "\t";
     653//  cout << "B " << B->Projection(this) << "\t";
     654//  cout << "C " << C->Projection(this) << "\t";
     655//  cout << endl;
     656  return A->Projection(this);
    657657};
    658658
    … …  
    664664bool Vector::LSQdistance(Vector **vectors, int num)
    665665{
    666         int j;
    667 
    668         for (j=0;j<num;j++) {
    669                 cout << Verbose(1) << j << "th atom's vector: ";
    670                 (vectors[j])->Output((ofstream *)&cout);
    671                 cout << endl;
    672         }
    673 
    674         int np = 3;
    675         struct LSQ_params par;
    676 
    677         const gsl_multimin_fminimizer_type *T =
    678                 gsl_multimin_fminimizer_nmsimplex;
    679         gsl_multimin_fminimizer *s = NULL;
    680         gsl_vector *ss, *y;
    681         gsl_multimin_function minex_func;
    682 
    683         size_t iter = 0, i;
    684         int status;
    685         double size;
    686 
    687         /* Initial vertex size vector */
    688         ss = gsl_vector_alloc (np);
    689         y = gsl_vector_alloc (np);
    690 
    691         /* Set all step sizes to 1 */
    692         gsl_vector_set_all (ss, 1.0);
    693 
    694         /* Starting point */
    695         par.vectors = vectors;
    696         par.num = num;
    697 
    698         for (i=NDIM;i--;)
    699                 gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.);
    700 
    701         /* Initialize method and iterate */
    702         minex_func.f = &LSQ;
    703         minex_func.n = np;
    704         minex_func.params = (void *)&par;
    705 
    706         s = gsl_multimin_fminimizer_alloc (T, np);
    707         gsl_multimin_fminimizer_set (s, &minex_func, y, ss);
    708 
    709         do
    710                 {
    711                         iter++;
    712                         status = gsl_multimin_fminimizer_iterate(s);
    713 
    714                         if (status)
    715                                 break;
    716 
    717                         size = gsl_multimin_fminimizer_size (s);
    718                         status = gsl_multimin_test_size (size, 1e-2);
    719 
    720                         if (status == GSL_SUCCESS)
    721                                 {
    722                                         printf ("converged to minimum at\n");
    723                                 }
    724 
    725                         printf ("%5d ", (int)iter);
    726                         for (i = 0; i < (size_t)np; i++)
    727                                 {
    728                                         printf ("%10.3e ", gsl_vector_get (s->x, i));
    729                                 }
    730                         printf ("f() = %7.3f size = %.3f\n", s->fval, size);
    731                 }
    732         while (status == GSL_CONTINUE && iter < 100);
    733 
    734         for (i=(size_t)np;i--;)
    735                 this->x[i] = gsl_vector_get(s->x, i);
    736         gsl_vector_free(y);
    737         gsl_vector_free(ss);
    738         gsl_multimin_fminimizer_free (s);
    739 
    740         return true;
     666  int j;
     667
     668  for (j=0;j<num;j++) {
     669    cout << Verbose(1) << j << "th atom's vector: ";
     670    (vectors[j])->Output((ofstream *)&cout);
     671    cout << endl;
     672  }
     673
     674  int np = 3;
     675  struct LSQ_params par;
     676
     677  const gsl_multimin_fminimizer_type *T =
     678    gsl_multimin_fminimizer_nmsimplex;
     679  gsl_multimin_fminimizer *s = NULL;
     680  gsl_vector *ss, *y;
     681  gsl_multimin_function minex_func;
     682
     683  size_t iter = 0, i;
     684  int status;
     685  double size;
     686
     687  /* Initial vertex size vector */
     688  ss = gsl_vector_alloc (np);
     689  y = gsl_vector_alloc (np);
     690
     691  /* Set all step sizes to 1 */
     692  gsl_vector_set_all (ss, 1.0);
     693
     694  /* Starting point */
     695  par.vectors = vectors;
     696  par.num = num;
     697
     698  for (i=NDIM;i--;)
     699    gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.);
     700
     701  /* Initialize method and iterate */
     702  minex_func.f = &LSQ;
     703  minex_func.n = np;
     704  minex_func.params = (void *)&par;
     705
     706  s = gsl_multimin_fminimizer_alloc (T, np);
     707  gsl_multimin_fminimizer_set (s, &minex_func, y, ss);
     708
     709  do
     710    {
     711      iter++;
     712      status = gsl_multimin_fminimizer_iterate(s);
     713
     714      if (status)
     715        break;
     716
     717      size = gsl_multimin_fminimizer_size (s);
     718      status = gsl_multimin_test_size (size, 1e-2);
     719
     720      if (status == GSL_SUCCESS)
     721        {
     722          printf ("converged to minimum at\n");
     723        }
     724
     725      printf ("%5d ", (int)iter);
     726      for (i = 0; i < (size_t)np; i++)
     727        {
     728          printf ("%10.3e ", gsl_vector_get (s->x, i));
     729        }
     730      printf ("f() = %7.3f size = %.3f\n", s->fval, size);
     731    }
     732  while (status == GSL_CONTINUE && iter < 100);
     733
     734  for (i=(size_t)np;i--;)
     735    this->x[i] = gsl_vector_get(s->x, i);
     736  gsl_vector_free(y);
     737  gsl_vector_free(ss);
     738  gsl_multimin_fminimizer_free (s);
     739
     740  return true;
    741741};
    742742
    … …  
    746746void Vector::AddVector(const Vector *y)
    747747{
    748         for (int i=NDIM;i--;)
    749                 this->x[i] += y->x[i];
     748  for (int i=NDIM;i--;)
     749    this->x[i] += y->x[i];
    750750}
    751751
    … …  
    755755void Vector::SubtractVector(const Vector *y)
    756756{
    757         for (int i=NDIM;i--;)
    758                 this->x[i] -= y->x[i];
     757  for (int i=NDIM;i--;)
     758    this->x[i] -= y->x[i];
    759759}
    760760
    … …  
    764764void Vector::CopyVector(const Vector *y)
    765765{
    766         for (int i=NDIM;i--;)
    767                 this->x[i] = y->x[i];
     766  for (int i=NDIM;i--;)
     767    this->x[i] = y->x[i];
    768768}
    769769
    … …  
    775775void Vector::AskPosition(double *cell_size, bool check)
    776776{
    777         char coords[3] = {'x','y','z'};
    778         int j = -1;
    779         for (int i=0;i<3;i++) {
    780                 j += i+1;
    781                 do {
    782                         cout << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: ";
    783                         cin >> x[i];
    784                 } while (((x[i] < 0) || (x[i] >= cell_size[j])) && (check));
    785         }
     777  char coords[3] = {'x','y','z'};
     778  int j = -1;
     779  for (int i=0;i<3;i++) {
     780    j += i+1;
     781    do {
     782      cout << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: ";
     783      cin >> x[i];
     784    } while (((x[i] < 0) || (x[i] >= cell_size[j])) && (check));
     785  }
    786786};
    787787
    … …  
    803803bool Vector::SolveSystem(Vector *x1, Vector *x2, Vector *y, double alpha, double beta, double c)
    804804{
    805         double D1,D2,D3,E1,E2,F1,F2,F3,p,q=0., A, B1, B2, C;
    806         double ang; // angle on testing
    807         double sign[3];
    808         int i,j,k;
    809         A = cos(alpha) * x1->Norm() * c;
    810         B1 = cos(beta + M_PI/2.) * y->Norm() * c;
    811         B2 = cos(beta) * x2->Norm() * c;
    812         C = c * c;
    813         cout << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl;
    814         int flag = 0;
    815         if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping
    816                 if (fabs(x1->x[1]) > MYEPSILON) {
    817                         flag = 1;
    818                 } else if (fabs(x1->x[2]) > MYEPSILON) {
    819                         flag = 2;
    820                 } else {
    821                         return false;
    822                 }
    823         }
    824         switch (flag) {
    825                 default:
    826                 case 0:
    827                         break;
    828                 case 2:
    829                         flip(&x1->x[0],&x1->x[1]);
    830                         flip(&x2->x[0],&x2->x[1]);
    831                         flip(&y->x[0],&y->x[1]);
    832                         //flip(&x[0],&x[1]);
    833                         flip(&x1->x[1],&x1->x[2]);
    834                         flip(&x2->x[1],&x2->x[2]);
    835                         flip(&y->x[1],&y->x[2]);
    836                         //flip(&x[1],&x[2]);
    837                 case 1:
    838                         flip(&x1->x[0],&x1->x[1]);
    839                         flip(&x2->x[0],&x2->x[1]);
    840                         flip(&y->x[0],&y->x[1]);
    841                         //flip(&x[0],&x[1]);
    842                         flip(&x1->x[1],&x1->x[2]);
    843                         flip(&x2->x[1],&x2->x[2]);
    844                         flip(&y->x[1],&y->x[2]);
    845                         //flip(&x[1],&x[2]);
    846                         break;
    847         }
    848         // now comes the case system
    849         D1 = -y->x[0]/x1->x[0]*x1->x[1]+y->x[1];
    850         D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2];
    851         D3 = y->x[0]/x1->x[0]*A-B1;
    852         cout << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n";
    853         if (fabs(D1) < MYEPSILON) {
    854                 cout << Verbose(2) << "D1 == 0!\n";
    855                 if (fabs(D2) > MYEPSILON) {
    856                         cout << Verbose(3) << "D2 != 0!\n";
    857                         x[2] = -D3/D2;
    858                         E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2;
    859                         E2 = -x1->x[1]/x1->x[0];
    860                         cout << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n";
    861                         F1 = E1*E1 + 1.;
    862                         F2 = -E1*E2;
    863                         F3 = E1*E1 + D3*D3/(D2*D2) - C;
    864                         cout << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
    865                         if (fabs(F1) < MYEPSILON) {
    866                                 cout << Verbose(4) << "F1 == 0!\n";
    867                                 cout << Verbose(4) << "Gleichungssystem linear\n";
    868                                 x[1] = F3/(2.*F2);
    869                         } else {
    870                                 p = F2/F1;
    871                                 q = p*p - F3/F1;
    872                                 cout << Verbose(4) << "p " << p << "\tq " << q << endl;
    873                                 if (q < 0) {
    874                                         cout << Verbose(4) << "q < 0" << endl;
    875                                         return false;
    876                                 }
    877                                 x[1] = p + sqrt(q);
    878                         }
    879                         x[0] =  A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
    880                 } else {
    881                         cout << Verbose(2) << "Gleichungssystem unterbestimmt\n";
    882                         return false;
    883                 }
    884         } else {
    885                 E1 = A/x1->x[0]+x1->x[1]/x1->x[0]*D3/D1;
    886                 E2 = x1->x[1]/x1->x[0]*D2/D1 - x1->x[2];
    887                 cout << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n";
    888                 F1 = E2*E2 + D2*D2/(D1*D1) + 1.;
    889                 F2 = -(E1*E2 + D2*D3/(D1*D1));
    890                 F3 = E1*E1 + D3*D3/(D1*D1) - C;
    891                 cout << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
    892                 if (fabs(F1) < MYEPSILON) {
    893                         cout << Verbose(3) << "F1 == 0!\n";
    894                         cout << Verbose(3) << "Gleichungssystem linear\n";
    895                         x[2] = F3/(2.*F2);
    896                 } else {
    897                         p = F2/F1;
    898                         q = p*p - F3/F1;
    899                         cout << Verbose(3) << "p " << p << "\tq " << q << endl;
    900                         if (q < 0) {
    901                                 cout << Verbose(3) << "q < 0" << endl;
    902                                 return false;
    903                         }
    904                         x[2] = p + sqrt(q);
    905                 }
    906                 x[1] = (-D2 * x[2] - D3)/D1;
    907                 x[0] = A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
    908         }
    909         switch (flag) { // back-flipping
    910                 default:
    911                 case 0:
    912                         break;
    913                 case 2:
    914                         flip(&x1->x[0],&x1->x[1]);
    915                         flip(&x2->x[0],&x2->x[1]);
    916                         flip(&y->x[0],&y->x[1]);
    917                         flip(&x[0],&x[1]);
    918                         flip(&x1->x[1],&x1->x[2]);
    919                         flip(&x2->x[1],&x2->x[2]);
    920                         flip(&y->x[1],&y->x[2]);
    921                         flip(&x[1],&x[2]);
    922                 case 1:
    923                         flip(&x1->x[0],&x1->x[1]);
    924                         flip(&x2->x[0],&x2->x[1]);
    925                         flip(&y->x[0],&y->x[1]);
    926                         //flip(&x[0],&x[1]);
    927                         flip(&x1->x[1],&x1->x[2]);
    928                         flip(&x2->x[1],&x2->x[2]);
    929                         flip(&y->x[1],&y->x[2]);
    930                         flip(&x[1],&x[2]);
    931                         break;
    932         }
    933         // one z component is only determined by its radius (without sign)
    934         // thus check eight possible sign flips and determine by checking angle with second vector
    935         for (i=0;i<8;i++) {
    936                 // set sign vector accordingly
    937                 for (j=2;j>=0;j--) {
    938                         k = (i & pot(2,j)) << j;
    939                         cout << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl;
    940                         sign[j] = (k == 0) ? 1. : -1.;
    941                 }
    942                 cout << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n";
    943                 // apply sign matrix
    944                 for (j=NDIM;j--;)
    945                         x[j] *= sign[j];
    946                 // calculate angle and check
    947                 ang = x2->Angle (this);
    948                 cout << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t";
    949                 if (fabs(ang - cos(beta)) < MYEPSILON) {
    950                         break;
    951                 }
    952                 // unapply sign matrix (is its own inverse)
    953                 for (j=NDIM;j--;)
    954                         x[j] *= sign[j];
    955         }
    956         return true;
    957 };
     805  double D1,D2,D3,E1,E2,F1,F2,F3,p,q=0., A, B1, B2, C;
     806  double ang; // angle on testing
     807  double sign[3];
     808  int i,j,k;
     809  A = cos(alpha) * x1->Norm() * c;
     810  B1 = cos(beta + M_PI/2.) * y->Norm() * c;
     811  B2 = cos(beta) * x2->Norm() * c;
     812  C = c * c;
     813  cout << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl;
     814  int flag = 0;
     815  if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping
     816    if (fabs(x1->x[1]) > MYEPSILON) {
     817      flag = 1;
     818    } else if (fabs(x1->x[2]) > MYEPSILON) {
     819      flag = 2;
     820    } else {
     821      return false;
     822    }
     823  }
     824  switch (flag) {
     825    default:
     826    case 0:
     827      break;
     828    case 2:
     829      flip(&x1->x[0],&x1->x[1]);
     830      flip(&x2->x[0],&x2->x[1]);
     831      flip(&y->x[0],&y->x[1]);
     832      //flip(&x[0],&x[1]);
     833      flip(&x1->x[1],&x1->x[2]);
     834      flip(&x2->x[1],&x2->x[2]);
     835      flip(&y->x[1],&y->x[2]);
     836      //flip(&x[1],&x[2]);
     837    case 1:
     838      flip(&x1->x[0],&x1->x[1]);
     839      flip(&x2->x[0],&x2->x[1]);
     840      flip(&y->x[0],&y->x[1]);
     841      //flip(&x[0],&x[1]);
     842      flip(&x1->x[1],&x1->x[2]);
     843      flip(&x2->x[1],&x2->x[2]);
     844      flip(&y->x[1],&y->x[2]);
     845      //flip(&x[1],&x[2]);
     846      break;
     847  }
     848  // now comes the case system
     849  D1 = -y->x[0]/x1->x[0]*x1->x[1]+y->x[1];
     850  D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2];
     851  D3 = y->x[0]/x1->x[0]*A-B1;
     852  cout << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n";
     853  if (fabs(D1) < MYEPSILON) {
     854    cout << Verbose(2) << "D1 == 0!\n";
     855    if (fabs(D2) > MYEPSILON) {
     856      cout << Verbose(3) << "D2 != 0!\n";
     857      x[2] = -D3/D2;
     858      E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2;
     859      E2 = -x1->x[1]/x1->x[0];
     860      cout << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n";
     861      F1 = E1*E1 + 1.;
     862      F2 = -E1*E2;
     863      F3 = E1*E1 + D3*D3/(D2*D2) - C;
     864      cout << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
     865      if (fabs(F1) < MYEPSILON) {
     866        cout << Verbose(4) << "F1 == 0!\n";
     867        cout << Verbose(4) << "Gleichungssystem linear\n";
     868        x[1] = F3/(2.*F2);
     869      } else {
     870        p = F2/F1;
     871        q = p*p - F3/F1;
     872        cout << Verbose(4) << "p " << p << "\tq " << q << endl;
     873        if (q < 0) {
     874          cout << Verbose(4) << "q < 0" << endl;
     875          return false;
     876        }
     877        x[1] = p + sqrt(q);
     878      }
     879      x[0] =  A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
     880    } else {
     881      cout << Verbose(2) << "Gleichungssystem unterbestimmt\n";
     882      return false;
     883    }
     884  } else {
     885    E1 = A/x1->x[0]+x1->x[1]/x1->x[0]*D3/D1;
     886    E2 = x1->x[1]/x1->x[0]*D2/D1 - x1->x[2];
     887    cout << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n";
     888    F1 = E2*E2 + D2*D2/(D1*D1) + 1.;
     889    F2 = -(E1*E2 + D2*D3/(D1*D1));
     890    F3 = E1*E1 + D3*D3/(D1*D1) - C;
     891    cout << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
     892    if (fabs(F1) < MYEPSILON) {
     893      cout << Verbose(3) << "F1 == 0!\n";
     894      cout << Verbose(3) << "Gleichungssystem linear\n";
     895      x[2] = F3/(2.*F2);
     896    } else {
     897      p = F2/F1;
     898      q = p*p - F3/F1;
     899      cout << Verbose(3) << "p " << p << "\tq " << q << endl;
     900      if (q < 0) {
     901        cout << Verbose(3) << "q < 0" << endl;
     902        return false;
     903      }
     904      x[2] = p + sqrt(q);
     905    }
     906    x[1] = (-D2 * x[2] - D3)/D1;
     907    x[0] = A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
     908  }
     909  switch (flag) { // back-flipping
     910    default:
     911    case 0:
     912      break;
     913    case 2:
     914      flip(&x1->x[0],&x1->x[1]);
     915      flip(&x2->x[0],&x2->x[1]);
     916      flip(&y->x[0],&y->x[1]);
     917      flip(&x[0],&x[1]);
     918      flip(&x1->x[1],&x1->x[2]);
     919      flip(&x2->x[1],&x2->x[2]);
     920      flip(&y->x[1],&y->x[2]);
     921      flip(&x[1],&x[2]);
     922    case 1:
     923      flip(&x1->x[0],&x1->x[1]);
     924      flip(&x2->x[0],&x2->x[1]);
     925      flip(&y->x[0],&y->x[1]);
     926      //flip(&x[0],&x[1]);
     927      flip(&x1->x[1],&x1->x[2]);
     928      flip(&x2->x[1],&x2->x[2]);
     929      flip(&y->x[1],&y->x[2]);
     930      flip(&x[1],&x[2]);
     931      break;
     932  }
     933  // one z component is only determined by its radius (without sign)
     934  // thus check eight possible sign flips and determine by checking angle with second vector
     935  for (i=0;i<8;i++) {
     936    // set sign vector accordingly
     937    for (j=2;j>=0;j--) {
     938      k = (i & pot(2,j)) << j;
     939      cout << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl;
     940      sign[j] = (k == 0) ? 1. : -1.;
     941    }
     942    cout << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n";
     943    // apply sign matrix
     944    for (j=NDIM;j--;)
     945      x[j] *= sign[j];
     946    // calculate angle and check
     947    ang = x2->Angle (this);
     948    cout << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t";
     949    if (fabs(ang - cos(beta)) < MYEPSILON) {
     950      break;
     951    }
     952    // unapply sign matrix (is its own inverse)
     953    for (j=NDIM;j--;)
     954      x[j] *= sign[j];
     955  }
     956  return true;
     957};
Note: See TracChangeset for help on using the changeset viewer.