C++ UnitCell类代码示例

static int l_uc_C(lua_State* L)
{
	UnitCell* uc = lua_tounitcell(L, 1);
	if(!uc) return 0;

	lua_pushvector(L, new Vector(uc->C()));
	return 1;
}

void Crystal::generate_reflections(double min_d)
{
    bp.clear();
    UnitCell reciprocal = uc->get_reciprocal();

    // set upper limit for iteration of Miller indices
    // TODO: smarter algorithm, like in uctbx::unit_cell::max_miller_indices()
    int max_h = 20;
    int max_k = 20;
    int max_l = 20;
    if (fabs(uc->alpha - M_PI/2) < 1e-9 && fabs(uc->beta - M_PI/2) < 1e-9 &&
                                          fabs(uc->gamma - M_PI/2) < 1e-9) {
        max_h = (int) (uc->a / min_d);
        max_k = (int) (uc->b / min_d);
        max_l = (int) (uc->c / min_d);
    }
    // Don't generate too many reflections (it could happen
    // when user chooses Q instead of 2T, or puts wrong wavelength)
    if (max_h * max_k * max_l > 8000)
        max_h = max_k = max_l = 20;

    for (int h = 0; h != max_h+1; h = inc_neg(h))
        for (int k = 0; k != max_k+1; k = inc_neg(k))
            for (int l = (h==0 && k==0 ? 1 : 0); l != max_l+1; l = inc_neg(l)) {
                double d = 1 / reciprocal.calculate_distance(h, k, l);
                //double d = uc->calculate_d(h, k, l); // the same
                if (d < min_d)
                    continue;

                // check for systematic absence
                if (is_sys_absent(sg_ops, h, k, l))
                    continue;

                Miller hkl = { h, k, l };
                bool found = false;
                for (vector<PlanesWithSameD>::iterator i = bp.begin();
                        i != bp.end(); ++i) {
                    if (fabs(d - i->d) < 1e-9) {
                        i->add(hkl, sg_ops);
                        found = true;
                        break;
                    }
                }
                if (!found) {
                    PlanesWithSameD new_p;
                    new_p.planes.push_back(Plane(hkl));
                    new_p.d = d;
                    new_p.lpf = 0.;
                    new_p.intensity = 0.;
                    new_p.enabled = true;
                    bp.push_back(new_p);
                }
            }
    sort(bp.begin(), bp.end());
    old_min_d = min_d;
}

/// Returns the lowest and highest d-Value in the list. Uses UnitCell and HKL
/// for calculation to prevent problems with potentially inconsistent d-Values
/// in Peak.
std::pair<double, double> SortHKL::getDLimits(const std::vector<Peak> &peaks,
                                              const UnitCell &cell) const {
  auto dLimitIterators = std::minmax_element(
      peaks.begin(), peaks.end(), [cell](const Peak &lhs, const Peak &rhs) {
        return cell.d(lhs.getHKL()) < cell.d(rhs.getHKL());
      });

  return std::make_pair(cell.d((*dLimitIterators.first).getHKL()),
                        cell.d((*dLimitIterators.second).getHKL()));
}

static int l_uc_addatomreduced(lua_State* L)
{
	UnitCell* uc = lua_tounitcell(L, 1);
	if(!uc) return 0;
	
	Atom* a = lua_toatom(L, 2);
	if(!a) return 0;
	
	uc->addAtomReducedCoorinates(a);
	return 0;
}

// Adjust matrix of current model
void CellAnglesPopup::adjustCurrentMatrix(int angleIndex, double value)
{
	// Get current model and set new angle in cell
	Model* model = parent_.aten().currentModelOrFrame();
	if (model)
	{
		UnitCell cell = model->cell();
		cell.setAngle(angleIndex, value);
		CommandNode::run(Commands::CellAxes, "ddddddddd", cell.parameter(UnitCell::CellAX), cell.parameter(UnitCell::CellAY), cell.parameter(UnitCell::CellAZ), cell.parameter(UnitCell::CellBX), cell.parameter(UnitCell::CellBY), cell.parameter(UnitCell::CellBZ),  cell.parameter(UnitCell::CellCX), cell.parameter(UnitCell::CellCY), cell.parameter(UnitCell::CellCZ)); 
	}
}

static int l_uc_translate(lua_State* L)
{
	UnitCell* uc = lua_tounitcell(L, 1);
	if(!uc) return 0;
	
	int x = lua_tointeger(L, 2);
	int y = lua_tointeger(L, 3);
	int z = lua_tointeger(L, 4);

	uc->translate(x, y, z);
	return 0;
}

static int l_uc_g2r(lua_State* L)
{
	UnitCell* uc = lua_tounitcell(L, 1);
	if(!uc) return 0;
	
	Vector v;
	lua_makevector(L, 2, v);

	Vector b = uc->globalToReduced(v);
	lua_pushvector(L, new Vector(b));
	return 1;
}

static int l_uc_applyoperator(lua_State* L)
{
	UnitCell* uc = lua_tounitcell(L, 1);
	if(!uc) return 0;

	if(lua_isfunction(L, 2))
	{
		uc->applyOperator(L, 2);
		return 0;
	}
	if(lua_isstring(L, 2))
	{
		if(!uc->applyOperator(lua_tostring(L, 2)))
			return luaL_error(L, "Failed to apply `%s'\n", lua_tostring(L, 2));
	}
	return 0;
}

