1D interpolation of 3D data where ndgrid is not possible to use

Hello all,

I have a program which simulates fields on a surface as a function of time. For example, there may be a sphere and an electric field at points on the sphere as a function of time. It results in 3D data with dimensions of Etot = [ detectorpoint X spatial direction X timepoints]. The main issue is that the time points are not on a world time line. For example, the times where I have data for detector point 1, Etot(1,1,:), and detector point 2, Etot(2,1,:), are not the same although I do know the times for each respective data point. I would like to interpolate along the 3rd dimension but due to the irregularity of the time points for each detector point ndgrid is not applicable.

The reason for looking into this method of using interpn instead of looping through Etot is because loops can be slow and I am expecting large arrays. I have a code that will run until it gets stuck on interpn. It works fine if you set the variable xyznum = 2 but anything above that is won't. Also you will note that if you run it for xyznum =2 the second way of interpolating is much faster. I think that explains my problem.

%%Constants
cvel = 1;
k = 1;
tau = 2*pi; 
Enot = 1;
ntcut = 1/2; 
%%Detector Setup

R = 1000; 
%!!! number of detector points
xyznum =3;
theta = linspace(0,pi,xyznum);
phi = linspace(0,2*pi,xyznum+1);
xrad = zeros(3,xyznum^2);
dtheta = abs(theta(2)-theta(1));
dphi = abs(phi(2)-phi(1));
m = 1;
for i = 1:xyznum
  for j = 1:xyznum
      xrad(:,m) = R*[sin(theta(i))*cos(phi(j));...
          sin(theta(i))*sin(phi(j));...
          cos(theta(i))];
      m = m+1;
  end
end
Detectpts = xyznum^2;
Rmin = R;
Rmax = R;
%!!! number of time points in world time
tpts = 2^nextpow2(10);
tworld_start = (Rmin)-2*ntcut*cvel*k*tau;
tworld_end = (Rmax)+2*ntcut*cvel*k*tau;
tworld = linspace(tworld_start,tworld_end,tpts);

%%Traj setup

 %!! number of time points in "experiment"
nt = 5;
t = linspace(-ntcut*tau,ntcut*tau,nt);
r = Enot*(cos(t)-1);
vx = Enot*sin(t);
ax = -Enot*cos(t);

rtot = [r;zeros(2,nt)];
Vtot = [vx;zeros(2,nt)];
Atot = [ax;zeros(2,nt)];

%%Field set up

Eworld = zeros(Detectpts,3,tpts);
Bworld = zeros(Detectpts,3,tpts);

cvel = 1;

% remove duplicates of xrad and find uniques indices
[~,iunique,~] = unique(xrad','rows');

%%Vectorized formula
% This segment calculates Etot and Btot.

xrad = repmat(xrad',[1,1,nt]);
rtot = permute(repmat(rtot,[1,1,Detectpts]),[3,1,2]);
beta = permute(repmat(Vtot,[1,1,Detectpts]),[3,1,2]);
betadot = permute(repmat(Atot,[1,1,Detectpts]),[3,1,2]);

R = xrad - rtot;
Rmag = repmat((sum(R.^2,2)).^(0.5),[1,3,1]);
n = R./Rmag;
OneMinusNB = ones(Detectpts,3,nt) - repmat(sum(n.*beta,2),[1,3,1]);
numer1 = repmat(sum(n.*betadot,2),[1,3,1]).*(n-beta);
numer2 = OneMinusNB.*betadot;
denom = Rmag.*OneMinusNB.*OneMinusNB.*OneMinusNB;
Etot = (numer1 - numer2)./denom; %dimensions [Detectpts X 3 X nt ] 
Btot = cross(n,Etot,2); %dimensions [Detectpts X 3 X nt ] 
% !!!! Here is were the time is calculated  which I will want to
% interpolate onto a world time
time = permute(repmat(t,[3,1,Detectpts]),[3,1,2]) + Rmag; % dimensions [Detectpts X 3 X nt ] 

%clear large variables
clear R Rmag n OneMinusNB numer1 numer2 denom

%%Interpolate Vectorized Formula

% Here is one method I have used but it takes rather long due to the
% slicing of Etot
EtotWorld = zeros(Detectpts,3,tpts);
BtotWorld = zeros(Detectpts,3,tpts);
tic
for i = 1:Detectpts
  for k = 1:3
      EE = Etot(i,k,:);
      BB = Btot(i,k,:);
      tt = time(i,k,:);
      A = interpn(tt,EE,tworld);
      EtotWorld(i,k,:) = A;
      A = interpn(tt,BB,tworld);
      BtotWorld(i,k,:) = A;
  end
end
toc

tic
x = 1:1:Detectpts;
X = permute(repmat(x,[3,1,nt]),[2,1,3]); % original
XX = permute(repmat(x,[3,1,tpts]),[2,1,3]); % interpolated
y = 1:1:3;
Y = repmat(y,[Detectpts,1,nt]); % original
YY = repmat(y,[Detectpts,1,tpts]); % interpolated
TWORLD = permute(repmat(tworld,[3,1,Detectpts]),[3,1,2]); % interpolated.
EInterpTest = interpn(X,Y,time,Etot,XX,YY,TWORLD);
toc

[KSSV]

YOu need to make X,Y,time,Etot into 4D matrices to make it work.

[Colton Fruhling]

Why?