from Structure analysis by Yasser
computation of deflections in various structures such as Frame, Truss, Plate,...

Frame.m

clear; clc; close all; 
disp ('                                                                                  ')
disp ('                                                                                  ')
disp ('                                                                                  ')
disp ('=================================================================================')
disp ('=================================================================================')
disp ('                                                                                  ')
disp ('                                                                                  ')
disp ('                                                                                  ')


disp (' ***:*** planner truss=t2,space truss=t3,grid,planner frame=f2,space frame=f3 *** ')
disp ('                 By: Yasser Bigdeli (2012) Email: yasse.bigdeli@gmail.com                         ')
disp ('                                                                                  ')
disp ('                                                                                  ')
disp ('                                                                                  ')
disp ('=================================================================================')
disp ('=================================================================================')
syms t2 t3 grid f2 f3
S=input('*** Name structure ***:*** planner truss=t2,space truss=t3,grid,planner frame=f2,space frame=f3 *** ');
format short 
N=input('*** How many nodes have in this structure( N = ?).  ? *** ');
for I=1:N
    I=I
    X(I)=input('X(I) =');
    Y(I)=input('Y(I) =');
    if S==t3 || S==f3 
        Z(I)=input('Z(I) =');
    else
        Z(I)=0;
    end
end
t=0;l=0;D=0;A=0;E=0;clear Iz(t);clear Iy(t);clear G(t);clear Jo(t);clear Ang(t);clear y;clear n
syms y n no
c=input('*** Dose materials specifications like together ? " yes=y  or no=n " ***');
if S==f3
    e=input('*** Have structural sections Angle with principle axis in local coordinate ? " yes=y  or no=n " ***');
    dd=input('*** Are all element Angle same ? " yes=y  or no=no " ***');
end
for I=1:N
    for J=I+1:N
        I,J
        Q=input('*** if is not element between in the nodse I and J that plase input zero else input every number ? ***');
        if Q~=0
            t=t+1
            B=t;
            u(t)=I;
            v(t)=J;
            if S~=f3
                Ang(t)=0;
            end
            Iz(t)=0;Iy(t)=0;G(t)=0;Jo(t)=0;Ang(t)=0;
            if t==1
                E(t)=input('E(t) =');
                Ang(t)=0;
                if S==f2
                    A(t)=0;
                end
                if S~=grid & S~=f2
                    A(t)=input('A(t) =');
                else
                    A(t)=0;
                end
                if S==f3 || S==grid
                    Iy(t)=input('Iy(t) =');
                    G(t)=input('G(t) =');
                    Jo(t)=input('Jo(t) =');
                    if S==f3
                        Iz(t)=input('Iz(t) =');
                    end
                end
                if S==f2
                    Iz(t)=input('Iz(t) =');
                    Iy(t)=0;
                    G(t)=0;
                    Jo(t)=0;
                end
                if S==f3 && e==y
                    Ang(t)=input('Angle(t) =');
                else 
                    Ang(t)=0;
                end
            elseif c==y 
                E(t)=E(t-1);
                A(t)=A(t-1);
                Iy(t)=Iy(t-1);
                Iz(t)=Iz(t-1);
                G(t)=G(t-1);
                Jo(t)=Jo(t-1);
                if S==f3 && dd==no  
                    Ang(t)=input('Angle(t) =');
                elseif S==f3 && dd==y
                        Ang(t)=Ang(t-1);
                end
            else
                E(t)=input('E(t) =');
                if S==f2
                    A(t)=0;
                end
                if S~=grid & S~=f2
                    A(t)=input('A(t) =');
                end
                if S==f3 || S==grid
                    Iy(t)=input('Iy(t) =');
                    G(t)=input('G(t) =');
                    Jo(t)=input('Jo(t) =');
                    if S==f3
                        Iz(t)=input('Iz(t) =');
                    end
                end
                if S==f2
                    Iz(t)=input('Iz(t) =');
                    Iy(t)=0;
                    G(t)=0;
                    Jo(t)=0;
                end
                if S==f3 && e==y && dd==no
                    Ang(t)=input('Angle(t) =');
                else
                    Ang(t)=0;
                end
            end
        L(t)=(((X(I)-X(J))^2+(Y(I)-Y(J))^2+(Z(I)-Z(J))^2)^(1/2));
        l(t)=(X(I)-X(J))/L(t);
        m(t)=(Y(I)-Y(J))/L(t);
        n(t)=(Z(I)-Z(J))/L(t);
        D(t)=((l(t))^2+(m(t))^2)^(1/2);
        end
    end
