Category : C Source Code
Archive : GRAFGEM3.ZIP
Filename : SQFINAL.C
Output of file : SQFINAL.C contained in archive : GRAFGEM3.ZIP
* This code calculates the volume, mass, and inertia tensors of
* superquadric ellipsoids and toroids. The code includes methods
* to numerically compute gamma functions and beta functions
*/
#include
#include
/*
* The following function, abgam() is based on a continued fraction numerical
* method found in Abremowitz and Stegun, Handbook of Mathematical Functions
*
*/
double abgam (x)
double x;
{
double gam[10],
result,
temp;
gam[0] = 1./ 12.;
gam[1] = 1./ 30.;
gam[2] = 53./ 210.;
gam[3] = 195./ 371.;
gam[4] = 22999./ 22737.;
gam[5] = 29944523./ 19733142.;
gam[6] = 109535241009./ 48264275462.;
temp = 0.5*log (2*M_PI) - x + (x - 0.5)*log (x)
+ gam[0]/(x + gam[1]/(x + gam[2]/(x + gam[3]/(x + gam[4] /
(x + gam[5]/(x + gam[6]/x))))));
return temp;
}
/*
* A method to compute the gamma() function.
*
*/
double gamma (x)
double x;
{
double result,
abgam ();
result = exp (abgam (x + 5))/(x*(x + 1)*(x + 2)*(x + 3)*(x + 4));
return result;
}
/*
* A method to compute the beta() function.
*/
double beta (m, n)
double m,
n;
{
double gamma ();
return (gamma (m)*gamma (n)/gamma (m + n));
}
/*
* A method to compute the volume of a superquadric ellipsoid
* with axis lengths a1, a2, a3, north-south exponent n
* and east-west exponent e
*/
double sqellipvol (a1, a2, a3, n, e)
double a1, /* x radius */
a2, /* y radius */
a3, /* z radius */
n, /* north-south param */
e; /* east-west param */
{
double beta ();
return ((2./ 3.)*a1*a2*a3*e*n*beta (e/2., e/2.)*beta (n, n/2.));
}
/*
* A method to compute the volume of a superquadric toroid
* with axis lengths a1, a2, a3, north-south exponent n
* east-west exponent e, and hole parameter alpha
*/
double sqtoroidvol (a1, a2, a3, n, e, alpha)
double a1, /* x radius */
a2, /* y radius */
a3, /* z radius */
n, /* north-south param */
e, /* east-west param */
alpha; /* torus hole-size parameter */
{
return (2.*a1*a2*a3*alpha*e*n*beta (e/2., e/2.)*beta (n/2., n/2.));
}
/*
* A procedure to print the inertia tensor of a canonical superquadric
* ellipsoid in its body coordinate system.
*/
void sq_ellipsoid_tensor (rho, a1, a2, a3, e, n)
double a1,
a2,
a3,
e,
n,
rho; /* density of material */
{
double iellip[3][3],
i1E,
i2E,
i3E,
iworld[3][3],
R[3][3];
i1E = (2./ 5.)* a1*a1*a1*a2*a3*e*n*beta(3.* e/2., e/2.)*beta(2.* n, n/2.);
i2E = (2./ 5.)* a1*a2*a2*a2*a3*e*n*beta(e/2., 3.* e/2.)*beta(2.* n, n/2.);
i3E = (2./ 5.)* a1*a2*a3*a3*a3*e*n*beta(e/2., e/2.)*beta(n, 3.* n/2.);
iellip[0][0] = 0;
iellip[1][0] = 0;
iellip[2][0] = 0;
iellip[0][1] = 0;
iellip[1][1] = 0;
iellip[2][1] = 0;
iellip[0][2] = 0;
iellip[1][2] = 0;
iellip[2][2] = 0;
iellip[0][0] = i2E + i3E;
iellip[1][1] = i1E + i3E;
iellip[2][2] = i1E + i2E;
printf ("ellipsoid inertia tensor in body coordinates\n");
printf (" iellip1 = %f %f %f \n", iellip[0][0], iellip[1][0], iellip[2][0]);
printf (" iellip2 = %f %f %f \n", iellip[0][1], iellip[1][1], iellip[2][1]);
printf (" iellip3 = %f %f %f \n", iellip[0][2], iellip[1][2], iellip[2][2]);
}
/*
* A procedure to print the inertia tensor of a canonical superquadric
* toroid in its body coordinate system.
*/
void sq_toroid_tensor (rho, a1, a2, a3, e, n, alpha)
double rho,
a1,
a2,
a3,
e,
n,
alpha;
{
double itor[3][3],
i1T,
i2T,
i3T;
i1T = a1*a1*a1*a2*a3*alpha*e*n*beta(3*e/2,e/2)*(2*alpha*alpha*beta(n/2,n/2)+
3*beta(3*n/2,n/2));
i2T = a1*a2*a2*a2*a3*alpha*e*n*beta(e/2,3*e/2)*(2*alpha*alpha*beta(n/2,n/2)+
3*beta(3*n/2,n/2));
i3T = a1*a2*a3*a3*a3*alpha*e*n*beta(e/2,e/2)*beta(n/2,3*n/2);
itor[0][0] = 0;
itor[1][0] = 0;
itor[2][0] = 0;
itor[0][1] = 0;
itor[1][1] = 0;
itor[2][1] = 0;
itor[0][2] = 0;
itor[1][2] = 0;
itor[2][2] = 0;
itor[0][0] = i2T + i3T;
itor[1][1] = i1T + i3T;
itor[2][2] = i1T + i2T;
printf ("toroid inertia tensor in body coordinates\n");
printf (" itor1 = %f %f %f \n", itor[0][0], itor[1][0], itor[2][0]);
printf (" itor2 = %f %f %f \n", itor[0][1], itor[1][1], itor[2][1]);
printf (" itor3 = %f %f %f \n", itor[0][2], itor[1][2], itor[2][2]);
}
/*
* A procedure to print the inertia tensor components in world coordinates,
* given an inertia tensor in body coordinates, and the 3x3 rotation matrix
* which rotates body vectors into world coordinates.
*/
void iworld (Ibody, R)
double Ibody[3][3],
R[3][3];
{
double Iworld[3][3];
Iworld[0][0] =
(R[0][0]*Ibody[0][0]+R[0][1]*Ibody[1][0]+R[0][2]*Ibody[2][0])*R[0][0] +
(R[0][0]*Ibody[0][1]+R[0][1]*Ibody[1][1]+R[0][2]*Ibody[2][1])*R[0][1] +
(R[0][0]*Ibody[0][2]+R[0][1]*Ibody[1][2]+R[0][2]*Ibody[2][2])*R[0][2];
Iworld[1][0] =
(R[1][0]*Ibody[0][0]+R[1][1]*Ibody[1][0]+R[1][2]*Ibody[2][0])*R[0][0]+
(R[1][0]*Ibody[0][1]+R[1][1]*Ibody[1][1]+R[1][2]*Ibody[2][1])*R[0][1]+
(R[1][0]*Ibody[0][2]+R[1][1]*Ibody[1][2]+R[1][2]*Ibody[2][2])*R[0][2];
Iworld[2][0] =
(R[2][0]*Ibody[0][0]+R[2][1]*Ibody[1][0]+R[2][2]*Ibody[2][0])*R[0][0]+
(R[2][0]*Ibody[0][1]+R[2][1]*Ibody[1][1]+R[2][2]*Ibody[2][1])*R[0][1]+
(R[2][0]*Ibody[0][2]+R[2][1]*Ibody[1][2]+R[2][2]*Ibody[2][2])*R[0][2];
Iworld[0][1] =
(R[0][0]*Ibody[0][0]+R[0][1]*Ibody[1][0]+R[0][2]*Ibody[2][0])*R[1][0]+
(R[0][0]*Ibody[0][1]+R[0][1]*Ibody[1][1]+R[0][2]*Ibody[2][1])*R[1][1]+
(R[0][0]*Ibody[0][2]+R[0][1]*Ibody[1][2]+R[0][2]*Ibody[2][2])*R[1][2];
Iworld[1][1] =
(R[1][0]*Ibody[0][0]+R[1][1]*Ibody[1][0]+R[1][2]*Ibody[2][0])*R[1][0]+
(R[1][0]*Ibody[0][1]+R[1][1]*Ibody[1][1]+R[1][2]*Ibody[2][1])*R[1][1]+
(R[1][0]*Ibody[0][2]+R[1][1]*Ibody[1][2]+R[1][2]*Ibody[2][2])*R[1][2];
Iworld[2][1] =
(R[2][0]*Ibody[0][0]+R[2][1]*Ibody[1][0]+R[2][2]*Ibody[2][0])*R[1][0]+
(R[2][0]*Ibody[0][1]+R[2][1]*Ibody[1][1]+R[2][2]*Ibody[2][1])*R[1][1]+
(R[2][0]*Ibody[0][2]+R[2][1]*Ibody[1][2]+R[2][2]*Ibody[2][2])*R[1][2];
Iworld[0][2] =
(R[0][0]*Ibody[0][0]+R[0][1]*Ibody[1][0]+R[0][2]*Ibody[2][0])*R[2][0]+
(R[0][0]*Ibody[0][1]+R[0][1]*Ibody[1][1]+R[0][2]*Ibody[2][1])*R[2][1]+
(R[0][0]*Ibody[0][2]+R[0][1]*Ibody[1][2]+R[0][2]*Ibody[2][2])*R[2][2];
Iworld[1][2] =
(R[1][0]*Ibody[0][0]+R[1][1]*Ibody[1][0]+R[1][2]*Ibody[2][0])*R[2][0]+
(R[1][0]*Ibody[0][1]+R[1][1]*Ibody[1][1]+R[1][2]*Ibody[2][1])*R[2][1]+
(R[1][0]*Ibody[0][2]+R[1][1]*Ibody[1][2]+R[1][2]*Ibody[2][2])*R[2][2];
Iworld[2][2] =
(R[2][0]*Ibody[0][0]+R[2][1]*Ibody[1][0]+R[2][2]*Ibody[2][0])*R[2][0]+
(R[2][0]*Ibody[0][1]+R[2][1]*Ibody[1][1]+R[2][2]*Ibody[2][1])*R[2][1]+
(R[2][0]*Ibody[0][2]+R[2][1]*Ibody[1][2]+R[2][2]*Ibody[2][2])*R[2][2];
printf ("toroid inertia tensor in body coordinates\n");
printf (" Iworld1 = %f %f %f \n", Iworld[0][0], Iworld[1][0], Iworld[2][0]);
printf (" Iworld2 = %f %f %f \n", Iworld[0][1], Iworld[1][1], Iworld[2][1]);
printf (" Iworld3 = %f %f %f \n", Iworld[0][2], Iworld[1][2], Iworld[2][2]);
}
/*
* sgn(x) returns -1.0 or 1.0 for speed.
* sgn(x) can return 0.0 on zero if you wish, depending on
* your convention
*/
double sgn(x)
double x;
{
if (x <= 0.0) return (-1.0);
else return(1.0);
}
/*
* computes position on the surface of a superquadric ellipsoid
* v goes from -Pi/2 to Pi/2; u goes from -Pi to Pi.
*
*/
void sqellipsoidposn(a1,a2,a3,n,e,alpha,u,v)
double a1,a2,a3,n,e,alpha,u,v;
{
double cu, su, cv, sv, x,y,z;
cu = cos(u);
su = sin(u);
cv = cos(v);
sv = cos(v);
x = a1*(alpha + pow(cv,n))*pow(cu,e)*sgn(cu)*sgn(cv);
y = a2*(alpha + pow(cv,n))*pow(su,e)*sgn(su)*sgn(cv);
z = a3*pow(sv,n)*sgn(sv);
}
/*
* computes position on the surface of a superquadric toroid
* u and v go from -Pi to Pi
*/
void sqtoroidposn(a1,a2,a3,n,e,u,v)
double a1,a2,a3,n,e,u,v;
{
double cu, su, cv, sv, x,y,z;
cu = cos(u);
su = sin(u);
cv = cos(v);
sv = cos(v);
x = a1*pow(cv,n)*pow(cu,e)*sgn(cu)*sgn(cv);
y = a2*pow(cv,n)*pow(su,e)*sgn(su)*sgn(cv);
z = a3*pow(sv,n)*sgn(sv);
}
/*
* A procedure to test some of the above code
*
*/
main () {
double x;
printf (" gamma(1)= 1.0 = %12.10lf\n", gamma (1.0));
printf (" gamma(1/2)^2= Pi =%12.10lf\n", gamma (0.5)*gamma (0.5));
printf (" gamma(2)= 1.0 = %12.10lf\n", gamma (2.0));
printf (" gamma(3)= 2.0 = %12.10lf\n", gamma (3.0));
printf (" gamma(4)= 6.0 = %12.10lf\n", gamma (4.0));
printf("\n");
printf ("beta(1,1)= 1.0 = %12.10lf\n", beta (1.0, 1.0));
printf ("beta(1,1/2)= 2.0 = %12.10lf\n", beta (1.0, 0.5));
printf ("beta(1/2,1/2)= Pi = %12.10lf\n", beta (0.5, 0.5));
printf("\n");
printf ("sq ellipsoid volume/pi= 4/3 = %12.10lf\n",
sqellipvol(1., 1., 1., 1., 1.)/M_PI);
printf ("sq toroid volume/pi^2 = 2.0 = %12.10lf\n",
sqtoroidvol (1., 1., 1., 1., 1., 1.)/M_PI/M_PI);
printf("\n");
sq_ellipsoid_tensor (1., 1., 1., 1., 1., 1., 1.);
sq_toroid_tensor (1., 1., 1., 1., 1., 1., 1.);
}
Very nice! Thank you for this wonderful archive. I wonder why I found it only now. Long live the BBS file archives!
This is so awesome! 😀 I’d be cool if you could download an entire archive of this at once, though.
But one thing that puzzles me is the “mtswslnkmcjklsdlsbdmMICROSOFT” string. There is an article about it here. It is definitely worth a read: http://www.os2museum.com/wp/mtswslnk/