c_code/main.c

//Compile as gcc -o bec main.c -lm
//Execute as ./bec in Linux

//--------------Preamble------------
#include <stdlib.h>
#include <stdio.h>
#include <complex.h>
#include <time.h>

//Datatype Aliases
#define TRUE (0 == 0)
#define FALSE (1 == 0)
typedef double complex fieldType;
typedef double parameterType;

//Global Constants
const int Size = 180;
const int maxTimesteps = 1000;
const int shotImageAt = 500;
const fieldType xy0 = 20;
const fieldType dt = 1e-3;
const parameterType tolerance = 1e-10;
const parameterType windingNumber = 1.0;
const fieldType complexZero = 0.0;

//Global Parameters
const fieldType g = 1.0;
const fieldType hbar = 1.0;
const fieldType smallOmegaSq = 0.01;
const parameterType mu = 1.0;
const fieldType m = 1.0;
const fieldType Omega = 0.0;
const parameterType phaseStepWidth = 10;
const parameterType phaseStepHeight = 2.5;
const parameterType phaseStepPosition = 80;

//Function Prototypes
void init(fieldType field[][Size], fieldType value);
void init1D(parameterType F[], parameterType value);
void copy(fieldType source[][Size], fieldType dest[][Size]);
void loadParameters(const char *filename); //Not Functional ATM
void steadyState(fieldType field[][Size], parameterType F[], parameterType dr);
void steadyVortexState(fieldType field[][Size], parameterType F[], parameterType dr);
void transferRadial(parameterType F[], fieldType field[][Size], parameterType xOffset, parameterType yOffset, parameterType dx, parameterType dy, parameterType dr);
void setupUniformPhase(fieldType field[][Size], parameterType phase, parameterType xOffset, parameterType yOffset, parameterType dx, parameterType dy);
void setupVortexPhase(fieldType field[][Size], parameterType xOffset, parameterType yOffset, parameterType dx, parameterType dy);
void applyPhaseStep(fieldType field[][Size]);
void evolveUsingTDGP(fieldType field[][Size], fieldType oldfield[][Size], fieldType xOffset, fieldType yOffset, fieldType h_t, fieldType h_x, fieldType h_y);
void setBoundaryFree(fieldType field[][Size]);
void setBoundaryVanish(fieldType field[][Size]);
void solveTridiagonal(fieldType a[], fieldType b[], fieldType c[], fieldType unknowns[], fieldType knowns[]);
void solveTridiagonal2(fieldType a[], fieldType b[], fieldType c[], fieldType unknowns[], fieldType knowns[]);
void writeImage(fieldType field[][Size]); //Not Functional ATM
void writeRadialSolution(parameterType F[], parameterType dr);
void writeResultsCSV(fieldType field[][Size]);
void writeResultsGNU(fieldType field[][Size]);
//----------------------------------

//======================Main=========================
int main()
{
  fieldType Field[Size][Size], OldField[Size][Size];
  parameterType radialSolution[Size];
  parameterType spatialStep = 2.0*(double)xy0/(double)Size, timeDiff = 0.0;
  parameterType g, mu, smallOmegaSq;
  int p = 0;
  time_t start,end;

  //Initialize.
  time(&start);
  init(Field,0.0);
  init(OldField,0.0);
  init1D(radialSolution,1.0);

  /*for(p = 0; p < Size; p ++)
    fprintf(stderr,"%e, ",radialSolution[p]);
    fprintf(stderr,"\n");*/

  printf("---------------------------------------------------------\n");
  printf("||\t BEC Simulation\n");

  //Calculate and Setup Steady State.
  //steadyState(Field,radialSolution,spatialStep);
  steadyVortexState(Field,radialSolution,spatialStep);
  transferRadial(radialSolution,Field,xy0,xy0,spatialStep,spatialStep,spatialStep);

  //Setup Initial Phase
  //setupUniformPhase(Field,1.0,xy0,xy0,spatialStep,spatialStep);
  setupVortexPhase(Field,xy0,xy0,spatialStep,spatialStep);
  applyPhaseStep(Field);

  //Evolve Field according TDGP Equation using parameters provided.
  fprintf(stderr,"||\t Evolving... ");
  for(p = 0; p < maxTimesteps; p ++)
    evolveUsingTDGP(Field,OldField,xy0,xy0,dt,spatialStep,spatialStep);
  fprintf(stderr,"Done\n");

  //Write Results of Simulation.
  writeResultsCSV(Field);
  writeRadialSolution(radialSolution,spatialStep);

  //Calculate Time Taken and Close Application.
  time(&end);
  timeDiff = difftime(end,start);
  printf("||\t The Process took %f minutes (%f hours) to Complete.\n", (double)(timeDiff/60), (double)(timeDiff/3600));
  printf("---------------------------------------------------------\n");
  return 0;
}
//===================================================

