Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Some minor bug fixes and some cleaning #96

Merged
merged 3 commits into from
Nov 12, 2020
Merged

Some minor bug fixes and some cleaning #96

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5cd0145
Now using the Rotations package instead of manually building the rota…
schustermartin Nov 12, 2020
c5e25c8
Skipping unnecesary translations and in the case of CylindricalPoints…
schustermartin Nov 12, 2020
77ebac9
Added another round to safely return an Integer. Int( Some float ) wo…
schustermartin Nov 12, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 3 additions & 10 deletions src/Geometries/VolumePrimitives/Cone.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,26 +8,19 @@ struct Cone{T} <: AbstractVolumePrimitive{T, 3} ## Only upright Cones at the mom
zStart::T
zStop::T
translate::Union{CartesianVector{T},Missing}
rotate::SMatrix{3,3,T}
rotate::Rotations.RotXYZ{T}
end



function Cone{T}( rStart1::T, rStop1::T, rStart2::T, rStop2::T, φStart::T, φStop::T, zStart::T, zStop::T, translate::Union{CartesianVector{T},Missing}, rotX::T, rotY::T, rotZ::T) where {T}
rotationMatrix::SMatrix{3,3,T} = SMatrix{3,3,T}([1 0 0;
0 cos(rotX) -sin(rotX);
0 sin(rotX) cos(rotX)]) *
SMatrix{3,3,T}([cos(rotY) 0 sin(rotY);
0 1 0;
-sin(rotY) 0 cos(rotY)]) *
SMatrix{3,3,T}([cos(rotZ) -sin(rotZ) 0;
sin(rotZ) cos(rotZ) 0;
0 0 1])
rotationMatrix::Rotations.RotXYZ{T} = RotXYZ(rotX, rotY, rotZ)

Cone{T}(rStart1, rStop1, rStart2, rStop2, φStart, φStop, zStart, zStop, translate, rotationMatrix)
end

function in(point::CylindricalPoint{T}, cone::Cone{T}) where T
(ismissing(cone.translate) || cone.translate == CartesianVector{T}(0.0, 0.0, 0.0)) ? nothing : point = CylindricalPoint(CartesianPoint(point) - cone.translate)
if ( point.z in ClosedInterval{T}(cone.zStart,cone.zStop) ) && ( point.φ in ClosedInterval{T}(cone.φStart,cone.φStop) ) && ( point.r in ClosedInterval{T}(minimum([cone.rStart1,cone.rStart2,cone.rStop1,cone.rStop2]), maximum([cone.rStart1,cone.rStart2,cone.rStop1,cone.rStop2])) )

