function asfs_matlab2nc_datareports(fname,begDir)
disp('asfs_matlab2nc_datareports ver. 1, 2018 August 27'); close all;

% C. Cox
% 8/24/2019
% 
%
% PURPOSE:
%
% This code reads files generated by the ASFS stations and creates a report
% on the contents of the files, data quality, plots, warnings, etc. The 
% program begins with "fname", which is the filename of the slow data file 
% you want a report on; 'fast' and 'met' data are also read and analyzed.
% The report is an html file that embeds stanalone figures.
%
% asfs_matlab2nc_datareports(fname)
% fname = filename of slow data in quotes, e.g., 'ASFS_30_slow_07202019_0500.dat' 
%
% Requires 
%   - Matlab R2019a (tested on 9.6.0.1072779)
%   - Matlab Statistics and Machine Learning Toolbox  
%
%
%                  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% CONTENTS:
%
% (1) Read the files
%       "slow", "met" and "fast" .dat files are read using readtable
%       - calls subfunction "cr1000x_reader"
%
% (2) Wite netCDF files
%       If flag "makenc" = 1, converts the slow/fast/met .dat files to nc
%       - calls subfunction "cdfPutlocal"
% 
% (3) Write report
%       Collects summary statistics and performs a series of simple tests
%       on instrument data quality, housekeeping, system health and file
%       contents. The results are written to an html file.
%       - calls subfunctions "colorizedstring", "despik_pct",
%         "solarPosition", and "hywex_rhi"
%
% (4) Figures
%       A series of .png files are created as quicklooks. These are
%       stand-alone figure files, but are also embedded into the html file.
%       The figures are not part of the html file, but are displayed within
%       the html file when the file is opened. The html code has paths that
%       point to the location of the figures. Thus, if the figures are
%       moved it won't fine them!
%
% (end) Close out
%       Writes a closing message
%
% (func) Local functions
%       These are the subfunctions that are listed in (1)-(4).
%
%
% 
% INSTRUCTIONS:
%
% You can set this up to run on a cron if you like. Perhaps with a bash
% script containing the command line statements to execute the script. This
% works on MacOSX where pathtomfile is the path to the directory where 
% asfs_matlab2nc_datareports.m is stored. You can add this path
% permanently (and then no need for addpath(genpath('pathtomfile')); in the
% statement below) by running addpath(genpath('pathtomfile')); and thens
% savepath from a matlab session once or adding the line to the startup.m
% file.
%
% matlab -nosplash -nodesktop -r "addpath(genpath('pathtomfile'));asfs_matlab2nc_datareports('ASFS_30_slow_07202019_0500.dat');exit"
%
% on a Windows maching it probably looks like this:
%
% matlab.exe -nosplash -nodesktop -r "addpath(genpath('pathtomfile'));asfs_matlab2nc_datareports('ASFS_30_slow_07202019_0500.dat');exit"
% or
% "C:\<a long path here>\matlab.exe" -nodisplay -nosplash -nodesktop -r "run('C:\<a long path here>\asfs_matlab2nc_datareports.m');"
%
%
% NOTES:
%
% - Matlab generates image files from figures in a manner similar to a screenshot so this won't ever truly run in the background. 
%
% - A day of data from one station takes about 36 sec to complete.
%
% - Each station must be run separately.
%
% - Self-contained. All ancillary functions not part the required software above are found in (func) Local Functions.
%
%
%
% ~~~~~~~~~~~~~~~   SETTINGS. You need to fill this out. ~~~~~~~~~~~~~~~~~~

%fname         = 'ASFS_30_slow_07202019_0500.dat'                                 ; % this is the filename of the slow data table file that you want a report on. met and fast data table files with the same date stamp are automatically.
%begDir        = '/Volumes/arguk/MOSAiC/MOSAiC_Documentation/Scripts/Quicklooks/' ; % this is the path to where fname is located
reDir         = begDir                                                           ; % ...or wherever you want the report to be stored
ncDir         = begDir                                                           ; % ...or wherever you want the netcdf files to be stored
splotDir      = begDir                                                           ; % ...or wherever you want the plots to be stored. 
fplotDir      = ''                                                               ; % where the html should find the plots. could be fplotDir = splotDir, but e.g., if set to '' makes it a relative path with figs sharing the directory with the html report
imtype        = 'png'                                                            ; % what's your favorite figure type? accepts 'png', 'jpeg', and 'tiff'
makenc        = 1                                                                ; % 1 = yes, convert the .dat files to .nc, anything else means skip that
metadata      = 'University of Colorado, MOSAiC. matthew.shupe@colorado.edu, PI' ; % this is a global attribute for the nc file. write whatever you like
nrec_exp      = 30                                                               ; % Probably should be 1440 (daily files): the number of records in the slow files (e.g., 30 = 30 min files). this is how frequently in minutes the Windows Task Manager executes the asfs_rename.py file
count_exp     = 12                                                               ; % 60s / 5s scan time. This shouldn't change.

% Configurable Physical Limits for color-coding and flagging data quality. 
%            ~~~~ Make adjustments if you like ~~~~
irt_targ_hi  = 5    ; % Maximum IRT surface brightness temperature [C]
irt_targ_lo  = -80  ; % Minimum IRT surface brightness temperature [C]
sr50d_hi     = 2.5  ; % Maximum distance from SR50 to surface [m]: make this value installation height + 0.5? 
sr50d_lo     = 1    ; % Minimum distance from SR50 to surface [m]: make this value installation height - 1? 
sr30_hi      = 1000 ; % SWD & SWU max [Wm^2]
sr30_lo      = -4   ; % SWD & SWU min [Wm^2]
ir20_hi      = 400  ; % LWD & LWU max [Wm^2] (FYI, 315 W/m2 ~ 0 C)
ir20_lo      = 50   ; % LWD & LWU min [Wm^2] (Note: LWD under clear skies can be as low as ~100 Wm2)
ir20_case_hi = 10   ; % Maximum IR20 instrument case temperature [C]
ir20_case_lo = -80  ; % Minimum IR20 instrument case temperature [C]
met_t_hi     = 10   ; % Maximum Vaisala air temperature [C]
met_t_lo     = -80  ; % Minimum Vaisala air temperature [C]
met_rh_hi    = 101  ; % Maximum Vaisala relative humidity [%]
met_rh_lo    = 5    ; % Minimum Vaisala relative humidity [%]
met_p_hi     = 1100 ; % Maximum Vaisala air pressure [hPa or ~mb]
met_p_lo     = 900  ; % Minimum Vaisala air pressure [hPa or ~mb]
lic_h2o_hi   = 25   ; % Maximum Licor h2o [mg/m3]
lic_h2o_lo   = 0    ; % Minimum Licor h2o [mg/m3]
lic_co2_hi   = 1000 ; % Maximum Licor co2 [g/m3]
lic_co2_lo   = 0    ; % Minimum Licor co2 [g/m3]

% Configurable limites for warnings found at the end of 
% Section 3 (Write Reports)
%            ~~~~ Make adjustments if you like ~~~~
licor_ice_limit = 50  ; % Licor signal strength that will be a threshold for ice. 50% is arbitraty and low.
dirty_licor     = 90  ; % Licor signal strength that will be a threshold for needing cleaning. 100% is clean. values below 100% are normal. licor provides little advice on limits. Based on testing experience in Boulder, 90% is probably ok. ASFS40 fell from ~99% to 88% form April to July without differing from the other two systems in July.
sr50_ice_limit  = 0.1 ; % Distance in meters that is the limit for saying its iced based on the assumption that if the sr50 distance is very small there must be soemthing in the way close to the source. 10 cm is arbitrary.
albedo_limit    = 0.5 ; % How low of an albedo before you get worried? 0.75-0.85 is expected for snow., 0.5-0.7 for bare ice. anything less than 0.5 and there is a problem: could be and ice SWU and ice-free SWD but differentially iced radiometers are probably the other way around (albedo >> 1) so probably open water/melt pond nearby
%hdg_limit       = 0.5; % (This is never going to work). Maximum rotation rate per day in degrees. Only tested if you have at least a day of data because the heading can be a bit erratic over shorter time scales

% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%% (1) Read the files
tic;
cd(begDir);

% Message local Matlab command prompt. disp is found throughout 
disp('Reading the files...');

% Set up the expected number of records (rows) in the fast and met files
nrec_exp_fast = nrec_exp*60*20;  % number of records in the fast files (30 min x 60 s/min x 20 Hz)
nrec_exp_met  = nrec_exp*12;     % number of records in the met files (30 min x 60 s/min * 5/60 Hz)

% Read the 'slow' data table.
% cr1000x_reader relies on matlab function readtable, which nicely handles
% unexpected data types that occasionally appear in Campbell files, most
% notably "NAN", which importdata cannot interpret.
%
% hdr1 = the first Campbell file header, which is metadata (logger serial number, compile results, etc)
% hdr2 = the column header for the variable names
% hdr3 = the column header for the variable units
% hdr4 = (if exists) the column header for the variable computation (e.g., sample, average, or standard devation).

try
    [hdr1,hdr2,hdr3,hdr4,file,nvars] = cr1000x_reader(fname,begDir);
    disp(['File ',fname,' successfully read']);
    readmessage{1} = {['File ',fname,' successfully read']};
    abort_flag = 0;
catch
    readmessage{1} = {['Read Failure: ',fname,' not found. Abort.']};
    disp('Read Failure: No slow file found.');
    abort_flag = 1; % if we can't find the file, we will print some feedback to the html and exit down a little farther
end

% Now load 'fast' and 'met' files in the same way as the slow file
fnamemet = fname;
fnamemet = strrep(fnamemet,'slow','met');
try
    [hdr1_m,hdr2_m,hdr3_m,hdr4_m,file_m,nvars_m] = cr1000x_reader(fnamemet,begDir);
    disp(['File ',fname,' successfully read']);
    readmessage{2} = {['File ',fnamemet,' successfully read']};
    abort_flag_met = 0;
