Hello Andrii,
I used the phrase “The old model, ported from an earlier code, without any physical model to back it up” because that’s the unfortunate case - I tried to find a source for this equation but didn’t, as it was ported from an in-house ray tracing code. That’s why the new reflection model was developed in the first place, because even though the old gives valid-looking results, it’s not acceptable to have a simulator with a key component not backed up by verified paper(s).
Sure, my Phd thesis is public, feel free to use those approximations. As you can see in Fig. 2.32 and Fig.2.36 right side, for very low grazing angles (~1 deg), the out-of-plane angle approximation is not that good. I paste at the end the C++ routines if they help.
- Just one question: In your graphs, and in your messages, you refer often to blue dashed curves, but I don’t see any. Is it because it’s always perfectly overlapping the green curve?
Thanks, Marton
Old reflection model surface perturbation:
std::tuple<Vector3d, Vector3d, Vector3d>
ParticleTracer::PerturbateSurface(const SimulationFacet &collidedFacet,
const double &sigmaRatio) {
double rnd1 = randomGenerator.rnd();
double rnd2 = randomGenerator.rnd(); // for debug
// Saturate(rnd1, 0.01, 0.99); //Otherwise thetaOffset would go to +/-
// infinity Saturate(rnd2, 0.01, 0.99); //Otherwise phiOffset would go to +/-
// infinity
double thetaOffset = atan(sigmaRatio * tan(PI * (rnd1 - 0.5)));
double phiOffset = atan(sigmaRatio * tan(PI * (rnd2 - 0.5)));
// Make a local copy
Vector3d nU_facet = collidedFacet.sh.nU;
Vector3d nV_facet = collidedFacet.sh.nV;
Vector3d N_facet = collidedFacet.sh.N;
// Random rotation around N (to discard U orientation thus make scattering
// isotropic)
double randomAngle = randomGenerator.rnd() * 2 * PI;
nU_facet = Rotate(nU_facet, Vector3d(0, 0, 0), N_facet, randomAngle);
nV_facet = Rotate(nV_facet, Vector3d(0, 0, 0), N_facet, randomAngle);
// Bending surface base vectors
Vector3d nU_rotated =
Rotate(nU_facet, Vector3d(0, 0, 0), nV_facet, thetaOffset);
Vector3d N_rotated =
Rotate(N_facet, Vector3d(0, 0, 0), nV_facet, thetaOffset);
nU_rotated = Rotate(nU_rotated, Vector3d(0, 0, 0), nU_facet,
phiOffset); // check if correct
Vector3d nV_rotated =
Rotate(nV_facet, Vector3d(0, 0, 0), nU_facet, phiOffset);
N_rotated = Rotate(N_rotated, Vector3d(0, 0, 0), nU_facet, phiOffset);
return std::make_tuple(nU_rotated, nV_rotated, N_rotated);
}
New reflection:
void ParticleTracer::PerformBounce_new(SimulationFacet &collidedFacet,
const int &reflType, const double &inTheta,
const double &inPhi) {
double outTheta, outPhi; // perform bounce without scattering, will perturbate
// these angles later if it's a rough surface
if (collidedFacet.sh.reflectType == REFLECTION_DIFFUSE) {
outTheta = acos(sqrt(randomGenerator.rnd()));
outPhi = randomGenerator.rnd() * 2.0 * PI;
} else if (collidedFacet.sh.reflectType == REFLECTION_SPECULAR) {
outTheta = PI - inTheta;
outPhi = inPhi;
} else { // material reflection, might have backscattering
switch (reflType) {
case REFL_FORWARD: // forward scattering
outTheta = PI - inTheta;
outPhi = inPhi;
break;
case REFL_DIFFUSE: // diffuse scattering
outTheta = acos(sqrt(randomGenerator.rnd()));
outPhi = randomGenerator.rnd() * 2.0 * PI;
break;
case REFL_BACK: // back scattering
outTheta = PI - inTheta;
outPhi = PI + inPhi;
break;
} // end switch (transparent pass treated at the Intersect() routine
} // end material reflection
if (collidedFacet.sh.doScattering) {
double incidentAngle = abs(inTheta);
if (incidentAngle > PI / 2)
incidentAngle = PI - incidentAngle; // coming from the normal side
double y = cos(incidentAngle);
double wavelength =
3E8 * 6.626E-34 / (energy * 1.6E-19); // energy[eV] to wavelength[m]
double specularReflProbability =
exp(-Sqr(4 * PI * collidedFacet.sh.rmsRoughness * y /
wavelength)); // Debye-Wallers factor, See "Measurements of
// x-ray scattering..." by Dugan, Sonnad, Cimino,
// Ishibashi, Scafers, eq.2
bool specularReflection = randomGenerator.rnd() < specularReflProbability;
if (!specularReflection) {
// Smooth surface reflection performed, now let's perturbate the angles
// Using Gaussian approximated distributions of eq.14. of the above
// article
double onePerTau =
collidedFacet.sh.rmsRoughness / collidedFacet.sh.autoCorrLength;
// Old acceptance-rejection algorithm
/*
size_t nbTries = 0; double outThetaPerturbated;
do {
double dTheta = Gaussian(2.9264*onePerTau); //Grazing angle
perturbation, depends only on roughness, must assure it doesn't go against
the surface outThetaPerturbated = outTheta + dTheta; nbTries++; } while
(((outTheta < PI / 2) != (outThetaPerturbated < PI / 2)) && nbTries <
10);*/
// New truncated Gaussian algorithm, see N. Chopin: Fast simulation of
// truncated Gaussian distributions, DOI: 10.1007/s11222-009-9168-1
double lowerBound = 0.0;
double upperBound = PI / 2;
if (outTheta > (PI / 2)) { // Limits: PI/2 .. PI instead of 0..PI/2
lowerBound += PI / 2;
upperBound += PI / 2;
}
double outThetaPerturbated = truncatedGaussian.GetGaussian(
outTheta, 2.9264 * onePerTau, lowerBound, upperBound);
double dPhi = randomGenerator.Gaussian(
(2.80657 * pow(incidentAngle, -1.00238) -
1.00293 * pow(incidentAngle, 1.22425)) *
onePerTau); // Out-of-plane angle perturbation, depends on roughness
// and incident angle
outTheta = outThetaPerturbated;
outPhi += dPhi;
}
}
particle.direction = PolarToCartesian(collidedFacet.sh.nU, collidedFacet.sh.nV,
collidedFacet.sh.N, outTheta, outPhi, false);
RecordHit(HIT_REF);
lastHitFacet = &collidedFacet;
if (collidedFacet.sh.countRefl) {
RecordHitOnTexture(&collidedFacet);
}
}
The full file is here: https://gitlab.cern.ch/molflow_synrad/synrad/-/blob/master_cli/src/Simulation/Particle.cpp?ref_type=heads
I want to make my code reliable. So, currently, I’m working on implementing the new Synrad+ model into my custom-built Geant4 code. All previous options worked well, so I hope this new option goes smoothly, too.
