5A1335 Symmetries in Physics

Program files

 
To use the Matlab files, cut along the indicated lines, save them as filename.m 
in the directory where you run Matlab.  

Use the 'help' facility to find out how to use them. 

The function files given here are:

cg.m -  calculates SU(2) Clebsch-Gordan coefficients
dj.m -  calculates reduced SU(2) matrix elements
mult.m  -  multiplicities in reduction of products of SU(2) irreps
sundim.m  - dimension of an irrep of SU(N) defined by a Young diagram
symdim.m -  dimension of an irrep of the permutation group S(N) 
terms.m - terms allowed by Pauli principle for n electrons

There are some additional auxiliary functions (separated by %%--------)
included inside the function files.  You can separate them into m-files if
you prefer, e.g. if running  an older version of Matlab. 

The files have been tested using Matlab 5.2.1. and 6.5
There is no guarantee that they give the correct answer in all cases. 

Updated 23 Oct 2004

-------------------------- cut here ------------------------------

function f=cg(j1,j2,m1,m2);
%> cg.m  calculates numerical SU(2) Clebsch-Gordan coefficients, 
%> 
%> call cg(j1,j2,m1,m2); 
%> Input: a consistent set of integers or half-integers (j1,j2,m1,m2)
%> Output: the first column = values of j = |j1 - j2|, ... , j1 + j2 .
%> the second column is the C-G coefficient 

jmin=max(abs(m1+m2),abs(j1-j2));
jmax=j1+j2;
jj=[jmin:jmax]';
mm=m1+m2;

x1=fact(j1+j2-jj).*fact(j2+jj-j1).*fact(jj+j1-j2)./fact(jj+j1+j2+1);
x2=fact(j1+m1)*fact(j1-m1)*fact(j2+m2)*fact(j2-m2).*fact(jj+mm).*fact(jj-mm);
x3=zeros(size(jj));

for t=0:jmax
x3=[x3, ((-1)^t)./(fact(t).*fact(jj-j2+m1+t).*fact(jj-j1-m2+t).*fact(jmax-jj-t).*fact(j1-m1-t).*fact(j2+m2-t))];
end
x4=sum(x3')';
f=[jj,sqrt(2*jj+1).*sqrt(x1.*x2).*x4];

%%--------------------------------------------------------------------
function f=fact(x);
%>  calculates the factorials fact(x) of 
%> a matrix with integer entries, sets infinity for negative values.
%> Call: fact(X)
%> Input: X = vector or matrix, integer values
%> Output: a vector or matrix of same size.
%> Compare the standard gamma function of MATLAB. 
%> We want to avoid all NaN outputs!!!

maxx=max(max(x));	y=zeros(size(x));

% Iterate the faculty for elements of values 1,2,...,maxx.

	for n=0:maxx-1

	z=~(x - maxx + n);	y=y+z; 	y=(maxx-n)*y;

	end 

z=~x; % Now deal with the zero arguments, answer is 1.

y=y+z;

% Now deal with all the remaining matrix elements  -
% the arguments are negative, the answer should be infinity -
% set each matrix element separately

w=~y; [w1,w2]=find(w); nn=length(w1);

	for m=1:nn

	y(w1(m),w2(m))=inf;

	end

f=y;

-------------------------- cut here ------------------------------
 
function f=dj(j,x);
%> dj.m defines a (2j+1)x(2j+1) reduced rotation matrix
%> called  d^j_m1_m2, defined as in Edmonds by the Jacobi polynomial.
%> Call  dj(j,x), 
%> Input: j positive integer or half-integer
%> x = beta in [0,pi] is the angle.  
%> Output: a (2j+1)x(2j+1) real matrix
%> Uses evalpol, jacob 

% GL 961123

%% First deal with half-integer case, renormalize to integers
q=ceil(j-0.1)-floor(j+0.1); 
jj=round(j+0.5*q);

x1=sin(0.5*x);
x2=cos(0.5*x);

A=[];

for m1=q:jj

m2=[-m1+q:m1]';
w0=x2.^(m1+m2-q).*x1.^(m1-m2);
w1=sqrt(fact(jj+m1-q).*fact(jj-m1)./fact(jj+m2-q)./fact(jj-m2));
w=evalpol(jacob(m1-m2,m1+m2-q,jj-m1),x1^2);
w=w0.*w1.*w;
ch=cos((m2-m1)*pi);

A(jj+m1+1-q,jj+1-m1:jj+m1+1-q)=w';
A(jj-m1+1:jj+m1+1-q,jj-m1+1)=flipud(w);
A(jj+1-m1:jj+m1+1-q,jj+m1+1-q)=ch.*w;
A(jj-m1+1,jj-m1+1:jj+m1+1-q)=flipud(ch.*w)';

end

f=A';

%%------------------------------------------------------------------
function f=fact(x);
%> calculates the factorials fact(x) of 
%> a matrix with integer entries, sets infinity for negative values.
%> Call: fact(X)
%> Input: X = vector or matrix, integer values
%> Output: a vector or matrix of same size.
%> Compare the standard gamma function of MATLAB. 
%> We want to avoid all NaN outputs!!!

maxx=max(max(x));	y=zeros(size(x));

% Iterate the faculty for elements of values 1,2,...,maxx.

	for n=0:maxx-1

	z=~(x - maxx + n);	y=y+z; 	y=(maxx-n)*y;

	end 

z=~x; % Now deal with the zero arguments, answer is 1.

y=y+z;

% Now deal with all the remaining matrix elements  -
% the arguments are negative, the answer should be infinity -
% set each matrix element separately

w=~y; [w1,w2]=find(w); nn=length(w1);

	for m=1:nn

	y(w1(m),w2(m))=inf;

	end

f=y;

%%-------------------------------------------------------------------
function f=jacob(a,b,n);
%>  calculates the coefficients of the powers of (1-X)/2 in 
%> the Jacobi polynomials of order n and integer parameters.
%> Call: jacob(a,b,n) where n = non-neg integer, 
%> a, b = COLUMN vectors with INTEGER components, a must be non-negative.
%> One of a,b may be scalar, otherwise vectors of same length.
%> It returns a matrix of dim(a) rows, each defining a polynomial.
%> The order of coefficients is from n to 0.
%> For evaluation write something like z=jacob(a,b,n)
%> y=evalpol(z,0.5*(1-X)); to give a vector of the size of a.
 
% GL 961123

len=max(length(a),length(b));	llen=[1:len]';

aa=-n; %Paras in hypergeom series = aa,bb,cc
bb=a+b+1+n;
cc=a+1+zeros(size(llen));
f=[]; %  Start on function vector
g=ones(size(llen));	h=ones(size(0:n));
a1=g;	a2=1;

	if n==0
	f=[g,f]; 
	else
	f=[g,f];

		for r=1:n
		a1=a1.*(a+r);
		a2=a2*r;
		g=g.*(aa+r-1).*(bb+r-1)./(cc+r-1)/r;
		f=[g,f];

		end

	f=(a1*h).*f/a2;
	f=round(f);
	end

%%----------------------------------------------------------------------
function f=evalpol(p,x);
%> The function file  does polynomial evaluation for a number 
%> of polynomials given by a matrix p, where each row represents the 
%> coefficients of a polynomial, in descending powers.  
%> Call evalpol(p,x); 
%> Input: p a matrix, x a row vector
%> Output = matrix, # rows given by p, # columns = length of x
%> (Compare MATLAB routine polyval) 
%>
%>  Goran Lindblad

pps=size(p); pp1=pps(1); nc=pps(2); nn=[0:nc-1]'; x1=size(x);

f=[];

if sum(x1)==2

	for n=1:pp1

		y = filter(1,[1,-x],p(n,:));
		y = y(nc);
		f=[f;y];
	
	end

elseif x1(1) == 1

xx = (ones(nc,1)*x).^(nn*ones(size(x)));

f=fliplr(p)*xx;

else

disp('> ERROR: the second argument in  must be a row vector!');

f=[];

end

-------------------------- cut here ------------------------------

function f=mult(x);
%> mult.m is a program which calculates the multiplicities of angular 
%> momenta
%> Call: mult(x),
%> Input:  row vector x integer or half-integer elements >=0.
%> Output: matrix of 3 columns, the first gives the possible angular
%> momenta, the second the number of number of times each value appears,
%> the third is the total multiplicity. 
%> The last line gives the number of angular momenta multiplets 
%> and the total number of states. 

N=length(x);

	if N==1	
	f=[x(1);1];	
	end

jtop=x(1)+x(2);
j0=abs(x(1)-x(2)); %This is zero or 0.5
jj=j0:jtop;
array=[jj;ones(size(jj))];

	if N==2
	f=array';
	end

	for p=3:N
	array=concat(array,x(p));
	end
	
	array=[array;(2*array(1,:)+1).*array(2,:)];
	
f=[array';[ 0 ,sum(array(2,:)),sum(array(3,:))]];

%%-----------------------------------------------------------------
function f = concat(x,y);
%>  concat.m defines concat(x,j), where x is an array
%>  [j1,...,jm;n1,...,nm]; jk and j are angular momenta,
%>  where the jk must be all integer or all half-integer.
%>  The nk must be non-negative integers
x1=x(1,:);	x2=x(2,:);

m=length(x1);	jtop=max(x1)+y;
j0=jtop-floor(jtop); % This is zero or 0.5

jj=j0:jtop; % This is the range of allowed values
w0=zeros(size(jj));
w=w0;

	for p=1:m;
	string=[abs(x1(p)-y):x1(p)+y]-j0+1;
	w1=w0;
	w1(string)=x2(p)*ones(size(string));
	w=w+w1;
	end

f=[jj;w];

--------------------- cut here -------------------------------------------

function f=sundim(N,Y);
%> This function calculates the dimension of the irrep of SU(N)
%> associated with a standard Young diagram Y
%> Call: sundim(N,Y),
%> Input: N = positive integer
%> Y= [y1, y2,..,yM]; positive integers forming a standard Young diagram
%> i.e. y1 >= y2 >=... >= yM, M <= N,
%> Output: an integer = the dimension of the irrep.
%> There is an error message if input is not correct.
%> Noninteger entries are rounded.

nn=round(N);

yy=Y; % must be integer and positive

yy=round(yy);

dd=max(diff(yy)); % must be negative

m1=max(yy); m2=min(yy); 

mm=length(yy); mxx=mm+m1+1;

prod1=1;prod2=1;

if nn < mm | m2 < 1 | dd > 0 

disp('> Disallowed values!');

else

%% calculate associated diagram

y2=ones(1,yy(1));

if mm==1

y2=y2;

else

for iter=2:mm

y2=mad(y2, ones(1,yy(iter)));

end 

y2=sum(y2);

end

mm1=ones(mm,1)*y2 - [0:mm-1]'*ones(1,m1);
mm1=crop(mm1, mxx,0);
mm2=yy'*ones(1,m1) - ones(mm,1)*[1:m1];
mm2=crop(mm2,mxx,0);
mm3=crop(mm1+mm2,3*mxx,1);
prod2=prod(prod(mm3));

for iter=1:mm

xx=1:yy(iter); xx=xx+nn-iter;
prod1=prod1*prod(xx);

end

end

x=prod1/prod2; error=abs(x-round(x));

if error > .1
disp('> Numerical errors!');
end

f=round(x);

%%------------------------------------------------------------
function f = crop(y, ymax, ymin);
%> The function file  crops a real vector or matrix: 
%> it replaces values y > ymax by ymax, values < ymin by ymin. 
%> Call: crop(y, ymax, ymin);
%>
%>  Goran Lindblad

% First make sure ymin < ymax

if ymax < ymin

z=[ymin,ymax];

ymax=z(1); ymin=z(2);

disp('> Interchanging ymax and ymin!');

elseif ymax == ymin

disp('> You have chosen ymax = ymin, a constant function is returned!');
 
end

% Now do the cropping 

wmax= 0.5*(sign(y-ymax) + 1); % 1 if y > ymax, 0 otherwise.
wmin= 0.5*(-sign(y-ymin) + 1); % 1 if y < ymin, 0 otherwise.
f=y.*(1-wmax).*(1-wmin) + ymax*wmax + ymin*wmin;

%%-----------------------------------------------------------------
function f=mad(A,x);
%> The function file  adds a row vector to an existing matrix
%> Call: mad(a,x)
%> Input: a is a matrix, x a row vector.
%> Output: matrix with x as a new last row to a, 
%> expanding the dimensions as needed.
%>
%>  Goran Lindblad 

if isempty(A)

f=x;

elseif isempty(x)

f=A;

else
sA=size(A);
sx=size(x);
if sx(1) > 1
disp('The vector x should be a row vector!');
return
end 
if sx(1)==0
return
end
if sA(2)==sx(2)
A=[A;x];
end
if sA(2) > sx(2)
x=[x,zeros(1,sA(2)-sx(2))];
A=[A;x];
end
if sA(2) < sx(2)
A=[A,zeros(sA(1),sx(2)-sA(2))];
A=[A;x];
end

f=A;

end

--------------------- cut here -------------------------------------------

function f=symdim(Y);
%>-----------------------------------------------------------
%> The file symdim calculates the dimension of the irrep of S(N)
%> (the permutation group of N objects)
%> defined by a standard Young diagram Y, and also the number
%> of group elements in the conjugacy class defined by Y
%> Call: symdim(Y),
%> Y= [y_1, y_2,.. y_r];
%> is a vector of positive integers forming a standard Young diagram
%> i.e. y_1 >= y_2 >=  ... >= y_r,
%>     sum_j y_j = N,
%> Output: the triple [N,D,M] 
%> D = the dimension of the irrep.
%> M = the number elements in the class 
%> There is an error message if input is not correct.
%> Noninteger entries are rounded.
%> G Lindblad 2004
%>------------------------------------------------------------

yy=Y; % components must be integer and positive

yy=round(yy);

nn=sum(yy);

dd=max(diff(yy)); % must be negative or zero

m1=max(yy); m2=min(yy); 

mm=length(yy);

if  m2 < 1 | dd > 0 

disp('> Disallowed values!');

elseif mm==1|m1==1

dd=1;

else 

kk=yy+mm-[1:mm];

p1=detpro(kk);

q1=prod(fact(kk));

dd=fact(nn)*p1/q1;

end

f=[nn,dd,cc(yy)];

%%--------------------------------------------------------
function f=cc(yy); 
%> calculates the number of elements in the class
%> formula (3.28) in Ledermann - Intro to the theory of finite groups
order=sum(yy);
lmax = max(yy);
nn=zeros(1,lmax);
mm=[1:lmax];
for iter=1:lmax
nn(iter) =  length(find(yy-iter == 0));
end 
%> nn(p) is the number of cycles of length p
%> mm is the vector of the length of possible cycles

f = fact(order)/prod(mm.^nn)/prod(fact(nn));

%%---------------------------------------------------
function f=detpro(X);
%> This function calculates the product 
%> product(j< k)[X(j)-X(k)]
%> Call: detpro(X)
%> Input: X a row or column vector, real or complex
%> Output: A real or complex number.

len=length(X);

if len==1

f = X;

else

yy=1; 

for iter=1:len-1

yy=yy*prod(X(iter)-X(iter+1:len));

end

end

f=yy;

%%------------------------------------------------------------------

function f=fact(x);
%>  calculates the factorials fact(x) of 
%> a matrix with integer entries, sets infinity for negative values.
%> Call: fact(X)
%> Input: X = vector or matrix, integer values
%> Output: a vector or matrix of same size.
%> Compare the standard gamma function of MATLAB. 
%> We want to avoid all NaN outputs!!!

maxx=max(max(x));       y=zeros(size(x));

% Iterate the faculty for elements of values 1,2,...,maxx.

        for n=0:maxx-1

        z=~(x - maxx + n);      y=y+z;  y=(maxx-n)*y;

        end 

z=~x; % Now deal with the zero arguments, answer is 1.

y=y+z;

% Now deal with all the remaining matrix elements  -
% the arguments are negative, the answer should be infinity -
% set each matrix element separately

w=~y; [w1,w2]=find(w); nn=length(w1);

        for m=1:nn

        y(w1(m),w2(m))=inf;

        end

f=y;

%  G Lindblad 2004
 
--------------------- cut here -------------------------------------------

function f=terms(N,L);
%> The function f= terms(N,L) finds the terms allowed by Pauli principle. 
%> Input: N = number of electrons, 
%> L = angular momentum of one-electron orbital, same for all.
%> Output: a matrix defining the multiplicity of the allowed terms 
%> of the N-electron state. S = column index, L = row index.
%> Reference: 
%> Curl and Kilpatrick, "Atomic term symbols by group theory",
%> Am J Phys 28(1960), 357--365, in particular eqn (16).

Smax=.5*N; % Maximal spin value, integer or half-integer.
S1=floor(Smax);
S0= Smax-S1;

ww=[]; % starting matrix for solution 

for iter=0:S1 % iteration over the allowed values of total spin

S=S0+iter;
nu=2*L+1;
aa=round(.5*N-S);
bb=round(2*S); 
cc=round(N*L+1 - S^2 - N*(N-2)/4);

% iteration of nominator in eqn (16)
tfactor=inline('[-1,zeros(1,x)]+[zeros(1,x),1]');
p1=1; 
p2=conv(tfactor(nu+1),tfactor(nu+1-aa));
p3=1;

	if aa==0
			p4=1;
	elseif aa==1
			p4=p2;
	elseif aa==2
			p4=conv(tfactor(nu),tfactor(nu));
			p4=conv(p2,p4);
	else
		for n=nu+2-aa:nu
% 		disp(n);
		p1=conv(p1,tfactor(n));
		end
		p4=conv(p1,p1);
		p4=conv(p4,p2);
	end

	if bb==0
			p3=1;
	else
		for n=nu+1-aa-bb:nu-aa
		p3=conv(p3,tfactor(n));
		end

	end

p4=conv(p4,p3); % this is P(q) in (16)
p4=conv(p4,tfactor(1)); % this is the nominator (1-q)P(q) of (16)

% now iteration of the denominator

q1=1; q2=1;

	if aa==0
	q1=1; 
		if bb==1
		q2=tfactor(1);
		else
			for n=1:bb
			q2=conv(q2,tfactor(n));
			end
		end

	elseif aa==1
	q1=tfactor(1);
		if bb==0
		q2=tfactor(2);
		else
			for n=1:aa+bb+1
			q2=conv(q2,tfactor(n));
			end
		q2=deconv(q2,tfactor(bb+1));
		end
		
	else 
		for n=1:aa
		q1=conv(q1,tfactor(n));
		end
		
		if bb==0
			for n=2:aa+1
			q2=conv(q2,tfactor(n));
			end
		else
		q2=q1;
			for n=aa+1:aa+bb+1
			q2=conv(q2,tfactor(n));
			end
		q2=deconv(q2,tfactor(bb+1));
			
		end

	end

q3=conv(q1,q2);
q3=[q3,zeros(1,cc)];
gen=-deconv(p4,q3);
ww=mad(ww,fliplr(gen));

end

f=ww';

%%-----------------------------------------------------------------
function f=mad(A,x);
%> The function file  adds a row vector to an existing matrix
%> Call: mad(a,x)
%> Input: a is a matrix, x a row vector.
%> Output: matrix with x as a new last row to a, 
%> expanding the dimensions as needed.
%>
%>  Goran Lindblad

%  GL 961125

if isempty(A)

f=x;

elseif isempty(x)

f=A;

else
sA=size(A);
sx=size(x);
if sx(1) > 1
disp('The vector x should be a row vector!');
return
end 
if sx(1)==0
return
end
if sA(2)==sx(2)
A=[A;x];
end
if sA(2) > sx(2)
x=[x,zeros(1,sA(2)-sx(2))];
A=[A;x];
end
if sA(2) < sx(2)
A=[A,zeros(sA(1),sx(2)-sA(2))];
A=[A;x];
end

f=A;

end

--------------------- cut here -------------------------------------------
UP