how plot amplitude of mode from eigenvector?
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Hi every body.
For solving an integral (this integral is Huygens Fresnel integral and it is useful for find modes of a laser resonator), I used Gauss-Legendre quadrature method and in final I got its eigenvectors and eigenvalues. but I couldn't plot the mode's amplitude and phase. If anyone has a experience about it, I hope help?
these codes are what is that I wrote.
clc
clear all
close all
a=0;
b=1;
N=200;
N=N-1;
N1=N+1; N2=N+2;
xu=linspace(a,b,N1)';
% Initial guess
y=cos((2*(0:N)'+1)*pi/(2*N+2))+(0.27/N1)*sin(pi*xu*N/N2);
% Legendre-Gauss Vandermonde Matrix
L=zeros(N1,N2);
% Derivative of LGVM
Lp=zeros(N1,N2);
% Compute the zeros of the N+1 Legendre Polynomial
% using the recursion relation and the Newton-Raphson method
y0=2;
% Iterate until new points are uniformly within epsilon of old points
while max(abs(y-y0))>eps
L(:,1)=1;
Lp(:,1)=0;
L(:,2)=y;
Lp(:,2)=1;
for k=2:N1
L(:,k+1)=( (2*k-1)*y.*L(:,k)-(k-1)*L(:,k-1) )/k;
end
Lp=(N2)*( L(:,N1)-y.*L(:,N2) )./(1-y.^2);
y0=y;
y=y0-L(:,N2)./Lp;
end
% Linear map from[-1,1] to [a,b]
x=(a*(1-y)+b*(1+y))/2;
% Compute the weights
w=(b-a)./((1-y.^2).*Lp.^2)*(N2/N1)^2;
LL=40.92;
A=1; B=LL; C=0;
n=N+1; m=n; s=n;
landa=10.6e-4;
k=2*pi/landa;
r1=linspace(a,b,n);
r2=linspace(a,b,m);
kernel1=zeros(n,m);
kernel2=zeros(n,m);
H1=zeros(n,m);
H2=zeros(n,m);
T1=zeros(1,m);
T2=zeros(1,m);
R=100*LL;
teta0=12.22e-3*pi/360;
T1(:)=exp(-2*1i*k*teta0*r1(:));
T2(:)=exp(-1i*k*r1(:).^2./R);
for s=1:n
kernel1(:,s)=-1i*(k/B)*r2(s).*besselj(0,r2(s).*r1(:)*k/B).*exp(1i*k*(A*r2(s).^2+B*r1(:).^2)/B/2);
H1(:,s)=kernel1(:,s).*T1(:);
end
for s=1:m
kernel2(:,s)=-1i*(k/B)*r1(s).*besselj(0,r1(s).*r2(:)*k/B).*exp(1i*k*(A*r1(s).^2+B*r2(:).^2)/B/2);
H2(:,s)=kernel2(:,s).*T2(:);
end
W=diag(w);
P1=H1*W;
P2=H2*W;
HH=P1*P2;
[V,D]=eig(HH);
I want to plot below pictures:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/152090/image.jpeg)
0 Kommentare
Antworten (0)
Siehe auch
Kategorien
Mehr zu Interpolation finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!