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));
Egap(x)=(1.0/sqrt(N))
4)Then I mask Epix(x) over Egap(x)'s domain, and vice-versa
Epix(x)=Maskgap(Epix); and Egap=maskPixel(Egap);
I combine both to get ETotal=Epix + Egap;
I get back a potential using a defined function
V(x)=FFT(F(k)*FFT(Etotal(x)));
I set some boundary conditions on the potential (example potential over the pixel =1000V)
V(x)=boundarycondition(V(x));[/quote]
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 :
[quote]
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.
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?
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:
Fourier transform to rho to obtain F(rho) (kz,kx)
Calculate the fourier transform of phi by means of the Fourier transformed Gauss equation:
F(phi) = - F(rho) / (kz^2 + kx^2)
Fourier transforming back to phi(z,x): phi = F^-1(F(phi))
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).
Cheers!