aboutsummaryrefslogtreecommitdiff
path: root/julia/CrystFEL/src
diff options
context:
space:
mode:
Diffstat (limited to 'julia/CrystFEL/src')
-rw-r--r--julia/CrystFEL/src/CrystFEL.jl94
-rw-r--r--julia/CrystFEL/src/cell.jl430
-rw-r--r--julia/CrystFEL/src/crystal.jl121
-rw-r--r--julia/CrystFEL/src/datatemplates.jl127
-rw-r--r--julia/CrystFEL/src/detgeom.jl47
-rw-r--r--julia/CrystFEL/src/diffcalc.jl44
-rw-r--r--julia/CrystFEL/src/image.jl356
-rw-r--r--julia/CrystFEL/src/indexing.jl151
-rw-r--r--julia/CrystFEL/src/mergeutils.jl78
-rw-r--r--julia/CrystFEL/src/millepede.jl42
-rw-r--r--julia/CrystFEL/src/peaklist.jl117
-rw-r--r--julia/CrystFEL/src/peaksearch.jl74
-rw-r--r--julia/CrystFEL/src/reflists.jl411
-rw-r--r--julia/CrystFEL/src/streams.jl138
-rw-r--r--julia/CrystFEL/src/symmetry.jl134
15 files changed, 2364 insertions, 0 deletions
diff --git a/julia/CrystFEL/src/CrystFEL.jl b/julia/CrystFEL/src/CrystFEL.jl
new file mode 100644
index 00000000..897d2614
--- /dev/null
+++ b/julia/CrystFEL/src/CrystFEL.jl
@@ -0,0 +1,94 @@
+"""
+ CrystFEL
+
+Julia bindings for CrystFEL data structures and routines
+
+## Quick start
+```julia
+ using CrystFEL
+ ...
+```
+"""
+module CrystFEL
+
+libcrystfel = "libcrystfel.so"
+
+# Configure libcrystfel to use Julia's memory management. This is needed so
+# that the Julia GC knows about the memory we allocate via libcrystfel
+# routines. Otherwise, potentially very large objects will be kept hanging
+# around in memory because Julia thinks it's using a very small amount of
+# memory, and rarely runs the GC. In the case of image structures, the
+# difference between apparent and true memory use can be a factor of a million!
+function __init__()
+ @ccall libcrystfel.set_mm_funcs(cglobal(:jl_malloc)::Ptr{Cvoid},
+ cglobal(:jl_free)::Ptr{Cvoid},
+ cglobal(:jl_calloc)::Ptr{Cvoid},
+ cglobal(:jl_realloc)::Ptr{Cvoid})::Cint
+end
+
+include("cell.jl")
+using .UnitCells
+export UnitCell, LatticeType, CenteringType, UniqueAxis
+export TriclinicLattice, MonoclinicLattice, OrthorhombicLattice
+export TetragonalLattice, HexagonalLattice, RhombohedralLattice, CubicLattice
+export PrimitiveCell, ACenteredCell, BCenteredCell, CCenteredCell
+export BodyCenteredCell, FaceCenteredCell, RhombohedralCell, RhombohedralCellOnHexagonalAxes
+export NoUniqueAxis, UnknownUniqueAxis, UniqueAxisA, UniqueAxisB, UniqueAxisC
+export rotatecell
+
+include("detgeom.jl")
+using .DetGeoms
+export Panel, DetGeom
+
+include("symmetry.jl")
+using .Symmetry
+export SymOpList, asymmetricindices
+
+include("datatemplates.jl")
+using .DataTemplates
+export DataTemplate, loaddatatemplate, wavelength, cameralength
+export translategroup!, rotategroup!
+
+include("peaklist.jl")
+using .PeakLists
+export PeakList
+
+include("reflists.jl")
+using .RefLists
+export RefList, loadreflist, savereflist!
+export Reflection, MergedReflection, UnmergedReflection
+
+include("crystal.jl")
+using .Crystals
+export Crystal, InternalCrystal
+
+include("image.jl")
+using .Images
+export Image
+
+include("diffcalc.jl")
+using .DiffractionCalculations
+export predictreflections, calculatepartialities!
+export PartialityModel, UnityModel, XSphereModel, OffsetModel, RandomModel, GeneralGaussianModel
+
+include("indexing.jl")
+using .Indexing
+export Indexer, index
+
+include("streams.jl")
+using .Streams
+export Stream, chunkwrite, chunkread, allcrystals
+
+include("millepede.jl")
+using .Millepede
+export Mille
+
+include("peaksearch.jl")
+using .PeakSearch
+export zaefpeaks, peakfinder8, peakfinder9
+
+include("mergeutils.jl")
+using .MergeUtils
+export @addmeasurement, cstddev, mergereflections
+
+end # of module
diff --git a/julia/CrystFEL/src/cell.jl b/julia/CrystFEL/src/cell.jl
new file mode 100644
index 00000000..a18810bb
--- /dev/null
+++ b/julia/CrystFEL/src/cell.jl
@@ -0,0 +1,430 @@
+module UnitCells
+
+using Random
+using Printf
+
+import ..CrystFEL: libcrystfel
+export UnitCell, LatticeType, CenteringType, UniqueAxis
+export TriclinicLattice, MonoclinicLattice, OrthorhombicLattice
+export TetragonalLattice, HexagonalLattice, RhombohedralLattice, CubicLattice
+export PrimitiveCell, ACenteredCell, BCenteredCell, CCenteredCell
+export BodyCenteredCell, FaceCenteredCell, RhombohedralCell, RhombohedralCellOnHexagonalAxes
+export NoUniqueAxis, UnknownUniqueAxis, UniqueAxisA, UniqueAxisB, UniqueAxisC
+export rotatecell
+
+
+# Represents the real C-side (opaque) structure.
+mutable struct InternalUnitCell end
+
+# The Julia-side structure, needed to house the pointer to the C structure
+# Without this, we would only ever have a Ptr{DataTemplate}, not a DataTemplate.
+mutable struct UnitCell
+ internalptr::Ptr{InternalUnitCell}
+end
+
+
+"""
+Enumeration of the seven Bravais lattice types: `TriclinicLattice`,
+`MonoclinicLattice`, `OrthorhombicLattice`, `TetragonalLattice`,
+`RhombohedralLattice`, `HexagonalLattice`, `CubicLattice`.
+"""
+@enum LatticeType begin
+ TriclinicLattice
+ MonoclinicLattice
+ OrthorhombicLattice
+ TetragonalLattice
+ RhombohedralLattice
+ HexagonalLattice
+ CubicLattice
+end
+
+
+"""
+Enumeration of unit cell centering possibilities: `PrimitiveCell`,
+`ACenteredCell`, `BCenteredCell`, `CCenteredCell`, `BodyCenteredCell`
+(I-centering), `FaceCenteredCell` (F-centering), and `RhombohedralCell`
+(R-centering, primitive rhombohedral cell). `RhombohedralCellOnHexagonalAxes`
+indicates "H-centering" as used by the protein data bank, which is different
+from the "triple hexagonal cell" described in the International Tables.
+"""
+@enum CenteringType begin
+ PrimitiveCell = Int('P')
+ ACenteredCell = Int('A')
+ BCenteredCell = Int('B')
+ CCenteredCell = Int('C')
+ BodyCenteredCell = Int('I')
+ FaceCenteredCell = Int('F')
+ RhombohedralCell = Int('R')
+ RhombohedralCellOnHexagonalAxes = Int('H')
+end
+
+
+"""
+Enumeration of unique axis possibilities. The possibilities are `UniqueAxisA`,
+`UniqueAxisB` and `UniqueAxisC`. Alternatively, use `NoUniqueAxis` if the type
+of unit cell does not have a unique axis (triclinic, orthorhombic, cubic or
+rhombohedral). `UnknownUniqueAxis` means that the unique axis is not known.
+"""
+@enum UniqueAxis begin
+ NoUniqueAxis = Int('*')
+ UnknownUniqueAxis = Int('?')
+ UniqueAxisA = Int('a')
+ UniqueAxisB = Int('b')
+ UniqueAxisC = Int('c')
+end
+
+
+"""
+ UnitCell(latticetype, centering, uniqueaxis, a, b, c, α, β, γ)
+
+Creates a CrystFEL UnitCell, in an undefined orientation, from the given
+parameters. The axis lengths (a, b, c) should be in *Ångstroms*, and the
+angles (α, β, γ) should be in *degrees* - note that this is different to the
+equivalent function in CrystFEL's C API.
+
+See the documentation for `LatticeType`, `CenteringType` and `UniqueAxis` for
+possible values. You can also use the characters `'a'`, `'b'` and `'c'` for
+`uniqueaxis`, or `'P'`, `'A'`, `'B'`, `'C'`, `'I'`, `'F'`, `'R'` and `'H`' for
+`centering.
+
+Corresponds to CrystFEL C API function `cell_new_from_parameters` with follow-up
+calls to `cell_set_centering`, `cell_set_lattice_type` and `cell_set_unique_axis`.
+"""
+function UnitCell(latticetype, centering, uniqueaxis, a, b, c, α, β, γ)
+
+ ninety(a) = isapprox(a, 90, atol=0.1)
+
+ if latticetype == OrthorhombicLattice
+ if !(ninety(α) && ninety(β) && ninety(γ))
+ throw(ArgumentError("All angles must be 90°"))
+ end
+
+ elseif latticetype == TetragonalLattice
+ if !(ninety(α) && ninety(β) && ninety(γ))
+ throw(ArgumentError("All angles must be 90°"))
+ end
+
+ elseif latticetype == CubicLattice
+ if !(ninety(α) && ninety(β) && ninety(γ))
+ throw(ArgumentError("All angles must be 90°"))
+ end
+
+ elseif latticetype == HexagonalLattice
+ if uniqueaxis == UniqueAxisA
+ if !isapprox(b, c, rtol=0.01)
+ throw(ArgumentError("b and c lengths should be equal"))
+ end
+ elseif uniqueaxis == UniqueAxisB
+ if !isapprox(a, c, rtol=0.01)
+ throw(ArgumentError("a and c lengths should be equal"))
+ end
+ elseif uniqueaxis == UniqueAxisC
+ if !isapprox(a, b, rtol=0.01)
+ throw(ArgumentError("a and b lengths should be equal"))
+ end
+ else
+ throw(ArgumentError("Hexagonal cell requires a unique axis"))
+ end
+ end
+
+ out = ccall((:cell_new_from_parameters, libcrystfel),
+ Ptr{InternalUnitCell},
+ (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble,Cdouble),
+ a/1e10, b/1e10, c/1e10, deg2rad(α), deg2rad(β), deg2rad(γ))
+ if out == C_NULL
+ throw(ArgumentError("Failed to create unit cell"))
+ end
+
+ ccall((:cell_set_centering, libcrystfel),
+ Cvoid, (Ptr{InternalUnitCell},Cchar),
+ out, centering)
+
+ ccall((:cell_set_unique_axis, libcrystfel),
+ Cvoid, (Ptr{InternalUnitCell},Cchar),
+ out, uniqueaxis)
+
+ ccall((:cell_set_lattice_type, libcrystfel),
+ Cvoid, (Ptr{InternalUnitCell},Cint),
+ out, latticetype)
+
+ cell = UnitCell(out)
+
+ finalizer(cell) do x
+ ccall((:cell_free, libcrystfel),
+ Cvoid, (Ptr{InternalUnitCell},), x.internalptr)
+ end
+
+ return cell
+end
+
+
+
+"""
+ UnitCell(latticetype, centering, a, b, c, α, β, γ)
+
+A convenience constructor which attempts to determine the unique axis
+automatically from the cell parameters. If the unique axis is not obvious,
+an `ArgumentError` will be thrown.
+"""
+function UnitCell(latticetype, centering, a, b, c, α, β, γ)
+
+ notninety(a) = !isapprox(a, 90, atol=0.5)
+ ninety(a) = isapprox(a, 90, atol=0.1)
+ onetwenty(a) = isapprox(a, 120, atol=0.1)
+
+ if latticetype == TriclinicLattice
+ ua = NoUniqueAxis
+ elseif latticetype == OrthorhombicLattice
+ ua = NoUniqueAxis
+ elseif latticetype == RhombohedralLattice
+ ua = NoUniqueAxis
+ elseif latticetype == CubicLattice
+ ua = NoUniqueAxis
+
+ elseif latticetype == MonoclinicLattice
+ if notninety(α) && ninety(β) && ninety(γ)
+ ua = UniqueAxisA
+ elseif ninety(α) && notninety(β) && ninety(γ)
+ ua = UniqueAxisB
+ elseif ninety(α) && ninety(β) && notninety(γ)
+ ua = UniqueAxisC
+ else
+ throw(ArgumentError("Can't determine unique axis"))
+ end
+
+ elseif latticetype == TetragonalLattice
+ if isapprox(b, c, rtol=0.01) && !isapprox(a, b, rtol=0.05)
+ ua = UniqueAxisA
+ elseif isapprox(a, c, rtol=0.01) && !isapprox(a, b, rtol=0.05)
+ ua = UniqueAxisB
+ elseif isapprox(a, b, rtol=0.01) && !isapprox(c, b, rtol=0.05)
+ ua = UniqueAxisC
+ else
+ throw(ArgumentError("Can't determine unique axis"))
+ end
+
+ elseif latticetype == HexagonalLattice
+ if onetwenty(α) && ninety(β) && ninety(γ)
+ ua = UniqueAxisA
+ elseif ninety(α) && onetwenty(β) && ninety(γ)
+ ua = UniqueAxisB
+ elseif ninety(α) && ninety(β) && onetwenty(γ)
+ ua = UniqueAxisC
+ else
+ throw(ArgumentError("Can't determine unique axis"))
+ end
+ end
+
+ UnitCell(latticetype, centering, ua, a, b, c, α, β, γ)
+
+end
+
+
+"""
+ UnitCell(latticetype, centering, a, b, c)
+
+Construct a `UnitCell` for an `OrthorhombicLattice`, `TetragonalLattice` or
+`CubicLattice`.
+"""
+function UnitCell(latticetype::LatticeType, centering::CenteringType, a::Real, b::Real, c::Real)
+ if latticetype in (OrthorhombicLattice, TetragonalLattice, CubicLattice)
+ UnitCell(latticetype, centering, a, b, c, 90, 90, 90)
+ else
+ throw(ArgumentError("More parameters needed for this type of lattice"))
+ end
+end
+
+
+"""
+ UnitCell(CubicLattice, centering, a)
+
+Construct a `UnitCell` for a `CubicLattice`.
+"""
+function UnitCell(latticetype::LatticeType, centering::CenteringType, a::Real)
+ if latticetype == CubicLattice
+ UnitCell(latticetype, centering, a, a, a, 90, 90, 90)
+ else
+ throw(ArgumentError("More parameters needed for this type of lattice"))
+ end
+end
+
+
+"""
+ UnitCell(RhombohedralLattice, a, α)
+
+Construct a `UnitCell` for a `RhombohedralLattice`.
+"""
+function UnitCell(latticetype::LatticeType, a::Real, α::Real)
+ if latticetype == RhombohedralLattice
+ UnitCell(latticetype, RhombohedralCell, a, a, a, α, α, α)
+ else
+ throw(ArgumentError("More parameters needed for this type of lattice"))
+ end
+end
+
+
+function getlatticetype(cell)
+ lt = ccall((:cell_get_lattice_type, libcrystfel),
+ Cint, (Ptr{InternalUnitCell},), cell.internalptr)
+ cen = ccall((:cell_get_centering, libcrystfel),
+ Cchar, (Ptr{InternalUnitCell},), cell.internalptr)
+ ua = ccall((:cell_get_unique_axis, libcrystfel),
+ Cchar, (Ptr{InternalUnitCell},), cell.internalptr)
+ return LatticeType(lt),CenteringType(cen),UniqueAxis(ua)
+end
+
+
+function getcellparams(cell)
+ let a=Ref{Cdouble}(0),
+ b=Ref{Cdouble}(0),
+ c=Ref{Cdouble}(0),
+ α=Ref{Cdouble}(0),
+ β=Ref{Cdouble}(0),
+ γ=Ref{Cdouble}(0)
+ ccall((:cell_get_parameters, libcrystfel),
+ Cvoid, (Ptr{InternalUnitCell},
+ Ref{Cdouble},Ref{Cdouble},Ref{Cdouble},
+ Ref{Cdouble},Ref{Cdouble},Ref{Cdouble}),
+ cell.internalptr, a, b, c, α, β, γ)
+ return (a=a[], b=b[], c=c[], α=α[], β=β[], γ=γ[])
+ end
+end
+
+
+# Returns the direct-space basis vectors as a Julia matrix
+# See matrix-notation.pdf for information. This returns an "M-matrix".
+function directcartesianmatrix(uc)
+ ax = Ref{Cdouble}(0)
+ ay = Ref{Cdouble}(0)
+ az = Ref{Cdouble}(0)
+ bx = Ref{Cdouble}(0)
+ by = Ref{Cdouble}(0)
+ bz = Ref{Cdouble}(0)
+ cx = Ref{Cdouble}(0)
+ cy = Ref{Cdouble}(0)
+ cz = Ref{Cdouble}(0)
+ out = @ccall libcrystfel.cell_get_cartesian(uc.internalptr::Ptr{InternalUnitCell},
+ ax::Ref{Cdouble}, ay::Ref{Cdouble}, az::Ref{Cdouble},
+ bx::Ref{Cdouble}, by::Ref{Cdouble}, bz::Ref{Cdouble},
+ cx::Ref{Cdouble}, cy::Ref{Cdouble}, cz::Ref{Cdouble})::Cint
+ if out != 0
+ throw(ErrorException("Failed to convert cell parameters"))
+ end
+ return [ax[] bx[] cx[]; ay[] by[] cy[]; az[] bz[] cz[]]
+end
+
+
+# Returns the reciprocal-space basis vectors as a Julia matrix
+# See matrix-notation.pdf for information. This returns an "R-matrix".
+function reciprocalcartesianmatrix(uc)
+ ax = Ref{Cdouble}(0)
+ ay = Ref{Cdouble}(0)
+ az = Ref{Cdouble}(0)
+ bx = Ref{Cdouble}(0)
+ by = Ref{Cdouble}(0)
+ bz = Ref{Cdouble}(0)
+ cx = Ref{Cdouble}(0)
+ cy = Ref{Cdouble}(0)
+ cz = Ref{Cdouble}(0)
+ out = @ccall libcrystfel.cell_get_reciprocal(uc.internalptr::Ptr{InternalUnitCell},
+ ax::Ref{Cdouble}, ay::Ref{Cdouble}, az::Ref{Cdouble},
+ bx::Ref{Cdouble}, by::Ref{Cdouble}, bz::Ref{Cdouble},
+ cx::Ref{Cdouble}, cy::Ref{Cdouble}, cz::Ref{Cdouble})::Cint
+ if out != 0
+ throw(ErrorException("Failed to convert cell parameters"))
+ end
+ return [ax[] ay[] az[]; bx[] by[] bz[]; cx[] cy[] cz[]]
+end
+
+
+function Base.propertynames(uc::UnitCell; private=false)
+ (:a, :b, :c, :α, :β, :γ, :latticetype, :cellparams,
+ :directcartesian, :reciprocalcartesian,
+ :internalptr)
+end
+
+
+function Base.getproperty(uc::UnitCell, name::Symbol)
+ if name === :internalptr
+ getfield(uc, :internalptr)
+ elseif name === :cellparams
+ return getcellparams(uc)
+ elseif name === :latticetype
+ return getlatticetype(uc)
+ elseif name === :a
+ return getcellparams(uc).a
+ elseif name === :b
+ return getcellparams(uc).b
+ elseif name === :c
+ return getcellparams(uc).c
+ elseif name === :α
+ return getcellparams(uc).α
+ elseif name === :β
+ return getcellparams(uc).β
+ elseif name === :γ
+ return getcellparams(uc).γ
+ elseif name === :directcartesian
+ return directcartesianmatrix(uc)
+ elseif name === :reciprocalcartesian
+ return reciprocalcartesianmatrix(uc)
+ end
+end
+
+
+function Base.show(io::IO, uc::UnitCell)
+ write(io, "UnitCell(")
+ lt,cen,ua = uc.latticetype
+ show(io, lt); write(io, ", ")
+ show(io, cen); write(io, ", ")
+ show(io, ua); write(io, ",\n ")
+ @printf(io, "%.3f Å, %.3f Å, %.3f Å, %.3f°, %.3f°, %.3f°",
+ uc.a*1e10, uc.b*1e10, uc.c*1e10,
+ rad2deg(uc.α), rad2deg(uc.β), rad2deg(uc.γ))
+ write(io, ")")
+end
+
+
+# This type is for talking to libcrystfel, and is named to avoid conflicting
+# with other quaternion libraries.
+mutable struct CrystFELQuaternion
+ w::Cdouble
+ x::Cdouble
+ y::Cdouble
+ z::Cdouble
+end
+
+function randomquat()
+ r = ()->2.0*rand(Float64)-1.0
+ q = [r(), r(), r(), r()]
+ q ./ √(sum(q.^2))
+end
+
+
+"""
+ rotatecell(uc::UnitCell, quaternion)
+
+Rotate a unit cell according to a quaternion (represented as a vector of 4 floats).
+"""
+function rotatecell(uc, quat)
+ q = CrystFELQuaternion(quat...)
+ out = @ccall libcrystfel.cell_rotate(uc.internalptr::Ptr{InternalUnitCell},
+ q::CrystFELQuaternion)::Ptr{InternalUnitCell}
+ if out == C_NULL
+ throw(ErrorException("Failed to rotate unit cell"))
+ end
+ UnitCell(out)
+end
+
+
+"""
+ rotatecell(uc::UnitCell)
+
+Rotate a unit cell at random in three dimensions. Use this routine for
+simulating serial crystallography datasets.
+
+Equivalent to CrystFEL routine `cell_rotate(uc, random_quaternion(<rng>))`.
+"""
+rotatecell(uc) = rotatecell(uc, randomquat())
+
+
+end # of module
diff --git a/julia/CrystFEL/src/crystal.jl b/julia/CrystFEL/src/crystal.jl
new file mode 100644
index 00000000..f05f651f
--- /dev/null
+++ b/julia/CrystFEL/src/crystal.jl
@@ -0,0 +1,121 @@
+module Crystals
+
+using Printf
+
+import ..CrystFEL: libcrystfel
+import ..CrystFEL.RefLists: RefList, InternalRefList, UnmergedReflection
+import ..CrystFEL.UnitCells: UnitCell, InternalUnitCell
+export Crystal, InternalCrystal
+
+# Represents the real C-side (opaque) structure.
+mutable struct InternalCrystal end
+
+mutable struct Crystal
+ internalptr::Ptr{InternalCrystal}
+ cell
+end
+
+
+function Crystal(cell::UnitCell; profileradius=2e6, mosaicity=0)
+
+ out = ccall((:crystal_new, libcrystfel),
+ Ptr{InternalCrystal}, ())
+
+ if out == C_NULL
+ throw(ArgumentError("Failed to create crystal"))
+ end
+
+ # We make a copy of the cell, to avoid memory model shenanigans
+ uccopy = ccall((:cell_new_from_cell, libcrystfel),
+ Ptr{InternalUnitCell}, (Ptr{InternalUnitCell},),
+ cell.internalptr)
+
+ ccall((:crystal_set_cell, libcrystfel),
+ Cvoid, (Ptr{InternalCrystal},Ptr{InternalUnitCell}),
+ out, uccopy)
+
+ ccall((:crystal_set_profile_radius, libcrystfel),
+ Cvoid, (Ptr{InternalCrystal},Cdouble),
+ out, profileradius)
+
+ ccall((:crystal_set_mosaicity, libcrystfel),
+ Cvoid, (Ptr{InternalCrystal},Cdouble),
+ out, mosaicity)
+
+ cr = Crystal(out, nothing)
+
+ finalizer(cr) do x
+ ccall((:crystal_free, libcrystfel), Cvoid, (Ptr{InternalCrystal},),
+ x.internalptr)
+ end
+
+ return cr
+end
+
+
+function Base.setproperty!(cr::Crystal, name::Symbol, val)
+ if name === :internalptr
+ setfield!(cr, :internalptr, val)
+ end
+end
+
+
+function getcell(cr)
+ out = @ccall libcrystfel.crystal_relinquish_cell(cr.internalptr::Ptr{InternalCrystal})::Ptr{InternalUnitCell}
+ if getfield(cr, :cell) === nothing || getfield(cr, :cell).internalptr != out
+ if out != C_NULL
+ setfield!(cr, :cell, UnitCell(out))
+ else
+ setfield!(cr, :cell, nothing)
+ end
+ end
+ return getfield(cr, :cell)
+end
+
+
+function Base.getproperty(cr::Crystal, name::Symbol)
+ if name === :internalptr
+ getfield(cr, :internalptr)
+ elseif name === :cell
+ return getcell(cr)
+ elseif name === :Bfac
+ return @ccall libcrystfel.crystal_get_Bfac(cr.internalptr::Ptr{InternalCrystal})::Cdouble
+ elseif name === :osf
+ return @ccall libcrystfel.crystal_get_osf(cr.internalptr::Ptr{InternalCrystal})::Cdouble
+ elseif name === :mos
+ return @ccall libcrystfel.crystal_get_mosaicity(cr.internalptr::Ptr{InternalCrystal})::Cdouble
+ elseif name === :r
+ return @ccall libcrystfel.crystal_get_profile_radius(cr.internalptr::Ptr{InternalCrystal})::Cdouble
+ elseif name === :resolution
+ return @ccall libcrystfel.crystal_get_resolution_limit(cr.internalptr::Ptr{InternalCrystal})::Cdouble
+ elseif name === :flag
+ return @ccall libcrystfel.crystal_get_user_flag(cr.internalptr::Ptr{InternalCrystal})::Cint
+ else
+ throw(ErrorException("Type Crystal has no field "*String(name)))
+ end
+end
+
+
+function Base.show(io::IO, mime::MIME"text/plain", cr::Crystal)
+ @printf(io, "CrystFEL.Crystal(%p):\n\n", cr.internalptr)
+ if cr.cell !== nothing
+ show(io, cr.cell)
+ write(io, "\n\n")
+ else
+ write(io, "Unit cell parameters not set\n\n")
+ end
+ println(io, " Linear scale factor: ", cr.osf)
+ println(io, " Debye-Walle factor: ", cr.Bfac)
+ println(io, " Mosaicity: ", cr.mos)
+ println(io, " Profile radius: ", cr.r/1e9, " nm⁻¹")
+ println(io, " Resolution limit: ", cr.resolution)
+ println(io, " Flag: ", cr.flag)
+end
+
+
+function Base.show(io::IO, cr::Crystal)
+ @printf(io, "CrystFEL.Crystal(%p)", cr.internalptr)
+end
+
+
+end # of module
diff --git a/julia/CrystFEL/src/datatemplates.jl b/julia/CrystFEL/src/datatemplates.jl
new file mode 100644
index 00000000..b063aa90
--- /dev/null
+++ b/julia/CrystFEL/src/datatemplates.jl
@@ -0,0 +1,127 @@
+module DataTemplates
+
+import ..CrystFEL: libcrystfel
+export DataTemplate, InternalDataTemplate, loaddatatemplate
+export wavelength, cameralength
+export translategroup!, rotategroup!
+
+# Represents the real C-side (opaque) structure.
+mutable struct InternalDataTemplate end
+
+# The Julia-side structure, needed to house the pointer to the C structure
+# Without this, we would only ever have a Ptr{DataTemplate}, not a DataTemplate.
+mutable struct DataTemplate
+ internalptr::Ptr{InternalDataTemplate}
+end
+
+"""
+ loaddatatemplate(filename)
+
+Creates a CrystFEL DataTemplate by loading a geometry file.
+
+Corresponds to CrystFEL C API function `data_template_new_from_file()`.
+"""
+function loaddatatemplate(filename::AbstractString)
+
+ out = ccall((:data_template_new_from_file, libcrystfel),
+ Ptr{InternalDataTemplate}, (Cstring,), filename)
+ if out == C_NULL
+ throw(ArgumentError("Failed to load geometry file"))
+ end
+
+ dt = DataTemplate(out)
+
+ finalizer(dt) do x
+ ccall((:data_template_free, libcrystfel),
+ Cvoid, (Ptr{InternalDataTemplate},), x.internalptr)
+ end
+
+ return dt
+end
+
+
+"""
+ wavelength(datatemplate)
+
+Retrieves the radiation wavelength from a `DataTemplate`, if possible.
+
+It's not always possible to do this. Some geometry files declare the
+wavelength as a parameter to be retrieved fom each image's metadata. A
+program using this routine should take this possibility into account.
+
+Corresponds to CrystFEL C API function `data_template_get_wavelength_if_possible`.
+"""
+function wavelength(dtempl::DataTemplate)
+ wl = ccall((:data_template_get_wavelength_if_possible, libcrystfel),
+ Cdouble, (Ptr{InternalDataTemplate},), dtempl.internalptr)
+ if isnan(wl)
+ return nothing
+ else
+ return wl
+ end
+end
+
+
+"""
+ cameralength(datatemplate)
+
+Retrieves the camera length from a `DataTemplate`, if possible.
+
+It's not always possible to do this. Some geometry files declare the
+detector position(s) as a parameter to be retrieved fom each image's metadata. A
+program using this routine should take this possibility into account.
+
+Corresponds to CrystFEL C API function `data_template_get_clen_if_possible`.
+"""
+function cameralength(dtempl::DataTemplate)
+ clen = ccall((:data_template_get_clen_if_possible, libcrystfel),
+ Cdouble, (Ptr{InternalDataTemplate},), dtempl.internalptr)
+ if isnan(clen)
+ return nothing
+ else
+ return clen
+ end
+end
+
+
+"""
+ translategroup!(datatemplate, groupname, xshift, yshift, zshift)
+
+Modifies `DataTemplate` by moving the specified panel group by the specified
+amount (in metres).
+
+Corresponds to CrystFEL C API function `data_template_translate_group_m`.
+"""
+function translategroup!(dtempl::DataTemplate, groupname, xshift, yshift, zshift)
+ r = @ccall libcrystfel.data_template_translate_group_m(dtempl.internalptr::Ptr{InternalDataTemplate},
+ groupname::Cstring,
+ xshift::Cdouble,
+ yshift::Cdouble,
+ zshift::Cdouble)::Cint
+ if r != 0
+ throw(ErrorException("Failed to shift DataTemplate"))
+ end
+
+end
+
+
+"""
+ rotategroup!(datatemplate, groupname, angle, axis)
+
+Modifies `DataTemplate` by rotating the specified panel group by the specified
+amount (in degrees) about the specified xaxis (:x, :y or :z).
+
+Corresponds to CrystFEL C API function `data_template_rotate_group`.
+"""
+function rotategroup!(dtempl::DataTemplate, groupname, angle, axis)
+ r = @ccall libcrystfel.data_template_rotate_group(dtempl.internalptr::Ptr{InternalDataTemplate},
+ groupname::Cstring,
+ deg2rad(angle)::Cdouble,
+ String(axis)[1]::Cchar)::Cint
+ if r != 0
+ throw(ErrorException("Failed to rotate DataTemplate"))
+ end
+
+end
+
+end # of module
diff --git a/julia/CrystFEL/src/detgeom.jl b/julia/CrystFEL/src/detgeom.jl
new file mode 100644
index 00000000..ba943d98
--- /dev/null
+++ b/julia/CrystFEL/src/detgeom.jl
@@ -0,0 +1,47 @@
+module DetGeoms
+export Panel
+
+mutable struct Panel
+ name::Cstring
+ cx::Cdouble
+ cy::Cdouble
+ cz::Cdouble
+ pixel_pitch::Cdouble
+ adu_per_photon::Cdouble
+ max_adu::Cdouble
+ fsx::Cdouble
+ fsy::Cdouble
+ fsz::Cdouble
+ ssx::Cdouble
+ ssy::Cdouble
+ ssz::Cdouble
+ w::Cint
+ h::Cint
+ group::Ptr{Cvoid}
+end
+
+
+mutable struct DetGeom
+ panels::Ptr{Panel}
+ n_panels::Cint
+ top_group::Ptr{Cvoid}
+end
+
+
+function Base.show(io::IO, p::Panel)
+ write(io, "Panel(")
+ write(io, "name=\"")
+ write(io, unsafe_string(p.name))
+ write(io, "\", center=(")
+ show(io, p.cx); write(io, ", "); show(io, p.cy); write(io, ", "); show(io, p.cz)
+ write(io, "), fs=(")
+ show(io, p.fsx); write(io, ", "); show(io, p.fsy); write(io, ", "); show(io, p.fsz)
+ write(io, "), ss=(")
+ show(io, p.ssx); write(io, ", "); show(io, p.ssy); write(io, ", "); show(io, p.ssz)
+ write(io, "), size=(")
+ show(io, p.w); write(io, ", "); show(io, p.h)
+ write(io, "))")
+end
+
+
+end # of module
diff --git a/julia/CrystFEL/src/diffcalc.jl b/julia/CrystFEL/src/diffcalc.jl
new file mode 100644
index 00000000..1c6aea9c
--- /dev/null
+++ b/julia/CrystFEL/src/diffcalc.jl
@@ -0,0 +1,44 @@
+module DiffractionCalculations
+
+import ..CrystFEL: libcrystfel
+import ..CrystFEL.Images: InternalImage, Image
+import ..CrystFEL.Crystals: InternalCrystal, Crystal
+import ..CrystFEL.RefLists: RefList, UnmergedReflection, InternalRefList
+import ..CrystFEL.Symmetry: SymOpList
+export predictreflections, calculatepartialities!
+export PartialityModel, UnityModel, XSphereModel, OffsetModel, RandomModel, GeneralGaussianModel
+
+
+"""
+Enumeration of the available partiality models.
+"""
+@enum PartialityModel begin
+ UnityModel
+ XSphereModel
+ OffsetModel
+ RandomModel
+ GeneralGaussianModel
+end
+
+
+function predictreflections(cr::Crystal, image::Image; maxres=1e10)
+
+ refls = @ccall libcrystfel.predict_to_res(cr.internalptr::Ptr{InternalCrystal},
+ image.internalptr::Ptr{InternalImage},
+ maxres::Cdouble)::Ptr{InternalRefList}
+ sym = SymOpList("1")
+ return RefList{UnmergedReflection}(refls, sym)
+end
+
+
+function calculatepartialities!(reflist::RefList{UnmergedReflection},
+ cr::Crystal, image::Image; model=XSphereModel, maxres=1e10)
+
+ @ccall libcrystfel.calculate_partialities(reflist.internalptr::Ptr{InternalRefList},
+ cr.internalptr::Ptr{InternalCrystal},
+ image.internalptr::Ptr{InternalImage},
+ model::Cint)::Cvoid
+end
+
+
+end # of module
diff --git a/julia/CrystFEL/src/image.jl b/julia/CrystFEL/src/image.jl
new file mode 100644
index 00000000..4090dc61
--- /dev/null
+++ b/julia/CrystFEL/src/image.jl
@@ -0,0 +1,356 @@
+module Images
+
+using Printf
+
+import ..CrystFEL: libcrystfel
+import ..CrystFEL.DataTemplates: DataTemplate, InternalDataTemplate
+import ..CrystFEL.DetGeoms: DetGeom
+import ..CrystFEL.PeakLists: PeakList, InternalPeakList
+import ..CrystFEL.Crystals: Crystal, InternalCrystal
+import ..CrystFEL.RefLists: RefList, InternalRefList, UnmergedReflection
+import ..CrystFEL.Symmetry: SymOpList
+export Image
+
+const HEADER_CACHE_SIZE = 128
+
+mutable struct CrystalRefListPair
+ crystal::Ptr{InternalCrystal}
+ reflist::Ptr{InternalRefList}
+ owns_crystal::Cint
+ owns_reflist::Cint
+end
+
+mutable struct InternalImage
+ dp::Ptr{Ptr{Cfloat}}
+ bad::Ptr{Ptr{Cint}}
+ sat::Ptr{Ptr{Cfloat}}
+ hit::Cint
+ crystals::Ptr{CrystalRefListPair}
+ n_crystals::Cint
+ indexed_by::Cint
+ n_indexing_tries::Cint
+ detgeom::Ptr{DetGeom}
+ data_source_type::Cint
+ filename::Cstring
+ ev::Cstring
+ data_block::Ptr{Cvoid}
+ data_block_size::Csize_t
+ meta_data::Cstring
+ header_cache::NTuple{HEADER_CACHE_SIZE, Ptr{Cvoid}}
+ n_cached_headers::Cint
+ id::Cint
+ serial::Cint
+ spectrum::Ptr{Cvoid}
+ lambda::Cdouble
+ div::Cdouble
+ bw::Cdouble
+ peak_resolution::Cdouble
+ peaklist::Ptr{InternalPeakList}
+ ida::Ptr{Cvoid}
+ owns_peaklist::Cint
+end
+
+
+mutable struct Image
+ internalptr::Ptr{InternalImage}
+ peaklist::Union{Nothing,PeakList}
+ crystals
+ reflists
+end
+
+
+function makecrystallist(image, listptr, n)
+
+ crystals = []
+
+ if listptr == C_NULL
+ return crystals
+ end
+
+ for i in 1:n
+ pairptr = unsafe_load(listptr, i)
+
+ # Re-use old Crystal if possible
+ n = findfirst(getfield(image, :crystals)) do x
+ x.internalptr == pairptr.crystal
+ end
+
+ if n !== nothing
+ cr = getfield(image, :crystals)[n]
+ else
+ cr = Crystal(pairptr.crystal, nothing)
+ end
+
+ if pairptr.reflist == C_NULL
+ reflist = nothing
+ else
+ reflist = RefList{UnmergedReflection}(pairptr.reflist, SymOpList("1"))
+ pairptr.owns_reflist = 0
+ finalizer(reflist) do x
+ @ccall libcrystfel.reflist_free(x.internalptr::Ptr{InternalRefList})::Cvoid
+ end
+ end
+ push!(crystals, (crystal=cr, reflections=reflist))
+ pairptr.owns_crystal = 0
+ unsafe_store!(listptr, pairptr, i)
+ finalizer(cr) do x
+ @ccall libcrystfel.crystal_free(x.internalptr::Ptr{InternalCrystal})::Cvoid
+ end
+ end
+
+ image.crystals = map(x->x.crystal, crystals)
+ image.reflists = map(x->x.reflections, crystals)
+ return crystals
+
+end
+
+
+function getpeaklist(image)
+ idata = unsafe_load(image.internalptr)
+ if (getfield(image, :peaklist) === nothing) ||
+ (idata.peaklist != getfield(image, :peaklist).internalptr)
+ if idata.peaklist != C_NULL
+ setfield!(image, :peaklist, PeakList(idata.peaklist))
+ # From now on, Julia is completely responsible for freeing the peaklist
+ idata.owns_peaklist = 0
+ unsafe_store!(image.internalptr, idata)
+ else
+ setfield!(image, :peaklist, nothing)
+ end
+ end
+ return getfield(image, :peaklist)
+end
+
+
+function Base.getproperty(image::Image, name::Symbol)
+ if name === :internalptr
+ getfield(image, :internalptr)
+ elseif name === :peaklist
+ getpeaklist(image)
+ else
+ idata = unsafe_load(image.internalptr)
+
+ if name === :crystals
+ return makecrystallist(image,
+ getfield(idata, :crystals),
+ getfield(idata, :n_crystals))
+
+ else
+ getfield(idata, name)
+ end
+ end
+end
+
+
+strdup(str) = @ccall libcrystfel.cfstrdup(str::Cstring)::Cstring
+
+
+function assert_type(val, type)
+ if !(val isa type)
+ throw(ArgumentError("Must be a "*string(type)*" (have "*string(typeof(val))*" instead)"))
+ end
+end
+
+
+function set_peaklist(image, new_pl)
+
+ assert_type(new_pl, PeakList)
+
+ idata = unsafe_load(image.internalptr)
+ if (idata.owns_peaklist == 1) && (idata.peaklist != C_NULL)
+ @ccall libcrystfel.image_feature_list_free(idata.peaklist::Ptr{InternalPeakList})::Cvoid
+ end
+ idata.peaklist = new_pl.internalptr
+ idata.owns_peaklist = 0
+ unsafe_store!(image.internalptr, idata)
+ setfield!(image, :peaklist, new_pl)
+
+end
+
+
+function Base.setproperty!(image::Image, name::Symbol, val)
+ if name === :internalptr
+ setfield!(image, :internalptr, val)
+ else
+
+ if name === :peaklist
+ return set_peaklist(image, val)
+
+ elseif name === :filename
+ assert_type(val, AbstractString)
+ val = strdup(val)
+
+ elseif name === :ev
+ assert_type(val, AbstractString)
+ val = strdup(val)
+
+ elseif name === :crystals
+ return setfield!(image, :crystals, val)
+
+ elseif name === :reflists
+ return setfield!(image, :reflists, val)
+
+ end
+
+ idata = unsafe_load(image.internalptr)
+ setproperty!(idata, name, val)
+ unsafe_store!(image.internalptr, idata)
+
+ end
+end
+
+
+function Base.propertynames(image::Image; private=false)
+ if private
+ fieldnames(InternalImage)
+ else
+ tuple(fieldnames(InternalImage)..., :internalptr)
+ end
+end
+
+
+function Base.push!(image::Image, cr::Crystal)
+ @ccall libcrystfel.image_add_crystal(image.internalptr::Ptr{InternalImage},
+ cr.internalptr::Ptr{InternalCrystal})::Cvoid
+
+ idata = unsafe_load(image.internalptr)
+ ncryst = idata.n_crystals
+ pairptr = unsafe_load(idata.crystal_refls, ncryst)
+ pairptr.owns_crystal = 0
+ unsafe_store!(idata.crystal_refls, pairptr, ncryst)
+ push!(getfield(image, :crystals), cr)
+end
+
+
+function Base.push!(image::Image, cr::Crystal, reflections::RefList{UnmergedReflection})
+ @ccall libcrystfel.image_add_crystal_refls(image.internalptr::Ptr{InternalImage},
+ cr.internalptr::Ptr{InternalCrystal},
+ reflections.internalptr::Ptr{InternalRefList})::Cvoid
+ idata = unsafe_load(image.internalptr)
+ ncryst = idata.n_crystals
+ pairptr = unsafe_load(idata.crystals, ncryst)
+ pairptr.owns_crystal = 0
+ pairptr.owns_reflist = 0
+ unsafe_store!(idata.crystals, pairptr, ncryst)
+ push!(getfield(image, :crystals), cr)
+ push!(getfield(image, :reflists), reflections)
+end
+
+
+"""
+ Image(dtempl::DataTemplate)
+
+Creates a CrystFEL image structure, not linked to any file or data block,
+i.e. for simulation purposes. This will fail if `dtempl` contains any
+references to metadata fields, e.g. `photon_energy = /LCLS/photon_energy eV`.
+
+Corresponds to CrystFEL C API function `image_create_for_simulation()`.
+"""
+function Image(dtempl::DataTemplate)
+
+ out = ccall((:image_create_for_simulation, libcrystfel),
+ Ptr{Image}, (Ref{InternalDataTemplate},), dtempl.internalptr)
+ if out == C_NULL
+ throw(ArgumentError("Failed to create image"))
+ end
+
+ image = Image(out, nothing, [], [])
+
+ finalizer(image) do x
+ ccall((:image_free, libcrystfel), Cvoid, (Ptr{InternalImage},), x.internalptr)
+ end
+
+ return image
+end
+
+
+"""
+ Image(dtempl::DataTemplate, filename::AbstractString, event::AbstractString,
+ no_image_data=false, no_mask_data=false)
+
+Loads an image from the filesystem.
+
+Corresponds to CrystFEL C API function `image_read()`.
+"""
+function Image(dtempl::DataTemplate,
+ filename::AbstractString,
+ event::AbstractString="//",
+ no_image_data=false,
+ no_mask_data=false)
+
+ out = @ccall libcrystfel.image_read(dtempl.internalptr::Ptr{InternalDataTemplate},
+ filename::Cstring, event::Cstring,
+ no_image_data::Cint, no_mask_data::Cint,
+ C_NULL::Ptr{Cvoid})::Ptr{Image}
+ if out == C_NULL
+ throw(ArgumentError("Failed to load image"))
+ end
+
+ image = Image(out, nothing, [], [])
+
+ finalizer(image) do x
+ ccall((:image_free, libcrystfel), Cvoid, (Ptr{InternalImage},), x.internalptr)
+ end
+
+ return image
+end
+
+
+function Base.show(io::IO, mime::MIME"text/plain", image::Image)
+
+ idata = unsafe_load(image.internalptr)
+ @printf(io, "CrystFEL.Image(%p):\n\n", image.internalptr)
+
+ println(io, " Serial number: ", idata.serial)
+ write(io, " Filename: ")
+ if idata.filename == C_NULL
+ write(io, "<not set>")
+ else
+ write(io, unsafe_string(idata.filename))
+ end
+ write(io, "\n")
+
+ write(io, " Frame ID: ")
+ if idata.ev == C_NULL
+ write(io, "<not set>")
+ else
+ write(io, unsafe_string(idata.ev))
+ end
+ write(io, "\n")
+
+ write(io, "\n")
+ println(io, " Wavelength: ", idata.lambda*1e10, " Å")
+ println(io, " Bandwidth: ", idata.bw*100, " %")
+ println(io, " Divergence: ", idata.div*1e3, " mrad")
+
+ write(io, "\n")
+ if idata.peaklist != C_NULL
+ let npk = @ccall libcrystfel.image_feature_count(idata.peaklist::Ptr{InternalPeakList})::Cint
+ println(io, " Number of peaks: ", npk)
+ end
+ else
+ println(io, " Number of peaks: 0 (no peak list)")
+ end
+
+ println(io, " Estimated resolution: ", 1e10/idata.peak_resolution, " Å")
+ write(io, " Hit flag: ")
+ if idata.hit != 0
+ write(io, "set")
+ else
+ write(io, "not set")
+ end
+ write(io, "\n")
+
+ write(io, "\n")
+ println(io, " Number of crystals: ", idata.n_crystals)
+ println(io, " Number of indexing attempts made: ", idata.n_crystals)
+ println(io, " Indexed by algorithm: ", idata.indexed_by)
+end
+
+
+function Base.show(io::IO, image::Image)
+ @printf(io, "CrystFEL.Image(%p)", image.internalptr)
+end
+
+
+end # of module
diff --git a/julia/CrystFEL/src/indexing.jl b/julia/CrystFEL/src/indexing.jl
new file mode 100644
index 00000000..543246ec
--- /dev/null
+++ b/julia/CrystFEL/src/indexing.jl
@@ -0,0 +1,151 @@
+module Indexing
+
+import ..CrystFEL: libcrystfel
+import ..CrystFEL.UnitCells: UnitCell, InternalUnitCell
+import ..CrystFEL.Images: Image, InternalImage
+import ..CrystFEL.DataTemplates: wavelength, cameralength
+export Indexer, index
+
+mutable struct IndexingPriv end
+
+mutable struct Indexer
+ indexingpriv::Ptr{IndexingPriv}
+end
+
+
+function indexflags(retry, multilattice, refine, peakcheck, cellcheck)
+ flags = 0
+ if retry
+ flags |= 1
+ end
+ if multilattice
+ flags |= 2
+ end
+ if refine
+ flags |= 4
+ end
+ if peakcheck
+ flags |= 32
+ end
+ if cellcheck
+ flags |= 64
+ end
+ return flags
+end
+
+
+"""
+ Indexer(methods, dtempl, cell;
+ tolerances=(0.05,0.05,0.05,1.5,1.5,1.5),
+ retry=true, multilattice=false, refine=true,
+ peakcheck=true, cellcheck=true,
+ wavelength_estimate=nothing, clen_estimate=missing,
+ n_threads=1)
+
+Creates a new CrystFEL indexing engine, which you can later apply to CrystFEL
+`Image` structures using `index()`.
+"""
+function Indexer(methods, dtempl, cell; tolerances=(0.05,0.05,0.05,1.5,1.5,1.5),
+ retry=true, multilattice=false, refine=true, peakcheck=true, cellcheck=true,
+ wavelength_estimate=nothing, clen_estimate=missing, n_threads=1)
+
+ taketwoopts = Ref{Ptr{Cvoid}}(C_NULL)
+ xgandalfopts = Ref{Ptr{Cvoid}}(C_NULL)
+ pinkindexeropts = Ref{Ptr{Cvoid}}(C_NULL)
+ felixopts = Ref{Ptr{Cvoid}}(C_NULL)
+ fromfileopts = Ref{Ptr{Cvoid}}(C_NULL)
+ asdfopts = Ref{Ptr{Cvoid}}(C_NULL)
+
+ @ccall libcrystfel.default_method_options(taketwoopts::Ref{Ptr{Cvoid}},
+ xgandalfopts::Ref{Ptr{Cvoid}},
+ pinkindexeropts::Ref{Ptr{Cvoid}},
+ felixopts::Ref{Ptr{Cvoid}},
+ fromfileopts::Ref{Ptr{Cvoid}},
+ asdfopts::Ref{Ptr{Cvoid}})::Cvoid
+
+ flags = indexflags(retry, multilattice, refine, peakcheck, cellcheck)
+
+ let wlfromdtempl = wavelength(dtempl)
+ if !isnothing(wlfromdtempl)
+ wavelength_estimate = wlfromdtempl
+ else
+ if isnothing(wavelength_estimate)
+ throw(ArgumentError("Wavelength cannot be determined from data template. "*
+ "Use Indexer(wavelength_estimate=...)"))
+ end
+ end
+ end
+
+ let clenfromdtempl = cameralength(dtempl)
+ if !isnothing(clenfromdtempl)
+ clen_estimate = clenfromdtempl
+ else
+ if isnothing(clen_estimate)
+ throw(ArgumentError("Camera length cannot be determined from data template. "*
+ "Use Indexer(clen_estimate=...)"))
+ end
+ end
+ end
+
+ tols = Vector{Cfloat}(undef, 6)
+ tols[1] = tolerances[1]
+ tols[2] = tolerances[2]
+ tols[3] = tolerances[3]
+ tols[4] = deg2rad(tolerances[4])
+ tols[5] = deg2rad(tolerances[5])
+ tols[6] = deg2rad(tolerances[6])
+
+ out = @ccall libcrystfel.setup_indexing(methods::Cstring,
+ cell.internalptr::Ptr{InternalUnitCell},
+ tols::Ptr{Cfloat},
+ flags::Cint,
+ wavelength_estimate::Cdouble,
+ clen_estimate::Cdouble,
+ n_threads::Cint,
+ taketwoopts[]::Ptr{Cvoid},
+ xgandalfopts[]::Ptr{Cvoid},
+ pinkindexeropts[]::Ptr{Cvoid},
+ felixopts[]::Ptr{Cvoid},
+ fromfileopts[]::Ptr{Cvoid},
+ asdfopts[]::Ptr{Cvoid})::Ptr{IndexingPriv}
+
+ if out == C_NULL
+ throw(ErrorException("Indexing setup failed"))
+ end
+
+ indexer = Indexer(out)
+
+ finalizer(indexer) do x
+ @ccall libcrystfel.cleanup_indexing(x.indexingpriv::Ptr{IndexingPriv})::Cvoid
+ end
+
+ return indexer
+
+end
+
+
+"""
+ index(image::Image, indexer::Indexer; mille=nothing)
+
+Index `image` using `indexer`.
+
+If `mille` is a valid `Mille` object, detector geometry alignment data will
+be written.
+"""
+function index(image::Image, idxr::Indexer; mille=nothing)
+
+ if mille === nothing
+ imille = C_NULL
+ else
+ imille = mille.internalptr
+ end
+
+ @ccall libcrystfel.index_pattern_4(image.internalptr::Ptr{InternalImage},
+ idxr.indexingpriv::Ptr{IndexingPriv},
+ C_NULL::Ptr{Cvoid},
+ C_NULL::Ptr{Cvoid},
+ imille::Ptr{Cvoid})::Cvoid
+end
+
+
+end # of module
diff --git a/julia/CrystFEL/src/mergeutils.jl b/julia/CrystFEL/src/mergeutils.jl
new file mode 100644
index 00000000..701992db
--- /dev/null
+++ b/julia/CrystFEL/src/mergeutils.jl
@@ -0,0 +1,78 @@
+module MergeUtils
+
+import ..CrystFEL: libcrystfel
+using ..CrystFEL.RefLists
+using ..CrystFEL.Symmetry
+import ..CrystFEL.UnitCells: InternalUnitCell
+export @addmeasurement, cstddev, mergereflections
+
+
+macro addmeasurement(measurement, weight,
+ mean, sumweight, wksp)
+ return quote
+ delta = $(esc(measurement)) - $(esc(mean))
+ newsumweight = $(esc(sumweight)) + $(esc(weight))
+ R = delta * $(esc(weight)) / newsumweight
+ $(esc(mean)) += R
+ $(esc(wksp)) += $(esc(sumweight)) * delta * R
+ $(esc(sumweight)) = newsumweight
+ end
+end
+
+cstddev(nmeas, work1, work2) = sqrt(work2/work1)/sqrt(nmeas)
+
+
+struct Polarisation
+ fraction::Cdouble
+ angle::Cdouble
+ disable::Cint
+end
+
+function polarisation_correction!(reflist, cell, polfrac, polangle)
+ pol = Polarisation(polfrac, rad2deg(polangle), 0)
+ @ccall libcrystfel.polarisation_correction(reflist.internalptr::Ptr{InternalRefList},
+ cell.internalptr::Ptr{InternalUnitCell},
+ pol::Ref{Polarisation})::Cvoid
+end
+
+
+function mergereflections(correction, crystalrefls, sym)
+
+ merged = RefList{MergedReflection}(sym)
+
+ for (cr,reflections) in crystalrefls
+
+ polarisation_correction!(reflections, cr.cell, 1.0, 0.0)
+
+ for refl in reflections
+
+ indices = asymmetricindices(sym, refl.indices)
+ model_version = get!(merged, indices)
+ @addmeasurement(correction(refl, cr), 1.0,
+ model_version.intensity,
+ model_version.temp1,
+ model_version.temp2)
+ model_version.nmeasurements += 1
+
+ end
+ end
+
+ for refl in merged
+ if refl.nmeasurements > 1
+ refl.sigintensity = cstddev(refl.nmeasurements, refl.temp1, refl.temp2)
+ else
+ # Cannot delete a reflection from a list (especially not
+ # while iterating), but setting nmeasurements to zero
+ # prevents it from being written to the output file.
+ refl.nmeasurements = 0
+ end
+ end
+
+ return merged
+
+end
+
+mergereflections(crystalrefls, sym) = mergereflections((refl,cr)->refl.intensity, crystalrefls, sym)
+
+
+end # of module
diff --git a/julia/CrystFEL/src/millepede.jl b/julia/CrystFEL/src/millepede.jl
new file mode 100644
index 00000000..a65256dc
--- /dev/null
+++ b/julia/CrystFEL/src/millepede.jl
@@ -0,0 +1,42 @@
+module Millepede
+
+import ..CrystFEL: libcrystfel
+export Mille
+
+mutable struct InternalMille end
+
+mutable struct Mille
+ internalptr::Ptr{InternalMille}
+end
+
+
+"""
+ Mille(filename)
+
+Creates a new "Mille" object, which can be passed to `CrystFEL.Indexing.index()`
+to accumulate detector geometry alignment information in the specified file.
+
+When you've finished adding data, call `close()` on the Mille object. This will
+be done automatically when the object is freed by the garbage collector, but
+that might not happen until Julia exits.
+
+Corresponds to CrystFEL C API routine `crystfel_mille_new`.
+"""
+function Mille(filename::AbstractString)
+ out = @ccall libcrystfel.crystfel_mille_new(filename::Cstring)::Ptr{InternalMille}
+ if out == C_NULL
+ throw(ArgumentError("Failed to create Millepede object"))
+ end
+ finalizer(close, Mille(out))
+end
+
+
+function Base.close(x::Mille)
+ if x.internalptr != C_NULL
+ @ccall libcrystfel.crystfel_mille_free(x.internalptr::Ptr{InternalMille})::Cvoid
+ x.internalptr = C_NULL
+ end
+end
+
+
+end # of module
diff --git a/julia/CrystFEL/src/peaklist.jl b/julia/CrystFEL/src/peaklist.jl
new file mode 100644
index 00000000..9eedad27
--- /dev/null
+++ b/julia/CrystFEL/src/peaklist.jl
@@ -0,0 +1,117 @@
+module PeakLists
+
+using Printf
+import ..CrystFEL: libcrystfel
+export PeakList, InternalPeakList
+
+mutable struct InternalPeak
+ fs::Cdouble
+ ss::Cdouble
+ panelnumber::Cint
+ intensity::Cdouble
+ name::Cstring
+end
+
+mutable struct InternalPeakList end
+
+mutable struct PeakList
+ internalptr::Ptr{InternalPeakList}
+end
+
+
+function PeakList()
+ out = @ccall libcrystfel.image_feature_list_new()::Ptr{InternalPeakList}
+ if out == C_NULL
+ throw(ArgumentError("Failed to create peak list"))
+ end
+ finalizer(PeakList(out)) do pl
+ @ccall libcrystfel.image_feature_list_free(pl.internalptr::Ptr{InternalPeakList})::Cvoid
+ end
+end
+
+
+Base.copy(::PeakList) = throw(ErrorException("Cannot copy(::PeakList), must deepcopy"))
+
+function Base.deepcopy(peaklist::PeakList)
+ out = @ccall libcrystfel.image_feature_list_copy(peaklist.internalptr::Ptr{InternalPeakList})::Ptr{InternalPeakList}
+ if out == C_NULL
+ throw(ErrorException("Failed to copy peak list"))
+ end
+ PeakList(out)
+end
+
+function Base.length(peaklist::PeakList)
+ @ccall libcrystfel.image_feature_count(peaklist.internalptr::Ptr{InternalPeakList})::Cint
+end
+
+Base.firstindex(peaklist::PeakList) = 1
+Base.lastindex(peaklist::PeakList) = length(peaklist)
+
+
+function Base.push!(peaklist::PeakList, fs, ss, panelnumber, intensity, name=nothing)
+
+ rname = isnothing(name) ? C_NULL : @ccall strdup(name::Cstring)::Cstring
+
+ @ccall libcrystfel.image_add_feature(peaklist.internalptr::Ptr{InternalPeakList},
+ fs::Cdouble, ss::Cdouble, panelnumber::Cint,
+ intensity::Cdouble, rname::Cstring)::Cvoid
+end
+
+
+function Base.getindex(peaklist::PeakList, n)
+ out = @ccall(libcrystfel.image_get_feature(peaklist.internalptr::Ptr{InternalPeakList},
+ (n-1)::Cint)::Ptr{InternalPeak})
+ if out == C_NULL
+ throw(BoundsError(peaklist, n))
+ end
+ pdata = unsafe_load(out)
+ if pdata.name == C_NULL
+ nname = nothing
+ else
+ nname = unsafe_string(pdata.name)
+ end
+ return (fs=pdata.fs, ss=pdata.ss, panelnumber=pdata.panelnumber,
+ intensity=pdata.intensity, name=nname)
+end
+
+
+function Base.iterate(peaklist::PeakList)
+ if length(peaklist) > 0
+ return peaklist[1],(1,length(peaklist))
+ else
+ return nothing
+ end
+end
+
+
+function Base.iterate(peaklist::PeakList, state)
+ let nxt = state[1]+1
+ len = state[2]
+ if nxt == len+1
+ return nothing
+ else
+ return peaklist[nxt],(nxt,state[2])
+ end
+ end
+end
+
+
+function Base.show(io::IO, ::MIME"text/plain", peaklist::PeakList)
+ println(io, "Peak list with ", length(peaklist), " peaks")
+ print(io, " fs ss panel intensity name")
+ let n = 0
+ for pk in Iterators.take(peaklist, 11)
+ if n == 10
+ # We have printed 10 already, and are here again. Truncate...
+ print(io, "\n ⋮ ⋮ ⋮ ⋮ ⋮")
+ break
+ end
+ write(io, "\n")
+ @printf(io, "%7.2f %7.2f %6i %10.2f %s",
+ pk.fs, pk.ss, pk.panelnumber, pk.intensity, pk.name)
+ n += 1
+ end
+ end
+end
+
+end # of module
diff --git a/julia/CrystFEL/src/peaksearch.jl b/julia/CrystFEL/src/peaksearch.jl
new file mode 100644
index 00000000..aaaee16c
--- /dev/null
+++ b/julia/CrystFEL/src/peaksearch.jl
@@ -0,0 +1,74 @@
+module PeakSearch
+
+import ..CrystFEL: libcrystfel
+import ..CrystFEL.Images: InternalImage
+import ..CrystFEL.PeakLists: PeakList, InternalPeakList
+
+export zaefpeaks, peakfinder8, peakfinder9
+
+
+function tf10(val)
+ if val
+ return 1
+ else
+ return 0
+ end
+end
+
+
+function zaefpeaks(image; threshold=100, mingrad=100000, minsnr=5,
+ radiusinn=4, radiusmid=5, radiusout=7, usesaturated=true)
+ out = @ccall libcrystfel.search_peaks(image.internalptr::Ptr{InternalImage},
+ threshold::Cfloat,
+ mingrad::Cfloat,
+ minsnr::Cfloat,
+ radiusinn::Cdouble,
+ radiusmid::Cdouble,
+ radiusout::Cdouble,
+ tf10(usesaturated)::Cint)::Ptr{InternalPeakList}
+ if out == C_NULL
+ throw(ErrorException("Peak search failed"))
+ end
+ PeakList(out)
+end
+
+
+
+function peakfinder8(image; threshold=100, minsnr=5, minpix=2, maxpix=200,
+ localbg=3, minres=0, maxres=5000, usesaturated=true, maxpeaks=2000)
+ out = @ccall libcrystfel.peakfinder8(image.internalptr::Ptr{InternalImage},
+ maxpeaks::Cint,
+ threshold::Cfloat,
+ minsnr::Cfloat,
+ minpix::Cint,
+ maxpix::Cint,
+ localbg::Cint,
+ minres::Cint,
+ maxres::Cint,
+ tf10(usesaturated)::Cint,
+ 0::Cint,
+ C_NULL::Ptr{Cvoid})::Ptr{InternalPeakList}
+ if out == C_NULL
+ throw(ErrorException("Peak search failed"))
+ end
+ PeakList(out)
+end
+
+
+function peakfinder9(image; minsnrbig=7, minsnrpeak=6, minsnrwhole=5, minbgsig=11,
+ brightpxcut=-Inf, window=5)
+ out = @ccall libcrystfel.search_peaks_peakfinder9(image.internalptr::Ptr{InternalImage},
+ minsnrbig::Cfloat,
+ minsnrpeak::Cfloat,
+ minsnrwhole::Cfloat,
+ minbgsig::Cfloat,
+ brightpxcut::Cfloat,
+ window::Cint)::Ptr{InternalPeakList}
+ if out == C_NULL
+ throw(ErrorException("Peak search failed"))
+ end
+ PeakList(out)
+end
+
+
+end # of module
diff --git a/julia/CrystFEL/src/reflists.jl b/julia/CrystFEL/src/reflists.jl
new file mode 100644
index 00000000..ba7311da
--- /dev/null
+++ b/julia/CrystFEL/src/reflists.jl
@@ -0,0 +1,411 @@
+module RefLists
+
+using Printf
+import ..CrystFEL: libcrystfel
+import ..CrystFEL.Symmetry: SymOpList, InternalSymOpList, symmetry_name
+export RefList, loadreflist, savereflist!
+export Reflection, UnmergedReflection, MergedReflection
+export InternalRefList
+
+
+# The internal libcrystfel structures, not exposed directly
+# We only ever have e.g. a Ptr{InternalRefList}, never a real InternalRefList
+mutable struct InternalRefList end
+mutable struct InternalReflection end
+mutable struct InternalRefListIterator end
+
+# The Julian exposed types
+abstract type Reflection end
+
+mutable struct RefList{T<:Reflection}
+ internalptr::Ptr{InternalRefList}
+ symmetry::SymOpList
+end
+
+mutable struct RefListIterator
+ lastrefl::Ptr{InternalReflection}
+ internalptr::Ptr{InternalRefListIterator}
+end
+
+mutable struct MergedReflection <: Reflection
+ internalptr::Ptr{InternalReflection}
+end
+
+mutable struct UnmergedReflection <: Reflection
+ internalptr::Ptr{InternalReflection}
+end
+
+
+function RefList{MergedReflection}(sym::SymOpList)
+ out = @ccall libcrystfel.reflist_new()::Ptr{InternalRefList}
+ if out == C_NULL
+ throw(ErrorException("Failed to create RefList"))
+ end
+ finalizer(RefList{MergedReflection}(out, sym)) do x
+ @ccall libcrystfel.reflist_free(x.internalptr::Ptr{InternalRefList})::Cvoid
+ end
+end
+
+
+function Base.iterate(reflist::RefList{T}) where T
+
+ rli = Ref{Ptr{InternalRefListIterator}}(C_NULL)
+ refl = ccall((:first_refl, libcrystfel),
+ Ptr{InternalReflection}, (Ptr{InternalRefList},Ref{Ptr{InternalRefListIterator}}),
+ reflist.internalptr, rli)
+
+ if refl == C_NULL
+ return nothing # no reflections!
+ end
+
+ iter = RefListIterator(refl,rli[])
+ finalizer(iter) do x
+ ccall((:free_reflistiterator, libcrystfel),
+ Cvoid, (Ptr{InternalRefListIterator},), x.internalptr)
+ end
+
+ return T(refl),iter
+
+end
+
+
+function Base.iterate(::RefList{T}, iter) where T
+
+ refl = ccall((:next_refl, libcrystfel),
+ Ptr{InternalReflection}, (Ptr{InternalReflection},Ptr{InternalRefListIterator}),
+ iter.lastrefl, iter.internalptr)
+
+ if refl == C_NULL
+ iter.internalptr = C_NULL # libcrystfel already freed it
+ return nothing
+ end
+
+ iter.lastrefl = refl
+
+ return T(refl),iter
+
+end
+
+
+Base.IteratorEltype(::RefList{T}) where T = T
+Base.isdone(iter::RefListIterator) = ((iter.internalptr == C_NULL) && (iter.lastrefl != C_NULL))
+Base.length(reflist::RefList) = ccall((:num_reflections, libcrystfel),
+ Cint, (Ptr{InternalRefList},), reflist.internalptr)
+
+
+function Base.getindex(reflist::RefList{T}, indices) where T
+
+ refl = @ccall libcrystfel.find_refl(reflist.internalptr::Ptr{InternalRefList},
+ indices[1]::Cint,
+ indices[2]::Cint,
+ indices[3]::Cint)::Ptr{InternalReflection}
+
+ if refl == C_NULL
+ return nothing
+ else
+ return T(refl)
+ end
+end
+
+
+function Base.get!(reflist::RefList{T}, indices) where T
+
+ refl = @ccall libcrystfel.find_refl(reflist.internalptr::Ptr{InternalRefList},
+ indices[1]::Cint,
+ indices[2]::Cint,
+ indices[3]::Cint)::Ptr{InternalReflection}
+
+ if refl == C_NULL
+ return push!(reflist, indices)
+ else
+ return T(refl)
+ end
+end
+
+
+function Base.push!(reflist::RefList{T}, indices) where T
+
+ refl = @ccall libcrystfel.add_refl(reflist.internalptr::Ptr{InternalRefList},
+ indices[1]::Cint,
+ indices[2]::Cint,
+ indices[3]::Cint)::Ptr{InternalReflection}
+
+ if refl == C_NULL
+ return nothing
+ else
+ return T(refl)
+ end
+end
+
+
+function loadreflist(filename::AbstractString)
+
+ psym = Ref{Cstring}()
+ out = ccall((:read_reflections_2, libcrystfel),
+ Ptr{InternalRefList}, (Cstring,Ref{Cstring}),
+ filename, psym)
+ if out == C_NULL
+ throw(ArgumentError("Failed to load reflection list"))
+ end
+
+ symmetryname = unsafe_string(psym[])
+ return RefList{MergedReflection}(out, SymOpList(symmetryname))
+
+end
+
+
+function savereflist!(reflist::RefList{MergedReflection}, filename::AbstractString)
+
+ r = @ccall libcrystfel.write_reflist_2(filename::Cstring,
+ reflist.internalptr::Ptr{InternalRefList},
+ reflist.symmetry.internalptr::Ptr{InternalSymOpList})::Cint
+
+ if r != 0
+ throw(ErrorException("Failed to save reflection list"))
+ end
+
+end
+
+
+function Base.show(io::IO, ::MIME"text/plain", reflist::RefList{MergedReflection})
+ println(io, "Merged reflection list in point group ", symmetry_name(reflist.symmetry))
+ print(io, " h k l intensity σ(intens) nmeas")
+ let n = 0
+ for refl in Iterators.take(reflist, 11)
+ if n == 10
+ # We have printed 10 already, and are here again. Truncate...
+ print(io, "\n ⋮ ⋮ ⋮ ⋮ ⋮ ⋮")
+ break
+ end
+ let ind = refl.indices
+ write(io, "\n")
+ @printf(io, "%4i %4i %4i %10.2f %10.2f %5i", ind[1], ind[2], ind[3],
+ refl.intensity, refl.sigintensity, refl.nmeasurements)
+ n += 1
+ end
+ end
+ end
+end
+
+
+function Base.show(io::IO, ::MIME"text/plain", reflist::RefList{UnmergedReflection})
+ println(io, "Unmerged reflection list in point group ", symmetry_name(reflist.symmetry))
+ print(io, " h k l intensity σ(intens) fs ss panel")
+ let n = 0
+ for refl in Iterators.take(reflist, 11)
+ if n == 10
+ # We have printed 10 already, and are here again. Truncate...
+ print(io, "\n ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮")
+ break
+ end
+ let ind = refl.indices,
+ pos = refl.detectorposition
+ write(io, "\n")
+ @printf(io, "%4i %4i %4i %10.2f %10.2f %7.2f %7.2f %6i",
+ ind[1], ind[2], ind[3], refl.intensity, refl.sigintensity,
+ pos.fs, pos.ss, pos.panelnumber)
+ n += 1
+ end
+ end
+ end
+end
+
+
+function detectorpos(refl::Reflection)
+ pfs = Ref{Cdouble}(0)
+ pss = Ref{Cdouble}(0)
+ pn = ccall((:get_panel_number, libcrystfel),
+ Cint, (Ptr{InternalReflection},), refl.internalptr)
+ ccall((:get_detector_pos, libcrystfel),
+ Cint, (Ptr{InternalReflection},Ref{Cdouble},Ref{Cdouble}),
+ refl.internalptr, pfs, pss)
+ (fs=pfs[], ss=pss[], panelnumber=pn)
+end
+
+
+function indices(refl::Reflection)
+ h = Ref{Cint}(0)
+ k = Ref{Cint}(0)
+ l = Ref{Cint}(0)
+ ccall((:get_indices, libcrystfel),
+ Cint, (Ptr{InternalReflection},Ref{Cint},Ref{Cint},Ref{Cint}),
+ refl.internalptr, h, k, l)
+ [h[], k[], l[]]
+end
+
+
+function symmetricindices(refl::Reflection)
+ h = Ref{Cint}(0)
+ k = Ref{Cint}(0)
+ l = Ref{Cint}(0)
+ ccall((:get_symmetric_indices, libcrystfel),
+ Cint, (Ptr{InternalReflection},Ref{Cint},Ref{Cint},Ref{Cint}),
+ refl.internalptr, h, k, l)
+ [h[], k[], l[]]
+end
+
+
+function Base.getproperty(refl::Reflection, name::Symbol)
+ if name === :intensity
+ ccall((:get_intensity, libcrystfel),
+ Cdouble, (Ptr{InternalReflection},), refl.internalptr)
+ elseif name === :sigintensity
+ ccall((:get_esd_intensity, libcrystfel),
+ Cdouble, (Ptr{InternalReflection},), refl.internalptr)
+ elseif name === :partiality
+ ccall((:get_partiality, libcrystfel),
+ Cdouble, (Ptr{InternalReflection},), refl.internalptr)
+ elseif name === :khalf
+ ccall((:get_khalf, libcrystfel),
+ Cdouble, (Ptr{InternalReflection},), refl.internalptr)
+ elseif name === :kpred
+ ccall((:get_kpred, libcrystfel),
+ Cdouble, (Ptr{InternalReflection},), refl.internalptr)
+ elseif name === :excitationerror
+ ccall((:get_exerr, libcrystfel),
+ Cdouble, (Ptr{InternalReflection},), refl.internalptr)
+ elseif name === :lorentzfactor
+ ccall((:get_lorentz, libcrystfel),
+ Cdouble, (Ptr{InternalReflection},), refl.internalptr)
+ elseif name === :phase
+ ccall((:get_phase, libcrystfel),
+ Cdouble, (Ptr{InternalReflection},), refl.internalptr)
+ elseif name === :peak
+ ccall((:get_peak, libcrystfel),
+ Cdouble, (Ptr{InternalReflection},), refl.internalptr)
+ elseif name === :meanbackground
+ ccall((:get_mean_bg, libcrystfel),
+ Cdouble, (Ptr{InternalReflection},), refl.internalptr)
+ elseif name === :temp1
+ ccall((:get_temp1, libcrystfel),
+ Cdouble, (Ptr{InternalReflection},), refl.internalptr)
+ elseif name === :temp2
+ ccall((:get_temp2, libcrystfel),
+ Cdouble, (Ptr{InternalReflection},), refl.internalptr)
+ elseif name === :nmeasurements
+ ccall((:get_redundancy, libcrystfel),
+ Cint, (Ptr{InternalReflection},), refl.internalptr)
+ elseif name === :flag
+ ccall((:get_flag, libcrystfel),
+ Cint, (Ptr{InternalReflection},), refl.internalptr)
+ elseif name === :detectorposition
+ detectorpos(refl)
+ elseif name === :indices
+ indices(refl)
+ elseif name === :symmetricindices
+ symmetricindices(refl)
+ else
+ getfield(refl, name)
+ end
+end
+
+
+function Base.setproperty!(refl::Reflection, name::Symbol, val)
+ if name === :intensity
+ @ccall libcrystfel.set_intensity(refl.internalptr::Ptr{InternalReflection},
+ val::Cdouble)::Cvoid
+ elseif name === :sigintensity
+ @ccall libcrystfel.set_esd_intensity(refl.internalptr::Ptr{InternalReflection},
+ val::Cdouble)::Cvoid
+ elseif name === :partiality
+ @ccall libcrystfel.set_partiality(refl.internalptr::Ptr{InternalReflection},
+ val::Cdouble)::Cvoid
+ elseif name === :khalf
+ @ccall libcrystfel.set_khalf(refl.internalptr::Ptr{InternalReflection},
+ val::Cdouble)::Cvoid
+ elseif name === :kpred
+ @ccall libcrystfel.set_kpred(refl.internalptr::Ptr{InternalReflection},
+ val::Cdouble)::Cvoid
+ elseif name === :excitationerror
+ @ccall libcrystfel.set_exerr(refl.internalptr::Ptr{InternalReflection},
+ val::Cdouble)::Cvoid
+ elseif name === :lorentzfactor
+ @ccall libcrystfel.set_lorentz(refl.internalptr::Ptr{InternalReflection},
+ val::Cdouble)::Cvoid
+ elseif name === :phase
+ @ccall libcrystfel.set_phase(refl.internalptr::Ptr{InternalReflection},
+ val::Cdouble)::Cvoid
+ elseif name === :peak
+ @ccall libcrystfel.set_peak(refl.internalptr::Ptr{InternalReflection},
+ val::Cdouble)::Cvoid
+ elseif name === :meanbackground
+ @ccall libcrystfel.set_mean_bg(refl.internalptr::Ptr{InternalReflection},
+ val::Cdouble)::Cvoid
+ elseif name === :temp1
+ @ccall libcrystfel.set_temp1(refl.internalptr::Ptr{InternalReflection},
+ val::Cdouble)::Cvoid
+ elseif name === :temp2
+ @ccall libcrystfel.set_temp2(refl.internalptr::Ptr{InternalReflection},
+ val::Cdouble)::Cvoid
+ elseif name === :nmeasurements
+ @ccall libcrystfel.set_redundancy(refl.internalptr::Ptr{InternalReflection},
+ val::Cint)::Cvoid
+ elseif name === :flag
+ @ccall libcrystfel.set_flag(refl.internalptr::Ptr{InternalReflection},
+ val::Cint)::Cvoid
+ elseif name === :detectorposition
+ @ccall libcrystfel.set_detector_pos(refl.internalptr::Ptr{InternalReflection},
+ val.fs::Cdouble, val.ss::Cdouble)::Cvoid
+ @ccall libcrystfel.set_panel_number(refl.internalptr::Ptr{InternalReflection},
+ val.panelnumber::Cint)::Cvoid
+ elseif name === :indices
+ throw(ErrorException("Cannot set the indices of a Reflection"))
+ elseif name === :symmetricindices
+ @ccall libcrystfel.set_symmetric_indices(refl.internalptr::Ptr{InternalReflection},
+ val[1]::Cint, val[2]::Cint, val[3]::Cint)::Cvoid
+ else
+ setfield!(refl, name, val)
+ end
+end
+
+
+function Base.propertynames(::UnmergedReflection; private=false)
+ names = (:intensity,:sigintensity,:partiality,:khalf,:kpred,:lorentzfactor,
+ :excitationerror,:phase,:peak,:meanbackground,:temp1,:temp2,
+ :nmeasurements,:flag,:detectorposition,:indices,:symmetricindices)
+ if private
+ tuple(names..., :internalptr)
+ else
+ names
+ end
+end
+
+
+function Base.propertynames(::MergedReflection; private=false)
+ names = (:intensity,:sigintensity,:phase,:temp1,:temp2,
+ :nmeasurements,:flag,:indices,:symmetricindices)
+ if private
+ tuple(names..., :internalptr)
+ else
+ names
+ end
+end
+
+
+function Base.show(io::IO, refl::MergedReflection)
+ write(io, "MergedReflection(")
+ show(io, refl.indices)
+ write(io, ", intensity=")
+ show(io, refl.intensity)
+ write(io, ", σ(intensity)=")
+ show(io, refl.sigintensity)
+ write(io, ", nmeasurements=")
+ show(io, refl.nmeasurements)
+ write(io, ")")
+end
+
+
+function Base.show(io::IO, refl::UnmergedReflection)
+ write(io, "UnmergedReflection(")
+ show(io, refl.indices)
+ write(io, " at ");
+ show(io, refl.detectorposition)
+ write(io, ", intensity=")
+ show(io, refl.intensity)
+ write(io, ", σ(intensity)=")
+ show(io, refl.sigintensity)
+ write(io, ")")
+end
+
+
+end # of module
diff --git a/julia/CrystFEL/src/streams.jl b/julia/CrystFEL/src/streams.jl
new file mode 100644
index 00000000..ce75c906
--- /dev/null
+++ b/julia/CrystFEL/src/streams.jl
@@ -0,0 +1,138 @@
+module Streams
+
+import ..CrystFEL: libcrystfel
+import ..CrystFEL.DataTemplates: DataTemplate, InternalDataTemplate
+import ..CrystFEL.Images: Image, InternalImage
+export Stream, chunkwrite, chunkread, allcrystals
+
+# Represents the real C-side (opaque) structure.
+mutable struct InternalStream end
+
+# The Julia-side structure, needed to house the pointer to the C structure
+mutable struct Stream
+ internalptr::Ptr{InternalStream}
+end
+
+
+"""
+ Stream(filename, "w", dtempl)
+
+Opens a CrystFEL stream for writing. Note that you must provide a `DataTemplate`,
+which is needed to translate "panel coordinates" to "file coordinates".
+
+Corresponds to CrystFEL C API routine `stream_open_for_write`.
+"""
+function Stream(filename, mode::AbstractString, dtempl::DataTemplate)
+
+ if mode == "w"
+ out = @ccall libcrystfel.stream_open_for_write(filename::Cstring,
+ dtempl.internalptr::Ptr{InternalDataTemplate})::Ptr{InternalStream}
+ if out == C_NULL
+ throw(ErrorException("Failed to open stream for reading"))
+ end
+
+ @ccall libcrystfel.stream_write_data_template(out::Ptr{InternalStream},
+ dtempl.internalptr::Ptr{InternalDataTemplate})::Cvoid
+
+ finalizer(close, Stream(out))
+
+ elseif mode =="r"
+ throw(ArgumentError("To open a stream for reading, don't provide the DataTemplate"))
+
+ else
+ throw(ArgumentError("Unrecognised CrystFEL stream mode"*mode))
+ end
+end
+
+
+"""
+ Stream(filename, "r")
+
+Opens a CrystFEL stream for reading.
+
+Close the stream with `close` when you've finished (this will happen
+automatically when the `Stream` object is finalized).
+
+Corresponds to CrystFEL C API routine `stream_open_for_read`.
+"""
+function Stream(filename, mode::AbstractString)
+
+ if mode == "r"
+ out = @ccall libcrystfel.stream_open_for_read(filename::Cstring)::Ptr{InternalStream}
+ if out == C_NULL
+ throw(ErrorException("Failed to open stream for reading"))
+ end
+ finalizer(close, Stream(out))
+
+ elseif mode == "w"
+ throw(ArgumentError("To open a stream for writing, you must provide "
+ *"a DataTemplate: use Stream(filename, \"w\", dtempl)"))
+
+ else
+ throw(ArgumentError("Unrecognised CrystFEL stream mode"*mode))
+ end
+end
+
+
+function Base.close(st::Stream)
+ if st.internalptr != C_NULL
+ @ccall libcrystfel.stream_close(st.internalptr::Ptr{InternalStream})::Cvoid
+ st.internalptr = C_NULL
+ end
+end
+
+
+function streamflags(peaks, reflections, imagedata)
+ flags = 0
+ if reflections
+ flags |= 2
+ end
+ if peaks
+ flags |= 4
+ end
+ if imagedata
+ flags |= 8
+ end
+ return flags
+end
+
+
+function chunkwrite(st::Stream, image::Image; peaks=true, reflections=true)
+ st.internalptr == C_NULL && throw(ErrorException("Stream is closed"))
+ flags = streamflags(peaks, reflections, false)
+ @ccall libcrystfel.stream_write_chunk(st.internalptr::Ptr{InternalStream},
+ image.internalptr::Ptr{InternalImage},
+ flags::Cint)::Cvoid
+end
+
+
+function chunkread(st::Stream; peaks=true, reflections=true, datageom=true)
+
+ st.internalptr == C_NULL && throw(ErrorException("Stream is closed"))
+
+ flags = streamflags(peaks, reflections, datageom)
+ out = @ccall libcrystfel.stream_read_chunk(st.internalptr::Ptr{InternalStream},
+ flags::Cint)::Ptr{InternalImage}
+ out == C_NULL && return nothing
+
+ finalizer(Image(out, nothing, [], [])) do x
+ ccall((:image_free, libcrystfel), Cvoid, (Ptr{InternalImage},), x.internalptr)
+ end
+
+end
+
+
+function allcrystals(st)
+ Channel() do ch
+ while true
+ image = chunkread(st, peaks=false, reflections=true, datageom=false)
+ image === nothing && break
+ for cr in image.crystals
+ put!(ch, (cr.crystal, cr.reflections))
+ end
+ end
+ end
+end
+
+
+end # of module
diff --git a/julia/CrystFEL/src/symmetry.jl b/julia/CrystFEL/src/symmetry.jl
new file mode 100644
index 00000000..4b4caa4d
--- /dev/null
+++ b/julia/CrystFEL/src/symmetry.jl
@@ -0,0 +1,134 @@
+module Symmetry
+
+import ..CrystFEL: libcrystfel
+export SymOpList, InternalSymOpList, InternalIntegerMatrix
+export symmetry_name, asymmetricindices
+
+
+# Types for pointers returned by libcrystfel
+mutable struct InternalSymOpList end
+mutable struct InternalIntegerMatrix end
+
+
+# Exposed types
+mutable struct SymOpList
+ internalptr::Ptr{InternalSymOpList}
+end
+
+
+mutable struct SymOp
+ internalptr::Ptr{InternalIntegerMatrix}
+end
+
+
+function SymOpList(pointgroup::AbstractString)
+
+ out = ccall((:get_pointgroup, libcrystfel),
+ Ptr{InternalSymOpList}, (Cstring,), pointgroup)
+ if out == C_NULL
+ throw(OutOfMemoryError())
+ end
+
+ sym = SymOpList(out)
+
+ finalizer(sym) do x
+ ccall((:free_symoplist, libcrystfel),
+ Cvoid, (Ptr{InternalSymOpList},), x.internalptr)
+ end
+
+ return sym
+
+end
+
+
+function Base.getindex(sym::SymOpList, i::Int)
+
+ if (i<1) || (i > length(sym))
+ throw(BoundsError())
+ end
+
+ out = ccall((:get_symop, libcrystfel),
+ Ptr{InternalIntegerMatrix},
+ (Ptr{InternalSymOpList},Ptr{Cvoid},Cint),
+ sym.internalptr,C_NULL,i-1)
+
+ if out == C_NULL
+ throw(OutOfMemoryError())
+ end
+
+ return SymOp(out)
+
+end
+
+
+function Base.length(sym::SymOpList)
+ return ccall((:num_equivs, libcrystfel),
+ Cint, (Ptr{InternalSymOpList},Ptr{Cvoid}),
+ sym.internalptr, C_NULL)
+end
+
+
+function hkl_op(op::SymOp)
+ s = ccall((:name_equiv, libcrystfel),
+ Cstring,
+ (Ptr{InternalIntegerMatrix},),
+ op.internalptr)
+ return unsafe_string(s)
+end
+
+
+function symmetry_name(sym::SymOpList)
+ s = ccall((:symmetry_name, libcrystfel),
+ Cstring,
+ (Ptr{InternalSymOpList},),
+ sym.internalptr)
+ return unsafe_string(s)
+end
+
+
+function Base.iterate(sym::SymOpList)
+ return (sym[1], 2)
+end
+
+
+function Base.iterate(sym::SymOpList, i)
+ if i > length(sym)
+ return nothing
+ else
+ return (sym[i], i+1)
+ end
+end
+
+
+function asymmetricindices(sym::SymOpList, indices)
+ ho = Ref{Cint}(0)
+ ko = Ref{Cint}(0)
+ lo = Ref{Cint}(0)
+ @ccall libcrystfel.get_asymm(sym.internalptr::Ptr{InternalSymOpList},
+ indices[1]::Cint, indices[2]::Cint, indices[3]::Cint,
+ ho::Ref{Cint}, ko::Ref{Cint}, lo::Ref{Cint})::Cvoid
+ (ho[], ko[], lo[])
+end
+
+
+function Base.show(io::IO, ::MIME"text/plain", sym::SymOpList)
+ print(io, length(sym), "-element SymOpList (\"", symmetry_name(sym), "\")")
+ for op in sym
+ print(io, "\n", hkl_op(op))
+ end
+end
+
+
+function Base.show(io::IO, sym::SymOpList)
+ write(io, "SymOpList(\"", symmetry_name(sym), "\")")
+end
+
+
+function Base.show(io::IO, op::SymOp)
+ write(io, "SymOp(")
+ write(io, "\"")
+ write(io, hkl_op(op))
+ write(io, "\")")
+end
+
+end # of module
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-01 09:18:30 +0000