Merge remote-tracking branch 'origin/master' into develop

.gitattributes

+*.nc filter=lfs diff=lfs merge=lfs -text
+*.mat filter=lfs diff=lfs merge=lfs -text

RU/.gitattributes

+*.nc filter=lfs diff=lfs merge=lfs -text
+*.mat filter=lfs diff=lfs merge=lfs -text
+*.pdf filter=lfs diff=lfs merge=lfs -text

RU/Ep_approx.m

+function E_p = Ep_approx(prpt,alph,delta_ext,L_cbb,L_hbb,B_ssm,L,d_cbb,d_hbb)
+%
+% function E_p = Ep_approx(prpt,alph,delta_ext,L_cbb,L_hbb,B_ssm,L,d_cbb,d_hbb)
+% 
+% The SSM and sensor act like a polarizer-analyzer pair.  This function calculates
+% the approximate radiance error generated by the polarizer-analyzer pair
+%
+% Inputs
+% 	prpt:		p_r * p_t for sensor + SSM 		(unitless)	[nchan x 1]
+% 	alph:		analyzer angle alpha for sensor 	(radians)	[scalar]
+% 	delta_ext:	polarizer angle	for N scene views	(radians) 	[1 x N]
+% 	L_cbb:		Radiance for cold cal reference 	(Rad units))	[nChan x N]
+% 	L_hbb:		Radiance for hot cal reference 		(Rad units))	[nChan x 1]
+% 	B_ssm:		Planck radiance for scene mirror temp 	(Rad units))	[nChan x 1]
+% 	L:		Scene radiance for N scene views 	(Rad units))	[nChan x 1] or [nChan x N]
+% 	d_cbb:		polarizer angle for CBB view 		(radians)	[scalar]
+% 	d_hbb:		polarizer angle for HBB view 		(radians)	[scalar]
+%
+% Outputs
+% 	E_p		calibration bias due to polarization	(Rad units)	[nChan x N]
+%
+% 2020-08-19	JKT	
+% University of Wisconsin-Madison Space Science and Engineering Center (UW-SSEC)
+
+if size(L,2) == 1
+	L = repmat(L,1,length(delta_ext));
+elseif size(L,2) ~= length(delta_ext)
+	fprintf(1, 'Ep_approx:  ERROR Invalid scene radiance dimensions\n');
+end
+if abs(d_cbb) > 2*pi
+	fprintf(1, 'Ep_approx:  ERROR delta CBB out of range, should be in units of radians\n')
+	return
+end
+if abs(d_hbb) > 2*pi
+	fprintf('Ep_approx:  ERROR delta HBB out of range, should be in units of radians\n')
+	return
+end
+E_p = NaN*ones(length(L_cbb),length(delta_ext));                   % N_{d.wn} x N_{d.delta_ext} 
+
+for ii = 1:length(delta_ext)
+	d_ext = delta_ext(ii);
+	L_ext = L(:,ii);
+	% approximate equation
+	term1 = prpt.*(L_ext*cos(2*(d_ext - alph)));
+	term2 = -prpt.*(L_hbb.*((L_ext - L_cbb)./(L_hbb - L_cbb))*cos(2*(d_hbb - alph)));
+	term3 = -prpt.*(L_cbb.*((L_hbb - L_ext)./(L_hbb - L_cbb))*cos(2*(d_cbb - alph)));
+	term4 = -B_ssm.*prpt.*(cos(2*(d_ext - alph)));
+	term5 = +B_ssm.*prpt.*(((L_ext - L_cbb)./(L_hbb - L_cbb))*cos(2*(d_hbb - alph)));
+	term6 = +B_ssm.*prpt.*(((L_hbb - L_ext)./(L_hbb - L_cbb))*cos(2*(d_cbb - alph)));
+	E_p(:,ii) = term1 + term2 + term3 + term4 + term5 + term6;
+end
+return

RU/ICTradModelRU_L1b_aux.m

