%
% sample matlab script for reading high resolution wind fields
% BYU June 2002
% revised for new latitude scaling
%
% how to extract the data from the file is show, though not all
% available data is extracted or lotted.

%fid1=fopen('/auto/cers1/users/long/src/linux/Seawinds/hdf/c/L1B/QS_S1B01201.20001321038.avewr.WRave2','r','ieee-be');
fid1=fopen('out.WRave2','r','ieee-be');

% plot retrieved winds

dumb=fread(fid1,1,'int');
a=fread(fid1,[5],'int')';
ngrid=a(1);
wctgrid=a(2);
watgrid=a(3);
ctgrid=a(4);
atgrid=a(5);
dumb=fread(fid1,1,'int');

dumb=fread(fid1,1,'int');
a=fread(fid1,[4],'int')';
natstart=a(1);
natend=a(2);
nctstart=a(3);
nctend=a(4);
dumb=fread(fid1,1,'int');

dumb=fread(fid1,1,'int');
a=fread(fid1,[4],'int')'
iastart=a(1);
iaend=a(2);
icstart=a(3);
icend=a(4);
dumb=fread(fid1,1,'int');

dumb=fread(fid1,1,'int');
ascnode=fread(fid1,1,'float');
a=fread(fid1,[3],'int')'
irev=a(1);
year=a(2);
day=a(3);
dumb=fread(fid1,1,'int');
if day < 0
  disp('Rain-containing file');
end

fa=iastart-iastart+1;
la=iaend-iastart+1;
fc=icstart-icstart+1;
lc=icend-icstart+1;
wid=min(lc-fc+1,la);
ind=[1:lc];

nchoice=zeros([icend-icstart+1,iaend-iastart+1]);
speed=zeros([icend-icstart+1,iaend-iastart+1]);
dir=zeros([icend-icstart+1,iaend-iastart+1]);
lat=zeros([icend-icstart+1,iaend-iastart+1]);
lon=zeros([icend-icstart+1,iaend-iastart+1]);
rain=zeros([icend-icstart+1,iaend-iastart+1]);
ic=zeros([icend-icstart+1,iaend-iastart+1]);
ia=zeros([1,iaend-iastart+1]);
wid=700;


for i=iastart:iaend
  dumb=fread(fid1,1,'int');
  b=fread(fid1,9*lc,'int16');
  n=b(9*(ind-1)+1);
  nchoice(ind,i-iastart+1)=n;
  n(find(n<=0 | n>4))=1;
  ind2=2*(n-1)+1+9*(ind'-1)+1;
  speed(ind,i-iastart+1)=b(ind2)/200.0; 
  dirt=b(ind2+1)/180;
  dirt(find(dirt<0))=dirt(find(dirt<0))+360.0;
  dir(ind,i-iastart+1)=dirt;
  
  a=fread(fid1,[2*lc],'int16');
  land=a(2*ind-1);
  land(find(land>0))=1;
  land(find(land<0))=0;
  lat(ind,i-iastart+1)=(abs(a(2*ind-1))-16201)/180.0;
  lon(ind,i-iastart+1)=a(2*ind)/180.0;

  c=fread(fid1,4*lc,'int16');
  ind3=n+4*(ind'-1)+1;
  rain(ind,i-iastart+1)=c(ind3)/100.0;

  ia(i-iastart+1)=fread(fid1,1,'int16');
  ct=fread(fid1,lc,'int16');
  ic(ind,i-iastart+1)=ct;
  
% apply land mask
  speed(ind,i-iastart+1)=speed(ind,i-iastart+1).*land+(1-land)*200;
    
  dumb=fread(fid1,1,'int');

end
fclose(fid1);


if (1==0)

myfigure(1)
image(nchoice*15);
txt = sprintf('title(''Choice*15   Alongtrack [%d to %d]  Crosstrack [%d to %d]'')',(i+iastart-1)/ngrid, (i+wid-1+iastart-1)/ngrid,(fc+icstart-1)/ngrid,(lc+icstart-1)/ngrid);
eval(txt)
colormap('gray')
colorbar;

end

myfigure(2)
colormap('gray');
%image(dir*64/360);
imagesc(dir,[0 360]);

cm=gray;
cm2=cm;
[ncm k]=size(cm);
for k=1:ncm/2
  cm2(k,:)=1-cm(2*k-1,:);
  cm2(ncm-k+1,:)=1-cm(2*k-1,:);
end
colormap(cm2);
colorbar;
txt = sprintf('title(''Direction (deg)   Alongtrack [%d to %d]  Crosstrack [%d to %d]'')',(iastart-1)/ngrid, (wid-1+iastart-1)/ngrid,(fc+icstart-1)/ngrid,(lc+icstart-1)/ngrid);
eval(txt)
%axis([0 700 0 600])
axis off


myfigure(3)
imagesc(speed,[0 50]);
%colormap('gray');
colormap('jet');
colorbar;
txt = sprintf('title(''Speed (m/s)   Alongtrack [%d to %d]  Crosstrack [%d to %d]'')',(iastart-1)/ngrid, (wid-1+iastart-1)/ngrid,(fc+icstart-1)/ngrid,(lc+icstart-1)/ngrid);
eval(txt)
%axis([0 700 0 600])
axis off

myfigure(4)
[m n]=size(speed)
[N X]=hist(reshape(speed,[1 m*n]),0.5+[0:49]);
plot(X,log(N))
xlabel('Wind speed (m/s)');
ylabel('log count')
title('Wind speed histogram');

myfigure(5)
imagesc(rain,[0 100]);
colormap('jet');
colorbar;
txt = sprintf('title(''rain (km-mm/hr)   Alongtrack [%d to %d]  Crosstrack [%d to %d]'')',(iastart-1)/ngrid, (wid-1+iastart-1)/ngrid,(fc+icstart-1)/ngrid,(lc+icstart-1)/ngrid);
eval(txt)
axis off
















