Simulation of VIRGO payload

Version 2.0

VERSIONS LOG BOOK

Generalities

General constants and settings

Needs["Graphics`Graphics`"];
Off[General::spell1];
Off[General::spell];

precis = 16;
SP  = SetPrecision[#, precis] &;
g = 980665/100000;
SteelYoungModulus = 186 * 10^9 *0.5; (*!!!*)
SteelPoissonRatio = 314/1000 *0.5;   (*!!!*)
SteelDensity    = 7874 *2300/7874;   (*!!!*)
AluminumDensity = 2698;
TitanumDensity  = 4540;
SeismicSpeed = 500;
PhiWire = 3 10^(-3) /1000;
PhiBlade  = 3 10^(-3);

FarTower = True;
BeamFiniteElements = 4;
BeamReductionFactor = 4;

BladeFiniteElements = 2;
BladeReductionFactor = 1;

Utility functions

Translation and rotation matrices at linear order

TranslateRotate[x_List, x0_List, t0_List] :=
   x0 + {{1,      -t0[[3]], +t0[[2]]},
         {+t0[[3]],       1, -t0[[1]]},
         {-t0[[2]], +t0[[1]],        1}} . x;

Fast Derivative

MyOuterD[expr_, sinistra_List, destra_List] :=
  Module[{ris, tmp},
         tmp = Expand[expr];
         ris = (D[tmp, #] &) /@ sinistra;
         Return[Outer[D[#1, #2] &, ris, destra]];
        ];

Numerical Values

Wires

L[6] = 110/100; 				r[6] =  1 1/1000;
L[7] = 73/100;  				r[7] =  3/20000;
L[8] = 73/100;
If[FarTower == True,
   r[8] = 150 10^(-6),
   r[8] = 150 10^(-6)];

Masses

M[6] = 80;
M[7] = 45;
If[FarTower == True,
   M[8] = 4233/100,
   M[8] = 2117/100];

Tensions

T[8] = g*M[8]/4;
T[7] = g*M[7]/4;
T[6] = g*(M[6] + M[7] + M[8]);

Momenta of Inertia

Jxx[6] = 2.5;
Jxy[6] = 0.01;
Jxz[6] = 0.0;
Jyy[6] = 0.95;
Jyz[6] = 0.0;
Jzz[6] = 2.96;

Jx[7] =        1.78134;
Jy[7] = Jz[7] = 1.25233;

Jx[8] =         1/2 M[8] (350/1000/2)^2;
Jy[8] = Jz[8] = M[8]/12*(3*(350/1000/2)^2 + (100/1000)^2);

Clamp positions (respect to CM)

Marionetta

clampNum[6]   = {0, 8};
clampEq[6, 0] = {0, 0, 0}; (* wire to the Cross of "Seven" *)

(* Wires holding the reference mass *)
clampEq[6, 1] = {+50/1000/2, +400/1000/2, 0}; 
clampEq[6, 2] = {-50/1000/2, +400/1000/2, 0};
clampEq[6, 3] = {-50/1000/2, -400/1000/2, 0};
clampEq[6, 4] = {+50/1000/2, -400/1000/2, 0};

(* Wires holding the mirror *)
clampEq[6, 5] = {+50/1000/2, +350/1000/2, 0};
clampEq[6, 6] = {-50/1000/2, +350/1000/2, 0};
clampEq[6, 7] = {-50/1000/2, -350/1000/2, 0};
clampEq[6, 8] = {+50/1000/2, -350/1000/2, 0};

Reference Mass

clampNum[7]   = {1,4};
clampEq[7, 1] = {+50/1000/2, +400/1000/2, 0};
clampEq[7, 2] = {-50/1000/2, +400/1000/2, 0};
clampEq[7, 3] = {-50/1000/2, -400/1000/2, 0};
clampEq[7, 4] = {+50/1000/2, -400/1000/2, 0};

Mirror

clampNum[8]   = {1,4};
clampEq[8, 1] = {+50/1000/2, +350/1000/2, 0};
clampEq[8, 2] = {-50/1000/2, +350/1000/2, 0};
clampEq[8, 3] = {-50/1000/2, -350/1000/2, 0};
clampEq[8, 4] = {+50/1000/2, -350/1000/2, 0};

Clamp position after CM rotations and translations

Module[{i,j},
       Do[clamp[i, j] =
              TranslateRotate[clampEq[i, j],
                              {x[i], y[i], z[i]},
                              {tx[i], ty[i], tz[i]}] - 
                              clampEq[i, j],
          {i,6,8}, {j,clampNum[i][[1]],clampNum[i][[2]]}]];

Formulas for wires

Wire Finite Element approximation

Mass matrix

MfemBeam[ρ_,  h_] := ((h * ρ)/420) (({{156, 22, 54, -13}, {22, 4, 13, -3}, {54, 13, 156, -22}, {-13, -3, -22, 4}}));

Stiffness matrix

KfemBeam[EI_,  h_,  T_] := (EI/(h^3))(({{12, 6, -12, 6}, {6, 4, -6, 2}, {-12, -6, 12, -6}, {6, 2, -6, 4}}))+(T/(30 h))(({{36, 3, -36, 3}, {3, 4, -3, -1}, {-36, -3, 36, -3}, {3, -1, -3, 4}}))

Generation of element lengths

\(\(BeamElementLengths[n_Integer, \ m_Integer]\  := \ 

        Module[{hn, \ nodes}, 
		If[OddQ[n], \ 
            Print["ElementLengthsBeam::error: ", \ n, \ "is odd"]; \ 
            Return[{}]]; 
		hn\  = \ 
            n/2\  - 1; 
		nodes\  = \ {1/2}; 
		If[m == 1, 
            Return[Table[1/n, {n}]]]; 
		While[
            hn\  > \ 0, 
			PrependTo[nodes, 
              nodes[\([1]\)]*1/m]; 
			\(hn--\)]; 
		PrependTo[nodes, 
            0]; 
		nodes\  = \ 
            Join[nodes, \ 1\  - \ Rest[Reverse[nodes]]]; 
		Return[
            Drop[nodes, 1]\  - \ Drop[nodes, \(-1\)]];];\)\)

Generation of potential and kinetic energy matrices

\(\(BeamMatrices[Xd_List, \ Td_List, \ Xu_List, \ Tu_List, ρ_, \ 
          L_, \ r_, \ T_, \ n_Integer, \ m_Integer]\  := 
	
        Module[{EI, \ h\ , hList, j\ , 
            regola, 
					\ \ \ \ \ \ \ \ \ \ \ \ Ukin = 0, Upot = 0, \ 
            xBeam, \ yBeam, \ txBeam, \ tyBeam, Kbase, \ Mbase, \ X, \ Y, \ 
            Z, \ KTorque, \ KTraction}, 
		EI\  = \ 
            SteelYoungModulus\ *\ \((1\  + \ 
                  I\ *\ PhiWire)\)*\((Pi\ *\ r^4)\)/4; 
		
		hList\  = 
            L*BeamElementLengths[n, m]; \ 
	
		 (*\ 
            notice\ that\ it\ is\ the\ linear\ density\ which\ counts\ *) 
		Kbase\  = \ \(KfemBeam[EI, #, T]\  &\)\  /@ \ 
              hList; 
	\ \ Mbase\  = \ \(MfemBeam[ρ*\((Pi*
                        r^2)\), #]\  &\)\  /@ \ hList; 
	
		Do[
            h\  = \ hList[\([j]\)]; 
			\ \ \ \ Y = {yBeam[
                  j - 1], \ \(-h\)*txBeam[j - 1], \ 
                yBeam[j], \ \(-h\)*txBeam[j]}; 
			\ \ \ \ X\  = \ {xBeam[
                  j - 1], \ h*tyBeam[j - 1], \ xBeam[j], \ 
                h*tyBeam[j]}; 
				\ \ \ \ Ukin\  = \ 
              Ukin\  + \ 1/2\ *\ Y\  . \ Mbase[\([j]\)]\  . \ Y\  + \ 
                1/2\ *\ X\  . \ Mbase[\([j]\)]\  . \ 
                    X; 
				\ \ \ \ Upot\  = \ 
              Upot\  + \ 1/2\ *\ Y\  . \ Kbase[\([j]\)]\  . \ Y\  + \ 
                1/2\ *\ X\  . \ Kbase[\([j]\)]\  . \ X, 
			\ \ \ \ {j, 
              1, n}]; 
		
		regola\  = \ {xBeam[0]\  -> \ Xu[\([1]\)], 
              tyBeam[0]\  -> \ Tu[\([2]\)], \ xBeam[n]\  -> \ Xd[\([1]\)], \ 
              tyBeam[n]\  -> \ Td[\([2]\)], 
				yBeam[0]\  -> \ 
                Xu[\([2]\)], \ txBeam[0]\  -> \ Tu[\([1]\)], \ 
              yBeam[n]\  -> \ Xd[\([2]\)], \ 
              txBeam[n]\  -> \ Td[\([1]\)]}; 
		
		varList\  = \ 
            Table[{xBeam[j], \ tyBeam[j], \ yBeam[j], \ txBeam[j]}, {j, 1, 
                  n - 1}] // Flatten; 
		
		Ukin\  = \ 
            Ukin /. regola; 
		Upot\  = \ 
            Upot /. \ regola; 
		
		 (*\ 
            add\ torsional\ and\ traction\ stiffness\ *) 
	
          KTorque\  = \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ SteelYoungModulus\ \((1\  + \ I\ PhiWire)\)/\((4\ Pi\ \((1\  + \ 
                        SteelPoissonRatio)\))\)\ *\ \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \((\((\ Pi\ *r^2\ )\)^2)\)/L; \ 
	
          Upot\  = \ 
            Upot\  + \ 
              1/2\ KTorque\ *\ \((\ 
                    Td[\([3]\)]\  - \ Tu[\([3]\)]\ )\)^2; \ 
		
	
          KTraction\  = \ 
            SteelYoungModulus\ *\ \((1\  + \ I\ PhiWire)\)\ *\ \((Pi\ *r^2)\)/
                L; \ 
\ \ Upot\  = \ 
            Upot\  + \ 
              1/2\ KTraction\ *\ \((\ 
                    Xd[\([3]\)]\  - \ 
                      Xu[\([3]\)]\ )\)^2; 
		
		Return[{Ukin, \ Upot, \ 
              varList}]];\)\)

Wire Potential Energy Formula

UWire[L_, r_, T_, Xd_List, Td_List, Xu_List, Tu_List] :=
  Module[{EI, k, mat, X, Y, KTorque, KTraction, ris},
         (* Torsional Stiffness *)
         KTorque =
             SteelYoungModulus (1 + I PhiWire)/(4 Pi (1 + SteelPoissonRatio)) * 
             (( Pi r^2 )^2)/L;
         ris = 1/2 KTorque * ( Td[[3]] - Tu[[3]] )^2;
         (* Traction Stiffness *)
         KTraction = SteelYoungModulus * (1 + I PhiWire) *
                     (Pi r^2)/L;
         ris += 1/2 KTraction * ( Xd[[3]] - Xu[[3]] )^2;
         (* Transverse stiffnesses *) 
         (* matrix between 

			X = x[i], L_i ty[i], x[i-1], L_i ty[i-1] 

			  or between

			Y =y[i], -L_i tx[i], y[i-1], -L_i tx[i-1]
		 *)
         EI = SteelYoungModulus * (1 + I PhiWire) *
               (Pi r^4) / 4;
         k = Sqrt[T L^2/(EI)];
         mat = T * L;
         mat = mat * k / ( L^2 * (k  - 2 Tanh[k/2]) );
         mat = mat *
            {{1, Tanh[k/2]/k, - 1, Tanh[k/2]/k},
             {Tanh[k/2]/k, 1/k (1/Tanh[k] - 1/k), 
            - Tanh[k/2]/k, 1/k (Tanh[k]/k - 1/Cosh[k])},
            {-1, - Tanh[k/2]/k, 1, - Tanh[k/2]/k},
            {Tanh[k/2]/k, 1/k (Tanh[k/2]/k - 1/Cosh[k]), 
             - Tanh[k/2]/k, 1/k (1/Tanh[k] - 1/k)}};
          mat = SP[mat];
          X = {Xd[[1]], L*Td[[2]],Xu[[1]], L*Tu[[2]]};
          Y = {Xd[[2]],-L*Td[[1]],Xu[[2]],-L*Tu[[1]]};
          ris += 1/2 X . mat . X + 1/2 Y . mat . Y;
          Return[ris];];

Lagrangian

Variables vectors

Q  =  {};
Q0  =  {x[5],  y[5],  z[5],  tx[5],  ty[5],  
          tz[5]};
Qphys  =  {};

Marionetta Energy

Q  =  
        Join[Q,  {x[6],  y[6],  z[6],  tx[6],  ty[6],  
            tz[6]}];
Pot[6]  =  Kin[6]  =  0;

Wire contribute

Module[{tmp},
		tmp = BeamMatrices[clamp[6,0], {tx[6], ty[6], tz[6]},
               				{x[5], y[5], z[5]},
               				{tx[5], ty[5], tz[5]},
		   					SteelDensity, L[6], r[6], T[6],
		   					BeamFiniteElements, BeamReductionFactor];
		Kin[6] += tmp[[1]];
		Pot[6] += tmp[[2]];
		Q = Join[Q, tmp[[3]]];];

CM and suspended raising

Pot[6] += 1/2 (tx[6]^2 + ty[6]^2)*
			(T[6] clamp[6,0][[3]] - 
			 T[7]*(clampEq[6,1] + clampEq[6,2]+
				  clampEq[6,3] + clampEq[6,4])[[3]] -
			 T[8]*(clampEq[6,5] + clampEq[6,6]+
				   clampEq[6,7] + clampEq[6,8])[[3]]);

Kinetic Energy

Kin[6] += 1/2 M[6] (x[6]^2 + y[6]^2 + z[6]^2) + 
	   1/2 (Jxx[6] tx[6]^2 + Jyy[6] ty[6]^2 + Jzz[6] tz[6]^2 +
	   	    2 Jxy[6] tx[6] ty[6] + 2 Jxz[6] tx[6] tz[6] + 
	   	    2 Jyz[6] ty[6] tz[6]);

Reference Mass Energy

Q  =  
        Join[Q,  {x[7],  y[7],  z[7],  tx[7],  ty[7],  
            tz[7]}];
Kin[7]  =  Pot[7]  =  0;

Wires contribution

Module[{tmp,i},
		Do[tmp = BeamMatrices[clamp[7,i], {tx[7], ty[7], tz[7]},
               				clamp[6,i], {tx[6], ty[6], tz[6]},
		   					SteelDensity, L[7], r[7], T[7],
		   					BeamFiniteElements, BeamReductionFactor];
		   Kin[7] += tmp[[1]];
		   Pot[7] += tmp[[2]];
		   Q = Join[Q, tmp[[3]]],
		   {i,1,4}];];

CM raising

Pot[7] += 1/2 (tx[7]^2 + ty[7]^2) T[7] *
		(clampEq[7,1] + clampEq[7,2] + 
		 clampEq[7, 3] + clampEq[7, 4])[[3]];

Kinetic energy

Kin[7] += 1/2 M[7] (x[7]^2 + y[7]^2 + z[7]^2) + 
	   1/2 (Jx[7] tx[7]^2 + Jy[7] ty[7]^2 + Jz[7] tz[7]^2);

Test Mass energy

Q  =  
        Join[Q,  {x[8],  y[8],  z[8],  tx[8],  ty[8],  
            tz[8]}];
Kin[8]  =  Pot[8]  =  0;

Wires contribution

Module[{tmp,i},
		Do[tmp = BeamMatrices[clamp[8,i], {tx[8], ty[8], tz[8]},
               				clamp[6,i+4], {tx[6], ty[6], tz[6]},
		   					SteelDensity, L[8], r[8], T[8],
		   					BeamFiniteElements, BeamReductionFactor];
		   Kin[8] += tmp[[1]];
		   Pot[8] += tmp[[2]];
		   Q = Join[Q, tmp[[3]]],
		   {i,1,4}];];

CM raising

Pot[8] += 1/2 (tx[8]^2 + ty[8]^2) T[8] *
		(clampEq[8,1] + clampEq[8,2] + 
		 clampEq[8, 3] + clampEq[8, 4])[[3]];

Kinetic energy

Kin[8] += 1/2 M[8] (x[8]^2 + y[8]^2 + z[8]^2) + 
	   1/2 (Jx[8] tx[8]^2 + Jy[8] ty[8]^2 + Jz[8] tz[8]^2);

Quadratic Lagrangian

Total energy

KinTot = Module[{i},Sum[Kin[i],{i,6,8}]];
PotTot = Module[{i},Sum[Pot[i],{i,6,8}]];

Vectors of variables for matrix formulation

Q = Union[Q];
Q0 = Union[Q0];
zeroRule = Join[Map[Rule[#, 0] &, Join[Q, Q0]], {force[i_] -> 0}];

Matrix computation

KinMat = SP[MyOuterD[KinTot, Q, Q]/. zeroRule];

PotMat = SP[MyOuterD[PotTot, Q, Q]/. zeroRule];

CoupMat= SP[-MyOuterD[PotTot,Q, Q0]/.zeroRule];

Transfer Matrix (direct solution)

TransferMatrix[nu_?NumberQ] :=
    Module[{ris,om},
           om = SP[2 Pi nu];
           ris = 
             LUBackSubstitution[LUDecomposition[-om^2 KinMat + PotMat],
           						CoupMat];
           Return[Transpose[ris]];
           ];

TransferMatrix::usage="TransferMatrix[nu_] builds numerically
the transfer matrix between the suspension point and all the
other degrees of freedom of the chain.

The rows of the matrix address the degrees of freedom of
the suspension point, in the same order as in the vector

	Q0 = {xSoil, tySoil, ySoil, txSoil, zSoil, tzSoil}

while the columns of the matrix address the degrees of
freedom of the chain, as in the vector Q. ";

Coils

Coils on the "Seventh" filter and magnets on Marionetta

Position of coils with respect to filter "seven" CM

(* RICONTROLLARE! *)
(* La posizione verticale e' data dalla lunghezza del
   filo di sospensione, diminuita della posizione verticale
   del clamp rispetto al CM e corretta per la spaziatura
   fra marionetta e gambe per il caso dei verticali *)
zCoilsFilterSeven = - L[6] + clampEq[5,1][[3]];
coilNum[5] = {1,4};
coil[5,1]  = {280/1000,0,zCoilsFilterSeven + 60/1000};
coil[5,2]  = {-280/1000,0,zCoilsFilterSeven + 60/1000};
coil[5,3]  = {-50/1000,391/1000,zCoilsFilterSeven};
coil[5,4]  = {50/1000,-391/1000,zCoilsFilterSeven};

Position of magnets with respect to Marionetta CM

(* RICONTROLLARE SE LA MARIONETTA E' DESTRORSA O
   SINISTRORSA! INFLUISCE SUL SEGNO DELLA COORD. X per i
   COILS *)
magnetNum[6] = {1,4};
magnet[6,1] = {280/1000,0,0};
magnet[6,2] = {-280/1000,0,0};
magnet[6,3] = {0,391/1000,0};
magnet[6,4] = {0,-391/1000,0};

Coils on the Reference Mass and magnets on the Test Mass

Position of coils with respect to Reference Mass CM

(* Assolutamente da RICONTROLLARE!! *)
coilNum[7] = {1,4};
coil[7,1] = {-60/1000, +350/2/1000, 0};
coil[7,2] = {-60/1000, -350/2/1000, 0};
coil[7,3] = {-60/1000, 0, +350/2/1000};
coil[7,4] = {-60/1000, 0, -350/2/1000};

Position of magnets with respect to Test Mass CM

(* Parecchio IPOTETICI !!! *)
magnetNum[8] = {1,4};
magnet[8,1] = {-50/1000, +350/2/1000, 0};
magnet[8,2] = {-50/1000, -350/2/1000, 0};
magnet[8,3] = {-50/1000, 0, +350/2/1000};
magnet[8,4] = {-50/1000, 0, -350/2/1000};

Coupling between coils ad magnets

Coil Forces

Naming of Forces for FilterSeven - Marionetta

Naming of Forces for ReferenceMass - TestMass

Forces acting vertically between filter "seven" and Marionetta

forceVector[5,1] = {0,0,-1/2(force[1] + force[2])};
forceVector[5,2] = {0,0,-1/2(force[1] - force[2])};

Forces acting horizontally between filter "seven" and Marionetta

forceVector[5,3] = { 1/2( force[3] - force[4]),0,0};
forceVector[5,4] = {-1/2(-force[3] - force[4]),0,0};

Forces acting horizontally between Reference and Test Mass

forceVector[7,1] = {1/2( 1/2(force[5] + force[7]) - force[6]),0,0};
forceVector[7,2] = {1/2( 1/2(force[5] + force[7]) + force[6]),0,0};
forceVector[7,3] = {1/2( 1/2(force[5] - force[7]) + force[8]),0,0};
forceVector[7,4] = {1/2( 1/2(force[5] - force[7]) - force[8]),0,0};

Definition of "Potential energy" from coil forces

Reaction on Filter Seven

Ext[5] = Module[{i},
				  Sum[-TranslateRotate[coil[5,i],
				  		{x[5], y[5], zVessel[5]},
				  		{tx[5], ty[5], tz[5]}] . (-forceVector[5,i]), 
				  	  {i,coilNum[5][[1]],coilNum[5][[2]]}]];

Action on Marionetta

Ext[6] = Module[{i},
				  Sum[-TranslateRotate[magnet[6,i],
				  		{x[6], y[6], z[6]},
				  		{tx[6], ty[6], tz[6]}] . forceVector[5,i], 
				  	  {i,magnetNum[6][[1]],magnetNum[6][[2]]}]];

Reaction on Reference Mass

Ext[7] = Module[{i},
				  Sum[-TranslateRotate[coil[7,i],
				  		{x[7], y[7], z[7]},
				  		{tx[7], ty[7], tz[7]}] . (-forceVector[7,i]), 
				  	  {i,coilNum[7][[1]],coilNum[7][[2]]}]];

Action on Test Mass

Ext[8] = Module[{i},
				Sum[-TranslateRotate[magnet[8,i],
				  	{x[8], y[8], z[8]},
				  	{tx[8], ty[8], tz[8]}] . forceVector[7,i], 
				  	{i,magnetNum[8][[1]],magnetNum[8][[2]]}]];

Total "External" potential energy (Dropping Filter Seven)

ExtTot = Module[{i},Sum[Ext[i],{i,6,8}]];

Coupling Matrix of forces

Forces = Module[{i}, Table[force[i],{i,1,8}]];
ExtMat = SP[-MyOuterD[ExtTot,Q, Forces]/.zeroRule];

Transfer Matrix Forces (direct solution)

TransferMatrixForces[nu_?NumberQ] :=
    Module[{ris, om},
           om = SP[2 Pi nu];
           tmpA = LUBackSubstitution[LUDecomposition[SP[-om^2 KinMat + PotMat]],
									 SP[ExtMat]];
           Return[Transpose[ris]];
           ];

TransferMatrixForces::usage="TransferMatrixForces[nu_] builds numerically
the transfer matrix between the internal forces in the chain and all the
other degrees of freedom of the chain.

The rows of the matrix address the internal forces,
in the same order as in the vector

	Forces = {force[1], ... force[8]}

while the columns of the matrix address the degrees of
freedom of the chain, as in the vector Q. ";

SaveTransferForces[nuMin_?NumberQ, nuMax_?NumberQ, nuPoints_Integer] :=
	Module[{nu, logNu, logNuMin, logNuStep,
          f1File, f2File, f3File, f4File, f5File, f6File, f7File,
          f1Stream, f2Stream, f3Stream, f4Stream, f5Stream, f6Stream, f7Stream,
          ris,
          xIndex, tyIndex, zIndex, tzIndex},
         logNuMin = Log[nuMin];
         logNuMax = Log[nuMax];
         logNuStep = (logNuMax - logNuMin)/(nuPoints - 1);
         
         (* define File Names *)
         f1File = StringJoin["f1.dat"];
         f2File = StringJoin["f2.dat"];
         f3File = StringJoin["f3.dat"];
         f4File = StringJoin["f4.dat"];
         f5File = StringJoin["f5.dat"];
         f6File = StringJoin["f6.dat"];
         f7File = StringJoin["f7.dat"];
         f8File = StringJoin["f8.dat"];
         
         (* open files *)
         f1Stream = OpenWrite[f1File];
         f2Stream = OpenWrite[f2File];
         f3Stream = OpenWrite[f3File];
         f4Stream = OpenWrite[f4File];
         f5Stream = OpenWrite[f5File];
         f6Stream = OpenWrite[f6File];
         (* f7Stream = OpenWrite[f7File];
            f8Stream = OpenWrite[f8File]; *)
         
         (* headings *)
         WriteString[f1Stream, "# Motions from sinusoidal Unit f1 (along z) force\n"];
         WriteString[f1Stream, "#f\tz[4]\tz[5]\tz[6]\tz[7]\tz[8]\n"];
         
         WriteString[f2Stream, "# Motions from sinusoidal Unit f2 (along ty) force\n"];
         WriteString[f2Stream, "#f\tty[4]\tx[4]\tty[5]\tx[5]\tty[6]\tx[6]\tty[7]\tx[7]\tty[8]\tx[8]\n"];
         
         WriteString[f3Stream, "# Motions from sinusoidal Unit f3 (along x) force\n"];
         WriteString[f3Stream, "#f\tx[4]\tty[4]\tx[5]\tty[5]\tx[6]\tty[6]\tx[7]\tty[7]\tx[8]\tty[8]\n"];
         
         WriteString[f4Stream, "# Motions from sinusoidal Unit f4 (along tz) force\n"];
         WriteString[f4Stream, "#f\ttz[4]\ttz[5]\ttz[6]\ttz[7]\ttz[8]\n"];
         
         WriteString[f5Stream, "# Motions from sinusoidal Unit f5 (along x) force\n"];
         WriteString[f5Stream, "#f\tx[4]\tty[4]\tx[5]\tty[5]\tx[6]\tty[6]\tx[7]\tty[7]\tx[8]\tty[8]\n"];
         
         WriteString[f6Stream, "# Motions from sinusoidal Unit f6 (along tz) force\n"];
         WriteString[f6Stream, "#f\ttz[4]\ttz[5]\ttz[6]\ttz[7]\ttz[8]\n"];
         
         
         (* Output Indices *)
         Do[xIndex[i] = Position[Q, x[i]][[1]][[1]];
            tyIndex[i] = Position[Q, ty[i]][[1]][[1]];
            zIndex[i] = If[i > 5, Position[Q, z[i]][[1]][[1]], Position[Q, zVessel[i]][[1]][[1]]];
            tzIndex[i] = Position[Q, tz[i]][[1]][[1]],
            {i,4,8}];
                     
         (* Cycle over frequencies *)
         Do[nu = N[Exp[logNu]];\
         	TiltTransNu = N[TiltTrans[nu]];
         	PosSeismNu  = N[PosSeism[nu]];
         	TiltSeismNu = N[TiltSeism[nu]];
            (* generate a part of the full transfer matrix *)
            ris = Abs /@ N[TransferMatrixForces[nu]];
            WriteString[f1Stream, CForm[nu],
                          "\t", CForm[ris[[1]][[zIndex[4]]]],
                          "\t", CForm[ris[[1]][[zIndex[5]]]],
                          "\t", CForm[ris[[1]][[zIndex[6]]]],
                          "\t", CForm[ris[[1]][[zIndex[7]]]],
                          "\t", CForm[ris[[1]][[zIndex[8]]]],
                          "\n"];
            WriteString[f2Stream, CForm[nu],
                          "\t", CForm[ris[[2]][[tyIndex[4]]]],
                          "\t", CForm[ris[[2]][[xIndex[4]]]],
                          "\t", CForm[ris[[2]][[tyIndex[5]]]],
                          "\t", CForm[ris[[2]][[xIndex[5]]]],
                          "\t", CForm[ris[[2]][[tyIndex[6]]]],
                          "\t", CForm[ris[[2]][[xIndex[6]]]],
                          "\t", CForm[ris[[2]][[tyIndex[7]]]],
                          "\t", CForm[ris[[2]][[xIndex[7]]]],
                          "\t", CForm[ris[[2]][[tyIndex[8]]]],
                          "\t", CForm[ris[[2]][[xIndex[8]]]],
                          "\n"];
            WriteString[f3Stream, CForm[nu],
                          "\t", CForm[ris[[3]][[xIndex[4]]]],
                          "\t", CForm[ris[[3]][[tyIndex[4]]]],
                          "\t", CForm[ris[[3]][[xIndex[5]]]],
                          "\t", CForm[ris[[3]][[tyIndex[5]]]],
                          "\t", CForm[ris[[3]][[xIndex[6]]]],
                          "\t", CForm[ris[[3]][[tyIndex[6]]]],
                          "\t", CForm[ris[[3]][[xIndex[7]]]],
                          "\t", CForm[ris[[3]][[tyIndex[7]]]],
                          "\t", CForm[ris[[3]][[xIndex[8]]]],
                          "\t", CForm[ris[[3]][[tyIndex[8]]]],
                          "\n"];
            WriteString[f4Stream, CForm[nu],
                          "\t", CForm[ris[[4]][[tzIndex[4]]]],                          
                          "\t", CForm[ris[[4]][[tzIndex[5]]]],
                          "\t", CForm[ris[[4]][[tzIndex[6]]]],
                          "\t", CForm[ris[[4]][[tzIndex[7]]]],
                          "\t", CForm[ris[[4]][[tzIndex[8]]]],
                          "\n"];
            WriteString[f5Stream, CForm[nu],
                          "\t", CForm[ris[[5]][[xIndex[4]]]],
                          "\t", CForm[ris[[5]][[tyIndex[4]]]],
                          "\t", CForm[ris[[5]][[xIndex[5]]]],
                          "\t", CForm[ris[[5]][[tyIndex[5]]]],
                          "\t", CForm[ris[[5]][[xIndex[6]]]],
                          "\t", CForm[ris[[5]][[tyIndex[6]]]],
                          "\t", CForm[ris[[5]][[xIndex[7]]]],
                          "\t", CForm[ris[[5]][[tyIndex[7]]]],
                          "\t", CForm[ris[[5]][[xIndex[8]]]],
                          "\t", CForm[ris[[5]][[tyIndex[8]]]],
                          "\n"];
            WriteString[f6Stream, CForm[nu],
                          "\t", CForm[ris[[6]][[tzIndex[4]]]],                          
                          "\t", CForm[ris[[6]][[tzIndex[5]]]],
                          "\t", CForm[ris[[6]][[tzIndex[6]]]],
                          "\t", CForm[ris[[6]][[tzIndex[7]]]],
                          "\t", CForm[ris[[6]][[tzIndex[8]]]],
                          "\n"],
            {logNu, logNuMin, logNuMax, logNuStep}];
        Close /@ {f1Stream, f2Stream, f3Stream, f4Stream, f5Stream, f6Stream};
         ];
         
SaveTransferForces::usage="SaveTransferForces[nuMin_, nuMax_, nuPoints_] computes
	some transfer functions between the soil motion and parts of the S.A.,
	using LU decomposition and saving the results in files.";

Transfer Functions

Eigenvalues and Eigenvectors

SAEigenVec::usage="Eigenvectors of the SA";
SAEigenVal::usage="Eigenvalues of the SA";

Module[{eigenMat, signs, i,om},
		eigenMat = Inverse[KinMat] . PotMat;
		Print["KinMat inverted."];
		SAEigenVec = Eigenvectors[eigenMat];
		Print["Eigenvalues computed."];
		SAEigenVal = Eigenvalues[eigenMat];
		(* check signs *)
		signs = Sign /@ (Abs /@ SAEigenVal);
		Print["signs = ", signs];
		Print["Eigenvectors found."];
		(* Normalization of vectors *)
		Do[SAEigenVec[[i]] = SAEigenVec[[i]] /
			Sqrt[SAEigenVec[[i]] . KinMat . SAEigenVec[[i]]],
		   {i,1,Length[SAEigenVec]}];
		Print["Eigenvectors normalized."];
		(* General solution *)
		TransferSolution = Sum[1/(SAEigenVal[[i]] - om^2)*
						SAEigenVec[[i]] *
						(SAEigenVec[[i]] . CoupMat . Q0),
					   {i,1,Length[SAEigenVec]}];
		Print["General solution found."];
		(* General solution for Forces *)
		TransferSolutionForces = 
		           Sum[1/(SAEigenVal[[i]] - om^2)*
						SAEigenVec[[i]] *
						(SAEigenVec[[i]] . ExtMat . Forces),
					   {i,1,Length[SAEigenVec]}];
		Print["General solutionForces found."];
		(* Different transfer functions *)
		TransferFunction[i_, j_, nu_] :=
			(TransferSolution[[j]]/.(Q0[[i]] -> 1))/. zeroRule/.(om -> 2 Pi nu);
		Print["TransferFunction defined"];
		TransferFunctionForces[i_, j_, nu_] :=
			(TransferSolutionForces[[j]]/.(Forces[[i]] -> 1))/. zeroRule/.(om -> 2 Pi nu);
		Print["TransferFunctionForces defined"];
		];

TransferSolution::usage="TransferSolution is a vector containing the general motion
of the elements of the SA due the motion of the soil";

TransferSolutionForces::usage="TransferSolutionForces is a vector containing the general
motion of the elements of the SA under the internal forces";

TransferFunction::usage="TransferFunction[input_, output_, nu_]
returns the value of the complex transfer functions between
points input and output of the chain.
The input index needs to be in the range of the length of
Q0; the output index in the range of the length of Q.";

TransferFunctionForces::usage=
"TransferFunctionForces[input_, output_, nu_]
returns the value of the complex transfer functions from
one of the forces in the force vector Forces to
one of the elements of the position vector Q";