% 本程序用于MTF计算

% 输出
% MTF_onorbit_all_Np [垂轨，沿轨]动态MTF
% MTF_onorbit_all 二维动态MTF
% Pic_blurred 模糊图像

% 输入
% ratio = 1.0829; % 其物理意义是截止频率,lambda*F# = ratio*d
% MTF_static = 0.16; % 静态传函
% MTF_defocus = 0.975;  % 离焦
% MTF_zg = 0.98; % 杂光
% MTF_hp = 0.990; % 行频失配
% MTF_plj = 0.964; % 偏流角
% MTF_hot = 0.94; % 力热变形
% MTF_gas = 0.8; % 大气
% dp_push = 1; % 推扫方向的像移

function [MTF_onorbit_all_Np, MTF_onorbit_all, Pic_blurred] = MTF_simulation(ratio,MTF_static,MTF_defocus,MTF_zg,MTF_hp,MTF_plj,MTF_hot,MTF_gas,dp_push,Pic_origin)
%% 理想点扩散函数,Airy Disk, 直径为2.44*lambda*F#
d = 1; % 像元尺寸
Np = 1/2/d; % 奈奎斯特频率
dx = 0.1;
x = -5:dx:(5-dx);
y = x;
fx = -1/2/dx:1/2/5:(1/2/dx-1/2/5);
M = length(x);
[xx,yy] = meshgrid(x,y);
zz = sqrt(xx.^2+yy.^2);
psf_Airy = (2*besselj(1,pi*1/ratio*zz)./(pi*1/ratio*zz)).^2;
psf_Airy(M/2+1,M/2+1) = 1;
MTF_Airy = abs(fftshift(fft2(psf_Airy/sum(psf_Airy(:)))));

Np_site = round(0.5/dx)+M/2+1; % 奈奎斯特频率所在位置
MTF_Airy_Np = MTF_Airy(M/2+1,Np_site); % 奈奎斯特频率处MTF

% figure()
% mesh(xx,yy,psf)
% 
% figure()
% imshow(uint8(sqrt(psf)*255))

%% 静态传函
% 矩形孔径
psf_rect = rect(1/d*xx).*rect(1/d*yy);
MTF_rect = abs(fftshift(fft2(psf_rect/sum(psf_rect(:)))));
MTF_rect_Np = MTF_rect(M/2+1,Np_site);  % 奈奎斯特频率处MTF,理论值是0.637，因为采样问题，可能在0.6314
% figure()
% mesh(xx,yy,psf_rect);
% psf_rect = conv2(psf,psf_rect,'same');


ratio_static = MTF_static/MTF_rect_Np/MTF_Airy_Np; % 理想psf+探测器采样，与静态MTF还会有差距，
if ratio_static > 1
    disp('MTF_static is wrong,check it please!')
    return;
end

if ratio_static <= 1
    guass_sigma = sqrt(-log(ratio_static)*2)/2/pi/Np;
    psf_gauss = gaussmf(zz,[guass_sigma  0]);
    MTF_gauss = abs(fftshift(fft2(psf_gauss/sum(psf_gauss(:)))));
    MTF_gauss_Np =  MTF_gauss(M/2+1,Np_site);  % 奈奎斯特频率处MTF
end

psf_Airy_Guass = conv2(psf_Airy,psf_gauss,'same'); % 该函数用于生成图片时需要用到
psf_static = conv2(psf_Airy_Guass,psf_rect,'same'); % 静态psf
psf_static = psf_static/max(psf_static(:)); % 静态psf,最大值归一化
MTF_static_all = MTF_Airy.*MTF_rect.*MTF_gauss; % 全频段的MTF
MTF_static_NP = MTF_static_all(M/2+1,Np_site); % 奈奎斯特频率处的传函

%% 动态传函
ratio_orbit = MTF_defocus*MTF_zg*MTF_hp*MTF_plj*MTF_hot*MTF_gas; % 图省事
guass_sigma2 = sqrt(-log(ratio_orbit)*2)/2/pi/Np;
psf_gauss2 = gaussmf(zz,[guass_sigma2  0]);
MTF_gauss2 = abs(fftshift(fft2(psf_gauss2/sum(psf_gauss2(:)))));
MTF_gauss2_Np =  MTF_gauss2(M/2+1,Np_site);  % 奈奎斯特频率处MTF

psf_moving = rect(1/dp_push*yy).*rect(1/0.0001*xx); % 推扫PSF
MTF_moving = abs(fftshift(fft2(psf_moving/sum(psf_moving(:))))); % 推扫MTF

temp = conv2(psf_static,psf_gauss2,'same');
temp = temp/max(temp(:));
psf_onorbit_all = conv2(temp,psf_moving,'same');  % 动态PSF
psf_onorbit_all = psf_onorbit_all/max(psf_onorbit_all(:)); % 归一化
MTF_onorbit_all = MTF_moving.*MTF_gauss2.*MTF_static_all;  % 动态传函
MTF_onorbit_all_Np = [MTF_onorbit_all(M/2+1,Np_site);MTF_onorbit_all(Np_site,M/2+1)];  % 奈奎斯特频率处MTF

psf_onorbit_norect = conv2(conv2(psf_Airy_Guass,psf_gauss2,'same'),psf_moving,'same'); % 未算矩形采样psf，用于生成图
psf_onorbit_norect = psf_onorbit_norect/max(psf_onorbit_norect(:)); 

% figure()
% imshow(uint8(sqrt(psf_onorbit_all)*255))
% axis equal

% 保存曲线数据
dlmwrite('x.txt',x,'delimiter','\t');
dlmwrite('fx.txt',fx,'delimiter','\t');
dlmwrite('psf_Airy.txt',psf_Airy,'delimiter','\t');
dlmwrite('psf_rect.txt',psf_rect,'delimiter','\t');
dlmwrite('psf_static.txt',psf_static,'delimiter','\t');
dlmwrite('psf_onorbit_all.txt',psf_onorbit_all,'delimiter','\t');
dlmwrite('MTF_Airy.txt',MTF_Airy,'delimiter','\t');
dlmwrite('MTF_rect.txt',MTF_rect,'delimiter','\t');
dlmwrite('MTF_static_all.txt',MTF_static_all,'delimiter','\t');
dlmwrite('MTF_onorbit_all.txt',MTF_onorbit_all,'delimiter','\t');


%% 生成blurred图像
psf_onorbit_all_rect_CMOS  = psf_onorbit_norect(1:round(1/dx):end,1:round(1/dx):end);
psf_onorbit_all_rect_CMOS  = psf_onorbit_all_rect_CMOS/sum(psf_onorbit_all_rect_CMOS(:));

Pic_blurred = conv2(Pic_origin,psf_onorbit_all_rect_CMOS,'same');
end