function magn=MfromT(B0,T,wn,wm,w1,wmas,Azz,DD,T1ea,T2ea,T1eb,T2eb,T1n,T2n,dt,N,Ia,Ib,Ig,Na,Nb,Ng,Je)
%   Here dt is the time step; N1 = number of rotor cycles with HFI off;
%   N2 = number of rotor cycles with HFI on.

B0=str2num(B0);
T=str2num(T);
wn=str2num(wn);
wm=str2num(wm);
w1=str2num(w1);
wmas=str2num(wmas);
Azz=str2num(Azz);
DD=str2num(DD);
T1ea=str2num(T1ea);
T2ea=str2num(T2ea);
T1eb=str2num(T1eb);
T2eb=str2num(T2eb);
T1n=str2num(T1n);
T2n=str2num(T2n);
dt=str2num(dt);
N=str2num(N);
Ia=str2num(Ia);
Ib=str2num(Ib);
Ig=str2num(Ig);
Na=str2num(Na);
Nb=str2num(Nb);
Ng=str2num(Ng);
Je=str2num(Je);



format long e
% spin operators, parameters of the system and things like that

ai=sqrt(-1);

sx=[0 0.5; 0.5 0];
sy=[0 -0.5*sqrt(-1); 0.5*sqrt(-1) 0];
sz=[0.5 0; 0 -0.5];
odin=[1 0; 0 1];

sp=[0 1; 0 0];
sm=[0 0; 1 0];

Hz1=zeros(8,8);
Hz2=zeros(8,8);
Hzn=zeros(8,8);

Hx1=zeros(8,8);
Hx2=zeros(8,8);
Hy1=zeros(8,8);
Hy2=zeros(8,8);
Hxn=zeros(8,8);
Hyn=zeros(8,8);
bigone=kron(odin,kron(odin,odin));


Hz1=kron(sz,kron(odin,odin));
Hz2=kron(odin,kron(sz,odin));
Hzn=kron(odin,kron(odin,sz));

Hx1=kron(sx,kron(odin,odin));
Hx2=kron(odin,kron(sx,odin));
Hxn=kron(odin,kron(odin,sx));
Hy1=kron(sy,kron(odin,odin));
Hy2=kron(odin,kron(sy,odin));
Hyn=kron(odin,kron(odin,sy));

Hz=Hz1+Hz2;

gammae=1.76086*1e11;
gxx=2.0032;
gyy=2.0026;
gzz=2.0032;             % g-tensor in MS

gxx2=2.0061;
gyy2=2.0022;
gzz2=2.0088;             

Bf=1.76*1.054*B0/(1.38*T);    % E/kT

Dpas=[DD/2 0 0; 0 DD/2 0; 0 0 -DD];
Apas=[Azz/2 0 0; 0 Azz/2 0; 0 0 -Azz];


magn=zeros(N/100,4);
%magn=zeros(N,4);
weicoe=0;
iiA=0;
for zIa=1:Na
    for zIb=1:Nb
        for zIg=1:Ng
            zbb1=0.5*(pi/Nb)+(zIb-1)*(pi)/Nb;
            weicoe=weicoe+sin(zbb1);
       end
    end
end

%                        for Initr=1:3
iiA=iiA+1
%   Hamiltonian, spin 1 is "e", spin 2 is "N"
    H0=zeros(8,8);
    
%    Nz=Initr-2;
Nz=0;
    
fx=1/(2.0023)*gxx*B0*gammae+2*pi*18.8*1e6*Nz;
fy=1/(2.0023)*gyy*B0*gammae+2*pi*92.4*1e6*Nz;
fz=1/(2.0023)*gzz*B0*gammae+2*pi*18.2*1e6*Nz;

fx2=1/(2.0023)*gxx2*B0*gammae+2*pi*18.8*1e6*Nz;
fy2=1/(2.0023)*gyy2*B0*gammae+2*pi*92.4*1e6*Nz;
fz2=1/(2.0023)*gzz2*B0*gammae+2*pi*18.2*1e6*Nz;

gpas=[fx 0 0; 0 fy 0; 0 0 fz];
gpas2=[fx2 0 0; 0 fy2 0; 0 0 fz2];

% Orientation of g-tensor 1
aa1=0.5*(2*pi/Na)+(Ia-1)*(2*pi)/Na-pi;
bb1=0.5*(pi/Nb)+(Ib-1)*(pi)/Nb;
gg1=0.5*(2*pi/Ng)+(Ig-1)*(2*pi)/Ng-pi;



%aa1=2*pi/3;
%bb1=pi/6;
%gg1=0;
% Orientation of g-tensor 2 with respect to g-1
aa2=pi*102/180;
bb2=pi*108/180;
gg2=pi*124/180;
% Orientation of DDI with respect to g-1
bbe=pi*33/180;
gge=pi/4;
% Orientation of HFI with respect to g-1
bbf=0;
ggf=0;

thetaM=0.9553166181245093; %magic angle!!!


%initial populations
ro=zeros(8,8);
Zs=exp(Bf/2)+exp(-Bf/2);
for i=1:8
ro(i,i)=0.5*exp(-Bf*Hz(i,i))/Zs;
end

%ro(1,1)=0;
   aaa=trace(ro);
   ro=ro/aaa;



rho=zeros(8,8);

Mr=round((2*pi/wmas)/dt);

tic
% LOOP OVER ROTOR CYCLE; HFI is on
Lpr2=eye(64,64);

for ifi=1:Mr

    time=(ifi-1)*dt;
    wtr=wmas*time;
%   Euler matrix for g-tensor 1

R11=eurot(aa1,bb1,gg1);
R22=eurot(aa2,bb2,gg2);
Rd=eurot(0,bbe,gge);
Ra=eurot(0,bbf,ggf);

Rlf=eurot(wtr,thetaM,0);
Rlf1=inv(Rlf);

g1=Rlf*R11*gpas*inv(R11)*Rlf1;
g2=Rlf*R11*(R22)*gpas2*inv(R22)*inv(R11)*Rlf1;

w1e=g1(3,3);
w2e=g2(3,3);

Dlf=Rlf*R11*inv(Rd)*Dpas*(Rd)*inv(R11)*Rlf1;
DDI=Dlf(3,3);

Alf=Rlf*Ra*Apas*inv(Ra)*Rlf1;
Azz1=Alf(3,3);
Azz2=Alf(3,1);
Azz3=Alf(3,2);


% Liouville space calcn
    H00=(w1e-wm)*Hz1+(w2e-wm)*Hz2-wn*Hzn+Azz1*Hz1*Hzn+Azz2*Hz1*Hxn+Azz3*Hz1*Hyn+DDI*(2*Hz1*Hz2-Hx1*Hx2-Hy1*Hy2)+Je*(Hz1*Hz2+Hx1*Hx2+Hy1*Hy2)+0.5*w1*(Hx1+Hx2);
    
[ev,vv]=eig(H00);
        ev1=inv(ev);


% now relaxation
Hz1a=ev1*Hz1*ev;
Hz2a=ev1*Hz2*ev;
Hzna=ev1*Hzn*ev;
Hza=Hz1a+Hz2a;

Hx1a=ev1*Hx1*ev;
Hx2a=ev1*Hx2*ev;
Hxna=ev1*Hxn*ev;

Hy1a=ev1*Hy1*ev;
Hy2a=ev1*Hy2*ev;
Hyna=ev1*Hyn*ev;


R1=zeros(8,8);
R2=zeros(8,8);
W2=zeros(8,8);
%   Evaluating matrix of relaxation of populations in the EB

for i=1:8
    for j=1:8
        ep=Hza(i,i)-Hza(j,j);
        Ef=exp(-ep*Bf/2)/(exp(ep*Bf/2)+exp(-ep*Bf/2));
        
        R1(i,j)=(1/T1ea)*(Hx1a(i,j)*Hx1a(j,i)+Hy1a(i,j)*Hy1a(j,i)+Hz1a(i,j)*Hz1a(j,i))+(1/T1eb)*(Hx2a(i,j)*Hx2a(j,i)+Hy2a(i,j)*Hy2a(j,i)+Hz2a(i,j)*Hz2a(j,i))+(1/T1n)*(Hxna(i,j)*Hxna(j,i)+Hyna(i,j)*Hyna(j,i)+Hzna(i,j)*Hzna(j,i));
        R1(i,j)=R1(i,j)*Ef;
        R2(i,j)=(1/T2ea)*abs(Hz1a(i,i)-Hz1a(j,j))^2+(1/T2eb)*abs(Hz2a(i,i)-Hz2a(j,j))^2+(1/T2n)*abs(Hzna(i,i)-Hzna(j,j))^2;
        W2(i,j)=exp(-dt*R2(i,j)-dt*ai*vv(i,i)+dt*ai*vv(j,j));
    end
end

for i=1:8
    R1(i,i)=0;
    W2(i,i)=0;
end

for i=1:8
    for j=1:8
    if (abs(i-j)>0);
        R1(i,i)=R1(i,i)-R1(j,i);
    end
    end
end
W1=expm(R1*dt);

% now super-matrices as 8*8*8*8 and then as 64*64

Rfull=zeros(64,64);
% Relaxation: only transitions between levels and diagonal phase relaxation
for ib=1:8
    for jb=1:8
        
        II=(ib-1)*8+ib;
        JJ=(jb-1)*8+jb;
        IJ=(ib-1)*8+jb;
        
    Rfull(II,JJ)=W1(ib,jb);
    if(abs(ib-jb)>0);
        Rfull(IJ,IJ)=W2(ib,jb);
    end
    end
end

SA=kron(ev,transpose(ev1));
SB=transpose(kron(ev,transpose(ev1)));

  Lpr2=SA*Rfull*SB*Lpr2;
        
end
toc

%   initial super-vector
mp0=zeros(64);
for ib=1:8
    for jb=1:8
        iii=(ib-1)*8+jb;
mp0(iii)=ro(ib,jb);
    end
end

mp=mp0;

% Loop for rotor cycles with HFI on
for in2=100:100:N
    in22=in2/100;
        time=(2*pi/wmas)*(in2);
    magn(in22,1)=time;
    mp=Lpr2^100*mp;
    DM=zeros(8,8);
for ib=1:8
    for jb=1:8
        iii=(ib-1)*8+jb;
DM(ib,jb)=mp(iii);
    end
end    
        magn(in22,2)=(magn(in22,2)+real(trace(Hz1*DM))*sin(bb1))/weicoe;
        magn(in22,3)=(magn(in22,3)+real(trace(Hz2*DM))*sin(bb1))/weicoe;
        magn(in22,4)=(magn(in22,4)+real(trace(Hzn*DM))*sin(bb1))/weicoe;
end


%             end


%for in2=1:N
%        magn(in2,2)=magn(in2,2)/weicoe;
%        magn(in2,3)=magn(in2,3)/weicoe;
%        magn(in2,4)=magn(in2,4)/weicoe;
%end

%save levrel.dat levs -ASCII
save magnrel.dat magn -ASCII
%    plot(levs(:,1),levs(:,2),levs(:,1),levs(:,3),levs(:,1),levs(:,4),levs(:,1),levs(:,5),levs(:,1),levs(:,6),levs(:,1),levs(:,7),levs(:,1),levs(:,8),levs(:,1),levs(:,9));


%save adiRD2s.dat Tad -ASCII
