% 本程序用于信噪比计算

function [SNR_output, mmm_pho_SNR, Picture_output] = SNR_simulation(Sensor_input,Optic_input,T_imaging_input,Light_input,pho_input,pho_L_infor,site_str)

%% Sensor_input 探测器性能参数输入
I_dark = Sensor_input(1); % 暗电流电子数，单位e-/s
Fn = Sensor_input(2); % 噪声因子，无量纲
QE = Sensor_input(3); % 平均量子效率，无量纲
Gain = Sensor_input(4); % 增益系数，无量纲
Noise_read = Sensor_input(5); % 读出噪声，单位e-
FWC = Sensor_input(6); % 满阱电子数，单位e-
ADC = Sensor_input(7); % 量化位数，单位bit

% 如果未指定探测器性能参数，则采用默认参数（HR4256性能参数）
if isnan(I_dark)
    I_dark  = 100*1000; % 暗电流噪声电子数，单位e-/s
end
if isnan(Fn)
    Fn = 1; % 噪声因子，无量纲
end
if isnan(QE)
    QE = 0.6; % 平均量子效率，无量纲
end
if isnan(Gain)
    Gain = 2.0; % 增益系数，无量纲
end
if isnan(Noise_read)
    Noise_read = 120; % 读出噪声，单位e-
end
if isnan(FWC)
    FWC = 140*1000; % 满阱电子数，单位e-
end
if isnan(ADC)
    ADC = 12; % 量化位数，单位bit
end



%% Optic_input 光学系统性能参数
Lambda_center = Optic_input(1); % 中心波长,单位nm
D_aperture = Optic_input(2); % 光学系统口径，单位m
GSD = Optic_input(3); % 星下点分辨率，单位m
H_orbit = Optic_input(4); % 轨道高度，单位km
Tao0 = Optic_input(5); % 光学系统透过率，无量纲
Epslong = Optic_input(6); % 光学系统面遮拦比，无量纲

% 如果未指定光学系统性能参数，则采用默认参数（商业星性能参数）
if isnan(Lambda_center)
    Lambda_center = 650; % 中心波长,单位nm
end
if isnan(D_aperture)
    D_aperture = 1.2; % 光学系统口径，单位m
end
if isnan(GSD)
    GSD = 0.25; % 星下点分辨率，单位m
end
if isnan(H_orbit)
    H_orbit = 500; % 轨道高度，单位km
end
if isnan(Tao0)
    Tao0 = 0.8; % 光学系统透过率，无量纲
end
if isnan(Epslong)
    Epslong = 0.15; % 光学系统面遮拦比，无量纲
end

%% 成像参数输入
t_TDI = T_imaging_input(1); % 积分时间，单位s，可以不输入，改输入积分级数M
M_TDI = T_imaging_input(2); % 积分级数，无量纲，整数

% 如果未指定成像参数，则采用默认参数（商业星性能参数） 
if isnan(M_TDI)
    M_TDI = 64; % 积分级数，无量纲，整数
end
if isnan(t_TDI)
    GM = 398603*10^9; % G:万有引力常数，M：地球质量
    R = 6371*1000; % R: 地球半径
    H = H_orbit*1000; % H: 轨道高度,单位m
    rs = R + H;% 半长轴a，单位km
    Vx = sqrt(GM/rs); % 飞行速度
    ws = Vx/rs; % 地面飞行速度
    GSD_Vx = R*ws; % 角速度 度/每秒
    t_TDI_single = GSD/GSD_Vx; % 单行积分时间
    t_TDI = M_TDI*t_TDI_single ; % 积分时间，单位s
end

%% 入瞳辐亮度输入
% 必须输入
% 若不输入，则默认高中低三个参数
if isnan(Light_input)
    Light_input = [4.35,13.36,55.7]; % 500-900nm,北纬夏季，10°@0.05，30°@0.2， 60°@0.65，单位W/m2/sr
end
h_planck = 6.63*1E-34; % 普朗克常数，单位Js
c_light = 3.0*1E8; % 光速，常量，单位m/s

%% 计算信噪比
S_signal_with_L = QE*(GSD/H_orbit/1000)^2*pi*(1-Epslong)/4*Tao0*(Lambda_center*1E-9)/(h_planck*c_light)*D_aperture^2*t_TDI; % 信号电子数/辐射亮度 的系数
S_signal = S_signal_with_L*Light_input; % 信号电子数
S_signal(S_signal > FWC/Gain) = FWC/Gain; % 信号电子数大于满阱，表示探测器饱和了
SNR_ratio = Gain*S_signal./sqrt(Fn^2*Gain^2*S_signal + Gain^2*I_dark*t_TDI+ Gain^2*Noise_read^2);  % 信噪比，倍数
SNR_output = 20*log10(SNR_ratio); % 信噪比，dB


