# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

import scipy.optimize
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

data = pd.read_csv('data.csv')
data['psi_up'] = data['180 + psi (up)'] - 180
data['psi_down'] = data['180 + psi (down)'] - 180
data['phi_1'] = data['2 phi_1'] / 2
data.head()


# $\varphi_2 = \psi + A - \varphi_1$
#
# $\sin\varphi_1 = n \sin\beta_1$
#
# $\sin\varphi_2 = n \sin\beta_2$
#
# $\beta_1+\beta_2=A$

# +
def get_error_function(phi_1, psi, A):
    A *= np.pi / 180
    phi_1 *= np.pi / 180
    psi *= np.pi / 180
    
    phi_2 = psi + A - phi_1
    
    def error_function(n):
        beta_1 = np.arcsin(np.sin(phi_1) / n)
        beta_2 = np.arcsin(np.sin(phi_2) / n)
        return beta_1 + beta_2 - A
    
    return error_function

def solve_n(phi_1, psi, A):
    if np.isnan(psi):
        return np.nan
    else:
        return scipy.optimize.newton_krylov(get_error_function(phi_1, psi, A), 1.55, f_tol=1e-12)


# +
A = 38
EPSILON = 1e-4

for direction in ['up', 'down']:
    data[f'n_{direction}']   = data.apply(lambda x: solve_n(x['phi_1'], x[f'psi_{direction}'], A), axis=1)
    data[f'd_phi_1_n_{direction}'] = (data.apply(lambda x: solve_n(x['phi_1'] + EPSILON, x[f'psi_{direction}'], A), axis=1)
                            - data[f'n_{direction}']) / EPSILON
    data[f'd_psi_{direction}_n_{direction}'] = (data.apply(lambda x: solve_n(x['phi_1'], x[f'psi_{direction}'] + EPSILON, A), axis=1)
                            - data[f'n_{direction}']) / EPSILON
# data['n_down'] = data.apply(lambda x: solve_n(x['phi_1'], x['psi_down'], A), axis=1)?
# -

data.head()

# +
d_angle = 0.5
d_phi_1 = d_angle / 2

for direction in ['up', 'down']:
    data[f'd_n_{direction}'] = ((data[f'd_phi_1_n_{direction}'] * d_phi_1) ** 2 + 
    (data[f'd_psi_{direction}_n_{direction}'] * d_angle) ** 2) ** 0.5

# +
data['cos2_theta_down'] = (np.sin(data.phi_1 / 180 * np.pi) / data.n_down) ** 2
data['d_cos2_theta_down'] = data['cos2_theta_down'] * (
    np.abs(np.cos(data.phi_1 / 180 * np.pi) * d_phi_1 / 180 * np.pi / np.sin(data.phi_1 / 180 * np.pi)) +
    np.abs(data.d_n_down / data.n_down)) * 2

data['cos2_theta_up'] = (np.sin(data.phi_1 / 180 * np.pi) / data.n_up) ** 2
data['d_cos2_theta_up'] = data['cos2_theta_up'] * (
    np.abs(np.cos(data.phi_1 / 180 * np.pi) * d_phi_1 / 180 * np.pi / np.sin(data.phi_1 / 180 * np.pi)) +
    np.abs(data.d_n_up / data.n_up)) * 2

data['d_cos2_theta_up'] = data['d_cos2_theta_up'].fillna(0)
data['d_cos2_theta_down'] = data['d_cos2_theta_down'].fillna(0)

data.head(1)
# -

data[['cos2_theta_down', 'd_cos2_theta_down', 'n_down', 'd_n_down']].to_csv('down.data', sep=' ', header=False, index=False)
data[['cos2_theta_up', 'd_cos2_theta_up', 'n_up', 'd_n_up']].to_csv('up.data', sep=' ', header=False, index=False)

np.sin((data.psi_up.min() + A + 1) / 2 / 180 * np.pi) / np.sin((A + 0.5) / 2 / 180 * np.pi)

data.psi_up.min()

phi = -4.75 / 180 * np.pi
A2 = A - 1
1/np.sin(A2 / 180 * np.pi) * (np.sin(phi) ** 2 + 1 + 2 * np.cos(A2 / 180 * np.pi) * np.sin(phi))

# plt.scatter((np.sin(data.phi_1 / 180 * np.pi) / data.n_down) ** 2, data.n_down)
plt.scatter((np.sin(data.phi_1 / 180 * np.pi) / data.n_down) ** 2, data.d_cos2_theta_up)
plt.show()

data.head(1)

columns = ['2 phi_1', '180 + psi (down)', '180 + psi (up)', 'n_down', 'd_n_down', 'n_up', 'd_n_up']
data[columns].to_csv('table.csv', sep=' ', index=False)

for direction in ['down', 'up']:
    columns = ['phi_1', f'psi_{direction}']
    plot2 = data[columns]
    plot2['d_phi_1'] = np.ones(plot2['phi_1'].shape[0]) * 0.25
    plot2['d_psi'] = np.ones((plot2.shape[0],)) * 0.5
    plot2.to_csv(f'plot2_{direction}.csv', sep=' ', index=False, header=False)

StarcoderdataJupyterScriptsDedupFiltered