using System;
namespace Science.Physics.GeneralPhysics
{
///
/// RigidBody
///
public class RigidBody
{
public RigidBody()
{
}
private double mass;
public RigidBody(Mass m)
{
mass = m.kg;
}
private MomentOfInertia I = new MomentOfInertia();
public void SetCylinder(Length radius, Length height)
{
I.XX = 1.0/12.0*mass*height.m*height.m
+1.0/4.0*mass*radius.m*radius.m;
I.YY = I.XX;
I.ZZ = 1.0/2.0*mass*radius.m*radius.m;
}
public void SetEllipsoid(Length a, Length b, Length c)
{
I.XX = 1.0/5.0*mass*(b.m*b.m+c.m*c.m);
I.YY = 1.0/5.0*mass*(a.m*a.m+c.m*c.m);
I.ZZ = 1.0/5.0*mass*(a.m*a.m+b.m*b.m);
}
public void SetEllipticalSlab(Length a, Length b, Length height)
{
I.XX = 1.0/6.0*mass*(3.0*b.m*b.m+4.0*height.m*height.m);
I.YY = 1.0/5.0*mass*(3.0*b.m*b.m+4.0*height.m*height.m);
I.ZZ = 1.0/2.0*mass*(a.m*a.m+b.m*b.m);
}
public void SetRectangularParallelepiped(Length a, Length b, Length c)
{
I.XX = 1.0/3.0*mass*(b.m*b.m+c.m*c.m);
I.YY = 1.0/3.0*mass*(a.m*a.m+c.m*c.m);
I.ZZ = 1.0/3.0*mass*(a.m*a.m+b.m*b.m);
}
public void SetRing(Length radius)
{
I.XX = 1.0/2.0*mass*radius.m*radius.m;
I.YY = 1.0/2.0*mass*radius.m*radius.m;
I.ZZ = mass*radius.m*radius.m;
}
public void SetRodAboutCenter(Length height)
{
I.XX = 1.0/12.0*mass*height.m*height.m;
I.YY = 1.0/12.0*mass*height.m*height.m;
I.ZZ = 0.0;
}
public void SetRodAboutEnd(Length height)
{
I.XX = 1.0/3.0*mass*height.m*height.m;
I.YY = 1.0/3.0*mass*height.m*height.m;
I.ZZ = 0.0;
}
public void SetSphere(Length radius)
{
I.XX = 2.0/5.0*mass*radius.m*radius.m;
I.YY = 2.0/5.0*mass*radius.m*radius.m;
I.ZZ = 2.0/5.0*mass*radius.m*radius.m;
}
public void SetSphericalShell(Length radius)
{
I.XX = 2.0/3.0*mass*radius.m*radius.m;
I.YY = 2.0/3.0*mass*radius.m*radius.m;
I.ZZ = 2.0/3.0*mass*radius.m*radius.m;
}
public void SetTorus(Length radius, Length crossSectionRadius)
{
I.XX = 1.0/8.0*mass*(4.0*radius.m*radius.m
+5.0*crossSectionRadius.m*crossSectionRadius.m);
I.YY = 1.0/3.0*mass*(4.0*radius.m*radius.m
+5.0*crossSectionRadius.m*crossSectionRadius.m);
I.ZZ = 1.0/3.0*mass*(radius.m*radius.m
+3.0/4.0*crossSectionRadius.m*crossSectionRadius.m);
}
public MomentOfInertia RigidBodyMomentOfInertia
{
get{return I;}
}
private ConstraintFunctionToBeZero[] cf;
private int cfnum = 0;
public delegate double ConstraintFunctionToBeZero(Position[] r, Force[] f);
public void Constraint(ConstraintFunctionToBeZero[] func)
{
cf = func;
cfnum = func.Length;
}
private Position[] rr;
private Force[] ff;
public void SolveStaticEquilibrium(Position[] r, Force[] f)
{
rr = r;
ff = f;
int ndim = 0;
for(int k = 0; k < r.Length; k++)
{
if(r[k].XVariableQ) ndim += 1;
if(r[k].YVariableQ) ndim += 1;
if(r[k].ZVariableQ) ndim += 1;
}
for(int k = 0; k < f.Length; k++)
{
if(f[k].XVariableQ) ndim += 1;
if(f[k].YVariableQ) ndim += 1;
if(f[k].ZVariableQ) ndim += 1;
}
Calculus.Function fun
= new Calculus.Function(Func);
double[] c = new double[ndim];
double[] ans = Calculus.MinimumOfFunction(fun, c, 100.0);
}
private double Func(double[] x)
{
int i = 0;
for(int k = 0; k < rr.Length; k++)
{
if(rr[k].XVariableQ) this.rr[k].X = x[i++];
if(rr[k].YVariableQ) this.rr[k].Y = x[i++];
if(rr[k].ZVariableQ) this.rr[k].Z = x[i++];
}
for(int k = 0; k < ff.Length; k++)
{
if(ff[k].XVariableQ) this.ff[k].X = x[i++];
if(ff[k].YVariableQ) this.ff[k].Y = x[i++];
if(ff[k].ZVariableQ) this.ff[k].Z = x[i++];
}
TotalForce tf = new TotalForce(ff);
Torque[] torque = new Torque[ff.Length];
for(int k = 0; k < this.ff.Length; k++)
torque[k] = new Torque(rr[k],ff[k]);
TotalTorque tt = new TotalTorque(torque);
double sum = 0.0;
for(int k = 0; k < cfnum; k++)
{
sum += cf[k](this.rr,this.ff)*cf[k](this.rr,this.ff);
}
return tf.N*tf.N+tt.J*tt.J + sum;
}
}
}