# -*- coding: utf-8 -*-
"""
Created on Fri Aug 27 15:14:48 2021

@author: Joschka
"""
from src import physical_constants as cs
import numpy as np


mu = 9.9* cs.mu_B/8

Grad_Bz = cs.m_Dy_164 * 9.81/(8*mu)

print("For levitation:")
print(f"dBz/dz = {Grad_Bz*1e4*1e-2:.4f} G/cm")
print("")

T = 1e-6
dt = 10e-3
sigma = np.sqrt(cs.k_B*T/cs.m_Dy_164) * dt

dz = 2 * sigma
print(f"cloud size: {sigma*1e3:.3f} mm")

#dz = 250e-6


Grad_Bz = 2 * dz * cs.m_Dy_164/(dt**2 * mu)


print(" ")
print("For Stern-Gerlach separation:")
print(f"dBz/dz = {Grad_Bz*1e4*1e-2:.4f} G/cm")
print(" ")

def cloud_pos(mF, dt = 10e-3):
    a = mF*mu*Grad_Bz/cs.m_Dy_164 - 9.81
    s = 0.5 * a *dt**2
    return s
a8 = -8*mu*Grad_Bz/cs.m_Dy_164 - 9.81
a8_neg = 8*mu*Grad_Bz/cs.m_Dy_164 - 9.81
a7 = -7*mu*Grad_Bz/cs.m_Dy_164 - 9.81
s8 = 0.5 * a8 * dt**2
s8_neg = 0.5 * a8_neg * dt**2
s7 = 0.5 * a7 * dt**2

for mF in np.arange(-8,9,1):
    print(f"for mF = {mF}: pos = {cloud_pos(mF) * 1e3} mm")


print(f"Total fall for m8 for acceleration due to magnetic field in same direction than gravity = {s8 * 1e3}mm")
print(f"Total lift for m8 for acceleration due to magnetic field in other direction than gravity = {s8_neg * 1e3}mm")
print(f"Fall without magnetic field = {-0.5 * 9.81 * dt**2 *1e3} mm")


print(f"s7 = {s7 * 1e3}mm")
print(f"diff = {(s8 - s7) * 1e3} mm")

# %%
# For thermal cloud size from paper

dt = 15e-3
sigma = 0.5e-3

T = 3.2e-6

sigma = np.sqrt(cs.k_B*T/cs.m_Dy_164) * dt


dz = 2 * sigma
print(f"cloud size: {sigma*1e3:.3f} mm")

#dz = 250e-6


Grad_Bz = 2 * dz * cs.m_Dy_164/(dt**2 * 1.241 * cs.mu_B)
print(" ")
print("For Stern-Gerlach separation:")
print(f"dBz/dz = {Grad_Bz*1e4*1e-2:.4f} G/cm")
print(" ")

# %%
dt1 = 10e-3
dt2 = 10e-3
dt = dt1 + dt2
sigma = np.sqrt(0.3e-3**2 + (np.sqrt(cs.k_B*T/cs.m_Dy_164) * dt)**2)

T = 3.2e-6

#sigma = np.sqrt(cs.k_B*T/cs.m_Dy_164) * dt

def cloud_sep(Grad_Bz, dt1 = 10e-3, dt2 = 0):
    Grad_Bz*=1e-2
    a = cs.mu_B * 1.241*Grad_Bz/cs.m_Dy_164
    s1 = 0.5 * a *dt1**2
    v1 = a * dt1
    s2 = s1 + v1 * dt2
    s3 = 0.5 * 9.81 *(dt1+dt2)**2
    #print(f"fall:{s3*1e3}mm")
    return s2*1e3

print(f"2 x sigma = {2*sigma*1e3}")
print(cloud_sep(8, dt1 = dt1, dt2 = dt2))

print(cloud_sep(30, dt1 = dt1, dt2 = dt2))