# -*- coding: utf-8 -*-
"""
Created on Mon Aug 16 11:49:41 2021

@author: Joschka
"""
import matplotlib.pyplot as plt
import numpy as np
from src import B_field_calculation as bf

from src import coil_class as BC


from IPython import get_ipython
get_ipython().run_line_magic('matplotlib', 'qt')
#get_ipython().run_line_magic('matplotlib', 'inline')

#set up axis
x_m = np.linspace(-0.05, 0.05, 51)
z_m = np.linspace(-0.05, 0.05, 201)


z = z_m*1e3
x = x_m*1e3 #for plotting in mm

#Import Values from simulation
B_z_sim = np.loadtxt('data/B_z_HH2.txt')
B_x_sim = np.loadtxt('data/B_x_HH2.txt')


################# My simulation #########################
I = 5
HH = 1
d_coils = 44
R_inner = 44-3*1.7

layers = 6
windings = 2
wire_width = 1.7
wire_height = 2.6




B_z, B_x = bf.B_multiple_raster_test(I,HH,R_inner,d_coils,layers,windings,wire_width, wire_height, x_m,z_m)
#B_test = B_field_ideal_AHH(layers*windings,I,R_inner*1e-3,d_coils*1e-3,z_m)

#B_x = np.concatenate((-np.flip(B_r),B_r[1:len(B_r)]))









#Calculate gradients/curvature
B_z_sim_grad = np.gradient(np.gradient(B_z_sim,z_m),z_m)/1e4
B_x_sim_grad = np.gradient(B_x_sim,x_m)/100
B_z_grad = np.gradient(np.gradient(B_z,z_m),z_m)/1e4
B_x_grad = np.gradient(B_x,x_m)/100


#Calculate relative differences in permille
rel_diff_Bz = (B_z-B_z_sim)/B_z
rel_diff_Bx = (B_x-B_x_sim)/B_x
rel_diff_Bz_grad = (B_z_grad-B_z_sim_grad)/B_z_grad

rel_diff_Bz_grad_mean = (B_z_grad-B_z_sim_grad)/np.mean(B_z_grad)
rel_diff_Bx_grad = (B_x_grad-B_x_sim_grad)/B_x_grad

#Plotting
plt.figure(1,figsize=(20,18))

plt.rcParams.update({'font.size': 15})
plt.suptitle("Helmholtz coil field Bz along z-axis, comparison of simulations", fontsize=30)


#Field plot
##########################
plt.subplot(3,2,1)
plt.plot(z,B_z,linestyle = "solid", label = r"$Bz$: Result via elliptic integrals")
plt.plot(z,B_z_sim,linestyle = "dashdot", label = r"$B_{z, sim}$: Numerical Matlab simulation")
plt.plot(z,(B_z-B_z_sim), label = r"$Bz - B_{z, sim}$")
#plt.xlim(-0.01,0.01)
plt.title("B-field" ,fontsize = 30)

plt.ylabel(r"$Bz$ [G]")
plt.xlabel("z-axis [mm]")
plt.legend()


#############################
plt.subplot(3,2,3)
plt.plot(z,(B_z-B_z_sim), label = r"$Bz - B_{z, sim}$")
plt.ylabel("absolute deviation [G]")
plt.xlabel("z-axis [mm]")
plt.legend()

#############################
plt.subplot(3,2,5)
plt.plot(z,1000*rel_diff_Bz, label = "$(Bz - B_{z, sim}) / Bz$")
plt.ylabel("relative deviation [‰]")
plt.xlabel("z-axis [mm]")
plt.legend()


######################Gradient plot############################

################
plt.subplot(3,2,2)
plt.plot(z,B_z_grad,linestyle = "solid", label = r"$\nabla_z^2 Bz$: Result via elliptic integrals")
plt.plot(z,B_z_sim_grad,linestyle = "dashdot", label = r"$\nabla_z^2 B_{z, sim}$: Numerical Matlab sim.")
plt.plot(z,(B_z_grad-B_z_sim_grad), label = r"$\nabla_z^2 Bz - \nabla_z^2 B_{z, sim}$")

plt.ylabel(r"$\nabla_z^2 Bz [G/cm^2]$")
plt.xlabel("z-axis [mm]")
plt.title("Curvature of B-field",fontsize = 30)
plt.legend(loc='lower right')


#################

plt.subplot(3,2,4)
plt.plot(z,(B_z_grad-B_z_sim_grad), label = r"$\nabla_z^2 Bz - \nabla_z^2 B_{z, sim}$")
plt.ylabel(r"absolute deviation $[G/cm^2]$")
plt.xlabel("z-axis [mm]")
plt.legend()

#####################
plt.subplot(3,2,6)
plt.plot(z,1000*rel_diff_Bz_grad, label = r"$(\nabla_z^2 Bz - \nabla_z^2 B_{z, sim}) / \nabla_z^2 Bz$")
#plt.ylim(-57,10)
plt.ylabel("relative deviation [‰]")
plt.xlabel("z-axis [mm]")
plt.legend()


plt.savefig("output/HH_benchmark_5A_6x2.pdf")
plt.show()

############### relative deviation with averaging by the mean not the individual value ########################################
plt.figure(2)

plt.plot(z,1000*rel_diff_Bz_grad_mean, label = r"$(\nabla_z^2 Bz - \nabla_z^2 B_{z, sim}) / mean(\nabla_z^2 Bz)$")
#plt.ylim(-57,10)
plt.ylabel("relative deviation [‰]")
plt.xlabel("z-axis [mm]")
plt.legend()
plt.savefig("output/HH_benchmark_5A_6x2_rel_deviation_via_mean.pdf")
plt.show()

##################### x-Axis #########################################################

plt.figure(3)

plt.rcParams.update({'font.size': 15})
plt.suptitle("Helmholtz coil field B_x along x-axis, comparison of simulations", fontsize=30)


#Field plot
##########################

plt.plot(x,B_x,linestyle = "solid", label = r"$B_x$: Result via elliptic integrals")
plt.plot(x,B_x_sim,linestyle = "dashdot", label = r"$B_{x, sim}$: Numerical Matlab simulation")
plt.plot(x,(B_x-B_x_sim), label = r"$B_x - B_{x, sim}$")
#plt.xlim(-0.01,0.01)
plt.title("B-field" ,fontsize = 30)

plt.ylabel(r"$B_x$ [G]")
plt.xlabel("x-axis [mm]")
plt.legend()




#################


plt.savefig("output/HH_benchmark_5A_6x2_x-axis.pdf")
plt.show()