BUG: forward() does not give same answer as first solution in forward_solutions_table()

[prjemian]

During testing to hold a rotation axis constant when computing the forward(h,k,l) axis positions, observed that the position returned by forward((1,1,1)) is not the same as the first position shown by forward_solutions_table([(1,1,1)], full=True). The problem observed is that phi is not held to the expected value by forward().

code

e4cv.reset_constraints()

print(f"{e4cv.calc.engine.modes=}")
e4cv.calc.engine.mode = "constant_phi"

e4cv.apply_constraints({"phi": Constraint(-180, 180, 21.3456, True)})
e4cv.show_constraints()
print(e4cv.forward_solutions_table([refl], full=True, digits=4))
print(e4cv.forward(refl))
e4cv.reset_constraints()

output

e4cv.calc.engine.modes=['bissector', 'constant_omega', 'constant_chi', 'constant_phi', 'double_diffraction', 'psi_constant']
===== ========= ========== ======= ====
axis  low_limit high_limit value   fit 
===== ========= ========== ======= ====
omega -180.0    180.0      0.0     True
chi   -180.0    180.0      0.0     True
phi   -180.0    180.0      21.3456 True
tth   -180.0    180.0      0.0     True
===== ========= ========== ======= ====

self._solutions=[PosCalcE4CV(omega=4.92922853140664, chi=37.66038381906326, phi=21.3456, tth=-28.430681927403345), PosCalcE4CV(omega=-4.92922853140664, chi=-142.33961618093676, phi=21.3456, tth=28.430681927403345)]
========= ======== ======= ========= ======= ========
(hkl)     solution omega   chi       phi     tth     
========= ======== ======= ========= ======= ========
(1, 1, 1) 0        4.9292  37.6604   21.3456 -28.4307
(1, 1, 1) 1        -4.9292 -142.3396 21.3456 28.4307 
========= ======== ======= ========= ======= ========

self._solutions=[PosCalcE4CV(omega=21.07167613962381, chi=45.00000223615021, phi=0.0, tth=-28.430681910610566), PosCalcE4CV(omega=-21.07167613962381, chi=-134.99999776384982, phi=0.0, tth=28.430681910610566)]
self._solutions=[PosCalcE4CV(omega=21.071676144724915, chi=45.000002234645585, phi=0.0, tth=-28.4306819088791), PosCalcE4CV(omega=-21.071676144724915, chi=-134.99999776535444, phi=0.0, tth=28.4306819088791)]
PosCalcE4CV(omega=21.071676144724915, chi=45.000002234645585, phi=0.0, tth=-28.4306819088791)

Note, a diagnostic print was added just before this line: https://github.com/bluesky/hklpy/blob/bb4ac84a4ee7747200d342c56da0719d09ae7c1c/hkl/engine.py#L258

Thanks, @strempfer, for identifying this situation.

[strempfer]

A general freeze command that freezes the phi in constant_phi mode and so on but also uses the psi engine for calculation of psi in constant psi mode would be great, similar to:

def freeze(*args):
    _geom_ = current_diffractometer()
    _check_geom_selected()
    mode = _geom_.calc.engine.mode
    if "constant" in mode.strip(" "):
        print(mode.strip(""))
        if "phi" in mode.strip(""):
            axis = "phi"
            value = _geom_.get_axis_constraints(axis).value
        elif "chi" in mode.strip(""):
            axis = "chi"
            value = _geom_.get_axis_constraints(axis).value
        elif "omega" in mode.strip(""):
            axis = "omega"
            value = _geom_.get_axis_constraints(axis).value
        elif "psi" in mode.strip(""):
            axis = "psi"
            (
                _h2,
                _k2,
                _l2,
                value,
            ) = _geom_.calc._engine.engine.parameters_values_get(
                Hkl.UnitEnum.USER
            )

        if len(args) == 0:
            value = int(
                (input("Freeze {} to ({})? ".format(axis, value))) or value
            )
        elif len(args) == 1:
            value = args[0]
        else:
            raise ValueError(
                "either no argument or angle needs to be provided."
            )
        if axis == "psi":
            _geom_.calc._engine.engine.parameters_values_set(
                [_h2, _k2, _l2, value], Hkl.UnitEnum.USER
            )
        else:
            _geom_.apply_constraints({axis: Constraint(0, 120, value, True)})
        print("{} frozen to {}".format(axis, value))
    else:
        raise ValueError("Function not available for mode '{}'".format(mode))

[strempfer]

It would be nice to have a function that returns for a given mode if there is an axis kept constant and the name of the axis. This would simplify the freeze function considerably.

[prjemian]

@strempfer - Like this:

In [1]: from hkl import SimulatedE6C

In [2]: sim6c = SimulatedE6C("", name="sim6c")

In [3]: [ax for ax in sim6c.calc.engine._engine.axis_names_get(0) if ax not in sim6c.calc.engine._engine.axis_names_get(1)]
Out[3]: ['mu', 'gamma']

In [4]: sim6c.calc.engine.mode = "constant_chi_vertical"

In [5]: [ax for ax in sim6c.calc.engine._engine.axis_names_get(0) if ax not in sim6c.calc.engine._engine.axis_names_get(1)]
Out[5]: ['mu', 'chi', 'gamma']

Thinking that such a function should also show what additional parameters are required by the mode (instead of making yet another function to remember).

[strempfer]

Perfect! This allows also to determine if in vertical (fixed mu and gamma) or horizontal (fixed omega and delta) mode.

From: Pete R Jemian @.> Date: Thursday, November 30, 2023 at 4:57 PM To: bluesky/hklpy @.> Cc: Strempfer, Joerg @.>, Mention @.> Subject: Re: [bluesky/hklpy] BUG: forward() does not give same answer as first solution in forward_solutions_table() (Issue #303)

@strempferhttps://github.com/strempfer - Like this:

In [1]: from hkl import SimulatedE6C

In [2]: sim6c = SimulatedE6C("", name="sim6c")

In [3]: [ax for ax in sim6c.calc.engine._engine.axis_names_get(0) if ax not in sim6c.calc.engine._engine.axis_names_get(1)]

Out[3]: ['mu', 'gamma']

In [4]: sim6c.calc.engine.mode = "constant_chi_vertical"

In [5]: [ax for ax in sim6c.calc.engine._engine.axis_names_get(0) if ax not in sim6c.calc.engine._engine.axis_names_get(1)]

Out[5]: ['mu', 'chi', 'gamma']

Thinking that such a function should also show what additional parameters are required by the mode (instead of making yet another function to remember).

— Reply to this email directly, view it on GitHubhttps://github.com/bluesky/hklpy/issues/303#issuecomment-1834688690, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AOAHWJ7AETNZWVKCTKIHVMLYHEFM3AVCNFSM6AAAAAA7TKZ5ZWVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTQMZUGY4DQNRZGA. You are receiving this because you were mentioned.Message ID: @.***>

[prjemian]

First, though, let's resolve why the two solutions are different.