Hitting a specified target region on a sphere

I am attempting to create a program that will count the number of hits to a specific rectangular area on the surface of a sphere. How the programs is supposed to work, is random lines are generated and if one of those line hits in the target area the count goes up one. My problem is I do not think I am generating the lines correctly and I really have know idea how to correctly set the count parameters. This is the code I have so far and how I think the lines should be generated and what the count parameter might be.

import numpy as np
import random as rand
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


fig = plt.figure()

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

ax.set_aspect("equal")  

#rough model of the earth

theta, phi =  np.mgrid[0:2*np.pi : 20j ,0:np.pi : 20j]

r = 6.3

x = r * np.cos(phi)*np.sin(theta)
y = r * np.sin(phi)*np.sin(theta)
z = r * np.cos(theta)

ax.plot_wireframe(x,y,z, color = "k")

#target area
lat1x = 46.49913179 * (2*np.pi/360)
lat2x = 46.4423682 * (2*np.pi/360)
long1y = -119.4049072 * (2*np.pi/360)
long2y = -119.5048141 * (2*np.pi/360)

lat3x = 46.3973998 * (2*np.pi/360)
lat4x = 46.4532495 * (2*np.pi/360)
long3y = -119.4495392 * (2*np.pi/360)
long4y = -119.3492884 * (2*np.pi/360)

def to_cartesian(lat,lon):
    x = r * np.cos(lon)*np.cos(lat)
    y = r * np.sin(lon)*np.cos(lat)
    z = r * np.sin(lat)
    return [x,y,z]

p1 = to_cartesian(lat1x,long1y)
p2 = to_cartesian(lat2x,long2y)
p3 = to_cartesian(lat3x,long3y)
p4 = to_cartesian(lat4x,long4y)

X = np.array([p1,p2,p3,p4])
ax.scatter(X[:,0],X[:,1],X[:,2], color = "r")

#random line path

n = 500

x0 = np.zeros(n)
y0 = np.zeros(n)
z0 = np.zeros(n)

x1 = np.zeros(n)
y1 = np.zeros(n)
z1 = np.zeros(n)

for k in range (n):
     theta = rand.uniform(0.0, np.pi)
     phi = rand.uniform(0, (2 * np.pi)) 
     x0[k] = 100 * np.sin(phi) * np.cos(theta)
     y0[k] = 100 * np.sin(phi) * np.sin(theta)
     z0[k] = 100 * np.cos(theta)

for j in range (n):
    theta = rand.uniform(0.0, np.pi)
    phi = rand.uniform(0, (2 * np.pi))
    x1[j] = 100 * np.sin(phi) * np.cos(theta)
    y1[j] = 100 * np.sin(phi) * np.sin(theta)
    z1[j] = 100 * np.cos(theta)
#ax.plot_wireframe([x0[k],x1[j]],[y0[k],y1[j]],[z0[k],z1[j]], color="g")

# count if hit target area

count = 0

for i in range (n):
    if np.any([x1[i]<=X[:0]<=x0[i]]) and np.any([y1[i]<=X[:1]<=y0[i]]) and 
    np.any([z1[i]<=X[:2]<=z0[i]]):
         count =+1
 print (count)

 plt.show()