close all

clear all

format short

format compact

% Check orthonormaality

fun00 = @(z) sqrt(2)*sin(1*pi*z)*sqrt(2).*sin(1*pi*z);

fun01 = @(z) sqrt(2)*sin(1*pi*z)*sqrt(2).*sin(2*pi*z);

fun02 = @(z) sqrt(2)*sin(1*pi*z)*sqrt(2).*sin(3*pi*z);

fun03 = @(z) sqrt(2)*sin(1*pi*z)*sqrt(2).*sin(4*pi*z);

fun04 = @(z) sqrt(2)*sin(1*pi*z)*sqrt(2).*sin(5*pi*z);

fun05 = @(z) sqrt(2)*sin(1*pi*z)*sqrt(2).*sin(6*pi*z);

r00 = integral(fun00,0,1);

r01 = integral(fun01,0,1);

r02 = integral(fun02,0,1);

r03 = integral(fun03,0,1);

r04 = integral(fun04,0,1);

r05 = integral(fun05,0,1);

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

fun11 = @(z) sqrt(2)*sin(2*pi*z)*sqrt(2).*sin(2*pi*z);

fun12 = @(z) sqrt(2)*sin(2*pi*z)*sqrt(2).*sin(3*pi*z);

fun13 = @(z) sqrt(2)*sin(2*pi*z)*sqrt(2).*sin(4*pi*z);

fun14 = @(z) sqrt(2)*sin(2*pi*z)*sqrt(2).*sin(5*pi*z);

fun15 = @(z) sqrt(2)*sin(2*pi*z)*sqrt(2).*sin(6*pi*z);

r11 = integral(fun11,0,1);

r12 = integral(fun12,0,1);

r13 = integral(fun13,0,1);

r14 = integral(fun14,0,1);

r15 = integral(fun15,0,1);

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

fun22 = @(z) sqrt(2)*sin(3*pi*z)*sqrt(2).*sin(3*pi*z);

fun23 = @(z) sqrt(2)*sin(3*pi*z)*sqrt(2).*sin(4*pi*z);

fun24 = @(z) sqrt(2)*sin(3*pi*z)*sqrt(2).*sin(5*pi*z);

fun25 = @(z) sqrt(2)*sin(3*pi*z)*sqrt(2).*sin(6*pi*z);

r22 = integral(fun22,0,1);

r23 = integral(fun23,0,1);

r24 = integral(fun24,0,1);

r25 = integral(fun25,0,1);

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

fun33 = @(z) sqrt(2)*sin(4*pi*z)*sqrt(2).*sin(4*pi*z);

fun34 = @(z) sqrt(2)*sin(4*pi*z)*sqrt(2).*sin(5*pi*z);

fun35 = @(z) sqrt(2)*sin(4*pi*z)*sqrt(2).*sin(6*pi*z);

r33 = integral(fun33,0,1);

r34 = integral(fun34,0,1);

r35 = integral(fun35,0,1);

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

fun44 = @(z) sqrt(2)*sin(5*pi*z)*sqrt(2).*sin(5*pi*z);

fun45 = @(z) sqrt(2)*sin(5*pi*z)*sqrt(2).*sin(6*pi*z);

r44= integral(fun44,0,1);

r45= integral(fun45,0,1);

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

fun55 = @(z) sqrt(2)*sin(6*pi*z)*sqrt(2).*sin(6*pi*z);

r55 = integral(fun55,0,1);

for n = 0:5

�� for m = 0:5

����� %n,m

����� fun = @(z,n,m) sqrt(2)*sin((n+1)*pi*z)*sqrt(2).*sin((m+1)*pi*z);

����� over(n+1,m+1) = integral(@(z)fun(z,n,m),0,1);

�� end

end

for n = 0:5

�� for m = 0:5

����� %n,m

����� fun = @(z,n,m) sqrt(2)*sin((n+1)*pi*z).*(z - 1/2)*sqrt(2).*sin((m+1)*pi*z);

����� dipole(n+1,m+1) = (integral(@(z)fun(z,n,m),0,1))^2;

�� end

end

over;

dipole

% All check out to double precision float tolerances

z = linspace(-1/2,1/2,1001);

zeros = 0*z;

% plot overlaps with the dipole operator

for n = 0:5

�� for m = n:5

����� if (abs(dipole(m+1,n+1)) < 1e-16) dipole(m+1,n+1) = 0; end

����� y = sin(pi*(n+1)*(z+0.5)).*z.*sin(pi*(m+1)*(z+0.5));

����� ymax = max(abs(y));

����� plot(z,y,'k-', z,zeros,'k--')

� ����axis([-1/2,1/2,-ymax,ymax])

����� s = sprintf('Dipole operator: |<%g|x|%g>|^2 = %g',n,m,dipole(m+1,n+1));

����� title(s,'fontsize',14)

����� xlabel('x/L','fontsize',18)

����� pause

����� close

�� end

end

%z = linspace(-1/2,1/2,1000001);

%zeros = 0*z;

%n = 100; m = 101;

%dipole = (integral(@(z)fun(z,n,m),0,1))^2;

%����� y = sin(pi*(n+1)*(z+0.5)).*z.*sin(pi*(m+1)*(z+0.5));

%����� ymax = max(abs(y));

%����� plot(z,y,'k-', z,zeros,'k--')

%����� axis([-1/2,1/2,-ymax,ymax])

%����� s = sprintf('Dipole operator: |<%g|x|%g>|^2 = %g',n,m,dipole);

%����� title(s,'fontsize',14)

%����� xlabel('x/L','fontsize',18)

�����

��