clear all;clc ;close all;warning off all
syms a z k r Fir n aB AL Ds as as1 as2
%%%% Field of a Gaussian beam wave after propagating a distance of z, excluding exp(jkz)
% Ur = 1/(1 + 2*j*a*z)*exp(-k*a*r^2/(1 + 2*j*a*z));   
% %%% derivatives of U_r wrt r (first and second orders) and wrt z
% Ur1 = diff(Ur,'r');Ur2 = diff(Ur1,'r');Urz = diff(Ur,'z'); 
% Ur1 = diff(Ur,'r');Ur2 = diff(Ur,'r',2);Urz = diff(Ur,'z'); %%% Alternative second derivative
% Urz2 = diff(Ur,'z',2); %%% Extra term form Helmholtz Eq
% PWE = Ur2 + Ur1/r + 2*j*k*Urz;simplify(PWE)
% HE = Ur2 + Ur1/r + Urz2 + 2*j*k*Urz;simplify(HE)
% lamda = 1.55e-6;k = 2*pi/lamda;alfas = 1e-2;a = 1/(k*alfas^2);  %%%% Inserting numeric values
% z = 1e3;r = 1e-2;n = 3;Fir = pi/4;aB = 100;%%%% Inserting numeric values
% PWE = eval(Ur2 + Ur1/r + 2*j*k*Urz) %%%%% Testing PWE numerically
% %Ur2s = num2str(eval(Ur2));
% ['Ur2 = ' num2str(eval(Ur2))]
% ['Ur1/r = ' num2str(eval(Ur1/r))]
% ['Urz2 = ' num2str(eval(Urz2))]
% ['2jkUrz = ' num2str(eval(2*j*k*Urz))]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Bessel Gaussian beam
Ur = -exp(-j*n*Fir)/(1 + 2*j*a*z)*exp(-(j*aB^2*z + 2*a*k^2*r^2)/(2*k*(1 + 2*j*a*z)))*besselj(n,aB*r/(1 + 2*j*a*z)); 
Ur1 = diff(Ur,'r');Ur2 = diff(Ur1,'r');UrFi2 = diff(Ur,'Fir',2);Urz = diff(Ur,'z');
PWE = Ur1 + Ur2/r + UrFi2/r^2 + 2*j*k*Urz;
simplify(PWE)
lamda = 1.55e-6;k = 2*pi/lamda;alfas = 1e-2;a = 1/(k*alfas^2);
z = 1e3;r = 1e-2;n = 6;Fir = pi/5;aB = 150;
PWE = eval(Ur2 + Ur1/r + UrFi2/r^2 + 2*j*k*Urz)
%%%% Case of Sino / Hyperbolic beam Q1 of ECE MT 22.11.2011
% ALarr = [AL 0];Dsarr = [Ds -Ds];asarr = [as as];aLarr = 1./(k*asarr.^2);Ur = 0; %%% For Gaussian beam
% % ALarr = [AL AL];Dsarr = [Ds -Ds];asarr = [as as];aLarr = 1./(k*asarr.^2);Ur = 0; %%% For cosh Gaussian beam
% ALarr = [AL AL];Dsarr = [j*Ds -j*Ds];asarr = [as as];aLarr = 1./(k*asarr.^2);Ur = 0; %%% For cos Gaussian beam
% ALarr = [AL -AL];Dsarr = [Ds -Ds];asarr = [as as];aLarr = 1./(k*asarr.^2);Ur = 0; %%% For sinh Gaussian beam
% ALarr = [j*AL -j*AL];Dsarr = [j*Ds -j*Ds];asarr = [as as];aLarr = 1./(k*asarr.^2);Ur = 0; %%% For sine Gaussian beam
% ALarr = [AL -0.5*AL];Dsarr = [0*Ds -0*Ds];asarr = [as1 as2];aLarr = 1./(k*asarr.^2);Ur = 0; %%% For annular Gaussian beam
% for L = 1:2
% Al = ALarr(L);al = aLarr(L);Dl = Dsarr(L);
% Ur = Ur + Al/(1 + 2*j*al*z)*exp(-(k*al*r^2 - j*(cos(Fir) + sin(Fir))*r*Dl + j/k*Dl^2*z)/(1 + 2*j*al*z));
% end
% Urr = diff(Ur,'r');Urrr = diff(Urr,'r');UrFi2 = diff(Ur,'Fir',2);Urz = diff(Ur,'z');
% PWE = Urrr + Urr/r + UrFi2/r^2 + 2*j*k*Urz;simplify(PWE)