%% Computing products of Fourier-Bessel functions for triangles. clear all clf %% first set up angles, coordinates, orthocenter, sectors for the triangle. Check. A(1) = (1/4)*pi; A(2) = (1/4)*pi; A(3) = pi-A(1)-A(2); lambda=2;sqrl=sqrt(lambda); M=1; N=1;Nmax=5,Mmax=5% to be used for plotting %Triangle x(1) = 0; x(2) = 1; x(3) = 1/2*(1-(sin(A(1))/sin(A(3)))^2+(sin(A(2))/sin(A(3)))^2); y(1) = 0; y(2) = 0; y(3) =(sin(A(2))*sin(A(1))/sin(A(3))); [x' y'] figure(1) subplot(1,2,1) hold on patch(x,y,ones(size(x))); tex(x(2),y(2),'B','Color', 'b'); text(x(3),y(3),'C','Color', 'b'); % plots the triangle % Orthocenter calculations. First the triangle sides: c = 1, b = y(3)/sin(A(1)), a = y(3)/sin(A(2)), dd = a*cos(A(2))*cos(A(3))+b*cos(A(3))*cos(A(1))+c*cos(A(1))*cos(A(2)); ox(1) = (b*cos(A(3))*cos(A(1))+c*cos(A(1))*cos(A(2))*x(3))/dd; ox(2) = c*cos(A(1))*cos(A(2))*y(3)/dd line( [x(3),ox(1)],[y(3),ox(2)]) line( [ox(1),ox(1)],[0,ox(2)]) line( [0 ox(1)],[0 ox(2)]) plot(ox(1),ox(2),'r*'); text(ox(1),ox(2),'O'); % Ra = sqrt(ox(1)^2+ox(2)^2); Rb = sqrt((1-ox(1))^2+ox(2)^2); t=linspace(0,1,30); [xa,ya]=pol2cart(t*A(1),Ra); plot(xa,ya,'g'); [xb,yb]=pol2cart(pi-t*A(2),Rb); plot(1+xb,yb,'r'); title(['Triangle, sectors, orthocenter, for triangle with base angles alpha,beta = ',num2str(A(1)),',' num2str(A(2))]); axis equal subplot(1,2,2) %% defining basis functions % PhiA PhiA = @(x,y,n) gamma(n*pi/A(1)+1)*besselj(n*pi/A(1), sqrl*sqrt(x.^2+y.^2)).*cos(n*pi*atan2(y, x)./A(1)); [th,rh]=meshgrid(t*A(1),t*Ra);[Xa,Ya]=pol2cart(th,rh);Za= PhiA(Xa,Ya,N); hold on patch(x,y,ones(size(x))); line( [x(3),ox(1)],[y(3),ox(2)]);line( [ox(1),ox(1)],[0,ox(2)]); line( [0 ox(1)],[0 ox(2)]) plot(ox(1),ox(2),'r*'); surf(Xa,Ya,Za);view(0,90),shading flat axis equal % PhiB PhiB = @( x, y, n) gamma(n*pi/A(2)+1)*besselj(n*pi/A(2), sqrl*sqrt((x-1).^2+y.^2)).*cos(n*pi*atan2(y, 1-x)./A(2)); [th,rh]=meshgrid(t*A(2),t*Rb);[Xb,Yb]=pol2cart(th,rh);Xb=1-Xb;Zb= PhiB(Xb,Yb,M); surf(Xb,Yb,Zb);view(0,90),colorbar, shading flat title(['Basis functions in sectors, n,m=',num2str(N),',',num2str(M)]) axis([-0.1 1.1 0 1.1]) %% product of basis functions ProdPhiAB=@(x,y,n,m)PhiA(x,y,n).*PhiB(x,y,m); RestrictAB1=@(x,y,n,m)(x>ox(1)).*(x.^2+y.^2ox(1)).*(x.^2+y.^2