Source Jam - aeneas-1.1

00036 {
00037  int i,j;
00038  double *v;
00039  double *nabla2R;
00040  double *R;
00041  double *ac[4],*bc[4],*cc[4],*dc[4];
00042  double *xi[3], l[4];
00043  double dx,dy,dz;
00044  
00045 // infinitesimal increments = 5 Angstrom
00046  dx=0.5e-9; dy=0.5e-9; dz=0.5e-9;
00047  
00048  v=malloc((Ng+1)*sizeof(double));
00049  nabla2R=malloc((Ng+1)*sizeof(double));
00050  R=malloc((Ng+1)*sizeof(double));
00051  for(i=0;i<3;i++) xi[i]=malloc((Ne+1)*sizeof(double));
00052  for(i=0;i<4;i++){
00053    ac[i]=malloc((Ne+1)*sizeof(double));
00054    bc[i]=malloc((Ne+1)*sizeof(double));
00055    cc[i]=malloc((Ne+1)*sizeof(double));
00056    dc[i]=malloc((Ne+1)*sizeof(double));
00057  }
00058 
00059 // barycentric coordinates of the Guass point
00060  l[0]=0.25;  l[1]=0.25;  l[2]=0.25;  l[3]=0.25;
00061 
00062 // calculation of the center of mass coordinates for every element
00063  for(i=0;i<Ne;i++){
00064       xi[0][i]=l[0]*coord[0][noeud_geo[0][i]-1]
00065               +l[1]*coord[0][noeud_geo[1][i]-1]
00066               +l[2]*coord[0][noeud_geo[2][i]-1]
00067               +l[3]*coord[0][noeud_geo[3][i]-1];
00068       xi[1][i]=l[0]*coord[1][noeud_geo[0][i]-1]
00069               +l[1]*coord[1][noeud_geo[1][i]-1]
00070               +l[2]*coord[1][noeud_geo[2][i]-1]
00071               +l[3]*coord[1][noeud_geo[3][i]-1];
00072       xi[2][i]=l[0]*coord[2][noeud_geo[0][i]-1]
00073               +l[1]*coord[2][noeud_geo[1][i]-1]
00074               +l[2]*coord[2][noeud_geo[2][i]-1]
00075               +l[3]*coord[2][noeud_geo[3][i]-1];
00076 //    printf("%g %g %g\n",xi[0][j][i],xi[1][j][i],xi[2][j][i]);
00077     }
00078 
00079 // calculation of the theta functions for every element
00080  for(i=0;i<Ne;i++){
00081     int l;
00082     int n=4;
00083     int indx[5];
00084     double d;
00085     double a[5][5];
00086     double b[5];
00087 
00088 // definition of the a matrix for the linear system
00089     {int k;
00090      for(k=0;k<4;k++){
00091        a[k+1][1]=coord[0][noeud_geo[k][i]-1];
00092        a[k+1][2]=coord[1][noeud_geo[k][i]-1];
00093        a[k+1][3]=coord[2][noeud_geo[k][i]-1];
00094        a[k+1][4]=1.;
00095      }}
00096 // LU decomposition of the A matrix
00097 #include "../linear_solver/ludcmp.h"
00098     for(l=0;l<4;l++){
00099       if(l==0){b[1]=1.;b[2]=0.;b[3]=0.;b[4]=0.;}
00100       if(l==1){b[1]=0.;b[2]=1.;b[3]=0.;b[4]=0.;}
00101       if(l==2){b[1]=0.;b[2]=0.;b[3]=1.;b[4]=0.;}
00102       if(l==3){b[1]=0.;b[2]=0.;b[3]=0.;b[4]=1.;}
00103 #include "../linear_solver/lubksb.h"
00104       ac[l][i]=b[1]; bc[l][i]=b[2]; cc[l][i]=b[3]; dc[l][i]=b[4];
00105 // printf("%g %g %g %g\n",b[1],b[2],b[3],b[4]);
00106     } // end of l-cycle
00107 
00108 // definition of R field on each node
00109     for(j=0;j<4;j++) R[noeud_geo[j][i]-1]=sqrt(NE[noeud_geo[j][i]-1]);
00110  } // end of i-cycle
00111 
00112 // just in case...
00113  for(i=0;i<Ng;i++) nabla2R[i]=0.;
00114 
00115 // definition of laplacian of R on each node
00116  for(i=0;i<Ne;i++){
00117    double h1,h2,h3;
00118 // second order finite difference scheme
00119    h1=R[noeud_geo[0][i]-1]*(ac[0][i]*(xi[0][i]+dx)+bc[0][i]*xi[1][i]+cc[0][i]*xi[2][i]+dc[0][i])
00120      +R[noeud_geo[1][i]-1]*(ac[1][i]*(xi[0][i]+dx)+bc[1][i]*xi[1][i]+cc[1][i]*xi[2][i]+dc[1][i])
00121      +R[noeud_geo[2][i]-1]*(ac[2][i]*(xi[0][i]+dx)+bc[2][i]*xi[1][i]+cc[2][i]*xi[2][i]+dc[2][i])
00122      +R[noeud_geo[3][i]-1]*(ac[3][i]*(xi[0][i]+dx)+bc[3][i]*xi[1][i]+cc[3][i]*xi[2][i]+dc[3][i]);
00123    h2=R[noeud_geo[0][i]-1]*(ac[0][i]*xi[0][i]+bc[0][i]*xi[1][i]+cc[0][i]*xi[2][i]+dc[0][i])
00124      +R[noeud_geo[1][i]-1]*(ac[1][i]*xi[0][i]+bc[1][i]*xi[1][i]+cc[1][i]*xi[2][i]+dc[1][i])
00125      +R[noeud_geo[2][i]-1]*(ac[2][i]*xi[0][i]+bc[2][i]*xi[1][i]+cc[2][i]*xi[2][i]+dc[2][i])
00126      +R[noeud_geo[3][i]-1]*(ac[3][i]*xi[0][i]+bc[3][i]*xi[1][i]+cc[3][i]*xi[2][i]+dc[3][i]);
00127    h3=R[noeud_geo[0][i]-1]*(ac[0][i]*(xi[0][i]-dx)+bc[0][i]*xi[1][i]+cc[0][i]*xi[2][i]+dc[0][i])
00128      +R[noeud_geo[1][i]-1]*(ac[1][i]*(xi[0][i]-dx)+bc[1][i]*xi[1][i]+cc[1][i]*xi[2][i]+dc[1][i])
00129      +R[noeud_geo[2][i]-1]*(ac[2][i]*(xi[0][i]-dx)+bc[2][i]*xi[1][i]+cc[2][i]*xi[2][i]+dc[2][i])
00130      +R[noeud_geo[3][i]-1]*(ac[3][i]*(xi[0][i]-dx)+bc[3][i]*xi[1][i]+cc[3][i]*xi[2][i]+dc[3][i]);
00131    for(j=0;j<4;j++) nabla2R[noeud_geo[j][i]-1]+=(h1-2.*h2+h3)/(dx*dx);
00132    h1=R[noeud_geo[0][i]-1]*(ac[0][i]*xi[0][i]+bc[0][i]*(xi[1][i]+dy)+cc[0][i]*xi[2][i]+dc[0][i])
00133      +R[noeud_geo[1][i]-1]*(ac[1][i]*xi[0][i]+bc[1][i]*(xi[1][i]+dy)+cc[1][i]*xi[2][i]+dc[1][i])
00134      +R[noeud_geo[2][i]-1]*(ac[2][i]*xi[0][i]+bc[2][i]*(xi[1][i]+dy)+cc[2][i]*xi[2][i]+dc[2][i])
00135      +R[noeud_geo[3][i]-1]*(ac[3][i]*xi[0][i]+bc[3][i]*(xi[1][i]+dy)+cc[3][i]*xi[2][i]+dc[3][i]);
00136    h2=R[noeud_geo[0][i]-1]*(ac[0][i]*xi[0][i]+bc[0][i]*xi[1][i]+cc[0][i]*xi[2][i]+dc[0][i])
00137      +R[noeud_geo[1][i]-1]*(ac[1][i]*xi[0][i]+bc[1][i]*xi[1][i]+cc[1][i]*xi[2][i]+dc[1][i])
00138      +R[noeud_geo[2][i]-1]*(ac[2][i]*xi[0][i]+bc[2][i]*xi[1][i]+cc[2][i]*xi[2][i]+dc[2][i])
00139      +R[noeud_geo[3][i]-1]*(ac[3][i]*xi[0][i]+bc[3][i]*xi[1][i]+cc[3][i]*xi[2][i]+dc[3][i]);
00140    h3=R[noeud_geo[0][i]-1]*(ac[0][i]*xi[0][i]+bc[0][i]*(xi[1][i]-dy)+cc[0][i]*xi[2][i]+dc[0][i])
00141      +R[noeud_geo[1][i]-1]*(ac[1][i]*xi[0][i]+bc[1][i]*(xi[1][i]-dy)+cc[1][i]*xi[2][i]+dc[1][i])
00142      +R[noeud_geo[2][i]-1]*(ac[2][i]*xi[0][i]+bc[2][i]*(xi[1][i]-dy)+cc[2][i]*xi[2][i]+dc[2][i])
00143      +R[noeud_geo[3][i]-1]*(ac[3][i]*xi[0][i]+bc[3][i]*(xi[1][i]-dy)+cc[3][i]*xi[2][i]+dc[3][i]);
00144    for(j=0;j<4;j++) nabla2R[noeud_geo[j][i]-1]+=(h1-2.*h2+h3)/(dy*dy);
00145    h1=R[noeud_geo[0][i]-1]*(ac[0][i]*xi[0][i]+bc[0][i]*xi[1][i]+cc[0][i]*(xi[2][i]+dz)+dc[0][i])
00146      +R[noeud_geo[1][i]-1]*(ac[1][i]*xi[0][i]+bc[1][i]*xi[1][i]+cc[1][i]*(xi[2][i]+dz)+dc[1][i])
00147      +R[noeud_geo[2][i]-1]*(ac[2][i]*xi[0][i]+bc[2][i]*xi[1][i]+cc[2][i]*(xi[2][i]+dz)+dc[2][i])
00148      +R[noeud_geo[3][i]-1]*(ac[3][i]*xi[0][i]+bc[3][i]*xi[1][i]+cc[3][i]*(xi[2][i]+dz)+dc[3][i]);
00149    h2=R[noeud_geo[0][i]-1]*(ac[0][i]*xi[0][i]+bc[0][i]*xi[1][i]+cc[0][i]*xi[2][i]+dc[0][i])
00150      +R[noeud_geo[1][i]-1]*(ac[1][i]*xi[0][i]+bc[1][i]*xi[1][i]+cc[1][i]*xi[2][i]+dc[1][i])
00151      +R[noeud_geo[2][i]-1]*(ac[2][i]*xi[0][i]+bc[2][i]*xi[1][i]+cc[2][i]*xi[2][i]+dc[2][i])
00152      +R[noeud_geo[3][i]-1]*(ac[3][i]*xi[0][i]+bc[3][i]*xi[1][i]+cc[3][i]*xi[2][i]+dc[3][i]);
00153    h3=R[noeud_geo[0][i]-1]*(ac[0][i]*xi[0][i]+bc[0][i]*xi[1][i]+cc[0][i]*(xi[2][i]-dz)+dc[0][i])
00154      +R[noeud_geo[1][i]-1]*(ac[1][i]*xi[0][i]+bc[1][i]*xi[1][i]+cc[1][i]*(xi[2][i]-dz)+dc[1][i])
00155      +R[noeud_geo[2][i]-1]*(ac[2][i]*xi[0][i]+bc[2][i]*xi[1][i]+cc[2][i]*(xi[2][i]-dz)+dc[2][i])
00156      +R[noeud_geo[3][i]-1]*(ac[3][i]*xi[0][i]+bc[3][i]*xi[1][i]+cc[3][i]*(xi[2][i]-dz)+dc[3][i]);
00157    for(j=0;j<4;j++) nabla2R[noeud_geo[j][i]-1]+=(h1-2.*h2+h3)/(dz*dz);
00158  }
00159 
00160  for(i=0;i<Ng;i++) v[i]=0.; // just in case...
00161  
00162 // updates the potential by means of Bohm potential
00163  for(i=0;i<Ne;i++){
00164   for(j=0;j<4;j++)
00165    if(R[noeud_geo[j][i]-1]!=0.)
00166      V[noeud_geo[j][i]-1]+=-(HBAR*HBAR/(2.*MSTAR[i_dom[i]][1]*M*Q))
00167                          *nabla2R[noeud_geo[j][i]-1]/R[noeud_geo[j][i]-1];
00168 // printf("QV = %g\n",-(HBAR*HBAR/(2.*MSTAR[i_dom[i]][1]*M*Q))*nabla2R[noeud_geo[0][i]-1]/R[noeud_geo[0][i]-1]);
00169 //printf("QV = %g\n",V[noeud_geo[0][i]-1]);
00170  }
00171  
00172 // free precedently allocated memory
00173  free(v);
00174  free(R);
00175  free(nabla2R);
00176  for(i=0;i<3;i++) free(xi[i]);
00177  for(i=0;i<4;i++){
00178    free(ac[i]);
00179    free(bc[i]);
00180    free(cc[i]);
00181    free(dc[i]);
00182  }
00183 }
Bohm_potential.h File Reference

Functions

Function Documentation


Functions
void	Bohm_potential (void)