c_code/main.c

00001 //Compile as gcc -o bec main.c -lm
00002 //Execute as ./bec in Linux
00003 
00004 //--------------Preamble------------
00005 #include <stdlib.h>
00006 #include <stdio.h>
00007 #include <complex.h>
00008 #include <time.h>
00009 
00010 //Datatype Aliases
00011 #define TRUE (0 == 0)
00012 #define FALSE (1 == 0)
00013 typedef double complex fieldType;
00014 typedef double parameterType;
00015 
00016 //Global Constants
00017 const int Size = 180;
00018 const int maxTimesteps = 1000;
00019 const int shotImageAt = 500;
00020 const fieldType xy0 = 20;
00021 const fieldType dt = 1e-3;
00022 const parameterType tolerance = 1e-10;
00023 const parameterType windingNumber = 1.0;
00024 const fieldType complexZero = 0.0;
00025 
00026 //Global Parameters
00027 const fieldType g = 1.0;
00028 const fieldType hbar = 1.0;
00029 const fieldType smallOmegaSq = 0.01;
00030 const parameterType mu = 1.0;
00031 const fieldType m = 1.0;
00032 const fieldType Omega = 0.0;
00033 const parameterType phaseStepWidth = 10;
00034 const parameterType phaseStepHeight = 2.5;
00035 const parameterType phaseStepPosition = 80;
00036 
00037 //Function Prototypes
00038 void init(fieldType field[][Size], fieldType value);
00039 void init1D(parameterType F[], parameterType value);
00040 void copy(fieldType source[][Size], fieldType dest[][Size]);
00041 void loadParameters(const char *filename); //Not Functional ATM
00042 void steadyState(fieldType field[][Size], parameterType F[], parameterType dr);
00043 void steadyVortexState(fieldType field[][Size], parameterType F[], parameterType dr);
00044 void transferRadial(parameterType F[], fieldType field[][Size], parameterType xOffset, parameterType yOffset, parameterType dx, parameterType dy, parameterType dr);
00045 void setupUniformPhase(fieldType field[][Size], parameterType phase, parameterType xOffset, parameterType yOffset, parameterType dx, parameterType dy);
00046 void setupVortexPhase(fieldType field[][Size], parameterType xOffset, parameterType yOffset, parameterType dx, parameterType dy);
00047 void applyPhaseStep(fieldType field[][Size]);
00048 void evolveUsingTDGP(fieldType field[][Size], fieldType oldfield[][Size], fieldType xOffset, fieldType yOffset, fieldType h_t, fieldType h_x, fieldType h_y);
00049 void setBoundaryFree(fieldType field[][Size]);
00050 void setBoundaryVanish(fieldType field[][Size]);
00051 void solveTridiagonal(fieldType a[], fieldType b[], fieldType c[], fieldType unknowns[], fieldType knowns[]);
00052 void solveTridiagonal2(fieldType a[], fieldType b[], fieldType c[], fieldType unknowns[], fieldType knowns[]);
00053 void writeImage(fieldType field[][Size]); //Not Functional ATM
00054 void writeRadialSolution(parameterType F[], parameterType dr);
00055 void writeResultsCSV(fieldType field[][Size]);
00056 void writeResultsGNU(fieldType field[][Size]);
00057 //----------------------------------
00058 
00059 //======================Main=========================
00060 int main()
00061 {
00062   fieldType Field[Size][Size], OldField[Size][Size];
00063   parameterType radialSolution[Size];
00064   parameterType spatialStep = 2.0*(double)xy0/(double)Size, timeDiff = 0.0;
00065   parameterType g, mu, smallOmegaSq;
00066   int p = 0;
00067   time_t start,end;
00068 
00069   //Initialize.
00070   time(&start);
00071   init(Field,0.0);
00072   init(OldField,0.0);
00073   init1D(radialSolution,1.0);
00074 
00075   /*for(p = 0; p < Size; p ++)
00076     fprintf(stderr,"%e, ",radialSolution[p]);
00077     fprintf(stderr,"\n");*/
00078 
00079   printf("---------------------------------------------------------\n");
00080   printf("||\t BEC Simulation\n");
00081 
00082   //Calculate and Setup Steady State.
00083   //steadyState(Field,radialSolution,spatialStep);
00084   steadyVortexState(Field,radialSolution,spatialStep);
00085   transferRadial(radialSolution,Field,xy0,xy0,spatialStep,spatialStep,spatialStep);
00086 
00087   //Setup Initial Phase
00088   //setupUniformPhase(Field,1.0,xy0,xy0,spatialStep,spatialStep);
00089   setupVortexPhase(Field,xy0,xy0,spatialStep,spatialStep);
00090   applyPhaseStep(Field);
00091 
00092   //Evolve Field according TDGP Equation using parameters provided.
00093   fprintf(stderr,"||\t Evolving... ");
00094   for(p = 0; p < maxTimesteps; p ++)
00095     evolveUsingTDGP(Field,OldField,xy0,xy0,dt,spatialStep,spatialStep);
00096   fprintf(stderr,"Done\n");
00097 
00098   //Write Results of Simulation.
00099   writeResultsCSV(Field);
00100   writeRadialSolution(radialSolution,spatialStep);
00101 
00102   //Calculate Time Taken and Close Application.
00103   time(&end);
00104   timeDiff = difftime(end,start);
00105   printf("||\t The Process took %f minutes (%f hours) to Complete.\n", (double)(timeDiff/60), (double)(timeDiff/3600));
00106   printf("---------------------------------------------------------\n");
00107   return 0;
00108 }
00109 //===================================================
00110 
00111 //Function Bodies
00112 void init(fieldType field[][Size], fieldType value)
00113 {
00114   int k,j;
00115 
00116   for(k = 0; k < Size; k ++)
00117     for(j = 0; j < Size; j ++)
00118       field[j][k] = value;
00119 }
00120 
00121 void init1D(parameterType F[], parameterType value)
00122 {
00123   int j;
00124 
00125   for(j = 0; j < Size; j ++)
00126     F[j] = value;
00127 }
00128 
00129 void copy(fieldType source[][Size], fieldType dest[][Size])
00130 {
00131   int j, k;
00132 
00133   for(k = 0; k < Size; k++)
00134     for(j = 0; j < Size; j++)
00135       dest[j][k] = source[j][k];
00136 }
00137 
00138 void loadParameters(const char *filename)
00139 {
00140   FILE *outFile;
00141   int j,k;
00142 
00143   outFile = fopen(filename,"r");
00144 
00145   fclose(outFile);
00146 }
00147 
00148 void steadyState(fieldType field[][Size], parameterType F[], parameterType dr)
00149 {
00150   parameterType sum = 0.0, oldSum = 1.0, r = 0.0;
00151   int j, radialExtent = Size-1;
00152 
00153   while(fabs(sum - oldSum) > tolerance && radialExtent > 0)
00154   {
00155     oldSum = sum;
00156     sum = 0.0;
00157     for(j = 1; j < radialExtent; j ++)
00158       {
00159         r = j*dr;
00160         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);
00161         sum += F[j];
00162       }
00163     F[radialExtent] = 0.0;
00164     F[0] = 2*F[1] - F[1];
00165   }
00166 }
00167 
00168 void steadyVortexState(fieldType field[][Size], parameterType F[], parameterType dr)
00169 {
00170   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;
00171   int j, radialExtent = Size-1, count = 0;
00172 
00173   if(radialExtent > 0)
00174     F[0] = 0.0;
00175 
00176   while(fabs(sum - oldSum) > tolerance && radialExtent > 0)
00177   {
00178     oldSum = sum;
00179     sum = 0.0;
00180     for(j = 1; j < radialExtent; j ++)
00181       {
00182         r = j*dr;
00183         internal = g*F[j]*F[j] + 0.5*smallOmegaSq*r*r - mu + n*Omega;
00184         jFrac = (F[j+1]-F[j-1])/j;
00185         frac1 = 0.5*dr*jFrac;
00186         jSq = (parameterType)j*j;
00187         frac2 = (n*n)/(jSq)*F[j];
00188         extra = F[j+1]+F[j-1];
00189         F[j] = 0.5*( extra+frac1-frac2 ) - dr*dr*F[j]*internal;
00190         sum += F[j];
00191         //fprintf(stderr,"%e \n", F[j]);
00192       }
00193     //exit(0);
00194     //fprintf(stderr,"%d ", count);
00195     //count ++;
00196     //F[radialExtent] = F[radialExtent-1];
00197     F[radialExtent] = 0.0;
00198   }
00199   //fprintf(stderr,"Tolerance Met!\n");
00200 }
00201 
00202 void transferRadial(parameterType F[], fieldType field[][Size], parameterType xOffset, parameterType yOffset, parameterType dx, parameterType dy, parameterType dr)
00203 {
00204   parameterType x, y;
00205   parameterType l, rSquared;
00206   parameterType magnitude;
00207   int j, k, rInteger;
00208 
00209   for(k = 0; k < Size; k++)
00210     {
00211       y = (k)*dy - yOffset;
00212       for(j = 0; j < Size; j++)
00213         {
00214           x = (j)*dx - xOffset;
00215           rSquared = x*x + y*y;
00216           l = sqrt(rSquared)/dr;
00217           rInteger = (int)(l);
00218           magnitude = F[rInteger] + (l-rInteger)*( F[rInteger+1] - F[rInteger] );
00219           field[j][k] = magnitude;
00220         }
00221     }
00222 }
00223 
00224 void setupUniformPhase(fieldType field[][Size], parameterType phase, parameterType xOffset, parameterType yOffset, parameterType dx, parameterType dy)
00225 {
00226   parameterType x, y;
00227   int j, k;
00228 
00229   for(k = 0; k < Size; k++)
00230     {
00231       y = k*dy - yOffset;
00232       for(j = 0; j < Size; j++)
00233         {
00234           x = j*dx - xOffset;
00235           field[j][k] = cabs(field[j][k])*cexp(I*phase);
00236         }
00237     }
00238 }
00239 
00240 void setupVortexPhase(fieldType field[][Size], parameterType xOffset, parameterType yOffset, parameterType dx, parameterType dy)
00241 {
00242   parameterType x, y;
00243   int j, k;
00244 
00245   for(k = 0; k < Size; k++)
00246     {
00247       y = k*dy - yOffset;
00248       for(j = 0; j < Size; j++)
00249         {
00250           x = j*dx - xOffset;
00251           field[j][k] = cabs(field[j][k])*cexp(I*atan2(y,x));
00252         }
00253     }
00254 }
00255 
00256 void applyPhaseStep(fieldType field[][Size])
00257 {
00258   int j, k;
00259   parameterType endCounter = (double)(phaseStepWidth), counter = 1.0, gradientFactor = 0.0;
00260 
00261   for(k = phaseStepPosition; k < Size; k ++)
00262     {
00263       for(j = 0; j < Size; j ++)
00264         {
00265           gradientFactor = counter/endCounter;
00266           field[j][k] = field[j][k]*cexp(I*gradientFactor*phaseStepHeight);
00267           //field[j][k] = cabs(field[j][k])*cexp(I*carg(field[j][k]))*cexp(I*gradientFactor*phaseStepHeight);
00268         }
00269       if(counter < endCounter)
00270         counter ++;
00271     }
00272 }
00273 
00274 void evolveUsingTDGP(fieldType currentField[][Size], fieldType oldField[][Size], fieldType xOffset, fieldType yOffset, fieldType h_t, fieldType h_x, fieldType h_y)
00275 {
00276   int maxFSIIterations = 4, xExtent = Size-1, yExtent = Size-1, j, k, p;
00277   fieldType a[xExtent], b[xExtent], c[xExtent], d[xExtent], X[xExtent];
00278   fieldType a1[yExtent], b1[yExtent], c1[yExtent], d1[yExtent], X1[yExtent];
00279   fieldType x, y;
00280   x = 4.0*h_x*h_x*m; //x is used as tmp storage
00281   y = 4.0*h_y*h_y*m; //y is used as tmp storage
00282   fieldType     alpha_x = I*((h_t*hbar)/x),
00283     alpha_y = I*((h_t*hbar)/y),
00284     Psi, oldPsi;
00285 
00286     fieldType dtCoEff1, dtCoEff2;
00287 
00288   copy(currentField,oldField);
00289 
00290   for(p = 0; p < maxFSIIterations; p ++)
00291     {
00292       for(k = 1; k < yExtent; k ++)
00293         {
00294           y = (k)*h_y - yOffset;
00295           for(j = 1; j < xExtent; j ++)
00296             {
00297               x = (j)*h_x - xOffset;
00298               Psi = currentField[j][k];
00299               oldPsi = oldField[j][k];
00300 
00301               dtCoEff1 = h_t*0.5*Psi*I;
00302               dtCoEff2 = h_t*0.5*oldPsi*I;
00303 
00304               a[j] = -alpha_x;
00305               b[j] = 1.0 + 2.0*alpha_x;
00306               c[j] = -alpha_x;
00307               //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));
00308               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));
00309               //fprintf(stderr,"%e, %e, oldField: %e\n", cabs(a[j]), cabs(d[j]), cabs(oldField[j][k]));
00310             }
00311         //exit(0);
00312           //Free Boundary Conditions
00313           b[1] = b[1] + 2.0*a[1];
00314           c[1] = c[1] - a[1];
00315           a[xExtent-1] = a[xExtent-1] - c[xExtent-1];
00316           b[xExtent-1] = b[xExtent-1] + 2.0*c[xExtent-1];
00317           solveTridiagonal2(a,b,c,X,d);
00318           for(j = 1; j < xExtent; j ++)
00319           {
00320             currentField[j][k] = X[j];
00321             //if(p == 3) fprintf(stderr,"%e %e \n", cabs(X[j]), carg(X[j]));
00322           }
00323           //if(p == 3) exit(0);
00324         }
00325       for(j = 1; j < xExtent; j ++)
00326         {
00327           for(k = 1; k < yExtent; k ++)
00328             {
00329               //a[k] = -alpha_y;
00330               //b[k] = 1.0 + 2.0*alpha_y;
00331               //c[k] = -alpha_y;
00332               d[k] = currentField[j][k];
00333             }
00334           //Free Boundary Conditions
00335           b[1] = b[1] + 2.0*a[1]; //DO NOT USE CONTRACTED OPERATORS HERE (like =- etc.)
00336           c[1] = c[1] - a[1];
00337           a[yExtent-1] = a[yExtent-1] - c[yExtent-1];
00338           b[yExtent-1] = b[yExtent-1] + 2.0*c[yExtent-1];
00339           //Solve Tridiagonal
00340           solveTridiagonal2(a,b,c,X,d);
00341           for(k = 1; k < yExtent; k ++)
00342           {
00343             currentField[j][k] = X[k];
00344         //fprintf(stderr,"%e, %e \n", cabs(X[k]), carg(X[k]));
00345           }
00346           //exit(0);
00347         }
00348       //Boundary Conditions for Edges
00349       setBoundaryFree(currentField);
00350     }
00351 }
00352 
00353 void setBoundaryFree(fieldType field[][Size])
00354 {
00355   int xExtent = Size-1, yExtent = Size-1, j, k;
00356 
00357   for(k = 1; k < xExtent; k ++)
00358     {
00359       field[k][0] = 2.0*field[k][1] - field[k][2];
00360       field[k][yExtent] = 2.0*field[k][yExtent-1] - field[k][yExtent-2];
00361     }
00362   for(j = 0; j <= yExtent; j ++)
00363     {
00364       field[0][j] = 2.0*field[1][j] - field[2][j];
00365       field[xExtent][j] = 2.0*field[xExtent-1][j] - field[xExtent-2][j];
00366     }
00367   field[0][0] = 2.0*field[1][1] - field[2][2];
00368   field[0][yExtent] = 2.0*field[1][yExtent-1] - field[2][yExtent-2];
00369   field[xExtent][0] = 2.0*field[xExtent-1][1] - field[xExtent-2][2];
00370   field[xExtent][yExtent] = 2.0*field[xExtent-1][yExtent-1] - field[xExtent-2][yExtent-2];
00371 }
00372 
00373 void setBoundaryVanish(fieldType field[][Size])
00374 {
00375   int xExtent = Size-1, yExtent = Size-1, j, k;
00376 
00377   //Phase Gradient Zero at Boundary
00378   for(j = 1; j < xExtent; j ++)
00379     {
00380       field[j][0] = complexZero*cexp(I*carg(field[j][0]));
00381       field[j][yExtent] = complexZero*cexp(I*carg(field[j][yExtent]));
00382     }
00383   for(k = 0; k <= yExtent; k ++)
00384     {
00385       field[0][k] = complexZero*cexp(I*carg(field[0][k]));
00386       field[xExtent][k] = complexZero*cexp(I*carg(field[xExtent][k]));
00387     }
00388 }
00389 
00390 void solveTridiagonal(fieldType a[], fieldType b[], fieldType c[], fieldType unknowns[], fieldType knowns[])
00391 {
00392   int n = Size-1, startElement = 1, j;
00393   fieldType tmp[n], bet;
00394 
00395   if(b[startElement] == complexZero)
00396     {
00397       fprintf(stderr,"Tridiagonal Error - First Element of Diagonal is Zero\n");
00398       exit(0);
00399     }
00400 
00401   unknowns[startElement] = knowns[startElement]/(bet = b[startElement]);
00402   for(j = startElement+1; j < n; j ++)
00403     {
00404       tmp[j] = c[j-1]/bet;
00405       bet = b[j] - a[j]*tmp[j];
00406       if(bet == complexZero)
00407         {
00408           fprintf(stderr,"Tridiagonal Error - Zero pivot found (check if Matrix A is Diagonally Dominant\n)");
00409           exit(0);
00410         }
00411       else
00412         unknowns[j] = (knowns[j] - a[j]*unknowns[j-1])/bet;
00413     }
00414   //Backsubstitution
00415   for(j = n-2; j >= startElement; j --)
00416     unknowns[j] = unknowns[j] - tmp[j+1]*unknowns[j+1];
00417 }
00418 
00419 void solveTridiagonal2(fieldType a[], fieldType b[], fieldType c[], fieldType unknowns[], fieldType knowns[])
00420 {
00421     fieldType tmp[Size];
00422         int n = Size, startElement = 1, j;
00423 
00424         for(j = startElement+1; j < n; j++)
00425         {
00426                 tmp[j] = a[j]/b[j-1];
00427                 b[j] = b[j] - tmp[j]*c[j-1];
00428                 knowns[j] = knowns[j] - tmp[j]*knowns[j-1];
00429         }
00430         unknowns[n-1] = knowns[n-1]/b[n-1];
00431         for(j = n-2; j >= startElement; j--)
00432                 unknowns[j] = (knowns[j] - c[j]*unknowns[j+1])/b[j];
00433 }
00434 
00435 void writeImage(fieldType field[][Size])
00436 {
00437   FILE *outFile;
00438   int j,k;
00439 
00440   outFile = fopen("data.dat","w");
00441 
00442   for(k = 0; k < Size; k ++)
00443     {
00444       for(j = 0; j < Size; j ++)
00445         fprintf(outFile,"%e , ", field[j][k]);
00446       fprintf(outFile,"\n");
00447     }
00448 
00449   fclose(outFile);
00450 }
00451 
00452 void writeRadialSolution(parameterType F[], parameterType dr)
00453 {
00454   FILE *outFile;
00455   int j;
00456 
00457   outFile = fopen("cradial.csv","w");
00458 
00459   //fprintf(outFile,"%d\n", Size);
00460   for(j = 0; j < Size; j ++)
00461     {
00462       fprintf(outFile,"%f, %e", j*dr, F[j]);
00463       fprintf(outFile,"\n");
00464     }
00465 
00466   fclose(outFile);
00467 }
00468 
00469 void writeResultsCSV(fieldType field[][Size])
00470 {
00471   //Writes Output of the Norm and the Phase of the Result in Excel Comma Separated Values (CSV) Format.
00472   FILE *outFile;
00473   int j,k;
00474 
00475   outFile = fopen("norm.csv","w");
00476 
00477   for(k = 0; k < Size; k ++)
00478     {
00479       for(j = 0; j < Size; j ++)
00480         fprintf(outFile,"%e , ",cabs(field[j][k]));
00481       fprintf(outFile,"\n");
00482     }
00483 
00484   fclose(outFile);
00485   outFile = fopen("phase.csv","w");
00486 
00487   for(k = 0; k < Size; k ++)
00488     {
00489       for(j = 0; j < Size; j ++)
00490         fprintf(outFile,"%e , ",carg(field[j][k]));
00491       fprintf(outFile,"\n");
00492     }
00493 
00494   fclose(outFile);
00495 }
00496 
00497 void writeResultsGNU(fieldType field[][Size])
00498 {
00499   //Writes Output of the Norm and the Phase of the Result in GNU Plot Format.
00500   FILE *outFile;
00501   int j,k;
00502 
00503   outFile = fopen("norm.dat","w");
00504 
00505   //fprintf(outFile,"%d\n", Size*Size);
00506   for(k = 0; k < Size; k ++)
00507       for(j = 0; j < Size; j ++)
00508       {
00509         fprintf(outFile,"%d\t%d\t%e", j, k, cabs(field[j][k]));
00510         fprintf(outFile,"\n");
00511       }
00512 
00513   fclose(outFile);
00514   outFile = fopen("phase.dat","w");
00515 
00516   //fprintf(outFile,"%d\n", Size*Size);
00517   for(k = 0; k < Size; k ++)
00518       for(j = 0; j < Size; j ++)
00519       {
00520         fprintf(outFile,"%d\t%d\t%e", j, k, carg(field[j][k]));
00521         fprintf(outFile,"\n");
00522       }
00523 
00524   fclose(outFile);
00525 }

Generated on Sat May 13 13:22:48 2006 for Quantum Construct (qC++) by  doxygen 1.4.6-NO