%READSCAT.M
%Driver to generate figure for scatterometer winds near hurricane
%INPUT: fname - filename of scatterometer data
%       storm - storm name
%       cntr - [latitude, longitude] of storm center, longitude from 0-360
%       pass - 1 for ascending pass, 2 for descending pass (for plotting)
%       pad - degrees from storm center to extract data
%This file is configured to read either quickscat or windsat data
%from data.remss.com. The routine for decoding binary files distinguishes
%between quickscat and windsat data based on the file name: quickscat files
%should begin with 'qscat', and windsat files should begin with 'wsat'.

clear all; close all; clc

%USER DEFINED INPUT VARIABLES
%------------------------------------------
fname = 'qscat_20080711v4'; % filename of scatterometer Data
storm = 'Bertha'; % storm name 
cntr = [29.0 297.7]; % [latitude, longitude] of storm center, longitude from 0-360
pass = 2;% 1 for ascending pass, 2 for descending pass (for plotting)
pad = 5; % degrees from storm center to extract data
%------------------------------------------

% CONSTANTS
Re = 6.371e6; % Mean Earth radius (m)

% EXTRACTING DATA
% mingmt: gmt time in minutes of day 
% windspd: wind speed magnitude (m/s)
% winddir: wind direction (relative to vertical, positive
%          anticlockwise, from 0-360 degrees)
% scatflag: (0=no rain, 1=rain)
% radrain: rain in km*mm/hr 

[mingmt,windspd,winddir,lat,lon] = get_scat_daily(fname,cntr,pad,pass);
nx = size(lat,1); % Number of points in x-dir (same as y-dir)

% CALCULATING ZONAL AND MERIDIONAL WINDS
u = windspd.*sind(winddir); % zonal wind speed (m/s)
v = windspd.*cosd(winddir); % meridional wind speed (m/s)

% CALCULATING POSITION ANGLE RELATIVE TO VERTICAL
pos_angle = zeros(nx,nx); % position angle (relative to vertical, positive
%                           anticlockwise, from 0-360 degrees). initializing with zeros
x = (lon - cntr(2))*2*pi*Re*cosd(mean(lat(:)))/360; % x-coordinate relative to center (m)
y = (lat - cntr(1))*2*pi*Re/360; % % y-coordinate relative to center (m)
radius = sqrt(x.^2 + y.^2); % Radius relative to center (m)

for i=1:nx
    for j=1:nx
        my_x = x(i,j);
        my_y = y(i,j);
        my_angle = atan(abs(my_x/my_y))*180/pi;
        
        % Ensuring that pos_angle varies from 0 to 360
        if my_x >= 0 && my_y >= 0 % first quadrant (0 to 90)
            pos_angle(i,j) = my_angle;
        elseif my_x > 0 && my_y < 0 % fourth quadrant (90 to 180)
            pos_angle(i,j) = 180 - my_angle;
        elseif my_x <= 0 && my_y < 0 % third quadrant (180 to 270)
            pos_angle(i,j) = 180 + my_angle;
        else % second quadrant (270 to 360)
            pos_angle(i,j) = 360 - my_angle;
        end
    end
end

%CALCULATE VR AND V_THETA 
vtheta = windspd .* sind(pos_angle - winddir); % Azimuthal wind speed (m/s)
vr = windspd .* cosd(pos_angle - winddir); % Radial wind speed (m/s)

% PLOTTING
    % WIND SPEED AS A FUNCTION OF LAT AND LON
figure
set(gcf,'color','w')
[c,h] = contour(lon,lat,windspd);
clabel(c,h);
hold on;
thet = 0:.05:2*pi;
for r = 1:pad
    rad = r*ones(1,length(thet));
    [xc,yc] = pol2cart(thet,rad);
    x_line = xc + cntr(2);
    y_line = yc + cntr(1);
    plot(x_line,y_line,'k:'); % distance rings (approximate)
end
plot(cntr(2),cntr(1),'kx');
axis square;
hold off;
title('Wind speed [m/s]');
xlabel('Longitude');
ylabel('Latitude');

    % WIND SPEED AS A FUNCTION OF RADIUS
figure 
set(gcf,'color','w')
scatter(reshape(radius,nx*nx,1)/1000,reshape(windspd,nx*nx,1),'+k');
title('Wind Speed')
xlabel('Radius [km]');
ylabel('Wind speed [m/s]')
set(gca,'xlim',[0 pad*110])

    % STREAMLINES
