Source Jam - aeneas-1.1

00030 {
00031  register int j;
00032  int n=4;
00033  int p=0;
00034  int ifr[3]={0.,0.,0.};
00035  int indx[5];
00036  double d;
00037  double a[5][5],b[5];
00038  double dkx,dky,dkz;
00039  double hmt, x, y, z;
00040  double ksquared,gk,thesquareroot;
00041 
00042  if(IV==9) return;
00043 // if(EL==-1) return; // just in case...
00044  if(i_dom[EL]==SIO2) return;
00045 
00046 // Valid for every material
00047  dkx=-QH*0.25*(Ex[noeud_geo[0][EL]-1]+Ex[noeud_geo[1][EL]-1]
00048               +Ex[noeud_geo[2][EL]-1]+Ex[noeud_geo[3][EL]-1])*tau;
00049  dky=-QH*0.25*(Ey[noeud_geo[0][EL]-1]+Ey[noeud_geo[1][EL]-1]
00050               +Ey[noeud_geo[2][EL]-1]+Ey[noeud_geo[3][EL]-1])*tau;
00051  dkz=-QH*0.25*(Ez[noeud_geo[0][EL]-1]+Ez[noeud_geo[1][EL]-1]
00052               +Ez[noeud_geo[2][EL]-1]+Ez[noeud_geo[3][EL]-1])*tau;
00053 // printf("KX=%g KY=%g KZ=%g\n",KX,KY,KZ);
00054 // printf("dkx=%g dky=%g dkz=%g\n",dkx,dky,dkz);
00055  
00056  x=y=z=0.;
00057  for(j=0;j<4;j++){
00058   x+=L[j]*coord[0][noeud_geo[j][EL]-1];
00059   y+=L[j]*coord[1][noeud_geo[j][EL]-1];
00060   z+=L[j]*coord[2][noeud_geo[j][EL]-1];
00061  }
00062  x/=4.; y/=4.; z/=4.;
00063 //if(z<-1.){
00064 //   printf("%g %g %g %g %g\n",L[0],L[1],L[2],L[3],L[0]+L[1]+L[2]+L[3]); 
00065 //   for(j=0;j<4;j++) printf("%d %g %g %g\n",j,coord[0][noeud_geo[j][EL]-1],coord[1][noeud_geo[j][EL]-1],coord[2][noeud_geo[j][EL]-1]);
00066 //   printf("%d %d x = %g y = %g z = %g\n",EL,noeud_geo[3][EL]-1,x,y,z);
00067 //}
00068 
00069 // x*=1.e-6; y*=1.e-6; z*=1.e-6;
00070   
00071  if(i_dom[EL]!=SIO2){
00072 // Electron drift process
00073    hmt=HM[i_dom[EL]][IV]*tau;
00074    ksquared=KX*KX+KY*KY+KZ*KZ;
00075    gk=HHM[i_dom[EL]][IV]*ksquared;
00076    thesquareroot=sqrt(1.+4.*alphaK[i_dom[EL]][IV]*gk);
00077 //   printf("!!! %g\n",thesquareroot);
00078  }
00079 // updating the position and pseudo-wave vector
00080 // (second order Runge-Kutta method)
00081 //printf("hmt = %g\n",hmt);
00082 // printf("X=%g Y=%g Z=%g\n",x,y,z);
00083 // printf("KX=%g KY=%g KZ=%g\n",KX,KY,KZ);
00084 // printf("dkx=%g dky=%g dkz=%g\n",dkx,dky,dkz);
00085 // printf("hmt = %g sq= %g\n",hmt,thesquareroot);
00086 // printf("%g %g %g\n",hmt*(KX+0.5*dkx)/thesquareroot,hmt*(KY+0.5*dky)/thesquareroot,
00087 //                     hmt*(KZ+0.5*dkz)/thesquareroot);
00088    x+=hmt*(KX+0.5*dkx)/thesquareroot;
00089    y+=hmt*(KY+0.5*dky)/thesquareroot;
00090    z+=hmt*(KZ+0.5*dkz)/thesquareroot;
00091    KX+=dkx;
00092    KY+=dky;
00093    KZ+=dkz;
00094 // printf("X=%g Y=%g Z=%g\n\n",x,y,z);
00095 
00096 // updating the barycentric coordinates
00097 //   x*=1.e6; y*=1.e6; z*=1.e6;
00098 //if(z>1.5) printf("x = %g y = %g z = %g\n",x,y,z);
00099 //   memset(&a,0,sizeof(a));
00100 //   memset(&b,0,sizeof(b));
00101    {int k;
00102     for(k=0;k<4;k++){
00103       a[1][k+1]=coord[0][noeud_geo[k][EL]-1];
00104       a[2][k+1]=coord[1][noeud_geo[k][EL]-1];
00105       a[3][k+1]=coord[2][noeud_geo[k][EL]-1];
00106       a[4][k+1]=1.;
00107    }}
00108    b[1]=4.*x;
00109    b[2]=4.*y;
00110    b[3]=4.*z;
00111    b[4]=4.;
00112 //printf("%g %g %g\n",coord[0][noeud_geo[1][EL]-1],coord[1][noeud_geo[1][EL]-1],
00113 //                    coord[2][noeud_geo[1][EL]-1]);
00114 //printf("%g %g %g\n",x,y,z);
00115 //printf("EL = %d\n",EL);
00116 //     {int ii,jj;
00117 //     for(ii=1;ii<=4;ii++){
00118 //        for(jj=1;jj<=4;jj++) printf("a[%d][%d] = %g   ",ii,jj,a[ii][jj]);
00119 //        printf("\n");
00120 //     }}
00121 // LU decomposition solver A.x=b
00122 #include "../linear_solver/ludcmp.h"
00123 #include "../linear_solver/lubksb.h"
00124 //   printf("before %g %g %g %g\n",L[0],L[1],L[2],L[3]);
00125    for(j=0;j<4;j++) L[j]=b[j+1];
00126 //   printf("** %g %g %g %g %g\n",L[0],L[1],L[2],L[3],L[0]+L[1]+L[2]+L[3]);
00127    if(L[0]<0. || L[1]<0. || L[2]<0. || L[3]<0.) goto label1;
00128    else return;
00129    label1:
00130 //    printf("EL = %d %g %g %g %g\n",EL,L[0],L[1],L[2],L[3]);
00131 //   printf("star %g %g %g %g %g\n",L[0],L[1],L[2],L[3],L[0]+L[1]+L[2]+L[3]);
00132     for(j=0;j<4;j++) if(L[j]<0.) p=j;
00133 // on the boundaries
00134 // =================
00135     if(p==0){
00136       ifr[0]=i_front[noeud_geo[1][EL]-1];
00137       ifr[1]=i_front[noeud_geo[2][EL]-1];
00138       ifr[2]=i_front[noeud_geo[3][EL]-1];
00139     }
00140     if(p==1){
00141       ifr[0]=i_front[noeud_geo[0][EL]-1];
00142       ifr[1]=i_front[noeud_geo[2][EL]-1];
00143       ifr[2]=i_front[noeud_geo[3][EL]-1];
00144     }
00145     if(p==2){
00146       ifr[0]=i_front[noeud_geo[0][EL]-1];
00147       ifr[1]=i_front[noeud_geo[1][EL]-1];
00148       ifr[2]=i_front[noeud_geo[3][EL]-1];
00149     }
00150     if(p==3){
00151       ifr[0]=i_front[noeud_geo[0][EL]-1];
00152       ifr[1]=i_front[noeud_geo[1][EL]-1];
00153       ifr[2]=i_front[noeud_geo[2][EL]-1];
00154     }
00155 // Schottky contact
00156     if(vois[p][EL]==-1 &&
00157        (ifr[0]==SCHOTTKY &&
00158         ifr[1]==SCHOTTKY &&
00159         ifr[2]==SCHOTTKY)){
00160             IV=9;
00161             return;
00162         }
00163 // ohmic contacts
00164     if(vois[p][EL]==-1 &&
00165        (ifr[0]==OHMIC &&
00166         ifr[1]==OHMIC &&
00167         ifr[2]==OHMIC)){
00168             IV=9;
00169             return;
00170         }
00171 // insulator contacts
00172     if(vois[p][EL]==-1 || i_dom[vois[p][EL]]==SIO2/* &&
00173        (ifr[0]==INSULATOR ||
00174         ifr[1]==INSULATOR ||
00175         ifr[2]==INSULATOR)*/){
00176         reflection:;
00177         double max=-4.;
00178         double x[2],y[2],z[2];
00179 // initial position
00180         x[0]=(L[0]*coord[0][noeud_geo[0][EL]-1]+
00181               L[1]*coord[0][noeud_geo[1][EL]-1]+
00182               L[2]*coord[0][noeud_geo[2][EL]-1]+
00183               L[3]*coord[0][noeud_geo[3][EL]-1])/(L[0]+L[1]+L[2]+L[3]);
00184         y[0]=(L[0]*coord[1][noeud_geo[0][EL]-1]+
00185               L[1]*coord[1][noeud_geo[1][EL]-1]+
00186               L[2]*coord[1][noeud_geo[2][EL]-1]+
00187               L[3]*coord[1][noeud_geo[3][EL]-1])/(L[0]+L[1]+L[2]+L[3]);
00188         z[0]=(L[0]*coord[2][noeud_geo[0][EL]-1]+
00189               L[1]*coord[2][noeud_geo[1][EL]-1]+
00190               L[2]*coord[2][noeud_geo[2][EL]-1]+
00191               L[3]*coord[2][noeud_geo[3][EL]-1])/(L[0]+L[1]+L[2]+L[3]);
00192         if(L[0]>L[1] && L[0]>L[2]) max=L[0];
00193         if(L[1]>L[0] && L[1]>L[2]) max=L[1];
00194         if(L[2]>L[1] && L[2]>L[0]) max=L[2];
00195         for(j=0;j<3;j++){
00196           L[j]/=max;
00197           if(L[j]<0.) L[j]=-L[j];
00198 //          L[j]=rnd();
00199         }
00200 // this is an elastic scattering
00201        L[3]=4.-(L[0]+L[1]+L[2]);
00202 // final position
00203         x[1]=(L[0]*coord[0][noeud_geo[0][EL]-1]+
00204               L[1]*coord[0][noeud_geo[1][EL]-1]+
00205               L[2]*coord[0][noeud_geo[2][EL]-1]+
00206               L[3]*coord[0][noeud_geo[3][EL]-1])/(L[0]+L[1]+L[2]+L[3]);
00207         y[1]=(L[0]*coord[1][noeud_geo[0][EL]-1]+
00208               L[1]*coord[1][noeud_geo[1][EL]-1]+
00209               L[2]*coord[1][noeud_geo[2][EL]-1]+
00210               L[3]*coord[1][noeud_geo[3][EL]-1])/(L[0]+L[1]+L[2]+L[3]);
00211         z[1]=(L[0]*coord[2][noeud_geo[0][EL]-1]+
00212               L[1]*coord[2][noeud_geo[1][EL]-1]+
00213               L[2]*coord[2][noeud_geo[2][EL]-1]+
00214               L[3]*coord[2][noeud_geo[3][EL]-1])/(L[0]+L[1]+L[2]+L[3]);
00215 // update the pseudo-wave vector
00216        KX=(x[1]-x[0])/tau;
00217        KY=(y[1]-y[0])/tau;
00218        KZ=(z[1]-z[0])/tau;
00219        return;
00220     }
00221 // scattering on SiO2
00222 /*    if(i_dom[vois[p][EL]]==SIO2){
00223         double max;
00224         if(L[0]>L[1] && L[0]>L[2]) max=L[0];
00225         if(L[1]>L[0] && L[1]>L[2]) max=L[1];
00226         if(L[2]>L[1] && L[2]>L[0]) max=L[2];
00227         for(j=0;j<3;j++){
00228           L[j]/=max;
00229           if(L[j]<0.) L[j]=-L[j];
00230         }
00231 // this is an elastic scattering
00232        L[3]=4.-(L[0]+L[1]+L[2]);
00233        return;
00234     }*/
00235 // in case of abrupt heterojunction, one has
00236 // to calculate the reflection and transmission
00237 // probability in order to see jow the particle
00238 // behaves.
00239    if(i_dom[EL]!=i_dom[vois[p][EL]] && i_dom[vois[p][EL]]!=SIO2
00240                                     && i_dom[vois[p][EL]]!=-1){
00241 //     printf("############### %d %d\n",EL,vois[p][EL]);
00242 //     printf("*************** %d %d\n\n",i_dom[EL],i_dom[vois[p][EL]]);
00243      double k1,k2,T,R;
00244      double ksquared,thesquareroot,energy;
00245 // calculates the particle energy
00246      ksquared=P[n][1]*P[n][1]+P[n][2]*P[n][2]+P[n][3]*P[n][3];
00247      thesquareroot=sqrt(1.+4.*alphaK[i_dom[EL]][IV]*HHM[i_dom[EL]][IV]*ksquared);
00248      energy=(thesquareroot-1.)/(2.*alphaK[i_dom[EL]][IV]);
00249      if(IV==2) energy+=EG[i_dom[EL]];
00250      energy*=Q; // <--- energy in Joule
00251 // calculates the Transmission probability
00252      k1=sqrt(2.*MSTAR[i_dom[EL]][IV]*M*energy)/HBAR;
00253      if(IV==2 && NOVALLEY[i_dom[vois[p][EL]]]==1){
00254       if((energy+(EMIN[i_dom[EL]][IV]-EMIN[i_dom[vois[p][EL]]][1])*Q)<0.) goto reflection;
00255       k2=sqrt(2.*MSTAR[i_dom[vois[p][EL]]][1]*M*
00256              (energy-(EMIN[i_dom[EL]][IV]-EMIN[i_dom[vois[p][EL]]][IV])*Q))/HBAR;
00257      }
00258      else{
00259       if((energy+(EMIN[i_dom[EL]][IV]-EMIN[i_dom[vois[p][EL]]][IV])*Q)<0.) goto reflection;
00260       k2=sqrt(2.*MSTAR[i_dom[vois[p][EL]]][IV]*M*
00261              (energy+(EMIN[i_dom[EL]][IV]-EMIN[i_dom[vois[p][EL]]][IV])*Q))/HBAR;
00262      }
00263      if(IV==2 && NOVALLEY[i_dom[vois[p][EL]]]==1){
00264        T=pow(2.*(k1/(MSTAR[i_dom[EL]][IV]*M))
00265          /(k1/(MSTAR[i_dom[EL]][IV]*M)+k2/(MSTAR[i_dom[vois[p][EL]]][1]*M)),2.);
00266        R=pow((k1/(MSTAR[i_dom[EL]][IV]*M)-k2/(MSTAR[i_dom[vois[p][EL]]][1]*M))
00267          /(k1/(MSTAR[i_dom[EL]][IV]*M)+k2/(MSTAR[i_dom[vois[p][EL]]][1]*M)),2.);
00268      }
00269      else{
00270        T=pow(2.*(k1/(MSTAR[i_dom[EL]][IV]*M))
00271          /(k1/(MSTAR[i_dom[EL]][IV]*M)+k2/(MSTAR[i_dom[vois[p][EL]]][IV]*M)),2.);
00272        R=pow((k1/(MSTAR[i_dom[EL]][IV]*M)-k2/(MSTAR[i_dom[vois[p][EL]]][IV]*M))
00273             /(k1/(MSTAR[i_dom[EL]][IV]*M)+k2/(MSTAR[i_dom[vois[p][EL]]][IV]*M)),2.);
00274      }
00275 //     printf("k1 = %g     k2 = %g   R = %g  T = %g\n",k1,k2,R*R/(R*R+T*T),T*T/(R*R+T*T)); 
00276      if(rnd()<R*R/(R*R+T*T)) goto reflection;
00277    }
00278 
00279 //   transmission:;
00280 // inside the device (between two different tetrahedra of the same material)
00281    EL=vois[p][EL];
00282 //   printf("new EL=%d\n",EL);
00283 //   printf("%g %g %g %g %g\n",L[0],L[1],L[2],L[3],L[0]+L[1]+L[2]+L[3]);
00284       memset(&a,0,sizeof(a));
00285       memset(&b,0,sizeof(b));
00286       {int k;
00287        for(k=0;k<4;k++){
00288          a[1][k+1]=coord[0][noeud_geo[k][EL]-1];
00289          a[2][k+1]=coord[1][noeud_geo[k][EL]-1];
00290          a[3][k+1]=coord[2][noeud_geo[k][EL]-1];
00291          a[4][k+1]=1.;
00292       }}
00293       b[1]=4.*x;
00294       b[2]=4.*y;
00295       b[3]=4.*z;
00296       b[4]=4.;
00297 // LU decomposition solver A.x=b
00298 #include "../linear_solver/ludcmp.h"
00299 #include "../linear_solver/lubksb.h"
00300       for(j=0;j<4;j++) L[j]=b[j+1];
00301       if(L[0]>=0. && L[1]>=0. && L[2]>=0. && L[3]>=0.){
00302 //         printf("final %g %g %g %g %g\n",L[0],L[1],L[2],L[3],L[0]+L[1]+L[2]+L[3]);
00303          return;
00304       }
00305       else goto label1;
00306 }
drift.h File Reference

Functions

Function Documentation


Functions
void	drift (double tau)