157 lines
4.8 KiB
Python
157 lines
4.8 KiB
Python
|
__author__ = 'markusro'
|
||
|
|
||
|
from numpy import *
|
||
|
import mathlib
|
||
|
|
||
|
import scipy.special
|
||
|
import scipy.integrate
|
||
|
|
||
|
numpat = re.compile('\d+\.\d+')
|
||
|
NUM_PROCESSES = 8
|
||
|
|
||
|
#define norm1(a,b) exp(gammln(b/a))/(a*pow(b/a,b/a))
|
||
|
|
||
|
#define glntau1(x,t,a,b) exp( -exp( (log(x)-log(t))*a )*b/a) * pow(x/t,b)
|
||
|
|
||
|
|
||
|
def filon(oms, x,y):
|
||
|
amps = N.zeros(len(oms), dtype='complex')
|
||
|
for i,om in enumerate(oms):
|
||
|
amps[i] = sum(diff(y)/diff(x)*(cos(om*x[1:])-cos(om*x[:-1])))/om**2
|
||
|
amps[i] += 1j*(y[0]/om + N.sum(diff(y)/diff(x)*(sin(om*x[1:])-sin(om*x[:-1])))/om**2)
|
||
|
return amps
|
||
|
|
||
|
|
||
|
class Yaff:
|
||
|
def __init__(self, dist_type=0):
|
||
|
self.dist_x = N.logspace(-10,4,512)
|
||
|
self.dist_y = N.zeros(self.dist_x.size)
|
||
|
self.dist_type = [self.yaff,
|
||
|
self.yaff_gg,
|
||
|
][dist_type]
|
||
|
|
||
|
|
||
|
|
||
|
def gg(self, p, tau):
|
||
|
tau0, a, b = p
|
||
|
"""
|
||
|
Generalized Gamma function
|
||
|
"""
|
||
|
NGG = a * (b/a)**(b/a) / scipy.special.gamma(b/a)
|
||
|
#g = exp(-b/a*exp( (log(tau)-log(tau0))*a))*(tau/tau0)**b
|
||
|
g = exp(-b/a*exp( (log(tau)-log(tau0))*a))*(tau/tau0)**b
|
||
|
return g*NGG
|
||
|
|
||
|
|
||
|
def ggb(self, p, tau):
|
||
|
tau0,a,b = p
|
||
|
norm_ggb = a*(1+b)/pi *b**(b/(1+b)) * sin(pi*b/(1+b))
|
||
|
g = (b*exp( (log(tau)-log(tau0)) *a) + exp( (log(tau)-log(tau0))*(-a*b)))**(-1)
|
||
|
return norm_ggb*g
|
||
|
|
||
|
|
||
|
def gg_hw(self, p ,tau):
|
||
|
tau0, a, b, g, s = p
|
||
|
"""
|
||
|
Generalized Gamma function with HF wing.
|
||
|
s gives the onset of the wing in units of tau0, i.e. 100 means 100*tau0 or two decades later
|
||
|
g is the slope of the wing, smaller than b; if not the slope will be b
|
||
|
"""
|
||
|
NGG = a * (b/a)**(b/a) / ( scipy.special.gamma(b/a) + s**(g-b)*(a/b)**((g-b)/a)* scipy.special.gamma(g/a) )
|
||
|
g = exp(-b/a*exp((log(tau) - log(tau0))*a))*(tau/tau0)**b * (1 + (tau*s/tau0)**(g-b))
|
||
|
return g*NGG
|
||
|
|
||
|
|
||
|
def dist_to_eps(omegas,dist_x,dist_y):
|
||
|
epp = zeros(len(omegas), dtype='complex')
|
||
|
for i,om in enumerate(omegas):
|
||
|
kern_imag = dist_y*om*dist_x/(1+(om*dist_x)**2)
|
||
|
kern_real = dist_y/(1+(om*dist_x)**2)
|
||
|
epp[i] = scipy.integrate.simps(kern_real,log(dist_x))
|
||
|
epp[i] += 1j*scipy.integrate.simps(kern_imag,log(dist_x))
|
||
|
return epp
|
||
|
|
||
|
def dist_to_relax(self, ts,dist_x,dist_y):
|
||
|
phi = zeros(len(ts))
|
||
|
for i,t in enumerate(ts):
|
||
|
phi[i] = phi_t(t,dist_x,dist_y)
|
||
|
return phi
|
||
|
|
||
|
|
||
|
def phi_t(self, t, dist_x, dist_y):
|
||
|
kern = dist_y*exp(-t/dist_x)
|
||
|
phi = scipy.integrate.simps(kern,log(dist_x))
|
||
|
return phi
|
||
|
|
||
|
|
||
|
def williams_watts(self, phi_1, phi_2, lambd):
|
||
|
return phi_1 * ( (1-lambd) + lambd*phi_2)
|
||
|
|
||
|
def yaff2(self, p, x):
|
||
|
"""
|
||
|
GGe+G b ---> GG+GB - lambda determined by p
|
||
|
"""
|
||
|
delta_eps, s_c, tau1,a1,b1, g,s, tau2, a2, b2 = p
|
||
|
p_ggb = tau2,a2,b2
|
||
|
p_gghw=tau1,a1,b1,g,s
|
||
|
lambd = 2*g/pi*a1*(b1/a1)**(b1/a1)
|
||
|
lambd /= scipy.special.gamma(b1/a1)*s**(b1-g) + (a1/b1)**((g-b1)/a1)*scipy.special.gamma(g/a1)
|
||
|
print "LAMBDA:",lambd,a1,b1,g,s
|
||
|
dist_ggb = self.ggb(p_ggb, dist_t)
|
||
|
phi_beta = self.dist_to_relax(ts, dist_t, dist_ggb ).real
|
||
|
phi_alpha = self.dist_to_relax(ts, dist_t, self.gg_hw(p_gghw, dist_t)).real
|
||
|
phi_tot = (1-lambd) + lambd*phi_beta
|
||
|
phi_tot *= phi_alpha
|
||
|
epp = delta_eps*2*pi*x*filon(2*pi*x, ts, phi_tot ).real
|
||
|
epp += s_c/(2*pi*x)
|
||
|
return epp
|
||
|
|
||
|
def yaff(self, p,x, cb=None):
|
||
|
"""
|
||
|
GG + GB
|
||
|
"""
|
||
|
delta_eps, tau1, a1, b1, lambd, tau2, a2, b2 = p
|
||
|
|
||
|
dist_ggb = ggb(self.dist_x,tau2,a2,b2)
|
||
|
dist_gg = gg(self.dist_x,tau1,a1,b1)
|
||
|
|
||
|
ts = logspace(-10,5,1024)
|
||
|
phi_beta = self.dist_to_relax(ts, self.dist_x, dist_ggb).real
|
||
|
phi_alpha = self.dist_to_relax(ts, self.dist_x, dist_gg).real
|
||
|
# William-Watts-Ansatz
|
||
|
phi_tot = (1-lambd) + lambd*phi_beta
|
||
|
phi_tot *= phi_alpha
|
||
|
epp = delta_eps*2*pi*x*filon(2*pi*x, ts, phi_tot ).real
|
||
|
|
||
|
if cb != None:
|
||
|
cb.next(p,x,epp)
|
||
|
return epp
|
||
|
|
||
|
def yaff_gg(self, p,x, cb=None):
|
||
|
"""
|
||
|
GG
|
||
|
"""
|
||
|
delta_eps, s_c, tau1,a1,b1 = p
|
||
|
|
||
|
dist_gg = gg(dist_t,tau1,a1,b1)
|
||
|
|
||
|
phi_alpha = dist_to_relax(ts, dist_t, dist_gg).real
|
||
|
epp = delta_eps*2*pi*x*filon(2*pi*x, ts, phi_alpha ).real
|
||
|
# epp = dist_to_eps(2*pi*x, dist_t, dist_gg ).real
|
||
|
epp += s_c/(2*pi*x)
|
||
|
if cb != None:
|
||
|
cb(p,x,epp)
|
||
|
return epp
|
||
|
|
||
|
|
||
|
def show_correlations(cov, ax):
|
||
|
cor = zeros(cov.shape)
|
||
|
for i in range(cor.shape[0]):
|
||
|
for j in range(cor.shape[1]):
|
||
|
if cov[i,i] == 0 or cov[j,j]==0:
|
||
|
cor[i,j] = 0
|
||
|
else:
|
||
|
cor[i,j] = cov[i,j]/sqrt(cov[i,i]*cov[j,j])
|
||
|
ax.matshow(abs(cor), cmap=cm.RdYlGn)
|
||
|
ax.colorbar()
|