%'Pragram for Newmark explicit integration method for MDOF systems%
%Special application: Discontinous loading with linear variations%
%follow the instructions of the program carfully!%
% when you need to enter the values of a matrix enter the values row by row%
%Written by Morteza Razi%

close all;
clc;
clear all;
nf= input('whats the degree of freedom?  ');
for i=1:nf
    disp('enter the values for degree:');
    disp(i);
    M(i)=input('input the mass  value:  ');
    d(i,1)=input('enter u0: ');
    v(i,1)=input('enter v0: ');
    for j=1:nf 
    disp([i,j]);
        k(i,j)=input('input stiffness next value---  ')
    m(i,j)=0;
    end;
    m(i,i)=M(i);
end;
gama=input('gama?   ');
betta=input('beeta=?  ');
answer1=input('Do you have the value of damping matrix components? y/n ==> ','s');
if answer1=='y' 
disp('enter the valuse of matrix c respectively ==>  ');
    for x=1:nf
        for y=1:nf
            c(i,j)=input('next value of c ==> ');
        end;
    end;
else
kesai=input('input the kesai value(damping ratio)==>  ');
omega=input('enter the circular frequency==>  ');    
    c=m*2*kesai*omega;
        end;
n=input('enter the number of steps  ');
duration=input('enter duration time  ');
delta_t=0.000001*fix(1000000*duration/n);
for i=1:n+1
    t(i)=(i-1)*delta_t;
    end;
span_num=input('input the number of time spans (discontinuty points minus 1) ? => -');
   for i=1:span_num
   disp('for span number');
   disp(i);
   span(i)=input('enter the timespan duration');
   end;
    timeendspan(1)=span(1);
    timestartspan(1)=0;
    for i=2:span_num
    timeendspan(i)=timeendspan(i-1)+span(i)
    timestartspan(i)=timeendspan(i-1);
    end;
        for i=1:span_num
        startstep(i)=fix(timestartspan(i)/delta_t)+1;
        endstep(i)=fix(timeendspan(i)/delta_t);
        spanstep(i)=fix((timeendspan(i)-timestartspan(i))/delta_t)+1;
        end;
    for i=1:span_num
         for j=1:nf
disp ('enter the values for span and degree of');
disp([i,j]);
         jump(i,j)=input('enter the first value of A');
         delta_A_lin=input('value of a for formula: delta_A=(a*t+b) for this step ');
         a_var=input('enter the a for formula: delta_A=(a*t+b) for this step');        
         p=startstep(i);
         delta_A(j,p)=jump(i,j); 
             for q=(startstep(i)+1):endstep(i)
                t_in_span=(q-startstep(i))*delta_t
                delta_A(j,q)=a_var*t_in_span+delta_A_lin;
             end;
        end;
    end;
    A=jump;
    a=m^(-1)*(-k*d-c*v+A(:,1));
    for dd=1:nf
    displac(dd,1)=d(dd);
    velo(dd,1)=v(dd);
    accel(dd,1)=a(dd);
    end;
for j=1:n
q_hat=1/betta/delta_t*v+a/2/betta;
r_hat=gama/betta*v+(gama/2/betta-1)*delta_t*a;
delta_A_hat=delta_A(:,j)+m*q_hat+c*r_hat;
s_hat=k+1/(betta*delta_t^2)*m+gama/betta/delta_t*c;
delta_d=inv(s_hat)*delta_A_hat;
delta_a=1/(betta*delta_t^2)*delta_d-q_hat;
delta_v=gama/(betta*delta_t)*(delta_d)-(r_hat);
disp(i);
d=d+delta_d;
v=v+delta_v;
a=a+delta_a;
for kk=1:nf
    displac(kk,j+1)=d(kk);
    velo(kk,j+1)=v(kk);
    accel(kk,j+1)=a(kk);
end;
end;
t(n+1)=t(n)+delta_t;
for i=1:nf
plot(t,displac(i,:));
end;
        for j=1:(n+1)
        javab(1,j)=t(j);
        ind=1;
        for jj=1:nf;
            ind=ind+1;
            javab(ind,j)=displac((ind-1),j);
            end;
       end;
fid = fopen('Displacements.txt', 'wt');
fprintf(fid, '%8.4f%8.4f%8.4f\n', javab);
fclose(fid);
disp('end');