% C:\Data\WSRA\WSRA_L4_ATOMIC.m      24feb2020

%path(path,'c:\data\WSRA')
path(path,'E:DATA\ATOMIC\P3\wsra')

filename = 'WSRA-L4-20200210I.nc';  

ncid = netcdf.open(filename,'NOWRITE');
directional_wave_spectrum = netcdf.getVar(ncid,0);
directional_wave_spectrum_180 = netcdf.getVar(ncid,1);
dominant_to_secondary_partition_angle = netcdf.getVar(ncid,2);
dominant_wave_direction = netcdf.getVar(ncid,3);
dominant_wave_height = netcdf.getVar(ncid,4);
dominant_wave_length = netcdf.getVar(ncid,5);
hurricane_eye_distance_east = netcdf.getVar(ncid,6);
hurricane_eye_distance_north = netcdf.getVar(ncid,7);
latitude = netcdf.getVar(ncid,8);
longitude = netcdf.getVar(ncid,9);
peak_spectral_variance = netcdf.getVar(ncid,10);
platform_course = netcdf.getVar(ncid,11);
platform_orientation = netcdf.getVar(ncid,12);
wsra_computed_roll = netcdf.getVar(ncid,13);
platform_radar_altitude = netcdf.getVar(ncid,14);
platform_speed_wrt_ground = netcdf.getVar(ncid,15);
rainfall_rate = netcdf.getVar(ncid,16);
rainfall_rate_median = netcdf.getVar(ncid,17);
sea_surface_mean_square_slope = netcdf.getVar(ncid,18);
sea_surface_mean_square_slope_median = netcdf.getVar(ncid,19);
sea_surface_wave_significant_height = netcdf.getVar(ncid,20);
secondary_wave_direction = netcdf.getVar(ncid,21);
secondary_wave_height = netcdf.getVar(ncid,22);
secondary_wave_length = netcdf.getVar(ncid,23);
timeWSRA = netcdf.getVar(ncid,24);
trajectory = netcdf.getVar(ncid,25);
wavenumber_east = netcdf.getVar(ncid,26);
wavenumber_north = netcdf.getVar(ncid,27);
wind_direction = netcdf.getVar(ncid,28);
wind_speed = netcdf.getVar(ncid,29);
delete_half_plane_angles = netcdf.getVar(ncid,30);
roll=ncread(filename,'wsra_computed_roll');
mss = sea_surface_mean_square_slope;



netcdf.close(ncid);

dtr = pi/180.;
dhpa = dtr*delete_half_plane_angles; 
dhpk = 0.0025*[28 25 22 19 16 13 10 7];
wh = sea_surface_wave_significant_height;
dwd=dominant_wave_direction;
dwh=dominant_wave_height;
dwl=dominant_wave_length;
swd=secondary_wave_direction;
swh=secondary_wave_height;
swl=secondary_wave_length;
nkm = hurricane_eye_distance_north/1000.;
ekm = hurricane_eye_distance_east/1000.;
gs = platform_speed_wrt_ground;
lt = latitude;
ln = longitude;
alt = platform_radar_altitude;
thacircl = (0:360)*dtr;
xcircl = sin(thacircl);
ycircl = cos(thacircl);
v = (0.1:0.1:0.9);
X = (-0.07854:0.00245437:0.07854);
Y = X;
dws = directional_wave_spectrum;
dws180 = directional_wave_spectrum_180;
mss = sea_surface_mean_square_slope;
croll = wsra_computed_roll;
rr = rainfall_rate;
wd = wind_direction + 180.;
ws = wind_speed;
trk = platform_course;
hdg = platform_orientation;
[n m]=size(roll);rollm=[];rollx=[];mssm=[];rrm=[];
for i=1:m;
    rollm(i)=nanmedianx(roll(:,i));
    rollx(i)=max(abs(roll(:,i)));
    mssm(i)=nanmedianx(mss(:,i));
    rrm(i)= nanmedianx(rr(:,i));
