% This is Transmittance formulation of mine 24.12.2006 
% Note that the notation sigma is used for rho_s
% Note that this programme generates intensity as well as tau (transmission coefficient)
% This way it can also be used for partially coherent general beam propagation
% Note that free space intensity is calculated concurrently by inserting Cn2 = 1e-40
% Errored terms are corrected on 9.03.2007
% This is the same as Transmittance, but modified for non-Kol paper 2011
function transmit635
clear;clc;clf reset;warning off MATLAB:divideByZero;close all
lamda = 1.55e-6;k = 2*pi/lamda;sig2 = 1e+10;L = 2e3;
Cn = 1e-15;Cn_fs = 1e-40;Vxmat = [0 4.7 25 100 200 250];
%%%%%%%%%%%% Source Parameters %%%%%%%%%%%%%%%%
ALarr = [1 0;1 1;1 1;1 -1];tetaarr = [0 0;0 0;0 0;0 0];
narr = [0 0;0 0;0 0;0 0];marr = [0 0;0 0;0 0;0 0];
alfasxarr = [1e-2 1e-2;1e-2 1e-2;1e-2 1e-2;1e-2 0.8e-2];alfasxarr = alfasxarr/sqrt(2); 
alfasyarr = [1e-2 1e-2;1e-2 1e-2;1e-2 1e-2;1e-2 0.8e-2];alfasyarr = alfasyarr/sqrt(2); 
Vxarr = [1*Vxmat(1) -1*Vxmat(1);1*Vxmat(4) -1*Vxmat(4);j*Vxmat(4) -j*Vxmat(4);
         j*Vxmat(1) -j*Vxmat(1)];Vyarr = Vxarr;
Fsxarr = [5e15 5e15;5e15 5e15;5e15 5e15;5e15 5e15];
Fsyarr = [5e15 5e15;5e15 5e15;5e15 5e15;5e15 5e15]; 

sxsetn = -1.e-1:0.2e-2:1e-1;sysetn = sxsetn;

sxset = sxsetn*0.3;syset = sysetn*0.3; 

%sxset = sxsetn*0.2;syset = sysetn*0.2;
%sxset = sxsetn*0.4;syset = sysetn*0.4; 
p1x = 0;p2x = 0;p1y = 0;p2y = 0;
% -------- Formulation ------------
axarr =1./alfasxarr;ayarr = 1.00./alfasyarr;bxarr = zeros(size(axarr)); byarr = zeros(size(ayarr));
alfaxarr = 1./(k*alfasxarr.^2) + j./Fsxarr;alfayarr = 1./(k*alfasyarr.^2) + j./Fsyarr;
Cstr = [' \itC_n^2 = \rm\bf1x10^ ^-^1^5 m^-^2^/^3'];
lamstr = [' \it\lambda = \rm\bf' num2str(lamda/1e-6) ' \mum'];
Lstr = [' \itL = \rm\bf' num2str(L/1e3) ' km'];
sigstr = [' \fontsize {11}\it\rho_s \rm\bf\fontsize {10}= ' num2str(sqrt(sig2)/1e-3,'%2.0f') ' mm '];
legc = {Cstr;' ';lamstr;' ';Lstr};
legtitle = [' Gaussian beam             ';
            ' cos-Gaussian beam         ';
            ' cosh-Gaussian beam        ';
            ' annular-Gaussian beam     '];
Vxstr = [' \itD_s = \rm\bf' num2str(abs(Vxarr(2,:))) ' m^-^1 '];
Vystr = [' \itD_s = \rm\bf' num2str(abs(Vyarr(2,:))) ' m^-^1 '];
alfastrsx = [' \it\alpha_s = \rm\bf' num2str(alfasxarr(1,:)*sqrt(2)/1e-2,'%4.1f') ' cm '];
alfastrsy = [' \it\alpha_s_y = \rm\bf' num2str(alfasyarr(1,:)/1e-2,'%4.1f') ' cm '];

