Changeset 6ac7ee for src/vector.cpp


Ignore:
Timestamp:
Feb 9, 2009, 5:24:10 PM (17 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:
d8b94a
Parents:
124df1
git-author:
Frederik Heber <heber@…> (02/09/09 15:55:37)
git-committer:
Frederik Heber <heber@…> (02/09/09 17:24:10)
Message:

Merge branch 'ConcaveHull' of ../espack2 into ConcaveHull

Conflicts:

molecuilder/src/boundary.cpp
molecuilder/src/boundary.hpp
molecuilder/src/builder.cpp
molecuilder/src/linkedcell.cpp
molecuilder/src/linkedcell.hpp
molecuilder/src/vector.cpp
molecuilder/src/vector.hpp
util/src/NanoCreator.c

Basically, this resulted from a lot of conversions two from spaces to one tab, which is my standard indentation. The mess was caused by eclipse auto-indenting. And in espack2:ConcaveHull was the new stuff, so all from ConcaveHull was replaced in case of doubt.
Additionally, vector had ofstream << operator instead ostream << ...

File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/vector.cpp

    • Property mode changed from 100644 to 100755
    r124df1 r6ac7ee  
    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));
     46};
     47
     48/** Calculates distance between this and another vector in a periodic cell.
     49 * \param *y array to second vector
     50 * \param *cell_size 6-dimensional array with (xx, xy, yy, xz, yz, zz) entries specifying the periodic cell
     51 * \return \f$| x - y |\f$
     52 */
     53double Vector::PeriodicDistance(const Vector *y, const double *cell_size) const
     54{
     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);
    4684};
    4785
    … …  
    5189 * \return \f$| x - y |^2\f$
    5290 */
    53 double Vector::PeriodicDistance(const Vector *y, const double *cell_size) const
    54 {
    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);
     91double Vector::PeriodicDistanceSquared(const Vector *y, const double *cell_size) const
     92{
     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);
    84122};
    85123
    … …  
    90128void Vector::KeepPeriodic(ofstream *out, double *matrix)
    91129{
    92 //  int N[NDIM];
    93 //  bool flag = false;
    94   //vector Shifted, TranslationVector;
    95   Vector TestVector;
    96 //  *out << Verbose(1) << "Begin of KeepPeriodic." << endl;
    97 //  *out << Verbose(2) << "Vector is: ";
    98 //  Output(out);
    99 //  *out << endl;
    100   TestVector.CopyVector(this);
    101   TestVector.InverseMatrixMultiplication(matrix);
    102   for(int i=NDIM;i--;) { // correct periodically
    103     if (TestVector.x[i] < 0) {  // get every coefficient into the interval [0,1)
    104       TestVector.x[i] += ceil(TestVector.x[i]);
    105     } else {
    106       TestVector.x[i] -= floor(TestVector.x[i]);
    107     }
    108   }
    109   TestVector.MatrixMultiplication(matrix);
    110   CopyVector(&TestVector);
    111 //  *out << Verbose(2) << "New corrected vector is: ";
    112 //  Output(out);
    113 //  *out << endl;
    114 //  *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;
    115153};
    116154
    … …  
    121159double Vector::ScalarProduct(const Vector *y) const
    122160{
    123   double res = 0.;
    124   for (int i=NDIM;i--;)
    125     res += x[i]*y->x[i];
    126   return (res);
     161        double res = 0.;
     162        for (int i=NDIM;i--;)
     163                res += x[i]*y->x[i];
     164        return (res);
    127165};
    128166
    129167
    130168/** Calculates VectorProduct between this and another vector.
    131  *  -# returns the Product in place of vector from which it was initiated
    132  *  -# ATTENTION: Only three dim.
    133  *  \param *y array to vector with which to calculate crossproduct
    134  *  \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&
    135173 */
    136174void Vector::VectorProduct(const Vector *y)
    137175{
    138   Vector tmp;
    139   tmp.x[0] = x[1]* (y->x[2]) - x[2]* (y->x[1]);
    140   tmp.x[1] = x[2]* (y->x[0]) - x[0]* (y->x[2]);
    141   tmp.x[2] = x[0]* (y->x[1]) - x[1]* (y->x[0]);
    142   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);
    143181
    144182};
    … …  
    151189void Vector::ProjectOntoPlane(const Vector *y)
    152190{
    153   Vector tmp;
    154   tmp.CopyVector(y);
    155   tmp.Normalize();
    156   tmp.Scale(ScalarProduct(&tmp));
    157   this->SubtractVector(&tmp);
     191        Vector tmp;
     192        tmp.CopyVector(y);
     193        tmp.Normalize();
     194        tmp.Scale(ScalarProduct(&tmp));
     195        this->SubtractVector(&tmp);
    158196};
    159197
    … …  
    164202double Vector::Projection(const Vector *y) const
    165203{
    166   return (ScalarProduct(y));
     204        return (ScalarProduct(y));
    167205};
    168206
    … …  
    172210double Vector::Norm() const
    173211{
    174   double res = 0.;
    175   for (int i=NDIM;i--;)
    176     res += this->x[i]*this->x[i];
    177   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));
    178216};
    179217
    … …  
    182220void Vector::Normalize()
    183221{
    184   double res = 0.;
    185   for (int i=NDIM;i--;)
    186     res += this->x[i]*this->x[i];
    187   if (fabs(res) > MYEPSILON)
    188     res = 1./sqrt(res);
    189   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);
    190228};
    191229
    … …  
    194232void Vector::Zero()
    195233{
    196   for (int i=NDIM;i--;)
    197     this->x[i] = 0.;
     234        for (int i=NDIM;i--;)
     235                this->x[i] = 0.;
    198236};
    199237
    … …  
    202240void Vector::One(double one)
    203241{
    204   for (int i=NDIM;i--;)
    205     this->x[i] = one;
     242        for (int i=NDIM;i--;)
     243                this->x[i] = one;
    206244};
    207245
    … …  
    210248void Vector::Init(double x1, double x2, double x3)
    211249{
    212   x[0] = x1;
    213   x[1] = x2;
    214   x[2] = x3;
     250        x[0] = x1;
     251        x[1] = x2;
     252        x[2] = x3;
    215253};
    216254
    … …  
    219257 * \return \f$\acos\bigl(frac{\langle x, y \rangle}{|x||y|}\bigr)\f$
    220258 */
    221 double Vector::Angle(Vector *y) const
    222 {
    223   return acos(this->ScalarProduct(y)/Norm()/y->Norm());
     259double Vector::Angle(const Vector *y) const
     260{
     261        return acos(this->ScalarProduct(y)/Norm()/y->Norm());
    224262};
    225263
    … …  
    230268void Vector::RotateVector(const Vector *axis, const double alpha)
    231269{
    232   Vector a,y;
    233   // normalise this vector with respect to axis
    234   a.CopyVector(this);
    235   a.Scale(Projection(axis));
    236   SubtractVector(&a);
    237   // construct normal vector
    238   y.MakeNormalVector(axis,this);
    239   y.Scale(Norm());
    240   // scale normal vector by sine and this vector by cosine
    241   y.Scale(sin(alpha));
    242   Scale(cos(alpha));
    243   // add scaled normal vector onto this vector
    244   AddVector(&y);
    245   // add part in axis direction
    246   AddVector(&a);
     270        Vector a,y;
     271        // normalise this vector with respect to axis
     272        a.CopyVector(this);
     273        a.Scale(Projection(axis));
     274        SubtractVector(&a);
     275        // construct normal vector
     276        y.MakeNormalVector(axis,this);
     277        y.Scale(Norm());
     278        // scale normal vector by sine and this vector by cosine
     279        y.Scale(sin(alpha));
     280        Scale(cos(alpha));
     281        // add scaled normal vector onto this vector
     282        AddVector(&y);
     283        // add part in axis direction
     284        AddVector(&a);
    247285};
    248286
    … …  
    254292Vector& operator+=(Vector& a, const Vector& b)
    255293{
    256   a.AddVector(&b);
    257   return a;
     294        a.AddVector(&b);
     295        return a;
    258296};
    259297/** factor each component of \a a times a double \a m.
    … …  
    264302Vector& operator*=(Vector& a, const double m)
    265303{
    266   a.Scale(m);
    267   return a;
    268 };
    269 
    270 /** Sums two vectors \a  and \b component-wise.
     304        a.Scale(m);
     305        return a;
     306};
     307
     308/** Sums two vectors \a and \b component-wise.
    271309 * \param a first vector
    272310 * \param b second vector
    … …  
    275313Vector& operator+(const Vector& a, const Vector& b)
    276314{
    277   Vector *x = new Vector;
    278   x->CopyVector(&a);
    279   x->AddVector(&b);
    280   return *x;
     315        Vector *x = new Vector;
     316        x->CopyVector(&a);
     317        x->AddVector(&b);
     318        return *x;
    281319};
    282320
    … …  
    288326Vector& operator*(const Vector& a, const double m)
    289327{
    290   Vector *x = new Vector;
    291   x->CopyVector(&a);
    292   x->Scale(m);
    293   return *x;
     328        Vector *x = new Vector;
     329        x->CopyVector(&a);
     330        x->Scale(m);
     331        return *x;
    294332};
    295333
    … …  
    300338bool Vector::Output(ofstream *out) const
    301339{
    302   if (out != NULL) {
    303     *out << "(";
    304     for (int i=0;i<NDIM;i++) {
    305       *out << x[i];
    306       if (i != 2)
    307         *out << ",";
    308     }
    309     *out << ")";
    310     return true;
    311   } else
    312     return false;
    313 };
    314 
    315 ofstream& operator<<(ofstream& ost,Vector& m)
    316 {
    317         m.Output(&ost);
     340        if (out != NULL) {
     341                *out << "(";
     342                for (int i=0;i<NDIM;i++) {
     343                        *out << x[i];
     344                        if (i != 2)
     345                                *out << ",";
     346                }
     347                *out << ")";
     348                return true;
     349        } else
     350                return false;
     351};
     352
     353ostream& operator<<(ostream& ost,Vector& m)
     354{
     355        ost << "(";
     356        for (int i=0;i<NDIM;i++) {
     357                ost << m.x[i];
     358                if (i != 2)
     359                        ost << ",";
     360        }
     361        ost << ")";
    318362        return ost;
    319363};
    … …  
    324368void Vector::Scale(double **factor)
    325369{
    326   for (int i=NDIM;i--;)
    327     x[i] *= (*factor)[i];
     370        for (int i=NDIM;i--;)
     371                x[i] *= (*factor)[i];
    328372};
    329373
    330374void Vector::Scale(double *factor)
    331375{
    332   for (int i=NDIM;i--;)
    333     x[i] *= *factor;
     376        for (int i=NDIM;i--;)
     377                x[i] *= *factor;
    334378};
    335379
    336380void Vector::Scale(double factor)
    337381{
    338   for (int i=NDIM;i--;)
    339     x[i] *= factor;
     382        for (int i=NDIM;i--;)
     383                x[i] *= factor;
    340384};
    341385
    … …  
    345389void Vector::Translate(const Vector *trans)
    346390{
    347   for (int i=NDIM;i--;)
    348     x[i] += trans->x[i];
     391        for (int i=NDIM;i--;)
     392                x[i] += trans->x[i];
    349393};
    350394
    … …  
    354398void Vector::MatrixMultiplication(double *M)
    355399{
    356   Vector C;
    357   // do the matrix multiplication
    358   C.x[0] = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
    359   C.x[1] = M[1]*x[0]+M[4]*x[1]+M[7]*x[2];
    360   C.x[2] = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
    361   // transfer the result into this
    362   for (int i=NDIM;i--;)
    363     x[i] = C.x[i];
     400        Vector C;
     401        // do the matrix multiplication
     402        C.x[0] = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
     403        C.x[1] = M[1]*x[0]+M[4]*x[1]+M[7]*x[2];
     404        C.x[2] = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
     405        // transfer the result into this
     406        for (int i=NDIM;i--;)
     407                x[i] = C.x[i];
    364408};
    365409
    … …  
    369413void Vector::InverseMatrixMultiplication(double *A)
    370414{
    371   Vector C;
    372   double B[NDIM*NDIM];
    373   double detA = RDET3(A);
    374   double detAReci;
    375 
    376   // calculate the inverse B
    377   if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
    378     detAReci = 1./detA;
    379     B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
    380     B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
    381     B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
    382     B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
    383     B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
    384     B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
    385     B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
    386     B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
    387     B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
    388 
    389     // do the matrix multiplication
    390     C.x[0] = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
    391     C.x[1] = B[1]*x[0]+B[4]*x[1]+B[7]*x[2];
    392     C.x[2] = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
    393     // transfer the result into this
    394     for (int i=NDIM;i--;)
    395       x[i] = C.x[i];
    396   } else {
    397     cerr << "ERROR: inverse of matrix does not exists!" << endl;
    398   }
     415        Vector C;
     416        double B[NDIM*NDIM];
     417        double detA = RDET3(A);
     418        double detAReci;
     419
     420        // calculate the inverse B
     421        if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
     422                detAReci = 1./detA;
     423                B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);            // A_11
     424                B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);            // A_12
     425                B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);            // A_13
     426                B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);            // A_21
     427                B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);            // A_22
     428                B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);            // A_23
     429                B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);            // A_31
     430                B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);            // A_32
     431                B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);            // A_33
     432
     433                // do the matrix multiplication
     434                C.x[0] = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
     435                C.x[1] = B[1]*x[0]+B[4]*x[1]+B[7]*x[2];
     436                C.x[2] = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
     437                // transfer the result into this
     438                for (int i=NDIM;i--;)
     439                        x[i] = C.x[i];
     440        } else {
     441                cerr << "ERROR: inverse of matrix does not exists!" << endl;
     442        }
    399443};
    400444
    … …  
    409453void Vector::LinearCombinationOfVectors(const Vector *x1, const Vector *x2, const Vector *x3, double *factors)
    410454{
    411   for(int i=NDIM;i--;)
    412     x[i] = factors[0]*x1->x[i] + factors[1]*x2->x[i] + factors[2]*x3->x[i];
     455        for(int i=NDIM;i--;)
     456                x[i] = factors[0]*x1->x[i] + factors[1]*x2->x[i] + factors[2]*x3->x[i];
    413457};
    414458
    … …  
    418462void Vector::Mirror(const Vector *n)
    419463{
    420   double projection;
    421   projection = ScalarProduct(n)/n->ScalarProduct(n);    // remove constancy from n (keep as logical one)
    422   // withdraw projected vector twice from original one
    423   cout << Verbose(1) << "Vector: ";
    424   Output((ofstream *)&cout);
    425   cout << "\t";
    426   for (int i=NDIM;i--;)
    427     x[i] -= 2.*projection*n->x[i];
    428   cout << "Projected vector: ";
    429   Output((ofstream *)&cout);
    430   cout << endl;
     464        double projection;
     465        projection = ScalarProduct(n)/n->ScalarProduct(n);              // remove constancy from n (keep as logical one)
     466        // withdraw projected vector twice from original one
     467        cout << Verbose(1) << "Vector: ";
     468        Output((ofstream *)&cout);
     469        cout << "\t";
     470        for (int i=NDIM;i--;)
     471                x[i] -= 2.*projection*n->x[i];
     472        cout << "Projected vector: ";
     473        Output((ofstream *)&cout);
     474        cout << endl;
    431475};
    432476
    … …  
    440484bool Vector::MakeNormalVector(const Vector *y1, const Vector *y2, const Vector *y3)
    441485{
    442   Vector x1, x2;
    443 
    444   x1.CopyVector(y1);
    445   x1.SubtractVector(y2);
    446   x2.CopyVector(y3);
    447   x2.SubtractVector(y2);
    448   if ((x1.Norm()==0) || (x2.Norm()==0)) {
    449     cout << Verbose(4) << "Given vectors are linear dependent." << endl;
    450     return false;
    451   }
    452 //  cout << Verbose(4) << "relative, first plane coordinates:";
    453 //  x1.Output((ofstream *)&cout);
    454 //  cout << endl;
    455 //  cout << Verbose(4) << "second plane coordinates:";
    456 //  x2.Output((ofstream *)&cout);
    457 //  cout << endl;
    458 
    459   this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
    460   this->x[1] = (x1.x[2]*x2.x[0] - x1.x[0]*x2.x[2]);
    461   this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
    462   Normalize();
    463 
    464   return true;
     486        Vector x1, x2;
     487
     488        x1.CopyVector(y1);
     489        x1.SubtractVector(y2);
     490        x2.CopyVector(y3);
     491        x2.SubtractVector(y2);
     492        if ((fabs(x1.Norm()) < MYEPSILON) || (fabs(x2.Norm()) < MYEPSILON) || (fabs(x1.Angle(&x2)) < MYEPSILON)) {
     493                cout << Verbose(4) << "Given vectors are linear dependent." << endl;
     494                return false;
     495        }
     496//      cout << Verbose(4) << "relative, first plane coordinates:";
     497//      x1.Output((ofstream *)&cout);
     498//      cout << endl;
     499//      cout << Verbose(4) << "second plane coordinates:";
     500//      x2.Output((ofstream *)&cout);
     501//      cout << endl;
     502
     503        this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
     504        this->x[1] = (x1.x[2]*x2.x[0] - x1.x[0]*x2.x[2]);
     505        this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
     506        Normalize();
     507
     508        return true;
    465509};
    466510
    … …  
    476520bool Vector::MakeNormalVector(const Vector *y1, const Vector *y2)
    477521{
    478   Vector x1,x2;
    479   x1.CopyVector(y1);
    480   x2.CopyVector(y2);
    481   Zero();
    482   if ((x1.Norm()==0) || (x2.Norm()==0)) {
    483     cout << Verbose(4) << "Given vectors are linear dependent." << endl;
    484     return false;
    485   }
    486 //  cout << Verbose(4) << "relative, first plane coordinates:";
    487 //  x1.Output((ofstream *)&cout);
    488 //  cout << endl;
    489 //  cout << Verbose(4) << "second plane coordinates:";
    490 //  x2.Output((ofstream *)&cout);
    491 //  cout << endl;
    492 
    493   this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
    494   this->x[1] = (x1.x[2]*x2.x[0] - x1.x[0]*x2.x[2]);
    495   this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
    496   Normalize();
    497 
    498   return true;
     522        Vector x1,x2;
     523        x1.CopyVector(y1);
     524        x2.CopyVector(y2);
     525        Zero();
     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;
    499543};
    500544
    … …  
    507551{
    508552        bool result = false;
    509   Vector x1;
    510   x1.CopyVector(y1);
    511   x1.Scale(x1.Projection(this));
    512   SubtractVector(&x1);
    513   for (int i=NDIM;i--;)
    514           result = result || (fabs(x[i]) > MYEPSILON);
    515 
    516   return result;
     553        Vector x1;
     554        x1.CopyVector(y1);
     555        x1.Scale(x1.Projection(this));
     556        SubtractVector(&x1);
     557        for (int i=NDIM;i--;)
     558                result = result || (fabs(x[i]) > MYEPSILON);
     559
     560        return result;
    517561};
    518562
    … …  
    525569bool Vector::GetOneNormalVector(const Vector *GivenVector)
    526570{
    527   int Components[NDIM]; // contains indices of non-zero components
    528   int Last = 0;  // count the number of non-zero entries in vector
    529   int j;  // loop variables
    530   double norm;
    531 
    532   cout << Verbose(4);
    533   GivenVector->Output((ofstream *)&cout);
    534   cout << endl;
    535   for (j=NDIM;j--;)
    536     Components[j] = -1;
    537   // find two components != 0
    538   for (j=0;j<NDIM;j++)
    539     if (fabs(GivenVector->x[j]) > MYEPSILON)
    540       Components[Last++] = j;
    541   cout << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl;
    542 
    543   switch(Last) {
    544     case 3:  // threecomponent system
    545     case 2:  // two component system
    546       norm = sqrt(1./(GivenVector->x[Components[1]]*GivenVector->x[Components[1]]) + 1./(GivenVector->x[Components[0]]*GivenVector->x[Components[0]]));
    547       x[Components[2]] = 0.;
    548       // in skp both remaining parts shall become zero but with opposite sign and third is zero
    549       x[Components[1]] = -1./GivenVector->x[Components[1]] / norm;
    550       x[Components[0]] = 1./GivenVector->x[Components[0]] / norm;
    551       return true;
    552       break;
    553     case 1: // one component system
    554       // set sole non-zero component to 0, and one of the other zero component pendants to 1
    555       x[(Components[0]+2)%NDIM] = 0.;
    556       x[(Components[0]+1)%NDIM] = 1.;
    557       x[Components[0]] = 0.;
    558       return true;
    559       break;
    560     default:
    561       return false;
    562   }
     571        int Components[NDIM]; // contains indices of non-zero components
     572        int Last = 0;    // count the number of non-zero entries in vector
     573        int j;  // loop variables
     574        double norm;
     575
     576        cout << Verbose(4);
     577        GivenVector->Output((ofstream *)&cout);
     578        cout << endl;
     579        for (j=NDIM;j--;)
     580                Components[j] = -1;
     581        // find two components != 0
     582        for (j=0;j<NDIM;j++)
     583                if (fabs(GivenVector->x[j]) > MYEPSILON)
     584                        Components[Last++] = j;
     585        cout << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl;
     586
     587        switch(Last) {
     588                case 3: // threecomponent system
     589                case 2: // two component system
     590                        norm = sqrt(1./(GivenVector->x[Components[1]]*GivenVector->x[Components[1]]) + 1./(GivenVector->x[Components[0]]*GivenVector->x[Components[0]]));
     591                        x[Components[2]] = 0.;
     592                        // in skp both remaining parts shall become zero but with opposite sign and third is zero
     593                        x[Components[1]] = -1./GivenVector->x[Components[1]] / norm;
     594                        x[Components[0]] = 1./GivenVector->x[Components[0]] / norm;
     595                        return true;
     596                        break;
     597                case 1: // one component system
     598                        // set sole non-zero component to 0, and one of the other zero component pendants to 1
     599                        x[(Components[0]+2)%NDIM] = 0.;
     600                        x[(Components[0]+1)%NDIM] = 1.;
     601                        x[Components[0]] = 0.;
     602                        return true;
     603                        break;
     604                default:
     605                        return false;
     606        }
    563607};
    564608
    … …  
    571615double Vector::CutsPlaneAt(Vector *A, Vector *B, Vector *C)
    572616{
    573 //  cout << Verbose(3) << "For comparison: ";
    574 //  cout << "A " << A->Projection(this) << "\t";
    575 //  cout << "B " << B->Projection(this) << "\t";
    576 //  cout << "C " << C->Projection(this) << "\t";
    577 //  cout << endl;
    578   return A->Projection(this);
     617//      cout << Verbose(3) << "For comparison: ";
     618//      cout << "A " << A->Projection(this) << "\t";
     619//      cout << "B " << B->Projection(this) << "\t";
     620//      cout << "C " << C->Projection(this) << "\t";
     621//      cout << endl;
     622        return A->Projection(this);
    579623};
    580624
    … …  
    594638        }
    595639
    596   int np = 3;
     640        int np = 3;
    597641        struct LSQ_params par;
    598642
    599   const gsl_multimin_fminimizer_type *T =
    600     gsl_multimin_fminimizer_nmsimplex;
    601   gsl_multimin_fminimizer *s = NULL;
    602   gsl_vector *ss, *y;
    603   gsl_multimin_function minex_func;
    604 
    605   size_t iter = 0, i;
    606   int status;
    607   double size;
    608 
    609   /* Initial vertex size vector */
    610   ss = gsl_vector_alloc (np);
    611   y = gsl_vector_alloc (np);
    612 
    613   /* Set all step sizes to 1 */
    614   gsl_vector_set_all (ss, 1.0);
    615 
    616   /* Starting point */
    617   par.vectors = vectors;
     643        const gsl_multimin_fminimizer_type *T =
     644                gsl_multimin_fminimizer_nmsimplex;
     645        gsl_multimin_fminimizer *s = NULL;
     646        gsl_vector *ss, *y;
     647        gsl_multimin_function minex_func;
     648
     649        size_t iter = 0, i;
     650        int status;
     651        double size;
     652
     653        /* Initial vertex size vector */
     654        ss = gsl_vector_alloc (np);
     655        y = gsl_vector_alloc (np);
     656
     657        /* Set all step sizes to 1 */
     658        gsl_vector_set_all (ss, 1.0);
     659
     660        /* Starting point */
     661        par.vectors = vectors;
    618662         par.num = num;
    619663
    … …  
    621665                gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.);
    622666
    623   /* Initialize method and iterate */
    624   minex_func.f = &LSQ;
    625   minex_func.n = np;
    626   minex_func.params = (void *)&par;
    627 
    628   s = gsl_multimin_fminimizer_alloc (T, np);
    629   gsl_multimin_fminimizer_set (s, &minex_func, y, ss);
    630 
    631   do
    632     {
    633       iter++;
    634       status = gsl_multimin_fminimizer_iterate(s);
    635 
    636       if (status)
    637         break;
    638 
    639       size = gsl_multimin_fminimizer_size (s);
    640       status = gsl_multimin_test_size (size, 1e-2);
    641 
    642       if (status == GSL_SUCCESS)
    643         {
    644           printf ("converged to minimum at\n");
    645         }
    646 
    647       printf ("%5d ", (int)iter);
    648       for (i = 0; i < (size_t)np; i++)
    649         {
    650           printf ("%10.3e ", gsl_vector_get (s->x, i));
    651         }
    652       printf ("f() = %7.3f size = %.3f\n", s->fval, size);
    653     }
    654   while (status == GSL_CONTINUE && iter < 100);
    655 
    656   for (i=(size_t)np;i--;)
    657     this->x[i] = gsl_vector_get(s->x, i);
    658   gsl_vector_free(y);
    659   gsl_vector_free(ss);
    660   gsl_multimin_fminimizer_free (s);
     667        /* Initialize method and iterate */
     668        minex_func.f = &LSQ;
     669        minex_func.n = np;
     670        minex_func.params = (void *)&par;
     671
     672        s = gsl_multimin_fminimizer_alloc (T, np);
     673        gsl_multimin_fminimizer_set (s, &minex_func, y, ss);
     674
     675        do
     676                {
     677                        iter++;
     678                        status = gsl_multimin_fminimizer_iterate(s);
     679
     680                        if (status)
     681                                break;
     682
     683                        size = gsl_multimin_fminimizer_size (s);
     684                        status = gsl_multimin_test_size (size, 1e-2);
     685
     686                        if (status == GSL_SUCCESS)
     687                                {
     688                                        printf ("converged to minimum at\n");
     689                                }
     690
     691                        printf ("%5d ", (int)iter);
     692                        for (i = 0; i < (size_t)np; i++)
     693                                {
     694                                        printf ("%10.3e ", gsl_vector_get (s->x, i));
     695                                }
     696                        printf ("f() = %7.3f size = %.3f\n", s->fval, size);
     697                }
     698        while (status == GSL_CONTINUE && iter < 100);
     699
     700        for (i=(size_t)np;i--;)
     701                this->x[i] = gsl_vector_get(s->x, i);
     702        gsl_vector_free(y);
     703        gsl_vector_free(ss);
     704        gsl_multimin_fminimizer_free (s);
    661705
    662706        return true;
    … …  
    668712void Vector::AddVector(const Vector *y)
    669713{
    670   for (int i=NDIM;i--;)
    671     this->x[i] += y->x[i];
     714        for (int i=NDIM;i--;)
     715                this->x[i] += y->x[i];
    672716}
    673717
    … …  
    677721void Vector::SubtractVector(const Vector *y)
    678722{
    679   for (int i=NDIM;i--;)
    680     this->x[i] -= y->x[i];
     723        for (int i=NDIM;i--;)
     724                this->x[i] -= y->x[i];
    681725}
    682726
    … …  
    686730void Vector::CopyVector(const Vector *y)
    687731{
    688   for (int i=NDIM;i--;)
    689     this->x[i] = y->x[i];
     732        for (int i=NDIM;i--;)
     733                this->x[i] = y->x[i];
    690734}
    691735
    … …  
    697741void Vector::AskPosition(double *cell_size, bool check)
    698742{
    699   char coords[3] = {'x','y','z'};
    700   int j = -1;
    701   for (int i=0;i<3;i++) {
    702     j += i+1;
    703     do {
    704       cout << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: ";
    705       cin >> x[i];
    706     } while (((x[i] < 0) || (x[i] >= cell_size[j])) && (check));
    707   }
     743        char coords[3] = {'x','y','z'};
     744        int j = -1;
     745        for (int i=0;i<3;i++) {
     746                j += i+1;
     747                do {
     748                        cout << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: ";
     749                        cin >> x[i];
     750                } while (((x[i] < 0) || (x[i] >= cell_size[j])) && (check));
     751        }
    708752};
    709753
    … …  
    725769bool Vector::SolveSystem(Vector *x1, Vector *x2, Vector *y, double alpha, double beta, double c)
    726770{
    727   double D1,D2,D3,E1,E2,F1,F2,F3,p,q=0., A, B1, B2, C;
    728   double ang; // angle on testing
    729   double sign[3];
    730   int i,j,k;
    731   A = cos(alpha) * x1->Norm() * c;
    732   B1 = cos(beta + M_PI/2.) * y->Norm() * c;
    733   B2 = cos(beta) * x2->Norm() * c;
    734   C = c * c;
    735   cout << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl;
    736   int flag = 0;
    737   if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping
    738     if (fabs(x1->x[1]) > MYEPSILON) {
    739       flag = 1;
    740     } else if (fabs(x1->x[2]) > MYEPSILON) {
    741       flag = 2;
    742     } else {
    743       return false;
    744     }
    745   }
    746   switch (flag) {
    747     default:
    748     case 0:
    749       break;
    750     case 2:
    751       flip(&x1->x[0],&x1->x[1]);
    752       flip(&x2->x[0],&x2->x[1]);
    753       flip(&y->x[0],&y->x[1]);
    754       //flip(&x[0],&x[1]);
    755       flip(&x1->x[1],&x1->x[2]);
    756       flip(&x2->x[1],&x2->x[2]);
    757       flip(&y->x[1],&y->x[2]);
    758       //flip(&x[1],&x[2]);
    759     case 1:
    760       flip(&x1->x[0],&x1->x[1]);
    761       flip(&x2->x[0],&x2->x[1]);
    762       flip(&y->x[0],&y->x[1]);
    763       //flip(&x[0],&x[1]);
    764       flip(&x1->x[1],&x1->x[2]);
    765       flip(&x2->x[1],&x2->x[2]);
    766       flip(&y->x[1],&y->x[2]);
    767       //flip(&x[1],&x[2]);
    768       break;
    769   }
    770   // now comes the case system
    771   D1 = -y->x[0]/x1->x[0]*x1->x[1]+y->x[1];
    772   D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2];
    773   D3 = y->x[0]/x1->x[0]*A-B1;
    774   cout << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n";
    775   if (fabs(D1) < MYEPSILON) {
    776     cout << Verbose(2) << "D1 == 0!\n";
    777     if (fabs(D2) > MYEPSILON) {
    778       cout << Verbose(3) << "D2 != 0!\n";
    779       x[2] = -D3/D2;
    780       E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2;
    781       E2 = -x1->x[1]/x1->x[0];
    782       cout << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n";
    783       F1 = E1*E1 + 1.;
    784       F2 = -E1*E2;
    785       F3 = E1*E1 + D3*D3/(D2*D2) - C;
    786       cout << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
    787       if (fabs(F1) < MYEPSILON) {
    788         cout << Verbose(4) << "F1 == 0!\n";
    789         cout << Verbose(4) << "Gleichungssystem linear\n";
    790         x[1] = F3/(2.*F2);
    791       } else {
    792         p = F2/F1;
    793         q = p*p - F3/F1;
    794         cout << Verbose(4) << "p " << p << "\tq " << q << endl;
    795         if (q < 0) {
    796           cout << Verbose(4) << "q < 0" << endl;
    797           return false;
    798         }
    799         x[1] = p + sqrt(q);
    800       }
    801       x[0] =  A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
    802     } else {
    803       cout << Verbose(2) << "Gleichungssystem unterbestimmt\n";
    804       return false;
    805     }
    806   } else {
    807     E1 = A/x1->x[0]+x1->x[1]/x1->x[0]*D3/D1;
    808     E2 = x1->x[1]/x1->x[0]*D2/D1 - x1->x[2];
    809     cout << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n";
    810     F1 = E2*E2 + D2*D2/(D1*D1) + 1.;
    811     F2 = -(E1*E2 + D2*D3/(D1*D1));
    812     F3 = E1*E1 + D3*D3/(D1*D1) - C;
    813     cout << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
    814     if (fabs(F1) < MYEPSILON) {
    815       cout << Verbose(3) << "F1 == 0!\n";
    816       cout << Verbose(3) << "Gleichungssystem linear\n";
    817       x[2] = F3/(2.*F2);
    818     } else {
    819       p = F2/F1;
    820       q = p*p - F3/F1;
    821       cout << Verbose(3) << "p " << p << "\tq " << q << endl;
    822       if (q < 0) {
    823         cout << Verbose(3) << "q < 0" << endl;
    824         return false;
    825       }
    826       x[2] = p + sqrt(q);
    827     }
    828     x[1] = (-D2 * x[2] - D3)/D1;
    829     x[0] = A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
    830   }
    831   switch (flag) { // back-flipping
    832     default:
    833     case 0:
    834       break;
    835     case 2:
    836       flip(&x1->x[0],&x1->x[1]);
    837       flip(&x2->x[0],&x2->x[1]);
    838       flip(&y->x[0],&y->x[1]);
    839       flip(&x[0],&x[1]);
    840       flip(&x1->x[1],&x1->x[2]);
    841       flip(&x2->x[1],&x2->x[2]);
    842       flip(&y->x[1],&y->x[2]);
    843       flip(&x[1],&x[2]);
    844     case 1:
    845       flip(&x1->x[0],&x1->x[1]);
    846       flip(&x2->x[0],&x2->x[1]);
    847       flip(&y->x[0],&y->x[1]);
    848       //flip(&x[0],&x[1]);
    849       flip(&x1->x[1],&x1->x[2]);
    850       flip(&x2->x[1],&x2->x[2]);
    851       flip(&y->x[1],&y->x[2]);
    852       flip(&x[1],&x[2]);
    853       break;
    854   }
    855   // one z component is only determined by its radius (without sign)
    856   // thus check eight possible sign flips and determine by checking angle with second vector
    857   for (i=0;i<8;i++) {
    858     // set sign vector accordingly
    859     for (j=2;j>=0;j--) {
    860       k = (i & pot(2,j)) << j;
    861       cout << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl;
    862       sign[j] = (k == 0) ? 1. : -1.;
    863     }
    864     cout << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n";
    865     // apply sign matrix
    866     for (j=NDIM;j--;)
    867       x[j] *= sign[j];
    868     // calculate angle and check
    869     ang = x2->Angle (this);
    870     cout << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t";
    871     if (fabs(ang - cos(beta)) < MYEPSILON) {
    872       break;
    873     }
    874     // unapply sign matrix (is its own inverse)
    875     for (j=NDIM;j--;)
    876       x[j] *= sign[j];
    877   }
    878   return true;
    879 };
     771        double D1,D2,D3,E1,E2,F1,F2,F3,p,q=0., A, B1, B2, C;
     772        double ang; // angle on testing
     773        double sign[3];
     774        int i,j,k;
     775        A = cos(alpha) * x1->Norm() * c;
     776        B1 = cos(beta + M_PI/2.) * y->Norm() * c;
     777        B2 = cos(beta) * x2->Norm() * c;
     778        C = c * c;
     779        cout << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl;
     780        int flag = 0;
     781        if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping
     782                if (fabs(x1->x[1]) > MYEPSILON) {
     783                        flag = 1;
     784                } else if (fabs(x1->x[2]) > MYEPSILON) {
     785                        flag = 2;
     786                } else {
     787                        return false;
     788                }
     789        }
     790        switch (flag) {
     791                default:
     792                case 0:
     793                        break;
     794                case 2:
     795                        flip(&x1->x[0],&x1->x[1]);
     796                        flip(&x2->x[0],&x2->x[1]);
     797                        flip(&y->x[0],&y->x[1]);
     798                        //flip(&x[0],&x[1]);
     799                        flip(&x1->x[1],&x1->x[2]);
     800                        flip(&x2->x[1],&x2->x[2]);
     801                        flip(&y->x[1],&y->x[2]);
     802                        //flip(&x[1],&x[2]);
     803                case 1:
     804                        flip(&x1->x[0],&x1->x[1]);
     805                        flip(&x2->x[0],&x2->x[1]);
     806                        flip(&y->x[0],&y->x[1]);
     807                        //flip(&x[0],&x[1]);
     808                        flip(&x1->x[1],&x1->x[2]);
     809                        flip(&x2->x[1],&x2->x[2]);
     810                        flip(&y->x[1],&y->x[2]);
     811                        //flip(&x[1],&x[2]);
     812                        break;
     813        }
     814        // now comes the case system
     815        D1 = -y->x[0]/x1->x[0]*x1->x[1]+y->x[1];
     816        D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2];
     817        D3 = y->x[0]/x1->x[0]*A-B1;
     818        cout << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n";
     819        if (fabs(D1) < MYEPSILON) {
     820                cout << Verbose(2) << "D1 == 0!\n";
     821                if (fabs(D2) > MYEPSILON) {
     822                        cout << Verbose(3) << "D2 != 0!\n";
     823                        x[2] = -D3/D2;
     824                        E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2;
     825                        E2 = -x1->x[1]/x1->x[0];
     826                        cout << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n";
     827                        F1 = E1*E1 + 1.;
     828                        F2 = -E1*E2;
     829                        F3 = E1*E1 + D3*D3/(D2*D2) - C;
     830                        cout << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
     831                        if (fabs(F1) < MYEPSILON) {
     832                                cout << Verbose(4) << "F1 == 0!\n";
     833                                cout << Verbose(4) << "Gleichungssystem linear\n";
     834                                x[1] = F3/(2.*F2);
     835                        } else {
     836                                p = F2/F1;
     837                                q = p*p - F3/F1;
     838                                cout << Verbose(4) << "p " << p << "\tq " << q << endl;
     839                                if (q < 0) {
     840                                        cout << Verbose(4) << "q < 0" << endl;
     841                                        return false;
     842                                }
     843                                x[1] = p + sqrt(q);
     844                        }
     845                        x[0] =  A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
     846                } else {
     847                        cout << Verbose(2) << "Gleichungssystem unterbestimmt\n";
     848                        return false;
     849                }
     850        } else {
     851                E1 = A/x1->x[0]+x1->x[1]/x1->x[0]*D3/D1;
     852                E2 = x1->x[1]/x1->x[0]*D2/D1 - x1->x[2];
     853                cout << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n";
     854                F1 = E2*E2 + D2*D2/(D1*D1) + 1.;
     855                F2 = -(E1*E2 + D2*D3/(D1*D1));
     856                F3 = E1*E1 + D3*D3/(D1*D1) - C;
     857                cout << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
     858                if (fabs(F1) < MYEPSILON) {
     859                        cout << Verbose(3) << "F1 == 0!\n";
     860                        cout << Verbose(3) << "Gleichungssystem linear\n";
     861                        x[2] = F3/(2.*F2);
     862                } else {
     863                        p = F2/F1;
     864                        q = p*p - F3/F1;
     865                        cout << Verbose(3) << "p " << p << "\tq " << q << endl;
     866                        if (q < 0) {
     867                                cout << Verbose(3) << "q < 0" << endl;
     868                                return false;
     869                        }
     870                        x[2] = p + sqrt(q);
     871                }
     872                x[1] = (-D2 * x[2] - D3)/D1;
     873                x[0] = A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
     874        }
     875        switch (flag) { // back-flipping
     876                default:
     877                case 0:
     878                        break;
     879                case 2:
     880                        flip(&x1->x[0],&x1->x[1]);
     881                        flip(&x2->x[0],&x2->x[1]);
     882                        flip(&y->x[0],&y->x[1]);
     883                        flip(&x[0],&x[1]);
     884                        flip(&x1->x[1],&x1->x[2]);
     885                        flip(&x2->x[1],&x2->x[2]);
     886                        flip(&y->x[1],&y->x[2]);
     887                        flip(&x[1],&x[2]);
     888                case 1:
     889                        flip(&x1->x[0],&x1->x[1]);
     890                        flip(&x2->x[0],&x2->x[1]);
     891                        flip(&y->x[0],&y->x[1]);
     892                        //flip(&x[0],&x[1]);
     893                        flip(&x1->x[1],&x1->x[2]);
     894                        flip(&x2->x[1],&x2->x[2]);
     895                        flip(&y->x[1],&y->x[2]);
     896                        flip(&x[1],&x[2]);
     897                        break;
     898        }
     899        // one z component is only determined by its radius (without sign)
     900        // thus check eight possible sign flips and determine by checking angle with second vector
     901        for (i=0;i<8;i++) {
     902                // set sign vector accordingly
     903                for (j=2;j>=0;j--) {
     904                        k = (i & pot(2,j)) << j;
     905                        cout << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl;
     906                        sign[j] = (k == 0) ? 1. : -1.;
     907                }
     908                cout << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n";
     909                // apply sign matrix
     910                for (j=NDIM;j--;)
     911                        x[j] *= sign[j];
     912                // calculate angle and check
     913                ang = x2->Angle (this);
     914                cout << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t";
     915                if (fabs(ang - cos(beta)) < MYEPSILON) {
     916                        break;
     917                }
     918                // unapply sign matrix (is its own inverse)
     919                for (j=NDIM;j--;)
     920                        x[j] *= sign[j];
     921        }
     922        return true;
     923};
Note: See TracChangeset for help on using the changeset viewer.