## Some notes on scattering 01

— 1. Back! —

After a very long hiatus I’m definitely going back to doing research in Physics and blogging. Yes, I know that I promised this a lot of times in the past but this time it is for real.

In December I went to Portugal and talked with my future advisor and we decided this time it’d be for real. We talked about the work he has done lately, namely with his latest PhD. student, and we tried to figure out what would be my research program.

We talked about a few possibilities and narrowed it down to basically two or three options. To tell you the truth I still don’t know which one of the options I’ll end up choosing and this is so because all of the time I was away from Physics.

Thus from the time being I’ll mainly do some catching up and I really hope that in one or two months time I’ll be able to look into an article and figure out what’s going on.

— 2. Reading plan —

In order to do the catching up I mentioned previously I’ll need to read texts mainly on Quantum Mechanics and Field Theory. The texts I’ll read are:

— 3. Scattering Theory —

This set of class notes aims to be a concise introduction to the theory of quantum scattering in the non-relativistic and relativistic formalism.

The second section of the class notes is about the mathematical description of particles in quantum mechanics under the guise of waves. First we describe one dimensional waves and later on the formalism is expanded to include three dimensional waves.

— 3.1. One dimensional wave packet —

In the instant ${t=0}$, a particle can be described in quantum mechanics by the use of a wave packet:

$\displaystyle \psi(x,t=0)=\dfrac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}dk\varphi(k)e^{ikx} \ \ \ \ \ (1)$

Where ${\varphi (k)}$ is the Fourier transform. As an example in the lectures the Fourier transform is taken to be

$\displaystyle \varphi (k) = \begin{cases} \dfrac{1}{\sqrt{2\Delta k}} \quad |k-\bar{k}|\leq 0 \Delta k \\ 0 \quad\quad\quad |k-\bar{k}|>0 \end{cases} \ \ \ \ \ (2)$

Now for that choice of the Fourier transform the wave function becomes

{\begin{aligned} \phi(x) &= \dfrac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}dk\varphi(k)e^{ikx} \\ &= \dfrac{1}{\sqrt{2\pi}}\dfrac{1}{\sqrt{2\Delta k}}\int_{\bar{k}-\Delta k}^{\bar{k}+\Delta k}dke^{ikx} \\ &= \dfrac{1}{2\sqrt{\pi\Delta k}ix}\left( e^{i(\bar{k}+\Delta k)x}-e^{i(\bar{k}-\Delta k)x} \right) \\ &=\dfrac{1}{2\sqrt{\pi\Delta k}ix}e^{ikx}\left( e^{i\Delta kx}-e^{-i\Delta kx} \right) \\ &= \dfrac{1}{2\sqrt{\pi\Delta k}ix}e^{ikx}\left(\cos(\Delta kx)+i\sin(\Delta kx)-\cos(\Delta kx)+i\sin(\Delta kx) \right)\\ &= \dfrac{e^{ikx}}{\sqrt{\pi\Delta k}}\dfrac{\sin(\Delta kx)}{x} \end{aligned}}

— 3.2. Momentum —

The average momentum of a particle is defined in quantum mechanics by the following relationship

$\displaystyle =\int_{-\infty}^{+\infty}dx\psi^*(x)\left\lbrace-i\dfrac{\partial }{\partial x}\right\rbrace\psi(x) \ \ \ \ \ (3)$

For the wave packet expression 1 and the Fourier transform example 2 that was chosen it follows

