clear all
close all

% chapter 3 

% figure 3.1
normspec([.5,1],0,1)

% figure 3.2
figure(2)
x=-3:.01:3;
f=cdf('norm',x,0,1);
plot(x,f)

%read arrays of yearly lake level
% data has been preprocessed for 1895-2010 period

%lake level data with year,number of values, lake level (ft)
load gsl_yr.txt;
time = gsl_yr(:,1);
%convert to meters
lev = gsl_yr(:,3)*.3048;

%use for figures 3.3 and 3.4
%dfittool(lev)

% load collins air temperature data (as well as other stations)
  [year,mo,dy,hr,atb,agd,cln,amb,alt] = textread('alta_stns.dat','%f %f %f %f %f %f %f %f %f');

%use for figure 3.5 
%comment out first call to dfittool above
%dfittool(cln)

% for figure 3.6
p = normspec([-1,1],0,1);

p = normspec([-2,2],0,1);

p= normspec ([2,Inf],0,1);

xq = norminv([0.05 0.95],0,1)

figure(6)
x=0:300; 
y=geocdf(x,0.01); 
stairs(x,y);

figure(7)
x=0:300; 
y=geocdf(x,0.02); 
stairs(x,y);

%figure 3.9
load klamath.txt
figure(8)
yr = klamath(:,1);
flow = klamath(:,2);
bar(yr,flow)
%dfittool(flow)

%figure 3.10
%dfittool(cln)
mu=-5.1
sig=5.9
normspec([mu-1.96*sig, mu+1.96*sig], mu, sig);

%figure 3.11

[nd,nn] = size(year);
for i=1:nd
    if(year(i)<50)
        year(i)=year(i)+2000;
    else
        year(i)=year(i)+1900;
    end
end
dates = datenum(year,mo,dy,hr,0,0);
figure(10)
hold on
x1 = [dates(1), dates(nd)];
y1 = [ 6, 6];
y2 = [-16,-16];
plot(x1,y1,'r');
plot(x1,y2,'g')
datetick('x',25)
bar(dates,cln)

% annual total utah precip (cm)
load  ut_ppt_yr.txt;
ppt = ut_ppt_yr(:,2);
meanp = mean(ppt);
sdp = std(ppt,1);
anomp = ppt - meanp;

p = normspec([-1.96*sdp,1.96*sdp],0,sdp);

sd5 = sdp/sqrt(5);
p = normspec([-1.96*sd5,1.96*sd5],0,sd5);

% now do for all 5 year periods in the record
h05 = zeros(112,1);
p05 = zeros(112,1);

for iyr=1897:2008
    ii = iyr-1896;
    valy = anomp(ii:ii+4);
[h,p,ci,stat]= ttest(valy,0,.05,'left');
h05(ii) = h;
% store the values so that the cases where the null hypothesis can be
% rejected appear near 1
p05(ii) = 1- p;
end
figure(13)
x = 1897:1:2008;
bar(x,h05,'r');
figure(14)
bar(x,p05,'b');

% a simple little comparison of how the t test values
% are distributed for a sample of 5 vs. a samplof 30, 
% which is almost the same as a gaussian
% i.e. if your sample size is greater than 30
% the t test values are distributed like a gaussian
figure(15);
x = -5:.1:5;
y = tpdf(x,5);
z = normpdf(x,0,1);
a = tpdf(x,30);
plot(x,y,'r',x,a,'g',x,z,'b');



