Hi all,

I’m working on some updates to my yt-based pyXSIM package and I’m stuck and I thought I would email this list for some guidance. 

I’m thinking of a situation where we have a data object (sphere, box, whatever) and it straddles a periodic boundary. I want to convert the coordinates such that the coordinates are translated with respect to some arbitrary origin (say, the center of a sphere but in theory it could be anywhere), but are also continuous, i.e, they do not wrap at the boundary. 

I’ve looked at the functions get_radius and get_periodic_rvec in yt/fields/field_functions.py, and based on that I have come up with the code below, but it doesn’t quite work for any arbitrary value of the “ctr” argument (the origin). I think it’s because the functions I mentioned only need to calculate absolute value differences in order to compute a radius. 

The x, y, and z arguments to the function below are the input coordinate arrays. Could anyone who is more familiar with this sort of thing point out what I should be doing? 

Best,

John Z

def get_periodic_coords(ds, ctr, x, y, z):
    coords = ds.arr(np.zeros((3, x.size)), "kpc")
    coords[0, :] = (x - ctr[0]).to("kpc")
    coords[1, :] = (y - ctr[1]).to("kpc")
    coords[2, :] = (z - ctr[2]).to("kpc")
    if sum(ds.periodicity) == 0:
        return coords
    dw = ds.domain_width.to("kpc")
    for i in range(3):
        if ds.periodicity[i]: 
            c = coords[i, :]
            cdw = c - dw[i]
            mins = np.argmin([np.abs(c), np.abs(cdw)], axis=0)
            ii = np.where(mins == 1)[0]
            coords[i, ii] = coords[i, ii] - dw[i]
    return coords