packages/MC/geometry/RTK_Cell.cc
//---------------------------------------------------------------------------//
/*!
 * \brief Find the region for a point.
 */
\brief Find the region for a point.
\ brief Find the region for a point .
int RTK_Cell::region(double x,
                     double y) const
{
    using def::X; using def::Y;
    REQUIRE(d_num_regions > 0);
    REQUIRE(x > d_extent[0][LO] || soft_equiv(x, d_extent[0][LO]));
    REQUIRE(x < d_extent[0][HI] || soft_equiv(x, d_extent[0][HI]));
    REQUIRE(y > d_extent[1][LO] || soft_equiv(y, d_extent[1][LO]));
    REQUIRE(y < d_extent[1][HI] || soft_equiv(y, d_extent[1][HI]));
    if (!d_num_shells)
    {
        REQUIRE(d_num_regions == 1);
        if (d_vessel)
        {
            double gx = this->l2g(x, X);
            double gy = this->l2g(y, Y);
            double R2 = gx*gx + gy*gy;
            if (d_inner && d_outer)
            {
                if (R2 < d_R0*d_R0)
                    return Geo_State_t::MODERATOR;
                else if (R2 > d_R1*d_R1)
                    return Geo_State_t::MODERATOR;
                else
                    return Geo_State_t::VESSEL;
            }
            else if (d_inner)
            {
                if (R2 > d_R0*d_R0)
                    return Geo_State_t::VESSEL;
                else
                    return  Geo_State_t::MODERATOR;
            }
            else if (d_outer)
            {
                if (R2 < d_R1*d_R1)
                    return Geo_State_t::VESSEL;
                else
                    return Geo_State_t::MODERATOR;
            }
            else
            {
                VALIDATE(false, "Vessel without bisection.");
            }
        }
        else
        {
            return Geo_State_t::MODERATOR;
        }
    }
    double rp2 = x*x + y*y;
    for (int n = 0; n < d_num_shells; ++n)
    {
        if (rp2 <= d_r[n] * d_r[n]) return n;
    }
    ENSURE(rp2 > d_r[d_num_shells-1]*d_r[d_num_shells-1]);
    return d_num_regions - 1;
}
int RTK_Cell :: region ( double x , double y ) const { using def :: X ; using def :: Y ; REQUIRE ( d_num_regions > 0 ) ; REQUIRE ( x > d_extent [ 0 ] [ LO ] || soft_equiv ( x , d_extent [ 0 ] [ LO ] ) ) ; REQUIRE ( x < d_extent [ 0 ] [ HI ] || soft_equiv ( x , d_extent [ 0 ] [ HI ] ) ) ; REQUIRE ( y > d_extent [ 1 ] [ LO ] || soft_equiv ( y , d_extent [ 1 ] [ LO ] ) ) ; REQUIRE ( y < d_extent [ 1 ] [ HI ] || soft_equiv ( y , d_extent [ 1 ] [ HI ] ) ) ; if ( ! d_num_shells ) { REQUIRE ( d_num_regions == 1 ) ; if ( d_vessel ) { double gx = this -> l2g ( x , X ) ; double gy = this -> l2g ( y , Y ) ; double R2 = gx * gx + gy * gy ; if ( d_inner && d_outer ) { if ( R2 < d_R0 * d_R0 ) return Geo_State_t :: MODERATOR ; else if ( R2 > d_R1 * d_R1 ) return Geo_State_t :: MODERATOR ; else return Geo_State_t :: VESSEL ; } else if ( d_inner ) { if ( R2 > d_R0 * d_R0 ) return Geo_State_t :: VESSEL ; else return Geo_State_t :: MODERATOR ; } else if ( d_outer ) { if ( R2 < d_R1 * d_R1 ) return Geo_State_t :: VESSEL ; else return Geo_State_t :: MODERATOR ; } else { VALIDATE ( false , " " ) ; } } else { return Geo_State_t :: MODERATOR ; } } double rp2 = x * x + y * y ; for ( int n = 0 ; n < d_num_shells ; ++ n ) { if ( rp2 <= d_r [ n ] * d_r [ n ] ) return n ; } ENSURE ( rp2 > d_r [ d_num_shells - 1 ] * d_r [ d_num_shells - 1 ] ) ; return d_num_regions - 1 ; }
\brief Find the region for a point.
\ brief Find the region for a point .
// if this is an empty box, return the moderator // check to see if it is in a vessel region // get coordinates of point in vessel coordinate system // calculate rp // loop through shells to find the containing shell // point is in the moderator region