catch
    readmessage{2} = {'Read Failure: No met file found.'};
    disp('Read Failure: No met file found.');
    abort_flag_met = 1;
end

fnamefast = fname;
fnamefast = strrep(fnamefast,'slow','fast');
try
    [hdr1_f,hdr2_f,hdr3_f,hdr4_f,file_f,nvars_f] = cr1000x_reader(fnamefast,begDir);
    disp(['File ',fname,' successfully read']);
    readmessage{3} = {['File ',fnamefast,' successfully read']};
    abort_flag_fast = 0;
catch
    readmessage{3} = {'Read Failure: No fast file found.'};
    disp('Read Failure: No fast file found.');
    abort_flag_fast = 1;
end

disp(['Finished reading files at ',num2str(round(toc*10)/10),'sec.']);

%% (2) Write a netCDFs

% if the flag is set to make nc files and we successfully read the data, then proceed
if makenc == 1 && abort_flag == 0
    
    disp ('Writing netCDF files.');
    
    % only cycle through files you found!
    js = 1:3;
    js = js.*(~[abort_flag abort_flag_fast abort_flag_met]);
    js(js==0) = [];

    % now cycle through the files and convert to netcdf
    for j = js
        
        if j == 1
            fname_tmp = fname; 
            hdr1_tmp = hdr1;
            hdr2_tmp = hdr2;
            hdr3_tmp = hdr3;
            hdr4_tmp = hdr4;
            file_tmp = file;
            nvar_tmp = nvars;
        elseif j == 2
            fname_tmp = fnamefast;
            hdr1_tmp = hdr1_f;
            hdr2_tmp = hdr2_f;
            hdr3_tmp = hdr3_f;
            hdr4_tmp = hdr4_f;
            file_tmp = file_f;
            nvar_tmp = nvars_f;
        elseif j == 3
            fname_tmp = fnamemet;
            hdr1_tmp = hdr1_m;
            hdr2_tmp = hdr2_m;
            hdr3_tmp = hdr3_m;
            hdr4_tmp = hdr4_m;
            file_tmp = file_m;
            nvar_tmp = nvars_m;
        end
        
        
        try
            % new nc file name will be...
            ncfname = [fname_tmp(1:end-4),'.nc'];

            % open nc file pointer
            ncid = netcdf.create(ncfname,'NC_WRITE');

            % add some global attributes

              % date created
            Date_Created = num2str(date);
            varid = netcdf.getConstant('NC_GLOBAL');
            netcdf.putAtt(ncid,varid,'Date_Created',Date_Created);

              % file meta data
            Source = char(hdr1_tmp);
            netcdf.putAtt(ncid,varid,'Source',char(cellstr(Source)));  
            netcdf.putAtt(ncid,varid,'Contact',metadata);  


            % write some data

              %  this just extracts years, months, days etc from the time stamp and
              %  puts them in the netcdf
            [yyyy mm dd hh mi ss] = datevec(datenum(file_tmp.TIMESTAMP));
            theatt.long_name = 'year';
            cdfPutlocal(yyyy,'yyyy','double',theatt,length(yyyy),ncid);  
            theatt.long_name = 'month';
            cdfPutlocal(mm,'mo','NC_SHORT',theatt,length(mm),ncid);  
            theatt.long_name = 'day';
            cdfPutlocal(dd,'day','NC_SHORT',theatt,length(dd),ncid);  
            theatt.long_name = 'hour UTC';
            cdfPutlocal(hh,'hh','NC_SHORT',theatt,length(hh),ncid);  
            theatt.long_name = 'minute';
            cdfPutlocal(mi,'mi','NC_SHORT',theatt,length(mi),ncid);  
            theatt.long_name = 'second';
            cdfPutlocal(ss,'ss','NC_SHORT',theatt,length(ss),ncid);  

            % insert data into the file beginning with the 2nd variable
            % (i.e., skip TIMESTAMP) and iterating through all columns

            for h = 2:nvar_tmp
                data = table2array(file_tmp(:,h));
                theatt.long_name = char(hdr2_tmp(h));
                theatt.unit      = char(hdr3_tmp(h));
                cdfPutlocal(data,char(hdr2_tmp(h)),'NC_DOUBLE',theatt,length(data),ncid);  
            end

            % close out and print some messages
            
            netcdf.close(ncid);

            ncmessage{j} = {['netCDF successfully created: ',ncfname]};
            disp(['netCDF successfully created: ',ncfname,' at ',num2str(round(toc*10)/10),'sec.']);
        catch
            ncmessage{j} = {['netCDF write failed: ',ncfname]};
            disp(['netCDF write failed: ',ncfname,' at ',num2str(round(toc*10)/10),'sec.']);
        end
    
    end
    
else
    ncmessage = {'netCDF flag turned off. No .nc files were created.'};
    disp('netCDF flag turned off. Skipping.');
end

%% (3) Write Report

% Create the html report file name and open a file point (fileID)
rpfname = [fname(1:end-4),'_datareport.html'];
fileID = fopen(rpfname,'wt');

