diff --git a/+quantity/Discrete.m b/+quantity/Discrete.m
index 88c6d9b9190c0581da7544404c2d418f1391e43b..6b1085d4a6bf9f781032c2eb3d92b8a294198410 100644
--- a/+quantity/Discrete.m
+++ b/+quantity/Discrete.m
@@ -47,7 +47,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
% size(valueOriginal{it}) == gridSize
% OR
% 2) adouble-array with
- % size(valueOriginal) == [gridSize, size(quantity)]
+ % size(valueOriginal) == [gridSize, size(quantity)]
% Furthermore, 'gridName' must be part of the name-value-pairs
% in varargin. Additional parameters can be specified using
% name-value-pair-syntax in varargin.
@@ -155,10 +155,10 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
obj = reshape(obj, objSize);
end
end% Discrete() constructor
-
+
%---------------------------
% --- getter and setters ---
- %---------------------------
+ %---------------------------
function i = isConstant(obj)
% the quantity is interpreted as constant if it has no grid or
% it has a grid that is only one point.
@@ -190,7 +190,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
end
obj.gridName = name;
end
-
+
function set.grid(obj, grid)
if ~iscell(grid)
grid = {grid};
@@ -285,7 +285,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
if isempty(a)
if nargin > 1
varargin{1}.assertSameGrid(varargin{2:end});
- end
+ end
return;
else
referenceGridName = a(1).gridName;
@@ -328,9 +328,9 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
else
referenceGridName = a(1).gridName;
referenceGrid = a(1).grid;
- referenceGridSize = a(1).gridSize(referenceGridName);
+ referenceGridSize = a(1).gridSize(referenceGridName);
end
-
+
for it = 1 : numel(varargin)
if isempty(varargin{it})
continue;
@@ -407,7 +407,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
% X3 ...]' when any of X1, X2, X3, etc. is an object.
%
% See also HORZCAT, CAT.
- c = cat(2, a, varargin{:});
+ c = cat(2, a, varargin{:});
end
function c = vertcat(a, varargin)
%VERTCAT Vertical concatenation.
@@ -448,27 +448,27 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
% B = CAT(DIM,A1,A2,A3,A4,...) concatenates the input
% arrays A1, A2, etc. along the dimension DIM.
%
- % When used with comma separated list syntax, CAT(DIM,C{:}) or
+ % When used with comma separated list syntax, CAT(DIM,C{:}) or
% CAT(DIM,C.FIELD) is a convenient way to concatenate a cell or
% structure array containing numeric matrices into a single matrix.
%
% Examples:
- % a = magic(3); b = pascal(3);
+ % a = magic(3); b = pascal(3);
% c = cat(4,a,b)
% produces a 3-by-3-by-1-by-2 result and
% s = {a b};
- % for i=1:length(s),
+ % for i=1:length(s),
% siz{i} = size(s{i});
% end
% sizes = cat(1,siz{:})
% produces a 2-by-2 array of size vectors.
- %
+ %
% See also NUM2CELL.
-
+
% Copyright 1984-2005 The MathWorks, Inc.
% Built-in function.
if nargin == 1
- objCell = {a};
+ objCell = {a};
else
objCell = [{a}, varargin(:)'];
@@ -490,7 +490,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
c = quantity.Discrete.empty(S);
return
else
- obj = objCell{objIdx};
+ obj = objCell{objIdx};
end
for k = 1:numel(objCell(~isEmpty))
@@ -505,7 +505,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
if isempty(value)
M = zeros([prod(obj(1).gridSize), size(value(l))]);
end
- M = reshape(M, [obj(1).gridSize, size(value)]);
+ M = reshape(M, [obj(1).gridSize, size(value)]);
o = quantity.Discrete( M, ...
'size', size(value), ...
'gridName', obj(1).gridName, ...
@@ -520,7 +520,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
[fineGrid, fineGridName] = getFinestGrid(objCell{~isEmpty});
for it = 1 : (numel(varargin) + 1) % +1 because the first entry is a
objCell{it} = objCell{it}.changeGrid(fineGrid, fineGridName);
- end
+ end
assertSameGrid(objCell{:});
argin = [{dim}, objCell(:)'];
c = builtin('cat', argin{:});
@@ -569,13 +569,13 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
end
solution = objInverseTemp.on(rhs);
-
-% solution = zeros(numel(obj), 1);
-% for it = 1 : numel(obj)
-% objInverseTemp = obj(it).invert(gridName);
-% solution(it) = objInverseTemp.on(rhs(it));
-% end
-% solution = reshape(solution, size(obj));
+
+ % solution = zeros(numel(obj), 1);
+ % for it = 1 : numel(obj)
+ % objInverseTemp = obj(it).invert(gridName);
+ % solution(it) = objInverseTemp.on(rhs(it));
+ % end
+ % solution = reshape(solution, size(obj));
end
function inverse = invert(obj, gridName)
@@ -611,7 +611,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
myParser.parse(varargin{:});
variableGrid = myParser.Results.variableGrid;
- myGridSize = [numel(variableGrid), ...
+ myGridSize = [numel(variableGrid), ...
numel(myParser.Results.initialValueGrid)];
% the time (s) vector has to start at 0, to ensure the IC.
@@ -632,14 +632,14 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
if numel(positiveVariableGrid) > 1
[resultGridPositive, odeSolutionPositive] = ...
ode45(@(y, z) obj(it).on(z), ...
- positiveVariableGrid, ...
- myParser.Results.initialValueGrid(icIdx), options);
+ positiveVariableGrid, ...
+ myParser.Results.initialValueGrid(icIdx), options);
end
if numel(negativeVariableGrid) >1
[resultGridNegative, odeSolutionNegative] = ...
ode45(@(y, z) obj(it).on(z), ...
- negativeVariableGrid, ...
- myParser.Results.initialValueGrid(icIdx), options);
+ negativeVariableGrid, ...
+ myParser.Results.initialValueGrid(icIdx), options);
end
if any(variableGrid == 0)
resultGrid = [flip(resultGridNegative(2:end)); 0 ; resultGridPositive(2:end)];
@@ -682,8 +682,8 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
assert(numel(values) == numel(gridName2Replace), ...
'gridName2Replace and values must be of same size');
- % here substitution starts:
- % The first (gridName2Replace{1}, values{1})-pair is
+ % here substitution starts:
+ % The first (gridName2Replace{1}, values{1})-pair is
% replaced. If there are more cell-elements in those inputs
% then subs() is called again for the remaining pairs
% (gridName2Replace{2:end}, values{2:end}).
@@ -725,10 +725,10 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
newValue = misc.diagNd(obj.on(newGridForOn), gridIndices);
newGrid = {newGridForOn{gridIndices(1)}, ...
newGridForOn{1:1:numel(newGridForOn) ~= gridIndices(1) ...
- & 1:1:numel(newGridForOn) ~= gridIndices(2)}};
+ & 1:1:numel(newGridForOn) ~= gridIndices(2)}};
newGridName = {values{1}, ...
obj(1).gridName{1:1:numel(obj(1).gridName) ~= gridIndices(1) ...
- & 1:1:numel(obj(1).gridName) ~= gridIndices(2)}};
+ & 1:1:numel(obj(1).gridName) ~= gridIndices(2)}};
else
% this is the default case. just grid name is
@@ -790,7 +790,12 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
% sich ungewöhnlich verhält, wenn man sie ohne idx argument
% aufruft.
if nargin == 1
- value = 1:obj.gridSize;
+
+ if numel(obj.gridSize) == 1
+ value = zeros(obj.gridSize, 1);
+ else
+ value = zeros(obj.gridSize, 1);
+ end
if isempty(value)
value = 0;
end
@@ -798,12 +803,15 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
if ~iscell(varargin)
varargin = {varargin};
end
- value = cellfun(@(v) v(varargin{:}), {obj.valueDiscrete});
- value = reshape(value, size(obj));
+ value = cellfun(@(v) v(varargin{:}), {obj.valueDiscrete}, ...
+ 'UniformOutpu', false);
+
+ valueSize = size(value{1});
+
+ value = reshape([value{:}], [valueSize(1:obj(1).nargin), size(obj)]);
end
end
-
function n = nargin(obj)
% FIXME: check if all funtions in this object have the same
% number of input values.
@@ -910,7 +918,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
end
end
-
+
function myLabel = labelHelper(gridNumber)
myLabel = ['$$', greek2tex(obj(rowIdx, columnIdx, figureIdx).gridName{gridNumber}), '$$'];
end
@@ -999,6 +1007,26 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
[newObj.grid] = deal(myGrid);
end
+ function [lowerBounds, upperBounds] = getBoundaryValues(obj, varargin)
+ % GETBOUNDARYVALUES returns the boundary values of the grid
+ % [lowerBounds, upperBounds] = getBoundaryValues(obj, gridName)
+ % returns the boundary values of the grid, specified by
+ % gridName
+ %
+ % [lowerBounds, upperBounds] = getBoundaryValues(obj) returns
+ % the boudary values of all grids defined for the object.
+
+ grids = obj(1).grid(obj.gridIndex(varargin{:}));
+ lowerBounds = zeros(1, numel(grids));
+ upperBounds = zeros(1, numel(grids));
+
+ for k = 1:numel(grids)
+ lowerBounds(k) = grids{k}(1);
+ upperBounds(k) = grids{k}(end);
+ end
+
+ end
+
end
%% math
@@ -1103,7 +1131,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
% misc.multArray is very efficient, but requires that the
% multiple dimensions of the input are correctly arranged.
[idx, permuteGrid] = computePermutationVectors(a, b);
-
+
parameters = struct(a(1));
parameters.name = [a(1).name, ' ', b(1).name];
parameters.grid = [...
@@ -1130,7 +1158,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
if b.gridIndex(parameters.gridName{it})
bGrid{b.gridIndex(parameters.gridName{it})} = ...
parameters.grid{it};
- end
+ end
end
parameters = misc.struct2namevaluepair(parameters);
@@ -1164,11 +1192,11 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
objDiscreteReshaped = permute(reshape(objDiscreteOriginal, ...
[prod(obj(1).gridSize), size(obj)]), [2, 3, 1]);
invDiscrete = zeros([prod(obj(1).gridSize), size(obj)]);
-
+
parfor it = 1 : size(invDiscrete, 1)
invDiscrete(it, :, :) = inv(objDiscreteReshaped(:, :, it));
end
-
+
y = quantity.Discrete(reshape(invDiscrete, size(objDiscreteOriginal)),...
'size', size(obj), ...
'name', ['(', obj(1).name, ')^{-1}'], ...
@@ -1208,11 +1236,11 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
function x = mldivide(A, B)
% mldivide see doc mldivide of matlab:
- % "x = A\B is the solution to the equation Ax = B. Matrices
+ % "x = A\B is the solution to the equation Ax = B. Matrices
% A and B must have the same number of rows."
x = inv(A)*B;
end
-
+
function x = mrdivide(B, A)
% mRdivide see doc mrdivide of matlab:
% "x = B/A solves the system of linear equations x*A = B for x.
@@ -1319,6 +1347,10 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
% gridName properties of obj for the "names" and returns its
% index as "idx" and its logical index as "log"
+ if nargin == 1
+ names = obj(1).gridName();
+ end
+
if ~iscell(names)
names = {names};
end
@@ -1349,9 +1381,9 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
if obj.isConstant && isempty(obj(1).gridName)
result = quantity.Discrete(zeros(size(obj)), ...
- 'size', size(obj), 'grid', obj(1).grid, ...
- 'gridName', obj(1).gridName, ...
- 'name', ['(d_{.}', obj(1).name, ')']);
+ 'size', size(obj), 'grid', obj(1).grid, ...
+ 'gridName', obj(1).gridName, ...
+ 'name', ['(d_{.}', obj(1).name, ')']);
return
end
@@ -1393,45 +1425,72 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
'name', ['(d_{', gridName, '}', obj(1).name, ')']);
end
- function I = int(obj, intGridName)
+ function I = int(obj, varargin)
% INT integrate
% INT(obj) is the definite integral of obj with respect to
% all of its independent variables over the grid obj(1).grid.
- % INT(obj, intGridName) is the definite integral of obj with
- % respect to the spatial coordinate intGridName, that has to
+ % INT(obj, gridName) is the definite integral of obj with
+ % respect to the spatial coordinate grid name, that has to
% be one of obj(1).gridName. The names must be a string or a
% array of strings.
+ % INT(obj, a, b)
+
+ intGridName = obj(1).gridName;
+ [a, b] = obj.getBoundaryValues;
if nargin == 1
- intGridName = obj(1).gridName;
+ % see default settings.
+ elseif nargin == 2 || nargin == 4
+ % (obj, z) OR (obj, z, a, b)
+ intGridName = varargin{1};
+ elseif nargin == 3
+ % (obj, a, b)
+ a = varargin{1};
+ b = varargin{2};
end
- if nargin > 2
- error('Not yet implemented');
+ if numel(intGridName) > 1
+ obj = obj.int(intGridName{1}, a(1), b(1));
+ I = obj.int(intGridName{2:end}, a(2:end), b(2:end));
+ return
end
+
+ %% do the integration over one dimension
I = obj.copy();
- [idxName, isGrid] = obj.gridIndex(intGridName); % get index of the variables that should be considered for integration
+ % get index of the variables that should be considered for
+ % integration
+ [idxGrid, isGrid] = obj.gridIndex(intGridName);
+ [domainIdx{1:obj(1).nargin}] = deal(':');
+ currentGrid = obj(1).grid{ idxGrid };
+ domainIdx{idxGrid} = find(currentGrid >= a(idxGrid) & currentGrid <= b(idxGrid));
+
+ myGrid = currentGrid(domainIdx{idxGrid});
+
for kObj = 1:numel(obj)
- J = obj(kObj).on(); % get numerical values
- for kName = 1:length(idxName) % iterate over all the integration dimensions
- currentGrid = obj(kObj).grid{idxName(kName)};% currentGrid that should be used for the integration.
- J = numeric.trapz_fast_nDim(currentGrid, J, idxName(kName));
+ J = numeric.trapz_fast_nDim(myGrid, obj(kObj).atIndex(domainIdx{:}), idxGrid);
+
+ % If the quantity is integrated, so that the first
+ % dimension collapses, this dimension has to be shifted
+ % away, so that the valueDiscrete can be set correctly.
+ if size(J,1) == 1
+ J = shiftdim(J);
end
if all(isGrid)
grdName = '';
- grd = {};
+ newGrid = {};
else
grdName = obj(kObj).gridName{~isGrid};
- grd = obj(kObj).grid{~isGrid};
+ newGrid = obj(kObj).grid{~isGrid};
end
-
+
% create result:
I(kObj).valueDiscrete = J;
I(kObj).gridName = grdName;
- I(kObj).grid = grd;
+ I(kObj).grid = newGrid;
end
+
end
function result = cumInt(obj, varargin)
@@ -1442,7 +1501,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
% result = grid{1}^boundaryGridName ...
% f(..., zeta, ... ) dintGrid
- % two input parser since some default values depend on
+ % two input parser since some default values depend on
% intGridName.
myParser1 = misc.Parser;
myParser1.addParameter('intGridName', obj(1).gridName{1});
@@ -1464,7 +1523,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
result = quantity.Discrete(F, 'size', size(obj), 'grid', myGrid, ...
'gridName', obj(1).gridName, 'name', ['int(', obj(1).name, ')']);
- % substitute gridName from original gridName (= intGridName) to
+ % substitute gridName from original gridName (= intGridName) to
% boundaryGridName
if ~iscell(boundaryGridName)
boundaryGridName = {boundaryGridName};
@@ -1557,7 +1616,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
function maxValue = max(obj)
% max returns the maximal value of all elements of A over all
% variables as a double array.
- maxValue = reshape(max(obj.on(), [], 1:obj(1).nargin), [size(obj)]);
+ maxValue = reshape(max(obj.on(), [], 1:obj(1).nargin), [size(obj)]);
end
function maxValue = MAX(obj)
@@ -1568,20 +1627,20 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
function minValue = min(obj)
% max returns the minimum value of all elements of A over all
% variables as a double array.
- minValue = reshape(min(obj.on(), [], 1:obj(1).nargin), [size(obj)]);
+ minValue = reshape(min(obj.on(), [], 1:obj(1).nargin), [size(obj)]);
end
function absQuantity = abs(obj)
% abs returns the absolut value of the quantity as a quantity
absQuantity = quantity.Discrete(abs(obj.on()), ...
'grid', obj(1).grid, 'gridName', obj(1).gridName, ...
- 'size', size(obj), 'name', ['|', obj(1).name, '|']);
+ 'size', size(obj), 'name', ['|', obj(1).name, '|']);
end
function meanValue = mean(obj)
% max returns the mean value of all elements of A over all
% variables as a double array.
- meanValue = reshape(mean(obj.on(), 1:obj(1).nargin), [size(obj)]);
+ meanValue = reshape(mean(obj.on(), 1:obj(1).nargin), [size(obj)]);
end
function f = volterra(K, x, varargin)
@@ -1610,13 +1669,13 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
ip.parse2ws(varargin{:});
if ip.isDefault('variable') && x.nargin > 1
- error(['For x quantities wrt. two independent' ...
+ error(['For x quantities wrt. two independent' ...
' variables, the name of the variable has to be' ...
'parametrized explicitly']);
end
- %% prepare numerical data
+ %% prepare numerical data
% do the multiplication under the integration and safe the
% information for the grid:
v = K*x;
@@ -1626,7 +1685,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
t.log = ~s.log;
t.grid = v(1).grid{t.log};
t.name = v(1).gridName{t.log};
- V = v.on();
+ V = v.on();
% check if s is the second independent variable in V
if s.idx == 1 % permutation needed in this case!
@@ -1659,17 +1718,17 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
- numeric.trapz_fast_nDim(s.grid(1:tIdx), V_t(1:tIdx, :, :)));
end
end
-
+
if fixedLowerBound
f = quantity.Discrete(F, 'size', size(v), 'grid', {t.grid}, ...
- 'gridName', {t.name}, 'name', ['volterra(', K(1).name, x(1).name, ')']);
+ 'gridName', {t.name}, 'name', ['volterra(', K(1).name, x(1).name, ')']);
else
f = quantity.Discrete(F, 'size', size(v), 'grid', {t.grid, s.grid}, ...
'gridName', {t.name, s.name}, 'name', ['volterra(', K(1).name, x(1).name, ')']);
end
end
-
+
function q = obj2value(obj)
% for a case with a 3-dimensional grid, the common cell2mat
% with reshape did not work. So it has to be done manually:
@@ -1812,7 +1871,7 @@ classdef (InferiorClasses = {?quantity.Symbolic, ?quantity.Operator}) Discrete
function [idx, permuteGrid] = computePermutationVectors(a, b)
% Computes the required permutation vectors to use
% misc.multArray for multiplication of matrices
-
+
uA = unique(a(1).gridName);
assert(numel(uA) == numel(a(1).gridName), 'Gridnames have to be unique!');
uB = unique(b(1).gridName);
diff --git a/+quantity/Symbolic.m b/+quantity/Symbolic.m
index 33d381a48155a91d48e3ee9b623dafb0a84ec58f..f0e3d75485b3bfa81a2b8675a3c561eadaf97ca6 100644
--- a/+quantity/Symbolic.m
+++ b/+quantity/Symbolic.m
@@ -493,10 +493,15 @@ classdef Symbolic < quantity.Function
'name', ['sqrtm(', x(1).name, ')']);
end
- function D = diff(obj, k)
+ function D = diff(obj, k, gridName)
if nargin == 1
k = 1;
end
+
+ if nargin == 3
+ error('Not yet implemented')
+ end
+
% TODO: specify for which variable it should be differentiated
D = obj.copy();
[D.name] = deal(['(d_{' char(obj(1).variable) '} ' D(1).name, ')']);
@@ -656,45 +661,12 @@ classdef Symbolic < quantity.Function
'variable', obj(1).variable);
end
-% function result = cumInt(obj, varargin)
-% % cumInt performes the integral
-% % result = int_boundaryGridName(1)^boundaryGridName(2) ...
-% % f(..., zeta, ... ) dintGrid
-% % or if numel(boundaryGridName) == 1
-% % result = grid{1}^boundaryGridName ...
-% % f(..., zeta, ... ) dintGrid
-%
-% % two input parser since some default values depend on
-% % intGridName.
-% myParser1 = misc.Parser;
-% myParser1.addParameter('intGridName', obj(1).gridName{1});
-% myParser1.parse(varargin{:});
-% intGridName = myParser1.Results.intGridName;
-% intGridIdx = obj.gridIndex(intGridName);
-% intVariable = obj(1).variable(intGridIdx);
-%
-% myParser2 = misc.Parser;
-% myParser2.addParameter('boundaryGridName', {intGridName});
-% myParser2.parse(varargin{:});
-% boundaryGridName = myParser2.Results.boundaryGridName;
-% if ~iscell(boundaryGridName)
-% boundaryGridName = {0, boundaryGridName};
-% end
-% if numel(boundaryGridName) == 1
-% boundaryGridName = {0, boundaryGridName{1}};
-% end
-% for it = 1 : numel(boundaryGridName)
-% if ischar(boundaryGridName{it})
-% boundaryValue{it} = obj(1).variable(obj.gridIndex(boundaryGridName{it}));
-% else
-% boundaryValue{it} = boundaryGridName{it};
-% end
-% end
-%
-% % integrate
-% resultSymbolic = int(obj.sym(), intVariable, boundaryValue{1}, boundaryValue{it});
-% result = quantity.Symbolic(resultSymbolic, 'name', ['int(', obj(1).name, ')'], ...
-% 'grid', obj(1).grid, 'variabe', obj(1).variable);
+
+
+ %% TODO das verhalten von int und cumint an das von quantity.Discrete anpassen!
+% function C = cumint(obj, z, a, b)
+% z = obj(1).grid{1};
+% C = obj.int(z, z(1), obj.variable);
% end
function C = int(obj, z, a, b)
diff --git a/+unittests/+quantity/testDiscrete.m b/+unittests/+quantity/testDiscrete.m
index 4dd258591393fbefc4f9285cd5f1b9552dfef59f..925b8f57e0e5c348a266169ebb7be28069069473 100644
--- a/+unittests/+quantity/testDiscrete.m
+++ b/+unittests/+quantity/testDiscrete.m
@@ -582,6 +582,20 @@ testCase.verifyEqual(a(end), A.atIndex(end));
testCase.verifyEqual(a(1), A.atIndex(1));
testCase.verifyEqual(a(23), A.atIndex(23));
+y = linspace(0,2,51);
+b1 = sin(z * pi * y);
+b2 = cos(z * pi * y);
+B = quantity.Discrete({b1; b2}, 'grid', {z, y}, 'gridName', {'z', 'y'});
+
+B_1_y = B.atIndex(1,1);
+b_1_y = [sin(0); cos(0)];
+testCase.verifyTrue(numeric.near(squeeze(B_1_y(1,1,:)), b_1_y));
+
+B_z_1 = B.atIndex(':',1);
+
+testCase.verifyTrue(all(B_z_1(:,1,1) == 0));
+testCase.verifyTrue(all(B_z_1(:,:,2) == 1));
+
end
function testGridJoin(testCase)