clear all; clc; %close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% perfect!

lamda=300/183/1000;
%[a Rh]=meshgrid(7:.02:7.8,22:.5:35);
a=3.5/1000;%(7.8-lamda/4)/2;
wa=0.6435*a;         % wahich represents the beam raidius of apertue;
Rh=25.72/1000;
w0=wa./sqrt((1+(pi.*wa.^2/lamda./Rh).^2));
za=Rh./(1+(lamda.*Rh/pi./wa./wa).^2);
F=w0*pi/2/lamda;
fi0=atan(lamda.*za/pi./w0./w0);


zc1=pi.*w0.^2/lamda;     f=165/1000;  z01=(f+sqrt(f.*f-4.*zc1.*zc1))/2;  
z=z01;
wz=w0.*sqrt(1+(z./zc1).^2);     Te=(70/1000./wz/0.3393).^2

fprintf('\n The flare angle of corrugated horn is %g°. \n\n The focal ratio of this beam is %g.\n',asin(a/Rh)*180/pi,w0*pi/2/lamda);
fprintf('\n 近軸近似度  %g ',w0/lamda);
mode=4;
res=ones(mode,mode)*1000;
fi=zeros(mode,mode);

intv1=200;
intv2=360;
[r theta]=meshgrid(0:1*a/intv1:1*a,0:2*pi/intv2:2*pi);
x=r.*cos(theta);
y=r.*sin(theta);
x1=sqrt(2)*x/wa;
y1=sqrt(2)*y/wa;
p1=besselj(0,2.405/a*r);    %figure; surf(x,y,p1); title('bessel口徑場');
for d=0:size(res,1)-1
    for s=0:size(res,2)-1
        p2=sqrt(1/(pi*wa*wa*2^(d+s-1)*prod(1:d)*prod(1:s))) .*polyval(HermitePoly(d),x1) .*polyval(HermitePoly(s),y1) .*exp(-(x.^2+y.^2)/wa/wa);
        p4=p2.*p1.*r;
        res(d+1,s+1)=dblsim(r,theta,p4);;
        fi(d+1,s+1)=exp(-j*(d+s+1)*fi0) *exp(j*2*pi/lamda*za);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
    end
end
A=res.*fi;
% horn=load('data/para_A.mat');   A1=horn.A;
% surf(abs(A-A1));

save('data/para_A.mat','lamda','a','w0','za','F','A','mode');
% 
% sar2=p1.*p1.*r;
% bar2=dblsim(r,theta,sar2);
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % %this is another part
% % 
% 
% r0=a;
% [r theta]=meshgrid(0:r0/60:r0,0:2*pi/intv2:2*pi);
% x=r.*cos(theta);
% y=r.*sin(theta);
% x1=sqrt(2)*x/wa;
% y1=sqrt(2)*y/wa;
% sz=mode-1;
% sar1=zeros(size(x,1),size(x,2));
% for m=1:sz
%     for n=1:sz
%         sar1=sar1+res(m+1,n+1) .*  sqrt(1/(pi*wa*wa*2^(m+n-1)*prod(1:m)*prod(1:n))) .*Hermite(m,x1) .*Hermite(n,y1) .*exp(-(x.^2+y.^2)/wa/wa); 
%     end
% end
% figure;surf(x,y,abs(sar1)); title('高階模重建bessel口徑場');
% figure;surf(x,y,abs(sar1-p1));title('高階模對bessel口徑場的近似程度');
%         
%         
% sar1=sar1.^2.*r;        
% bar1=dblsim(r,theta,sar1);
% 
% fprintf('\n 高斯模主模效率為 %g.\n',bar1/bar2);