import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import argrelextrema
import matplotlib.patheffects as PathEffects

# Constants---------------------------------------------------------------------
G = 6.6741*10**-11                                                  #m3/kgs2
Msun = 2*10**30                                                     #kg
Mcharon = 1.6e21                                                    #kg
Mpluto = 1.3e22                                                     #kg
rpluto = 1.188e6                                                    #m
rcharon = 6.06e5                                                    #m
AU = 1496*10**8                                                     #m
Mstar = Msun                                                        #kg
r = 39*AU                                                           #m
rho = 1778.                                                         #kg/m3
v_hw = -32.                                                         #m/s
Omega0 = np.sqrt(G*Mstar/r**3)                                      #1/s
xco = 2*v_hw/(3.*Omega0)                                            #m
vkepler = np.sqrt(G*Mstar/r)                                        #m/s
period = 2 * np.pi * r / vkepler                                    #s
variablessaved = 16


# Functions-------------------------------------------------------------------------------

# Makingconstants gives a complete set of new constants dependent on the system and mass-ratio
def makingconstants(ratio, tau_s):

    # Subsequent values--------------------------------------------------
    Mtot = Mpluto + Mcharon
    M1 = ratio * Mtot / (ratio + 1)
    M2 = Mtot / (ratio + 1)
    rHill = r * ((M1 + M2) / (3.*Mstar))**(1/3.)                        #m
    a_b = 0.0025                                                        # rHill
    distance = a_b * rHill                                              #m
    rplanet1 = ((3 * M1) / (4 * np.pi * rho))**(1/3.)                   #m
    rplanet2 = ((3 * M2) / (4 * np.pi * rho))**(1/3.)                   #m
    r1 = M2 / (M1 + M2) * distance                                      #m
    r2 = M1 / (M1 + M2) * distance                                      #m
    Omega1 = np.sqrt(G * Mtot / (distance)**3)                          #1/s
    Omega2 = np.sqrt(G * Mtot / (distance)**3)                          #1/s
    vk1 = Omega1 * r1                                                   #m/s
    vk2 = Omega1 * r2                                                   #m/s
    t_s = tau_s/Omega0                                                  #s

    return M1, M2, distance, tau_s, rplanet1, rplanet2, r1, r2, rHill, vk1, vk2, Omega1, Omega2, t_s

# MakingR0 gives a starting array for R0 given a certain starting phi and x_coordinate
def makingR0(x_0, phistart):

    # Making R0 the correct dimensions, then filling in all the values    
    R0 = np.zeros((11))

    R0[0] = x_0                                                                               # x_0
    R0[1] = max([50*rplanet1, 2e2*np.sqrt(G*(M1 + M2)*t_s/abs(v_hw))])                        # y_0
    R0[2] = 0.                                                                                # z_0
    R0[3] = 2*v_hw*tau_s/(1+tau_s**2)                                                         # vx_0
    R0[4] = 1.*v_hw/(1+tau_s**2) - 3./2*Omega0*x_0                                            # vy_0
    R0[5] = 0.                                                                                # vz_0
    R0[6] = phistart                                                                          # phi_0
    R0[7] = 0.                                                                                # t_0
    R0[8] = 0                                                                                 # hit
    R0[9] = ((x_0 - r1 * np.cos(phistart))**2 + (R0[1] - r1 * np.sin(phistart))**2)**0.5      # Distance1
    R0[10] = ((x_0 + r2 * np.cos(phistart))**2 + (R0[1] + r2 * np.sin(phistart))**2)**0.5     # Distance2
    
    return R0

# Singlesimulation makes a complete simulation for multiple x_starting coordinates given one phi-value
# NOTE: this function is not used in the paper, but can be used when it is desirable to simulate multiple pebbles at once
def singlesimulationold(R0, iterations, fast = True, printinginfo=True):

    '''Dimensions of R0: (subscript denotes which pebble)
    ---------------------------------------
    [[ x_0,  y_0,  z_0,  vx_0,  vy_0,  vz_0,  phi,  t, hit],
     [ x_1,  y_1,              ...,                 t, hit],
     [ x_2,  y_2,              ...,                 t, hit],
                               ...,
     [ x_n,  y_n,              ...,                 t, hit]]
    -------------------------------------
    '''

    # Making one array to encapture all pebble trajectories
    R = np.zeros((iterations, len(R0), len(R0[0])))
    R[0] = R0
    
    # Give the starting time step
    h = 100000
    percentage = 0
    
    # Spirals says how many spirals we count as a "hit", spiraliterations gives the amount of iterations we measure
    spirals = 8
    minimaarray = np.zeros((len(R0), spirals))
    
    # For loop going through the simulation
    for i, Rn in enumerate(R):

        # Skip the starting values
        if i>0:
            
            if i % (len(R) / 100) ==0 and printinginfo == True:
                percentage += 1
                print("########################################################################")
                print("Sub-Simulation at {}%, there are {} pebbles left".format(percentage,len(np.where(R[i-1,:,8] == 0)[0])))
                print("Variable time step is {}s".format(h))
            
            # Keep score of the hits, if every pebble is hit or out of bounds you can stop the simulation
            hits = R[i-1,:,8]

            # Checking if every pebble has already been done simulating
            if min(hits) == 0.:

                #calculating the next iteration + the variable time step
                R[i], h = RKF45(R[i-1], h)

            # Else we give every iteration the last known value
            else:
                R[i:len(R)], h = R[i-1], h
                print("Sub-Simulation done prematurely at iteration {}".format(i))
                break
        
        # We check every 500 iterations if there are orbits found to minimize computation time
        if fast == True and i % 500 == 0 and i > 0:

            # Make two lists of distances between pebbles and both planetesimals
            distance1 = R[0:i,:,-2]
            distance2 = R[0:i,:,-1]
            
            # Make lists of all the minima of pebbles (so we look for eliptical orbits)
            minima1, pebbleminima1 = argrelextrema(distance1, np.less)
            minima2, pebbleminima2 = argrelextrema(distance2, np.less)
            minima1counts = np.histogram(pebbleminima1, bins=np.arange(0,len(R[0])+1))[0]
            minima2counts = np.histogram(pebbleminima2, bins=np.arange(0,len(R[0])+1))[0]

            nohitpebbles = np.where(hits == 0)[0]

            # Going through all the pebbles that are still simulating
            for ipebble in nohitpebbles:

                # These are the last minima found on planetesimal 1 and 2
                pebbleminima1last = minima1[np.where(pebbleminima1==ipebble)][-spirals:]
                pebbleminima2last = minima2[np.where(pebbleminima2==ipebble)][-spirals:]

                # Make sure we do not divide by zero
                if len(pebbleminima1last) > 0 and len(pebbleminima2last) > 0:

                    # Convert these iterations to actual distances
                    lastdistances2 = distance2[pebbleminima2last, ipebble]
                    lastdistances1 = distance1[pebbleminima1last, ipebble]
                    
                    # If there are minima found, and they are all within one orbital distance (such that it isn't spiralling around both, just one) they will be accreted
                    if len(pebbleminima1last) == spirals:

                        # We also look if they are further away than the other planetesimal
                        if max(lastdistances1) < r1 and max(lastdistances1) < min(lastdistances2):
                            R[i,ipebble,8] = 1

                    # Same for planetesimal 2:
                    if len(pebbleminima2last) == spirals:
                        
                        if max(lastdistances2) < r2 and max(lastdistances2) < min(lastdistances1):
                            R[i,ipebble,8] = 2
                
    return R

