% WBloop_2021_(WMO, WBANo, startyear) fixed re_DD_file and re_WBGT_file 13 Jan 2021
% REALIZED re_PC_file columns can add up to more thna days in the year when
% there have been NaN between classes, particularly patchy arctic icy or heat
% adjust icy:hea:hum:nat:ven:eva:DXX by factor=(total days/checkhours) returned in PC(2:8) 
% WBloop2020d % Basis [mean(stPa) Thumid BASEenthalpy]: Given mean station pressure (Pascals) and Thumid (wb �C), Basis of cooling enthalpy (kJ/kg-da) 
% ad to basis not enough dewpoint data
% humidification dehumidification
function [ASHRAEcz, timeframe, HDDstats, CDDstats, ChDstats, WBGTstatsW, WBGTstatsY, DXstats, PC, savings, Basis] = WBloop_2021_u(WMO,WBANo, startyear)
statyears=14; % minimum sample of years to draw rates
significant=0.05; % threshold to reject null hypothesis (no slope)
% initialize what is to be returned
    savings(1:5)=NaN;
    HDDstats(1:10)=NaN;
    CDDstats(1:10)=NaN;
    ASHRAEcz=NaN;
    MMM(1:12)=NaN;
    DXstats(1:10)=NaN;
    WBGTstatsW(1:10)=NaN;
    WBGTstatsY(1:10)=NaN;
    ChDstats(1:10)=NaN;
    PC(1:10)=NaN;
%%% past WB16Blocal.m below here - then remove WMO & startyear,
% also uncomment % close(gcf); % close current figure (two occurances)

% WB16B_local.m 24th March 2016
% to do yet:
% X. scale y-axis days per year [0 61 122 183 244 305 366]
% X. ensure dpC < dbC  on psych chart (wb_psych.m & H2O_psych.m)
% 5. check start end Nobs

% v16b add maximum monthly mean (MMM) over the most recent years
% v16b also average DX days per year

% www1/ftp, but also needs archive drive (USB D on windows)
%
% Be sure to include function psych.m 
% [wbC, H2O] = psych(dbC, dpC, stPa)

% NOW USE PEAK DATA ON FIGURE 2
% WBGT now on basis of maxC and MAXwbC
%%%% AND WRITE OUTPUT FILES AS GRAPHICS & TABLES
% checking db cooling degree days 
% substitution DX on enthalpy rather than wb DD

%%%% to do list
% ground temperatures
% design temperatures
% psych chart MIN-MEAN-MAX
%%
unix=1; % switch to 1 for unix, zero for windows
% unix=1; % for unix
if unix==1
    slash='/';
else
    slash='\';
end
%%
% psycrometrics from ASHRAE Fundamentals 1993
C(1) = -5.6745359E+03;
% C(2) = -5.1523058E-01;
C(2) = +6.3925247;
C(3) = -9.6778430E-03;
C(4) = +6.2215701E-07;
C(5) = +2.0747825E-09;
C(6) = -9.4840240E-13;
C(7) = +4.1635019;
C(8) = -5.8002206E+03;
% C(9) = +5.5162560;
C(9) = 1.3914993;
C(10)= -4.8640239E-02;
C(11)= +4.1764768E-05;
C(12)= -1.4452093E-08;
C(13)= +6.5459673;
C(14)= 6.54;
C(15)=14.526;
C(16)= 0.7389;
C(17)= 0.09486;
C(18)= 0.04569;
%%
ZoneBase=[18 10]; % ASHRAE Table A3 from 2013 Standard 169 base of heating and cooling degree days
dbcomfort=[18.3 24 27]; % �C dry bulb temperature bounds for comfort %%%%% adding third stage for ventilation
wbcomfort=[9.4 19.4 24]; % �C wet bulb temperature for wetbulb degree bounds for comfort %%%%% adding third stage for ventilation
TVhot=dbcomfort(3); TVhumid=wbcomfort(3); %%%%% 27�C dry bulb and 24�C wet bulb ventilation comfort limits (V)
Tcold=dbcomfort(1); Thot=dbcomfort(2); % Tcold=18.3; Thot=24; % dry bulb C limits of natural comfort
Thumid=wbcomfort(2); Tdry=wbcomfort(1); % Thumid=19.4 wet bulb C limits of comfort % Tdry = 9.4
Ticing=5; deltaTice=5; % Climate Zone Map "ice-build-up" not above 5C and dewpoint not more than 5C below drybulb

wbgtcolour= [15 20 26.6 27.8 29.4 31.1 32.2]; % �C base temperatures for "blue", "clear", "white", "green", "yellow", "red", "black"
% global wet bulb temperature (WBGT) classifications

% startyear= 1990; % rather than whole timeseries
stopyear=2020; % 2016; % stop time series


% 
if unix==1
    localpath='/home/home01/fbsepet/nb/restore.1605057541/ftp.ncdc.noaa.gov/pub/data/gsod/'; % arc3 nobackup
	% localpath='/RDS/Q0079/pub/NOAAdata/www1.ncdc.noaa.gov/pub/data/gsod/'; % running unix
    keypath='/home/home01/fbsepet/nb/WB/'; % running unix on arc 3
else
    localpath='D:\gsod\'; % running windows with removable drive
    keypath='D:\WB\'; % running windows
end 
listfile='isd-history_2021.csv';
filekey = strcat(keypath,listfile);
delimiter = ',';
startRow = 2;
%"USAF","WBAN","STATION NAME","CTRY","STATE","ICAO","LAT","LON","ELEV(M)","BEGIN","END"
formatSpec = '%f%f%s%s%s%s%f%f%f%f%f%[^\n\r]';
fileID = fopen(filekey,'r');
dataArray0 = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'EmptyValue' ,NaN,'HeaderLines' ,startRow-1, 'ReturnOnError', false);
fclose(fileID);

% startline=1000; stepper=3;
% for WMO=27260:10:27280
% WMO=27270;
    %j=startline:startline+stepper % initial test of looping structure down the list
    
goodyears=0; % count how many of the timeseries actually work
USAF = dataArray0{:, 1}; % WMO

% WMO = USAF(j); % work down the list based on index j

WBAN = dataArray0{:, 2};
STATIONNAME = dataArray0{:, 3};
CTRY = dataArray0{:, 4};
STATE = dataArray0{:, 5};
ICAO = dataArray0{:, 6};
LAT = dataArray0{:, 7};
LON = dataArray0{:, 8};
ELEVM = dataArray0{:, 9};
BEGIN = dataArray0{:, 10};
END = dataArray0{:, 11};
clearvars filekey delimiter startRow formatSpec fileID ans;
%% 
noaaRow=find(USAF==WMO & WBAN==WBANo);
    WBANstring=num2str(WBANo);
    switch length(WBANstring) % how long is this string of digits?
        case 1 % one digit ... to -five-
            WBANstring=strcat('0000',WBANstring); 
        case 2 % two digits ... to -five-
            WBANstring=strcat('000',WBANstring); 
        case 3 % three digits ... to -five-
            WBANstring=strcat('00',WBANstring); 
        case 4 % four digits ... to -five-
            WBANstring=strcat('0',WBANstring); 
        case 5 % five digits ... to -five-
            % as is 
        otherwise % zero digits ... to 99999
            WBANstring='-99999-'; 
    end

WMOstring=num2str(WMO);
switch length(WMOstring) % how long is this string of digits?
    case 1
        WMOstring=strcat('00000',WMOstring); 
    case 2
        WMOstring=strcat('0000',WMOstring); 
    case 3
        WMOstring=strcat('000',WMOstring); 
    case 4
        WMOstring=strcat('00',WMOstring); 
    case 5
        WMOstring=strcat('0',WMOstring); 
    case 6
        % as is  
    otherwise
        WMOstring=strcat('999999'); 
end
%%
outpath=strcat(keypath,'output',slash,WMOstring,'-',WBANstring,slash);
if ~exist(outpath,'dir'); mkdir(outpath); end

% output graphics filenames
fig1PNG=strcat(outpath,WMOstring,'-',WBANstring,'_fig1.png');
fig2PNG=strcat(outpath,WMOstring,'-',WBANstring,'_fig2.png');
% fig1EPS=strcat(outpath,WMOstring,'_fig1.eps');fig2EPS=strcat(outpath,WMOstring,'_fig2.eps');
% fig1WMF=strcat(outpath,WMOstring,'_fig1.emf');fig2WMF=strcat(outpath,WMOstring,'_fig2.emf');

% output data filenames
QAfile=strcat(outpath,WMOstring,'-',WBANstring,'_QA.txt');
re_DD_file=strcat(outpath,WMOstring,'-',WBANstring,'_DD.txt');
re_WBGT_file=strcat(outpath,WMOstring,'-',WBANstring,'_WBGT.txt');
re_PC_file=strcat(outpath,WMOstring,'-',WBANstring,'_PC.txt');

gdays=0; % log of all days processed

starter=max(startyear,min(ceil(BEGIN(noaaRow)/10000))) % after this point use startyear to set minumum goodyears
    startyear=starter; 
stopper=min(stopyear,max(floor(END(noaaRow)/10000))) % after this point use stopyear to set maximum goodyears
    stopyear=stopper;
yrcnt=0;
% annum=2000
QA(1:1+stopper-starter,1)=WMO; % identity of world meteorological organization site (filling first column)
QA(1:1+stopper-starter,2)=WBANo; % US station id
QA(1:1+stopper-starter,3:10)=NaN; % record of quality assurance
timeset=0; % switched to one after recording the first occurence
MM(starter:stopper,1:12,1:3)=NaN; % monthly mean v16b (20Mar2017)
% re_WBGT(1:(1+stopper-starter),1:10)=NaN;
% re_PC(1:(1+stopper-starter),1:10)=NaN;
% re_DD(1:(1+stopper-starter),1:13)=NaN;
for annum=starter:stopper % change for year of interest
yrcnt=yrcnt+1;
QA(yrcnt,3)=annum; % quality assurance record year numbered down third column
% https://en.wikipedia.org/wiki/Atmospheric_pressure
SSLP=101325; % standard sea level pressure (Pascals)
L=0.0065; % K/m temperature lapse rate, = g/cp for dry air
Cp=1007; % J/(kg�K)% cp	constant pressure specific heat	~ 
UTP=22632; % standard pressure at top of troposphere (Pascals)
T0=288.15; % K sea level standard temperature
g=9.80665; % m/s2 Earth-surface gravitational acceleration	
molar=0.0289644; % kg/mol molar mass of dry air
R0=8.31447; % J/(mol�K) Universal Gas Counstant
TROPO=11019; % standard thickness of troposphere (meters)
BGAMSL=-611; % Base Geometric Altitude above MSL % lookup station data
M=csvread(strcat(keypath,'stations.txt'));
[~,idx]=ismember(WMO,M(:,1));
idx(idx==0)=[];
data=M(idx,:);
if isempty(idx)
    LATI=mean(LAT(noaaRow));
    LONG=mean(LON(noaaRow));
    ELEV=mean(ELEVM(noaaRow));
    disp('not in ASHRAE station dataset ... use standard atmosphere model of STDP');
    % STDP=(SSLP+(UTP-SSLP)*(ELEV-BGAMSL)/TROPO)/1000 % kPa default station pressure
    STDP=(SSLP/1000)*(1-L*ELEV/T0)^((g*molar)/(R0*L)); % kPa default station pressure
else
    LATI=data(2); % degrees north
    LONG=data(3); % degrees east
    ELEV=data(4); % meters
    STDP=data(5); % kPa
end
QA(yrcnt,4)=STDP; % To check the nominal station pressure in case files do not have such
clear M data;
placelist = strcat(keypath,'placenames.txt');
delimiter = ''; formatSpec = '%s%[^\n\r]'; fileID = fopen(placelist,'r');
dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'EmptyValue' ,NaN, 'ReturnOnError', false);
fclose(fileID);
PLACENAME = dataArray{:, 1};
if isempty(idx)
    place=STATIONNAME(min(noaaRow));
else
    place=PLACENAME{idx};
end
version = strcat(place,' (WMO ',WMOstring,') year ',num2str(annum)) % for example 2010 Brisbane Aero WMO 945780
clearvars placelist delimiter formatSpec fileID dataArray ans PLACENAME idx;
%% filenames and paths
filename=strcat(WMOstring,'-',WBANstring,'-',num2str(annum),'.op'); % input datafile
zfilename=strcat(filename,'.gz'); % zipped input datafile

SORTEDnoaapath=strcat(keypath,'NOAAdata',slash,WMOstring,'-',WBANstring,slash); % sorted NOAA data
outpath=strcat(keypath,'output',slash,WMOstring,'-',WBANstring,slash); % output of this script
wetpath=strcat(keypath,'wetput',slash,WMOstring,'-',WBANstring,slash); % wetbulb transformation of NOAA data

minWBfile=strcat(wetpath,filename,'-MIN-WB-W.txt'); % output MIN wetbulb datafile
WBfile=strcat(wetpath,filename,'-WB-W.txt'); % output wetbulb datafile
maxWBfile=strcat(wetpath,filename,'-MAX-WB-W.txt'); % output MAX wetbulb datafile


if ~exist(SORTEDnoaapath,'dir'); mkdir(SORTEDnoaapath); end
if ~exist(outpath,'dir'); mkdir(outpath); end
if ~exist(wetpath,'dir'); mkdir(wetpath); end

skip=0; % first assume you will find it either locally or on the server
% no need to download if available on local disk\
test1=exist(strcat(localpath,num2str(annum),slash,zfilename),'file')==0; % test 1 if it's missing
test2=0; % assume it's been copied over, because it's not ever there in localpath

if or(test1 , test2) % then a file is NOT already copied onto localpath
	% in which case requires a download of the zipped file
	try
        noaaFTP='ftp.ncdc.noaa.gov';
        NCDC=ftp(noaaFTP); % open ftp connection
        cd(NCDC,'pub/data/gsod'); % change directory to GSOD global daily
        cd(NCDC,num2str(annum)); % change directory to the year of interest
        details=dir(NCDC,strcat(WMOstring,'*'));
        if isempty(details)
            skip =1; % this year is not covered on their server
        else
            if ~exist(strcat(localpath,num2str(annum)),'dir')
                mkdir(strcat(localpath,num2str(annum)))
            end
            % download the zipped WMO station data for that particular year
            % there is only one case
            details=dir(NCDC, strcat(WMOstring,'-',WBANstring,'*'));
            if isempty(details)
            else
                QA(yrcnt,2)=WBAN(noaaRow);
                mget(NCDC,zfilename);
                gunzip(zfilename,'.');
                movefile(zfilename,strcat(localpath,num2str(annum),slash,zfilename));
            end
        end
        close(NCDC);
    catch
        formatSpeca = string('There is either nothing found or no FTP connection to %s so try WWW');
        str = sprintf(formatSpeca,strcat(noaaFTP,'/',num2str(annum),'/',zfilename))
        try
            noaaWWW='https://www1.ncdc.noaa.gov/pub/data/gsod/';
            fullURLpath=strcat(noaaWWW,num2str(annum)); % drill down to the year
            if ~exist(strcat(localpath,num2str(annum)),'dir'); mkdir(strcat(localpath,num2str(annum))); end
            % download the zipped WMO station data for that particular year
            % there is only one case
            QA(yrcnt,2)=WBAN(noaaRow);
            urlwrite(strcat(fullURLpath,'/',zfilename),strcat(localpath,num2str(annum),slash,zfilename));
        catch
            formatSpecb = string('There is either no WWW connection or nothing found at %s');
            str = sprintf(formatSpecb,strcat(fullURLpath,'/',zfilename))
            skip=1; % there is no internet connection to noaa servers so let's try with what we have locally
        end
	end
end
% if skip has value of one then there are results this year
if skip==0 % go ahead and process this year 
    %% following from import wizard
startRow = 2;
% Read columns of data as strings:
% For more information, see the TEXTSCAN documentation.
formatSpec = '%6s%6s%6s%2s%2s%8s%3s%8s%3s%8s%3s%8s%3s%7s%3s%7s%3s%7s%7s%8s%8s%1s%7s%6s%s%[^\n\r]';
%%
testthere=0; % skeptically preset that one of the files is not found locally
if exist(strcat(localpath,num2str(annum),slash,zfilename))==2
    testthere=1; % switch if standard filename found locally
    try %%%%%%
        gunzip(strcat(localpath,num2str(annum),slash,zfilename),'.'); %%%%%% catch if file is corrupt
    catch
        formatSpecc = string('The file %s may be corrupted so try a fresh FTP download.');
        noaaFTP='ftp.ncdc.noaa.gov';
        str = sprintf(formatSpecc,strcat(noaaFTP,'/',num2str(annum),'/',zfilename))
    	try
            noaaFTP='ftp.ncdc.noaa.gov';
            NCDC=ftp(noaaFTP); % open ftp connection
            cd(NCDC,'pub/data/gsod'); % change directory to GSOD global daily
            cd(NCDC,num2str(annum)); % change directory to the year of interest
            details=dir(NCDC,strcat(WMOstring,'*'));
            if isempty(details)
                skip =1; % this year is not covered on their server
            else
                if ~exist(strcat(localpath,num2str(annum)),'dir'); mkdir(strcat(localpath,num2str(annum))); end
                % download the zipped WMO station data for that particular year
                % there is only one case
                details=dir(NCDC, strcat(WMOstring,'-',WBANstring,'*'));
                if isempty(details)
                else
                    QA(yrcnt,2)=WBAN(noaaRow);
                    mget(NCDC,zfilename);
                    gunzip(zfilename,'.');
                    movefile(zfilename,strcat(localpath,num2str(annum),slash,zfilename));
                end
            end
            close(NCDC);
        catch
            formatSpeca = string('There is either nothing found or no FTP connection to %s so try WWW');
            str = sprintf(formatSpeca,strcat(noaaFTP,'/',num2str(annum),'/',zfilename))
            try
                noaaWWW='https://www1.ncdc.noaa.gov/pub/data/gsod/';
                fullURLpath=strcat(noaaWWW,num2str(annum)); % drill down to the year
                if ~exist(strcat(localpath,num2str(annum)),'dir'); mkdir(strcat(localpath,num2str(annum))); end
                % download the zipped WMO station data for that particular year
                % there is only one case
                QA(yrcnt,2)=WBAN(noaaRow);
                urlwrite(strcat(fullURLpath,'/',zfilename),strcat(localpath,num2str(annum),slash,zfilename));
            catch
                formatSpecb = string('There is either no WWW connection or nothing found at %s');
                str = sprintf(formatSpecb,strcat(fullURLpath,'/',zfilename))
                skip=1; % there is no internet connection to noaa servers so let's try with what we have locally
            end
        end
    end %%%%%%
    fileID = fopen(filename,'r');
    dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'HeaderLines' ,startRow-1, 'ReturnOnError', false);
    fclose(fileID);
    delete(filename);
end
if testthere==1
    startyear=max(annum,startyear); % possibly increase the startyear = annum
    stopyear=min(annum,stopyear); % possibly decrease the stopyear = annum
end
% puts a copy of filename in local working directory
%% Read columns of data according to format string.
% This call is based on the structure of the file used to generate this
% code. If an error occurs for a different file, try regenerating the code
% from the Import Tool.
% dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '', 'HeaderLines' ,startRow-1, 'ReturnOnError', false);

%% Close the text file.
% fclose(fileID);
% delete(filename);
%% Convert the contents of columns containing numeric strings to numbers.
% Replace non-numeric strings with NaN.
raw = repmat({''},length(dataArray{1}),length(dataArray)-1);
for col=1:length(dataArray)-1
    raw(1:length(dataArray{col}),col) = dataArray{col};
end
numericData = NaN(size(dataArray{1},1),size(dataArray,2));

for col=[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,23,24,25]
    % Converts strings in the input cell array to numbers. Replaced non-numeric
    % strings with NaN.
    rawData = dataArray{col};
    for row=1:size(rawData, 1);
        % Create a regular expression to detect and remove non-numeric prefixes and
        % suffixes.
        regexstr = '(?<prefix>.*?)(?<numbers>([-]*(\d+[\,]*)+[\.]{0,1}\d*[eEdD]{0,1}[-+]*\d*[i]{0,1})|([-]*(\d+[\,]*)*[\.]{1,1}\d+[eEdD]{0,1}[-+]*\d*[i]{0,1}))(?<suffix>.*)';
        try
            result = regexp(rawData{row}, regexstr, 'names');
            numbers = result.numbers;
            
            % Detected commas in non-thousand locations.
            invalidThousandsSeparator = false;
            if any(numbers==',');
                thousandsRegExp = '^\d+?(\,\d{3})*\.{0,1}\d*$';
                if isempty(regexp(thousandsRegExp, ',', 'once'));
                    numbers = NaN;
                    invalidThousandsSeparator = true;
                end
            end
            % Convert numeric strings to numbers.
            if ~invalidThousandsSeparator;
                numbers = textscan(strrep(numbers, ',', ''), '%f');
                numericData(row, col) = numbers{1};
                raw{row, col} = numbers{1};
            end
        catch me
        end
    end
end

%% Split data into numeric and cell columns.
rawNumericColumns = raw(:, [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,23,24,25]);
rawCellColumns = raw(:, 22);


%% Allocate imported array to column variable names
YEAR = cell2mat(rawNumericColumns(:, 3));
MO = cell2mat(rawNumericColumns(:, 4));
DA = cell2mat(rawNumericColumns(:, 5));
TEMP = cell2mat(rawNumericColumns(:, 6));
DEWP = cell2mat(rawNumericColumns(:, 8));
STP = cell2mat(rawNumericColumns(:, 12));
MAX = cell2mat(rawNumericColumns(:, 20));
MIN = cell2mat(rawNumericColumns(:, 21));
PRCP = cell2mat(rawNumericColumns(:, 22));

%% fill in isolated one or two missing station temperature readings
QA(yrcnt,5)=length(YEAR); % number of records that year
Tfill=0;
if length(TEMP)>3
if and(and(TEMP(1)==9999.9,TEMP(2)==9999.9),and(not(TEMP(3)==9999.9),not(TEMP(4)==9999.9)))
    TEMP(2)= TEMP(3)-(TEMP(4)-TEMP(3));
    Tfill=Tfill+1;
end
if and(TEMP(1)==9999.9,and(not(TEMP(2)==9999.9),not(TEMP(3)==9999.9)))
    TEMP(1)= TEMP(2)-(TEMP(3)-TEMP(2));
    Tfill=Tfill+1;
end
for tex=2:length(TEMP)-2
    if and(and(TEMP(tex)==9999.9,TEMP(tex+1)==9999.9),and(not(TEMP(tex-1)==9999.9),not(TEMP(tex+2)==9999.9)))
        TEMP(tex)=(2*TEMP(tex-1)+TEMP(tex+2))/3;
        Tfill=Tfill+1;
    end
end
for tex=2:length(TEMP)-1
    if and(TEMP(tex)==9999.9,and(not(TEMP(tex-1)==9999.9),not(TEMP(tex+1)==9999.9)))
        TEMP(tex)=(TEMP(tex-1)+TEMP(tex+1))/2;
        Tfill=Tfill+1;
    end
end
if and(and(TEMP(end-1)==9999.9,TEMP(end)==9999.9),and(not(TEMP(end-3)==9999.9),not(TEMP(end-2)==9999.9)))
    TEMP(end-1)= TEMP(end-2)-(TEMP(end-3)-TEMP(end-2));
    Tfill=Tfill+1;
end
if and(TEMP(end)==9999.9,and(not(TEMP(end-2)==9999.9),not(TEMP(end-1)==9999.9)))
    TEMP(end)= TEMP(end-1)-(TEMP(end-2)-TEMP(end-1));
    Tfill=Tfill+1;
end
if Tfill>0
    QA(yrcnt,6)=Tfill; % output
else
    QA(yrcnt,6)=0;
end
if max(TEMP)==9999.9
    disp('not enough temperature data ...');
    QA(yrcnt,6)=NaN;
end
for dexy=1:length(TEMP)
    if TEMP(dexy)==9999.9
        TEMP(dexy)=NaN;
    end
end
else
    QA(yrcnt,6)=NaN;
end
%% fill in isolated one or two missing station dewpoint readings
Dfill=0;
if length(DEWP)>3
if and(and(DEWP(1)==9999.9,DEWP(2)==9999.9),and(not(DEWP(3)==9999.9),not(DEWP(4)==9999.9)))
    DEWP(2)= DEWP(3)-(DEWP(4)-DEWP(3));
    Dfill=Dfill+1;
end
if and(DEWP(1)==9999.9,and(not(DEWP(2)==9999.9),not(DEWP(3)==9999.9)))
    DEWP(1)= DEWP(2)-(DEWP(3)-DEWP(2));
    Dfill=Dfill+1;
end
for tex=2:length(DEWP)-2
    if and(and(DEWP(tex)==9999.9,DEWP(tex+1)==9999.9),and(not(DEWP(tex-1)==9999.9),not(DEWP(tex+2)==9999.9)))
        DEWP(tex)=(2*DEWP(tex-1)+DEWP(tex+2))/3;
        Dfill=Dfill+1;
    end
end
for tex=2:length(DEWP)-1
    if and(DEWP(tex)==9999.9,and(not(DEWP(tex-1)==9999.9),not(DEWP(tex+1)==9999.9)))
        DEWP(tex)=(DEWP(tex-1)+DEWP(tex+1))/2;
        Dfill=Dfill+1;
    end
end
if and(and(DEWP(end-1)==9999.9,DEWP(end)==9999.9),and(not(DEWP(end-3)==9999.9),not(DEWP(end-2)==9999.9)))
    DEWP(end-1)= DEWP(end-2)-(DEWP(end-3)-DEWP(end-2));
    Dfill=Dfill+1;
end
if and(DEWP(end)==9999.9,and(not(DEWP(end-2)==9999.9),not(DEWP(end-1)==9999.9)))
    DEWP(end)= DEWP(end-1)-(DEWP(end-2)-DEWP(end-1));
    Dfill=Dfill+1;
end
if Dfill>0
    QA(yrcnt,7)=Dfill; % output
else
    QA(yrcnt,7)=0;
end
if max(DEWP)==9999.9
    disp('not enough dewpoint data ...');
    QA(yrcnt,7)=NaN;
end
for dexy=1:length(DEWP)
    if DEWP(dexy)==9999.9
        DEWP(dexy)=NaN;
    end
end
else
    QA(yrcnt,7)=NaN;
end
%% fill in isolated one or two missing station pressure readings

Pfill=0;
%{
if length(STP)>3
if and(and(STP(1)==9999.9,STP(2)==9999.9),and(not(STP(3)==9999.9),not(STP(4)==9999.9)))
    STP(2)= STP(3)-(STP(4)-STP(3));
    Pfill=Pfill+1;
end
if and(STP(1)==9999.9,and(not(STP(2)==9999.9),not(STP(3)==9999.9)))
    STP(1)= STP(2)-(STP(3)-STP(2));
    Pfill=Pfill+1;
end
for tex=2:length(STP)-2
    if and(and(STP(tex)==9999.9,STP(tex+1)==9999.9),and(not(STP(tex-1)==9999.9),not(STP(tex+2)==9999.9)))
        STP(tex)=(2*STP(tex-1)+STP(tex+2))/3;
        Pfill=Pfill+1;
    end
end
for tex=2:length(STP)-1
    if and(STP(tex)==9999.9,and(not(STP(tex-1)==9999.9),not(STP(tex+1)==9999.9)))
        STP(tex)=(STP(tex-1)+STP(tex+1))/2;
        Pfill=Pfill+1;
    end
end
if and(and(STP(end-1)==9999.9,STP(end)==9999.9),and(not(STP(end-3)==9999.9),not(STP(end-2)==9999.9)))
    STP(end-1)= STP(end-2)-(STP(end-3)-STP(end-2));
    Pfill=Pfill+1;
end
if and(STP(end)==9999.9,and(not(STP(end-2)==9999.9),not(STP(end-1)==9999.9)))
    STP(end)= STP(end-1)-(STP(end-2)-STP(end-1));
    Pfill=Pfill+1;
end
%}
for tex=1:length(STP)
    if STP(tex)==9999.9
        STP(tex)=STDP*10; % 1000 hectopascals = 100 kPa
        Pfill=Pfill+1;
    end
end
if and(Pfill>0, Pfill<=366)
    QA(yrcnt,8)=Pfill; % output
else
    QA(yrcnt,8)=0;
end

%% fill in isolated one or two missing station MAX daily temperature readings
maxTfill=0;
if length(MAX)>3
if and(and(MAX(1)==9999.9,MAX(2)==9999.9),and(not(MAX(3)==9999.9),not(MAX(4)==9999.9)))
    MAX(2)= MAX(3)-(MAX(4)-MAX(3));
    maxTfill=maxTfill+1;
end
if and(MAX(1)==9999.9,and(not(MAX(2)==9999.9),not(MAX(3)==9999.9)))
    MAX(1)= MAX(2)-(MAX(3)-MAX(2));
    maxTfill=maxTfill+1;
end
for tex=2:length(MAX)-2
    if and(and(MAX(tex)==9999.9,MAX(tex+1)==9999.9),and(not(MAX(tex-1)==9999.9),not(MAX(tex+2)==9999.9)))
        MAX(tex)=(2*MAX(tex-1)+MAX(tex+2))/3;
        maxTfill=maxTfill+1;
    end
end
for tex=2:length(MAX)-1
    if and(MAX(tex)==9999.9,and(not(MAX(tex-1)==9999.9),not(MAX(tex+1)==9999.9)))
        MAX(tex)=(MAX(tex-1)+MAX(tex+1))/2;
        maxTfill=maxTfill+1;
    end
end
if and(and(MAX(end-1)==9999.9,MAX(end)==9999.9),and(not(MAX(end-3)==9999.9),not(MAX(end-2)==9999.9)))
    MAX(end-1)= MAX(end-2)-(MAX(end-3)-MAX(end-2));
    maxTfill=maxTfill+1;
