Go to the source code of this file.
Functions | |
| void | creation (int i, double t) |
|
||||||||||||
|
Definition at line 29 of file creation.h. References alphaK, BKTQ, EL, GM, i_dom, IV, KX, KY, KZ, L, NOVALLEY, PI, rnd(), SIO2, SMH, and TS. Referenced by ensemble_montecarlo(). 00030 { 00031 double c1,c2,c3,c4,c5,c6,c7; 00032 00033 // the following rows are just for avoiding warning messages during compilation 00034 c1=0.; 00035 c2=0.; 00036 c3=0.; 00037 c4=0.; 00038 c5=0.; 00039 c6=0.; 00040 c7=0.; 00041 00042 L[0]=rnd(); 00043 L[1]=rnd(); 00044 L[2]=rnd(); 00045 L[3]=4.-(L[0]+L[1]+L[2]); 00046 EL=i; 00047 00048 // We assume that the particles are initially 00049 // at near thermal equilibrium 00050 if(i_dom[i]!=SIO2){ 00051 IV=1; 00052 if(NOVALLEY[i_dom[i]]==2 && rnd()>0.8) IV=2; 00053 c1=log(rnd()); 00054 c2=SMH[i_dom[i]][IV]*sqrt(-1.5*BKTQ*c1*(1.-alphaK[i_dom[i]][IV]*1.5*BKTQ*c1)); 00055 c3=rnd(); 00056 c4=sqrt(1.-c3*c3); 00057 c5=2.*PI*rnd(); 00058 c6=sin(c5); 00059 c7=cos(c5); 00060 } 00061 00062 // definition of the incoming pseudo-vector 00063 // double modulus; 00064 00065 // calculates the MODULUS OF THE PSEUDO-VECTOR 00066 KX=c2*c3; 00067 KY=c2*c4*c6; 00068 KZ=c2*c4*c7; 00069 /* 00070 // The following rows are needed if one wants to insert incoming 00071 // particles with a normal velocity, normal to the tetrahedra surface. 00072 00073 modulus=sqrt(pow(KX,2.)+pow(KY,2.)+pow(KZ,2.)); 00074 00075 // calculates the internal unit vector 00076 // i.e. the unit vector which is perpendicular to the "OHMIC face" 00077 // and which points to the internal side of the tetrahedra 00078 int h,l=0; 00079 double x[4],y[4],z[4]; 00080 double nx,ny,nz; 00081 double x21,y21,z21; 00082 double x23,y23,z23; 00083 double nmodulus; 00084 00085 for(h=0;h<4;h++){ 00086 if(i_front[noeud_geo[h][i]-1]==OHMIC){ 00087 x[l]=coord[0][noeud_geo[h][i]-1]; 00088 y[l]=coord[1][noeud_geo[h][i]-1]; 00089 z[l]=coord[2][noeud_geo[h][i]-1]; 00090 l++; 00091 } 00092 if(i_front[noeud_geo[h][i]-1]!=OHMIC){ 00093 x[3]=coord[0][noeud_geo[h][i]-1]; 00094 y[3]=coord[1][noeud_geo[h][i]-1]; 00095 z[3]=coord[2][noeud_geo[h][i]-1]; 00096 } 00097 } 00098 //printf("** %g %g %g\n",x[0],y[0],z[0]); 00099 //printf("** %g %g %g\n",x[1],y[1],z[1]); 00100 //printf("** %g %g %g\n",x[2],y[2],z[2]); 00101 //printf("** %g %g %g\n",x[3],y[3],z[3]); 00102 00103 x21=x[1]-x[0]; 00104 y21=y[1]-y[0]; 00105 z21=z[1]-z[0]; 00106 x23=x[1]-x[2]; 00107 y23=y[1]-y[2]; 00108 z23=z[1]-z[2]; 00109 nx=(y21*z23-y23*z21); 00110 ny=(z21*x23-z23*x21); 00111 nz=(x21*y23-x23*y21); 00112 nmodulus=sqrt(pow(nx,2.)+pow(ny,2.)+pow(nz,2.)); 00113 nx/=nmodulus; 00114 ny/=nmodulus; 00115 nz/=nmodulus; 00116 // check if the n vector points to the internal of the tetrahedra 00117 double scalar; 00118 scalar=nx*(x[3]-x[1])+ny*(y[3]-y[1])+nz*(z[3]-z[1]); 00119 scalar/=sqrt(pow(nx*(x[3]-x[1]),2.) 00120 +pow(ny*(y[3]-y[1]),2.) 00121 +pow(nz*(z[3]-z[1]),2.)); 00122 //printf("theta = %g %g\n",acos(scalar),scalar); 00123 if(!(acos(scalar)>=0. && acos(scalar)<=0.5*PI)){ 00124 nx*=-1.; 00125 ny*=-1.; 00126 nz*=-1.; 00127 } 00128 //printf("%g %g %g %g\n",sqrt(pow(nx,2.)+pow(ny,2.)+pow(nz,2.)),nx,ny,nz); 00129 // final definition of the pseudo-wave vector 00130 KX=nx*modulus; 00131 KY=ny*modulus; 00132 KZ=nz*modulus; 00133 */ 00134 TS=t-log(rnd())/GM[i_dom[i]]; 00135 }
|