{\begin{aligned} &=\int_{-\infty}^{+\infty}dx\psi^*(x)\left\lbrace-i\dfrac{\partial }{\partial x}\right\rbrace\psi(x) \\ &= \int_{-\infty}^{+\infty}dx \left\lbrace \dfrac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}dk\varphi^*(k)e^{-ikx} \right\rbrace\left\lbrace-i\dfrac{\partial }{\partial x}\right\rbrace \left\lbrace\dfrac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}dk'\varphi(k')e^{ik'x}\right\rbrace \\ &=\int_{-\infty}^{+\infty}dx \left\lbrace \dfrac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}dk\varphi^*(k)e^{-ikx} \right\rbrace \left\lbrace\dfrac{1}{\sqrt{2\pi}}\int_{-\infty}^{+\infty}dk'\varphi(k')k'e^{ik'x}\right\rbrace\\ &= \int_{-\infty}^{+\infty}dk\int_{-\infty}^{+\infty}dk'\varphi^*(k)\varphi(k')k'\dfrac{1}{2\pi}\int_{-\infty}^{+\infty}dx e^{i(k'-k)x}\\ &=\int_{-\infty}^{+\infty}dk\int_{-\infty}^{+\infty}dk'\varphi^*(k)\varphi(k')k'\delta(k'-k)\\ &= \int_{-\infty}^{+\infty}dk|\varphi(k)|^2k\\ &= \int_{\bar{k}-\Delta k}^{\bar{k}+\Delta k}dk\dfrac{k}{2\Delta k}\\ &= \dfrac{1}{2\Delta k}\left( (\bar{k}+\Delta k)^2-(\bar{k}-\Delta k)^2 \right)\\ &= \dfrac{1}{4\Delta k}(\bar{k}^2+\Delta k^2+2\bar{k}\Delta k-\bar{k}^2-\Delta k^2+2\bar{k}\Delta k)\\ &= \bar{k} \end{aligned}}

The previous derivation shows that the average value of the momentum of a particle is ${\bar{k}}$ which is the central value of the ${k}$ distribution.

As I was studying the lecture notes a part about the previous derivation bothered me. I really didn’t like that we had to use the function given as an example of the Fourier transform to achieve the result. This result is patently valid for a more general class of Fourier transforms and a more general proof must exist.

The facts that make this result stand are:

1. The Fourier transform must be normalizable ${\left(\displaystyle\int_{-\infty}^{+\infty}|\varphi(k)|^2=1\right)}$.
2. The Fourier transform must be symmetric around ${\bar{k}}$.

Introducing the change of variable ${u=k-\bar{k}}$ it is ${du=dk}$ and ${\varphi (u)=\varphi (-u)}$ which means that ${\varphi (u)}$ is an even function.

Therefore our derivation can be:

{\begin{aligned} &=\int_{-\infty}^{+\infty}dk|\varphi(k)|^2k\\ &= \int_{-\infty}^{+\infty}du|\varphi(u+\bar{k})|^2(u+\bar{k})\\ &= \bar{k}\int_{-\infty}^{+\infty}|\varphi(u+\bar{k})|^2+\int_{-\infty}^{+\infty}du |\varphi(u+\bar{k})|^2u\\ &= \bar{k}+0\\ &= \bar{k} \end{aligned}}

Where ${\displaystyle\bar{k}\int_{-\infty}^{+\infty}|\varphi(u+\bar{k})|^2=\bar{k}}$ follows from the fact that ${\varphi (k)}$ is normalized (and a shift doesn’t alter the nature of an integral) and ${\int_{-\infty}^{+\infty}du |\varphi(u+\bar{k})|^2u=0}$ follows from the fact that ${|\varphi(u+\bar{k})|^2u}$ is an odd function.

## Source code in C to simulate a Speckle Potential

Addendum: After writing this post I had a nagging feeling. Some kind of unease about the codes I provided. After taking a quick read of my second edition of Numerical Recipes in C I realized what might have caused that nagging feeling. Copyright! I realized that by posting portions of their code as it was I was violating the terms of he license.

I’m not the biggest fan of the concept of copyright, but I have to respect those that are, and I have to respect them even more when I use their work in ways that are positive to my line of work. Thus I’ve decided to redact everything that is from the Numerical Recipes but I’ll leave the lines of code that belong to me.

And just to clarify things: everything that comes from me on this blog can be used by others as long as they respect the terms put forward in this license: Attribution-NonCommercial-ShareAlike 3.0 Unported (CC BY-NC-SA 3.0)

In this post I provided a code that enable one to compute some properties of fermions when subjected to a cigar shaped trap potential.

Now I’ll provide the code I did back in the day in order to simulate a speckle potential. This code uses some of the routines available in the Numerical Recipes in C and like in the previous posts I’ll also provide the graphs I obtained back then.

I used the files nr.h, nrutil.h and nrutil.c from Numerical Recipes in C.

And this is my code:


#include<stdio.h>
#include<math.h>
#include <stdlib.h>
#include<string.h>
#include<malloc.h>

#include"nr.h"
#include"nrutil.h"

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

/*
Replaces data by its ndim-dimensional discrete Fourier transform, if isign is input as 1.
nn[1..ndim] is an integer array containing the lengths of each dimension (number of complex
values), which MUST all be powers of 2. data is a real array of length twice the product of
these lengths, in which the data are stored as in a multidimensional complex array: real and
imaginary parts of each element are in consecutive locations, and the rightmost index of the
array increases most rapidly as one proceeds along data. For a two-dimensional array, this is
equivalent to storing the array by rows. If isign is input as -1, data is replaced by its inverse
transform times the product of the lengths of all dimensions.
*/

void fourn(float data[], unsigned long nn[], int ndim, int isign);

float ran2(long *idum);
float gasdev(long *idum);

void histogram(float vector[],int num,unsigned long int dimension, char name[]);
void average_intensity(float vector[],int dim);
void probability_density_function(int vector_count[],float vector_classes[], float bin,int nbin,int dim,char file_name[]);

int main()
{
unsigned long int dim=512;
int ndim=2;
int i,j,k,num_aux,l;
unsigned long *nn;
float **grid,**matrix_aux;
/*
**grid is the matrix were we are going to save the values of the field intensities at each point of our grid
**matrix_aux is the matrix that is going to be used on the definition of the window function
*/
float *data; // The values of the electric fields are going to be saved on this vector
float vector_histogram[dim*dim],aux;
long int iseed=-1;
FILE *f_out,*f_intensity;

char nome[35]="Speckle_Potential_";
char nome_2[20]="Intensity_";
char ext[]=".dat";
char num[5];

sprintf(num,"%d",(int)dim);

strcat(nome,num);
strcat(nome,ext);

strcat(nome_2,num);
strcat(nome_2,ext);

num_aux=2*dim*dim;
data=vector(1,num_aux);
nn=lvector(1,ndim);
grid=matrix(1,dim,1,dim);
matrix_aux=matrix(1,dim,1,2*dim);

f_out=fopen(nome,"w+");
f_intensity=fopen(nome_2,"w+");

/*
Initialization of the nn vector with the number of points of each dimension.
*/

for(i=1;i<=ndim;i++)
nn[i]=dim;

/*
Initialization of the elements of the grid
*/

/*
Real part of the electric field initialization
*/

for(i=0; i<=(num_aux-2)/2; i++)
data[i*2+1]=gasdev(&iseed);

/*
Imaginary part of the electric field initialization
*/

for(i=1;i<=num_aux/2;i++)
data[2*i]=gasdev(&iseed);

/*
End of the initialization of the grid
*/

/*
Fast fourier Transform algorithm
*/

fourn(data,nn,ndim,1);

/*
End of the Fourier Transform
*/

/*
Window function definition
*/

for(i=1,k=1;i<=dim;i++)
for(j=1;j<=2*dim;j++)
{
matrix_aux[i][j]=data[k];
k=k+1;
}

for(i=1;i<=dim;i++)
for(j=1;j<=2*dim;j++)
{
if( (i>1 && i<dim) && (j>2 && j<2*dim-1) )
matrix_aux[i][j]=matrix_aux[i][j];
else
matrix_aux[i][j]=0;
}

for(i=1,k=1;i<=dim;i++)
for(j=1;j<=2*dim;j++)
{
data[k]=matrix_aux[i][j];
k=k+1;
}

/*
Inverse Fourier Transform in order to create E'(x,y)
*/

fourn(data,nn,ndim,-1);

/*
E'(x,y) field has been created but lacks the normalization procedure
so we have to normalize it by hand.
*/

for(i=1;i<=num_aux;i++)
data[i]=data[i]/(dim*dim);

for(i=1,k=1,l=0;i<=dim;i++)
for(j=1;j<=dim;j++)
{
grid[i][j]=data[k]*data[k]+data[k+1]*data[k+1];
k=k+2;
fprintf(f_out,"%d\t%d\t%f\n",i,j,grid[i][j]);
fprintf(f_intensity,"%f\n",grid[i][j]);
vector_histogram[l]=grid[i][j];
l=l+1;
}

histogram(vector_histogram,dim,dim*dim,nome_2);
average_intensity(vector_histogram,dim*dim);

free_vector(data,1,num_aux);
free_lvector(nn,1,ndim);
free_matrix(grid,1,dim,1,dim);
free_matrix(matrix_aux,1,dim,1,2*dim);

fclose(f_out);
fclose(f_intensity);

return 0;
}

void histogram(float vector[],int num,unsigned long int dimension, char name[])
{
float class;
float num_max,num_min,aux,sum_aux=0;
FILE *f_hist;
int interval=100,i,j;
float hist[interval];
int count[interval];
char file_name[30]="Histogram_";

strcat(file_name,name);

for(i=0;i<interval;i++)
count[i]=0;

f_hist=fopen(file_name,"w+");

num_max=vector[0];
num_min=vector[0];

for(i=0;i<dimension;i++)
{
aux=vector[i];
if(aux>num_max)
num_max=aux;
}

for(i=0;i<dimension;i++)
{
aux=vector[i];
if(aux<num_min)
num_min=aux;
}
num_min=0;
class=(num_max-num_min)/(interval);

for(i=0;i<interval;i++)
hist[i]=num_min+i*class;

for(i=0;i<dimension;i++)
{
aux=vector[i];
for(j=0;j<interval;j++)
if((aux>=num_min+j*class) && (aux<=num_min+(j+1)*class))
count[j]=count[j]+1;
}

for(i=0;i<interval;i++)
{
sum_aux=0;
for(j=i;j<interval-1;j++)
sum_aux=sum_aux+count[j];

fprintf(f_hist,"%f\t%f\n",(hist[i]+hist[i+1])/2,sum_aux/dimension);
}

probability_density_function(count,hist,class,interval,dimension,name);

fclose(f_hist);
}

void average_intensity(float vector[], int num)
{
int i;
float average=0;
FILE *f_average;

f_average=fopen("Average_Intensity.dat","a+");

for(i=0;i<num;i++)
average=average+vector[i];

fprintf(f_average,"Average intensity in %d*%d grid is %f\n",(int)sqrt(num),(int)sqrt(num),average/num);
}

void probability_density_function(int vector_count[],float vector_classes[], float bin,int nbin,int dim,char file_name[])
{
float function_density[nbin],data_frequency[nbin];
FILE *f_density;
int i;
char name[30]="Density_Function_";

strcat(name,file_name);

f_density=fopen(name,"w+");

for(i=0;i<nbin;i++)
data_frequency[i]=vector_count[i]/(float)dim;

for(i=0;i<nbin;i++)
function_density[i]=data_frequency[i]/bin;

for(i=0;i<nbin;i++)
fprintf(f_density,"%f\t%f\n",vector_classes[i],function_density[i]);

fclose(f_density);

}

void fourn(float data[], unsigned long nn[], int ndim, int isign)
{
This code belongs to Numerical Recipes in C so I've redacted it.
}

float ran2(long *idum)
/*
Long period (> 2^1018) random number generator of L'Ecuyer with Bays-Durham shue
and added safeguards. Returns a uniform random deviate between 0.0 and 1.0 (exclusive of
the endpoint values). Call with idum a negative integer to initialize; thereafter, do not alter
idum between successive deviates in a sequence. RNMX should approximate the largest floating
value that is less than 1.
*/

{
This code belongs to Numerical Recipes in C so I've redacted it.
}

float gasdev(long *idum)
/*Returns a normally distributed deviate with zero mean and unit variance, using ran1(idum)
as the source of uniform deviates.*/
{

This code belongs to Numerical Recipes in C so I've redacted it.
}



We can change the size of the grid that we use to simulate the formation of the speckle pattern. The sizes that I used were 64×64, 128×128, 256×256 and 512×512. This progression allowed me to realize that Statistical Properties of the simulated speckle pattern got closer and closer to the statistical properties of a real speckle pattern.

As can be seen from the graphs the agreement gets better and better as the size of the grid one uses gets bigger.

References:

1. M. Modugno, cond-mat/0509807.
2. The chapter “Statistical Properties of Laser Speckle Patterns” by J. W. Goodman on Laser Speckle and Related Phenomena (Topics in Applied Physics).
Posted in 00 General, 03 Programming | | Leave a comment

## Trick with partial derivatives in Statistical Physics

— 1. The problem —

Some time ago I was studying Statistical Physics and I ended up solving an interesting exercise. The exercise in question is from the book Statistical Physics by Mandl and it is on page 66:

A system possesses three energy levels ${E_1=\epsilon}$, ${E_2=2 \epsilon}$ and ${E_3=3 \epsilon}$, with degeneracies ${g(E_1)=g(E_3)=1}$, ${g(E_2)=2}$. Find the heat capacity of the system.

The first time I solved this problem I did it in a normal way, but then I realized something that allowed to simplify my calculations and solve the problem in a more elegant way.

First we have to calculate the partition function, ${Z}$, for this system.

— 1.1. First Solution —

Since what matters in this problem is the difference of energies we let ${E_2=0}$.

{\begin{aligned} Z &= \displaystyle \sum_{E_r}g(E_r)e^{-\beta E_r} \\ &= e^{\beta \epsilon}+2+2e^{-\beta \epsilon} \\ &= 2(1+ \cosh (\beta \epsilon)) \end{aligned}}

After calculating the partition function we have to calculate the mean energy, ${\bar{E}}$, of this system. By definition it is:

{\begin{aligned} \bar{E} &= -\dfrac{\partial}{\partial \beta} \ln Z \\ &= -\dfrac{\partial}{\partial \beta} \ln 2(1+ \cosh (\beta \epsilon)) \\ &= -\dfrac{2 \epsilon \sinh (\beta \epsilon)}{2(1+\cosh (\beta \epsilon))} \\ &= -\dfrac{\epsilon \sinh (\beta \epsilon)}{1+\cosh (\beta \epsilon)} \end{aligned}}

By definition ${c}$ is

{\begin{aligned} c &= \dfrac{\partial \bar{E}}{\partial T} \\ &= -\dfrac{\partial}{\partial T} \left( \dfrac{\epsilon \sinh (\beta \epsilon)}{1+\cosh (\beta \epsilon)}\right) \end{aligned}}

To differentiate the last expression with respect to ${T}$ isn’t that hard, but it is bothersome and if one isn’t careful the possibility of making a mistake shouldn’t be discarded.

Keeping in mind that in Statistical Physics it is ${\beta = 1/(kT)}$, it follows that (this is the point were I start to use the mathematical trick):

{\begin{aligned} c &= \dfrac{\partial \bar{E}}{\partial T} \\ &= k \dfrac{\partial \bar{E}}{\partial (kT)} \\ &= k \dfrac{\partial \bar{E}}{\partial (1/\beta)} \\ &= -\beta ^2 k \dfrac{\partial \bar{E}}{\partial \beta}\\ &= \beta ^2 k \dfrac{\partial}{\partial \beta}\left( \dfrac{\epsilon \sinh (\beta \epsilon)}{1+\cosh (\beta \epsilon)}\right) \end{aligned}}

Please notice that I haven’t differentiated anything so far. All I’ve done is to change variables in order to calculate the derivative in an easier way.

This last expression is already easy to differentiate but we are caught up in the moment and we’ll just make one more change of variable.

{\begin{aligned} c &= \beta ^2 k \dfrac{\partial}{\partial \beta}\left( \dfrac{\epsilon \sinh (\beta \epsilon)}{1+\cosh (\beta \epsilon)}\right) \\ &= \beta ^2 \epsilon ^2 k \dfrac{\partial}{\partial (\beta \epsilon)}\left( \dfrac{ \sinh (\beta \epsilon)}{1+\cosh (\beta \epsilon)}\right) \\ &= x ^2 k \dfrac{\partial}{\partial x}\left( \dfrac{ \sinh x}{1+\cosh x}\right) \\ &= x ^2 k \dfrac{\cosh x + \cosh ^2 x - \sinh ^2 x}{(1+ \cosh x)^2} \\ &= x ^2 k \dfrac{\cosh x + 1}{(1+ \cosh x)^2} \\ &= \dfrac{x ^2 k}{1+ \cosh x} \end{aligned}}

Yes, in this case the simplification of the hard work involved wasn’t that big, but I think that we should keep in mind this type of reasoning. This way we can apply it to similar scenarios that appear countless times in Physics and Mathematics.

— 1.2. Second Solution —

The previous steps could be even more simplified if we let ${E_1=0}$ and remember an elementary algebraic identity.

{\begin{aligned} Z &= 1+2e^{-\beta \epsilon}+e^{-\beta 2\epsilon} \\ &= (1+e^{-\beta \epsilon})^2 \\ \end{aligned}}

With this expression the calculus of ${\bar{E} }$ and ${c}$ is a lot easier.

{\begin{aligned} \bar{E} &= -\dfrac{\partial}{\partial \beta} \ln Z \\ &= -\dfrac{\partial}{\partial \beta} \ln(1+e^{-\beta \epsilon})^2 \\ &= -2\dfrac{\partial}{\partial \beta} \ln(1+e^{-\beta \epsilon}) \\ &= \dfrac{2\epsilon e^{-\beta\epsilon}}{1+e^{-\beta \epsilon}} \\ &= \dfrac{2\epsilon}{1+e^{\beta \epsilon}} \end{aligned}}

Now for ${c}$ it is

{\begin{aligned} c &= \dfrac{\partial \bar{E}}{\partial T} \\ &= \dfrac{\partial}{\partial T}\left( \dfrac{2\epsilon}{1+e^{\beta \epsilon}}\right) \end{aligned}}

Even though this expression is relatively easy to differentiate with respect to ${T}$ than ${-\dfrac{\epsilon \sinh (\beta \epsilon)}{1+\cosh (\beta \epsilon)}}$ we’ll use the same technique of change of variables that we already know:

{\begin{aligned} c &= \dfrac{\partial \bar{E}}{\partial T} \\ &= k \dfrac{\partial \bar{E}}{\partial (kT)} \\ &= k \dfrac{\partial \bar{E}}{\partial (1/\beta)} \\ &= -\beta ^2 k \dfrac{\partial \bar{E}}{\partial \beta}\\ &= -2\beta^2 \epsilon k \dfrac{\partial}{\partial \beta}\left( \dfrac{1}{1+e^{\beta \epsilon}}\right)\\ &= -2x^2 k \dfrac{\partial}{\partial x}\left( \dfrac{1}{1+e^{x}}\right)\\ &= 2k\dfrac{x^2e^x}{(1+e^{x})^2} \end{aligned}}

An expression that is, apparently, different from ${\dfrac{x ^2 k}{1+ \cosh x}}$.

The proof that these two expressions for ${c}$ are indeed the same is left as an exercise for the reader.

— 2. Appendixes —

In this section we will show that the value of ${c}$ doesn’t depend on what level we consider to be the zero level of energy. We will also show some elementary properties of derivatives taking into account the change of variables.

— 2.1. Derivatives and changes of variables —

Throughout this text we used the properties of changes of variables in derivatives. The goal of this section is to give some quick proofs of the said properties.

Let ${u=u(x)}$. It is an elementary result that ${\dfrac{\partial}{\partial u}=\dfrac{\partial x}{\partial u}\dfrac{\partial}{\partial x}}$.

If it is ${u=kx}$ then

$\displaystyle \dfrac{\partial}{\partial u}=\dfrac{1}{k}\dfrac{\partial}{\partial x} \ \ \ \ \ (1)$

If it is ${u=1/x}$ then

$\displaystyle \dfrac{\partial}{\partial u}=-\dfrac{1}{u^2}\dfrac{\partial}{\partial x} \ \ \ \ \ (2)$

Rewriting equation 1 and equation 2 in order to ${\dfrac{\partial}{\partial x}}$ we get the equalities used in the text.

— 2.2. The invariance of ${c}$

In the two deductions of ${c}$ we took the the ${0}$ level of energy the level that was most convenient. In the interest of having a self contained post we’ll give a justification of why this can be so.

The partition function of the system is ${Z=e^{-\beta \epsilon}+2e^{-\beta 2\epsilon}+e^{-\beta 3\epsilon}}$.

First we let ${E_2=0}$, and the partition function was ${Z_2=e^{\beta \epsilon}+2+e^{-\beta \epsilon}}$; while in the other resolution we let ${E_1=0}$, obtaining the partition function ${Z_1=1+2e^{-\beta \epsilon}+e^{-\beta 2\epsilon}}$.

Hence ${Z_2=e^{\beta E_2} Z}$ and ${Z_1=e^{\beta E_1} Z}$. This method is readily generalized to ${p}$ levels of energy where one has ${Z_n=e^{\beta E_n} Z}$, with ${n\leq p}$.

What we want to show is that for ${Z}$ and ${Z_n}$ it follows ${c_n=c}$.

{\begin{aligned} c_n &= \dfrac{\partial\bar{E}_n}{\partial T} \\ &= -\dfrac{\partial}{\partial T}\dfrac{\partial}{\partial \beta}\ln Z_n \\ &= -\dfrac{\partial}{\partial T}\dfrac{\partial}{\partial \beta}\ln \left( e^{-\beta E_n}Z\right) \\ &= -\dfrac{\partial}{\partial T}\dfrac{\partial}{\partial \beta}\ln e^{-\beta E_n} -\dfrac{\partial}{\partial T}\dfrac{\partial}{\partial \beta}\ln Z \\ &= \dfrac{\partial}{\partial T}\dfrac{\partial}{\partial \beta}(\beta E_n)+c\\ &= \dfrac{\partial}{\partial T} E_n+c\\ &= 0+c\\ &= c \end{aligned}}

## Mathematica Resources

Well just because I suck at Mathematica programming it doesn’t mean you should suck at it too. Nor does it mean that I should suck at it forever.

With that in mind I want to bring to your attention a couple of resources on Mathematica that are certainly very useful.

The first one is An Introduction to Programming with Mathematica which a very nice book with a lot of examples. In Wolfram Library Archive you can find some additional files for the book.

The other book I’ll mention is freely available and quite frankly is very, very good. You can get it here with a lot of additional resources.

I’ll give myself two or three months to read these books and in the meantime I hope to become better at Mathematica.

Posted in 00 General, 02 Literature | Tagged | 2 Comments

## Some more jazz

In a previous post I’ve already mentioned a few books that are useful for people that want to get into Quantum Mechanics and Field Theory.

This time I want to add two more suggestions:

Yep, and this is it. In the nex week or so a few developments on this issue will start to appear in here.

In another news, the semester is finally over (which means that the pains of grading exams are about to come) and that just means more free time for me to get back to shape.

Don’t get me wrong, though! I do love the interaction with all of my students but after a while I just love it a tiny little bit less…

Ps: This website by Gerard ‘t Hooft also has a lot of goodies: HOW to BECOME a GOOD THEORETICAL PHYSICIST.

Posted in 00 General, 02 Literature | Leave a comment

## Some resources on Field Theory and “all” that jazz

Here are some online resources on Field Theory related issues:

Though I’m hardly a specialist in this (yet!) I can see that all of these books are great and would be a great help for anyone wanting to enter this exciting field of research.

Kleinert also has this page were his lecture notes are available. My recommendation is for you to download his Particle and Quantum Fields files.

Posted in 00 General, 02 Literature | 1 Comment

## A quickie on Special Relativity

This post is intended to be a very concise rebuttal to one of the most popular misconceptions about the theory of Special Relativity.

— 1. Brief Historical Review —

For the ones that really want to have a better grasp of this complicated issue you can go to:

But what I really want to say is that we can basically say that the catalyst of Relativity was the Electromagnetic Theory of Maxwell. After deriving the full set Maxwell’s equations Maxwell was able to show that the Electromagnetic Field propagated itself like a wave whose speed was ${c=1/\sqrt{\epsilon_0 \mu_0}}$.

This simple fact had at least two problems to it:

1. According to Galileo’s transformations any speed was relative to the frame were it was being measured. But this value of ${c}$ was an invariant.
2. According to classical Physical ideas if some phenomenon was described by a wave equation than something had to be waving.

In the case of the electromagnetic field that question of what exactly was waving was a very though nut to crack. Just remember the simple fact that we can see the Sun even though the space between us and the Sun is essentially void.

To solve these two problems with one go the Physicists proposed that there was a substance that existed all over the space that waved whenever the electromagnetic radiation was being carried. And the speed of light was indeed an invariant, but it was an invariant relatively to the aether…

But that of course begged the question of why didn’t we find any variation in the speed of light since the Earth’s movement relative to the aether was changing all the time?

The answer was that the methods weren’t good enough to detect such a little variation. And the hunt was on to devise a method with enough sensibility to detect the more than obvious changes.

A lot of things were tried, all of them failed, a lot of partial answers were proposed to explain this failure but a coherent picture was missing!

— 2. Enter Einstein —

This of course changed in 1905 when a young man named Einstein publish some of the best articles in the history of our business.

One of these articles was a major rethinking of the concepts of space and time and really is a beauty to be read for any self-respecting physicist. But I digress…

With just two simple (and seemingly contradictory axioms) Einstein was able to show why no variation in the value of the speed of light could be detected and he also predicted many more phenomena that initiated a whole new cultural Zeitgeist for the years to come.

His two principles were:

1. The same laws of electrodynamics and optics will be valid for all frames of reference for which the equations of mechanics hold good.
2. Light is always propagated in empty space with a definite velocity c which is independent of the state of motion of the emitting body

— 3. The misconception —

After this somewhat lengthy introduction is time for me to state what I consider to be the popular misconception: “The theory of relativity says that nothing can travel faster than light.”

Well, we surely don’t see that being said by Einstein in the previous quote! What he says is that the speed of light is an invariant (even though even in University textbooks we can read the first postulate as saying that the speed of light is the limit of all speeds in the Universe and that it is the same for all observers).

— 3.1. Moving Frames —

Well if it isn’t a postulate (if it were a postulate it would always be up for experimental scrutiny) maybe it is something that we can derive in Special Relativity?

Well, sort of…

Imagine that you have two inertial frames ${S}$ and ${S'}$ and you consider ${S}$ to be a stationary frame and that ${S'}$ has speed ${v}$ relative to ${S}$. If you know that the speed of a given body to be ${u'_x}$ in ${S'}$ than you can show that the velocity of the same body in the inertial frame ${S}$, ${u_x}$, to be:

$\displaystyle u_x=\frac{u'_x+v}{1+\frac{u'_xv}{c^2}}$

Now if you admit that ${u'_x < c}$ and that ${v < c}$ you can show that ${u_x < c}$. But if you don’t admit the initial hypothesis than the value of ${u_x}$ isn’t bounded above by ${c}$.

And now for a proof of the previous assertion:

If ${u'x < c}$ and ${v < c}$, that for some ${\alpha, \beta > 0}$ it is ${u'_x=c-\alpha}$ and ${v =c-\beta}$.

Now ${u_x}$ is

{\begin{aligned} u_x &= \dfrac{c-\alpha+c-\beta}{1+(c-\alpha)(c-\beta)/c^2} \\ &= \dfrac{2c-(\alpha+\beta)}{1+(c^2-(\alpha+\beta)+\alpha\beta)/c^2} \\ &= \dfrac{2c-(\alpha+\beta)}{2c^2-(\alpha+\beta)c+\alpha\beta}c^2 \\ &< \dfrac{2c-(\alpha+\beta)}{2c^2-(\alpha+\beta)c}c^2 \\ &= \dfrac{2c-(\alpha+\beta)}{2c-(\alpha+\beta)}c \\ &= c \end{aligned}}

Which is ${u_x < c}$ just like we said it would be.

Again all of that this shows is that one can’t boost a body to have a speed bigger than ${c}$ by using an inertial frame whose speed is smaller than ${c}$ if the speed of the body in ${S'}$ is smaller than ${c}$ too.

— 3.2. Accelerating bodies —

But we still have another chance, though. What if we don’t use boosts? What if we accelerate a particle during a sufficient long time interval? Doesn’t a speed bigger than ${c}$ happens eventually?

Let’s see.

First we have to define the linear momentum of a particle in Special Relativity. If we want the linear momentum to be conserved in relativistic collisions the definition that makes sense is

$\displaystyle \vec{p}=\gamma m \vec{u}$

Where ${\gamma}$ has its usual meaning and ${\vec{u}}$ is the speed of the body.

The definition of force still is ${\vec{F}=\dfrac{d\vec{p}}{dt}}$. Thus ${\vec{F}=\dfrac{d}{dt}(\gamma m \vec{u})}$.

It is possible to show that for a body of constant mass (which is always the case if something doesn’t get aggregated to the body or if the body doesn’t get broken up – summing up: we aren’t falling for the “relativistic mass” erroneous idea) one has

$\displaystyle a=\frac{F}{m}(1-u^2/c^2)^{3/2}$

Now imagine that you are pushing on a body and its speed is getting bigger and bigger. Than what we have is ${u \rightarrow c}$ but in this case it is ${a \rightarrow 0}$. Thus the closer the speed of the body gets to ${c}$ the less it is accelerating until that finally its acceleration is ${0}$.

The conclusion is that if a body starts with ${u < c}$ than it would take an infinite amount of time to get to a speed of ${c}$ and from then on its speed would be ${c}$ for ever.

— 4. Conclusion —

Contrary to what is normally said the fact that the speed of light is limit of all speeds on the real world isn’t a postulate.

It does hold with the caveat that the bodies have to have speeds less than ${c}$ to begin with, but nothing in Special Relativity forces the bodies to have speeds less than ${c}$ when they start moving.

— 5. Appendix —

I think that in order to finish this post I just have to mention that it is possible to derive the theory of Special Relativity with just the first postulate and one doesn’t need to make any mention to the speed of light.

The idea is the following one: imagine that you have three inertial frames ${S_1}$, ${S_2}$ and ${S_3}$ whose axes are parallel to each other.

Then you boost ${S_2}$ relatively to ${S_1}$ with a velocity ${v_{21}}$.

You also boost ${S_3}$ relatively to ${S_2}$ with a velocity ${v_{32}}$.

The question now is: Are the axes of ${S_1}$ and ${S_3}$ still parallel?

If you answered the previous question with a resounding yes, than what you’re doing is saying that we always live an Euclidean timespace and that the dynamics that rules us are the classic ones.

If on the other hand you just said: Gee, I don’t really know… than you’re saying that is none of our business to tell the spacetime continuum how it ought to behave and you are willing to study it and only then conclude what its true nature is.

This is just what Feigenbaum did in this great article: The Theory of Relativity – Galileo’s Child whose reading I highly reccomend and whose abstract I just have to quote (no bold in the original):

We determine the Lorentz transformations and the kinematic content and dynamical framework of special relativity as purely an extension of Galileo’s thoughts. No reference to light is ever required: The theories of relativity are logically independent of any properties of light. The thoughts of Galileo are fully realized in a system of Lorentz transformations with a parameter ${1/c^2}$, some undetermined, universal constant of nature; and are realizable in no other. Isotropy of space plays a deep and pivotal role in all of this, since here three-dimensional space appears at first blush, and persists until the conclusion: Relativity can never correctly be fully developed in just one spatial dimension.

## New avenues into Quantum Mechanics

In recent times two articles that can do the seemingly impossible in Quantum Mechanics have been published and they generated some buzz on the interweb.

This first article Observing the Average Trajectories of Single Photons in a Two-Slit Interferometer is notable because the authors show how they were able to observe the trajectories of photons in a double slit experiment and still managed to observe a clear interference pattern. A thing that is impossible to do according to the Complementarity principle.

In the second article, Direct measurement of the quantum wavefunction, the authors computed the the transverse spatial wavefunction of a single photon by means that they consider to be direct. For me, that haven’t read the article, so far it doesn’t seem to be a so direct method as claimed, but nevertheless the level of experimental expertise even to get an indirect computation of the wave function is certainly worthy of respect.

These spectacular achievements were possible because these two teams used weak measurement techniques, together with statistical ensembles of photons and non simultaneous measurements of the complementary variables they set out to determine.

The two abstracts are here (the bold isn’t in the original):

• Direct measurement of the quantum wavefunction

The wavefunction is the complex distribution used to completely describe a quantum system, and is central to quantum theory. But despite its fundamental role, it is typically introduced as an abstract element of the theory with no explicit definition. Rather, physicists come to a working understanding of the wavefunction through its use to calculate measurement outcome probabilities by way of the Born rule. At present, the wavefunction is determined through tomographic methods which estimate the wavefunction most consistent with a diverse collection of measurements. The indirectness of these methods compounds the problem of defining the wavefunction. Here we show that the wavefunction can be measured directly by the sequential measurement of two complementary variables of the system. The crux of our method is that the first measurement is performed in a gentle way through weak measurement so as not to invalidate the second. The result is that the real and imaginary components of the wavefunction appear directly on our measurement apparatus. We give an experimental example by directly measuring the transverse spatial wavefunction of a single photon, a task not previously realized by any method. We show that the concept is universal, being applicable to other degrees of freedom of the photon, such as polarization or frequency, and to other quantum systems ? for example, electron spins, SQUIDs (superconducting quantum interference devices) and trapped ions. Consequently, this method gives the wavefunction a straightforward and general definition in terms of a specific set of experimental operations. We expect it to expand the range of quantum systems that can be characterized and to initiate new avenues in fundamental quantum theory.

• Observing the Average Trajectories of Single Photons in a Two-Slit Interferometer

A consequence of the quantum mechanical uncertainty principle is that one may not discuss the path or ?trajectory? that a quantum particle takes, because any measurement of position irrevocably disturbs the momentum, and vice versa. Using weak measurements, however, it is possible to operationally define a set of trajectories for an ensemble of quantum particles. We sent single photons emitted by a quantum dot through a double-slit interferometer and reconstructed these trajectories by performing a weak measurement of the photon momentum, postselected according to the result of a strong measurement of photon position in a series of planes. The results provide an observationally grounded description of the propagation of subensembles of quantum particles in a two-slit interferometer.

And an excellent explanation of why this was possible can be found in this blog post Watching Photons Interfere: “Observing the Average Trajectories of Single Photons in a Two-Slit Interferometer”

Posted in 03 Physics, 03 Quantum Mechanics | 1 Comment

## Cigar-shaped trap potential

At more or less at the same time that I did the Mathematica code to simulate a speckle pattern producing potential I also had to create a code that recreated a cigar-shaped trap potential. This was done with the goal of calculating a few properties of fermions.

The idea was to first start with a cigar shaped potential and then superimpose a random potential on it (the strength of this random potential would be variable) so that we could see what would happen to the properties of the fermions (the sample that we used was of 10000 fermions) as the randomness would get stronger and stronger.

This post is a preamble and the first in a series that will attempt to derive and explain a few things about this line of work that I tried to follow in my past.

For now I’ll only post the C code I wrote back then, and the graphics I obtained but in future posts I’ll talk more about this, because, even though it isn’t related with my current line of work, it is an interesting topic.

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

/*********************************************
* Introduction of experimental parameters   *
*********************************************/

float h=1.00;
float Mf=1.00; 	/* This is the fermion mass. */
float wlong=80.00;
float wtrans=200.00; /* wlong is the frequency on the x-axis and wtrans
the frequency on the y-z plane. */
float delta_x=1.0;
float delta_r=0.1;
float L;
float l;

int n_sample=10000; /* This is the number of fermions given in the sample. */

double Ef,E0,Smax,sqrtPi;

/*********************************************************
* End of the introduction the experimental parameters   *
*********************************************************/

const int Num=2600; /* Num will be the variable that will permit us
to calculate the number of fermions in function
of the fermi energy. */

double hermitepoly(int n, double x); /* Function that will calculate the Hermite
polynomial for a given degree and given point.*/

unsigned long int N_fermiEnergy(double a); /*Function that will calculate the number of fermions
for each fermi energy. */

double fermienergy(void); /* By examining where this functions changes signs we will calculate
the fermi energy of a sample given an initial number of fermions. */

double factorial(int n); /*Calculates the factorial of a given number. */

/* The three wave functions that we obtain after separation of variables. */
/* And th~ese are the properly squared nondimensional functions already. */

double psi1_x(int mx, double x);
double psi2_y(int my, double y);
double psi3_z(int mz, double z);

/* Functions that will allow us to plot the fermion density by calculating
their contribution in the expression of the fermion density*/

double Ns(int n);

double density(double x, double y, double z);

void plot_distribution(void); /* Function that will "plot" the fermion density. */

void integral(void); /* Functional that will calculate the number of fermions acording
to function distribution to certify that the program really works. */

int main(void)
{
E0 = h*(wlong/2 + wtrans);
Ef=fermienergy();
printf("Ef=%lf\n",Ef);
Smax = (Ef-E0)/(h*wtrans)+1;

sqrtPi=sqrt(3.1415927);

L=sqrt(h/(Mf*wlong));
l=sqrt(h/(Mf*wtrans));

plot_distribution();
integral();

return 0;
}

double hermitepoly (int n, double x)
{
double a;
double b;
double c;
int i;

if (n == 0)
return 1;
if (n == 1)
return 2 * x;
a = 1.0;
b = 2.0 * x;
for (i = 2; i <= n; i++)
{
c = b;
b = 2.0 * x * b - 2.0 * (i - 1.0) * a;
a = c;
}

return b;

}

unsigned long int N_fermiEnergy(double a)
{
unsigned long int *Num_niveis;
unsigned long int Ntot=0;
double smax;
int i;

smax=(int)((a-E0)/(h*wtrans)+1);

Num_niveis = (unsigned long int*) malloc((int)smax);

for (i=0; i < smax; i++)
{

Num_niveis[i]=(i+1)*((int)((a-(E0+i*h*wtrans))/(h*wlong))+1);

Ntot = Num_niveis[i] + Ntot;
}

free(Num_niveis);

return Ntot;

}

double fermienergy(void)
{

double aux,x1=0,x2=0;
long int k;
unsigned long int n_fermions;

for(aux=h*(wlong/2 + wtrans); aux<Num; aux=aux+wlong/12)
{
n_fermions=N_fermiEnergy(aux);

k = n_fermions-n_sample;
if(k==0)
return aux;
else
{
if(k<0)
{
x1=aux;
x2=aux+wlong/12;
}

}
}

return (x1+x2+wlong/6)/2;

}

double factorial(int n)
{
double res=1;
int i;

for(i=1;i<=n;i++)
res=res*i;

return res;
}

double psi1_x(int mx, double x)
{
return (l/L)*pow(hermitepoly(mx,x*l/L),2)*(1/(pow(2,mx)*factorial(mx)*sqrtPi))*exp(-(pow(x,2)*pow(l,2)/pow(L,2)));
}

double psi2_y(int my, double y)
{
return pow(hermitepoly(my,y),2)*((1/(pow(2,my)*factorial(my)*sqrtPi))*exp(-pow(y,2)));
}

double psi3_z(int mz, double z)
{
return pow(hermitepoly(mz,z),2)*((1/(pow(2,mz)*factorial(mz)*sqrtPi))*exp(-(pow(z,2))));
}

double Ns(int s)
{
return (((Ef-(E0+s*h*wtrans))/(h*wlong))+1);
}

double density(double x, double y, double z)
{
int s, t, p,n;
double squared_yz,squared_x, squared=0;

for(s=0;s<=Smax;s++)
{
n=(int)Ns(s);
squared_yz=0;
squared_x=0;

for(p=0;p<n;p++)
squared_x = squared_x + psi1_x(p,x);

for(t=0;t<=s;t++)
squared_yz = squared_yz + psi2_y(t,y)*psi3_z(s-t,z);

squared = squared + squared_x*squared_yz;
}
return squared;
}

void plot_distribution()
{
double i,j; /*i is the variable that will provide the steps in the x direction
and j the one that will provide the steps in the "r direction". */
FILE *fp;
char X='x';
char R='r';
char str[]= "Density";
int dim_x=14;
int dim_r=7;

fp=fopen("fermion_y_10000.dat","w+");

if(fp==NULL)
printf("File impossible to open!!!\n");
else
{

fprintf(fp,"%c\t%c\t%s\n",X,R,str);
for(i=-dim_x;i<=dim_x;i=i+delta_x)
for(j=0;j<=dim_r;j=j+delta_r)
{
fprintf(fp,"%lf\t%lf\t%lf\n",i,j,density(i,j,0.0));
}
}

fclose(fp);
}

void integral(void)
{
double i,j; /*i is the variable that will provide the steps in the x direction
and j the one the will provide the steps in the "r direction". */

double fermion_number=0;
double base_area=delta_r*delta_x;
int dim_x=14;
int dim_r=7;

for(i=-dim_x;i<=dim_x;i=i+delta_x)
for(j=0;j<=dim_r;j=j+delta_r)
fermion_number = fermion_number + j*density(i+delta_x/2,j+delta_r/2,0.0)*base_area;

printf("The number of fermions is %lf.\n",2*3.1415927*fermion_number);
}



To finalize this post I’ll just show a few graphics that were produced:

Posted in 00 General, 03 Programming | 1 Comment

## 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…

Posted in 00 General, 03 Programming | Tagged | 3 Comments