C++ Direction::cartesian方法代码示例

Position Trust2Geometry::generatePosition() const
{
    double X = _random->uniform();
    if (X<_M0)
    {
        double dd0 = 0.0;
        double x, y, z;
        while (dd0<_R0*_R0)
        {
            x = (2.0*_random->uniform()-1.0) * _L0;
            y = (2.0*_random->uniform()-1.0) * _L0;
            z = (2.0*_random->uniform()-1.0) * _L0;
            dd0 = (x+_L0)*(x+_L0) + (y+_L0)*(y+_L0) + (z+_L0)*(z+_L0);
        }
        return Position(x,y,z);
    }
    else if (X<_M0+_M1)
    {
        Direction bfk = _random->direction();
        double kx, ky, kz;
        bfk.cartesian(kx,ky,kz);
        double r = _R1*pow(_random->uniform(),1.0/3.0);
        double x1, y1, z1;
        _bfr1.cartesian(x1,y1,z1);
        double x = x1 + r*kx;
        double y = y1 + r*ky;
        double z = z1 + r*kz;
        return Position(x,y,z);
    }
    else
    {
        Direction bfk = _random->direction();
        double kx, ky, kz;
        bfk.cartesian(kx,ky,kz);
        double r = _R2*pow(_random->uniform(),1.0/3.0);
        double x2, y2, z2;
        _bfr2.cartesian(x2,y2,z2);
        double x = x2 + r*kx;
        double y = y2 + r*ky;
        double z = z2 + r*kz;
        return Position(x,y,z);
    }
}

DustGridPath ParticleTreeDustGridStructure::path(Position bfr, Direction bfk) const
{
    // Initialize the path
    DustGridPath path(bfr, bfk, 3000);

    // If the photon package starts outside the dust grid, move it into the first grid cell that it will pass
    Position r = path.moveinside(extent(), _eps);

    // Get the node containing the current location;
    // if the position is not inside the grid, return an empty path
    const TreeNode* node = root()->whichnode(r);
    if (!node) return path.clear();

    // Start the loop over nodes/path segments until we leave the grid.
    // Use a different code segment depending on the search method.
    double x,y,z;
    r.cartesian(x,y,z);
    double kx,ky,kz;
    bfk.cartesian(kx,ky,kz);

    while (node)
    {
        // calculate the distances to the planes of the node walls that are forward of this point
        double xnext = (kx<0.0) ? node->xmin() : node->xmax();
        double ynext = (ky<0.0) ? node->ymin() : node->ymax();
        double znext = (kz<0.0) ? node->zmin() : node->zmax();
        double dsx = (xnext-x)/kx;
        double dsy = (ynext-y)/ky;
        double dsz = (znext-z)/kz;

        // find the distance to the nearest wall (avoiding the wall containing the entry point)
        // and add the corresponding path segment (unless there is no exit point due to rounding errors)
        double ds = DBL_MAX;  // very large, but not infinity (so that infinite dsi values are discarded)
        if (dsx>0 && dsx<ds) ds = dsx;
        if (dsy>0 && dsy<ds) ds = dsy;
        if (dsz>0 && dsz<ds) ds = dsz;
        if (ds<DBL_MAX)
            path.addsegment(cellnumber(node), ds);
        else
            ds = 0;

        // advance the current point
        x += (ds+_eps)*kx;
        y += (ds+_eps)*ky;
        z += (ds+_eps)*kz;

        // always search from the root node down
        node = root()->whichnode(Vec(x,y,z));
    }

    return path;
}

DustGridPath SpheDustGridStructure::path(Position bfr, Direction bfk) const
{
    // Determination of the initial position and direction of the path,
    // and calculation of some initial values

    DustGridPath path(bfr, bfk, 2*_Nr + 2);
    double x,y,z;
    bfr.cartesian(x,y,z);
    double kx,ky,kz;
    bfk.cartesian(kx,ky,kz);

    // Move the photon package to the first grid cell that it will pass.
    // If it does not pass any grid cell, return an empty path.

    double r = bfr.radius();
    double q = x*kx + y*ky + z*kz;
    double p = sqrt((r-q)*(r+q));
    if (r>_rmax)
    {
        if (q>0.0 || p>_rmax) return path.clear();
        else
        {
            r = _rmax - 1e-8*(_rv[_Nr]-_rv[_Nr-1]);
            double qmax = sqrt((_rmax-p)*(_rmax+p));
            double ds = (qmax-q);
            path.addsegment(-1,ds);
            q = qmax;
        }
    }

    // Determination of the initial grid cell

    int i = NR::locate_clip(_rv, r);

    // And here we go...

    double rN, qN;

    // Inward movement (not always...)

    if (q<0.0)
    {
        int imin = NR::locate_clip(_rv,p);
        rN = _rv[i];
        qN = -sqrt((rN-p)*(rN+p));
        while (i>imin)
        {
            int m = i;
            double ds = qN-q;
            path.addsegment(m, ds);
            i--;
            q = qN;
            rN = _rv[i];
            qN = -sqrt((rN-p)*(rN+p));
        }
    }

    // Outward movement

    rN = _rv[i+1];
    qN = sqrt((rN-p)*(rN+p));
    while (true)
    {
        int m = i;
        double ds = qN-q;
        path.addsegment(m, ds);
        i++;
        if (i>=_Nr-1) return path;
        else
        {
            q = qN;
            rN = _rv[i+1];
            qN = sqrt((rN-p)*(rN+p));
        }
    }
}