In [7]:
#spiral arm parameters
function logisticFunction(x,x0, w)
return 1. / (1. + exp(-2. * (abs(x) - x0) / w));
end
function getRegularField(x,y,z)
kpc=1.0;muG=1.0;
# spiral arm parameters
pitch = 11.5 * pi*1.0 / 180;
sinPitch = sin(pitch);
cosPitch = cos(pitch);
tan90MinusPitch = tan(pi / 2 - pitch);
rArms=zeros(8)
rArms[1] = 5.1 * kpc;
rArms[2] = 6.3 * kpc;
rArms[3] = 7.1 * kpc;
rArms[4] = 8.3 * kpc;
rArms[5] = 9.8 * kpc;
rArms[6] = 11.4 * kpc;
rArms[7] = 12.7 * kpc;
rArms[8] = 15.5 * kpc;
# regular field parameters
bRing = 0.1 * muG;
hDisk = 0.40 * kpc;
wDisk = 0.27 * kpc;
bDisk=zeros(8)
bDisk[1] = 0.1 * muG;
bDisk[2] = 3.0 * muG;
bDisk[3] = -0.9 * muG;
bDisk[4] = -0.8 * muG;
bDisk[5] = -2.0 * muG;
bDisk[6] = -4.2 * muG;
bDisk[7] = 0.0 * muG;
bDisk[8] = 2.7 * muG;
bNorth = 1.4 * muG;
bSouth = -1.1 * muG;
rNorth = 9.22 * kpc;
rSouth = 17 * kpc;
wHalo = 0.20 * kpc;
z0 = 5.3 * kpc;
bX = 4.6 * muG;
thetaX0 = 49.0 * pi*1.0 / 180;
sinThetaX0 = sin(thetaX0);
cosThetaX0 = cos(thetaX0);
tanThetaX0 = tan(thetaX0);
rXc = 4.8 * kpc;
rX = 2.9 * kpc;
# striated field parameter
sqrtbeta = sqrt(1.36);
# turbulent field parameters
bDiskTurb=zeros(8)
bDiskTurb[1] = 10.81 * muG;
bDiskTurb[2] = 6.96 * muG;
bDiskTurb[3] = 9.59 * muG;
bDiskTurb[4] = 6.96 * muG;
bDiskTurb[5] = 1.96 * muG;
bDiskTurb[6] = 16.34 * muG;
bDiskTurb[7] = 37.29 * muG;
bDiskTurb[8] = 10.35 * muG;
bDiskTurb5 = 7.63 * muG;
zDiskTurb = 0.61 * kpc;
bHaloTurb = 4.68 * muG;
rHaloTurb = 10.97 * kpc;
zHaloTurb = 2.84 * kpc;
bx=0;by=0;bz=0;
r = sqrt(x * x + y * y); # in-plane radius
d = sqrt(x * x + y * y + z*z); # distance to galactic center
if (d < 1 ) || (d > 20 )
return [bx,by,bz]; # 0 field for d < 1 kpc or d > 20 kpc
end
phi=atan2(y,x); # phi de -pi a pi
sinPhi = sin(phi);
cosPhi = cos(phi);
lfDisk = logisticFunction(z, hDisk, wDisk);
# disk field
if r > 3
if r < 5
# molecular ring
bMag = bRing * (5 * kpc / r) * (1 - lfDisk);
bx += -bMag * sinPhi;
by += bMag * cosPhi;
else
# spiral region
i_0=7
for i=1:7
r11=rArms[i] *exp((phi*pitch-pi*pitch)) ;
r12=rArms[i+1] *exp((phi*pitch-pi*pitch)) ;
r21=rArms[i] *exp(((phi+2*pi)*pitch-pi*pitch));
r22=rArms[i+1] *exp(((phi+2*pi)*pitch-pi*pitch));
r31=rArms[i] *exp(((phi-2*pi)*pitch-pi*pitch));
r32=rArms[i+1] *exp(((phi-2*pi)*pitch-pi*pitch));
if (r>=r11 && r<r12) || (r>=r21 && r<r22) || (r>=r31 && r<r32)
i_0=i
end
end
bMag = bDisk[i_0];
bMag *= (5 * kpc / r) * (1 - lfDisk);
bx += bMag * (sinPitch * cos(phi) - cosPitch * sin(phi));
by += bMag * (sinPitch * sin(phi) + cosPitch * cos(phi));
end
end
# toroidal halo field
bMagH = exp(-abs(z) / z0) * lfDisk;
if (z >= 0)
bMagH *= bNorth * (1 - logisticFunction(r, rNorth, wHalo));
else
bMagH *= bSouth * (1 - logisticFunction(r, rSouth, wHalo));
end
bx += -bMagH * sinPhi;
by += bMagH * cosPhi;
# poloidal halo field
rc = rXc + abs(z) / tanThetaX0;
if (r < rc)
# varying elevation region
rp = r * rXc / rc;
bMagX = bX * exp(-1 * rp / rX) * (rp / r)^ 2.;
thetaX = atan2(abs(z), (r - rp));
if (z == 0)
thetaX = pi / 2.;
end
sinThetaX = sin(thetaX);
cosThetaX = cos(thetaX);
else
# constant elevation region
rp = r - abs(z) / tanThetaX0;
bMagX = bX * exp(-rp / rX) * (rp / r);
sinThetaX = sinThetaX0;
cosThetaX = cosThetaX0;
end
bx += sign(z) * bMagX * cosThetaX * cosPhi;
by += sign(z) * bMagX * cosThetaX * sinPhi;
bz += bMagX * sinThetaX;
#print ("JF2012 x=$x y=$y z=$z \n");
#print ("JF2012 Bx=$bx By=$by Bz=$bz \n");
return [bx,by,bz];
end
Out[7]:
In [8]:
function RG4_D(f,XX0)
N=50; # Cambiar este no. de iteraciones si se requiere
c=2.998*10^(8.); # Usamos MKS
kpc=3.0857*10^19.;
kyear=3.154*10^10.;
t=0;niter=0;
h=20*kpc/c/N/2.;
#Esta es la función Runge-Kutta con sus coeficientes
D=0;
while (D<20.0*kpc )
D=sqrt(XX0[1]*XX0[1]+XX0[2]*XX0[2]+XX0[3]*XX0[3]);
k1=f(XX0,t)
k2=f(XX0+(0.5*k1*h),t+.5*h)
k3=f(XX0+(0.5*k2*h),t+.5*h)
k4=f(XX0+(k3*h),t+h)
XX0=XX0 + ((h/6) * (k1 + 2 * (k2 + k3) + k4))
t += h;niter += 1;
end
return XX0,t,niter #Regresa un arreglo donde están
end
Out[8]:
In [9]:
function Propa_JF2012_simple(Charge,Mass,Energy, galactic_l, galactic_b,my_print)
# Entradas: Charge,Mass,Energy, galactic_l, galactic_b, my_print in degrees
######## Constantes
EeV=1.602*10^(-1.);
Mp=1.67*10^(-27.);
Qp=1.602*10^(-19.)
c=2.998*10^(8.);
kpc=3.0857*10^19.;
muGauss=10.^(-10.); # Tesla
kyear=3.154*10^10.;
############## Cond. iniciales
E0=Energy*EeV;# p=Ev/c^2
q=-Charge*Qp; # units of p charge, negative foe backpropagation
m0=Mass*Mp # units of p mass
############## Cond. iniciales
galactic_l_rad=deg2rad(galactic_l);
galactic_b_rad=deg2rad(galactic_b);
nx=cos(galactic_b_rad)*cos(galactic_l_rad);
ny=cos(galactic_b_rad)*sin(galactic_l_rad);
nz=sin(galactic_b_rad);
n_inicial=[nx,ny,nz]
n_inicial=n_inicial/norm(n_inicial);
v0=c*sqrt(1-m0^2c^4/E0^2)*n_inicial;
x0=-8.5 * kpc; # Sun position
y0=0;
z0=0;
#X0_noRel=[x0,y0,z0,v0[1],v0[2],v0[3]];
X0_Rel=[x0,y0,z0,E0*v0[1]/c^2,E0*v0[2]/c^2,E0*v0[3]/c^2];
function f_Rel(X,t)
# X es [x,y,z,px,py,pz]
# Distancia en m y tiempo en s
c=2.998*10^(8.);
M=m0*sqrt(1+(X[4]^2+X[5]^2+X[6]^2)/(m0*c)^2); # masa relativista
#M=E0/c^2 # masa constante
(Bx,By,Bz)=getRegularField(X[1]/kpc,X[2]/kpc,X[3]/kpc)*muGauss; # B en Teslas
return [X[4]/M,X[5]/M,X[6]/M,q/M*(X[5]*Bz-X[6]*By),q/M*(X[6]*Bx-X[4]*Bz),q/M*(X[4]*By-X[5]*Bx)]
end
xt, t_f, niter=RG4_D(f_Rel,X0_Rel);
x_f=xt[1]/kpc+8.5;
y_f=xt[2]/kpc;
z_f=xt[3]/kpc;
px=xt[4];
py=xt[5];
pz=xt[6];
p_f=[px,py,pz];
M=m0*sqrt(1+(px^2+py^2+pz^2)/(m0*c)^2);
E_f=M*c^2/EeV;
t_f=t_f/kyear;
R_f=sqrt(x_f^2+y_f^2+z_f^2);
alpha=rad2deg(acos(dot(n_inicial,p_f)/norm(p_f)));
l_f=rad2deg(atan2(py,px));
if l_f <0; l_f += 360.;end
b_f=rad2deg(atan2(pz,sqrt(px*px+py*py)));
if my_print ==1
@printf("niter \tE_f\tx_f\ty_f\tz_f\tl_f \tb_f \talpha \tt_f \n");
@printf("%d \t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\n",niter,E_f,x_f,y_f,z_f,l_f, b_f,alpha,t_f);
end
alpha=@sprintf("%.2f",alpha)
return alpha
end
Out[9]:
In [ ]:
#Using libraries
using PyPlot
#Include some code
include("JF2012_simple.jl")
include("RG4_D.jl")
include("Propa_JF2012_simple.jl")
#Features constantes
Charge=1 # 1 para p, 8 para O, 26 para Fe
Mass=1 # 1 para p, 16 para O, 56 para Fe
Energy=60 #Eev
N_lon=360
N_lat=180
granularity=1
latitude=N_lat*granularity
longitude=N_lon*granularity
Deviation=zeros(latitude,longitude)
print("Theta=")
#Propagation
for i= 1:latitude
if i%granularity==0
print(i/granularity,"-")
end
for j= 1:longitude
lat=90-(i-0.5)/granularity
lon=(j-0.5)/granularity-180.0
alpha = round(float(Propa_JF2012_simple(Charge,Mass,Energy,lon,lat,0)),12);
Deviation[i,j]=alpha
end
end
#Plot
PyPlot.figure(figsize=(12,7))
PyPlot.subplot(111,projection="hammer")
x=linspace(-pi,pi,longitude)
y=linspace(pi/2,-pi/2,latitude)
PyPlot.pcolormesh(x,y,Deviation,cmap="jet",shading="interp")
PyPlot.colorbar()
PyPlot.grid(true)
PyPlot.show()
In [ ]: