Vec3r rotateVector(const Matrix44r& mat, const Vec3r& v) {
    Vec3r res;
    for (unsigned i = 0; i < 3; i++) {
      res[i] = 0;
      for (unsigned j = 0; j < 3; j++)
	res[i] += mat(i, j) * v[j];
    }
    res.normalize();
    return res;
  }

bool Line::intersect(const SphereVolume &sphere,
                           Real         &enter, 
                           Real         &exit  ) const
{
    Vec3r v;
    Pnt3r center;

    sphere.getCenter(center);

    Real radius;
    Real h;
    Real b;
    Real d;
    Real t1;
    Real t2;

    radius = sphere.getRadius();

    v = center - _pos;

    h = (v.dot(v))-(radius * radius);
    b = (v.dot(_dir));

    if(h >= 0.f && b <= 0.f)
        return false;

    d = b * b - h;

    if(d < 0.f)
        return false;

    d  = osgSqrt(d);
    t1 = b - d;

//    if (t1 > 1)
//        return false;

    t2 = b + d;

    if( t1 < TypeTraits<Real>::getDefaultEps() )
    {
        if( t2 < TypeTraits<Real>::getDefaultEps() /*|| t2 > 1*/)
        {
            return false;
        }
    }

    enter = t1;
    exit  = t2;

    return true;
}

void FEdgeXDetector::preProcessFace(WXFace *iFace){
    Vec3r firstPoint = iFace->GetVertex(0)->GetVertex();
    Vec3r N = iFace->GetNormal();

    // Compute the dot product between V (=_Viewpoint - firstPoint) and N:
    Vec3r V(_Viewpoint - firstPoint);
    N.normalize();
    V.normalize();
    iFace->SetDotP(N * V);

    // compute the distance between the face center and the viewpoint:
    Vec3r dist_vec(iFace->center() - _Viewpoint);
    iFace->SetZ(dist_vec.norm());
}

inline static real angle(WOEdge *h)
{
	const Vec3r& n1 = h->GetbFace()->GetNormal();
	const Vec3r& n2 = h->GetaFace()->GetNormal();
	const Vec3r v = h->GetVec();
	real sine = (n1 ^ n2) * v / v.norm();
	if (sine >= 1.0) {
		return M_PI / 2.0;
	}
	if (sine <= -1.0) {
		return -M_PI / 2.0;
	}
	return ::asin(sine);
}

void FEdgeXDetector::ProcessSilhouetteFace(WXFace *iFace, bool meshSilhouettes)
{

    real NdotVepsilonHack = 0;// 0.05; //0.1; //0.01; //0.1;

    // SILHOUETTE LAYER
    // Compute the dot products between View direction and N at each vertex
    // of the face:
    Vec3r point;
    int closestPointId = 0;
    real dist, minDist = FLT_MAX;
    int numVertices = iFace->numberOfVertices();
    WXFaceLayer * faceLayer = new WXFaceLayer(iFace, Nature::SILHOUETTE, true);

    Vec3r normal;
    if(meshSilhouettes){
        // Use per face normal
        normal = (iFace->GetVertex(2)->GetVertex() - iFace->GetVertex(0)->GetVertex()) ^ (iFace->GetVertex(1)->GetVertex() - iFace->GetVertex(0)->GetVertex());
        normal.normalize();
    }

    for(int i=0; i<numVertices; i++){
        point = iFace->GetVertex(i)->GetVertex();
        if(!meshSilhouettes){
            // Use per vertex normal
            normal = iFace->GetVertexNormal(i);
            normal.normalize();
        }
        Vec3r V(_Viewpoint - point);
        V.normalize();
        real d = normal * V + NdotVepsilonHack;
        faceLayer->PushDotP(d);
        // Find the point the closest to the viewpoint
        Vec3r dist_vec(point - _Viewpoint);
        dist = dist_vec.norm();
        if(dist < minDist) {
            minDist = dist;
            closestPointId = i;
        }

        // store ndotv at the vertex for use in the region-based visibility
        //    assert(dynamic_cast<WXVertex*>(iFace->GetVertex(i))!=NULL);
        ((WXVertex*)iFace->GetVertex(i))->setNdotV(d);
    }
    // Set the closest point id:
    faceLayer->SetClosestPointIndex(closestPointId);
    // Add this layer to the face:
    iFace->AddSmoothLayer(faceLayer);
}

