Advertisement
kestutisma

mc36sm_ss.jl

Apr 10th, 2019
1,343
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Julia 5.66 KB | None | 0 0
  1. export MC36SM_Mindaugo_SS
  2. function MC36SM_Mindaugo_SS(vj::Float64, par::Array{Float64,2}, p::Array{Float64}, pc1c2::Float64, pc2c1::Float64)
  3.     PPP = zeros(36,36);
  4.     p1 = 0;
  5.     p2 = 0;
  6.     p3 = 0;
  7.     p4 = 0;
  8.  
  9.     A = par[:,1]';
  10.    v0 = par[:,2]';
  11.     Go = par[:,3]';
  12.    Gc = par[:,4]';
  13.     Ro = par[:,5]';
  14.    Rc = par[:,6]';
  15.     Pt = par[:,7]';
  16.    P = par[:,8]';;
  17.  
  18.    g = zeros(36,4);
  19.    v = zeros(36,4);
  20.    k = zeros(36,4);
  21.    R = zeros(36,4);
  22.  
  23.    # //states conductances
  24.    for i1=1:2
  25.        for i2=1:3
  26.            for i3=1:3
  27.                for i4=1:2
  28.                    if (i1==1)
  29.                        g[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,1]=Go[1];
  30.                        R[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,1]=Ro[1];
  31.                    else
  32.                        g[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,1]=Gc[1];
  33.                        R[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,1]=Rc[1];
  34.                    end
  35.                    if (i2==1)
  36.                        g[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,2]=Go[2];
  37.                        R[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,2]=Ro[2];
  38.                    else
  39.                        g[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,2]=Gc[2];
  40.                        R[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,2]=Rc[2];
  41.                    end
  42.                    if (i3==1)
  43.                        g[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,3]=Go[3];
  44.                        R[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,3]=Ro[3];
  45.                    else
  46.                        g[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,3]=Gc[3];
  47.                        R[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,3]=Rc[3];
  48.                    end
  49.                    if (i4==1)
  50.                        g[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,4]=Go[4];
  51.                        R[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,4]=Ro[4];
  52.                    else
  53.                        g[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,4]=Gc[4];
  54.                        R[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,4]=Rc[4];
  55.                    end
  56.                end
  57.            end
  58.        end
  59.    end
  60.  
  61.    h=g;
  62.  
  63.    for a = 1:4 #// voltages and rectification
  64.     #    @bp
  65.        gg = 1 ./ sum( (1 ./ g), dims = 2 ); #sums rows
  66.        v=vj*[gg gg gg gg]./g;
  67.        g=g .* exp.(v./R);
  68.    end
  69.  
  70.    for i = 1:36 #//k=exp(repmat(A,16,1).*(v.*repmat(P,16,1)-repmat(v0,16,1)));
  71.        for j = 1:4
  72.             # @bp
  73.            k[i,j] = exp(A[j] * (v[i,j] * P[j] - v0[j]));
  74.        end
  75.    end
  76.  
  77.    K = Pt[1];
  78.  
  79.    # //------------Matrica P-------------------------------------------------
  80.    for i1=1:2 for i2=1:3 for i3=1:3 for i4=1:2
  81.     for j1=1:2 for j2=1:3 for j3=1:3 for j4=1:2
  82.       i=(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1;
  83.       j=(j1-1)*18+(j2-1)*6+(j3-1)*2+(j4-1)*1+1;
  84.  
  85.  
  86.       # %-------be Pt---------------------
  87.       # % p1
  88.        if (i1==1)
  89.            if (j1==1)
  90.                p1=(1-K*k[i,1]/(1+k[i,1]));
  91.            else
  92.                p1=K*k[i,1]/(1+k[i,1]);
  93.            end
  94.        else
  95.            if (j1==2)
  96.                p1=(1-K/(1+k[i,1]));
  97.            else
  98.                p1=K/(1+k[i,1]);
  99.            end
  100.        end
  101.    # %------------------------------
  102.    #   % p2
  103.        if (i2==1)
  104.            if (j2==1)
  105.                p2=(1-K*k[i,2]/(1+k[i,2]));
  106.            end
  107.            if (j2==2)
  108.                p2=K*k[i,2]/(1+k[i,2]);
  109.            end
  110.            if (j2==3)
  111.                p2=0;
  112.            end
  113.        end
  114.  
  115.        if (i2==2)
  116.            if (j2==1)
  117.                p2=(K/(1+k[i,2]))*(1-pc1c2);
  118.            end
  119.            if (j2==2)
  120.                p2=(1-K/(1+k[i,2]))*(1-pc1c2);
  121.            end
  122.            if (j2==3)
  123.                p2=pc1c2;
  124.            end
  125.        end
  126.  
  127.        if (i2==3)
  128.            if (j2==1)
  129.                p2=0;
  130.            end
  131.            if (j2==2)
  132.                p2=pc2c1;
  133.            end
  134.            if (j2==3)
  135.                p2=1-pc2c1;
  136.            end
  137.        end
  138.      # %----------------------------------------
  139.      # % p3
  140.        if (i3==1)
  141.            if (j3==1)
  142.                p3=(1-K*k[i,3]/(1+k[i,3]));
  143.            end
  144.            if (j3==2)
  145.                p3=K*k[i,3]/(1+k[i,3]);
  146.            end
  147.            if (j3==3)
  148.                p3=0;
  149.            end
  150.        end
  151.  
  152.        if (i3==2)
  153.            if (j3==1)
  154.                p3=(K/(1+k[i,3]))*(1-pc1c2);
  155.            end
  156.            if (j3==2)
  157.                p3=(1-K/(1+k[i,3]))*(1-pc1c2);
  158.            end
  159.            if (j3==3)
  160.                p3=pc1c2;
  161.            end
  162.        end
  163.  
  164.        if (i3==3)
  165.            if (j3==1)
  166.                p3=0;
  167.            end
  168.            if (j3==2)
  169.                p3=pc2c1;
  170.            end
  171.            if (j3==3)
  172.                p3=1-pc2c1;
  173.            end
  174.        end
  175.      # %----------------------------
  176.      # %p4
  177.        if (i4==1)
  178.            if (j4==1)
  179.                p4=(1-K*k[i,4]/(1+k[i,4]));
  180.            else
  181.                p4=K*k[i,4]/(1+k[i,4]);
  182.            end
  183.        else
  184.            if (j4==2)
  185.                p4=(1-K/(1+k[i,4]));
  186.            else
  187.                p4=K/(1+k[i,4]);
  188.            end
  189.        end
  190.  
  191.     PPP[i,j]=p1*p2*p3*p4;
  192.  
  193.    end
  194.    end
  195.    end
  196.    end
  197.    end
  198.    end
  199.    end
  200.    end
  201.  
  202.    # // -------------- SS skaiciavimas
  203.    d_g = 100000;
  204.    g_old = 100000;
  205.    g_final = 0;
  206.    i = 0;
  207.    while (d_g / (g_old + g_final) > .001e-5)
  208.        i = i + 1;
  209.        q=p*PPP;
  210.        ggg = 1 ./ sum( (1 ./ g), dims=2 ); #%sums rows
  211.        g_final = (q*ggg)[1]; #[1] nes ggg array2d tipo, nors turi tik 1 reiksme
  212.  
  213.        # for i = 0; i < 36; i++)    p[i] = q[i]; //grazina p
  214.        p .= q;
  215.  
  216.        d_g = abs(g_old - g_final);
  217.        g_old = g_final;
  218.    end
  219.    return g_final;
  220. end
  221. # end
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement