2d double integral with trapz
4 次查看(过去 30 天)
显示 更早的评论
Hi: I am using Comsol to connect to matlab and run the following code. For some reason, the function trapz crashes. It was working before. Can any one help me with this?. what is wrong with this function ?. Was there any change in the use of this function in recent versions of matlab? please help. Naema
%Set parameters to compute fieber field.
lambda0 = 830e-9;
x0 = 0E-6; y0 = 0E-6; %x and y fiber offset
%%a - core radius;
%for SMF - 28: a = 4.1e-6; NA = .12 @ 1550nm;
a = 2e-6; NA = .12; V = 2*pi*a/lambda0*NA;
r_f = (5.3/2)*1E-6; %a*(.65 + 1.619/V^(3/2) + 2.879/V^6);
%Extract SPP field from Femlab.
fem.xmesh=meshextend(fem);
gridNumIntervals = 800;
Ymax = 5e-6;
Xmax = 5e-6;
[x, y] = meshgrid(-Xmax:2*Xmax/gridNumIntervals:Xmax, -Ymax:2*Ymax/gridNumIntervals:Ymax);
p=[x(:)'; y(:)'];
Espp=postinterp(fem,'Ey',p);
Espp=reshape(Espp, size(x));
Efib = exp(-((x-x0).^2 + (y-y0).^2)./r_f^2);
if 0
%Some times Femlab does not want to plot, so better save the fields and
%plot them in matlab.
figure(1);
subplot(211);
hold on;
surf(x*1e6, y*1e6, abs(Espp));
colorbar;
shading flat;
axis image;
xlabel('x(\mum)', 'fontsize', 12);
ylabel('y(\mum)', 'fontsize', 12);
% colormap(gray) %good colours 'bone' 'cool' 'winter'
subplot(212);
hold on;
surf(x*1e6, y*1e6, abs(Efib));
colorbar;
shading flat;
axis image;
xlabel('x(\mum)', 'fontsize', 12);
ylabel('y(\mum)', 'fontsize', 12);
%colormap(gray) %good colours 'bone' 'cool' 'winter'
end;
%Save the SPP and fiber fields and the coordinates.
% save('C:\Working folder\Theoretical\spp coupling eff\CEFF.mat','x','y','Espp','Efib');
save('C:\kholoud desk_2\Haneen\Israel prob\CEFF.mat','x','y','Espp','Efib');
%Compute coupling efficeincy.
fn=Efib.*conj(Espp);
fd=sqrt(trapz2D(x,y,Espp.*conj(Espp))*trapz2D(x,y,Efib.*conj(Efib)));
C = trapz2D(x,y,fn./fd);
absC = abs(C);
Ceff = absC;
CdB = -10*log10(Ceff);
beta0=2*pi/lambda0;
beta = sqrt(-fem.sol.lambda); %from the fem.sol.lambda
neff = beta/beta0; %beta*lambdaSweep(fInt)/2/pi;
MPA=20*log10(exp(1))*imag(beta/100)
%Display results.
disp('--------------------------------------------');
disp('');
disp(['Grid size in x: ', num2str(2*Xmax*1E6/gridNumIntervals),' um']);
disp(['Grid size in y: ', num2str(2*Ymax*1E6/gridNumIntervals),' um']);
disp('');
disp(['Fiber offset in x:', num2str(x0*1E6),' um']);
disp(['Fiber offset in y:', num2str(y0*1E6),' um']);
disp(['Complex coupling factor: ', num2str(C)]);
disp(['Coupling efficiency: ', num2str(Ceff)]);
disp(['Coupling loss per facet: ', num2str(CdB), ' dB']);
##########################################
here is the function associated with the code: % % A double integration of the tabulated data. % f(x,y) must be genarated in terms of meshgrid % Integral computes by approximation of the integral of f(x,y) via the trapezoidal method. % Firstly, an integration loop within inner integral over X is % performed obtaining the function of area A(y). Secondly, outer loop % integration is performed for A(y) over y limits % % Example: % Q = dblquad(inline('y*sin(x)+x*cos(y)'), pi, 2*pi, 0, pi); % Q = -9.8696 % [X,Y] = meshgrid(pi:pi/100:2*pi, 0:pi/100:pi ); % Z=Y.*sin(X)+X.*cos(Y); % int_2D_tabulated(X,Y, Z ) % ans = -9.8688
function V = trapz2D(X, Y, Z )
V = trapz( Y(:,1), (trapz(X(1,:),Z, 2)) ); return
0 个评论
回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!