static real angle_from_cotan(WVertex *vo, WVertex *v1, WVertex *v2)
{
	/* cf. Appendix B and the caption of Table 1 from [Meyer et al 2002] */
	real udotv, denom;

	Vec3r u (v1->GetVertex() - vo->GetVertex());
	Vec3r v(v2->GetVertex() - vo->GetVertex());

	udotv = u * v;
	denom = sqrt(u.squareNorm() * v.squareNorm() - udotv * udotv);

	/* Note: I assume this is what they mean by using atan2(). -Ray Jones */

	/* tan = denom/udotv = y/x (see man page for atan2) */
	return (fabs(atan2(denom, udotv)));
}

void MonaghanKernel::monaghanGradient( const Vec3r& r, Vec3r& gradient )
{
    HReal dist = r.length();
    HReal q = dist*m_invH;
    gradient.fill(0.0);
    if( q >= 0 && q < 1 )
    {
        HReal scalar = -3.0f*(2-q)*(2-q);
        scalar += 12.0f*(1-q)*(1-q);
        gradient = (m_g*m_invH*scalar/dist)*r;
    }
    else if ( q >=1 && q < 2 )
    {
        HReal scalar = -3.0f*(2-q)*(2-q);
        gradient = (m_g*scalar*m_invH/dist)*r;
    }
}

void Grid::castInfiniteRay(const Vec3r& orig,
			   const Vec3r& dir,
			   OccludersSet& occluders,
			   unsigned timestamp) {
  Vec3r end = Vec3r(orig + FLT_MAX * dir / dir.norm());
  bool inter = initInfiniteRay(orig, dir, timestamp);
  if(!inter)
      return;
  allOccludersGridVisitor visitor(occluders);
  castRayInternal(visitor);
}

// precondition1: P is inside the sphere
// precondition2: P,V points to the outside of the sphere (i.e. OP.V > 0)
static bool sphere_clip_vector(const Vec3r& O, real r, const Vec3r& P, Vec3r& V)
{
	Vec3r W = P - O;
	real a = V.squareNorm();
	real b = 2.0 * V * W;
	real c = W.squareNorm() - r * r;
	real delta = b * b - 4 * a * c;
	if (delta < 0) {
		// Should not happen, but happens sometimes (numerical precision)
		return true;
	}
	real t = - b + ::sqrt(delta) / (2.0 * a);
	if (t < 0.0) {
		// Should not happen, but happens sometimes (numerical precision)
		return true;
	}
	if (t >= 1.0) {
		// Inside the sphere
		return false;
	}

	V[0] = (t * V.x());
	V[1] = (t * V.y());
	V[2] = (t * V.z());

	return true;
}

void FEdgeXDetector::preProcessFace(WXFace *iFace)
{
	Vec3r firstPoint = iFace->GetVertex(0)->GetVertex();
	Vec3r N = iFace->GetNormal();

	// Compute the dot product between V (=_Viewpoint - firstPoint) and N:
	Vec3r V;
	if (_orthographicProjection) {
		V = Vec3r(0.0, 0.0, _Viewpoint.z() - firstPoint.z());
	}
	else {
		V = Vec3r(_Viewpoint - firstPoint);
	}
	N.normalize();
	V.normalize();
	iFace->setDotP(N * V);

	// compute the distance between the face center and the viewpoint:
	if (_orthographicProjection) {
		iFace->setZ(iFace->center().z() - _Viewpoint.z());
	}
	else {
		Vec3r dist_vec(iFace->center() - _Viewpoint);
		iFace->setZ(dist_vec.norm());
	}
}

    void VelocityMotor::internal_update()
    {
        // check if we have a solid
        if (mData.solid == NULL || isEnabled() == false) return ;

        Vec3r targetVelocity = mData.velocity;
        Solid * solid = mData.solid;

        Vec3r currentAchievedVelocity = solid->getGlobalLinearVel();

        if (doesGravityAffectSolid())
        {
            Vec3r gravity = mSimulator->getGravity();
            if (gravity.length() > 0)
            {
                Vec3r gravity_velocity = project(gravity, currentAchievedVelocity);
                currentAchievedVelocity -= gravity_velocity;
            }
        }

        Vec3r deltaVelocity = targetVelocity - currentAchievedVelocity;

        Vec3r forceVector = deltaVelocity / mSimulator->getStepSize() * solid->getMass();
        if (!doesGravityAffectSolid())
            forceVector -= mSimulator->getGravity() * solid->getMass();

        if (forceVector.length() > getMaximumForce())
        {
            forceVector.normalize();
            forceVector *= getMaximumForce();
        }

        Force controllingForce;
        controllingForce.duration = 0;
        controllingForce.singleStep = true;
        controllingForce.type = GLOBAL_FORCE;
        controllingForce.vec = forceVector;

        solid->addForce(controllingForce);
    }

