https://hal.archives-ouvertes.fr/hal-02535772
Revision 5070d8d89ae3101cc81fd626379cd43f94843f2a authored by Software Heritage on 07 April 2020, 15:31 UTC, committed by Software Heritage on 07 April 2020, 15:31 UTC
0 parent
Tip revision: 5070d8d
gii-a_G1986.py
``````import numpy as np
import matplotlib.pyplot as plt

def g0(x):
""" Eq.(6) of 1986A&A...160..195G """
return 1./( x + np.sqrt(x*x + 4./np.pi) )

def gc(xp,x):
""" Eq.(5) of 1986A&A...160..195G """
# symmetry relation
if (xp <= 0.):
xp=abs(xp)
if (x <=0.):
x=abs(x)
#----------
if (xp <= x):
return g0(x)
else:
return np.exp(x*x - xp*xp)*g0(xp)

def ierfc(x):
""" Eq.(9) of 1986A&A...160..195G """
return np.exp(-x*x)*(1.-2*x*g0(x))/np.sqrt(np.pi)

def gw0(xp,x):
""" Eq.(8a) of 1986A&A...160..195G """
return ierfc(0.5*abs(xp-x))

def gw1(xp,x):
""" Eq.(8B) of 1986A&A...160..195G """
return (2.-xp/x)*ierfc(0.5*abs(xp-x))

def psic(x):
""" Doppler profile, after Eq.(10) of 1986A&A...160..195G """
return np.exp(-x*x)/np.sqrt(np.pi)

def psiw(a,x):
""" Lorentzian profile, after Eq.(10) of 1986A&A...160..195G """
return a/(a*a+x*x)/np.pi

def gtrans(apr,a,xp,x):
""" Eq.(10) of 1986A&A...160..195G """
rpsi=psic(x)/(psic(x)+psiw(a,x))
if (apr==0):
return rpsi*gc(xp,x) + (1.-rpsi)*gw0(xp,x)
if (apr==1):
return rpsi*gc(xp,x) + (1.-rpsi)*gw1(xp,x)

def giia(apr,a,xp,x):
if (abs(x) < 2):
return gc(xp,x)
if (2 <= abs(x) <= 4):
return gtrans(apr,a,xp,x)
if (abs(x) > 4):
if (apr==0):
return gw0(xp,x)
if (apr==1):
return gw1(xp,x)

dx=0.1
xp=np.arange(-10.,10.+dx,dx)

# may be 0 or 1 (see Gouttebroze 1986)
apr=1

a=0.001
# for G86 figures' reproduction

plt.figure(1)
for x in [-5.,0., 2., 4., 6., 8.]:
g2=np.zeros((len(xp)),'Float64')
for k in np.arange(len(xp)):
xx=xp[k]
# how to make this vector, passing all xp's once?
g2[k]=giia(apr,a,xx,x)
# --------------------
plt.plot(xp,g2)
plt.xlabel("x'",fontsize=18)
plt.ylabel("g_II-A(x',x)",fontsize=18)
plt.xlim(left=-7.5,right=10.)
plt.ylim(top=1.0,bottom=0.)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)

plt.figure(2)
for x in [0., 1., 1.5, 2.,2.5,3.,3.5, 4.]:
#for x in [0., -1., -1.5, -2.,-2.5,-3.,-3.5, -4.]:
g2=np.zeros((len(xp)),'Float64')
for k in np.arange(len(xp)):
xx=xp[k]
# how to make this vector, passing all xp's once?
g2[k]=giia(apr,a,xx,x)
# --------------------
plt.plot(xp,g2)
plt.xlabel("x'",fontsize=18)
plt.ylabel("g_II-A(x',x)",fontsize=18)
plt.xlim(left=-5.,right=5.)
plt.ylim(top=1.0,bottom=0.)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)

plt.show()

``````

Computing file changes ...