Advertisement
kestutisma

test

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