How to perform integral ?

clear all;

clc;

close all;

dim=64;

L=531*10^-9;

k=2*pi/L;

[x,y]=meshgrid(-L:2*L/(dim-1):L);

r=sqrt(x.^2+y.^2);

alp=asin(0.95);

h=(alp)/64;

z=-32*L:1*L:31*L;

c=numel(z);

Era=zeros([dim dim c]);

Ela=zeros([dim dim c]);

Er=0;

Ez=0;

Hp=0;

for ii=1:c;

   % electric field components

f=@(theta) -(cos(theta))^0.5*sin(2*theta)*bessel(1,k*r*sin(theta));

 g=@(theta) -2*1i*(cos(theta))^0.5*(sin(theta))^2*bessel(0,k*r*sin(theta));%*exp(1i*k*z(ii)*cos(theta));

for jj=1:1:63

theta=(jj*h); Er=Er+(f(0)+2*f(jj*h)+f(alp));

Era(:,:,ii)=Er.*exp(1i*k*z(ii)*cos(theta));

end

 for jj=1:1:63

theta=(jj*h);

    Ez=Ez+(g(0)+2*g(jj*h)+g(alp));

Ela(:,:,ii)=Ez.*exp(1i*k*z(ii)*cos(theta));

 end

end

Er1=(h/2)*(Era);

 Ez1=(h/2)*(Ela);

Ir=Er1.*conj(Er1);

Iz=Ez1.*conj(Ez1);

Itot=Ir+Iz;

Itot=squeeze(Itot(dim/2,:,:)).';

figure(1)

imagesc(Itot),colormap gray

Hi,

By using this Program I am calculating intensity of radially polarized beam

in YZ plane.

By using this program for z=0, that is in XY plane I am getting correct

 intensity pattern, but when I am varying z from -32:1:31 in YZ plane I not

 getting correct intensity pattern that I am getting half part only other

half part I am missing from figure(1);

Here I am performing Integration of following function by using trepazoidal

rule

f=@(theta)

-(cos(theta))^0.5*sin(2*theta)*bessel(1,k*r*sin(theta))*exp(1i*k*z(ii)*cos(theta));

 g=@(theta) -2*1i*(cos(theta))^0.5*

(sin(theta))^2*bessel(0,k*r*sin(theta))*exp(1i*k*z(ii)*cos(theta));

I have problem in Integration of this to functions

we will get intensity as a function of r and z

can any one help me in this regard,it is great help for me

thank you