Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- module testMC36SM
- 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(1, 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;
- while (d_g / (g_old + g_final) > .001e-5)
- i = i + 1;
- p=p*PPP;
- ggg = 1 ./ sum( (1 ./ g), dims=2 ); #%sums rows
- g_final = (p*ggg)[1];
- 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