end
w=zeros(12,12);
k=zeros(N,N);r=cell(B,1);R=cell(B,1);
if S==t2
    i=2;
elseif S==t3 || S==grid || S==f2
    i=3;
elseif S==f3
    i=6;
end
j=i-1;    
K=zeros(j*N,j*N);
P=zeros(i*N,1);
p=zeros(i,1);
Ppij=zeros(i,1);
Ppji=zeros(i,1);
U=zeros(i*N,1);
U2=zeros(i*N,1);
if S==grid || S==f2 || S==f3
    q1=zeros(i*N,i*N);
    q2=zeros(i*N,i*N);
    q3=zeros(i*N,i*N);
    q4=zeros(i*N,i*N);
    q5=zeros(i*N,i*N);
    q6=zeros(i*N,i*N);
    z1=zeros(i*N,i*N);
    z2=zeros(i*N,i*N);
    r=zeros(i*N,i*N);
    R=zeros(i*N,i*N);
    r1=zeros(i*N,i*N);r2=zeros(i*N,i*N);r3=zeros(i*N,i*N);r4=zeros(i*N,i*N);
    q11=zeros(i*N,i*N);
    q44=zeros(i*N,i*N);
end
aa=zeros(2*i*N,1);
s=zeros(i,1);
zzz=zeros(i,1);
clear WI(t,I);
clear WI(t,J);
for t=1:B
    I=u(t);
    J=v(t);
    if S==grid || S==f2 || S==f3
        [I J];
        k(6*I-5:6*I,6*I-5:6*I)=[(E(t)*A(t))/L(t) 0 0 0 0 0 
            0 (12*(E(t)*Iz(t))/((L(t))^3)) 0 0 0 (-6*(E(t)*Iz(t))/((L(t))^2))
            0 0 (12*(E(t)*Iy(t))/((L(t))^3)) 0 (6*(E(t)*Iy(t))/((L(t))^2)) 0
            0 0 0 ((G(t)*Jo(t))/L(t)) 0 0
            0 0 (6*(E(t)*Iy(t))/((L(t))^2)) 0 (4*(E(t)*Iy(t))/L(t)) 0 
            0 (-6*(E(t)*Iz(t))/((L(t))^2)) 0 0 0 (4*(E(t)*Iz(t))/L(t))];
        k(6*I-5:6*I,6*J-5:6*J)=[(E(t)*A(t))/L(t) 0 0 0 0 0
            0 (12*(E(t)*Iz(t))/((L(t))^3)) 0 0 0 (-6*(E(t)*Iz(t))/((L(t))^2))
            0 0 (-12*(E(t)*Iy(t))/((L(t))^3)) 0 (-6*(E(t)*Iy(t))/((L(t))^2)) 0
            0 0 0 ((G(t)*Jo(t))/L(t)) 0 0
            0 0 (-6*(E(t)*Iy(t))/((L(t))^2)) 0 (-2*(E(t)*Iy(t))/L(t)) 0
            0 (-6*(E(t)*Iz(t))/((L(t))^2)) 0 0 0 (2*(E(t)*Iz(t))/L(t))];
        k(6*J-5:6*J,6*I-5:6*I)=k(6*I-5:6*I,6*J-5:6*J);
        k(6*J-5:6*J,6*J-5:6*J)=k(6*I-5:6*I,6*I-5:6*I);
        if S==f3
            WI(t,I)={k(6*I-5:6*I,6*I-5:6*I)};
            WI(t,J)={k(6*I-5:6*I,6*J-5:6*J)};
            WJ(t,I)={k(6*J-5:6*J,6*I-5:6*I)};
            WJ(t,J)={k(6*J-5:6*J,6*J-5:6*J)};
        end
        WII=k(6*I-5:6*I,6*I-5:6*I);
        WIJ=k(6*I-5:6*I,6*J-5:6*J);
        WJI=WIJ;
        WJJ=WII;
        i=6;j=5;
        q1(i*I-j:6*I,i*J-j:i*J)=[1 0 0 0 0 0
            0 cos(Ang(t)) sin(Ang(t)) 0 0 0
            0 -sin(Ang(t)) cos(Ang(t)) 0 0 0
            0 0 0 1 0 0 
            0 0 0 0 cos(Ang(t)) sin(Ang(t)) 
            0 0 0 0 -sin(Ang(t)) cos(Ang(t))];
        q11(6*I-5:6*I,6*J-5:6*J)=q1(6*I-5:6*I,6*J-5:6*J);
        q1=q1(6*I-5:6*I,6*J-5:6*J);
        if D(t)~=0
            q2(i*I-j:i*I,i*J-j:i*J)=[l(t) m(t) n(t) 0 0 0
                -m(t)/D(t) l(t)/D(t) 0 0 0 0
                -l(t)*n(t)/D(t) -m(t)*n(t)/D(t) D(t) 0 0 0
                0 0 0 l(t) m(t) n(t)
                0 0 0 -m(t)/D(t) l(t)/D(t) 0
                0 0 0 -l(t)*n(t)/D(t) -m(t)*n(t)/D(t) D(t)];
            q2=q2(6*I-5:6*I,6*J-5:6*J);
        elseif D(t)==0
            q3(i*I-j:i*I,i*J-j:i*J)=[0 0 1 0 0 0
                0 1 0 0 0 0
                -1 0 0 0 0 0
                0 0 0 0 0 1
                0 0 0 0 1 0
                0 0 0 -1 0 0];
            q3=q3(6*I-5:6*I,6*J-5:6*J);
        end
        q4(i*I-j:i*I,i*J-j:i*J)=[1 0 0 0 0 0
            0 -cos(Ang(t)) -sin(Ang(t)) 0 0 0
            0 sin(Ang(t)) -cos(Ang(t)) 0 0 0
            0 0 0 1 0 0 
            0 0 0 0 -cos(Ang(t)) -sin(Ang(t)) 
            0 0 0 0 sin(Ang(t)) -cos(Ang(t))];
        q44(6*I-5:6*I,6*J-5:6*J)=q4(6*I-5:6*I,6*J-5:6*J);
        q4=q4(6*I-5:6*I,6*J-5:6*J);
        if S==grid
            i=3;j=2;
            [I J];
            k(i*I-j:i*I,i*I-j:i*I)=WII(3:5,3:5);
            k(i*I-j:i*I,i*J-j:i*J)=WIJ(3:5,3:5);
            k(i*J-j:i*J,i*I-j:i*I)=WJI(3:5,3:5);
            k(i*J-j:i*J,i*J-j:i*J)=WJJ(3:5,3:5);
            r(i*I-j:i*I,i*J-j:i*J)=q2(3:5,3:5);
            R(i*I-j:i*I,i*J-j:i*J)=[1 0 0 
                                    0 -l(t) -m(t)
                                    0 m(t) -l(t)];
            r2(i*I-j:i*I,i*J-j:i*J)=r(i*I-j:i*I,i*J-j:i*J);
            r1(i*I-j:i*I,i*J-j:i*J)=r(i*I-j:i*I,i*J-j:i*J)';
            r3(i*I-j:i*I,i*J-j:i*J)=R(i*I-j:i*I,i*J-j:i*J);
            r4(i*I-j:i*I,i*J-j:i*J)=R(i*I-j:i*I,i*J-j:i*J)';
        end
        if S==f2
            i=3;j=2;
            k(i*I-j:i*I,i*I-j:i*I)=[WII(1,1) WII(1,2) WII(1,6)
                                    WII(2,1) WII(2,2) WII(2,6)
                                    WII(6,1) WII(6,2) WII(6,6)];
            k(i*I-j:i*I,i*J-j:i*J)=[WIJ(1,1) WIJ(1,2) WIJ(1,6)
                                    WIJ(2,1) WIJ(2,2) WIJ(2,6)
                                    WIJ(6,1) WIJ(6,2) WIJ(6,6)];
            k(i*J-j:i*J,i*I-j:i*I)=[WJI(1,1) WJI(1,2) WJI(1,6)
                                    WJI(2,1) WJI(2,2) WJI(2,6)
                                    WJI(6,1) WJI(6,2) WJI(6,6)];
            k(i*J-j:i*J,i*J-j:i*J)=[WJJ(1,1) WJJ(1,2) WJJ(1,6)
                                    WJJ(2,1) WJJ(2,2) WJJ(2,6)
                                    WJJ(6,1) WJJ(6,2) WJJ(6,6)];
            r(i*I-j:i*I,i*J-j:i*J)=[q2(1,1) q2(1,2) q2(1,6)
                                    q2(2,1) q2(2,2) q2(2,6)
                                    q2(6,1) q2(6,2) q2(6,6)];
            R(i*I-j:i*I,i*J-j:i*J)=[-l(t) -m(t) 0 
                                    m(t) -l(t) 0
                                    0 0 1];
            r2(i*I-j:i*I,i*J-j:i*J)=r(i*I-j:i*I,i*J-j:i*J);
            r1(i*I-j:i*I,i*J-j:i*J)=r(i*I-j:i*I,i*J-j:i*J)';
            r3(i*I-j:i*I,i*J-j:i*J)=R(i*I-j:i*I,i*J-j:i*J);
            r4(i*I-j:i*I,i*J-j:i*J)=R(i*I-j:i*I,i*J-j:i*J)';
        end
        if S==f3       
            if D(t)~=0  
                r1(i*I-j:i*I,i*J-j:i*J)=q2'*q1';
                r2(i*I-j:i*I,i*J-j:i*J)=q1*q2;
                q5(i*I-j:i*I,i*J-j:i*J)=[-l(t) -m(t) -n(t) 0 0 0
                    m(t)/D(t) -l(t)/D(t) 0 0 0 0
                    l(t)*n(t)/D(t) m(t)*n(t)/D(t) D(t) 0 0 0
                    0 0 0 -l(t) -m(t) -n(t)
                    0 0 0 m(t)/D(t) -l(t)/D(t) 0
                    0 0 0 l(t)*n(t)/D(t) m(t)*n(t)/D(t) D(t)];
                q5=q5(6*I-5:6*I,6*J-5:6*J);
                r3(i*I-j:i*I,i*J-j:i*J)=q4*q5;
                r4(i*I-j:i*I,i*J-j:i*J)=q5'*q4';
            else
                r1(i*I-j:i*I,i*J-j:i*J)=q3'*q1';
                r2(i*I-j:i*I,i*J-j:i*J)=q1*q3;
                q6(i*I-j:i*I,i*J-j:i*J)=[0 0 -1 0 0 0
                    0 -1 0 0 0 0
                    -1 0 0 0 0 0
                    0 0 0 0 0 -1
                    0 0 0 0 -1 0
                    0 0 0 -1 0 0];
                q6=q6(6*I-5:6*I,6*J-5:6*J);
                r3(i*I-j:i*I,i*J-j:i*J)=q4*q6;
                r4(i*I-j:i*I,i*J-j:i*J)=q6'*q4';
            end
        end
        [I J];
        WI(t,I)={k(i*I-j:i*I,i*I-j:i*I)};
        WI(t,J)={k(i*I-j:i*I,i*J-j:i*J)};
        WJ(t,I)={k(i*J-j:i*J,i*I-j:i*I)};
        WJ(t,J)={k(i*J-j:i*J,i*J-j:i*J)};
        K(i*I-j:i*I,i*I-j:i*I)=r1(i*I-j:i*I,i*J-j:i*J)*k(i*I-j:i*I,i*I-j:i*I)*r2(i*I-j:i*I,i*J-j:i*J)+K(i*I-j:i*I,i*I-j:i*I);
        K(i*I-j:i*I,i*J-j:i*J)=r1(i*I-j:i*I,i*J-j:i*J)*k(i*I-j:i*I,i*J-j:i*J)*r3(i*I-j:i*I,i*J-j:i*J);
        K(i*J-j:i*J,i*I-j:i*I)=r4(i*I-j:i*I,i*J-j:i*J)*k(i*J-j:i*J,i*I-j:i*I)*r2(i*I-j:i*I,i*J-j:i*J);
        K(i*J-j:i*J,i*J-j:i*J)=r4(i*I-j:i*I,i*J-j:i*J)*k(i*J-j:i*J,i*J-j:i*J)*r3(i*I-j:i*I,i*J-j:i*J)+K(i*J-j:i*J,i*J-j:i*J);
    elseif S==t2 || S==t3
        [t I J];
        W=[(E(t)*A(t))/L(t)];
        WI(t,I)={W};
        WI(t,J)={W};
        WJ(t,I)={W};
        WJ(t,J)={W};
        k(I,I)=W;
        k(I,J)=W;
        k(J,I)=W';
        k(J,J)=W;
        if S==t2
            r=[l(t) m(t)];
            R=-r;
            K(2*I-1:2*I,2*I-1:2*I)=r'*k(I,I)*r+K(2*I-1:2*I,2*I-1:2*I);
            K(2*I-1:2*I,2*J-1:2*J)=r'*k(I,J)*R;
            K(2*J-1:2*J,2*I-1:2*I)=R'*k(J,I)*r;
            K(2*J-1:2*J,2*J-1:2*J)=R'*k(J,J)*R+K(2*J-1:2*J,2*J-1:2*J);
        else
            r=[l(t) m(t) n(t)];
            R=-r;
            K(3*I-2:3*I,3*I-2:3*I)=r'*k(I,I)*r+K(3*I-2:3*I,3*I-2:3*I);
            K(3*I-2:3*I,3*J-2:3*J)=r'*k(I,J)*R;
            K(3*J-2:3*J,3*I-2:3*I)=R'*k(J,I)*r;
            K(3*J-2:3*J,3*J-2:3*J)=R'*k(J,J)*R+K(3*J-2:3*J,3*J-2:3*J);
        end    
    end
