TTensor

• centerMass
• density
• density:
• mass
• mass:
• t0
• t1
• t2
• tensor
• tp

• calcPlane:
• faceIntegrate:tri:
• projectionIntegrate:
• result
• volumeIntegrate

• addFaces:vertices:
• initialize
• initializeDensity:
• initializeMass:

• initialize
• initializeDensity:
• initializeMass:
• testCube
• testCube:position:
• testSphere
• testSphere:
• testSphere:position:

	notDone ifTrue:[self result.].
	^ centerMass.

	^ density

	density _ d.
	mass _ nil.
	notDone _ true.

	^ mass.

	mass _ ms.
	density _ nil.
	notDone _ true.

	^ t0.

	^ t1.

	^ t2.

	notDone ifTrue:[ ^ self result.].
	^ tensor.

	^ tp.

	| abc d abcd |
	abc _ (((tri at: 2)-(tri at:1)) cross: ((tri at: 3)-(tri at: 1))).
	abc squaredLength = 0 ifTrue:[^ nil ].
	abc normalize.
	d _ (abc dot: (tri at: 1)) negated.
	abcd _ abc asB3DVector4.
	abcd w: d.
	^ abcd.

	| fArray pArray k1 k2 k3 k4 a1 a2 a3 b1 b2 b3 ab w p1 pa paa paaa pb pbb pbbb pab paab pabb|

	pArray _ self projectionIntegrate: tri.
	fArray _ FloatArray ofSize: 12.
	p1 _ pArray at:1.
	pa _ pArray at:2.
	paa _ pArray at:3.
	paaa _ pArray at:4.
	pb _ pArray at: 5.
	pbb _ pArray at: 6.
	pbbb _ pArray at: 7.
	pab _ pArray at: 8.
	paab _ pArray at: 9.
	pabb _ pArray at: 10.

	k1 _ 1.0/(plane at: c). k2 _ k1*k1. k3 _ k2*k1. k4 _ k3*k1.
	
	a1 _ plane at: a.
	a2 _ a1*a1.
	a3 _ a1*a2.
	b1 _ plane at: b.
	b2 _ b1*b1.
	b3 _ b1*b2.
	ab _ a1*b1.
	w _ plane at: 4.

	fArray at: 1 put: k1 * pa.
	fArray at: 2 put: k1 * pb.
	fArray at: 3 put: k2 *((a1*pa) + (b1*pb) + (w * p1)) negated .
	
	fArray at: 4 put: k1 * paa.
	fArray at: 5 put: k1 * pbb.
	fArray at: 6 put: k3 * ((a2*paa) + (2*ab*pab) + (b2*pbb) +
		(w*((2*((a1*pa)+(b1*pb))) + (w*p1)))).

	fArray at: 7 put: k1 * paaa.
	fArray at: 8 put: k1 * pbbb.
	fArray at: 9 put: k4 *((a3*paaa) + (3*a2*b1*paab)
		+ (3*a1*b2*pabb) + (b3*pbbb))negated.
	fArray at: 9 put: (fArray at: 9) +
		(k4*3*w*((a2*paa) + (2*a1*b1*pab) + (b2*pbb))negated).
	fArray at: 9 put: (fArray at: 9) +
		(k4*w*w*((3*((a1*pa) + (b1*pb))) + (w*p1)) negated).

	fArray at: 10 put: k1 * paab.
	fArray at: 11 put: k2  * ((a1*pabb) + (b1*pbbb) + (w*pbb)) negated.
	fArray at: 12 put: k3 * ((a2*paaa) + (2*a1*b1*paab) + (b2*pabb)).
	fArray at: 12 put: (fArray at: 12) + (k3 * (w*((2*((a1*paa) + (b1*pab))) + (w*pa)))).

	^  fArray.
	

