Skip to content

@dharmatech /pse-pro-4.17.cs
Created

Embed URL

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
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();
}
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.