Defining a function (including vector dot product) for all the points in 3D

I am trying to build the following function in a three dimensional domain.

where k is a constant vector, X is the position vector, c is a constant number, and t is time.

k is a vector of size [1 3], X is an array of size [NX*NY*NZ 3] that represents the points in the three-dimensional domain, c is a constant, and t is an array of size [1 NT].

The following is the setup of the problem.

    dx = 0.1;
    dy = 0.5;
    dz = 0.1;

    [x, y, z] = meshgrid( (1:100)*dx, (1:100)*dy, (1:100)*dz );
    X = [x(:) y(:) z(:)];

    k = [1 2 3];
    c = 0.5;    
    t = 0:0.1:1;

I thought about using arrayfun and repeating the vector k using repmat and dot it with X in the second dimension but I don't know what I should do for the multiplication of c and t.

In fact, the following loop works but it is very slow (takes 200 seconds on my machine).

    f = zeros(numel(X)/3, numel(t));
    for n = 1:numel(t)
        for i = 1:numel(X)/3
            f(i, n) = tan(dot(k, X(i,:)+c*t(n)));
        end
    end

What would be an efficient way of defining the function for all the points and all the times? The output of this function, for example, looks like an array of size [NX*NY*NZ NT].

[Udit Gupta]

Can you give an example with a few points in the X array?

Answer 1

[Matt J]

 f = tan(  bsxfun(@plus, X*k(:), c*t)  );

[AP]

I would be very thankful if you could let me know how this statement works so that I learn it for future use.

[Matt J]

X*k(:) computes all the dot(X(i,:),k)'s

For column vector u and row vector v,

 Y=bsxfun(@plus, u,v)

computes a matrix Y with Y(i,j)=u(i)+v(j). See "doc bsxfun".

[Matt J]

This is probably a bit faster

    k = [1 2 3];
    c = 0.5;    
    t = 0:0.1:1;

    [x, y, z] = meshgrid( (1:100)*(k(1)*dx), 
                          (1:100)*(k(2)*dy), 
                          (1:100)*(k(3)*dz);

   f=tan( bsxfun(@plus, x(:)+y(:)+z(:), ct) );

Answer 2

[George Papazafeiropoulos]

% data
dx = 0.1;
dy = 0.5;
dz = 0.1;
[x, y, z] = meshgrid( (1:100)*dx, (1:100)*dy, (1:100)*dz );
X = [x(:) y(:) z(:)];
k = [1 2 3];
c = 0.5;
t = 0:0.1:1;
lt=length(t);
% engine
u=numel(X)/3;
t=t(ones(u,1),:);
X=repmat(X,lt,1);
t=t(:);
t=t(:,ones(1,3));
u1=sum(k(ones(numel(X)/3,1),:).*(X+c*t),2);
ff=tan(u1);
% result
ff=reshape(ff,u,[])