"void
InertiaTensor::FaceIntegrals(TVVector3 norm, float d, TVVector3 tri[3])
{
	double w;
	double k1, k2, k3, k4;

	ProjectionIntegrals(tri);

	w = d;
	k1 = 1 / norm[C]; k2 = k1 * k1; k3 = k2 * k1; k4 = k3 * k1;

1	Fa = k1 * Pa;
2	Fb = k1 * Pb;
3	Fc = -k2 * (norm[A]*Pa + norm[B]*Pb + w*P1);

4	Faa = k1 * Paa;
5	Fbb = k1 * Pbb;
6	Fcc = k3 * (sqr(norm[A])*Paa + 2*norm[A]*norm[B]*Pab + sqr(norm[B])*Pbb
		+ w*(2*(norm[A]*Pa + norm[B]*Pb) + w*P1));

7	Faaa = k1 * Paaa;
8	Fbbb = k1 * Pbbb;
9	Fccc = -k4 * (cube(norm[A])*Paaa + 3*sqr(norm[A])*norm[B]*Paab 
		+ 3*norm[A]*sqr(norm[B])*Pabb + cube(norm[B])*Pbbb
		+ 3*w*(sqr(norm[A])*Paa + 2*norm[A]*norm[B]*Pab + sqr(norm[B])*Pbb)
		+ w*w*(3*(norm[A]*Pa + norm[B]*Pb) + w*P1));

10	Faab = k1 * Paab;
11	Fbbc = -k2 * (norm[A]*Pabb + norm[B]*Pbbb + w*Pbb);
12	Fcca = k3 * (sqr(norm[A])*Paaa + 2*norm[A]*norm[B]*Paab + sqr(norm[B])*Pabb
		+ w*(2*(norm[A]*Paa + norm[B]*Pab) + w*Pa));
}

"

	| pArray a0 a1 b0 b1 da db a02 a03 a04 b02 b03 b04 a12 a13 b12 b13 c1 ca caa caaa cb cbb cbbb cab kab caab kaab cabb kabb |
	
	pArray _ FloatArray ofSize:10.

	1 to: 3 do:[ :i || i2 |
		i2 _ 1+(i\\3).
		a0 _ (tri at: i) at: a.
		b0 _ (tri at: i) at: b.
		a1 _ (tri at: i2) at: a.
		b1 _ (tri at: i2) at: b.
		da _ a1 - a0.
		db _ b1 - b0.
		a02 _ a0*a0. a03 _ a02*a0. a04 _ a03*a0.
		b02 _ b0*b0. b03 _ b02*b0. b04 _ b03*b0. 
		a12 _ a1*a1. a13 _ a12*a1.
		b12 _ b1*b1. b13 _ b12*b1.

		c1 _ a1+a0.
		ca _ (a1*c1) + a02. 		caa _ (a1*ca)+a03. caaa _ (a1*caa)+a04.
		cb _ (b1*(b1+b0)) + b02. cbb _ (b1*cb)+b03. cbbb _ (b1*cbb)+b04.
		cab _ (3*a12)+(2*a1*a0)+a02. kab _ a12+(2*a1*a0)+(3*a02).
		caab _ (a0*cab)+(4*a13). kaab _ (a1*kab)+(4*a03).
		cabb _ (4*b13)+(3*b12*b0)+(2*b1*b02)+b03.
		kabb _ b13+(2*b12*b0)+(3*b1*b02)+(4*b03).

		pArray at: 1 put: ((pArray at: 1) + (db*c1)).
		pArray at: 2 put: ((pArray at: 2) + (db*ca)).
		pArray at: 3 put: ((pArray at: 3) + (db*caa)).
		pArray at: 4 put: ((pArray at: 4) + (db*caaa)).
		pArray at: 5 put: ((pArray at: 5) + (da*cb)).
		pArray at: 6 put: ((pArray at: 6) + (da* cbb)).
		pArray at: 7 put: ((pArray at: 7) + (da*cbbb)).
		pArray at: 8 put: ((pArray at: 8) + (db*((b1*cab)+(b0*kab)))).
		pArray at: 9 put: ((pArray at: 9) + (db*((b1*caab)+(b0*kaab)))).
		pArray at: 10 put: ((pArray at: 10) + (da*((a1*cabb)+(a0*kabb)))).
		].

	pArray at: 1 put: (pArray at: 1)/2.0.
	pArray at: 2 put: (pArray at: 2)/6.0.
	pArray at: 3 put: (pArray at: 3)/12.0.
	pArray at: 4 put: (pArray at: 4)/20.0.
	pArray at: 5 put: (pArray at: 5)/-6.0.
	pArray at: 6 put: (pArray at: 6)/-12.0.
	pArray at: 7 put: (pArray at: 7)/-20.0.
	pArray at: 8 put: (pArray at: 8)/24.0.
	pArray at: 9 put: (pArray at: 9)/60.0.
	pArray at:10 put: (pArray at: 10)/-60.0.

	^ pArray.
