diff options
| -rw-r--r-- | CMakeLists.txt | 4 | ||||
| -rw-r--r-- | README.md | 9 | ||||
| -rw-r--r-- | angleclass.cpp | 42 | ||||
| -rw-r--r-- | angleclass.h | 22 | ||||
| -rw-r--r-- | anglerange.cpp | 227 | ||||
| -rw-r--r-- | anglerange.h | 31 | ||||
| -rw-r--r-- | angleset.cpp | 188 | ||||
| -rw-r--r-- | angleset.h | 91 | ||||
| -rw-r--r-- | main.cpp | 393 |
9 files changed, 745 insertions, 262 deletions
diff --git a/CMakeLists.txt b/CMakeLists.txt index 888ad9f..5add3ab 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,8 +3,8 @@ cmake_minimum_required(VERSION 2.8) aux_source_directory(. SRC_LIST) add_executable(${PROJECT_NAME} ${SRC_LIST}) -SET(CMAKE_CXX_FLAGS "-Wall ${CMAKE_CXX_FLAGS} -std=c++11") -SET(CMAKE_CXX_FLAGS_DEBUG "-Wall ${CMAKE_CXX_FLAGS_DEBUG} -std=c++11") +SET(CMAKE_CXX_FLAGS "-Wall ${CMAKE_CXX_FLAGS} -std=c++11 -DNDEBUG") +SET(CMAKE_CXX_FLAGS_DEBUG "-Wall ${CMAKE_CXX_FLAGS_DEBUG} -std=c++11 -O0") IF(UNIX) TARGET_LINK_LIBRARIES(${PROJECT_NAME} m) @@ -62,6 +62,15 @@ Both points together made writing the correct conditions for upper and lower bou of the result ranges a pain, and I'm pretty certain there might still be the one or the other bug hidden. +##Compilation +The easiest way to compile this program is by using the included cmake file. Simply run cmake, +followed by make. The default build should use sane CFLAGS, and have the NDEBUG macro defined. + +If you, for whatever reason, prefer to compile manually, just link all the +.cpp files together in one executable. If you don't want to waste a lot of processor cycles, +make sure to define the NDEBUG macro, as otherwise some quite CPU-heavy assertion checks are in +place. For most compilers this can be done by passing the "-DNDEBUG" command line argument. + ##Usage Using this program is quite easy: Just supply the input as command line parameters in this order: a1, a2, alpha, b1min, b1max, b2min, b2max, betamin, betamax diff --git a/angleclass.cpp b/angleclass.cpp index e1981d1..8f30d69 100644 --- a/angleclass.cpp +++ b/angleclass.cpp @@ -52,58 +52,58 @@ void angleclass::setval(const double setme) value=shiftinrange(setme); } -double angleclass::getval() +double angleclass::getval() const { return(value); } -angleclass angleclass::operator *(angleclass other) +angleclass angleclass::operator *(const angleclass other) const { - return(angleclass(value*other.getval())); + return(angleclass(value*other.value)); } -angleclass angleclass::operator-(angleclass other) +angleclass angleclass::operator-(const angleclass other) const { - return(angleclass(value-other.getval())); + return(angleclass(value-other.value)); } -angleclass angleclass::operator+(angleclass other) +angleclass angleclass::operator+(const angleclass other) const { - return(angleclass(value+other.getval())); + return(angleclass(value+other.value)); } -angleclass angleclass::operator/(angleclass other) +angleclass angleclass::operator/(const angleclass other) const { - return(angleclass(value/other.getval())); + return(angleclass(value/other.value)); } -bool angleclass::operator<(angleclass other) +bool angleclass::operator<(const angleclass other) const { - return(value<other.getval()); + return(value<other.value); } -bool angleclass::operator>(angleclass other) +bool angleclass::operator>(const angleclass other) const { - return(value>other.getval()); + return(value>other.value); } -bool angleclass::operator<=(angleclass other) +bool angleclass::operator<=(const angleclass other) const { - return(value<=other.getval()); + return(value<=other.value); } -bool angleclass::operator>=(angleclass other) +bool angleclass::operator>=(const angleclass other) const { - return(value>=other.getval()); + return(value>=other.value); } -bool angleclass::operator==(angleclass other) +bool angleclass::operator==(const angleclass other) const { - return(value==other.getval()); + return(value==other.value); } -bool angleclass::operator!=(angleclass other) +bool angleclass::operator!=(const angleclass other) const { - return(value!=other.getval()); + return(value!=other.value); } diff --git a/angleclass.h b/angleclass.h index 5b2e99e..8135c0a 100644 --- a/angleclass.h +++ b/angleclass.h @@ -39,19 +39,19 @@ public: angleclass(const double setme); void setval(const double setme); - double getval(); + double getval() const; - angleclass operator+(angleclass other); - angleclass operator-(angleclass other); - angleclass operator*(angleclass other); - angleclass operator/(angleclass other); + angleclass operator+(const angleclass other) const; + angleclass operator-(const angleclass other) const; + angleclass operator*(const angleclass other) const; + angleclass operator/(const angleclass other) const; //beware: these only compare the numeric value of the angle. - bool operator<(angleclass other); - bool operator>(angleclass other); - bool operator<=(angleclass other); - bool operator>=(angleclass other); - bool operator==(angleclass other); - bool operator!=(angleclass other); + bool operator<(const angleclass other) const; + bool operator>(const angleclass other) const; + bool operator<=(const angleclass other) const; + bool operator>=(const angleclass other) const; + bool operator==(const angleclass other) const; + bool operator!=(const angleclass other) const; }; #endif // ANGLECLASS_H diff --git a/anglerange.cpp b/anglerange.cpp index 7dda2e9..c62f760 100644 --- a/anglerange.cpp +++ b/anglerange.cpp @@ -43,6 +43,8 @@ anglerange::anglerange() upperborder = angleclass(0.0); upperset = false; lowerset = false; + fullcircle= false; + sortby = SRT_SIZE; } anglerange::anglerange(const angleclass &lower, const angleclass &upper) @@ -51,6 +53,8 @@ anglerange::anglerange(const angleclass &lower, const angleclass &upper) upperborder = upper; upperset = true; lowerset = true; + fullcircle = false; + sortby = SRT_SIZE; } bool anglerange::isempty() const @@ -84,17 +88,81 @@ void anglerange::setempty() { upperset = false; lowerset = false; - lowerborder = 0.0; - upperborder = 0.0; + fullcircle = false; +// lowerborder = 0.0; +// upperborder = 0.0; +//commented out, as users (I) should not rely on this. +} + +void anglerange::setsorttype(const anglerange::sorttype value) +{ + sortby=value; +} + +const anglerange::sorttype anglerange::getsorttype() +{ + return(sortby); +} + +bool anglerange::iscircle() const +{ + return(fullcircle && !isempty()); +} + +void anglerange::setcircle(bool value) +{ + if(value) + { + lowerset=true; + upperset=true; + lowerborder=0.0; + upperborder=0.0; + } + fullcircle=value; +} + +bool anglerange::isinside(const angleclass &val) const +{ + if(isempty()) + { + return(false); + } + else + { + if(fullcircle) + { + return(true); + } + else if(lowerborder>upperborder) + { + return((val>=lowerborder || val<=upperborder)); + } + else + { + return(val>=lowerborder && val<=upperborder); + } + } } anglerange anglerange::overlap(const anglerange &other) { //storage for the return value: anglerange retval; //anglerange constructor without arguments: emtpy range [0:0] + retval.setsorttype(sortby); //return value inherits current sorttype. //first: if one of the input ranges is empty, don't do anything, leave retval empty! + if(!(isempty()) && !(other.isempty())) { + //special treatment if one of the ranges is a full circle: + if(fullcircle) + { + retval = other; + retval.setsorttype(sortby); + } + else if(other.fullcircle) + { + retval = *this; + } //Here we need to take care to always follow the counter-clockwise convention. //This leaves us with four possibilities: //a) Both ranges contain the zero-line @@ -103,7 +171,7 @@ anglerange anglerange::overlap(const anglerange &other) //d) Both ranges are regular. //first check if this range contains the zero-line - if(lowerborder>upperborder) + else if(lowerborder>upperborder) { //ok, this range is around zero. Let's check the other one too: if(other.getlower()>other.getupper()) @@ -174,8 +242,19 @@ anglerange anglerange::overlap(const anglerange &other) anglerange anglerange::combine(const anglerange &other) { anglerange retval; + retval.setsorttype(sortby); //return value inherits current sorttype. //check if any of the ranges is empty. If yes: return an empty range as well. if(!(isempty()) && !(other.isempty())){ + //again: special treatment if one of the ranges is a full circle: + if(fullcircle) + { + retval = *this; + } + if(other.fullcircle) + { + retval = other; + retval.setsorttype(sortby); + } //Here we need to take care to always follow the counter-clockwise convention. //This leaves us with four possibilities: //a) Both ranges contain the zero-line @@ -184,14 +263,23 @@ anglerange anglerange::combine(const anglerange &other) //d) Both ranges are regular. //first check if this range contains the zero-line - if(lowerborder>upperborder) + else if(lowerborder>upperborder) { //ok, this range is around zero. Let's check the other one too: if(other.getlower()>other.getupper()) { //ok, both ranges contain the zero-line. - retval.setlower(fmin(lowerborder.getval(),other.getlower().getval())); - retval.setupper(fmax(upperborder.getval(),other.getupper().getval())); + //let's see if the result is a full circle + if(lowerborder<=other.upperborder || other.lowerborder <=upperborder) + { + //full circle + retval.setcircle(true); + } + else + { + retval.setlower(fmin(lowerborder.getval(),other.getlower().getval())); + retval.setupper(fmax(upperborder.getval(),other.getupper().getval())); + } } else { @@ -203,6 +291,11 @@ anglerange anglerange::combine(const anglerange &other) retval.setlower(lowerborder); retval.setupper(upperborder); } + //again: check if result is a full circle + else if(other.lowerborder<=upperborder && other.upperborder>=lowerborder) + { + retval.setcircle(true); + } else if(other.getupper()>=lowerborder) //not >, as per definition upperborder is inside range { retval.setlower(other.getlower()); @@ -226,6 +319,11 @@ anglerange anglerange::combine(const anglerange &other) retval.setlower(other.getlower()); retval.setupper(other.getupper()); } + //full circle? + else if(lowerborder <= other.getupper() && upperborder >= other.getlower()) + { + retval.setcircle(true); + } else if(lowerborder <= other.getupper()) { retval.setlower(other.getlower()); @@ -239,6 +337,7 @@ anglerange anglerange::combine(const anglerange &other) } else //both ranges regular { + //this also means: none of them contains zero, the result cannot be a full circle //check if there's an overlap double curmax = fmin(other.getupper().getval(),upperborder.getval()); double curmin = fmax(other.getlower().getval(),lowerborder.getval()); @@ -252,3 +351,119 @@ anglerange anglerange::combine(const anglerange &other) } return(retval); } + +bool anglerange::operator<(const anglerange &other) const +{ + if(isempty() || other.isempty()) + { + return(false); + } + else + { + switch(sortby) + { + case SRT_LOWER: + return(lowerborder < other.lowerborder); + break; + case SRT_UPPER: + return(upperborder < other.upperborder); + break; + case SRT_SIZE: + //I know that this if elseif thing can be written as a single statement. This is for readability. + if(other.fullcircle) + { + return(!fullcircle); + } + else if(fullcircle) + { + return(false); + } + else + { + //angleclass takes care of zero crossing. C++ is awesome. + return(upperborder - lowerborder < other.upperborder - other.lowerborder); + } + break; + default: + throw std::out_of_range("Invalid sort field set for anglerange! Go, debug.\n"); + } + } +} + +bool anglerange::operator>(const anglerange &other) const +{ + //as == doesn't care about sort order but fully compares, we have to implement two sort operators fully... + if(isempty() || other.isempty()) + { + return(false); + } + else + { + switch(sortby) + { + case SRT_LOWER: + return(lowerborder > other.lowerborder); + break; + case SRT_UPPER: + return(upperborder > other.upperborder); + break; + case SRT_SIZE: + if(fullcircle) + { + return(!other.fullcircle); + } + else if(other.fullcircle) + { + return(false); + } + else + { + //angleclass takes care of zero crossing. C++ is awesome. + return(upperborder - lowerborder > other.upperborder - other.lowerborder); + } + break; + default: + throw std::out_of_range("Invalid sort field set for anglerange! Go, debug.\n"); + } + } +} + +bool anglerange::operator<=(const anglerange &other) const +{ + if(isempty() || other.isempty()) + { + return(false); + } + else + { + return(!(operator>(other))); + } +} +bool anglerange::operator>=(const anglerange &other) const +{ + if(isempty() || other.isempty()) + { + return(false); + } + else + { + return(!(operator<(other))); + } +} + + +bool anglerange::operator==(const anglerange &other) const +{ + if(!(isempty()) && !(other.isempty())) + { + //both ranges set, we need to compare field by field + return( (lowerborder == other.lowerborder) && (upperborder == other.upperborder) ); + } + else + return(isempty() && other.isempty()); +} +bool anglerange::operator!=(const anglerange &other) const +{ + //needn't check for isempty here, as == already checks for it. + return(!operator==(other)); +} diff --git a/anglerange.h b/anglerange.h index ab98cbc..6d6d286 100644 --- a/anglerange.h +++ b/anglerange.h @@ -38,14 +38,22 @@ #define ANGLERANGE_H #include "angleclass.h" +#include <stdexcept> class anglerange { +public: + typedef enum sorttype {SRT_LOWER,SRT_UPPER,SRT_SIZE} sorttype; private: angleclass lowerborder; angleclass upperborder; bool upperset; bool lowerset; + bool fullcircle; + sorttype sortby; //for comparison operators + //default value: SRT_SIZE + //lower means sort by lower border, upper sort by upper border, and size sort by (upper-lower) + //only size cares about zero-crossing, the others just compare numeric values. public: //Default constructor: Both limits 0.0, marked as empty anglerange(); @@ -54,9 +62,18 @@ public: void setupper(const angleclass &newupper); void setempty(); //marks the range as empty. + //to change sort field + void setsorttype(const sorttype value); + const sorttype getsorttype(); + bool isempty() const; //returns true when the range is empty (or when one limit isn't set) - angleclass getlower() const; //warning: Does not check if empty - angleclass getupper() const; //warning: Does not check if empty + angleclass getlower() const; //warning: Does not check if empty, does not check if circle + angleclass getupper() const; //warning: Does not check if empty, does not check if circle + bool isinside(const angleclass &val) const; //check if angleclas val is in the range. + + //to deal with full circles: + void setcircle(bool value); + bool iscircle() const; //this function calculates the overlap between this anglerange and another one, returning it as //a new anglerange. @@ -66,6 +83,16 @@ public: //if they are disjoint, an empty range is given back. anglerange combine(const anglerange &other); + bool operator<(const anglerange &other) const; + bool operator>(const anglerange &other) const; + bool operator<=(const anglerange &other) const; + bool operator>=(const anglerange &other) const; + + //these operators don't use sort order. They really compare by element! + //also: Two empty ranges are considered equal! + bool operator==(const anglerange &other) const; + bool operator!=(const anglerange &other) const; + }; #endif // ANGLERANGE_H diff --git a/angleset.cpp b/angleset.cpp new file mode 100644 index 0000000..fc759b8 --- /dev/null +++ b/angleset.cpp @@ -0,0 +1,188 @@ +/* + * LatticeMatch calculator - class used for angles + * Copyright (C) 2015 Andreas Grois + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + * + * + * This class extends anglerange with the ability to deal with disjoint ranges. + * + * The internal storage isn't exposed, and there is no way to address an individual sub-range. + * The reason for this is that the functions to add or remove ranges will resort the contents + * of the internal storage arbitrarily. + * + * This class is part of the LatticeMatch program. + * + * To contact the author either use electronic mail: Andreas Grois <andreas.grois@jku.at> + * or write to: + * Andreas Grois, Institute for Semiconductor and SolidState physics, Johannes Kepler University, + * Altenbergerstraße 69, 4040 Linz, AUSTRIA + */ + +#include "angleset.h" + +angleset::angleset() +{ + //storage should be initialized as an empty vector, all we need here is: + consistent=true; +} + +angleset::angleset(const anglerange &firstrange) +{ + assert(storage.size()==0); + if(!firstrange.isempty()) + { + storage.push_back(firstrange); + storage[0].setsorttype(anglerange::SRT_LOWER); + } + consistent=true; +} + +angleset::angleset(const angleclass &firstlower, const angleclass &firstupper) +{ + assert(storage.size()==0); + storage.push_back(anglerange(firstlower,firstupper)); + storage[0].setsorttype(anglerange::SRT_LOWER); + consistent=true; +} + +void angleset::reserve(size_t n) +{ + storage.reserve(n); +} + +void angleset::combine() +{ + //the following should hopefully also work when storage is empty. + if(storage.size()>=2){ //only do something if there are at least two elements + for(std::vector<anglerange>::iterator curpos=storage.end()-2; //second last element + curpos>=storage.begin(); //learned the hard way, that for an empty vector .begin() returns (unsigned)0, so .begin()-2>.begin()... + curpos--){ + std::vector<anglerange>::iterator compare = storage.end()-1; + while(compare!=curpos) + { + anglerange tmp = curpos->combine(*compare); + if(!(tmp.isempty())) + { + *curpos=tmp; + compare--; //I'm not certain if it's necessary to first decrement before deleting, + //but better safe than sorry. + storage.erase(compare+1); + } + else + compare--; + } + } + } + consistent=true; +} + +bool angleset::isempty() const +{ + return(storage.empty()); +} + +bool angleset::iscircle() +{ + if(!consistent) + combine(); + //if only one element in storage should be a full circle, combine will reduce storage to a single element vector. + if(!storage.empty()) + { + return(storage.front().iscircle()); + } + else{ + return(false); + } +} + +void angleset::add(const anglerange &value) +{ + //quick and dirty: We add the range, sort through, and combine if necessarry. + //this means: we *should not reuse* this function for adding anglesets. + if(!value.isempty()) + { + storage.push_back(value); + storage.back().setsorttype(anglerange::SRT_LOWER); + combine(); + } +} + +void angleset::add(const angleset &value) +{ + storage.reserve(storage.size()+value.storage.size()); + storage.insert(storage.end(),value.storage.begin(),value.storage.end()); + //if value is a valid angleset, sort order is set to SRT_LOWER for all fields. + combine(); +} + +void angleset::add(const double &lower, const double &upper) +{ + storage.push_back(anglerange(lower,upper)); + storage.back().setsorttype(anglerange::SRT_LOWER); + combine(); +} + + +angleset angleset::overlap(const anglerange &other) +{ + if(!consistent) + combine(); + angleset retval; + for_each( + storage.begin(), + storage.end(), + [&](anglerange i) + { + retval.add(i.overlap(other)); //given that add checks if the value passed to it is empty: this should work?!? + } + ); + return(retval); +} + +angleset angleset::overlap(const angleset &other) +{ + //just add the overlap of each anglerange in other ;-) + angleset retval; + for_each( + other.storage.begin(), + other.storage.end(), + [&](anglerange i) + { + retval.add(overlap(i)); + } + ); + return(retval); +} + +std::vector<anglerange> angleset::getranges() +{ + return storage; +} + +const std::vector<anglerange>& angleset::getrangesref() +{ + return storage; +} + +void angleset::clear() +{ + storage.clear(); +} + +void angleset::sort() +{ + std::sort(storage.begin(),storage.end()); +} diff --git a/angleset.h b/angleset.h new file mode 100644 index 0000000..781a499 --- /dev/null +++ b/angleset.h @@ -0,0 +1,91 @@ +/* + * LatticeMatch calculator - class used for angles + * Copyright (C) 2015 Andreas Grois + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + * + * + * This class extends anglerange with the ability to deal with disjoint ranges. + * + * The internal storage isn't exposed, and there is no way to address an individual sub-range. + * The reason for this is that the functions to add or remove ranges will resort the contents + * of the internal storage arbitrarily. + * + * This class is part of the LatticeMatch program. + * + * To contact the author either use electronic mail: Andreas Grois <andreas.grois@jku.at> + * or write to: + * Andreas Grois, Institute for Semiconductor and SolidState physics, Johannes Kepler University, + * Altenbergerstraße 69, 4040 Linz, AUSTRIA + */ + +#ifndef ANGLESET_H +#define ANGLESET_H +#include <vector> +#include <cassert> +#include <algorithm> +#include "anglerange.h" + +class angleset +{ +private: + std::vector<anglerange> storage; + bool consistent; + void combine(); +public: + angleset(); + angleset(const anglerange &firstrange); + angleset(const angleclass &firstlower, const angleclass &firstupper); + + //returns if the set is empty. Currently this means: internal storage is an empty vector. + bool isempty() const; + bool iscircle();//returns true when the whole angleset is a circle. + + //add or remove single ranges + void add(const anglerange &value); + //void remove(const anglerange &value); //not implemented yet + void add(const double &lower, const double &upper); + //void remove(const double &lower, const double &upper); //not implemented yet + + //here it gets interesting: add or remove complete sets. + void add(const angleset &value); + //void remove(const angleset &value); //not implemented yet + + //in planning: add_deferred functions, that set consistent to false, and do not call combine() at the end. + //but for this of course each and every other function has to check for consistency + + //having this return a new angleset is not consistent with the above add/remove functions, but it + //is consistent to older code in anglerange.h + angleset overlap(const anglerange &other); + angleset overlap(const angleset &other); + + void reserve(size_t n); //just forwards storage's reserve function. + + //empties the range + void clear(); + + //sorts the range + void sort(); + + //*generates* a new vector consisting of non-overlapping, unique angleranges. + std::vector<anglerange> getranges(); + //gives you a reference to the internal storage, but const + //reference, to keep this relatively fast, const so you don't accidentally modify it. + //If the internal storage format changes, this might be changed to give a copy instead. + const std::vector<anglerange>& getrangesref(); //I hope this does what I think it does... + +}; + +#endif // ANGLESET_H @@ -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"; }); } } |