%% Computatiton of dewpoint Temperature
tempVals=cell2mat(struct2cell(load('Temp.mat')))
humiVals=cell2mat(struct2cell(load('Humi.mat')))

%% Means
avgTemp=mean(tempVals);
avgHumi=mean(humiVals);


%%
alpha=6.1121;
beta=18.678;
lambda=257.14;
theta=234.5;

P_s=alpha*exp(beta-avgTemp/theta)*(avgTemp/(lambda+avgTemp)); % saturation water vapour pressure
Gamma=log(avgHumi/100*P_s/alpha);

T_dpC=lambda*Gamma/(beta-Gamma); % dewpoint temperature

avgTempF=9/5*(avgTemp)+32 % temp in farenheit (do not use in computation, for visualization only)
T_dpF= 9/5*(T_dpC)+32    % dewpoint temp in farenheit (do not use in computation, for visualization only)

% [avgTemp,T_dpC] % for side-by-side visualization
% [avgTempF,T_dpF]% for side-by-side visualization

%% Uncertainty Computation
noOfMeasurements=25;

stdTemp=std(tempVals); % standard deviation of temp data
stdHumi=std(humiVals); % standard deviation of humidity data

tScore=1.711; % t-score for 95% probability and dof=24

randUncenTemp=tScore*stdTemp/sqrt(noOfMeasurements); % standard uncertainty in temp
randUncenHumi=tScore*stdHumi/sqrt(noOfMeasurements); % standard uncertainty in humidity

resUncenTemp=1/2*1; % uncertainty in temp due to resolution
resUncenHumi=1/2*.01; % uncertainty in humidity due to resolution

tempUncen= sqrt(randUncenTemp^2+resUncenTemp^2); % uncertainty in temperature
humiUncen=sqrt(randUncenHumi^2+resUncenHumi^2); % uncertainty in humidity

uncenAlpha=.01*alpha;   % uncertainty in alpha
uncenBeta=.01*beta;     % uncertainty in beta
uncenLambda=.01*lambda; % uncertainty in lambda
uncenTheta=.01*theta;   % uncertainty in theta

syms a b l t T H

P_ssym=a*exp((b-T/t)*(T/(l+T)));   % saturation water vapour pressure as a symbolic function using a=alpha, b=beta, l=lambda, t=theta, T=temp, H=humidity as symbolic variables
gammasym=log(H/100*P_ssym/a);      % gamma as a symbolic function using a=alpha, b=beta, l=lambda, t=theta, T=temp, H=humidity as symbolic variables
T_dpsym=(l*gammasym/(b-gammasym)); % dewpoint temperature as a symbolic function using a=alpha, b=beta, l=lambda, t=theta, T=temp, H=humidity as symbolic variables

dTdpdT=double(simplify(subs(diff(T_dpsym,T),[a b l t T H], [alpha beta lambda theta avgTemp avgHumi])));      % propagation coefficient d(T_dp)/dT evaluated at averages
dTdpdH=double(simplify(subs(diff(T_dpsym,H),[a b l t T H], [alpha beta lambda theta avgTemp avgHumi])));      % propagation coefficient d(T_dp)/dH evaluated at averages
dTdpdalpha=double(simplify(subs(diff(T_dpsym,a),[a b l t T H], [alpha beta lambda theta avgTemp avgHumi])));  % propagation coefficient d(T_dp)/d(alpha) evaluated at averages
dTdpdbeta=double(simplify(subs(diff(T_dpsym,b),[a b l t T H], [alpha beta lambda theta avgTemp avgHumi])));   % propagation coefficient d(T_dp)/d(beta) evaluated at averages
dTdpdlambda=double(simplify(subs(diff(T_dpsym,l),[a b l t T H], [alpha beta lambda theta avgTemp avgHumi]))); % propagation coefficient d(T_dp)/d(lambda) evaluated at averages
dTdpdtheta=double(simplify(subs(diff(T_dpsym,t),[a b l t T H], [alpha beta lambda theta avgTemp avgHumi])));  % propagation coefficient d(T_dp)/d(theta) evaluated at averages

uncenDewpoint= sqrt((dTdpdT*tempUncen)^2+(dTdpdH*humiUncen)^2+(dTdpdalpha*uncenAlpha)^2+(dTdpdbeta*uncenBeta)^2+(dTdpdlambda*uncenLambda)^2+(dTdpdtheta*uncenTheta)^2);% propagate uncertainties to dewpoint temp using propagation formula

dewpointLowerBound=double(simplify(subs(T_dpsym,[a b l t T H],[alpha beta lambda theta avgTemp avgHumi])))-uncenDewpoint % compute upperbound in dewpoint temp from uncertainty
dewpointUpperBound= double(simplify(subs(T_dpsym,[a b l t T H],[alpha beta lambda theta avgTemp avgHumi])))+uncenDewpoint % compute lowerbound in dewpoint temp from uncertainty