# Acceleration calculates a velocity and acceleration vector for a given position and velocity
def acceleration(Rpebble, kvalue = np.zeros([6])):

    # Making useable arrays---------------------------------------------------------------
    Rpebble2 = Rpebble[0:6] + kvalue
    x = Rpebble2[0]
    y = Rpebble2[1]
    z = Rpebble2[2]
    vx = Rpebble2[3]
    vy = Rpebble2[4]
    vz = Rpebble2[5]
    phi = Rpebble[6]
    t = Rpebble[7]
    hit = Rpebble[8]

    # Calculate the locations and distances of the planetesimals w/ respect to pebble---------------------------------
    planetesimal1x = r1 * np.cos(phi)
    planetesimal1y = r1 * np.sin(phi)
    planetesimal1z = 0.
    planetesimal2x = r2 * -1 * np.cos(phi)
    planetesimal2y = r2 * -1 * np.sin(phi)
    planetesimal2z = 0.
    distance1 = ((x - planetesimal1x)**2 + (y - planetesimal1y)**2 + (z - planetesimal1z)**2)**0.5
    distance2 = ((x - planetesimal2x)**2 + (y - planetesimal2y)**2 + (z - planetesimal2z)**2)**0.5

    # Calculate the difference forces----------------------------------------------------------------------------------
    a_cor = np.array([2 * Omega0 * vy, -2*Omega0*vx, 0])
    a_tidal = np.array([3*x*Omega0**2, 0, 0])
    a_drag = -1./t_s * np.array([vx, vy - v_hw + 3./2*Omega0*x, vz])
    a_grav = -1.*G*M1/distance1**3 * np.array([x - planetesimal1x, y - planetesimal1y, z - planetesimal1z]) \
    -1.*G*M2/distance2**3 * np.array([x - planetesimal2x, y - planetesimal2y, z - planetesimal2z]) 

    newacc = a_cor + a_tidal + a_drag + a_grav

    # If statements to look for hits:
    if distance1 < rplanet1:
        hit = 1
    elif distance2 < rplanet2:
        hit = 2
    
    if True:
        # If statements to see if it is out of bounds
        if x < -1 * rHill:
            hit = 3
        elif y > rHill and vy > 0.:
            hit = 3
        elif xco < -1*rHill and y < -1 * rHill:
            hit = 3    

    # Update all the new values ---------------------------
    Rpebblenew = np.array([vx, vy, vz, newacc[0], newacc[1], newacc[2]])

    return Rpebblenew, hit, distance1, distance2

# RKF45 calculates every next iteration and defines next time-step
def RKF45(Rpebble, h):

    # Epsilon sets error-margin (needs to be at least 1e-6)
    epsilon = 1e-8
    accepted = False

    # Only return a new iteration when error is accepted, otherwhise calculate iteration again with new stepsize
    while not accepted:
        
        # First, all the k-values (with the first iteration we also calculate possible hits)-------------------------------
        k1, hit, distance1, distance2 = acceleration(Rpebble)
        k1 = h * k1
        k2 = h * acceleration(Rpebble,       1./4 * k1)[0]
        k3 = h * acceleration(Rpebble,      3./32 * k1 +      9./32 * k2)[0]
        k4 = h * acceleration(Rpebble, 1932./2197 * k1 - 7200./2197 * k2 + 7296./2197 * k3)[0]
        k5 = h * acceleration(Rpebble,   439./216 * k1 -         8. * k2 +  3680./513 * k3 -  845./4104 * k4)[0]
        k6 = h * acceleration(Rpebble,    -8./27  * k1 +         2. * k2 - 3544./2565 * k3 + 1859./4104 * k4 - 11./40 * k5)[0]

        # Then, calculate both the 5th, and 6th order values---------------------------------------------------------------
        vector = Rpebble[0:6]
        vector1 = (25./216 * k1 + 1408./2565  * k3 + 2197./4104   * k4 - 1./5  * k5)                 #5th order
        vector1_ = (16./135 * k1 + 6656./12825 * k3 + 28561./56430 * k4 - 9./50 * k5 + 2./55 * k6)   #6th order
        
        # Updating phi and t-----------------------------------------------------------------------------------------------
        phi = Rpebble[6]
        t = Rpebble[7]
        phinew = (phi + h*Omega1) % (2 * np.pi)
        tnew = t + h

        # Calculating delta and new variable time-step according to Es-Hagh (2005)-----------------------------------------
        rn1 = np.linalg.norm((vector1 - vector1_)/h)
        if rn1 > 0.:
            delta = .84 * (epsilon / rn1) ** 0.25
        else:
            delta = 5.

        if delta <= 0.1:
            h = 0.1 * h
        elif delta > 4.:
            h = 4. * h
        elif delta > 1.:
            h = delta * h
        
        # Accept if difference is less then accepted error. If still not satisfied, change h regardless
        if rn1 < epsilon:
            accepted = True
        elif accepted == False and delta > 0.1 and delta <= 1.:
            h = delta * h

    # Now to return the values in the right dimensions (so adding phi, t and hits)-------------------------------------
    Rpebblenew = np.zeros((len(Rpebble)))
    Rpebblenew[0:6] = vector + vector1_
    Rpebblenew[6] = phinew
    Rpebblenew[7] = tnew
    Rpebblenew[8] = hit
    Rpebblenew[9] = distance1
    Rpebblenew[10] = distance2

    return Rpebblenew, h

