Mathematica code

A long time ago I had to make a program that would simulate a speckle pattern. I ended up doing it in C and in Mathematica.

Now I’ll just post in here my Mathematica code to see how it looks like in wordpress.com style. I’ll be following the guidelines in this page.

Without further ado here is the code:

```ClearAll["Global`*"]

dim = 64

grid64 = Table[
RandomReal[NormalDistribution[0, 1]] +
I*RandomReal[NormalDistribution[0, 1]], {i, dim}, {j, dim}];

fourier64 = Fourier[grid64];

(*This part of the code will define the window function. Notice that the \
values of i and j can be changed and this shouldn't have any effect on the \
final result*)

For[i = 1, i <= dim, i++,
For[j = 1, j <= dim, j++,
If[(i > 1 && i < dim) && (j > 1 &&  j < dim),
fourier64[[i, j]] = fourier64[[i, j]], fourier64[[i, j]] = 0]]]

finalfield64 = InverseFourier[fourier64];

gridintensity64 = Table[, {i, dim}, {j, dim}];

(*This part of the code calculates the speckle field intensity in each point \
of the grid.*)

For[i = 1, i <= dim, i++,
For[j = 1, j <= dim, j++,
gridintensity64[[i, j]] = Abs[finalfield64[[i, j]]]^2   ]]

list64 = Table[0, {i, dim*dim}];

(* In this part of the code we will save the values of the various \
intensities at each point of the grid by the following order: Running every \
column for each row as the indices go from 1 to dim! *)

For[i = 1; k = 1, i <= dim, i++,
For[j = 1, j <= dim, j++,
list64[[k]] = gridintensity64[[i, j]]; k = k + 1
]]

Export["Intensity_64.dat", list64]

(*The inverse of average64 should appear as the coefficient in the negative \
exponential that expresses what's the probability that a given the intensity \
excedes a certain value of I*)

average64 = Mean[list64]

(*Now will start the study of the statistical properties of the intensities*)

MinIntensity64 = Min[list64]

MaxIntensity64 = Max[list64]

interval = 100

bin = (MaxIntensity64 - MinIntensity64)/interval

(*HistogramX will be the array that will keep the values of the intervals for \
the intensities. HistogramY will be the array that will keep the relevant \
probability of a given intensity interval*)

HistogramX = Table[0, {i, interval}];

HistogramY = Table[0, {i, interval}];

For[i = 1, i <= interval, i++,
HistogramX[[i]] = MinIntensity64 + (i - 1)*bin]

HistogramY = BinCounts[list64, {MinIntensity64, MaxIntensity64, bin}];

(*Let us normalize the previous results so that we can have probabilities*)

HistogramY = 1.0*HistogramY/(dim*dim);

IntensityHistogram64 = Table[0, {i, 2*interval}];

For[i = 1, i <= interval, i++,
IntensityHistogram64[[2*(i - 1) + 1]] = HistogramX[[i]];
IntensityHistogram64[[2*i]] = HistogramY[[i]] ]

IntensityHistogram64 = Partition[IntensityHistogram64, 2];

graphIntensityHistogram64 =
ListPlot[IntensityHistogram64, PlotLabel -> "Intensity Histogram 64*64"]

Export["Intensity_Histogram_64.dat", IntensityHistogram64]

Export["Intensity_Histogram_64.gif", graphIntensityHistogram64]

(* Calculus of the Probability Density Function of our generated intensity \
values*)

h = bin

AcumulatedIntensities64 = Table[0, {i, interval}];

ExcedeIntensity64 = Table[0, {i, 2*interval}];

For[i = 1, i <= interval, i++,
For[k = i, k <= interval, k++,
AcumulatedIntensities64[[i]] =
AcumulatedIntensities64[[i]] + HistogramY[[k]]  ]]

For[i = 1, i <= interval, i++,
ExcedeIntensity64[[2*(i - 1) + 1]] = HistogramX[[i]];
ExcedeIntensity64[[2*i]] = AcumulatedIntensities64[[i]] ]

ExcedeIntensity64 = Partition[ExcedeIntensity64, 2];

graphExcedeIntensity64=ListPlot[ExcedeIntensity64, PlotLabel ->
"Probability that the intensity exceeds I",PlotRange->{{0,MaxIntensity64},{0,1}}]

Export["Excede_intensity_64.dat", ExcedeIntensity64]

Export["Excede_intensity_64.gif", graphExcedeIntensity64]

grafcontour64 =
ListContourPlot[gridintensity64, PlotLabel -> "Contour Plot 64x64 grid"]

Export["Contour_Plot_64.gif", grafcontour64]

Correlation64 =
ListCorrelate[gridintensity64, gridintensity64, {-1, 1}, 0]/(dim*dim);

Correlation64Plot =
ListContourPlot[Correlation64, PlotLabel -> "Autocorrelation Function"]

Export["Correlation_Plot_64.gif", Correlation64Plot]

(*It looks like the expected autocorrelation but it is much more angulous \
than the real one and also the maximum is much more evident. This is due to \
size of the grid and the spacing we took between points. Nevertheless it \
seems to be in good agreement. *)
```

Please don’t make fun of my code as this was was written a long time ago and my Mathematica skills weren’t that good back then (they still aren’t but let’s just ignore that for the time being).

And in case you want to see the graphs I managed to produce here they are: 