• API
• FAQ
• Tools
• Archive
SHARE
TWEET # Untitled a guest Jun 16th, 2019 54 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. cosϕSol[θS_, sinθh_, cosθh_, a_] =
2.   cosϕ /. Solve[cosϕ*Sin[θS] sinθh + Cos[θS]*cosθh == a, cosϕ][];
3. θhSol[θS_, a_] =
4.   θh /. Solve[cosϕSol[θS, Sin[θh], Cos[θh], a] == 1, θh] /. {C -> 0};
5. A[ES_, mS_] = Sqrt[4*ES^2 - 125^2]/Sqrt[4*ES^2 - 4*mS^2];
6. integral1[ES_, θS_, mS_] :=
7.  NIntegrate[1,
8.   {θh,
9.    Max[Min[θhSol[θS, A[ES, mS]][], θhSol[θS, A[ES, mS]][]], 0],
10.    Min[Max[θhSol[θS, A[ES, mS]][], θhSol[θS, A[ES, mS]][]], Pi]},
11.   {ϕh, 0,
12.    Min[Pi, ArcCos[cosϕSol[θS, Sin[θh], Cos[θh], A[ES, mS]]]]}]
13.
14. Clear[cosϕSol, θhSol, A, integral1]
15.
16. cosϕSol[θS_, sinθh_, cosθh_, a_] :=
17.   cosϕ /.
18.    Solve[cosϕ*Sin[θS] sinθh +
19.        Cos[θS]*cosθh == a, cosϕ][];
20.
21. θhSol[θS_,
22.   a_] := θh /.
23.    Solve[cosϕSol[θS, Sin[θh], Cos[θh], a] ==
24.       1, θh] /. {C -> 0}
25.
26. A[ES_, mS_] := Sqrt[4*ES^2 - 125^2]/Sqrt[4*ES^2 - 4*mS^2];
27.
28. integral1[ES_, θS_, mS_] :=
29.  NIntegrate[1,
30.   {θh,
31.    Max[Min[θhSol[θS,
32.        A[ES, mS]][], θhSol[θS, A[ES, mS]][]], 0],
33.    Min[Max[θhSol[θS,
34.        A[ES, mS]][], θhSol[θS, A[ES, mS]][]], Pi]},
35.   {ϕh, 0,
36.    Min[Pi, ArcCos[
37.      cosϕSol[θS, Sin[θh], Cos[θh],
38.       A[ES, mS]]]]},
39.   Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}]
40.
41. integral1[1500, 0.01, 40]
42.
43. (* During evaluation of In:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.
44.
45. During evaluation of In:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.
46.
47. During evaluation of In:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.
48.
49. During evaluation of In:= General::stop: Further output of Solve::ifun will be suppressed during this calculation.
50.
51. During evaluation of In:= NIntegrate::nlim: ϕh = Min[3.14159,3.14159 -0.545349 I] is not a valid limit of integration. *)
52.
53. (* NIntegrate[1, {θh,
54.   Max[Min[θhSol[0.01, A[1500, 40]][[
55.      1]], θhSol[0.01, A[1500, 40]][]], 0],
56.   Min[Max[θhSol[0.01, A[1500, 40]][[
57.      2]], θhSol[0.01, A[1500, 40]][]], π]}, {ϕh, 0,
58.    Min[π,
59.    ArcCos[cosϕSol[0.01, Sin[θh], Cos[θh],
60.      A[1500, 40]]]]},
61.  Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}] *)
62.
63. integral1[ES_, θS_, mS_] :=
64.  NIntegrate[1,
65.   {θh,
66.    Max[Min[θhSol[θS,
67.        A[ES, mS]][], θhSol[θS, A[ES, mS]][]], 0],
68.    Min[Max[θhSol[θS,
69.        A[ES, mS]][], θhSol[θS, A[ES, mS]][]], Pi]},
70.   {ϕh, 0,
71.    Min[Pi, Re@
72.      ArcCos[cosϕSol[θS, Sin[θh], Cos[θh],
73.        A[ES, mS]]]]},
74.   Method -> {"GlobalAdaptive", "SymbolicProcessing" -> 0}]
75.
76. integral1[1500, 0.01, 40]
77.
78. (* During evaluation of In:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.
79.
80. During evaluation of In:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.
81.
82. During evaluation of In:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information.
83.
84. During evaluation of In:= General::stop: Further output of Solve::ifun will be suppressed during this calculation. *)
85.
86. (* 0.0981353 *)
87.
88. integral1[ES_, θS_, mS_] :=
89.   NIntegrate[1,
90.    {θh,
91.     Max[Min[θhSol[θS, A[ES, mS]][], θhSol[θS, A[ES, mS]][]], 0],
92.     Min[Max[θhSol[θS, A[ES, mS]][], θhSol[θS, A[ES, mS]][]], Pi]},
93.    {ϕh, 0,
94.     Min[Pi, ArcCos[Piecewise[{{#, -1 <= # <= 1}}, 1.] &@
95.        cosϕSol[θS, Sin[θh], Cos[θh], A[ES, mS]]]
96.      ]}];
97.
98. integral1[1500, 0.01, 40]
99.
100. (*  0.0289182  *)
101.
102. integral1[ES_, θS_, mS_] :=
103.   With[{θh0 = Sort@Clip[θhSol[θS, A[ES, mS]], {0, Pi}]},
104.    NIntegrate[1,
105.     {θh, θh0[], θh0[]},
106.     {ϕh,0,
107.      ArcCos[Piecewise[{{#, -1 <= # <= 1}}, 1.] &@
108.        cosϕSol[θS, Sin[θh], Cos[θh], A[ES, mS]]
109.       ]},
110.     MinRecursion -> 2]
111.    ];
112.
113. integral1[1500, 0.01, 40]
114. (*  0.0289182  *)
115.
116. ClearAll[realACos];  (* another way to code a real arc cosine *)
117. realACos[x_?NumericQ /; -1 <= x <= 1] := ArcCos[x];
118. realACos[_?NumericQ] := 0;
119.
120. integral2[ES_, θS_, mS_] :=
121.   With[{θh0 = Sort@Clip[θhSol[θS, A[ES, mS]], {0, Pi}]},
122.    NIntegrate[
123.     realACos[cosϕSol[θS, Sin[θh], Cos[θh], A[ES, mS]]],
124.     {θh, θh0[], θh0[]},
125.     MinRecursion -> 5, MaxRecursion -> 20,
126.     Method -> "GaussKronrodRule"]
127.    ];
128.
129. integral2[1500, 0.01, 40]
130. (*  0.0289182  *)
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.
Not a member of Pastebin yet?