//Function Bodies
void init(fieldType field[][Size], fieldType value)
{
  int k,j;

  for(k = 0; k < Size; k ++)
    for(j = 0; j < Size; j ++)
      field[j][k] = value;
}

void init1D(parameterType F[], parameterType value)
{
  int j;

  for(j = 0; j < Size; j ++)
    F[j] = value;
}

void copy(fieldType source[][Size], fieldType dest[][Size])
{
  int j, k;

  for(k = 0; k < Size; k++)
    for(j = 0; j < Size; j++)
      dest[j][k] = source[j][k];
}

void loadParameters(const char *filename)
{
  FILE *outFile;
  int j,k;

  outFile = fopen(filename,"r");

  fclose(outFile);
}

void steadyState(fieldType field[][Size], parameterType F[], parameterType dr)
{
  parameterType sum = 0.0, oldSum = 1.0, r = 0.0;
  int j, radialExtent = Size-1;

  while(fabs(sum - oldSum) > tolerance && radialExtent > 0)
  {
    oldSum = sum;
    sum = 0.0;
    for(j = 1; j < radialExtent; j ++)
      {
        r = j*dr;
        F[j] = 0.5*(F[j+1]+F[j-1]+0.5*dr*(F[j+1]-F[j-1])/j) - dr*dr*F[j]*(g*F[j]*F[j] + 0.5*smallOmegaSq*r*r - mu);
        sum += F[j];
      }
    F[radialExtent] = 0.0;
    F[0] = 2*F[1] - F[1];
  }
}

void steadyVortexState(fieldType field[][Size], parameterType F[], parameterType dr)
{
  parameterType sum = 0.0, oldSum = 1.0, r = 0.0, n = windingNumber, internal = 0.0, frac1 = 0.0, frac2 = 0.0, jSq = 0.0, extra = 0.0, jFrac = 0.0;
  int j, radialExtent = Size-1, count = 0;

  if(radialExtent > 0)
    F[0] = 0.0;

  while(fabs(sum - oldSum) > tolerance && radialExtent > 0)
  {
    oldSum = sum;
    sum = 0.0;
    for(j = 1; j < radialExtent; j ++)
      {
        r = j*dr;
        internal = g*F[j]*F[j] + 0.5*smallOmegaSq*r*r - mu + n*Omega;
        jFrac = (F[j+1]-F[j-1])/j;
        frac1 = 0.5*dr*jFrac;
        jSq = (parameterType)j*j;
        frac2 = (n*n)/(jSq)*F[j];
        extra = F[j+1]+F[j-1];
        F[j] = 0.5*( extra+frac1-frac2 ) - dr*dr*F[j]*internal;
        sum += F[j];
        //fprintf(stderr,"%e \n", F[j]);
      }
    //exit(0);
    //fprintf(stderr,"%d ", count);
    //count ++;
    //F[radialExtent] = F[radialExtent-1];
    F[radialExtent] = 0.0;
  }
  //fprintf(stderr,"Tolerance Met!\n");
}

