PCDSystem1D.m

Go to the documentation of this file.

namespace models{
namespace pcd{


/* (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 PCDSystem1D
  :public models.pcd.BasePCDSystem {
    public:

        PCDSystem1D(models.BaseFullModel model) {
            this = this@models.pcd.BasePCDSystem(model);

            /*  Set core function */
            this.f= models.pcd.CoreFun1D(this);

            /*  Spatial resolution (in real sizes) */
            this.Omega= [0 1] * this.Model.L;
            this.h= this.Model.L/24;

            /*  Remove params set in BasePCDSystem (simpler for 1D) */
            this.Params= data.ModelParam.empty;

            /*  Add input param (is getting inserted after the BasePCDSystem
             * constructor params, so number 9!) */
            this.addParam(" U ", 1e-3, " Range ", [1e-5, 1e-2], " Desired ", 50);
        }


        function  plot(models.BaseFullModel model,double t,matrix<double> y) {

            if length(t) > 700
                idx = round(linspace(1,length(t),700));
                t = t(idx);
                y = y(:,idx);
            end
            states = [" alive "," unstable "," dead "];

            figure;
            X = t;
            Y = this.Omega(1):this.h:this.Omega(2);
            m = model.System.Dims(1);
            doplot(y(1:m,:)," Caspase-8 (x_a) ",1);
            doplot(y(m+1:2*m,:)," Caspase-3 (y_a) ",2);
            doplot(y(2*m+1:3*m,:)," Pro-Caspase-8 (x_i) ",3);
            doplot(y(3*m+1:end,:)," Pro-Caspase-3 (y_i) ",4);

            function doplot(y,thetitle,pnr)
                subplot(2,2,pnr);
                mesh(X,Y,y);
                /* surf(X,Y,y,'EdgeColor','none'); */
                xlabel(" Time [s] ");
                ylabel(sprintf(" %.2e to %.2e: Cell core to hull [m] ",this.Omega(1),this.Omega(2)));
                grid off;
                mi = min(y(:));
                Ma = max(y(:));
                if abs((mi-Ma) / mi) < 1e-14
                    mi = .999*mi; Ma=1.001*Ma;
                end
                axis([0 max(t) this.Omega mi Ma]);
                di = abs(this.Model.SteadyStates(:,pnr)-y(end));
                reldi = di ./ (this.Model.SteadyStates(:,pnr)+eps);
                reldistr = Utils.implode(reldi," ,  "," %2.3e ");
                if any(reldi > .1) || any(reldi < 10)
                    [~, id] = min(di);
                    title(sprintf(" Model '%s', %s concentrations\nCell state at T=%d: %s\n%s ", model.Name, thetitle,...
                    max(t),states[id],reldistr));
                else
                    title(sprintf(" Model '%s', %s concentrations\n%s ", model.Name, thetitle,reldistr));
                end
            end
        }
    protected:

        function  newSysDimension() {
            m = prod(this.Dims);

            x0 = zeros(4*m,1);
            x0(1:2*m) = 1e-16;
            x0(2*m+1:end) = 1e-9;
            this.x0= dscomponents.ConstInitialValue(x0);

            e = ones(m,1);
            A = spdiags([e -2*e e],-1:1,m,m)/this.hs^2;
            A(1) = -1/this.hs^2;
            A(end) = -1/this.hs^2;
            A = blkdiag(A,this.Diff(1)*A,this.Diff(2)*A,this.Diff(3)*A);
            this.A= dscomponents.LinearCoreFun(A);

            /*  Extracts the caspase-3 concentrations from the result
 *             C = zeros(m,4*m);
 *             C(:,m+1:2*m) = diag(ones(m,1));
             *this.C = dscomponents.LinearInputConv(sparse(C)); */
        }
};
}
}