# Singlesimulation makes a complete simulation for a single pebble given R0
def singlesimulation(R0, iterations, fast = True):
    '''Dimensions of R: (subscript denotes which pebble)
    ---------------------------------------
    [[ x_0,  y_0,  z_0,  vx_0,  vy_0,  vz_0,  phi_0,  t_0, hit, Distance_1, Distance_2],
     [ x_1,  y_1,                           ...,                Distance_1, Distance_2],
     [ x_2,  y_2,                           ...,                Distance_1, Distance_2],
                                            ...,
     [ x_n,  y_n,                           ...,                Distance_1, Distance_2]]
    -------------------------------------
    '''

    # Make starting array
    R = np.zeros((iterations, len(R0)))
    R[0] = R0

    # Starting value for h (time step)
    h = 1000000.
    spirals = 8

    # Going through the simulation
    for i, Rn in enumerate(R):
        
        # Skip first iteration
        if i > 0:

            # Print info
            if False and i % 100 == 0:
                print(f'#################### iteration {i} #######################')
                print(f'Distances: {int(R[i-1,-2] / rplanet1)} & {int(R[i-1,-1] / rplanet2)}')
                print(f'Time step h: {h}')
                print("")

            # Calculate every next iteration    
            R[i], h = RKF45(R[i-1], h)
            
            # Stop simulation when there is a hit/out of bounds
            if R[i,8] != 0:
                R[i+1:len(R)] = R[i]
                #print(f"simulation stopped at iteration {i}")
                break
            
            # Check every 500 iterations if enough spirals are found
            if fast == True and i % 100 == 0:
                
                # Make lists with minimal extrema
                distances1 = R[:i+1,9]
                distances2 = R[:i+1,10]
                minima1iterations = argrelextrema(distances1, np.less)[0]
                minima2iterations = argrelextrema(distances2, np.less)[0]

                # There needs to be at least a "spirals" amount of minima on both planetesimals
                if len(minima1iterations) >= spirals and len(minima2iterations) >= spirals:

                    # Only take the last values
                    lastiterations1 = minima1iterations[-spirals:]
                    lastiterations2 = minima2iterations[-spirals:]
                    minima1 = distances1[lastiterations1]
                    minima2 = distances2[lastiterations2]

                    # Check if all the spirals are on one of the two planetesimals 
                    # (also make sure the minima arent further away than the orbital distance itself)
                    if min(minima2) > max(minima1) and max(minima1) < distance:
                        R[i,8] = 1
                        #print("It happened with planetesimal 1!")
                    elif min(minima1) > max(minima2) and max(minima2) < distance:
                        R[i,8] = 2
                        #print("It happened with planetesimal 2!")

    return R

