<encoding="UTF-8">
磁気形状因子の計算（双極子近似ではなく6次の項までちゃんと考える）

references:
[1] E. Balcar and S. W. Lovesey, "Theory of Magnetic Neutron and Photon Scattering" (Oxford).
[2] S. W. Lovesey, "Theory of Neutron Scattering from Condensed Matter" vol. 2, (Oxford).
[3] International Tables for Crystallography, Vol. C, ed. A. J. C. Wilson and E. Prince, (1999).


(1) 
球ベッセル関数J_kの展開係数を入力し，<j_0(kappa)>, <j_2(kappa)>, <j_4(kappa)>, <j_6(kappa)>を定義．
(* Coefficients for Ce2+ by P. J. Brown. [3] *)
rjcoef[0] = {0.2953, 17.6846, 0.2923, 6.7329, 0.4313, 5.3827, -0.0194};
rjcoef[2] = {0.9809, 18.063, 1.8413, 7.7688, 0.9905, 2.8452, 0.012};
rjcoef[4] = {-0.6468, 10.5331, 0.4052, 5.6243, 0.3412, 1.5346, 0.008};
rjcoef[6] = {-0.1212, 7.994, -0.0639, 4.0244, 0.1519, 1.0957, 0.0078};

frj[l_, s_] := ((Sum[ rjcoef[l][[i]]*Exp[-rjcoef[l][[i + 1]]*s^2], {i, 1, 5, 2} ] + rjcoef[l][[7]])*If[l == 0, 1, s^2]);

Plot[{frj[0, s], frj[2, s], frj[4, s], frj[6, s]}, {s, 0, 1.2}, PlotRange -> All]


(2) いろいろな定義．
(* 2.1 Definition (3.4.9) and (3.4.10) *)

AA[K_, Kd_, l_] := ((-1)^(l + 1)*Sqrt[(2*l + 3)*(l + 1)/(2*l + 1)]
               *ThreeJSymbol[{l, 0}, {K, 0}, {l + 1, 0}]*SixJSymbol[{1, Kd, K}, {l, l + 1, l}])
BB[K_, Kd_, l_] := ((-1)^l*Sqrt[(2*l - 1)*l/(2*l + 1)]
               *ThreeJSymbol[{l, 0}, {K, 0}, {l - 1, 0}]*SixJSymbol[{1, Kd, K}, {l, l - 1, l}])