void general_ewald_sum(UnitCell &uc, const Vector3d &a1, const Vector3d &a2, const Vector3d &a3, const Vector3i &Rcutoff)
{
    int N = uc.size();
    const double sigma = M_PI;

    // renormalize electron occupation
    double unrenorm_Ne = 0, Ne = 0;
    for (UnitCell::iterator it = uc.begin(); it < uc.end(); ++it) {
        unrenorm_Ne += it->ne;
        Ne += it->C;
    }
    for (UnitCell::iterator it = uc.begin(); it < uc.end(); ++it) 
        it->ne *= Ne/unrenorm_Ne;

    // get reciprocal vectors
    Vector3d b1, b2, b3;
    double uc_vol = a1.dot(a2.cross(a3));
    b1 = 2*M_PI*a2.cross(a3)/uc_vol;
    b2 = 2*M_PI*a3.cross(a1)/uc_vol;
    b3 = 2*M_PI*a1.cross(a2)/uc_vol;
    
    // Ewald sum
    double Vs, Vl;
    Vector3d k;
    for (int id = 0; id < N; ++id) {
        Vs = 0; Vl = 0;
        for (int m = -Rcutoff[0]; m < Rcutoff[0]; ++m)
            for (int n = -Rcutoff[1]; n < Rcutoff[1]; ++n)
                for (int l = -Rcutoff[2]; l < Rcutoff[2]; ++l) {
                    k = n*b1 + m*b2 + l*b3;
                    double ksquare = k.dot(k);
                    for (int i = 0; i < N; ++i) {
                        double dR = (uc[id].r - (uc[i].r + m*a1 + n*a2 + l*a3)).norm();
                        // for long-range terms
                        if ( !(m == 0 && n == 0 && l == 0) )
                            Vl += 4*M_PI/uc_vol*(uc[i].C-uc[i].ne)/ksquare * cos(k.dot(uc[id].r-uc[i].r)) * exp(-pow(sigma,2)*ksquare/2.);
                        
                        // for short-range terms
                        if ( !(m == 0 && n == 0 && l == 0 && i == id) )
                            Vs += (uc[i].C - uc[i].ne) / dR * erfc(dR/sqrt(2)/sigma);
                    }
                }
        uc[id].V = Vs + Vl - (uc[id].C - uc[id].ne)*sqrt(2/M_PI)/sigma;
    }
}

static int l_uc_al2bv(lua_State* L)
{
	UnitCell* uc = lua_tounitcell(L, 1);
	if(!uc) return 0;
	
	double v[6];
	for(int i=0; i<6; i++)
	{
		v[i] = lua_tonumber(L, i+2);
		if(v[i] == 0)
			return luaL_error(L, "UnitCell:anglesLengthsToBasisVectors requires 6 numbers");
	}
	
	uc->anglesLengthsToBasisVectors(
		v[0], v[1], v[2],
		v[3], v[4], v[5]);
	return 0;
}

void PawleyFunction::setPeakPositions(std::string centreName, double zeroShift,
                                      const UnitCell &cell) const {
  for (size_t i = 0; i < m_hkls.size(); ++i) {
    double centre = getTransformedCenter(cell.d(m_hkls[i]));

    m_peakProfileComposite->getFunction(i)
        ->setParameter(centreName, centre + zeroShift);
  }
}

void Crystal::toggleUnitCell()
{
  if (m_molecule->unitCell()) {
    m_molecule->setUnitCell(NULL);
    m_molecule->emitChanged(Molecule::UnitCell | Molecule::Removed);
  }
  else {
    UnitCell *cell = new UnitCell;
    cell->setCellParameters(static_cast<Real>(3.0),
                            static_cast<Real>(3.0),
                            static_cast<Real>(3.0),
                            static_cast<Real>(90.0) * DEG_TO_RAD,
                            static_cast<Real>(90.0) * DEG_TO_RAD,
                            static_cast<Real>(90.0) * DEG_TO_RAD);
    m_molecule->setUnitCell(cell);
    m_molecule->emitChanged(Molecule::UnitCell | Molecule::Added);
    editUnitCell();
  }
}

static int l_uc_copy(lua_State* L)
{
	UnitCell* uc = lua_tounitcell(L, 1);
	if(!uc) return 0;
	
	UnitCell* uc2 = new UnitCell();//uc->A(), uc->B(), uc->C());
	
	uc2->r2g = uc->r2g;
	uc2->g2r = uc->g2r;
	
	uc2->setA(uc->A());
	uc2->setB(uc->B());
	uc2->setC(uc->C());
	
	for(unsigned int i=0; i<uc->atoms.size(); i++)
	{
		uc2->addAtomGlobalCoorinates(new Atom(*uc->atoms[i]));
	}
	lua_pushunitcell(L, uc2);
	return 1;
}

// Wrap a vector in an Orthorombic UnitCell
static Vector3D wrap_orthorombic(const UnitCell& cell, const Vector3D& vect) {
    Vector3D res;
    res[0] = static_cast<float>(vect[0] - round(vect[0]/cell.a())*cell.a());
    res[1] = static_cast<float>(vect[1] - round(vect[1]/cell.b())*cell.b());
    res[2] = static_cast<float>(vect[2] - round(vect[2]/cell.c())*cell.c());
    return res;
}

// Wrap a vector in an Orthorombic UnitCell
static Vector3D wrap_triclinic(const UnitCell& cell, const Vector3D& vect) {
    Vector3D res = vect;

    auto mat = cell.matricial();
    for (size_t i=2 ; i != static_cast<size_t>(-1) ; i--) {
        while (fabs(res[i]) > mat[i][i]/2) {
            if (res[i] < 0) {
                res[0] += static_cast<float>(mat[i][0]);
                res[1] += static_cast<float>(mat[i][1]);
                res[2] += static_cast<float>(mat[i][2]);
            } else {
                res[0] -= static_cast<float>(mat[i][0]);
                res[1] -= static_cast<float>(mat[i][1]);
                res[2] -= static_cast<float>(mat[i][2]);
            }
        }
    }
    return res;
}