Add angle sets (disjoint ranges) - initial support

By adding a new angleset class, things get much easier to read. Also, the part of the previous commit that dealt with combining the individual ranges was kind of stupid.
Things missing:
o) functions to remove angle ranges from a set
o) deferred add functions, that allow adding multiple ranges without calling the O(n²) function angleset::combine() in between.

1 files changed, 173 insertions, 220 deletions

diff --git a/main.cpp b/main.cpp
index 51e9506..0f3f62e 100644
--- a/main.cpp
+++ b/main.cpp
@@ -68,7 +68,7 @@
 #include <cmath>
 #include <cstdio>
 #include <algorithm>
-#include "anglerange.h"
+#include "angleset.h"
 
 using namespace std;
 
@@ -87,7 +87,7 @@ int main(int argc, char* argv[])
     {
         //read in command line arguments
         double commargs[9], swapper;
-        unsigned int i,j;
+        unsigned int i;
         for(i=1;i<(unsigned)argc;i++)
         {
             //what is this stringstreams stuff everyone is hyped about?!?
@@ -139,16 +139,16 @@ int main(int argc, char* argv[])
         //solutions run from -maxn to maxn, and from -maxm to maxm.
         //Due to the ambiguity of asin, two solutions exist for each value of n and m
         //This means, that (2*maxn+1)*2 solutions exist, the same for m.
-        anglerange pxranges[4*maxn+2];
-        anglerange qxranges[4*maxm+2];
+        angleset pxranges;
+        pxranges.reserve(4*maxn+2); //reserve memory, so adding stuff is faster...
+        angleset qxranges;
+        qxranges.reserve(4*maxm+2);
 
         //up to now it was more or less dull C code. Now comes the first real difference:
         //As asin changes sign together with its argument, one has to treat positive and negative n differently
         //first the special case n=0:
-        pxranges[0].setlower(commargs[ALPHA]);
-        pxranges[0].setupper(commargs[ALPHA]);
-        pxranges[1].setlower(commargs[ALPHA] - M_PI);
-        pxranges[1].setupper(commargs[ALPHA] - M_PI);
+        pxranges.add(commargs[ALPHA],commargs[ALPHA]);
+        pxranges.add(commargs[ALPHA] - M_PI, commargs[ALPHA] - M_PI);
         //now the slightly more difficult case: n>0
         for(i=1;i<=maxn;i++){
             //while factoring out the asin doesn't improve performance much - it's only used twice, it improves readability, as the important thing in the formulas below
@@ -159,37 +159,57 @@ int main(int argc, char* argv[])
             //we need to consider that sin(alpha) can be negative. In that case the sign of the asin will change as well.
             if(asina1b1min>=0)
             {
-                pxranges[2*i].setlower(commargs[ALPHA] - asina1b1min);
-                pxranges[2*i].setupper(commargs[ALPHA] - asina1b1max);
-                pxranges[2*i+1].setlower(commargs[ALPHA] - M_PI + asina1b1max);
-                pxranges[2*i+1].setupper(commargs[ALPHA] - M_PI + asina1b1min);
+                pxranges.add(
+                            commargs[ALPHA] - asina1b1min,
+                            commargs[ALPHA] - asina1b1max
+                            );
+                pxranges.add(
+                            commargs[ALPHA] - M_PI + asina1b1max,
+                            commargs[ALPHA] - M_PI + asina1b1min
+                            );
                 //and the most difficult case: n<0
                 //here the arcsin is negative. as i is positive, I'll just change the sign in front of the arcsin.
-                pxranges[2*i+2*maxn].setlower(commargs[ALPHA] + asina1b1max);
-                pxranges[2*i+2*maxn].setupper(commargs[ALPHA] + asina1b1min);
-                pxranges[2*i+2*maxn+1].setlower(commargs[ALPHA] - M_PI - asina1b1min);
-                pxranges[2*i+2*maxn+1].setupper(commargs[ALPHA] - M_PI - asina1b1max);
+                pxranges.add(
+                            commargs[ALPHA] + asina1b1max,
+                            commargs[ALPHA] + asina1b1min
+                            );
+                pxranges.add(
+                            commargs[ALPHA] - M_PI - asina1b1min,
+                            commargs[ALPHA] - M_PI - asina1b1max
+                            );
             }
             else
             {
                 //just as above, but with upper and lower limits switched
-                pxranges[2*i].setlower(commargs[ALPHA] - asina1b1max);
-                pxranges[2*i].setupper(commargs[ALPHA] - asina1b1min);
-                pxranges[2*i+1].setlower(commargs[ALPHA] - M_PI + asina1b1min);
-                pxranges[2*i+1].setupper(commargs[ALPHA] - M_PI + asina1b1max);
+                pxranges.add(
+                            commargs[ALPHA] - asina1b1max,
+                            commargs[ALPHA] - asina1b1min
+                            );
+                pxranges.add(
+                            commargs[ALPHA] - M_PI + asina1b1min,
+                            commargs[ALPHA] - M_PI + asina1b1max
+                            );
                 //n<0
-                pxranges[2*i+2*maxn].setlower(commargs[ALPHA] + asina1b1min);
-                pxranges[2*i+2*maxn].setupper(commargs[ALPHA] + asina1b1max);
-                pxranges[2*i+2*maxn+1].setlower(commargs[ALPHA] - M_PI - asina1b1max);
-                pxranges[2*i+2*maxn+1].setupper(commargs[ALPHA] - M_PI - asina1b1min);
+                pxranges.add(
+                            commargs[ALPHA] + asina1b1min,
+                            commargs[ALPHA] + asina1b1max
+                            );
+                pxranges.add(
+                            commargs[ALPHA] - M_PI - asina1b1max,
+                            commargs[ALPHA] - M_PI - asina1b1min
+                            );
             }
         }
         //same nonsense for qxranges
         //first the easy part: m=0;
-        qxranges[0].setlower(commargs[ALPHA] - commargs[BETAMAX]);
-        qxranges[0].setupper(commargs[ALPHA] - commargs[BETAMIN]);
-        qxranges[1].setlower(commargs[ALPHA] - commargs[BETAMAX] - M_PI);
-        qxranges[1].setupper(commargs[ALPHA] - commargs[BETAMIN] - M_PI);
+        qxranges.add(
+                    commargs[ALPHA] - commargs[BETAMAX],
+                    commargs[ALPHA] - commargs[BETAMIN]
+                    );
+        qxranges.add(
+                    commargs[ALPHA] - commargs[BETAMAX] - M_PI,
+                    commargs[ALPHA] - commargs[BETAMIN] - M_PI
+                    );
         //now the slightly more difficult case: m>0;
         for(i=1;i<=maxm;i++){
             //also here, the asin values are factored out for improved readability.
@@ -198,60 +218,63 @@ int main(int argc, char* argv[])
             //same here: keep in mind that sin(alpha) can be negative:
             if(asina1b2min>=0)
             {
-                qxranges[2*i].setlower(commargs[ALPHA] - commargs[BETAMAX] - asina1b2min);
-                qxranges[2*i].setupper(commargs[ALPHA] - commargs[BETAMIN] - asina1b2max);
-                qxranges[2*i+1].setlower(commargs[ALPHA] - commargs[BETAMAX] - M_PI + asina1b2max);
-                qxranges[2*i+1].setupper(commargs[ALPHA] - commargs[BETAMIN] - M_PI + asina1b2min);
+                qxranges.add(
+                            commargs[ALPHA] - commargs[BETAMAX] - asina1b2min,
+                            commargs[ALPHA] - commargs[BETAMIN] - asina1b2max
+                            );
+                qxranges.add(
+                            commargs[ALPHA] - commargs[BETAMAX] - M_PI + asina1b2max,
+                            commargs[ALPHA] - commargs[BETAMIN] - M_PI + asina1b2min
+                            );
                 //and last, but not leasst, the most difficult, m<0 - here the arcsin is negative;
                 //as i is positive, I'll just change the sign in front of the arcsin.
-                qxranges[2*i+2*maxm].setlower(commargs[ALPHA] - commargs[BETAMAX] + asina1b2max);
-                qxranges[2*i+2*maxm].setupper(commargs[ALPHA] - commargs[BETAMIN] + asina1b2min);
-                qxranges[2*i+1+2*maxm].setlower(commargs[ALPHA] - commargs[BETAMAX] - M_PI - asina1b2min);
-                qxranges[2*i+1+2*maxm].setupper(commargs[ALPHA] - commargs[BETAMIN] - M_PI - asina1b2max);
+                qxranges.add(
+                            commargs[ALPHA] - commargs[BETAMAX] + asina1b2max,
+                            commargs[ALPHA] - commargs[BETAMIN] + asina1b2min
+                            );
+                qxranges.add(
+                            commargs[ALPHA] - commargs[BETAMAX] - M_PI - asina1b2min,
+                            commargs[ALPHA] - commargs[BETAMIN] - M_PI - asina1b2max
+                            );
             }
             else
             {
-                qxranges[2*i].setlower(commargs[ALPHA] - commargs[BETAMAX] - asina1b2max);
-                qxranges[2*i].setupper(commargs[ALPHA] - commargs[BETAMIN] - asina1b2min);
-                qxranges[2*i+1].setlower(commargs[ALPHA] - commargs[BETAMAX] - M_PI + asina1b2min);
-                qxranges[2*i+1].setupper(commargs[ALPHA] - commargs[BETAMIN] - M_PI + asina1b2max);
+                qxranges.add(
+                            commargs[ALPHA] - commargs[BETAMAX] - asina1b2max,
+                            commargs[ALPHA] - commargs[BETAMIN] - asina1b2min
+                            );
+                qxranges.add(
+                            commargs[ALPHA] - commargs[BETAMAX] - M_PI + asina1b2min,
+                            commargs[ALPHA] - commargs[BETAMIN] - M_PI + asina1b2max
+                            );
                 //m<0
-                qxranges[2*i+2*maxm].setlower(commargs[ALPHA] - commargs[BETAMAX] + asina1b2min);
-                qxranges[2*i+2*maxm].setupper(commargs[ALPHA] - commargs[BETAMIN] + asina1b2max);
-                qxranges[2*i+1+2*maxm].setlower(commargs[ALPHA] - commargs[BETAMAX] - M_PI - asina1b2max);
-                qxranges[2*i+1+2*maxm].setupper(commargs[ALPHA] - commargs[BETAMIN] - M_PI - asina1b2min);
-            }
-        }
-        //Now the real big change from the plain C version: Using the anglerange class to generate a list
-        //of overlaps instead of just outputting all of them.
-
-        std::vector<anglerange> xoverlaps; //keep this, as the overlap between this and yoverlaps gives coincident results
-        for(i=0;i<4*maxn+2;i++)
-        {
-            for(j=0;j<4*maxm+2;j++)
-            {
-                anglerange tmp = pxranges[i].overlap(qxranges[j]);
-                if(!(tmp.isempty()))
-                {
-                    xoverlaps.push_back(tmp);
-                }
+                qxranges.add(
+                            commargs[ALPHA] - commargs[BETAMAX] + asina1b2min,
+                            commargs[ALPHA] - commargs[BETAMIN] + asina1b2max
+                            );
+                qxranges.add(
+                            commargs[ALPHA] - commargs[BETAMAX] - M_PI - asina1b2max,
+                            commargs[ALPHA] - commargs[BETAMIN] - M_PI - asina1b2min
+                            );
             }
         }
+        //Calculate the overlap between these two:
+        angleset xoverlaps=pxranges.overlap(qxranges);
 
         //ok, same thing for qy, py:
         //I know that this is wasteful regarding RAM, but well, we're still talking about bytes, not about megabytes ;-)
         unsigned int maxo, maxp;
         maxo=fabs(commargs[B1MAX]/(commargs[A2]*sin(commargs[ALPHA])));
         maxp=fabs(commargs[B2MAX]/(commargs[A2]*sin(commargs[ALPHA])));
-        anglerange qyranges[4*maxo+2];
-        anglerange pyranges[4*maxp+2];
+        angleset qyranges;
+        qyranges.reserve(4*maxo+2);
+        angleset pyranges;
+        pyranges.reserve(4*maxp+2);
 
         //now let's start with qyranges. As previously we need to consider the "sign" of o,p, and sin(alpha)
         //first: o=0
-        qyranges[0].setlower(0.0);
-        qyranges[0].setupper(0.0);
-        qyranges[1].setlower(M_PI);
-        qyranges[1].setupper(M_PI);
+        qyranges.add(0.0,0.0);
+        qyranges.add(M_PI,M_PI);
         for(i=1;i<=maxo;i++)
         {
             //also here: factor out the asin for improved readability.
@@ -261,36 +284,56 @@ int main(int argc, char* argv[])
             if(asina2b1max>=0)
             {
                 //case: o>0
-                qyranges[2*i].setlower(asina2b1max);
-                qyranges[2*i].setupper(asina2b1min);
-                qyranges[2*i+1].setlower(M_PI - asina2b1min);
-                qyranges[2*i+1].setupper(M_PI - asina2b1max);
+                qyranges.add(
+                            asina2b1max,
+                            asina2b1min
+                            );
+                qyranges.add(
+                            M_PI - asina2b1min,
+                            M_PI - asina2b1max
+                            );
                 //case: o<0
-                qyranges[2*i+2*maxo].setlower( -asina2b1min);
-                qyranges[2*i+2*maxo].setupper( -asina2b1max);
-                qyranges[2*i+1+2*maxo].setlower(M_PI + asina2b1max);
-                qyranges[2*i+1+2*maxo].setupper(M_PI + asina2b1min);
+                qyranges.add(
+                            -asina2b1min,
+                            -asina2b1max
+                            );
+                qyranges.add(
+                            M_PI + asina2b1max,
+                            M_PI + asina2b1min
+                            );
             }
             else
             {
                 //case: o>0
-                qyranges[2*i].setlower(asina2b1min);
-                qyranges[2*i].setupper(asina2b1max);
-                qyranges[2*i+1].setlower(M_PI - asina2b1max);
-                qyranges[2*i+1].setupper(M_PI - asina2b1min);
+                qyranges.add(
+                            asina2b1min,
+                            asina2b1max
+                            );
+                qyranges.add(
+                            M_PI - asina2b1max,
+                            M_PI - asina2b1min
+                            );
                 //case: o<0
-                qyranges[2*i+2*maxo].setlower( -asina2b1max);
-                qyranges[2*i+2*maxo].setupper( -asina2b1min);
-                qyranges[2*i+1+2*maxo].setlower(M_PI + asina2b1min);
-                qyranges[2*i+1+2*maxo].setupper(M_PI + asina2b1max);
+                qyranges.add(
+                            -asina2b1max,
+                            -asina2b1min
+                            );
+                qyranges.add(
+                            M_PI + asina2b1min,
+                            M_PI + asina2b1max
+                            );
             }
         }
         //that was too easy. Probably it's buggy as hell...
         //now to py
-        pyranges[0].setlower( -commargs[BETAMAX]);
-        pyranges[0].setupper( -commargs[BETAMIN]);
-        pyranges[1].setlower(M_PI - commargs[BETAMAX]);
-        pyranges[1].setupper(M_PI - commargs[BETAMIN]);
+        pyranges.add(
+                    -commargs[BETAMAX],
+                    -commargs[BETAMIN]
+                    );
+        pyranges.add(
+                    M_PI - commargs[BETAMAX],
+                    M_PI - commargs[BETAMIN]
+                    );
         for(i=1;i<=maxp;i++)
         {
             //and again: readability
@@ -299,156 +342,66 @@ int main(int argc, char* argv[])
             if(asina2b2max>=0)
             {
                 //case: p>0
-                pyranges[2*i].setlower(asina2b2max - commargs[BETAMAX]);
-                pyranges[2*i].setupper(asina2b2min - commargs[BETAMIN]);
-                pyranges[2*i+1].setlower(M_PI - asina2b2min - commargs[BETAMAX]);
-                pyranges[2*i+1].setupper(M_PI - asina2b2max - commargs[BETAMIN]);
+                pyranges.add(
+                            asina2b2max - commargs[BETAMAX],
+                            asina2b2min - commargs[BETAMIN]
+                            );
+                pyranges.add(
+                            M_PI - asina2b2min - commargs[BETAMAX],
+                            M_PI - asina2b2max - commargs[BETAMIN]
+                            );
                 //case: p<0
-                pyranges[2*i+2*maxp].setlower( -asina2b2min - commargs[BETAMAX]);
-                pyranges[2*i+2*maxp].setupper( -asina2b2max - commargs[BETAMIN]);
-                pyranges[2*i+1+2*maxp].setlower(M_PI + asina2b2max - commargs[BETAMAX]);
-                pyranges[2*i+1+2*maxp].setupper(M_PI + asina2b2min - commargs[BETAMIN]);
+                pyranges.add(
+                            -asina2b2min - commargs[BETAMAX],
+                            -asina2b2max - commargs[BETAMIN]
+                            );
+                pyranges.add(
+                            M_PI + asina2b2max - commargs[BETAMAX],
+                            M_PI + asina2b2min - commargs[BETAMIN]
+                            );
             }
             else
             {
                 //ok, here the asin is of opposite sign!
                 //case p>0
-                pyranges[2*i].setlower(asina2b2min - commargs[BETAMAX]);
-                pyranges[2*i].setupper(asina2b2max - commargs[BETAMIN]);
-                pyranges[2*i+1].setlower(M_PI - asina2b2max - commargs[BETAMAX]);
-                pyranges[2*i+1].setupper(M_PI - asina2b2min - commargs[BETAMIN]);
+                pyranges.add(
+                            asina2b2min - commargs[BETAMAX],
+                            asina2b2max - commargs[BETAMIN]
+                            );
+                pyranges.add(
+                            M_PI - asina2b2max - commargs[BETAMAX],
+                            M_PI - asina2b2min - commargs[BETAMIN]
+                            );
                 //case: p<0
-                pyranges[2*i+2*maxp].setlower( -asina2b2max - commargs[BETAMAX]);
-                pyranges[2*i+2*maxp].setupper( -asina2b2min - commargs[BETAMIN]);
-                pyranges[2*i+1+2*maxp].setlower(M_PI + asina2b2min - commargs[BETAMAX]);
-                pyranges[2*i+1+2*maxp].setupper(M_PI + asina2b2max - commargs[BETAMIN]);
+                pyranges.add(
+                            -asina2b2max - commargs[BETAMAX],
+                            -asina2b2min - commargs[BETAMIN]
+                            );
+                pyranges.add(
+                            M_PI + asina2b2min - commargs[BETAMAX],
+                            M_PI + asina2b2max - commargs[BETAMIN]
+                            );
             }
         }
         //99 bottles of bugs on the wall, 99 bottles of bugs. You get one down and fix it up, 99 bottles of bugs...
         //100 bottles of bugs on the wall, 100 bottles of bugs....
 
-        std::vector<anglerange> yoverlaps;
-        for(i=0;i<4*maxo+2;i++)
-        {
-            for(j=0;j<4*maxp+2;j++)
-            {
-                anglerange tmp = qyranges[i].overlap(pyranges[j]);
-                if(!(tmp.isempty()))
-                {
-                    yoverlaps.push_back(tmp);
-                }
-            }
-        }
-
-        //to be a commensurate match, an angle has to be in both, x- and yoverlaps
-        std::vector<anglerange> commensurate;
-        for(i=0;i<xoverlaps.size();i++)
-        {
-            for(j=0;j<yoverlaps.size();j++)
-            {
-                anglerange tmp = xoverlaps[i].overlap(yoverlaps[j]);
-                if(!(tmp.isempty()))
-                {
-                    commensurate.push_back(tmp);
-                }
-            }
-        }
-/*
-        cout << "Possible coincident lattice matches with px, qx:" << std::endl;
-        for(i=0;i<xoverlaps.size();i++)
-        {
-            cout << (xoverlaps.at(i)).getlower().getval()*180/M_PI << " " << (xoverlaps.at(i)).getupper().getval()*180/M_PI << std::endl;
-        }
-        cout << "Possible coincident lattice matches with qy, py:" << std::endl;
-        for(i=0;i<yoverlaps.size();i++)
-        {
-            cout << (yoverlaps.at(i)).getlower().getval()*180/M_PI << " " << (yoverlaps.at(i)).getupper().getval()*180/M_PI << std::endl;
-        }
-*/
-        //to make the output more readable, combine the ranges.
-        //first: copy xranges and yranges together:
-        std::vector<anglerange> ranges;
-        ranges.reserve(xoverlaps.size()+yoverlaps.size()+1);
-        ranges.insert(ranges.end(),xoverlaps.begin(),xoverlaps.end());
-        ranges.insert(ranges.end(),yoverlaps.begin(),yoverlaps.end());
-        if(ranges.size()>1) //combining elements only makes sense, if there are at least 2 elements to combine...
-        {
-            bool stop;
-            do
-            {
-                i=0;
-                stop=false;
-                do
-                {
-                    j=i+1;
-                    do
-                    {
-                        anglerange tmp = ranges.at(i).combine(ranges.at(j));
-                        if(!(tmp.isempty()))
-                        {
-                            stop=true;
-                            ranges.push_back(tmp);
-                            ranges.erase(ranges.begin()+j);
-                            ranges.erase(ranges.begin()+i);
-                        }
-                        j++;
-                    }
-                    while(!stop && j<ranges.size());
-                    i++;
-                }
-                while(!stop && i<ranges.size()-1); //the -1 is there, as for i=ranges.size() there is no valid j any more.
-            }
-            while(stop); //yes, here is stop, not !stop, as stop means: stop to iterate over elements, go back to main loop, which should only terminate, if the internal loops finished.
-        }
-
-        //the same thing for commensurate...
-        if(commensurate.size()>1)
-        {
-            bool stop;
-            do
-            {
-                i=0;
-                stop=false;
-                do
-                {
-                    j=i+1;
-                    do
-                    {
-                        anglerange tmp = commensurate.at(i).combine(commensurate.at(j));
-                        if(!(tmp.isempty()))
-                        {
-                            stop=true;
-                            commensurate.push_back(tmp);
-                            commensurate.erase(commensurate.begin()+j);
-                            commensurate.erase(commensurate.begin()+i);
-                        }
-                        j++;
-                    }
-                    while(!stop && j<commensurate.size());
-                    i++;
-                }
-                while(!stop && i<commensurate.size()-1); //the -1 is there, as for i=commensurate.size() there is no valid j any more.
-            }
-            while(stop); //yes, here is stop, not !stop, as stop means: stop to iterate over elements, go back to main loop, which should only terminate, if the internal loops finished.
-        }
+        angleset yoverlaps = pyranges.overlap(qyranges);
 
-        //again,for readability: sort the ranges by their lower limit:
-        {
-            std::sort(ranges.begin(),ranges.end(),[](const anglerange &a,const anglerange &b){return(a.getlower()<b.getlower());});
-            std::sort(commensurate.begin(),commensurate.end(),[](const anglerange &a,const anglerange &b){return(a.getlower()<b.getlower());});
-        }
+        angleset coincident;
+        coincident=xoverlaps;
+        coincident.add(yoverlaps);
 
+        //to be a commensurate match, an angle has to be in both, x- and yoverlaps
+        angleset commensurate = xoverlaps.overlap(yoverlaps);
 
-        cout << "Possible coincident lattice matches:" << std::endl;
-        for(i=0;i<ranges.size();i++)
-        {
-            cout << (ranges.at(i)).getlower().getval()*180/M_PI << " " << (ranges.at(i)).getupper().getval()*180/M_PI << std::endl;
-        }
-        cout << "Possible commensurate lattice matches:" << std::endl;
-        for(i=0;i<commensurate.size();i++)
-        {
-            cout << (commensurate.at(i)).getlower().getval()*180/M_PI << " " << (commensurate.at(i)).getupper().getval()*180/M_PI << std::endl;
-        }
+        //quick and dirrrty
+        coincident.sort();
+        commensurate.sort();
+        cout << "Coincident Matches:\n";
+        for_each(coincident.getrangesref().begin(),coincident.getrangesref().end(),[](anglerange i) {cout << i.getlower().getval()*180/M_PI << " " << i.getupper().getval()*180/M_PI << "\n"; });
+        cout << "Commensurate Matches:\n";
+        for_each(commensurate.getrangesref().begin(),commensurate.getrangesref().end(),[](anglerange i) {cout << i.getlower().getval()*180/M_PI << " " << i.getupper().getval()*180/M_PI << "\n"; });
     }
 }