"
void 
InertiaTensor::ProjectionIntegrals(TVVector3 tri[3])
{
	double a0, a1, da;
	double b0, b1, db;
	double a0_2, a0_3, a0_4, b0_2, b0_3, b0_4;
	double a1_2, a1_3, b1_2, b1_3;
	double C1, Ca, Caa, Caaa, Cb, Cbb, Cbbb;
	double Cab, Kab, Caab, Kaab, Cabb, Kabb;

	P1 = Pa = Pb = Paa = Pab = Pbb = Paaa = Paab = Pabb = Pbbb = 0.0;

	for (int i = 0; i < 3; i++) {
		a0 = tri[i][A];
		b0 = tri[i][B];
		a1 = tri[(i+1) % 3][A];
		b1 = tri[(i+1) % 3][B];
		da = a1 - a0;
		db = b1 - b0;
		a0_2 = a0 * a0; a0_3 = a0_2 * a0; a0_4 = a0_3 * a0;
		b0_2 = b0 * b0; b0_3 = b0_2 * b0; b0_4 = b0_3 * b0;
		a1_2 = a1 * a1; a1_3 = a1_2 * a1; 
		b1_2 = b1 * b1; b1_3 = b1_2 * b1;

		C1 = a1 + a0;
		Ca = a1*C1 + a0_2; Caa = a1*Ca + a0_3; Caaa = a1*Caa + a0_4;
		Cb = b1*(b1 + b0) + b0_2; Cbb = b1*Cb + b0_3; Cbbb = b1*Cbb + b0_4;
		Cab = 3*a1_2 + 2*a1*a0 + a0_2; Kab = a1_2 + 2*a1*a0 + 3*a0_2;
		Caab = a0*Cab + 4*a1_3; Kaab = a1*Kab + 4*a0_3;
		Cabb = 4*b1_3 + 3*b1_2*b0 + 2*b1*b0_2 + b0_3;
		Kabb = b1_3 + 2*b1_2*b0 + 3*b1*b0_2 + 4*b0_3;

		P1 += db*C1;
		Pa += db*Ca;
		Paa += db*Caa;
		Paaa += db*Caaa;
		Pb += da*Cb;
		Pbb += da*Cbb;
		Pbbb += da*Cbbb;
		Pab += db*(b1*Cab + b0*Kab);
		Paab += db*(b1*Caab + b0*Kaab);
		Pabb += da*(a1*Cabb + a0*Kabb);
	}

1	P1 /= 2.0;
2	Pa /= 6.0;
3	Paa /= 12.0;
4	Paaa /= 20.0;
5	Pb /= -6.0;
6	Pbb /= -12.0;
7	Pbbb /= -20.0;
8	Pab /= 24.0;
9	Paab /= 60.0;
10	Pabb /= -60.0;
}"

	| cms tt1 tt2 ttp|

	tt1 _ t1/2.0.
	tt2 _ t2/3.0.
	ttp _ tp/2.0.
	
	mass ifNotNil:[
		density _ mass/t0.] ifNil:[
		mass _ density*t0.].

	centerMass _ tt1/t0.
	cms _ centerMass * centerMass.
