/*
 * This code is generated by BioUML FrameWork 
 * for BIOMD0000000128.xml diagram  at 2008.03.20 15:09:05
 */
import biouml.plugins.simulation.ae.NewtonSolver;
import biouml.plugins.simulation.java.JavaBaseModel;
import ru.biosoft.math.MathRoutines;

public class BIOMD0000000128 extends JavaBaseModel
{

/*
 * Write rules to calculate equation parameters
 */
    private void __internalVarInitFunc_0(double time, double[] x)
    {
        minf = 1/(1 + Math.exp((vm - x[3])/sm));
        ik = gk*x[6]*(x[3] - vk);
        ainf = 1/(1 + dact/x[0]);
        hinfer = 1/(1 + x[0]/dinh);
        jerp = kserca*x[0];
        binf = IP3/(IP3 + dip3);
        o = Math.pow(ainf, 3)*Math.pow(binf, 3)*Math.pow(hinfer, 3);
        jmemtot = -(alpha*ica + kc*x[0]);
        jerleak = perl*(x[2] - x[0]);
        jerip3 = o*(x[2] - x[0]);
        jertot = jerleak + jerip3 - jerp;
        ninf = 1/(1 + Math.exp((vn - x[3])/sn));
        perl_inf = x[5]*x[1]*Math.pow(x[0], 4)/(Math.pow(ki, 4) + Math.pow(x[0], 4));
        hinf = 1/(1 + Math.exp((vh - x[3])/sh));
        ica = gca*minf*(x[3] - vca);
        igirk = girk*x[4]*(x[3] - vk);
    }


/*
 * Write rules to calculate equation parameters excluding internal variables.
 */
    public void __internalRateVarInitFunc_0(double time, double[] x)
    {
        rate_reaction_0000001 = cell*f*(jertot + jmemtot);
        rate_reaction_0000002 = -fer*sigmav*jertot*cell;
        rate_reaction_000003 = cell*ETswitch*((cAMPlow - x[1])/taudir);
    }

    public void Init()
    {
        initialValues = getInitialValues();
/*
 * Initialize variables
 */
        cell = 1.0; // initial value of $cell
        alpha = 4.5E-6; // initial value of alpha
        cAMPlow = 0.2; // initial value of cAMPlow
        cm = 5300.0; // initial value of cm
        dact = 0.35; // initial value of dact
        dinh = 0.4; // initial value of dinh
        dip3 = 0.5; // initial value of dip3
        f = 0.01; // initial value of f
        fer = 0.01; // initial value of fer
        gca = 2000.0; // initial value of gca
        girk = 1000.0; // initial value of girk
        gk = 3500.0; // initial value of gk
        kc = 0.15; // initial value of kc
        ki = 0.5; // initial value of ki
        kserca = 0.4; // initial value of kserca
        lambda = 1.25; // initial value of lambda
        perl = 5.0E-4; // initial value of perl
        sh = 70.0; // initial value of sh
        sigmav = 10.0; // initial value of sigmav
        sm = 12.0; // initial value of sm
        sn = 5.0; // initial value of sn
        taudir = 20000.0; // initial value of taudir
        tauh = 20.0; // initial value of tauh
        taun = 20.0; // initial value of taun
        vca = 25.0; // initial value of vca
        vh = -20.0; // initial value of vh
        vk = -75.0; // initial value of vk
        vm = -20.0; // initial value of vm
        vn = -16.0; // initial value of vn
    }

