I’m studying a scalar field u(x,t) on a 2D box LxL, with periodic boundary conditions.
What I would like to do is to calculate the cross-correlation between u(x,t1) and u(x+r,t1) at varying of r, where r is the euclidean distance.

After discretized the domain (then now it could be see as a lattice), my problem is about how to reach, for each node of the lattice, the points on the circle of radius r.

I have implemented a simple way to do it, but it doesn’t work with the periodic boundary conditions:

r=4 #for instance
dd=1:N
kk=[Int64[] for i in dd]
for i in dd
for j in dd
if norm(A[:,i]-A[:,j], 2)==r
append!(kk[i],[j])
end
end
end

Where A is the 2xN matrix containing the coordinates of each node of the lattice. Here the kk vector contains for each node his “neighbors” at the euclidean distance r, but it not consider the bc.

I’m sure that exist a better way to solve this probem, every help is welcome, thanks!

Ah, I think I misunderstood. You’re only looking for 1-dimensional distance? If so, disregard my previous comment. But CircularArrays.jl might still be useful. Or just use the mod1() function (which is part of base julia).

I only want to know what points are at distance r from a fixed point (x1,y1), and my troubles are at the boundary of my domain, here the norm can’t see all the points at the correct distance.
Now I will take a look at CircularArrays.jl, but how do you think it could be used?

@Frazze do you need to preserve the flat coordinates with spherical boundary conditions or you are after an actual sphere geometry? Notice that you can perform the variogram calculations linked above on spheres just fine. Try to define your domain as domain = triangulate(Sphere((0.,0.), 1.0)) in GeoStats.jl and see if that addresses the problem.

We are still polishing the features of the project related to manifolds and non-Euclidean geometry, but maybe you can cook something up already. Try to read the GeoStats.jl and Meshes.jl docs to see if you find something useful.

The coordinates are only a “trick” for me to try to solve the problem. I have a 2D domain where some ODEs are calculated. Then I want to calculate the cross correlation of the solution between all the points of the box, but I wish to maintain the regular square-lattice discretization.

I have decided to implement this approach because the issue seems to pop up and the details (esp normalization and scaling) tend to be quite tricky. You can find the prototype here

I’m sorry @josuagrw, I take this post to do another question, still related to the topic.

What if I would like to calculate the power spectra (then the ft of spatial correlation) of mi scalar field at one fixed time t just having the solution of the field u(x,y,t)?

There is a way to do it without calculate before the correlation?

Thank you so much, I’m really new with this stuff…

First off: With periodic boundary conditions, the power spectrum is discrete because only specific wavelengths fulfill that boundary condition.

Once you realize that, it actually makes everything easier: The spectral density (power spectrum) is just the square absolute value of the fourier transform, so

there are different ways to do it. The naive way is

using Plots
scatter(
reshape(norm.(ks), :),
reshape(power_spectrum, :)
)

which typically looks quite messy (depending on how isotropic your data really is). If you want / need to clean that up more, it really depends on what kind of data you are analysing and how you want to present it.

Hi @Frazze sorry I forgot to tell you about fftshift(): Basically the modes are not ordered sequentially in memory, so plot(power_spectrum) will give you a wrong picture.

This is why it is important to use fftshiftorfftfreq (or their equivalentes rfftshift and rfftfreq)

Ok thank you,
Just another thing that it’s not so clear for me: What about the array dimension of power_spectrum?
I can’t understand why, starting from a data array of dim= 16384 (a domain 128x128 reshaped) with rfft I reach a vector of dim=8193, in particular this vector dim is different than the k_norm array dim, then I don’t know how to represent it…

Basically rfft halves the memory because the Fourier transform of a real input has a certain symmetry. Look here to see how you can calculate the modes in the 2D case