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: