00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 void device_setup(void)
00030 {
00031 long int n=0;
00032 int i,np,m;
00033 double c1,c2=0.,c3,c4,c5,c6,c7;
00034
00035
00036 for(i=0;i<Ne;i++){
00037 double xa,ya,za;
00038 double xb,yb,zb;
00039 double xc,yc,zc;
00040 double xd,yd,zd;
00041 double ambx, amby, ambz;
00042 double bmcx, bmcy, bmcz;
00043 double cmdx, cmdy, cmdz;
00044
00045 xa=coord[0][noeud_geo[0][i]-1];
00046 ya=coord[1][noeud_geo[0][i]-1];
00047 za=coord[2][noeud_geo[0][i]-1];
00048
00049 xb=coord[0][noeud_geo[1][i]-1];
00050 yb=coord[1][noeud_geo[1][i]-1];
00051 zb=coord[2][noeud_geo[1][i]-1];
00052
00053 xc=coord[0][noeud_geo[2][i]-1];
00054 yc=coord[1][noeud_geo[2][i]-1];
00055 zc=coord[2][noeud_geo[2][i]-1];
00056
00057 xd=coord[0][noeud_geo[3][i]-1];
00058 yd=coord[1][noeud_geo[3][i]-1];
00059 zd=coord[2][noeud_geo[3][i]-1];
00060
00061 ambx=xa-xb; amby=ya-yb; ambz=za-zb;
00062 bmcx=xb-xc; bmcy=yb-yc; bmcz=zb-zc;
00063 cmdx=xc-xd; cmdy=yc-yd; cmdz=zc-zd;
00064
00065 VOLUME[i]=fabs(ambx*bmcy*cmdz+amby*bmcz*cmdx+ambz*bmcx*cmdy
00066 -cmdx*bmcy*ambz-cmdy*bmcz*ambx-cmdz*bmcx*amby)/6.;
00067
00068
00069
00070
00071
00072
00073
00074 }
00075
00076
00077 for(i=0;i<Ne;i++) EPP[i]=DDMAX*VOLUME[i]/NP1;
00078
00079 for(i=0;i<Ne;i++){
00080
00081
00082 np=(int)(ND[i]*VOLUME[i]/EPP[i]+0.5);
00083 if(ND[i]==0.) np=0;
00084
00085 if(vois[0][i]==-1 || vois[1][i]==-1 ||
00086 vois[2][i]==-1 || vois[3][i]==-1) np/=2;
00087 if(np>0){
00088 for(m=1;m<=np;m++){
00089 n++;
00090 if(n>NPMAX){
00091 printf("device setup error : too much particles n = %ld > %d!\n",n,NPMAX);
00092 exit(EXIT_FAILURE);
00093 }
00094
00095
00096
00097
00098
00099 IV=1;
00100 c1=log(rnd());
00101 if(i_dom[i]!=SIO2 && NOVALLEY[i_dom[i]]==1)
00102 c2=SMH[i_dom[i]][IV]*sqrt(-1.5*BKTQ*c1*(1.-alphaK[i_dom[i]][1]*1.5*BKTQ*c1));
00103 if(i_dom[i]!=SIO2 && NOVALLEY[i_dom[i]]==2){
00104 IV=1;
00105 c2=SMH[i_dom[i]][IV]*sqrt(-1.5*BKTQ*c1*(1.-alphaK[i_dom[i]][1]*1.5*BKTQ*c1));
00106 if(rnd()>0.8){
00107 IV=2;
00108 c2=SMH[i_dom[i]][IV]*sqrt(-1.5*BKTQ*c1*(1.-alphaK[i_dom[i]][2]*1.5*BKTQ*c1));
00109 }
00110 }
00111 c3=1.-2.*rnd();
00112 c4=sqrt(1.-c3*c3);
00113 c5=2.*PI*rnd();
00114 c6=sin(c5);
00115 c7=cos(c5);
00116
00117 P[n][0]=IV;
00118
00119 P[n][1]=c2*c3*c6;
00120 P[n][2]=c2*c4*c6;
00121 P[n][3]=c2*c7;
00122
00123 P[n][4]=-log(rnd())/GM[i_dom[i]];
00124
00125 P[n][5]=rnd();
00126 P[n][6]=rnd();
00127 P[n][7]=rnd();
00128 P[n][8]=4.-(P[n][5]+P[n][6]+P[n][7]);
00129
00130
00131 P[n][9]=i;
00132 if(i_dom[i]==SIO2) n--;
00133 }
00134 }
00135 }
00136 INUM = n;
00137 NPMAX0 = INUM;
00138 printf("Initial Number of Electron Super-particles = %d\n", INUM);
00139 }