using CuArrays,CUDAnative,CUDAdrv,BenchmarkTools
N = 3
k = ones(N,N)*5
M = 100
rho = zeros(N,N)
aa = range(-5,5,length=N)
K = zeros(N,N)
A = repeat(-aa’,N,1)
B = repeat(aa,1,N)
C=ones(N,N)-10
phi = repeat(range(0,2pi,length=M),1,M)
theta = repeat(range(0,pi/2,length=M),1,M)’
Rho = zeros(M,M)
Ave_Energy = ones(N,N)*pi/N^2
junyundu = 0
f = zeros(N,N)
function getnum(rho,K,f,A,B,C,phi,theta)
while junyundu < 0.9
goal_E = ones(N,N)/N^2
number_E = zeros(N,N)
E = zeros(N,N)
Threads.@threads for i = 1:M
Threads.@threads for j = 1:M
f .= @views (A[:, :].^2 .+ B[:,:].^2 .+ C[:, :] .^2).^0.5
K .= @views 2 .* k[:, :] .+ f[:,:]
rho .= @views (K[:, :] .^2 .- f[:, :].^2) ./ (2 .* (K[:, :] .- A[:, :] .* sin(theta[i, j]) .* cos(phi[i, j]) .-
B[:, :] .*sin(theta[i, j]) .* sin(phi[i, j]) .- C[:, :] .* cos(theta[i, j])))
Rho[i, j] = maximum(maximum(rho[:]))
FF = findall(x->x .== maximum(maximum(rho)),rho)
number_E[FF] = number_E[FF] .+ 1/length(FF)
end
end
E = number_E/(M*M)*pi
Var = goal_E - E
da = findall(x->x > 0, Var)
xiao = findall(x->x < 0, Var)
Var = Ave_Energy - E
fangcha = sum(sum(Var.^2))
if Var == 0
k .= @views k[:, :] .+ 0
else
k .= @views k[:, :] .+ Var[:, :] .* 0.1
end
qa = maximum(maximum(E))
qb = minimum(minimum(E))
qc = sum(sum(E))/(N*N)
global junyundu = 1-(qa-qb)/qc
println(meannum)
end
end