%%%%%%%%%%%%% Beam at Source Plane %%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% Note this is used for establishing equvalence only %%%%%%
[sx,sy]=meshgrid(sxset,syset);sxn = sx / 1e-2;syn = sy / 1e-2;tic
I_sarr = [];maxgenarr = [];
for ru = 1:length(ALarr(:,1))
%for ru = 1:1    
u_s = zeros(size(sx));
for Ls = 1:length(ALarr(1,:))
    AL = ALarr(ru,Ls);teta = tetaarr(ru,Ls);n = narr(ru,Ls);m = marr(ru,Ls);
    ax = axarr(ru,Ls);bx = bxarr(ru,Ls);ay = ayarr(ru,Ls);by = byarr(ru,Ls);
    alfax = alfaxarr(ru,Ls);alfay = alfayarr(ru,Ls);Vx = Vxarr(ru,Ls);Vy = Vyarr(ru,Ls);
    u_sarr = AL.*exp(-j*teta).*H(n,ax*sx+bx).*H(m,ay*sy+by).* .... 
    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));%maxgen = 1; %%% This setting is used for area scintillation
I_s = I_s / maxgen;maxgenarr = [maxgenarr maxgen];
I_sarr = cat(3,I_sarr,I_s);end;
%iatla = 1;if iatla == 2  % For skipping receiver plane plots
%%%%%%%%%%%%%%%%% Intensity Profile at Receiver Plane %%%%%%%%%%%
% Beam after Propogation
pxsetarr = [5*sxset;5*sxset;5*sxset;5*sxset];pxsetarr = 1*pxsetarr;pysetarr = pxsetarr;
%pxsetarr = [1.0*sxset;1.8*sxset;4*sxset;60*sxset];pysetarr = pxsetarr;
%pxsetarr = [1.0*sxset;1*sxset;1*sxset;2*sxset];pysetarr = pxsetarr;
%pxsetarr = [1.8*sxset;1.8*sxset;1.8*sxset];pysetarr = pxsetarr;
I_arr = [];I_arr_fs = [];tet_arr = ones(size(sx));
%%%%% Loop for beam type %%%%%%
for ru = 1:length(ALarr(:,1))
%for ru = 1:1  
    
%    I_arr = cat(3,I_arr,I_sarr(:,:,ru));I_arr_fs = cat(3,I_arr_fs,I_sarr(:,:,ru)); 
%%%%%% Loop for propagation lengths %%%%%%

    pxset = pxsetarr(ru,:);pyset = pysetarr(ru,:);
    [px,py]=meshgrid(pxset,pyset);pxn = px / 1e-2;pyn = py / 1e-2;
    ro2 = (0.546*Cn*k^2*L)^(-1.2);b = k/(2*L);ro2_fs = (0.546*Cn_fs*k^2*L)^(-1.2);
    I_r = zeros(size(sx));I_r_fs = zeros(size(sx));

%%%%% Loop for summation in each beam type %%%%%%    
for L1 = 1:length(ALarr(1,:))
for L2 = 1:length(ALarr(1,:)) 
AL1 = ALarr(ru,L1);AL2 = ALarr(ru,L2);tetaL1 = tetaarr(ru,L1);tetaL2 = tetaarr(ru,L2); 
nL1 = narr(ru,L1);mL1 = marr(ru,L1);nL2 = narr(ru,L2);mL2 = marr(ru,L2);
alfaxL1 = alfaxarr(ru,L1);alfaxL2 = alfaxarr(ru,L2);alfayL1 = alfayarr(ru,L1);alfayL2 = alfayarr(ru,L2);
VxL1 = Vxarr(ru,L1);VxL2 = Vxarr(ru,L2);VyL1 = Vyarr(ru,L1);VyL2 = Vyarr(ru,L2);
axL1 = axarr(ru,L1);axL2 = axarr(ru,L2);bxL1 = bxarr(ru,L1);bxL2 = bxarr(ru,L2);
ayL1 = ayarr(ru,L1);ayL2 = ayarr(ru,L2);byL1 = byarr(ru,L1);byL2 = byarr(ru,L2);

