function urm6_QFM_example
%
% QFM = Quark Flavour Model
%
% This function is a full example of URMT under Pythagoras conditions,
% i.e. n=2 exponent, starting with the the trivial 'Pythaogrean double'
% solution (0,1,1) for X+, i.e. x=0, y^2=z^2.
%
% The solution is the same as the 'Numeric Example' given in Chapter 14
% of the fifth book 'A Quark Flavour Model'.
%
% The initial solution is time-evolved by 2s to give X-
% as the ubiquitous (3,4,5) Pythag. triple, albeit in the
% (4,-3,-5) form. The X+ is invariant to time-evolution and
% remains as (0,1,1).
%
% The solution is then further frequency-evolved by 2Hz (f3=2) to give
% the non-trivial X+ = (-12,5,13) Pythsgorean triple.
%
% The time, frequency-evolved solution solution is then lifted to
% a full URM6 (6x6) six-quark flavour model, using evolutionary
% parameters t3=0,t4=1,t5=2,t6=3;
%
% Of course, the numbers given are just one example, and the initial
% solution and the time and frequency domain values are completely arbitrary.
% So too, the 'lifted' URM6 solution is parameterised by four
% arbitrary parameters 't3,t4,t5,t6' all of which can be varied to give
% different numerical solutions to the same quark flavour model.
%
% Output to the screen is only for check values that should be zero in all
% cases if no errors. 'zero' here means both exactly zero or
% numerically close to zero, i.e. of order < 1e-10,
% allowing for machine numerical error. Indeed, most URM6 values shown
% are of order 1e-15 and not exactly zero.
%
% There is also output of selected eigenvectors and matrices to the file
% urm6example.txt. Note that several eigenvectors and matrices are scaled by
% a factor of 2 or 3, e.g. quark charge Q is printed as 3*Q, to remove the 1/3
% charge and make all output integer, not rational.
%
% Only the standard set of eigenvectors and matrices is shown, not
% the sum and difference sets, this is purely to keep the script
% relatively simple and short (still over 1200 lines).
%
% Modification History
% 04/11/2024 V1.0 comment tidy, f2 renamed f3 in line with Book 5
% 27/09/2024 As per urmt_QFM_example but with k=0, l=1, s=0, u=-1
% to give the X+=(0,1,1) trival Pythag solution and as used in
% the Numeric Example, ch14 of Book 5, p289.
%
% www.urmt.org
% 04/11/2024 RJM V1.0

urmt_filename='urm6example.txt'; % print output filename

% start with the URM3 parametric solution in temrs of just the two integers k,l
% using the standard parameterisation of a Pythagorean triple, as below.
% x=12;y=5;z=13;  % coordinates x,y,z, Pythag triple (12,5,13)
% trivial X+ solution (0,1,1)
l = 1;
k = 0;
% l = 3; % earlier test values, not used
% k = 2; % ditto
x = 2*k*l;
y = l^2 - k^2;
z = l^2 + k^2;

C = 1; % unity eigenvalue C to enforce unity roots, e.g. u^2==C^2 mod p

% assign s,u as a particular solution to the Diophantine equation C=ks - lu
% trivial X+ solution (0,1,1)
s = 0;
u = -1;
% s = 2; % earlier test values, not used
% u = 1; % ditto
C - (k*s - l*u) % check the URM3 Pythag solution for s,u given k,l

% P=-7;Q=+4;R=-8; % unity roots - dynamical variables P,Q,R
P = -(k*s + l*u);
Q = (l*s - k*u);
R = -(l*s+k*u);

% assign conjugate dynamical variables using Pythagoras conditions
PB=P;QB=Q;RB=-R;

% the URM3 unity root matrix A0
A0=[0 R QB;RB 0 P;Q PB 0];

K=P*PB+Q*QB+R*RB; % Kinetic term K
V=(P*Q*R + PB*QB*RB)/C; % Potential term V = 0 under Pythagoras conditions

% Determine the invariant eigenvalue C from the dynamical conservation equation
C2=K+V; % C^2 = kinetic plus potential energy
Cb=sqrt(C2);
Cb - C % verify it is still unity

% a few useful constants
twoc=2*C;
twoc2=2*C2;
root2 = sqrt(2);
root2c = root2*C;

% Determine the scale factors (alpha,beta,gamma)
% from the divisibilty relations. Don't assume initial (x,y,z)
% solution is zero in any coordinate

if x==0, alpha=0; else
    alpha = (C2 - P*PB)/x;
    % it is possible for the 'unity' root P to be zero for trivial
    % Pythag solutions, i.e. when the Pythag triple is really a double of
    % the form (x,0,x) or (0,y,y) so only do this check when non-trivial
    if P ~= 0
      mod(P^2,x) - C2 % verify unity root P
    end
end;

if y==0, beta=0; else
    beta = (C2 - Q*QB)/y;
    if Q ~= 0
      mod(Q^2,y) - C2 % verify unity root Q for non-trivial Q, see P above
    end
end;

if z==0, gamma = 0; else
    gamma = (C2 - R*RB)/z;
    if R ~= 0
    tmp = mod(R^2,z); % verify unity root R for non-trivial R, see P above
    if tmp > 0, tmp=tmp-z; end;
    tmp + C2 % = 0
    end
end

% check alpha, beta, gamma agree with the URM3 parametric solution
% Note that this (alpha,beta,gamma) is the dual Pythag triple (-4,-3,5), i.e.
% the classic (3,4,5) to within sign and swap of x and y
alpha - -2*2*u
beta - (u^2 - s^2)
gamma - (u^2 + s^2)

% construct the three column eigenvectors
Xp=[x;y;z];
Xm=[alpha;beta;-gamma];
X0=[P;-Q;R];

% form the 'T operator' (Minkowski form)
I3=eye(3);
T3=I3;
T3(3,3)=-1;

% form the URMT conjugate, row eigenvectors
mX=urmt_conj(Xp);
pX=urmt_conj(Xm);
oX=urmt_conj(X0);

% check non-zero eigenvector inner products
pX*Xp - twoc2
oX*X0 - C2 % the dynamical conservation equation (DCE)
mX*Xm - twoc2

% orthogonal eigenvector inner products, should all be zero by definition
pX*Xm % Pythagoras
mX*Xp % Pythagoras
oX*Xm
oX*Xp
pX*X0
mX*X0

% eigenvector equations
A0*Xp - C*Xp % plus  eigenvector Xp, eigenvalue +C
A0*X0        % zero  eigenvector X0, eigenvalue  0
A0*Xm + C*Xm % minus eigenvector Xm, eigenvalue -C

% conjugate eigenvector equations
pX*A0 - C*pX
oX*A0
mX*A0 + C*mX

A0b = (Xp*pX - Xm*mX)/twoc; % spectral decomposition form of A0
A0b - A0
A02 = A0*A0; % A0^2
A0*(A02 - C2*I3) % Cayley-Hamilton polynomial, eigenvalues +C,0,-C

Ap=[0 z -y;-z 0 x;-y x 0]; % plus matrix A+
Apb = (X0*mX - Xp*oX)/C;   % outer product construction of A+
Apb - Ap % confirm the above two matrix constructions are the same
Ap*Ap*Ap % Ap^3 = 0 => three zero eigenvalues by Cayley-Hamilton

% A+ eigenvector equations
Ap*Xp % = 0 annihilator A+ * X+ = 0
Ap*Xm - 2*C*X0 % = A+ * X- = 2C*X0
Ap*X0 + C*Xp % = 0 because A+ * X0 = -C*X+

% A+ conjugate eigenvector equations
pX*Ap + twoc*oX
mX*Ap
oX*Ap - C*mX

% minus matrix, A- = dual of A+
Am=[0 -gamma -beta;gamma 0 alpha;-beta alpha 0];
Amb = -(X0*pX - Xm*oX)/C; % outer product construction of A-
Amb - Am % confirm the above two are the same
Am*Am*Am % A-^3 = 0 => three zero eigenvalues by Cayley-Hamilton

% Am eigenvector equations
Am*Xm % = 0 annihilator
Am*Xp + 2*C*X0 % = A0+ * X- = 2C*X0
Am*X0 - C*Xm % = 0 because A0+ * X0 = -C*X+

% A- conjugate eigenvector equations
pX*Am
mX*Am - twoc*oX
oX*Am + C*pX

% commutators - same as angular momentum, three state +,0,-
commutator(A0,Ap) - C*Ap % [A0,A+] = C.A+
commutator(A0,Am) + C*Am % [A0,A-] = -C.A-
commutator(Ap,Am) - twoc*A0 % [A+,A-] = 2C.A0

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% print the eigenvectors and matrices
% for the intial solution, See (14.3), (14.4) p290 and (14.8) p291
%
fout = fopen( urmt_filename, 'wt' ); % open output file

fprintf(fout,'\nURM3 Initial solution for k=%d, l=%d, s=%d, u=%d\n\n',...
  k,l,s,u);
  % Eigenvectors - note they are the trivial forms
print_nx1vector( Xp,X0,Xm,    '       X+        X0        X-', fout,3);

% A matrices
print_nxnarray( Ap, 'A+ =', fout,3);
print_nxnarray( A0, 'A0 =', fout,3);
print_nxnarray( Am, 'A- =', fout,3);

fclose(fout);
% End print
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Eigenvector temporal (time-domain) evolution
%
% Time evolve the URMT solution and show it also has an
% exponentiated form, but without the trig and with a finite
% cutoff.
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t3 = 2.0; % evolve 2 time units
% t3 = 0.0; % no evolution, zero test value

Xpt = Xp; % X+ time invariant, i.e. no time evolution
X0t = -t3*Xp + X0; % evolve X0
Xmt = -t3*t3*Xp + 2*t3*X0 + Xm; % evolve X-

% evolve A0
A0t = A0 - t3*Ap;

% evolved eigenvector equations - same form as unevolved
A0t*X0t
A0t*Xpt - C*Xpt % OK
A0t*Xmt + C*Xmt % OK

% Exponentiated evolution operator
Et = expm(t3*Ap/C);

% Inverse, time-reversed evolution operator
Eti = expm(-t3*Ap/C);
Et*Eti - I3 % check inverse - sanity.

% Exponential operator evolution of unity root matrix A0
Et*A0*Eti - A0t

% series expansion form of evolution operator with finite cut-off at (A+)^2 term
% since (A+)^n = 0 for n>2
Et - (I3 + t3*Ap/C + (t3*Ap/C)^2/2 )

% Verify individual eigenvector evolution equations
Et*Xp - Xpt
Et*X0 - X0t
Et*Xm - Xmt

% evolution in eigenvector matrix form
XX = [Xp X0 Xm];     % inital eigenvector matrix
XXt = [Xpt X0t Xmt]; % evolved eigenvector matrix
XXt - Et*XX

% obtain evolved eigenvector conjugates
pXt = transpose(T3*Xmt);
oXt = transpose(T3*X0t);
mXt = transpose(T3*Xpt);

