"Fossies" - the Fresh Open Source Software Archive

Member "fityk-1.3.1/wxgui/ceria.cpp" (12 Aug 2016, 24737 Bytes) of package /linux/misc/fityk-1.3.1.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "ceria.cpp" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 1.3.0_vs_1.3.1.

    1 // Author: Marcin Wojdyr
    2 // Licence: GNU General Public License ver. 2+
    3 
    4 #include "ceria.h"
    5 
    6 #ifdef _MSC_VER
    7 #define _USE_MATH_DEFINES
    8 #endif
    9 #include <stdio.h>
   10 #include <ctype.h>
   11 #include <cmath>
   12 #include <assert.h>
   13 #include <string.h>
   14 #include <algorithm>
   15 
   16 using namespace std;
   17 
   18 static
   19 Pos apply_seitz(Pos const& p0, SeitzMatrix const& s)
   20 {
   21     Pos p;
   22     p.x = s.R[0]*p0.x + s.R[1]*p0.y + s.R[2]*p0.z + s.T[0]/12.;
   23     p.y = s.R[3]*p0.x + s.R[4]*p0.y + s.R[5]*p0.z + s.T[1]/12.;
   24     p.z = s.R[6]*p0.x + s.R[7]*p0.y + s.R[8]*p0.z + s.T[2]/12.;
   25     return p;
   26 }
   27 
   28 
   29 string fullHM(const SpaceGroupSetting *sgs)
   30 {
   31     if (sgs == NULL)
   32         return "";
   33     else if (sgs->ext == 0)
   34         return sgs->HM;
   35     else
   36         return string(sgs->HM) + ":" + sgs->ext;
   37 }
   38 
   39 static
   40 bool check_symmetric_hkl(const SgOps &sg_ops, const Miller &p1,
   41                                               const Miller &p2)
   42 {
   43     for (size_t i = 0; i < sg_ops.seitz.size(); ++i)
   44     {
   45         const int* R = sg_ops.seitz[i].R;
   46         int h = R[0] * p1.h + R[3] * p1.k + R[6] * p1.l;
   47         int k = R[1] * p1.h + R[4] * p1.k + R[7] * p1.l;
   48         int l = R[2] * p1.h + R[5] * p1.k + R[8] * p1.l;
   49         if (h == p2.h && k == p2.k && l == p2.l)
   50             return  true;
   51         // we assume Friedel symmetry here
   52         if (h == -p2.h && k == -p2.k && l == -p2.l)
   53             return true;
   54     }
   55     return false;
   56 }
   57 
   58 static
   59 bool is_position_empty(const vector<Pos>& pp, const Pos& p)
   60 {
   61     const double eps = 0.01;
   62     for (vector<Pos>::const_iterator j = pp.begin(); j != pp.end(); ++j) {
   63         if ((fabs(p.x - j->x) < eps || fabs(p.x - j->x) > 1 - eps) &&
   64             (fabs(p.y - j->y) < eps || fabs(p.y - j->y) > 1 - eps) &&
   65             (fabs(p.z - j->z) < eps || fabs(p.z - j->z) > 1 - eps))
   66             return false;
   67     }
   68     return true;
   69 }
   70 
   71 static
   72 double mod1(double x) { return x - floor(x); }
   73 
   74 void add_symmetric_images(Atom& a, const SgOps& sg_ops)
   75 {
   76     assert(a.pos.size() == 1);
   77     // iterate over translation, Seitz matrices and inversion
   78     for (size_t nt = 0; nt != sg_ops.tr.size(); ++nt) {
   79         const TransVec& t = sg_ops.tr[nt];
   80         for (size_t ns = 0; ns != sg_ops.seitz.size(); ++ns) {
   81             Pos ps = apply_seitz(a.pos[0], sg_ops.seitz[ns]);
   82             Pos p = { mod1(ps.x + t.x/12.),
   83                       mod1(ps.y + t.y/12.),
   84                       mod1(ps.z + t.z/12.) };
   85             if (is_position_empty(a.pos, p)) {
   86                 a.pos.push_back(p);
   87             }
   88             if (sg_ops.inv) {
   89                 Pos p2 = { mod1(-ps.x + t.x/12.),
   90                            mod1(-ps.y + t.y/12.),
   91                            mod1(-ps.z + t.z/12.) };
   92                 if (is_position_empty(a.pos, p2)) {
   93                     a.pos.push_back(p2);
   94                 }
   95             }
   96         }
   97     }
   98 }
   99 
  100 int parse_atoms(const char* s, Crystal& cr)
  101 {
  102     int line_with_error = -1;
  103 
  104     int atom_count = 0;
  105     // to avoid not necessary copying, don't use atoms.clear()
  106     for (int line_nr = 1; ; ++line_nr) {
  107         // skip whitespace
  108         while(isspace(*s) && *s != '\n')
  109             ++s;
  110         if (*s == '\n') {
  111             ++s;
  112             continue;
  113         }
  114         if (*s == '\0')
  115             break;
  116 
  117         // usually the atom is not changed, so we first parse data
  118         // into a new struct Atom, and if it is different we copy it
  119         // and do calculations
  120         Atom a;
  121 
  122         // parse symbol
  123         const char* word_end = s;
  124         while (isalnum(*word_end))
  125             ++word_end;
  126         int symbol_len = word_end - s;
  127         if (symbol_len == 0 || symbol_len >= 8) {
  128             line_with_error = line_nr;
  129             break;
  130         }
  131         memcpy(a.symbol, s, symbol_len);
  132         a.symbol[symbol_len] = '\0';
  133         s = word_end;
  134 
  135         // parse x, y, z
  136         char *endptr;
  137         a.pos[0].x = strtod(s, &endptr);
  138         s = endptr;
  139         a.pos[0].y = strtod(s, &endptr);
  140         s = endptr;
  141         a.pos[0].z = strtod(s, &endptr);
  142 
  143         // check if the numbers were parsed
  144         if (endptr == s || // one or more numbers were not parsed
  145             (!isspace(*endptr) && *endptr != '\0')) // e.g. "Si 0 0 0foo"
  146         {
  147             line_with_error = line_nr;
  148             break;
  149         }
  150 
  151         s = endptr;
  152 
  153         // if there is more than 4 words in the line, ignore the extra words
  154         while (*s != '\n' && *s != '\0')
  155             ++s;
  156         if (*s == '\n')
  157             ++s;
  158 
  159         // check if the atom needs to be copied, and copy it if necessary
  160         if (atom_count == (int) cr.atoms.size()) {
  161             a.xray_sf = find_in_it92(a.symbol);
  162             a.neutron_sf = find_in_nn92(a.symbol);
  163             cr.atoms.push_back(a);
  164             add_symmetric_images(cr.atoms[atom_count], cr.sg_ops);
  165         } else {
  166             Atom &b = cr.atoms[atom_count];
  167             if (strcmp(a.symbol, b.symbol) != 0) {
  168                 memcpy(b.symbol, a.symbol, 8);
  169                 b.xray_sf = find_in_it92(b.symbol);
  170                 b.neutron_sf = find_in_nn92(b.symbol);
  171             }
  172             if (a.pos[0].x != b.pos[0].x
  173                     || a.pos[0].y != b.pos[0].y
  174                     || a.pos[0].z != b.pos[0].z) {
  175                 b.pos.resize(1);
  176                 b.pos[0] = a.pos[0];
  177                 add_symmetric_images(b, cr.sg_ops);
  178             }
  179         }
  180         ++atom_count;
  181     }
  182     // vector cr.atoms may be longer than necessary
  183     cr.atoms.resize(atom_count);
  184     return line_with_error;
  185 }
  186 
  187 
  188 void PlanesWithSameD::add(Miller const& hkl, const SgOps& sg_ops)
  189 {
  190     for (vector<Plane>::iterator i = planes.begin(); i != planes.end(); ++i) {
  191         if (check_symmetric_hkl(sg_ops, *i, hkl)) {
  192             i->multiplicity++;
  193             return;
  194         }
  195     }
  196     // equivalent plane not found
  197     planes.push_back(Plane(hkl));
  198 }
  199 
  200 Crystal::Crystal()
  201     : uc(NULL)
  202 {
  203     atoms.reserve(16);
  204     sg_ops.tr.push_back(TransVec(0, 0, 0)); // always keep trivial tr
  205     // add unit seitz matrix
  206     SeitzMatrix sm = { { 1, 0, 0,
  207                          0, 1, 0,
  208                          0, 0, 1, },
  209                        { 0, 0, 0 } };
  210     sg_ops.seitz.push_back(sm);
  211 }
  212 
  213 Crystal::~Crystal()
  214 {
  215     delete uc;
  216 }
  217 
  218 // in the same order as in enum CrystalSystem
  219 const char *CrystalSystemNames[] = {
  220     "Undefined", // 0
  221     NULL,
  222     "Triclinic", // 2
  223     "Monoclinic",
  224     "Orthorhombic",
  225     "Tetragonal",
  226     "Trigonal",
  227     "Hexagonal",
  228     "Cubic" // 8
  229 };
  230 
  231 const char* get_crystal_system_name(CrystalSystem xs)
  232 {
  233     return CrystalSystemNames[xs];
  234 }
  235 
  236 CrystalSystem get_crystal_system(int space_group)
  237 {
  238     if (space_group <= 2)
  239         return TriclinicSystem;
  240     else if (space_group <= 15)
  241         return MonoclinicSystem;
  242     else if (space_group <= 74)
  243         return OrthorhombicSystem;
  244     else if (space_group <= 142)
  245         return TetragonalSystem;
  246     else if (space_group <= 167)
  247         return TrigonalSystem;
  248     else if (space_group <= 194)
  249         return HexagonalSystem;
  250     else
  251         return CubicSystem;
  252 }
  253 
  254 int get_sg_order(const SgOps& sg_ops)
  255 {
  256     return sg_ops.tr.size() * sg_ops.seitz.size() * (sg_ops.inv ? 2 : 1);
  257 }
  258 
  259 static
  260 char parse_sg_extension(const char *symbol, char *qualif)
  261 {
  262     if (symbol == NULL || *symbol == '\0') {
  263         if (qualif != NULL)
  264             qualif[0] = '\0';
  265         return 0;
  266     }
  267     char ext = 0;
  268     while (isspace(*symbol) || *symbol == ':')
  269         ++symbol;
  270     if (isdigit(*symbol) || *symbol == 'R' || *symbol == 'H') {
  271         ext = *symbol;
  272         ++symbol;
  273         while (isspace(*symbol) || *symbol == ':')
  274             ++symbol;
  275     }
  276     if (qualif != NULL) {
  277         strncpy(qualif, symbol, 4);
  278         qualif[4] = '\0';
  279     }
  280     return ext;
  281 }
  282 
  283 static
  284 const SpaceGroupSetting* parse_hm_or_hall(const char *symbol)
  285 {
  286     // copy and 'normalize' symbol (up to ':') to table s
  287     char s[32];
  288     for (int i = 0; i < 32; ++i) {
  289         if (*symbol == '\0' || *symbol == ':') {
  290             s[i] = '\0';
  291             break;
  292         } else if (isspace(*symbol)) {
  293             s[i] = ' ';
  294             ++symbol;
  295             while (isspace(*symbol))
  296                 ++symbol;
  297         } else {
  298             // In HM symbols, first character is upper case letter.
  299             // In Hall symbols, the second character is upper case.
  300             // The first and second chars are either upper or ' ' or '-'.
  301             // The rest of alpha chars is lower case. 
  302             s[i] = (i < 2 ? toupper(*symbol) : tolower(*symbol));
  303             ++symbol;
  304         }
  305     }
  306     // now *symbol is either '\0' or ':'
  307     for (const SpaceGroupSetting *p = space_group_settings;
  308                                                 p->sgnumber != 0; ++p) {
  309         if (strcmp(p->HM, s) == 0) {
  310             if (*symbol == ':') {
  311                 char ext = parse_sg_extension(symbol+1, NULL);
  312                 while (p->ext != ext) {
  313                     ++p;
  314                     if (strcmp(p->HM, s) != 0) // full match not found
  315                         return NULL;
  316                 }
  317             }
  318             return p;
  319         } else if (strcmp(p->Hall + (p->Hall[0] == ' ' ? 1 : 0), s) == 0) {
  320             return p;
  321         }
  322     }
  323     return NULL;
  324 }
  325 
  326 static
  327 const SpaceGroupSetting* find_space_group_setting(int sgn, const char *setting)
  328 {
  329     char qualif[5];
  330     char ext = parse_sg_extension(setting, qualif);
  331     const SpaceGroupSetting *p = find_first_sg_with_number(sgn);
  332     if (p == NULL)
  333         return NULL;
  334     while (p->ext != ext || strcmp(p->qualif, qualif) != 0) {
  335         ++p;
  336         if (p->sgnumber != sgn) // not found
  337             return NULL;
  338     }
  339     return p;
  340 }
  341 
  342 const SpaceGroupSetting* parse_any_sg_symbol(const char *symbol)
  343 {
  344     if (symbol == NULL)
  345         return NULL;
  346     while (isspace(*symbol))
  347         ++symbol;
  348     if (isdigit(*symbol)) {
  349         const char* colon = strchr(symbol, ':');
  350         int sgn = strtol(symbol, NULL, 10);
  351         return find_space_group_setting(sgn, colon);
  352     } else {
  353         return parse_hm_or_hall(symbol);
  354     }
  355 }
  356 
  357 void Crystal::set_space_group(const SpaceGroupSetting* sgs_)
  358 {
  359     sgs = sgs_;
  360     sg_ops.tr.resize(1); // leave only the trivial tr (0, 0, 0)
  361     sg_ops.inv = false;
  362     sg_ops.inv_t.x = sg_ops.inv_t.y = sg_ops.inv_t.z = 0;
  363     if (sgs == NULL)
  364         return;
  365 
  366     int n_seitz = get_seitz_mx_count(sgs);
  367     // the first seitz matrix is unit matrix
  368     sg_ops.seitz.resize(n_seitz + 1);
  369     for (int i = 0; i != n_seitz; ++i)
  370         get_seitz_mx(sgs, i, &sg_ops.seitz[i+1]);
  371 
  372     switch (sgs->Hall[1])
  373     {
  374         case 'A': sg_ops.tr.push_back(TransVec(0,6,6)); break;
  375         case 'B': sg_ops.tr.push_back(TransVec(6,0,6)); break;
  376         case 'C': sg_ops.tr.push_back(TransVec(6,6,0)); break;
  377         case 'I': sg_ops.tr.push_back(TransVec(6,6,6)); break;
  378         case 'P': break;
  379         case 'R': sg_ops.tr.push_back(TransVec(8,4,4));
  380                   sg_ops.tr.push_back(TransVec(4,8,8));
  381                   break;
  382         case 'F': sg_ops.tr.push_back(TransVec(0,6,6));
  383                   sg_ops.tr.push_back(TransVec(6,0,6));
  384                   sg_ops.tr.push_back(TransVec(6,6,0));
  385                   break;
  386         default: assert(0);
  387     }
  388     if (sgs->Hall[0] == '-') {
  389         sg_ops.inv = true;
  390     } else {
  391         const char* t = strstr(sgs->Hall, " -1");
  392         if (t != NULL) {
  393             sg_ops.inv = true;
  394             t += 3;
  395             if      (strcmp(t, "ab") == 0) sg_ops.inv_t = TransVec(6, 6, 0);
  396             else if (strcmp(t, "ac") == 0) sg_ops.inv_t = TransVec(6, 0, 6);
  397             else if (strcmp(t, "bc") == 0) sg_ops.inv_t = TransVec(0, 6, 6);
  398             else if (strcmp(t, "ad") == 0) sg_ops.inv_t = TransVec(9, 3, 3);
  399             else if (strcmp(t, "bw") == 0) sg_ops.inv_t = TransVec(0, 6, 3);
  400             else if (strcmp(t, "d" ) == 0) sg_ops.inv_t = TransVec(3, 3, 3);
  401             else if (strcmp(t, "n" ) == 0) sg_ops.inv_t = TransVec(6, 6, 6);
  402             else assert(0);
  403         }
  404     }
  405 }
  406 
  407 // returns true if exists t in sg_ops.tr, such that: h*(t+T) != n
  408 // used by is_sys_absent()
  409 static
  410 bool has_nonunit_tr(const SgOps& sg_ops, const int* T, int h, int k, int l)
  411 {
  412     for (vector<TransVec>::const_iterator t = sg_ops.tr.begin();
  413                                                 t != sg_ops.tr.end(); ++t)
  414         if (((T[0]+t->x) * h + (T[1]+t->y) * k + (T[2]+t->z) * l) % 12 != 0)
  415             return true;
  416     return false;
  417 }
  418 
  419 static
  420 bool is_sys_absent(const SgOps& sg_ops, int h, int k, int l)
  421 {
  422     for (size_t i = 0; i < sg_ops.seitz.size(); ++i) {
  423         const int* R = sg_ops.seitz[i].R;
  424         const int* T = sg_ops.seitz[i].T;
  425         int M[3] = { h * R[0] + k * R[3] + l * R[6],
  426                      h * R[1] + k * R[4] + l * R[7],
  427                      h * R[2] + k * R[5] + l * R[8] };
  428         if (h == M[0] && k == M[1] && l == M[2]) {
  429             if (has_nonunit_tr(sg_ops, T, h, k, l))
  430                 return true;
  431         } else if (h == -M[0] && k == -M[1] && l == -M[2] && sg_ops.inv) {
  432             int ts[3] = { sg_ops.inv_t.x - T[0],
  433                           sg_ops.inv_t.y - T[1],
  434                           sg_ops.inv_t.z - T[2] };
  435             if (has_nonunit_tr(sg_ops, ts, h, k, l))
  436                 return true;
  437         }
  438     }
  439     return false;
  440 }
  441 
  442 
  443 // helper to generate sequence 0, 1, -1, 2, -2, 3, ...
  444 static int inc_neg(int h) { return h > 0 ? -h : -h+1; }
  445 
  446 void Crystal::generate_reflections(double min_d)
  447 {
  448     bp.clear();
  449     UnitCell reciprocal = uc->get_reciprocal();
  450 
  451     // set upper limit for iteration of Miller indices
  452     // TODO: smarter algorithm, like in uctbx::unit_cell::max_miller_indices()
  453     int max_h = 20;
  454     int max_k = 20;
  455     int max_l = 20;
  456     if (fabs(uc->alpha - M_PI/2) < 1e-9 && fabs(uc->beta - M_PI/2) < 1e-9 &&
  457                                           fabs(uc->gamma - M_PI/2) < 1e-9) {
  458         max_h = (int) (uc->a / min_d);
  459         max_k = (int) (uc->b / min_d);
  460         max_l = (int) (uc->c / min_d);
  461     }
  462     // Don't generate too many reflections (it could happen
  463     // when user chooses Q instead of 2T, or puts wrong wavelength)
  464     if (max_h * max_k * max_l > 8000)
  465         max_h = max_k = max_l = 20;
  466 
  467     for (int h = 0; h != max_h+1; h = inc_neg(h))
  468         for (int k = 0; k != max_k+1; k = inc_neg(k))
  469             for (int l = (h==0 && k==0 ? 1 : 0); l != max_l+1; l = inc_neg(l)) {
  470                 double d = 1 / reciprocal.calculate_distance(h, k, l);
  471                 //double d = uc->calculate_d(h, k, l); // the same
  472                 if (d < min_d)
  473                     continue;
  474 
  475                 // check for systematic absence
  476                 if (is_sys_absent(sg_ops, h, k, l))
  477                     continue;
  478 
  479                 Miller hkl = { h, k, l };
  480                 bool found = false;
  481                 for (vector<PlanesWithSameD>::iterator i = bp.begin();
  482                         i != bp.end(); ++i) {
  483                     if (fabs(d - i->d) < 1e-9) {
  484                         i->add(hkl, sg_ops);
  485                         found = true;
  486                         break;
  487                     }
  488                 }
  489                 if (!found) {
  490                     PlanesWithSameD new_p;
  491                     new_p.planes.push_back(Plane(hkl));
  492                     new_p.d = d;
  493                     new_p.lpf = 0.;
  494                     new_p.intensity = 0.;
  495                     new_p.enabled = true;
  496                     bp.push_back(new_p);
  497                 }
  498             }
  499     sort(bp.begin(), bp.end());
  500     old_min_d = min_d;
  501 }
  502 
  503 
  504 // stol = sin(theta)/lambda
  505 static
  506 void set_F2(Plane& p, const vector<Atom>& atoms,
  507                  RadiationType radiation, double stol)
  508 {
  509     // calculating F_hkl, (Pecharsky & Zavalij, eq. (2.89) and (2.103))
  510     // assuming population g = 1
  511     // assuming temperature factor t = 1
  512     double F_real = 0.;
  513     double F_img = 0;
  514     for (vector<Atom>::const_iterator i = atoms.begin(); i != atoms.end(); ++i){
  515         double f = 1.;
  516         if (radiation == kXRay && i->xray_sf)
  517             f = calculate_it92_factor(i->xray_sf, stol*stol);
  518         else if (radiation == kNeutron && i->neutron_sf)
  519             f = i->neutron_sf->bond_coh_scatt_length;
  520         for (vector<Pos>::const_iterator j = i->pos.begin();
  521                                                 j != i->pos.end(); ++j) {
  522             double hx = p.h * j->x + p.k * j->y + p.l * j->z;
  523             F_real += f * cos(2*M_PI * hx);
  524             F_img += f * sin(2*M_PI * hx);
  525         }
  526     }
  527     p.F2 = F_real*F_real + F_img*F_img;
  528     //printf("hkl=(%d %d %d) F=(%g, %g) F2=%g\n",
  529     //       p.h, p.k, p.l, F_real, F_img, p.F2);
  530 }
  531 
  532 static
  533 void set_lpf(PlanesWithSameD &bp, RadiationType radiation, double lambda)
  534 {
  535     if (lambda == 0)
  536         bp.lpf = 1.;
  537     else {
  538         double T = asin(bp.stol() * lambda); // theta
  539         if (radiation == kXRay) {
  540             // for x-rays, we assume K=0.5 and
  541             // LP = (1 + cos(2T)^2) / (cos(T) sin(T)^2)
  542             //  (Pecharsky & Zavalij, eq. (2.70), p. 192)
  543             bp.lpf = (1 + cos(2*T)*cos(2*T)) / (cos(T)*sin(T)*sin(T));
  544         } else if (radiation == kNeutron) {
  545             // Kisi & Howard, Applications of Neutron Powder Diffraction (2.38)
  546             // no polarization only the Lorentz factor:
  547             // L = 1 / (4 sin^2(T) cos(T))
  548             bp.lpf = 1 / (4 * sin(T)*sin(T)*cos(T));
  549             // for TOF diffractometers with a fixed diffraction angle:
  550             // L = d^4 sin(theta)
  551         }
  552     }
  553 }
  554 
  555 void Crystal::update_intensities(RadiationType r, double lambda)
  556 {
  557     if (atoms.empty())
  558         return;
  559     for (vector<PlanesWithSameD>::iterator i = bp.begin(); i != bp.end(); ++i) {
  560         set_lpf(*i, r, lambda);
  561         double t = 0;
  562         for (vector<Plane>::iterator j = i->planes.begin();
  563                                      j != i->planes.end(); ++j) {
  564             set_F2(*j, atoms, r, i->stol());
  565             t += j->multiplicity * j->F2;
  566         }
  567         i->intensity = i->lpf * t;
  568     }
  569 }
  570 
  571 double UnitCell::calculate_V() const
  572 {
  573     // Giacovazzo p.62
  574     double cosA = cos(alpha), cosB = cos(beta), cosG = cos(gamma);
  575     double t = 1 - cosA*cosA - cosB*cosB - cosG*cosG + 2*cosA*cosB*cosG;
  576     return a*b*c * sqrt(t);
  577 }
  578 
  579 /*
  580 double UnitCell::calculate_d(int h, int k, int l) const
  581 {
  582     double sinA=sin(alpha), sinB=sin(beta), sinG=sin(gamma),
  583            cosA=cos(alpha), cosB=cos(beta), cosG=cos(gamma);
  584     return
  585       sqrt((1 - cosA*cosA - cosB*cosB - cosG*cosG + 2*cosA*cosB*cosG)
  586               / (  (h/a*sinA) * (h/a*sinA)  //(h/a*sinA)^2
  587                  + (k/b*sinB) * (k/b*sinB)  //(k/b*sinB)^2
  588                  + (l/c*sinG) * (l/c*sinG)  //(l/c*sinG)^2
  589                  + 2*h*l/a/c*(cosA*cosG-cosB)
  590                  + 2*h*k/a/b*(cosA*cosB-cosG)
  591                  + 2*k*l/b/c*(cosB*cosG-cosA)
  592                 )
  593           );
  594 }
  595 */
  596 
  597 //     [       1/a               0        0   ]
  598 // M = [  -cosG/(a sinG)    1/(b sinG)    0   ]
  599 //     [     a*cosB*         b*cosA*      c*  ]
  600 //
  601 // where A is alpha, B is beta, G is gamma, * means reciprocal space.
  602 // (Giacovazzo, p.68)
  603 void UnitCell::set_M()
  604 {
  605     double sinA=sin(alpha), sinB=sin(beta), sinG=sin(gamma),
  606            cosA=cos(alpha), cosB=cos(beta), cosG=cos(gamma);
  607     M[0][0] = 1/a;
  608     M[0][1] = 0.;
  609     M[0][2] = 0.;
  610     M[1][0] = -cosG/(a*sinG);
  611     M[1][1] = 1/(b*sinG);
  612     M[1][2] = 0.;
  613     M[2][0] = b * c * sinA / V // a*
  614               * (cosA*cosG-cosB) / (sinA*sinG); //cosB*
  615     M[2][1] = a * c * sinB / V // b*
  616               * (cosB*cosG-cosA) / (sinB*sinG); //cosA*
  617     M[2][2] = a * b * sinG / V; // c*
  618 
  619     M_1[0][0] = a;
  620     M_1[0][1] = 0;
  621     M_1[0][2] = 0;
  622     M_1[1][0] = b * cosG;
  623     M_1[1][1] = b * sinG;
  624     M_1[1][2] = 0;
  625     M_1[2][0] = c * cosB;
  626     M_1[2][1] = -c * sinB * (cosB*cosG-cosA) / (sinB*sinG);
  627     M_1[2][2] = 1 / M[2][2];
  628 }
  629 
  630 // returns UnitCell reciprocal to this, i.e. that has parameters a*, b*, ...
  631 // (Giacovazzo, p. 64)
  632 UnitCell UnitCell::get_reciprocal() const
  633 {
  634     double ar = b * c * sin(alpha) / V;
  635     double br = a * c * sin(beta) / V;
  636     double cr = a * b * sin(gamma) / V;
  637     double cosAr = (cos(beta)*cos(gamma)-cos(alpha)) / (sin(beta)*sin(gamma));
  638     double cosBr = (cos(alpha)*cos(gamma)-cos(beta)) / (sin(alpha)*sin(gamma));
  639     double cosGr = (cos(alpha)*cos(beta)-cos(gamma)) / (sin(alpha)*sin(beta));
  640     return UnitCell(ar, br, cr, acos(cosAr), acos(cosBr), acos(cosGr));
  641 }
  642 
  643 // returns 1/|v|, where v = M * [h k l]; 
  644 double UnitCell::calculate_distance(double h, double k, double l) const
  645 {
  646     double v2 = 0.;
  647     for (int i = 0; i != 3; ++i) {
  648         double t = h * M_1[0][i] + k * M_1[1][i] + l * M_1[2][i];
  649         v2 += t*t;
  650     }
  651     return sqrt(v2);
  652 }
  653 
  654 const Anode anodes[] = {
  655     { "Cu", 1.54056, 1.54439 },
  656     { "Cr", 2.28970, 2.29361 },
  657     { "Fe", 1.93604, 1.93998 },
  658     { "Mo", 0.70930, 0.71359 },
  659     { "Ag", 0.55941, 0.56380 },
  660     { "Co", 1.78901, 1.79290 },
  661     { NULL, 0, 0 }
  662 };
  663 
  664 
  665 static
  666 const char* default_cel_files[][2] = {
  667 
  668 {"bSiC",
  669 "cell  4.358 4.358 4.358  90 90 90\n"
  670 "Si   14  0 0 0\n"
  671 "C     6  0.25 0.25 0.25\n"
  672 "rgnr 216"
  673 },
  674 
  675 {"aSiC",
  676 "cell  3.082 3.082 15.123  90 90 120\n"
  677 "SI1  14  0       0       0\n"
  678 "SI2  14  0.3333  0.6667  0.1667\n"
  679 "SI3  14  0.3333  0.6667  0.8333\n"
  680 "C1    6  0.0     0.0     0.125\n"
  681 "C2    6  0.3333  0.6667  0.2917\n"
  682 "C3    6  0.3333  0.6667  0.9583\n"
  683 "rgnr 186",
  684 },
  685 
  686 {"NaCl",
  687 "cell  5.64009 5.64009 5.64009  90 90 90\n"
  688 "Na   11   0   0   0\n"
  689 "Cl   17   0.5 0   0\n"
  690 "rgnr 225"
  691 },
  692 
  693 {"diamond",
  694 "cell  3.5595 3.5595 3.5595  90 90 90\n"
  695 "C     6  0 0 0\n"
  696 "rgnr 227"
  697 },
  698 
  699 {"Si",
  700 "cell  5.4309 5.4309 5.4309  90 90 90\n"
  701 "Si   14  0 0 0\n"
  702 "rgnr 227"
  703 },
  704 
  705 {"CeO2",
  706 "cell  5.41 5.41 5.41  90 90 90\n"
  707 "Ce   58  0   0   0\n"
  708 "O     8  0.25 0.25 0.25\n"
  709 "rgnr 225"
  710 },
  711 
  712 {"Zn",
  713 "cell  2.665 2.665 4.946  90 90 120\n"
  714 "Zn   30  0.33333 0.66667 0.25\n"
  715 "rgnr 194"
  716 },
  717 
  718 {NULL, NULL}
  719 };
  720 
  721 
  722 CelFile read_cel_file(FILE *f)
  723 {
  724     CelFile cel = { 0., 0., 0., 0., 0., 0., NULL, vector<AtomInCell>() };
  725     if (!f)
  726         return cel;
  727     char s[20];
  728     int r = fscanf(f, "%4s %lf %lf %lf %lf %lf %lf",
  729                s, &cel.a, &cel.b, &cel.c, &cel.alpha, &cel.beta, &cel.gamma);
  730     if (r != 7)
  731         return cel;
  732     if (strcmp(s, "cell") != 0 && strcmp(s, "CELL") != 0
  733             && strcmp(s, "Cell") != 0) {
  734         return cel;
  735     }
  736     while (1) {
  737         r = fscanf(f, "%12s", s);
  738         if (r != 1)
  739             return cel;
  740         if (strcmp(s, "RGNR") == 0 || strcmp(s, "rgnr") == 0
  741                 || strcmp(s, "Rgnr") == 0)
  742             break;
  743         AtomInCell atom;
  744         r = fscanf(f, "%d %lf %lf %lf", &atom.Z, &atom.x, &atom.y, &atom.z);
  745         if (r != 4) {
  746             return cel;
  747         }
  748         cel.atoms.push_back(atom);
  749         // skip the rest of the line
  750         for (int c = fgetc(f); c != '\n' && c != EOF; c = fgetc(f))
  751             ;
  752     }
  753     int sgn;
  754     r = fscanf(f, "%d", &sgn);
  755     if (r != 1)
  756         return cel;
  757     if (sgn < 1 || sgn > 230)
  758         return cel;
  759     cel.sgs = find_first_sg_with_number(sgn);
  760     for (int c = fgetc(f); c != '\n' && c != EOF; c = fgetc(f)) {
  761         if (c == ':') {
  762             r = fscanf(f, "%8s", s);
  763             if (r == 1)
  764                 cel.sgs = find_space_group_setting(sgn, s);
  765             break;
  766         } else if (isdigit(c)) {
  767             ungetc(c, f);
  768             int pc_setting;
  769             r = fscanf(f, "%d", &pc_setting);
  770             if (r == 1)
  771                 cel.sgs = get_sg_from_powdercell_rgnr(sgn, pc_setting);
  772             break;
  773         } else if (!isspace(c))
  774             break;
  775     }
  776     return cel;
  777 }
  778 
  779 void write_cel_file(CelFile const& cel, FILE *f)
  780 {
  781     //for (int i = 1; i != 104; ++i)
  782     //    assert(find_Z_in_pse(i)->Z == i);
  783     fprintf(f, "cell  %g %g %g  %g %g %g\n", cel.a, cel.b, cel.c,
  784                                              cel.alpha, cel.beta, cel.gamma);
  785     for (vector<AtomInCell>::const_iterator i = cel.atoms.begin();
  786                                                    i != cel.atoms.end(); ++i) {
  787         const t_pse* pse = find_Z_in_pse(i->Z);
  788         fprintf(f, "%-2s   %2d  %g %g %g\n",
  789                 pse->symbol, i->Z, i->x, i->y, i->z);
  790     }
  791     int sgn = cel.sgs->sgnumber;
  792     fprintf(f, "rgnr %d", sgn);
  793     if (sgn != 1 && (cel.sgs-1)->sgnumber == sgn) {
  794         fprintf(f, " :");
  795         if (cel.sgs->ext != 0)
  796             fprintf(f, "%c", cel.sgs->ext);
  797         fprintf(f, "%s", cel.sgs->qualif);
  798     }
  799     fprintf(f, "\n");
  800 }
  801 
  802 void write_default_cel_files(const char* path_prefix)
  803 {
  804     for (const char*(*s)[2] = default_cel_files; (*s)[0] != NULL; ++s) {
  805         string filename = string(path_prefix) + (*s)[0] + ".cel";
  806         FILE *f = fopen(filename.c_str(), "w");
  807         if (!f)
  808             continue;
  809         fprintf(f, "%s\n", (*s)[1]);
  810         fclose(f);
  811     }
  812 }
  813 
  814