%%%%%% Loop for transmittance ratio %%%%%%
Dx = 0.5*k*alfaxL1*sig2*ro2 + sig2 + 0.25*ro2 -j*b*sig2*ro2;
Dx_fs = 0.5*k*alfaxL1*sig2*ro2_fs + sig2 + 0.25*ro2_fs -j*b*sig2*ro2_fs;
Dy = 0.5*k*alfayL1*sig2*ro2 + sig2 + 0.25*ro2 -j*b*sig2*ro2;
Dy_fs = 0.5*k*alfayL1*sig2*ro2_fs + sig2 + 0.25*ro2_fs -j*b*sig2*ro2_fs;
Px = 0.5*k*(0.5*k*alfaxL1*conj(alfaxL2)*sig2*ro2 + (alfaxL1 + conj(alfaxL2))*(sig2 + 0.25*ro2) + ....
    j*b*sig2*ro2*(alfaxL1 - conj(alfaxL2)) + 2*b^2*sig2*ro2/k);
Px_fs = 0.5*k*(0.5*k*alfaxL1*conj(alfaxL2)*sig2*ro2_fs + (alfaxL1 + conj(alfaxL2))*(sig2 + 0.25*ro2_fs) + ....
    j*b*sig2*ro2_fs*(alfaxL1 - conj(alfaxL2)) + 2*b^2*sig2*ro2_fs/k);
Py = 0.5*k*(0.5*k*alfayL1*conj(alfayL2)*sig2*ro2 + (alfayL1 + conj(alfayL2))*(sig2 + 0.25*ro2) + ....
    j*b*sig2*ro2*(alfayL1 - conj(alfayL2)) + 2*b^2*sig2*ro2/k);
Py_fs = 0.5*k*(0.5*k*alfayL1*conj(alfayL2)*sig2*ro2_fs + (alfayL1 + conj(alfayL2))*(sig2 + 0.25*ro2_fs) + ....
    j*b*sig2*ro2_fs*(alfayL1 - conj(alfayL2)) + 2*b^2*sig2*ro2_fs/k);
Qx = (p1x - p2x)*( - 3*j*b*sig2*ro2 - 0.5*j*b*ro2^2 + 0.5*k*alfaxL1*sig2*ro2) + ....
    j*b*sig2*ro2^2*(k*alfaxL1 - 2*j*b)*px + j*ro2*conj(VxL2)*Dx -j*ro2*VxL1*(sig2 + 0.25*ro2);
Qx_fs = (p1x - p2x)*( - 3*j*b*sig2*ro2_fs - 0.5*j*b*ro2_fs^2 + 0.5*k*alfaxL1*sig2*ro2_fs) + ....
    j*b*sig2*ro2_fs^2*(k*alfaxL1 - 2*j*b)*px + j*ro2_fs*conj(VxL2)*Dx_fs -j*ro2*VxL1*(sig2 + 0.25*ro2_fs);
Qy = (p1y - p2y)*( - 3*j*b*sig2*ro2 - 0.5*j*b*ro2^2 + 0.5*k*alfayL1*sig2*ro2) + ....
    j*b*sig2*ro2^2*(k*alfayL1 - 2*j*b)*py + j*ro2*conj(VyL2)*Dy -j*ro2*VyL1*(sig2 + 0.25*ro2);
Qy_fs = (p1y - p2y)*( - 3*j*b*sig2*ro2_fs - 0.5*j*b*ro2_fs^2 + 0.5*k*alfayL1*sig2*ro2_fs) + ....
    j*b*sig2*ro2_fs^2*(k*alfayL1 - 2*j*b)*py + j*ro2_fs*conj(VyL2)*Dy_fs -j*ro2*VyL1*(sig2 + 0.25*ro2_fs);

Ex = (sig2*ro2*Px*(-p1x + p2x - 2*j*b*ro2*px - j*ro2*VxL1).^2 + Qx.^2)/(4*ro2^2*Dx*Px);
Ey = (sig2*ro2*Py*(-p1y + p2y - 2*j*b*ro2*py - j*ro2*VyL1).^2 + Qy.^2)/(4*ro2^2*Dy*Py);
Ex_fs = (sig2*ro2_fs*Px_fs*(-p1x + p2x - 2*j*b*ro2_fs*px - j*ro2_fs*VxL1).^2 + Qx_fs.^2)/(4*ro2_fs^2*Dx_fs*Px_fs);
Ey_fs = (sig2*ro2_fs*Py_fs*(-p1y + p2y - 2*j*b*ro2_fs*py - j*ro2_fs*VyL1).^2 + Qy_fs.^2)/(4*ro2_fs^2*Dy_fs*Py_fs);

