<jupyter_start><jupyter_text>Project to tangent plane<jupyter_code>galcen_frame = coord.Galactocentric(z_sun=0*u.kpc, galcen_distance=8.1*u.kpc)
sim_c = stream_reg.to_coord_frame(coord.ICRS, galactocentric_frame=galcen_frame)
ra_reg =  sim_c.ra
dec_reg = sim_c.dec
from astropy.wcs import WCS

def get_projected_coords(c, pal5_c):
    """
    Parameters
    ----------
    c : astropy SkyCoord object
        The coordinates that you want to project
    pal5_c : astropy SkyCoord object
        The coordinates of Pal 5 the cluster
    """
    offset_fr = coord.SkyOffsetFrame(origin=pal5_c)
    c2 = c.transform_to(offset_fr)

    wcs = WCS(naxis=2)
    wcs.wcs.crpix = [0., 0.]
    wcs.wcs.cdelt = [1., 1.]
    wcs.wcs.crval = [0, 0]
    wcs.wcs.ctype = ["RA---ZEA", "DEC--ZEA"]

    return wcs.wcs_world2pix(c2.lon.degree, c2.lat.degree, 0)
stream_proj = get_projected_coords(sim_c, pal5_c)
stream_proj_x = stream_proj[0]
stream_proj_y = stream_proj[1]
fig,axes = plt.subplots(1,2,figsize=(14,7))#, sharey=True, sharex=True)
matplotlib.rcParams.update({'font.size': 22})

axes[0].set_aspect('equal')
axes[1].set_aspect('equal')
axes[0].set_xlim(300,200)
axes[1].set_xlim(80,-20)
axes[0].set_ylim(-80,20)
axes[1].set_ylim(-80,20)

axes[0].set_title('SkyCords', fontsize=18)
axes[0].scatter(ra_reg, dec_reg, s=1)
axes[0].set_xlabel('RA [deg]')
axes[0].set_ylabel('dec [deg]')

axes[1].set_title('Area perserving projection', fontsize=18)
axes[1].scatter(stream_proj_x,stream_proj_y, s=1)
axes[1].set_xlabel('$X$ [deg]')
axes[1].set_ylabel('$Y$ [deg]')

fig.tight_layout()<jupyter_output><empty_output><jupyter_text>Fit polynomial to leading and trailing armUsing new are perserving coordinate system<jupyter_code>#First sort data pased on "X" coord. 

stream_proj_x_sort = np.sort(stream_proj_x)
#stream_proj_y_sort = np.sort(stream_proj_y)
inds = stream_proj_x.argsort()
stream_proj_y_sort = stream_proj_y[inds]
#Gotta sort sim data. and separate into leading trailing particles based on prog. position 

sim_lead_index = np.where((stream_proj_y_sort < -1) & (stream_proj_x_sort < 20)& (stream_proj_y_sort > -40))
sim_trail_index = np.where((stream_proj_y_sort > -1))


polyfit_lead = np.polyfit(stream_proj_x_sort[sim_lead_index] , stream_proj_y_sort[sim_lead_index] , deg=3)
polyfit_trail = np.polyfit(stream_proj_x_sort[sim_trail_index] , stream_proj_y_sort[sim_trail_index] , deg=3)
stream_poly_y_lead = stream_proj_x_sort[sim_lead_index]**3*polyfit_lead[0] + stream_proj_x_sort[sim_lead_index]**2*polyfit_lead[1] + stream_proj_x_sort[sim_lead_index]*polyfit_lead[2] + polyfit_lead[3] 
stream_poly_y_trail = stream_proj_x_sort[sim_trail_index]**3*polyfit_trail[0] + stream_proj_x_sort[sim_trail_index]**2*polyfit_trail[1] + stream_proj_x_sort[sim_trail_index]*polyfit_trail[2] + polyfit_trail[3] 
#Have x-values for poly fit that are not the sim. points
#poly_x_lead = np.linspace(-10,0,len(stream_proj_x_sort[sim_lead_index]))
#poly_x_trail = np.linspace(0,60,len(stream_proj_x_sort[sim_trail_index]))
#poly_x_trail.shape 
fig,axes = plt.subplots(1,1,figsize=(7,7))#, sharey=True, sharex=True)
matplotlib.rcParams.update({'font.size': 22})


axes.set_aspect('equal')
axes.set_xlim(80,-20)
axes.set_ylim(-80,20)

axes.set_title('Area perserving projection', fontsize=18)
axes.scatter(stream_proj_x_sort[sim_lead_index],stream_proj_y_sort[sim_lead_index], s=1, color='orange')
axes.scatter(stream_proj_x_sort[sim_trail_index],stream_proj_y_sort[sim_trail_index], s=1, color='steelblue')

axes.plot(stream_proj_x_sort[sim_lead_index],stream_poly_y_lead, color="maroon") #poly fit to leading arm
axes.plot(stream_proj_x_sort[sim_trail_index],stream_poly_y_trail, color="black") #poly fit to trailing arm

axes.set_xlabel('$X$ [deg]')
axes.set_ylabel('$Y$ [deg]')

fig.tight_layout()<jupyter_output><empty_output><jupyter_text>Find normal,v, and tangential, u, vectors to points along polynomialFor u-vector (tangent)For any given point,k, along the polynomial, use the values of polynomial surrounding k to find the tangent to the point (i.e. derivative).For v-vector (normal)For any given point,k, along the polynomial, use the values of polynomial surrounding k to find the curvature in this point (double derivative) to know the "sign" of curvature and therefore how to place v.<jupyter_code>#Looking at one point, k, along the trailing poly-fit
k_neigh = 1 #look at indeces next surrounding k 
k_index = 900 # will make loop

k_x = stream_proj_x_sort[sim_trail_index][k_index]
k_y = stream_poly_y_trail[k_index]

k_x_neigh = stream_proj_x_sort[sim_trail_index][k_index - k_neigh]
k_y_neigh = stream_poly_y_trail[k_index - k_neigh]

k_x_neigh_neg = stream_proj_x_sort[sim_trail_index][k_index + k_neigh]
k_y_neigh_neg = stream_poly_y_trail[k_index + k_neigh]

#Part of poly fit line to use:
k_poly_x = stream_proj_x_sort[sim_trail_index][k_index-k_neigh:k_index+k_neigh]# part of polynomial that surrounds node k 
k_poly_y = stream_poly_y_trail[k_index-k_neigh:k_index+k_neigh]

fig,axes = plt.subplots(1,1,figsize=(5,5))#, sharey=True, sharex=True)
matplotlib.rcParams.update({'font.size': 16})


axes.set_aspect('equal')
axes.set_xlim(20,0)
axes.set_ylim(-10,10)

axes.set_title('Area perserving projection', fontsize=18)
axes.scatter(stream_proj_x_sort[sim_trail_index],stream_proj_y_sort[sim_trail_index], s=1, color='steelblue', label="sim. points, trail")

#axes.plot(poly_trail_x ,stream_poly_y_trail, color="black", label="polyfit, trail") #poly fit to trailing arm
axes.plot(stream_proj_x_sort[sim_trail_index],stream_poly_y_trail, color="black", label="polyfit, trail") #poly fit to trailing arm
axes.plot(k_poly_x,k_poly_y,color="black", label="polyfit, trail") #poly fit to trailing arm



axes.plot(k_x,k_y, marker="o", markersize=5, color='black', linestyle='none', label="node, k") #poly fit to trailing arm
axes.plot(k_x_neigh,k_y_neigh, marker="o", markersize=5, color='red', linestyle='none', label="k, neigh") #poly fit to trailing arm
axes.plot(k_x_neigh_neg,k_y_neigh_neg, marker="o", markersize=5, color='magenta', linestyle='none', label="k, neigh") #poly fit to trailing arm

axes.set_xlabel('$X$ [deg]')
axes.set_ylabel('$Y$ [deg]')

