Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- 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
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement