function [ reflect ] = ls8_dn2refl(img,meta,band)
%ls8_dn2refl Convert landsat 8 digital numbers reflectance
%   Uses formulae at http://landsat.usgs.gov/Landsat8_Using_Product.php
% reflect=output reflectance
% img = landsat 8 band image
% meta = scene metadata variable (create with ls8_meta_parser)
% band = image band
% function written by Neil Arnold 

% Get relevant metadata
switch band
    case 1
        band_mult=meta.RADIOMETRIC_RESCALING.REFLECTANCE_MULT_BAND_1;
        band_add=meta.RADIOMETRIC_RESCALING.REFLECTANCE_ADD_BAND_1;
    case 2
        band_mult=meta.RADIOMETRIC_RESCALING.REFLECTANCE_MULT_BAND_2;
        band_add=meta.RADIOMETRIC_RESCALING.REFLECTANCE_ADD_BAND_2;
    case 3
        band_mult=meta.RADIOMETRIC_RESCALING.REFLECTANCE_MULT_BAND_3;
        band_add=meta.RADIOMETRIC_RESCALING.REFLECTANCE_ADD_BAND_3;
    case 4
        band_mult=meta.RADIOMETRIC_RESCALING.REFLECTANCE_MULT_BAND_4;
        band_add=meta.RADIOMETRIC_RESCALING.REFLECTANCE_ADD_BAND_4;
    case 5
        band_mult=meta.RADIOMETRIC_RESCALING.REFLECTANCE_MULT_BAND_5;
        band_add=meta.RADIOMETRIC_RESCALING.REFLECTANCE_ADD_BAND_5;
    case 6
        band_mult=meta.RADIOMETRIC_RESCALING.REFLECTANCE_MULT_BAND_6;
        band_add=meta.RADIOMETRIC_RESCALING.REFLECTANCE_ADD_BAND_6;
    case 7
        band_mult=meta.RADIOMETRIC_RESCALING.REFLECTANCE_MULT_BAND_7;
        band_add=meta.RADIOMETRIC_RESCALING.REFLECTANCE_ADD_BAND_7;
    case 8
        band_mult=meta.RADIOMETRIC_RESCALING.REFLECTANCE_MULT_BAND_8;
        band_add=meta.RADIOMETRIC_RESCALING.REFLECTANCE_ADD_BAND_8;
    case 9
        band_mult=meta.RADIOMETRIC_RESCALING.REFLECTANCE_MULT_BAND_8;
        band_add=meta.RADIOMETRIC_RESCALING.REFLECTANCE_ADD_BAND_8;
end
        
% TOA reflectance (no solar correction)
toa=band_mult*double(img)+band_add;

% Get solar zenith angle from metadata
solzen=90-meta.IMAGE_ATTRIBUTES.SUN_ELEVATION;

% Calculate reflectance with solar correction
reflect=toa/cosd(solzen);

end

