pse-pro-4.19.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 Time { get { return time.Value; } }

        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 (NotNull(a.position.x, b.position.x, a.velocity.x) && a.acceleration.x == 0.0)
                return (b.position.X - a.position.X) / a.velocity.X;

            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 placekicker must kick a football from a point 36.0 m
        // (about 40 yards) from the goal, and half the crowd
        // hopes the ball will clear the crossbar, which is 3.05 m
        // high. When kicked, the ball leaves the ground with a
        // speed of 20.0 m/s at an angle of 53.0° to the horizontal.
        //
        // (a) By how much does the ball clear or fall short of
        // clearing the crossbar?
        //
        // (b) Does the ball approach the crossbar while 
        // still rising or while falling?

        // PSE PRO 4.19

        // xi = 0   
        // yi = 0

        // vix = 20 cos(53 deg)
        // viy = 20 sin(53 deg)

        // ax = 0
        // ay = -9.8

        // xf = 36
        // yf =  ?

        // t = ?

        // (2.3)        xf = xi + vxi t + 1/2 ax t^2
        // ax=0         xf = xi + vxi t
        //              t = (xf - xi) / vxi

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

            obj_A.position = new Point(0, 0);
            obj_A.velocity = Point.FromAngle((53.0).Radians(), 20);
            obj_A.acceleration = new Point(0, -9.8);
            obj_A.time = 0;

            var obj_B = new Obj();

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

            obj_B.time = Obj.CalcTime(obj_A, obj_B);

            obj_B = obj_A.AtTime(obj_B.Time);

            Console.WriteLine("object at A: " + obj_A);
            Console.WriteLine("object at B: " + obj_B);

            var field_goal_height = 3.05;

            Console.WriteLine("ball clears goal by {0:F3} meters",
                (obj_B.position.y - field_goal_height));

            Console.WriteLine("ball approaches the goal while {0}",
                obj_B.velocity.y > 0 ? "rising" : "falling");

            Console.ReadLine();
        }
    }
}