0)
{
l_2 = v*t2i0; /* l_3 is the length or trace in the entering wall of the sample can, after scattering */
l_1-= v*t2i0; /* l_1 should not contain the part of the trace in the sample can if scattering happens in the entering wall */
}
else l_2=0.0; /* l_3 should be zero if the powder itself scatters, or the exiting wall of the sample can */
}
if(!cylinder_intersect(&t20, &t21, x, y, z,vout_x, vout_y, vout_z, radius+d_V, h)) /* a modifier ... */
{
if (ti0 < t0) printf("\n---> %12.6lg / %12.6lg - %12.6lg - %12.6lg / %12.6lg ---> %8.6lg %8.6lg %8.6lg - %8.6lg\n",t0,ti0,dt,ti1,t1,x,y,z,l_powder);
printf("\nCoherent: %d, powder: %d\n",coherent,powder);
printf("%lg %lg %lg (%lg %lg) %lg %lg\n",t0,ti0,dt,ti1,t1,t21,t2i1);
printf("%6.1lf deg scattering angle (%lg, %d)\n", asin(q_v/(2.0*v))*2.0*RAD2DEG,q[i]*K2V,i);
printf("FATAL ERROR: Did not hit outer cylinder from inside.\n");
exit(1); /* This can really not happen ... at least, theoretically ... (TH 10/10/99) */
}
l_2 += v*(t21-t2i1);
if (powder==0)
{
my_s = V_my_s;
l_full=l_can;
}
else
{
l_full=l_powder;
if (coherent==1) my_s = my_s_v2/(v*v); /* my_s_v2 only known for coherent case */
else my_s=my_i; /* of course ... even if it does not sound logic at the first glimpse */
}
/* l_full is now the full path without scattering through the finally scattering material, powder or can (TH 10/10/99) */
my_total_powder=my_i+my_a_v/v+my_s_v2_total/(v*v);
my_total_can =V_my_s+V_my_a_v/v;
p *= fabs(l_full*my_s*exp(-my_total_powder*(l_p+l_1)-my_total_can*(l_c+l_2)))/(1.-(double)transmission);
/*printf("%10.2lg, %10.2lg, p=%10.2lg\n",my_total_powder,my_s_v2_total,p);*/
end_counter++;
/* if (p>1) exit(1); causes trouble with new source and weight>1*/
/*
if ((unsigned long)fmod(end_counter,50) == 0)
{
if (powder==0) printf("-");
if (coherent==1) printf("*");
else if (powder==1) printf("+");
fflush(stdout);
#if __dest_os == __mac_os
Event_loop();
#endif
}
*/
}
else ABSORB; /* ... if the neutron does not hit the sample ... */
%}
FINALLY
%{
fclose(outfile);
/*
printf("\n");
for (i=1;i<=nbInt;i++)
{
printf("i=%4d N=%8d =%18.15f \n",i,Nsum[i],Psum[i]/Nsum[i]);
}
*/
%}
MCDISPLAY
%{
magnify("xyz");
circle("xz", 0, h/2.0, 0, radius);
circle("xz", 0, -h/2.0, 0, radius);
circle("xz", 0, h/2.0, 0, radius+d_V);
circle("xz", 0, -h/2.0, 0, radius+d_V);
line(-radius, -h/2.0, 0, -radius, +h/2.0, 0);
line(+radius, -h/2.0, 0, +radius, +h/2.0, 0);
line(0, -h/2.0, -radius, 0, +h/2.0, -radius);
line(0, -h/2.0, +radius, 0, +h/2.0, +radius);
line(-radius-d_V, -h/2.0, 0, -radius-d_V, +h/2.0, 0);
line(+radius+d_V, -h/2.0, 0, +radius+d_V, +h/2.0, 0);
line(0, -h/2.0, -radius-d_V, 0, +h/2.0, -radius-d_V);
line(0, -h/2.0, +radius+d_V, 0, +h/2.0, +radius+d_V);
%}
END