MRdata=read_hwd_data_trd(2013,1,27,2013,2,26);

% plot the time series of zonal wind at 88 km altitude (height bin 3)
% during the data period

x1 = MRdata.time(3,:,1);
x2 = MRdata.time(3,:,2)+ 24*ones(size(x1));
x3 = MRdata.time(3,:,3)+ 48*ones(size(x2));
x4 = MRdata.time(3,:,4)+ 72*ones(size(x3));

y1 = MRdata.u(3,:,1);
y2 = MRdata.u(3,:,2);
y3 = MRdata.u(3,:,3);
y4 = MRdata.u(3,:,4);

data = [x1 x2 x3 x4; y1 y2 y3 y4]';

datax = data(:,1);
datay = data(:,2);

Datafitting = testfit(datax,datay);
ci = confint(Datafitting);