Questions about the sampling function

[MohamedNedal]

Hello, I was checking the function sample_at_coords and I have a couple of questions about this part here.

1. What should be the heliospheric frame for the input coordinates (lon, lat, distance)?
2. The input coordinates (lon, lat, distance), should their units be (deg, deg, Rsun) or (rad, rad, Rsun)?

[thepeteriley]

Hi Mohamed,
The frame is Carrington coordinates. This is how the model output is produced (if you're using a different routine to pull the in situ data, you'll need to make sure they're in the same frame). The units are as follows:

(lon, lat, radius) in (deg, deg, m)

However, bear in mind that the raw units of the HDF files are in radians and solar radii. Depending on how you're interacting with the model results, it's a good idea to check them. But if you're using the PsiPy routines and variables like

psp_coords.lon, psp_coords.lat, psp_coords.radius

then these are in the aforementioned units (deg, deg, m).

You can check this all out by running plot_in_situ_comparison.py, which can be downloaded from:

https://psipy.readthedocs.io/en/stable/auto_examples/sampling/plot_in_situ_comparison.html#sphx-glr-auto-examples-sampling-plot-in-situ-comparison-py

It already includes a print statement with the units described.

Let me know if you need more clarification. We should probably add a sentence or two in the ReadMeDoc to explain this :)

[MohamedNedal]

Thank you, Dr. Riley @thepeteriley , I appreciate that :)
I would like to sample density from the MAS datacube along a trajectory defined by these coordinates:
(x, y) --> (<Longitude -1529.53270443 arcsec>, <Latitude -755.9071465 arcsec>)
and z is defined as a 1D array Z = np.linspace(-3, 3, 50)*u.solRad
I got a density profile and I need confirmation that it's correct/incorrect. Here's what I did:

# extract MAS density data 
mas_output = MASOutput('/path/to/mas/data')
rho = mas_output['rho']

# solar radius in arcsec 
Rsun_arcsec = 959.798683

# define lon and lat of a point 
x = -1529.53*u.arcsec
y = -755.907*u.arcsec

# convert the lon and lat to m unit 
X = (x.value/Rsun_arcsec)*u.solRad
Y = (y.value/Rsun_arcsec)*u.solRad
X = X.to('m')
Y = Y.to('m')

# make a 1D array of points on the z-axis to sample the data at 
npoints = 50
Z = np.linspace(-3, 3, npoints)*u.solRad
Z = Z.to('m')

mas_sampled_rho = []
# go over every value in Z array and use it with the aformentioned (X,Y) coords 
for z in Z:
    # make a skycoord object with heliocentric frame 
    c_heliocentric = SkyCoord(X, Y, z, frame=frames.Heliocentric, obstime='2019-04-03T12:00:04.840', observer='Earth')

    # transform to heliospheric Carrington frame 
    c_heliocarrington = c_heliocentric.transform_to(frames.HeliographicCarrington)
    
    # sample that point from the MAS MHD datacube 
    sample = rho.sample_at_coords(c_heliocarrington.lon, c_heliocarrington.lat, c_heliocarrington.radius)
    mas_sampled_rho.append(sample)

# plot the sampled density profile 
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(Z.to('solRad'), mas_sampled_rho, label='MAS')
ax.set_xlabel('Distance [$R_\odot$]')
ax.set_ylabel(r'$\rho$ [cm$^{-3}$]')
ax.set_yscale('log')
plt.show()

This is the output density profile:

I tried to visualize that trajectory and this is what I got:

# init a 3D plot 
fig = plt.figure(figsize=(7,7))
ax = fig.add_subplot(111, projection='3d')

# draw the Sun 
phi_sun, theta_sun = np.mgrid[0:2*np.pi:100j, 0:np.pi:100j]
x_sun = np.sin(theta_sun)*np.cos(phi_sun)
y_sun = np.sin(theta_sun)*np.sin(phi_sun)
z_sun = np.cos(theta_sun)
sun_surf = ax.plot_surface(x_sun, y_sun, z_sun, color='orange', alpha=1, label='Sun')
sun_surf._facecolors2d = sun_surf._facecolor3d
sun_surf._edgecolors2d = sun_surf._edgecolor3d

# draw a big sphere that represents the inner part of the MAS domain 
phi_mas, theta_mas = np.mgrid[0:2*np.pi:200j, 0:np.pi:200j]
big_radius = 6
x_mas = np.sin(theta_mas)*np.cos(phi_mas)*big_radius
y_mas = np.sin(theta_mas)*np.sin(phi_mas)*big_radius
z_mas = np.cos(theta_mas)*big_radius
mas_surf = ax.plot_surface(x_mas, y_mas, z_mas, color='dodgerblue', alpha=0.1, label='MAS')
mas_surf._facecolors2d = mas_surf._facecolor3d
mas_surf._edgecolors2d = mas_surf._edgecolor3d

# draw Sun-Earth LOS 
x_pos = 2
y_pos = 0
z_pos = 0
x_direct = 1
y_direct = 0
z_direct = 0
ax.quiver(x_pos, y_pos, z_pos, x_direct, y_direct, z_direct, length=2, arrow_length_ratio=0.3, color='black')
ax.text(x_pos, y_pos, z_pos, 'Earth LOS', size=10, zorder=1, color='black')

# draw coord axes at the origin point -- the center of the Sun 
limit = 4
# // to Y 
ax.plot([0, 0], [-limit, limit], [0, 0], '--', color='gray')
# // to X 
ax.plot([-limit, limit], [0, 0], [0, 0], '--', color='gray')
# // to Z 
ax.plot([0, 0], [0, 0], [-limit, limit], '--', color='gray')

# define lon and lat of a point 
x = -1529.53*u.arcsec
y = -755.907*u.arcsec

# convert the lon and lat to Rsun unit 
x = (x.value/Rsun_arcsec)*u.solRad
y = (y.value/Rsun_arcsec)*u.solRad

# make a 1D array of points on the z-axis to sample the data at 
npoints = 50
Z = np.linspace(-3, 3, npoints)*u.solRad

# draw a line passing through the point across the MAS domain 
line_label = 'Sampling path'
rad_sample, theta_sample, phi_sample = [], [], []
mas_dist, mas_lon, mas_lat = [], [], []

for z in Z:
    # Convert the point from cartesian to spherical coordinates 
    r = np.sqrt(x.value**2 + y.value**2 + z.value**2)
    theta = np.arctan2(np.sqrt(x.value**2 + y.value**2), z.value)
    phi = np.arctan2(y.value, x.value)

    # Plot the point 
    ax.scatter(r*np.cos(theta), r*np.sin(theta)*np.cos(phi), r*np.sin(theta)*np.sin(phi), 
               c='crimson', marker='o', s=7, label=line_label)
    line_label = '_nolegend_'
    
    rad_sample.append(r)
    theta_sample.append(theta)
    phi_sample.append(phi)
    
    # dist --> X = r*np.cos(theta) 
    # lon  --> Y = r*np.sin(theta)*np.cos(phi) 
    # lat  --> Z = r*np.sin(theta)*np.sin(phi)
    mas_dist.append(r*np.cos(theta))
    mas_lon.append(r*np.sin(theta)*np.cos(phi))
    mas_lat.append(r*np.sin(theta)*np.sin(phi))

# plot settings 
ax.set_xlim3d(left=-limit, right=limit)
ax.set_ylim3d(bottom=-limit, top=limit)
ax.set_zlim3d(bottom=-limit, top=limit)
ax.set_box_aspect([limit, limit, limit])
for axis in 'xyz':
    getattr(ax, 'set_{}label'.format(axis))(axis+' [R$_\odot]$')
ax.legend(loc='upper left', prop={'size':10}, bbox_to_anchor=(-0.15, 1.1))
ax.grid(False)
plt.show()