function fig8_20

%  DATE  TMAX  TMIN  TOBS
%  temperature values are in Celsius degrees to tenths
%  temp=-9999 means not recorded (need to skip these)
%  data source: http://www.ncdc.noaa.gov
%  date(1)=1 is 2014/01/01  and  date(end)=2015/12/16

load tempss.txt
AA=sortrows(tempss);
[n1 n2 n3 n4]=size(AA);

ic=0;
for i=1:n1
    if AA(i,2) ~= -9999
        year=floor(AA(i,1)/10000);
        day= round((AA(i,1)/100 - floor(AA(i,1)/100))*100);
        m=floor(AA(i,1)/100);
        month= round((m/100 - floor(m/100))*100);
        ic=ic+1;
        datez=datenum(year,month,day);
        if ic==1
            date0=datez;
        end
        date(ic)=datez-date0+1;
        temp(ic)=AA(i,2)*0.1;
    end
end
nd=length(date)

%%% spline fit
m=5;
ptsm=linspace(min(date),max(date),m);
dp=ptsm(2)-ptsm(1);
pts=[ ptsm(1)-dp ptsm ptsm(m)+dp];
B=zeros(m+2,nd);
for id=1:nd
    for im=1:m+2
        x=(date(id)-pts(im))/dp;
        B(im,id)=bbspline(x);
    end
end
C=B*B';
cond(C,inf)
d=B*temp';
a=C\d;
% evaluate interpolation function s(x)
nx=200;
dat=linspace(min(date),max(date),nx);
for ix=1:nx
    sum=0;
    for k=1:m+2
        xx=(dat(ix)-pts(k))/dp;
        sum=sum+a(k)*bbspline(xx);
    end
    s(ix)=sum;
end

clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
subaxis(1,1,1,1,'MT',0.01,'MB',0.16,'MR',-0.02,'ML',0.05,'P',0.04)
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
hold on

%%% plot entire data set
%plot(date,temp,'or','MarkerSize',8,'LineWidth',2)

%%% plot nx (randomly choosen) data points
nx=50;
r=randi(nd,1,nx);
for i=1:nx
    d(i)=date(r(i));
    t(i)=temp(r(i));
end

plot(d,t,'or','MarkerSize',8,'LineWidth',2)
plot(dat,s,'b','LineWidth',2)
grid on
box on
xlabel('Date')
ylabel('Temperature')
tmax=datenum(2015,12,31)-date0+1
t2=datenum(2014,7,1)-date0+1;
t3=datenum(2015,1,1)-date0+1;
t4=datenum(2015,7,1)-date0+1;
axis([0 tmax -20 40])
set(gca,'xtick',[1 t2 t3 t4 tmax])
set(gca,'XTickLabel',{'Jan';'July';'Jan';'July';'Jan'})
set(gca,'FontSize',14,'FontWeight','bold')


%%% generate data for exercise in text
nxx=15;
dd(1)=date(1);
tt(1)=temp(1);
dd(nxx)=date(nd);
tt(nxx)=temp(nd);
nd
nn=round(nd/(nxx-1))
i1=1;
for i=nn:nn:nd
    i1=i1+1;
    iz=i+randi(16)-8;
    dd(i1)=date(iz);
    tt(i1)=temp(iz);
end
dd
tt

%exportgraphics(gcf,'/Users/mark/Desktop/temp.eps')

% figure(2)
% clf
% % get(gcf)
% set(gcf,'Position', [23 813 658 230])
% hold on
% axis([0 tmax -20 40])
% plot(dd,tt,'ro','LineWidth',2)


function y=bbspline(x)
% Calculate the value of a cubic B-spline at point x
x=abs(x) ;
if x>2,
    y=0 ;
else
    if x>1,
        y=(2-x)^3/6 ;
    else
        y=2/3-x^2*(1-x/2) ;
    end
end