figure
set(gcf,'color','w')
h = streamslice(lon(:,1)',lat(1,:),u',v');
set(h,'linewidth',1.5)
axis square;
title(sprintf('Streamlines (%s)', storm));
xlim([min(lon(:)) max(lon(:))])
ylim([min(lat(:)) max(lat(:))])
xlabel('Longitude (°)')
ylabel('Latitude (°)')

    % AZIMUTHAL WIND SPEED AS A FUNCTION OF RADIUS
figure
set(gcf,'color','w')
scatter(reshape(radius,1,nx*nx)/1000,reshape(abs(vtheta),1,nx*nx),'+k');
xlabel('Radius [km]')
ylabel('v_{\theta} [m/s]')
title('Magnitude of azimuthal wind speed')
% set(gca,'xlim',[0 pad*110])

    % CONTOUR OF RADIAL VELOCITY (vr)
figure
set(gcf,'color','w')
cmin = -20; cint = 0.5; cmax = 20;
contourf(lon,lat,vr)
colormap(gca,bluepalered(cmin,cmax,cint))
colorbar
xlabel('Longitude (°)')
ylabel('Latitude (°)')
title('Radial wind speed v_r (m/s)')
caxis([cmin cmax])

    % CONTOUR OF AZIMUTHAL VELOCITY (vtheta)
figure
set(gcf,'color','w')
cmin = -20; cint = 0.5; cmax = 20;
contourf(lon,lat,vtheta)
colormap(gca,bluepalered(cmin,cmax,cint))
colorbar
xlabel('Longitude (°)')
ylabel('Latitude (°)')
title('Azimuthal wind speed v_{\theta} (m/s)')
caxis([cmin cmax])

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% FUNCTIONS

% DISPATCH BASED ON FILENAME
function [mingmt,windspd,winddir,lat,lon] = ...
    get_scat_daily(filename,cntr,pad,pass)

    if startsWith(filename, 'qscat')
        [mingmt,windspd,winddir,lat,lon] = ...
            get_qscat_daily(filename,cntr,pad,pass);
    elseif startsWith(filename, 'wsat')
        [mingmt,windspd,winddir,lat,lon] = ...
            get_windsat_daily(filename,cntr,pad,pass);
    else
        error('File name does not start with qscat or wsat')
    end

end

% READ QUICKSCAT FILES
function [mingmt,windspd,winddir,lat,lon] = ...
    get_qscat_daily(filename,cntr,pad,pass)   

    xscale=[6.,.2,1.5];

    if (exist(filename)==2)
        fid=fopen(filename,'r');
        for iasc=1:2
            for ivar=1:4
                dat=fread(fid,[1440 720],'uchar');
                switch ivar
                    case 1
                        dat =dat*xscale(ivar);
                        mingmt_raw(:,:,iasc)=dat;
                    case 2
                        dat=dat*xscale(ivar);
                        windspd_raw(:,:,iasc)=dat;
                    case 3
                        dat=dat*xscale(ivar);
                        winddir_raw(:,:,iasc)=dat;
                end   % switch
            end	  % ivar loop
        end    % iasc loop
        fclose(fid);

        index = find(mingmt_raw > 1440);
        windspd_raw(index) = -999.;
        winddir_raw(index) = -999.;
        
        lon_raw = 0.125:.25:359.875;
        lat_raw = -89.875:.25:89.875;
        elim = find(windspd_raw==-999);
        windspd_raw(elim) = NaN;
        winddir_raw(elim) = NaN;

        %FIND INDICES OF LAT/LON BOX
        lt1 = find(lat_raw <= (cntr(1)-pad),1,'last');
        lt2 = find(lat_raw >= (cntr(1)+pad),1,'first');
        ln1 = find(lon_raw <= (cntr(2)-pad),1,'last');
        ln2 = find(lon_raw >= (cntr(2)+pad),1,'first');

        % CALCULATE (APPROXIMATE) RADIUS FROM STORM CENTER 
        % FOR ALL POINTS WITHIN LAT/LON BOX
        [lat,lon]=meshgrid(lat_raw(lt1:lt2),lon_raw(ln1:ln2));
        mingmt = mingmt_raw(ln1:ln2,lt1:lt2,pass);
        windspd = windspd_raw(ln1:ln2,lt1:lt2,pass);
        winddir = winddir_raw(ln1:ln2,lt1:lt2,pass);

    else
        mingmt=[];windspd=[];winddir=[];
    end    % file exists
end

% READ WINDSAT FILES
function [mingmt,windspd,winddir,lat,lon] = ...
    get_windsat_daily(filename,cntr,pad,pass)   

    xscale=[6.,.2,1.5];

    if (exist(filename)==2)
        fid=fopen(filename,'r');
        for iasc=1:2
            for ivar=1:9
                dat=fread(fid,[1440 720],'uchar');
                switch ivar
                    case 1
                        dat =dat*xscale(1);
                        mingmt_raw(:,:,iasc)=dat;
                    case 8
                        dat=dat*xscale(2);
                        windspd_raw(:,:,iasc)=dat;
                    case 9
                        dat=dat*xscale(3);
                        winddir_raw(:,:,iasc)=dat;
                end   % switch
            end	  % ivar loop
        end    % iasc loop
        fclose(fid);

        index = find(mingmt_raw > 1440);
        windspd_raw(index) = -999.;
        winddir_raw(index) = -999.;
        
        lon_raw = 0.125:.25:359.875;
        lat_raw = -89.875:.25:89.875;
        elim = find(windspd_raw==-999);
        windspd_raw(elim) = NaN;
        winddir_raw(elim) = NaN;

        %FIND INDICES OF LAT/LON BOX
        lt1 = find(lat_raw <= (cntr(1)-pad),1,'last');
        lt2 = find(lat_raw >= (cntr(1)+pad),1,'first');
        ln1 = find(lon_raw <= (cntr(2)-pad),1,'last');
        ln2 = find(lon_raw >= (cntr(2)+pad),1,'first');

        % CALCULATE (APPROXIMATE) RADIUS FROM STORM CENTER 
        % FOR ALL POINTS WITHIN LAT/LON BOX
        [lat,lon]=meshgrid(lat_raw(lt1:lt2),lon_raw(ln1:ln2));
        mingmt = mingmt_raw(ln1:ln2,lt1:lt2,pass);
        windspd = windspd_raw(ln1:ln2,lt1:lt2,pass);
        winddir = winddir_raw(ln1:ln2,lt1:lt2,pass);

    else
        mingmt=[];windspd=[];winddir=[];
    end    % file exists
end

% COLORMAP FUNCTION 
function [ map ] = bluepalered(cmin,cmax,cint)
    n=floor(abs(cmin/cint));
    m=floor(abs(cmax/cint));
    r = [linspace(0,0.8,n) linspace(1,1,m)];
    g = [linspace(0.3,1,n) linspace(1,0.3,m)];
    b = [linspace(1,1,n) linspace(0.8,0,m)];
    map = [r(:), g(:), b(:)];
end
