I am currently working on simulating random points on a circle based on a non-uniform distribution. Initially, I successfully generated these random points using a set of code lines. However, some of the simulated random points fell outside the circle when wrapping these lines within a function and calling it. Despite ensuring the correctness of the input variables for the function, my guess is that when the code executes within the function, the vector containing the probabilities associated with specific polar coordinates becomes disordered, leading to points being generated outside the circle.

What could be the issue? Sorry but the code is too long to be pasted here.

If the code is too long to share, please share an example of an input/output pair (or more than one). There is a good chance this will be enough to generate some working code, which might be even shorter.

struct containing the domain
circ = Circle(0., 0., R)

#generate random number
Np = rand(Poisson(Npar))

#create a matrix containing the cartesian position of the points
X = Array{Float64}(undef, 0, Nd);

function simulate_model(particle::Particle, circ::Circle, Np::Int64, Rk::Float64, D::Float64, X::Matrix{Float64})

for i in 1:Np
#this first block returns the correct results
x, y = GenerateHit(circ, Rk);
track = Track(x, y, Rk)
#integral contains the probability of generate a point at radius Ļ, theta the porbability of generating a point at angle Īø along radius Ļ
#Dd is a variable needed for other purposes (but it is calculated correctly)
#radius is a vector of possible radii
integral, theta, Dd, radius = distribute_D(particle, circ, track);
#next line returns points outside the circle
X_, Y_ = calculate_position(particle, circ, integral, theta, Dd, radius);
#I check if the position is within a radius circ.r, if outside is wrong
dist = sqrt.(X_[:,1].*X_[:,1] .+ X_[:,2].*X_[:,2])
if size(dist[dist .> circ.r],1) != 0
println("Error")
return X_
end
#I concatenate the matrix of positions
X = vcat(X,X_)
D += Dd
end
return X

end

The problem is the function calculate_position that gives different answers, whether executed inside the function or as a single line of code. I noted that changing variables from global to local will give different results without converging into a stable solution.

The input seems correct. I am mostly wondering how it is possible that if I run the function I get correct positions but if I run a function that calls calculate_position I get wrong positions.

To quote code enter three backquotes (`) on a separate line, to start code block, and end with another line with three backquotes.

It would be much easier to work with the code for the problematic calculate_position function, or at the minimum the parameters causing the output to go outside circle.

How are you doing your conversion from polar to cartesian coordinates and how are you sampling your \theta values? My guess is that youāre either sampling \theta outside of the [0,2\pi] range, or thereās a problem when converting to cartesian coordinates. I just played around with the example below and it seems to work fine:

using Distributions
using StatsPlots
function circle_shape(h, k, r)
Īø = LinRange(0, 2 * Ļ, 500)
h .+ r * sin.(Īø), k .+ r * cos.(Īø)
end
r = 10
circ = circle_shape(0,0,r)
d = truncated(Normal(5), 0, 2Ļ)
Īø = rand(d, 10)
x = r .* cos.(Īø)
y = r .* sin.(Īø)
plot(circ, aspect_ratio=1, legend=false)
scatter!(x,y)

It sounds like the issue youāre seeing is particular to the function definition that you havenāt shared, so Iām not sure thereās much we can do to help you debug that blindly.

I will add, though, that you might be interested in Meshes.jl which has a lot of built-in tooling for geometric definitions, parametric functions for generating points, and even random sampling tools. Hereās a quick example of

using Distributions
using Meshes
# Create a 2D surface bounded by a unit circle
center = Meshes.Point(0, 0)
radius = 1.0
circle = Meshes.Ball(center, radius)
# Some distributions from Distributions.jl with supports in [0,1]
dist_r = Normal(0.5, 0.1)
dist_phi = Normal(0.25, 0.1) # centered on pi/2
# Get a collection of N sample Meshes.Point's on the surface of the circle
# Using your geometric variable name as a function turns it into a
# parametric function, e.g. circle(r, phi) -> the Point at that location
N = 10
sample_rs = rand(dist_r, N)
sample_phis = rand(dist_phis, N)
samples = circle.(sample_rs, sample_phis)