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 montecarlo_setup(int material)
00031 {
00032 double wo,no,dos,aco,oge[7],oga[7];
00033
00034 double dos1,dos2,am1,am2,cl,deq,dij;
00035 double hwe,hwij,wij,we,ne,nij;
00036 double poe,poa,ope,opa,eqe,eqa,qmin,qmax;
00037 double initialenergy,sei,finalenergy,sef;
00038 double eps,epf,ep,bimp,cimp,qd;
00039 double ak,qq,wk;
00040 int ie,i;
00041 int z2=4;
00042
00043
00044 BKTQ=KB*TL/Q;
00045 QH=Q/HBAR;
00046
00047
00048
00049 if(NOVALLEY[material]==2){
00050
00051 am1=MSTAR[material][1]*M;
00052 am2=MSTAR[material][2]*M;
00053
00054 EC[material][1]=EMIN[material][1];
00055 EC[material][2]=EMIN[material][2];
00056
00057
00058 eps=EPSR[material]*EPS0;
00059 epf=EPF[material]*EPS0;
00060 ep=1./(1./epf-1./eps);
00061
00062
00063 cl=RHO[material]*pow(UL[material],2.);
00064 deq=dij=DTK[material][0]*Q;
00065 hwe=hwij=HWO[material][0];
00066
00067 AF[material][1]=alphaK[material][1];
00068 AF[material][2]=alphaK[material][2];
00069
00070 SMH[material][1]=sqrt(2.*am1*Q)/HBAR;
00071 SMH[material][2]=sqrt(2.*am2*Q)/HBAR;
00072 HHM[material][1]=HBAR*HBAR/(2.*am1*Q);
00073 HHM[material][2]=HBAR*HBAR/(2.*am2*Q);
00074 HM[material][1]=HBAR/am1;
00075 HM[material][2]=HBAR/am2;
00076
00077 wo=HWO[material][0]*Q/HBAR;
00078 wij=hwij*Q/HBAR;
00079 we=hwe*Q/HBAR;
00080
00081 no=1./(exp(HWO[material][0]/BKTQ)-1.);
00082 nij=1./(exp(hwij/BKTQ)-1.);
00083 ne=1./(exp(hwe/BKTQ)-1.);
00084
00085 dos1=pow(sqrt(2.*am1*Q)/HBAR,3.)/pow(2.*PI,2.);
00086 dos2=pow(sqrt(2.*am2*Q)/HBAR,3.)/pow(2.*PI,2.);
00087
00088 poe=Q/8./PI/ep*Q*wo*(no+1.);
00089 poa=poe*no/(1.+no);
00090 aco=2.*PI*DA[material]/Q*DA[material]*BKTQ/HBAR*Q/cl;
00091 ope=PI*dij/wij*dij/RHO[material]/Q*(nij+1.);
00092 opa=ope*nij/(1.+nij);
00093 eqe=PI*deq/we*deq/RHO[material]/Q*(ne+1.);
00094 eqa=eqe*ne/(1.+ne);
00095
00096
00097 cimp=CIMP;
00098 qd=sqrt(Q*cimp/BKTQ/eps);
00099 QD2=qd*qd;
00100 bimp=2.*PI*cimp*Q*Q/HBAR*Q/eps/eps;
00101
00102
00103 for(ie=1;ie<=DIME;ie++){
00104 initialenergy=DE*((double)(ie));
00105 sei=sqrt(initialenergy);
00106
00107
00108
00109
00110 finalenergy=initialenergy-HWO[material][0];
00111 if(finalenergy>0.){
00112 sef=sqrt(finalenergy);
00113 qmax=sef+sei;
00114 qmin=sei-sef;
00115 SWK[material][1][1][ie]=poe*SMH[material][1]*sei/initialenergy/Q*log(qmax/qmin);
00116 }
00117 else SWK[material][1][1][ie]=0.;
00118
00119 finalenergy=initialenergy+HWO[material][0];
00120 sef=sqrt(finalenergy);
00121 qmax=sef+sei;
00122 qmin=sef-sei;
00123 SWK[material][1][2][ie]=SWK[material][1][1][ie]
00124 +poa*SMH[material][1]*sei/initialenergy/Q*log(qmax/qmin);
00125
00126 finalenergy=initialenergy-hwij+EC[material][1]-EC[material][2];
00127 if(finalenergy>0.){
00128 sef=sqrt(finalenergy*(1.+AF[material][2]*finalenergy));
00129 SWK[material][1][3][ie]=SWK[material][1][2][ie]+z2*ope*sef*dos2*
00130 (1.+2.*AF[material][2]*finalenergy);
00131 }
00132 else SWK[material][1][3][ie]=SWK[material][1][2][ie];
00133
00134 finalenergy=initialenergy+hwij+EC[material][1]-EC[material][2];
00135 if(finalenergy>0.){
00136 sef=sqrt(finalenergy*(1.+AF[material][2]*finalenergy));
00137 SWK[material][1][4][ie]=SWK[material][1][3][ie]+z2*opa*sef*dos2*
00138 (1.+2.*AF[material][2]*finalenergy);
00139 }
00140 else SWK[material][1][4][ie]=SWK[material][1][3][ie];
00141
00142 finalenergy=initialenergy;
00143 sef=sqrt(finalenergy*(1.+AF[material][1]*finalenergy));
00144 SWK[material][1][5][ie]=SWK[material][1][4][ie]+aco*sef*dos1*
00145 (1.+2.*AF[material][1]*finalenergy);
00146
00147 finalenergy=initialenergy;
00148 sef=sqrt(finalenergy*(1.+AF[material][1]*finalenergy));
00149 ak=SMH[material][1]*sef;
00150 qq=QD2*(4.*ak*ak+QD2);
00151 wk=bimp/qq*sef*dos1*(1.+2.*AF[material][1]*finalenergy);
00152 if(wk>1.e14) wk=1.e14;
00153 SWK[material][1][6][ie]=SWK[material][1][5][ie]+wk;
00154
00155
00156 wk=PII[material]*pow((initialenergy-ETH[material])/ETH[material],2.);
00157 if(IIFLAG==YES) SWK[material][1][7][ie]=SWK[material][1][6][ie]+wk;
00158 if(IIFLAG==NO) SWK[material][1][7][ie]=SWK[material][1][6][ie];
00159
00160
00161
00162
00163
00164 finalenergy=initialenergy-HWO[material][0];
00165 if(finalenergy>0.){
00166 sef=sqrt(finalenergy);
00167 qmax=sef+sei;
00168 qmin=sei-sef;
00169 SWK[material][2][1][ie]=poe*SMH[material][2]*sei/initialenergy/Q*log(qmax/qmin);
00170 }
00171 else SWK[material][2][1][ie]=0.;
00172
00173 finalenergy=initialenergy+HWO[material][0];
00174 sef=sqrt(finalenergy);
00175 qmax=sef+sei;
00176 qmin=sef-sei;
00177 SWK[material][2][2][ie]=SWK[material][2][1][ie]
00178 +poa*SMH[material][2]*sei/initialenergy/Q*log(qmax/qmin);
00179
00180
00181 finalenergy=initialenergy-hwe;
00182 if(finalenergy>0.){
00183 sef=sqrt(finalenergy*(1.+AF[material][2]*finalenergy));
00184 SWK[material][2][3][ie]=SWK[material][2][2][ie]
00185 +(z2-1.)*eqe*sef*dos2*(1.+2.*AF[material][2]*finalenergy);
00186 }
00187 else SWK[material][2][3][ie]=SWK[material][2][2][ie];
00188
00189 finalenergy=initialenergy+hwe;
00190 sef=sqrt(finalenergy*(1.+AF[material][2]*finalenergy));
00191 SWK[material][2][4][ie]=SWK[material][2][3][ie]
00192 +(z2-1.)*eqa*sef*dos2*(1.+2.*AF[material][2]*finalenergy);
00193
00194
00195 finalenergy=initialenergy-hwij+EC[material][2]-EC[material][1];
00196 if(finalenergy>0.){
00197 sef=sqrt(finalenergy*(1.+AF[material][1]*finalenergy));
00198 SWK[material][2][5][ie]=SWK[material][2][4][ie]+ope*sef*dos1*
00199 (1.+2.*AF[material][1]*finalenergy);
00200 }
00201 else SWK[material][2][5][ie]=SWK[material][2][4][ie];
00202
00203 finalenergy=initialenergy+hwij+EC[material][2]-EC[material][1];
00204 if(finalenergy>0.){
00205 sef=sqrt(finalenergy*(1.+AF[material][1]*finalenergy));
00206 SWK[material][2][6][ie]=SWK[material][2][5][ie]+opa*sef*dos1*
00207 (1.+2.*AF[material][1]*finalenergy);
00208 }
00209 else SWK[material][2][6][ie]=SWK[material][2][5][ie];
00210
00211 finalenergy=initialenergy;
00212 sef=sqrt(finalenergy*(1.+AF[material][2]*finalenergy));
00213 SWK[material][2][7][ie]=SWK[material][2][6][ie]+aco*sef*dos2*
00214 (1.+2.*AF[material][2]*finalenergy);
00215
00216 finalenergy=initialenergy;
00217 sef=sqrt(finalenergy*(1.+AF[material][2]*finalenergy));
00218 ak=SMH[material][2]*sef;
00219 qq=QD2*(4.*ak*ak+QD2);
00220 wk=bimp/qq*sef*dos2*(1.+2.*AF[material][2]*finalenergy);
00221 if(wk>1.e14) wk=1.e14;
00222 SWK[material][2][8][ie]=SWK[material][2][7][ie]+wk;
00223 }
00224
00225
00226 wk=PII[material]*pow((initialenergy-ETH[material])/ETH[material],2.);
00227 if(IIFLAG==YES) SWK[material][2][9][ie]=SWK[material][2][8][ie]+wk;
00228 if(IIFLAG==NO) SWK[material][2][9][ie]=SWK[material][2][8][ie];
00229
00230
00231 GM[material]=SWK[material][1][7][1];
00232 for(ie=1;ie<=DIME;ie++){
00233 if(SWK[material][1][7][ie]>GM[material]) GM[material]=SWK[material][1][7][ie];
00234 if(SWK[material][2][9][ie]>GM[material]) GM[material]=SWK[material][2][9][ie];
00235 }
00236 printf("GAMMA[%d] = %g \n",material, GM[material]);
00237 for(i=1;i<=7;i++)
00238 for(ie=1;ie<=DIME;ie++)
00239 SWK[material][1][i][ie]/=GM[material];
00240 for(i=1;i<=9;i++)
00241 for(ie=1;ie<=DIME;ie++)
00242 SWK[material][2][i][ie]/=GM[material];
00243 }
00244
00245
00246
00247
00248
00249 if(NOVALLEY[material]==1){
00250 SMH[material][1]=sqrt(2.*MSTAR[material][1]*M*Q)/HBAR;
00251 HHM[material][1]=HBAR*HBAR/(2.*MSTAR[material][1]*M*Q);
00252 HM[material][1]=HBAR/(MSTAR[material][1]*M);
00253
00254 dos=pow((sqrt(2.*MSTAR[material][1]*M)*sqrt(Q)/HBAR),3.)/(4.*PI*PI);
00255
00256 aco=2.*PI*(DA[material]/Q)
00257 *DA[material]*(BKTQ/HBAR)*(Q/(RHO[material]*UL[material]*UL[material]));
00258
00259 for(i=1;i<=6;i++){
00260
00261 oge[i]=0.;
00262 oga[i]=0.;
00263 if(ZF[material][i-1]!=0.){
00264 wo=HWO[material][i-1]*Q/HBAR;
00265 no=1./(exp(HWO[material][i-1]/BKTQ) - 1.);
00266 oge[i]=ZF[material][i-1]*PI*(DTK[material][i-1]*Q/wo)
00267 *((DTK[material][i-1]*Q/RHO[material])/Q)*(no+1.);
00268 oga[i]=oge[i]*no/(1.+no);
00269 }
00270 }
00271
00272 for(ie=1; ie<=DIME; ++ie) SWK[material][0][0][ie]=0.;
00273 for(ie=1; ie<=DIME; ++ie){
00274 initialenergy=DE*((double) ie);
00275 sei=sqrt(initialenergy);
00276
00277 for(i=1;i<=6;i++){
00278 finalenergy=initialenergy-HWO[material][i-1];
00279 SWK[material][0][i*2-1][ie]=SWK[material][0][i*2-2][ie];
00280 if(finalenergy>0.){
00281 sef=sqrt(finalenergy*(1.+alphaK[material][1]*finalenergy));
00282 SWK[material][0][i*2-1][ie]=SWK[material][0][i*2-2][ie]
00283 +oge[i]*sef*dos*(1.+2.*alphaK[material][1]*finalenergy);
00284 }
00285 finalenergy=initialenergy+HWO[material][i-1];
00286 SWK[material][0][i*2][ie]=SWK[material][0][i*2-1][ie];
00287 if(finalenergy>0.){
00288 sef=sqrt(finalenergy*(1.+alphaK[material][1]*finalenergy));
00289 SWK[material][0][i*2][ie]=SWK[material][0][i*2-1][ie]
00290 +oga[i]*sef*dos*(1.+2.*alphaK[material][1]*finalenergy);
00291 }
00292 }
00293
00294 finalenergy=initialenergy;
00295 sef=sqrt(finalenergy*(1.+alphaK[material][1]*finalenergy));
00296 SWK[material][0][13][ie]=SWK[material][0][12][ie]
00297 +aco*sef*dos*(1.+2.*alphaK[material][1]*finalenergy);
00298
00299
00300 wk=PII[material]*pow((initialenergy-ETH[material])/ETH[material],2.);
00301 if(IIFLAG==YES) SWK[material][0][14][ie]=SWK[material][0][13][ie]+wk;
00302 if(IIFLAG==NO) SWK[material][0][14][ie]=SWK[material][0][13][ie];
00303 }
00304
00305
00306 GM[material]=SWK[material][0][14][1];
00307 for(ie=1;ie<=DIME;++ie)
00308 if(SWK[material][0][14][ie]>GM[material]) GM[material]=SWK[material][0][14][ie];
00309 printf("GAMMA[%d] = %g \n", material, GM[material]);
00310 for(ie=1;ie<=DIME;ie++)
00311 for(i=1;i<=14;i++)
00312 SWK[material][0][i][ie]/=GM[material];
00313 }
00314
00315 }