srf_generate(center, fwhm, shape, n, yedge)

function [wl, y, wn] = srf_generate(center, fwhm, shape, n, yedge)
% SRF_GENERATE Generate a MODTRAN-compatible spectral response function 
%
% [wl, y, wn] = srf_generate(center, fwhm, shape, n, yedge)
%
% Inputs:
%   center    center wavelength
%     fwhm    desired full width half max
%    shape    'gauss', 'bartlett', 'welch', 'cosine' (default: 'gauss')
%        n    number of samples                      (default: 65)
%    yedge    edge values                            (default: 0.001)
%
% Outputs:
%       wl    vector of wavelengths (nm)
%        y    vector of values
%       wn    vector of wavenumbers (cm^-1)
%
% Example:
%
%    [wl, y, wn] = srf_generate(380.10, 9.9);
%    plot(wl, y);
%    plot(wn, y);
%
% The latest version should be available via:
% wget \
% http://apex-esa.cvs.sourceforge.net/viewvc/*checkout*/apex-esa/tools/srf_generate.m
%
% inspired by Armin Gnter's freely distributable gauss.m
% http://ftp.physik.uni-greifswald.de/~guenter/matlab/signalproc/gauss.m
% and http://mathworld.wolfram.com/FullWidthatHalfMaximum.html
%
    if nargin < 5 yedge = 0.001;   end
    if nargin < 4     n = 65;      end
    if nargin < 3 shape = 'gauss'; end
    if nargin < 2 error('Give at least band center and fwhm.'); end

    if n < 3 error('Need at least 3 samples'); end

    n2 = floor((n-1)/2); 			% mid point
    x = [-n2:n2+mod(n-1,2)]';			% calculate sample indexes

    switch lower(shape)
    case 'gauss'
        sigma = sqrt(-n2.^2/2/log(yedge));	% calculate sigma
        nfwhm = 2 * sqrt(2 * log(2)) * sigma;	% normlized fwhm 
        y = exp(-x.^2/2/sigma^2);		% gaussian
    case 'bartlett'
        fudge = n2 * (1+yedge+(yedge/100));	
        nfwhm = n2 + (yedge * fudge);		% normlized fwhm 
        y = 1 - (abs(x) / nfwhm);		% bartlett
        y(find(y<0))=realmin;			% in case of odd samples
    case 'welch'
        fudge = n2 * (.5+(yedge*.378));
        alph = n2 + (yedge * fudge);		% calculate alpha
        nfwhm = sqrt(2) * alph;			% normalized fwhm 
        y = 1 - ((x .^ 2) / (alph ^ 2));	% welch
        y(find(y<0))=realmin;			% in case of odd samples
    case 'cosine'
        fudge = n2 * (.6365+(yedge*.5));
        alph = n2 + (yedge * fudge);		% calculate alpha
        nfwhm = (4/3) * alph;			% normalized fwhm 
        y = cos((pi .* x) / (2 * alph));	% cos
        y(find(y<0))=realmin;			% in case of odd samples
    case 'box'
        y = ones(size(x));
        y(1) = realmin;
        y(2*n2+1:end) = realmin;
        nfwhm = trapz(y);
    otherwise
        error('Shape must be one of: gauss, bartlett, welch, cosine');
    end

    delta = fwhm / nfwhm;			% determine sample delta

    wl = (ones(n, 1) .* center) + (delta .* x);	% wavelengths (nm)
    wn = 1e7 ./ wl;				% wavenumbers (cm^-1)
end