void transferRadial(parameterType F[], fieldType field[][Size], parameterType xOffset, parameterType yOffset, parameterType dx, parameterType dy, parameterType dr)
{
  parameterType x, y;
  parameterType l, rSquared;
  parameterType magnitude;
  int j, k, rInteger;

  for(k = 0; k < Size; k++)
    {
      y = (k)*dy - yOffset;
      for(j = 0; j < Size; j++)
        {
          x = (j)*dx - xOffset;
          rSquared = x*x + y*y;
          l = sqrt(rSquared)/dr;
          rInteger = (int)(l);
          magnitude = F[rInteger] + (l-rInteger)*( F[rInteger+1] - F[rInteger] );
          field[j][k] = magnitude;
        }
    }
}

void setupUniformPhase(fieldType field[][Size], parameterType phase, parameterType xOffset, parameterType yOffset, parameterType dx, parameterType dy)
{
  parameterType x, y;
  int j, k;

  for(k = 0; k < Size; k++)
    {
      y = k*dy - yOffset;
      for(j = 0; j < Size; j++)
        {
          x = j*dx - xOffset;
          field[j][k] = cabs(field[j][k])*cexp(I*phase);
        }
    }
}

void setupVortexPhase(fieldType field[][Size], parameterType xOffset, parameterType yOffset, parameterType dx, parameterType dy)
{
  parameterType x, y;
  int j, k;

  for(k = 0; k < Size; k++)
    {
      y = k*dy - yOffset;
      for(j = 0; j < Size; j++)
        {
          x = j*dx - xOffset;
          field[j][k] = cabs(field[j][k])*cexp(I*atan2(y,x));
        }
    }
}

void applyPhaseStep(fieldType field[][Size])
{
  int j, k;
  parameterType endCounter = (double)(phaseStepWidth), counter = 1.0, gradientFactor = 0.0;

  for(k = phaseStepPosition; k < Size; k ++)
    {
      for(j = 0; j < Size; j ++)
        {
          gradientFactor = counter/endCounter;
          field[j][k] = field[j][k]*cexp(I*gradientFactor*phaseStepHeight);
          //field[j][k] = cabs(field[j][k])*cexp(I*carg(field[j][k]))*cexp(I*gradientFactor*phaseStepHeight);
        }
      if(counter < endCounter)
        counter ++;
    }
}