"Compute the inertia tensor."
	tensor _ B3DMatrix4x4 new.
	tensor a11: (density *((tt2 at:2) + (tt2 at:3))) - (mass * ((cms at: 2) + (cms at: 3))).
	tensor a22: (density *((tt2 at:3) + (tt2 at:1))) - (mass * ((cms at: 3) + (cms at: 1))).
	tensor a33: (density *((tt2 at:1) + (tt2 at:2))) - (mass * ((cms at: 1) + (cms at: 2))).
	tensor a12: (mass * (centerMass at:1) * (centerMass at: 2)) - (density * (ttp at: 1)). 
	tensor a21: tensor a12.
	tensor a23:  (mass * (centerMass at:2) * (centerMass at: 3)) - (density * (ttp at: 2)). 
	tensor a32: tensor a23.
	tensor a13:  (mass * (centerMass at:1) * (centerMass at: 3)) - (density * (ttp at: 3)) . 
	tensor a31: tensor a13.
	tensor a44: 1.0.
	notDone _ false.
	^ tensor.


"




 
"

"InertiaTensor::InertiaTensor(RMFace* faces, int numFaces, OSPoint3D* points, float denseMass, ITFlags flags):
mFaces(faces), mNumFaces(numFaces), mPoints(points)
{
	TVMatrix33 it;

	VolumeIntegrals();

	if(flags == ITDENSITY)
	{
		mDensity = denseMass;
		mMass = mDensity * T0;
	}
	else if(flags == ITMASS)
	{
		mMass = denseMass;
		mDensity = mMass/T0;
	}

	//TVVector3 mCM;
	//TVMatrix33 mJ;
	 /* compute center of mass */
	mCM.x = T1[X] / T0;
	mCM.y = T1[Y] / T0;
	mCM.z = T1[Z] / T0;

	/* compute inertia tensor */
	mJ[X][X] = mDensity * (T2[Y] + T2[Z]);
	mJ[Y][Y] = mDensity * (T2[Z] + T2[X]);
	mJ[Z][Z] = mDensity * (T2[X] + T2[Y]);
	mJ[X][Y] = mJ[Y][X] = - mDensity * TP[X];
	mJ[Y][Z] = mJ[Z][Y] = - mDensity * TP[Y];
	mJ[Z][X] = mJ[X][Z] = - mDensity * TP[Z];

	/* translate inertia tensor to center of mass */
	mJ[X][X] -= mMass * (mCM.y*mCM.y + mCM.z*mCM.z);
	mJ[Y][Y] -= mMass * (mCM.z*mCM.z + mCM.x*mCM.x);
	mJ[Z][Z] -= mMass * (mCM.x*mCM.x + mCM.y*mCM.y);
	mJ[X][Y] = mJ[Y][X] += mMass * mCM.x * mCM.y; 
	mJ[Y][Z] = mJ[Z][Y] += mMass * mCM.y * mCM.z; 
	mJ[Z][X] = mJ[X][Z] += mMass * mCM.z * mCM.x; 
}
"

	| tri plane aPlane fArray t |
	tri _ B3DVector3Array ofSize:3.
	1 to: faces size by: 3 do:[ :i |
		tri at: 1 put: (vertices at: 1+(faces at: i)).
		tri at: 2 put: (vertices at: 1+(faces at: i+1)).
		tri at: 3 put: (vertices at: 1+(faces at: i+2)).
		plane _ self calcPlane: tri.
		plane ifNotNil:[
			aPlane _ plane abs.
	" a,b, and c are the indices into the plane equation s.t. c is the maximal length component."
			c _ 1.
			(aPlane at: 2) > (aPlane at: c) ifTrue:[ c _ 2. ].
			(aPlane at: 3) > (aPlane at: c) ifTrue:[ c _ 3 ].

			a _ 1+(c\\3).
			b _ 1+(a\\3).
			fArray _ self faceIntegrate: plane tri: tri.
			a == 1 ifTrue:[ t _ fArray at: 1.].
			b == 1 ifTrue:[ t _ fArray at: 2.].
			c == 1 ifTrue:[ t _ fArray at: 3.].
			t0 _ t0 + ((plane at: 1) * t).
			t1 at: a put: (t1 at: a) + ((plane at:a)*(fArray at:4)).
			t1 at: b put: (t1 at: b) + ((plane at:b)*(fArray at:5)).
			t1 at: c put: (t1 at: c) + ((plane at:c)*(fArray at:6)).
			t2 at: a put: (t2 at: a) + ((plane at:a)*(fArray at:7)).
			t2 at: b put: (t2 at: b) + ((plane at:b)*(fArray at:8)).
			t2 at: c put: (t2 at: c) + ((plane at:c)*(fArray at:9)).
			tp at: a put: (tp at: a) + ((plane at:a)*(fArray at:10)).
			tp at: b put: (tp at: b) + ((plane at:b)*(fArray at:11)).
			tp at: c put: (tp at: c) + ((plane at:c)*(fArray at: 12)).
			].
		].

