clc;
clear;

%Wyliczone wcześniej wartośći w1 i w2
w1=6.719; %wartość w1 w rad/s
w2=14.884; %wartość w2 w rad/s

%Dane układu
m1=2; %masa 1 w kg
m2=3; %masa 2 w kg
c1=200; %współczynnik sprężystości 1 w N/m
c2=300; %współczynnik sprężystości 2 w N/m
c3=400; %współczynnik sprężystości 3 w N/m
F=10; %siła działająca na układ w N

%Dane cewki
%I=0.87173; %natężenie prądu w A
N=100; %ilość uzwojeń
r_r=3*10^(-3); %promień rdzenia w m
ob_r=2*pi*r_r; %obwód rdzenia m
sr_dr=1.2*10^(-3);  %średnica drutu m
l_uzw=ob_r; %długość zwoju m
L1=N*l_uzw; %długość przewodnika m
L2=N*sr_dr; %długość cewki m
u=4000; %względna przenikalność magnetyczna
u0=4*pi*10^(-7); %przenikalność magnetyczna próżni (stała) w H/m
%H=N*I/L2; %natężenie pola magnetycznego w A/m
%B=u*u0*H; %indukcja pola magnetycznego w N/A*m
%G=B*I*L1; %siła z jaką działa cewka w N

%Obliczenia szukanego natężenia prądu
%Wprowadzenie minimalnych aplitud
A_1=0.001;
A_2=0.002;

Y1=[-m1*w1.^2+c1 -c1; -c1 -m2*w1.^2+c1+c2]; %macierz dla w1
Y2=[-m1*w2^2+c1 -c1; -c1 -m2*w2.^2+c1+c2]; %macierz dla w2
Z=[F; 0];
B=[A_1;A_2];
G1=Y1*B;
G2=Y2*B;
G1_szuk=abs(G1+Z)
G2_szuk=abs(G2+Z)

I1_szuk=sqrt((G1_szuk*L2)/(N*L1*u*u0))
I2_szuk=sqrt((G2_szuk*L2)/(N*L1*u*u0))

%Obliczenie apmplitud dla układu po doprowadzeniu wyliczonego natęzenia
%prądu do cewek
w=0:0.001:20;
x=0;
A1_k=zeros(1,20000);
A2_k=zeros(1,20000);
for w=0:001:20
    if w==w1 && w==w2
    W=abs(c1*c2-c1*m1*w.^2-c1*m2*w.^2-c2*m1*w.^2+m1*m2*w.^4);
    A1=abs((F*(-m2*w.^2+c1+c2))./W)
    A2=(c1*F)./W
    elseif w=w1
    G=(I1_szuk(1,1).^2*L1*N*u*u0)/L2;
    G2=(I1_szuk(1,2).^2*L1*N*u*u0)/L2;
    W=abs(c1*c2-c1*m1*w.^2-c1*m2*w.^2-c2*m1*w.^2+m1*m2*w.^4); %wyznacznik macierzy
    A1=abs(((F-G)*(-m2*w.^2+c1+c2))/W);
    A2=(c1*(F-G2))./W;
    elseif w=w2
    G=(I2_szuk(1,1).^2*L1*N*u*u0)/L2;
    G2=(I2_szuk(1,2).^2*L1*N*u*u0)/L2;
    W=abs(c1*c2-c1*m1*w.^2-c1*m2*w.^2-c2*m1*w.^2+m1*m2*w.^4); %wyznacznik macierzy
    A1=abs(((F-G)*(-m2*w.^2+c1+c2))/W);
    A2=(c1*(F-G2))./W;
    end
end

subplot(2,1,1)
plot(w,A1_k)
subplot(2,1,2)
plot(w,A2_k)