end
A=K;
a=input('*** How many nodes have externed force and moment   ? *** ');
x=0;
for x=1:a
    I=input('*** What node is that (I) ? *** ');
    if S==t2
        P(2*I-1,1)=input('PX(I) =');
        P(2*I,1)=input('PY(I) =');
    elseif S==t3
        P(3*I-2,1)=input('PX(I) =');
        P(3*I-1,1)=input('PY(I) =');
        P(3*I,1)=input('PZ(I) =');
    elseif S==grid
        P(3*I-2,1)=input('PZ(I) =');
        P(3*I-1,1)=input('MX(I) =');
        P(3*I,1)=input('MY(I) =');
    elseif S==f2
        P(3*I-2,1)=input('PX(I) =');
        P(3*I-1,1)=input('PY(I) =');
        P(3*I,1)=input('MZ(I) =');
    elseif S==f3
        P(6*I-5,1)=input('PX(I) =');
        P(6*I-4,1)=input('PY(I) =');
        P(6*I-3,1)=input('PZ(I) =');
        P(6*I-2,1)=input('MX(I) =');
        P(6*I-1,1)=input('MY(I) =');
        P(6*I,1)=input('MZ(I) =');
    end
end
P2=P;
h=0;
K2=A;
C=zeros(j*N,1);
aa=zeros(12*N,1);
Ismo=input('*** How many nodes have the boundry condition   ? *** ');
if Ismo>0
    for x=1:Ismo
        F=input('*** Which node have boundry condition (I) ? *** ');
        C(x,1)=F;
        smij=input('*** How many  boundry condition is in this node  ? *** ');
        f(F)=smij;
        if smij>0
            j=0;i=0;
            if S==f3 && smij==6
                for g=1:smij
                    s(g,1)=g;
                    aa((6*F-6)+s(g,1),1)=s(g);
                    i=6*F-6+g-h;
                    K(i,:)=[];
                    K(:,i)=[];
                    P(i,:)=[];
                    U(i,:)=[];
                    h=h+1;
                end
            else
                if S==t3 || S==grid || S==f2 
                    if smij==3
                        for g=1:smij
                            s(g,1)=g;
                            aa((3*F-3)+s(g,1),1)=s(g,1);
                            i=(3*F)+g-3-h;
                            K(i,:)=[];
                            K(:,i)=[];
                            P(i,:)=[];
                            U(i,:)=[];
                            h=h+1;
                        end
                    else
                        for g=1:smij
                            s(g,1)=input('s(g,1) =');
                            aa((3*F-3)+s(g,1),1)=s(g,1);
                            i=(3*F)+s(g,1)-3-h;
                            K(i,:)=[];
                            K(:,i)=[];
                            P(i,:)=[];
                            U(i,:)=[];
                            h=h+1;
                        end
                    end
                else
                    if S==t2 && smij==2
                            for g=1:smij
                                s(g,1)=g;
                                aa((2*F-2)+s(g,1),1)=s(g,1);
                                i=(2*F)+g-2-h;
                                K(i,:)=[];
                                K(:,i)=[];
                                P(i,:)=[];
                                U(i,:)=[];
                                h=h+1;
                            end
                    else
                        if S~=t3 || S~=grid || S~=f2
                            for g=1:smij
                                s(g,1)=input('*** What node is that (i) ? *** ');         
                                if S==t2
                                    j=2;ii=2;
                                end
                                if S==f3
                                    j=6;ii=6;
                                end
                                aa((ii*F-ii)+s(g,1),1)=s(g,1);
                                i=(j*F)+s(g,1)-j-h;
                                K(i,:)=[];
                                K(:,i)=[];
                                P(i,:)=[];
                                U(i,:)=[];
                                h=h+1;
                            end
                        end
                    end
                end
            end
        end
    end