"void 
InertiaTensor::VolumeIntegrals()
{
	double nx, ny, nz;

	T0 = T1[X] = T1[Y] = T1[Z] 
		= T2[X] = T2[Y] = T2[Z] 
		= TP[X] = TP[Y] = TP[Z] = 0;

	for (int i = 0; i < mNumFaces; i++) 
	{
		TVVector3 tri[3];
		TVVector3 norm;
		float d;
		tri[0] = mPoints[mFaces[i].P1.nVIndex];
		tri[1] = mPoints[mFaces[i].P2.nVIndex];
		tri[2] = mPoints[mFaces[i].P3.nVIndex];

		if(CalcPlaneNorm(tri, norm, d))
		{
			nx = fabs(norm[X]);
			ny = fabs(norm[Y]);
			nz = fabs(norm[Z]);
			if (nx > ny && nx > nz) C = X;
			else C = (ny > nz) ? Y : Z;
			A = (C + 1) % 3;
			B = (A + 1) % 3;
			FaceIntegrals(norm, d, tri);

			T0 += norm[X] * ((A == X) ? Fa : ((B == X) ? Fb : Fc));

			T1[A] += norm[A] * Faa;
			T1[B] += norm[B] * Fbb;
			T1[C] += norm[C] * Fcc;
			T2[A] += norm[A] * Faaa;
			T2[B] += norm[B] * Fbbb;
			T2[C] += norm[C] * Fccc;
			TP[A] += norm[A] * Faab;
			TP[B] += norm[B] * Fbbc;
			TP[C] += norm[C] * Fcca;
		}
	}

	T1[X] /= 2; T1[Y] /= 2; T1[Z] /= 2;
	T2[X] /= 3; T2[Y] /= 3; T2[Z] /= 3;
	TP[X] /= 2; TP[Y] /= 2; TP[Z] /= 2;
}
"

	faces _ fcs.
	vertices _ verts.
	self volumeIntegrate.
	notDone _ true.

	mass _ 1.0.
	t0 _ 0.0.
	t1 _ B3DVector3 new.
	t2 _ B3DVector3 new.
	tp _ B3DVector3 new.
	notDone _ true.

	density _ d.
	t0 _ 0.0.
	t1 _ B3DVector3 new.
	t2 _ B3DVector3 new.
	tp _ B3DVector3 new.
	notDone _ true.

	mass _ mss.
	t0 _ 0.0.
	t1 _ B3DVector3 new.
	t2 _ B3DVector3 new.
	tp _ B3DVector3 new.
	notDone _ true.

	^ self new initialize.

	^ self new initializeDensity: d.

	^ self new initializeMass: m.

	
	| cube quadFaces triFaces it counter |

	cube _ TCube new.
	cube location: (B3DVector3 x:1 y:2 z: 3).
	cube update.
	quadFaces _ cube quadFaces.
	triFaces _ IntegerArray ofSize: (3*2*6).
	counter _ 1.
	1 to: quadFaces size by: 4 do:[ :i |
		triFaces at: counter put: (quadFaces at: i)-1.
		triFaces at: counter + 1 put: (quadFaces at: i+1)-1.
		triFaces at: counter + 2 put: (quadFaces at: i+2)-1.
		triFaces at: counter + 3 put: (quadFaces at: i+2)-1.
		triFaces at: counter + 4 put: (quadFaces at: i+3)-1.
		triFaces at: counter + 5 put: (quadFaces at: i)-1.
		counter _ counter+6.].
	it _ TTensor createNewWithMass: 1.0.
	it addFaces: triFaces vertices: cube vertices.
	it result.
	Transcript show: 'tensor: '; cr; show: it tensor; cr.
	Transcript show: 'center of mass: '; cr; show: it centerMass; cr.
	Transcript show: 'density: '; show: it density; cr.
	^ it.
	

	

	
	| cube quadFaces triFaces it counter |

	cube _ TCube new.
	cube location: position.
	cube extent: scale.
	cube update.
	quadFaces _ cube quadFaces.
	triFaces _ IntegerArray ofSize: (3*2*6).
	counter _ 1.
	1 to: quadFaces size by: 4 do:[ :i |
		triFaces at: counter put: (quadFaces at: i)-1.
		triFaces at: counter + 1 put: (quadFaces at: i+1)-1.
		triFaces at: counter + 2 put: (quadFaces at: i+2)-1.
		triFaces at: counter + 3 put: (quadFaces at: i+2)-1.
		triFaces at: counter + 4 put: (quadFaces at: i+3)-1.
		triFaces at: counter + 5 put: (quadFaces at: i)-1.
		counter _ counter+6.].
	it _ TTensor createNewWithMass: 1.0.
	it addFaces: triFaces vertices: cube vertices.
	it result.
	Transcript show: 'tensor: '; cr; show: it tensor; cr.
	Transcript show: 'center of mass: '; cr; show: it centerMass; cr.
	Transcript show: 'density: '; show: it density; cr.
	^ it.
	

	

	
	| v f vert radius segments position seg2 ringSin ringCos rts rtc rbs rbc pi2 ax it |

	radius _ 1.0.
	segments _30.
	seg2 _ 1+ (segments * 2) .
	ringSin _ FloatArray ofSize: seg2.
	ringCos _ FloatArray ofSize: seg2.
	vert _ B3DVector3 new.
	pi2 _ Float pi *2.0.
	position _ B3DVector3 x: 10 y:10 z: 10.
	"position _ B3DVector3 new."
	v _ OrderedCollection new.
	f _ OrderedCollection new.
	1 to: seg2-1 do:[ :index | 
		ax _ ((index-1) * pi2)/ (seg2-1).
		ringSin at:index put: ax sin.
		ringCos at:index put: ax cos.].
	ringSin at: seg2 put: (ringSin at: 1).
	ringCos at: seg2 put: (ringCos at: 1).

	rts _ 0.0.
	rtc _ 1.0.
	rbs _ ringSin at: 2.
	rbc _ ringCos at: 2.
	1 to: segments do:[ :iv |
		1 to: seg2 do: [ :ih | 
			vert x: rts*(ringSin at: ih) y: rtc z: rts*(ringCos at: ih).
			v add:  (vert * radius)+ position.
			vert x: rbs*(ringSin at: ih) y: rbc z: rbs*(ringCos at: ih).
			v add: (vert * radius)+ position.
			ih  = 1 ifFalse:[
				f add: v size - 4.
				f add: v size - 3.
				f add: v size - 2.
				f add: v size - 3.
				f add: v size - 1.
				f add: v size - 2.
			].].
		rts _ rbs.
		rtc _ rbc.
		rbs _ ringSin at: iv+2.
		rbc _ ringCos at: iv+2.
		].
