(*-----------------------------------------------------------------------*) 
  
 (*-----------------------------------------------------------------------*) 

  
f[order_, QCC_, ωRF_,  Δt_,  n_]:=(
For [m = 1, m ≤ n, m++,
ωQ=(QCC/(2 (Sqrt[ 6])))(
           (Subscript[d, 1  C])Cos[(m - 
                                        1) Δt*ωrot]       +     (Subscript[d, 2  C])Cos[(m - 
                                        1)*2  Δt* ωrot]               
    + (Subscript[d, 1  S])Sin[(m - 
                                        1) Δt*ωrot]       +     (Subscript[d, 2  S])Sin[(m - 
                                        1)*2  Δt* ωrot]  ); 

ωQ21 = 0;     ωQ22 = 0 ;
If [order == 2,{
(Subscript[W, 20])=  (Subscript[d, 21  C])Cos[(m - 
                                        1) Δt*ωrot]       +    (Subscript[d, 22  C])Cos[(m - 
                                        
                                        1)*2 Δt* ωrot]               
    +(Subscript[d, 21  S])Sin[(m - 
                                        1) Δt*ωrot]       +    (Subscript[d, 22  S])Sin[(m - 
                                        1)*2 Δt* ωrot];
(Subscript[W, 40])= (Subscript[a, 40])
   +(Subscript[d, 44  S])Sin[(m - 
                                        1)*4  Δt*ωrot]  +   (Subscript[d, 44  C])Cos[(m - 
                                        1)*4  Δt*ωrot]
   +(Subscript[d, 43  S])Sin[(m - 
                                        1)*3  Δt*ωrot]  +   (Subscript[d, 43  C])Cos[(m - 
                                        1)*3  Δt*ωrot]
   +(Subscript[d, 42  S])Sin[(m - 
                                        1)*2  Δt*ωrot]  +   (Subscript[d, 42  C])Cos[(m - 
                                        1)*2  Δt*ωrot] 
   +(Subscript[d, 41  S])Sin[(m - 
                                        1)*   Δt*ωrot]  +   (Subscript[d, 41  C])Cos[(m - 
                                        1)*    Δt*ωrot];
ωQ21 = ((-1 )/(Subscript[ω, 0])) ((QCC^2)/36) (((-21)/(2 (Sqrt[ 70]))) 
                                      (Subscript[W, 40])   + (9/(2 (Sqrt[ 5]))) 
                                      (Subscript[W, 00])); 
ωQ22 = ((-1 )/(Subscript[ω, 0])) ((QCC^2)/36) ((27/(2 (Sqrt[ 70]))) 
                                      (Subscript[W, 40])  + (6/(Sqrt[14])) 
                                      (Subscript[W, 20])  - (3/(2 (Sqrt[ 5]))) (Subscript[W, 00]));
}];  
Ha=(({{ωQ + ωQ21, -((Sqrt[3])/2) ωRF, 0, 0}, {-((Sqrt[3])/2) ωRF, -ωQ + ωQ22, -ωRF, 0}, {0, -ωRF, -ωQ - ωQ22, -((Sqrt[3])/2) ωRF}, {0, 0, -((Sqrt[3])/2) ωRF, ωQ - ωQ21}}));    
{HT, Tp} = Eigensystem[N[Ha]];                                       T = 
                      Transpose[Tp];
(Subscript[n, 1]) = 
                      DiagonalMatrix[
                        Exp[-I*Δt*HT]] ;  
(Subscript[ρ, 1]) = T . (Subscript[n, 1]) . Tp;           (Subscript[ρ, 2]) = T . Conjugate[(Subscript[n, 1])] . Tp;              (Subscript[ρ, 0]) = (Subscript[ρ, 1]) . (Subscript[ρ, 0]) . (Subscript[ρ, 2]);        
s[m] = (Subscript[ρ, 0]);
] ;   
);     
  