sumx = 0;sumy =0;sumx_fs = 0;sumy_fs =0;
for Lx1 = 0:fix(nL1/2);for Lnx1 = 0:nL1-2*Lx1;for kx1 = 0:fix(Lnx1/2);for Lnx11 = 0:Lnx1-2*kx1 
for Lx2 = 0:fix(nL2/2);for Lnx2 = 0:nL2-2*Lx2;for kx2 = 0:fix((Lnx11+Lnx2)/2)    
    TLx1 = funT(Lx1);TLx2 = funT(Lx2);facLx1 = funb(nL1,2*Lx1);facLx2 = funb(nL2,2*Lx2);facLnx1 = funb(nL1-2*Lx1,Lnx1);
    facLnx11 = funb(Lnx1-2*kx1,Lnx11);facLnx2 = funb(nL2-2*Lx2,Lnx2);
    
    sumx = sumx + (-1)^(Lx1+Lx2)*2^(nL1+nL2-Lx1-Lx2-Lnx1-Lnx2)*TLx1*TLx2*facLx1*facLx2*facLnx1*facLnx11*facLnx2* ....
    factorial(Lnx1)/(factorial(Lnx1-2*kx1)*factorial(kx1))*factorial(Lnx11+Lnx2)/(factorial(Lnx11+Lnx2-2*kx2)*factorial(kx2))* ....
    axL1^Lnx1*(conj(axL2))^Lnx2*bxL1^(nL1-2*Lx1-Lnx1)*conj(bxL2)^(nL2-2*Lx2-Lnx2)*sig2^(Lnx1-kx1-Lnx11)*Dx^(-Lnx1+kx1+kx2)* ....
    (-p1x+p2x-2*j*b*ro2*px-j*ro2*VxL1).^(Lnx1-2*kx1-Lnx11)*(sig2+0.25*ro2)^Lnx11*ro2^(kx1-Lnx11-Lnx2+2*kx2)* ....
    Px^(-Lnx11-Lnx2+kx2).*Qx.^(Lnx11+Lnx2-2*kx2);

    sumx_fs = sumx_fs + (-1)^(Lx1+Lx2)*2^(nL1+nL2-Lx1-Lx2-Lnx1-Lnx2)*TLx1*TLx2*facLx1*facLx2*facLnx1*facLnx11*facLnx2* ....
    factorial(Lnx1)/(factorial(Lnx1-2*kx1)*factorial(kx1))*factorial(Lnx11+Lnx2)/(factorial(Lnx11+Lnx2-2*kx2)*factorial(kx2))* ....
    axL1^Lnx1*(conj(axL2))^Lnx2*bxL1^(nL1-2*Lx1-Lnx1)*conj(bxL2)^(nL2-2*Lx2-Lnx2)*sig2^(Lnx1-kx1-Lnx11)*Dx_fs^(-Lnx1+kx1+kx2)* ....
    (-p1x+p2x-2*j*b*ro2_fs*px-j*ro2_fs*VxL1).^(Lnx1-2*kx1-Lnx11)*(sig2+0.25*ro2_fs)^Lnx11*ro2_fs^(kx1-Lnx11-Lnx2+2*kx2)* ....
    Px_fs^(-Lnx11-Lnx2+kx2).*Qx_fs.^(Lnx11+Lnx2-2*kx2);
end;end;end;end;end;end;end

