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ϕ][[1]];
  3. θhSol[θS_, a_] =
  4.   θh /. Solve[cosϕSol[θS, Sin[θh], Cos[θh], a] == 1, θh] /. {C[1] -> 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]][[1]], θhSol[θS, A[ES, mS]][[2]]], 0],
  10.    Min[Max[θhSol[θS, A[ES, mS]][[2]], θhSol[θS, A[ES, mS]][[1]]], 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ϕ][[1]];
  20.  
  21. θhSol[θS_,
  22.   a_] := θh /.
  23.    Solve[cosϕSol[θS, Sin[θh], Cos[θh], a] ==
  24.       1, θh] /. {C[1] -> 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]][[1]], θhSol[θS, A[ES, mS]][[2]]], 0],
  33.    Min[Max[θhSol[θS,
  34.        A[ES, mS]][[2]], θhSol[θS, A[ES, mS]][[1]]], 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[132]:= 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[132]:= 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[132]:= 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[132]:= General::stop: Further output of Solve::ifun will be suppressed during this calculation.
  50.  
  51. During evaluation of In[132]:= 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]][[2]]], 0],
  56.   Min[Max[θhSol[0.01, A[1500, 40]][[
  57.      2]], θhSol[0.01, A[1500, 40]][[1]]], π]}, {ϕ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]][[1]], θhSol[θS, A[ES, mS]][[2]]], 0],
  68.    Min[Max[θhSol[θS,
  69.        A[ES, mS]][[2]], θhSol[θS, A[ES, mS]][[1]]], 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[134]:= 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[134]:= 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[134]:= 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[134]:= 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]][[1]], θhSol[θS, A[ES, mS]][[2]]], 0],
  92.     Min[Max[θhSol[θS, A[ES, mS]][[2]], θhSol[θS, A[ES, mS]][[1]]], 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[[1]], θh0[[2]]},
  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[[1]], θh0[[2]]},
  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. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top