void ThirdPersonCamera::LookAt(const Vec3r &eye, const Vec3r &target, const Vec3r &up)
{
    m_eye = eye;
    m_target = target;
    m_targetYAxis = up;

    m_zAxis = eye - target;
    m_zAxis.Normalise();

    m_viewDir = m_zAxis * -1;

    m_xAxis = Vec3r::CrossProduct(up, m_zAxis);
    m_xAxis.Normalise();

    m_yAxis = Vec3r::CrossProduct(m_zAxis, m_xAxis);
    m_yAxis.Normalise();
   // m_xAxis.Normalise();

    m_viewMatrix.Access(0,0) = m_xAxis.X;
    m_viewMatrix.Access(1,0) = m_xAxis.Y;
    m_viewMatrix.Access(2,0) = m_xAxis.Z;
	m_viewMatrix.Access(3,0) = -Vec3r::DotProduct(m_xAxis, eye);

    m_viewMatrix.Access(0,1) = m_yAxis.X;
    m_viewMatrix.Access(1,1) = m_yAxis.Y;
    m_viewMatrix.Access(2,1) = m_yAxis.Z;
	m_viewMatrix.Access(3,1) = -Vec3r::DotProduct(m_yAxis, eye);
				
    m_viewMatrix.Access(0,2) = m_zAxis.X;
    m_viewMatrix.Access(1,2) = m_zAxis.Y;
    m_viewMatrix.Access(2,2) = m_zAxis.Z;    
	m_viewMatrix.Access(3,2) = -Vec3r::DotProduct(m_zAxis, eye);

	m_orientation.FromMatrix(m_viewMatrix.m_matrix);

    Vec3r offset = m_target - m_eye;

	m_offsetDistance = offset.Length();
}

HReal MonaghanKernel::monaghanValue( const Vec3r & r )
{
    HReal value = 0.0;
    HReal q = r.length()*m_invH;
    if( q >= 0 && q < 1 )
    {
        value = m_v*( (2-q)*(2-q)*(2-q) - 4.0f*(1-q)*(1-q)*(1-q));
    }
    else if ( q >=1 && q < 2 )
    {
        value = m_v*( (2-q)*(2-q)*(2-q) );
    }
    else
    {
        value = 0.0f;
    }
    return value;
}

void FEdgeXDetector::ProcessSilhouetteFace(WXFace *iFace)
{
	// SILHOUETTE LAYER
	Vec3r normal;
	// Compute the dot products between View direction and N at each vertex of the face:
	Vec3r point;
	int closestPointId = 0;
	real dist, minDist = FLT_MAX;
	int numVertices = iFace->numberOfVertices();
	WXFaceLayer *faceLayer = new WXFaceLayer(iFace, Nature::SILHOUETTE, true);
	for (int i = 0; i < numVertices; i++) {
		point = iFace->GetVertex(i)->GetVertex();
		normal = iFace->GetVertexNormal(i);
		normal.normalize();
		Vec3r V;
		if (_orthographicProjection) {
			V = Vec3r(0.0, 0.0, _Viewpoint.z() - point.z());
		}
		else {
			V = Vec3r(_Viewpoint - point);
		}
		V.normalize();
		real d = normal * V;
		faceLayer->PushDotP(d);
		// Find the point the closest to the viewpoint
		if (_orthographicProjection) {
			dist = point.z() - _Viewpoint.z();
		}
		else {
			Vec3r dist_vec(point - _Viewpoint);
			dist = dist_vec.norm();
		}
		if (dist < minDist) {
			minDist = dist;
			closestPointId = i;
		}
	}
	// Set the closest point id:
	faceLayer->setClosestPointIndex(closestPointId);
	// Add this layer to the face:
	iFace->AddSmoothLayer(faceLayer);
}