end
K
U=inv(K)*P
h=0;ee=0;ss=0;
for I=1:N
    [I];
    if S==t2
        i=2;j=1;
    elseif S==t3 || S==grid || S==f2
        i=3;j=2;
    elseif S==f3
        i=6;j=5;
    end
    if Ismo==0
        U2(i*I-j:i*I,1)=U(i*I-j:i*I,1);
    else   
        for x=1:I
            F=C(x,1);
            if I==F
                ee=F;
            end
        end
        if I==ee
            for g=1:i
                ss=aa((i*ee-i)+g,1);
                z=i*I-i;
                if ss~=0
                    U2(z+g,1)=0;
                    h=h+1;
                end
                if ss==0
                    oo=U(z+g-h,1);
                    U2(z+g,1)=oo;
                end
            end
        elseif I~=ee
            U2(i*I-j:i*I,1)=U(i*I-j-h:i*I-h,1);
        end
    end
end
for t=1:B
    I=u(t);J=v(t);
    nodes.IJ=[I,J]
    if S==grid || S==f2 || S==f3
        cell2mat(WI(t,I));
        r3(i*I-j:i*I,i*J-j:i*J);
        r1(i*I-j:i*I,i*J-j:i*J)'*cell2mat(WI(t,I))*r2(i*I-j:i*I,i*J-j:i*J);
        U2(i*I-j:i*I,1);
        cell2mat(WI(t,J));
        r2(i*I-j:i*I,i*J-j:i*J); 
        r4(i*I-j:i*I,i*J-j:i*J)'*cell2mat(WI(t,J))*r2(i*I-j:i*I,i*J-j:i*J);
        U2(i*J-j:i*J,1);
        if S==f3 && e==y
            Ppij=q11(i*I-j:i*I,i*J-j:i*J)*cell2mat(WI(t,I))*r2(i*I-j:i*I,i*J-j:i*J)*U2(i*I-j:i*I,1)+q11(i*I-j:i*I,i*J-j:i*J)*cell2mat(WI(t,J))*r3(i*I-j:i*I,i*J-j:i*J)*U2(i*J-j:i*J,1)
            Ppji=q44(i*I-j:i*I,i*J-j:i*J)*cell2mat(WJ(t,I))*r2(i*I-j:i*I,i*J-j:i*J)*U2(i*I-j:i*I,1)+q44(i*I-j:i*I,i*J-j:i*J)*cell2mat(WJ(t,J))*r3(i*I-j:i*I,i*J-j:i*J)*U2(i*J-j:i*J,1)
        end
        pij=cell2mat(WI(t,I))*r2(i*I-j:i*I,i*J-j:i*J)*U2(i*I-j:i*I,1)+cell2mat(WI(t,J))*r3(i*I-j:i*I,i*J-j:i*J)*U2(i*J-j:i*J,1)
        pji=cell2mat(WJ(t,I))*r2(i*I-j:i*I,i*J-j:i*J)*U2(i*I-j:i*I,1)+cell2mat(WJ(t,J))*r3(i*I-j:i*I,i*J-j:i*J)*U2(i*J-j:i*J,1)
        PIJ=r1(i*I-j:i*I,i*J-j:i*J)*cell2mat(WI(t,I))*r2(i*I-j:i*I,i*J-j:i*J)*U2(i*I-j:i*I,1)+r1(i*I-j:i*I,i*J-j:i*J)*cell2mat(WI(t,J))*r3(i*I-j:i*I,i*J-j:i*J)*U2(i*J-j:i*J,1)
        PJI=r4(i*I-j:i*I,i*J-j:i*J)*cell2mat(WJ(t,I))*r2(i*I-j:i*I,i*J-j:i*J)*U2(i*I-j:i*I,1)+r4(i*I-j:i*I,i*J-j:i*J)*cell2mat(WJ(t,J))*r3(i*I-j:i*I,i*J-j:i*J)*U2(i*J-j:i*J,1)
    else
        if S==t2
            r=[l(t) m(t)];
            R=-r;
        elseif S==t3
            r=[l(t) m(t) n(t)];
            R=-r;
        end
        cell2mat(WI(t,I));  
        U2(i*I-j:i*I,1);
        cell2mat(WI(t,J));
        U2(i*J-j:i*J,1);
        pij=cell2mat(WI(t,I))*r*U2(i*I-j:i*I,1)+cell2mat(WI(t,J))*R*U2(i*J-j:i*J,1)  
        pji=cell2mat(WJ(t,I))*r*U2(i*I-j:i*I,1)+cell2mat(WJ(t,J))*R*U2(i*J-j:i*J,1)
        PIJ=r'*cell2mat(WI(t,I))*r*U2(i*I-j:i*I,1)+r'*cell2mat(WI(t,J))*R*U2(i*J-j:i*J,1)
        PJI=R'*cell2mat(WI(t,I))*r*U2(i*I-j:i*I,1)+R'*cell2mat(WJ(t,J))*R*U2(i*J-j:i*J,1)
    end
end

Contact us