% This is Genbeam intensities using ES ExpPolyHerm4 quadrapule routine - 13.07.2009
% Note that cartesian coordinates are used
function Genbeam_rec_Lag
clear;clc;clf;warning off MATLAB:divideByZero;close all
lamda = 1.55e-6; k = 2*pi/lamda;Cn = 1e-15;sig2 = 2e+2;
narr = [1 0 3 0] ;narr = 0;Vxmat = [0 4.7 25 100 200 250]; 
sxset = -1.e-1:0.2e-2:1e-1;syset = sxset;
%sxset = sxset*0.8;syset = sxset; 
sxset = sxset*0.3;syset = sxset; 
%%%%%%%%%%%%%%%%%%%%%%%%%%      Genbeam    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sx,sy]=meshgrid(sxset,syset);tic
%%%%% Source settings %%%%%%%%%%

alfas = 1e-2;alfa = 1/(k*alfas^2);alfasx = alfas/sqrt(1);alfasy = alfas/sqrt(1);  
alfasarr = [1e-2 1e-2 1e-2 1e-2];alfasarr = alfasarr*1;Larr = 0.01e3:0.8e3:3.6e3;
ALarr = [0.5 0.5];narr = [0 0];marr = [0 0];
alfasxarr = [1e-2 1e-2];%alfasxarr = alfasxarr*3; 
alfasyarr = [1e-2 1e-2];%alfasyarr = alfasyarr*3; 
Vxarr = [1*Vxmat(1) -1*Vxmat(1)];Vyarr = Vxarr;Fsxarr = [5e15 5e15];Fsyarr = [5e15 5e15];
axarr =1./alfasxarr;ayarr = 1.00./alfasyarr;
alfaxarr = 1./(k*alfasxarr.^2) + j./Fsxarr;alfayarr = 1./(k*alfasyarr.^2) + j./Fsyarr;
%%%%%%%%%%%%% Source Beam %%%%%%%%%%%%%%%%%%
I_sarr = [];maxgenarr = [];u_s = zeros(size(sx));
for is = 1:length(narr)
    AL = ALarr(is);n = narr(is);m = marr(is);ax = axarr(is);ay = ayarr(is);
    alfax = alfaxarr(is);alfay = alfayarr(is);Vx = Vxarr(is);Vy = Vyarr(is);
%%% Source in cartesian coordinates
u_sarr = AL*Her(n,ax*sx).*Her(m,ay*sy).*exp(-(0.5*k*alfax.*sx.^2 +j*Vx*sx)).*exp(-(0.5*k*alfay.*sy.^2 +j*Vy*sy));
u_s = u_s + u_sarr;end 
I_s = u_s.*conj(u_s);maxgen = max(max(I_s));I_s = I_s / max(max(I_s));maxgenarr = [maxgenarr maxgen];
sxn = sx / 1e-2;syn = sy / 1e-2;I_sarr = cat(3,I_sarr,I_s); 
N = [0 0 0 0];LG = [0 0 0 0;0 0 0 0];ib = 1;
R = sqrt(sxset.^2+sxset.^2);RS = [-R(1:round(numel(R)/2)) R(round(numel(R)/2)+1:end)];RSn = RS / 1e-2;
rxarr = 2*sxset;Larr = 1.0e3;ry = 0;I_rmat = [];
R = sqrt(rxarr.^2+rxarr.^2);RP = [-R(1:round(numel(R)/2)) R(round(numel(R)/2)+1:end)];RPn = RP / 1e-2;
[rxm,rym] = meshgrid(rxarr,rxarr);
%%%%%%%%%%%%% Receiver Beam - Finding <I(r,L)> %%%%%%%%%%%%%%%%%
for iL = 1:length(Larr)
    L = Larr(iL);I_rarr = [];b = k/(2*L);
    ro2 = (0.546*Cn*k^2*L)^(-1.2);
for ir = 1:length(RP)
    rx = rxarr(ir);ry = rxarr(ir);I_r = 0;%[iL ir]
for L1 = 1:length(narr);for L2 = 1:length(narr);AL1 = ALarr(L1);AL2 = ALarr(L2);
    n1 = narr(L1);n2 = narr(L2);m1 = marr(L1);m2 = marr(L2);axL1 = axarr(L1);axL2 = axarr(L2);ayL1 = ayarr(L1);ayL2 = ayarr(L2);
    alfasxL1 = alfasxarr(L1);alfasxL2 = alfasxarr(L2);alfasyL1 = alfasyarr(L1);alfasyL2 = alfasyarr(L2);
    VxL1 = Vxarr(L1);VxL2 = Vxarr(L2);VyL1 = Vyarr(L1);VyL2 = Vyarr(L2);
    q11 = 0.5/alfasxL1^2 + 1/ro2 + 0.5/sig2 - j*b;q22 = 0.5/alfasxL2^2 + 1/ro2 + 0.5/sig2 + j*b;
    q33 = 0.5/alfasyL1^2 + 1/ro2 + 0.5/sig2 - j*b;q44 = 0.5/alfasyL2^2 + 1/ro2 + 0.5/sig2 + j*b;
    q12 = -1/ro2 - 0.5/sig2;q34 = q12;q13 = 0;q14 = 0;q23 = 0;q24 = 0;
    Q = [q11 q12 q13 q14;q12 q22 q23 q24;q13 q23 q33 q34;q14 q24 q34 q44];
    R = [j*b*rx+0.5*j*VxL1 -j*b*rx-0.5*j*conj(VxL2) j*b*ry+0.5*j*VyL1 -j*b*ry-0.5*j*conj(VyL2)];
    H = [n1 n2 m1 m2;axL1 axL2 ayL1 ayL2];
    I_r = I_r + AL1*conj(AL2)*ExpPolyHermLagu401(Q,R,N,H,LG);
    end;end;
I_rarr = [I_rarr I_r];end
I_rmat = [I_rmat;I_rarr];end
I_rmat = I_rmat*k^2/(4*pi^2*L^2)/maxgenarr(ib);%I_rmat
figure(2)
I_s = I_sarr(:,:,ib);plot(RSn,diag(I_s),'-k','LineWidth',2);%I_s
axis ([min(RSn)*1.2 max(RSn)*1.2 0 1.1])
set(gcf,'Color',[1 1 1]);set(gca,'FontSize',16)
xlabel('\its_x\rm\bf  axis in cm','FontSize',14,'FontWeight','bold');
zlabel('\fontname{Verdana}\fontsize{12}\it I_s_N ','FontSize',12,'FontWeight','bold');grid on;
figure(3);%I_rmat
I_s = I_rmat;plot(RPn,I_s,'--k','LineWidth',3)
set(gcf,'Color',[1 1 1]);set(gca,'FontSize',16)
xlabel('\its_x\rm\bf  axis in cm','FontSize',14,'FontWeight','bold');
zlabel('\fontname{Verdana}\fontsize{12}\it I_s_N ','FontSize',12,'FontWeight','bold');grid on;
axis ([min(RPn)*1 max(RPn)*1 0 max(max(I_s))]);toc;I_s = I_s';
save('C:\Documents and Settings\HTERemote\My Documents\MATLAB\Trans\Data\Cosh_fromlagu4','I_s','-ASCII')
%%%%%%%%%%%%%%%%% Two Intensity Profiles at Source Plane %%%%%%%%%%%
function yout = Her(nn,xx)
mat = cat(4,ones(size(xx)),2*xx,4*xx.^2-2);
for say = 4:nn+1
    mat(:,:,:,say) = 2*xx.*mat(:,:,:,say-1) - 2*(say-2)*mat(:,:,:,say-2);
end
yout = mat(:,:,:,nn+1);