end
if and(MAX(end)==9999.9,and(not(MAX(end-2)==9999.9),not(MAX(end-1)==9999.9)))
    MAX(end)= MAX(end-1)-(MAX(end-2)-MAX(end-1));
    maxTfill=maxTfill+1;
end
if maxTfill>0
    QA(yrcnt,9)=maxTfill; % output
else
    QA(yrcnt,9)=0;
end
if max(MAX)==9999.9
    disp('not enough MAX daily temperature data ...');
    QA(yrcnt,9)=NaN;
end
for dexy=1:length(MAX)
    if MAX(dexy)==9999.9
        MAX(dexy)=NaN;
    end
end
else
    QA(yrcnt,9)=NaN;
end
%% fill in isolated one or two missing station MIN daily temperature readings
minTfill=0;
if length(MIN)>3
if and(and(MIN(1)==9999.9,MIN(2)==9999.9),and(not(MIN(3)==9999.9),not(MIN(4)==9999.9)))
    MIN(2)= MIN(3)-(MIN(4)-MIN(3));
    minTfill=minTfill+1;
end
if and(MIN(1)==9999.9,and(not(MIN(2)==9999.9),not(MIN(3)==9999.9)))
    MIN(1)= MIN(2)-(MIN(3)-MIN(2));
    minTfill=minTfill+1;
end
for tex=2:length(MIN)-2
    if and(and(MIN(tex)==9999.9,MIN(tex+1)==9999.9),and(not(MIN(tex-1)==9999.9),not(MIN(tex+2)==9999.9)))
        MIN(tex)=(2*MIN(tex-1)+MIN(tex+2))/3;
        minTfill=minTfill+1;
    end
end
for tex=2:length(MIN)-1
    if and(MIN(tex)==9999.9,and(not(MIN(tex-1)==9999.9),not(MIN(tex+1)==9999.9)))
        MIN(tex)=(MIN(tex-1)+MIN(tex+1))/2;
        minTfill=minTfill+1;
    end
end
if and(and(MIN(end-1)==9999.9,MIN(end)==9999.9),and(not(MIN(end-3)==9999.9),not(MIN(end-2)==9999.9)))
    MIN(end-1)= MIN(end-2)-(MIN(end-3)-MIN(end-2));
    minTfill=minTfill+1;
end
if and(MIN(end)==9999.9,and(not(MIN(end-2)==9999.9),not(MIN(end-1)==9999.9)))
    MIN(end)= MIN(end-1)-(MIN(end-2)-MIN(end-1));
    minTfill=minTfill+1;
end
if minTfill>0
    QA(yrcnt,10)=minTfill; % output
else
    QA(yrcnt,10)=0;
end
if max(MIN)==9999.9
    disp('not enough MIN daily temperature data ...');
    QA(yrcnt,10)=NaN;
end
for dexy=1:length(MIN)
    if MIN(dexy)==9999.9
        MIN(dexy)=NaN;
    end
end
else
    QA(yrcnt,10)=NaN;
end

%% fill in isolated one or two missing station daily PRCP readings
PRCPfill=0;
if length(PRCP)>3
if and(and(PRCP(1)==9999.9,PRCP(2)==9999.9),and(not(PRCP(3)==9999.9),not(PRCP(4)==9999.9)))
    PRCP(2)= PRCP(3)-(PRCP(4)-PRCP(3));
    PRCPfill=PRCPfill+1;
else
    PRCP(2)=NaN;
end
if and(PRCP(1)==9999.9,and(not(PRCP(2)==9999.9),not(PRCP(3)==9999.9)))
    PRCP(1)= PRCP(2)-(PRCP(3)-PRCP(2));
    PRCPfill=PRCPfill+1;
end
for tex=2:length(PRCP)-2
    if and(and(PRCP(tex)==9999.9,PRCP(tex+1)==9999.9),and(not(PRCP(tex-1)==9999.9),not(PRCP(tex+2)==9999.9)))
        PRCP(tex)=(2*PRCP(tex-1)+PRCP(tex+2))/3;
        PRCPfill=PRCPfill+1;
    end
end
for tex=2:length(PRCP)-1
    if and(PRCP(tex)==9999.9,and(not(PRCP(tex-1)==9999.9),not(PRCP(tex+1)==9999.9)))
        PRCP(tex)=(PRCP(tex-1)+PRCP(tex+1))/2;
        PRCPfill=PRCPfill+1;
    end
end
if and(and(PRCP(end-1)==9999.9,PRCP(end)==9999.9),and(not(PRCP(end-3)==9999.9),not(PRCP(end-2)==9999.9)))
    PRCP(end-1)= PRCP(end-2)-(PRCP(end-3)-PRCP(end-2));
    PRCPfill=PRCPfill+1;
end
if and(PRCP(end)==9999.9,and(not(PRCP(end-2)==9999.9),not(PRCP(end-1)==9999.9)))
    PRCP(end)= PRCP(end-1)-(PRCP(end-2)-PRCP(end-1));
    PRCPfill=PRCPfill+1;
end
if PRCPfill>0
    QA(yrcnt,11)=PRCPfill; % output
else
    QA(yrcnt,11)=0;
end
if max(PRCP)==9999.9
    disp('not enough PRCP daily temperature data ...');
    QA(yrcnt,11)=NaN;
end
for dexy=1:length(PRCP)
    if PRCP(dexy)==9999.9
        PRCP(dexy)=NaN;
    end
end
else
    QA(yrcnt,11)=NaN;
end
%% Clear temporary variables
clearvars zfilename startRow formatSpec fileID dataArray ans raw col numericData rawData row regexstr result numbers invalidThousandsSeparator thousandsRegExp me rawNumericColumns rawCellColumns;

dateline=datenum(YEAR,MO,DA);
dbC=round((TEMP-32)*50/9)/10; 
dpC=round((DEWP-32)*50/9)/10;
stPa=round(STP*100); % sea level typically 101300 Pascals or 1013 hectopascals
maxC=round((MAX-32)*50/9)/10;
minC=round((MIN-32)*50/9)/10;
mm=single(int8(round(PRCP*254)))/10;
clear TEMP DEWP MAX MIN PRCP STP
%%
if length(dbC)>=300 % 346 = 95% of 365 reasonable threshold to consider the year
    goodyears=goodyears+1; % count one more goodyear found on server
    yrcnter(goodyears)=yrcnt; % growing vector of the completed records

%% check if WB is on local disk
if (exist(strcat(wetpath,minWBfile),'file')==0 || exist(strcat(wetpath,WBfile),'file')==0 || exist(strcat(wetpath,maxWBfile),'file')==0) % then at least one of the wetbulb datafiles is not already on local disk

    %% applied to MIN daily temperature (minC)
    MINwbC(1:length(minC))=NaN; minH2O(1:length(minC))=NaN;
    for i=1:length(minC)
        MINwbC(i) = psych(minC(i), dpC(i), stPa(i));
        minH2O(i) = H2O_psych(minC(i), dpC(i), stPa(i));
        WB_H2O(i,1)=YEAR(i); % year
        WB_H2O(i,2)=MO(i); % month
        WB_H2O(i,3)=DA(i); % day           
        WB_H2O(i,4)=double(int32(minC(i)*10))/10; % min drybulb
        WB_H2O(i,5)=double(int32(MINwbC(i)*10))/10; % MIN wetbulb
        WB_H2O(i,6)=double(int32(dpC(i)*10))/10; % dewpoint
        WB_H2O(i,7)=double(int32(minH2O(i)*10))/10; % absolute humidity
        WB_H2O(i,8)=round(stPa(i)); % pressure
    end
    MINwbC=MINwbC'; % vertical array of MAX wetbulb temperature
    minH2O=minH2O'; % verical array of Humidity expressed in g H20 per kg dry air 
    Wmin=minH2O/1000; % non-dimensional humidity
    % output coincidence of (db wb dp H2O stPa) in matrix WB_H2O
    fid = fopen(minWBfile,'w'); % opens file for writing 'w'
    fprintf(fid,'year mo da minDB minWB dewpt H2O_W StnPa\n');
    for dealer=1:size(WB_H2O,1)     
        fprintf(fid,'%4g %2g %2g %5.1f %5.1f %5.1f %5.1f %5.0f\n',WB_H2O(dealer,:));
    end
    fclose(fid); clear WB_H2O dealer
    %% applied to MEAN daily temperature (dbC)    
    wbC(1:length(dbC))=NaN; H2O(1:length(dbC))=NaN;
    for i=1:length(dbC)
        wbC(i) = psych(dbC(i), dpC(i), stPa(i));
        H2O(i) = H2O_psych(dbC(i), dpC(i), stPa(i));
        WB_H2O(i,1)=YEAR(i); % year
        WB_H2O(i,2)=MO(i); % month
        WB_H2O(i,3)=DA(i); % day           
        WB_H2O(i,4)=double(int32(dbC(i)*10))/10; % drybulb
        WB_H2O(i,5)=double(int32(wbC(i)*10))/10; % wetbulb
        WB_H2O(i,6)=double(int32(dpC(i)*10))/10; % dewpoint
        WB_H2O(i,7)=double(int32(H2O(i)*10))/10; % absolute humidity
        WB_H2O(i,8)=round(stPa(i)); % pressure
    end
    wbC=wbC'; % vertical array of MEAN wetbulb temperature
    H2O=H2O'; % verical array of Humidity expressed in g H20 per kg dry air 
    W=H2O/1000; % non-dimensional humidity
    % output coincidence of (db wb dp H2O stPa) in matrix WB_H2O
    fid = fopen(WBfile,'w'); % opens file for writing 'w'
    fprintf(fid,'year mo da aveDB aveWB dewpt H2O_W StnPa\n');
    for dealer=1:size(WB_H2O,1)     
        fprintf(fid,'%4g %2g %2g %5.1f %5.1f %5.1f %5.1f %5.0f\n',WB_H2O(dealer,:));
    end
    fclose(fid); clear WB_H2O dealer
    %% applied to MAX daily temperature (maxC)
    MAXwbC(1:length(maxC))=NaN; maxH2O(1:length(maxC))=NaN;
    for i=1:length(maxC)
        MAXwbC(i) = psych(maxC(i), dpC(i), stPa(i));
        maxH2O(i) = H2O_psych(maxC(i), dpC(i), stPa(i));
        WB_H2O(i,1)=YEAR(i); % year
        WB_H2O(i,2)=MO(i); % month
        WB_H2O(i,3)=DA(i); % day           
        WB_H2O(i,4)=double(int32(maxC(i)*10))/10; % max drybulb
        WB_H2O(i,5)=double(int32(MAXwbC(i)*10))/10; % MAX wetbulb
        WB_H2O(i,6)=double(int32(dpC(i)*10))/10; % dewpoint
        WB_H2O(i,7)=double(int32(maxH2O(i)*10))/10; % absolute humidity
        WB_H2O(i,8)=round(stPa(i)); % pressure
    end
    MAXwbC=MAXwbC'; % vertical array of MAX wetbulb temperature
    maxH2O=maxH2O'; % verical array of Humidity expressed in g H20 per kg dry air 
    Wmax=maxH2O/1000; % non-dimensional humidity
    % output coincidence of (db wb dp H2O stPa) in matrix WB_H2O
    fid = fopen(maxWBfile,'w'); % opens file for writing 'w'
    fprintf(fid,'year mo da maxDB maxWB dewpt H2O_W StnPa\n');
    for dealer=1:size(WB_H2O,1)     
        fprintf(fid,'%4g %2g %2g %5.1f %5.1f %5.1f %5.1f %5.0f\n',WB_H2O(dealer,:));
    end
    fclose(fid); clear WB_H2O dealer
else % WB_H2O datafile exists in wetpath so no need to do the previous analysis again
    
    % Open the minWB file
    formatSpec = '%4f%3f%3f%6f%6f%6f%6f%f%[^\n\r]';
    fileID = fopen(minWBfile,'r');
    dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '',  'ReturnOnError', false);
    fclose(fileID);
    % yyyy = dataArray{:, 1};
    % mm = dataArray{:, 2};
    % diem = dataArray{:, 3};
    % db_C = dataArray{:, 4};
    MINwbC = dataArray{:, 5};
    % dpC = dataArray{:, 6};
    Wmin = dataArray{:, 7};
    % Pa_ = dataArray{:, 8};
    clearvars fileID dataArray ans;
    
    
    fileID = fopen(WBfile,'r');
    dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '',  'ReturnOnError', false);
    fclose(fileID);
    % yyyy = dataArray{:, 1};
    % mm = dataArray{:, 2};
    % diem = dataArray{:, 3};
    % db_C = dataArray{:, 4};
    wbC = dataArray{:, 5};
    dpC = dataArray{:, 6};
    W = dataArray{:, 7};
    % Pa_ = dataArray{:, 8};
    clearvars fileID dataArray ans;
    
    % Open the maxWB file
    fileID = fopen(maxWBfile,'r');
    dataArray = textscan(fileID, formatSpec, 'Delimiter', '', 'WhiteSpace', '',  'ReturnOnError', false);
    fclose(fileID);
    % yyyy = dataArray{:, 1};
    % mm = dataArray{:, 2};
    % diem = dataArray{:, 3};
    % db_C = dataArray{:, 4};
    MAXwbC = dataArray{:, 5};
    % dpC = dataArray{:, 6};
    Wmax = dataArray{:, 7};
    % Pa_ = dataArray{:, 8};
    clearvars formatSpec fileID dataArray ans;
    
end % now you have wetbulb temperatures
%%
minDB=min(1,min(floor(dbC))); 
if minDB<1 
    minDB=1; % don't consider DD base temperatures under freezing conditions
end
maxDB=max(ceil(dbC));
minWB=min(1,min(floor(wbC))); 
if minWB<1 
    minWB=1; % don't consider DD base temperatures under freezing conditions
end
maxWB=max(ceil(wbC));

% wet bulb degree days base temperature Tbase (C) varied minWB:maxWB
if isnan(maxDB)
    DBCDD(minDB)=NaN;
else
    DBCDD(minDB:maxDB)=NaN; % not zeroing
end
if isnan(maxDB)
    DBHDD(minDB)=NaN;
else
    DBHDD(minDB:maxDB)=NaN; % not zeroing
end
if isnan(maxWB)
    WBCDD(minWB)=NaN;
else
    WBCDD(minWB:maxWB)=NaN; % not zeroing
end
if isnan(maxWB)
    WBHDD(minWB)=NaN
else
    WBHDD(minWB:maxWB)=NaN; % not zeroing
end
for doy=1:i % day of the year
    for TbaseWB=minWB:maxWB;
        if wbC(doy) > TbaseWB 
            WBCDD(TbaseWB) = WBCDD(TbaseWB) + (int32(wbC(doy)*100)/100 - TbaseWB); 
        end
        if wbC(doy) < TbaseWB 
            WBHDD(TbaseWB) = WBHDD(TbaseWB) + (TbaseWB - int32(wbC(doy)*100)/100); 
        end
    end
    for TbaseDB=minDB:maxDB;
        if dbC(doy) > TbaseDB
            DBCDD(TbaseDB) = DBCDD(TbaseDB) + (int32(dbC(doy)*100)/100 - TbaseDB);
        end
        if dbC(doy) < TbaseDB
            DBHDD(TbaseDB) = DBHDD(TbaseDB) + (TbaseDB - int32(dbC(doy)*100)/100);
        end
    end 
% invert new j-vectors from row to be column-wise
end
WBCDD=WBCDD';  WBHDD=WBHDD';  DBCDD=DBCDD'; DBHDD=DBHDD';


%% capture data for psychrometric chart plot from all years

% stPa = > AstPa
% dbC = > AdbC
% wbC = > AwbC
% H2O = > AH2O
% dpC = > AdpC

AstPa(gdays+1:gdays+length(dbC))=stPa;
AdbC(gdays+1:gdays+length(dbC))=dbC;
AwbC(gdays+1:gdays+length(dbC))=wbC;
AH2O(gdays+1:gdays+length(dbC))=H2O;
AdpC(gdays+1:gdays+length(dbC))=dpC;


AmaxC(gdays+1:gdays+length(maxC))=maxC;
AMAXwbC(gdays+1:gdays+length(MAXwbC))=MAXwbC;
AminC(gdays+1:gdays+length(minC))=minC;
AMINwbC(gdays+1:gdays+length(MINwbC))=MINwbC;

gdays=gdays+length(dbC); %

%% days (mean conditions) naturally comfortable, requiring heating, evap cooling, and active cooling

ASHP_days=0; heating_degree_days=0; icing_days=0; icing_degree_days=0;
evapworks_days=0; DXcooling_days=0; DXXcooling_days=0; DXcoolDehum_days=0; wb_cooling_degree_days(1:3)=0; db_cooling_degree_days=0; 
    CwbDays=0; cooling_enthalpy_days(1:3)=0; ASHRAEHDD=0; ASHRAECDD=0;
wbgtdd=0; wsCDD(1:3)=0; wsChD(1:3)=0;% wet bulb cooling degree days and enthalpy-days that can be avoided from DX with evap condensers
DXcoolDehumCdbDD=0; avoided_dbcdd=0; % standard cooling degree days, and those avoided with evap cooling
humidif_days=0; CdbDays=0; CwbDays=0; natural_days=0; wbgt_days(1:9)=0; ventn_days=0;
substitution(1:length(dbC),1:3)=NaN; % savings by wetting DX vapour compression cooling condenser
shed(1:length(dbC),1:3)=NaN; % shed load by wetting DX vapour compression cooling condenser
%% base enthalpy for substitution of evaporative condensers
%       ASHRAE Fundamentals 1993 SI edition Chapter 6, equation 30
%       enthalpy (kJ/kg) 1.006*Tdb(C)+W(kg/kg)*(2501+1.805*Tdb(C))
% adjusting units H2O/1000 to obtain dimensionless W
% for the purposes for establishing a base, use comfort wetbulb threshold Thumid
% Thumid=19.4 wbC limits of comfort (adjustable at start of this script)
% and furthermore assume saturation (set dbC the same)
% and so Ws (H20/1000) at saturation must be locked down at this point
%       HOF SI chapter 6 equation 21
%       Ws = 0.62198 * pw / (p - pw)

WBA=Thumid +273.15; % absolute temperature Kelvins of wet bulb condition
if Thumid<0 % below freezing
	lnpwsx = C(1)/WBA + C(2) + C(3)*WBA + C(4)*WBA^2 + C(5)*WBA^3 + C(6) * WBA^4 + C(7)*log(WBA);
	pwsx=exp(lnpwsx); % Pascals vapour pressure at saturation with WB
else % above freezing
	lnpwsx = C(8)/WBA + C(9) + C(10)*WBA + C(11)*WBA^2 + C(12)*WBA^3 + C(13)*log(WBA);
	pwsx=exp(lnpwsx); % Pascals vapour pressure at saturation with WB
end % now we have vapour pressure at WB saturation
% estimate pressure from average of the vector stPa
Wsatx=0.62198*pwsx/(mean(stPa)-pwsx); % g-H2O/g-air at saturation of WB
%
BASEenthalpy=1.006*Thumid+Wsatx*(2501+1.805*Thumid);
Basis(1:3)=[mean(stPa) Thumid BASEenthalpy]; % Given mean station pressure (Pascals) and Thumid (wb �C), Basis of cooling enthalpy (kJ/kg-da)
Basis(4)=NaN; % to be completed later after dewpoint data availability is known
%%

% WBGT now on basis of maxC and MAXwbC instead of dbC and wbC
for ldx=1:length(dbC)
    
    % degree days above wbgtcolour(3) "yellow", "red" AND "black"
    if 0.3*maxC(ldx)+0.7*MAXwbC(ldx)>=wbgtcolour(3)
        wbgtdd=wbgtdd+0.3*maxC(ldx)+0.7*MAXwbC(ldx)-wbgtcolour(3);
    end
    % classification of days by seven colours of wet bulb global temperature
    
        % v16b count days that year as wbgt_days(9)
   if ~isnan(MAXwbC(ldx))
       if isnan(wbgt_days(9))
           wbgt_days(9)=1;
       else
            wbgt_days(9)=wbgt_days(9)+1;
       end
   end
    
    if 0.3*maxC(ldx)+0.7*MAXwbC(ldx)<wbgtcolour(1)
        % below all wbgt scheme colours "deep purple"
        if isnan(wbgt_days(1)); wbgt_days(1)=1; else
        wbgt_days(1)=wbgt_days(1)+1; end
    end
    if 0.3*maxC(ldx)+0.7*MAXwbC(ldx)>=wbgtcolour(1) && 0.3*maxC(ldx)+0.7*MAXwbC(ldx)<wbgtcolour(2)
        % "deep blue"
        if isnan(wbgt_days(2)); wbgt_days(2)=1; else
        wbgt_days(2)=wbgt_days(2)+1; end
    else
        if 0.3*maxC(ldx)+0.7*MAXwbC(ldx)>=wbgtcolour(2) && 0.3*maxC(ldx)+0.7*MAXwbC(ldx)<wbgtcolour(3)
            % "ligth blue"
            if isnan(wbgt_days(3)); wbgt_days(3)=1; else
            wbgt_days(3)=wbgt_days(3)+1; end
        else
            if 0.3*maxC(ldx)+0.7*MAXwbC(ldx)>=wbgtcolour(3) && 0.3*maxC(ldx)+0.7*MAXwbC(ldx)<wbgtcolour(4)
                % white
                if isnan(wbgt_days(4)); wbgt_days(4)=1; else
                wbgt_days(4)=wbgt_days(4)+1; end
            else
                if 0.3*maxC(ldx)+0.7*MAXwbC(ldx)>=wbgtcolour(4) && 0.3*maxC(ldx)+0.7*MAXwbC(ldx)<wbgtcolour(5)
                    % green
                    if isnan(wbgt_days(5)); wbgt_days(5)=1; else
                    wbgt_days(5)=wbgt_days(5)+1; end
                else
                    if 0.3*maxC(ldx)+0.7*MAXwbC(ldx)>=wbgtcolour(5) && 0.3*maxC(ldx)+0.7*MAXwbC(ldx)<wbgtcolour(6)
                        % yellow
                        if isnan(wbgt_days(6)); wbgt_days(6)=1; else
                        wbgt_days(6)=wbgt_days(6)+1; end
                    else
                        if 0.3*maxC(ldx)+0.7*MAXwbC(ldx)>=wbgtcolour(6) && 0.3*maxC(ldx)+0.7*MAXwbC(ldx)<wbgtcolour(7)
                            % red
                            if isnan(wbgt_days(7)); wbgt_days(7)=1; else
                            wbgt_days(7)=wbgt_days(7)+1; end
                        else
                            if 0.3*maxC(ldx)+0.7*MAXwbC(ldx)>=wbgtcolour(7)
                                % black
                                if isnan(wbgt_days(8)); wbgt_days(8)=1; else
                                wbgt_days(8)=wbgt_days(8)+1; end
                            end
                        end
                    end
                end
            end
        end
    end
    if ~isnan(wbgt_days(9))
    wbgt_pc=100.*wbgt_days(1:8)/wbgt_days(9); end
    
    
    % v16b (20Mar2017) find monthly mean MM(YEAR,MO)
    if isnan(MM(annum,MO(ldx),2))
        MM(annum,MO(ldx),2)=dbC(ldx); % sum
        MM(annum,MO(ldx),3)=1; % count
        MM(annum,MO(ldx),1)=dbC(ldx); % mean
    else
        MM(annum,MO(ldx),2)=MM(annum,MO(ldx),2) + dbC(ldx); % sum
        MM(annum,MO(ldx),3)=MM(annum,MO(ldx),3) + 1; % count
        MM(annum,MO(ldx),1)=MM(annum,MO(ldx),2)/MM(annum,MO(ldx),3); % mean
    end
    % Should have set NaN outside loop
    
    if dbC(ldx)>ZoneBase(2) % greater than 10C (cooling desired in summer)
        ASHRAECDD=ASHRAECDD+dbC(ldx)-ZoneBase(2);
    end
    if dbC(ldx)<ZoneBase(1) % less than 18C (heating desired in winter)
        ASHRAEHDD=ASHRAEHDD+ZoneBase(1)-dbC(ldx);
    end
    
    % re_PC(10)
    % all potential cooling equal or above drybulb base temperature Thot 24C
    if dbC(ldx)>=Thot
        if isnan(CdbDays)
            CdbDays=1;
        else
            CdbDays=CdbDays+1;
        end
        db_cooling_degree_days=db_cooling_degree_days+dbC(ldx)-Thot;
    end
    
   % cooling AND dehumidification that calls for vapour compression work
   if wbC(ldx)>=Thumid % wetbulb >19.4C % re_PC(11)
       if isnan(CwbDays)
           CwbDays=1;
       else
           CwbDays=CwbDays+1;
       end 
       if isnan(DXcoolDehum_days)
           DXcoolDehum_days=1;
       else
           DXcoolDehum_days=DXcoolDehum_days+1;
       end
       DXcoolDehumCdbDD=DXcoolDehumCdbDD+dbC(ldx)-Thot;
       % carnot efficiency gain if evaporative condenser is
       % COPc = Qc/(Qh-Qc)
       % COPc = Tl/(Th-T1)
       % COPd = ADP_K/(db_K - ADP_K)
       % COPw = ADP_K/(wb_K - ADP_K)
       % ratio (d/w) = (wb-ADP)/(db-ADP)
       % reducing work into motor with latent heat of water
        ADP=8; % nominal A/C evaporator temperature
       % ADP=18; % nominal chilled beam evaporator temperature
       % ?T=8 K for split units and ?T=16 K for rooftop
       % overall ?T=10 K for the analysis of air-cooled dry heat exchangers
        deltaTdry=10; deltaTwet=7.5;
       % For cooling towers, range= 5.5 K and approach=2.5 K
       % based on discussions with industry representatives
       % For water based and ground heat exchangers systems approach = 2 K
       % and range=4K.  Additionally, all systems except dry air-cooled systems
       % additional temperature increase of 2K for additional heat exchanger
       % between the refrigeration liquid and water-based heat rejection circuit.
       % reference:-
       % Bruelisauera, Meggers, and Leibundgut (2013) 
       % Frontiers of Architectural Research 2:292�300
       % Choosing heat sinks for cooling in tropical climates
        substitution(ldx,1)=(1-(MINwbC(ldx)+deltaTwet-ADP)/(minC(ldx)+deltaTdry-ADP));
            wsCDD(1)=wsCDD(1)+(MINwbC(ldx)-Thumid)*substitution(ldx,1);
            wb_cooling_degree_days(1)=wb_cooling_degree_days(1)+MINwbC(ldx)-Thumid;
            enthalpy(ldx,1)=1.006*minC(ldx)+Wmin(ldx)*(2501+1.805*minC(ldx));
            if enthalpy(ldx,1)>BASEenthalpy
                cooling_enthalpy_days(1)=cooling_enthalpy_days(1)+enthalpy(ldx,1)-BASEenthalpy;
                wsChD(1)=wsChD(1)+(enthalpy(ldx,1)-BASEenthalpy)*substitution(ldx,1);
            end
            shed(ldx,1)=100*(enthalpy(ldx,1)-BASEenthalpy)*substitution(ldx,1)/(enthalpy(ldx,1)-BASEenthalpy);
            
        substitution(ldx,2)=(1-(wbC(ldx)+deltaTwet-ADP)/(dbC(ldx)+deltaTdry-ADP));
            wsCDD(2)=wsCDD(2)+(wbC(ldx)-Thumid)*substitution(ldx,2);
            wb_cooling_degree_days(2)=wb_cooling_degree_days(2)+wbC(ldx)-Thumid;
            enthalpy(ldx,2)=1.006*dbC(ldx)+W(ldx)*(2501+1.805*dbC(ldx));
            if enthalpy(ldx,2)>BASEenthalpy
                cooling_enthalpy_days(2)=cooling_enthalpy_days(2)+enthalpy(ldx,2)-BASEenthalpy;
                wsChD(2)=wsChD(2)+(enthalpy(ldx,2)-BASEenthalpy)*substitution(ldx,2);
            end
            shed(ldx,2)=100*(enthalpy(ldx,2)-BASEenthalpy)*substitution(ldx,2)/(enthalpy(ldx,2)-BASEenthalpy);
            
        substitution(ldx,3)=(1-(MAXwbC(ldx)+deltaTwet-ADP)/(maxC(ldx)+deltaTdry-ADP));
            wsCDD(3)=wsCDD(3)+(MAXwbC(ldx)-Thumid)*substitution(ldx,3);
            wb_cooling_degree_days(3)=wb_cooling_degree_days(3)+MAXwbC(ldx)-Thumid;
            enthalpy(ldx,3)=1.006*maxC(ldx)+Wmax(ldx)*(2501+1.805*maxC(ldx));
            if enthalpy(ldx,3)>BASEenthalpy
                cooling_enthalpy_days(3)=cooling_enthalpy_days(3)+enthalpy(ldx,3)-BASEenthalpy;
                wsChD(3)=wsChD(3)+(enthalpy(ldx,3)-BASEenthalpy)*substitution(ldx,3);
            end
            shed(ldx,3)=100*(enthalpy(ldx,3)-BASEenthalpy)*substitution(ldx,3)/(enthalpy(ldx,3)-BASEenthalpy);
   end

    % re_PC(6)
    %%%%% ventilation offers an alternative to dehumidifcation above Thumid
    %%%%% but ventilation only covers dehumidification as far as TVhot 
    %%%%% and covers part of direct expansion cooling as far as TVhumid
    % re_PC(7) and(and(Tcold<=dbC(ldx) , dbC(ldx)<=TVhot) , and(Thumid<wbC(ldx) , wbC(ldx)<=TVhumid)) 
    block1=(Thumid<wbC(ldx) && wbC(ldx)<=TVhumid && Tcold<=dbC(ldx) && dbC(ldx)<=TVhot);
    % ... OR and(and(Thot<dbC(ldx) , dbC(ldx)<=TVhot) , and(Thumid<=wbC(ldx) , wbC(ldx)<=TVhumid)); 
    block2=(Thot<dbC(ldx) && dbC(ldx)<=TVhot && Tdry<=wbC(ldx) && wbC(ldx)<=TVhumid);
    % where last block of dehumification by ventilation instead of evap cooling %%%%%
    if block1 || block2
        if isnan(ventn_days)
            ventn_days=1;
        else
            ventn_days=ventn_days+1;
        end
        avoided_dbcdd=avoided_dbcdd+dbC(ldx)-Thot;
    end
    % Tcold=18.3db and   (Tdry=9.4db   or   Thot=24db) and Tdry=9.4db
    % humidification BELOW lower wet bound of comfort (drybulb>Tcold & wetbulb<Tdry)
    if and(Tcold<=dbC(ldx) , wbC(ldx)<Tdry) % || and(Thot<=dbC(ldx) , wbC(ldx)<=Tdry)
        if isnan(humidif_days)
            humidif_days=1;
        else
            humidif_days=humidif_days+1;
        end
        % re_PC(4) and(and(Tcold<=dbC(ldx) , dbC(ldx)<=Thot) , wbC(ldx)<Tdry); % humidification
    end
    
    % note natural comfort is inclusive of bounding temperatures (dry and wet)
    if (Tcold<=dbC(ldx) && dbC(ldx)<=Thot)
        if (Tdry<=wbC(ldx) && wbC(ldx)<=Thumid)
            if isnan(natural_days)
                natural_days=1;
            else
                natural_days=natural_days+1;
            end
            % re_PC(5) and(and(Tcold<=dbC(ldx) , dbC(ldx)<=Thot) , and(Thumid>=wbC(ldx) , wbC(ldx)>=Tdry)); % comfort
        end
    end
    % heating required BELOW lower dry bound of comfort
    if dbC(ldx)<Tcold
        heating_degree_days=heating_degree_days+Tcold-dbC(ldx);
        % note icing added to non-icing heating days as another HVAC catagory
        if dbC(ldx)<Ticing && (dbC(ldx)-dpC(ldx))<deltaTice
            % re_PC(3) and(dbC(ldx)<Ticing , (dbC(ldx)-dpC(ldx))<deltaTice) ; % icing of heatpumps
            % icing is a problem for air-cooled heatpumps 
            % So consider ground source heat pumps or wood pellets
            if isnan(icing_days)
                icing_days=1;
            else
                icing_days=icing_days+1;
            end
            icing_degree_days=icing_degree_days+deltaTice+dpC(ldx)-dbC(ldx);
        else
            % re_PC(2) and(dbC(ldx)<Tcold , not(re_PC(3))); % heating without icing of heatpumps
            % heating without icing problems
            if isnan(ASHP_days)
                ASHP_days=1;
            else
                ASHP_days=ASHP_days+1;
            end
        end
    end
    % evaporative cooling separated from vapour compression and/or ventilation cooling
    % re_PC(7) and(and(Tdry<wbC(ldx) , wbC(ldx)<=Thumid) , TVhot<dbC(ldx)); % evaporative cooling
    if Tdry<wbC(ldx) && wbC(ldx)<=Thumid && TVhot<dbC(ldx)
        if isnan(evapworks_days)
            evapworks_days=1;
        else
            evapworks_days=evapworks_days+1;
        end
        avoided_dbcdd=avoided_dbcdd+dbC(ldx)-Thot;
    end
    % cooling that needs DX vapour compression work
    if TVhot<dbC(ldx) || Thumid<wbC(ldx)
        if isnan(DXcooling_days)
            DXcooling_days=1;
        else
            DXcooling_days=DXcooling_days+1; % DX required eXclusive of direct evaporative cooling opportunity
        end
    end
    if (TVhot<dbC(ldx) && Thumid<wbC(ldx)) || TVhumid<wbC(ldx)
        % cooling demand eXXclusive of both ventilation and direct evaporative cooling opportunities 
        % re_PC(8) and(and(Thot<dbC(ldx) , wbC(ldx)>Thumid),~and(dbC(ldx)<=TVhot , wbC(ldx)<=TVhumid)); %%%%%
        if isnan(DXXcooling_days)
            DXXcooling_days=1;
        else
            DXXcooling_days=DXXcooling_days+1; % eXXcluding ventilation free cooling opportunity
        end
    end