axes.legend()
fig.tight_layout()
#Find u - derivative of the part of poly surrounding k
deriv = np.polyder(polyfit_trail)   
print(deriv)
#I already have a point to evaluate in: k_x, and I have it's value k_y
slope_at_node = np.polyval(deriv, k_x)
print(slope_at_node)
ylow = (k_x_neigh_neg - k_x )  * slope_at_node + k_y
yhigh = (k_x_neigh- k_x) * slope_at_node + k_y
#polyfit_trail
polyfit_trail
deriv_2 = np.polyder(deriv) #second derivative 
second_der_at_node = np.polyval(deriv_2, k_x)
xmin = np.round(k_x) - 3
xmax = np.round(k_x) + 3
ymin = np.round(k_y) - 3
ymax = np.round(k_y) + 3
fig,axes = plt.subplots(1,1,figsize=(5,5))#, sharey=True, sharex=True)
matplotlib.rcParams.update({'font.size': 16})


axes.set_aspect('equal')
axes.set_xlim(xmax,xmin)
axes.set_ylim(ymin,ymax)

axes.set_title('Area perserving projection', fontsize=18)
axes.plot(stream_proj_x_sort[sim_trail_index],stream_poly_y_trail, color="black", label="polyfit, trail", linestyle='--') #poly fit to trailing arm


axes.plot(k_x,k_y, marker="o", markersize=5, color='black', linestyle='none', label="node, k") #poly fit to trailing arm
axes.plot(k_x_neigh,k_y_neigh, marker="o", markersize=5, color='red', linestyle='none', label="k, neigh") #poly fit to trailing arm
axes.plot(k_x_neigh_neg,k_y_neigh_neg, marker="o", markersize=5, color='red', linestyle='none', label="k, neigh") #poly fit to trailing arm

#axes.plot([k_x,k_x_neigh],[k_y,yhigh], label="tangent") #from node, k
axes.plot([k_x_neigh_neg,k_x_neigh],[ylow,yhigh], label="tangent")


axes.set_xlabel('$X$ [deg]')
axes.set_ylabel('$Y$ [deg]')

axes.legend(loc="lower left", fontsize=14)
fig.tight_layout()<jupyter_output><empty_output><jupyter_text>I now have the tangent and sign of the secon derivative: compute u, v unit vectors from node k.<jupyter_code>u = np.array([[0,k_x-k_x_neigh],[0,k_y-yhigh]])
u_norm = u/np.sqrt((np.sum(u**2)))
u_x =  k_x - u_norm[0] 
u_y =  k_y - u_norm[1] 
import math

def rotate(origin, point, angle):
    """
    Rotate a point counterclockwise by a given angle around a given origin.

    The angle should be given in radians.
    """
    ox, oy = origin
    px, py = point

    qx = ox + math.cos(angle) * (px - ox) - math.sin(angle) * (py - oy)
    qy = oy + math.sin(angle) * (px - ox) + math.cos(angle) * (py - oy)
    return qx, qy
#Second derivative at node determines which way the v-vector is pointing
sign = 1
if second_der_at_node < 0:
    sign = -1
v_x2, v_y2 = rotate(origin = (u_x[0],u_y[0]), point = (u_x[1],u_y[1]), angle = sign*np.pi/2)
v_x = np.array([u_x[0], v_x2])
v_y = np.array([u_y[0], v_y2])
fig,axes = plt.subplots(1,1,figsize=(5,5))#, sharey=True, sharex=True)
matplotlib.rcParams.update({'font.size': 16})

axes.set_aspect('equal')
axes.set_xlim(xmax,xmin)
axes.set_ylim(ymin,ymax)

axes.set_title('Area perserving projection', fontsize=18)

axes.plot(k_x,k_y, marker="o", markersize=5, color='black', linestyle='none', label="node, k") #poly fit to trailing arm
axes.plot(stream_proj_x_sort[sim_trail_index],stream_poly_y_trail, color="black", label="polyfit, trail", linestyle='--') #poly fit to trailing arm
axes.plot(u_x,u_y, label='u, unit-vec') #from node, k
axes.plot(v_x,v_y, label='v, unit-vec') #from node, k


axes.set_xlabel('$X$ [deg]')
axes.set_ylabel('$Y$ [deg]')

axes.legend(loc="lower left", fontsize=14)
fig.tight_layout()<jupyter_output><empty_output><jupyter_text>Function that: given "node", k on polynomial fit, returns u,v-vectiors<jupyter_code>def u_v_unitvectors(node_k, polyfit_x, polyfit_y, fit):
    """Given an index at node, k, for a polynimial fit, the function returns a tangential, u, and 
    norm, v, unitvectors at node, k.
    
    "node_k" is and index included in the polyfit_arr. 
    
    "fit" is an array of coecfficients of the polynomial fit.
    "polyfit_x" is the array associated with the polynomial fit. 
    "polyfit_y" is the polynomial fit y-values 


    The function returns nromalized normal, v, and tangent, u, vectors from point at node_k along polyfit.
    
    """

    #Along the  poly-fit, k_index is the index of the node of interest, 
    #and k_neigh is it's nearest point on poly
    k_neigh = 1 
    k_index = node_k 


    #Computes the derivative from the polynomial fit coefficients
    deriv = np.polyder(fit)   
    slope_at_node = np.polyval(deriv, polyfit_x[k_index])
    
    #Sec. deriv determines which way the v-vector is pointing
    #Sign: is fed into the rotation-function below
    deriv_2 = np.polyder(deriv) 
    second_der_at_node = np.polyval(deriv_2, polyfit_x[k_index])
    sign = 1
    if second_der_at_node < 0:
        sign = -1


    #The x and y-value of the the tangent line for the u-vector 
    x_neigh = polyfit_x[k_index] - polyfit_x[k_index - k_neigh]
    y_neigh = polyfit_y[k_index] - ((polyfit_x[k_index - k_neigh]\
                                     - polyfit_x[k_index]) * slope_at_node + polyfit_y[k_index] )

 
    #tangential norm vec. 
    u = np.array([[0, x_neigh],[0, y_neigh]])
    u_norm = u/np.sqrt((np.sum(u**2)))
    u_x =  polyfit_x[k_index] - u_norm[0] 
    u_y =  polyfit_y[k_index] - u_norm[1] 
    u = np.array([u_x,u_y]) #normalized
    
    
    #normal-vector to node, k normalized 
    #calls "rotate function"
    v_x2, v_y2 = rotate(origin = (u_x[0],u_y[0]), point = (u_x[1],u_y[1]), angle = sign*np.pi/2)
    v_x = np.array([u_x[0], v_x2]) #same origin as u-vev
    v_y = np.array([u_y[0], v_y2]) #same origin as u-vev
    v = np.array([v_x,v_y]) 
    
    return u, v<jupyter_output><empty_output><jupyter_text>Gaussians along stream Vary "size" (width), amoplitude (density) and maybe position as<jupyter_code>def covarience(u, h, s_k):
    """Given tangent normalized vector at a given node k,
    this function takes the "bandwidth", h, and the function, s_k, 
    and returns the covarience matrix"""
    Ck_prim = np.zeros([2,2])
    Ck_prim[0,0] = h**2
    Ck_prim[1,1] = s_k**2
    
 #   print(u)
    theta = np.arctan2(u[1,0] -u[1,1] ,u[0,0] - u[0,1]) #might need to flip this when not in "flipped coord. sys"
  #  print(theta)

    R = np.zeros([2,2]) #rotation from u,v to x,y (sky-tangent-plane, area conserved )
    R[0,0] = np.cos(theta)
    R[0,1] = - np.sin(theta)
    R[1,0] = np.sin(theta)
    R[1,1] = np.cos(theta)
    
    Ck = R @ Ck_prim @ R.T
   # print(Ck.shape)
    return Ck
from scipy.special import logsumexp
from scipy.stats import multivariate_normal
#u[:,:,0]
def ln_density(xy, nodes, s_k, a_k, h, u):
    # TODO: compute C at location of each node, store as "Cs" array or list
    
    ln_dens = np.zeros([len(nodes),len(xy)])
    for k in range(len(nodes)):
      #  print(h)
        C = covarience(u[:,:,k], h, s_k[k])
        ln_dens[k] = multivariate_normal.logpdf(xy,nodes[k],C) + np.log(a_k[k]) 
        
                       
    return logsumexp(ln_dens,axis=0)<jupyter_output><empty_output><jupyter_text>Call functions for Gaussian mixture model<jupyter_code>print(polyfit_lead)
