run:=Module[{QCC,  q, amplitude},
(Subscript[ω, 0]) = 2  π*larmorFrequencyMhz;
dim = 
                        IntegerPart[2  quadrupoleSpin + 1];
ms = 
                        Range[quadrupoleSpin, -quadrupoleSpin, -1]; 
Iz = DiagonalMatrix[ms];    
ωrot = 
                        2  π*spinRatekHz*(10^(-3));
QCC = 2  π*QCCMHz;      
q = 
                        2  quadrupoleSpin (2  quadrupoleSpin - 1);
A = ((Sqrt[6])/(2  q)) QCC*
                          DiagonalMatrix[(3  (ms^2) - 
                                  quadrupoleSpin (quadrupoleSpin + 1))/
                              3];
If[
                        quadrupoleOrder == 
                          2, {
  W0 = (Sqrt[5]) (3 + (η^2))/
                                10; 
  amplitude = ((-1 )/(Subscript[ω, 0])) 
                              ((QCC^2)/(q^2)); 
  A4 = (amplitude/(2 (Sqrt[ 70]))) (18  quadrupoleSpin (quadrupoleSpin + 1) - 
                                  34  (ms^2) - 5) 
                              ms; 
  A2 = (amplitude/(2 (Sqrt[ 14]))) (8  quadrupoleSpin (quadrupoleSpin + 1) - 
                                  12  (ms^2) - 3) 
                              ms; 
  A0 = (amplitude/(Sqrt[5])) ( quadrupoleSpin (quadrupoleSpin + 1) -   3 
                                     (ms^2)) ms;
}];
If[
                        numberOfGammaAngles == 1,  
                        crystal[(Subscript[α, PR]), (Subscript[β, PR]),  (Subscript[γ, PR])], 
    powder[powderFile,  numberOfGammaAngles]];
];
  
  
radiofreq[ωRF_]:=Module[{RF,  ymax},
HRF = Table[0, {dim}, {dim}];
RF = - π*ωRF*(10^(-3)) (Sqrt[quadrupoleSpin (quadrupoleSpin + 1) - ms (ms - 1)]);
ymax = 
                        IntegerPart[quadrupoleSpin + 0.5];
For [y = 1,  
                        y ≤ 
                          ymax,  y++, 
    HRF[[y, y + 1]] = HRF[[y + 1, 
                                  y]] = HRF[[dim - y, 
                                    dim - y + 1]] = HRF[[dim - y + 1, 
                                      dim - y]] = 
                                  RF[[y]];
]; 
];
  
  
pulse[
                  tau_, ωRF_] := (
radiofreq[ωRF]; 
nbPulseIncrement = 
                    IntegerPart[tau/Δt]; 
If[
                    nbPulseIncrement == 1,  pulsevar, 
                    pulsefixe];
);
  
  
  
pulseFunction[poub_]:=
   Module[{V20,W20,W40,frequency,HQ1,HQ2,Ha,HT,Tp,T,mat,ρ11,ρ22},
 frequency = (poub - 
                              1)*Δt*ωrot;
 V20=  d1C*Cos[frequency]    +  d2C*Cos[2  frequency]
     +d1S*Sin[frequency]    +  d2S*Sin[2  frequency];

 HQ1 = V20*A;       HQ2 = 0;

 If[quadrupoleOrder == 2,{
  W20=d21C*Cos[frequency] + d22C*Cos[2  frequency]
     +d21S*
                                    Sin[frequency] + d22S*Sin[
                                      2  frequency];
  W40= a40
      +d44S*Sin[4  frequency]  +   d44C*Cos[
                                      4  frequency]
      +d43S*
                                    Sin[3  frequency]  +   d43C*Cos[
                                      3  frequency]
      +d42S*Sin[2  frequency]  +   d42C*Cos[
                                      2  frequency]
      +d41S*Sin[   frequency]  +   d41C*Cos[   frequency];

  HQ2=DiagonalMatrix[A4*W40 + A2*W20 -A0*W0 ]; 
}];

Ha = 
                        HRF + HQ1 + HQ2;
{HT, Tp} = 
                        Eigensystem[N[Ha]];        T = Transpose[Tp];
mat = 
                        DiagonalMatrix[
                          Exp[-I*Δt*HT]];
ρ11 = T . mat . Tp;   ρ22 = T . Conjugate[mat] . Tp;   (Subscript[ρ, 0]) = ρ11 . (Subscript[ρ, 0]) . ρ22; 
];
  
