Quit[];

En este notebook se calcula numéricamente la sección eficaz diferencial
para el pozo esférico de potencial.

Definición de los ángulos que aparecen en el problema.

(*Ángulo de incidencia.*)
α[b_] :=ArcSin[b];

(*Ángulo de transmisión.*)
β[b_, n_]:= ArcSin[(b/n)];

(*Ángulo de dispersión.*)
αinf[b_, n_] := π + 2α[b] - 2β[b, n];

1) Resultado de simular las posiciones finales de las Npar partículas
 sobre una pantalla esférica de radio R para un índice de refracción npr.

(*Radio de la pantalla. Experimenten con distintos valores de R
 para ver el efecto de no estar en la región asintótica de las trayectorias.*)
R = 2;

(*Índice de refracción.*)
npr = 1.5;

(*Número de partículas.*)
Npar = 50000;

(*Se elijen al azar los ángulos azimutales*)
ϕpr = 2. π RandomReal[{0, 1}, Npar];

(*Esta fórmula cálcula el ángulo polar sobre la pantalla. Tiene en cuenta que la pantalla no está en el infinito. 
Cuanto menor sea R, mayores serán los efectos de proyectar sobre una pantalla a distancia finita.*)
αR = αinf[bes = (Sqrt[ RandomReal[{0, 1}, Npar]]), npr] - ArcSin[bes/R];

(*Esta función calcula la dirección en coordenadas cartesianas en función de los ángulos esféricos.*)
coord[θ_, ϕ_]:= ({Cos[θ],0,0} + Sin[θ]{0,Sin[ϕ], Cos[ϕ]})
SetAttributes[coord, Listable]; 

(*Se calculan las posiciones de las partículas dispersadas sobre la pantalla.*)
puntos = R coord[αR, ϕpr];

(*La pantalla es una esfera de radio R (en realidad 0.99 R, para que los puntos no queden justo sobre la esfera.*)
esferaR = Graphics3D[{Orange,Opacity[1], Specularity[White, 2000], Sphere[{0,0, 0},0.99 R]}];

(*Se grafican los puntos junto con la pantalla*)
Show[
esferaR,
Graphics3D[{White, Opacity[1],AbsolutePointSize[0],Point[puntos]}], 

ViewVertical-> {0,0,1}, 
ImageSize-> 500, 
PlotRange-> All, 
Lighting-> {{"Point", GrayLevel[0.4], {-10R, 0, 0}}, {"Ambient", GrayLevel[0.4]}},
ViewVector-> {10{-4, -3, 2}, {0, 0, 0}}, 
ViewVertical-> {0,0,1},
Background-> Black]

2) Cálculo numérico de ((dσ)/(dΩ)).

(*Se simula el resultado de arrojar npar partículas. En Theta se guardan los npar ángulos de dispersión.*)
npr =1.5;
Npar = 1000000;
Theta = αinf[(Sqrt[RandomReal[{0,1}, Npar]]), npr]-π;

(*Se cuentan las partículas dispersadas dentro de cada intervalo angular de ancho Δ, entre 0 y el valor máximo en Theta.*)
tbin = BinCounts[Theta, {0, Max[Theta], Δ = 0.01}];

(*Se calcula la expresión aproximada de dσ/dΩ. Se toma como punto de evaluación el centro de cada intervalo.*)
σExp = Table[{θ = i Δ-(Δ/2), ((tbin[[i]])/(2π(Npar/π) Δ Sin[θ]))}, {i, 1, Length@tbin}];

(*Gráfico de σdata.*)
g1 = ListPlot[σExp, PlotRange-> All, PlotStyle->Orange]

(*Gráfico de la sección eficaz calculada analíticamente.*)
g2 = Plot[(((npr^2)(npr - Cos[θ/2])(npr Cos[θ/2]-1))/(4Cos[(θ/2)](((1 + (npr^2)-2npr Cos[(θ/2)]))^2))), {θ, 0, π- 2 ArcSin[(1/npr)]}, PlotRange-> All]

(*Comparación del resultado numérico con la curva teórica.*)
graf = Show[g2, g1,  AxesOrigin-> {0,0}, ImageSize-> 600, PlotRange-> All]