print(polyfit_trail)
 ### Leading arm #####
polyfit_lead_x = np.linspace(-20,0,1000)
polyfit_lead_y = np.poly1d(polyfit_lead)(polyfit_lead_x)

#### Trailing Arm ########
#Redefine my polyfit_x, polyfit_y
polyfit_trail_x = np.linspace(0,80,1000)
polyfit_trail_y = np.poly1d(polyfit_trail)(polyfit_trail_x)

np.linspace(1,len(polyfit_lead_x)-1,10, dtype=int)
np.arange(1,len(polyfit_lead_x)-1,np.int(len(polyfit_lead_x)/10.))
#z = np.array([ 1, 101, 201, 301, 401, 501,  651, 901])
#z= np.geomspace(50,len(polyfit_lead_x)-1,100, dtype=int)
#z
from utilities import get_uniform_idx
polyfit_lead_x = np.linspace(-20,0,10000)
polyfit_lead_y = np.poly1d(polyfit_lead)(polyfit_lead_x)
polyfit_trail_x = np.linspace(0,80,10000)
polyfit_trail_y = np.poly1d(polyfit_trail)(polyfit_trail_x)
nodes_arr = get_uniform_idx(polyfit_lead_x, polyfit_lead_y, spacing=5)
nodes_arr_trail = get_uniform_idx(polyfit_trail_x, polyfit_trail_y, spacing=5)
## Here's the array of "indices"/node placements
# Space out various nodes and save their positions in "nodes"-array 
i = 0
# nodes_arr = np.arange(1,len(polyfit_lead_x)-1,np.int(len(polyfit_lead_x)/10.))
#nodes_arr = z
nodes = np.zeros([len(nodes_arr),2])
for k in nodes_arr:
    nodes[i,:] = np.array([polyfit_lead_x[k],polyfit_lead_y[k]])
    i = i + 1
    
    
####TRAILING####
i = 0
#nodes_arr_trail = np.arange(1,len(polyfit_trail_x)-1,np.int(len(polyfit_trail_x)/10.))
nodes_trail = np.zeros([len(nodes_arr_trail),2])
for k in nodes_arr_trail:
    nodes_trail[i,:] = np.array([polyfit_trail_x[k],polyfit_trail_y[k]])
    i = i + 1
#nodes
### Call u/v-function to make u-vector of length based on number of nodes/indices to feed ln_density
u = np.zeros([2,2,len(nodes_arr)])
v = np.zeros([2,2,len(nodes_arr)])
for k in range(len(nodes_arr)):
    u[:,:,k],v[:,:,k] = u_v_unitvectors(node_k = nodes_arr[k] , \
                        polyfit_x = polyfit_lead_x, polyfit_y = polyfit_lead_y, fit = polyfit_lead)
    
    
### Call u/v-function to make u-vector of length based on number of nodes/indices to feed ln_density
u_trail = np.zeros([2,2,len(nodes_arr_trail)])
v_trail = np.zeros([2,2,len(nodes_arr_trail)])
for k in range(len(nodes_arr_trail)):
    u_trail[:,:,k],v_trail[:,:,k] = u_v_unitvectors(node_k = nodes_arr_trail[k] , \
                        polyfit_x = polyfit_trail_x, polyfit_y = polyfit_trail_y, fit = polyfit_trail)    
x, y = np.meshgrid(np.linspace(-30, 80, 128),
                   np.linspace(-80, 30, 128))
xy = np.stack((x.ravel(), y.ravel())).T
# this below should be executable 
K = len(nodes)
ak = np.ones(K)
# ak = np.random.uniform(size=K)
# ak /= ak.sum()
sk = np.random.uniform(0.5, 2, size=K)
derp = ln_density(xy, nodes = nodes, s_k=sk, a_k=ak, h = 3, u = u)


