module testMC36SM
using LinearAlgebra
export run
function run()
    println("starting..")
    vj = LinRange(0, 120, 100);
    gj = zeros(length(vj))
    parAllFast = [ .15 28 20e-9 20e-9/10 800 50 .5 -1;       .15 28 20e-9 20e-9/10 800 50 .5 1 ];
    parAllSlow = [ .1 36 100e-9 1e-15 1e10 1e10 .5 -1;       .1 36 100e-9 1e-15 1e10 1e10 .5 1 ];
    par = [parAllFast; parAllSlow];
    pc1c2 = 0;
    pc2c1 = 0;
    pc1c2 = 0.01;
    pc2c1 = 0.001;
    ppp = zeros(36);
    ppp[1] = 1;
    for i = 1:length(vj)
        gj[i] = MC36SM_Mindaugo_SS(vj[i], par, ppp, pc1c2, pc2c1);
    end
    println("done.")
    return vj, gj
end

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]';;

   K = Pt[1];
   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

   for a = 1:4 #// voltages and rectification
       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
       for j = 1:4
           k[i,j] = exp(A[j] * (v[i,j] * P[j] - v0[j]));
       end
   end

   # //------------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

   d_g = 100000;
   g_old = 100000;
   g_final = 0;
   i = 0;
   p_ = copy(p) #36-vec
   ggg = 1 ./ sum( (1 ./ g), dims=2 ); #%sums rows
   while (d_g / (g_old + g_final) > .001e-5)
       i = i + 1;
       #p=p*PPP;
       mul!(p_, PPP', p)
       p, p_ = p_, p

       g_final = dot(ggg, p) #(p*ggg)[1];
       d_g = abs(g_old - g_final);
       g_old = g_final;
   end
   return g_final;
end
end