Skip to content
This repository has been archived by the owner on Aug 6, 2024. It is now read-only.

[BUG] consistent offset and lobes in residual #30

Open
4 tasks done
nstarman opened this issue Jan 22, 2021 · 0 comments
Open
4 tasks done

[BUG] consistent offset and lobes in residual #30

nstarman opened this issue Jan 22, 2021 · 0 comments
Assignees
Labels
bug Something isn't working
Milestone

Comments

@nstarman
Copy link
Collaborator

Description

Figure out why and solve!

Unknown-6

Expected behavior

  • dependent on l, m max
  • dependent on the grid

Actual behavior

  • dependent on l, m max : it's not
  • dependent on the grid : maybe it is?

Steps to Reproduce

# THIRD PARTY
import agama
import astropy.coordinates as coord
import astropy.units as u
import numpy as np
from astropy.visualization import quantity_support
from galpy import df as gdf
from galpy import potential as gpot

# PROJECT-SPECIFIC
from discO import (
    GaussianMeasurementErrorSampler,
    PotentialFitter,
    PotentialSampler,
    conf,
)
from discO.plugin.agama.fitter import (
    AGAMAMultipolePotentialFitter,
    AGAMAPotentialFitter,
)
from discO.core.pipeline import Pipeline
from discO.core.residual import GridResidual

quantity_support();

mass = 1e12 * u.solMass
r0 = 10 * u.kpc  # scale factor

agama.setUnits(mass=1, length=1, velocity=1)

hernquist_pot = gpot.HernquistPotential(amp=2 * mass, a=r0)
sampler = PotentialSampler(gdf.isotropicHernquistdf(hernquist_pot))
measurer = GaussianMeasurementErrorSampler(c_err=1*u.percent)  # 10% error
fitter = PotentialFitter(None, key="agama", pot_type="multipole", symmetry="s", gridsizeR=50, lmax=0, mmax=0)

_x = np.linspace(-5, 5, num=50)
x, y, z = np.meshgrid(_x, _x, _x) * u.kpc
points = coord.CartesianRepresentation(x, y, z)
residualer = GridResidual(points, original_pot=hernquist_pot)

pipe = Pipeline(sampler=sampler, fitter=fitter, residualer=residualer)
pipe

resid = pipe(n=10000, frame=coord.Galactocentric(), random=0)


from mpl_toolkits.mplot3d import axes3d
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np

fig = plt.figure()
ax = fig.gca(projection='3d')

rar = np.abs(resid.represent_as(coord.SphericalRepresentation).vf_distance.flatten().value)
norm = mpl.colors.Normalize(vmin=min(rar), vmax=max(rar), clip=True)
mapper = cm.ScalarMappable(norm=norm, cmap=cm.inferno_r)

uvw = resid.vf_xyz.value * 1e13

qv = ax.quiver(
    x, y, z, *uvw,
    color=mapper.to_rgba(rar, alpha=0.1),
    length=0.1, )
ax.set_xlim(-8, 8)
ax.set_ylim(-8, 8)
ax.set_zlim(-8, 8)

plt.colorbar(qv)

plt.show()

System Details

macOS-10.16-x86_64-i386-64bit
Python 3.8.2 | packaged by conda-forge | (default, Apr 24 2020, 07:56:27)
[Clang 9.0.1 ]
astropy 4.3.dev310+g59d625df4.d20201220
astroquery 0.4.2.dev6458
Matplotlib 3.2.1
Numpy 1.18.4
Scipy 1.4.1

Checklist

  • Check out our contributing guidelines.
  • Check out our code of conduct.
  • Search the GitHub repository to see if a similar issue has already been posted. If a similar issue is closed, have a quick look to see if you are satisfied by the resolution. If not please go ahead and open an issue!
  • Please check that the development version still produces the same bug. You can install development version with pip install git+https://github.com/astropy/astropy command.
@nstarman nstarman added the bug Something isn't working label Jan 22, 2021
@nstarman nstarman added this to the v0.1 milestone Jan 22, 2021
@nstarman nstarman self-assigned this Jan 22, 2021
Sign up for free to subscribe to this conversation on GitHub. Already have an account? Sign in.
Labels
bug Something isn't working
Projects
None yet
Development

No branches or pull requests

2 participants