C++ MatrixPtr::multiplyVector方法代码示例

double Miller::expectedRadius(double spotSize, double mosaicity, vec *hkl)
{
    vec usedHKL;
    
    if (hkl == NULL)
    {
        MatrixPtr newMatrix = MatrixPtr();
        rotateMatrixHKL(latestHRot, latestKRot, 0, matrix, &newMatrix);
        
        usedHKL = new_vector(h, k, l);
        hkl = &usedHKL;
        
        newMatrix->multiplyVector(hkl);
    }
    
    
    spotSize = fabs(spotSize);
    double radMos = fabs(mosaicity) * M_PI / 180;
    
    double distanceFromOrigin = length_of_vector(*hkl);
    
    double spotSizeIncrease = fabs(radMos * distanceFromOrigin);
    
    double radius = (spotSize + spotSizeIncrease);
    
    return radius;
}

/**
 Generates twinned/untwinned reflection IDs for easy searching.
 */
void Reflection::generateReflectionIds()
{
    if (millerCount() == 0)
    {
        std::cout << "Warning! Miller count is 0" << std::endl;
    }
    
    int h = miller(0)->getH();
    int k = miller(0)->getK();
    int l = miller(0)->getL();
    
    vec hkl = new_vector(h, k, l);
    
    for (int i = 0; i < ambiguityCount(); i++)
    {
        MatrixPtr ambiguityMat = matrixForAmbiguity(i);
        
        int asuH, asuK, asuL;
        
        ccp4spg_put_in_asu(spaceGroup, hkl.h, hkl.k, hkl.l, &asuH, &asuK, &asuL);
        vec hklAsu = new_vector(asuH, asuK, asuL);
        
        ambiguityMat->multiplyVector(&hklAsu);
       
        int newId = reflectionIdForMiller(hklAsu);
        
        reflectionIds.push_back(newId);
    }
}

vec Miller::getTransformedHKL(MatrixPtr matrix)
{
    vec hkl = new_vector(h, k, l);
    
    matrix->multiplyVector(&hkl);
    
    return hkl;
}

void Miller::recalculatePartiality(MatrixPtr rotatedMatrix, double mosaicity,
                                   double spotSize, double wavelength, double bandwidth, double exponent, bool binary)
{
    if (model == PartialityModelFixed)
    {
        return;
    }
    
    if (model == PartialityModelBinary)
    {
  //      binary = true;
    }
    
    vec hkl = new_vector(h, k, l);
    
    rotatedMatrix->multiplyVector(&hkl);
    
    this->wavelength = getEwaldSphereNoMatrix(hkl);
    
    double tempPartiality = partialityForHKL(hkl, mosaicity,
                                            spotSize, wavelength, bandwidth, exponent, binary);
    
    if (binary && tempPartiality == 1.0)
        return;
    
    double normPartiality = 1;
    
    if (normalised && tempPartiality > 0)
    {
        normPartiality = calculateNormPartiality(rotatedMatrix, mosaicity, spotSize, wavelength, bandwidth, exponent);
    }
    else
    {
        partiality = 0;
        return;
    }
    
    partiality = tempPartiality / normPartiality;
    
    if (partiality > 1.2)
        partiality = 0;
    if (partiality > 1.0)
        partiality = 1;
    
    if (partiality > 1 && std::isfinite(partiality))
    {
        std::cout << "Partiality: " << partiality << std::endl;
    }
        
    if ((!std::isfinite(partiality)) || (partiality != partiality))
    {
        partiality = 0;
    }
}

void Shoebox::complexShoebox(double wavelength, double bandwidth, double radius)
{
    MillerPtr tempMiller = getMiller();
    
    if (!tempMiller)
    {
        return;
    }
    
    MatrixPtr mat;
    mat = tempMiller->getMatrix();
    
    vec hkl = new_vector(tempMiller->getH(), tempMiller->getK(), tempMiller->getL());
    mat->multiplyVector(&hkl);
    
    double hk = sqrt(hkl.h * hkl.h + hkl.k * hkl.k);
    
    double maxHK = hk + radius;
    double minHK = hk - radius;
    
    double minWave = wavelength - bandwidth * bandwidthMultiplier;
    double maxWave = wavelength + bandwidth * bandwidthMultiplier;
    
    double max_mm = (maxHK * detectorDistance / (1 / maxWave + hkl.l));
    double min_mm = (minHK * detectorDistance / (1 / minWave + hkl.l));
    
    double width_mm = (radius * detectorDistance / (1 / wavelength));
    
    double mmDiff = max_mm - min_mm;
    
    double widthLeak = 1 * pixelLeak;

    double pixelHeight = mmDiff / mmPerPixel + pixelLeak;
    double pixelWidth = 2 * width_mm / mmPerPixel + widthLeak;
    
    double angle = cartesian_to_angle(hkl.h, hkl.k);
    angle += M_PI / 2;
    
    Box smallBox;
    Box newShoebox;
    compressBigShoebox(pixelWidth, pixelHeight, angle, smallBox);
    chopBox(smallBox);
    centreShoebox(smallBox);
    putBoxOnPaddedBackground(smallBox, newShoebox);

    shoebox = newShoebox;
    printShoebox();
}

