mc36smLib.cpp

#include "GjLib.hpp"
#include <math.h>

using namespace std;

int add(int a, int b)
{
    return (a + b);
}

double MC36SM_Mindaugo_SS(double vj, double *par, double *p, double Pt, double pc1c2, double pc2c1)
{
    double *A = par;
    double *v0 = par + 1 * 4;
    double *Go = par + 2 * 4;
    double *Gc = par + 3 * 4;
    double *Ro = par + 4 * 4;
    double *Rc = par + 5 * 4;
    double *P = par + 6 * 4;

    double g[36][4] = {};// initialize to 0
    double v[36][4] = {};
    double k[36][4] = {};
    double R[36][4] = {};

    //states conductances
    for (int i1 = 0; i1 < 2; i1++)
        for (int i2 = 0; i2 < 3; i2++)
            for (int i3 = 0; i3 < 3; i3++)
                for (int i4 = 0; i4 < 2; i4++)
                {
                    if (i1 == 0)
                    {
                        g[i1 * 18 + i2 * 6 + i3 * 2 + i4 * 1][0] = Go[0];
                        R[i1 * 18 + i2 * 6 + i3 * 2 + i4 * 1][0] = Ro[0];
                    }
                    else
                    {
                        g[i1 * 18 + i2 * 6 + i3 * 2 + i4 * 1][0] = Gc[0];
                        R[i1 * 18 + i2 * 6 + i3 * 2 + i4 * 1][0] = Rc[0];
                    }
                    if (i2 == 0)
                    {
                        g[i1 * 18 + i2 * 6 + i3 * 2 + i4 * 1][1] = Go[1];
                        R[i1 * 18 + i2 * 6 + i3 * 2 + i4 * 1][1] = Ro[1];
                    }
                    else
                    {
                        g[i1 * 18 + i2 * 6 + i3 * 2 + i4 * 1][1] = Gc[1];
                        R[i1 * 18 + i2 * 6 + i3 * 2 + i4 * 1][1] = Rc[1];
                    }
                    if (i3 == 0)
                    {
                        g[i1 * 18 + i2 * 6 + i3 * 2 + i4 * 1][2] = Go[2];
                        R[i1 * 18 + i2 * 6 + i3 * 2 + i4 * 1][2] = Ro[2];
                    }
                    else
                    {
                        g[i1 * 18 + i2 * 6 + i3 * 2 + i4 * 1][2] = Gc[2];
                        R[i1 * 18 + i2 * 6 + i3 * 2 + i4 * 1][2] = Rc[2];
                    }
                    if (i4 == 0)
                    {
                        g[i1 * 18 + i2 * 6 + i3 * 2 + i4 * 1][3] = Go[3];
                        R[i1 * 18 + i2 * 6 + i3 * 2 + i4 * 1][3] = Ro[3];
                    }
                    else
                    {
                        g[i1 * 18 + i2 * 6 + i3 * 2 + i4 * 1][3] = Gc[3];
                        R[i1 * 18 + i2 * 6 + i3 * 2 + i4 * 1][3] = Rc[3];
                    }

                }

    double h[36][4] = {};
    for (int i = 0; i < 4; i++)     for (int j = 0; j < 36; j++)        h[j][i] = g[j][i];

    for (int a = 0; a < 4; a++) // voltages and rectification
    {
        double gg[36] = {};
        for (int i = 0; i < 36; i++)        //gg = 1. / sum(1. / g, 2); %sums rows
        {
            for (int j = 0; j < 4; j++)
                gg[i] += 1 / g[i][j];
            gg[i] = 1 / gg[i];
        }
        for (int i = 0; i < 4; i++) //v=vj*[gg gg gg gg]./g;
            for (int j = 0; j < 36; j++)
                v[j][i] = vj * gg[j] / g[j][i];
        for (int i = 0; i < 4; i++) //g=g.*exp(v./R);
            for (int j = 0; j < 36; j++)
                g[j][i] = h[j][i] * exp(v[j][i] / R[j][i]);
    }

    for (int i = 0; i < 36; i++) //k=exp(repmat(A,16,1).*(v.*repmat(P,16,1)-repmat(v0,16,1)));
        for (int j = 0; j < 4; j++)
            k[i][j] = exp(A[j] * (v[i][j] * P[j] - v0[j]));

    double K = Pt;
    double PPP[36][36] = {};
    //------------Matrica P-------------------------------------------------
    double l = pc1c2 / pc2c1;
    double p1, p2, p3, p4;
    for (int i1 = 0; i1 < 2; i1++)
        for (int i2 = 0; i2 < 3; i2++)
            for (int i3 = 0; i3 < 3; i3++)
                for (int i4 = 0; i4 < 2; i4++)
                    for (int j1 = 0; j1 < 2; j1++)
                        for (int j2 = 0; j2 < 3; j2++)
                            for (int j3 = 0; j3 < 3; j3++)
                                for (int j4 = 0; j4 < 2; j4++)
                                {
                                    int i = i1 * 18 + i2 * 6 + i3 * 2 + i4 * 1;
                                    int j = j1 * 18 + j2 * 6 + j3 * 2 + j4 * 1;
                                    ////////////// p1
                                    if (i1 == 0)
                                    {
                                        if (j1 == 0)
                                            p1 = (1 - K * k[i][0] / (1 + k[i][0]));
                                        else
                                            p1 = K * k[i][0] / (1 + k[i][0]);
                                    }
                                    else
                                    {
                                        if (j1 == 1)
                                            p1 = (1 - K / (1 + k[i][0]));
                                        else
                                            p1 = K / (1 + k[i][0]);
                                    }
                                    ///////////////p2
                                    if (i2 == 0)
                                    {
                                        if (j2 == 0)
                                            p2 = (1 - K * k[i][1] / (1 + k[i][1]));
                                        if (j2 == 1)
                                            p2 = K * k[i][1] / (1 + k[i][1]);
                                        if (j2 == 2)
                                            p2 = 0;
                                    }
                                    else if (i2 == 1)
                                    {
                                        if (j2 == 0)
                                            p2 = (K / (1 + k[i][1]))*(1 - pc1c2);
                                        if (j2 == 1)
                                            p2 = (1 - K / (1 + k[i][1]))*(1 - pc1c2);
                                        if (j2 == 2)
                                            p2 = pc1c2;
                                    }
                                    else if (i2 == 2)
                                    {
                                        if (j2 == 0)
                                            p2 = 0;
                                        if (j2 == 1)
                                            p2 = pc2c1;
                                        if (j2 == 2)
                                            p2 = 1 - pc2c1;
                                    }
                                    ////////////////p3
                                    if (i3 == 0)
                                    {
                                        if (j3 == 0)
                                            p3 = (1 - K * k[i][2] / (1 + k[i][2]));
                                        if (j3 == 1)
                                            p3 = K * k[i][2] / (1 + k[i][2]);
                                        if (j3 == 2)
                                            p3 = 0;
                                    }
                                    else if (i3 == 1)
                                    {
                                        if (j3 == 0) //c1->o
                                            p3 = (K / (1 + k[i][2]))*(1 - pc1c2);
                                        if (j3 == 1)
                                            p3 = (1 - K / (1 + k[i][2]))*(1 - pc1c2);
                                        if (j3 == 2)
                                            p3 = pc1c2;
                                    }
                                    else if (i3 == 2)
                                    {
                                        if (j3 == 0)
                                            p3 = 0;
                                        if (j3 == 1)
                                            p3 = pc2c1;
                                        if (j3 == 2)
                                            p3 = 1 - pc2c1;
                                    }
                                    /////////////////p4
                                    if (i4 == 0)
                                    {
                                        if (j4 == 0)
                                            p4 = (1 - K * k[i][3] / (1 + k[i][3]));
                                        else
                                            p4 = K * k[i][3] / (1 + k[i][3]);
                                    }
                                    else
                                    {
                                        if (j4 == 1)
                                            p4 = (1 - K / (1 + k[i][3]));
                                        else
                                            p4 = K / (1 + k[i][3]);
                                    }
                                    PPP[i][j] = p1 * p2*p3*p4;
                                }

    // -------------- SS skaiciavimas
    double d_g = 100000;
    double g_old = 100000;
    double g_final = 0;
    int i = 0;
    while (d_g / (g_old + g_final) > .001e-5)
    {
        i++;
        double q[36] = {};
        for (int j = 0; j < 36; j++) for (int i = 0; i < 36; i++)    q[j] += p[i] * PPP[i][j]; // q=p*PPP;
        double ggg[36] = {};
        for (int i = 0; i < 36; i++)        //ggg = 1. / sum(1. / g, 2); %sums rows
        {
            for (int j = 0; j < 4; j++) ggg[i] += 1 / g[i][j];
            ggg[i] = 1 / ggg[i];
        }

        g_final = 0;
        for (int j = 0; j < 36; j++) //g_final = q*ggg';
            g_final += q[j] * ggg[j];

        for (int i = 0; i < 36; i++)    p[i] = q[i]; //grazina p

        d_g = fabs(g_old - g_final);
        g_old = g_final;
    }
    return g_final;
}

// Returns gj
double* MC36SM_Mindaugo_SS_t(int n, double *vj, double *par, double *p, double Pt, double pc1c2, double pc2c1) {

}