for Ly1 = 0:fix(mL1/2);for Lmy1 = 0:mL1-2*Ly1;for ky1 = 0:fix(Lmy1/2);for Lmy11 = 0:Lmy1-2*ky1 
for Ly2 = 0:fix(mL2/2);for Lmy2 = 0:mL2-2*Ly2;for ky2 = 0:fix((Lmy11+Lmy2)/2)      
    TLy1 = funT(Ly1);TLy2 = funT(Ly2);facLy1 = funb(mL1,2*Ly1);facLy2 = funb(mL2,2*Ly2);facLmy1 = funb(mL1-2*Ly1,Lmy1);
    facLmy11 = funb(Lmy1-2*ky1,Lmy11);facLmy2 = funb(mL2-2*Ly2,Lmy2);
    
    sumy = sumy + (-1)^(Ly1+Ly2)*2^(mL1+mL2-Ly1-Ly2-Lmy1-Lmy2)*TLy1*TLy2*facLy1*facLy2*facLmy1*facLmy11*facLmy2* ....
    factorial(Lmy1)/(factorial(Lmy1-2*ky1)*factorial(ky1))*factorial(Lmy11+Lmy2)/(factorial(Lmy11+Lmy2-2*ky2)*factorial(ky2))* ....
    ayL1^Lmy1*(conj(ayL2))^Lmy2*byL1^(mL1-2*Ly1-Lmy1)*conj(byL1)^(mL2-2*Ly2-Lmy2)*sig2^(Lmy1-ky1-Lmy11)*Dy^(-Lmy1+ky1+ky2)* ....
    (-p1y+p2y-2*j*b*ro2*py-j*ro2*VyL1).^(Lmy1-2*ky1-Lmy11)*(sig2+0.25*ro2)^Lmy11*ro2^(ky1-Lmy11-Lmy2+2*ky2)* ....
    Py^(-Lmy11-Lmy2+ky2).*Qy.^(Lmy11+Lmy2-2*ky2);
    
    sumy_fs = sumy_fs + (-1)^(Ly1+Ly2)*2^(mL1+mL2-Ly1-Ly2-Lmy1-Lmy2)*TLy1*TLy2*facLy1*facLy2*facLmy1*facLmy11*facLmy2* ....
    factorial(Lmy1)/(factorial(Lmy1-2*ky1)*factorial(ky1))*factorial(Lmy11+Lmy2)/(factorial(Lmy11+Lmy2-2*ky2)*factorial(ky2))* ....
    ayL1^Lmy1*(conj(ayL2))^Lmy2*by^(mL1-2*Ly1-Lmy1)*conj(by)^(mL2-2*Ly2-Lmy2)*sig2^(Lmy1-ky1-Lmy11)*Dy_fs^(-Lmy1+ky1+ky2)* ....
    (-p1y+p2y-2*j*b*ro2_fs*py-j*ro2_fs*VyL1).^(Lmy1-2*ky1-Lmy11)*(sig2+0.25*ro2_fs)^Lmy11*ro2_fs^(ky1-Lmy11-Lmy2+2*ky2)* ....
    Py_fs^(-Lmy11-Lmy2+ky2).*Qy_fs.^(Lmy11+Lmy2-2*ky2);
end;end;end;end;end;end;end;
I_r = I_r + AL1*conj(AL2)*exp(-j*(tetaL1 + tetaL2))*sig2*ro2/sqrt(Px*Py)*exp(Ex).*exp(Ey).*sumx.*sumy;
I_r_fs = I_r_fs + AL1*conj(AL2)*exp(-j*(tetaL1 + tetaL2))*sig2*ro2_fs/sqrt(Px_fs*Py_fs)*exp(Ex_fs).*exp(Ey_fs).*sumx_fs.*sumy_fs;

end;end  
I_r = I_r*b^2*exp(j*b*(p1x^2-p2x^2+p1y^2-p2y^2))*exp(-((p1x-p2x)^2+(p1y-p2y)^2)/ro2);
I_r_fs = I_r_fs*b^2*exp(j*b*(p1x^2-p2x^2+p1y^2-p2y^2))*exp(-((p1x-p2x)^2+(p1y-p2y)^2)/ro2_fs);
tet = I_r./I_r_fs;I_r = I_r/maxgenarr(ru);I_r_fs = I_r_fs/maxgenarr(ru);
I_arr = cat(3,I_arr,I_r);I_arr_fs = cat(3,I_arr_fs,I_r_fs);%iL

tet_arr = cat(3,tet_arr,tet);end; % End of loop for one beam