"	1 to: f size by: 3 do:[ : i | |p1 p2 center scale aForm|
		aForm _ Form extent: 1@1.
		aForm fillBlack.
		center _ 500@500.
		scale _ 100@100.
		p1 _ ((v at: 1+(f at: i)) at:2)@((v at: 1+(f at: i))at:3).
		p2 _ ((v at: 1+(f at: i+1)) at:2)@((v at: 1+(f at: i+1))at:3).
		(Line from: (p1*scale)+center to: (p2*scale)+center withForm: aForm) display.
		p1 _ p2.
		p2 _ ((v at: 1+(f at: i+2)) at:2)@((v at: 1+(f at: i+2))at: 3).
		(Line from: (p1*scale)+center to: (p2*scale)+center withForm: aForm) display.
		p1 _ p2.
		p2 _((v at: 1+(f at: i)) at:2)@((v at: 1+(f at: i))at:3)..
		(Line from: (p1*scale)+center to: (p2*scale)+center withForm: aForm)display.
		]."
	it _ TTensor initializeMass: 1.0.
	it addFaces: f vertices: v.
	it result.

	Transcript show:'Sphere tensor:'; cr.
	Transcript show: it tensor; cr.
	Transcript show:'Center of mass:'; cr.
	Transcript show: it centerMass; cr.
	Transcript show: 'Density: '; cr; show: it density; cr.
	^ it.
	

	

	
	| v f vert radius position seg2 ringSin ringCos rts rtc rbs rbc pi2 ax it |

	radius _ 1.0.
	seg2 _ 1+ (segments * 2) .
	ringSin _ FloatArray ofSize: seg2.
	ringCos _ FloatArray ofSize: seg2.
	vert _ B3DVector3 new.
	pi2 _ Float pi *2.0.
	position _ B3DVector3 x: 10 y:10 z: 10.
	"position _ B3DVector3 new."
	v _ OrderedCollection new.
	f _ OrderedCollection new.
	1 to: seg2-1 do:[ :index | 
		ax _ ((index-1) * pi2)/ (seg2-1).
		ringSin at:index put: ax sin.
		ringCos at:index put: ax cos.].
	ringSin at: seg2 put: (ringSin at: 1).
	ringCos at: seg2 put: (ringCos at: 1).

	rts _ 0.0.
	rtc _ 1.0.
	rbs _ ringSin at: 2.
	rbc _ ringCos at: 2.
	1 to: segments do:[ :iv |
		1 to: seg2 do: [ :ih | 
			vert x: rts*(ringSin at: ih) y: rtc z: rts*(ringCos at: ih).
			v add:  (vert * radius)+ position.
			vert x: rbs*(ringSin at: ih) y: rbc z: rbs*(ringCos at: ih).
			v add: (vert * radius)+ position.
			ih  = 1 ifFalse:[
				f add: v size - 4.
				f add: v size - 3.
				f add: v size - 2.
				f add: v size - 3.
				f add: v size - 1.
				f add: v size - 2.
			].].
		rts _ rbs.
		rtc _ rbc.
		rbs _ ringSin at: iv+2.
		rbc _ ringCos at: iv+2.
		].