r1,r2 = get_intersection_rs_for_given_z(point.z,cone)
Expand Down
22 changes: 8 additions & 14 deletions src/Geometries/VolumePrimitives/Tube.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,14 +30,14 @@ function Tube{T}(dict::Dict{Any, Any}, inputunit_dict::Dict{String,Unitful.Units
else
@warn "please specify a height of the Tube 'h'"
end
phi_interval = if haskey(dict, "phi")

phi_interval = if haskey(dict, "phi")
Interval(geom_round(T(deg2rad(dict["phi"]["from"]))), geom_round(T(deg2rad(dict["phi"]["to"]))))
else
Interval(geom_round(T(deg2rad(0))), geom_round(T(deg2rad(360))))
end

r_interval = if haskey(dict["r"], "from")
r_interval = if haskey(dict["r"], "from")
Interval(geom_round(ustrip(uconvert(u"m", T(dict["r"]["from"]) * inputunit_dict["length"] ))), geom_round(ustrip(uconvert(u"m", T(dict["r"]["to"]) * inputunit_dict["length"]))))
else
Interval(geom_round(0), geom_round(ustrip(uconvert(u"m", T(dict["r"]) * inputunit_dict["length"]))))
Expand All @@ -47,7 +47,7 @@ function Tube{T}(dict::Dict{Any, Any}, inputunit_dict::Dict{String,Unitful.Units
r_interval,
phi_interval,
Interval(zStart, zStop),
translate
translate
)
end

Expand All @@ -56,20 +56,14 @@ function Geometry(T::DataType, t::Val{:tube}, dict::Dict{Union{Any,String}, Any}
end

function in(point::CartesianPoint{T}, tube::Tube{T}) where T
ismissing(tube.translate) ? nothing : point -= tube.translate
(ismissing(tube.translate) || tube.translate == CartesianVector{T}(0.0,0.0,0.0)) ? nothing : point -= tube.translate
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

tube.translate should never be zero

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Well, we have to change the general struct of our primitives anyhow when we add Rotations.
This if ismissing should never appear. I guess we will use parametric types to handle this on compilation level.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For both, translation and rotation.

point = convert(CylindricalPoint,point)
if point.r in tube.r_interval && point.φ in tube.φ_interval && point.z in tube.z_interval
return true
end
return false
return (point.r in tube.r_interval && point.φ in tube.φ_interval && point.z in tube.z_interval)
end

function in(point::CylindricalPoint{T}, tube::Tube{T}) where T
ismissing(tube.translate) ? nothing : point = CylindricalPoint(CartesianPoint(point)-tube.translate)
if point.r in tube.r_interval && point.φ in tube.φ_interval && point.z in tube.z_interval
return true
end
return false
(ismissing(tube.translate) || tube.translate == CartesianVector{T}(0.0,0.0,0.0)) ? nothing : point = CylindricalPoint(CartesianPoint(point)-tube.translate)
return (point.r in tube.r_interval && point.φ in tube.φ_interval && point.z in tube.z_interval)
end


Expand Down
18 changes: 9 additions & 9 deletions src/Simulation/Simulation.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ function NamedTuple(sim::Simulation{T}) where {T <: SSDFloat}
electron_drift_field = NamedTuple(sim.electron_drift_field),
hole_drift_field = NamedTuple(sim.hole_drift_field)
)
if length(wps_strings) > 0
if length(wps_strings) > 0
nt = merge(nt, (weighting_potentials = NamedTuple{Tuple(wps_syms)}(NamedTuple.( skipmissing(sim.weighting_potentials))),))
end
return nt
Expand Down Expand Up @@ -398,7 +398,7 @@ function _calculate_potential!( sim::Simulation{T}, potential_type::UnionAll, co
only_2d::Bool = length(grid.axes[2]) == 1 ? true : false
if CS == :cylindrical
cyclic::T = grid.axes[2].interval.right - grid.axes[2].interval.left
n_φ_sym::Int = only_2d ? 1 : Int(round(T(2π) / cyclic, digits = 3))
n_φ_sym::Int = only_2d ? 1 : round(Int, round(T(2π) / cyclic, digits = 3))
n_φ_sym_info_txt = if only_2d
"φ symmetry: Detector is φ-symmetric -> 2D computation."
else
Expand Down Expand Up @@ -703,12 +703,12 @@ end



calculate_stored_energy(sim::Simulation) =
calculate_stored_energy(sim::Simulation) =
calculate_stored_energy(sim.electric_field, sim.ϵ)

function calculate_stored_energy(ef::ElectricField{T,3,S}, ϵ::DielectricDistribution{T,3,S}) where {T <: SSDFloat, S}
W::T = 0

cylindric::Bool = S == :cylindrical
cartesian::Bool = !cylindric
r0_handling::Bool = typeof(ef.grid.axes[1]).parameters[2] == :r0
Expand All @@ -725,7 +725,7 @@ function calculate_stored_energy(ef::ElectricField{T,3,S}, ϵ::DielectricDistrib
Δmp1::Vector{T} = diff(mp1)
Δmp2::Vector{T} = diff(mp2)
Δmp3::Vector{T} = diff(mp3)

w1r::Vector{T} = inv.(Δmp1) .* (mp1[2:end] .- ax1)
w1l::Vector{T} = inv.(Δmp1) .* (ax1 - mp1[1:end-1])
w2r::Vector{T} = inv.(Δmp2) .* (mp2[2:end] .- ax2)
Expand All @@ -751,11 +751,11 @@ function calculate_stored_energy(ef::ElectricField{T,3,S}, ϵ::DielectricDistrib
ev::SArray{Tuple{3},Float32,1,3} = ef.data[i1, i2, i3]
dV::T = _Δmp3 * _Δmp2 * _Δmp1
_ϵ::T = sum([
ϵ[i1, i2, i3] * w1l[i1] * w2l[i2] * w3l[i3],
ϵ[i1, i2, i3] * w1l[i1] * w2l[i2] * w3l[i3],
ϵ[i1 + 1, i2, i3] * w1r[i1] * w2l[i2] * w3l[i3],
ϵ[i1, i2 + 1, i3] * w1l[i1] * w2r[i2] * w3l[i3],
ϵ[i1, i2 + 1, i3] * w1l[i1] * w2r[i2] * w3l[i3],
ϵ[i1, i2, i3 + 1] * w1l[i1] * w2l[i2] * w3r[i3],
ϵ[i1 + 1, i2 + 1, i3] * w1r[i1] * w2r[i2] * w3l[i3],
ϵ[i1 + 1, i2 + 1, i3] * w1r[i1] * w2r[i2] * w3l[i3],
ϵ[i1 + 1, i2, i3 + 1] * w1r[i1] * w2l[i2] * w3r[i3],
ϵ[i1, i2 + 1, i3 + 1] * w1l[i1] * w2r[i2] * w3r[i3],
ϵ[i1 + 1, i2 + 1, i3 + 1] * w1r[i1] * w2r[i2] * w3r[i3]
Expand All @@ -781,4 +781,4 @@ Calculates the capacitance of an detector in Farad.
function calculate_capacitance(sim::Simulation{T}) where {T <: SSDFloat}
W = calculate_stored_energy(sim)
return uconvert(u"pF", 2 * W / (_get_abs_bias_voltage(sim.detector)^2))
end
end
12 changes: 6 additions & 6 deletions test/test_script.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
#=
This is a script to test the hole simulation chain with all standard detector types.
This is a script to test the whole simulation chain with all standard detector types.
The script also produces some output plots but only very basic ones, since it
is only to test the core functionality of the package. (So not detector specific plots)
=#
Expand All @@ -10,7 +10,7 @@ mkpath(outputdir)

@info "Loading packages"
using Plots; pyplot()
using SolidStateDetectors
using SolidStateDetectors; SSD = SolidStateDetectors

T = Float32
@info "Testing now for Float32:"
Expand Down Expand Up @@ -73,7 +73,7 @@ for key in [:InvertedCoax, :BEGe, :Coax, :CGD, :Spherical]
)
end
savefig(joinpath(outputdir, "$(key)_1_Electric_Potential_$(nref)_refinements"))
if nref != nrefs[end]
if nref != nrefs[end]
refine!(simulation, ElectricPotential, (1e-5, 1e-5, 1e-5), (1e-5, 1e-5, 1e-5), update_other_fields = true)
end
@show size(simulation.electric_potential.grid)
Expand Down Expand Up @@ -111,11 +111,11 @@ for key in [:InvertedCoax, :BEGe, :Coax, :CGD, :Spherical]
pos = if key == :InvertedCoax
CylindricalPoint{T}[ CylindricalPoint{T}( 0.02, deg2rad(10), 0.025 ) ]
elseif key == :CGD
CartesianPoint{T}[ CartesianPoint{T}( 0.006, 0.005, 0.005 ) ]
CartesianPoint{T}[ CartesianPoint{T}( 0.006, 0.005, 0.005 ) ]
elseif key == :BEGe
CylindricalPoint{T}[ CylindricalPoint{T}( 0.016, deg2rad(10), 0.015 ) ]
CylindricalPoint{T}[ CylindricalPoint{T}( 0.016, deg2rad(10), 0.015 ) ]
elseif key == :Coax
CylindricalPoint{T}[ CylindricalPoint{T}( 0.016, deg2rad(10), 0.005 ) ]
CylindricalPoint{T}[ CylindricalPoint{T}( 0.016, deg2rad(10), 0.005 ) ]
elseif key == :Spherical
CylindricalPoint{T}[ CylindricalPoint{T}( 0.00, deg2rad(0), 0.0 ) ]
end
Expand Down