I can't figure out what is wrong with my C code for finding the second derivative of a function in 2d (Laplacian) using FFTW. I have taken a small bit of code from my main code that deals with finding the derivative below. Can anyone help please as I have tried everything I can think of, but nothing works? It is meant to be a an easy thing to do as the derivative of a function when you Fourier transform it turns into a multiplication in Fourier space.
/*
* File: main.c
* Author: stephen
*
* Created on 19 September 2022, 17:37
*/
#include <stdio.h>
#include <stdlib.h>
#include <tgmath.h>
#include <complex.h>
#include <fftw3.h>
/*
*
*/
int main(int argc, char** argv) {
int nx = 256;
int ny = 256;
fftw_complex *gauss = (fftw_complex *) fftw_malloc(sizeof (fftw_complex)
* nx * ny);
fftw_complex *Fgauss = (fftw_complex *) fftw_malloc(sizeof (fftw_complex)
* nx * ny);
fftw_complex *Fgaussderiv =
(fftw_complex *) fftw_malloc(sizeof (fftw_complex) * nx * ny);
fftw_complex *derivgauss =
(fftw_complex *) fftw_malloc(sizeof (fftw_complex) * nx * ny);
double *dderivgauss = (double *) fftw_malloc(sizeof ( double) * nx * ny);
fftw_complex *derivgaussact =
(fftw_complex *) fftw_malloc(sizeof (fftw_complex) * nx * ny);
fftw_plan plan_backward_gauss = fftw_plan_dft_2d(nx, ny, Fgaussderiv,
derivgauss, FFTW_BACKWARD, FFTW_EXHAUSTIVE);
fftw_plan plan_forward_gauss = fftw_plan_dft_2d(nx, ny, gauss,
Fgauss, FFTW_FORWARD, FFTW_EXHAUSTIVE);
double sig = 256.0;
double a = 1.0/sig*sqrt(2.0*M_PI);
double width = 20.0-E-3;
double height = 20.0-E-3;
double dx = width/(double)nx;
double dy = height/(double)ny;
for (int i = -nx/2; i <= nx/2 - 1; ++i) {
for (int j = -ny/2; j <= ny/2 - 1; ++j) {
//Gaussian in 2d
gauss[(nx/2 + i) + (ny/2 + j)*nx] =
CMPLX(a*exp(-((double)i * (double)i*dx*dx + (double)j *
(double)j*dy*dy)/2.0*sig*sig), 0.0);
//Laplacian of a Gaussian in 2d
derivgaussact[(nx/2 + i) + (ny/2 + j)*nx] = CMPLX((((-sig*sig +
(double)i*(double)i*dx*dx)+(-sig*sig +
(double)j*(double)j*dy*dy))/sig*sig*sig*sig)*a*exp(
-((double)i * (double)i*dx*dx + (double)j *
(double)j*dy*dy)/2.0*sig*sig), 0.0);
}
}
//Take the 2D Fourier transform of the input Gaussian
fftw_execute(plan_forward_gauss);
//Shift the Fourier transform so the low frequencies are at the centre
//Tested and works as expected.
fftshift2D(Fgauss, nx, ny);
double *KxKxKyKy = (double*) malloc(sizeof (double)*nx * ny);
for (int j = -ny / 2; j <= (ny / 2) - 1; j++) {
for (int i = -nx / 2; i <= (nx / 2) - 1; i++) {
//Multiply by the square of the wavenumber in Fourier space
KxKxKyKy[i + nx / 2 + ((ny / 2)+j)*nx] =
-((4.0*M_PI*M_PI/width*width)*(double)i*(double)i +
(4.0*M_PI*M_PI/height*height)*(double)j*(double)j);
}
}
for (int j = 0; j < ny; ++j) {
for (int i = 0; i < nx; ++i) {
Fgaussderiv[i + nx * j] = Fgauss[i + nx * j]*KxKxKyKy[i + nx * j];
}
}
//Inverse shift of the frequencies.
//Tested and works as expected.
ifftshift2D(Fgaussderiv, nx, ny);
//Fourier transform back
fftw_execute(plan_backward_gauss);
for (int i = 0; i < nx; ++i) {
for (int j = 0; j < ny; ++j) {
//FFTW scales result.
derivgauss[i + nx*j] = derivgauss[i + nx*j]/(double)nx*(double)ny;
}
}
//Output the data.
char filename[80];
FILE * fp;
sprintf(filename, "/home/stephen/imagedata/%d%09d.txt", 1, 5);
fp = fopen(filename, "w+");
for (i = 0; i < nx* ny; i++) {
fprintf(fp,"%0.15E, %0.15E, %d\n",
creal(derivgauss[i]),creal(derivgaussact[i]),i);
}
fclose(fp);
//Clean up.
free(KxKxKyKy);
free(dderivgauss);
fftw_free(derivgauss);
fftw_free(gauss);
fftw_free(Fgauss);
fftw_free(Fgaussderiv);
fftw_destroy_plan(plan_backward_gauss);
fftw_destroy_plan(plan_forward_gauss);
return (EXIT_SUCCESS);
}