%% pho_input 反射率矩阵或者图片数据
% 反射率分布，必须输入
if isnan(pho_input)
    Picture_output = nan;
    mmm_pho_SNR = nan;
else
    if length(size(pho_input)) >= 3 % 表示输入为rgb图像
        pho_input = rgb2gray(pho_input);
        pho_input = im2double(pho_input);
    else
        %         length(size(pho_input)) <= 2 % 表示输入为图像或者矩阵
        if max(pho_input)<=1.0 % 表述输入为0-1的矩阵,已经是double型
            ;
        else
            pho_input = im2double(pho_input);
        end
    end
    
    %% pho_L_infor 最大反射率，辐射亮度L与反射率pho的k,b输入
    pho_max = pho_L_infor(1); % 最大的反射率
    if isnan(pho_max)
        pho_max = 0.9; % 如果不输入，那就按0.9，云的反射率
    end
    if ~isnan(pho_max)
        pho_input = pho_max*pho_input/max(max(pho_input)); % 认为最大反射率pho_max
    end    
    
    %% 仿真输出  图片为Picture_output.tiff
    [sizeM,sizeN] = size(pho_input); %求尺寸
    pho_L_k = pho_L_infor(2);  % pho到辐射亮度L的斜率k
    pho_L_b = pho_L_infor(3);  % pho到辐射亮度L的本底b
    
    % 计算最大/平均/最小反射率，最大/平均/最小信噪比
    mmm_pho = [max(pho_input(:)),mean(pho_input(:)),min(pho_input(:))]; % 平均反射率
    S_signal_mmm = S_signal_with_L*(pho_L_k*mmm_pho + pho_L_b); % 信号电子数
    S_signal_mmm(S_signal_mmm  > FWC/Gain) = FWC/Gain; % 信号电子数大于满阱，表示探测器饱和了
    SNR_ratio_mmm  = Gain*S_signal_mmm./sqrt(Fn^2*Gain^2*S_signal_mmm  + Gain^2*I_dark*t_TDI+ Gain^2*Noise_read^2);  % 信噪比，倍数
    SNR_output_mmm  = 20*log10(SNR_ratio_mmm); % 信噪比，dB
    mmm_pho_SNR = [mmm_pho; SNR_output_mmm];
    
    % 图像生产：输出信号→DN
    L_light = pho_L_k*pho_input + pho_L_b; % 辐射亮度
    S_signal = S_signal_with_L*L_light; % 信号电子数
    % noise_white = 1/0.289*(rand(sizeM,sizeN)-0.5); %  均值为0，rms为1
    noise_S_signal = Fn*sqrt(S_signal).*(1/0.289*(rand(sizeM,sizeN)-0.5)); % 光子噪声
    noise_S_dark = sqrt(I_dark*t_TDI).*(1/0.289*(rand(sizeM,sizeN)-0.5));  % 暗电流噪声
    noise_S_read = Noise_read*(1/0.289*(rand(sizeM,sizeN)-0.5));  % 读出噪声
    S_signal_output = S_signal + I_dark*t_TDI + noise_S_signal + noise_S_dark + noise_S_read; % 总电子数（信号+暗电流+噪声）
    S_signal_output(S_signal_output>FWC/Gain) = FWC/Gain; % 饱和了
    DN_signal_output = round(S_signal_output/(FWC/Gain)*(2^ADC -1)); %
    
    %imwrite(uint8(DN_signal_output/max(DN_signal_output(:))*255),[site_str,'Picture_output_LaShen.png']);
    imwrite(uint8(DN_signal_output/max(DN_signal_output(:))*255),site_str+'Picture_output_LaShen.png');
    if ADC > 16;
        disp('不支持16bit以上量化！！！')
    elseif (ADC >8)&&(ADC<=16)
        Picture_output = uint16(DN_signal_output); % 输出的
        %imwrite(Picture_output,[site_str,'Picture_output.tiff']);
        imwrite(Picture_output,site_str+'Picture_output.tiff');
    else
        Picture_output = uint8(DN_signal_output); % 输出的
        %imwrite(Picture_output,[site_str,'Picture_output.tiff']);
        imwrite(Picture_output,site_str+'Picture_output.tiff');
    end
end