public void AddForce(Vector worldForce, Vector worldOffset)
 {
     //add linar force
     m_forces += worldForce;
     //and it's associated torque
     m_torque += worldOffset % worldForce;
 }

        private static double Zoom(CompiledFunc func, Vector x, Vector dir, double aLow, double aHigh, double fZero, double diffZero)
        {
            var normDir = dir.Normalize();
            double aMid = 0;
            double fValue = 0;
            double diff = 0;

            while (Math.Abs(aLow - aHigh) > EPS)
            {
                aMid = aLow + (aHigh - aLow)/2;
                fValue = func.Eval(x + aMid*dir);

                if (fValue > fZero + C1*aMid*diffZero || fValue >= func.Eval(x + aLow*dir))
                    aHigh = aMid;
                else
                {
                    diff = (func.Differentiate(x + aMid*dir)*normDir).Sum();

                    if (Math.Abs(diff) <= -C2*diffZero)
                        return aMid;

                    if (diff*(aHigh - aLow) >= 0)
                        aHigh = aLow;

                    aLow = aMid;
                }
            }

            return aMid;
        }

 /// <summary>
 /// Solves a linear system through inversing a matrix
 /// </summary>
 /// <param name="paramM"></param>
 /// <param name="yV"></param>
 /// <returns></returns>
 public Vector LinearSolve(Matrix paramM, Vector yV)
 {
     if (!paramM.IsSquare() || paramM.Height != yV.Dim)
         throw new Exception("Dim Error");
     Vector xV = paramM.Inverse() * yV;
     return xV;
 }

 public static void Write(this BinaryWriter writer, Vector vector)
 {
     writer.Write(vector.X);
     writer.Write(vector.Y);
     writer.Write(vector.Z);
     writer.Write(vector.W);
 }

 public void Move(Vector offset)
 {
     foreach(Node node in Nodes) {
         node.X += offset.X;
         node.Y += offset.Y;
     }
 }

 public override float D(ref Vector wh)
 {
     float costhetah = Math.Abs(BrdfBase.CosTheta(ref wh));
     return (exponent + 2f) *
            MathLab.INVTWOPI *
            (float)Math.Pow(Math.Max(0f, costhetah), exponent);
 }

      public GaussianQuadrature(int n, GaussianOrthogonalPolynomial orthPoly)
      {
         x_ = new Vector(n);
         w_ = new Vector(n);

        // set-up matrix to compute the roots and the weights
        Vector e = new Vector(n-1);

        int i;
        for (i=1; i < n; ++i) 
        {
            x_[i] = orthPoly.alpha(i);
            e[i-1] = Math.Sqrt(orthPoly.beta(i));
        }
        x_[0] = orthPoly.alpha(0);

        TqrEigenDecomposition tqr = new TqrEigenDecomposition( x_, e,
                               TqrEigenDecomposition.EigenVectorCalculation.OnlyFirstRowEigenVector,
                               TqrEigenDecomposition.ShiftStrategy.Overrelaxation);

        x_ = tqr.eigenvalues();
        Matrix ev = tqr.eigenvectors();

        double mu_0 = orthPoly.mu_0();
        for (i=0; i<n; ++i) {
            w_[i] = mu_0*ev[0,i]*ev[0,i] / orthPoly.w(x_[i]);
        }
      }

 public PinholeCamera(Point eye, Point lat, Vector up, float fov)
 {
     this.eye = eye;
     this.lat = lat;
     this.up = up;
     this.fov = fov;
 }

 public override void DisturbVector(Vector<double> vecToBeDisturbed)
 {
     double[] samples = new double[vecToBeDisturbed.Count];
     _gauss.Samples(samples);
     for(int i = 0; i < vecToBeDisturbed.Count; ++i)
         vecToBeDisturbed.At(i, vecToBeDisturbed.At(i) + (double)samples[i]);
 }

	private Notes _notes; // usermodel needs to Set this

	/**
	 * Constructs a Slide from the Slide record, and the SlideAtomsSet
	 *  Containing the text.
	 * Initialises TextRuns, to provide easier access to the text
	 *
	 * @param slide the Slide record we're based on
	 * @param notes the Notes sheet attached to us
	 * @param atomSet the SlideAtomsSet to Get the text from
	 */
	public Slide(NPOI.HSLF.record.Slide slide, Notes notes, SlideAtomsSet atomSet, int slideIdentifier, int slideNumber) {
        base(slide, slideIdentifier);

		_notes = notes;
		_atomSet = atomSet;
		_slideNo = slideNumber;

 		// Grab the TextRuns from the PPDrawing
		TextRun[] _otherRuns = FindTextRuns(getPPDrawing());

		// For the text coming in from the SlideAtomsSet:
		// Build up TextRuns from pairs of TextHeaderAtom and
		//  one of TextBytesAtom or TextCharsAtom
		Vector textRuns = new Vector();
		if(_atomSet != null) {
			FindTextRuns(_atomSet.GetSlideRecords(),textRuns);
		} else {
			// No text on the slide, must just be pictures
		}

		// Build an array, more useful than a vector
		_Runs = new TextRun[textRuns.Count+_otherRuns.Length];
		// Grab text from SlideListWithTexts entries
		int i=0;
		for(i=0; i<textRuns.Count; i++) {
			_Runs[i] = (TextRun)textRuns.Get(i);
            _Runs[i].SetSheet(this);
		}
		// Grab text from slide's PPDrawing
		for(int k=0; k<_otherRuns.Length; i++, k++) {
			_Runs[i] = _otherRuns[k];
            _Runs[i].SetSheet(this);
		}
	}

 public static void Write(this BinaryWriter bw, Vector vector)
 {
     bw.Write(vector.X);
     bw.Write(vector.Y);
     bw.Write(vector.Z);
     bw.Write(vector.W);
 }

 /// <summary>
 /// VMP message to 'A'.
 /// </summary>
 /// <param name="X">Incoming message from 'X'. Must be a proper distribution.  If any element is uniform, the result will be uniform.</param>
 /// <param name="A">Incoming message from 'A'.</param>
 /// <param name="result">Modified to contain the outgoing message.</param>
 /// <returns><paramref name="result"/></returns>
 /// <remarks><para>
 /// The outgoing message is the exponential of the integral of the log-factor times incoming messages, over all arguments except 'A'.
 /// The formula is <c>int log(f(A,x)) q(x) dx</c> where <c>x = (X,B)</c>.
 /// </para></remarks>
 /// <exception cref="ImproperMessageException"><paramref name="X"/> is not a proper distribution</exception>
 /// <exception cref="ImproperMessageException"><paramref name="B"/> is not a proper distribution</exception>
 public static VectorGaussianArray AAverageLogarithm([SkipIfAllUniform] GaussianArray2D X, [SkipIfAllUniform] VectorGaussianArray B, [SkipIfAllUniform] VectorGaussianArray result)
 {
     int I = X.GetLength(0), J = X.GetLength(1), K = B[0].Dimension;
     var Ebbp = new PositiveDefiniteMatrix[J];
     var mb = new Vector[J];
     for (int j = 0; j < J; j++)
     {
         Ebbp[j] = new PositiveDefiniteMatrix(K, K);
         mb[j] = Vector.Zero(K);
         B[j].GetMeanAndVariance(mb[j], Ebbp[j]);
         Ebbp[j].SetToSumWithOuter(Ebbp[j], 1, mb[j], mb[j]);
     }
     for (int i = 0; i < I; i++)
     {
         result[i].Precision.SetAllElementsTo(0);
         result[i].MeanTimesPrecision.SetAllElementsTo(0);
         for (int j = 0; j < J; j++)
         {
             // nb: would be more memory efficient to have a SetToAPlusCB routine
             result[i].Precision.SetToSum(result[i].Precision, Ebbp[j] * X[i, j].Precision);
             result[i].MeanTimesPrecision.SetToSum(result[i].MeanTimesPrecision, mb[j] * X[i, j].MeanTimesPrecision);
         }
     }
     return result;
 }

 public UserCmd()
 {
     ViewAngles = new Vector();
     HeadAngles = new Vector();
     HeadOffset = new Vector();
     CrosshairTrace = new Vector();
 }

        // f' by x  (df[i] by dx[j])
        //пользователь вводит производную каждой функции фи по каждому х
        public double dfx(Vector xt, double t, ref Vector z, int i, int j)
        {
            double[] ut = new double[dimU];
            Program.control(t, ref z, ref ut);

            return evaluate(ref xdtdx[i, j], ref xt, ut[0], 0, t);
        }

    public Vector CalculateForce(Vector relativeGroundSpeed, float timeStep)
    {
        //calculate speed of tire patch at ground
        Vector patchSpeed = -m_forwardAxis * m_wheelSpeed * m_wheelRadius;

        //get velocity difference between ground and patch
        Vector velDifference = relativeGroundSpeed + patchSpeed;

        //project ground speed onto side axis
        float forwardMag = 0;
        Vector sideVel = velDifference.Project(m_sideAxis);
        Vector forwardVel = velDifference.Project(m_forwardAxis, out forwardMag);

        //calculate super fake friction forces
        //calculate response force
        Vector responseForce = -sideVel * 2.0f;
        responseForce -= forwardVel;

        //calculate torque on wheel
        m_wheelTorque += forwardMag * m_wheelRadius;

        //integrate total torque into wheel
        m_wheelSpeed += m_wheelTorque / m_wheelInertia * timeStep;

        //clear our transmission torque accumulator
        m_wheelTorque = 0;

        //return force acting on body
        return responseForce;
    }

    private static int VectorIntEquals()
    {
        const int Pass = 100;
        const int Fail = -1;

        Vector<int> A = new Vector<int>(3);
        Vector<int> B = new Vector<int>(3);
        Vector<int> C = new Vector<int>(5);

        bool result = A.Equals(B);
        if (!result)
        {
            return Fail;
        }

        result = A.Equals(C);
        if (result)
        {
            return Fail;
        }

        if (A.Equals(Vector<int>.Zero))
        {
            return Fail;
        }

        if (!Vector<int>.Zero.Equals(Vector<int>.Zero))
        {
            return Fail;
        }

        return Pass;
    }

        public void Test()
        {
            double k = 1.0;
            double theta = 0.02;
            double sigma = 0.08;
            double r0 = 0.01;
            Vector par = new Vector(3);
            par[0] = k;
            par[1] = theta;
            par[2] = sigma;

            double strike = 0.98;
            double noz = 100.0;

            double T = 1.0;
            double S = 2.0;

            double callFairmat = noz * CIRCap.BondCall(r0, 0.0, T, S, strike, par);
            double callBenchmark = 0.317304505282290;

            Console.WriteLine("CallFairmat   = " + callFairmat);
            Console.WriteLine("CallBenchmark = " + callBenchmark);

            double maxError = 1e-10;
            Assert.Less(Math.Abs(callFairmat - callBenchmark), maxError);
        }

        public Line(Point P, Point Q)
        {
            if(P.X == Q.X) // se for vertical
            {
                _a = new Point((P.Y < Q.Y) ? P : Q);
                _b = new Point((P.Y < Q.Y) ? Q : P);
                _center = new Point(_a.X, 0);
                _radius = 0;
                _alfa = new Point (Center.X,0);
                _beta = new Point (Center.X,-1);
                _director = new Vector(new Point(0,1));

            }
            else
            {
                _a = new Point((P.X < Q.X) ? P : Q);
                _b = new Point((P.X < Q.X) ? Q : P);
                if (A.Y == B.Y)
                    _center = new Point((A.X + B.X) / 2, 0);
                else
                    _center = new Point(((A.X * A.X) - (B.X * B.X) + (A.Y * A.Y) * (B.Y * B.Y)) / (2 * (A.X - B.X)), 0);
                _radius = (float)A.EuclidianDistance(Center);
                _alfa = new Point(Center.X-Radius, 0);
                _beta = new Point(Center.X+Radius, 0);
                _director = new Vector(Center, A).Normal;

            }
        }

 /// <summary>
 ///   Initializes a new instance of the <see cref = "ItemAutoSubscription" /> class.
 /// </summary>
 /// <param name = "item">
 ///   The subscribed item.
 /// </param>
 /// <param name = "itemPosition">
 ///   The item position.
 /// </param>
 /// <param name = "itemRegion">
 ///   The item Region.
 /// </param>
 /// <param name = "subscription">
 ///   The subscription.
 /// </param>
 public ItemAutoSubscription(Item item, Vector itemPosition, Region itemRegion, IDisposable subscription)
 {
     this.ItemPosition = itemPosition;
     this.item = item;
     this.subscription = subscription;
     this.WorldRegion = itemRegion;
 }

 public override float Pdf(ref Vector wo, ref Vector wi)
 {
     Vector H = (wo + wi).Normalize();
     float costheta = Math.Abs(H.z);
     float blinn_pdf = ((exponent + 2f) * (float)Math.Pow(costheta, exponent)) / (2f * MathLab.M_PI * 4f * (wo & H));
     return blinn_pdf;
 }