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.
Then I posted a Mathematica code that allows one one to simulate a speckle 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:
, 
and 
, with degeneracies 
, 
. Find the heat capacity of the system.
, for this system.
.
, of this system. By definition it is:
is
isn’t that hard, but it is bothersome and if one isn’t careful the possibility of making a mistake shouldn’t be discarded.
, it follows that (this is the point were I start to use the mathematical trick):
and remember an elementary algebraic identity.
and 
we’ll use the same technique of change of variables that we already know:
.
. It is an elementary result that 
.
then
then
we get the equalities used in the text.
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.
.
; while in the other resolution we let 
.
and 
. This method is readily generalized to 
levels of energy where one has 
, with 
.
it follows 
.