+function B = ICTradModelRU_L1b_aux(sensor,T_ICT,T_CRIS,ICT_Param,ru_term,ru_value)
+%
+% function B = ICTradModelRU_L1b_aux(sensor,T_ICT,T_CRIS,ICT_Param,ru_term,ru_value)
+%
+% Compute ICT predicted radiance using ICT radiometric model prescribed by ATBD
+% with perturbation to a single term
+%
+% Inputs
+%	sensor		sparse spectral params structure
+%	  .band
+%	  .v
+%	  .output_range
+%	T_ICT		ICT Temperature
+%	T_CRIS		CrIS temperature struct
+%	ICT_Param	ICT radiometric model parameters
+%	ru_term		string indicating ru_term to perturb
+%	ru_value	size of perturbation
+% 
+% Outputs
+%	B		Structure of predicted ICT radiance variables
+%	.emitted	emitted radiance
+%	.wn		wavenumber scale
+%	.SSM_Baffle	SSM Baffle term reflected radiance
+%	.ICT_baffle	ICT baffle term reflected radiance
+% 	.OMA		OMA term reflected radiance
+%	.BS_warm	BS_warm term reflected radiance
+% 	.BS_cold	BS_cold term reflected radiance
+%	.Space		Earth/Space term reflected radiance
+%	.reflected	total reflected radiance
+%	.total		total radiance (emitted and reflected)
+%	.V		view factors (struct)
+%	.e		emissivities (struct)
+%	.r		reflectivies (.SSM)
+%	.T		temperatures (struct)
+%	.e_ICT		ICT emissivity
+%	.T_ICT		ICT temperature
+%	.ru_term	RU term to perturb
+%	.ru_value	value of perturbation
+%	.no_perturb_val
+%	.perturbed_val
+%
+% 2020-09-13	JKT
+% University of Wisconsin-Madison Space Science and Engineering Center (UW-SSEC)
+
+wn = sensor.v_onaxis;
+
+ICT_wnum = ICT_Param.(sprintf('%s_ICT_wnum',lower(sensor.band)));
+ICT_emiss = ICT_Param.(sprintf('%s_ICT_emiss',lower(sensor.band)));
+ICT_Baffle_emiss = ICT_Param.(sprintf('%s_ICT_Baffle_emiss',lower(sensor.band)));
+Housing_emiss = ICT_Param.(sprintf('%s_Housing_emiss',lower(sensor.band)));
+ScanBaffle_emiss = ICT_Param.(sprintf('%s_ScanBaffle_emiss',lower(sensor.band)));
+ScanMirror_emiss = ICT_Param.(sprintf('%s_ScanMirror_emiss',lower(sensor.band)));
+Earth_emiss = ICT_Param.(sprintf('%s_Earth_emiss',lower(sensor.band)));
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% ICT emissivity interpolation to wavenumber grid
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Interpolate 4-min-EP emissivity to the on-axis sensor grid. 
+%% The wavenumber values corresponding to the emissivity values in the engineering packet are on the normal spectral resolution and on user grid
+e_ICT = interp1(ICT_wnum,ICT_emiss,wn,'nearest','extrap');
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Temperature, emissivity and view factor init
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%% SSM
+%%----------
+%r.SSM = 1-ICT_Param.Band(iband).Emissivity.ScanMirrorPts;				
+r.SSM = 1-ScanMirror_emiss;
+
+%% SSM_Baffle
+%%----------
+%e.SSM_Baffle = ICT_Param.Band(iband).Emissivity.ScanBafflePts;				
+e.SSM_Baffle = ScanBaffle_emiss;
+T.SSM_Baffle = T_CRIS.ssm_baffle_temp + T_CRIS.ssm_baffle_corr_temp;
+T.SSM_Baffle(isnan(T.SSM_Baffle)) = 0;
+V.SSM_Baffle = ICT_Param.View_Factor.ScanBaffle;                			
+
+%% ICT_baffle 
+%%----------
+%e.ICT_baffle = ICT_Param.Band(iband).Emissivity.ICT_BafflePts;  			
+e.ICT_baffle = ICT_Baffle_emiss;
+T.ICT_baffle = T_ICT;
+V.ICT_baffle = ICT_Param.View_Factor.ICTBaffle;       					
+
+%% OMA + Frame + Other
+%%----------
+%e.OMA = ICT_Param.Band(iband).Emissivity.HousingPts;					
+e.OMA = Housing_emiss;
+T.OMA = mean(0.5*T_CRIS.oma_struct_temp_1 + 0.5*T_CRIS.oma_struct_temp_2);
+V.OMA = ICT_Param.View_Factor.OMAandFrame;
+
+%% BS_warm
+%%----------
+e.BS_warm = e.OMA;									
+T.BS_warm = T.OMA;
+V.BS_warm = ICT_Param.View_Factor.BeamSplitterWarm;  					
+
+%% BS_cold
+%%----------
+T.BS_cold = T_ICT;
+V.BS_cold = ICT_Param.View_Factor.BeamSplitterCold;    					
+
+%% Space/Earth Reflected
+%%----------
+e.Space = Earth_emiss;
+T.Space = ICT_Param.Earth_Temperature;      						
+V.Space = ICT_Param.View_Factor.Space;   						
+
+no_perturb_val = [];	% initialize, will remain at this value for no perturbation
+perturbed_val = [];	% initialize, will remain at this value for no perturbation
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Perturbation 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+switch ru_term
+	case 'e_ICT'
+		no_perturb_val = e_ICT;
+		e_ICT = e_ICT*(1-ru_value);
+		e_ICT(e_ICT > 1) = 1;
+		e_ICT(e_ICT < 0) = 0;
+		perturbed_val = e_ICT;
+
+	case 'T_ICT'
+		no_perturb_val = T_ICT;
+		T_ICT = T_ICT + ru_value;
+		perturbed_val = T_ICT;
+
+	case 'SSM_scan_mirror_baffle'
+		no_perturb_val = T.SSM_Baffle;
+		T.SSM_Baffle = T.SSM_Baffle + ru_value;
+		perturbed_val = T.SSM_Baffle;
+
+	case 'T_ScanBaffleCorrection'
+		no_perturb_val = T.SSM_Baffle;
+		T.SSM_Baffle = T.SSM_Baffle + ru_value;
+		perturbed_val = T.SSM_Baffle;
+
+	case 'T_ICT_baffle'
+		no_perturb_val = T.ICT_baffle;
+		T.ICT_baffle = T.ICT_baffle + ru_value;
+		perturbed_val = T.ICT_baffle;
+
+	case 'T_OMA'					% OMA + Frame + Other
+		no_perturb_val = T.OMA;
+		T.OMA = T.OMA + ru_value;
+		perturbed_val = T.OMA;
+
+	case 'T_Frame'					% no longer used
+		no_perturb_val = T.Frame;
+		T.Frame = T.Frame + ru_value;
+		perturbed_val = T.Frame;
+
+	case 'T_BS_warm'
+		no_perturb_val = T.BS_warm;
+		T.BS_warm = T.BS_warm + ru_value;
+		perturbed_val = T.BS_warm;
+
+	case 'T_BS_cold'
+		no_perturb_val = T.BS_cold;
+		T.BS_cold = T.BS_cold + ru_value;
+		perturbed_val = T.BS_cold;
+
+	case 'T_Space'
+		no_perturb_val = T.Space;
+		T.Space = T.Space + ru_value;
+		perturbed_val = T.Space;
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% ICT emission
+%  (ref:  eq 75 in 474-00032, Rev C)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+B.emitted = e_ICT .* bt2rad(wn,T_ICT);
+B.wn = wn;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%  Reflected term 1: SSM Baffle 
+%  (ref:  eq 76 in 474-00032, Rev C)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if V.SSM_Baffle == 0
+	B.SSM_Baffle = zeros(size(B.wn));
+else
+	B.SSM_Baffle = (1 - e_ICT) .* e.SSM_Baffle * V.SSM_Baffle .* bt2rad(wn,T.SSM_Baffle);
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%  Reflected term 2:  ICT Baffle  
+%  (ref:  eq 78a in 474-00032, Rev C)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if V.ICT_baffle == 0
+	B.ICT_baffle = zeros(size(B.wn));
+else
+	B.ICT_baffle = (1 - e_ICT) .* e.ICT_baffle * V.ICT_baffle .* bt2rad(wn,T.ICT_baffle);
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%  Reflected term 3a,b,c:  OMA, frame, other, Beamsplitter Warm 
+% (ref:  eq 77 in 474-00032, Rev C)
+% V = 0.214 for OMA, Other, Frame terms in ATBD (= 0.0590 + 0.0100 + 0.1450)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% OMA
+%%----------
+if V.OMA == 0
+	B.OMA = zeros(size(B.wn));
+else
+	B.OMA = (1 - e_ICT) .* e.OMA * r.SSM * V.OMA .* bt2rad(wn,T.OMA);
+end
+
+%% BS_Warm
+%%----------
+if V.BS_warm == 0
+	B.BS_warm = zeros(size(B.wn));;
+else
+	B.BS_warm = (1 - e_ICT) .* e.BS_warm * r.SSM * V.BS_warm .* bt2rad(wn,T.BS_warm);
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%  Reflected term 4:  Beamsplitter Cold 
+%  (ref:  eq 78b in 474-00032, Rev C)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if V.BS_cold == 1
+	B.BS_cold = (1 - e_ICT) .* r.SSM^2 .* 0.5 * V.BS_cold .* bt2rad(wn,T.BS_cold);			% Consistent with ATBD
+else
+	B.BS_cold = (1 - e_ICT) .* r.SSM^2 .* 0.5 * V.BS_cold .* bt2rad(wn,T.BS_cold);			% SNPP
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%  Reflected term 5:  Space/Earth View  
+%  (ref:  eq 78c in 474-00032, Rev C)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if V.Space == 0
+	B.Space = zeros(size(B.wn));;
+else
+	B.Space = (1 - e_ICT) .* V.Space .* e.Space .* bt2rad(wn,T.Space);
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% total reflected radiance into ICT cavity
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+V.Total = V.SSM_Baffle + V.ICT_baffle + V.OMA + V.BS_warm + V.BS_cold + V.Space; 	% should total 1, for check only
+B.reflected = B.SSM_Baffle + B.ICT_baffle + B.OMA + B.BS_warm + B.BS_cold + B.Space;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%  Total radiance from ICT
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+B.total = B.emitted + B.reflected;
+
+B.wn = wn;
+B.V = V;
+B.e = e;
+B.r = r;
+B.T = T;
+B.e_ICT = e_ICT;
+B.T_ICT = T_ICT;
+B.ru_term = ru_term;
+B.ru_value = ru_value;
+B.no_perturb_val = no_perturb_val;
+B.perturbed_val = perturbed_val;
+
+return