void evolveUsingTDGP(fieldType currentField[][Size], fieldType oldField[][Size], fieldType xOffset, fieldType yOffset, fieldType h_t, fieldType h_x, fieldType h_y)
{
  int maxFSIIterations = 4, xExtent = Size-1, yExtent = Size-1, j, k, p;
  fieldType a[xExtent], b[xExtent], c[xExtent], d[xExtent], X[xExtent];
  fieldType a1[yExtent], b1[yExtent], c1[yExtent], d1[yExtent], X1[yExtent];
  fieldType x, y;
  x = 4.0*h_x*h_x*m; //x is used as tmp storage
  y = 4.0*h_y*h_y*m; //y is used as tmp storage
  fieldType     alpha_x = I*((h_t*hbar)/x),
    alpha_y = I*((h_t*hbar)/y),
    Psi, oldPsi;

    fieldType dtCoEff1, dtCoEff2;

  copy(currentField,oldField);

  for(p = 0; p < maxFSIIterations; p ++)
    {
      for(k = 1; k < yExtent; k ++)
        {
          y = (k)*h_y - yOffset;
          for(j = 1; j < xExtent; j ++)
            {
              x = (j)*h_x - xOffset;
              Psi = currentField[j][k];
              oldPsi = oldField[j][k];

              dtCoEff1 = h_t*0.5*Psi*I;
              dtCoEff2 = h_t*0.5*oldPsi*I;

              a[j] = -alpha_x;
              b[j] = 1.0 + 2.0*alpha_x;
              c[j] = -alpha_x;
              //d[j] = oldPsi + alpha_x*( oldField[j+1][k] + oldField[j-1][k] + oldField[j][k+1] + oldField[j][k-1] - 4.0*oldPsi) - h_t*0.5*I/hbar*Psi*(0.5*m*smallOmegaSq*(x*x + y*y) + g*cabs(Psi)*cabs(Psi)) - h_t*0.5*I/hbar*oldPsi*(0.5*m*smallOmegaSq*(x*x + y*y) + g*cabs(oldPsi)*cabs(oldPsi));
              d[j] = oldPsi + alpha_x*( oldField[j+1][k] + oldField[j-1][k] + oldField[j][k+1] + oldField[j][k-1] - 4.0*oldPsi) - dtCoEff1/hbar*(0.5*m*smallOmegaSq*(x*x + y*y) + g*cabs(Psi)*cabs(Psi)) - dtCoEff2/hbar*(0.5*m*smallOmegaSq*(x*x + y*y) + g*cabs(oldPsi)*cabs(oldPsi));
              //fprintf(stderr,"%e, %e, oldField: %e\n", cabs(a[j]), cabs(d[j]), cabs(oldField[j][k]));
            }
        //exit(0);
          //Free Boundary Conditions
          b[1] = b[1] + 2.0*a[1];
          c[1] = c[1] - a[1];
          a[xExtent-1] = a[xExtent-1] - c[xExtent-1];
          b[xExtent-1] = b[xExtent-1] + 2.0*c[xExtent-1];
          solveTridiagonal2(a,b,c,X,d);
          for(j = 1; j < xExtent; j ++)
          {
            currentField[j][k] = X[j];
            //if(p == 3) fprintf(stderr,"%e %e \n", cabs(X[j]), carg(X[j]));
          }
          //if(p == 3) exit(0);
        }
      for(j = 1; j < xExtent; j ++)
        {
          for(k = 1; k < yExtent; k ++)
            {
              //a[k] = -alpha_y;
              //b[k] = 1.0 + 2.0*alpha_y;
              //c[k] = -alpha_y;
              d[k] = currentField[j][k];
            }
          //Free Boundary Conditions
          b[1] = b[1] + 2.0*a[1]; //DO NOT USE CONTRACTED OPERATORS HERE (like =- etc.)
          c[1] = c[1] - a[1];
          a[yExtent-1] = a[yExtent-1] - c[yExtent-1];
          b[yExtent-1] = b[yExtent-1] + 2.0*c[yExtent-1];
          //Solve Tridiagonal
          solveTridiagonal2(a,b,c,X,d);
          for(k = 1; k < yExtent; k ++)
          {
            currentField[j][k] = X[k];
        //fprintf(stderr,"%e, %e \n", cabs(X[k]), carg(X[k]));
          }
          //exit(0);
        }
      //Boundary Conditions for Edges
      setBoundaryFree(currentField);
    }
}

void setBoundaryFree(fieldType field[][Size])
{
  int xExtent = Size-1, yExtent = Size-1, j, k;

  for(k = 1; k < xExtent; k ++)
    {
      field[k][0] = 2.0*field[k][1] - field[k][2];
      field[k][yExtent] = 2.0*field[k][yExtent-1] - field[k][yExtent-2];
    }
  for(j = 0; j <= yExtent; j ++)
    {
      field[0][j] = 2.0*field[1][j] - field[2][j];
      field[xExtent][j] = 2.0*field[xExtent-1][j] - field[xExtent-2][j];
    }
  field[0][0] = 2.0*field[1][1] - field[2][2];
  field[0][yExtent] = 2.0*field[1][yExtent-1] - field[2][yExtent-2];
  field[xExtent][0] = 2.0*field[xExtent-1][1] - field[xExtent-2][2];
  field[xExtent][yExtent] = 2.0*field[xExtent-1][yExtent-1] - field[xExtent-2][yExtent-2];
}

void setBoundaryVanish(fieldType field[][Size])
{
  int xExtent = Size-1, yExtent = Size-1, j, k;

  //Phase Gradient Zero at Boundary
  for(j = 1; j < xExtent; j ++)
    {
      field[j][0] = complexZero*cexp(I*carg(field[j][0]));
      field[j][yExtent] = complexZero*cexp(I*carg(field[j][yExtent]));
    }
  for(k = 0; k <= yExtent; k ++)
    {
      field[0][k] = complexZero*cexp(I*carg(field[0][k]));
      field[xExtent][k] = complexZero*cexp(I*carg(field[xExtent][k]));
    }
}

