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
00030 void scat(int Material)
00031 {
00032 int j=0;
00033 register int i,ie=0;
00034 double ksquared,thesquareroot,superparticle_energy;
00035 double r1,finalenergy=0.,finalk,cosinus,sinus,fai;
00036 double f,cb,cf,sf,skk,a11,a12,a13,a21,a22,a23,a32,a33,x1,x2,x3,sb,r2;
00037 double ki,kf,cs,sn;
00038 double aaa=1.999,bbb=0.333;
00039 double ra1,ra2;
00040
00041 if(IV==9) return;
00042 if(Material==SIO2) return;
00043
00044
00045
00046 if(Material!=SIO2 && NOVALLEY[Material]==1){
00047 ksquared=KX*KX+KY*KY+KZ*KZ;
00048 thesquareroot=sqrt(1.+4.*alphaK[Material][1]*HHM[Material][1]*ksquared);
00049 superparticle_energy=(thesquareroot-1.)/(2.*alphaK[Material][1]);
00050 if(superparticle_energy<=0.) return;
00051 ie=((int)(superparticle_energy/DE))+1;
00052 if(ie>DIME) ie=DIME;
00053
00054
00055 r1 = rnd();
00056
00057 for(i=1;i<=6;i++){
00058
00059 if(r1<=SWK[Material][0][i*2-1][ie] && j==0){
00060 finalenergy=superparticle_energy-HWO[Material][i-1];
00061 j=1;
00062 }
00063
00064 if((r1<=SWK[Material][0][i*2][ie]) && j==0){
00065 finalenergy=superparticle_energy+HWO[Material][i-1];
00066 j=1;
00067 }
00068 }
00069
00070 if((r1<=SWK[Material][0][13][ie]) && j==0){
00071 finalenergy=superparticle_energy;
00072 j=1;
00073 }
00074
00075 if((r1<=SWK[Material][0][14][ie]) && j==0){
00076 ra1=0.5-rnd();
00077 ra2=pow(rnd(),2.);
00078 finalenergy=bbb*superparticle_energy*(1.0+aaa*ra1*ra2);
00079 finalenergy=superparticle_energy-finalenergy;
00080 j=1;
00081 }
00082 if((finalenergy<=0.) || j==0) return;
00083
00084
00085 finalk = SMH[Material][1]*sqrt(finalenergy*(1.+alphaK[Material][1]*finalenergy));
00086 cosinus = 1.-2.*rnd();
00087 sinus = sqrt(1.-cosinus*cosinus);
00088 fai = 2.*PI*rnd();
00089 KX = finalk*cosinus;
00090 KY = finalk*sinus*cos(fai);
00091 KZ = finalk*sinus*sin(fai);
00092 return;
00093 }
00094
00095
00096
00097 if(Material!=SIO2 && NOVALLEY[Material]==2){
00098 ksquared=KX*KX+KY*KY+KZ*KZ;
00099 ki=sqrt(ksquared);
00100 thesquareroot=sqrt(1.+4.*alphaK[Material][IV]*HHM[Material][IV]*ksquared);
00101 superparticle_energy=(thesquareroot-1.)/(2.*alphaK[Material][IV]);
00102 if(superparticle_energy<=0.) return;
00103 ie=((int)(superparticle_energy/DE))+1;
00104 if(ie>DIME) ie=DIME;
00105
00106
00107
00108 if(IV==1){
00109 r1 = rnd();
00110
00111
00112 if(r1<=SWK[Material][1][1][ie] && j==0){
00113 finalenergy=superparticle_energy-HWO[Material][0];
00114 if(finalenergy<0.) return;
00115 j=1;
00116
00117 kf=SMH[Material][IV]*sqrt(finalenergy*(1.+alphaK[Material][IV]*finalenergy));
00118 f=2.*ki*kf/(ki-kf)/(ki-kf);
00119 if(f<=0.) return;
00120 cb=(1.+f-pow(1.+2.*f,rnd()))/f;
00121
00122 sb=sqrt(1.-cb*cb);
00123 fai=2.*PI*rnd();
00124 cf=cos(fai);
00125 sf=sin(fai);
00126 skk=sqrt(KX*KX+KY*KY);
00127 a11=KY/skk;
00128 a12=KX*KZ/skk/ki;
00129 a13=KX/ki;
00130 a21=-KX/skk;
00131 a22=KY*KZ/skk/ki;
00132 a23=KY/ki;
00133 a32=-skk/ki;
00134 a33=KZ/ki;
00135 x1=kf*sb*cf;
00136 x2=kf*sb*sf;
00137 x3=kf*cb;
00138 KX=a11*x1+a12*x2+a13*x3;
00139 KY=a21*x1+a22*x2+a23*x3;
00140 KZ=a32*x2+a33*x3;
00141 return;
00142 }
00143
00144
00145 if((r1<=SWK[Material][1][2][ie]) && j==0){
00146 finalenergy=superparticle_energy+HWO[Material][0];
00147 if(finalenergy<0.) return;
00148 j=1;
00149
00150 kf=SMH[Material][IV]*sqrt(finalenergy*(1.+alphaK[Material][IV]*finalenergy));
00151 f=2.*ki*kf/(ki-kf)/(ki-kf);
00152 if(f<=0.) return;
00153 cb=(1.+f-pow((1.+2.*f),rnd()))/f;
00154
00155 sb=sqrt(1.-cb*cb);
00156 fai=2.*PI*rnd();
00157 cf=cos(fai);
00158 sf=sin(fai);
00159 skk=sqrt(KX*KX+KY*KY);
00160 a11=KY/skk;
00161 a12=KX*KZ/skk/ki;
00162 a13=KX/ki;
00163 a21=-KX/skk;
00164 a22=KY*KZ/skk/ki;
00165 a23=KY/ki;
00166 a32=-skk/ki;
00167 a33=KZ/ki;
00168 x1=kf*sb*cf;
00169 x2=kf*sb*sf;
00170 x3=kf*cb;
00171 KX=a11*x1+a12*x2+a13*x3;
00172 KY=a21*x1+a22*x2+a23*x3;
00173 KZ=a32*x2+a33*x3;
00174 return;
00175 }
00176
00177 if((r1<=SWK[Material][1][3][ie]) && j==0){
00178 finalenergy=superparticle_energy-HWO[Material][0]+EC[Material][1]-EC[Material][2];
00179 if(finalenergy<0.) return;
00180 IV=2;
00181 j=1;
00182
00183 kf = SMH[Material][IV]*sqrt(finalenergy*(1.+alphaK[Material][IV]*finalenergy));
00184 cs = 1.-2.*rnd();
00185 sn = sqrt(1.-cs*cs);
00186 fai = 2.*PI*rnd();
00187 KX = kf*cs;
00188 KY = kf*sn*cos(fai);
00189 KZ = kf*sn*sin(fai);
00190 return;
00191 }
00192
00193 if((r1<=SWK[Material][1][4][ie]) && j==0){
00194 finalenergy=superparticle_energy+HWO[Material][0]+EC[Material][1]-EC[Material][2];
00195 if(finalenergy<0.) return;
00196 IV=2;
00197 j=1;
00198
00199 kf = SMH[Material][IV]*sqrt(finalenergy*(1.+alphaK[Material][IV]*finalenergy));
00200 cs = 1.-2.*rnd();
00201 sn = sqrt(1.-cs*cs);
00202 fai = 2.*PI*rnd();
00203 KX = kf*cs;
00204 KY = kf*sn*cos(fai);
00205 KZ = kf*sn*sin(fai);
00206 return;
00207 }
00208
00209 if((r1<=SWK[Material][1][5][ie]) && j==0){
00210 finalenergy=superparticle_energy;
00211 if(finalenergy<0.) return;
00212 finalk=sqrt(ksquared);
00213 j=1;
00214
00215 kf = SMH[Material][IV]*sqrt(finalenergy*(1.+alphaK[Material][IV]*finalenergy));
00216 cs = 1.-2.*rnd();
00217 sn = sqrt(1.-cs*cs);
00218 fai = 2.*PI*rnd();
00219 KX = kf*cs;
00220 KY = kf*sn*cos(fai);
00221 KZ = kf*sn*sin(fai);
00222 return;
00223 }
00224
00225 if((r1<=SWK[Material][1][6][ie]) && j==0){
00226 finalenergy=superparticle_energy;
00227 if(finalenergy<0.) return;
00228 r2=rnd();
00229 cb=1.-r2/(0.5+(1.-r2)*ksquared/QD2);
00230 kf=ki;
00231
00232 sb=sqrt(1.-cb*cb);
00233 fai=2.*PI*rnd();
00234 cf=cos(fai);
00235 sf=sin(fai);
00236 skk=sqrt(KX*KX+KY*KY);
00237 a11=KY/skk;
00238 a12=KX*KZ/skk/ki;
00239 a13=KX/ki;
00240 a21=-KX/skk;
00241 a22=KY*KZ/skk/ki;
00242 a23=KY/ki;
00243 a32=-skk/ki;
00244 a33=KZ/ki;
00245 x1=kf*sb*cf;
00246 x2=kf*sb*sf;
00247 x3=kf*cb;
00248 KX=a11*x1+a12*x2+a13*x3;
00249 KY=a21*x1+a22*x2+a23*x3;
00250 KZ=a32*x2+a33*x3;
00251 return;
00252 }
00253
00254 if((r1<=SWK[Material][1][7][ie]) && j==0){
00255 ra1=0.5-rnd();
00256 ra2=pow(rnd(),2.);
00257 finalenergy=bbb*superparticle_energy*(1.0+aaa*ra1*ra2);
00258 finalenergy=superparticle_energy-finalenergy;
00259 j=1;
00260 }
00261 if((finalenergy<=0.) || j==0) return;
00262
00263
00264 finalk = SMH[Material][1]*sqrt(finalenergy*(1.+alphaK[Material][1]*finalenergy));
00265 cosinus = 1.-2.*rnd();
00266 sinus = sqrt(1.-cosinus*cosinus);
00267 fai = 2.*PI*rnd();
00268 KX = finalk*cosinus;
00269 KY = finalk*sinus*cos(fai);
00270 KZ = finalk*sinus*sin(fai);
00271 return;
00272 }
00273
00274
00275 if(IV==2){
00276 r1 = rnd();
00277
00278
00279 if(r1<=SWK[Material][2][1][ie] && j==0){
00280 finalenergy=superparticle_energy-HWO[Material][0];
00281 if(finalenergy<0.) return;
00282 j=1;
00283
00284 kf=SMH[Material][IV]*sqrt(finalenergy*(1.+alphaK[Material][IV]*finalenergy));
00285 f=2.*ki*kf/(ki-kf)/(ki-kf);
00286 if(f<=0.) return;
00287 cb=(1.+f-pow((1.+2.*f),rnd()))/f;
00288
00289 sb=sqrt(1.-cb*cb);
00290 fai=2.*PI*rnd();
00291 cf=cos(fai);
00292 sf=sin(fai);
00293 skk=sqrt(KX*KX+KY*KY);
00294 a11=KY/skk;
00295 a12=KX*KZ/skk/ki;
00296 a13=KX/ki;
00297 a21=-KX/skk;
00298 a22=KY*KZ/skk/ki;
00299 a23=KY/ki;
00300 a32=-skk/ki;
00301 a33=KZ/ki;
00302 x1=kf*sb*cf;
00303 x2=kf*sb*sf;
00304 x3=kf*cb;
00305 KX=a11*x1+a12*x2+a13*x3;
00306 KY=a21*x1+a22*x2+a23*x3;
00307 KZ=a32*x2+a33*x3;
00308 return;
00309 }
00310
00311 if((r1<=SWK[Material][2][2][ie]) && j==0){
00312 finalenergy=superparticle_energy+HWO[Material][0];
00313 if(finalenergy<0.) return;
00314 j=1;
00315
00316 kf=SMH[Material][IV]*sqrt(finalenergy*(1.+alphaK[Material][IV]*finalenergy));
00317 f=2.*ki*kf/(ki-kf)/(ki-kf);
00318 if(f<=0.) return;
00319 cb=(1.+f-pow((1.+2.*f),rnd()))/f;
00320
00321 sb=sqrt(1.-cb*cb);
00322 fai=2.*PI*rnd();
00323 cf=cos(fai);
00324 sf=sin(fai);
00325 skk=sqrt(KX*KX+KY*KY);
00326 a11=KY/skk;
00327 a12=KX*KZ/skk/ki;
00328 a13=KX/ki;
00329 a21=-KX/skk;
00330 a22=KY*KZ/skk/ki;
00331 a23=KY/ki;
00332 a32=-skk/ki;
00333 a33=KZ/ki;
00334 x1=kf*sb*cf;
00335 x2=kf*sb*sf;
00336 x3=kf*cb;
00337 KX=a11*x1+a12*x2+a13*x3;
00338 KY=a21*x1+a22*x2+a23*x3;
00339 KZ=a32*x2+a33*x3;
00340 return;
00341 }
00342
00343 if((r1<=SWK[Material][2][3][ie]) && j==0){
00344 finalenergy=superparticle_energy-HWO[Material][0];
00345 if(finalenergy<0.) return;
00346 j=1;
00347
00348 kf = SMH[Material][IV]*sqrt(finalenergy*(1.+alphaK[Material][IV]*finalenergy));
00349 cs = 1.-2.*rnd();
00350 sn = sqrt(1.-cs*cs);
00351 fai = 2.*PI*rnd();
00352 KX = kf*cs;
00353 KY = kf*sn*cos(fai);
00354 KZ = kf*sn*sin(fai);
00355 return;
00356 }
00357
00358 if((r1<=SWK[Material][2][4][ie]) && j==0){
00359 finalenergy=superparticle_energy+HWO[Material][0];
00360 if(finalenergy<0.) return;
00361 j=1;
00362
00363 kf = SMH[Material][IV]*sqrt(finalenergy*(1.+alphaK[Material][IV]*finalenergy));
00364 cs = 1.-2.*rnd();
00365 sn = sqrt(1.-cs*cs);
00366 fai = 2.*PI*rnd();
00367 KX = kf*cs;
00368 KY = kf*sn*cos(fai);
00369 KZ = kf*sn*sin(fai);
00370 return;
00371 }
00372
00373 if((r1<=SWK[Material][2][5][ie]) && j==0){
00374 finalenergy=superparticle_energy-HWO[Material][0]+EC[Material][2]-EC[Material][1];
00375 if(finalenergy<0.) return;
00376 IV=1;
00377 j=1;
00378
00379 kf = SMH[Material][IV]*sqrt(finalenergy*(1.+alphaK[Material][IV]*finalenergy));
00380 cs = 1.-2.*rnd();
00381 sn = sqrt(1.-cs*cs);
00382 fai = 2.*PI*rnd();
00383 KX = kf*cs;
00384 KY = kf*sn*cos(fai);
00385 KZ = kf*sn*sin(fai);
00386 return;
00387 }
00388
00389 if((r1<=SWK[Material][2][6][ie]) && j==0){
00390 finalenergy=superparticle_energy+HWO[Material][0]+EC[Material][2]-EC[Material][1];
00391 if(finalenergy<0.) return;
00392 IV=1;
00393 j=1;
00394
00395 kf = SMH[Material][IV]*sqrt(finalenergy*(1.+alphaK[Material][IV]*finalenergy));
00396 cs = 1.-2.*rnd();
00397 sn = sqrt(1.-cs*cs);
00398 fai = 2.*PI*rnd();
00399 KX = kf*cs;
00400 KY = kf*sn*cos(fai);
00401 KZ = kf*sn*sin(fai);
00402 return;
00403 }
00404
00405 if((r1<=SWK[Material][2][7][ie]) && j==0){
00406 finalenergy=superparticle_energy;
00407 if(finalenergy<0.) return;
00408 finalk=sqrt(ksquared);
00409 j=1;
00410
00411 kf = SMH[Material][IV]*sqrt(finalenergy*(1.+alphaK[Material][IV]*finalenergy));
00412 cs = 1.-2.*rnd();
00413 sn = sqrt(1.-cs*cs);
00414 fai = 2.*PI*rnd();
00415 KX = kf*cs;
00416 KY = kf*sn*cos(fai);
00417 KZ = kf*sn*sin(fai);
00418 return;
00419 }
00420
00421 if((r1<=SWK[Material][2][8][ie]) && j==0){
00422 finalenergy=superparticle_energy;
00423 if(finalenergy<0.) return;
00424 r2=rnd();
00425 cb=1.-r2/(0.5+(1.-r2)*ksquared/QD2);
00426 kf=ki;
00427
00428 sb=sqrt(1.-cb*cb);
00429 fai=2.*PI*rnd();
00430 cf=cos(fai);
00431 sf=sin(fai);
00432 skk=sqrt(KX*KX+KY*KY);
00433 a11=KY/skk;
00434 a12=KX*KZ/skk/ki;
00435 a13=KX/ki;
00436 a21=-KX/skk;
00437 a22=KY*KZ/skk/ki;
00438 a23=KY/ki;
00439 a32=-skk/ki;
00440 a33=KZ/ki;
00441 x1=kf*sb*cf;
00442 x2=kf*sb*sf;
00443 x3=kf*cb;
00444 KX=a11*x1+a12*x2+a13*x3;
00445 KY=a21*x1+a22*x2+a23*x3;
00446 KZ=a32*x2+a33*x3;
00447 return;
00448 }
00449
00450
00451 if((r1<=SWK[Material][2][9][ie]) && j==0){
00452 ra1=0.5-rnd();
00453 ra2=pow(rnd(),2.);
00454 finalenergy=bbb*superparticle_energy*(1.0+aaa*ra1*ra2);
00455 finalenergy=superparticle_energy-finalenergy;
00456 j=1;
00457 }
00458 if((finalenergy<=0.) || j==0) return;
00459
00460
00461 finalk = SMH[Material][2]*sqrt(finalenergy*(1.+alphaK[Material][2]*finalenergy));
00462 cosinus = 1.-2.*rnd();
00463 sinus = sqrt(1.-cosinus*cosinus);
00464 fai = 2.*PI*rnd();
00465 KX = finalk*cosinus;
00466 KY = finalk*sinus*cos(fai);
00467 KZ = finalk*sinus*sin(fai);
00468 return;
00469 }
00470 }
00471 }