% begin by collecting some information
if abort_flag == 0
    
    disp('Writing report...');
    
    % Learn some things about the data
    hdr1_split = strsplit(hdr1,',');
    ftype = char(hdr1_split(end)); % slow, met or fast
    npts = height(file)*width(file); % total number of values in the matrix (height x width) = (number rows x number cols)
    if abort_flag_fast == 0; npts_f = height(file_f)*width(file_f); end
    if abort_flag_met == 0; npts_m = height(file_m)*width(file_m); end

    % Matlab serial date!
    % dn = number of days since 1/1/0000
    %   this a standard matlab date format and is convenient because many
    %   native functions expect it.
    dn = datenum(file.TIMESTAMP);
    
    % Ephemeris calculation
    % I doubt this is very accuarate at low sun angles (use Reda's sun_position.m instead)
    % but good enough for what we need here
    latitude = file.gps_lat_deg_Avg + file.gps_lat_min_Avg./60;
    longitude = file.gps_lon_deg_Avg + file.gps_lon_min_Avg./60;
    altitude = file.gps_alt_Avg./1000;
    szas = NaN.*dn;
    for k = 1:length(latitude)
        location.latitude = latitude(k);
        location.longitude = longitude(k);
        location.altitude = altitude(k);
        [angles,projection] = solarPosition(dn(k),location.latitude,location.longitude,0,0,0);
        szas(k) = angles(1);
    end
else
    hdr1_split = [{''} {''} {''} {''} {''} {''}]; % hold places with empties if we are going to abort
end

% =========================================================================
% Header
%
%   Write the header of the html file including metadata about the file and
%   status of reads and writes.
%
%   The approach to write information to the html file is the print to file
%  (using fprintf) to the html file (fileID) a string that is written in
%  html. The subfunction colorizestring allows the html text to be colored.


fprintf(fileID,'%s', ...
    ['<b>=========================================================='] ...
    );
fprintf(fileID,'%s', ...
    ['<br/>Data report for ',fname,' (inc. slow/fast/met)<br/><br/>'] ...
    );
fprintf(fileID,'%s', ...
    ['Station ', char(hdr1_split(2)),'<br/><br/>'] ...
    );
fprintf(fileID,'%s', ...
    ['Logger Program Version ', char(hdr1_split(6)),'</b><br/><br/>'] ...
    );
for k = 1: length(readmessage)
    fprintf(fileID,'%s', ...
        [char(readmessage{k}),'</b><br/>'] ...
        );
end
for k = 1: length(ncmessage)
    fprintf(fileID,'%s', ...
        [char(ncmessage{k}),'</b><br/>'] ...
        );
end
fprintf(fileID,'%s','<br/>');
fprintf(fileID,'%s', ...
    ['Legend: <br/>'] ...
    );
fprintf(fileID,'%s', ...
    [colorizestring('black','Black = FYI'),'<br/>'] ...
    );
fprintf(fileID,'%s', ...
    [colorizestring('green','Green = Data within limits'),'<br/>'] ...
    );
fprintf(fileID,'%s', ...
    [colorizestring('orange','Orange = Caution, data approaching limits'),'<br/>'] ...
    );
fprintf(fileID,'%s', ...
    [colorizestring('red','Red = Data out of limits'),'<br/>'] ...
    );
fprintf(fileID,'%s', ...
    ['<b>==========================================================</b><br/><br/><br/>'] ...
    );

% ok, now that we have written and html header, we can exit if we failed to
% read the slow file.
if abort_flag == 1; return; end

% =========================================================================





% =========================================================================
% Positioning
%
%   Now we will write a section on the geographic location data for the
%   ASFS: lat/lon/altitude/heading. We will also report on the quality of
%   the gps data and on the station tilt recorded by the metek's
%   inclinometer.

fprintf(fileID,'%s', ...
    ['<b>Positioning:</b><br/><br/>'] ...
    );


fprintf(fileID,'%s', ...
    ['Station begin position ', num2str(file.gps_lat_deg_Avg(1)), '&deg ', num2str(file.gps_lat_min_Avg(1),'%.8f'),''' N, ' ...
                                num2str(file.gps_lon_deg_Avg(1)), '&deg ', num2str(file.gps_lon_min_Avg(1),'%.8f'),''' E<br/>'] ...
    );
fprintf(fileID,'%s', ...
    ['Station end position ',   num2str(file.gps_lat_deg_Avg(end)), '&deg ', num2str(file.gps_lat_min_Avg(end),'%.8f'),''' N, ' ...
                                num2str(file.gps_lon_deg_Avg(end)), '&deg ', num2str(file.gps_lon_min_Avg(end),'%.8f'),''' E<br/>'] ...
    );

fprintf(fileID,'%s', ...
    ['Station begin heading ',   num2str(file.gps_hdg_Avg(1)./100),'&deg<br/>'] ...
    );
fprintf(fileID,'%s', ...
    ['Station end heading ',   num2str(file.gps_hdg_Avg(1)./100),'&deg<br/>'] ...
    );

fprintf(fileID,'%s', ...
    ['Station begin altitude ',   num2str(file.gps_alt_Avg(1)./1000),'m<br/>'] ...
    );
fprintf(fileID,'%s', ...
    ['Station end altitude ',   num2str(file.gps_alt_Avg(end)./1000),'m<br/>'] ...
    );

if mode(file.gps_qc) == 0
    gpqc = '0 = invalid';
    tcolor = 'red';
elseif mode(file.gps_qc) == 1 
    gpqc = '1 = GPS fix (SPS)';
    tcolor = 'green';   
elseif mode(file.gps_qc) == 2
    gpqc = '2 = DGPS fix';
    tcolor = 'green';   
elseif mode(file.gps_qc) == 3
    gpqc = '3 = PPS fix';
    tcolor = 'orange';   
elseif mode(file.gps_qc) == 4
    gpqc = '4 = Real Time Kinematic';
    tcolor = 'orange';   
elseif mode(file.gps_qc) == 5
    gpqc = '5 = Float RTK';
    tcolor = 'orange';   
elseif mode(file.gps_qc) == 6
    gpqc = '6 = estimated (dead reckoning) (2.3 feature)';
    tcolor = 'orange';   
elseif mode(file.gps_qc) == 7
    gpqc = '7 = Manual input mode';
    tcolor = 'orange';   
elseif mode(file.gps_qc) == 7
    gpqc = '8 = Simulation mode';
    tcolor = 'orange';   
end
fprintf(fileID,'%s', ...
    colorizestring(tcolor,['Modal value for GPS fix = ',num2str(mode(file.gps_qc)),' with a minimum of ',num2str(min(file.gps_nsat_Avg)),' satellites.<br/>']) ...
    );

% inclinometer
maxIncl = max([nanmean(file.metek_InclX_Avg) nanmean(file.metek_InclY_Avg)]);

if maxIncl < 1
    tcolor = 'green';
elseif 1 <= maxIncl <= 5
    tcolor = 'orange';   
elseif maxIncl > 5
    tcolor = 'red';  
end


fprintf(fileID,'%s', ...
    colorizestring(tcolor,['Maximum (X or Y) value observed of the Metek Inclinometer was ',num2str(maxIncl),'&deg.<br/><br/>']) ...
    );

% =========================================================================





% =========================================================================
% File Contents
%
%   Basic counts of the contents of the files. Does the number of rows
%   match the expected number? How many missing values and were can they
%   be found?

fprintf(fileID,'%s', ...
    ['<b>File Contents:</b><br/><br/>'] ...
    );


% How many records?
if height(file) < nrec_exp
    tcolor = 'red';
else
    tcolor = 'green';
end
fprintf(fileID,'%s', ...
    ['There are ',num2str(npts),' observations in the slow file: '] ...
    );
fprintf(fileID,'%s', ...
    colorizestring(tcolor,['Number records (rows): ',num2str(height(file)),'. Expected: ',num2str(nrec_exp),'<br/>']) ...
    );
% Number of missing data
nmd  = sum(ismissing(file),'all');
[r c] = find(ismissing(file));
c = unique(c);
nanvars = [];
for k = 1:length(c)
   nanvars = [nanvars,char(hdr2(c(k))),'   ']; 
end

if nmd > 0
    fprintf(fileID,'%s', ...
        colorizestring('orange',['Found ',num2str(nmd),' NaNs (',num2str(round(nmd/npts*10000)/100),'%) in the following Variables: ',nanvars,'<br/><br/>']) ...
        );
else
    fprintf(fileID,'%s', ...
        colorizestring('green',['No NaNs were found.<br/><br/>']) ...
        );
end


if abort_flag_met == 0
    if height(file_m) < nrec_exp_met
        tcolor = 'red';
    else
        tcolor = 'green';
    end
    fprintf(fileID,'%s', ...
        ['There are ',num2str(npts_m),' observations in the met file: '] ...
        );
    fprintf(fileID,'%s', ...
        colorizestring(tcolor,['Number records (rows): ',num2str(height(file_m)),'. Expected: ',num2str(nrec_exp_met),'<br/>']) ...
        );
    % Number of missing data
    nmd  = sum(ismissing(file_m),'all');
    [r c] = find(ismissing(file_m));
    c = unique(c);
    nanvars = [];
    for k = 1:length(c)
       nanvars = [nanvars,char(hdr2_m(c(k))),'   ']; 
    end

    if nmd > 0
        fprintf(fileID,'%s', ...
            colorizestring('orange',['Found ',num2str(nmd),' NaNs (',num2str(round(nmd/npts_m*10000)/100),'%) in the following Variables: ',nanvars,'<br/><br/>']) ...
            );
    else
        fprintf(fileID,'%s', ...
            colorizestring('green',['No NaNs were found.<br/><br/>']) ...
            );
    end
end

if abort_flag_fast == 0
    if height(file_f) < nrec_exp_fast
        tcolor = 'red';
    else
        tcolor = 'green';
    end
    fprintf(fileID,'%s', ...
        ['There are ',num2str(npts_f),' observations in the fast file: '] ...
        );
    fprintf(fileID,'%s', ...
        colorizestring(tcolor,['Number records (rows): ',num2str(height(file_f)),'. Expected: ',num2str(nrec_exp_fast),'<br/>']) ...
        );
    % Number of missing data
    nmd  = sum(ismissing(file_f),'all');
    [r c] = find(ismissing(file_f));
    c = unique(c);
    nanvars = [];
    for k = 1:length(c)
       nanvars = [nanvars,char(hdr2_f(c(k))),'   ']; 
    end

    if nmd > 0
        fprintf(fileID,'%s', ...
            colorizestring('orange',['Found ',num2str(nmd),' NaNs (',num2str(round(nmd/npts_f*10000)/100),'%) in the following Variables: ',nanvars,'<br/><br/>']) ...
            );
    else
        fprintf(fileID,'%s', ...
            colorizestring('green',['No NaNs were found.<br/><br/>']) ...
            );
    end
end

% =========================================================================





% =========================================================================
% Housekeeping
%
%   Temperatures of box/electronics, power information, etc.

fprintf(fileID,'%s', ...
    ['<b>Housekeeping:</b><br/><br/>'] ...
    );

fprintf(fileID,'%s',['Temperatures:<br/>']);

% Panel temeperature. This is the temperature of the terminal block on the
% CR1000X and is a reasonable estimate of the interface box temperature. We
% report the max and min temps in the file.
ptmpmax = max(file.PTemp_Avg);
ptmpmin = min(file.PTemp_Avg);

% dx/dy = change in panel temperature / day. We want to know the maximum
% absolute value of this metric. If it gets too large, the change in
% temperature may exceed the rate that the logger calibrates its voltage
% measurement, which temporarily increases uncertainty in IR20 and flux
% plates.
dxdy = abs(diff(file.PTemp_Avg)./diff(dn)); % C/day
dxdy(isinf(dxdy)) = NaN; % no change
maxdxdy=max(dxdy)/1440; % C/min

if ptmpmax > 50
    tcolor = 'red';
elseif (ptmpmax > 0) && (ptmpmax < 50)
    tcolor = 'orange';
else 
    tcolor = 'green';
end
fprintf(fileID,'%s', ...
    colorizestring(tcolor,['The maximum temperature of the logger panel was ',num2str(ptmpmax),' C<br/>']) ...
    );
if ptmpmin < -40
    tcolor = 'red'; 
elseif (ptmpmax < -25) && (ptmpmax > -40)
    tcolor = 'orange';
else 
    tcolor = 'green';
end
fprintf(fileID,'%s', ...
    colorizestring(tcolor,['The minimum temperature of the logger panel was ',num2str(ptmpmin),' C<br/>']) ...
    );
if maxdxdy > 1
    tcolor = 'red';
elseif (maxdxdy > 0.1) && (maxdxdy < 1)
    tcolor = 'orange';
else 
    tcolor = 'green';
end
fprintf(fileID,'%s', ...
    colorizestring(tcolor,['The maximum abs(temperature change) was ',num2str(maxdxdy),' C/min<br/>']) ...
    );
efoytintmax = max(file.efoy_Tint_Avg);
if efoytintmax > 40
    tcolor = 'red';
elseif (efoytintmax >= 0) && (efoytintmax < 40)
    tcolor = 'orange';
else 
    tcolor = 'green';
end
fprintf(fileID,'%s', ...
    colorizestring(tcolor,['Maximum EFOY internal temperature: ',num2str(efoytintmax),' C<br/>']) ...
    );

efoytintmin = min(file.efoy_Tint_Avg);
if efoytintmin < -40
    tcolor = 'red';
elseif (efoytintmin >= -40) && (efoytintmin < -30)
    tcolor = 'orange';
else 
    tcolor = 'green';
end
fprintf(fileID,'%s', ...
    colorizestring(tcolor,['Minimum EFOY internal temperature: ',num2str(efoytintmin),' C<br/>']) ...
    );

fprintf(fileID,'%s','<br/>');


fprintf(fileID,'%s',['Power:<br/>']);

% Battery voltage: actually the voltage supply to the logger. should be 12V
bvm = nanmean(file.batt_volt_Avg);
if (bvm > 12.5) || (bvm < 11.5)
    tcolor = 'red';
else 
    tcolor = 'green';
end
fprintf(fileID,'%s', ...
    colorizestring(tcolor,['Voltage supplied to the logger was ',num2str(bvm),' V<br/>']) ...
    );

evm = nanmean(file.efoy_Ubat_Avg);
if (evm < 12.3*2) || (evm > 14.2*2)
    tcolor = 'red';
    fprintf(fileID,'%s', ...
        colorizestring(tcolor,['EFOY battery voltages are out of range: ',num2str(evm),' V<br/>']) ...
        );
else 
    tcolor = 'green';
    fprintf(fileID,'%s', ...
        colorizestring(tcolor,['EFOY battery voltages are in range: ',num2str(evm),' V<br/>']) ...
        );
end

fprintf(fileID,'%s','<br/>');


fprintf(fileID,'%s',['EFOY:<br/>']);

efoyfuell = nanmean(file.efoy_FuellSt_Avg);
if efoyfuell < 10
    tcolor = 'red';
elseif (efoyfuell >= 10) && (efoyfuell < 20)
    tcolor = 'orange';
else 
    tcolor = 'green';
end
fprintf(fileID,'%s', ...
    colorizestring(tcolor,['EFOY remaining fuel estimate: ',num2str(efoyfuell),'%<br/>']) ...
    );
efoyerror = file.efoy_Error_Max;
if any(efoyerror) ~= 0
    tcolor = 'red';
    ee = unique(efoyerror); ee(ee == 0) = [];
    eestr = [];
    for k = 1:length(ee)
       eestr = [eestr num2str(ee(k)) '  '];
    end
    fprintf(fileID,'%s', ...
        colorizestring(tcolor,['The EFOY had errors: ',num2str(eestr),'<br/>']) ...
        );
else 
    tcolor = 'green';
    fprintf(fileID,'%s', ...
        colorizestring(tcolor,['There were no EFOY errors.<br/>']) ...
        );
end

fprintf(fileID,'%s','<br/>');



fprintf(fileID,'%s',['Data Logging:<br/>']);
% Counts
if strcmp(ftype,'slow') 
    skpscns = count_exp*nrec_exp - sum(file.counts_main_Tot);
    if skpscns >= count_exp*nrec_exp.*0.01
        tcolor = 'red';
    elseif (skpscns > 0) && (skpscns < count_exp*nrec_exp.*0.01)
        tcolor = 'orange';
    else 
        tcolor = 'green';  
    end
    fprintf(fileID,'%s', ...
        colorizestring(tcolor,['There were ',num2str(skpscns),' Skipped Scans<br/>']) ...
        );  
    
    
    maxst  = max(file.call_time_mainscan_Max);
    if maxst > 5000
        tcolor = 'red';
    elseif (maxst > 4000) && (maxst < 5000)
        tcolor = 'orange';
    else 
        tcolor = 'green';  
    end   
    fprintf(fileID,'%s', ...
        colorizestring(tcolor,['The maximum scan time was ',num2str(maxst),' mSec<br/>']) ...
        );  
    stmean = nanmean(file.call_time_mainscan_Max);
    if stmean > 5000
        tcolor = 'red';
    elseif (stmean > 4000) && (stmean < 5000)
        tcolor = 'orange';
    else 
        tcolor = 'green';  
    end   
    fprintf(fileID,'%s', ...
        colorizestring(tcolor,['The mean scan time was ',num2str(stmean),' mSec<br/>']) ...
        );  
    
end

fprintf(fileID,'%s','<br/>');
% =========================================================================





% =========================================================================
% Instrumentation
%
%   Each instrument is infdiviudally checked for (configurable) physically
%   possible limits, which are set by the user at the top of the code. The
%   amount of missing data is also reported. There are a few other quality 
%   indicators hard-coded in here too based on previous experience and 
%   manufacturer recommendations.

fprintf(fileID,'%s', ...
    '<b>Instrumentation:</b><br/><br/>' ...
    );

% IRT
irtb = file.apogee_body_T_Avg;
irtt = file.apogee_targ_T_Avg;
mv = length(find(isnan(irtt)));
bv = length(find(irtt > irt_targ_hi | irtt < irt_targ_lo));
fprintf(fileID,'%s',['IRT:<br/>']);
if mv == 0 && bv == 0
    fprintf(fileID,'%s',colorizestring('green','OK<br/><br/>'));
else
    fprintf(fileID,'%s', ...
        colorizestring('orange',['There are ',num2str(mv),' missing records and ',num2str(bv),' records out-of-range. <br/><br/>']) ...
        );
end

% SR50A
sr50d = file.sr50_dist_Avg;
sr50q = file.sr50_qc_Avg;
mv = length(find(isnan(sr50d)));
bv = length(find(sr50d > sr50d_hi | sr50d < sr50d_lo));

fprintf(fileID,'%s','SR50A:<br/>');
if mv == 0 && bv == 0 && nanmean(file.sr50_qc_Avg) >= 152 && nanmean(file.sr50_qc_Avg <= 210)
    fprintf(fileID,'%s',colorizestring('green','OK<br/><br/>'));
else
    if mv > 0
        tcolor_m = 'orange';
    elseif mv > nrec_exp.*0.1
        tcolor_m = 'red';
    else
        tcolor_m = 'green';
    end
    
    if bv > 0
        tcolor_b = 'orange';
    elseif mv > nrec_exp.*0.1
        tcolor_b = 'red';
    else
        tcolor_b = 'green';
     end
    
    if nanmean(file.sr50_qc_Avg) < 152 || nanmean(file.sr50_qc_Avg) > 300
        tcolor_q = 'red'; 
    elseif nanmean(file.sr50_qc_Avg) > 210 && nanmean(file.sr50_qc_Avg) <= 300 
        tcolor_q = 'orange';
    else
        tcolor_q = 'green';
    end
       
    fprintf(fileID,'%s',colorizestring(tcolor_m,['There are ',num2str(mv),' missing records and ']));
    fprintf(fileID,'%s',colorizestring(tcolor_b,['there are ',num2str(bv),' records out-of-range and ']));
    fprintf(fileID,'%s',colorizestring(tcolor_q,['the data quality is ',num2str(nanmean(file.sr50_qc_Avg)),'.<br/><br/>']));

end

% Vaisala
met_t = file.vaisala_T_Avg;
met_p = file.vaisala_P_Avg;
met_rh = file.vaisala_RH_Avg;
mv = length(find(isnan(met_t)));
bv = length(find(met_t > met_t_hi | met_t < met_t_lo));
fprintf(fileID,'%s','Vaisala:<br/>');
if mv == 0 && bv == 0
    fprintf(fileID,'%s',colorizestring('green','Temperature OK<br/>'));
else
    fprintf(fileID,'%s', ...
        colorizestring('orange',['There are ',num2str(mv),' missing temperature records and ',num2str(bv),' records out-of-range. <br/>']) ...
        );
end
mv = length(find(isnan(met_rh)));
bv = length(find(met_rh > met_rh_hi | met_rh < met_rh_lo));
if mv == 0 && bv == 0
    fprintf(fileID,'%s',colorizestring('green','Relative Humidity OK<br/>'));
else
    fprintf(fileID,'%s', ...
        colorizestring('orange',['There are ',num2str(mv),' missing humidity records and ',num2str(bv),' records out-of-range. <br/>']) ...
        );
end
mv = length(find(isnan(met_p)));
bv = length(find(met_p > met_p_hi | met_p < met_p_lo));
if mv == 0 && bv == 0
    fprintf(fileID,'%s',colorizestring('green','Pressure OK<br/>'));
else
    fprintf(fileID,'%s', ...
        colorizestring('orange',['There are ',num2str(mv),' missing pressure records and ',num2str(bv),' records out-of-range. <br/><br/>']) ...
        );
end


% SR30
mv = length(find(isnan(file.sr30_swd_Irr_Avg)));
bv = length(find(file.sr30_swd_Irr_Avg > sr30_hi | file.sr30_swd_Irr_Avg < sr30_lo));
fantach=nanmean(file.sr30_swd_fantach_Avg);
heater=nanmean(file.sr30_swd_heatA_Avg);
fprintf(fileID,'%s','SR30 SWD (upward facing):<br/>');
if mv == 0 && bv == 0 && fantach > 50 && heater > 50
    fprintf(fileID,'%s',colorizestring('green','OK<br/>'));
else
    fprintf(fileID,'%s', ...
        colorizestring('orange',['There are ',num2str(mv),' missing records and ',num2str(bv),' records out-of-range. <br/>']) ...
        );
end
if fantach > 50 
    fprintf(fileID,'%s', ...
        colorizestring('green',['The SWD fan is on (',num2str(fantach),' Hz).<br/>']) ...
        );
else
     fprintf(fileID,'%s', ...
        colorizestring('red',['The SWD fan is dead (',num2str(fantach),' Hz).<br/><br/>']) ...
        );   
end
if heater > 50
    fprintf(fileID,'%s', ...
        colorizestring('green',['The SWD heater is on (',num2str(heater),' mA).<br/><br/>']) ...
        );
else
     fprintf(fileID,'%s', ...
        colorizestring('red',['The SWD heater is off (',num2str(heater),' mA).<br/><br/>']) ...
        );   
end


mv = length(find(isnan(file.sr30_swu_Irr_Avg)));
bv = length(find(file.sr30_swu_Irr_Avg > sr30_hi | file.sr30_swu_Irr_Avg < sr30_lo));
fantach=nanmean(file.sr30_swu_fantach_Avg);
heater=nanmean(file.sr30_swu_heatA_Avg);
fprintf(fileID,'%s','SR30 SWU (downward facing):<br/>');
if mv == 0 && bv == 0 && fantach > 50 && heater > 50
    fprintf(fileID,'%s',colorizestring('green','OK<br/>'));
else
    fprintf(fileID,'%s', ...
        colorizestring('orange',['There are ',num2str(mv),' missing records and ',num2str(bv),' records out-of-range. <br/>']) ...
        );
end
if fantach > 50 
    fprintf(fileID,'%s', ...
        colorizestring('green',['The SWU fan is on (',num2str(fantach),' Hz).<br/>']) ...
        );
else
     fprintf(fileID,'%s', ...
        colorizestring('red',['The SWU fan is dead (',num2str(fantach),' Hz).<br/>']) ...
        );   
end
if heater > 50
    fprintf(fileID,'%s', ...
        colorizestring('green',['The SWU heater is on (',num2str(heater),' mA).<br/>']) ...
        );
else
     fprintf(fileID,'%s', ...
        colorizestring('red',['The SWU heater is off (',num2str(heater),' mA).<br/>']) ...
        );   
end

% Albedo
alb = nanmean(file.sr30_swu_IrrC_Avg) ./ nanmean(file.sr30_swd_IrrC_Avg);
nightind = find(szas > 80); 
alb(nightind) = NaN;
if length(nightind) == length(alb)
    fprintf(fileID,'%s',colorizestring('black',['Too dark to calculate albedo.<br/><br/>']));
else
    fprintf(fileID,'%s',colorizestring('black',['Mean Albedo: ',num2str(nanmean(alb)),'.<br/><br/>']));
end



% IR20
mv = length(find(isnan(file.ir20_lwd_Wm2_Avg)));
bv = length(find(file.ir20_lwd_Wm2_Avg > sr30_hi | file.ir20_lwd_Wm2_Avg < sr30_lo));
ir20case=nanmean(file.ir20_lwd_DegC_Avg);
fantach=nanmean(file.ir20_lwd_fan_Avg);
fprintf(fileID,'%s','IR20 LWD (upward facing):<br/>');
if mv == 0 && bv == 0 && fantach > 2000 && ir20case < ir20_case_hi && ir20case > ir20_case_lo
    fprintf(fileID,'%s',colorizestring('green','OK<br/>'));
else
    fprintf(fileID,'%s', ...
        colorizestring('orange',['There are ',num2str(mv),' missing records and ',num2str(bv),' records out-of-range. <br/>']) ...
        );
end

if ir20case < ir20_case_hi && ir20case > ir20_case_lo
    fprintf(fileID,'%s', ...
        colorizestring('green','The LWD case temperature is in range.<br/>') ...
        );
else
     fprintf(fileID,'%s', ...
        colorizestring('red','The LWD case temperature is out of range. This affects the flux. Check wires and resistor on CR1000X. Could also be caused by ground loop or unexpected voltage source.<br/>') ...
        );   
end

if fantach > 2000
    fprintf(fileID,'%s', ...
        colorizestring('green',['The LWD fan is on (',num2str(fantach),' mV).<br/>']) ...
        );
elseif fantach < 400
     fprintf(fileID,'%s', ...
        colorizestring('red',['The LWD fan is dead (',num2str(fantach),' mV).<br/>']) ...
        ); 
else
     fprintf(fileID,'%s', ...
        colorizestring('orange',['The LWD fan is struggling (',num2str(fantach),' mV). Possibly clogged with snow/ice.<br/>']) ...
        ); 
end

fprintf(fileID,'%s','<br/>');

mv = length(find(isnan(file.ir20_lwu_Wm2_Avg)));
bv = length(find(file.ir20_lwu_Wm2_Avg > sr30_hi | file.ir20_lwu_Wm2_Avg < sr30_lo));
ir20case=nanmean(file.ir20_lwu_DegC_Avg);
fantach=nanmean(file.ir20_lwu_fan_Avg);
fprintf(fileID,'%s','IR20 LWU (downward facing):<br/>');
if mv == 0 && bv == 0 && fantach > 2000 && ir20case < ir20_case_hi && ir20case > ir20_case_lo
    fprintf(fileID,'%s',colorizestring('green','OK<br/>'));
else
    fprintf(fileID,'%s', ...
        colorizestring('orange',['There are ',num2str(mv),' missing records and ',num2str(bv),' records out-of-range. <br/>']) ...
        );
end
if ir20case < ir20_case_hi && ir20case > ir20_case_lo 
    fprintf(fileID,'%s', ...
        colorizestring('green','The LWD case temperature is in range.<br/>') ...
        );
else
     fprintf(fileID,'%s', ...
        colorizestring('red','The LWD case temperature is out of range. This affects the flux. Check wires and resistor on CR1000X. Could also be caused by ground loop or unexpected voltage source.<br/>') ...
        );   
end
if fantach > 2000
    fprintf(fileID,'%s', ...
        colorizestring('green',['The LWD fan is on (',num2str(fantach),' mV).<br/><br/>']) ...
        );
elseif fantach < 400
     fprintf(fileID,'%s', ...
        colorizestring('red',['The LWD fan is dead (',num2str(fantach),' mV).<br/><br/>']) ...
        ); 
else
     fprintf(fileID,'%s', ...
        colorizestring('orange',['The LWD fan is struggling (',num2str(fantach),' mV). Possibly clogged with snow/ice.<br/><br/>']) ...
        ); 
end



% Licor
fprintf(fileID,'%s','Licor:<br/>');
if abort_flag_fast == 0
    
    licor_diag = (file_f.licor_diag_out);
    licor_diag(isnan(licor_diag)) = 0;
    licor_diag_bin = dec2bin(licor_diag);
    b1 = str2num(licor_diag_bin(:,1));
    b2 = str2num(licor_diag_bin(:,2));
    b3 = str2num(licor_diag_bin(:,3));
    b4 = str2num(licor_diag_bin(:,4));
    b5 = str2num(licor_diag_bin(:,5));
    b6 = str2num(licor_diag_bin(:,6));
    b7 = str2num(licor_diag_bin(:,7));
    b8 = str2num(licor_diag_bin(:,8));  
    ss = NaN.*licor_diag;
    pll = NaN.*licor_diag;
    detector_temp = NaN.*licor_diag;
    chopper_temp = NaN.*licor_diag;
    for k = 1:length(licor_diag)
        ss(k)            = 6.67*(b5(k)*2^3+b6(k)*2^2+b2(k)*2^1+b8(k)*2^0); % in percent, the signal strength. we will ignore and report the value 
        pll(k)           = b3(k); % PLL: 1 = good, 0 = bad
        detector_temp(k) = b2(k); % temp near set point = 1, not near set point = 0
        chopper_temp(k)  = b1(k); % temp near set point = 1, not near set point = 0
    end
    h2o_outliers = despik_pct(file_f.licor_h2o_out);
    h2o_outliers = h2o_outliers + length(find(isnan(file_f.licor_h2o_out)))/length(file_f.licor_h2o_out)*100;
    co2_outliers = despik_pct(file_f.licor_co2_out);
    co2_outliers = co2_outliers + length(find(isnan(file_f.licor_co2_out)))/length(file_f.licor_co2_out)*100;

    if h2o_outliers > 1 && h2o_outliers < 5
        tcolor = 'orange';
    elseif h2o_outliers >= 5
        tcolor = 'red';
    else 
        tcolor = 'green';
    end
    fprintf(fileID,'%s', ...
        colorizestring(tcolor,['H2O outliers & missing values: ',num2str(h2o_outliers),'%.<br/>']) ...
        );
    if co2_outliers > 1 && co2_outliers < 5
        tcolor = 'orange';
    elseif co2_outliers >= 5
        tcolor = 'red';
    else 
        tcolor = 'green';
    end
    fprintf(fileID,'%s', ...
        colorizestring(tcolor,['CO2 outliers & missing values: ',num2str(co2_outliers),'%.<br/>']) ...
        );


    pct_detector_pass = length(find(detector_temp == 1)) / length(detector_temp) * 100;
    if pct_detector_pass > 95 && pct_detector_pass < 99
        tcolor = 'orange';
    elseif pct_detector_pass <= 95
        tcolor = 'red';
    else 
        tcolor = 'green';
    end
    fprintf(fileID,'%s', ...
        colorizestring(tcolor,['Percent of detector temperatures that were good: ',num2str(pct_detector_pass),'%.<br/>']) ...
        );
    pct_chopper_pass = length(find(chopper_temp == 1)) / length(chopper_temp) * 100;
    if pct_chopper_pass > 95 && pct_chopper_pass < 99
        tcolor = 'orange';
    elseif pct_chopper_pass <= 95
        tcolor = 'red';
    else 
        tcolor = 'green';
    end
    fprintf(fileID,'%s', ...
        colorizestring(tcolor,['Percent of chopper temperatures that were good: ',num2str(pct_chopper_pass),'%.<br/>']) ...
        );
else
     fprintf(fileID,'%s', ...
        colorizestring('orange','Diagnostics Limited: No fast file available.<br/>') ...
        );
end

signal_str = nanmean(file.licor_co2_str_out_Avg);
if signal_str > 90 && signal_str < 99
    tcolor = 'orange';
elseif signal_str < 90
    tcolor = 'red';
else 
    tcolor = 'green';
end
fprintf(fileID,'%s', ...
    colorizestring(tcolor,['Licor Signal Strength: ',num2str(signal_str),'%.<br/>']) ...
    );

fprintf(fileID,'%s','<br/>');


% Metek
fprintf(fileID,'%s','Metek:<br/>');
fprintf(fileID,'%s',colorizestring('blue','Metek heater & pathstate disgnostics not yet implemented...<br/>'));
if abort_flag_fast == 0
    metek_diag = (file_f.metek_diag_out);
    metek_x_outliers = despik_pct(file_f.metek_x_out);
    metek_x_outliers = metek_x_outliers + length(find(isnan(file_f.metek_x_out)))/length(file_f.metek_x_out)*100;
    metek_y_outliers = despik_pct(file_f.metek_y_out);
    metek_y_outliers = metek_y_outliers + length(find(isnan(file_f.metek_y_out)))/length(file_f.metek_y_out)*100;
    metek_z_outliers = despik_pct(file_f.metek_z_out);
    metek_z_outliers = metek_z_outliers + length(find(isnan(file_f.metek_z_out)))/length(file_f.metek_z_out)*100;
    metek_T_outliers = despik_pct(file_f.metek_T_out);
    metek_T_outliers = metek_T_outliers + length(find(isnan(file_f.metek_T_out)))/length(file_f.metek_T_out)*100;
      
    if metek_x_outliers > 1 && metek_x_outliers < 5
        tcolor = 'orange';
    elseif metek_x_outliers >= 5
        tcolor = 'red';
    else 
        tcolor = 'green';
    end
    fprintf(fileID,'%s', ...
        colorizestring(tcolor,['Metek x outliers (based on Grachev despik.m) & missing values: ',num2str(metek_x_outliers),'%.<br/>']) ...
        );
    
    if metek_y_outliers > 1 && metek_y_outliers < 5
        tcolor = 'orange';
    elseif metek_y_outliers >= 5
        tcolor = 'red';
    else 
        tcolor = 'green';
    end
    fprintf(fileID,'%s', ...
        colorizestring(tcolor,['Metek y outliers (based on Grachev despik.m) & missing values: ',num2str(metek_y_outliers),'%.<br/>']) ...
        ); 
    
    if metek_z_outliers > 1 && metek_z_outliers < 5
        tcolor = 'orange';
    elseif metek_z_outliers >= 5
        tcolor = 'red';
    else 
        tcolor = 'green';
    end
    fprintf(fileID,'%s', ...
        colorizestring(tcolor,['Metek z outliers (based on Grachev despik.m) & missing values: ',num2str(metek_z_outliers),'%.<br/>']) ...
        );
    
    if metek_T_outliers > 1 && metek_T_outliers < 5
        tcolor = 'orange';
    elseif metek_T_outliers >= 5
        tcolor = 'red';
    else 
        tcolor = 'green';
    end
    fprintf(fileID,'%s', ...
        colorizestring(tcolor,['Metek T outliers (based on Grachev despik.m) & missing values: ',num2str(metek_T_outliers),'%.<br/>']) ...
        );
end

fprintf(fileID,'%s','<br/>');

% =========================================================================





% =========================================================================
% Warnings
%
%   Messages are provided when the data suggests the station may need
%   attention. Note that the thresholds are somewhat arbitrary 

fprintf(fileID,'%s', ...
    '<b>Station Warnings:</b><br/><br/>' ...
    );

warnflag = 0;

if nanmean(sr50d) > sr50d_hi
    fprintf(fileID,'%s',colorizestring('red','WARNING: SR50 distance is expectedly large. ASFS may be tilting or SR50 may be off level.<br/><br/>'));
    warnflag = 1;
end  
if nanmean(sr50d) < sr50d_lo
    fprintf(fileID,'%s',colorizestring('red','WARNING: SR50 distance is expectedly small. ASFS may be drifted in or SR50 may be persistently iced.<br/><br/>'));
    warnflag = 1;
end
if nanmean(alb) < albedo_limit
    fprintf(fileID,'%s',colorizestring('red','WARNING: Albedo <',num2str(albedo_limit),': Melt ponds or open water may be in FOV; else SWU iced more than SWD.<br/><br/>'));
    warnflag = 1;
end
if nanmean(signal_str) < dirty_licor && nanmean(signal_str) > licor_ice_limit
    fprintf(fileID,'%s',colorizestring('red','WARNING: Licor optics and radiometer domes may need cleaning of salt or light icing.<br/><br/>'));
    warnflag = 1;
end
if efoyfuell < 10
    fprintf(fileID,'%s',colorizestring('red','WARNING: ASFS is almost out of gas.<br/><br/>'));
    warnflag = 1;
end

% Trying to make a flag for changes in heading. Based on the Boulder data, this might need a long period to be confident.
% linear fit to heading
% a = polyfit(dn,file.gps_hdg_Avg./100,1); slp = a(1); % linear regress in deg/day
% if abs(slp) > hdg_limit && (dn(end)-dn(1)) > 0.96 % if you have at least a day of data. short durations will show a lot of variability
%     fprintf(fileID,'%s',colorizestring('red','WARNING: ASFS heading is changing fast!.<br/><br/>'));
%     warnflag = 1;
% end    


if warnflag == 0
    fprintf(fileID,'%s',colorizestring('green',['There are no station warnings to report.']));
end

fprintf(fileID,'%s','<br/><br/>');

% =========================================================================





% =========================================================================
% Icing indicators
%
%   Using data available to estimate the liklihood that ASFS is icing.

fprintf(fileID,'%s', ...
    ['<b>Icing Indicators:</b><br/><br/>'] ...
    );

% Relative Humidity w.r.t. ice
rhi = hywex_rhi(file.vaisala_RH_Avg,file.vaisala_T_Avg+273.15);

if length(find(rhi >= 100)) > 0 && length(find(rhi >= 100)) < length(rhi)
    fprintf(fileID,'%s',colorizestring('orange',['Possible Icing: RHi > 100% some of the time.<br/>']));
elseif length(find(rhi >= 100)) > 0 && length(find(rhi >= 100)) == length(rhi)
    fprintf(fileID,'%s',colorizestring('red',['Likely Icing: RHi persistently higher than 100%.<br/>']));
else 
    fprintf(fileID,'%s',colorizestring('green',['Icing unlikely: Air subsaturated w.r.t. ice.<br/>']));
end

% SR50
if nanmean(sr50d) < sr50_ice_limit
    fprintf(fileID,'%s',colorizestring('red',['SR50 distance is <',num2str(sr50_ice_limit.*10),' cm: Icing likely.<br/>']));
else
    fprintf(fileID,'%s',colorizestring('green',['SR50 not attenuated by ice.<br/>']));
end

% Optical icing detector
iced_times = find(file.licor_co2_str_out_Avg < licor_ice_limit);
iced_perc = length(iced_times)/length(file.licor_co2_str_out_Avg)*100;
if ~isempty(iced_times)
    fprintf(fileID,'%s',colorizestring('orange',['Optical Icing Flag: Licor optics obstructed ',num2str(round(iced_perc*100)/100),'% of the time. ASFS may be affected by icing. <br/>']));
else
    fprintf(fileID,'%s',colorizestring('green',['Licor optics are ice-free.<br/>']));
end  

disp(['Finished writing report at ',num2str(round(toc*10)/10),'sec.']);

% =========================================================================

%% (4) Figures

disp('Making figures...');

% Abort if we don't have a ong, jpeg or tiff selected
if sum(strcmp('png',imtype)+strcmp('jpeg',imtype)+strcmp('tiff',imtype)) == 0
    disp(['Please provide a valid image type for the figure files.']);
    return
end

fprintf(fileID,'%s','<br/><br/>');

fprintf(fileID,'%s', ...
    ['<b>Quicklooks:</b><br/><br/>'] ...
    );

cd(splotDir)

% XL = x-limits of the plots set to the begin and end times of the file.
XL = [dn(1) dn(end)];
% XT is the spacing of the x-ticks
if nrec_exp == 30 
    XT = floor(XL(1)):5/1440:ceil(XL(end));               % x-tick every 5 min for a 30 min file
elseif nrec_exp == 60
    XT = floor(XL(1)):10/1440:ceil(XL(end));              % x-tick every 10 min for a 60 min file
elseif nrec_exp == 1440
    XT = floor(XL(1)):1/24:ceil(XL(end));                 % x-tick every hour for a daily file
else
    XT = floor(XL(1)):round(length(dn)/10):ceil(XL(end)); % x-tick every 10th of the span for anything else
end

% plot sizes
plot_ht = 14;                        % plot width of the figure file in cm
plot_wd = 24;                        % plot height of the figure file in cm
html_ht = 500;                       % plot height in the html file. no idea what the units are, but 500 looks good.
html_wd = html_ht * plot_wd/plot_ht; % maintain the aspect ratio

f1 = [fname(1:7),'_',fname(end-16:end-4),'_wind_vel.',imtype];             % make a file name for the figure
v = sqrt(file.metek_x_Avg.^2+file.metek_y_Avg.^2);                         % horizontal wind velocity calc
clf;                                                                       % clear current figure
plot(dn,v,'.-','markersize',15); hold on;                                  % plot(x,y,'linestyle options')
ylim([0 25]); xlim(XL); grid on;                                           % x and y limits of the plot, overprint a grid at ticks
set(gca,'xtick',XT);                                                       % set the ticks of the current axis (gca = get current axis)
datetick('x','HH:MM','keepticks','keeplimits');                            % since x is a matlab serial date we can format it how we like, here hour:minute
xlabel('HH:MM UTC');                                                       % string for the x label
ylabel('Metek Hor Vel [m/s]');                                             % string for the y label
title([{'Metek Horizontal Wind Velocity'};{strrep(f1(1:end-4),'_','')}]);  % plot title, here a main title and subtitle
set(gca,'fontname','times','fontsize',12);                                 % fonts
set(gcf,'PaperPositionMode','manual');                                     % gcf = get current figure. this says we will defien the size ourselves
set(gcf,'PaperUnits','centimeters');                                       % and the units will be cm
set(gcf,'PaperPosition',[0 0 plot_wd plot_ht]);                            % and the size will be [left bottom with height]
print(['-d',imtype],'-r300',f1);                                           % save the figure with switch -dimtype (e.g., -dpng) at 300 dpi using the f1 filename
fprintf(fileID,'%s',['<img src="',fplotDir,f1,'" width="',num2str(html_wd),'" height="',num2str(html_ht),'">']); % embed this image into the html file

fprintf(fileID,'%s','<br/><br/>');
f2 = [fname(1:7),'_',fname(end-16:end-4),'_temps.',imtype];
clf;
plot(dn,file.vaisala_T_Avg,'.-','markersize',15); hold on;
plot(dn,file.apogee_targ_T_Avg,'.-','markersize',15);
plot(dn,((file.ir20_lwu_Wm2_Avg./(0.985*5.67*10^-8)).^0.25)-273.15,'.-','markersize',15);
xlim(XL); grid on;
ylabel('Temperature [C]');
set(gca,'xtick',XT);
datetick('x','HH:MM','keepticks','keeplimits');
xlabel('HH:MM UTC');
ylabel('Vaisala RH [%]');
legend('Vaisala Air Temp.','IRT Sfc Brightness Temp.','LWU Skin Temp (\epsilon=0.985)');
title([{'Temperature & Humidity'};{strrep(f2(1:end-4),'_','')}]);
set(gca,'fontname','times','fontsize',12);
set(gcf,'PaperPositionMode','manual');
set(gcf,'PaperUnits','centimeters');
set(gcf,'PaperPosition',[0 0 plot_wd plot_ht]);
print(['-d',imtype],'-r300',f2);
fprintf(fileID,'%s',['<img src="',fplotDir,f2,'" width="',num2str(html_wd),'" height="',num2str(html_ht),'">']);

fprintf(fileID,'%s','<br/><br/>');
f3 = [fname(1:7),'_',fname(end-16:end-4),'_humid.',imtype];
clf;
h1=plot(dn,file.vaisala_RH_Avg,'.-','markersize',15); hold on;
h2=plot(dn,rhi,'.-','markersize',15);
ylim([0 105]); xlim(XL); grid on;
set(gca,'xtick',XT);
datetick('x','HH:MM','keepticks','keeplimits');
xlabel('HH:MM UTC');
ylabel('RH [%]');
legend('RH_w','RH_i');
title([{'Relative Humidity'};{strrep(f3(1:end-4),'_','')}]);
set(gca,'fontname','times','fontsize',12);
set(gcf,'PaperPositionMode','manual');
set(gcf,'PaperUnits','centimeters');
set(gcf,'PaperPosition',[0 0 plot_wd plot_ht]);
print(['-d',imtype],'-r300',f3);
fprintf(fileID,'%s',['<img src="',fplotDir,f3,'" width="',num2str(html_wd),'" height="',num2str(html_ht),'">']);

fprintf(fileID,'%s','<br/><br/>');
f4 = [fname(1:7),'_',fname(end-16:end-4),'_longwave.',imtype];
clf;
plot(dn,file.ir20_lwd_Wm2_Avg,'.-','markersize',15); hold on;
plot(dn,file.ir20_lwu_Wm2_Avg,'.-','markersize',15);
xlim(XL); grid on;
set(gca,'xtick',XT);
datetick('x','HH:MM','keepticks','keeplimits');
xlabel('HH:MM UTC');
ylabel('Longwave Radiation [Wm^-^2]');
legend('LWD','LWU');
title([{'Longwave Radiation'};{strrep(f4(1:end-4),'_','')}]);
set(gca,'fontname','times','fontsize',12);
set(gcf,'PaperPositionMode','manual');
set(gcf,'PaperUnits','centimeters');
set(gcf,'PaperPosition',[0 0 plot_wd plot_ht]);
print(['-d',imtype],'-r300',f4);
fprintf(fileID,'%s',['<img src="',fplotDir,f4,'" width="',num2str(html_wd),'" height="',num2str(html_ht),'">']);

fprintf(fileID,'%s','<br/><br/>');
f5 = [fname(1:7),'_',fname(end-16:end-4),'_shortwave.',imtype];
clf;
plot(dn,file.sr30_swd_Irr_Avg,'.-','markersize',15); hold on;
plot(dn,file.sr30_swu_Irr_Avg,'.-','markersize',15);
xlim(XL); grid on;
set(gca,'xtick',XT);
datetick('x','HH:MM','keepticks','keeplimits');
xlabel('HH:MM UTC');
ylabel('Solar (Global) Radiation [Wm^-^2]');
legend('SWD','SWU');
title([{'Solar Radiation'};{strrep(f5(1:end-4),'_','')}]);
set(gca,'fontname','times','fontsize',12);
set(gcf,'PaperPositionMode','manual');
set(gcf,'PaperUnits','centimeters');
set(gcf,'PaperPosition',[0 0 plot_wd plot_ht]);
print(['-d',imtype],'-r300',f5);
fprintf(fileID,'%s',['<img src="',fplotDir,f5,'" width="',num2str(html_wd),'" height="',num2str(html_ht),'">']);

if any(~isnan(alb))
    fprintf(fileID,'%s','<br/><br/>');
    f6 = [fname(1:7),'_',fname(end-16:end-4),'_albedo.',imtype];
    clf;
    plot(dn,alb,'.-','markersize',15); hold on;
    ylim([0 1]); xlim(XL); grid on;
    set(gca,'xtick',XT);
    datetick('x','HH:MM','keepticks','keeplimits');
    xlabel('HH:MM UTC');
    ylabel('Albedo');
    title([{'Albedo'};{strrep(f6(1:end-4),'_','')}]);
    set(gca,'fontname','times','fontsize',12);
    set(gcf,'PaperPositionMode','manual');
    set(gcf,'PaperUnits','centimeters');
    set(gcf,'PaperPosition',[0 0 plot_wd plot_ht]);
    print(['-d',imtype],'-r300',f6);
    fprintf(fileID,'%s',['<img src="',fplotDir,f6,'" width="',num2str(html_wd),'" height="',num2str(html_ht),'">']);
end

fprintf(fileID,'%s','<br/><br/>');
f7 = [fname(1:7),'_',fname(end-16:end-4),'_fp.',imtype];
clf;
plot(dn,file.fp_A_Wm2_Avg,'.-','markersize',15); hold on;
plot(dn,file.fp_B_Wm2_Avg,'.-','markersize',15);
xlim(XL); grid on;
set(gca,'xtick',XT);
datetick('x','HH:MM','keepticks','keeplimits');
xlabel('HH:MM UTC');
ylabel('Conductive Flux [Wm^-^2]');
legend('Plate A','Plate B');
title([{'Flux Plates'};{strrep(f7(1:end-4),'_','')}]);
set(gca,'fontname','times','fontsize',12);
set(gcf,'PaperPositionMode','manual');
set(gcf,'PaperUnits','centimeters');
set(gcf,'PaperPosition',[0 0 plot_wd plot_ht]);
print(['-d',imtype],'-r300',f7);
fprintf(fileID,'%s',['<img src="',fplotDir,f7,'" width="',num2str(html_wd),'" height="',num2str(html_ht),'">']);

fprintf(fileID,'%s','<br/><br/>');
f8 = [fname(1:7),'_',fname(end-16:end-4),'_pres.',imtype];
clf;
plot(dn,file.vaisala_P_Avg,'.-','markersize',15); hold on;
xlim(XL); grid on;
set(gca,'xtick',XT);
datetick('x','HH:MM','keepticks','keeplimits');
xlabel('HH:MM UTC');
ylabel('Atm. Pres. [hPa]');
title([{'Air Pressure'};{strrep(f8(1:end-4),'_','')}]);
set(gca,'fontname','times','fontsize',12);
set(gcf,'PaperPositionMode','manual');
set(gcf,'PaperUnits','centimeters');
set(gcf,'PaperPosition',[0 0 plot_wd plot_ht]);
print(['-d',imtype],'-r300',f8);
fprintf(fileID,'%s',['<img src="',fplotDir,f8,'" width="',num2str(html_wd),'" height="',num2str(html_ht),'">']);

disp(['Figures finishd at ',num2str(round(toc*10)/10),'sec.']);

%% (end) Close out

% this doesn't do much. just a couple messages that say we are done now.
disp(['Report completed on ',datestr(now),'. Runtime: ',num2str(round(toc*10)/10),'sec.']);
fprintf(fileID,'%s',['<br/><br/><br/>Report completed on ',datestr(now),'. Runtime: ',num2str(round(toc*10)/10),'sec. Questions? christopher.j.cox@noaa.gov <br/>']);

end % match main function at line 1

%% (func) Local functions 


% =========================================================================
function [hdr1,hdr2,hdr3,hdr4,file,nvars] = cr1000x_reader(fname,begDir)
%
% C. Cox 
% August 2019
%
% cr1000x_reader reads .dat files form the loggers
%
% Relies on matlab function readtable, which nicely handles
% unexpected data types that occasionally appear in Campbell files, most
% notably "NAN", which importdata cannot interpret.
%
% Arguments: 
%
%   fname  = name of the file you want to read (string)
%   begDir = directory path to fname
%
% Returns:
%
%   hdr1  = the first Campbell file header, which is metadata (logger serial number, compile results, etc)
%   hdr2  = the column header for the variable names
%   hdr3  = the column header for the variable units
%	hdr4  = the column header for the variable computation (e.g., sample, average, or standard devation).
%   file  = a matlab table containing the data. Works like a structure (file.variable; e.g., file.TIMESTAMP or file.PTemp_Avg)
%   nvars = number of variables (i.e., columns)       

% grab the header using importdata and do some string manipulation because readtable is picky
tmp = importdata([begDir,fname]);
nhdrs = size(tmp.textdata,1)-size(tmp.data,1);

% these are the headers
hdr1 = tmp.textdata(1,1);
hdr2 = tmp.textdata(2,1);
hdr3 = tmp.textdata(3,1);
hdr4 = tmp.textdata(4,1);

% CR1000X information as char format sans double quotes
hdr1 = erase(char(hdr1),'"');

% Variable Names as char format sans double quotes
hdr2 = erase(char(hdr2),'"');
% string parse by comma delimeter. now we have an array of nvars strings
% with the variables name column headers
hdr2 = string(strsplit(hdr2,','));

% Variable Units. As with variable names, but for the units.
hdr3 = strsplit(char(hdr3),',');

% There might some variables without units, so find those and make them say 'N/A
for k = 1:length(hdr3)
    tmp = char(hdr3(k));
    tmp = strrep(tmp,'"','');
    if isempty(tmp); tmp = 'N/A'; end
    hdr3(k) = {tmp};
end
hdr3 = string(hdr3);

% Sample type (some files). As with variable names, but for the storage parameters 
hdr4 = strsplit(char(hdr4),',');
% There might some variables without a value, so find those and make them say 'N/A
for k = 1:length(hdr4)
    tmp = char(hdr4(k));
    tmp = strrep(tmp,'"','');
    if isempty(tmp); tmp = 'N/A'; end
    hdr4(k) = {tmp};
end
hdr4 = string(hdr4);

% number of variables/columns
nvars = length(hdr2);

% Set the options that will be passed to readtable. Some of these are not
% apparently available before Matlab 2019!
% Note that these strings do need to be " and not '. This is a new matlab
% thing. Strings ("..") are now distinct from chars ('..') and I think have
% some special powers more like object oriented programming than array
% programming. 

% Setup the Import Options
opts = delimitedTextImportOptions("NumVariables",nvars);
% Specify range and delimiter
opts.DataLines = [nhdrs+1, Inf];                                        % begin reading at the first data line and Inf(inite) denotes unknown file size
opts.Delimiter = ",";
% Specify column names and types
opts.VariableNames = hdr2;                                              % what are the columns I'm looking for?
opts.VariableTypes = cellstr(["datetime", repmat("double",1,nvars-1)]); % tell matlab that the first column is a datestring and that everything after that should be imported as a double
opts = setvaropts(opts, 1, "InputFormat", "yyyy-MM-dd HH:mm:ss");       % this is the format of the date used by Campbell
opts.ExtraColumnsRule = "ignore";                                       % ignore extra columns that were not specifid by hdr2 above
opts.EmptyLineRule = "read";                                            % for proper bean counting we should import NaNs if there is an empty line, but this probably won't happen 

% Import the data
% using readtable instead of importdata because importdata cannot deal with
%   some unexpected data types that appear (esp "NAN"). This is way faster 
%   than a line-by-line read too.
file = readtable([begDir,fname], opts);

end
% =========================================================================





% =========================================================================
function cdfPutlocal(Var,thename,prec,Att,Size,ncid)
%
% C. Cox 
% August 2019
%
% Inserting data into a netCDF requires a few steps such as defining dims 
% attributes and such. It is cleaner to put those steps in a function. 
%

% get a dimension id after telling the netcdf function the size of the data
dimid = [];
for hh = 1:length(Size)
    dimid = [dimid; netcdf.defDim(ncid,[thename,'d',num2str(hh)],Size(hh))];
end

% get a variable id after giving the netcdf function the dimid, precision and variable name
varid = netcdf.defVar(ncid,thename,prec,dimid);

% get out of (the default) define mode
netcdf.endDef(ncid);

% insert the variable in the the space allocated in the previous steps
netcdf.putVar(ncid,varid,Var);

% go back to define mode and add the attributes (if Att exists). Att is a
% structure so it is Att.att1, Att.att2, Att.att3, etc. where att# is the
% name of the attribute and it is equal to a string containing the
% metadata.
netcdf.reDef(ncid);
if ~isempty(Att)
    Attv = struct2cell(Att);
    Attn = fieldnames(Att);
    for hh = 1:length(Attn)
        netcdf.putAtt(ncid,varid,char(Attn(hh)),char(Attv(hh)));
    end
end

end
% =========================================================================





% =========================================================================
function colorizedstring = colorizestring(color, stringtocolorize)
%
% Some clever blogger showed how to print html text to a file with color
% https://blogs.mathworks.com/community/2009/05/25/colorizing-text-output/
%
% colorizestring wraps the given string in html that colors that text.

colorizedstring = ['<font color="', color, '">', stringtocolorize, '</font>'];

end
% =========================================================================





% =========================================================================
function pct=despik_pct(uz)
%
% A. Grachev/ C. Fairall
%
% Adapted by C. Cox
% August 2019 
% 
% This is the despiking code used by Fairall's group in EC calculations.
% I've lightly edited it and set it to only return the percent of flagged
% spikes. It is being used to identify data that would be excluded from an
% EC calculation.

[np,a]=size(uz);

uu2=sort(uz,1);
uu=uu2(1:np,1);
mp=floor(np/2);
mu=uu(mp);
sp=floor(.84*np);
sm=floor(.16*np);
sig=(uu(sp)-uu(sm))/2;
dsig=max(4*sig,.5);
im=1;
while abs(mu-uu(im))>dsig;
   im=im+1;
end;
ip=np;
while abs(uu(ip)-mu)>dsig;
   ip=ip-1;
end;
pct=(im+np-ip)/np*100;

end
% =========================================================================





% =========================================================================
function wvout = hywex_rhi(wvin,Tamb_K)
%
% P. Rowe
% May 1, 2002
%
% Adapted by C. Cox
% August 2019
%
% Calculated relative humidity with respect to ice following Hyland and
% Wexlar. This code block is lifted from humidRS80.m by P. Rowe.

rh_w = wvin;

c0 = -2.9912729e3;   c1 = -6.0170128e3;   c2 = 1.887643854e1;
c3 = -2.8354721e-2;   c4 = 1.7838301e-5;   c5 = -8.4150417e-10;
c6 = 4.4412543e-13;   D = 2.858487;
term = (c0*Tamb_K.^(-2)) + (c1*Tamb_K.^(-1)) + (c2*Tamb_K.^0) + ...
(c3*Tamb_K.^1) + (c4*Tamb_K.^2)+ (c5*Tamb_K.^3) + (c6*Tamb_K.^4);

Psat_w = exp(term + D*log(Tamb_K))  ; 	%Pa

Ph2o = 0.01 * rh_w .* Psat_w ;

td = Tamb_K;
c0 = -5.6745359e3 ;   c1 = 6.3925247;   c2 = -9.6778430e-3;
c3 = 6.2215701e-7 ;   c4 = 2.0747825e-9;   c5 = -9.4840240e-13;
D = 4.1635019;
% %
term = (c0*td.^(-1)) + (c1*td.^(0)) + (c2*td.^1) + ...
(c3*td.^2) + (c4*td.^3)+(c5*td.^4);
Psat_i = exp(term + (D*log(td)))  ; %Pa

rh_i = 100*Ph2o./Psat_i ;

wvout = rh_i ;

end
% =========================================================================





% =========================================================================
function [angles,projection] = solarPosition(datetime,latitude,longitude, ...
                                             time_zone,rotation,dst)
%
% Note by C. Cox August 2019: This code was downloaded from the Mathworks 
% File Exchange and chosen mostly because it is very compact and allows for
% a quick ephemeris calculation internally within this script. I tested it
% against Reda's sun_position program and it was close enough for what we
% need here. May not be great though and makes no atmospheric corrections.
                                         
%SOLARPOSITION Calculate solar position using most basic algorithm
%   This is the most basic algorithm. It is documented in Seinfeld &
%   Pandis, Duffie & Beckman and Wikipedia.
%
% [ANGLES,PROJECTION] = SOLARPOSITION(DATE,TIME,LATITUDE,LONGITUDE,TIME_ZONE)
% returns ZENITH & AZIMUTH for all DATE & TIME pairs at LATITUDE, LONGITUDE.
% ANGLES = [ZENITH,AZIMUTH] and PROJECTION = [PHI_X, PHI_Y]
% PHI_X is projection on x-z plane & PHI_Y is projection on y-z plane.
% DATETIME can be string, vector [YEAR, MONTH, DAY, HOURS, MINUTES, SECONDS],
%   cellstring or matrix N x [YEAR, MONTH, DAY, HOURS, MINUTES, SECONDS] for N
%   times.
% LATITUDE [degrees] and LONGITUDE [degrees] are the coordinates of the site.
% TIME_ZONE [hours] of the site.
% ROTATION [degrees] clockwise rotation of system relative to north.
% DST [logical] flag for daylight savings time, typ. from March to November
%   in the northern hemisphere.
%
% References:
% http://en.wikipedia.org/wiki/Solar_azimuth_angle
% http://en.wikipedia.org/wiki/Solar_elevation_angle
%
% Mark A. Mikofski
% Copyright (c) 2013
%
% datetime
if iscellstr(datetime) || ~isvector(datetime)
    datetime = datenum(datetime); % [days] dates & times
else
    datetime = datetime(:); % convert datenums to row
end
date = floor(datetime); % [days]
[year,~,~] = datevec(date);
time = datetime - date; % [days]
% constants
toRadians = @(x)x*pi/180; % convert degrees to radians
toDegrees = @(x)x*180/pi; % convert radians to degrees
% Equation of time
d_n = mod(date-datenum(year,1,1)+1,365); % day number
B = 2*pi*(d_n-81)/365; % ET parameter
ET = 9.87*sin(2*B)-7.53*cos(B)-1.5*sin(B); % [minutes] equation of time
% approximate solar time
solarTime = ((time*24-double(dst))*60+4*(longitude-time_zone*15)+ET)/60/24;
latitude_rad = toRadians(latitude); % [radians] latitude
rotation_rad = toRadians(rotation); % [radians] field rotation
t_h = (solarTime*24-12)*15; % [degrees] hour angle
t_h_rad = toRadians(t_h); % [radians]
delta = -23.45 * cos(2*pi*(d_n+10)/365); % [degrees] declination
delta_rad = toRadians(delta); % [radians]
theta_rad = acos(sin(latitude_rad)*sin(delta_rad)+ ...
    cos(latitude_rad)*cos(delta_rad).*cos(t_h_rad)); % [radians] zenith
theta = toDegrees(theta_rad); % [degrees] zenith
elevation = 90 - theta; % elevation
day = elevation>0; % day or night?
cos_phi = (cos(theta_rad)*sin(latitude_rad)- ...
    sin(delta_rad))./(sin(theta_rad)*cos(latitude_rad)); % cosine(azimuth)
% azimuth [0, 180], absolute value measured from due south, so east = west = 90,
% south = 0, north = 180
phi_south = acos(min(1,max(-1,cos_phi)));
% azimuth [0, 360], measured clockwise from due north, so east = 90,
% south = 180, and west = 270 degrees
phi_rad = NaN(size(phi_south)); % night azimuth is NaN
% shift from ATAN to ATAN2, IE: use domain from 0 to 360 degrees instead of
% from -180 to 180
phi_rad(day) = pi + sign(t_h(day)).*phi_south(day); % Shift domain to 0-360 deg
% projection of sun angle on x-z plane, measured from z-direction (up)
phi_x = toDegrees(atan2(sin(phi_rad-rotation_rad).*sin(theta_rad), ...
    cos(theta_rad))); % [degrees]
% projection of sun angle on y-z plane, measured from z-direction (up)
phi_y = toDegrees(atan2(cos(phi_rad-rotation_rad).*sin(theta_rad), ...
    cos(theta_rad))); % [degrees]
phi = toDegrees(phi_rad); % [degrees] azimuth
angles = [theta, phi]; % [degrees] zenith, azimuth
projection = [phi_x,phi_y]; % [degrees] x-z plane, y-z plane

end
% =========================================================================