clear all
fid = fopen('imet_14838.dat','r');

header = fgetl(fid);
count_line = 0
data = [];

while 1
count_line = count_line + 1;
line = fgetl(fid);
if ~isstr(line), break, end

line_cor = line;
line_cor = strrep(line_cor,'/',' ');
line_cor = strrep(line_cor,':',' ');

line_data = (sscanf(line_cor,'%f'))';

if length(line_data) == 22;
  data = [data;line_data];
else
  disp(['bad character in line ' num2str(count_line), ':']);
  disp(line_cor);
  disp('This part has been read successfully:');
  disp(num2str(line_data));
  return;
end

end
fclose(fid)

mon = data(:,1);
day = data(:,2);
year = data(:,3) + 2000;
hour = data(:,4);
mnt = data(:,5);
sec = data(:,6);

wind_sp = data(:,16); 
wind_dir = (data(:,17)/360)*2*pi;
wind_x = sin(wind_dir).*wind_sp; 
wind_y = cos(wind_dir).*wind_sp;

time = datenum(year,mon,day,hour,mnt,sec) - datenum(2004,1,1,0,0,0)

% calculate speed from difference in x,y abd delta t
lat = data(:,7);
lon = data(:,8);

[x,y] = ll2xy(lat,lon);
delta_xy = sqrt(diff(x).^2 + diff(y).^2);
delta_t = diff(time);
speed_time = time(1:(end-1)) + diff(time)/2;

% filter out some erratic points
aa = find((delta_xy < 800) & (delta_t ~= 0));
delta_xy_f = delta_xy(aa);
delta_t_f = delta_t(aa);
speed_time_f = speed_time(aa);

% speed in knots
speed = (delta_xy_f/1853.25)./(delta_t_f*24);

figure(1)
plot(speed_time_f,speed);
title('estimated ship velocity');
xlabel('yearday');
ylabel('knots')

figure(2)
y = zeros(size(time,1),1);
quiver(time,y,wind_x,wind_y,0);