end;
indey=1:m;

psv = peak_spectral_variance;
mx = max(size(trajectory));
flagr=ones(length(dwh),1);
for i=1:length(dwh)
    ij=find(abs(croll(:,i))>3);
    if(~isempty(ij));flagr(i)=0;end;
end;

clf
hold on
%axis([-153.5 -153.5+1.26*5/cos(14*dtr) 12 17])
for i=10:10:mx
plot(ln(i),lt(i),'bo','Linewidth',2.5)
end
for i=20:20:mx
plot(ln(i),lt(i),'ro','Linewidth',2.5)
end
plot(ln,lt,'k')
plot(ln,lt,'k.','markersize',10)
title(filename)
%plot([0 0],[-200 400],'k')
%plot([-300 400],[0 0],'k')
xlabel('LONGITUDE                   ')
ylabel('LATITUDE')
% print


pause

clf
plot(ws,alt,'k.','markersize',10)
title(filename)
xlabel('WIND SPEED (m/s)')
ylabel('AIRCRAFT ALTITUDE (m)')

pause

clf
hold on
for k=1:mx
  for j=1:5
    if(abs(croll(j,k))<4.) 
       plot(k,mss(3,k),'ko')
    end
  end
end
title([filename])
xlabel('SPECTRUM NUMBER')
ylabel('mean square slope (mss)')
axis([0 mx 0 0.08])

pause

clf
hold on
for k=1:mx
       plot(k,wh(k),'ko')
end
title([filename])
xlabel('SPECTRUM NUMBER')
ylabel('WAVE HEIGHT (m)')
axis([0 mx 0 5])

pause

for k=1:mx
clf
hold on
%scl = 0.08;
%scl = 0.0628;
%scl = 0.042;
scl = 0.05;
axis(scl*[-1 1.52 -1 1])   
xlabel('               east wavenumber (rad/m)')
ylabel('                                      north wavenumber (rad/m)')
       latfx =  fix(lt(k)*100)*0.01;
       lonfx = -fix(ln(k)*100)*0.01;
       nkmfx = fix(nkm(k)*10)*0.1;
       ekmfx = fix(ekm(k)*10)*0.1;
       trk = dtr*platform_course(k);
       hgt = fix(alt(k)+0.5);
for i=0:35
    tha = i*10*dtr; 
plot([0 scl*sin(tha)],[0 scl*cos(tha)],'g')    
end
for i=0:3:35
    tha = i*10*dtr; 
plot([0 scl*sin(tha)],[0 scl*cos(tha)],'b')    
end
for i=100:100:300
    wn = 2*pi/i;
    plot(wn*xcircl,wn*ycircl,'r')
end
for i=150:100:350
    wn = 2*pi/i;
    plot(wn*xcircl,wn*ycircl,'--r')
end
xws = -wind_speed(k)*sin(wind_direction(k)*dtr)/1000.;
yws = -wind_speed(k)*cos(wind_direction(k)*dtr)/1000.;
plot([0 xws],[0 yws],'b','Linewidth',4)
plot([0 0.01*sin(trk)],[0 0.01*cos(trk)],'r','Linewidth',3)

% directional wave spectrum
contour(X,Y,dws(:,:,k)/psv(k),v,'k')

% directional wave spectrum with artifact lobes retained
% contour(X,Y,dws180(:,:,k)/psv(k),v,'k')
title([filename,'   spec # ',num2str(k),',  ',num2str(latfx),'N, ',num2str(lonfx),'W                                         '])

text(-0.95*scl,0.9*scl,[num2str(fix(wh(k)*10)*0.1),' m SWH'])
text(-0.95*scl,0.8*scl,[num2str(fix(alt(k)+0.5)),' m hgt'])
text( 0.47*scl,0.9*scl,[num2str(psv(k)),' psv'])

pause

end

