clc;clear
%%
%计算Esun
%把原始太阳强度读进来。来源于马一川
load('sun')

%用GF5B的 中心波长和半宽算它的Esun
filename = 'G:\迅雷下载\假阳性\ZY1E\ZY1E_AHSI_E124.72_N46.77_20230602_019496_L1A0000610672\ZY1E_AHSI_SWIR_Spectralresponse.raw';
fid = fopen(filename,'r');
if fid<0
    warndlg('打开文件失败!');
    return;
else
        FormatString=repmat('%f ',1,2);%取前2列数据
        FWHM_SWIR = cell2mat(    textscan(fid,FormatString,180,'HeaderLines',0,'delimiter',',')   );
end
fclose(fid);
filename = 'G:\迅雷下载\假阳性\ZY1E\ZY1E_AHSI_E124.72_N46.77_20230602_019496_L1A0000610672\ZY1E_AHSI_VNIR_Spectralresponse.raw';
fid = fopen(filename,'r');
if fid<0
    warndlg('打开文件失败!');
    return;
else
        FormatString=repmat('%f ',1,2);%取前2列数据
        FWHM_VNIR = cell2mat(    textscan(fid,FormatString,150,'HeaderLines',0,'delimiter',',')   );
end
fclose(fid);

FWHM = [FWHM_VNIR;FWHM_SWIR]; % VNIR在前，SWIR在后

%公式参考HAPI.py
step = 1; %需要与sun的间隔保持一致,1nm
AF_Wing = 5;
x = -AF_Wing:step:AF_Wing;
for i = 1:length(FWHM)
    wavelength = FWHM(i,1);
    [~,index] = min(abs(sun(:,1) - wavelength));
    energy(i,:) = sun(index-AF_Wing/step : index+AF_Wing/step,2);
end
for i = 1:length(FWHM)
    g = FWHM(i,2); %g为半宽
    slit(i,:) = 2 / g * sqrt(log(2) / pi) * exp(-log(2)*(2*x/g).^2);  %该处做了改动
    slit(i,:) = slit(i,:) / (sum(slit(i,:))*step);
end
Esun = sum(energy .* slit,2) * step;

save('FWHM_SWIR.mat','FWHM_SWIR');
save('FWHM_VNIR.mat','FWHM_VNIR');
save('Esun.mat','Esun');
%%
%从txt里读 scale factor和offset

clc;clear
filename = 'G:\迅雷下载\假阳性\ZY1E\ZY1E_AHSI_E124.72_N46.77_20230602_019496_L1A0000610672\ZY1E_AHSI_SWIR_RadCal.raw';
fid = fopen(filename,'r');
if fid<0
    warndlg('打开文件失败!');
    return;
else
        FormatString=repmat('%f ',1,2);%取前2列数据
        swir_rad_cal = cell2mat(    textscan(fid,FormatString,180,'HeaderLines',0,'delimiter',',')   );
end
fclose(fid);
filename = 'G:\迅雷下载\假阳性\ZY1E\ZY1E_AHSI_E124.72_N46.77_20230602_019496_L1A0000610672\ZY1E_AHSI_VNIR_RadCal.raw';
fid = fopen(filename,'r');
if fid<0
    warndlg('打开文件失败!');
    return;
else
        FormatString=repmat('%f ',1,2);%取前2列数据
        vnir_rad_cal = cell2mat(    textscan(fid,FormatString,150,'HeaderLines',0,'delimiter',',')   );
end
fclose(fid);
save('swir_rad_cal.mat','swir_rad_cal');
save('vnir_rad_cal.mat','vnir_rad_cal');