Matlab实现热晕

clc;
clear;  


bc=1.06e-6;      % 波长
lamda=1.06e-6;    % 波长
 wo=0.3; % 光束束腰半径（常师兄设置为0.3）
 a=0.3;   % 设置束腰半径(在这里要注意对束腰半径的定义)
power=1e5;% 激光功率
L=0.6;% 网格大小
N=512;
k=2*pi/bc;% 波数
z=10000;%z距离  %1km-10km
 dz=500;
 bushu=z/dz;
 v=4;  % 风速  1.5以上
 
 zf=5000; % 描述激光波束的焦距
gamma=1.4; 
Cs=330; % 声速
A0=(2*power / (pi*wo.^2))^(1/2);
  
buchang=2*L/N;
[x,y]=meshgrid(-L:buchang:L-buchang,-L:buchang:L-buchang); %正负不等
xh=x(1,:);
yh=y(:,1);

  alph=4.06*10^(-7)+2.57*10^(-6)+3.63*10^(-5)+2.93*10^(-4); %消光系数  =吸收系数+散射系数       使用中纬度地区0-1km夏季 lamda=1.06μm
 xishou=2.57*10^(-6)+3.63*10^(-5); %吸收系数
 

 
%  Eu=A0.*(exp(-(x.^2+y.^2)/wo.^2));     %  高斯光束
  Eu=A0.*(exp(-(x.^2+y.^2)/wo.^2)).*exp(-1i*k*(x.^2+y.^2)/2/zf);%聚焦

delta=L/N;
[x1,y1]=meshgrid(-L/buchang:1:L/buchang-1,-L/buchang:1:L/buchang-1);
r = sqrt(x.^2 + y.^2) ;
fanghuqiang = exp(-(r/(0.95*L)).^10) ;   %超高斯分布吸收边界
zhenkong=zhenkongchuanshu(L,x1,y1,dz,k);
h=waitbar(0,'计算中，请等待...');

for l=1:bushu
     if l==1
            Eu1=myifft2(zhenkong.*myfft2(Eu)).*exp(-(alph)*dz/2)^0.5;
     else
            Eu1=myifft2(zhenkong.*myfft2(Eu2)).*exp(-(alph)*dz/2)^0.5.*fanghuqiang;

     end

      Ip=abs(Eu1).^2;

  
   

 I=(abs(Eu2)).^2;
 phase1=angle(Eu);%返回复数
 phase2=angle(Eu2);%返回 

  %pic=imagesc(xh,yh,I)%彩图or灰度图
  figure(1);
  %前两个参数用于设置figure的位置，后两个参数用于设置figure的大小（长和宽）；
  
  %picI=imshow(I,[ ]);
  picI=imshow(I,[ ],'border','tight','initialmagnification','fit');
  set(gcf,'Position',[0 30 1280 1024]);
  %set(gca,'Position',[0.01 0.01 2 2]);

%set (gcf,'Position',[0,0,500,500]);
axis normal;

saveas(gcf,'1.bmp','bmp');
%print(gcf,'-dbmp','bmp')
  
  %saveas(picI,'1.bmp');
  
  figure(2);
 
%   set(gcf,'Position',[0 0 1280 1024], 'color','w');
  pic=imagesc(xh,yh,I);
   axis('off');
   saveas(pic,'2.bmp');
  figure(3)
  picphase1=imshow(phase2,[ ]);%传输过后phase
  figure(4)
  picphase2=imshow(phase1,[ ],'border','tight','initialmagnification','fit');;%未传输phase
%    set(gcf,'Position',[0 30 1280 1024]);
   axis normal;
  saveas(gcf,'2.bmp','bmp');