Ignore:
Timestamp:
Jul 23, 2009, 1:45:24 PM (16 years ago)
Author:
Frederik Heber <heber@…>
Children:
47548d
Parents:
560995 (diff), 1b2aa1 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
git-author:
Frederik Heber <heber@…> (07/23/09 12:34:47)
git-committer:
Frederik Heber <heber@…> (07/23/09 13:45:24)
Message:

Merge branch 'MultipleMolecules'

Conflicts:

molecuilder/src/analyzer.cpp
molecuilder/src/atom.cpp
molecuilder/src/boundary.cpp
molecuilder/src/boundary.hpp
molecuilder/src/builder.cpp
molecuilder/src/config.cpp
molecuilder/src/datacreator.hpp
molecuilder/src/helpers.cpp
molecuilder/src/joiner.cpp
molecuilder/src/moleculelist.cpp
molecuilder/src/molecules.cpp
molecuilder/src/molecules.hpp
molecuilder/src/parser.cpp
molecuilder/src/parser.hpp
molecuilder/src/vector.cpp
molecuilder/src/verbose.cpp

merges:

compilation fixes:

File:
1 edited

Legend:

Unmodified
Added
Removed

  • molecuilder/src/vector.cpp

    r560995 rf39735  
    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
    … …  
    221259double Vector::Angle(const Vector *y) const
    222260{
    223   return acos(this->ScalarProduct(y)/Norm()/y->Norm());
     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 /** Prints a 3dim vector to a stream.
    316  * \param ost output stream
    317  * \param v Vector to be printed
    318  * \return output stream
    319  */
     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
    320353ostream& operator<<(ostream& ost,Vector& m)
    321354{
    322   ost << "(";
    323   for (int i=0;i<NDIM;i++) {
    324     ost << m.x[i];
    325     if (i != 2)
    326       ost << ",";
    327   }
    328   ost << ")";
    329   return ost;
     355        ost << "(";
     356        for (int i=0;i<NDIM;i++) {
     357                ost << m.x[i];
     358                if (i != 2)
     359                        ost << ",";
     360        }
     361        ost << ")";
     362        return ost;
    330363};
    331364
    … …  
    335368void Vector::Scale(double **factor)
    336369{
    337   for (int i=NDIM;i--;)
    338     x[i] *= (*factor)[i];
     370        for (int i=NDIM;i--;)
     371                x[i] *= (*factor)[i];
    339372};
    340373
    341374void Vector::Scale(double *factor)
    342375{
    343   for (int i=NDIM;i--;)
    344     x[i] *= *factor;
     376        for (int i=NDIM;i--;)
     377                x[i] *= *factor;
    345378};
    346379
    347380void Vector::Scale(double factor)
    348381{
    349   for (int i=NDIM;i--;)
    350     x[i] *= factor;
     382        for (int i=NDIM;i--;)
     383                x[i] *= factor;
    351384};
    352385
    … …  
    356389void Vector::Translate(const Vector *trans)
    357390{
    358   for (int i=NDIM;i--;)
    359     x[i] += trans->x[i];
     391        for (int i=NDIM;i--;)
     392                x[i] += trans->x[i];
    360393};
    361394
    … …  
    365398void Vector::MatrixMultiplication(double *M)
    366399{
    367   Vector C;
    368   // do the matrix multiplication
    369   C.x[0] = M[0]*x[0]+M[3]*x[1]+M[6]*x[2];
    370   C.x[1] = M[1]*x[0]+M[4]*x[1]+M[7]*x[2];
    371   C.x[2] = M[2]*x[0]+M[5]*x[1]+M[8]*x[2];
    372   // transfer the result into this
    373   for (int i=NDIM;i--;)
    374     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];
    375408};
    376409
    … …  
    380413void Vector::InverseMatrixMultiplication(double *A)
    381414{
    382   Vector C;
    383   double B[NDIM*NDIM];
    384   double detA = RDET3(A);
    385   double detAReci;
    386 
    387   // calculate the inverse B
    388   if (fabs(detA) > MYEPSILON) {;  // RDET3(A) yields precisely zero if A irregular
    389     detAReci = 1./detA;
    390     B[0] =  detAReci*RDET2(A[4],A[5],A[7],A[8]);    // A_11
    391     B[1] = -detAReci*RDET2(A[1],A[2],A[7],A[8]);    // A_12
    392     B[2] =  detAReci*RDET2(A[1],A[2],A[4],A[5]);    // A_13
    393     B[3] = -detAReci*RDET2(A[3],A[5],A[6],A[8]);    // A_21
    394     B[4] =  detAReci*RDET2(A[0],A[2],A[6],A[8]);    // A_22
    395     B[5] = -detAReci*RDET2(A[0],A[2],A[3],A[5]);    // A_23
    396     B[6] =  detAReci*RDET2(A[3],A[4],A[6],A[7]);    // A_31
    397     B[7] = -detAReci*RDET2(A[0],A[1],A[6],A[7]);    // A_32
    398     B[8] =  detAReci*RDET2(A[0],A[1],A[3],A[4]);    // A_33
    399 
    400     // do the matrix multiplication
    401     C.x[0] = B[0]*x[0]+B[3]*x[1]+B[6]*x[2];
    402     C.x[1] = B[1]*x[0]+B[4]*x[1]+B[7]*x[2];
    403     C.x[2] = B[2]*x[0]+B[5]*x[1]+B[8]*x[2];
    404     // transfer the result into this
    405     for (int i=NDIM;i--;)
    406       x[i] = C.x[i];
    407   } else {
    408     cerr << "ERROR: inverse of matrix does not exists!" << endl;
    409   }
     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        }
    410443};
    411444
    … …  
    420453void Vector::LinearCombinationOfVectors(const Vector *x1, const Vector *x2, const Vector *x3, double *factors)
    421454{
    422   for(int i=NDIM;i--;)
    423     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];
    424457};
    425458
    … …  
    429462void Vector::Mirror(const Vector *n)
    430463{
    431   double projection;
    432   projection = ScalarProduct(n)/n->ScalarProduct(n);    // remove constancy from n (keep as logical one)
    433   // withdraw projected vector twice from original one
    434   cout << Verbose(1) << "Vector: ";
    435   Output((ofstream *)&cout);
    436   cout << "\t";
    437   for (int i=NDIM;i--;)
    438     x[i] -= 2.*projection*n->x[i];
    439   cout << "Projected vector: ";
    440   Output((ofstream *)&cout);
    441   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;
    442475};
    443476
    … …  
    451484bool Vector::MakeNormalVector(const Vector *y1, const Vector *y2, const Vector *y3)
    452485{
    453   Vector x1, x2;
    454 
    455   x1.CopyVector(y1);
    456   x1.SubtractVector(y2);
    457   x2.CopyVector(y3);
    458   x2.SubtractVector(y2);
    459   if ((x1.Norm()==0) || (x2.Norm()==0)) {
    460     cout << Verbose(4) << "Given vectors are linear dependent." << endl;
    461     return false;
    462   }
    463 //  cout << Verbose(4) << "relative, first plane coordinates:";
    464 //  x1.Output((ofstream *)&cout);
    465 //  cout << endl;
    466 //  cout << Verbose(4) << "second plane coordinates:";
    467 //  x2.Output((ofstream *)&cout);
    468 //  cout << endl;
    469 
    470   this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
    471   this->x[1] = (x1.x[2]*x2.x[0] - x1.x[0]*x2.x[2]);
    472   this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
    473   Normalize();
    474 
    475   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;
    476509};
    477510
    … …  
    487520bool Vector::MakeNormalVector(const Vector *y1, const Vector *y2)
    488521{
    489   Vector x1,x2;
    490   x1.CopyVector(y1);
    491   x2.CopyVector(y2);
    492   Zero();
    493   if ((x1.Norm()==0) || (x2.Norm()==0)) {
    494     cout << Verbose(4) << "Given vectors are linear dependent." << endl;
    495     return false;
    496   }
    497 //  cout << Verbose(4) << "relative, first plane coordinates:";
    498 //  x1.Output((ofstream *)&cout);
    499 //  cout << endl;
    500 //  cout << Verbose(4) << "second plane coordinates:";
    501 //  x2.Output((ofstream *)&cout);
    502 //  cout << endl;
    503 
    504   this->x[0] = (x1.x[1]*x2.x[2] - x1.x[2]*x2.x[1]);
    505   this->x[1] = (x1.x[2]*x2.x[0] - x1.x[0]*x2.x[2]);
    506   this->x[2] = (x1.x[0]*x2.x[1] - x1.x[1]*x2.x[0]);
    507   Normalize();
    508 
    509   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;
    510543};
    511544
    … …  
    517550bool Vector::MakeNormalVector(const Vector *y1)
    518551{
    519   bool result = false;
    520   Vector x1;
    521   x1.CopyVector(y1);
    522   x1.Scale(x1.Projection(this));
    523   SubtractVector(&x1);
    524   for (int i=NDIM;i--;)
    525     result = result || (fabs(x[i]) > MYEPSILON);
    526 
    527   return result;
     552        bool result = false;
     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;
    528561};
    529562
    … …  
    536569bool Vector::GetOneNormalVector(const Vector *GivenVector)
    537570{
    538   int Components[NDIM]; // contains indices of non-zero components
    539   int Last = 0;  // count the number of non-zero entries in vector
    540   int j;  // loop variables
    541   double norm;
    542 
    543   cout << Verbose(4);
    544   GivenVector->Output((ofstream *)&cout);
    545   cout << endl;
    546   for (j=NDIM;j--;)
    547     Components[j] = -1;
    548   // find two components != 0
    549   for (j=0;j<NDIM;j++)
    550     if (fabs(GivenVector->x[j]) > MYEPSILON)
    551       Components[Last++] = j;
    552   cout << Verbose(4) << Last << " Components != 0: (" << Components[0] << "," << Components[1] << "," << Components[2] << ")" << endl;
    553 
    554   switch(Last) {
    555     case 3:  // threecomponent system
    556     case 2:  // two component system
    557       norm = sqrt(1./(GivenVector->x[Components[1]]*GivenVector->x[Components[1]]) + 1./(GivenVector->x[Components[0]]*GivenVector->x[Components[0]]));
    558       x[Components[2]] = 0.;
    559       // in skp both remaining parts shall become zero but with opposite sign and third is zero
    560       x[Components[1]] = -1./GivenVector->x[Components[1]] / norm;
    561       x[Components[0]] = 1./GivenVector->x[Components[0]] / norm;
    562       return true;
    563       break;
    564     case 1: // one component system
    565       // set sole non-zero component to 0, and one of the other zero component pendants to 1
    566       x[(Components[0]+2)%NDIM] = 0.;
    567       x[(Components[0]+1)%NDIM] = 1.;
    568       x[Components[0]] = 0.;
    569       return true;
    570       break;
    571     default:
    572       return false;
    573   }
     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        }
    574607};
    575608
    … …  
    582615double Vector::CutsPlaneAt(Vector *A, Vector *B, Vector *C)
    583616{
    584 //  cout << Verbose(3) << "For comparison: ";
    585 //  cout << "A " << A->Projection(this) << "\t";
    586 //  cout << "B " << B->Projection(this) << "\t";
    587 //  cout << "C " << C->Projection(this) << "\t";
    588 //  cout << endl;
    589   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);
    590623};
    591624
    … …  
    597630bool Vector::LSQdistance(Vector **vectors, int num)
    598631{
    599   int j;
    600 
    601   for (j=0;j<num;j++) {
    602     cout << Verbose(1) << j << "th atom's vector: ";
    603     (vectors[j])->Output((ofstream *)&cout);
    604     cout << endl;
    605   }
    606 
    607   int np = 3;
    608   struct LSQ_params par;
    609 
    610   const gsl_multimin_fminimizer_type *T =
    611     gsl_multimin_fminimizer_nmsimplex;
    612   gsl_multimin_fminimizer *s = NULL;
    613   gsl_vector *ss, *y;
    614   gsl_multimin_function minex_func;
    615 
    616   size_t iter = 0, i;
    617   int status;
    618   double size;
    619 
    620   /* Initial vertex size vector */
    621   ss = gsl_vector_alloc (np);
    622   y = gsl_vector_alloc (np);
    623 
    624   /* Set all step sizes to 1 */
    625   gsl_vector_set_all (ss, 1.0);
    626 
    627   /* Starting point */
    628   par.vectors = vectors;
    629   par.num = num;
    630 
    631   for (i=NDIM;i--;)
    632     gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.);
    633 
    634   /* Initialize method and iterate */
    635   minex_func.f = &LSQ;
    636   minex_func.n = np;
    637   minex_func.params = (void *)&par;
    638 
    639   s = gsl_multimin_fminimizer_alloc (T, np);
    640   gsl_multimin_fminimizer_set (s, &minex_func, y, ss);
    641 
    642   do
    643     {
    644       iter++;
    645       status = gsl_multimin_fminimizer_iterate(s);
    646 
    647       if (status)
    648         break;
    649 
    650       size = gsl_multimin_fminimizer_size (s);
    651       status = gsl_multimin_test_size (size, 1e-2);
    652 
    653       if (status == GSL_SUCCESS)
    654         {
    655           printf ("converged to minimum at\n");
    656         }
    657 
    658       printf ("%5d ", (int)iter);
    659       for (i = 0; i < (size_t)np; i++)
    660         {
    661           printf ("%10.3e ", gsl_vector_get (s->x, i));
    662         }
    663       printf ("f() = %7.3f size = %.3f\n", s->fval, size);
    664     }
    665   while (status == GSL_CONTINUE && iter < 100);
    666 
    667   for (i=(size_t)np;i--;)
    668     this->x[i] = gsl_vector_get(s->x, i);
    669   gsl_vector_free(y);
    670   gsl_vector_free(ss);
    671   gsl_multimin_fminimizer_free (s);
    672 
    673   return true;
     632        int j;
     633
     634        for (j=0;j<num;j++) {
     635                cout << Verbose(1) << j << "th atom's vector: ";
     636                (vectors[j])->Output((ofstream *)&cout);
     637                cout << endl;
     638        }
     639
     640        int np = 3;
     641        struct LSQ_params par;
     642
     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;
     662        par.num = num;
     663
     664        for (i=NDIM;i--;)
     665                gsl_vector_set(y, i, (vectors[0]->x[i] - vectors[1]->x[i])/2.);
     666
     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);
     705
     706        return true;
    674707};
    675708
    … …  
    679712void Vector::AddVector(const Vector *y)
    680713{
    681   for (int i=NDIM;i--;)
    682     this->x[i] += y->x[i];
     714        for (int i=NDIM;i--;)
     715                this->x[i] += y->x[i];
    683716}
    684717
    … …  
    688721void Vector::SubtractVector(const Vector *y)
    689722{
    690   for (int i=NDIM;i--;)
    691     this->x[i] -= y->x[i];
     723        for (int i=NDIM;i--;)
     724                this->x[i] -= y->x[i];
    692725}
    693726
    … …  
    697730void Vector::CopyVector(const Vector *y)
    698731{
    699   for (int i=NDIM;i--;)
    700     this->x[i] = y->x[i];
     732        for (int i=NDIM;i--;)
     733                this->x[i] = y->x[i];
    701734}
    702735
    … …  
    708741void Vector::AskPosition(double *cell_size, bool check)
    709742{
    710   char coords[3] = {'x','y','z'};
    711   int j = -1;
    712   for (int i=0;i<3;i++) {
    713     j += i+1;
    714     do {
    715       cout << Verbose(0) << coords[i] << "[0.." << cell_size[j] << "]: ";
    716       cin >> x[i];
    717     } while (((x[i] < 0) || (x[i] >= cell_size[j])) && (check));
    718   }
     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        }
    719752};
    720753
    … …  
    736769bool Vector::SolveSystem(Vector *x1, Vector *x2, Vector *y, double alpha, double beta, double c)
    737770{
    738   double D1,D2,D3,E1,E2,F1,F2,F3,p,q=0., A, B1, B2, C;
    739   double ang; // angle on testing
    740   double sign[3];
    741   int i,j,k;
    742   A = cos(alpha) * x1->Norm() * c;
    743   B1 = cos(beta + M_PI/2.) * y->Norm() * c;
    744   B2 = cos(beta) * x2->Norm() * c;
    745   C = c * c;
    746   cout << Verbose(2) << "A " << A << "\tB " << B1 << "\tC " << C << endl;
    747   int flag = 0;
    748   if (fabs(x1->x[0]) < MYEPSILON) { // check for zero components for the later flipping and back-flipping
    749     if (fabs(x1->x[1]) > MYEPSILON) {
    750       flag = 1;
    751     } else if (fabs(x1->x[2]) > MYEPSILON) {
    752       flag = 2;
    753     } else {
    754       return false;
    755     }
    756   }
    757   switch (flag) {
    758     default:
    759     case 0:
    760       break;
    761     case 2:
    762       flip(&x1->x[0],&x1->x[1]);
    763       flip(&x2->x[0],&x2->x[1]);
    764       flip(&y->x[0],&y->x[1]);
    765       //flip(&x[0],&x[1]);
    766       flip(&x1->x[1],&x1->x[2]);
    767       flip(&x2->x[1],&x2->x[2]);
    768       flip(&y->x[1],&y->x[2]);
    769       //flip(&x[1],&x[2]);
    770     case 1:
    771       flip(&x1->x[0],&x1->x[1]);
    772       flip(&x2->x[0],&x2->x[1]);
    773       flip(&y->x[0],&y->x[1]);
    774       //flip(&x[0],&x[1]);
    775       flip(&x1->x[1],&x1->x[2]);
    776       flip(&x2->x[1],&x2->x[2]);
    777       flip(&y->x[1],&y->x[2]);
    778       //flip(&x[1],&x[2]);
    779       break;
    780   }
    781   // now comes the case system
    782   D1 = -y->x[0]/x1->x[0]*x1->x[1]+y->x[1];
    783   D2 = -y->x[0]/x1->x[0]*x1->x[2]+y->x[2];
    784   D3 = y->x[0]/x1->x[0]*A-B1;
    785   cout << Verbose(2) << "D1 " << D1 << "\tD2 " << D2 << "\tD3 " << D3 << "\n";
    786   if (fabs(D1) < MYEPSILON) {
    787     cout << Verbose(2) << "D1 == 0!\n";
    788     if (fabs(D2) > MYEPSILON) {
    789       cout << Verbose(3) << "D2 != 0!\n";
    790       x[2] = -D3/D2;
    791       E1 = A/x1->x[0] + x1->x[2]/x1->x[0]*D3/D2;
    792       E2 = -x1->x[1]/x1->x[0];
    793       cout << Verbose(3) << "E1 " << E1 << "\tE2 " << E2 << "\n";
    794       F1 = E1*E1 + 1.;
    795       F2 = -E1*E2;
    796       F3 = E1*E1 + D3*D3/(D2*D2) - C;
    797       cout << Verbose(3) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
    798       if (fabs(F1) < MYEPSILON) {
    799         cout << Verbose(4) << "F1 == 0!\n";
    800         cout << Verbose(4) << "Gleichungssystem linear\n";
    801         x[1] = F3/(2.*F2);
    802       } else {
    803         p = F2/F1;
    804         q = p*p - F3/F1;
    805         cout << Verbose(4) << "p " << p << "\tq " << q << endl;
    806         if (q < 0) {
    807           cout << Verbose(4) << "q < 0" << endl;
    808           return false;
    809         }
    810         x[1] = p + sqrt(q);
    811       }
    812       x[0] =  A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
    813     } else {
    814       cout << Verbose(2) << "Gleichungssystem unterbestimmt\n";
    815       return false;
    816     }
    817   } else {
    818     E1 = A/x1->x[0]+x1->x[1]/x1->x[0]*D3/D1;
    819     E2 = x1->x[1]/x1->x[0]*D2/D1 - x1->x[2];
    820     cout << Verbose(2) << "E1 " << E1 << "\tE2 " << E2 << "\n";
    821     F1 = E2*E2 + D2*D2/(D1*D1) + 1.;
    822     F2 = -(E1*E2 + D2*D3/(D1*D1));
    823     F3 = E1*E1 + D3*D3/(D1*D1) - C;
    824     cout << Verbose(2) << "F1 " << F1 << "\tF2 " << F2 << "\tF3 " << F3 << "\n";
    825     if (fabs(F1) < MYEPSILON) {
    826       cout << Verbose(3) << "F1 == 0!\n";
    827       cout << Verbose(3) << "Gleichungssystem linear\n";
    828       x[2] = F3/(2.*F2);
    829     } else {
    830       p = F2/F1;
    831       q = p*p - F3/F1;
    832       cout << Verbose(3) << "p " << p << "\tq " << q << endl;
    833       if (q < 0) {
    834         cout << Verbose(3) << "q < 0" << endl;
    835         return false;
    836       }
    837       x[2] = p + sqrt(q);
    838     }
    839     x[1] = (-D2 * x[2] - D3)/D1;
    840     x[0] = A/x1->x[0] - x1->x[1]/x1->x[0]*x[1] + x1->x[2]/x1->x[0]*x[2];
    841   }
    842   switch (flag) { // back-flipping
    843     default:
    844     case 0:
    845       break;
    846     case 2:
    847       flip(&x1->x[0],&x1->x[1]);
    848       flip(&x2->x[0],&x2->x[1]);
    849       flip(&y->x[0],&y->x[1]);
    850       flip(&x[0],&x[1]);
    851       flip(&x1->x[1],&x1->x[2]);
    852       flip(&x2->x[1],&x2->x[2]);
    853       flip(&y->x[1],&y->x[2]);
    854       flip(&x[1],&x[2]);
    855     case 1:
    856       flip(&x1->x[0],&x1->x[1]);
    857       flip(&x2->x[0],&x2->x[1]);
    858       flip(&y->x[0],&y->x[1]);
    859       //flip(&x[0],&x[1]);
    860       flip(&x1->x[1],&x1->x[2]);
    861       flip(&x2->x[1],&x2->x[2]);
    862       flip(&y->x[1],&y->x[2]);
    863       flip(&x[1],&x[2]);
    864       break;
    865   }
    866   // one z component is only determined by its radius (without sign)
    867   // thus check eight possible sign flips and determine by checking angle with second vector
    868   for (i=0;i<8;i++) {
    869     // set sign vector accordingly
    870     for (j=2;j>=0;j--) {
    871       k = (i & pot(2,j)) << j;
    872       cout << Verbose(2) << "k " << k << "\tpot(2,j) " << pot(2,j) << endl;
    873       sign[j] = (k == 0) ? 1. : -1.;
    874     }
    875     cout << Verbose(2) << i << ": sign matrix is " << sign[0] << "\t" << sign[1] << "\t" << sign[2] << "\n";
    876     // apply sign matrix
    877     for (j=NDIM;j--;)
    878       x[j] *= sign[j];
    879     // calculate angle and check
    880     ang = x2->Angle (this);
    881     cout << Verbose(1) << i << "th angle " << ang << "\tbeta " << cos(beta) << " :\t";
    882     if (fabs(ang - cos(beta)) < MYEPSILON) {
    883       break;
    884     }
    885     // unapply sign matrix (is its own inverse)
    886     for (j=NDIM;j--;)
    887       x[j] *= sign[j];
    888   }
    889   return true;
    890 };
     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.