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:

Not much I know but it is a prelude for things to come…

Advertisements

About ateixeira

Physics, Math, Life and Poetry. Lately it has been more about Physics and Math though...
This entry was posted in 00 General, 03 Programming and tagged . Bookmark the permalink.

3 Responses to Mathematica code

  1. Pingback: Cigar-shaped trapped potential | ateixeira

  2. Pingback: Mathematica Resources | ateixeira

  3. Pingback: Source code in C to simulate a Speckle Potential | ateixeira

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s