end
if timeset==0; % only the first occurence
    timeframe(1)=annum; 
    timeset=1;
end 
timeframe(2)=annum; % increments to record the latest year
re_DD(yrcnt,1)=annum; % year
re_DD(yrcnt,2)=icing_degree_days; % output icing degree days
re_DD(yrcnt,3)=heating_degree_days; % heating degree days (base 18.3 db) ***used
re_DD(yrcnt,4)=db_cooling_degree_days; % dry bulb cooling degree days (base 24 C db) ***used
re_DD(yrcnt,5)=wb_cooling_degree_days(2); % MEAN wb degree days (base 19.4 wb) ***used
re_DD(yrcnt,6)=cooling_enthalpy_days(2); % MEAN enthalpy-days
% re_DD(yrcnt,7)=100*wsCDD(1)/wb_cooling_degree_days(1); % percent MIN avoided wb degree days (base 19.4 wb) with evap condenser ***used
% re_DD(yrcnt,8)=100*wsCDD(2)/wb_cooling_degree_days(2); % percent MEAN avoided wb degree days (base 19.4 wb) with evap condenser ***used
% re_DD(yrcnt,9)=100*wsCDD(3)/wb_cooling_degree_days(3); % percent MAX avoided wb degree days (base 19.4 wb) with evap condenser ***used 

re_DD(yrcnt,7)=100*wsChD(1)/cooling_enthalpy_days(1); % percent MIN avoided enthalpy days (base 19.4 wb) with evap condenser ***used
re_DD(yrcnt,8)=100*wsChD(2)/cooling_enthalpy_days(2); % percent MEAN avoided enthalpy days (base 19.4 wb) with evap condenser ***used
re_DD(yrcnt,9)=100*wsChD(3)/cooling_enthalpy_days(3); % percent MAX avoided enthalpy days (base 19.4 wb) with evap condenser ***used 
tom=max(shed(:,1)); dick=max(shed(:,2)); harry=max(shed(:,3));
re_DD(yrcnt,10)=max([tom dick harry]); % BEST sheded enthalpy days (base 19.4 wb) with evap condenser ***used
tom=max(min(shed(:,1)),0); dick=max(min(shed(:,2)),0); harry=max(min(shed(:,3)),0); 
re_DD(yrcnt,11)=min([tom dick harry]); % LEAST sheded enthalpy days (base 19.4 wb) with evap condenser ***used 
if isnan(QA(yrcnt,7)); % not enough dewpoint data this year
    re_DD(yrcnt,5:11)=NaN; 
end 
re_DD(yrcnt,12)=ASHRAEHDD; % ASHRAE Table A3 from 2013 Standard 169 base 64.4�F = 18�C
re_DD(yrcnt,13)=ASHRAECDD; % ASHRAE Table A3 from 2013 Standard 169 base 50�F = 10�C

re_PC(yrcnt,1)=annum; % year
    icing=icing_days; % *100/daysoyear; % output heatpump icing days
re_PC(yrcnt,2)=icing; % days subject to heatpump icing
    ASHP=ASHP_days; % *100/daysoyear; /daysoyear;
re_PC(yrcnt,3)=ASHP; % days calling for ASHP heating WITHOUT ICING
    humidif=humidif_days; % *100/daysoyear; /daysoyear;
re_PC(yrcnt,4)=humidif; % days that call for humidification
    natural=natural_days; % *100/daysoyear; /daysoyear;
re_PC(yrcnt,5)=natural; % days naturally comfortable
    ventn=ventn_days; % *100/daysoyear; /daysoyear; %%%%% using ventilation instead of dehumidification
re_PC(yrcnt,6)=ventn; % days that call for ventilation "dehumidification" %%%%%
    evapworks=evapworks_days; % *100/daysoyear; /daysoyear;
re_PC(yrcnt,7)=evapworks; % days that could be served with evap cooling
    DXXcooling=DXXcooling_days; % *100/daysoyear; /daysoyear;
re_PC(yrcnt,8)=DXXcooling; % days that really must have DX cooling (eXXcluding vent & evap)
if isnan(QA(yrcnt,7)); % not enough dewpoint data this year
    re_PC(yrcnt,[2 4 7 8])=NaN; 
end
%%%%% check totals
daysoyear=icing_days+ASHP_days+humidif_days+natural_days+evapworks_days+ventn_days+DXXcooling_days;
re_PC(yrcnt,9)=daysoyear;
% allcool days db_cool and h_dehum
overThot=CdbDays; 
overThumid=CwbDays;
re_PC(yrcnt,10)=overThot; % days that drybulb calls for all types of cooling Thot=24
re_PC(yrcnt,11)=overThumid; % days that drybulb calls for all types of cooling Thumid=19.4C 67F
%%
    HVAC='????'; % be sure to have exactly four characters
    winter='__'; summer='__'; % assume nothing
    floatdays=10; % nominally 10 days per year
    if re_PC(2)>=floatdays % Too Icy in winter for heat pump
    	winter='Gs'; % gas or groundsource winter
    else
        if (re_PC(2)+re_PC(3))>=floatdays % ASHP airsource heat pump will work
            winter='Hp'; % ASHP heat pump would do the job
        else
            if (re_PC(2)+re_PC(3)+re_PC(4))>=floatdays % humidification
                winter='Hu'; % heat pump could serve evaporation load
            end
        end
    end
    if re_PC(8)>=floatdays % DXX
        summer='Ac'; % dx air-conditioning for cooling
    else
        if (re_PC(7)+re_PC(8))>=floatdays % evap cooler
            summer='Ev'; % evaporative cool summer
        else
            if (re_PC(6)+re_PC(7)+re_PC(8))>=floatdays % ventilation
                summer='Vt';
            end
        end
    end
    HVAC=strcat(winter,summer);
    %%
re_WBGT(yrcnt,1)=annum; % year
re_WBGT(yrcnt,2)=wbgt_days(1); % days BELOW WBGT 15
re_WBGT(yrcnt,3)=wbgt_days(2); % days clear WBGT 21
re_WBGT(yrcnt,4)=wbgt_days(3); % days ????? WBGT 26.6
re_WBGT(yrcnt,5)=wbgt_days(4); % days white WBGT 26.6
re_WBGT(yrcnt,6)=wbgt_days(5); % days green WBGT 27.8
re_WBGT(yrcnt,7)=wbgt_days(6); % days yellow WBGT 29.4
re_WBGT(yrcnt,8)=wbgt_days(7); % days red WBGT 31.1
re_WBGT(yrcnt,9)=wbgt_days(8); % days black WBGT 32.2
re_WBGT(yrcnt,10)=wbgt_days(9); % total days that year

if isnan(QA(yrcnt,7)); % not enough dewpoint data this year
    re_WBGT(yrcnt,2:9)=NaN; 
end

else % no results this year %%%
    re_DD(yrcnt,1)=annum; % year
    re_DD(yrcnt,2)=NaN; % MEAN modelled icing degree days
    re_DD(yrcnt,3)=NaN; % MEAN heating degree days
    re_DD(yrcnt,4)=NaN; % MEAN dry bulb cooling degree days
    re_DD(yrcnt,5)=NaN; % MEAN wb degree days
    re_DD(yrcnt,6)=NaN; % MEAN cooling enthalpy days
    re_DD(yrcnt,7)=NaN; % percent MIN substituted enthalpyDays from Carnot evap condenser
    re_DD(yrcnt,8)=NaN; % percent MEAN substituted enthalpyDays from Carnot evap condenser
    re_DD(yrcnt,9)=NaN; % percent MAX substituted enthalpyDays from Carnot evap condenser
    re_DD(yrcnt,10)=NaN; % BEST percent shed from Carnot evap condenser
    re_DD(yrcnt,11)=NaN; % LEAST percent shed from Carnot evap condenser
    re_DD(yrcnt,12)=NaN; % ASHRAE Table A3 from 2013 Standard 169 HDD base 18�C (64.4�F)
    re_DD(yrcnt,13)=NaN; % ASHRAE Table A3 from 2013 Standard 169 CDD base 10�C (50�F) 
    
    re_PC(yrcnt,1)=annum; % year
    re_PC(yrcnt,2)=NaN; % days subject to heatpump icing
    re_PC(yrcnt,3)=NaN; % days calling for ASHP heating WITHOUT ICING
    re_PC(yrcnt,4)=NaN; % days that call for humidification
    re_PC(yrcnt,5)=NaN; % days naturally comfortable
    re_PC(yrcnt,6)=NaN; % days dehumid by ventilation
    re_PC(yrcnt,7)=NaN; % days that could be served with evap cooling
    re_PC(yrcnt,8)=NaN; % days that call for DXX cooling 
    re_PC(yrcnt,9)=NaN; % check total days that sum up all HVAC modes
    re_PC(yrcnt,10)=NaN; % days that call for all sorts of cooling equal or over Thot 24C dry bulb
    re_PC(yrcnt,11)=NaN; % days that call for all sorts of dehumidification equal or over Thumid 19.4C wet bulb
    
    re_WBGT(yrcnt,1)=annum; % year
    re_WBGT(yrcnt,2)=NaN; % days BELOW WBGT 15
    re_WBGT(yrcnt,3)=NaN; % days clear WBGT 21
    re_WBGT(yrcnt,4)=NaN; % days ????? WBGT 26.6
    re_WBGT(yrcnt,5)=NaN; % days white WBGT 26.6
    re_WBGT(yrcnt,6)=NaN; % days green WBGT 27.8
    re_WBGT(yrcnt,7)=NaN; % days yellow WBGT 29.4
    re_WBGT(yrcnt,8)=NaN; % days red WBGT 31.1
    re_WBGT(yrcnt,9)=NaN; % days black WBGT 32.2
    re_WBGT(yrcnt,10)=NaN; % total days that year
    
end % with regard to this year
clearvars -except HVAC Basis significant *stats PC savings WBGTstats* statyears MM timeframe timeset ZoneBase j slash unix keypath localpath startline stepper dataArray0 fig* re_* QAfile BASEenthalpy STP deltaTice Ticing Thumid Tdry Tcold Thot dbcomfort wbcomfort gdays AMINwbC AMAXwbC AmaxC AminC AstPa AdpC AdbC AwbC AH2O C QA goodyears yrcnter startyear stopyear mm shed substitution filename zfilename altfile altZfile noaaRow SORTEDnoaapath WBANstring WMO WMOstring starter stopper annum TVhot TVhumid dbcomfort wbcomfort wbgtcolour re_WBGT re_DD re_PC yrcnt USAF WBAN STATIONNAME CTRY STATE ICAO LAT LON ELEVM BEGIN END
else end % skipping if there was not a sufficient number of records that year
end % interannual loop
%%
if not(goodyears==0) % carry on with at least one good year
Basis(4)=NaN; % #4 dewpoint-datayears-register, added to the end of Basis
Basis(4)=nnz(~(isnan(QA(:,7)))); % number of non-NaNs , in other words dewpoint data years

blanks(1:goodyears)=0; % initially assume no gaps
for stone=1:length(re_WBGT(:,2)) % from the second good year until the second before the last
    if isnan(re_WBGT(stone,2)) % then there is a gap
        blanks(stone)=1;
    else
        blanks(stone)=0;
    end
