clear
close all


% Wind data source:
% https://www.ncdc.noaa.gov/cdo-web/datasets/GHCND/stations/GHCND:USW00094889/detail
% Have to add the thing to your cart and then order it. Comes via email. 

% In the .csv file, AWND = average wind speed.
% From the documentation at: https://www.ncei.noaa.gov/pub/data/cdo/documentation/GHCND_documentation.pdf
% the units for wind speed are m/s because this 
% download was set for metric units

% read in the data
data = readtable("WindData_AnnArborAirport_2013-2022.csv");

% pull out the wind speed
windSpeed = data.AWND;

% remove the NaNs
windSpeed(isnan(windSpeed)) = [];

% calculate mean and standard deviation
windSpeedMean = mean(windSpeed);
windSpeedStdDev = std(windSpeed);

disp('The mean value of the wind speeds is:');
disp(windSpeedMean);

disp('The standard deviation of the wind speeds is:');
disp(windSpeedStdDev);


% calculate the probability that the wind speed will be 
% between 1.25 m/s and 5.50 m/s. 

probLowerBound = normcdf(1.25,windSpeedMean,windSpeedStdDev);
probUpperBound = normcdf(5.50,windSpeedMean,windSpeedStdDev);
probInRange = probUpperBound - probLowerBound;

disp('The probability that the wind speed will be between 1.25 m/s and 5.5 m/s is:')
disp(probInRange);



% calculate how often the wind speed will be higher than a given value

windGustSpeed = 10.0; % 10 m/s = about 22 mph
numOccurrences = 1 ./ (1 - normcdf(windGustSpeed,windSpeedMean,windSpeedStdDev));

% the original measured data has an implicit time interval to it: the data
% was recorded once per day. So each "occurrence" corresponds to one day.

disp(['The wind speed will be around ' num2str(windGustSpeed) ' m/s about']);
disp(['once every ' num2str(numOccurrences) ' days']);
disp(['(' num2str(numOccurrences/365) ' years).' ]);
