pse-pro-4.17.cs

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace PSE_eq4_8_eq4_9
{
    public static class MathExtensionMethods
    {
        public static double Sqrt(this double n)
        { return Math.Sqrt(n); }

        public static double Sin(this double n)
        { return Math.Sin(n); }

        public static double Cos(this double n)
        { return Math.Cos(n); }

        public static double Asin(this double n)
        { return Math.Asin(n); }

        public static double Acos(this double n)
        { return Math.Acos(n); }

        public static double Pow(this double a, double b)
        { return Math.Pow(a, b); }

        public static double Degrees(this double n)
        { return 180 * n / Math.PI; }

        public static double Radians(this double n)
        { return n * Math.PI / 180; }
    }

    class Point
    {
        public double? x, y;

        public double X { get { return x.Value; } }
        public double Y { get { return y.Value; } }

        public override string ToString()
        { return string.Format("Point({0,7:F3}, {1,7:F3})", x, y); }

        public static Point FromAngle(double angle, double mag)
        { return new Point(angle.Cos() * mag, angle.Sin() * mag); }

        public Point(double? x_val, double? y_val) { x = x_val; y = y_val; }

        public static Point operator +(Point a, Point b)
        { return new Point(a.x + b.x, a.y + b.y); }

        public static Point operator -(Point a, Point b)
        { return new Point(a.x - b.x, a.y - b.y); }

        public static Point operator *(Point p, double n)
        {return new Point(p.x * n, p.y * n); }

        public static Point operator *(double n, Point p)
        { return new Point(p.x * n, p.y * n); }

        public static Point operator /(Point p, double n)
        {return new Point(p.x / n, p.y / n); }

        public static Point operator /(double n, Point p)
        { return new Point(p.x / n, p.y / n); }

        public double Norm()
        { return Math.Sqrt((double)(x * x + y * y)); }

        public double ToAngle()
        {
            return Math.Atan2((double) y, (double) x);
        }
    }

    // two-dimensional motion with constant acceleration

    //          vf_ = vi_ + a_ t        (1)

    // (1) vi   vi_ = vf_ - a_ t        (1.1)

    // (1) a    a_ = (vf_ - vi_) / t    (1.2)

    // (1) t    t = (vf_ - vi_) / a_    (1.3)

    // (1) x component:     vfx = vix + ax t      (1.x)
    // (1) y component:     vfy = viy + ay t      (1.y)

    // (1.x) /. (3)     vfx = vi cos(th) + ax t     (1.x.1)
    // (1.y) /. (4)     vfy = vi sin(th) + ay t     (1.y.1)

    //////////////////////////////////////////////////////////////////////
    
    //          rf_ = ri_ + vi_ t + 1/2 a_ t^2             (2)

    // (2) vi_   vi_ = (rf_ - ri_ - 1/2 a_ t^2) / t        (2.1)

    // (2) a_    a_ = (rf_ - ri_ - vi_ t) / (1/2 t^2)      (2.2)

    //////////////////////////////////////////////////////////////////////

    // (2) x components:    xf = xi + vxi t + 1/2 ax t^2   (2.3)

    //////////////////////////////////////////////////////////////////////

    // (2) y components:    yf = yi + vyi t + 1/2 ay t^2   (2.4)

    // (2.4) t:     0 = 1/2 ay t^2  +  vyi t  +  yi - yf

    //              t = (- vyi +- Sqrt(vyi^2 - 4 * 1/2 ay (yi - yf))) / (2 * 1/2 ay)

    //              t = (- vyi +- Sqrt(vyi^2 - 2 ay (yi - yf))) / ay    (2.5)

    //////////////////////////////////////////////////////////////////////

    #region // (2.3) and (2.4) in terms of speed and direction

    // vxi = vi cos(th)   (3)

    // vyi = vi sin(th)   (4)

    // (2.3) /. (3)     xf = xi + vi cos(th) t + 1/2 ax t^2     (2.3.1)

    // (2.4) /. (4)     yf = yi + vi sin(th) t + 1/2 ay t^2     (2.4.1)

    #endregion

    //////////////////////////////////////////////////////////////////////

    #region // projectile motion - level plane - initial angle or speed known
    
    // projectile motion (ax = 0)

    // (2.3.1):     xf = xi + vi cos(th) t          (2.3.2)

    //////////////////////////////////////////////////////////////////////
    
    // projectile motion : level plane (yi = yf)

    // (2.4.1):     0 = vi sin(th) t + 1/2 ay t^2

    //              - vi sin(th) t = 1/2 ay t^2

    //              - vi sin(th) = 1/2 ay t           (2.4.2)

    //////////////////////////////////////////////////////////////////////

    // (2.3.2) t:    t = (xf - xi) / (vi cos(th))      (2.3.2.1)

    // (2.4.2) /. (2.3.2.1)

    //      - vi sin(th) = 1/2 ay (xf - xi) / (vi cos(th))

    //      - vi sin(th) cos(th) = 1/2 ay (xf - xi) / vi

    //      - vi^2 sin(th) cos(th) = 1/2 ay (xf - xi)

    //      vi^2 sin(th) cos(th) = - 1/2 ay (xf - xi)

    //      vi^2 2 sin(th) cos(th) = - ay (xf - xi)

    //      vi^2 sin(2 th) = - ay (xf - xi)             (2.4.2.1)

    //////////////////////////////////////////////////////////////////////

    // (2.4.2.1) th:    sin(2 th) = - ay (xf - xi) / vi^2

    //                  2 th = arcsin(- ay (xf - xi) / vi^2)

    //                  th = 1/2 arcsin(- ay (xf - xi) / vi^2)      (2.4.2.1.1)

    //////////////////////////////////////////////////////////////////////

    // (2.4.2.1) vi:    vi^2 = - ay (xf - xi) / sin(2 th)

    //                  vi = +- sqrt(- ay (xf - xi) / sin(2 th))    (2.4.2.1.2)

    #endregion

    //////////////////////////////////////////////////////////////////////

    #region // projectile motion incline

    // xB = xA + vAx t + 1/2 ax t^2
    // yB = yA + vAy t + 1/2 ay t^2

    // xB - xA = d cos(phi)
    // yB - yA = d sin(phi)


    // ax = 0
    // xA = 0

    // xB = xA + vAx t

    // xB - xA = vAx t

    // d cos(phi) = vAx t

    // cos(phi) = vAx t / d

    // cos(phi) / vAx = t / d

    // t = d cos(phi) / vAx


    // yB - yA = vAy t + 1/2 ay t^2

    // d sin(phi) = vAy t + 1/2 ay t^2

    // sin(phi) = vAy t / d + 1/2 ay t^2 / d

    // sin(phi) = vAy cos(phi) / vAx + 1/2 ay t^2 / d

    // sin(phi) = vAy cos(phi) / vAx + 1/2 ay (d cos(phi) / vAx)^2 / d

    // sin(phi) = vAy cos(phi) / vAx + 1/2 ay d cos^2(phi) / vAx^2

    // d = {sin(phi) - vAy / vAx cos(phi)} / {1/2 ay cos^2(phi) / vAx^2} 

    #endregion

    //////////////////////////////////////////////////////////////////////

    class Obj
    {
        public static bool NotNull(params object[] objects)
        {
            for (int i = 0; i < objects.Length; i++)
                if (objects[i] == null) return false;

            return true;
        }

        static double Discriminant(double a, double b, double c)
        {
            return b*b - 4*a*c;
        }

        public static double QuadraticA(double a, double b, double c)
        {
            return (-b + Math.Sqrt(b * b - 4 * a * c)) / (2 * a);
        }

        public static double QuadraticB(double a, double b, double c)
        {
            return (-b - Math.Sqrt(b * b - 4 * a * c)) / (2 * a);
        }

        public Point position, velocity, acceleration;
        public double? time;

        public double? angle; // constraint

        public double? speed; // constraint

        public double Speed { get { return speed.Value; } }
        public double Angle { get { return angle.Value; } }

        public override string ToString()
        {
            return string.Format("Obj(pos={0}, vel={1}, acc={2}, time={3:F3})",
                position, velocity, acceleration, time);
        }

        public Obj AtTime(double t)
        {
            var obj = new Obj();

            obj.time = t;

            var dt = (double) (t - time);

            obj.acceleration = acceleration;
            obj.velocity = velocity + acceleration * dt;
            obj.position = position + velocity * dt + 0.5 * acceleration * dt * dt;

            return obj;
        }

        public static Point CalcAcceleration(Obj a, Obj b)
        {
            var dt = (double)(b.time - a.time);

            if (b.velocity != null && a.velocity != null)
                return (b.velocity - a.velocity) / dt;

            if (b.position != null && a.position != null && a.velocity != null)
                return (b.position - a.position - a.velocity * dt) / (1 / 2 * dt * dt);

            throw new Exception();
        }

        public static double CalcTime(Obj a, Obj b, int solution=0)
        {
            if (b.velocity.x != null && a.velocity.x != null && a.acceleration.x != null && a.acceleration.x != 0.0)
                return (double) ((b.velocity.x - a.velocity.x) / a.acceleration.x);

            if (b.velocity.y != null && a.velocity.y != null && a.acceleration.y != null && a.acceleration.y != 0.0)
                return (double) ((b.velocity.y - a.velocity.y) / a.acceleration.y);

            Func<double, double, double, double> quadratic;

            if (solution == 0)
                quadratic = QuadraticA;
            else
                quadratic = QuadraticB;

            if (a.position.x != null &&
                b.position.x != null &&
                a.velocity.x != null &&
                a.acceleration.x != null &&
                a.acceleration.x != 0.0)
                return
                    quadratic(
                        (double)(0.5 * a.acceleration.x),
                        (double)a.velocity.x,
                        (double)(a.position.x - b.position.x));

            if (a.position.y != null &&
                b.position.y != null &&
                a.velocity.y != null &&
                a.acceleration.y != null &&
                a.acceleration.y != 0.0)
                return 
                    quadratic(
                        (double) (0.5 * a.acceleration.y), 
                        (double) a.velocity.y, 
                        (double) (a.position.y - b.position.y));

            throw new Exception();
        }

        public static List<double> CalcTimes(Obj a, Obj b)
        {
            var result = new List<double>();

            if (a.position.y != null &&
                b.position.y != null &&
                a.velocity.y != null &&
                a.acceleration.y != null)
            {
                var discriminant =
                    Discriminant(
                        (double)(0.5 * a.acceleration.y),
                        (double)a.velocity.y,
                        (double)(a.position.y - b.position.y));

                if (discriminant > 0)
                {
                    result.Add(
                        QuadraticA(
                            (double)(0.5 * a.acceleration.y),
                            (double)a.velocity.y,
                            (double)(a.position.y - b.position.y)));

                    result.Add(
                        QuadraticB(
                            (double)(0.5 * a.acceleration.y),
                            (double)a.velocity.y,
                            (double)(a.position.y - b.position.y)));

                    return result;
                }

                if (discriminant == 0.0)
                {
                    result.Add(
                        QuadraticA(
                            (double)(0.5 * a.acceleration.y),
                            (double)a.velocity.y,
                            (double)(a.position.y - b.position.y)));
                    return result;
                }

                throw new Exception("no real values");
            }

            throw new Exception();
        }

        public static Point CalcInitialVelocity(Obj a, Obj b, int solution=0)
        {
            if (NotNull(a.time, b.time))
            {
                var dt = (double)(b.time - a.time);

                if (NotNull(b.velocity.x, b.velocity.y, b.acceleration.x, b.acceleration.y))
                    return b.velocity - b.acceleration * dt;

                if (NotNull(b.position.x, b.position.y, b.acceleration.x, b.acceleration.y))
                    return (b.position - a.position - 0.5 * b.acceleration * dt * dt) / dt;
            }

            if (a.acceleration.x == 0 &&
                b.acceleration.x == 0 &&
                a.position.y == b.position.y)
            {
                if (a.speed != null)
                {
                    a.angle = 0.5 * Math.Asin(-a.acceleration.Y * (b.position.X - a.position.X) / a.Speed.Pow(2));

                    return Point.FromAngle(a.Angle, a.Speed);
                }

                if (a.angle != null)
                {
                    a.speed = Math.Sqrt(-a.acceleration.Y * (b.position.X - a.position.X) / Math.Sin(2 * a.Angle));

                    return Point.FromAngle(a.Angle, a.Speed);
                }
            }

            if (a.position.x == 0 &&
                a.position.y == 0 &&
                b.position.x.HasValue &&
                b.position.y.HasValue &&
                a.acceleration.x == 0 &&
                b.acceleration.x == 0 &&
                a.acceleration.y.HasValue)
            {
                if (a.speed.HasValue)
                {
                    var xf = b.position.X;
                    var yf = b.position.Y;
                    var vi = a.Speed;
                    var ay = a.acceleration.Y;
                    
                    if (solution == 0)
                        a.angle = (Math.Asin((yf - ay * xf.Pow(2) / vi.Pow(2)) / Math.Sqrt(xf.Pow(2) + yf.Pow(2))) + Math.Atan(yf / xf)) / 2.0;

                    if (solution == 1)
                        a.angle = (Math.PI - Math.Asin((yf - ay * xf.Pow(2) / vi.Pow(2)) / Math.Sqrt(xf.Pow(2) + yf.Pow(2))) + Math.Atan(yf / xf)) / 2.0;

                    return Point.FromAngle(a.Angle, a.Speed);
                }
            }

            throw new Exception();
        }

        //public double CalcRange()
        //{
        //    return
        //        Math.Pow(velocity.Norm(), 2) *
        //        Math.Sin(2 * velocity.ToAngle()) /
        //        Math.Abs((double)acceleration.y);
        //}
    
        public Point ProjectileInclineIntersection(double phi)
        {
            if (NotNull(phi, velocity.x, velocity.y, acceleration.y)
                &&
                acceleration.x == 0.0)
            {
                var d =
                    (phi.Sin() - velocity.y / velocity.x * phi.Cos())
                    /
                    (0.5 * acceleration.y * phi.Cos().Pow(2) / ((double)velocity.x).Pow(2));

                return
                    new Point(
                        position.x + d * phi.Cos(),
                        position.y + d * phi.Sin());
            }

            throw new Exception();
        }

        // ar = v^2 / r

        public double CalcRadialAcceleration(double radius)
        { return velocity.Norm().Pow(2) / radius; }
    }

    class Program
    {
        // A cannon with a muzzle speed of 1 000 m/s is used to
        // start an avalanche on a mountain slope. The target is 
        // 2 000 m from the cannon horizontally and 800 m above
        // the cannon. At what angle, above the horizontal, should
        // the cannon be fired?

        // PSE PRO 4.17

        // xi =    0
        // yi =    0
        // xf = 2000
        // yf =  800

        // vi = 1000

        // ax =  0
        // ay = -9.8

        
        // (2.3.1): xf = xi + vi cos(th) t + 1/2 ax t^2
        //          xf = vi cos(th) t
        // t        t = xf / (vi cos(th))       (2.3.1.1)

        // (2.4.1):     yf = yi + vi sin(th) t + 1/2 ay t^2
        //              yf = vi sin(th) t + 1/2 ay t^2
        // /. (2.3.1.1) yf = vi sin(th) {xf / (vi cos(th))} + 1/2 ay {xf / (vi cos(th))}^2
        //              yf = sin(th) xf / cos(th) + 1/2 ay xf^2 / (vi^2 cos^2(th))
        // * 2 cos^2(th)    2 cos^2(th) yf = 2 cos^2(th) sin(th) xf / cos(th) + 2 cos^2(th) 1/2 ay xf^2 / (vi^2 cos^2(th))
        //                  2 cos^2(th) yf = 2 cos(th) sin(th) xf + ay xf^2 / vi^2
        // power-reducing / half angle identity     cos^2(x) = [1 + cos(2x)]/2 :
        //                  2 [1 + cos(2 th)]/2 yf = 2 cos(th) sin(th) xf + ay xf^2 / vi^2
        //                  [1 + cos(2 th)] yf = 2 cos(th) sin(th) xf + ay xf^2 / vi^2
        // double angle formula 2 sin(x) cos(x) = sin(2x) :
        //                  [1 + cos(2 th)] yf = sin(2 th) xf + ay xf^2 / vi^2
        //                  yf + yf cos(2 th) = sin(2 th) xf + ay xf^2 / vi^2
        //                  sin(2 th) xf - yf cos(2 th) = yf - ay xf^2 / vi^2
        // xf = r cos(phi)
        // yf = r sin(phi)
        //                  sin(2 th) r cos(phi) - r sin(phi) cos(2 th) = yf - ay xf^2 / vi^2
        //                  r [ sin(2 th) cos(phi) - sin(phi) cos(2 th) ] = yf - ay xf^2 / vi^2
        // sum/difference identity  sin(x) cos(y) - sin(y) cos(x) = sin(x - y) :
        //                  r sin(2 th - phi) = yf - ay xf^2 / vi^2
        // r = sqrt(xf^2 + yf^2)
        //                  sqrt(xf^2 + yf^2) sin(2 th - phi) = yf - ay xf^2 / vi^2
        // tan(phi) = yf / xf       phi = arctan(yf / xf):
        //                  sqrt(xf^2 + yf^2) sin(2 th - arctan(yf / xf)) = yf - ay xf^2 / vi^2
        // sin(2 th - arctan(yf / xf)) = [yf - ay xf^2 / vi^2] / sqrt(xf^2 + yf^2)
        //
        // arcsin 1: 2 th - arctan(yf / xf) = arcsin{[yf - ay xf^2 / vi^2] / sqrt(xf^2 + yf^2)}         (arcsin1)
        //
        // and
        //
        // arcsin 2: 2 th - arctan(yf / xf) = PI - arcsin{[yf - ay xf^2 / vi^2] / sqrt(xf^2 + yf^2)}    (arcsin2)

        // (arcsin1):
        //
        // 2 th - arctan(yf / xf) = arcsin{[yf - ay xf^2 / vi^2] / sqrt(xf^2 + yf^2)}
        // 2 th = arcsin{[yf - ay xf^2 / vi^2] / sqrt(xf^2 + yf^2)} + arctan(yf / xf)
        // th = {arcsin{[yf - ay xf^2 / vi^2] / sqrt(xf^2 + yf^2)} + arctan(yf / xf)} / 2

        // (arcsin2):
        //
        // 2 th - arctan(yf / xf) = PI - arcsin{[yf - ay xf^2 / vi^2] / sqrt(xf^2 + yf^2)}
        // 2 th = PI - arcsin{[yf - ay xf^2 / vi^2] / sqrt(xf^2 + yf^2)} + arctan(yf / xf)
        // th = [PI - arcsin{[yf - ay xf^2 / vi^2] / sqrt(xf^2 + yf^2)} + arctan(yf / xf)] / 2

        static void Main(string[] args)
        {
            var obj_A = new Obj();

            obj_A.position = new Point(0, 0);
            obj_A.velocity = new Point(null, null);
            obj_A.acceleration = new Point(0, -9.8);

            obj_A.speed = 1000.0;

            var obj_B = new Obj();

            obj_B.position = new Point(2000.0, 800.0);
            obj_B.velocity = new Point(null, null);
            obj_B.acceleration = new Point(0, -9.8);

            Console.WriteLine(Obj.CalcInitialVelocity(obj_A, obj_B, solution: 0).ToAngle().Degrees());
            Console.WriteLine(Obj.CalcInitialVelocity(obj_A, obj_B, solution: 1).ToAngle().Degrees());
                        
            Console.ReadLine();
        }
    }
}