end

figure('Name',['Plots of WMO station (', strcat(num2str(WMO),') ',STATIONNAME{noaaRow})],'NumberTitle','on','Position', [10 40 560*2 420*1.5])
subplot(2,2,1)
scaleax=max(re_DD(yrcnter,3)+re_DD(yrcnter,4)); % maximum sum of drybulb heating and cooling degree days
g(7)=area(re_DD(yrcnter,1),scaleax+0*re_DD(yrcnter,1)); % plot the maximum sum of drybulb heating and cooling degree days
    set(g(7),'FaceColor',[1,0,0]); % red above db  HDD
hold on; 
g(6)=area(re_DD(yrcnter,1),scaleax-re_DD(yrcnter,3)); % plot dry bulb heating degree days ***
    set(g(6),'FaceColor',[1,1,1]); % white below
g(5)=area(re_DD(yrcnter,1),re_DD(yrcnter,4)); % drybulb cooling DD ***
    set(g(5),'FaceColor',[0,0,1]); % blue below db CDD
legend('dbHDD','','dbCDD','Location','NorthEastOutside')
if scaleax>0 
    ylim([0,scaleax]);
end
ylabel({'     < heating degrees days (winter)';'(summer) cooling degree days >     '})
title(strcat('dry bulb degree days at WMO stn ',num2str(WMO)));
% blank over NaNs
g(9)=area(re_DD(:,1),scaleax*blanks); % Totals scaleax
	set(g(9),'FaceColor',[.75,.75,.75]); % grey
%   

subplot(2,2,2)

g(6)=area(re_DD(yrcnter,1),re_DD(yrcnter,10)); % BEST percent shed Carnot evap substituted wetbulb cooling DD
    set(g(6),'FaceColor',[1,0,0]); % red
hold on;
g(5)=area(re_DD(yrcnter,1),re_DD(yrcnter,9)); % percent MAX Carnot evap substituted wetbulb cooling DD
    set(g(5),'FaceColor',[1,1,0]); % yellow
g(4)=area(re_DD(yrcnter,1),re_DD(yrcnter,8)); % percent MEAN Carnot evap substituted wetbulb cooling DD
    set(g(4),'FaceColor',[0,1,0]); % green
g(3)=area(re_DD(yrcnter,1),re_DD(yrcnter,7)); % percent MIN Carnot evap substituted wetbulb cooling DD
    set(g(3),'FaceColor',[0,0,1]); % blue
g(6)=area(re_DD(yrcnter,1),re_DD(yrcnter,11)); % LEAST percent shed Carnot evap substituted wetbulb cooling DD
    set(g(6),'FaceColor',[1,1,1]); % white
legend('best shed','average peak','average mean','average low','low shed','Location','NorthWestOutside')
title('DX cooling replaceable by evaporative condenser')  
ylim([0,100]);
% blank over NaNs
g(1)=area(re_DD(:,1),max(re_DD(yrcnter,4))*blanks); % Totals scaleax
	set(g(1),'FaceColor',[.75,.75,.75]); % grey
%
% display hours of heating, comfort, cooling, evap cooling
subplot(2,2,3)
h(1)=area(re_PC(yrcnter,1),re_PC(yrcnter,8)+re_PC(yrcnter,7)+re_PC(yrcnter,6)+re_PC(yrcnter,5)+re_PC(yrcnter,4)+re_PC(yrcnter,3)+re_PC(yrcnter,2)); % total percent days incl icing
    set(h(1),'FaceColor',[1,0.5,0]); % orange icing
hold on;
h(2)=area(re_PC(yrcnter,1),re_PC(yrcnter,8)+re_PC(yrcnter,7)+re_PC(yrcnter,6)+re_PC(yrcnter,5)+re_PC(yrcnter,4)+re_PC(yrcnter,3)); % all HVAC except icing
    set(h(2),'FaceColor',[1,0,0]); % red ASHP heating
h(3)=area(re_PC(yrcnter,1),re_PC(yrcnter,8)+re_PC(yrcnter,7)+re_PC(yrcnter,6)+re_PC(yrcnter,5)+re_PC(yrcnter,4)); % days humification over heating
    set(h(3),'FaceColor',[1,1,0]); % ); % yellow humidifcation
h(4)=area(re_PC(yrcnter,1),re_PC(yrcnter,8)+re_PC(yrcnter,7)+re_PC(yrcnter,6)+re_PC(yrcnter,5)); % days passive over cooling
    set(h(4),'FaceColor',[1,1,1]); % white passive
h(5)=area(re_PC(yrcnter,1),re_PC(yrcnter,8)+re_PC(yrcnter,7)+re_PC(yrcnter,6)); % days ventilation
    set(h(5),'FaceColor',[0,1,0]); % green dehumification by ventilation
h(6)=area(re_PC(yrcnter,1),re_PC(yrcnter,8)+re_PC(yrcnter,7)); % days evap cooling
    set(h(6),'FaceColor',[0,0,1]); % blue evap cool
h(7)=area(re_PC(yrcnter,1),re_PC(yrcnter,8)); % days DX
    set(h(7),'FaceColor',[0.5,0,0.5]); % purple DX cooling

ylabel({'days per year'})
lefty=[0 61 122 183 244 305 366];
set(gca,'YLim',[min(lefty) max(lefty)])
set(gca,'YTick',lefty)
set(gca,'YTickLabel',lefty)
% ylim([0,366]);   
    
% legend('DX cool','evap cool','ventn','passive','humidif','ASHP','icing','Location','NorthEast')
legend('icing','ASHP','humidif','passive','ventn','evap cool','DX cool','Location','NorthEastOutside')
% title('percent of daily means recommending various HVAC modes')
title('daily means recommending various HVAC modes')
% blank over NaNs
h(9)=area(re_PC(:,1),366*blanks); % 366 days per year
	set(h(9),'FaceColor',[.75,.75,.75]); % grey
%
% display WBGT colours
subplot(2,2,4)
h(8)=area(re_WBGT(yrcnter,1),re_WBGT(yrcnter,2)+re_WBGT(yrcnter,3)+re_WBGT(yrcnter,4)+re_WBGT(yrcnter,5)+re_WBGT(yrcnter,6)+re_WBGT(yrcnter,7)+re_WBGT(yrcnter,8)+re_WBGT(yrcnter,9)); % "cold" WBGT 
    set(h(8),'FaceColor',[0.5,0,0.5]); % deep purple
hold on;
h(7)=area(re_WBGT(yrcnter,1),re_WBGT(yrcnter,3)+re_WBGT(yrcnter,4)+re_WBGT(yrcnter,5)+re_WBGT(yrcnter,6)+re_WBGT(yrcnter,7)+re_WBGT(yrcnter,8)+re_WBGT(yrcnter,9)); % "cool" WBGT
    set(h(7),'FaceColor',[0,0,0.5]); % deep blue
h(6)=area(re_WBGT(yrcnter,1),re_WBGT(yrcnter,4)+re_WBGT(yrcnter,5)+re_WBGT(yrcnter,6)+re_WBGT(yrcnter,7)+re_WBGT(yrcnter,8)+re_WBGT(yrcnter,9)); % "neutral" WBGT
    set(h(6),'FaceColor',[0,0,1]); % light blue
h(5)=area(re_WBGT(yrcnter,1),re_WBGT(yrcnter,5)+re_WBGT(yrcnter,6)+re_WBGT(yrcnter,7)+re_WBGT(yrcnter,8)+re_WBGT(yrcnter,9)); % "OK" WBGT
    set(h(5),'FaceColor',[1,1,1]); % white
h(4)=area(re_WBGT(yrcnter,1),re_WBGT(yrcnter,6)+re_WBGT(yrcnter,7)+re_WBGT(yrcnter,8)+re_WBGT(yrcnter,9)); % "warm" WBGT 
    set(h(4),'FaceColor',[0,1,0]); % green
h(3)=area(re_WBGT(yrcnter,1),re_WBGT(yrcnter,7)+re_WBGT(yrcnter,8)+re_WBGT(yrcnter,9)); % "stressed" WBGT
    set(h(3),'FaceColor',[1,1,0]); % yellow
h(2)=area(re_WBGT(yrcnter,1),re_WBGT(yrcnter,8)+re_WBGT(yrcnter,9)); % "bad" WBGT 
    set(h(2),'FaceColor',[1,0,0]); % red
h(1)=area(re_WBGT(yrcnter,1),re_WBGT(yrcnter,9)); % "worst" WBGT
    set(h(1),'FaceColor',[0,0,0]); % black
title('shaded WBGT conditions at daily maximum')

ylabel({'days per year'})
lefty=[0 61 122 183 244 305 366];
set(gca,'YLim',[min(lefty) max(lefty)])
set(gca,'YTick',lefty)
set(gca,'YTickLabel',lefty)
% ylim([0,366]);

str1=strcat(' >=',num2str(wbgtcolour(1)));
str2=strcat(' >=',num2str(wbgtcolour(2)));
str3=strcat(' >=',num2str(wbgtcolour(3)));
str4=strcat(' >=',num2str(wbgtcolour(4)));
str5=strcat(' >=',num2str(wbgtcolour(5)));
str6=strcat(' >=',num2str(wbgtcolour(6)));
str7=strcat(' >=',num2str(wbgtcolour(7)));
legend('colder ',str1,str2,str3,str4,str5,str6,str7,'Location','NorthWestOutside')
% blank over NaNs
h(9)=area(re_WBGT(:,1),366*blanks); % should total up to 366 days
	set(h(9),'FaceColor',[.75,.75,.75]); % grey

print(fig1PNG,'-dpng'); % print(fig1EPS,'-depsc'); print(fig1WMF,'-dmeta'); % save as png %, eps, and windows metafile
close(gcf); % close current figure
else
    goodyears
end
    %%
noaaRow
WBAN(noaaRow)
STATIONNAME{noaaRow};
% CTRY{noaaRow}
% STATE{noaaRow}
% LAT(noaaRow)
% LON(noaaRow)
% ELEVM(noaaRow)
% BEGIN(noaaRow)
% END(noaaRow)

%%
if not(gdays==0) % there was at least one day of data

% AstPa
% AdbC
% AH2O
% AwbC
% AdpC

gidx = and(and(Tcold<=AdbC , AdbC<=Thot) , and(Thumid>=AwbC , AwbC>=Tdry)); % comfort
nidx = and(Tcold<=AdbC , AwbC<Tdry); % humidification
% "block 1" (dehumification and cooling by ventilation) ... which has an enthalpy savings compared to vapour compression cooling
% Tcold<=AdbC<=TVhot AND Thumid<AwbC<=TVhumid
% OR "block 2" (cooling by ventilation) ... which is NOT an enthalpy-duty
% Tdry<=AwbC<=Thumid AND Thot<AdbC<=TVhot % dehumification by ventilation %%%%%
pidx = or(and(and(Tcold<=AdbC,AdbC<=TVhot),and(Thumid<AwbC,AwbC<=TVhumid)),and(and(Tdry<=AwbC,AwbC<=Thumid),and(Thot<AdbC,AdbC<=TVhot))); % dehumification AND/OR cooling by ventilation %%%%%
iidx = and(AdbC<Ticing , (AdbC-AdpC)<deltaTice) ; % icing of heatpumps
hidx = and(AdbC<Tcold , not(iidx)); % ASHP heating without icing of heatpumps
eidx = and(and(Tdry<=AwbC , AwbC<=Thumid) , AdbC>TVhot); % evaporative cooling
didx = and(and(Thot<AdbC , AwbC>Thumid),~and(AdbC<=TVhot , AwbC<=TVhumid)); %%%%%
%%%%% didx is now the vapour-compression cooling eXXcluding dehum by ventilation and evap cooling

% prepare saturation curve
% kilo Pascals average pressure of dataset
kPa_max=max(AstPa)/1000;
kPa_average=mean(AstPa)/1000;
kPa_min=min(AstPa)/1000;
kPa_q1=median(AstPa)/1000;
kPa_q1=quantile(AstPa,.25)/1000;
kPa_median=median(AstPa)/1000;
kPa_q3=quantile(AstPa,.75)/1000;
pout=[kPa_max kPa_average kPa_min; kPa_q1 kPa_median kPa_q3] % output
minx=floor(min(AdbC))-3; maxx=ceil(max(AdbC));xstep=1;
xsteps=maxx-minx; % number of one degree marks on x-axis
kount=0; % intialize going into loop
while kount<=xsteps
    % work out saturation curve versus dry bulb
    kount=kount+1;
    db_axis(kount)=minx+(kount-1)*xstep;
    dbA=db_axis(kount) +273.15; % absolute temperature Kelvins of x-axis increment
    if db_axis(kount)<0 % below freezing
        lnpwsdp = C(1)/dbA + C(2) + C(3)*dbA + C(4)*dbA^2 + C(5)*dbA^3 + C(6) * dbA^4 + C(7)*log(dbA);
        pwsxdp=exp(lnpwsdp); % Pascals vapour pressure at saturation with db
    else % above freezing
        lnpwsdp = C(8)/dbA + C(9) + C(10)*dbA + C(11)*dbA^2 + C(12)*dbA^3 + C(13)*log(dbA);
        pwsxdp=exp(lnpwsdp); % Pascals vapour pressure at saturation with db
    end % now we have vapour pressure at DP saturation
    
    AH2O_sat_average(kount)=0.62198*pwsxdp/(1000*kPa_average-pwsxdp);
    AH2O_sat_q1(kount)=0.62198*pwsxdp/(1000*kPa_q1-pwsxdp);
    AH2O_sat_median(kount)=0.62198*pwsxdp/(1000*kPa_median-pwsxdp);
    AH2O_sat_q3(kount)=0.62198*pwsxdp/(1000*kPa_q3-pwsxdp);
    db_desicated_average(kount) = AH2O_sat_average(kount)*(2501-2.326*db_axis(kount))/1.006+db_axis(kount);
    db_desicated_q1(kount) = AH2O_sat_q1(kount)*(2501-2.326*db_axis(kount))/1.006+db_axis(kount);
    db_desicated_median(kount) = AH2O_sat_median(kount)*(2501-2.326*db_axis(kount))/1.006+db_axis(kount);
    db_desicated_q3(kount) = AH2O_sat_q3(kount)*(2501-2.326*db_axis(kount))/1.006+db_axis(kount);
end
%%
if max(AH2O)>0
strF2=['Psychrometric Chart for WMO station (', strcat(num2str(WMO),') ',STATIONNAME{noaaRow})];
figure('Name',strF2,'NumberTitle','on','Position', [60 40 560*1.5 420*1.5])
% Add graphic of saturation curve alternatives
hgh(1)=plot(db_axis,1000*AH2O_sat_average);
hold on
hgh(2)=plot(db_axis,1000*AH2O_sat_q1);
hgh(3)=plot(db_axis,1000*AH2O_sat_median);
hgh(4)=plot(db_axis,1000*AH2O_sat_q3);
% add to the graphic the median wetbulb lines
for index =1:xsteps;
    % plot([db_axis(index) db_desicated_average(index)],[1000*AH2O_sat_average(index) 0]);
    hgh(index)=plot([db_axis(index) db_desicated_median(index)],[1000*AH2O_sat_median(index) 0], 'k');
