# Orbyts - UCL | Stellar International School

## Histogram

plt.hist(theta_Bnhat_list, bins=5, 
         facecolor='#313132', 
         edgecolor='#e0e0e0', 
         linewidth=0.5, 
         alpha=0.7)

## Coplanarity Normal

def coplanarity_normal(Bu, Bd):
    """ Coplanarity normal to the shock surface
        derived from coplanarity analysis formula.

    Args:
        Bu (list, array): Upstream magnetic field vector.
        Bd (list, array): Downstream magnetic field vector.

    Returns:
        n_cp (array): Coplanrity normal vector. 

    Example:
    >>> Bu = [-0.2, -0.3, 0.2]
    >>> Bd = [-0.5, -1, 0.4]
    >>> n_cp = coplanarity_normal(Bu, Bd)
    >>> print(n_cp, np.linalg.norm(n_cp))
    """
    
    Bu = np.array(Bu)
    Bd = np.array(Bd)
    dB = Bd-Bu
    BdXBu = np.cross(Bd, Bu)
    dBXBdXBu = np.cross(dB, BdXBu)
    n_cp = dBXBdXBu / np.linalg.norm(dBXBdXBu)
    return n_cp

## Calculate coplanarity normal

n_cp_samples = []
for i in range(N_CHUNKS):
    for j in range(N_CHUNKS):
        Bu = B_up_mean[i][0:3] # starrting from start of B_u
        Bd = B_down_mean[j][0:3] # starting from the end of B_d
        n_cp = coplanarity_normal(Bu, Bd)
        n_cp_samples.append(n_cp)
        #print(n_cp, np.linalg.norm(n_cp))
n_cp_samples = np.array(n_cp_samples)
print(n_cp_samples.shape)

## Calculate shock angle relative to the coplanarity shock normal

B_up_mean_averaged = B_up_mean[:,0:3].mean(axis=0)
theta_Bncp_samples = []

# loop through each coplanarity normal and do dot product with B_u
for i in range(len(n_cp_samples)):
    theta_Bncp = np.arccos(np.dot(B_up_mean_averaged, n_cp_samples[i])/np.linalg.norm(B_up_mean_averaged))*180/np.pi
    if theta_Bncp > 90:
        theta_Bncp = 180 - theta_Bncp
    theta_Bncp_samples.append(theta_Bncp)

theta_Bncp_mean = np.mean(theta_Bncp_samples)
theta_Bncp_std = np.std(theta_Bncp_samples)
theta_Bncp_mean, theta_Bncp_std

## Calculate shock angle relative to coplanarity normal in a loop:
The first part is the same as the code you have before when calculating the model normal.

theta_Bncp_samples_list = []
theta_Bncp_mean_list = []
theta_Bncp_std_list = []
ncp_samples_list = []

for i in range(len(df_BS_crossings)):
    print(f'BS crossing {i}:')

    # Get crossing time
    crossing_time = df_BS_crossings.iloc[i].TIME

    # Define interval
    time_start = crossing_time - pd.Timedelta(7, unit='m')
    time_end = crossing_time + pd.Timedelta(7, unit='m')
    df_crossing = df_1s.query("TIME > @time_start & TIME < @time_end")
    
    # Define B_up and B_down
    if df_BS_crossings.iloc[i].DIRECTION == 'I':
        B_up = df_crossing.query("TIME < @crossing_time")
        B_down = df_crossing.query("TIME > @crossing_time")
    elif df_BS_crossings.iloc[i].DIRECTION == 'O':
        B_up = df_crossing.query("TIME > @crossing_time")
        B_down = df_crossing.query("TIME < @crossing_time")

    # Split into 8 chunks
    N_CHUNKS = 8

    # Split B_up and B_down into 8 subintervals 
    B_up_chunks = np.array_split(B_up, N_CHUNKS)
    B_down_chunks = np.array_split(B_down, N_CHUNKS)
    
    # Calculate mean of each chunk giving 8 mean values for B_up and B_down
    B_up_mean = np.zeros((N_CHUNKS, 4))
    for j in range(N_CHUNKS):
        B_up_mean[j] = B_up_chunks[j].mean().values

    B_down_mean = np.zeros((N_CHUNKS, 4))
    for j in range(N_CHUNKS):
        B_down_mean[j] = B_down_chunks[j].mean().values

