% SRAsfmr9sep04.m
yeary=2004;yeardol=num2str(yeary);
cfdir=['c:\data\cwf\dropbox\data\cblast\SRA\' yeardol '\'];
jx=dir([cfdir 'hdr*.asc']);
[k,kx]=size(jx);
for i=1:k
    disp([num2str(i) '  ' jx(i).name])
end;
j=1;
j=input('File number?  ')
fldol=jx(j).name;
hdr=load([cfdir fldol]);

h   = hdr(:,1);               %aircraft height (m)
Ns  = hdr(:,2);               %starting spectrum number ofm5-specm average
lt  = hdr(:,3)/100.;          %latitude (deg)
ln  = hdr(:,4)/100.;          %longitude (deg)
wht = hdr(:,5)/100.;          %SWH (total wave ht, m)
Spk = hdr(:,6)/8.64e4;        %peak spectral density (m^4/rad^2/
gsn = hdr(:,7)/10.;           %N comp of ground speed (m/s)
gse = hdr(:,8)/10.;           %E comp of ground speed (m/s)
hed = hdr(:,9)/10.;           %aircraft heading (deg)
un  = hdr(:,10)/10.;          %N comp of wind speed (m/s) flt level
ue  = hdr(:,11)/10.;          %E comp of wind speed (m/s) flt level
sod = hdr(:,12);              %UTC seconds of day
Nspa = hdr(:,13);             %Actual number spec our of 5 used in the average
nkm0 = hdr(:,14)/10.;         %location in km N of eye   
ekm0 = hdr(:,15)/10.;         %location in km E of eye
wl1 = hdr(:,16)/10.;          %Wavelentth (m) dominant spec peak
wd1 = hdr(:,17)/10.;          %Propagation direction (Deg) dominant spec peak
wh1 = hdr(:,18)/100.;         %Significant wave height (m) dominant spec peak
wl2 = hdr(:,19)/10.;          %Wavelentth (m) secondary spec peak
wd2 = hdr(:,20)/10.;          %Propagation direction (Deg) secondary spec peak
wh2 = hdr(:,21)/100.;         %Significant wave height (m) secondary spec peak
ii= wh2==0;wh2(ii)=NaN;
ws1 = wh1./wl1;%wave slope
wst = wht./wl1;
hod = sod/3600.;              
cp1=sqrt(9.83/2/pi*wl1);%Phase speed (m/s), dominant wave peak
cp2=sqrt(9.83/2/pi*wl2);ii=find(cp2<1);cp2(ii)=NaN;%Phase speed (m/s), secondary wave peak

% %%%%   rotate to storm heading
 dayy=fldol(4:5);
 dateyy=fldol(4:10)
% head=-65;
% if dayy=='12';head=-75;end;
% if dayy=='14';head=-30;end;
    

nkm=nkm0;%nkm0*cos(head*pi/180)+ekm0*sin(head*pi/180);
ekm=ekm0;%-nkm0*sin(head*pi/180)+ekm0*cos(head*pi/180);
%%%%%

%plot(ln,lt,'.')
figure;
plot(ekm,nkm,'.')
axis([-250 250 -200 200])
xlabel('km east of eye')
ylabel('km north of eye')
title('SRA storm-relative wave spectrum locations on 9sep04')

figure;plot(wh1,ws1,'k.')
xlabel('SRA dominant wave height (m)')
ylabel('SRA wave steepness')
title(['Wave steepness dominant wave ' dateyy])
axis([0 14 0 0.07])


plot(wht,wst,'k.')
xlabel('SRA total wave height (m)')
ylabel('SRA wave steepness')
title(['Wave steepness total wave height  '  dateyy])
axis([0 14 0 0.07])

rdist=sqrt(ekm.^2+nkm.^2);speed=sqrt(un.^2+ue.^2);
figure;plot(rdist,wh1,'o',rdist,wh2,'.');xlabel('Radial distance (km)');ylabel('SRA total wave height (m)');title(['Wave height  '  dateyy]);
legend('Peak1','Peak2');
figure;plot(rdist,cp1,'o',rdist,cp2,'.');xlabel('Radial distance (km)');ylabel('SRA wave phase speed (m/s)');title(['Wave Phase speed  '  dateyy]);



figure;scatter(ekm,nkm,20,speed,'linewidth',3);colorbar;title(['Flt Level Wind speed (m/s)  ' dateyy]);
figure;scatter(ekm,nkm,20,cp1,'linewidth',3);colorbar;title(['Wave speed (m/s)  '   dateyy]);
figure;scatter(ekm,nkm,20,wht,'linewidth',3);colorbar;title(['Wave height (m)  '   dateyy]);
figure;scatter(ln,lt,20,wht,'linewidth',3);colorbar;title(['Wave height (m)  '   dateyy]);



