kill(all);

dati recuperati da modale FEM, 
modi propri normalizzati a massa modale unitaria.

m:5;

fmin:0$
fmax:1000$

zeta:0.01$

frequenze proprie

freq:[
83.37,
524.7,
828.9,
1488,
1608
];

spostamenti dell'unico g.d.l. caricato

disp:[
863.197,
-865.160,
0,
-870.81,
0
];

azione esterna sull'unico g.d.l. caricato

fr:1.2$
fi:0$

una quantità da ottenere in output

s11:[
-11673.2,
-65554.5,
0,
166612,
0
];

ciclo sui modi, calcolando per ogni l-esimo modo 
il fattore scalare  complesso
(xir[l] + %i xii[l])
che va a modularlo. 
Utilizzo nella sua definizione la risposta dell'oscillatore 1d 
equivalente a quel modo

xir_list:[]$
xii_list:[]$
for l : 1 thru m do (
    qr:fr*disp[l],
    qi:fi*disp[l],
    omega: 2*%pi*freq[l],
    r:f/freq[l],
    xir: (2*qi*r*zeta-qr*r^2+qr)/(4*omega^2*r^2*zeta^2+omega^2*r^4-2*omega^2*r^2+omega^2),
    xii:-(2*qr*r*zeta+qi*r^2-qi)/(4*omega^2*r^2*zeta^2+omega^2*r^4-2*omega^2*r^2+omega^2),
    xir_list:append(xir_list,[xir]),
    xii_list:append(xii_list,[xii])
);

in alternativa posso utilizzare le formule in fondo al 
paragrafo sulla risposta armonica per sovrapposizione modale
qui estraggo direttamente modulo e fase

xir_list2:[]$
xii_list2:[]$
for l : 1 thru m do (
    qr:fr*disp[l],
    qi:fi*disp[l],
    omega: 2*%pi*freq[l],
    r:f/freq[l],
    a:1-r^2,
    b:2*zeta*r,
    xir:(a*qr+b*qi)/omega^2/(a^2+b^2),
    xii:(a*qi-b*qr)/omega^2/(a^2+b^2),
    xir_list2:append(xir_list2,[xir]),
    xii_list2:append(xii_list2,[xii])
);

scostamento numerico tra le due soluzioni

wxplot2d(
append(
    (xir_list-xir_list2),
    (xii_list-xii_list2)
),
[f,fmin,fmax]
);

nullo...

mostro i fattori modulanti associati ai modi propri

for l : 1 thru m do(
    disp(
        [
            xi[l],
            ratsimp(
                xir_list[l]+%i*xii_list[l]
            )
        ]
    )
),ratprint=false;

procedo all'accumulo dei contributi tensionali associati ai vari modi

s11r:0$
s11i:0$
for l : 1 thru m do(
    disp(l),
    wxplot2d(
        s11[l]*[xir_list[l],xii_list[l]],
        [f,fmin,fmax],
        [legend,"Re","Im"]
    ),
    [s11r,s11i]:[s11r,s11i]+s11[l]*[xir_list[l],xii_list[l]]
);

valore statico, da FEM statico 39.2043

[s11r,s11i],f=0,numer;

plotto i risultati

wxplot2d(
[s11r,s11i,sqrt(s11r^2+s11i^2)],
[f,fmin,fmax],
[legend,"Re","Im","Mag"]
);

zoom nell'intorno della prima risonanza

wxplot2d(
[s11r,s11i,sqrt(s11r^2+s11i^2)],
[f,freq[1]-5,freq[1]+5],
[legend,"Re","Im","Mag"]
);

andamento al di fuori delle risonanze 
(zoom su valori limitati di amplificazione)

wxplot2d(
[s11r,s11i,sqrt(s11r^2+s11i^2)],
[f,fmin,fmax],[y,-100,100],
[legend,"Re","Im","Mag"]
);

magnitudo su scala logaritmica

wxplot2d(
[sqrt(s11r^2+s11i^2)],
[f,fmin,fmax],
[legend,"Mag"],
[logy]
);