real CurvePoint::curvature2d_as_angle() const
{
#if 0
	Vec3r edgeA = (_FEdges[0])->orientation2d();
	Vec3r edgeB = (_FEdges[1])->orientation2d();
	Vec2d N1(-edgeA.y(), edgeA.x());
	N1.normalize();
	Vec2d N2(-edgeB.y(), edgeB.x());
	N2.normalize();
	return acos((N1 * N2));
#endif
	if (__A == 0)
		return __B->curvature2d_as_angle();
	if (__B == 0)
		return __A->curvature2d_as_angle();
	return ((1 - _t2d) * __A->curvature2d_as_angle() + _t2d * __B->curvature2d_as_angle());
}

    bool intersectRayBBox(const Vec3r& orig, const Vec3r& dir,   // ray origin and direction
        const Vec3r& boxMin, const Vec3r& boxMax, // the bbox
        real t0, real t1,
        real& tmin, real& tmax,                  // I0=orig+tmin*dir is the first intersection, I1=orig+tmax*dir is the second intersection
        real epsilon){

            float tymin, tymax, tzmin, tzmax;
            Vec3r inv_direction(1.0/dir[0], 1.0/dir[1], 1.0/dir[2]);
            int sign[3];
            sign[0] = (inv_direction.x() < 0);
            sign[1] = (inv_direction.y() < 0);
            sign[2] = (inv_direction.z() < 0);

            Vec3r bounds[2];
            bounds[0] = boxMin;
            bounds[1] = boxMax;

            tmin = (bounds[sign[0]].x() - orig.x()) * inv_direction.x();
            tmax = (bounds[1-sign[0]].x() - orig.x()) * inv_direction.x();
            tymin = (bounds[sign[1]].y() - orig.y()) * inv_direction.y();
            tymax = (bounds[1-sign[1]].y() - orig.y()) * inv_direction.y();
            if ( (tmin > tymax) || (tymin > tmax) )
                return false;
            if (tymin > tmin)
                tmin = tymin;
            if (tymax < tmax)
                tmax = tymax;
            tzmin = (bounds[sign[2]].z() - orig.z()) * inv_direction.z();
            tzmax = (bounds[1-sign[2]].z() - orig.z()) * inv_direction.z();
            if ( (tmin > tzmax) || (tzmin > tmax) )
                return false;
            if (tzmin > tmin)
                tmin = tzmin;
            if (tzmax < tmax)
                tmax = tzmax;
            return ( (tmin < t1) && (tmax > t0) );
        }

/**
 * @brief 円柱を設定
 * @param [in]     R     　　 Controlクラスのポインタ
 * @param [in]     pos_x     中心座標(x,y) （無次元）
 * @param [in]     pos_y     中心座標(x,y) （無次元）
 * @param [in]     radius    半径 （無次元）
 * @param [in]     len_z     円柱のZ方向の長さ （無次元）
 * @param [in]     mid_solid 固体媒質ID
 * @param [out]    cut       カット情報
 * @param [out]    bid       境界ID
 */