void solveTridiagonal(fieldType a[], fieldType b[], fieldType c[], fieldType unknowns[], fieldType knowns[])
{
  int n = Size-1, startElement = 1, j;
  fieldType tmp[n], bet;

  if(b[startElement] == complexZero)
    {
      fprintf(stderr,"Tridiagonal Error - First Element of Diagonal is Zero\n");
      exit(0);
    }

  unknowns[startElement] = knowns[startElement]/(bet = b[startElement]);
  for(j = startElement+1; j < n; j ++)
    {
      tmp[j] = c[j-1]/bet;
      bet = b[j] - a[j]*tmp[j];
      if(bet == complexZero)
        {
          fprintf(stderr,"Tridiagonal Error - Zero pivot found (check if Matrix A is Diagonally Dominant\n)");
          exit(0);
        }
      else
        unknowns[j] = (knowns[j] - a[j]*unknowns[j-1])/bet;
    }
  //Backsubstitution
  for(j = n-2; j >= startElement; j --)
    unknowns[j] = unknowns[j] - tmp[j+1]*unknowns[j+1];
}

void solveTridiagonal2(fieldType a[], fieldType b[], fieldType c[], fieldType unknowns[], fieldType knowns[])
{
    fieldType tmp[Size];
        int n = Size, startElement = 1, j;

        for(j = startElement+1; j < n; j++)
        {
                tmp[j] = a[j]/b[j-1];
                b[j] = b[j] - tmp[j]*c[j-1];
                knowns[j] = knowns[j] - tmp[j]*knowns[j-1];
        }
        unknowns[n-1] = knowns[n-1]/b[n-1];
        for(j = n-2; j >= startElement; j--)
                unknowns[j] = (knowns[j] - c[j]*unknowns[j+1])/b[j];
}

void writeImage(fieldType field[][Size])
{
  FILE *outFile;
  int j,k;

  outFile = fopen("data.dat","w");

  for(k = 0; k < Size; k ++)
    {
      for(j = 0; j < Size; j ++)
        fprintf(outFile,"%e , ", field[j][k]);
      fprintf(outFile,"\n");
    }

  fclose(outFile);
}

void writeRadialSolution(parameterType F[], parameterType dr)
{
  FILE *outFile;
  int j;

  outFile = fopen("cradial.csv","w");

  //fprintf(outFile,"%d\n", Size);
  for(j = 0; j < Size; j ++)
    {
      fprintf(outFile,"%f, %e", j*dr, F[j]);
      fprintf(outFile,"\n");
    }

  fclose(outFile);
}

void writeResultsCSV(fieldType field[][Size])
{
  //Writes Output of the Norm and the Phase of the Result in Excel Comma Separated Values (CSV) Format.
  FILE *outFile;
  int j,k;

  outFile = fopen("norm.csv","w");

  for(k = 0; k < Size; k ++)
    {
      for(j = 0; j < Size; j ++)
        fprintf(outFile,"%e , ",cabs(field[j][k]));
      fprintf(outFile,"\n");
    }

  fclose(outFile);
  outFile = fopen("phase.csv","w");

  for(k = 0; k < Size; k ++)
    {
      for(j = 0; j < Size; j ++)
        fprintf(outFile,"%e , ",carg(field[j][k]));
      fprintf(outFile,"\n");
    }

  fclose(outFile);
}

void writeResultsGNU(fieldType field[][Size])
{
  //Writes Output of the Norm and the Phase of the Result in GNU Plot Format.
  FILE *outFile;
  int j,k;

  outFile = fopen("norm.dat","w");

  //fprintf(outFile,"%d\n", Size*Size);
  for(k = 0; k < Size; k ++)
      for(j = 0; j < Size; j ++)
      {
        fprintf(outFile,"%d\t%d\t%e", j, k, cabs(field[j][k]));
        fprintf(outFile,"\n");
      }

  fclose(outFile);
  outFile = fopen("phase.dat","w");

  //fprintf(outFile,"%d\n", Size*Size);
  for(k = 0; k < Size; k ++)
      for(j = 0; j < Size; j ++)
      {
        fprintf(outFile,"%d\t%d\t%e", j, k, carg(field[j][k]));
        fprintf(outFile,"\n");
      }

  fclose(outFile);
}

---