Instantly share code, notes, and snippets.

# dharmatech/pse-pro-4.17.cs Created Feb 7, 2012

What would you like to do?
 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 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 CalcTimes(Obj a, Obj b) { var result = new List(); 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(); } } }