Introduced back the normalization for the integration of the velocity distribution because average velocity was multiplied back in the calculation of the flux, added the missing end for parfor loop.

This commit is contained in:
Karthik 2021-07-15 16:47:13 +02:00
parent 504b7769f9
commit b29baf29d5

View File

@ -12,7 +12,7 @@ function ret = initialVelocitySampling(this, MOTObj)
* velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ... * velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ...
* this.OvenTemperatureinKelvin))); * this.OvenTemperatureinKelvin)));
c = integral(VelocityDistribution, 0, this.VelocityCutoff); c = integral(VelocityDistribution, 0, this.VelocityCutoff) / integral(VelocityDistribution, 0, inf);
this.ReducedFlux = c * this.ReducedClausingFactor * this.calculateFreeMolecularRegimeFlux(); this.ReducedFlux = c * this.ReducedClausingFactor * this.calculateFreeMolecularRegimeFlux();
@ -31,29 +31,30 @@ function ret = initialVelocitySampling(this, MOTObj)
NormalizationConstantForVelocityDistribution = this.velocityDistributionFunction(MaximumVelocityAllowed); NormalizationConstantForVelocityDistribution = this.velocityDistributionFunction(MaximumVelocityAllowed);
parfor i = 1:n parfor i = 1:n
% Rejection Sampling of speed % Rejection Sampling of speed
y = rand(1);
x = this.VelocityCutoff * rand(1);
while y > ((NormalizationConstantForVelocityDistribution)^-1 * this.velocityDistributionFunction(x)) %As long as this loop condition is satisfied, reject the corresponding x value
y = rand(1); y = rand(1);
x = this.VelocityCutoff * rand(1); x = this.VelocityCutoff * rand(1);
end while y > ((NormalizationConstantForVelocityDistribution)^-1 * this.velocityDistributionFunction(x)) %As long as this loop condition is satisfied, reject the corresponding x value
SampledVelocityMagnitude(i) = x; % When loop condition is not satisfied, accept x value and store as sample y = rand(1);
x = this.VelocityCutoff * rand(1);
end
SampledVelocityMagnitude(i) = x; % When loop condition is not satisfied, accept x value and store as sample
% Rejection Sampling of polar angle % Rejection Sampling of polar angle
w = rand(1);
z = this.ExitDivergence * rand(1);
while w > ((NormalizationConstantForAngularDistribution)^-1 * 2 * pi * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value
w = rand(1); w = rand(1);
z = this.ExitDivergence * rand(1); z = this.ExitDivergence * rand(1);
while w > ((NormalizationConstantForAngularDistribution)^-1 * 2 * pi * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value
w = rand(1);
z = this.ExitDivergence * rand(1);
end
SampledPolarAngle(i) = z; %When loop condition is not satisfied, accept x value and store as sample
% Sampling of azimuthal angle
SampledAzimuthalAngle(i)= 2 * pi * rand(1);
ret(i,:)=[SampledVelocityMagnitude(i)*cos(SampledPolarAngle(i)), SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*cos(SampledAzimuthalAngle(i)), ...
SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*sin(SampledAzimuthalAngle(i))];
end end
SampledPolarAngle(i) = z; %When loop condition is not satisfied, accept x value and store as sample
% Sampling of azimuthal angle
SampledAzimuthalAngle(i)= 2 * pi * rand(1);
ret(i,:)=[SampledVelocityMagnitude(i)*cos(SampledPolarAngle(i)), SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*cos(SampledAzimuthalAngle(i)), ...
SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*sin(SampledAzimuthalAngle(i))];
end end