#trailinf
K = len(nodes_trail)
ak = np.random.uniform(size=K)
ak /= ak.sum()
derp_trail = ln_density(xy, nodes = nodes_trail, s_k = np.ones(K), a_k = ak, h = 3, u = u_trail)
fig,axes = plt.subplots(1,2,figsize=(14,7))#, sharey=True, sharex=True)
matplotlib.rcParams.update({'font.size': 22})

axes[0].set_aspect('equal')
axes[1].set_aspect('equal')
axes[0].set_xlim(80,-30)
axes[1].set_xlim(80,-30)
axes[0].set_ylim(-80,30)
axes[1].set_ylim(-80,30)

axes[0].set_title('Leading', fontsize=18)
axes[0].pcolormesh(x, y, np.exp(derp.reshape(x.shape) - np.max(derp)),cmap ='OrRd')
#axes[0].plot(polyfit_lead_x[nodes_arr] ,polyfit_lead_y[nodes_arr] , linestyle='--',color='black')#,marker='o', label='node_k')

axes[0].set_xlabel('RA [deg]')
axes[0].set_ylabel('dec [deg]')

axes[1].set_title('Trailing arm', fontsize=18)
axes[1].pcolormesh(x, y, np.exp(derp_trail.reshape(x.shape) - np.max(derp_trail)),cmap ='Blues')
#axes[1].plot(polyfit_trail_x[nodes_arr_trail] ,polyfit_trail_y[nodes_arr_trail] , linestyle='--',color='black')#,marker='o', label='node_k')

axes[1].set_xlabel('$X$ [deg]')
axes[1].set_ylabel('$Y$ [deg]')

fig.tight_layout()

#Checking I populate the correct "gaussian"
fig,axes = plt.subplots(1,1,figsize=(7,7))#, sharey=True, sharex=True)
matplotlib.rcParams.update({'font.size': 16})

# xmin = np.round(polyfit_x[node_k]) - 3
# xmax = np.round(polyfit_x[node_k]) + 3
# ymin = np.round(polyfit_y[node_k]) - 3
# ymax = np.round(polyfit_y[node_k]) + 3

axes.set_xlim(80,-20)
axes.set_ylim(-60,60)

axes.set_aspect('equal')
# axes.set_xlim(xmax,xmin)
# axes.set_ylim(ymin,ymax)

#axes.scatter(test_Ck.T[0],test_Ck.T[1], s=0.002, zorder = -100, color='g')
axes.plot(polyfit_x ,polyfit_y , linestyle='--', label='polyfit', color='steelblue')
axes.plot(polyfit_lead_x ,polyfit_lead_y , linestyle='--', label='polyfit', color='orange')

axes.plot(u[0],u[1])#, label='u, unit-vec', c='red') #from node, k
axes.plot(v[0],v[1])#, label='v, unit-vec') #from node, k

axes.plot(u_trail[0],u_trail[1])#, label='u, unit-vec', c='red') #from node, k
axes.plot(v_trail[0],v_trail[1])#, label='v, unit-vec') #from node, k

axes.plot(polyfit_trail_x[nodes_arr_trail] ,polyfit_trail_y[nodes_arr_trail] , linestyle='none',color='black',marker='o', label='node_k', markersize=2)
axes.plot(polyfit_lead_x[nodes_arr] ,polyfit_lead_y[nodes_arr] , linestyle='none',color='black',marker='o', label='node_k', markersize=2)

axes.set_xlabel('$X$ [deg]')
axes.set_ylabel('$Y$ [deg]')

axes.legend(loc="lower left", fontsize=14)
fig.tight_layout()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1553408597481546, 'fraction_numerical': 0.023284845853234912, 'mean_word_length': 5.373227257004497, 'pattern_counts': {'":': 0, '<': 33, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '584432', 'license': 'MIT', 'n_tokens_mistral': 8076, 'n_tokens_neox': 7593, 'n_words': 1460, 'path': 'notebooks/archive/DensityModel.ipynb', 'repo_name': 'adrn/Pal5sBiggestFan', 'size': '18424'}