end
% plot data points on psychromtric chart now that we have added wetbulb lines
hg(1)=scatter(AdbC(gidx),AH2O(gidx), 'o'); % passive comfort black
    set(hg(1),'MarkerEdgeColor',[0,0,0]); % black 1
hg(2)=scatter(AdbC(nidx),AH2O(nidx), '.'); % humidification cyan
    set(hg(2),'MarkerEdgeColor',[1,1,0]); % light purple 2
hg(3)=scatter(AdbC(pidx),AH2O(pidx), '^'); % ventilation green
    set(hg(3),'MarkerEdgeColor',[0,1,0]); % dark green 3 
hg(4)=scatter(AdbC(hidx),AH2O(hidx), '>'); % ASHP heating red
    set(hg(4),'MarkerEdgeColor',[1,0,0]); % red 4
hg(5)=scatter(AdbC(iidx),AH2O(iidx), '<'); % icing orange
    set(hg(5),'MarkerEdgeColor',[1,0.5,0]); % orange 5
hg(6)=scatter(AdbC(eidx),AH2O(eidx), 'o');% evaporative cooling blue
    set(hg(6),'MarkerEdgeColor',[0,0,1]); % blue 6
hg(7)=scatter(AdbC(didx),AH2O(didx), '*'); % vapour-compression cooling purple
    set(hg(7),'MarkerEdgeColor',[0.5,0,0.5]); % purple 7
vectledge=[length(AdbC(gidx))>0 length(AdbC(nidx))>0 length(AdbC(pidx))>0 length(AdbC(hidx))>0 length(AdbC(iidx))>0 length(AdbC(eidx))>0 length(AdbC(didx))>0];
clear pvector
if isequal(vectledge,[    1     1     1     1     1     1     1]) % all modes
    pvector=         [hg(1) hg(2) hg(3) hg(4) hg(5) hg(6) hg(7)];
    legend(pvector, 'passive','humidification','ventilation','ASHP','icing','evaporative cool','vapour compression cooling','Location','northwest'); end
if isequal(vectledge,[    1     0     1     1     0     0     1]) % heat, passive, ventn, and DX
    pvector=         [hg(1)       hg(3) hg(4)             hg(7)];
    legend(pvector,'passive','ventilation','ASHP','vapour compression cooling','Location','northwest'); end
if isequal(vectledge,[    1     1     1     1     1     1     0]) % no DX cooling
    pvector=         [hg(1) hg(2) hg(3) hg(4) hg(5) hg(6)      ];
    legend(pvector,'passive','humidification','ventilation','ASHP','icing','evaporative cool','Location','northwest'); end
if isequal(vectledge,[    1     1     1     1     1     0     1]) % no evap cooling
    pvector=         [hg(1) hg(2) hg(3) hg(4) hg(5)        hg(7)];
    legend(pvector,'passive','humidification','ventilation','ASHP','icing','vapour compression cooling','Location','northwest'); end
if isequal(vectledge,[    1     1     1     1     1     0     0]) % no DX cooling, nor evap cooling
    pvector=         [hg(1) hg(2) hg(3) hg(4) hg(5)            ];
    legend(pvector,'passive','humidification','ventilation','ASHP','icing','Location','northwest'); end
if isequal(vectledge,[    1     1     1     1     0     0     1]) % all modes except icing, nor evaporative cool
    pvector=         [hg(1) hg(2) hg(3) hg(4)             hg(7)];
    legend(pvector,'passive','humidification','ventilation','ASHP','vapour compression cooling','Location','northwest'); end
if isequal(vectledge,[    1     0     1     1     0     1     1]) % all modes but icing nor humidification
    pvector=         [hg(1)       hg(3) hg(4)        hg(6) hg(7)];
    legend(pvector,'passive','ventilation','ASHP','evaporative cool','vapour compression cooling','Location','northwest'); end
if isequal(vectledge,[    1     1     1     1     0     1     1]) % all modes but icing
    pvector=         [hg(1) hg(2) hg(3) hg(4)       hg(6) hg(7)];
    legend(pvector,'passive','humidification','ventilation','ASHP','evaporative cool','vapour compression cooling','Location','northwest'); end
if isequal(vectledge,[    1     1     1     1     0     1     0]) % all modes but icing nor DX
    pvector=         [hg(1) hg(2) hg(3) hg(4)       hg(6)      ];
    legend(pvector,'passive','humidification','ventilation','ASHP','evaporative cool','Location','northwest'); end
if isequal(vectledge,[    1     1     1     1     0     0     0]) % all modes but icing, nor evaporative cool, nor DX cooling
    pvector=         [hg(1) hg(2) hg(3) hg(4)                  ];
    legend(pvector,'passive','humidification','ventilation','ASHP','Location','northwest'); end
if isequal(vectledge,[    1     1     1     0     0     0     1]) % no icing (nor heating), nor evaporative cool
    pvector=         [hg(1) hg(2) hg(3)                   hg(7)];
    legend(pvector,'passive','humidification','ventilation','vapour compression cooling','Location','northwest'); end
if isequal(vectledge,[    1     1     1     0     0     1     0]) % no icing (nor heating), nor DX cooling
    pvector=         [hg(1) hg(2) hg(3)             hg(6)      ];
    legend(pvector,'passive','humidification','ventilation','Location','northwest'); end
if isequal(vectledge,[    1     1     1     0     0     1     1]) % no icing (nor ASHP heating)
    pvector=         [hg(1) hg(2) hg(3)             hg(6) hg(7)];
    legend(pvector,'passive','humidification','ventilation','evaporative cool','vapour compression cooling','Location','northwest'); end
if isequal(vectledge,[    1     1     1     0     1     1     1]) % no heating
    pvector=         [hg(1) hg(2) hg(3)       hg(5) hg(6) hg(7)];
    legend(pvector,'passive','humidification','ventilation','icing','evaporative cool','vapour compression cooling','Location','northwest'); end
if isequal(vectledge,[    1     1     0     1     1     0     0]) % no DX cooling, nor evap cooling, nor dehumidification
    pvector=         [hg(1) hg(2)       hg(4) hg(5)            ];
    legend(pvector,'passive','humidification','ASHP','icing','Location','northwest'); end
if isequal(vectledge,[    1     0     1     0     0     1     1]) % no icing (nor ASHP heating), nor humidification
    pvector=         [hg(1)       hg(3)             hg(6) hg(7)];
    legend(pvector,'passive','ventilation','evaporative cool','vapour compression cooling','Location','northwest'); end
if isequal(vectledge,[    1     0     0     1     1     1     0]) % no humidification, dehum, nor vapour compression cooling
    pvector=         [hg(1)             hg(4) hg(5) hg(6)      ];
    legend(pvector,'passive','ASHP','icing','evaporative cool','Location','northwest'); end
if isequal(vectledge,[    1     0     1     1     1     1     1]) % no humidification
    pvector=         [hg(1)       hg(3) hg(4) hg(5) hg(6) hg(7)];
    legend(pvector,'passive','ventilation','ASHP','icing','evaporative cool','vapour compression cooling','Location','northwest'); end
if isequal(vectledge,[    1     0     1     1     1     0     0]) % no humidification, evap cooling, nor DX cooling
    pvector=         [hg(1)       hg(3) hg(4) hg(5)            ];
    legend(pvector,'passive','ventilation','ASHP','icing','Location','northwest'); end
if isequal(vectledge,[    1     0     0     1     1     0     0]) % only comfort, ASHP heating, and icing
    pvector=         [hg(1)             hg(4) hg(5)            ];
    legend(pvector,'passive','ASHP','icing','Location','northwest'); end
if isequal(vectledge,[    0     1     1     1     1     1     1]) % no comfort
    pvector=         [       hg(2) hg(3) hg(4) hg(5) hg(6) hg(7)];
    legend(pvector,'humidification','ventilation','ASHP','icing','evaporative cool','vapour compression cooling','Location','northwest'); end
if isequal(vectledge,[    0     0     0     1     1     0     0]) % only ASHP heating and icing
    pvector=         [                  hg(4) hg(5)            ];
    legend(pvector,'ASHP','icing','Location','northwest'); end
if isequal(vectledge,[    0     0     1     0     0     0     1]) % dehumidification and vapour compression cooling
    pvector=         [            hg(3)                   hg(7)];
    legend(pvector,'ventilation','vapour compression cooling','Location','northwest'); end
if isequal(vectledge,[    1     0     1     0     0     0     1]) % only passive, ventn & DX (Port Moresby)
    pvector=         [hg(1)       hg(3)                   hg(7)];
    legend(pvector,'passive','ventilation','vapour compression cooling','Location','northwest'); end

if isequal(vectledge,[    0     0     0     0     0     0     0]); end % no legend


axis([minx maxx 0 max(AH2O)]);
title(strF2)
xlabel('dry bulb (�C)') % x-axis label
ylabel('humidity (g H20/kg air)') % y-axis label

print(fig2PNG,'-dpng');% print(fig2EPS,'-depsc'); print(fig2WMF,'-dmeta'); % save as png %, eps, and windows metafile
end % selective print of figure 2
% output tabular data 
% quality assurance of years
fid = fopen(QAfile,'w') % opens file for writing 'w'
fprintf(fid,'%s\n','WMO_id _WBAN_ YEAR StnPres dys T_f D_f P_f min max ppt');
for dealer=1:size(QA,1)     
	fprintf(fid,'%6g %5g %4g %7.4f %3g %3g %3g %3g %3g %3g %3g\n',QA(dealer,:));
end
fclose(fid); clear dealer
%%
% WBGT
fid = fopen(re_WBGT_file,'w'); % opens file for writing 'w'
fprintf(fid,'YEAR  colder >=%4.1f >=%4.1f >=%4.1f >=%4.1f >=%4.1f >=%4.1f >=%4.1f totals\n',wbgtcolour);
for dealer=1:size(re_WBGT,1)     
	fprintf(fid,'%5g %6.0f %6.0f %6.0f %6.0f %6.0f %6.0f %6.0f %6.0f %6.0f\n',re_WBGT(dealer,:));
end
fclose(fid); clear dealer
%%
% Degree days DD
fid = fopen(re_DD_file,'w'); % opens file for writing 'w'
fprintf(fid,'YEAR  icyDD%1.0f HDD%3.1f CDD%3.1f wbDD%3.1f hD%4.1fJpg minEv  averEv   maxEv  bestEv leastEv HDD%4.1f CDD%4.1f\n',Ticing, Tcold, Thot, Thumid, BASEenthalpy, ZoneBase);
for dealer=1:size(re_DD,1)     
	fprintf(fid,'%4.0f %7.0f %7.0f %7.0f %7.0f %7.0f %7.2f %7.2f %7.2f %7.2f %7.2f %7.0f %7.0f\n',re_DD(dealer,:));
end
fclose(fid); clear dealer
%%
% percent modes HVAC
fid = fopen(re_PC_file,'w'); % opens file for writing 'w'
fprintf(fid,'%s \n','YEAR  icy heat humi comf vent evap  DXX days drycool dehum HVAC');
for dealer=1:size(re_PC,1)     
	fprintf(fid,'%4g %4g %4g %4g %4g %4g %4g %4g %4g    %4g  %4g %s\n',re_PC(dealer,:),HVAC);
end
fclose(fid); clear dealer

if ~isnan(max(re_DD(:,6)))
    ChDmdl=fitlm(re_DD(:,1),re_DD(:,6),'constant'); % cooling_enthalpy_days taken from the MEAN time of each day
    ChDstats(1)=double(int32(mean(ChDmdl.Fitted))); % constant estimate
    ChDstats(2)=double(int32(ChDmdl.RMSE)); % standard deviation
    ChDmdl = fitlm(re_DD(:,1),re_DD(:,6),'linear'); % assuming ChD changes with time
    ChDstats(3)=ChDmdl.Coefficients{2,1}; % slope per year
    ChDstats(4)=double(int32(ChDmdl.RMSE)); % Standard Estimate of Error a.k.a. root mean square error (STEYX)
    ChDstats(5)=ChDmdl.Rsquared.Adjusted; % R� adjusted correlation coefficient
    ChDstats(6)=ChDmdl.NumObservations; % n= years
    ChDstats(8)=ChDmdl.Coefficients.pValue(2);% values of pDW>0.05 suggests that the residuals are not autocorrelated.
    [ChDstats(10),ChDstats(9)] = dwtest(ChDmdl,'exact','both'); % two-sided Durbin-Watson test to determine if there is any autocorrelation among the residuals of the linear model
    if ChDstats(8)<=significant && ChDstats(6)>=statyears % slope regression <= 0.05 is significant (true or false) % subject to >=statyears
        ChDstats(7)=double(int32(predict(ChDmdl,max(re_DD(:,1))))); % +ChDmdl.RMSE)); % linear regression at end of timeline %%NOT+SEE%%
    else
        ChDstats(7)=NaN; % ChDstats(1)+ChDstats(2); % constant estimate + standard deviation as design guideline
    end
else
    ChDstats(1:10)=NaN;
end

if ~isnan(max(re_DD(:,11)))
    LEASTmdl=fitlm(re_DD(:,1),re_DD(:,11),'constant'); % assuming stationary percent LEAST avoided enthalpy days (base 19.4 wb) with evap condenser
    savings(1)=double(int32(mean(LEASTmdl.Fitted))); % constant estimate
else     savings(1)=NaN; end

if ~isnan(max(re_DD(:,7)))
    MINmdl=fitlm(re_DD(:,1),re_DD(:,7),'constant'); % assuming stationary percent MIN avoided enthalpy days (base 19.4 wb) with evap condenser
    savings(2)=double(int32(mean(MINmdl.Fitted))); % constant estimate
else     savings(2)=NaN; end

if ~isnan(max(re_DD(:,8)))
    MEANmdl=fitlm(re_DD(:,1),re_DD(:,8),'constant'); % assuming stationary percent MEAN avoided enthalpy days (base 19.4 wb) with evap condenser
    savings(3)=double(int32(mean(MEANmdl.Fitted))); % constant estimate
else     savings(3)=NaN; end

if ~isnan(max(re_DD(:,9)))
    MAXmdl=fitlm(re_DD(:,1),re_DD(:,9),'constant'); % assuming stationary percent MAX avoided enthalpy days (base 19.4 wb) with evap condenser
    savings(4)=double(int32(mean(MAXmdl.Fitted))); % constant estimate
else     savings(4)=NaN; end

if ~isnan(max(re_DD(:,10)))
    BESTmdl=fitlm(re_DD(:,1),re_DD(:,10),'constant'); % assuming stationary percent BEST avoided enthalpy days (base 19.4 wb) with evap condenser
    savings(5)=double(int32(mean(BESTmdl.Fitted))); % constant estimate
else     savings(5)=NaN; end

% assuming stationary HVAC days per year, estimate constant days per year assignment of PC(1:10)
if ~isnan(max(re_PC(:,9)))
    mdl=fitlm(re_DD(:,1),re_PC(:,9),'constant'); 
    PC(1)=double(int32(mean(mdl.Fitted))); % (1) check total days that sum up all HVAC modes