This next part is different, because we need to calculate coplanarity normal:

    # Calculate the coplanarity normal (64 samples)
    n_cp_samples = []
    # ENTER YOUR CODE HERE
    
    
    
    n_cp_samples = np.array(n_cp_samples)
    ncp_samples_list.append(n_cp_samples)

    # Calculate shock angle relative to the coplanarity shock normal
    # ENTER YOUR CODE HERE
    
    
    # Calculate mean
    theta_Bncp_mean = np.mean(theta_Bncp_samples)
    theta_Bncp_mean_list.append(theta_Bncp_mean)
    theta_Bncp_samples_list.append(theta_Bncp_samples)

    # Calculate errors: 1 standard deviation
    theta_Bncp_std = np.std(theta_Bncp_samples)
    theta_Bncp_std_list.append(theta_Bncp_std)

### The solution to Analysis 2

theta_Bncp_samples_list = []
theta_Bncp_mean_list = []
theta_Bncp_std_list = []
ncp_samples_list = []

for i in range(len(df_BS_crossings)):
    print(f'BS crossing {i}:')

    # Get crossing time
    crossing_time = df_BS_crossings.iloc[i].TIME

    # Define interval
    time_start = crossing_time - pd.Timedelta(7, unit='m')
    time_end = crossing_time + pd.Timedelta(7, unit='m')
    df_crossing = df_1s.query("TIME > @time_start & TIME < @time_end")
    
    # Define B_up and B_down
    if df_BS_crossings.iloc[i].DIRECTION == 'I':
        B_up = df_crossing.query("TIME < @crossing_time")
        B_down = df_crossing.query("TIME > @crossing_time")
    elif df_BS_crossings.iloc[i].DIRECTION == 'O':
        B_up = df_crossing.query("TIME > @crossing_time")
        B_down = df_crossing.query("TIME < @crossing_time")

    # Split into 8 chunks
    N_CHUNKS = 8

    # Split B_up and B_down into 8 subintervals 
    B_up_chunks = np.array_split(B_up, N_CHUNKS)
    B_down_chunks = np.array_split(B_down, N_CHUNKS)
    
    # Calculate mean of each chunk giving 8 mean values for B_up and B_down
    B_up_mean = np.zeros((N_CHUNKS, 4))
    for j in range(N_CHUNKS):
        B_up_mean[j] = B_up_chunks[j].mean().values

    B_down_mean = np.zeros((N_CHUNKS, 4))
    for j in range(N_CHUNKS):
        B_down_mean[j] = B_down_chunks[j].mean().values

    # Calculate the coplanarity shock normal (64 samples)
    n_cp_samples = []
    for j in range(N_CHUNKS):
        for k in range(N_CHUNKS):
            Bu = B_up_mean[j][0:3] # starrting from start of B_u
            Bd = B_down_mean[k][0:3] # starting from the end of B_d
            n_cp = coplanarity_normal(Bu, Bd)
            n_cp_samples.append(n_cp)
            #print(n_cp, np.linalg.norm(n_cp))
    n_cp_samples = np.array(n_cp_samples)
    ncp_samples_list.append(n_cp_samples)

    # Calculate shock angle relative to the coplanarity shock normal
    B_up_mean_averaged = B_up_mean[:,0:3].mean(axis=0)
    theta_Bncp_samples = []
    for j in range(len(n_cp_samples)):
        theta_Bncp = np.arccos(np.dot(B_up_mean_averaged, n_cp_samples[j])/np.linalg.norm(B_up_mean_averaged))*180/np.pi
        if theta_Bncp > 90:
            theta_Bncp = 180 - theta_Bncp
        theta_Bncp_samples.append(theta_Bncp)
    theta_Bncp_mean = np.mean(theta_Bncp_samples)
    theta_Bncp_mean_list.append(theta_Bncp_mean)
    theta_Bncp_samples_list.append(theta_Bncp_samples)

    # errors correspond to 1 standard deviation
    theta_Bncp_std = np.std(theta_Bncp_samples)
    theta_Bncp_std_list.append(theta_Bncp_std)

### Plot the results

plt.hist(theta_Bncp_mean_list, bins=5, facecolor='#002045', edgecolor='#e0e0e0', linewidth=0.5, alpha=0.7)
plt.xlim([0,90])
plt.xlabel(r'Shock angle $\theta_{Bncp}$(degree)')
plt.ylabel('Frequency')
#plt.savefig("./output/fig_shock_angle_coplanarity_distribution.jpg", dpi=300)

## :busts_in_silhouette: Group 1 (Yvonne, Dayana, Sachiin, Diviyesh)
You will be responsible for plotting skewness vs kurtosis of coplanarity shock angles.

The skewness tells whether the distribution of the shock angle for a particular crossing is 'normal' or not.

A skewness of 0 means normally distributed (bell curve).
Skewness > 0 means more weight in the left tail.
Skewness < 0 means more weight in the right tail.

A kurtosis of 3 means normally distributed (bell curve).
Kurtosis > 3 means heavier tail (less peaked).
Kurtosis < 3 means lighter tail (more peaked).

[Click here](https://scontent-hkg4-2.xx.fbcdn.net/v/t1.6435-9/96665492_125491382456362_2119748312570527744_n.jpg?_nc_cat=110&ccb=1-7&_nc_sid=8bfeb9&_nc_ohc=BHo_-DU9KowAX_5_YEq&_nc_ht=scontent-hkg4-2.xx&oh=00_AT-GOh1aHnR_RtiROkuztdiIrBpYW6UlBgugtyNBZeJ4Ig&oe=62E25E37) to see example distributions for different skewness and kurtosis.

# Calculate skewness and kurtosis 
from scipy.stats import skew
from scipy.stats import kurtosis

theta_Bncp_skewness = [skew(ENTER CODE HERE) for i in range(len(theta_Bncp_samples_list))]

theta_Bncp_kurtosis = [kurtosis(ENTER CODE HERE, fisher=False) for i in range(len(theta_Bncp_samples_list))]

### Plot the results
Make the plot and save it in your colab session.

plt.scatter(theta_Bncp_kurtosis, theta_Bncp_skewness, marker='+', color='k')
plt.axvline(x=3, color='k', linestyle='--', alpha=0.5)
plt.axhline(y=0, color='k', linestyle='--', alpha=0.5)
plt.xlabel(r'Kurtosis of $\theta_{BN}$')
plt.ylabel(r'Skewness of $\theta_{BN}$')
#plt.savefig("fig_shock_angle_skewness_kurtosis.jpg", dpi=300)

## :busts_in_silhouette: Group 2 (Belle, Junhao, Rui, Tingyu)
You will be responsible for plotting the angle between coplanarity normal and model normal ($\theta_{CPM}$) vs coplanarity shock angle:

You will need to store the $\hat{n}$ vectors from Analysis 1 where we calculated the shock angle using the BS model.

nhat_list = []
...

nhat_list.append(nhat)

Now to calculate the angle between the model normal and the coplanarity normal we do the following:

# loop through the nhat and ncp, and compute the dot product
theta_cpm_samples_list = []
theta_cpm_mean_list = []
theta_cpm_std_list = []
for i in range(len(ncp_samples_list)):
    nhat = nhat_list[i]
    ncp_samples = ncp_samples_list[i]
    theta_cpm_samples = []
    for j in range(len(ncp_samples)):
        theta_cpm = # ENTER CODE HERE
        if theta_cpm > 90:
            theta_cpm = 180 - theta_cpm
        theta_cpm_samples.append(theta_cpm)
    
    theta_cpm_samples_list.append(theta_cpm_samples)
    theta_cpm_mean_list.append(np.mean(theta_cpm_samples))
    theta_cpm_std_list.append(np.std(theta_cpm_samples))

### Plot results
Make the plot and save it in your colab session.

plt.errorbar(x=theta_Bncp_mean_list, y=theta_cpm_mean_list, 
            xerr=theta_Bncp_std_list, yerr=theta_cpm_std_list, 
            fmt='+', color='k', linewidth=1, alpha=0.5)
plt.axvline(x=45, color='k', linestyle='--', alpha=0.5)
plt.xlabel(r'Shock Angle $\theta_{BN}$ (degrees)')
plt.ylabel(r'Angle between model and coplanarity normals $\theta_{CPM}$ (degrees)')
#plt.savefig("fig_thetacpm_vs_shock_angle.jpg", dpi=300)

## :busts_in_silhouette: Group 3 (Angelina, Boston, Boon, Ian)
You will be responsible for plotting Shock angle deflection $\alpha$ vs field strength ratio $Bd/Bu$

Deflection angle is the angle between $B_d$ and $B_u$.

# shock deflection angle
# Do dot product of Bu vs Bd to get 64 deflection angles, calculate 1 standard deviation
# Calculate Bd/Bu to get 64 ratios, calculate 1 standard deviation
deflection_angle_samples_list = []
deflection_angle_mean_list = []
deflection_angle_std_list = []
BdToBu_samples_list = []
BdToBu_mean_list = []
BdToBu_std_list = []

for i in range(len(df_BS_crossings)):
    print(f'BS crossing {i}:')

    # Get crossing time
    crossing_time = df_BS_crossings.iloc[i].TIME

    # Define interval
    time_start = crossing_time - pd.Timedelta(7, unit='m')
    time_end = crossing_time + pd.Timedelta(7, unit='m')
    df_crossing = df_1s.query("TIME > @time_start & TIME < @time_end")
    
    # Define B_up and B_down
    if df_BS_crossings.iloc[i].DIRECTION == 'I':
        B_up = df_crossing.query("TIME < @crossing_time")
        B_down = df_crossing.query("TIME > @crossing_time")
    elif df_BS_crossings.iloc[i].DIRECTION == 'O':
        B_up = df_crossing.query("TIME > @crossing_time")
        B_down = df_crossing.query("TIME < @crossing_time")

    # Split B_up and B_down into 8 subintervals
    N_CHUNKS = 8
    B_up_chunks = np.array_split(B_up, N_CHUNKS)
    B_down_chunks = np.array_split(B_down, N_CHUNKS)
    
    # Calculate mean of each chunk giving 8 mean values for B_up and B_down
    B_up_mean = np.zeros((N_CHUNKS, 4))
    for j in range(N_CHUNKS):
        B_up_mean[j] = B_up_chunks[j].mean().values

    B_down_mean = np.zeros((N_CHUNKS, 4))
    for j in range(N_CHUNKS):
        B_down_mean[j] = B_down_chunks[j].mean().values
        
    # Deflection angle analysis STARTS HERE
    BdToBu_samples = []
    deflection_angle_samples = []
    for j in range(N_CHUNKS):
        for k in range(N_CHUNKS):
            Bu = B_up_mean[j][0:3] # starting from start of B_u
            Bd = B_down_mean[k][0:3] # starting from the end of B_d
            # Calculate deflection angle $\alpha$
            deflection_angle = #ENTER CODE HERE
            if deflection_angle > 90:
                deflection_angle = 180 - deflection_angle
            deflection_angle_samples.append(deflection_angle)

            # Calculate Bd/Bu ratio
            BdToBu_samples.append(np.linalg.norm(Bd)/np.linalg.norm(Bu))


    # store results in a list
    deflection_angle_samples_list.append(deflection_angle_samples)
    deflection_angle_mean_list.append(np.nanmean(deflection_angle_samples))
    deflection_angle_std_list.append(np.nanstd(deflection_angle_samples))

    BdToBu_samples_list.append(BdToBu_samples)
    BdToBu_mean_list.append(np.nanmean(BdToBu_samples))
    BdToBu_std_list.append(np.nanstd(BdToBu_samples))

### Plot results
Make the plot and save it in your colab session.

plt.errorbar(x=BdToBu_mean_list, y=deflection_angle_mean_list,
            xerr=BdToBu_std_list, yerr=deflection_angle_std_list, 
            fmt='+', color='k', linewidth=1, alpha=0.5)
plt.xlabel(r'Downstream / Upstream Field Ratio $B_d/B_u$')
plt.ylabel(r'Shock deflection angle $\alpha$ (degrees)')
#plt.savefig("fig_shock_deflection_angle_vs_BdToBu.jpg", dpi=300)

# Presentation
Intro and method sections:
* Slides 4-9: Angelina, Dayana, Belle, Boston
* Slides 10-13: Divyesh, Junhao, Boon, Rui, Ian, Tingyu

Results sections:
* Results 1: Rui
* Results 2: Boon
* Results 3: Group 1
* Results 4: Group 2
* Results 5: Group 3

Conclusion:
* Slide :