

% Read the IEEE_big_endian binary data, 6-hourly (4 times per day/file)
% The FORTRAN write Statment were:  
%		do it=1,itm
%			write(filename_unit) u
%			write(filename_unit) v
%                enddo ! it; itm=4 per day/file
%   where the arrays in FORTRAN are: u(im,jm), v(im,jm), 
% w/ im=1440 (for lon index),jm=719 (for lat index),
% w/ lon=0:.25:359.75; lat=-89.75:.25:89.75;

clear all	% house clearning for new data
close all

lonp = 360-124.3; %deg E
latp = 46.9; %deg N
tstart = datenum(2001,5,1);
tend = datenum(2001,8,6);
flagplot = true(1);
lim = [lonp-5,lonp+5,latp-5,latp+5];

%******** DO NOT CHANGE **************************************
undef=-9999.;   % bad vaule, will be assigned as NaN in MATLAB
xdim=1440;
ydim=719;
tdim=4;
mapsiz=xdim*ydim;
lon=0:0.25:359.75;
lat=-89.75:.25:89.75;
%******** DO NOT CHANGE **************************************

t = tstart:6/24:tend;
nt = length(t);

%Initialize data
up = zeros(nt,1);
vp = up;
w = zeros(ydim,xdim);
u = w; v = w;

if(flagplot)
    figure    
    ix = find(lim(1)<lon & lim(2)>lon);
    iy = find(lim(3)<lat & lim(4)>lat);
    h = imagesc(lon(ix),lat(iy),w(iy,ix));
    set(gca,'ydir','normal')
    hold on
    h3 = quiver(lon(ix),lat(iy),u(iy,ix),v(iy,ix),'k');
    axis equal
    set(gca,'xlim',lim(1:2),'ylim',lim(3:4))    
    caxis([0,20])
    colorbar    
    h2 = plot(lonp,latp,'sk','markersize',4);
    xlabel('Longitude, ºE')
    ylabel('Latitude, ºN')
    title('Surface Winds, m/s')
end

% return
ii=0;
ndays = ceil(nt/4);
for j=1:ndays
    str = datestr(tstart+j-1,'yyyymmdd');
    fname=['..\',str(1:4),'\uv',str,'.be'];
    disp(['Reading: ',fname])
    fid=fopen(fname,'r','b');	% 'r' = read only; 'b'=IEEE_big_endian
    if(fid==-1)
        error(['ERROR: Could not find file: ',fname])
    end
    fseek(fid,0,'bof');			% position at statrt of data file
    for it=1:tdim				% loop over the time records
        % now u first:
        fseek(fid,4,'cof');			% skip the 4-byte header in FORTRAN Binary write
        u=fread(fid,mapsiz,'float32');		% float32 = real*4 = single precision
        u(u==undef)=nan;			% assign bad values to NaN
        u=reshape(u,[xdim ydim])';		% re-shape 1-D u to 2-D, AND transpose x-y indexing: 1st for lat; 2nd for lon
        fseek(fid,4,'cof');			% skip 4-byte eof in FORTRAM binary write
        
        % now for v:
        fseek(fid,4,'cof');			% skip 4-byte header in FORTRAM binary write
        v=fread(fid,mapsiz,'float32');
        v(v==undef)=nan;
        v=reshape(v,[xdim ydim])';
        fseek(fid,4,'cof');
        
        %Plotting
        if(flagplot)
            w=(u.^2+v.^2).^0.5;	% speed/magnitude
            set(h,'cdata',w(iy,ix))
            set(h3,'Udata',u(iy,ix),'Vdata',v(iy,ix))
            drawnow
%             pause 
        end
        
        %Interpolation
        ii=ii+1;
        if(ii>nt)
            break
        end
        up(ii) = interp2(lon,lat,u,lonp,latp,'linear');
        vp(ii) = interp2(lon,lat,v,lonp,latp,'linear');
    end
    fclose(fid);
end

%Interpolate gaps
ind = isnan(up);
up(ind) = interp1(t(~ind),up(~ind),t(ind));
ind = isnan(vp);
vp(ind) = interp1(t(~ind),vp(~ind),t(ind));

%Speed and direction
wp = sqrt(up.^2+vp.^2);
dirp = mod(atan2(vp,up)*180/pi+180,360);

%Plotting
figure
plot(t,up,'b.-',t,vp,'r.-')
datetick('x',2)
ylabel('Wind Speed, m/s')
legend('U','V',0)

figure
subplot(2,1,1)
plot(t,wp,'b.-')
ylabel('Wind Speed, m/s')
datetick('x',2)
subplot(2,1,2)
plot(t,dirp,'r.-')
ylabel('Wind Direction, deg')
datetick('x',2)

%Write SMS files
tp = (t-tstart)'*24;
fid = fopen('Wind_Speed.xys','w');
fprintf(fid,'XYS 1 %d "%s"\r',nt,datestr(tstart,31));
fprintf(fid,'%f %f\r',[tp,wp]');
fclose(fid);
fid = fopen('Wind_Direction.xys','w');
fprintf(fid,'XYS 2 %d "%s"\r',nt,datestr(tstart,31));
fprintf(fid,'%f %f\r',[tp,dirp]');
fclose(fid);

return


