SHSystem.m

Go to the documentation of this file.

namespace models{
namespace motorunit{


/* (Autoinserted by mtoc++)
 * This source code has been filtered by the mtoc++ executable,
 * which generates code that can be processed by the doxygen documentation tool.
 *
 * On the other hand, it can neither be interpreted by MATLAB, nor can it be compiled with a C++ compiler.
 * Except for the comments, the function bodies of your M-file functions are untouched.
 * Consequently, the FILTER_SOURCE_FILES doxygen switch (default in our Doxyfile.template) will produce
 * attached source files that are highly readable by humans.
 *
 * Additionally, links in the doxygen generated documentation to the source code of functions and class members refer to
 * the correct locations in the source code browser.
 * However, the line numbers most likely do not correspond to the line numbers in the original MATLAB source files.
 */

class SHSystem
  :public models.BaseFirstOrderSystem {
    public: /* ( Transient ) */

        d = "[]";
    public: /* ( Constant ) */

        static const dm = 6;
        static const dsa = 56;
        static const .double C_m_slow = 0.58;
        static const C_m_fast = 1;


        static const .rowvec<double> ForceOutputScalingPolyCoeff = 1e5*"[0.072141692014344  -0.456658053996476   1.274105896595146  -2.041231359447655   2.058191740934228  -1.353159489096666   0.584711952400098  -0.164458365881107   0.029358733675881  -0.003174883073694   0.000165817182406    0.000048348565452   0.000012193246968]";
    public:

        noiseGen;
    public:

        upperlimit_poly;
    private: /* ( Transient ) */

        peakstart = false;

        
    public: /* ( Transient ) */

        SHSystem(models.BaseFullModel model) {
            this = this@models.BaseFirstOrderSystem(model);

            this.addParam(" fibre_type ", 0, " Range ", [0 1]);
            this.addParam(" mean_current_factor ", 3, " Range ", [0 9]);

            this.NumStateDofs= this.dm+this.dsa;
            this.updateDimensions;

            /* % Set system components
             * Core nonlinearity */
            this.f= models.motorunit.SHDynamics(this);

            /*  Load mean current limiting polynomial */
            s = load(models.motoneuron.Model.FILE_UPPERLIMITPOLY);
            this.upperlimit_poly= s.upperlimit_poly;

            /*  Setup noise input */
            ng = models.motoneuron.NoiseGenerator;
            this.noiseGen= ng;
            this.Inputs[1] = @ng.getInput;

            /*  Linear input B for motoneuron */
            this.B= this.assembleB;

            /*  Affine-Linear output C */
            this.C= this.assembleC;

            /*  Constant initial values
             * See also: models.motorunit.experiment.InitialConditions */
            if model.DynamicInitialConditions
                this.x0= this.assembleX0;
            else
                this.x0= dscomponents.ConstInitialValue(this.getConstInitialStates);
            end

            this.updateSparsityPattern;
        }
        function  prepareSimulation(colvec<double> mu,integer inputidx) {

            /*  Limit mean current depending on fibre type */
            mu(2) = min(polyval(this.upperlimit_poly,mu(1)),mu(2));

            prepareSimulation@models.BaseFirstOrderSystem(this, mu, inputidx);
        }
        function  setConfig(colvec<double> mu,integer inputidx) {
            setConfig@models.BaseFirstOrderSystem(this, mu, inputidx);
            this.noiseGen.setFibreType(mu(1));
            this.MaxTimestep= this.Model.dt*100;
            /*  Helper method to reset the "peak counter" of the dynamics. */
            if this.Model.SinglePeakMode
                this.peakstart= false;
            end
        }
        function status = singlePeakModeOutputFcn(unused1,matrix<double> y,flag) {
            if ~strcmp(flag," done ")
                if y(this.dm+1) > -20 && ~this.peakstart
                    this.peakstart= true;
                elseif y(this.dm+1) < -20 && this.peakstart
                    this.peakstart= false;
                    this.f.hadpeak= true;
                end
            end
            status = 0;
        }

    
    private: /* ( Transient ) */


        function x0 = assembleX0() {
            mc = metaclass(this);
            s = load(fullfile(fileparts(which(mc.Name))," x0coeff.mat "));
            x0 = dscomponents.AffineInitialValue;
            m = size(s.coeff,1);
            for k=1:m
                x0.addMatrix(sprintf(" polyval([%s],mu(1)) ",...
                    sprintf(" %g  ",s.coeff(k,:))),full(sparse(k,1,1,m,1)));
            end
        }
        function B = assembleB() {
            B = dscomponents.AffLinInputConv;
            /*  Base noise input mapping */
            B.addMatrix(" 1./(pi*(exp(log(100)*mu(1,:))*3.55e-05 + 77.5e-4).^2) ",...
                sparse(2,1,1,this.dm+this.dsa,2));
            /*  Independent noise input mapping with µ_2 as mean current factor */
            B.addMatrix(" mu(2,:)./(pi*(exp(log(100)*mu(1,:))*3.55e-05 + 77.5e-4).^2) ",...
                sparse(2,2,1,this.dm+this.dsa,2));
        }
        function C = assembleC() {
            C = dscomponents.AffLinOutputConv;
            /*  Extract V_s */
            C.addMatrix(" 1 ",sparse(1,7,1,2,this.dm+this.dsa));
            /*  Add scaling for force output A_s */
            str = " 1 ";
            if this.Model.SingleTwitchOutputForceScaling
                str = [" polyval([ " sprintf(" %g  ",this.ForceOutputScalingPolyCoeff) " ],mu(1)) "];
            end
            C.addMatrix(str, sparse(2,59,1,2,this.dm+this.dsa));
        }
        function x0 = getConstInitialStates() {
            x0_sarc = zeros(56,1);
            x0_sarc(1:14) = [-79.974; -80.2; 5.9; 150.9; 5.9; 12.7; 133;...
                133; 0.009466; 0.9952; 0.0358; 0.4981; 0.581; 0.009466];
            x0_sarc(15:37) = [0.9952; 0.0358; 0.4981; 0.581; 0; 0; 0; 0; 0; 1;...
                0; 0; 0; 0; 0.1; 1500; 0.1; 1500; 25; 615; 615; 811; 811];
            x0_sarc(38:45) = [16900; 16900; 0.4; 0.4; 7200; 7200; 799.6; 799.6];
            x0_sarc(46:54) = [1000; 1000; 3; 0.8; 1.2; 3; 0.3; 0.23; 0.23];
            x0_sarc(55:56) = [0.23; 0.23]; /* 0.0; 0.05 */

            x0 = [zeros(6,1); x0_sarc];
        }
};
}
}