"	1 to: f size by: 3 do:[ : i | |p1 p2 center scale aForm|
		aForm _ Form extent: 1@1.
		aForm fillBlack.
		center _ 500@500.
		scale _ 100@100.
		p1 _ ((v at: 1+(f at: i)) at:2)@((v at: 1+(f at: i))at:3).
		p2 _ ((v at: 1+(f at: i+1)) at:2)@((v at: 1+(f at: i+1))at:3).
		(Line from: (p1*scale)+center to: (p2*scale)+center withForm: aForm) display.
		p1 _ p2.
		p2 _ ((v at: 1+(f at: i+2)) at:2)@((v at: 1+(f at: i+2))at: 3).
		(Line from: (p1*scale)+center to: (p2*scale)+center withForm: aForm) display.
		p1 _ p2.
		p2 _((v at: 1+(f at: i)) at:2)@((v at: 1+(f at: i))at:3)..
		(Line from: (p1*scale)+center to: (p2*scale)+center withForm: aForm)display.
		]."
	it _ TTensor initializeMass: 1.0.
	it addFaces: f vertices: v.
	it result.

	Transcript show:'Sphere tensor:'; cr.
	Transcript show: it tensor; cr.
	Transcript show:'Center of mass:'; cr.
	Transcript show: it centerMass; cr.
	Transcript show: 'Density: '; cr; show: it density; cr.
	^ it.
	

	

	
	| v f vert radius seg2 ringSin ringCos rts rtc rbs rbc pi2 ax it |

	radius _ 1.0.
	seg2 _ 1+ (segments * 2) .
	ringSin _ FloatArray ofSize: seg2.
	ringCos _ FloatArray ofSize: seg2.
	vert _ B3DVector3 new.
	pi2 _ Float pi *2.0.
	v _ OrderedCollection new.
	f _ OrderedCollection new.
	1 to: seg2-1 do:[ :index | 
		ax _ ((index-1) * pi2)/ (seg2-1).
		ringSin at:index put: ax sin.
		ringCos at:index put: ax cos.].
	ringSin at: seg2 put: (ringSin at: 1).
	ringCos at: seg2 put: (ringCos at: 1).

	rts _ 0.0.
	rtc _ 1.0.
	rbs _ ringSin at: 2.
	rbc _ ringCos at: 2.
	1 to: segments do:[ :iv |
		1 to: seg2 do: [ :ih | 
			vert x: rts*(ringSin at: ih) y: rtc z: rts*(ringCos at: ih).
			v add:  (vert * radius)+ position.
			vert x: rbs*(ringSin at: ih) y: rbc z: rbs*(ringCos at: ih).
			v add: (vert * radius)+ position.
			ih  = 1 ifFalse:[
				f add: v size - 4.
				f add: v size - 3.
				f add: v size - 2.
				f add: v size - 3.
				f add: v size - 1.
				f add: v size - 2.
			].].
		rts _ rbs.
		rtc _ rbc.
		rbs _ ringSin at: iv+2.
		rbc _ ringCos at: iv+2.
		].
