Taking the curl of a vector field in matlab using Fourier transforms

[Biftheunderstudy]

I've been given a data set which is stored in a 3D array and I need to do some analysis of it. One part of the analysis is to take the curl, and I've been told that this will be easier to do if I take a Fourier transform first. This is because the curl in k-space is just ikxA.

Now, I've been attempting to use the example that is provided in the Matlab help documents to ensure that I can do the curl properly in k-space...but it's not working.

Here is my Matlab script, if someone can figure out why it is not working. (By not working, I mean I can't produce the same plots as the example)

Also, keep in mind that I've done a few iterations to work out bugs so maybe I've changed something for the worse. For instance, whether to fftshift or not, and the range of the k-vectors.

clear all;

figure
load wind
[curlx,curly,curlz,cav] = curl(x,y,z,u,v,w);
h = slice(x,y,z,cav,[90 134],59,0); 
shading interp
daspect([1 1 1]); 
axis tight
colormap hot(16)
camlight
set([h(1),h(2)],'ambientstrength',.6)

%Use ffts to solve for the curl now,
u_ =fftn(u);
v_ = fftn(v);
w_ =fftn(w);

%Find the ranges of the data
min_x = min(min(min(x)));
min_y = min(min(min(y)));
min_z = min(min(min(z)));
max_x = max(max(max(x)));
max_y = max(max(max(y)));
max_z = max(max(max(z)));

%length of the indeces
nx=length(u_(1,:,1)); ny=length(v_(:,1,1)); nz=length(w_(1,1,:));

%array spacing
dx = (max_x-min_x)/(nx-1);
dy = (max_y-min_y)/(ny-1);
dz = (max_z-min_z)/(nz-1);


%Build the k-vectors
% ordering is:
% 0, n/(N/2*dx), ... , (N/4-1)/(N/2*dx), +/- 1/(2*dx),-(N/4-1)/(N/2*dx),..
% or:
% 0, 1, ...., k_max, -k_max,...,-1
for j=1:round(nx/2);
    kx1(j) = (j-1.)/nx/2./dx;
end
for j=round(nx/2)+1:nx 
    kx1(j) = (j-nx-1)/nx/2./dx;
end

for j=1:round(ny/2);
    ky1(j) = (j-1.)/ny/2./dy;
end
for j=round(nx/2)+1:ny 
    ky1(j) = (j-ny-1)/ny/2./dy;
end

for j=1:round(nz/2);
    kz1(j) = (j-1.)/nz/2./dz;
end
for j=round(nz/2)+1:nz 
    kz1(j) = (j-nz-1)/nz/2./dz;
end

% kx1 = -(max_x-min_x)/2.:dx:(max_x-min_x)/2.;
% ky1 = -(max_y-min_y)/2.:dy:(max_y-min_y)/2.;
% kz1 = -(max_z-min_z)/2.:dz:(max_z-min_z)/2.;

[kx,ky,kz]=meshgrid(kx1,ky1,kz1);
%Take the cross product of (ikx,iky,ikz)x(u,v,w)
%ax = by*cz - bz*cy
%ay = bz*cx - bx*cz
%az = bx*cy - by*cx
curlx = 1i.*ky.*w_ - 1i.*kz.*v_;
curly = 1i.*kz.*u_ - 1i.*kx.*w_;
curlz = 1i.*kx.*v_ - 1i.*ky.*u_;
%This part is the same as in the curl function
nrm = sqrt(u_.*conj(u_) + v_.*conj(v_) + w_.*conj(w_));
cav2 = .5 * (curlx.*u_ + curly.*v_ + curlz.*w_) ./nrm;

%Inverse transform
curlx2 = abs(ifftn(curlx));
curly2 = abs(ifftn(curly));
curlz2 = abs(ifftn(curlz));
cav2 = abs(ifftn(cav2));

%Plot the data, see if it matches the first figure
figure
h = slice(x,y,z,cav2,[90 134],59,0); 
shading interp
daspect([1 1 1]); 
axis tight
colormap hot(16)
camlight
set([h(1),h(2)],'ambientstrength',.6)

%The second example in the matlab documents
figure
k = 4;
x = x(:,:,k); y = y(:,:,k); u = u(:,:,k); v = v(:,:,k); 
cav = curl(x,y,u,v);
pcolor(x,y,cav); shading interp
hold on;
quiver(x,y,u,v,'y')
hold off
colormap copper

w = w(:,:,k);
u_ = fftn(u)/dx;
v_ = fftn(v)/dx;
w_ = fftn(w)/dx;

[kx,ky]=meshgrid(kx1,ky1);

kz = zeros(size(u_));

curlx = zeros(size(u_));
curly = zeros(size(u_));
curlz = zeros(size(u_));

%Do the cross product again
%ax = by*cz - bz*cy
%ay = bz*cx - bx*cz
%az = bx*cy - by*cx
curlx = 0.;
curly = 0.;
curlz = 1i*kx.*v_ - 1i*ky.*u_;

curlz2 = real(ifftn(curlz));

nrm = sqrt(u.^2 + v.^2);
cav2 = .5 * curlz2;

figure
pcolor(x,y,cav2); shading interp
hold on;
quiver(x,y,u,v,'y')
hold off
colormap copper

 

[William Gaatjes]

Cyclowizard has profound knowledge of matlab...
Maybe he can help you ?

 

[TecHNooB]

If I still didn't have grad school apps to finish I'd help you :|

 

[Wizlem]

I'm not too familiar with fft but is it possible the fft and ifft cause enough calculation error to throw off the numbers? Everything looks fine to me but I don't know how to go about calculating the k vectors.

 

[Biftheunderstudy]

I would expect that if the fft were to cause error like that, I would at least get a similar looking plot...which I kind of do for the second example. The first is WAY off though...

 

[TecHNooB]

I havent looked at your code, but if you're modifying the fft, make sure the sum of the imaginary parts of your ifft are almost zero (should be very very small). If not, you have imaginary components that you probably shouldn't have which means you did something wrong in processing the fft. Also, use real(ifft(x)) instead of abs(ifft(x)).

 

[Biftheunderstudy]

What do you mean by modifying the fft? Do you mean by taking the fft, doing some work on the data, then ifft?

All the transforms were done using Matlab's built in fft and ifft functions (which are just fftw).

Also, I've tried taking the real/imaginary/abs of the output and it's still wrong, just different wrong.

 

[TecHNooB]

What do you mean by modifying the fft? Do you mean by taking the fft, doing some work on the data, then ifft?

All the transforms were done using Matlab's built in fft and ifft functions (which are just fftw).

Also, I've tried taking the real/imaginary/abs of the output and it's still wrong, just different wrong.

real and abs should look similar if not the same. if that's not the case, your data from the ifft is not all real and you messed up somewhere when processing the fft data.