void mcgb()
{
long i,j,k,ii,ir,idr1,idr2,odr1,odr2,nx,odtf1,odtf2,idf1,idf2;
double rxin,ryin,rzin,ron,rom,vn,vm,vim,vin,nxr;
double y1,y2,rr,dux,duy,duz,kro,lnew,drmm,dtfm;
int ok;
if (ronew!=ro) {
lnew=pow((n/ronew),(1.0/3.0));kro=lnew/l;
for (i=1;i<=n;i++){
rx[i]*=kro; ry[i]*=kro; rz[i]*=kro; };
ro=ronew;
l=lnew; l2=l/2; rc=0.45*l;drmax=pow(ro,-1.0/3.0);
if (drm1>=drmax) drm1=drmax; if (drm2>=drmax) drm2=drmax;
};
vn=ugb();
idr1=0;idf1=0;idr2=0;idf2=0;odr1=0;odtf1=0;odr2=0;odtf2=0;
for (i=1;i<=nt;i++){
for (ii=1;ii<=n;ii++) {
ir=(long)(n*drand48()+1);
if (ir<=n1) {
idr1++;
if (idr1>ndr) {
if ((odr1*2)>ndr) {
if (drm1<drmax) drm1*=1.05;
} else drm1*=0.95;
odr1=0;idr1=1;
};
drmm=drm1;
}
else {
idr2++;
if (idr2>ndr) {
if ((odr2*2)>ndr) {
if (drm2<drmax) drm2*=1.05;
}
else drm2*=0.95;
odr2=0;idr2=1;
};
drmm=drm2;
};
uxi=ux[ir];uyi=uy[ir]; uzi=uz[ir];
vin=uir(ir,rx[ir],ry[ir],rz[ir]);
rxin=rx[ir]+(2*drand48()-1)*drmm;
ryin=ry[ir]+(2*drand48()-1)*drmm;
rzin=rz[ir]+(2*drand48()-1)*drmm;
if (rxin>l) rxin-=l;if (rxin<0) rxin+=l;
if (ryin>l) ryin-=l;if (ryin<0) ryin+=l;
if (rzin>l) rzin-=l;if (rzin<0) rzin+=l;
vim=uir(ir,rxin,ryin,rzin);
ok=0;
if (gbt==1) {
vm=vn-vin+vim;
if (vm<vn) ok=1; else {if (((vn-vm)/t)>log(drand48())) ok=1;};
};
if (ok==1) {
if (ir<=n1) odr1++; else odr2++;
rx[ir]=rxin; ry[ir]=ryin; rz[ir]=rzin; vn=vm; };
ir=(long)(n*drand48()+1);
if (ir<=n1) {
idf1++;
if (idf1>ndr) {
if ((odtf1*2)>ndr) {
if (dtfm1<0.75) dtfm1*=1.05;
} else dtfm1*=0.95;
odtf1=0;idf1=1;
};
dtfm=dtfm1;
}
else {
idf2++;
if (idf2>ndr) {
if ((odtf2*2)>ndr) {
if (dtfm2<0.75) dtfm2*=1.05;
}
else dtfm2*=0.95;
odtf2=0;idf2=1;
};
dtfm=dtfm2;
};
uxi=ux[ir]; uyi=uy[ir]; uzi=uz[ir];
vin=uir(ir,rx[ir],ry[ir],rz[ir]);
do {
y1=2*drand48()-1; y2=2*drand48()-1; rr=y1*y1+y2*y2; }
while (rr>1);
dux=dtfm*(2*y1*sqrt(1-rr)); duy=dtfm*(2*y2*sqrt(1-rr));duz=dtfm*(1-2*rr);
uxi+=dux; uyi+=duy; uzi+=duz;
rr=sqrt(uxi*uxi+uyi*uyi+uzi*uzi);uxi/=rr; uyi/=rr; uzi/=rr;
vim=uir(ir,rx[ir],ry[ir],rz[ir]);
ok=0;
if (gbt==1) {
vm=vn-vin+vim;
if (vm<vn) ok=1;
else {
if (((vn-vm)/t)>log(drand48())) ok=1;};
};
if (ok==1) {
if (ir<=n1) odtf1++; else odtf2++;
ux[ir]=uxi; uy[ir]=uyi; uz[ir]=uzi; vn=vm;
};
};
};
u=vn/n;
}