    /*
     * Model variables initial values
     */
    protected double rate_reaction_0000001;
    protected double rate_reaction_0000002;
    protected double rate_reaction_000003;
    protected double cell;
    protected double ETswitch;
    protected double IP3;
    protected double ainf;
    protected double alpha;
    protected double binf;
    protected double cAMPlow;
    protected double cm;
    protected double dact;
    protected double dinh;
    protected double dip3;
    protected double f;
    protected double fer;
    protected double gca;
    protected double girk;
    protected double gk;
    protected double hinf;
    protected double hinfer;
    protected double ica;
    protected double igirk;
    protected double ik;
    protected double jerip3;
    protected double jerleak;
    protected double jerp;
    protected double jertot;
    protected double jmemtot;
    protected double kc;
    protected double ki;
    protected double kserca;
    protected double lambda;
    protected double minf;
    protected double ninf;
    protected double o;
    protected double perl;
    protected double perl_inf;
    protected double sh;
    protected double sigmav;
    protected double sm;
    protected double sn;
    protected double taudir;
    protected double tauh;
    protected double taun;
    protected double time;
    protected double vca;
    protected double vh;
    protected double vk;
    protected double vm;
    protected double vn;

    public double[] extendResult(double time,double [] x)
    {
        this.time = time;

        __internalVarInitFunc_0(time, x);

        double[] y = new double[26];
        y[0] = x[0];
        y[1] = x[1];
        y[2] = x[2];
        y[3] = ETswitch;
        y[4] = IP3;
        y[5] = x[3];
        y[6] = ainf;
        y[7] = binf;
        y[8] = girk;
        y[9] = x[4];
        y[10] = hinf;
        y[11] = hinfer;
        y[12] = ica;
        y[13] = igirk;
        y[14] = ik;
        y[15] = x[5];
        y[16] = jerip3;
        y[17] = jerleak;
        y[18] = jerp;
        y[19] = jertot;
        y[20] = jmemtot;
        y[21] = minf;
        y[22] = x[6];
        y[23] = ninf;
        y[24] = o;
        y[25] = perl_inf;
        return y;
    }
    public double[] getInitialValues()
    {
        double [] x = new double[7];
        this.time = 0.0;
        x[0] = 0.3; // - $"cell.c"
        x[1] = 1.0; // - $"cell.cAMP"
        x[2] = 260.0; // - $"cell.cer"
        x[3] = -60.0; // - V
        x[4] = 0.0; // - h
        x[5] = 1.0; // - inh
        x[6] = 0.0; // - n

        __internalVarInitFunc_0(time, x);
        __internalRateVarInitFunc_0(time, x);

        return x;
    }

/*
 * code for algebraic rules calculations
 */

/*
 * end of code for algebraic rules calculations
 */

    protected void calculateRates(double time, double[] x)
    {

        __internalVarInitFunc_0(time, x);
        __internalRateVarInitFunc_0(time, x);

    }

        /*
         * calculate dy/dt for 'BIOMD0000000128.xml' model
         */
    public void __internalDyDt_0(double time, double [] x, double[] result)
    {
        result[0] = +rate_reaction_0000001;
        result[1] = +rate_reaction_000003;
        result[2] = +rate_reaction_0000002;
        result[3] = +(-ica - ik - igirk)/cm;
        result[4] = +(hinf - x[4])/tauh;
        result[5] = +ETswitch*((0.2 - x[5])/taudir);
        result[6] = +lambda*(ninf - x[6])/taun;
    }
    protected double [] calculateResult(double time, double[] x)
    {
        double[] result = new double[7];
        __internalDyDt_0(time, x, result);
        return result;
    }
    public double[] dy_dt(double time, double[] x)
    {
        this.time = time;
        calculateRates( time,x );

        return calculateResult( time,x );
    }

    public double[] checkEvent(double time, double[] x)
    {
        this.time = time;

        __internalVarInitFunc_0(time, x);
        __internalRateVarInitFunc_0(time, x);

        double [] flags = new double[1];
        flags[0] = (time > 60000) ? +1 : -1; //event_0000001
        return flags;
    }

    public void processEvent(int __internalVar12060041455150, double time, double[] x)
    {
        this.time = time;
        if ( __internalVar12060041455150 == 0) //event_0000001
        {
            IP3 = 0.3;
            girk = 3000;
            ETswitch = 1;
        }
    }

} // class ...