# Plotonetrajectory gives a plot of a single pebble trajectory
def plotonetrajectory(R, moment = 0):

    # making the right lists for plotting
    ratio = rHill
    xlist = R[:,0] / ratio
    ylist = R[:,1] / ratio
    zlist = R[:,2] / ratio
    vxlist = R[:,3]
    vylist = R[:,4]
    vzlist = R[:,5]
    philist = R[:,6]
    tlist = R[:,7]
    result = R[-1,8]

    # Visualizing planetesimal orbit and hillsphere
    hillcirclex, hillcircley, circle1x, circle1y, circle2x, circle2y = np.zeros(100), np.zeros(100), np.zeros(100), np.zeros(100), np.zeros(100), np.zeros(100)
    planet1x, planet1y, planet2x, planet2y = np.zeros(100), np.zeros(100), np.zeros(100), np.zeros(100)

    # take the last value for phi and define the position of both planetesimals
    phiend = R[-1,6]
    p1x, p1y = r1*np.cos(phiend), r1*np.sin(phiend)
    p2x, p2y = -1*r2*np.cos(phiend), -1*r2*np.sin(phiend)
    for theta in range(100):
        hillcirclex[theta] = rHill*np.cos(theta/99.*2*np.pi) / ratio
        hillcircley[theta] = rHill*np.sin(theta/99.*2*np.pi) / ratio
        circle1x[theta] = r1*np.cos(theta/99.*2*np.pi) / ratio
        circle1y[theta] = r1*np.sin(theta/99.*2*np.pi) / ratio
        circle2x[theta] = r2*np.cos(theta/99.*2*np.pi) / ratio
        circle2y[theta] = r2*np.sin(theta/99.*2*np.pi) / ratio
        planet1x[theta] = (p1x + 1*rplanet1*np.cos(theta/99.*2*np.pi)) / ratio
        planet1y[theta] = (p1y + rplanet1*np.sin(theta/99.*2*np.pi)) / ratio
        planet2x[theta] = (p2x + 1*rplanet2*np.cos(theta/99.*2*np.pi)) / ratio
        planet2y[theta] = (p2y + 1*rplanet2*np.sin(theta/99.*2*np.pi)) / ratio


    fig = plt.figure(figsize=(4,4))
    plt.plot(circle1x, circle1y, color="black", linewidth = .5, linestyle = '--', label='Path of planetesimals')
    plt.plot(circle2x, circle2y, color="black", linewidth = .5, linestyle = '--', label='Path of planetesimals')
    plt.plot(xlist, ylist, color = "red", linewidth = 0.8, label = "Path of Pebble")
    print(len(xlist))

    if moment > 0:
        arrowlocation = arrowiterations[moment]
        grootte = 2
        plt.quiver(xlist[arrowlocation], ylist[arrowlocation], (xlist[arrowlocation]-xlist[arrowlocation-1])*grootte, \
            (ylist[arrowlocation]-ylist[arrowlocation-1])*grootte, color = 'red', headwidth = 7, minshaft = 2.5)
    plt.plot(planet1x, planet1y, color="black", linestyle = '-')#, label="Planet 1")
    plt.plot(planet2x, planet2y, color="black", linestyle = '-')#, label="Planet 2")
    plt.plot(hillcirclex, hillcircley, color="grey", linestyle = '--', linewidth = 1, label='Hill sphere')
    plt.axvline(x=xco/ratio, color="red", linestyle='--', label='x_co')
    plt.axis('square')
    plt.xlim(-.015,.015)
    plt.ylim(-.015,.015)
    plt.xticks([])
    plt.yticks([])  
    plt.xlabel('X')
    plt.ylabel('Y')
    #plt.savefig("figures/trajectoryD.pdf")
    plt.show()

# Fullsimulation does the full simulation (all different x- and phi-values) while saving inbetween
def fullsimulation(xlist, philist, loadingname = "0", saveas = "newfile.txt"):

    # Either loading in the file, or making a new list
    if loadingname == "0":
        hitlist = np.zeros((len(philist), len(xlist)))
    else:
        hitlist = np.loadtxt(loadingname)
    
    newhitlist = np.copy(hitlist)
    totaliterations = len(xlist) * len(philist)

    # Going through the list
    for iphi, row in enumerate(hitlist):

        # Only start simulating once a row is not simulated
        if sum(row) == 0:

            # Going through the various x-values
            for ix in range(len(row)):
                
                counter = len(xlist) * iphi + ix
                R0 = makingR0(xlist[ix], philist[iphi])
                R = singlesimulation(R0, iterations, fast = True)
                newhitlist[iphi, ix] = R[-1, 8]
                print(f'simulation {counter}/{totaliterations}, result = {int(R[-1,8])}')
            
            # Saving after every row
            np.savetxt(saveas, newhitlist)