void IP_Cylinder::setCircle(Control* R,
                            const REAL_TYPE pos_x,
                            const REAL_TYPE pos_y,
                            const REAL_TYPE radius,
                            const REAL_TYPE len_z,
                            const int mid_solid,
                            long long* cut,
                            int* bid)
{
  Vec3r pch(pitch);      ///< 無次元格子幅
  Vec3r org(origin);     ///< ノードローカル基点座標の無次元値

  
  REAL_TYPE dx = pitch[0];
  REAL_TYPE dy = pitch[1];
  REAL_TYPE dz = pitch[2];
  
  REAL_TYPE cx = pos_x;  ///< 円柱の中心座標（無次元）
  REAL_TYPE cy = pos_y;  ///< 円柱の中心座標（無次元）
  REAL_TYPE rs = radius; ///< 円柱の半径（無次元）
  
  int mid_s = mid_solid;
  
  // 球のbbox
  Vec3r box_min;  ///< Bounding boxの最小値
  Vec3r box_max;  ///< Bounding boxの最大値
  Vec3i box_st;   ///< Bounding boxの始点インデクス
  Vec3i box_ed;   ///< Bounding boxの終点インデクス
  box_min.assign( -(REAL_TYPE)rs+cx, -(REAL_TYPE)rs+cy, 0.0 ); // Z方向はダミー
  box_max.assign(  (REAL_TYPE)rs+cx,  (REAL_TYPE)rs+cy, 0.0 );
  box_st = find_index(box_min, org, pch);
  box_ed = find_index(box_max, org, pch);
  
  int ix = size[0];
  int jx = size[1];
  int kx = size[2];
  int gd = guide;
  
  
  // ローカルな無次元基点座標
  REAL_TYPE ox = origin[0];
  REAL_TYPE oy = origin[1];
  REAL_TYPE oz = origin[2];
  
  // グローバルな無次元座標
  REAL_TYPE zs = oz;
  REAL_TYPE ze;
  
  if ( mode == dim_2d )
  {
    ze = oz + region[2] + dz; // Z方向には全セルを対象, dzは安全係数
  }
  else
  {
    ze = oz + len_z; // Z-面からの距離
  }
  
  
  // カット情報
  Vec3r p[5];
  Vec3r base;  // セルのシフト量
  Vec3r b;     // セルセンタ座標
  REAL_TYPE lb[5]; // 内外判定フラグ
  
  for (int k=1; k<=kx; k++) {
    for (int j=box_st.y-1; j<=box_ed.y+1; j++) {
      for (int i=box_st.x-1; i<=box_ed.x+1; i++) {
        
        REAL_TYPE z = org.z + 0.5*dz + dz*(REAL_TYPE)(k-1);
        
        if ( z <= ze )
        {
          base.assign((REAL_TYPE)i-0.5, (REAL_TYPE)j-0.5, (REAL_TYPE)k-0.5);
          b = org + base*pch;
          
          p[0].assign(b.x   , b.y   , b.z   ); // p
          p[1].assign(b.x-dx, b.y   , b.z   ); // w
          p[2].assign(b.x+dx, b.y   , b.z   ); // e
          p[3].assign(b.x   , b.y-dy, b.z   ); // s
          p[4].assign(b.x   , b.y+dy, b.z   ); // n
          
          // (cx, cy, *)が球の中心
          for (int l=0; l<5; l++) {
            REAL_TYPE x = p[l].x - cx;
            REAL_TYPE y = p[l].y - cy;
            REAL_TYPE rr = sqrt(x*x + y*y);
            lb[l] = ( rr <= rs ) ? -1.0 : 1.0; // 内側がマイナス
          }
          
//.........这里部分代码省略.........

bool Line::intersect(const CylinderVolume &cyl, 
                           Real           &enter,  
                           Real           &exit ) const
{
    Real  radius = cyl.getRadius();

    Vec3r adir;
    Vec3r o_adir;
    Pnt3r apos;

    cyl.getAxis(apos, adir);

    o_adir = adir;
    adir.normalize();

    bool isect;

    Real  ln;
    Real  dl;
    Vec3r RC;
    Vec3r n;
    Vec3r D;

    RC = _pos - apos;

    n  = _dir.cross (adir);
    ln =  n  .length(    );

    if(ln == 0.f)    // IntersectionLine is parallel to CylinderAxis
    {
        D  = RC - (RC.dot(adir)) * adir;
        dl = D.length();

        if(dl <= radius)   // line lies in cylinder
        {
            enter = 0.f;
            exit  = Inf;
        }
        else
        {
            return false;
        }
    }
    else
    {
        n.normalize();

        dl    = osgAbs(RC.dot(n));        //shortest distance
        isect = (dl <= radius);

        if(isect)
        {                 // if ray hits cylinder
            Real  t;
            Real  s;
            Vec3r O;

            O = RC.cross(adir);
            t = - (O.dot(n)) / ln;
            O = n.cross(adir);

            O.normalize();

            s = osgAbs (
                (osgSqrt ((radius * radius) - (dl * dl))) / (_dir.dot(O)));

            exit = t + s;

            if(exit < 0.f)
                return false;

            enter = t - s;

            if(enter < 0.f)
                enter = 0.f;
        }
        else
        {
            return false;
        }
    }

    Real t;

    Plane bottom(-adir, apos);

    if(bottom.intersect(*this, t))
    {
        if(bottom.isInHalfSpace(_pos))
        {
            if(t > enter) 
                enter = t;
        }
        else
        {
            if(t < exit) 
                exit = t;
        }
    }
    else
    {
//.........这里部分代码省略.........

/*! Intersect the line with a triangle.
    \param [in] v0,v1,v2  Points definiting a triangle in CW orientation.
    \param [out] t  If hit, this returns the distance the hit is down the line.
    \param [in,out] norm If non-NULL, this is set to the normal at the point of
            intersection
    \returns True if there is an intersection.

    This algorithm is based on "Fast, Minimum Storage Ray/Triangle Instersection" by
    T. Moeller and B. Trumbore, with the addition of avoiding the computation of
    inv_det when no intersection happens.
 */
bool Line::intersect(const Pnt3r &v0, 
                     const Pnt3r &v1,
                     const Pnt3r &v2, 
                           Real  &t,
                           Vec3r *norm) const
{
    // Eps (1E-6f) didn't work with very small geometries!
    static const Real sEps = 1E-10f;

    // find vectors for two edges sharing v0.
    Vec3r edge1 = v1 - v0;
    Vec3r edge2 = v2 - v0;

    // begin calculating determinant - also used to calculate U parameter.
    Vec3r pvec = _dir.cross(edge2);

    // if determinant is near zero, ray lies in plane of triangle.
    Real det = edge1.dot(pvec);

    Vec3r qvec;

    if(det > sEps)
    {
        // calculate distance from v0 to ray origin.
        Vec3r tvec = _pos - v0;

        // calculate U parameter and test bounds.
        Real u = tvec.dot(pvec);

        if(u < 0.f || u > det)
        {
            return false;
        }

        // prepare to test V parameter.
        qvec = tvec.cross(edge1);

        // calculate V parameter and test bounds.
        Real v = _dir.dot(qvec);

        if(v < 0.f || u + v > det)
        {
            return false;
        }
    }
    else if(det < -sEps)
    {
        // calculate distance from v0 to ray origin.
        Vec3r tvec = _pos - v0;

        // calculate U parameter and test bounds.
        Real u = tvec.dot(pvec);

        if(u > 0.f || u < det)
        {
            return false;
        }

        // prepare to test V parameter.
        qvec = tvec.cross(edge1);

        // calculate V parameter and test bounds.
        Real v = _dir.dot(qvec);

        if(v > 0.f || u + v < det)
        {
            return false;
        }
    }
    else
    {
        return false;  // ray is parallel to the plane of the triangle.
    }

    Real inv_det = 1.0f / det;

    // calculate t, ray intersects triangle.
    t = edge2.dot(qvec) * inv_det;

    if(norm != NULL)
    {
        *norm = edge1.cross(edge2);
        norm->normalize();
    }

    return true;
}

// assume the point is inside the triangle
Vec3r ComputeBarycentricCoords(const Vec3r A,const Vec3r B,const Vec3r C,const Vec3r P)
{
  Vec3r BA = B-A;
  Vec3r CA = C-A;
  Vec3r N = BA^CA;
  
  Vec3r PA = A-P;
  Vec3r PB = B-P;
  Vec3r PC = C-P;

  double areaA = (PB^PC) * N;
  double areaB = (PC^PA) * N;
  double areaC = (PA^PB) * N;

#if 0
  // checking
  printf("Barycentric debug: ----------------------------------\n");
  printf("A = [%f %f %f], B = [%f %f %f], C = [%f %f %f]\n",
	 A.x(), A.y(), A.z(),	 B.x(), B.y(), B.z(),	 C.x(), C.y(), C.z());
  printf("P = [%f %f %f]\n", P.x(), P.y(), P.z());
  printf("areas = [%f %f %f]\n", areaA, areaB, areaC);
  printf("plot3([%f %f %f %f],[%f %f %f %f],[%f %f %f %f],'r-');\n",
	 A.x(), B.x(), C.x(), A.x(),
	 A.y(), B.y(), C.y(), A.y(),
	 A.z(), B.z(), C.z(), A.z());
  printf("hold on; plot3(%f,%f,%f,'bo');\n",
	 P.x(),P.y(),P.z());
#endif


  assert(areaA > -0.1 && areaB > -0.1 && areaC > -0.1);
  if (areaA < 0)
    areaA = 0;
  if (areaB < 0)
    areaB = 0;
  if (areaC < 0)
    areaC = 0;

  double totalArea = areaA + areaB + areaC;

  double a = areaA / totalArea;
  double b = areaB / totalArea;
  double c = areaC / totalArea;


  //  printf("c = [%f %f %f]\n", a,b,c);

  Vec3r result(a,b,c);

  return result;
}