85 lines
3.0 KiB
Python
85 lines
3.0 KiB
Python
import math
|
|
import numpy as np
|
|
from scipy.optimize import fsolve
|
|
|
|
|
|
def antiparallel2gap(K, phi, undudict):
|
|
gLH = K2gap(K, undudict['K-value_LH'])
|
|
if phi >= 0.0:
|
|
gLV = K2gap(K, undudict['K-value_apLV+'])
|
|
gC = K2gap(K, undudict['K-value_45+'])
|
|
dgLV = gLV - gLH
|
|
dgC = gC - gLH - dgLV/2
|
|
else:
|
|
gLV = K2gap(K, undudict['K-value_apLV-'])
|
|
gC = K2gap(K, undudict['K-value_45-'])
|
|
dgLV = gLV - gLH
|
|
dgC = gC - gLH - dgLV/2
|
|
return gLH + dgLV * np.sin(0.5 * phi)**2 + dgC * np.sin(phi)**2
|
|
|
|
|
|
def K2gap(Kval, fitparam):
|
|
g2K_func = np.poly1d(fitparam[::-1])
|
|
tau_init = 1.0
|
|
k_log = float(np.log(Kval))
|
|
return float(fsolve(k_log - g2K_func, tau_init))
|
|
|
|
|
|
|
|
## TODO: implement the proper model:
|
|
|
|
|
|
def antiparallel_rad2k_fit(radial, phi_in_degree, undudict):
|
|
key = "anti-parallel-" if phi_in_degree < 0 else "anti-parallel+"
|
|
undudict = undudict[key]["K-value"]
|
|
phi = math.radians(phi_in_degree)
|
|
amp1 = getfit(radial, undudict['fitpars_amp1'])
|
|
amp2 = getfit(radial, undudict['fitpars_amp2'])
|
|
amp3 = getfit(radial, undudict['fitpars_amp3'])
|
|
KLH = np.exp(getfit(radial, undudict['fitpars_KLH']))
|
|
KLV = np.exp(getfit(radial, undudict['fitpars_KLV']))
|
|
K45 = np.exp(getfit(radial, undudict['fitpars_K45']))
|
|
dKLV = KLV - KLH
|
|
dK45 = K45 - KLH - dKLV/2
|
|
#print(radial, phi_in_degree, amp1, amp2, amp3, KLH, dKLV, dK45)
|
|
if phi >= 0.0:
|
|
return KLH + dKLV * np.sin(0.5*phi)**2 + dK45 * np.sin(phi)**2 + amp1 * np.sin(2*phi) + amp2 * np.sin(2.0*phi)**2 + amp3 * np.cos(6*phi)
|
|
else:
|
|
return KLH + dKLV * np.sin(0.5*phi)**2 + dK45 * np.sin(phi)**2 + amp1 * np.sin(2*phi) + amp2 * np.sin(2.0*phi)**2 + amp3 * np.cos(8*phi)
|
|
|
|
|
|
def getfit(val, fitparam):
|
|
fit_func = np.poly1d(fitparam)
|
|
val = float(val)
|
|
return fit_func(val)
|
|
|
|
|
|
def antiparallel_k2rad_fit(Kval, phi_in_degree, undudict):
|
|
key = "anti-parallel-" if phi_in_degree < 0 else "anti-parallel+"
|
|
undudict = undudict[key]["K-value"]
|
|
phi = math.radians(phi_in_degree)
|
|
amp1 = getfit(Kval, undudict['fitpars_amp1'])
|
|
amp2 = getfit(Kval, undudict['fitpars_amp2'])
|
|
amp3 = getfit(Kval, undudict['fitpars_amp3'])
|
|
RadLH = getfit_inverse(np.log(Kval), undudict['fitpars_KLH'])
|
|
RadLV = getfit_inverse(np.log(Kval), undudict['fitpars_KLV'])
|
|
Rad45 = getfit_inverse(np.log(Kval), undudict['fitpars_K45'])
|
|
dRadLV = RadLV - RadLH
|
|
dRad45 = Rad45 - RadLH - dRadLV/2
|
|
#print(radial, phi_in_degree, amp1, amp2, amp3, RadLH, RadLV, Rad45)
|
|
if phi >= 0.0:
|
|
return RadLH + dRadLV * np.sin(0.5*phi)**2 + dRad45 * np.sin(phi)**2 + amp1 * np.sin(2*phi) + amp2 * np.sin(2.0*phi)**2 + amp3 * np.cos(6*phi)
|
|
else:
|
|
return RadLH + dRadLV * np.sin(0.5*phi)**2 + dRad45 * np.sin(phi)**2 + amp1 * np.sin(2*phi) + amp2 * np.sin(2.0*phi)**2 + amp3 * np.cos(8*phi)
|
|
|
|
|
|
def getfit_inverse(val, fitparam):
|
|
fit_func = np.poly1d(fitparam)
|
|
tau_init = 1.0
|
|
val = float(val)
|
|
return float(fsolve(val - fit_func, tau_init))
|
|
|
|
|
|
|
|
|