# Heatmapplot makes a heatmapplot
def heatmapplot(hitmap, xstart, xeind):

    Letters = ['A','B','C','D','E','F']

    #hitmap[np.where(hitmap == 0)] = 10
    hitmap[np.where(hitmap == 1)] = 0.5
    hitmap[np.where(hitmap == 2)] = 0.9
    hitmap[np.where(hitmap == 3)] = 0

    for n in moments:
        chosenphi = n[1]
        chosenx = n[0]


    aspectratio = 1.75

    plt.figure(figsize=(14,8))
    plt.imshow(hitmap, cmap='YlOrBr')
    plt.scatter(moments[:,0],moments[:,1], s=16, color = 'blue')
    plt.xlabel(r'$x_0[R_\mathrm{H}]$', fontsize = 16)
    plt.ylabel(r'$\varphi_0[\mathrm{rad}]$', fontsize = 16)
    plt.axhline(y=250, color="red", linestyle='--', linewidth = 2)
    plt.yticks([0,80,160,239,319,398,478], [0,1,2,3,4,5,6], fontsize = 15)
    xticks = np.linspace(xstart, xeind, 7)
    plt.xticks([100,300,500,700,900],np.round(xticks[1:6]/rHill,2), fontsize = 15)
    for i in range(len(moments)):
        plt.text(moments[i,0]+5, moments[i,1], Letters[i], fontsize=20,path_effects=[PathEffects.withStroke(linewidth=3,foreground="w")])

    plt.savefig('hitmap6.pdf',bbox_inches='tight')
    plt.show()

# Making a distribution corrected for the shearing field of co-rotating frame of reference
def cumdist(x1, x2, xamount):
    '''Input: inner release distance x1, outter release distance x2, amount of pebbles
	returns: shear corrected distribution over x1 and x2 with array of x-coordinates
	''' 

    x = np.linspace(0,1,xamount)
    term1 = -1*v_hw / (1 + tau_s**2)
    term2 = (x2 - x1)
    term3 = 0.75 * Omega0
    term4 = (x2**2 - x1**2)

    N = (term1 * term2) + (term3 * term4)

    part1 = 2 * N * (3 * Omega0)**-1
    part2 = (3 * Omega0 * x * N**-1) + (term1**2 * N**-2) + (3 * Omega0 * N**-2 * (term1 * x1 + 0.75 * Omega0 * x1**2))
    part3 = -1*term1 * N**-1

    dist = (part2**0.5 + part3) * part1

    return dist