double Miller::calculateNormPartiality(MatrixPtr rotatedMatrix, double mosaicity,
                                      double spotSize, double wavelength, double bandwidth, double exponent)
{
    vec hkl = new_vector(h, k, l);
    
    rotatedMatrix->multiplyVector(&hkl);
    
    double d = length_of_vector(hkl);
    
    double newH = 0;
    double newK = sqrt((4 * pow(d, 2) - pow(d, 4) * pow(wavelength, 2)) / 4);
    double newL = 0 - pow(d, 2) * wavelength / 2;
    
    vec newHKL = new_vector(newH, newK, newL);
    
    double normPartiality = partialityForHKL(newHKL, mosaicity,
                                            spotSize, wavelength, bandwidth, exponent);
    
    return normPartiality;
}

bool Miller::crossesBeamRoughly(MatrixPtr rotatedMatrix, double mosaicity,
                         double spotSize, double wavelength, double bandwidth)
{
    vec hkl = new_vector(h, k, l);
    
    rotatedMatrix->multiplyVector(&hkl);
    
    double inwards_bandwidth = 0; double outwards_bandwidth = 0;
    
    limitingEwaldWavelengths(hkl, mosaicity, spotSize, wavelength, &outwards_bandwidth, &inwards_bandwidth);
    
    double limit = 2 * bandwidth;
    double inwardsLimit = wavelength + limit;
    double outwardsLimit = wavelength - limit;
    
    if (outwards_bandwidth > inwardsLimit || inwards_bandwidth < outwardsLimit)
    {
        partiality = 0;
        return false;
    }
    
    partiality = 1;
    return true;
}

void Miller::positionOnDetector(MatrixPtr transformedMatrix, int *x,
                                int *y)
{
    double x_coord = 0;
    double y_coord = 0;
    
    vec hkl = new_vector(h, k, l);
    transformedMatrix->multiplyVector(&hkl);
    
    std::pair<double, double> coord = getImage()->reciprocalCoordinatesToPixels(hkl);
    x_coord = coord.first;
    y_coord = coord.second;
    
    bool even = shoebox->isEven();

    int intLastX = int(x_coord);
    int intLastY = int(y_coord);
    
    if (even)
    {
        intLastX = round(x_coord);
        intLastY = round(y_coord);
    }
    
    lastX = x_coord;
    lastY = y_coord;
    
    if (!Panel::shouldUsePanelInfo())
    {
        int search = indexer->getSearchSize();
        getImage()->focusOnAverageMax(&intLastX, &intLastY, search, peakSize, even);
        
        shift = std::make_pair(intLastX + 0.5 - x_coord, intLastY + 0.5 - y_coord);
        
        *x = intLastX;
        *y = intLastY;
    }
    else
    {
        int search = indexer->getSearchSize();
        Coord bestShift = Panel::shiftForMiller(this);
        
        if (bestShift.first != FLT_MAX)
        {
            double shiftedX = lastX + bestShift.first;
            double shiftedY = lastY + bestShift.second;
            
            int xInt = shiftedX;
            int yInt = shiftedY;
            
            getImage()->focusOnAverageMax(&xInt, &yInt, search, peakSize, even);
            
            shift = std::make_pair(xInt + 0.5 - x_coord, yInt + 0.5 - y_coord);
            
            *x = xInt;
            *y = yInt;
        }
        else
        {
            int search = indexer->getSearchSize();
            
            getImage()->focusOnAverageMax(&intLastX, &intLastY, search, peakSize, even);
            
            *x = intLastX;
            *y = intLastY;
        }
    }
}