proj4

function MN_PIASECKI

nr=input('Wymuszenie (1, 2, 3, 4, 5, 6, 7) ');
while (nr>7 | mod(nr,1)~=0 | nr<1)
    disp('Zly numer wymuszenia');
    nr=input('Wymuszenie (1, 2, 3, 4, 5, 6, 7) ');
end

t0=0;%s%
tk=0.06;%s%
R1=1000 %Ohm%
C1=0.000000032;%F%
R2=1000 %Ohm%
C2=0.000000032;%F%
RL=1000 %Ohm%

u10=0;
u20=0;
il0=0;

n=5000;
h=(tk-t0)/n;

u1=[];
u2=[];
il=[];

u1(1)=u10;
u2(1)=u20;
il(1)=il0;

t=t0:h:tk;

for i=1:n
    u1(i+1)=u1(i)+h*fu1(t(i),u1(i),u2(i),R1,C1,R2,C2,RL,nr);
    u2(i+1)=u2(i)+h*fu2(t(i),u1(i),u2(i),R1,C1,R2,C2,RL,nr);
    il(i+1)=il(i)+h*fi(t(i),u2(i),il(i),R1,C1,R2,C2,RL,nr);
end

s=0;
for i=0:n
    if i==0
        s=s+fc(il(i+1),u2(i+1),u1(i+1),RL);
    elseif i==n
        s=s+fc(il(i+1),u2(i+1),u1(i+1),RL);
    elseif mod(i,2)==0
        s=s+4*fc(il(i+1),u2(i+1),u1(i+1),RL);
    elseif mod(i,2)==1
        s=s+2*fc(il(i+1),u2(i+1),u1(i+1),RL);
end

end

P=1000*(h/3)*s;
wynik=['P = ',num2str(P)];
disp(wynik);

end


function [wymuszenie]=E(t,nr)
    if nr==1
        wymuszenie=2;
    elseif nr==2
        wymuszenie=sin(2*pi*50*t);
    elseif nr==3
        wymuszenie=sin(2*pi*600*t);
    elseif nr==4
        wymuszenie=sin(2*pi*1750*t);
    elseif nr==5
        wymuszenie=sin(2*pi*12000*t);
    elseif nr==6
        wymuszenie=sin(2*pi*21000*t);
    elseif nr==7
        if mod(10000*t,5)<(0.5/2)
            wymuszenie=1;
        elseif mod(10000*t,5)>=(0.5/2)
            wymuszenie=0;
        end
    end
end

function [du1]=fu1(t,u1,u2,R1,C1,R2,C2,RL,nr)
    du1=(1/(C1*R1))*(E(t,nr)-(u1-u2));
end

function [du2]=fu2(t,u1,u2,R1,C1,R2,C2,RL,nr)
    du2=(1/C2)*(((u1-u2)/R2)-(u2/RL));
end

function [dil]=fi(t,u1,u2,il,R1,C1,R2,C2,RL,nr)
    dil=u2/RL;
end

function [int]=fc(il,u1,u2,RL,R2)
    int=(u2*u2)/RL;
end