# Function to make 6 specific trajectoryplots:
def sixtrajectories():

    # Making the figure
    fig = plt.figure(figsize=(15,2.5))
    gs = fig.add_gridspec(1,6, hspace=0)
    axs = gs.subplots(sharex=True, sharey=True)
    Letters = ['A','B','C','D','E','F']

    # going through the 5 simulations:
    for i in range(6):

        # Doing the simulation
        [chosenx,chosenphi] = moments[i]
        arrowlocation = arrowiterations[i]
        R0 = makingR0(xlist[chosenx], philist[chosenphi])
        R = singlesimulation(R0, iterations)

        # Loading in info from simulation
        ratio = rHill
        X = R[:,0] / ratio
        Y = R[:,1] / ratio
        phiend = R[-1,6]

        # Making surroundings:
        if True:
            # Visualizing planetesimal orbit and hillsphere
            hillcirclex, hillcircley, circle1x, circle1y, circle2x, circle2y = np.zeros(100), np.zeros(100), np.zeros(100), np.zeros(100), np.zeros(100), np.zeros(100)
            planet1x, planet1y, planet2x, planet2y = np.zeros(100), np.zeros(100), np.zeros(100), np.zeros(100)

            # take the last value for phi and define the position of both planetesimals
            
            p1x, p1y = r1*np.cos(phiend), r1*np.sin(phiend)
            p2x, p2y = -1*r2*np.cos(phiend), -1*r2*np.sin(phiend)
            for theta in range(100):
                hillcirclex[theta] = rHill*np.cos(theta/99.*2*np.pi) / ratio
                hillcircley[theta] = rHill*np.sin(theta/99.*2*np.pi) / ratio
                circle1x[theta] = r1*np.cos(theta/99.*2*np.pi) / ratio
                circle1y[theta] = r1*np.sin(theta/99.*2*np.pi) / ratio
                circle2x[theta] = r2*np.cos(theta/99.*2*np.pi) / ratio
                circle2y[theta] = r2*np.sin(theta/99.*2*np.pi) / ratio
                planet1x[theta] = (p1x + 1*rplanet1*np.cos(theta/99.*2*np.pi)) / ratio
                planet1y[theta] = (p1y + rplanet1*np.sin(theta/99.*2*np.pi)) / ratio
                planet2x[theta] = (p2x + 1*rplanet2*np.cos(theta/99.*2*np.pi)) / ratio
                planet2y[theta] = (p2y + 1*rplanet2*np.sin(theta/99.*2*np.pi)) / ratio

        # Plotting:
        axs[i].plot(circle1x, circle1y, color="black", linewidth = .5, linestyle = '--', label='Path of planetesimals')
        axs[i].plot(circle2x, circle2y, color="black", linewidth = .5, linestyle = '--', label='Path of planetesimals')
        axs[i].plot(planet1x, planet1y, color="black", linestyle = '-', label="Planet 1")
        axs[i].plot(planet2x, planet2y, color="black", linestyle = '-', label="Planet 2")
        axs[i].plot(X, Y, color = "red", linewidth = 0.8, label = "Path of Pebble")
        grootte = 2
        axs[i].quiver(X[arrowlocation], Y[arrowlocation], (X[arrowlocation]-X[arrowlocation-1])*grootte, \
            (Y[arrowlocation]-Y[arrowlocation-1])*grootte, color = 'red', headwidth = 7, minshaft = 1.0)
        axs[i].set_xticks([])
        axs[i].set_yticks([])
        axs[i].text(-0.012, 0.010, Letters[i], fontsize=17)

        axs[i].set_xlim(-.013,.013)
        axs[i].set_ylim(-.013,.013)
        axs[i].set_aspect('equal')

    fig.tight_layout()
    plt.savefig('6trajectories(6x1).pdf')
    plt.show()
        
    
#----------------------- Running the code -----------------------#

Nx = 1000
Nphi = 500
xstart, xend = 1.25e9, 3.25e9
tau_s = 0.01
xlist = np.linspace(xstart, xend, Nx)
philist = np.linspace(0,2 *np.pi, Nphi+1)[:-1]
iterations = 5000

# Loading in our fiducial hitmap (Which you would get if you run fullsimulation() with the lists above):
ratio = 1
R1 = np.loadtxt('fiducialhitmap.txt')
M1, M2, distance, tau_s, rplanet1, rplanet2, r1, r2, rHill, vk1, vk2, Omega1, Omega2, t_s = makingconstants(ratio, tau_s)


# The "moments" and "arrowiteratoins"-arrays are just to help with the making of figure 2
moments = np.array([[259,323],[385,254],[530,360],[580,148],[790,224],[943,375]])
arrowiterations = [720,695,650,650,710,860]
heatmapplot(R1, xstart, xend)
sixtrajectories()


# For a simulation of a single pebble:
chosenx = 200
chosenphi = 150
R0 = makingR0(xlist[chosenx], philist[chosenphi])
R = singlesimulation(R0, iterations)
plotonetrajectory(R)