Scheduled maintenance on Monday 2019-06-24 between 10:00-11:00 CEST

Commit 3ac39d77 authored by Lukas Platz's avatar Lukas Platz

added Point Source model based on MfCorrelatedField

parent e64541c3
......@@ -82,7 +82,8 @@ from .library.dynamic_operator import (dynamic_operator,
from .library.light_cone_operator import LightConeOperator
from .library.wiener_filter_curvature import WienerFilterCurvature
from .library.correlated_fields import CorrelatedField, MfCorrelatedField
from .library.correlated_fields import (CorrelatedField, MfCorrelatedField,
from .library.adjust_variances import (make_adjust_variances_hamiltonian,
......@@ -141,3 +141,54 @@ def MfCorrelatedField(target, amplitudes, name='xi'):
vol = reduce(mul, [sp.scalar_dvol**-0.5 for sp in hsp])
return ht(vol*A*ducktape(hsp, None, name))
def MfPartiallyCorrelatedField(target, energy_amplitude, name='xi_u'):
"""Constructs an operator which turns white Gaussian excitation fields
into a correlated field defined on a DomainTuple with two entries.
On the first domain, a white correlation structure is assumed. On
the second domain, the correlation structures given by energy_amplitude
is used.
This operator may be used as a model for multi-frequency reconstructions
with correlation structure only in the energy direction.
target : Domain, DomainTuple or tuple of Domain
Target of the operator. Must contain exactly two spaces.
It is assumed that the second space is the energy domain.
energy_amplitude: Operator
amplitude operator for the energy correlation structure
name : string
:class:`MultiField` key for xi-field.
Correlated field
In NIFTy, non-harmonic RGSpaces are by definition periodic. Therefore
the operator constructed by this method will output a correlated field
with *periodic* boundary conditions. If a non-periodic field is needed,
one needs to combine this operator with a :class:`FieldZeroPadder` or even
two (one for the energy and one for the spatial subdomain)
tgt = DomainTuple.make(target)
if len(tgt) != 2:
raise ValueError
h_space = DomainTuple.make([dom.get_default_codomain() for dom in tgt])
ht1 = HarmonicTransformOperator(h_space, target=tgt[0], space=0)
ht2 = HarmonicTransformOperator(, target=tgt[1], space=1)
ht = ht2 @ ht1
p_space =[0]
power_distributor = PowerDistributor(h_space[-1], p_space)
A = power_distributor(energy_amplitude)
dd = ContractionOperator(h_space, 0).adjoint
return ht((dd @ A)*ducktape(h_space, None, name))
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment