clear all;clc ;close all;warning off all
syms A0 Ac Fi0 Fi alfas Fs alfa r k z kz
%%%% Field of a plane wave after propagating a distance of z
UR = A0*exp(j*Fi0 + j*k*z); 
UR = exp(j*k*z + j*k*r^2/(2*z))/(4*pi*z);
UR = sin(Fi/2)*sin(alfa*r)/sqrt(alfa*r)*exp(-j*kz*z);k2 = alfa^2 + kz^2; 
%%% Derivatives of U(R) wrt r Fi and z
URr = diff(UR,'r');URrr = diff(URr,'r');URFi = diff(UR,'Fi');URFiFi = diff(URFi,'Fi');URz = diff(UR,'z');URzz = diff(URz,'z'); 

% %%% Alternative expressions for derivatives
%URr = diff(UR,'r');URrr = diff(UR,2,'r');URFiFi = diff(UR,2,'Fi');URz = diff(UR,'z');URzz = diff(UR,2,'z'); 
% %%% Analytic expression of Helmholtz equation  
%HE = URrr + URr/r + URFiFi/r^2 + 1*URzz + k^2*UR;simplify(HE)
HE = r^2*URrr + r*URr + URFiFi + r^2*URzz + r^2*k2*UR;simplify(HE)
%simplify(r^2*URrr + r*URr)
% % %%% Numeric evaluation
% lamda = 1.55e-6;k = 2*pi/lamda;r = 1e-2;alfa = 10;Fi = pi/2.3;z = 1e3;
% HE = eval(URrr + URr/r + URFiFi/r^2 + 1*URzz + k^2*UR)

% %%% The case of a spherical wave, since spherical wave representation is
% %%% an approximation, it does not satisfy Helmholtz equation
% UR = exp(j*k*z + j*k*r^2/(2*z))/(4*pi*z);
% URr = diff(UR,'r');URrr = diff(UR,2,'r');URFiFi = diff(UR,2,'Fi');URzz = diff(UR,2,'z');;
% % %%% Analytic expression 
%HE = URrr + URr/r + URFiFi/r^2 + URzz + k^2*UR;simplify(HE)
% 
% %%% The case of a Gaussian beam wave, this does not satisfy Helmholtz equation, but PWE
% alfa = 1/(k*alfas^2) + 0.5*j/Fs;
UR = Ac*exp(-k*alfa*r^2/(1 + 2*j*alfa*z))/(1 + 2*j*alfa*z);%UR = Ac*exp(-k*alfa*r^2);
URr = diff(UR,'r');URrr = diff(UR,2,'r');URFiFi = diff(UR,2,'Fi');URz = diff(UR,'z');
% % %% Analytic expression 
PWE = 1.0*URrr + URr/r + URFiFi/r^2 + 2*j*k*URz; %['PWE result of Gaussian beam = ']
simplify(PWE)