Bug with Hamilton Potential and Frame Units

[mlucey]

Hey Adrian,

I just updated Gala so that I could use the orbit.to_galpy_orbit() function, but now I have a weird bug popping up. I am trying to rotate an SCF potential computed from an N-body snapshot, but I get this VaueError:

File "gala/potential/hamiltonian/chamiltonian.pyx", line 62, in gala.potential.hamiltonian.chamiltonian.Hamiltonian.init ValueError: Potential and Frame must have compatible unit systems (UnitSystem (kpc, Myr, solMass, rad) vs UnitSystem (kpc, Myr, solMass, rad))

The units are equivalent (at least according to the printed message) so I'm not sure why it thinks frame.units != potential.units.

As a quick hacky fix I tried to just go in to the source code and comment out the if --> raise error statements but I also have no idea how Cython works so my quick fix is not working for me. If you don't want to deal with this weird unit mismatch thing, would you possibly be able to just slack me the correct way to recompile with the lines commented out for now?

Thanks so much! Gala is my hero!

[adrn]

Weird! Can you share the code (here in a comment) where you instantiate the SCFPotential and Hamiltonian?

[mlucey]

Yeah, below are the snippets. Maybe it's something weird because I am loading a saved SCFpotential?

`from gala.units import galactic

S,T = scf.compute_coeffs_discrete(xyz, mass=mass*u.Msun/totalmass, nmax=8, lmax=8, r_s=1.,skip_m=False) SCFpotential = scf.SCFPotential(Snlm=S, Tnlm=T, m=totalmass, r_s=1.,units=galactic) SCFpotential.save('/home/mrl2968/Desktop/Bar_Dynamics/New390_pot.yml')

pot = load('/home/mrl2968/Desktop/Bar_Dynamics/New390_pot.yml',module=scf)

Om_bar = omega u.km/u.s/u.kpc frame = gp.ConstantRotatingFrame(Omega=[0,0,Om_bar.value]Om_bar.unit, units=galactic) H = gp.Hamiltonian(potential=pot, frame=frame) o=H.integrate_orbit(c,dt=1*u.Myr, n_steps=1000)`

[adrn]

That all looks right, so it could be a bug with the way unit systems are compared - oops!

Here's a hack that might fix it: Change the frame definition line to

frame = gp.ConstantRotatingFrame(Omega=[0,0,Om_bar.value] * Om_bar.unit,
                                 units=pot.units)

[mlucey]

That worked, thanks so much!