%QSCATHURR.M
%Driver to generate figure for QuikSCAT winds near hurricane
%INPUT: fname - filename of QuikSCAT 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

clear all; 
%close all; clc

%USER DEFINED INPUT VARIABLES
%------------------------------------------
fname = 'qscat_20080711v4'; % filename of QuikSCAT Data
storm = 'Bertha'; % storm name 
cntr = [29.0 297.7];

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,scatflag,radrain,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)

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% EXERCISE
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 1. Estimate the Rossby number (Ro) as a function of radius. 
% Remember that on Earth, Ro depends on the Earth's rotation rate 
% as well as the sine of the latitude. 
% You will find that variables named 'lat' and 'radius' useful. 
% Note that those are both matrices, hence Ro will also be a matrix.

% 2. Plot Ro as a function of radius
% You may find the 'scatter' function useful. 
% This function takes vectors as inputs, so you will need to reshape
% your matrices. Take a look at how that's done below. Make sure you
% understand what the reshape command does. 

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

% 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

% FUNCTION TO READ IN THE BIN FILE
function [mingmt,windspd,winddir,scatflag,radrain,lat,lon]=get_scat_daily(filename,cntr,pad,pass)    %this subroutine will read the RSS scatterometer daily bytemaps.
    %reads version-3a files released October 2006
    %
    %input argument is the full path file name:
    %   get_scat_daily_v03(filename)
    %
    %output arguments:
    %   [mingmt,windspd,winddir,scatflag,radrain]
    %   mingmt is gmt time in minutes of day
    %   windspd in m/s	(10 meter surface wind)
    %   winddir in degrees 	(oceanographic convention, blowing North = 0)
    %   scatflag		(0=no rain, 1=rain)
    %   radrain   in km*mm/hr  (-999.0 = no collocation available)
    %                       (  -1.0 = adjacent cells had rain, but not this cell)
    %                       (   0.0 = radiometer data exist and show no rain)
    %				(   0.5 - 31.0 = columnar rain rate in km*mm/hr)
    %
    %  The center of the first cell of the 1440 column and 720 row map is at 0.125 E longitude and -89.875 latitude.
    %  The center of the second cell is 0.375 E longitude, -89.875 latitude.
    % 		XLAT=0.25*ILAT-90.125
    %		XLON=0.25*ILON-0.125
    %
    %please read the description file on www.remss.com
    %for infomation on the various fields, or contact RSS support:
    % http://www.remss.com/support

    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;
                    case 4
                        scatflag_raw(:,:,iasc) = dat-2*fix(dat/2); % bit 1
                        rad_flag = fix((dat-4*fix(dat/4))/2); % bit 2
                        temp=fix(dat/4); % bits 3-8

                        temp2 = -999*ones(size(dat));

                        for i=1:prod(size(temp2)),
                            if rad_flag(i)==1
                                if temp(i)==0
                                    temp2(i)=0;
                                elseif temp(i)==1
                                    temp2(i)=-1;
                                else
                                    temp2(i)=temp(i)/2-0.5;
                                end
                            end

                        end  % i loop

                        radrain_raw(:,:,iasc)=temp2;

                end   % switch
            end	  % ivar loop
        end    % iasc loop
        fclose(fid);

        index = find(mingmt_raw > 1440);
        windspd_raw(index) = -999.;
        winddir_raw(index) = -999.;
        scatflag_raw(index)= -999.;
        radrain_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)); % note that the indices in *grid will be (lon,lat)!
        mingmt = mingmt_raw(ln1:ln2,lt1:lt2,pass);
        windspd = windspd_raw(ln1:ln2,lt1:lt2,pass);
        winddir = winddir_raw(ln1:ln2,lt1:lt2,pass);
        scatflag = scatflag_raw(ln1:ln2,lt1:lt2,pass);
        radrain = radrain_raw(ln1:ln2,lt1:lt2,pass);

    else
        mingmt=[];windspd=[];winddir=[];scatflag=[];radrain=[];
    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