pulsevar:=(
coef4 += 1;
pulseFunction[coef4];
);
  
pulsefixe:=Module[{pulseDuration},
 For [b = 1,  
                        b ≤ 
                          nbPulseIncrement, b++, {
   pulseDuration = b + coef4; 
   pulseFunction[
                            pulseDuration];
}];
 coef4 += 
                        nbPulseIncrement;
];
  
  
acq0 := (s[0] = (Subscript[ρ, 0]););
  
  
acq[
                  z_] := (s[z] = (Subscript[ρ, 0]););
  
  
store[sr_]:=(
d[sr] = (Subscript[ρ, 0]); 
indicRotor = coef4;
);
  
  
recall[r_]:=(
(Subscript[ρ, 0]) = 
                    d[r]; 
coef4 = indicRotor; 
);
  
  
filterElt[elements_]:=Module[{ρ6, long, elt1, elt2},
ρ6 = Table[0, {dim}, {dim}];
long = 
                        Length[elements];
For [i = 1, 
                        i ≤ long,  i++, {
 
                          elt1 = elements[[i, 1]]; 
 
                          elt2 = 
                            elements[[i, 
                                2]]; 
 ρ6[[elt1, 
                                elt2]] = (Subscript[ρ, 0])[[elt1, 
                                elt2]];   
}];
(Subscript[ρ, 0]) = ρ6;
];
  
  
filterCoh[coh_]:=Module[{ρ6, quanta, qq, long, x, elt1, elt2},
ρ6 = Table[0, {dim}, {dim}];

For [i = 1, 
                        i ≤ 
                          Length[
                            coh], i++, {
  quanta = 
                            coh[[i]]; 
  qq = 
                            Abs [quanta]; 
  long = 
                            dim - qq; 
  x = 
                            Range[1, long]; 
  For [
                            z = 1, 
                            z ≤ 
                              long, z++, {
    If[
                                quanta ≥ 0, {elt1 = x[[z]]; 
                                  elt2 = 
                                    x[[z]] + 
                                      qq},   
                             {elt2 = x[[z]]; 
                                  elt1 = 
                                    x[[z]] + 
                                      qq}]; 
    ρ6[[elt1, elt2]] = (Subscript[ρ, 0])[[elt1, 
                                    elt2]]; 
}]; 
}];

(Subscript[ρ, 0]) = ρ6;
];
  
  
  
