function [C]=WV2exoref(gy,gm,gd,sz)
%
% Calculates WV2 exoatmospheric (Planetary; Top Of Atmosphere) reflectance conversion factor for radiance to reflectance.
% Based on Chander & Markham (2003) and Updike & Comp (2010).
% http://landsathandbook.gsfc.nasa.gov/pdfs/L5TMLUTIEEE2003.pdf
% http://www.digitalglobe.com/downloads/
% Radiometric_Use_of_WorldView-2_Imagery.pdf
%
%                Small, 2012
%
%  Input:
%       gy,gm,gd = gregorian date [y,m,d]
%       sz = solar zenith angle (90-solar elevation) in degrees
%
%  Output:
%      C = conversion factors from spectral radiance (watts/m^2 * ster * um) to unitless reflectance for WV2 bands 1-8.
%
% No warrantee given or implied.  Use at your own risk.  Not for navigation

format long

% Earth-Sun distance vs julian day

esd=[1, 0.9832;  15, 0.9836;  32, 0.9853;  46, 0.9878;  60, 0.9909;  74, 0.9945; 91, 0.9993;  106, 1.0033; 121, 1.0076; 135, 1.0109; 152, 1.0140; 166, 1.0158; 182, 1.0167; 196, 1.0165; 213, 1.0149; 227, 1.0128; 242, 1.0092; 258, 1.0057; 274, 1.0011; 288, 0.9972; 305, 0.9925;319, 0.9892;335, 0.9860;349, 0.9843; 365, 0.9833]; 

% Gregorian to Julian;  Calls Rich Signell's julian.m function from
% http://woodshole.er.usgs.gov/staffpages/rsignell/rsignell.html

jd=julian(gy,gm,gd)-julian(gy,1,1)+1

% solar emission by L7 band
% esun=[1969.000,1840.000,1551.000,1044.000,225.700,82.07];

% solar emission by WV2 band
esun=[1758.2229,1974.2416,1856.4104,1738.4791,1559.4555,1342.0695,1069.7302,861.2866];

% Spectral Band Spectral Irradiance [W m^-2 um^-1]
% Panchromatic 1580.8140
% Coastal 1758.2229
% Blue 1974.2416
% Green 1856.4104
% Yellow 1738.4791
% Red 1559.4555
% Red Edge 1342.0695
% NIR1 1069.7302
% NIR2 861.2866


% seasonal Earth-Sun astronomical distance (from table)
ad=interp1q(esd(:,1),esd(:,2),jd);

% Radiance->ExoReflectance conversion factor
C=(pi*ones(size(esun))*ad^2)./(esun*cos(sz*pi/180));