RU/PLOT/cris_ru_plot.m

+function f_h = cris_ru_plot(RU,pname_out,fname_out);
+% function f_h = cris_ru_plot(RU,pname_out,fname_out);
+%
+% This is a function to plot the granule mean (all FOVs, all FORs, all scans) brightness temperature, 
+% and 3-sigma (k=3 coverage factor) RU in the same 9-panel format as figure 1 in the 
+% NASA Cross-track Infrared Sounder (CrIS) Level 1B Radiometric Uncertainty Description Document
+%
+% Inputs
+% RU		RU structure generated by cris_gran_RU_ncparam.m 
+% pname_out	full path to output directory where figure will be saved in .png format
+% fname_out	name of output file that figure will be saved to in .png format (.png extension will be added to fname_out)
+%
+% Outputs
+% f_h		handle to figure
+%
+% JKT, University of Wisconsin-Madison Space Science and Engineering Center (UW-SSEC)
+
+rucol(1,:) = [0.8431,0.1882,0.1529];
+rucol(2,:) = [0.9961,0.8784,0.5451];
+rucol(3,:) = [0.6510,0.8510,0.4157];
+rucol(4,:) = [0.1020,0.5961,0.3137];
+rucol(5,:) = [0.3294,0.1529,0.5608];
+rucol(6,:) = [0.8706,0.4667,0.6824];
+rucol(7,:) = [0,0,0];
+
+cmap_btscat = parula(128);
+
+%% Fontsize and Fontweight
+fs = 8;fw = 'n';
+
+ML = 0.1;		% margin left
+MR = 0.225;		% margin right
+MB = 0.2;		% margin bottom
+MT = 0.1;		% margin top
+HS = 0.03;		% horizontal spacing
+VS = 0.08;		% vertical spacing
+
+ylims_dbt = [-0.25 0.32];
+ylims_dbtscatter = [0 0.32];
+
+bands{1} = 'lw';
+bands{2} = 'mw';
+bands{3} = 'sw';
+
+w = 12; 
+h = 10;
+p = 0.01;
+
+f_h = figure;
+clf
+set(f_h,'Units','inches');
+screenpos = get(f_h,'Position');
+set(f_h,...
+	'Position',[screenpos(1:2) w h],...
+	'PaperUnits','inches',...
+	'PaperPosition',[p*w p*h w h],...
+	'PaperSize',[w*(1+2*p) h*(1+2*p)]);
+for ii = 1:length(bands)
+	clear leg_str
+	v = RU.(bands{ii}).v;
+	rad1 = squeeze(nanmean(RU.(bands{ii}).rad,2));
+	radu.T_ict = squeeze(nanmean(RU.(bands{ii}).ru_rad.dT_ict,2));
+	leg_str{1} = 'u(T_{ICT})';
+	radu.e_ict = squeeze(nanmean(RU.(bands{ii}).ru_rad.de_ict,2));
+	leg_str{2} = 'u(\epsilon_{ICT})';
+	radu.refl_meas = squeeze(nanmean(RU.(bands{ii}).ru_rad.refl_meas,2));
+	leg_str{3} = 'u(T_{Refl,meas})';
+	radu.refl_mod = squeeze(nanmean(RU.(bands{ii}).ru_rad.refl_mod,2));
+	leg_str{4} = 'u(T_{Refl,model})';
+	radu.polarization = squeeze(nanmean(RU.(bands{ii}).ru_rad.pol,2));
+	leg_str{5} = 'u(Pol Corr)';
+	radu.NLC = squeeze(nanmean(RU.(bands{ii}).ru_rad.dnlc,2));
+	leg_str{6} = 'u(NLC)';
+	radu.total = squeeze(nanmean(RU.(bands{ii}).ru_rad.total,2));
+	leg_str{7} = 'Total RU';
+	bt1 = rad2bt_l1b(v,rad1);
+	ru_flds = fieldnames(radu);
+	for iflds = 1:length(ru_flds)
+		btu.(ru_flds{iflds}) = rad2bt_l1b(v,rad1+radu.(ru_flds{iflds}))-bt1;
+	end
+
+	ru_ax(ii) = subplot(3,3,ii);
+	plot(v,nanmean(bt1(:,:),2),'color',[0.1294,0.4431,0.7098],'linewidth',1)
+	set(gca,'FontSize',fs,'FontWeight',fw,'XLim',[min(v)-20 max(v)+20],'YLim',[180 300])
+	if ii == 1
+		ylabel('BT (K)');
+	end
+	if ii == 2
+		title(sprintf('%s, All FOV mean',strrep(fname_out,'_','\_')));
+	end
+	grid on
+
+	ru_ax(ii+3) = subplot(3,3,ii+3);
+	hold on
+	for iflds = 1:length(ru_flds)
+		dy = btu.(ru_flds{iflds});
+		plot(v,nanmean(dy(:,:),2),'color',rucol(iflds,:),'linewidth',1)
+	end
+	set(gca,'FontSize',fs,'FontWeight',fw,'XLim',[min(v)-20 max(v)+20],'YLim',ylims_dbt)
+	xlabel('wavenumber');
+	if ii == 1
+		ylabel('3-sigma RU (K)')
+	end
+	grid on
+	if ii == 3
+		pos = get(gca,'position');
+		h_l = legend(leg_str);
+		h_pos = get(h_l,'Pos');
+		set(h_l,'Position',[1-h_pos(3),h_pos(2:4)]);
+	end
+
+	ru_ax(ii+6) = subplot(3,3,ii+6);
+	scatter(nanmean(bt1(:,:),2),nanmean(btu.total(:,:),2),5,v,'filled');
+	set(gca,'FontSize',fs,'FontWeight',fw,'YLim',[ylims_dbtscatter],'Box','on','XLim',[200 300])
+	xlabel('BT (K)');
+	if ii == 1
+		ylabel('3-sigma RU (K)');
+	end
+	grid on
+	pos = get(gca,'position');
+	hc = colorbar('hor');set(hc,'FontSize',8,'FontWeight',fw);
+	set(ru_ax(ii+6),'Position',pos);
+	set(get(hc,'XLabel'),'String','wavenumber','FontSize',8,'FontWeight','n')
+	cbpos = get(hc,'Pos');
+	set(hc,'Position',[cbpos(1) cbpos(2)+0.01 cbpos(3) cbpos(4)])
+	colormap(cmap_btscat)
+end
+print('-dpng',fullfile(pname_out,sprintf('RU_%s_%02d.png',fname_out,1)));
+
+return