% check conjugate eigenvector evolution equations
mXt - mX % invariant
pXt - pX*(I3 - t3*Ap/C + (t3*Ap/C)^2/2 )
oXt - (t3*mX + oX);

% check conjugate eigenvector evolution equations using inverse

% evolution operator
pX*Eti - pXt
mX*Eti - mXt
oX*Eti - oXt

% check conjugate eigenvector equations
pXt*A0t - C*pXt
oXt*A0t
mXt*A0t + C*mXt

% check inner products
pXt*Xpt - twoc2
mXt*Xmt - twoc2
pXt*Xmt
mXt*Xpt

% check evolved, outer product forms of A matrices
A0t - (Xpt*pXt - Xmt*mXt)/twoc
Amt = (Xmt*oXt - X0t*pXt)/C;
Apt = (X0t*mXt - Xpt*oXt)/C;

Apt - Ap % A+ time-invariant
Amt - (-t3*t3*Ap + 2*t3*A0 + Am)

% check operator evolution equations
Apt - Et*Ap*Eti % invariant
A0t - Et*A0*Eti
Amt - Et*Am*Eti

Apt*Xpt % = 0 annihilator
Apt*Xmt - 2*C*X0t % = A0+ * X- = 2C*X0
Apt*X0t + C*Xpt

Amt*Xmt % = 0 annihilator
Amt*Xpt + 2*C*X0t
Amt*X0t - C*Xmt % = 0 because A0+ * X0 = -C*X+

% conjugate evolved relations
pXt*Apt + twoc*oXt
mXt*Apt
oXt*Apt - C*mXt

pXt*Amt
mXt*Amt - 2*C*oXt
oXt*Amt + C*pXt

% check commutator relations of evolved A matrices
commutator(A0t,Apt) - C*Apt
commutator(A0t,Amt) + C*Amt
commutator(Apt,Amt) - twoc*A0t

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% test of unity root properties for evolved dynamical vars P,Q,R
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% extract dynamical vars from evolved A0
P = A0t(2,3);
Q = A0t(3,1);
R = A0t(1,2);

X0t - [P;-Q;R] % time-evolved X0 eigenvector

% extract coordinates from A+ (time-invariant as per X+)
x = Apt(2,3);
y = -Apt(1,3);
z = Apt(1,2);

Xpt - [x;y;z]

% extract scale factors from evolved A-
alpha = Amt(2,3);
beta = -Amt(1,3);
gamma = Amt(2,1);

Xmt - [alpha;beta;-gamma] % time-evolved X- eigenvector

% check divisibility equations
(C2 - P*P) - alpha*x
(C2 - Q*Q) - beta*y
(C2 + R*R) - gamma*z

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% print the eigenvectors and matrices
%
fout = fopen( urmt_filename, 'at' ); % open output file

fprintf(fout,'\nURM3 Time-domain Evolution time t3=%6.1f\n\n',t3);

print_nx1vector( Xpt,X0t,Xmt,    '       X+        X0        X-', fout,3 );

% A matrices
print_nxnarray( Apt, 'A+ =', fout,3);
print_nxnarray( A0t, 'A0 =', fout,3);
print_nxnarray( Amt, 'A- =', fout,3);

fclose(fout);
% End print
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% now evolve in the freq domain
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f3 = 2;

A0f = A0t - f3*Amt;

Xmf = Xmt; % invariant
X0f = -f3*Xmt + X0t;
Xpf = -f3*f3*Xmt + 2*f3*X0t + Xpt;

pXf = urmt_conj(Xmf);
mXf = urmt_conj(Xpf);
oXf = urmt_conj(X0f);

A0f*X0f
A0f*Xpf - C*Xpf
A0f*Xmf + C*Xmf

Amf = (Xmf*oXf - X0f*pXf)/C;
Apf = (X0f*mXf - Xpf*oXf)/C;
Amf - Amt % invariant

pXf*Amf
mXf*Amf - twoc*oXf
oXf*Amf + C*pXf

pXf*Apf + twoc*oXf
mXf*Apf
oXf*Apf - C*mXf

commutator(A0f,Apf) - C*Apf
commutator(A0f,Amf) + C*Amf
commutator(Apf,Amf) - twoc*A0f

Ef = expm(-f3*Amt/C); % !! NOTE the sign change from time domain !!
Efi = expm(f3*Amt/C); % !! ditto !!

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% print the eigenvectors and matrices
%
% fout = fopen( urmt_filename, 'at' ); % open output file

% x9 A matrices
% print_nxnarray( Ef,  'Ef =', fout,3);
% print_nxnarray( Efi, 'Efi =', fout,3);

% fclose(fout);
% End print
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Ef*Efi - I3 % sanity!

Ef*A0t*Efi - A0f
Ef*Amt*Efi - Amf
Ef*Apt*Efi - Apf

Amf - Amt % invariant
A0f - (-f3*Amt + A0t)
Apf - (-f3*f3*Amt + 2*f3*A0t + Apt)

% !! NOTE the sign change from time domain
% f3 is the negative of the time domain
Ef - (I3 - f3*Amt/C + (f3*Amt/C)^2/2 )
Ef*Xpt - Xpf
Ef*X0t - X0f
Ef*Xmt - Xmf

XX=[Xpt, X0t, Xmt];
XXf=[Xpf, X0f, Xmf];
Ef*XX - XXf

XX=[Xmt, X0t, Xpt]; % reverse ordering for STM
XXf=[Xmf, X0f, Xpf];
Ef*XX - XXf

XXfi=[mXf; 2*oXf; pXf]/twoc2; % inverse eigenvector matrix
XXfi*XXf - I3

% state transition matrix equivalents
Mf=[1 -f3 -f3*f3;0 1 2*f3;0 0 1];
XXf - XX*Mf
Ef - XXf*Mf*XXfi % equivalence of  STM and exponential
det(Ef) - 1
det(Mf) - 1
[vf df] = eigs(Ef);
[vfs dfs] = eigs(Mf);

Ef*Xmt - Xmt % invariant - sanity

% reverse ordering - no good for STM, but works for Et
% XX=[Xp, X0, Xm];
% XXf=[Xpf, X0f, Xmf];
% XXf - XX*Mf

% URMT Pauli equivalents
% these are the spin defns with complex Ad0, and give the standard
% commutation relations all with 'i' in the structure factors
% As0f = (Apf + Amf)/2;
% Ad0f = -(Apf - Amf)/2i;

% A1f = A0f;  % A1 == Ax == A0
% A2f = Ad0f; % A2 == Ay == -Ad0
% A3f = As0f; % A3 == Az == As0

% these are the spin defns with complex A2, and give the standard
% commutation relations all with 'i' in the structure factors

% commutator(A2f,A3f) - (C*1i*A1f)
% commutator(A1f,A2f) - (C*1i*A3f)
% commutator(A3f,A1f) - (C*1i*A2f)

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% print the eigenvectors and matrices
%
fout = fopen( urmt_filename, 'at' ); % open output file

fprintf(fout,'\nURM3 Frequency-domain Evolution time f3=%6.1f\n\n',f3);

print_nx1vector( Xpf,X0f,Xmf,    '       X+        X0        X-', fout,3 );

% A matrices
print_nxnarray( Apf, 'A+ =', fout,3);
print_nxnarray( A0f, 'A0 =', fout,3);
print_nxnarray( Amf, 'A- =', fout,3);

fclose(fout);
% End print
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% reset the solution for the next stage - lift URM3 to URM6
A0t = A0f;
Amt = Amf;  % invariant
Apt = Apf;

Xpt = Xpf;
X0t = X0f;
Xmt = Xmf;

pXt = pXf;
mXt = mXf;
oXt = oXf;

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Time-domain lifting of URM3 solution to URM6, i.e. 3x3 to 6x6
%
% Lift the URM3 solution to URM6 using four temporal params t3,t4,t5,t6
%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

A0_i = A0t; % save initial URM3 A matrices
Ap_i = Apt;

% initial URM6 unity root matrix
A0 = zeros(6,6);
A0(4:6,4:6) = A0_i; % A0 = x axis, keep in mind
A0t=zeros(6,6);

Ap = zeros(6,6); % URM6 Ap matrix
Ap(4:6,4:6) = Ap_i;

% the four URM6 evolution time parameters, example values
% t3=0; t4=1; t5=1; t6=1;
t3=0; t4=1; t5=2; t6=3; % as per ch14, book 5, see p305 for actual soln

% populate URM6 A0 matrix as per evolved form
A0t(4:6,4:6) = A0_i - t3*Ap_i; % A0 = x axis, keep in mind
A0t(1,4:6)=-t6*mXt;
A0t(4:6,1)=t6*Xpt;
A0t(2,4:6)=-t5*mXt;
A0t(4:6,2)=t5*Xpt;
A0t(3,4:6)=-t4*mXt;
A0t(4:6,3)=t4*Xpt;

% URM6 initial eigenvectors
Xp = [0;0;0;Xpt];
Xm = [0;0;0;Xmt];
X0A = [0;0;0;X0t];
X0B = [0;0;C;0;0;0];
X0C = [0;C;0;0;0;0];
X0D = [C;0;0;0;0;0];

% URM6 evolved eignevectors
Xpt = Xp; % invariant
X0At = -t3*Xp + X0A;
X0Bt = -t4*Xp + X0B;
X0Ct = -t5*Xp + X0C;
X0Dt = -t6*Xp + X0D;
X0t = (X0At + X0Bt + X0Ct + X0Dt)/2; % combined zero eignevector X60

% Full four temporal params, time-domain solution
Xmt = -(t3*t3+t4*t4+t5*t5+t6*t6)*Xp + 2*t3*X0A  + 2*t4*X0B + 2*t5*X0C + 2*t6*X0D + Xm;

% URM6 evolved eigenvector equations
A0t*X0At
A0t*X0Bt
A0t*X0Ct
A0t*X0Dt
A0t*Xpt - C*Xpt % OK
A0t*X0t
A0t*Xmt + C*Xmt % OK

% URM6 T operator (Minkowski form)
I6 = eye(6,6);
T6 = I6;
T6(6,6) = -1;

% conjugate URM6 eigenvectors
pXt = transpose(T6*Xmt);
mXt = transpose(T6*Xpt);
oXAt = transpose(T6*X0At);
oXBt = transpose(T6*X0Bt);
oXCt = transpose(T6*X0Ct);
oXDt = transpose(T6*X0Dt);
oXt = transpose(T6*X0t);
oXt - (oXAt + oXBt + oXCt + oXDt)/2 % combined conjugate eigenvector

A0tb = (Xpt*pXt - Xmt*mXt)/twoc; % outer product form of A0
A0tb - A0t % check identical to original

% check evolved conjugate eigenvector equations
pXt*A0t - C*pXt
mXt*A0t + C*mXt
oXt*A0t

% check evolved eigenvector inner products
pXt*Xpt - twoc2
mXt*Xmt - twoc2

oXAt*X0At - C2
oXBt*X0Bt - C2
oXCt*X0Ct - C2
oXDt*X0Dt - C2
oXt*X0t - C2 % normalised form

% check eigenvector orthogonality
pXt*Xmt
mXt*Xpt

oXAt*Xpt
oXBt*Xpt
oXCt*Xpt
oXDt*Xpt

oXAt*Xmt
oXBt*Xmt
oXCt*Xmt
oXDt*Xmt

oXAt*X0Bt
oXAt*X0Ct
oXAt*X0Dt

oXBt*X0Ct
oXBt*X0Dt
oXCt*X0Dt

oXt*Xpt
oXt*Xmt

Apt = (X0t*mXt - Xpt*oXt)/C; % raising operator A+ outer product form
Amt = -(X0t*pXt - Xmt*oXt)/C; % dual of A+ - the minus is because -C = dual(C)

% ********************************************************
% REDUNDANT - NOT REQUIRED
% evolve the URM6 solution a further two time units

if 0 % disable, not required

tn = 2;

Mtn = expm(tn*Apt/C); % exponentiated evolution matrix

XX = [Xpt X0At X0Bt X0Ct X0Dt Xmt]; % construct eigenvector matrix
XXn = Mtn*XX; % evolve the eigenvectors

% extract the four 'zero' eigenvectors and check inner products
X0Atn = XXn(:,2);
oXAtn = transpose(T6*X0Atn);
oXAtn*X0Atn - C2

X0Btn = XXn(:,3);
oXBtn = transpose(T6*X0Btn);
oXBtn*X0Btn - C2

X0Ctn = XXn(:,4);
oXCtn = transpose(T6*X0Ctn);
oXCtn*X0Ctn - C2

X0Dtn = XXn(:,5);
oXDtn = transpose(T6*X0Dtn);
oXDtn*X0Dtn - C2

% check the evolved eigenvectors against the evolution equations
X0Atn - (-(tn/2)*Xpt + X0At); % invariant
X0Btn - (-(tn/2)*Xpt + X0Bt);
X0Ctn - (-(tn/2)*Xpt + X0Ct);
X0Dtn - (-(tn/2)*Xpt + X0Dt);

end % disable, not required
% ****************************

% check standard evolved eigenvector relations - same form as unevolved
A0t*Xpt - C*Xpt
A0t*X0t
A0t*Xmt + C*Xmt

Apt*Xpt % = 0 annihilator
Apt*X0t + C*Xpt
Apt*Xmt - 2*C*X0t % = A0+ * X- = 2C*X0

Amt*Xmt % = 0 annihilator
Amt*Xpt + 2*C*X0t
Amt*X0t - C*Xmt % = 0 because A0+ * X0 = -C*X+

% conjugate evolved relations
pXt*A0t - C*pXt
oXt*A0t % = 0 annihilator
mXt*A0t + C*mXt

pXt*Apt + twoc*oXt
mXt*Apt % = 0 annihilator
oXt*Apt - C*mXt

pXt*Amt % = 0 annihilator
mXt*Amt - 2*C*oXt
oXt*Amt + C*pXt

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% print the URM6 evolved eigenvectors and matrices
%

fout = fopen( urmt_filename, 'at' ); % open output file

fprintf(fout,'\nURM6 time-domain lifting t3=%0.1f,t4=%0.1f,t5=%0.1f,t6=%0.1f\n\n',...
    t3,t4,t5,t6);

print_nx1vector( Xpt,2*X0t,Xmt,    '       X+      2*X0        X-', fout,6);
print_nx1vector_two( X0At,X0Bt,    '       X0A     X0B', fout,6);
print_nx1vector_two( X0Ct,X0Dt,    '       X0C     X0D', fout,6);

% A matrices
print_nxnarray( 2*Apt, '2*A+ =', fout,6);
print_nxnarray( 2*A0t, '2*A0 =', fout,6);
print_nxnarray( 2*Amt, '2*A- =', fout,6);

fclose(fout);

% End print
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% assign time-evolved eigenvector forms back to their std names,
% i.e. initialise the eignevector solution to the time-evolved values
Xp = Xpt;
Xm = Xmt;
X0A = X0At;
X0B = X0Bt;
X0C = X0Ct;
X0D = X0Dt;

% time-evolved conjugate eigenvectors
pX = pXt;
mX = mXt;
oXA = oXAt;
oXB = oXBt;
oXC = oXCt;
oXD = oXDt;

% now use this initial solution for the URM6 quark flavour model

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Quarks and charge operators

% For third component of isospin, operator Iz, these need to be the
% eigenvectors of the URM6 A0 (~A) matrix

u = Xpt/root2c; % up
d = Xmt/root2c; % down
s = X0At/C; % strange
c = X0Ct/C; % charm
b = X0Bt/C; % bottom
t = X0Dt/C; % top

% In this scheme, Izb = -Iz, anti-particles ub = URMT conj d,
% and db = URMT conj u

ub = urm6_conj(d); % ub ~ px
db = urm6_conj(u); % db ~ mx
sb = urm6_conj(s); % self-conjugate zero eigenvectors, as per X0A,B,C,B
cb = urm6_conj(c); % self-conjugate
bb = urm6_conj(b); % self-conjugate
tb = urm6_conj(t); % self-conjugate

% check anti-quark eigenvector defns versus equivalent
% URMT-conjugated eigenvectors
ub - pXt/root2c
db - mXt/root2c
sb - oXAt/C
cb - oXCt/C
bb - oXBt/C
tb - oXDt/C

% anti-quark, quark inner products

% up quark inner products
ub*u - 1
ub*d % orthogonal = 0
ub*s
ub*c
ub*b
ub*t

% down quark inner products
db*d - 1
db*s
db*c
db*b
db*t

% strange quark inner products
sb*s - 1
sb*c
sb*b
sb*t

% charm quark inner products
cb*c - 1
cb*b
cb*t

% bottom quark inner products
bb*b - 1
bb*t

% top quark inner products
tb*t - 1

Iz = (u*ub - d*db)/2; % operator third, z component of isospin
Izb = -Iz; % anti-particle form of Iz (negation of Iz in URMT)

Im = d*ub; % isospin lowering operator I-
Ip = u*db; % isospin raising operator I+

Ix = (Ip + Im)/2;  % x-axis spin operator
Iy = (Ip - Im)/2i; % y-axis spin operator

% check isospin x,y,z operator commutator relations
commutator(Ix,Iy) - 1i*Iz
commutator(Iy,Iz) - 1i*Ix
commutator(Iz,Ix) - 1i*Iy

Iz - A0t/twoc % check same as scaled URM6 matrix

% check isospin eigenvector equations
Iz*u - u/2 % +1/2 isospin eigenvalue up quark
Iz*s % zero isospin
Iz*d + d/2 % -1/2 isospin eigenvalue down quark

% check isospin raising, lowering operator commutator relations
commutator(Iz,Ip) - Ip
commutator(Iz,Im) + Im
commutator(Ip,Im) - 2*Iz

% check isopsin raising operator action
Ip*u % annihilates an up quark
Ip*s % annihilates a strange quark
Ip*d - u % raise a down to an up

% check isospin lowering operator action
Im*u - d % lower an up to a down
Im*s % annihilates a strange quark
Im*d % annihilates a down quark

% anti-particle raising, lowering operators - negation of std forms
Ipb = -Ip;
Imb = -Im;

% check anti-particle raising equations
ub*Ipb + db % anti-up to -ve anti-down
sb*Ipb % annihilates an anti-strange quark
db*Ipb % annihilates an anti-down quark

% check anti-particle lowering equations
ub*Imb % annihilates an anti-up quark
sb*Imb % annihilates an anti-strange quark
db*Imb + ub % anti-down to -ve anti-up

% construct 'charge' operators
S = -s*sb; % strangeness
C = c*cb; % charm
Bq = -b*bb; % bottom
Tq = t*tb; % top
B = (u*ub + d*db + s*sb + c*cb + b*bb + t*tb)/3; % Baryon
Y = (u*ub + d*db - 2*(s*sb + b*bb) + 4*(c*cb + t*tb))/3; % hypercharge
Q = (2*u*ub + 2*c*cb + 2*t*tb - d*db - s*sb - b*bb)/3;   % charge

Y - (B + S + C + Bq + Tq) % check consistent defn hypercharge operator
Q - (Iz + Y/2) % check consistent defn charge operator

% baryon number eigenvector equations - eigenvalue = quark baryon number = 1/3
B*u - u/3
B*d - d/3
B*s - s/3
B*c - c/3
B*b - b/3
B*t - t/3

% baryon number anti-quarks
Bb = -B;  % URMT anti-particle operator is negation of std form
% anti-quark baryon number = -1/3
ub*Bb + ub/3
db*Bb + db/3
sb*Bb + sb/3
cb*Bb + cb/3
bb*Bb + bb/3
tb*Bb + tb/3

% charge
Q*u - (2/3)*u % up charge = 2/3
Q*d + (1/3)*d % down charge = -1/3
Q*s + (1/3)*s % strange charge = -1/3
Q*c - (2/3)*c % charm charge = 2/3
Q*b + (1/3)*b % bottom charge = -1/3
Q*t - (2/3)*t % top charge = 2/3

% anti-quark charge = negation of above
Qb = -Q;  % URMT anti-particle operator is negation of std form
ub*Qb + (2/3)*ub
db*Qb - (1/3)*db
sb*Qb - (1/3)*sb
cb*Qb + (2/3)*cb
bb*Qb - (1/3)*bb
tb*Qb + (2/3)*tb

% hypercharge
Y*u - (1/3)*u
Y*d - (1/3)*d
Y*s + (2/3)*s
Y*c - (4/3)*c
Y*b + (2/3)*b
Y*t - (4/3)*t

% hypercharge anti-quarks negation of above
Yb = -Y; % URMT anti-particle operator is negation of std form
ub*Yb + (1/3)*ub
db*Yb + (1/3)*db
sb*Yb - (2/3)*sb
cb*Yb + (4/3)*cb
bb*Yb - (2/3)*bb
tb*Yb + (4/3)*tb

% URMT quark operator equivalents

% Baryon operator URMT eigenvector form
Burmt = (1/(6*C2))*(Xpt*pXt + Xmt*mXt) + ...
    (1/(3*C2))*(X0At*oXAt + X0Bt*oXBt) + ...
    (1/(3*C2))*(X0Ct*oXCt + X0Dt*oXDt);
B - Burmt % check consistent defn baryon number operator

Yurmt = (1/(3*C2))*(Xp*pX/2 + Xm*mX/2 - 2*(X0A*oXA + X0B*oXB) + ...
    4*(X0C*oXC + X0D*oXD) ); % URMT eigenvector defn of hypercharge operator
Y - Yurmt % check hypercharge

Qurmt = (1/(3*C2))*(Xp*pX + 2*X0C*oXC + 2*X0D*oXD - ...
(Xm*mX/2 + X0A*oXA + X0B*oXB) ); % URMT eigenvector defn of charge operator
Q - Qurmt % check charge

% alternate URMT eigenvector defn of hypercharge operator -
% should be the same as above
Qurmt = (1/(3*C2))*(Xp*pX + 2*X0C*oXC + 2*X0D*oXD) - ...
    (1/(6*C2))*(Xm*mX + 2*X0A*oXA + 2*X0B*oXB);
Q - Qurmt % check alternate charge def

% URMT eigenvector defn of third component of isospin
Izurmt = (Xp*pX - Xm*mX)/(4*C2);
Iz - Izurmt % check isospin

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% print the final URM6 representation
%
fout = fopen( urmt_filename, 'at' ); % open output file

fprintf(fout,'\nURM6 quark flavour representation\n\n');

% print quark eigenvectors with root2 factor and rounding - the latter
% only required for small numerical error, the root2 factor essentially
% converts them to integer to within small, machine numerical error
% uf = round(u*root2);
% df = round(d*root2);
% sf = round(s*root2);
% cf = round(c*root2);
% bf = round(b*root2);
% tf = round(t*root2);

% ubf = round(ub*root2);
% dbf = round(db*root2);
% sbf = round(sb*root2);
% cbf = round(cb*root2);
% bbf = round(bb*root2);
% tbf = round(tb*root2);

% print quark column vectors
print_nx1vector_two( u*root2,d*root2,    '   up*sqrt(2) , down*sqrt(2)', fout,6);
print_nx1vector_two( s,c,    ' strange  charm', fout,6);
print_nx1vector_two( t,b,    '      top    bottom', fout,6);
% print_nx1vector( uf,df,sf,    '     up       down     strange', fout,6);

% operators
print_nxnarray( 2*Iz, '2*Iz =', fout,6);
print_nxnarray( 2*Im, '2*Im =', fout,6);
print_nxnarray( 2*Ip, '2*Ip =', fout,6);

print_nxnarray( 3*Q, '3*Q =', fout,6);
print_nxnarray( 3*Y, '3*Y =', fout,6);
print_nxnarray( 3*B, '3*B =', fout,6);

% print_nxnarray( A0t, 'A0 =', fout,6);
% print_nxnarray( Amt, 'A- =', fout,6);

fclose(fout);

% End print
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('end');

end

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Auxiliary functions as used above
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Commutation
function C = commutator(A,B)
    C = A*B - B*A;
end

% URMT/URM3 conjugation
function B=urmt_conj(A)
    T3 = [1 0 0;0 1 0;0 0 -1];
    B=transpose(T3*A);
end

% URM6 conjugation
function B=urm6_conj(A)
    T6 = eye(6,6);
    T6(6,6) = -1;
    B=transpose(T6*A);
end

% Print an nxn array 'A' (size 'nn' x 'nn')
% with a message 'start_msg' that may be blank, e.g. start_msg='',
% and optional output to file handle 'fout', that must be opened in
% write/append mode beforehand.
function print_nxnarray( A, start_msg, fout, nn )
% determine if any element is imaginary

imagf = any( (abs( imag(A(:)) ) > 1e-7) > 0 );

fprintf(fout,'%s\n',start_msg);

if imagf > 0
    for ii=1:nn
    for jj=1:nn
        % zero the matrix entry if considered zero within numeric error 1e-7 here
        if abs( A(ii,jj) ) < 1e-7
          A(ii,jj) = 0;
        end
        % if abs( real( A(ii,jj) ) ) > 1e-7
            % dp = abs( imag(A(ii,jj)) - floor(imag(A(ii,jj))) );
            dp = abs( real(A(ii,jj)) - round(real(A(ii,jj))) );
            if dp > 1e-7
                fprintf(fout,'   %6.1f ',real(A(ii,jj)) );
            else
                fprintf(fout,'   %6.0f ',real(A(ii,jj)) );
            end
        % else
            fprintf(fout,'       0 ');
        % end

        if abs( imag( A(ii,jj) ) ) > 1e-7
            oneunit = abs( abs(imag(A(ii,jj))) - 1 );
            if oneunit < 1e-7
                if imag(A(ii,jj)) > 0
                    fprintf(fout,'+      i,' );
                else
                    fprintf(fout,'+     -i,' );
                end
            else
                dp = abs( imag(A(ii,jj)) - round(imag(A(ii,jj))) );
                if dp > 1e-7
                    fprintf(fout,'+ %6.1fi,',imag(A(ii,jj)) );
                else
                    fprintf(fout,'+ %6.0fi,',imag(A(ii,jj)) );
                end
            end
        else
            fprintf(fout,'        ,');
        end
    end% jj
    fprintf(fout,'\n');
    end% ii
else
    for ii=1:nn
        for jj=1:nn
          % zero the matrix entry if considered zero within numeric error 1e-7 here
          if abs( A(ii,jj) ) < 1e-7
            A(ii,jj) = 0;
          end
            % dp = abs( imag(A(ii,jj)) - floor(imag(A(ii,jj))) );
            dp = abs( real(A(ii,jj)) - round(real(A(ii,jj))) );
            if dp > 1e-7
                fprintf(fout,'   %6.1f ',real(A(ii,jj)) );
            else
                fprintf(fout,'   %6.0f ',real(A(ii,jj)) );
            end
        end% jj
        fprintf(fout,'\n');
    end% ii
end
fprintf(fout,'\n');
end

% print three vectors vp,v0,vm, of length 'nn'-element in a
% side-by-side column format, with a message 'start_msg' that
% may be blank, e.g. start_msg='',
% and optional output to file handle 'fout', that
% must be opened in write mode/append beforehand.
function print_nx1vector( vp,v0,vm, start_msg, fout,nn )
fprintf(fout,'%s\n\n',start_msg);
  % check for any fractional component
  % dp = 0.0;
  % dp = dp + abs( real(vp) - floor(real(vp)) );
  % dp = dp + abs( real(v0) - floor(real(v0)) );
  % dp = dp + abs( real(vm) - floor(real(vm)) );
  % dpf = sum(dp(:));
for ii=1:nn
  % check for any fractional component
  dp = 0.0;
  dp = dp + abs( real(vp(ii)) - round(real(vp(ii))) );
  dp = dp + abs( real(v0(ii)) - round(real(v0(ii))) );
  dp = dp + abs( real(vm(ii)) - round(real(vm(ii))) );
  if dp > 1e-7
  % if dpf > 1e-7
    fprintf(fout,'   %6.1f ',vp(ii),v0(ii),vm(ii) );
  else
    fprintf(fout,'   %6.0f ',vp(ii),v0(ii),vm(ii) );
  end
  fprintf(fout,'\n');
end
fprintf(fout,'\n');
end

% print two vectors vp,v0, of length 'nn'-element in a
% side-by-side column format, with a message 'start_msg' that
% may be blank, e.g. start_msg='',
% and optional output to file handle 'fout', that
% must be opened in write mode/append beforehand.
function print_nx1vector_two( vp,v0, start_msg, fout,nn )
fprintf(fout,'%s\n\n',start_msg);
for ii=1:nn
  % check for any fractional component
  dp = 0.0;
  dp = dp + abs( real(vp(ii)) - round(real(vp(ii))) );
  dp = dp + abs( real(v0(ii)) - round(real(v0(ii))) );
  if dp > 1e-7
    fprintf(fout,'   %6.1f ',vp(ii),v0(ii) );
  else
    fprintf(fout,'   %6.0f ',vp(ii),v0(ii) );
  end
  fprintf(fout,'\n');
end
fprintf(fout,'\n');
end

% print four vectors v0a,v0b,v0c,v0d, of length 'nn'-element in a
% side-by-side column format, with a message 'start_msg' that
% maybe be blank, e.g. start_msg='',
% and optional output to file handle 'fout', that
% must be opened in write mode/append beforehand.
function print_nx1_zerovectors( v0a,v0b,v0c, v0d, start_msg, fout,nn )
fprintf(fout,'%s\n\n',start_msg);
for ii=1:nn
  % check for any fractional component
  dp = 0.0;
  dp = dp + abs( real(v0a(ii)) - round(real(v0a(ii))) );
  dp = dp + abs( real(v0b(ii)) - round(real(v0b(ii))) );
  dp = dp + abs( real(v0c(ii)) - round(real(v0c(ii))) );
  dp = dp + abs( real(v0d(ii)) - round(real(v0d(ii))) );
  if dp > 1e-7
    fprintf(fout,'   %6.1f ',v0a(ii),v0b(ii),v0c(ii),v0d(ii) );
  else
    fprintf(fout,'   %6.0f ',v0a(ii),v0b(ii),v0c(ii),v0d(ii) );
  endif
  fprintf(fout,'\n');
end
fprintf(fout,'\n');
end

% print a single complex vector 'V', of length 'nn'-element in a
% side-by-side column format, with a message 'start_msg' that
% may be blank, e.g. start_msg='',
% and optional output to file handle 'fout', that
% must be opened in write mode/append beforehand.
function print_nx1vector_cmplx( V, start_msg, fout, nn )

% determine if any element is imaginary
imagf = any( (abs( imag(V(:)) ) > 1e-7) > 0 );

fprintf(fout,'%s\n',start_msg);

if imagf > 0
    for ii=1:nn
        if abs( real( V(ii) ) ) > 1e-7
            % dp = abs( imag(V(ii)) - floor(imag(V(ii))) );
            dp = abs( real(V(ii)) - round(real(V(ii))) );
            if dp > 1e-7
                fprintf(fout,'   %6.1f ',real(V(ii)) );
            else
                fprintf(fout,'   %6.0f ',real(V(ii)) );
            end
        else
            fprintf(fout,'       0 ');
        end

        if abs( imag( V(ii) ) ) > 1e-7
            oneunit = abs( abs(imag(V(ii))) - 1 );
            if oneunit < 1e-7
                if imag(V(ii)) > 0
                    fprintf(fout,'+      i,' );
                else
                    fprintf(fout,'+     -i,' );
                end
            else
                dp = abs( imag(V(ii)) - round(imag(V(ii))) );
                if dp > 1e-7
                    fprintf(fout,'+ %6.1fi,',imag(V(ii)) );
                else
                    fprintf(fout,'+ %6.0fi,',imag(V(ii)) );
                end
            end
        else
            fprintf(fout,'        ,');
        end
        fprintf(fout,'\n');
    end% ii
else
    for ii=1:nn
        % dp = abs( imag(V(ii)) - floor(imag(V(ii))) );
        dp = abs( real(V(ii)) - round(real(V(ii))) );
        if dp > 1e-7
            fprintf(fout,'   %6.1f ',real(V(ii)) );
        else
            fprintf(fout,'   %6.0f ',real(V(ii)) );
        end
        fprintf(fout,'\n');
    end% ii
end

fprintf(fout,'\n');

end