"	1 to: f size by: 3 do:[ : i | |p1 p2 center scale aForm|
		aForm _ Form extent: 1@1.
		aForm fillBlack.
		center _ 500@500.
		scale _ 100@100.
		p1 _ ((v at: 1+(f at: i)) at:2)@((v at: 1+(f at: i))at:3).
		p2 _ ((v at: 1+(f at: i+1)) at:2)@((v at: 1+(f at: i+1))at:3).
		(Line from: (p1*scale)+center to: (p2*scale)+center withForm: aForm) display.
		p1 _ p2.
		p2 _ ((v at: 1+(f at: i+2)) at:2)@((v at: 1+(f at: i+2))at: 3).
		(Line from: (p1*scale)+center to: (p2*scale)+center withForm: aForm) display.
		p1 _ p2.
		p2 _((v at: 1+(f at: i)) at:2)@((v at: 1+(f at: i))at:3)..
		(Line from: (p1*scale)+center to: (p2*scale)+center withForm: aForm)display.
		]."
	it _ TTensor initializeMass: 1.0.
	it addFaces: f vertices: v.
	it result.

	Transcript show:'Sphere tensor:'; cr.
	Transcript show: it tensor; cr.
	Transcript show:'Center of mass:'; cr.
	Transcript show: it centerMass; cr.
	Transcript show: 'Density: '; cr; show: it density; cr.
	^ it.