I was trying to do something for a project where I create two random 3D vectors, project them onto a 2D plane by removing their x-component, taking the dot product between the two to find their separation and then forcing their separation to be between 0º and 90º.

I'm then running into the issue where I can't plot my data how I want to. I want to create 3 CDF's onto the same plot where we specify angles between 0-90º. I want to show the CDF for the following regions 0-20, 70-90 and 0-90. The issue I'm running into is either: my CDF will either be a single linear line, a CDF that only works between the specified angles and completely disappears outside of that range, and or a CDF that works only within its specified region and then flatlines outside of my specified range.

I will include an example of what I want my plot to look like.

import numpy as np
import matplotlib.pyplot as plt

num_samples = 50000 #number of samples
theta_all = np.arccos(np.random.uniform(-1, 1, num_samples)) # Uniform in cos(theta) 0-180º
phi_all = np.random.uniform(0, 2 * np.pi, num_samples) # uniform vectors between 0-2pi or 0-360º

#create magnetic field vector components
Bx = np.sin(theta_all) * np.cos(phi_all) 
By = np.sin(theta_all) * np.sin(phi_all) 
Bz = np.cos(theta_all) 
B = np.vstack((Bx, By, Bz)).T # Magnetic field vector as a matrix

#create disk vectors components
d_theta = np.arccos(np.random.uniform(-1, 1, num_samples)) 
d_phi = np.random.uniform(0, 2 * np.pi, num_samples) 
dx = np.sin(d_theta) * np.cos(d_phi) 
dy = np.sin(d_theta) * np.sin(d_phi) 
dz = np.cos(d_theta)
dv = np.vstack((dx, dy, dz)).T

# --- Project onto y-z plane (ignore x component) ---
dv_proj = dv[:, 1:]  # (y, z)
B_proj = B[:, 1:]    # (y, z)

# --- Dot product in y-z plane ---
dot_product = np.sum(dv_proj * B_proj, axis=1)
mag_dv = np.linalg.norm(dv_proj, axis=1)
mag_B = np.linalg.norm(B_proj, axis=1)

# Avoid divide by zero
valid = (mag_dv > 1e-9) & (mag_B > 1e-9)
cos_gamma = np.ones(num_samples)
cos_gamma[valid] = dot_product[valid] / (mag_dv[valid] * mag_B[valid])
cos_gamma = np.clip(cos_gamma, -1.0, 1.0)

# --- Projected angles in degrees ---
gamma_deg = np.degrees(np.arccos(cos_gamma))
projected_angles = np.minimum(gamma_deg, 180 - gamma_deg)