tabgraph[nom_String]:=(  
Print["**********************************"];
For[a = 0,  
                      a ≤ 
                        np,   a++, pas[a] = a Δt;];

tableau = 
                      Table[{t, pas[t], NumberForm[h[t], 10]}, {t, 0, np}];  
\(Print[
                      TableForm[
                        tableau, \ 
\ \ TableHeadings -> {None, {"Rang", "t(μs)", "intensity"}}]]\);

Print["**********************************"];
\(ListPlot[
                      Table[{pas[t], h[t]}, {t, 0, 
                          np}], 
\ \ \ \ PlotJoined -> 
                        True, 
\ \ \ \ AxesLabel -> {"t (μs)", "Intensity (A.U.)"}, \ \ \ \ 
\ \ \ \ PlotStyle -> {Hue[
                            0.1]}, 
\ \ \ \ TextStyle -> {FontFamily -> "Times", FontSize -> 12}]\);
Print["**********************************"];

Clear[excel];
\(excel[
                        data_List]\  := 
\ \ \ \ \ Module[\ \ {\ file1\  = \ 
                            OpenWrite[
                              nom]\ }, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Scan[\((\ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ WriteString[file1, \ 
                                First[#]]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Scan[
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ WriteString[file1, \ "	", \ #] &, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Rest[#]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ WriteString[
                                file1, \ "
"]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ )\) &, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ data
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Close[file1]
\ \ \ \ ]\);
excel[tableau];
);
  
  
crystal[αd_, βd_,  γd_]:=Module[{αangle, βangle, γangle, 
                      elt11, elt22},
αangle = αd*π/180;       βangle = βd*π/180;       γangle = γd*π/180;

For[i = 0,  
                        i ≤ np,  i++, h[i] = 0];

elt11 = 
                        detectelt[[1, 1]];  elt22 = detectelt[[1, 2]];

(Subscript[ρ, 0]) = startOperator;  

coefv20Rot[αangle, βangle, γangle];
If[
                        quadrupoleOrder  ==  2,  
                        coefwx0Rot[αangle, βangle, γangle]];

coef4 = 0;    fsimulation;

For[i = 0,  
                        i ≤ 
                          np,  i++, 
    h[i] = 
                          N[Im[s[i][[elt11, elt22]] ] ] /
                            numberOfGammaAngles ];
];
  
  
powder[rep_, maxγ_]:=
  Module[{αangle, βangle, γangle, 
                      proba, elt11, elt22, fileSize, xtalFile},
For[i = 0,  i ≤ np,  i++, 
                        h[i] = 0];

elt11 = 
                        detectelt[[1, 1]];  elt22 = detectelt[[1, 2]];

xtalFile  = 
                        ReadList[rep ,  {Number,  Number,  Number}];
fileSize =  
                        Length[xtalFile];
Print["powderFile: rep", 
                        fileSize, "_simp"];

For[j = 1, j ≤ fileSize, j++,{
  Print[j,  "/",  fileSize];
  αangle = 
                                xtalFile[[j, 1]]*π/180;  
  βangle = 
                                xtalFile[[j, 2]]*π/180;
  proba = xtalFile[[j, 3]];

  For [g = 0, g <  maxγ, g++,{
    (Subscript[ρ, 0]) = startOperator;  
    γangle = ((2 g π)/max)γ;
    coefv20Rot[αangle, βangle, γangle];
    If[quadrupoleOrder == 2, 
                                        coefwx0Rot[αangle, βangle, γangle]];

    coef4 = 0;    fsimulation;

    For[k = 0, k ≤ np, k++,h[k]=h[k]+ proba*s[k]];
 }];
 }];

For[i = 0,  
                        i ≤ 
                          np,  i++, 
    h[i] = 
                          N[Im[h[i][[elt11, elt22]] ] ] /
                            numberOfGammaAngles ];
];
  
  
coefv20Rot[αeuler_, βeuler_, γeuler_] := 
              Module[{a1, a2, b1, 
                  b2}, 
(Subscript[c, 2  α]) = 
                  Cos[2  αeuler];       (Subscript[s, 2  α]) =  
                  Sin[2  αeuler];  
(Subscript[c, β]) = 
                  Cos[βeuler];             (Subscript[s, β]) =  
                  Sin[βeuler];              
(Subscript[c, 2  β]) = 
                  Cos[2  βeuler];       (Subscript[s, 2  β]) = 
                  Sin[2  βeuler]; 
(Subscript[c, γ]) = 
                  Cos[γeuler];             (Subscript[s, γ]) =  
                  Sin[γeuler];             
(Subscript[c, 2  γ]) = 
                  Cos[2  γeuler];       (Subscript[s, 2  γ]) =  
                  Sin[2  γeuler];                                               
a1 = -η*(Subscript[s, 2  α])*
                    (Subscript[s, β])/(Sqrt[3]);                   b1 = -(-3 + η *(Subscript[c, 2  α])) 
                    (Subscript[s, 2  β])/(2 (Sqrt[ 3])); 
a2 = -η*(Subscript[c, β])*
                    (Subscript[s, 2  α])/(Sqrt[6]);                   b2 = -(η*(Subscript[c, 2  α])* (3 + (Subscript[c, 2  β])) + 
                          6  ((Subscript[s, β])^2))/(4 (Sqrt[6])); 
d2S = 
                  a2 *(Subscript[c, 2  γ]) +  
                    b2*(Subscript[s, 2  γ]);              d1S = 
                  a1*(Subscript[c, γ])  + b1*(Subscript[s, γ]); 
d2C = 
                  a2 *(Subscript[s, 2  γ])  -  
                    b2 *(Subscript[c, 2  γ]);            d1C = 
                  a1*(Subscript[s, γ])  - b1*(Subscript[c, γ]);
];
  
 (*-- 
                provides the parameters d4xC and d4xS of W40 --*) 
  
coefwx0Rot[αeuler_, βeuler_, γeuler_] := 
  Module[{a22, b22, a21, b21, a41, a42, a43, a44, 
                b41, b42, b43, b44}, 
(Subscript[c, 4  α]) = 
                Cos[4  αeuler];             
(Subscript[c, 2  α]) = 
                Cos[2  αeuler];       (Subscript[s, 2  α]) =  
                Sin[2  αeuler]; 
        
(Subscript[c, 4  β]) = 
                Cos[4  βeuler];       (Subscript[s, 4  β]) = 
                Sin[4  βeuler]; 
(Subscript[c, 2  β]) = 
                Cos[2  βeuler];       (Subscript[s, 2  β]) = 
                Sin[2  βeuler]; 
(Subscript[c, β]) = 
                Cos[βeuler];             (Subscript[s, β]) =  
                Sin[βeuler]; 
 
(Subscript[c, γ]) = 
                Cos[γeuler];            (Subscript[s, γ]) =  
                Sin[γeuler];             
(Subscript[c, 2  γ]) = 
                Cos[2  γeuler];      (Subscript[s, 2  γ]) =  
                Sin[2  γeuler];                                             
(Subscript[c, 3  γ]) = Cos[3  γeuler];      (Subscript[s, 3  γ]) =  
                Sin[3  γeuler];            
(Subscript[c, 4  γ]) = 
                Cos[4  γeuler];      (Subscript[s, 4  γ]) =  
                Sin[4  γeuler];    
 
a22 = -(Sqrt[2/7])  η*(Subscript[c, β])* 
                  (Subscript[s, 2  α]) ;         b22 = -(η*
                          (Subscript[c, 2  α]) (3 + (Subscript[c, 2  β])) + 
                        ((Subscript[s, β])^2) (-3 + (η^2)))/(2 (Sqrt[14])); 
a21 = -(2/(Sqrt[7])) η *(Subscript[s, β])* 
                  (Subscript[s, 2  α]) ;     b21 = (-3 - 
                      2  (Subscript[c, 2  α]) *η + (η^2)) 
                  (Subscript[s, 2  β])/(2 (Sqrt[7])); 
d22S = 
                a22*(Subscript[c, 2  γ]) +  
                  b22*(Subscript[s, 2  γ]);       d21S = 
                a21 *(Subscript[c, γ])  +  
                  b21 *(Subscript[s, γ]); 
d22C = 
                a22 *(Subscript[s, 2  γ])  -  
                  b22 *(Subscript[c, 2  γ]);    d21C = 
                a21* (Subscript[s, γ])  -  
                  b21 *(Subscript[c, γ]); 
a40 = ((-(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)); 
a41 =  ((Sqrt[5/7])/72) η*
                  (Subscript[s, 2  α])*  (Subscript[s, β]) (15 + 
                      21  (Subscript[c, 2  β]) + 
                      14 η *(Subscript[c, 2  α])* 
                        ((Subscript[s, β])^2)); 
b41 = ((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  β])); 
a42 = - ((Sqrt[5/14]) /
                        18) η * 
                  (Subscript[c, β])*(Subscript[s, 2  α]) (-9 + 
                      21  (Subscript[c, 2  β]) + 
                      14 η *(Subscript[c, 2  α])*
                        ((Subscript[s, β])^2)); 
b42 = ((-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)); 
a43 = -((Sqrt[35])/
                        72) η (-3 - 
                      9  (Subscript[c, 2  β]) + η *
                        (Subscript[c, 2  α])*(5 + 3 *(Subscript[c, 2  β])))*
                  (Subscript[s, 2  α])*  (Subscript[s, β]); 
b43 = -((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  β]); 
a44 = -((Sqrt[35/2])/72) η *(Subscript[c, β])* 
                  (Subscript[s, 2  α]) (η *
                        (Subscript[c, 2  α])* (3 + (Subscript[c, 2  β])) + 
                      6 ((Subscript[s, β])^2)); 
b44 = -((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)); 
d41S = 
                a41*(Subscript[c, γ])   + b41*(Subscript[s, γ]);        d41C = 
                a41*(Subscript[s, γ])  - b41*(Subscript[c, γ]) ;       d42S = 
                a42*(Subscript[c, 2  γ]) + b42*(Subscript[s, 2  γ]);    d42C = 
                a42*(Subscript[s, 2  γ]) - 
                  b42*(Subscript[c, 2  γ]) ;   
d43S = 
                a43*(Subscript[c, 3  γ]) + b43*(Subscript[s, 3  γ]);    d43C = 
                a43*(Subscript[s, 3  γ]) - 
                  b43*(Subscript[c, 3  γ]) ; 
d44S = 
                a44*(Subscript[c, 4  γ]) + b44*(Subscript[s, 4  γ]);    d44C = 
                a44*(Subscript[s, 4  γ]) - 
                  b44*(Subscript[c, 4  γ]);
];
 (* 
            Save these functions in a notebook called QUADRUPOLE *) 
Save["QUADRUPOLE", run, radiofreq, pulse, pulseFunction, 
            pulsevar, pulsefixe, acq0, acq, store, recall, filterElt, 
            filterCoh, tabgraph, crystal, powder, coefv20Rot, 
            coefwx0Rot];