Guest User

Untitled

a guest
Nov 20th, 2017
105
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 2.21 KB | None | 0 0
  1. ne = 10;
  2. nl = 20;
  3. nL = 2*nl;
  4. nst = 2 (2 nl + 1);
  5. mvec1 = Table[-(nl - i), {i, 0, nL}]
  6. mvec = Join[mvec1, mvec1]
  7.  
  8. Lfxn[[Mu]_] := Mod[[Mu] - 1, nL] + 1
  9. sfxn[[Mu]_] := Sign[[Mu] - ( nst/2 + 0.5)]
  10. kdfxn[i_, j_] := If[i == j, 1, 0]
  11. Lvec = Table[i, {i, 1, nL}];
  12. svec = Table[sfxn[i[Mu]], {i[Mu], 1, nst}];
  13. avec = Table[kdfxn[i, j], {i, 1, ne}, {j, 1, nst}];
  14.  
  15. Vc = Block[{l, L, m, m1, p},
  16. With[{code =
  17. N[0.5 ThreeJSymbol[{l, 0}, {l, 0}, {L, 0}]^2 ThreeJSymbol[{l,
  18. m}, {l, m1}, {L, -(m + m1)}] ThreeJSymbol[{l, -p}, {l,
  19. p - (m + m1)}, {L, (m + m1)}] (2 l + 1)^2 (-1)^(
  20. p - (m + m1))]},
  21. Compile[{{l, _Integer}, {L, _Integer}, {m, _Integer}, {m1,
  22. _Integer}, {p, _Integer}}, code, CompilationTarget -> "C"]]];
  23.  
  24. Vex = Block[{l, L, m, m1, p},
  25. With[{code =
  26. N[0.5 ThreeJSymbol[{l, 0}, {l, 0}, {L,
  27. 0}]^2 ThreeJSymbol[{l, (-p)}, {l, m}, {l,
  28. p - m}] ThreeJSymbol[{l, m1 - m + p}, {l, -m1}, {L,
  29. m - p}] (2 l + 1)^2 (-1)^(m + m1)]},
  30. Compile[{{l, _Integer}, {L, _Integer}, {m, _Integer}, {m1,
  31. _Integer}, {p, _Integer}}, code, CompilationTarget -> "C"]]];
  32.  
  33. chmat = With[{Vcc = Vc, Vexx = Vex, kkdfxn = (If[# == #2, 1, 0] &)},
  34. Compile[{{nl, _Integer}, {nL, _Integer}, {nst, _Integer}, {ne,
  35. _Integer}, {mvec, _Real}, {Lvec, _Real}, {avec, _Real}, {svec, _Real}},
  36. Block[{ms, ms1, kf01, L0, tmp1, tmp2, m0, m10, p0, l0, tmp}, Table[
  37. ms = Compile`GetElement[svec, nms];
  38. ms1 = Compile`GetElement[svec, nms1];
  39. L0 = Compile`GetElement[Lvec, nL0];
  40. l0 = nl;
  41.  
  42. tmp = 0.0;
  43. Do[
  44. m0 = Compile`GetElement[mvec, nm];
  45. m10 = Compile`GetElement[mvec, nm1];
  46. p0 = Compile`GetElement[mvec, nm3];
  47.  
  48. tmp1 = Vcc[l0, L0, m0, m10, p0];
  49. tmp2 = Vexx[l0, L0, m0, m10, p0];
  50. kf01 = kkdfxn[ms, ms1];
  51.  
  52. Do[
  53. tmp += (tmp1 kf01) Compile`GetElement[avec, j,
  54. nm3] Compile`GetElement[avec, j, nm3 + nm1 - nm] -
  55. tmp2 Compile`GetElement[avec, j, nm3] Compile`GetElement[
  56. avec, j, nm1 - nm + nm3], {j, 1, ne}], {nm1, 1, nst}, {nm,
  57. 1, nst}, {nm3, 1, nst}];
  58.  
  59. , {nms, 1, nst}, {nms1, 1, nst}, {nL0, 1, nL}]]
  60.  
  61. , CompilationTarget -> "C",
  62. CompilationOptions -> {"InlineCompiledFunctions" -> True},
  63. RuntimeOptions -> "Speed"]];
  64.  
  65. chmat[nl, nL, nst, ne, mvec, Lvec, avec, svec]
Add Comment
Please, Sign In to add comment