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);
    }
}