%%%%%%%% Four Individual Intensity Plots %%%%%%%%%%%%%%%%%
for iL = 1:length(ALarr(:,1))
[px,py]=meshgrid(pxsetarr(iL,:),pysetarr(iL,:));pxn = px/1e-2;pyn = py/1e-2;
Lstr = [' \itL = \rm\bf' num2str(L/1e3) ' km'];figure(1);
I_r = I_arr(:,:,iL);subplot(2,2,iL);meshc(pxn,pyn,real(I_r));view([35 45]);colormap([0 0 0])
title(legtitle(iL,:),'FontSize',12,'FontWeight','bold');
xlabel('\itr_x\rm\bf axis in cm','FontSize',10,'FontWeight','bold','Rotation',-14);
ylabel('\itr_y\rm\bf axis in cm','FontSize',10,'FontWeight','bold','Rotation',25)
zlabel('\fontname{Verdana}\fontsize{12}\itI_n ','FontSize',12,'FontWeight','bold')    
set(gcf,'Renderer','Zbuffer');set(gcf,'Color','White');set(gca,'FontSize',11);
axis ([min(min(pxn))*1 max(max(pxn))*1 min(min(pyn))*1 max(max(pyn))*1 0 max(max(I_r))]);

legs1 = 'Slanted axis 1';legs2 = 'Slanted axis 2';
leg1 = {alfastrsx;' ';Vxstr};
if iL == 4
text(min(min(pxn))*0.95,max(max(pyn))*0.1,max(max(I_r))*(+1.3),leg1,'FontSize',11,'FontWeight','bold','Edgecolor','black','BackgroundColor','white');
text(min(min(pxn))*0.85,max(max(pyn))*0.5,max(max(I_r))*(+0.4),legc,'FontSize',11,'FontWeight','bold','Edgecolor','black','BackgroundColor','white','Rotation',0);
end;figure(2);subplot(2,2,iL);
[C,h] = contour(pxn,pyn,real(I_r),'-k');axis square;fax = 1;
line([min(min(pxn))*fax max(max(pxn))*fax],[0 0],'LineWidth',1.8,'color','black','LineStyle','--');
line([0 0],[min(min(pyn))*fax max(max(pyn))*fax],'LineWidth',1.8,'color','black','LineStyle','--');
title(legtitle(iL,:),'FontSize',12,'FontWeight','bold');
if iL == 2
line([min(min(pxn))*fax max(max(pxn))*fax],[min(min(pyn))*fax max(max(pyn))*fax],'LineWidth',2,'color','black','LineStyle','-');
text(min(min(pxn))*0.75,max(max(pyn))*0.8,max(max(I_r))*(+1.3),legs1,'FontSize',11,'FontWeight','bold','BackgroundColor','white');
end
if iL == 3
line([min(min(pxn))*fax max(max(pxn))*fax],[max(max(pyn))*fax min(min(pyn))*fax],'LineWidth',2,'color','black','LineStyle','-');
text(min(min(pxn))*0.75,max(max(pyn))*0.8,max(max(I_r))*(+1.3),legs2,'FontSize',11,'FontWeight','bold','BackgroundColor','white');
end
axis ([min(min(pxn))*fax max(max(pxn))*fax min(min(pyn))*fax max(max(pyn))*fax])
set(h,'LineWidth',1.2);H = clabel(C,h);set(H,'FontSize',6,'LineWidth',3)
set(gca,'FontSize',11);set(gcf,'Color',[1 1 1])
xlabel('\itr_x\rm\bf axis in cm','FontSize',12,'FontWeight','bold');
ylabel('\itr_y\rm\bf axis in cm','FontSize',12,'FontWeight','bold');
if iL == 4
text(min(min(pxn))*0.95,max(max(pyn))*0.1,max(max(I_r))*(+1.3),leg1,'FontSize',11,'FontWeight','bold','Edgecolor','black','BackgroundColor','white');
text(min(min(pxn))*0.85,max(max(pyn))*0.5,max(max(I_r))*(+0.4),legc,'FontSize',11,'FontWeight','bold','Edgecolor','black','BackgroundColor','white','Rotation',0);

end;
end;
%%%%%%% Functions Used %%%%%%%%%%%%%%%

function y = funb(a,b)
y = factorial(a)/factorial(b)/factorial(a-b);
function y = funT(el)
T = 1;
for j = 1:el
    T = (2*j-1)*T;
end
y = T;
function yout = H(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);