RU/PLOT/rad2bt_l1b.m

+
+function bt = rad2bt_l1b(freq,radiance);
+%
+% function bt = rad2bt_l1b(freq,radiance);
+%
+% compute brightness temperature given wavenumbers and radiances.
+%
+% Inputs:
+%	freq	  wavenumbers 			[NCHAN x 1] 			units:  1/cm 
+%	radiance  radiances 			[NCHAN x NFOV x NFOR x NSCAN]	units:  mW/(m^2 sr. 1/cm)
+% Outputs:
+%	bt	  computed brightnes temp	[NCHAN x NFOV x NFOR x NSCAN]	units:  Kelvin
+%
+% JKT, University of Wisconsin-Madison Space Science and Engineering Center (UW-SSEC)
+% based on rad2bt by DCT 1999-11-11
+%
+% not required for RU calculation, only for cris_ru_plot.m which is provided as an example of
+% plotting the CrIS RU in brightness temperature units [K]
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% fundamental constants:
+%  (Cohen, E.R. and B.N. Taylor, The 1986 CODATA recommended values
+%  of the fundamental physical constants, Journal of Research of
+%  the National Bureau of Standard, 92(2), March-April 1987.)
+
+h = 6.6260755E-34;  % Planck constant in Js
+c = 2.99792458E8;   % photon speed in m/s
+k = 1.380658E-23;   % Boltzmann constant in J/K
+
+c1 = 2*h*c*c*1e8;
+c2 = h*c/k*1e2;
+
+freq = freq(:);		% make sure freq is N x 1 and not 1 x N
+[NCHAN,NFOV,NFOR,NSCAN] = size(radiance);
+[NCHAN_freq] = length(freq);
+if NCHAN_freq ~= NCHAN
+	error('incompatible wavenumber and radiance dimensions');
+end
+
+bt = nan*radiance;
+freq = repmat(freq,1,NFOV,NFOR,NSCAN);
+bt = c2.*freq./(log(1+(c1.*freq.*freq.*freq)./(radiance/1e3)));
+
+ind_neg = find(radiance <= 0);
+ind_imag = find(imag(radiance) ~= 0);
+
+if length(ind_neg) ~= 0
+	warning('WARNING: negative inputs to RAD2BT');
+	warning('Setting resulting BT values to NaN');
+	bt(ind_neg) = NaN;
+end
+if length(ind_imag) ~= 0
+	warning('WARNING: imaginary inputs to RAD2BT');
+	warning('Setting resulting BT values to NaN');
+	bt(ind_imag) = NaN;
+end
+return
+
+

