
% Script to calculate g value for Sentinel-2's red (or green) band using the
% Sentinel-2 spectral response function data for the red band and applying
% the wavelength parameters from Pope and Fry (1997) for a and from
% Buiteveld et al. (1994) for b.

% The relationship (following Pope et al., 2016) is: 
% g = 2a + b 
% where a is the absorption coefficient and b is the scattering
% coefficient. 

% For details of the method, see supplement to the paper on The Cryosphere
% (https://doi.org/10.5194/tc-2018-56). 

%% Enter Sentinel-2 spectral response function and corresponding wavelengths 

S2_wavelengths = [646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684]; 

S2_spec_response = [0.002584161
0.034529234
0.149970293
0.464825627
0.81774561
0.965323761
0.983869281
0.996900219
1
0.995449387
0.991334334
0.977215143
0.936802357
0.873776426
0.814165825
0.776668942
0.764863979
0.77509069
0.801358663
0.830827664
0.857111548
0.883580858
0.908949643
0.934759132
0.955930529
0.968110072
0.973218918
0.971572352
0.969002767
0.965711927
0.960481316
0.944810855
0.884152441
0.706167011
0.422966946
0.189853092
0.063171769
0.020615018
0.00203356]; 

%% Enter measured data for a values from Pope and Fry (1997) 

pope_fry_wavelength = [645.0
647.5
650.0
652.5
655.0
657.5
660.0
662.5
665.0
667.5
670.0
672.5
675.0
677.5
680.0
682.5
685.0]; 

pope_fry_a = [0.325
0.335
0.340
0.358
0.371
0.393
0.410
0.424
0.429
0.436
0.439
0.448
0.448
0.461
0.465
0.478
0.486]; 

data = [pope_fry_wavelength,pope_fry_a]; 
cbar_div_num = 81; 

siz_cmap2 = size(pope_fry_a,1); 
x2=linspace(-1,1,2);
y2=linspace(-1,1,siz_cmap2);
yq2=linspace(-1,1,cbar_div_num);
[X2,Y2]=meshgrid(x2,y2);
[Xq2,Yq2]=meshgrid(x2,yq2);
interpolated_a = interp2(X2,Y2,data,Xq2,Yq2);

extracted_a = interpolated_a(3:2:end-1,:); % get the a values for S2 wavelengths 

%% Enter measured data for b from Buiteveld et al. (1994) 

buiteveld_b = [0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0007
0.0006
0.0006
0.0006
0.0006
0.0006
0.0006
0.0006
0.0006
0.0006
0.0006
0.0006
0.0006
0.0006
0.0006
0.0006
0.0006
0.0006
0.0006
0.0006
0.0006
0.0006
0.0006]; 

%% Set up calculation array for ease of following 

S2_wavelengths = extracted_a(:,1); 
S2_a = extracted_a(:,2); 
S2_calcs = [S2_wavelengths,S2_spec_response,S2_a,buiteveld_b]; 

%% Do the calculations as per Pope et al. supplement 

S2_2a = S2_a .* 2; 
S2_2a_b = S2_2a + buiteveld_b;  

S2_calcs2 = [S2_calcs,S2_2a_b]; 

%% Weight g values according to spectral response function data 

sum_all_spec = sum(S2_calcs2(:,2)); 

weightings = zeros(size(S2_calcs2,1),1);

for i = 1:length(S2_calcs2) 
    
    spec_resp = S2_calcs2(i,2); 
    spec_resp_g = S2_calcs2(i,5); 
    
    weight = spec_resp_g * spec_resp; 
    weightings(i,1) = weight; 
    
end 

g_value_sum = sum(weightings(:,1)); 
g_value = g_value_sum / sum_all_spec; 