powderMAS[order_, ω0Mhz_, QCCMHz_, ωRFkHz_, 
            VrotkHz_, tf_, tau_, η_, maxα_, maxβ_, 
            maxγ_, quanta_]:=(
  (Subscript[ω, 0])  = ω0Mhz*2  π*(10^3);           QCCbis = QCCMHz*2  π*(10^3);        ωbRFbis = ωRFkHz*2  π;       
  ωrot = VrotkHz*2  π;            Δt = tau*(10^(-3));                                 ns = tf/tau;  
  (Subscript[W, 00]) = ((Sqrt[5]) /10) (3 + (η^2));

  
For [i = 0,  
                i ≤ ns,  i++, h[i] = 0;];   
 
For [j = 0, j <  maxα, j++,Print[j];   
For [k = 1, k <  maxβ, k++,                       
                                                                               
For [c = 0, c <  maxγ, c++,                    
{  
(Subscript[ρ, 0]) = 
                                    DiagonalMatrix[{(3/2), 
                                        
                                        (1/2), -(1/2), -(3/2)}];

(Subscript[c, α]) =  
                                    Cos[((2 π j)/max)α];     (Subscript[s, α]) =  
                                    Sin[((2 π j)/max)α];      (Subscript[c, 2  α]) = 
                                    2  ((Subscript[c, α])^2) - 1;       (Subscript[s, 2  α]) = 
                                    2  (Subscript[c, α]) (Subscript[s, α]); 
(Subscript[c, β]) =  Cos[((π k)/max)β];     (Subscript[s, β]) =  
                                    Sin[((π k)/max)β];      (Subscript[c, 2  β]) = 
                                    2  ((Subscript[c, β])^2) - 1;       (Subscript[s, 2  β]) = 
                                    2  (Subscript[c, β]) (Subscript[s, β]);
(Subscript[c, γ]) =  
                                    Cos[((2 π  c)/max)γ];    (Subscript[s, γ]) =  
                                    Sin[((2 π  c)/max)γ];     (Subscript[c, 2  γ]) = 
                                    2  ((Subscript[c, γ])^2) - 1;      (Subscript[s, 2  γ]) = 
                                    2  (Subscript[c, γ]) (Subscript[s, γ]);
                
                                               
(Subscript[a, 1]) = -η (Subscript[s, 2  α]) 
                                      (Subscript[s, β])/(Sqrt[3]);                    (Subscript[b, 1]) = -(-3 + η (Subscript[c, 2  α])) (Subscript[s, 2  β])/(2 (Sqrt[ 3]));
(Subscript[a, 2]) = -η (Subscript[c, β]) (Subscript[s, 2  α])/(Sqrt[6]);                    (Subscript[b, 2]) = -(η (Subscript[c, 2  α]) (3 + (Subscript[c, 2  β])) + 
                                        6  ((Subscript[s, β])^2))/(4 (Sqrt[6]));
(Subscript[d, 2  S]) = (Subscript[a, 2]) (Subscript[c, 2  γ]) +  (Subscript[b, 2]) (Subscript[s, 2  γ]);                  (Subscript[d, 1  S]) = (Subscript[a, 1]) (Subscript[c, γ])  +  (Subscript[b, 1]) (Subscript[s, γ]);
(Subscript[d, 2  C]) = (Subscript[a, 2]) (Subscript[s, 2  γ])  -  (Subscript[b, 2]) (Subscript[c, 2  γ]);                 (Subscript[d, 1  C]) = 
                                    (Subscript[a, 1]) (Subscript[s, γ])  - (Subscript[b, 1]) 
                                        (Subscript[c, γ]); 

If [order == 2,{
 (Subscript[c, 4  α]) = 
                                        2  ((Subscript[c, 2  α])^2) - 1; 
 (Subscript[c, 4  β]) = 
                                        1 - 8  ((Subscript[c, β])^2) (1 - 
                                        ((Subscript[c, β])^2));         (Subscript[s, 4  β]) = 
                                        4  
                                        (Subscript[c, β]) (Subscript[c, 2  β]) 
                                        (Subscript[s, β]);
 (Subscript[c, 3  γ]) = (Subscript[c, γ]) (4  ((Subscript[c, γ])^2) - 3);                (Subscript[s, 3  γ]) = (Subscript[s, γ]) (4  ((Subscript[c, γ])^2) - 1);  
 (Subscript[c, 4  γ]) = 
                                        2  ((Subscript[c, 2  γ])^2) - 1;                        (Subscript[s, 4  γ]) = 4  (Subscript[c, γ]) (Subscript[c, 2  γ]) (Subscript[s, γ]);
    
(Subscript[a, 22]) = -(Sqrt[2/7])  η (Subscript[c, β]) (Subscript[s, 2  α]) ;          (Subscript[b, 22]) = -(η (Subscript[c, 2  α]) (3 + (Subscript[c, 2  β])) + 
                                        ((Subscript[s, β])^2) (-3 + (η^2)))/(2 (Sqrt[14]));
(Subscript[a, 21]) = -(2/(Sqrt[7])) η (Subscript[s, β]) (Subscript[s, 2  α]) ;      (Subscript[b, 21]) = (-3 - 
                                        2  (Subscript[c, 2  α]) η + (η^2)) (Subscript[s, 2  β])/(2 (Sqrt[7]));
(Subscript[d, 22  S]) = (Subscript[a, 22]) (Subscript[c, 2  γ]) +  (Subscript[b, 22]) (Subscript[s, 2  γ]);          (Subscript[d, 21  S]) = (Subscript[a, 21]) (Subscript[c, γ])  +  (Subscript[b, 21]) (Subscript[s, γ]);
(Subscript[d, 22  C]) = (Subscript[a, 22]) (Subscript[s, 2  γ])  -  (Subscript[b, 22]) (Subscript[c, 2  γ]);         (Subscript[d, 21  C]) = (Subscript[a, 21]) (Subscript[s, γ])  -  (Subscript[b, 21]) (Subscript[c, γ]);
    
(Subscript[a, 40]) = ((-(Sqrt[7/10]))/2304) ((18 + (η^2)) (9 + 
                                        20  (Subscript[c, 2  β]) + 
                                        35 (Subscript[c, 4  β])) + 
                                        240 η  (Subscript[c, 2  α]) (5 + 7 (Subscript[c, 2  β])) ((Subscript[s, β])^2) + 
                                        280 (η^2) (Subscript[c, 4  α])  ((Subscript[s, β])^4));
(Subscript[a, 41]) =  ((Sqrt[5/7])/72) η (Subscript[s, 2  α]) (Subscript[s, β]) (15 + 
                                        21  (Subscript[c, 2  β]) + 
                                        14 η (Subscript[c, 2  α]) 
                                        ((Subscript[s, β])^2));
(Subscript[b, 41]) = ((Sqrt[5/7]) /288) ((-18 - (η^2) - 
                                        12 η (Subscript[c, 2  α]) + 
                                        7 (η^2) 
                                        (Subscript[c, 4  α])) 
                                        (Subscript[s, 2  β]) - 
                                        7 (((-3 + η (Subscript[c, 2  α])))^2) (Subscript[s, 4  β]));
(Subscript[a, 42]) = - ((Sqrt[5/14]) /
                                        18) η (Subscript[c, β]) (Subscript[s, 2  α]) (-9 + 21  (Subscript[c, 2  β]) + 
                                        14 η (Subscript[c, 2  α]) 
                                        ((Subscript[s, β])^2));
(Subscript[b, 42]) = ((-1)/72) (Sqrt[5/14]) (3 η (Subscript[c, 2  α]) (5 + 
                                        4 (Subscript[c, 2  β]) + 
                                        7 (Subscript[c, 4  β])) + (7 (η^2) (Subscript[c, 4  α]) (3 + 
                                        (Subscript[c, 2  β])) + (18 + (η^2)) (5 + 7  (Subscript[c, 2  β]))) ((Subscript[s, β])^2));
(Subscript[a, 43]) = -((Sqrt[35])/72) η (-3 - 9 
                                         (Subscript[c, 2  β]) + η (Subscript[c, 2  α]) (5 + 3 (Subscript[c, 2  β]))) (Subscript[s, 2  α]) 
                                        (Subscript[s, β]);
(Subscript[b, 43]) = -((Sqrt[35])/288) (-18 - (η^2) - 
                                        12 η (Subscript[c, 2  α]) + 
                                        7 (η^2) (Subscript[c, 4  α]) + 
                                        2 (((-3 + η (Subscript[c, 2  α])))^2) (Subscript[c, 2  β])) (Subscript[s, 2  β]);
(Subscript[a, 44]) = -((Sqrt[35/2])/72) η (Subscript[c, β]) 
                                        (Subscript[s, 2  α]) (η (Subscript[c, 2  α]) (3 + (Subscript[c, 2  β])) + 6 ((Subscript[s, β])^2));
(Subscript[b, 44]) = -((Sqrt[35/2])/2304) ((η^2) (Subscript[c, 4  α]) (35 + 
                                        28 (Subscript[c, 2  β]) + 
                                        (Subscript[c, 4  β])) + 
                                        48 η (Subscript[c, 2  α]) (3 + (Subscript[c, 2  β])) ((Subscript[s, β])^2) + 8 (18 + (η^2)) 
                                        ((Subscript[s, β])^4));
(Subscript[d, 41  S]) = (Subscript[a, 41]) (Subscript[c, γ])   + (Subscript[b, 41]) (Subscript[s, γ]);       (Subscript[d, 41  C]) = (Subscript[a, 41]) 
                                        (Subscript[s, γ])  - (Subscript[b, 41]) 
                                        (Subscript[c, γ]);      (Subscript[d, 42  S]) = (Subscript[a, 42]) 
                                        (Subscript[c, 2  γ]) + (Subscript[b, 42]) 
                                        (Subscript[s, 2  γ]);    
(Subscript[d, 42  C]) = (Subscript[a, 42]) (Subscript[s, 2  γ]) - (Subscript[b, 42]) (Subscript[c, 2  γ]) ;   (Subscript[d, 43  S]) = (Subscript[a, 43]) 
                                        (Subscript[c, 3  γ]) + (Subscript[b, 43]) 
                                        (Subscript[s, 3  γ]);   (Subscript[d, 43  C]) = (Subscript[a, 43]) 
                                        (Subscript[s, 3  γ]) - (Subscript[b, 43]) 
                                        (Subscript[c, 3  γ]);
(Subscript[d, 44  S]) = (Subscript[a, 44]) (Subscript[c, 4  γ]) + (Subscript[b, 44]) (Subscript[s, 4  γ]);    (Subscript[d, 44  C]) = (Subscript[a, 44]) 
                                        (Subscript[s, 4  γ]) - (Subscript[b, 44]) 
                                        (Subscript[c, 4  γ]);
}];  

f[order, 
                                    QCCbis, ωbRFbis, Δt, 
                                    ns];    

For [i = 1, i ≤ ns, i++,
  
If [quanta == 1,{s[i] = 
                                        Im[s[i][[3, 2]]], h[i]=h[i]+(Subscript[s, β])s[i]*2/5}];
    
 If [quanta == 3,{s[i] = 
                                        Im[s[i][[4, 1]]], h[i]=h[i]+(Subscript[s, β])s[i]}];
];    

}];    
 ];       
];           

For [i = 1,  
                i ≤ 
                  ns,  i++,  h[i] = 
                    h[i]/(maxα (maxβ - 1) 
                          maxγ); ];   

  
  
Print["**********************************"];
For[a = 0,  
                a ≤ ns,  a++, time[a] = a*tau;];
powderMAS = 
                Chop[Table[{tt, time[tt], NumberForm[h[tt], 10]}, {tt, 0, 
                      ns}]]; 
\(Print[
                TableForm[
                  powderMAS, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ TableHeadings -> {None, {"Rang", "t(μs)", "intensity"}}]]\);
  
Print["**********************************"];
\(ListPlot[
                Table[{\ tt*tau, h[tt]}, {tt, 0, 
                    ns}], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotJoined -> 
                  True, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotLabel -> "Int=f(t)", 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ AxesLabel -> {"t(μs)", "Int.(U.A.)"}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ PlotStyle -> {Hue[0.1]}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ TextStyle -> {FontFamily -> "Times", FontSize -> 12}]\);
);

powderMAS[    2        , 
            105.8731007,     8           ,        100    ,       15         , 20 ,  
            1.0    ,  -1 ,      2   ,     3     ,   3     ,     1         ];

(*-----------------------------------------------------------------------------*) 
 (*                   Table powderMAS . 
              m in Microsoft EXCEL format              *) 
 (*-----------------------------------------------------------------------------*) 
Clear[writeExcel];
writeExcel[filename_String,  data_List] :=
     Module[  { file  =  OpenWrite[filename] },
                      Scan[( 
                               WriteString[file,  First[#]];
                                \(Scan[
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ WriteString[
                                  file, \ "	", \ #] &, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Rest[#]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]\);    
                                WriteString[
                              file,  "\n"]
                           )&,
                      data
                     ];    
                       Close[file]
    ]  ;

writeExcel["powderMAS.m",  powderMAS];

Remove[order, ω0Mhz, QCCMHz, ωRFkHz, VrotkHz, tf, 
      tau, η, maxα, maxβ, maxγ, 
      QCCbis,  ωbRFbis, 
      ns, 
             α, β, γ, i, j, k, c, h, ρ0, f, a, powderMAS, n, s, m, ωQ, Ha, T, 
      Tp, HT, n1, ωQ21, ωQ22]