自学教程：C++ radToDeg函数代码示例

double calcSunsetUTC(double JD, double latitude, double longitude){    double t = calcTimeJulianCent(JD);    // *** First pass to approximate sunset    double  eqTime = calcEquationOfTime(t);    double  solarDec = calcSunDeclination(t);    double  hourAngle = calcHourAngleSunset(latitude, solarDec);    double  delta = longitude - radToDeg(hourAngle);    double  timeDiff = 4 * delta;	// in minutes of time	    double  timeUTC = 720 + timeDiff - eqTime;	// in minutes	    double  newt = calcTimeJulianCent(calcJDFromJulianCent(t) + timeUTC/1440.0);     eqTime = calcEquationOfTime(newt);    solarDec = calcSunDeclination(newt);    hourAngle = calcHourAngleSunset(latitude, solarDec);    delta = longitude - radToDeg(hourAngle);    timeDiff = 4 * delta;    timeUTC = 720 + timeDiff - eqTime; // in minutes    // printf("************ eqTime = %f  /nsolarDec = %f /ntimeUTC = %f/n/n",eqTime,solarDec,timeUTC);    return timeUTC;}

//Hakee yhden muuttujan uudet arvot __declspec( dllexport ) void pull( const void * _in, int in_size, void * _out, int out_sz ){	int index;	Variable *var;	cpBody *body;	cpShape *shape;	index = PEEKINT(INPUT_MEMBLOCK,0);	var = vhGetVariable(&mVariableHandler,index);	switch (var->mType)	{	case VarTypeBody:		body = (cpBody*)var->mPtr;				POKEFLOAT(OUTPUT_MEMBLOCK,0,radToDeg(cpBodyGetAngle(body)));		POKEVECT(OUTPUT_MEMBLOCK,4,cpBodyGetPos(body));		POKEFLOAT(OUTPUT_MEMBLOCK,12,radToDeg(cpBodyGetAngVel(body)));		POKEVECT(OUTPUT_MEMBLOCK,16,cpBodyGetVel(body));		break;	case VarTypeShape:		shape = (cpShape*)var->mPtr;		POKEINT(OUTPUT_MEMBLOCK,0,VarTypeShape);			default:		MessageBoxA(NULL,"cpPull: Invalid variable type","cpChipmunk error",MB_OK);		exit(0);		break;	}}

void SceneNode::getTransform( Vec3f &trans, Vec3f &rot, Vec3f &scale ){	_relTrans.decompose( trans, rot, scale );	rot.x = radToDeg( rot.x );	rot.y = radToDeg( rot.y );	rot.z = radToDeg( rot.z );}

void IMU::LogState(int channel, char *device, int value) {#ifdef _DEBUG	struct _timeb tstruct;	char *timeline;	char buffer[100];		char buffer1[100];		_ftime(&tstruct);	timeline = ctime(&(tstruct.time));	strcpy(buffer, timeline + 11);		intToBinaryString(buffer1, value);		fprintf(logFile, "%.8s.%03hu Ch %03o %s %s PIPA %o %o %o CDUCMD %o %o %o GYRO %o IMU %.2f %.2f %.2f/n", buffer, tstruct.millitm, channel, device, buffer1,		    agc.GetErasable(0, RegPIPAX), 		    agc.GetErasable(0, RegPIPAY), 			agc.GetErasable(0, RegPIPAZ), 			agc.GetErasable(0, RegCDUXCMD),			agc.GetErasable(0, RegCDUYCMD),			agc.GetErasable(0, RegCDUZCMD),			//state->Erasable[0][RegCDUX],			//state->Erasable[0][RegCDUY],			//state->Erasable[0][RegCDUZ],						agc.GetErasable(0, RegGYROCTR),			radToDeg(Gimbal.X),			radToDeg(Gimbal.Y),			radToDeg(Gimbal.Z));				fflush(logFile);#endif}

 bool WindowGroup::operator==(const WindowGroup& other) const{  if (!m_shadingControl && !other.shadingControl()){    if (m_space.handle() == other.space().handle()){      if (m_construction.handle() == other.construction().handle()){        double angle = std::abs(radToDeg(getAngle(m_outwardNormal, other.outwardNormal())));        const double tol = 1.0;        if (angle < tol){          return true;        }      }    }    return false;  }else if (m_shadingControl && !other.shadingControl()){    return false;  }else if (!m_shadingControl && other.shadingControl()){    return false;  }  if (m_space.handle() == other.space().handle()){    if (m_construction.handle() == other.construction().handle()){      if (m_shadingControl->handle() == other.shadingControl()->handle()){        double angle = std::abs(radToDeg(getAngle(m_outwardNormal, other.outwardNormal())));        const double tol = 1.0;        if (angle < tol){          return true;        }      }    }  }  return false;}

  SensorRunnable(const sensors_event_t& data, const sensor_t* sensors, ssize_t size)  {    mSensorData.sensor() = HardwareSensorToHalSensor(data.type);    mSensorData.accuracy() = HardwareStatusToHalAccuracy(SensorseventStatus(data));    mSensorData.timestamp() = data.timestamp;    if (mSensorData.sensor() == SENSOR_GYROSCOPE) {      // libhardware returns gyro as rad.  convert.      mSensorValues.AppendElement(radToDeg(data.data[0]));      mSensorValues.AppendElement(radToDeg(data.data[1]));      mSensorValues.AppendElement(radToDeg(data.data[2]));    } else if (mSensorData.sensor() == SENSOR_PROXIMITY) {      mSensorValues.AppendElement(data.data[0]);      mSensorValues.AppendElement(0);      // Determine the maxRange for this sensor.      for (ssize_t i = 0; i < size; i++) {        if (sensors[i].type == SENSOR_TYPE_PROXIMITY) {          mSensorValues.AppendElement(sensors[i].maxRange);        }      }    } else if (mSensorData.sensor() == SENSOR_LIGHT) {      mSensorValues.AppendElement(data.data[0]);    } else {      mSensorValues.AppendElement(data.data[0]);      mSensorValues.AppendElement(data.data[1]);      mSensorValues.AppendElement(data.data[2]);    }    mSensorData.values() = mSensorValues;  }

void InfoPanel::buildDSOPage(const DeepSkyObject* dso,                             const Universe* universe,                             QTextStream& stream){    string name = universe->getDSOCatalog()->getDSOName(dso);    stream << "<h1>" << name.c_str() << "</h1><br>/n";        Vec3d eqPos = astro::eclipticToEquatorial(celToJ2000Ecliptic(dso->getPosition()));    Vec3d sph = rectToSpherical(eqPos);    int hours = 0;    int minutes = 0;    double seconds = 0;    astro::decimalToHourMinSec(radToDeg(sph.x), hours, minutes, seconds);    stream << "RA: " << hours << "h " << abs(minutes) << "m " << abs(seconds) << "s<br>/n";    int degrees = 0;    astro::decimalToDegMinSec(radToDeg(sph.y), degrees, minutes, seconds);    stream << "Dec: " << degrees << QString::fromUtf8(UTF8_DEGREE_SIGN) << " " <<        abs(minutes) << "' " << abs(seconds) << "/"<br>/n";    Vec3d galPos = astro::equatorialToGalactic(eqPos);    sph = rectToSpherical(galPos);        astro::decimalToDegMinSec(radToDeg(sph.x), degrees, minutes, seconds);    stream << "L: " << degrees << QString::fromUtf8(UTF8_DEGREE_SIGN) << " " <<        abs(minutes) << "' " << abs(seconds) << "/"<br>/n";    astro::decimalToDegMinSec(radToDeg(sph.y), degrees, minutes, seconds);    stream << "B: " << degrees << QString::fromUtf8(UTF8_DEGREE_SIGN) << " " <<        abs(minutes) << "' " << abs(seconds) << "/"<br>/n";}

void Polar::bvmg(double bt_longitude,double bt_latitude, double wp_longitude, double wp_latitude,                  double w_angle, double w_speed,                  double *heading, double *wangle){    double maxwangle,twaOrtho;    double maxheading;    double wanted_heading;    w_angle=degToRad(w_angle);    Orthodromie orth(bt_longitude,bt_latitude,wp_longitude,wp_latitude);    wanted_heading=orth.getAzimutRad();    twaOrtho = w_angle - wanted_heading;    myBvmgWind(twaOrtho,w_speed,&maxwangle);    maxheading = fmod((wanted_heading-maxwangle+twaOrtho), TWO_PI);    if (maxheading < 0)    {      maxheading += TWO_PI;    }    maxwangle = fmod(maxwangle, TWO_PI);    if (maxwangle > PI)    {      maxwangle -= TWO_PI;    } else if (maxwangle < -PI)    {      maxwangle += TWO_PI;    }#if 0  qWarning() << "BVMG: Wind " << w_speed << "kts " << radToDeg(w_angle);  qWarning() << "BVMG Wind Angle : wanted_heading " << radToDeg(wanted_heading);  qWarning() << "BVMG Wind Angle : heading " << radToDeg(maxheading) << ", wind angle " << radToDeg(maxwangle);#endif /* DEBUG */    *heading = radToDeg(maxheading);    *wangle = radToDeg(maxwangle);}

/*** update the Horde tentacle according to Smr _tentacleIK*/static void updateTentacle(){  //set of floats storing Horde's tentacle joints state  float tx,ty,tz,rx,ry,rz,foo;  double rrx,rry,rrz;  // a Horde handle towards a tentacle Bone (joint)  H3DNode tentacleBone;  // and its corresponding node in Smr  SMRKinematicJoint* ikJoint;  //get the transformation parameters from the SmrTentacle  for( unsigned int i=0; i < _tentacleIK->getNumJoints()-1 ; i++)  {    ikJoint = _tentacleIK->getJoint(i);    //get Horde tentacle's corresponding joint joint (Horde)     h3dFindNodes(H3DRootNode,ikJoint->getName().c_str(), H3DNodeTypes::Joint);    tentacleBone = h3dGetNodeFindResult(0);    h3dGetNodeTransform(tentacleBone,&tx,&ty,&tz,&rx,&ry,&rz,&foo,&foo,&foo);    //put the rotation values into euler angles (in degrees)    ikJoint->getOrientation().toEulerAngles(rrx,rry,rrz);    // update tentacle's joint according to Smr kinematic chain equivalent    h3dSetNodeTransform( tentacleBone, tx, ty, tz, radToDeg( ((float)(rrx)) ), radToDeg( ((float)(rry)) ), radToDeg( ((float)(rrz)) ), foo, foo, foo );  }  // orientate Horde tentacles root joint correctly //No more in the new version !  ikJoint = _tentacleIK->getJoint(0);  h3dFindNodes(H3DRootNode,ikJoint->getName().c_str(), H3DNodeTypes::Joint);  tentacleBone = h3dGetNodeFindResult(0);  h3dGetNodeTransform(tentacleBone,&tx,&ty,&tz,&rx,&ry,&foo,&foo,&foo,&foo);  h3dSetNodeTransform(tentacleBone,tx,ty,tz,rx,ry,rz,foo,foo,foo);}

void Camera::lookAt(){	assert(m_target != m_pos);	glm::vec3 direction = glm::normalize(m_target - m_pos);	m_phi = radToDeg(asinf(-direction.y));	m_theta = -radToDeg(atan2f(-direction.x, -direction.z));	normalizeAngles();}

Rotation *getAngles(point3d *p1, point3d* p2){	point3d *d = new point3d(p2->x - p1->x, p2->y - p1->y, p2->z - p1->z);	float heading = radToDeg(atan2(d->x,d->z));	float pitch =  radToDeg(atan2(d->y, sqrt( (d->x * d->x) + (d->z * d->z) )));	if(d->x < 0)			pitch *= -1;	return new Rotation(heading, pitch, 0);}

void Polar::myBvmgWind(double w_angle, double w_speed, double *wangle){    /* FIXME, this can be optimized a lot */    double speed=0;    double t_heading=0;    double t=0;    double imax = 90;    double t_max = -100;    double t_max2 = -100;    int i=0;    /* -90 to +90 form desired diretion */    for (i=0; i<imax; i++)    {        t_heading = w_angle + degToRad(((double)i));        speed = this->getSpeed(w_speed, A180(radToDeg(t_heading)));        if (speed < 0.0)        {            continue;        }        t = speed * cos(degToRad(((double)i)));        if (t > t_max)        {            t_max = t;            *wangle = t_heading;        } else if ( t_max - t > (t_max/20.0))        {            break;  /* cut if lower enough from current maximum */        }    }    for (i=0; i<imax; i++)    {        t_heading = w_angle - degToRad(((double)i));        speed = this->getSpeed(w_speed,  A180(radToDeg(t_heading)));        if (speed < 0.0)        {          continue;        }        t = speed * cos(- degToRad(((double)i)));        if (t > t_max2)        {          t_max2 = t;          if (t > t_max)          {            t_max = t;            *wangle = t_heading;          }        } else if (t_max2 - t > (t_max2/20.0))        {          break;        }    }}

void Camera::lookAt(Vec3 position){	if (position == mPosition)	{		std::cout << "MEGA ERROR: You are trying to look at your origin" << std::endl;		return;	}	Vec3 direction = glm::normalize(position - mPosition);	mVerticalAngle = radToDeg(asinf(-direction.y));	mHorizontalAngle = -radToDeg(atan2f(-direction.x, -direction.z));	normalizeAngles();}

Matrix3x4 SceneNode::getLocalTransform(){	Matrix3x4 transform = Matrix3x4();	transform.translate(position_);	transform.rotateX(radToDeg(rotation_.x_));	transform.rotateY(radToDeg(rotation_.y_));	transform.rotateZ(radToDeg(rotation_.z_));	transform.scale(scale_);	return transform;}

  bool GlareSensor_Impl::setTransformation(const openstudio::Transformation& transformation)   {    Vector3d translation = transformation.translation();    this->setPositionXCoordinate(translation.x());    this->setPositionYCoordinate(translation.y());    this->setPositionZCoordinate(translation.z());    EulerAngles eulerAngles = transformation.eulerAngles();    setPsiRotationAroundXAxis(radToDeg(eulerAngles.psi()));    setThetaRotationAroundYAxis(radToDeg(eulerAngles.theta()));    setPhiRotationAroundZAxis(radToDeg(eulerAngles.phi()));        return true;    }

Foam::tmp<Foam::scalarField> Foam::cellQuality::nonOrthogonality() const{    tmp<scalarField> tresult    (        new scalarField        (            mesh_.nCells(), 0.0        )    );    scalarField& result = tresult();    scalarField sumArea(mesh_.nCells(), 0.0);    const vectorField& centres = mesh_.cellCentres();    const vectorField& areas = mesh_.faceAreas();    const labelList& own = mesh_.faceOwner();    const labelList& nei = mesh_.faceNeighbour();    forAll(nei, faceI)    {        vector d = centres[nei[faceI]] - centres[own[faceI]];        vector s = areas[faceI];        scalar magS = mag(s);        scalar cosDDotS =            radToDeg(Foam::acos(min(1.0, (d & s)/(mag(d)*magS + VSMALL))));        result[own[faceI]] = max(cosDDotS, result[own[faceI]]);        result[nei[faceI]] = max(cosDDotS, result[nei[faceI]]);    }

void GraphicsSubsystem::render(Renderable::Ptr renderable) {	Rect boundingRect = renderable->boundingRect();	Vector2 pos = renderable->worldPos();	if (renderable->drawLayer() == DrawLayer::World) {		if (m_camera) { // Parallax			Vector2 camPos = m_camera->worldPos();			pos += (camPos - pos) * renderable->parallax();		}		boundingRect.translate(pos);		if (!viewport().intersects(boundingRect))			return;	} else {		boundingRect.translate(pos);		if (!m_screen.intersects(boundingRect))			return;	}	// Set position and scale	glPushMatrix();	glTranslatef(pos.x, pos.y, 0);	glRotatef(radToDeg(renderable->worldAngle()), 0, 0, 1.0f);	glScalef(renderable->scale(), renderable->scale(), 1.0f);	renderable->render();	glPopMatrix();}

void Cursor::onProcess(const WidgetEvents & events){	if (!isVisible())		return;	sf::Vector2f pos = events.mousePosition - getPosition();	if (myCursorTransform.getPosition() != pos)	{		myCursorTransform.setPosition(pos);		invalidateVertices();	}	if (isRotatable())	{		sf::Vector2f deltaVec = myBackPosition - myCursorTransform.getPosition();		sf::Vector2f normVec = vectorLength(deltaVec) < 0.0001f ? sf::Vector2f(1, 0) : deltaVec / vectorLength(deltaVec);		myBackPosition = myCursorTransform.getPosition() + normVec * getRotationOffset();		float angle = radToDeg(angleBetween(myBackPosition, myCursorTransform.getPosition())) + getAngle();		myCursorTransform.setRotation(angle);		/*		if (myCursorTransform.getRotation() != angle)		{			myCursorTransform.setRotation(angle);			invalidateVertices();		}		*/	}}

Foam::cylindricalCS Foam::arcEdge::calcAngle(){    vector a = p2_ - p1_;    vector b = p3_ - p1_;    // find centre of arcEdge    scalar asqr = a & a;    scalar bsqr = b & b;    scalar adotb = a & b;    scalar denom = asqr*bsqr - adotb*adotb;    if (mag(denom) < VSMALL)    {        FatalErrorInFunction            << denom            << abort(FatalError);    }    scalar fact = 0.5*(bsqr - adotb)/denom;    point centre = 0.5*a + fact*((a ^ b) ^ a);    centre += p1_;    // find position vectors w.r.t. the arcEdge centre    vector r1(p1_ - centre);    vector r2(p2_ - centre);    vector r3(p3_ - centre);    // find angles    angle_ = radToDeg(acos((r3 & r1)/(mag(r3) * mag(r1))));    // check if the vectors define an exterior or an interior arcEdge    if (((r1 ^ r2) & (r1 ^ r3)) < 0.0)    {        angle_ = 360.0 - angle_;    }    vector tempAxis;    if (angle_ <= 180.0)    {        tempAxis = r1 ^ r3;        if (mag(tempAxis)/(mag(r1)*mag(r3)) < 0.001)        {            tempAxis = r1 ^ r2;        }    }    else    {        tempAxis = r3 ^ r1;    }    radius_ = mag(r3);    // set up and return the local coordinate system    return cylindricalCS("arcEdgeCS", centre, tempAxis, r1);}

void Foam::fv::LeishmanBeddoesSD::correct(    scalar magU,    scalar alphaDeg,    scalar& cl,    scalar& cd,    scalar& cm){    LeishmanBeddoesSGC::correct    (        magU,        alphaDeg,        cl,        cd,        cm    );    // CTCorrection is a correction to ensure that the model reduces to static    // values when pitch rate and angle of attack is low enough.    if (CTCorrection_)    {        const scalar CC_corr_const = 1;        const scalar scaleStart = CC_corr_const * 0.5 * alphaSS_;        const scalar scaleEnd = CC_corr_const * alphaSS_;        const scalar rScaleStart = 0.01;        const scalar rScaleEnd = 0.02;        scalar alphaDeg = radToDeg(alpha_);        scalar scaleFactorAlpha =            (scaleEnd - fabs(alphaDeg))/(scaleEnd-scaleStart);        if (scaleFactorAlpha < 0)        {            scaleFactorAlpha = 0;        }        if (scaleFactorAlpha > 1)        {            scaleFactorAlpha = 1;        }        scalar scaleFactor_r =            (rScaleEnd - abs(r_))/(rScaleEnd-rScaleStart);        if (scaleFactor_r < 0)        {            scaleFactor_r = 0;        }        if (scaleFactor_r > 1)        {            scaleFactor_r = 1;        }        scalar scaleFactor = scaleFactorAlpha*scaleFactor_r;        scalar CTProfileData = profileData_.chordwiseCoefficient(alphaDeg);        scalar CTVal = CT_ + (CTProfileData - CTStatic_)*scaleFactor;        cl = CN_*cos(alpha_) + CTVal*sin(alpha_);        cd = CN_*sin(alpha_) - CTVal*cos(alpha_) + CD0_*(1 - scaleFactor);    }}

float RB2DBody::angle() const{    if (mBody) {        mAngle = radToDeg(mBody->GetAngle());    }    return mAngle;}

float RCPBody::angle() const{    if (mType != Static && mBody) {        mAngle = radToDeg(cpBodyGetAngle(mBody));    }    return mAngle;}

//helper function - converts x pixel at a cetrain zoom  to a longitude double Calc::xAtZoomToLongitude(int x, int zoom ){    //OK    double arc = (double) EARTHCIRCUM / ((1 << zoom) * 256.0);    double lon = (( (double) x * arc ) - (EARTHCIRCUM/2) ) / (double) EARTHRADIUS;    return radToDeg(lon);}

double calcSunDeclination(double t){    double e = calcObliquityCorrection(t);    double lambda = calcSunApparentLong(t);    double sint = sin(degToRad(e)) * sin(degToRad(lambda));    double theta = radToDeg(asin(sint));    return theta; // in degrees}

void Affine3DQ::prettyPrint1(std::ostream &out)  const{    Quaternion o = rotor.normalised();    Vector3dd axis = o.getAxis();    double   angle = radToDeg(o.getAngle());    out << "Pos:" <<  shift << std::endl;    out << "Rotation around: " << axis << " angle " << angle << "deg" << std::endl;}

  bool IlluminanceMap_Impl::setTransformation(const openstudio::Transformation& transformation)   {    bool test;    Vector3d translation = transformation.translation();    this->setOriginXCoordinate(translation.x());    this->setOriginYCoordinate(translation.y());    this->setOriginZCoordinate(translation.z());    EulerAngles eulerAngles = transformation.eulerAngles();    test = this->setPsiRotationAroundXAxis(radToDeg(eulerAngles.psi()));    OS_ASSERT(test);    test = this->setThetaRotationAroundYAxis(radToDeg(eulerAngles.theta()));    OS_ASSERT(test);    test = this->setPhiRotationAroundZAxis(radToDeg(eulerAngles.phi()));    OS_ASSERT(test);    return true;  }

void draw_separator(t_meter *x,  t_object *view, t_rect *rect){    t_jgraphics *g = jbox_start_layer((t_object *)x, view, hoa_sym_separator_layer, rect->width, rect->height);    	if (g)	{        double rotateAngle, channelWidth;        t_jmatrix transform;        t_jrgba black = rgba_addContrast(x->f_color_mbg, -0.12);        		jgraphics_matrix_init(&transform, 1, 0, 0, -1, x->f_center, x->f_center);		jgraphics_set_matrix(g, &transform);				// skelton separators and leds bg:		for(int i=0; i < x->f_meter->getNumberOfPlanewaves(); i++)		{			channelWidth = radToDeg(x->f_meter->getPlanewaveWidth(i));            rotateAngle = radToDeg(x->f_meter->getPlanewaveAzimuthMapped(i) + x->f_meter->getPlanewavesRotationX()) - (channelWidth*0.5);			if (!x->f_rotation)			{				rotateAngle += channelWidth;				rotateAngle *= -1;			}			jgraphics_rotate(g, degToRad(rotateAngle));            // separator			if (x->f_meter->getNumberOfPlanewaves() > 1)			{                jgraphics_set_line_width(g, 1);                jgraphics_set_source_jrgba(g, &black);				jgraphics_move_to(g, 0., x->f_rayonInt);				jgraphics_line_to(g, 0, x->f_rayonExt);				jgraphics_stroke(g);			}						jgraphics_rotate(g, degToRad(-rotateAngle));		}				jbox_end_layer((t_object*)x, view, hoa_sym_separator_layer);	}	jbox_paint_layer((t_object *)x, view, hoa_sym_separator_layer, 0., 0.);}

/** * Given a single number between 0..1, generate a latitude, longitude (in degrees) and a 3D * (x, y, z) point on a sphere with a radius of 1. */static void unitToLatLonDeg(        const double unit1, const double unit2, double *latDeg, double *lonDeg) {    // Calculate uniformly distributed 3D point on sphere (radius = 1.0):    // http://mathproofs.blogspot.co.il/2005/04/uniform-random-distribution-on-sphere.html    const double theta0 = (2.0 * PI) * unit1;    const double theta1 = acos(1.0 - (2.0 * unit2));    double x = sin(theta0) * sin(theta1);    double y = cos(theta0) * sin(theta1);    double z = cos(theta1);    // Convert Carthesian 3D point into lat/lon (radius = 1.0):    // http://stackoverflow.com/questions/1185408/converting-from-longitude-latitude-to-cartesian-coordinates    const double latRad = asin(z);    const double lonRad = atan2(y, x);    // Convert radians to degrees.    *latDeg = my_isnan(latRad) ? 90.0 : radToDeg(latRad);    *lonDeg = my_isnan(lonRad) ? 180.0 : radToDeg(lonRad);}

        forAll(faceCentres, faceI)        {            vector d = faceCentres[faceI] - centres[faceCells[faceI]];            vector s = faceAreas[faceI];            scalar magS = mag(s);            scalar cosDDotS =                radToDeg(Foam::acos(min(1.0, (d & s)/(mag(d)*magS + VSMALL))));            result[faceCells[faceI]] = max(cosDDotS, result[faceCells[faceI]]);        }