else     PC(1)=NaN; end
if ~isnan(max(re_PC(:,2)))
    mdl=fitlm(re_DD(:,1),re_PC(:,2),'constant'); 
    PC(2)=double(int32(mean(mdl.Fitted))); % (2) days subject to heatpump icing
else     PC(2)=NaN; end
if ~isnan(max(re_PC(:,3)))
    mdl=fitlm(re_DD(:,1),re_PC(:,3),'constant'); 
    PC(3)=double(int32(mean(mdl.Fitted))); % (3) days calling for ASHP heating WITHOUT ICING
else     PC(3)=NaN; end
if ~isnan(max(re_PC(:,4)))
    mdl=fitlm(re_DD(:,1),re_PC(:,4),'constant'); 
    PC(4)=double(int32(mean(mdl.Fitted))); % (4) days that call for humidification
else     PC(4)=NaN; end %%%%%% humidification
if ~isnan(max(re_PC(:,5)))
    mdl=fitlm(re_DD(:,1),re_PC(:,5),'constant'); 
    PC(5)=double(int32(mean(mdl.Fitted))); % (5) days naturally comfortable
else     PC(5)=NaN; end
if ~isnan(max(re_PC(:,6)))
    mdl=fitlm(re_DD(:,1),re_PC(:,6),'constant'); 
    PC(6)=double(int32(mean(mdl.Fitted))); % (6) days dehumid by ventilation 
else     PC(6)=NaN; end
if ~isnan(max(re_PC(:,7)))
    mdl=fitlm(re_DD(:,1),re_PC(:,7),'constant'); 
    PC(7)=double(int32(mean(mdl.Fitted))); % (7) days that could be served with evap cooling
else     PC(7)=NaN; end
if ~isnan(max(re_PC(:,8)))
    mdl=fitlm(re_DD(:,1),re_PC(:,8),'constant'); 
    PC(8)=double(int32(mean(mdl.Fitted))); % (8) days that call for DXX cooling 
else     PC(8)=NaN; end

checkhours=NaN;
for checkstep=2:8 % check sum of icy:DXX
    if~isnan(PC(checkstep))
        if(isnan(checkhours))
            checkhours=PC(checkstep);
        else
            checkhours=checkhours+PC(checkstep);
        end
    end
end

for checkstep=2:8 % adjust icy:DXX by factor=(total days/checkhours)
    if~isnan(PC(checkstep))
        if and(checkhours>0,~isnan(PC(1)))
            PC(checkstep)=PC(checkstep)*PC(1)/checkhours;
        else
            PC(checkstep)=NaN;
        end
    end
end

if ~isnan(max(re_PC(:,10)))
    mdl=fitlm(re_DD(:,1),re_PC(:,10),'constant'); 
    PC(9)=double(int32(mean(mdl.Fitted))); % (9) days that call for all sorts of cooling over Thot 24 db
else     PC(9)=NaN; end
if ~isnan(max(re_PC(:,11)))
    mdl=fitlm(re_DD(:,1),re_PC(:,11),'constant'); 
    PC(10)=double(int32(mean(mdl.Fitted))); % (10) days that call for all sorts of dehumidification over Thumid 19.4 wb
else     PC(10)=NaN; end


if ~isnan(max(re_PC(:,8)))
    DXmdl = fitlm(re_PC(:,1),re_PC(:,8),'constant'); % assuming stationary DX days per year
    DXstats(1)=double(int32(mean(DXmdl.Fitted))); % constant estimate
    DXstats(2)=double(int32(DXmdl.RMSE)); % standard deviation
    DXmdl = fitlm(re_PC(:,1),re_PC(:,8),'linear'); % assuming DX changes with time
    DXstats(3)=DXmdl.Coefficients{2,1}; % slope per year
    DXstats(4)=double(int32(DXmdl.RMSE)); % Standard Estimate of Error
    DXstats(5)=DXmdl.Rsquared.Adjusted; % R� adjusted correlation coefficient
    DXstats(6)=DXmdl.NumObservations; % n= years
    DXstats(8)=DXmdl.Coefficients.pValue(2);% values of pDW>0.05 suggests that the residuals are not autocorrelated.
    [DXstats(10),DXstats(9)] = dwtest(DXmdl,'exact','both'); % two-sided Durbin-Watson test to determine if there is any autocorrelation among the residuals of the linear model
    if DXstats(8)<=significant && DXstats(6)>=statyears % slope regression <= 0.05 is significant (true or false) % subject to >=statyears
        DXstats(7)=double(int32(predict(DXmdl,max(re_PC(:,1))))); % +DXmdl.RMSE)); % linear regression at end of timeline %%NOT+SEE%%
    else
        DXstats(7)=NaN; % DXstats(1)+DXstats(2); % constant estimate + standard deviation as design guideline
    end
else
    DXstats(1:10)=NaN;
end
WBGTwhite = re_WBGT(:,5)+ re_WBGT(:,6)+ re_WBGT(:,7)+ re_WBGT(:,8)+re_WBGT(:,9); % WBGT (@ Tmax) > 29.4C sum of white, green, yellow, red and black WBGT days
if ~isnan(max(WBGTwhite))
    WBGTmdlW = fitlm(re_WBGT(:,1),WBGTwhite,'constant'); % assuming stationary WBGT stressed days per year
    WBGTstatsW(1)=double(int32(mean(WBGTmdlW.Fitted))); % constant estimate
    WBGTstatsW(2)=double(int32(WBGTmdlW.RMSE)); % standard deviation
    WBGTmdlW = fitlm(re_WBGT(:,1),WBGTwhite,'linear'); % assuming WBGT changes with time
    WBGTstatsW(3)=WBGTmdlW.Coefficients{2,1}; % slope per year
    WBGTstatsW(4)=double(int32(WBGTmdlW.RMSE)); % Standard Estimate of Error
    WBGTstatsW(5)=WBGTmdlW.Rsquared.Adjusted; % R� adjusted correlation coefficient
    WBGTstatsW(6)=WBGTmdlW.NumObservations; % n= years
    WBGTstatsW(8)=WBGTmdlW.Coefficients.pValue(2);% values of pDW>0.05 suggests that the residuals are not autocorrelated.
    [WBGTstatsW(10),WBGTstatsW(9)] = dwtest(WBGTmdlW,'exact','both'); % two-sided Durbin-Watson test to determine if there is any autocorrelation among the residuals of the linear model
    if WBGTstatsW(8)<=significant && WBGTstatsW(6)>=statyears % slope regression <= 0.05 is significant (true or false) % subject to >=statyears
        WBGTstatsW(7)=double(int32(predict(WBGTmdlW,max(re_PC(:,1))))); % +WBGTmdlW.RMSE)); % linear regression at end of timeline %%NOT+SEE%%
    else
        WBGTstatsW(7)=NaN; % WBGTstatsW(1)+WBGTstatsW(2); % constant estimate + standard deviation as design guideline
    end
else
    WBGTstatsW(1:10)=NaN;
end

WBGTyellow = re_WBGT(:,7)+ re_WBGT(:,8)+re_WBGT(:,9); % WBGT (@ Tmax) > 29.4C sum of yellow, red and black WBGT days
if ~isnan(max(WBGTyellow))
    WBGTmdlY = fitlm(re_WBGT(:,1),WBGTyellow,'constant'); % assuming stationary WBGT stressed days per year
    WBGTstatsY(1)=double(int32(mean(WBGTmdlY.Fitted))); % constant estimate
    WBGTstatsY(2)=double(int32(WBGTmdlY.RMSE)); % standard deviation
    WBGTmdlY = fitlm(re_WBGT(:,1),WBGTyellow,'linear'); % assuming WBGT changes with time
    WBGTstatsY(3)=WBGTmdlY.Coefficients{2,1}; % slope per year
    WBGTstatsY(4)=double(int32(WBGTmdlY.RMSE)); % Standard Estimate of Error
    WBGTstatsY(5)=WBGTmdlY.Rsquared.Adjusted; % R� adjusted correlation coefficient
    WBGTstatsY(6)=WBGTmdlY.NumObservations; % n= years
    WBGTstatsY(8)=WBGTmdlY.Coefficients.pValue(2);% values of pDW>0.05 suggests that the residuals are not autocorrelated.
    [WBGTstatsY(10),WBGTstatsY(9)] = dwtest(WBGTmdlY,'exact','both'); % two-sided Durbin-Watson test to determine if there is any autocorrelation among the residuals of the linear model
    if WBGTstatsY(8)<=significant && WBGTstatsY(6)>=statyears % slope regression <= 0.05 is significant (true or false) % subject to >=statyears
        WBGTstatsY(7)=double(int32(predict(WBGTmdlY,max(re_PC(:,1))))); % +WBGTmdlY.RMSE)); % linear regression at end of timeline %%NOT+SEE%%
    else
        WBGTstatsY(7)=NaN; % WBGTstatsY(1)+WBGTstatsY(2); % constant estimate + standard deviation as design guideline
    end
else
    WBGTstatsY(1:10)=NaN;
end
    
if ~isnan(max(re_DD(:,12)))
    ashraeHDDmdl = fitlm(re_DD(:,1),re_DD(:,12),'constant'); % ASHRAE stationary HDD18
    ashraeHDDstats(1)=int32(mean(ashraeHDDmdl.Fitted));
else
    ashraeHDDstats(1)=NaN;
end
if ~isnan(max(re_DD(:,3)))
    HDDmdl = fitlm(re_DD(:,1),re_DD(:,3),'constant'); % assuming stationary HDD18.3
    HDDstats(1)=double(int32(mean(HDDmdl.Fitted))); % constant estimate
    HDDstats(2)=double(int32(HDDmdl.RMSE)); % standard deviation
    HDDmdl = fitlm(re_DD(:,1),re_DD(:,3),'linear'); % assuming HDD changes with time
    HDDstats(3)=HDDmdl.Coefficients{2,1}; % slope per year
    HDDstats(4)=double(int32(HDDmdl.RMSE)); % Standard Estimate of Error
    HDDstats(5)=HDDmdl.Rsquared.Adjusted; % R� adjusted correlation coefficient
    HDDstats(6)=HDDmdl.NumObservations; % n= years 
    HDDstats(8)=HDDmdl.Coefficients.pValue(2);% values of pDW>0.05 suggests that the residuals are not autocorrelated.
    [HDDstats(10),HDDstats(9)] = dwtest(HDDmdl,'exact','both'); % two-sided Durbin-Watson test to determine if there is any autocorrelation among the residuals of the linear model
    if HDDstats(8)<=significant && HDDstats(6)>=statyears % slope regression <= 0.05 is significant (true or false) % subject to >=statyears
        HDDstats(7)=double(int32(predict(HDDmdl,max(re_DD(:,1))))); % +HDDmdl.RMSE)); % linear regression at end of timeline %%NOT+SEE%%
    else
        HDDstats(7)=NaN; % HDDstats(1)+HDDstats(2); % constant estimate + standard deviation as design guideline
    end
else
    HDDstats(1:10)=NaN;
end

if ~isnan(max(re_DD(:,13)))
    ashraeCDDmdl = fitlm(re_DD(:,1),re_DD(:,13),'constant'); % ASHRAE stationary CDD10
    ashraeCDDstats(1)=int32(mean(ashraeCDDmdl.Fitted));
else
    ashraeCDDstats(1)=NaN;
end
if ~isnan(max(re_DD(:,4)))
    CDDmdl = fitlm(re_DD(:,1),re_DD(:,4),'constant'); % assuming stationary CDD24
    CDDstats(1)=double(int32(mean(CDDmdl.Fitted))); %constant estimate 
    CDDstats(2)=double(int32(CDDmdl.RMSE)); % standard deviation
    CDDmdl = fitlm(re_DD(:,1),re_DD(:,4),'linear'); % assuming CDD changes with time
    CDDstats(3)=CDDmdl.Coefficients{2,1}; % slope per year
    CDDstats(4)=double(int32(CDDmdl.RMSE)); % Standard Estimate of Error
    CDDstats(5)=CDDmdl.Rsquared.Adjusted; % R� adjusted correlation coefficient
    CDDstats(6)=CDDmdl.NumObservations; % n= years
    CDDstats(8)=CDDmdl.Coefficients.pValue(2);% values of pDW>0.05 suggests that the residuals are not autocorrelated.
    [CDDstats(10),CDDstats(9)] = dwtest(CDDmdl,'exact','both'); % two-sided Durbin-Watson test to determine if there is any autocorrelation among the residuals of the linear model
    if CDDstats(8)<=significant && CDDstats(6)>=statyears % slope regression <= 0.05 is significant (true or false) % subject to >=statyears
        CDDstats(7)=double(int32(predict(CDDmdl,max(re_DD(:,1))))); % +CDDmdl.RMSE)); % linear regression at end of timeline %%NOT+SEE%%
    else
        CDDstats(7)=NaN; % CDDstats(1)+CDDstats(2); % constant estimate + standard deviation as design guideline
    end
else
    CDDstats(1:10)=NaN;
end

ASHRAEcz=NaN;
if ashraeHDDstats(1)<=2000
    if ashraeCDDstats(1)>6000
        ASHRAEcz=0;
    else
        if ashraeCDDstats(1)>5000
            ASHRAEcz=1;
        else
            if ashraeCDDstats(1)>3500
                ASHRAEcz=2;
            else
                if ashraeHDDstats(1)>3500
                    ASHRAEcz=2;
                else
                    ASHRAEcz=3;
                end
            end
        end
    end
else
    if ashraeHDDstats(1)<=3000 && ashraeCDDstats(1)<=3500
         ASHRAEcz=4;
    else
        if ashraeHDDstats(1)<=4000 && ashraeCDDstats(1)<=3500
            ASHRAEcz=5;
        else
            if ashraeHDDstats(1)<=5000 && ashraeCDDstats(1)<=2500
                ASHRAEcz=6;
            else
                if ashraeHDDstats(1)<=7000 && ashraeCDDstats(1)<=2500
                    ASHRAEcz=7;
                else
                    if ashraeCDDstats(1)<=2500
                        ASHRAEcz=8;
                    end
                end
            end
        end
    end
end
% v16b maximum mean monthly MMM vector over the years ... vector for 12 months
MMM=max(MM(starter:stopper,:,1));

else
    gdays % nothing to report
    savings(1:5)=NaN;
    HDDstats(1:10)=NaN;
    CDDstats(1:10)=NaN;
    ASHRAEcz=NaN;
    MMM(1:12)=NaN;
    DXstats(1:10)=NaN;
    WBGTstatsW(1:10)=NaN;
    WBGTstatsY(1:10)=NaN;
    ChDstats(1:10)=NaN;
    PC(1:10)=NaN;
end
close(gcf); % close current figure
% clear fig* re_* QAfile BASEenthalpy STP gdays AMINwbC AMAXwbC AmaxC AminC AstPa AdpC AdbC AwbC AH2O C QA goodyears yrcnter mm substitution filename zfilename altfile altZfile noaaRow SORTEDnoaapath cases WBANstring WMOstring starter stopper annum re_WBGTstatsY re_DD re_PC yrcnt USAF WBAN STATIONNAME CTRY STATE ICAO LAT LON ELEVM BEGIN END
end

