FFTW , histos and matrix

Hi everyone,

I’m currently trying to build a simple 2D Electrostatic field solver using Fourrier transform .

I begin with a TMatrixF That I fill with my initial Potential (V) guess .

From there, I want to (let’s call initial guess Vini) :

[quote]1) Vin(k)i=(1.0/sqrt(N))*FFT(Vini(x)); // I transform to the f-space
2) Epix = F(k)*V(k), Egap=F2(k)*V(k) // I use some formula to get the E field in differents domain from Vini(k);
3) Then , Epix(x)=(1.0/sqrt(N))*FFT(Epix(x));

4)Then I mask Epix(x) over Egap(x)'s domain, and vice-versa
Epix(x)=Maskgap(Epix); and Egap=maskPixel(Egap);

  1. I combine both to get ETotal=Epix + Egap;

  2. I get back a potential using a defined function

  3. I set some boundary conditions on the potential (example potential over the pixel =1000V)

And I repeat 1 to 7 until it converge.

Now my problem is that it seems FFTW works only for histograms ?
What if I want to do it with a TMatrix ?

Now, I try to fill an histo with my matrix, then do my operation on histos like if they were matrix. That is a pain in the ass and I’m stuck in a hell of error of the type :

Error: Can’t call TH1::operator in current scope solver.c:80:
Possible candidates are…
Error: Binary operator oprand missing solver.c:80:
Error: Binary operator oprand missing solver.c:80:
Error: Binary operator oprand missing solver.c:80:
Error: operator ‘/’ divided by zero solver.c:80:
Error: Illegal pointer operation (tovalue) solver.c:80:
Error: improper lvalue solver.c:80:[/quote]

Is there any elegant and symbolic way to do this without too much pain ?

Thank you for your sugestion. I can post the code I have for now if it can help.

Mathieu B.

Hi Mathieu,

Not sure what you exactly want / miss . but glancing at tutorials/FFT.C,
I see that one can use the FFT functionality with histograms and
just arrays .

The data transfer between array’s and matrix is easily done through
functionality like TMatrixD/F::GetMatrixArray() TMatrixD/F::Use…


Hello, everyone!
I want to operate with matrix M_ij as a function.
For example, I want to fill matrix with f(x) (x is not initialized), then to find trace of it and draw a function Tr(M_ij) as a function of x.
Can ROOT do it?

Thanks in advance.

SMatrix has the Trace method:

root.cern.ch/doc/master/classRO … 38d51b34fb

Dear rooters,

I would like to share here a root macro dealing with FFT’s in 2D geometry.
The physical problem is just a classic electrostatic problem, i.e.:
Calculating the electric fields for a given charge density distribution at rest.
The basic equation is Gauss law (gradient of the electric field equal to the charge density):
grad(E) = - grad(grad phi) = rho
Here phi is the electrostatic potential such that E = -grad(phi).

Making a Fourier transform in both sides of the Gauss equation leads to: -(kz^2 + kx^2) F(phi) = F(rho).
Where kz and kx are the corresponding frequencies in z and x direction respectively.

The root macro deals with a 2D Gaussian charge distribution with (\sigma_z = 1 and \sigma_x = 1),
but this can be easily change to whatever charge distribution.
The charge distribution is normalised such that the peak density at the origin is equal to -1 (rho(0,0) = -1).

The algorithm goes as follows:

  1. Fourier transform to rho to obtain F(rho) (kz,kx)

  2. Calculate the fourier transform of phi by means of the Fourier transformed Gauss equation:
    F(phi) = - F(rho) / (kz^2 + kx^2)

  3. Fourier transforming back to phi(z,x): phi = F^-1(F(phi))

  4. Take the gradient of phi for the electric fields.

The macro generates 2D histograms for all the quantities described here.
They can be easily drawn in a root interactive session.

here the macro in question:
Gauss2D.C (5.84 KB)
The result seems OK for me,
There is a simple check against the analytical solution by the end of the macro (commented out).