(* 2.2 Reduced matrix element (theta||W (0,K')||theta) and (theta||W (1,K)||theta). Table 3.1 and 3.2. *)
rmW0 = {
    {Sqrt[3], 0, Sqrt[7], 0, Sqrt[11], 0, 0},
    {3/2*Sqrt[55/7], 0, 0, 0, -Sqrt[429/14], 0, 0},
    {Sqrt[39], 0, -2*Sqrt[91/33], 0, Sqrt[442/33], 0, 0},
    {1/2*Sqrt[195], 0, -Sqrt[455/33], 0, Sqrt[1105/66], 0, 0},
    {3*Sqrt[55/14], 0, 0, 0, -Sqrt[429/7], 0, 0},
    {Sqrt[21/2], 0, 7*Sqrt[1/2], 0, Sqrt[77/2], 0, 0},
    {0, 0, 0, 0, 0, 0, 0},
    {Sqrt[21/2], 0, 7*Sqrt[1/2], 0, Sqrt[77/2], 0, 0},
    {3*Sqrt[55/14], 0, 0, 0, -Sqrt[429/7], 0, 0},
    {1/2*Sqrt[195], 0, -Sqrt[455/33], 0, Sqrt[1105/66], 0, 0},
    {Sqrt[39], 0, -2*Sqrt[91/33], 0, Sqrt[442/33], 0, 0},
    {3/2*Sqrt[55/7], 0, 0, 0, -Sqrt[429/14], 0, 0},
    {Sqrt[3], 0, Sqrt[7], 0, Sqrt[11], 0, 0}};

rmW1 = {
    {Sqrt[3], 0, Sqrt[15], 0, 3*Sqrt[3], 0, Sqrt[39], 0, 0},
    {2*Sqrt[33/7], 0, Sqrt[715/42], 0, -2*Sqrt[39/7], 0, -Sqrt[1105/21], 0, 0},
    {Sqrt[390/7], 0, 5/3*Sqrt[13/11], 0, -2/11*Sqrt[1105/3], 0, 5/33*Sqrt[41990/7], 0, 0},
    {2*Sqrt[195/7], 0, -5/2*Sqrt[13/22], 0, 1/11*Sqrt[3315/2], 0, -5/22*Sqrt[20995/7], 0, 0},
    {Sqrt[165], 0, -Sqrt[143/6], 0, 2*Sqrt[39/5], 0, Sqrt[221/3], 0, 0},
    {2*Sqrt[42], 0, -Sqrt[70/3], 0, -Sqrt[42], 0, -Sqrt[182/3], 0, 0},
    {6, 0, 0, 0, 0, 0, 0, 0},
    {2*Sqrt[42], 0, -Sqrt[70/3], 0, -Sqrt[42], 0, -Sqrt[182/3], 0, 0},
    {Sqrt[165], 0, -Sqrt[143/6], 0, 2*Sqrt[39/5], 0, Sqrt[221/3], 0, 0},
    {2*Sqrt[195/7], 0, -5/2*Sqrt[13/22], 0, 1/11*Sqrt[3315/2], 0, -5/22*Sqrt[20995/7], 0, 0},
    {Sqrt[390/7], 0, 5/3*Sqrt[13/11], 0, -2/11*Sqrt[1105/3], 0, 5/33*Sqrt[41990/7], 0, 0},
    {2*Sqrt[33/7], 0, Sqrt[715/42], 0, -2*Sqrt[39/7], 0, -Sqrt[1105/21], 0, 0},
    {Sqrt[3], 0, Sqrt[15], 0, 3*Sqrt[3], 0, Sqrt[39], 0, 0}};

(* 2.3 Definition (3.6.8) for S'=S, L'=L. *)
Ao[K_, Kd_, l_, n4f_, L_, S_, J_, Jd_] := N[I^(Kd - 1)*(2*l + 1)^2*(-1)^(Jd + S + L)
     *Sqrt[2*(2*Jd + 1)/3/(2*S + 1)/(2*Kd + 1)]
     *If[Abs[J - Jd] <= Kd <= (J + Jd) && Kd <= 2*L,SixJSymbol[{Jd, Kd, J}, {L, S, L}], 0]
     *If[Kd <= (2*l - 1), AA[Kd, Kd, l], 0]*rmW0[[n4f, Kd]]
     *Switch[K, Kd + 1, Sqrt[Kd], Kd - 1, Sqrt[Kd + 1]]]*(rj[Kd - 1] + rj[Kd + 1])

(* Sum of the c.f.p in the last line of (3.6.8), deduced by comparison with (3.6.9). *)
cfpsum = {
      {-0.142857, 0, -0.142857, 0, -0.142857, 0, 0},
      {-0.180187, 0, 0., 0, 0.123888, 0, 0},
      {0.196116, 0, -0.0682789, 0, 0.0600205, 0, 0},
      {0.196116, 0, -0.0682789, 0, 0.0600205, 0, 0},
      {-0.180187, 0, 0., 0, 0.123888, 0, 0},
      {0, 0, 0, 0, 0, 0, 0},
      {0, 0, 0, 0, 0, 0, 0},
      {-0.142857, 0, -0.142857, 0, -0.142857, 0, 0},
      {-0.180187, 0, 0., 0, 0.123888, 0, 0},
      {0.196116, 0, -0.0682789, 0, 0.0600205, 0, 0},
      {0.196116, 0, -0.0682789, 0, 0.0600205, 0, 0},
      {-0.180187, 0, 0., 0, 0.123888, 0, 0},
      {-0.142857, 0, -0.142857, 0, -0.142857, 0, 0}};

(*  Definition of E (K,K') (3.6.4) for S'=S, L'=L. *)
Eo[K_, Kd_, l_, n4f_, L_, S_, J_, Jd_] := N[I^(K)/2/Sqrt[3]*Sqrt[l*(l + 1)]
        *Sqrt[(2*l + 1)*(2*K + 1)*(2*Jd + 1)]*(2*l + 1)*(2*Kd + 1)*(2*L + 1)
        *(-1)^(S + L + L + Jd)*ThreeJSymbol[{l, 0}, {K, 0}, {l, 0}]
        *SixJSymbol[{1, Kd, K}, {l, l, l}]*SixJSymbol[{L, Kd, L}, {Jd, S, J}]
        *cfpsum[[n4f, Kd]]]*(rg[K])

(* 2.4 Definition of 9j-Symbol for S'=S, L'=L. Appendix D .3. 
  Nonvanishing only when |K-K'|<=1<=K+K', 0<=K<=2L, and |J-J'|<=K'<=J+J'. *)
nineJSymbol[K_, Kd_, S_, L_, J_, Jd_] := Sum[
     If[Abs[L - S] <= X <= (L + S) && Abs[Jd - K] <= X <= (Jd + K),
      (-1)^(2*X)*(2*X + 1)*SixJSymbol[{1, S, S}, {L, J, X}]*SixJSymbol[{K, L, L}, {S, X, Jd}]
      *SixJSymbol[{Kd, Jd, J}, {X, 1, K}] , 0], {X, Abs[J - 1], J + 1}]

(* 2.5 Definition (3.6.12) for S'=S, L'=L. *)
Cs[K_, Kd_, l_, n4f_, L_, S_, J_, Jd_] := N[I^(K + 2)*Sqrt[1/6]
        *Sqrt[(2*l + 1)*(2*l + 1)*(2*Jd + 1)*(2*Kd + 1)*(2*Kd + 1)]
        *(-1)^(L + L + l)*If[K <= 2*l, ThreeJSymbol[{l, 0}, {K, 0}, {l, 0}], 0]
        *If[Abs[K - Kd] <= 1 <= (K + Kd) && 0 <= K <= 2*L && Abs[J - Jd] <= Kd <= (J + Jd),
        nineJSymbol[K, Kd, S, L, J, Jd], 0]*rmW1[[n4f, K + 1]]]*rj[K]

(* 2.6 Definition (4.2.10), (4.2.11).  *)
Bs[K_, Kd_, l_, n4f_, L_, S_, J_, Jd_] := Switch[Kd, K,
      If[Max[2, Abs[J - Jd]] <= K <= Min[2*l, 2*L, J + Jd], Cs[K, K, l, n4f, L, S, J, Jd], 0]
    , K + 1, (K + 2)/(2*K + 3)*
    (If[Max[0, Abs[J - Jd] - 1] <= K <= Min[2*l, 2*L, J + Jd - 1], Cs[K, K + 1, l, n4f, L, S, J, Jd], 0]
   + If[Max[0, Abs[J - Jd] - 1] <= K <= Min[2*l - 2, 2*L - 2, J + Jd - 1],
      N[Sqrt[(K + 1)/(K + 2)]]*Cs[K + 2, K + 1, l, n4f, L, S, J, Jd], 0])
    , K - 1, (K - 1)/(2*K - 1)*
    (If[Max[2, Abs[J - Jd] + 1] <= K <= Min[2*l, 2*L, J + Jd + 1], Cs[K, K - 1, l, n4f, L, S, J, Jd], 0]
   + If[Max[2, Abs[J - Jd] + 1] <= K <= Min[2*l + 2, 2*L + 2, J + Jd + 1],
      N[Sqrt[K/(K - 1)]]*Cs[K - 2, K - 1, l, n4f, L, S, J, Jd], 0])   ]


(3) 自由イオンにおける多重項内および多重項間遷移
基底J多重項（始状態）は既にJとして入力されているので，終状態の多重項Jdを入力．
(* 3.1 structure facotr G (q,J,J') for elastic intramultiplet and inelastic inter-multiplet transition: eq. (5.2.20).  Appendix E.3. Input Jd*)
Jd = 5/2;

Gqj = Expand[Sum[3/(Kd + 1)*(Ao[Kd - 1, Kd, l4f, n4f, L, S, J, Jd]
    + Bs[Kd - 1, Kd, l4f, n4f, L, S, J, Jd])^2, {Kd, 1, 2*l4f + 1, 2}]
   + Sum[3/(2*K + 1)*Bs[K, K, l4f, n4f, L, S, J, Jd]^2, {K, 2, 2*l4f, 2}]]
coef = {
  Coefficient[Gqj, rj[0]^2], Coefficient[Gqj, rj[0]*rj[2]],
  Coefficient[Gqj, rj[2]^2], Coefficient[Gqj, rj[2]*rj[4]],
  Coefficient[Gqj, rj[4]^2], Coefficient[Gqj, rj[4]*rj[6]],
  Coefficient[Gqj, rj[6]^2]}

(* 3.2 For elastic scatterings, J'=J,  the coefficient of rj[0]^2 is given by g^2J (J+1)/6. *)
g = 1 + (J*(J + 1) + S*(S + 1) - L*(L + 1))/(2*J*(J + 1));
Print[N[g^2*J*(J + 1)/6]]
fac = Table[{s,
     (coef[[1]]*frj[0, s]^2 + coef[[2]]*frj[0, s]*frj[2, s]
       + coef[[3]]*frj[2, s]^2 + coef[[4]]*frj[2, s]*frj[4, s]
       + coef[[5]]*frj[4, s]^2 + coef[[6]]*frj[4, s]*frj[6, s]
       + coef[[7]]*frj[6, s]^2) }  , {s, 0, 1.2, 0.02}];
ffree = Interpolation[fac]

fig1=Plot[ffree[x],{x,0,1.2}]

(* 3.3 Dipole approximation for the intramultiplet transition *)

facdip = Table[{s,
     g^2*J*(J + 1)/6*(frj[0, s]^2 + 2*(2/g - 1)*frj[0, s]*frj[2, s]
        + (2/g - 1)^2*frj[2, s]^2) } , {s, 0, 1.2, 0.02}];
fdip = Interpolation[facdip]
Plot[{ffree[x], fdip[x]}, {x, 0, 1.2}, PlotStyle -> {RGBColor[0, 0, 0], RGBColor[1, 0, 0]}]


(4) 基底J多重項での磁気形状因子の計算．結晶場や分子場での固有関数は他で計算されて決まっているとする．

(* 4.1 Table of A (K,K') and C (K,K'). The outputs represent : {(K, K'), {A (j0), A (j2), A (j4), A (j6)}, {C (j0), C (j2), C (j4), C (j6)}} *)
ACtable = Flatten[Table[ {{K, (K + jj)},
      {f = If[1 <= (K + jj) <= (2*l4f + 1), Expand[Ao[K, K + jj, l4f, n4f, L, S, J, J]], 0];
       Coefficient[f, rj[0]], Coefficient[f, rj[2]], Coefficient[f, rj[4]], Coefficient[f, rj[6]]},
      {f = If[1 <= (K + jj) <= (2*l4f + 1), Expand[Cs[K, K + jj, l4f, n4f, L, S, J, J]], 0];
       Coefficient[f, rj[0]], Coefficient[f, rj[2]],Coefficient[f, rj[4]], Coefficient[f, rj[6]]}}
    , {K, 0, 2*l4f + 2, 2}, {jj, -1, 1}], 1]

(* 4.2 Table of E (K,K') and C (K,K'). The outputs represent : {(K, K'), {E (g0), E (g2), E (g4), E (g6)}, {C (j0), C (j2), C (j4), C (j6)}} *)
ECtable = Flatten[Table[ {{K, (K + jj)},
      {f = If[1 <= (K + jj) <= (2*l4f + 1), Expand[Eo[K, K + jj, l4f, n4f, L, S, J, J]], 0];
       Coefficient[f, rg[0]], Coefficient[f, rg[2]],Coefficient[f, rg[4]], Coefficient[f, rg[6]]},
      {f = If[1 <= (K + jj) <= (2*l4f + 1), Expand[Cs[K, K + jj, l4f, n4f, L, S, J, J]], 0];
       Coefficient[f, rj[0]], Coefficient[f, rj[2]],Coefficient[f, rj[4]], Coefficient[f, rj[6]]}}
    , {K, 0, 2*l4f + 2, 2}, {jj, -1, 1}], 1]

Do[
 K = ACtable[[i, 1, 1]]; If[ACtable[[i, 1, 2]] == ACtable[[i, 1, 1]] + 1,
  Print[ACtable[[i, 1]], Tab, 2/(K + 2)*ECtable[[i, 2]], Tab, (2*K + 3)/(K + 2)*ACtable[[i, 2]]]
  ], {i, 1, 15}]


(5) Qベクトル一般形の計算．
既に固有値と固有関数は決まっていて，{valsa, vecsa}であるとする．

zA = Sum[Exp[-valsa[[i]]/temp], {i, 1, 2*J + 1}];
tempfacA = Chop[Table[Exp[-valsa[[i]]/temp]/zA, {i, 1, 2*J + 1}]]

(* 5.1 初期状態iina~iinb，終状態ifna~ifnbという計算の空間を指定．無視してよい準位は無視する． *)
iina = 1; iinb = 4;  
ifna = 1; ifnb = 4;

(* 5.2 Make table of the matrix elements of Q between initial states and final states. eq. (3.6.10) and (3.6.14).*)

Clear[rj0, rj2, rj4, rj6, th, ph];
Module[{M, Md, K, Kd},
 Do[
   QelspinA[rK] = Table[Table[  Chop[ComplexExpand[Sum[M = J - j + 1; Md = J - k + 1;
     If[Conjugate[vecsa[[ii, j]]]*vecsa[[if, k]] == 0, {0, 0, 0, 0},
        Conjugate[vecsa[[ii, j]]]*vecsa[[if, k]]*Sum[ K = ACtable[[i, 1, 1]]; Kd = ACtable[[i, 1, 2]];
     If[1 <= Kd <= Min[(2*l4f + 1), (J + J)] && Abs[q - M + Md] <= K && Abs[M - Md] <= Kd,
         SphericalHarmonicY[K, q - M + Md, th, ph]*N[Sqrt[4*Pi]
          *(ACtable[[i, 3]])*ClebschGordan[{K, q - M + Md}, {Kd, M - Md}, {1, q}]
          *ClebschGordan[{Kd, M - Md}, {J, Md}, {J, M}] ]
          , {0, 0, 0, 0}]   , {i, rK/2*3 + 1, rK/2*3 + 3}] ]
       , {j, 1, 2*J + 1}, {k, 1, 2*J + 1}] ]]  , {q, -1, 1}], {ii,iina, iinb}, {if, ifna, ifnb}] ;
   QelorbtA[rK] = Table[Table[  Chop[ComplexExpand[Sum[M = J - j + 1; Md = J - k + 1;
     If[Conjugate[vecsa[[ii, j]]]*vecsa[[if, k]] == 0, {0, 0, 0, 0},
        Conjugate[vecsa[[ii, j]]]*vecsa[[if, k]]*Sum[ K = ACtable[[i, 1, 1]]; Kd = ACtable[[i, 1, 2]];
     If[1 <= Kd <= Min[(2*l4f + 1), (J + J)] && Abs[q - M + Md] <= K && Abs[M - Md] <= Kd && Kd == K + 1,
         SphericalHarmonicY[K, q - M + Md, th, ph]*N[Sqrt[4*Pi]*(2*K + 3)/(K + 2)
          *(ACtable[[i, 2]])*ClebschGordan[{K, q - M + Md}, {Kd, M - Md}, {1, q}]
          *ClebschGordan[{Kd, M - Md}, {J, Md}, {J, M}] ]
          , {0, 0, 0, 0}]      , {i, rK/2*3 + 1, rK/2*3 + 3}] ]
       , {j, 1, 2*J + 1}, {k, 1, 2*J + 1}] ]]  , {q, -1, 1}], {ii,iina, iinb}, {if, ifna, ifnb}]
    , {rK, 0, 6, 2}];
  ]

(* 5.3  Transform Q from spherical to Cartesian components. *)
Do[
  QxyzspinA[rK] = Table[Chop[{(QelspinA[rK][[ii, if, 1]] - QelspinA[rK][[ii, if, 3]])/Sqrt[2],
          (QelspinA[rK][[ii, if, 1]] + QelspinA[rK][[ii, if, 3]])*I/Sqrt[2],
           QelspinA[rK][[ii, if, 2]]}], {ii, 1, iinb - iina + 1}, {if, 1, ifnb - ifna + 1}];
  QxyzorbtA[rK] = Table[Chop[{(QelorbtA[rK][[ii, if, 1]] - QelorbtA[rK][[ii, if, 3]])/Sqrt[2],
          (QelorbtA[rK][[ii, if, 1]] + QelorbtA[rK][[ii, if, 3]])*I/Sqrt[2],
           QelorbtA[rK][[ii, if, 2]]}], {ii, 1, iinb - iina + 1}, {if, 1, ifnb - ifna + 1}];
   , {rK, 0, 6, 2}];


(* 5.4 Function of Qxyz table between initial and final states. *)
Do[
  QspinA[rK][th_, ph_, s_] := Evaluate[Table[Table[Chop[ ComplexExpand[
       Sum[QxyzspinA[rK][[ii, ii, alp, jK/2 + 1]]*frj[jK, s], {jK, 0, 6, 2}] 
        ] ], {alp, 1, 3}] , {ii, 1, iinb - iina + 1}]];
  QorbtA[rK][th_, ph_, s_] := Evaluate[Table[Table[Chop[ ComplexExpand[
       Sum[QxyzorbtA[rK][[ii, ii, alp, jK/2 + 1]]*frj[jK, s], {jK, 0, 6, 2}] 
        ] ], {alp, 1, 3}] , {ii, 1, iinb - iina + 1}]];
    , {rK, 0, 6, 2}]


(6) Bragg反射

(* scattering vector in 1/Angstrom *)
qvec={1,1,1}  

(* 6.1 Magnetic Form-Factor *)
qx = qvec[[1]]; qy = qvec[[2]]; qz = qvec[[3]];
tau = Sqrt[qvec[[1]]^2 + qvec[[2]]^2 + qvec[[3]]^2];
kapp = qvec/tau;  qval = tau/4/Pi;
Which[qx == 0 && qy > 0, phi = Pi/2, qx == 0 && qy < 0, phi = 3*Pi/2, qy == 0 && qx > 0, phi = 0, 
  qy == 0 && qx < 0, phi = Pi, qx == 0 && qy == 0, phi = 0, True, phi = N[ArcTan[qx, qy]]];
Which[qz > 0 && qx == 0 && qy == 0, theta = 0, qz < 0 && qx == 0 && qy == 0, theta = Pi, 
  True, theta = N[ArcCos[qz/tau]]];

ffspinA = Chop[Table[Sum[QspinA[rK][theta, phi, qval][[ii]], {rK, 0, 6, 2}], {ii, 1, iinb - iina + 1}] ];
fforbtA = Chop[Table[Sum[QorbtA[rK][theta, phi, qval][[ii]], {rK, 0, 6, 2}], {ii, 1, iinb - iina + 1}] ];
ffA = ffspinA + fforbtA;
qperpA = Table[Cross[kapp, Cross[ffA[[ii]], kapp]], {ii, 1, iinb - iina + 1}];

(* 6.2 Bragg reflection of Unpolarized Neutrons *)
qperp = Sum[qperpA[[j]]*tempfac[[iina + j - 1]], {j, 1, iinb - iina + 1}] ;
sqwBragg=Chop[ Conjugate[qperp].qperp];
Print["S(q)_Bragg=",sqwBragg]

(* 6.3 Bragg reflection of Polarized Neutrons *)
(* Definition of the Pauli matrix. *)
pauli = {{{0, 1}, {1, 0}}, {{0, -I}, {I, 0}}, {{1, 0}, {0, -1}}};
pauli2 = Table[(pauli[[i]].pauli[[j]]), {i, 1, 3}, {j, 1, 3}];

(* Input directions (theta,phi) of incident (1) and scattered (2) neutrons. check1 and 2 must be zero. *)
th1 = Pi/2; ph1 = -Pi/4;  th2 = 0; ph2 = 0;

nsdir1 = N[{Sin[th1] Cos[ph1], Sin[th1] Sin[ph1], Cos[th1]}];
nsdir2 = N[{Sin[th2] Cos[ph2], Sin[th2] Sin[ph2], Cos[th2]}];
Print["direction of incident neutron spin = ", nsdir1]
Print["direction of scattered neutron spin = ", nsdir2]
(* Calculate spin functions of incident and scattered neutrons. *)
sp1 = N[{Cos[th1/2]*Exp[-I*ph1], Sin[th1/2]}];
sp2 = N[{Cos[th2/2]*Exp[-I*ph2], Sin[th2/2]}];
sp1nsf = sp1;
sp1sf = N[{-Sin[th1/2]*Exp[-I*ph1], Cos[th1/2]}];
check1 = Chop[(Table[Conjugate[sp1].pauli[[i]].sp1, {i, 1, 3}] - nsdir1)];
check2 = Chop[(Table[Conjugate[sp2].pauli[[i]].sp2, {i, 1, 3}] - nsdir2)];

nsig = Table[Chop[Conjugate[sp2].pauli[[i]].sp1], {i, 1, 3}];
nsigsf = Table[Chop[Conjugate[sp1sf].pauli[[i]].sp1], {i, 1, 3}];
nsignsf = Table[Chop[Conjugate[sp1nsf].pauli[[i]].sp1], {i, 1, 3}];
nsig2 = Chop[Table[Conjugate[sp1].pauli2[[i, j]].sp1, {i, 1, 3}, {j, 1, 3}]];

Print["incident neutron spin state: chi=", Chop[sp1]]
Print["scattered neutron spin state: chi'=", Chop[sp2]]
Print["check1=", check1, Tab, "check2=", check2]

(* Sq_Bragg *)
sqwNSF = Abs[nsignsf.qperp]^2; (* Non-Spin-Flip (chi->chi) *)
sqwSF = Abs[nsigsf.qperp]^2;  (* Spin-Flip (chi->-chi) *)
sqw = Abs[nsig.qperp]^2;  (*  chi->chi' *)
sqwNPA = Sum[nsig2[[i, j]]*Conjugate[qperp[[i]]]*qperp[[j]], {i, 1, 3}, {j, 1, 3}]; (* without polarization analysis (chi->  )  *) 
Print["S(q, chi-NSF)_Bragg=",sqwNSF]
Print["S(q, chi-SF)_Bragg=",sqwSF]
Print["S(q, chi-chi')_Bragg=",sqw]
Print["S(q, chi- )_Bragg=",sqwNPA]


A格子とB格子からなる系で構造因子がAとBの差で表される超格子反射の場合：
(5)以下をB格子{valsb, vecsb}についても行い，
zB, tempfacB, QelspinB, QelorbtB, QxyzspinB, QxyzorbtB, QspinB, QorbtB, ffspinB, fforbtB, ffB, qperpB
を計算する．次に(6.2)以下で
qperp = (Sum[qperpA[[j]]*tempfac[[iina + j - 1]], {j, 1, iinb - iina + 1}] 
            - Sum[qperpB[[j]]*tempfac[[iina + j - 1]], {j, 1, iinb - iina + 1}]);
として計算すればよい．