mc36sm_ss.jl

export MC36SM_Mindaugo_SS
function MC36SM_Mindaugo_SS(vj::Float64, par::Array{Float64,2}, p::Array{Float64}, pc1c2::Float64, pc2c1::Float64)
    PPP = zeros(36,36);
    p1 = 0;
    p2 = 0;
    p3 = 0;
    p4 = 0;

    A = par[:,1]';
    v0 = par[:,2]';
    Go = par[:,3]';
    Gc = par[:,4]';
    Ro = par[:,5]';
    Rc = par[:,6]';
    Pt = par[:,7]';
    P = par[:,8]';;

   g = zeros(36,4);
   v = zeros(36,4);
   k = zeros(36,4);
   R = zeros(36,4);

   # //states conductances
   for i1=1:2
       for i2=1:3
           for i3=1:3
               for i4=1:2
                   if (i1==1)
                       g[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,1]=Go[1];
                       R[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,1]=Ro[1];
                   else
                       g[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,1]=Gc[1];
                       R[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,1]=Rc[1];
                   end
                   if (i2==1)
                       g[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,2]=Go[2];
                       R[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,2]=Ro[2];
                   else
                       g[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,2]=Gc[2];
                       R[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,2]=Rc[2];
                   end
                   if (i3==1)
                       g[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,3]=Go[3];
                       R[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,3]=Ro[3];
                   else
                       g[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,3]=Gc[3];
                       R[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,3]=Rc[3];
                   end
                   if (i4==1)
                       g[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,4]=Go[4];
                       R[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,4]=Ro[4];
                   else
                       g[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,4]=Gc[4];
                       R[(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1,4]=Rc[4];
                   end
               end
           end
       end
   end

   h=g;

   for a = 1:4 #// voltages and rectification
    #    @bp
       gg = 1 ./ sum( (1 ./ g), dims = 2 ); #sums rows
       v=vj*[gg gg gg gg]./g;
       g=g .* exp.(v./R);
   end

   for i = 1:36 #//k=exp(repmat(A,16,1).*(v.*repmat(P,16,1)-repmat(v0,16,1)));
       for j = 1:4
            # @bp
           k[i,j] = exp(A[j] * (v[i,j] * P[j] - v0[j]));
       end
   end

   K = Pt[1];

   # //------------Matrica P-------------------------------------------------
   for i1=1:2 for i2=1:3 for i3=1:3 for i4=1:2
    for j1=1:2 for j2=1:3 for j3=1:3 for j4=1:2
      i=(i1-1)*18+(i2-1)*6+(i3-1)*2+(i4-1)*1+1;
      j=(j1-1)*18+(j2-1)*6+(j3-1)*2+(j4-1)*1+1;


      # %-------be Pt---------------------
      # % p1
       if (i1==1)
           if (j1==1)
               p1=(1-K*k[i,1]/(1+k[i,1]));
           else
               p1=K*k[i,1]/(1+k[i,1]);
           end
       else
           if (j1==2)
               p1=(1-K/(1+k[i,1]));
           else
               p1=K/(1+k[i,1]);
           end
       end
   # %------------------------------
   #   % p2
       if (i2==1)
           if (j2==1)
               p2=(1-K*k[i,2]/(1+k[i,2]));
           end
           if (j2==2)
               p2=K*k[i,2]/(1+k[i,2]);
           end
           if (j2==3)
               p2=0;
           end
       end

       if (i2==2)
           if (j2==1)
               p2=(K/(1+k[i,2]))*(1-pc1c2);
           end
           if (j2==2)
               p2=(1-K/(1+k[i,2]))*(1-pc1c2);
           end
           if (j2==3)
               p2=pc1c2;
           end
       end

       if (i2==3)
           if (j2==1)
               p2=0;
           end
           if (j2==2)
               p2=pc2c1;
           end
           if (j2==3)
               p2=1-pc2c1;
           end
       end
     # %----------------------------------------
     # % p3
       if (i3==1)
           if (j3==1)
               p3=(1-K*k[i,3]/(1+k[i,3]));
           end
           if (j3==2)
               p3=K*k[i,3]/(1+k[i,3]);
           end
           if (j3==3)
               p3=0;
           end
       end

       if (i3==2)
           if (j3==1)
               p3=(K/(1+k[i,3]))*(1-pc1c2);
           end
           if (j3==2)
               p3=(1-K/(1+k[i,3]))*(1-pc1c2);
           end
           if (j3==3)
               p3=pc1c2;
           end
       end

       if (i3==3)
           if (j3==1)
               p3=0;
           end
           if (j3==2)
               p3=pc2c1;
           end
           if (j3==3)
               p3=1-pc2c1;
           end
       end
     # %----------------------------
     # %p4
       if (i4==1)
           if (j4==1)
               p4=(1-K*k[i,4]/(1+k[i,4]));
           else
               p4=K*k[i,4]/(1+k[i,4]);
           end
       else
           if (j4==2)
               p4=(1-K/(1+k[i,4]));
           else
               p4=K/(1+k[i,4]);
           end
       end

    PPP[i,j]=p1*p2*p3*p4;

   end
   end
   end
   end
   end
   end
   end
   end

   # // -------------- SS skaiciavimas
   d_g = 100000;
   g_old = 100000;
   g_final = 0;
   i = 0;
   while (d_g / (g_old + g_final) > .001e-5)
       i = i + 1;
       q=p*PPP;
       ggg = 1 ./ sum( (1 ./ g), dims=2 ); #%sums rows
       g_final = (q*ggg)[1]; #[1] nes ggg array2d tipo, nors turi tik 1 reiksme

       # for i = 0; i < 36; i++)    p[i] = q[i]; //grazina p
       p .= q;

       d_g = abs(g_old - g_final);
       g_old = g_final;
   end
   return g_final;
end
# end