
function [d,globalPC,localPC] = crisSDR_recrad(globalPC_filename,granule_filename,rflag);

%
% function [d,globalPC,localPC] = crisSDR_recrad(globalPC,filename,rflag);
%
% A function to produce a structure d with fields analogous to the output of concatsdr_pd.m
% but containing radiance variables containing reconstructed radiances.
%
% Inputs:
%    globalPC_filename:  Name of the NetCDF file containing the L1B gloabl PCs.  E.g.:
%                  globalPC_filename = 'globalPCs_L1B_version0_20220204.nc';
%
%    granule_filename:   Name of file of CrISPCA granule level data.  E.g.:
%                  granule_filename = 'CRISPCAvX_j01_d20220118_t0009439_e0010137_b21588_c20220118012935784364_ssec.nc';
%
%    rflag:   Flag to do the radiance reconstruction using only the Global PCs (2) or with Global and Local.
%
% Outputs:
%    d.ES_RealLW(717,9,30,4)   LW band reconstructed radiances
%    d.ES_RealMW(869,9,30,4)   MW band  "
%    d.ES_RealSW(637,9,30,4)   SW band  "
%
% Other:
%    - Uses of pca_expandM.m (appended at end)
%    - Execution time per granule is about 0.5s
%
% DCT, 25 Jul 2023
%

if nargin <= 2
  rflag = 0;
end

% Load global PC data
ncfile = globalPC_filename;
globalPC.U = ncread(ncfile,'U');  % PCs (hamming apodized)
globalPC.M = ncread(ncfile,'M');  % Mean (noise normalized) spectrum
globalPC.D = ncread(ncfile,'D');  % Eigenvalues
% Also load the wavenumbers (2223), the NEDN values used for noise 
% normalization
v = ncread(ncfile,'v');
global_nedn = ncread(ncfile,'nedn');
globalPC.nedn = global_nedn;

% Load the granule-specific data.  This includes the global 
% Coeficients and the local PCs and local coeficients.
ncfile = granule_filename;
localPC.U = ncread(ncfile,'U');
localPC.M = ncread(ncfile,'M');
localPC.D = ncread(ncfile,'D');
globalCoef_int16 = ncread(ncfile,'globalCoef_int16');
globalCoef_min = ncread(ncfile,'globalCoef_min');
globalCoef_rng = ncread(ncfile,'globalCoef_rng');
localCoef_int16 = ncread(ncfile,'localCoef_int16');
localCoef_min = ncread(ncfile,'localCoef_min');
localCoef_rng = ncread(ncfile,'localCoef_rng');

% Compute Coefs from integers and min + range
globalCoef = single(globalCoef_int16)*globalCoef_rng/(2^16-1)+globalCoef_min;
localCoef = single(localCoef_int16)*localCoef_rng/(2^16-1)+localCoef_min;

globalPC.Coef = globalCoef;
localPC.Coef = localCoef;

% For SDRs, get these from the GCRSO files
% Also load the other essential variables
%d.lat = ncread(ncfile,'lat');
%d.lon = ncread(ncfile,'lon');
%d.sol_zen = ncread(ncfile,'sol_zen');
%d.sat_zen = ncread(ncfile,'sat_zen');
%d.sat_azi = ncread(ncfile,'sat_azi');
%d.asc_flag = ncread(ncfile,'asc_flag');
%d.obs_time_tai93 = ncread(ncfile,'obs_time_tai93');


% Compute reconstructed radiances
sz = size(localCoef);
globalCoef = globalCoef(:,:);
localCoef = localCoef(:,:);
[~,nrec] = size(localCoef);
global_nednm = repmat(global_nedn,1,nrec);

if rflag == 0
  RecRad = pca_expandM(globalCoef,globalPC).*global_nednm + pca_expandM(localCoef,localPC).*global_nednm;
elseif rflag == 2
  RecRad = pca_expandM(globalCoef,globalPC).*global_nednm;  % global only
end
  
% Put the reconstructed radiances back into (nChannels,nFovs,nFORs,nScans) array
sz(1) = length(v);
RecRad = reshape(RecRad,sz);

% And store as normal L1B variable names/arrays
d.ES_RealLW = RecRad(1:717,:,:,:);
d.ES_RealMW = RecRad(718:1586,:,:,:);
d.ES_RealSW = RecRad(1587:2223,:,:,:);


function [x,rc] = pca_expandM( coef, PC )

[ntr,nX] = size(coef);
PCM = repmat(PC.M,1,nX);
col = PC.U(:,1:ntr) * coef + PCM;
x = col;
rc = 0;