RU/README.md

+1.  Download the code
+
+-  go to https://gitlab.ssec.wisc.edu/cris_l1b/user/-/tree/master/RU
+
+-  Click download and select preference for download format (.zip,
+    .tar.gz, .tar.bz2, .tar)
+
+-  Decompress the downloaded archive file on your local machine
+
+2.  Run the code
+
+To do a test run of the code in matlab, you will need to modify the User
+Input section in the wrapper script (cris_calbudget_1gran_nc_wrapper.m)
+to set
+
+-   the path to the input v3 L1b granule (two test input granules are
+    provided in the TEST_INPUT directory) (l1b_pname)
+
+-   the satellite name (sat_name, ‘SNPP’ or ‘J1’)
+
+-   the path to the radiometric uncertainty code that you have
+    downloaded or cloned from gitlab (ru_code_pname)
+
+-   the save RU data flag (saveTotalRU_FLAG, set to 1 to save RU file)
+
+-   the output path (pname_out)
+
+3.  Compare your output to sample output
+
+-   Sample radiometric uncertainty files that were generated from the
+    sample input granule files are provided in the TEST_OUTPUT directory

RU/bt2rad.m

+function radiance = bt2rad(freq,bt);
+%
+% function radiance = bt2rad(freq,bt);
+%
+% compute radiance given wavenumbers and brightness temperature.
+%
+% Inputs:
+%	freq	  wavenumbers (Nx1) in 1/cm 
+%	bt        brightnes temperature (Nx1) in Kelvin
+% Outputs:
+%	radiance  Planck radiances (Nx1) in mW/(m^2 sr. 1/cm)
+%
+% see also RADTOT.M, RAD2BT.M, TTORAD.M
+% 
+% DCT 11/11-99
+% University of Wisconsin-Madison Space Science and Engineering Center (UW-SSEC)
+%
+
+% fundamental constants:
+%  (Cohen, E.R. and B.N. Taylor, The 1986 CODATA recommended values
+%  of the fundamental physical constants, Journal of Research of
+%  the National Bureau of Standard, 92(2), March-April 1987.)
+
+h = 6.6260755E-34;  % Planck constant in Js
+c = 2.99792458E8;   % photon speed in m/s
+k = 1.380658E-23;   % Boltzmann constant in J/K
+
+c1 = 2*h*c*c*1e8;
+c2 = h*c/k*1e2;
+
+radiance = 1e3 * c1.*(freq.*freq.*freq)./(exp((c2.*freq)./bt)-1);
+
+ind = find(bt < 0);
+if length(ind) ~= 0
+	radiance = nan*real(radiance);
+	disp('WARNING: negative inputs to BT2RAD')
+end
+ind = find(~isreal(bt));
+if length(ind) ~= 0
+	radiance = nan*real(radiance);
+	disp('WARNING: imaginary inputs to BT2RAD')
+end
+ind = find(isnan(bt));
+if length(ind) ~= 0
+	radiance = nan*real(radiance);