function varargout = seems3m(varargin)
% WARING1a M-file for WARING1a.fig
%      Based on the SEEMS7 image analysis program
%
%      SEEMS3M, by itself, creates a new SEEMS3M or raises the existing
%      singleton*.
%
%      H = SEEMS3M returns the handle to a new SEEMS3M or the handle to
%      the existing singleton*.
%
%      SEEMS3M('CALLBACK',hObject,eventData,handles,...) calls the local
%      function named CALLBACK in SEEMS3M.M with the given input arguments.
%
%      SEEMS3M('Property','Value',...) creates a new SEEMS3M or raises the
%      existing singleton*.  Starting from the left, property value pairs are
%      applied to the GUI before seems3m_OpeningFunction gets called.  An
%      unrecognized property name or invalid value makes property application
%      stop.  All inputs are passed to seems3m_OpeningFcn via varargin.
%
%      *See GUI Options on GUIDE's Tools menu.  Choose "GUI allows only one
%      instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help seems3m

% Last Modified by GUIDE v2.5 06-Aug-2003 14:28:01

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name',       mfilename, ...
                   'gui_Singleton',  gui_Singleton, ...
                   'gui_OpeningFcn', @seems3m_OpeningFcn, ...
                   'gui_OutputFcn',  @seems3m_OutputFcn, ...
                   'gui_LayoutFcn',  [] , ...
                   'gui_Callback',   []);
if nargin & isstr(varargin{1})
    gui_State.gui_Callback = str2func(varargin{1});
end

if nargout
    [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
    gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT



% --- Executes just before seems3m is made visible.
function seems3m_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% varargin   command line arguments to seems3m (see VARARGIN)

%  *** Note that it appears that all variables passed through the handles
%  structure must be reserved in the open function
%  *** Note that it appears that both handles and appdata are write-once/read-many rather than write-many/read-many when applied to arrays.
%  Use global variables for the latter.

%  Version History
%   SEEMS6d reads, displays, and saves all 3 image types ok
%   SEEMS6e crops the ssh image and gets its color bar
%   SEEMS6f reads the ssh data to screen successfully
%   SEEMS6g fixes the color bar gaps
%   SEEMS6h distinguishes between out of bounds and overland pixes for ssh
%   SEEMS6i handles ssh transects including navigation correctly
%   SEEMS6m adds linear application of chl colorbar and not quite perfect chl navigation
%   SEEMS6n adds log applicaton of the chl colorbar and much improved chl navigation
%   SEEMS6o perfects the chl navigation to single-pixel accuracy
%   SEEMS6p reverses ssh scale (it's flipped on the image because it is a vertical scale
%   SEEMS6q cleans up the debugging statements
%   SEEMS7 is the operational version - identical to SEEMS6q
%   WARING1a is a SAR wind analysis program - based on SEEMS7
%   WARING1b is a cleaned up version eliminating extranious plots
%   WARING1c begins to add actual SAR analysis features - it updates the figure only when the SAR plot button is pushed
%   WARING1d adds a single plotted sar wind vector in response to user control
%   WARING1e partially adds archiving and multi-record editing - it is a running but not correct fallback snapshot
%   WARING1f completes this task
%   WARING1g adds automatic archive save and load for sar wind vectors
%   WARING1h adds radio buttons for the determination of which data source to edit and the ability to delete SAR winds
%   WARING1i adds the ability to grey out SAR winds when they aren't editable
%   WARING1j adds the ability to switch off SAR editing
%   WARING2a begins the process of adding QuikSCAT winds - does file load or default field generation
%   WARING2b adds the basic editor capbility for QuikSCAT winds
%   WARING2c adds the display of QuikSCAT winds
%   WARING2d gets the editor to actually change the color and status of the correct winds 
%   WARING2e adds a slider bar and text indicator for QuikSCAT and MM5 editing radius
%   WARING3a adds real data input but doesn't yet display correctly or handle image navigation correctly.
%   WARING3b gets the sar image to display as a greyscale
%   WARING3c gets the wind vector colors correct and the sar grayscale
%       adjusted so that weak wind features show up (a fourth root
%       enhancement curve standing in for cmod4)
%   WARING3d gets the lat lon text line reading correctly for the mouse
%       click and the lat lon range read from the input dataset.
%       The sar winds still don't plot and the quikscat winds have vanished
%   WARING3e gets the subjective sar wind direction plotter/editor/saver/reloader working again
%   WARING3f gets the scatterometer wind plotter/editor/saver/reloader working again
%   WARING3g redoes the image size to 600x660 with size and min latpixel
%       and min lon pixel determining the max lat and lon pixels rhather than
%       the other way around.
%   WARING3h add sliders to control the min lat and lon pixels
%   WARING3i allows those sliders to actually move the SAR image
%   WARING3j moves the Quikscat and subjective sar vectors similarly
%   WARING3k allows the edits SAR subjective to respond to the sliders as well.
%   WARING3l allows the QuikSCAT winds to be edited correctly and the lat lon reading at the left of the screen and
%       in the subjective SAR file to be correct
%   WARING3m incorporates my wind direction interpolator and has the
%       SAR image respond to the edited winds
%   WARING3n is identical except for improved inline documentation
%   WARING4a adds a 2x2 averaging to shrink array size from 1200 to 600.  Also adds loading of the the real cmod-calculated winds vro sigmatowind2a_init.m.
%   WARING4b should add real cmod-calculated winds in place of my cos function in the wind recalculator.
%   WARING4c is a bug-fix version wherein the blended wind directions overright the intial sar wind directions and so are used in the recomputation of sar wind speed
%   WARING4d adds file saving of all the original input variables
%   WARING4e adds u and v output and fixes a sign error in v computation from wind direction that was flipping that component of the blended
%       wind diretion.
%   WARING4f cleans up some of the comments and removes a few unnecessary print statements.
%   WARING4g is the operational CMOD4 version
%   WARING5a adds basic CMOD5 subroutine but does not call it
%   WARING5b Calls the CMOD subroutine - works fine if the CMOD 5 response
%       function has rcs as a monotonic function of wind speed but crashes the
%       interplation if it doesn't.
%   WARING5c goes over to a slow-to-fast search before implementing the
%       interpolation - included code to create an aritifical non-monotonic
%       trend for testing the the slow-to-fast search - it works
%   WARING5d - eliminates the test code - note that all versions up through
%       this use a buggy early version of the JHUAPL CMOD 5 subroutine
%   WARING5e - replaces the CMOD 5 subroutine with an updated and corrected
%       version from JHUAPL
%   WARING5f - adds a check for NaN winds (i.e. backscatter greater than
%       any theoretically possible) and replaces such values with the maximum
%       wind possible.
%   WARING5g - improves the check for backscatter above theory and adds a
%       check for backscatter below theory.
%   WARING5h - attempts to solve the CMOD 5 problem, among other things
%   sampling sardir instead of averaging it (a 0 / 360 issue)
%   WARING5i - a debug version in which sardir is passed through untouched.
%   WARING5j - the CMOD4 version of WARING5h - it also includes a crude
%       patch in cmod4hyb2a to treat all theta values smaller than that for
%       which the program was designed as the smallest value for which it was
%       designed.  A crude hack at best.  It also includes a
%       predimensioning of the speed array so that the ENVI sat images that
%       don't fill the 1200x1200 domain don't result in speed images that
%       are smaller than 600x600.  
%   waring5k - the CMOD4 version as in WARING5J, but with all variables
%       decimated instead of averaged.
%   WARING5L = uses new cmod5
%   WARING5L5 = starts to eliminate the excess code and output


%  Clear the screen
clc

%  Clear out the main image file - in case there was a previous run
global compositeim
clear compositeim

% Set variables used for image display croping and display 
handles.minLatPixel     =     1;
handles.maxLatPixel     =   600;
handles.minLonPixel     =     1;
handles.maxLonPixel     =   600;

% The lower left corner variables show where the image is placed within the window.
handles.lowerleftX      =   350;  
handles.lowerleftY      =    50;

% Set default values for the SAR image geography.  These will be
% overwritten later so this block of code is probably redundant
global maxLat minLat maxLon minLon
minLat          =  25.2;
maxLat          =  44.5;
minLon          = -81.5;
maxLon          = -57.5;

%  Specify the size of the part of the SAR image that is actually shown on
%  screen at any one time and compute where the upper right hand corner of
%  the image is in the window.
global imageWidth imageHeight
imageWidth = 600;
handles.maxLonPixel = imageWidth + handles.minLonPixel - 1;
imageHeight = 600;
handles.maxLatPixel = imageHeight + handles.minLatPixel - 1;

% Specify the defualt value of the image offset - i.e. the position of the
% image before you start moving it around with the sliders
global sarLatShift sarLonShift
 sarLatShift = 0;
 sarLonShift = 0;

% Reserve space for a click counter
global clickCount
clickCount = 0;

% Set default for the selection radius for existing wind values 
global clickRadius changeRadius
clickRadius = 10;  % For SAR in pixels
changeRadius = 1; % For MM5 and QuikSCAT - in degrees lat and lon

% Reserve space for the variables used in text output
global lat lon statusText

% Reserve space for the cropped SAR image
global sarImage
sarImage = zeros(handles.maxLatPixel-handles.minLatPixel+1,handles.maxLonPixel-handles.minLonPixel+1,3); 

% Specify the SAR wind vector colors - note that when using SAR data the
% image array is of type double so these colors need to be in the range
% 0-1.
global sarNewColor sarOldColor sarInactiveColor sctInactiveColor
sarNewColor = [60,120,240]/256;
sarOldColor = [150, 180, 240]/256;
sarInactiveColor = [240, 240, 240]/256;
sctInactiveColor = [240, 240, 240]/256;

% Choose default command line output for seems3n
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

%  Set up the SAR image via a function call so I can reuse code
%   First load the raw SAR and Scattermometer data from Pete
disp('about to loadInputFcn')
loadInputFcn(hObject, eventdata, handles);
%   Second see if there is an edited Scatterometer file available for use
%   instead
disp('about to loadSctFcn')
loadSctFcn(hObject, eventdata, handles);
%   Third see if there are SAR subjective winds to be loaded
disp('about to loadSarFcn')
loadSarFcn(hObject, eventdata, handles);
%  Compute the direction-dependent wind speed from the sar radar cross section
disp('about to calcSarWindFcn')
calcSarWindFcn(hObject, eventdata, handles);
%  Display the SAR wind speed with Scatterometer and subjective SAR wind
%  directions overlaid.
disp('about to showSarFcn')
showSarFcn(hObject, eventdata, handles);

disp('about to set up SAR direction slider')
% Slider setup for SAR subjective wind direction
slider_step(1) = 5/360;
slider_step(2) = 0.25;
sarWindDirection = 180;
set(handles.SarDirectionSlider,'sliderstep',slider_step,'max',360,'min',0,'Value',sarWindDirection);
%  Set the status text under the slider
sarDirectionText = strcat('SAR Wind Direction => ',int2str(round(sarWindDirection)));
set(handles.sarDirectionText,'string',sarDirectionText); 

disp('about to set up editor radiaus slider')
% Slider setup for edit brush radius
slider_step(1) = 0.01;
slider_step(2) = 0.1;
set(handles.editRadiusSlider,'sliderstep',slider_step,'max',2,'min',0,'Value',changeRadius);
%  Set the status text under the slider
numFigures = 2;
editRadiusText = strcat('Edit Radius => ',num2str(changeRadius,numFigures));
set(handles.editRadiusText,'string',editRadiusText); 

% Set defaults for showing the various wind vector types - including setting the radio buttons
global showSAR showMM5 showSCT
showSAR = 1;
showMM5 = 0;
showSCT = 0;
set(handles.SARbutton,'Value',showSAR);
set(handles.MM5button,'Value',showMM5);
set(handles.QuikSCATbutton,'Value',showSCT);

% Set up the speed array so it won't shrink if the image doesn't span it
global speed
speed = zeros(handles.maxLatPixel-handles.minLatPixel+1,handles.maxLonPixel-handles.minLonPixel+1); 

disp('Editor is ready for use')

% Choose default command line output for seems3m
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes seems3m wait for user response (see UIRESUME)
% uiwait(handles.figure1);



% --- Outputs from this function are returned to the command line.
function varargout = seems3m_OutputFcn(hObject, eventdata, handles)
% varargout  cell array for returning output args (see VARARGOUT);
% hObject    handle to figure
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% Get default command line output from handles structure
% ---- Don't mness with this subroutine ----
varargout{1} = handles.output;



% --- Executes on button press in closeButton.
function closeButton_Callback(hObject, eventdata, handles)
% hObject    handle to closeButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global sarWind sarCount
global sctSpd sctDir sctStatus sctLat sctLon sctCount
global qsdir qslat qslon qsspd sardir sarlat sarlon sarlookdir sarrcs sarspeed sartheta 
global speed
% Get the current position of the GUI from the handles structure
% to pass to the modal dialog.
pos_size = get(handles.figure1,'Position');
% Call modaldlg with the argument 'Position'.
user_response = closedlg('Title','Confirm Close');
switch user_response
case {'No'}
% take no action
case 'Yes'
% Prepare to close GUI application window
%                  .
%                  .
%                  .
%  Close the GUI application window and all its subwindows
close all
%  Compute the u and v components of the final wind for saving to disk.
    radianWindDirection = pi * (sardir) / 180; % Align so that the wind direction is FROM.
    u = speed .* (-sin(radianWindDirection));
    v = speed .* (-cos(radianWindDirection));
    disp('about to save output')
%  Save the subjective sar winds
save('sarwinds', 'sarWind', 'sarCount')
%  Save the edited Scatterometer winds
save('sctwinds', 'sctSpd', 'sctDir', 'sctLat', 'sctLon', 'sctStatus', 'sctCount')
%  Save the edited edited sardir and sarspd along with all the other variables in Pete's original SAR input file.
save('outputData', 'qsdir', 'qslat' ,'qslon', 'qsspd', 'sardir', 'sarlat', 'sarlon', 'sarlookdir', 'sarrcs', 'sarspeed', 'sartheta', 'speed', 'u', 'v')
end


% --- Used to set up the SAR direction slider
% --- Executes during object creation, after setting all properties.
function SarDirectionSlider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to SarDirectionSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called
% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



% --- Used to respond to the SAR direction slider
% --- Executes on slider movement.
function SarDirectionSlider_Callback(hObject, eventdata, handles)
global sarWindDirection irec sarWind
% hObject    handle to SarDirectionSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
%  Get slider value
sarWindDirection = get(handles.SarDirectionSlider, 'value');
%  Store the slider value in the wind direction element of the sawWind data
%  strucutre
sarWind(irec,5) = sarWindDirection;
%  Set the status text under the slider
sarDirectionText = strcat('SAR Wind Direction => ',int2str(round(sarWindDirection)));
set(handles.sarDirectionText,'string',sarDirectionText); 
%  Replot the wind records
showSarFcn(hObject, eventdata, handles);
%  Update the handles structure
guidata(hObject, handles);



% --- SAR subjective wind data loading function
function loadSarFcn(hObject, eventdata, handles)
global sarImage sarWind sarCount
% hObject    handle to closeButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
%
% Display an initial image as a banner - **** note that none of the handles
%    reset here seem to get OUT despite similar ones doing so in other
%    functions.  Therefore I use the application data approach.
%    It looks like the handles only work out of callback functions but that
%    they work in and out of all functions.  I may have this backwards.
%
if (exist('sarwinds.mat') == 2)  % If the SAR subjective winds file exists, load it
    display 'Loading an existing SAR subjective wind file'
    load('sarwinds')  
else % Note that there are no SAR subjective winds so far
    display 'SAR subjective wind file does not yet exist'
    sarCount = 0;
end
guidata(hObject, handles);



% --- raw data loading function
function loadInputFcn(hObject, eventdata, handles)
global qslon qslat qsspd qsdir qstheta sarrcs sardir sarspeed sarlon sarlat sarlookdir sartheta
global maxLat minLat maxLon minLon
global sarHeight sarWidth
global ni nj w nimage
% hObject    handle to closeButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
%
if (exist('inputData.mat') == 2)  % If the SAR backscatter file exists, load it
    display 'Loading an existing input data file file'
    load('inputData')
else % Note that there are SAR backscatter data available
    display 'Input data file does not exist'
    sarCount = 0;
end
set(hObject, 'Units', 'pixels'); 

%  Average the sar variables over a 2x2 array to get them to fit in memory better
[isize,jsize] = size(sardir)
size(sartheta)
sarlat = sarlat(1:2:isize-1);
sarlon = sarlon(1:2:isize-1);
% sardir might be better done as a sampling (i.e. just the first term) than
% as a mean to avoid 0 / 360 problems
sardir = sardir(1:2:isize-1,1:2:isize-1);
sarspeed = sarspeed(1:2:isize-1,1:2:isize-1);
sarrcs = sarrcs(1:2:isize-1,1:2:isize-1);
sartheta = sartheta(1:2:isize-1,1:2:isize-1);
[iblock,jblock] = size(sarspeed);
pack

%  Locate the non-zero pixels (i.e. those that are actually in the swath)
[ni,nj] = find(sarrcs > 0.0);
w = find(sarrcs > 0.0);
nimage=size(w);

%  Reset image lat lon range to match the current dataset
maxLat = max(sarlat)
minLat = min(sarlat)
maxLon = max(sarlon)
minLon = min(sarlon)
[sarHeight,sarWidth] = size(sarrcs)
% Update handles structure
guidata(hObject, handles);



% --- QuikSCAT loading fuction -- it uses the raw data if there are no edited data available
function loadSctFcn(hObject, eventdata, handles)
global sctDir sctSpd sctLat sctLon sctStatus sctCount
global maxLat minLat maxLon minLon
global qslon qslat qsspd qsdir qstheta sarrcs sardir sarspeed sarlon sarlat sarlookdir sartheta
if (exist('sctwinds.mat') == 2)  % If the subjectively analyzed scattermometer wind file exists, load it
    display 'Loading an existing QuikSCAT wind file'
    load('sctwinds')
else % Note that there are no SAR subjective winds so far
    display 'QuikSCAT wind file does not exist'
    % Figure out which hemisphere one is in and correct the Scatterometer
    % to +180/-180 range used in Meteorology and Geography
    qshemisphere = (sign(qslon-180)+1)/2;
    qslon = (qshemisphere.*(qslon-360))+((1-qshemisphere).*qslon);
    % Determine how many Scatterometer observations there are
    [ymax,xmax] = size(qslat)
    % Pick out just those Scatterometer observations within the SAR image
    sctCount = 0
    for i=1:ymax
        for j=1:xmax
            if ((qslat(i,j) < maxLat) && (qslat(i,j) > minLat) && (qslon(i,j) < maxLon) && (qslon(i,j) > minLon))
                sctCount = sctCount + 1;
                sctSpd(sctCount) = qsspd(i,j);
                sctDir(sctCount) = mod(qsdir(i,j)+180,360); % This 180 flip may be case specific if our test case just happens to have the ambiguity backwards.
                sctStatus(sctCount) = 1;  %  Set all scatterometer winds to be ON by default
                sctLat(sctCount) = qslat(i,j);
                sctLon(sctCount) = qslon(i,j);
            end
        end
    end
end
guidata(hObject, handles);


% --- Primary image display function
function showSarFcn(hObject, eventdata, handles)
global sarImage sarCount sarWind compositeim sarWindDirection clickX clickY clickCount 
global sarNewColor sarOldColor sarInactiveColor sctInactiveColor
global showSAR showMM5 showSCT irec imageWidth imageHeight
global sctLat sctLon sctDir sctSpd sctStatus sctCount
global maxLat minLat maxLon minLon
global sarHeight sarWidth sarLonShift sarLatShift sarImageAll
% hObject    handle to closeButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
%
% Display an initial image as a banner - **** note that none of the handles
%    reset here seem to get OUT despite similar ones doing so in other
%    functions.  Therefore I use the application data approach.
%    It looks like the handles only work out of callback functions but that
%    they work in and out of all functions.  I may have this backwards.
set(hObject, 'Units', 'pixels'); 

%  Crop the SAR image so that just the current part shows on screen - ****
%  this is the line that crashes, because my image got cropped for ENVISAT,
%  but not for the previous satelite.  Rats.  ****
sarImage = flipud(sarImageAll(handles.minLatPixel+sarLatShift:handles.maxLatPixel+sarLatShift,handles.minLonPixel+sarLonShift:handles.maxLonPixel+sarLonShift,:)); % Crop image
[cropheight, cropwidth, junk]=size(sarImage);
%  Add subjective SAR wind directions to the cropped SAR image
sarScale = (max(max(sarImage)));
compositeim(:,:,1) = (sarImage)/sarScale;
compositeim(:,:,2) = (sarImage)/sarScale;
compositeim(:,:,3) = (sarImage)/sarScale;
vecLength = 10;
headRadius = 1;
for jrec = 1:sarCount % Loop over all existing sar wind records
    radianWindDirection = pi * (sarWind(jrec,5)-180) / 180; % Align so vector points DOWN wind
    u = sin(radianWindDirection);
    v = -cos(radianWindDirection);
    if (showSAR == 1)  % Use shades of blue if SAR is being edited
        if (jrec == irec)  % Set the appropriate color for the current wind record
            sarColor = sarNewColor;
        else
            sarColor = sarOldColor;
        end
    else  % Otherwise use greyed out shades of blue for SAR if it isn't being edited
        sarColor = sarInactiveColor;
    end
    %  Plot vector head
    for i = -headRadius:headRadius
        for j = -headRadius:headRadius
            yPos = sarWind(jrec,2)+j + sarLatShift;
            xPos = sarWind(jrec,1)+i - sarLonShift;
            if ((xPos > 0) && (yPos > 0) && (xPos <= imageWidth) && (yPos <= imageHeight))
                compositeim(yPos,xPos,:)=sarColor;
            end
        end
    end
    %  Plot vector body trailing off DOWN wind
    for i = 0:vecLength
        x = round(i * u);
        y = round(i * v);
        yPos = sarWind(jrec,2)+y + sarLatShift;
        xPos = sarWind(jrec,1)+x - sarLonShift;
        if ((xPos > 0) && (yPos > 0) && (xPos <= imageWidth) && (yPos <= imageHeight))
            compositeim(yPos,xPos,:)=sarColor;
        end
    end
end
% Add QuikSCAT wind directions to the cropped SAR image
vecLength = 10;
headRadius = 1;
%  Determine what faction of the SAR image is visible
fracHeight = imageHeight/sarHeight;
fracWidth = imageWidth/sarWidth;
%  Compute the shifted range of lat and lon visible in the cropped image
deltaLat = fracHeight * (maxLat - minLat);
deltaLon = fracWidth * (maxLon - minLon);
%  Compute the clipped image's minmum lat and lon
minLatCrop = minLat + ((maxLat - minLat) * (sarLatShift/sarHeight));
minLonCrop = minLon + ((maxLon - minLon) * (sarLonShift/sarHeight));
for i=1:sctCount  % Loop over all existing QuikSCAT wind records
    radianWindDirection = pi * (sctDir(i)-180) / 180; % Align so vector points DOWN wind
    u = sin(radianWindDirection);
    v = -cos(radianWindDirection);
    if (showSCT == 1)  % Use shades of blue if QuikSCAT is being edited
        if (sctStatus(i) == 1)
            sarColor = sarNewColor;
        else
            sarColor = sarOldColor;
        end
    else  % Otherwise use greyed out shades of blue for QickSCAT if it isn't being edited
        sarColor = sctInactiveColor;
    end
    %  Figure out the pixel values for the lat and lon of each QuikSCAT point
    sctY = imageHeight - round(((sctLat(i)-minLatCrop)/deltaLat)*imageHeight);
    sctX = round(((sctLon(i)-minLonCrop)/deltaLon)*imageWidth);
    %  Plot vector head
    for ih = -headRadius:headRadius
        for jh = -headRadius:headRadius
            yPos = sctY+ih;
            xPos = sctX+jh;
            if ((xPos > 0) && (yPos > 0) && (xPos <= imageWidth) && (yPos <= imageHeight))
                compositeim(yPos,xPos,:)=sarColor;
            end
        end
    end
    %  Plot vector body trailing off DOWN wind
    for k = 0:vecLength
        x = round(k * u);
        y = round(k * v);
        yPos = sctY+y;
        xPos = sctX+x;
        if ((xPos > 0) && (yPos > 0) && (xPos <= imageWidth) && (yPos <= imageHeight))
            compositeim(yPos,xPos,:)=sarColor;
        end
    end
end
% Display the SAR image with all its overlays
axes(handles.axis1);
image(compositeim);
disp('Screen update completed')
set(handles.axis1,'Visible', 'off','Units', 'pixels','Position', [handles.lowerleftX handles.lowerleftY cropwidth cropheight]);
guidata(hObject, handles);


% --- This is the main user interaction fuction
%   - It controls the response when the user clicks on the SAR image
% --- Executes on mouse press over figure background, over a disabled or
% --- inactive control, or over an axes background.
function figure1_WindowButtonDownFcn(hObject, eventdata, handles)
global sarImage sarCount sarWind statusText  clickCount clickRadius irec imageWidth imageHeight showSAR showMM5 showSCT
global deltaLat deltaLon latVal lonVal clickX clickY
global maxLat minLat maxLon minLon
global sarLatShift sarLonShift sarHeight sarWidth
% hObject    handle to figure1 (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% Increment click counter
% Get data for current mouse click, clipping to image area)
clickPos = get(hObject,'CurrentPoint');
%  These lat/lons need to be scaled by the image lat lon range
%  SAR pixel to lat lon conversion
deltaLat = maxLat - minLat;
deltaLon = maxLon - minLon;
latVal = min(sarHeight,max(1,round(clickPos(1,2))-handles.lowerleftY+sarLatShift));  
lonVal = min(sarWidth,max(1,round(clickPos(1,1))-handles.lowerleftX+sarLonShift));  
clickX = lonVal;
clickY = imageHeight - latVal; % flip click because mouse coordinates are upside down form image coordinates in MATLAB
%  Convert lat and lon from pixels to degrees
lonVal = minLon + (lonVal*deltaLon/sarWidth);
latVal = minLat + (latVal*deltaLat/sarHeight);
handles.newLat = latVal;
handles.newLon = lonVal;
% Print out values for current point: Lat, Lon
statusText = [num2str(latVal,4), ',    ',num2str(lonVal,4)];
set(handles.latestPointText,'string',statusText); 
%  Edit the appropriate data type
if (showSAR == 1)    
    sar_EditFcn(hObject, eventdata, handles);
elseif (showSCT == 1)
    sct_EditFcn(hObject, eventdata, handles);
elseif (showMM5 == 1)
end
%  Update the handles structure
guidata(hObject, handles);



% --- Edit the subjective SAR winds
function sar_EditFcn(hObject, eventdata, handles)
global sarImage sarCount sarWind statusText clickCount clickRadius irec imageWidth imageHeight showSAR showMM5 showSCT
global deltaLat deltaLon latVal lonVal clickX clickY
clickCount = clickCount + 1;
% Open existing SAR wind record or create a new one as needed
irec = 0;
if (sarCount > 0)  % Check to see if there are any existing winds the user might wish to edit
    for i=1:sarCount  % See if the user has chosen to edit an existing wind
        clickDistance = sqrt((clickX-sarWind(i,1))^2+(clickY-sarWind(i,2))^2);
        if clickDistance < clickRadius
            irec = i;  % Note which record is being edited
        end
    end
end
if (irec == 0) % Set up values for the case where the user wants to create a new wind record
    sarCount = sarCount + 1;
    sarWind(sarCount,1) = clickX;
    sarWind(sarCount,2) = clickY;
    sarWind(sarCount,3) = lonVal;
    sarWind(sarCount,4) = latVal;
    sarWind(sarCount,5) = 180;  %  Set a default wind direction
    irec = sarCount; %  Set the current record pointer to the newly created record
end
%  In either case, set the slider to reflect the data record (old or new)
set(handles.SarDirectionSlider, 'value', sarWind(irec,5));
sarDirectionText = strcat('SAR Wind Direction =>',int2str(round(sarWind(irec,5))));
set(handles.sarDirectionText,'string',sarDirectionText); 
%  Plot the newly created or edited wind record
showSarFcn(hObject, eventdata, handles);
%  Update the handles structure
guidata(hObject, handles);



% --- Edit the Scatterometer winds
function sct_EditFcn(hObject, eventdata, handles)
%global sarImage sarCount sarWind statusText clickCount clickRadius irec imageWidth imageHeight showSAR showMM5 showSCT
global changeRadius sctDir sctSpd sctLat sctLon sctStatus sctCount
global deltaLat deltaLon latVal lonVal clickX clickY
% Check the QuikSCAT records to see which ones to plot
for i=1:sctCount
    % See if the wind is close enough to need changing - note that this is an oval with the same lat and lon radius in degrees rather than in km.
    clickDistance = sqrt((latVal-sctLat(i))^2+(lonVal-sctLon(i))^2);
    if (clickDistance < changeRadius)
        sctStatus(i) = 1 - sctStatus(i);  % toggle the status
    end
end
%  Plot the edited wind records
showSarFcn(hObject, eventdata, handles);
%  Update the handles structure
guidata(hObject, handles);



% --- Executes on button press in sarPlotButton.
function sarPlotButton_Callback(hObject, eventdata, handles)
% hObject    handle to sarPlotButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
%  Redisplay the SAR image and all its various overlays
showSarFcn(hObject, eventdata, handles);



% --- Executes on button press in cmodRecomputeButton.
function cmodRecomputeButton_Callback(hObject, eventdata, handles)
% hObject    handle to cmodRecomputeButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% --- Does the SAR wind computation via CMOD
calcSarWindFcn(hObject, eventdata, handles);
% --- Redisplays the updated SAR image with all its overlays
showSarFcn(hObject, eventdata, handles);


% --- Does the direction-dependent SAR wind computation via a call to CMOD
function calcSarWindFcn(hObject, eventdata, handles)
global sarrcs sarImageAll sarlat sarlon sardir sarspeed sarlookdir sartheta
global sctDir sctSpd sctLat sctLon sctStatus sctCount
global sarWind sarCount
global sigma
global wind_speed wind_dir theta rlook alpha sigmavec
global ni nj w nimage
global speed
% Merge scat and subjective sar data from speed, direction, lat and lon
dataCount = 0;
for i=1:sctCount  
    if (sctStatus(i) > 0) % Load the active Scatterometer winds into the analysis arrays
        dataCount = dataCount + 1;
        dir(dataCount) = sctDir(i);
        xi(dataCount) = sctLon(i);
        yi(dataCount) = sctLat(i);
    end
end
for i = 1:sarCount  % load the subjective sar winds into the analysis arrays
    dataCount = dataCount + 1;
    dir(dataCount) = sarWind(i,5);
    xi(dataCount) = sarWind(i,3);
    yi(dataCount) = sarWind(i,4);
end
% Compute u and v data 
%   - Note that this ignores wind speed, it merely computes the direction vector components.
%       Thus, a strong jet won't dominate the analysis in it's local area.
%       And the speedless SAR subjective winds can be used.
dir = dir*pi/180; 
u = -sin(dir);
v = -cos(dir);
%  tic and toc serve as timers so one can monitor the performance of these
%  compuationally intensive routines on your own system
tic
%  Create the output grid
[XI,YI] = meshgrid(sarlon,sarlat);
meshtime = toc
tic
%  Analyze the QuikSCAT u and v onto the grid
u_ZI_nn = griddata(xi,yi,u,XI,YI,'nearest');
v_ZI_nn = griddata(xi,yi,v,XI,YI,'nearest');
nearest_neighbor = toc
% %  Display the result
% figure
% image(flipud(u_ZI_nn))
% title('nearest neighbor u')
% %  Display the result
% figure
% image(flipud(v_ZI_nn))
% title('nearest neighbor v')
tic
%  Analyze the QuikSCAT wind speeds onto the grid
u_ZI_bc = griddata(xi,yi,u,XI,YI,'cubic');
v_ZI_bc = griddata(xi,yi,v,XI,YI,'cubic');
bi_cubic = toc
% %  Display the result
% figure
% image(flipud(u_ZI_bc))
% title('bi-cibic u')
% figure
% image(flipud(v_ZI_bc))
% title('bi-cibic v')
%  Create merged bc and nn product 
%   - tricky code to get around all the
%   NaNs inherent in bi-cubic interpolation from a grid that doesn't extend
%   to the edges of the domain
t3 = 1000;
u_nanflag = isnan(u_ZI_bc);                     % Determine which points are missing
u_ZI_bc = max(u_nanflag,u_ZI_bc+t3);          % Replace the missing points with 1s as place holders
u_ZI_comp = max(u_ZI_bc,u_nanflag.*(u_ZI_nn+t3));
u_ZI_comp = u_ZI_comp - t3;                 % Remove the offset 
%  Do some diagnostic output
% figure
% image(flipud(u_ZI_comp))
% title('composite of nn and bc u')
v_nanflag = isnan(v_ZI_bc);                     % Determine which points are missing
v_ZI_bc = max(v_nanflag,v_ZI_bc+t3);          % Replace the missing points with 1s as place holders
v_ZI_comp = max(v_ZI_bc,v_nanflag.*(v_ZI_nn+t3));
v_ZI_comp = v_ZI_comp - t3;                 % Remove the offset
% figure
% image(flipud(v_ZI_comp))
% title('composite of nn and bc v')
dir_ZI_comp = mod(180+atan2(u_ZI_comp,v_ZI_comp)*180/pi,360);
figure
image(flipud(dir_ZI_comp)/2)
title('composite wind direction / 2')
colormap(jet(180))
colorbar('horiz')
xlabel('Feel free to delete this figure when you no longer need it.')





% ****************** Begin CMOD calling code
%  Start by overwriting the SAR wind direction array with the composite direction array
disp('about to overwrite sardir')
size(sardir)
size(dir_ZI_comp)
sardir = dir_ZI_comp;

disp('about to begin new cmod code')
%
% For test run 1, a value of 35 m/s is good. A better 
% alternative is to come up with a smart land mask that
% would allow you to avoid calculating winds over land.
%
maxwind = 35.0;
%
% 0.2 m/s spacing for wind_speed is arbitrary and some
% time savings may be had by going coarser (at the 
% expense of accuracy). May want to play around with
% this some. 
% 1 m/s is more than enough to prevent degradations inimage quality.
%wind_speed = (0:0.2:maxwind);  % Finer than need be for any wind above 2 m/s
%wind_speed = (0:2.0:maxwind); % Perhaps marginally too course a resolution
wind_speed = (0:1.0:maxwind);  % Works well for any wind above 2 m/s
%
rlook = sarlookdir;
%
% I am using alpha = 0.6 because that is what we use
% at APL. Paris Vachon would argue for 1.0. We can 
% adapt to the science as Envisat's interleaved polarization
% mode answers the question of what it should really be.
%
alpha = 0.6;
nspeed = size(wind_speed);
nx = max(ni)-min(ni)+1;
ny = max(nj)-min(nj)+1;
sigma = zeros(nspeed(2),nx,ny);
%  Do some preliminary housekeeping to set up loop limits and such.
np=5;
n = size(sigma);
ny = n(2);
nz = n(3);
%
% The following for loop is critical. The objective is to
% build up a large sigma array that contains sigma values
% for each wind speed, theta, wind direction triad (holding
% sarlook and alpha constant). The
% product of this loop is an array of sigma = f(u,theta,phi)
% where u = the wind speed from 0:maxwind incremented by 
%           whatever increment is desired.
%   theta = the incidence angles at each point. IMPORTANT: make
%           sure that the theta values that you use match 
%           exactly with the image theta values. If the 
%           theta values get screwed up here, it is impossible
%           to fix later.
%     phi = The wind direction in meteorological coordinates.
%           Just like theta need to make sure these match up with
%           the image values.
%
% Once these are defined, cmod is run for each value to
% obtain a corresponding sigma value.
%   *** GSY vectorized the wind loopl.  Getting a factor of 20 to 40
%   speed-up.  Vectorizing the image loop might actually have yielded
%   even greater gains but would have required coming up with a new way
%   to deal with pixels with zero RCS. ***
%
tic
for j = 1:nimage(1)
    theta = sartheta(ni(j),nj(j));
    wind_dir = sardir(ni(j),nj(j));
    % Uncomment whichever CMOD version you want
%     sigmavec = cmod4hyb(wind_speed, wind_dir, theta, r_look, alpha);
    sigmavec = cmod5hyb1c(wind_speed, wind_dir, theta, rlook, alpha);
    %  Then store the sigma vector
    sigma(:,ni(j),nj(j)) = sigmavec;
end
cmod_time = toc
%
%  Use linear interpolation to look up the right wind in the sima table.
tic
speed = zeros(handles.maxLatPixel-handles.minLatPixel+1,handles.maxLonPixel-handles.minLonPixel+1); 
% isn't shrunk in processing. ****
gsyn = max(size(wind_speed));
for j = 1:nimage(1)
    temp = sigma(:,ni(j),nj(j));
    utemp = wind_speed(:);
    xi = sarrcs(ni(j),nj(j));
    % yi = interp1(x,Y,xi,method)
    % yi = interp1(x,Y,xi,'linear')  - note that 'spline' would also work
    % Interpolating speed in rcs space - but rcs need not be monotonic in
    % cmod 5!
    % Check for non-monotonic behavior in rcs(wind_speed) and truncate
    % series at first maxima if needed
%     size(temp)
    gsydif = temp(1:gsyn-1)-temp(2:gsyn);
%     size(gsydif)
    utemp2 = utemp;
    %  3-15-2006 - this could probably be replaced with a simple sort and
    %  pulling the index of the max value out.
    monotonic_check = abs(mean(sign(gsydif)));
    gsyend = length(utemp);
    if (monotonic_check < 1.0)
        gsyend = 1;
    %   Non-monotonic backscatter - wind relationship  
    %   Search through RCS list for a match, starting at low wind speed
        for i = 2:gsyn-1
            %  If trend still hasn't revesed store the result
            if (gsydif(i)*gsydif(i-1)) > 0
                gsyend = i;
            end
        end
        %  Truncate the series at the first maxima
        temp = temp(1:gsyend);
        utemp2 = utemp(1:gsyend);
    end
    %                            x    y     xi
    gsytempwind =        interp1(temp,utemp2,xi,'linear');
    %  Replace NaN with the saturation valueof wind speed in those cases where observed
    %  backscatter is greater than any value theoretically possible for any
    %  lesser speed  
    if (gsytempwind >= 0.0)
        speed(ni(j),nj(j)) = gsytempwind;
    elseif (xi < min(temp))
        %  Set wind to zero if backscatter is below theoretical minimum
        speed(ni(j),nj(j)) = 0.0;
    else
        %  There are two alternative responses to not finding the observed
        %  backscatter: 1)use the fastest wind on the scale, 2) use the
        %  wind corresponding to the biggest backscatter
        %speed(ni(j),nj(j)) = max(utemp); % Nominally 35 m/s although this constant can be changed above
        speed(ni(j),nj(j)) = utemp(gsyend); % Speed corresponding to the maximum backscatter.
    end
end
linearlookuptabletime = toc
%
% "Saturate" the image by taking those places where the
% RCS values are greater than the maximum value in sigma
% and setting them to maxspeed. This is effectively
% clipping the image at max wind.
%
% Again, this may not be the smartest thing in the world
% to do.
%
ww=find(sarrcs > max(max(max(sigma))));
speed(ww) = maxwind;
% ***************** End CMOD caalling code


sarImageAll = max(speed,0); %  Eliminate the slightly negative approximations to zero
%  Scale image to make the weak areas of the wind field visible by using an exponent less than 1
%  Nominal exponent is 1.
scaleFactor = 1;
sarImageAll = (sarImageAll).^scaleFactor;
sarMin = min(min(sarImageAll));
sarMax = max(max(sarImageAll));
sarRange = sarMax-sarMin;
sarImageAll = (sarImageAll-sarMin)/sarRange;

%  Update the handles structure
guidata(hObject, handles);



% --- Executes on button press in sarSave 
%   - Doesn't actually do anything in this version of the Blender
%   - There isn't a sarSave butten in this version so this lack of action is good
%   - Saving happens automatically at exit in this version of the Blender,
%   but future versions could allow useer requested saves by implementing
%   this button and copying the save statement and the corresponding global statement to this function.
function sarSave_Callback(hObject, eventdata, handles)
% hObject    handle to sarSave (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)



% --- Responds to user interaction with the SAR radio buttons that control what
% variable is being edited
% --- Executes on button press in SARbutton.
function SARbutton_Callback(hObject, eventdata, handles)
% hObject    handle to SARbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
%  Toggles the radio buttons
global showSAR showMM5 showSCT
% Set appropriate values of the radio buttons
showMM5 = 0;
showSCT = 0;
showSAR = 1;
set(handles.MM5button,'Value',showMM5);
set(handles.QuikSCATbutton,'Value',showSCT);
%  Redisplay the SAR image with the appropriate overlay colors
showSarFcn(hObject, eventdata, handles);
% Hint: get(hObject,'Value') returns toggle state of SARbutton



% --- Responds to user interaction with the Scatterometer radio buttons that control what
% variable is being edited
% --- Executes on button press in QuikSCATbutton.
function QuikSCATbutton_Callback(hObject, eventdata, handles)
% hObject    handle to QuikSCATbutton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global showSAR showMM5 showSCT
% Set appropriate values of the radio buttons
showMM5 = 0;
showSAR = 0;
showSCT = 1;
set(handles.MM5button,'Value',showMM5);
set(handles.SARbutton,'Value',showSAR);
%  Redisplay the SAR image with the appropriate overlay colors
showSarFcn(hObject, eventdata, handles);
% Hint: get(hObject,'Value') returns toggle state of QuikSCATbutton


% --- Executes on button press in MM5button.
function MM5button_Callback(hObject, eventdata, handles)
% hObject    handle to MM5button (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global showSAR showMM5 showSCT
showSAR = 0;
showSCT = 0;
showMM5 = 1;
set(handles.SARbutton,'Value',showSAR);
set(handles.QuikSCATbutton,'Value',showSCT);
showSarFcn(hObject, eventdata, handles);
% Hint: get(hObject,'Value') returns toggle state of MM5button



% --- If the SAR edit radio button is 'on', hitting the SAR delete button
% deletes the current subjective SAR wind record
%function SARdeleteButton_ButtonDownFcn(hObject, eventdata, handles)
function SARdeleteButton_Callback(hObject, eventdata, handles)
% hObject    handle to SARdeleteButton (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
global irec sarWind sarCount
%  Delete the current record
j = 0;
for i=1:sarCount
    if (i == irec)  % ignore the record to be deleted 
    else % recopy the other records
        j = j + 1;
        sarWind(j,:) = sarWind(i,:);
    end
end
%  Update the bookkeeping
sarCount = j;
sarWind = sarWind(1:sarCount,:);
irec = j;
%  Plot the SAR image and all its overlays using shortened subjective SAR wind record
showSarFcn(hObject, eventdata, handles);
guidata(hObject, handles);



% --- Creates the edit radius slider
% --- Executes during object creation, after setting all properties.
function editRadiusSlider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to editRadiusSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called
% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



% --- Responds to user interaction with the edit radius slider.
% --- Executes on slider movement.
function editRadiusSlider_Callback(hObject, eventdata, handles)
global changeRadius
% hObject    handle to editRadiusSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
%  Get slider value
changeRadius = get(handles.editRadiusSlider, 'value');
%  Set the status text under the slider
numFigures = 2
editRadiusText = strcat('Edit Radius => ',num2str(changeRadius,numFigures));
set(handles.editRadiusText,'string',editRadiusText); 
%  Update the handles structure
guidata(hObject, handles);



% --- Create the east-west slider for SAR image navigagion
% --- Executes during object creation, after setting all properties.
function lonSlider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to lonSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called
% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



% --- Executes on slider movement for east-west navigation of the SAR image
function lonSlider_Callback(hObject, eventdata, handles)
global imageHeight imageWidth sarHeight sarWidth sarLonShift sarLatShift
% hObject    handle to lonSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
%  Get slider value
sarLonShiftRange = sarWidth - imageWidth;
%  Compute how many pixels to offset the SAR and Quikscat data by
sarLonShift = get(handles.lonSlider, 'value');
sarLonShift = floor(sarLonShiftRange * sarLonShift);
%  Replot the wind records
showSarFcn(hObject, eventdata, handles);
%  Update the handles structure
guidata(hObject, handles);




% --- Create the north-south slider for SAR image navigagion
% --- Executes during object creation, after setting all properties.
function latSlider_CreateFcn(hObject, eventdata, handles)
% hObject    handle to latSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    empty - handles not created until after all CreateFcns called
% Hint: slider controls usually have a light gray background, change
%       'usewhitebg' to 0 to use default.  See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
    set(hObject,'BackgroundColor',[.9 .9 .9]);
else
    set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end



% --- Executes on slider movement for north-south SAR image navigation.
function latSlider_Callback(hObject, eventdata, handles)
global imageHeight imageWidth sarHeight sarWidth sarLonShift sarLatShift
% hObject    handle to latSlider (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'Value') returns position of slider
%        get(hObject,'Min') and get(hObject,'Max') to determine range of slider
%  Get slider value
sarLatShiftRange = sarHeight - imageHeight;
%  Compute how many pixels to offset the SAR and Quikscat data by
sarLatShift = get(handles.latSlider, 'value');
sarLatShift = floor(sarLatShiftRange * sarLatShift);
%  Replot the wind records
showSarFcn(hObject, eventdata, handles);
%  Update the handles structure
guidata(hObject, handles);



% --- George Young's bottom-up interpolation function.  A replacement for Matlab's interp1 for non-monotonic data.
function yi = interp1gsy(x,y,xi,method2use)
%  method2use is "n" for nearest neighbor, "l" for linear regression
%  See how many points there are
npts = length(x);
%  Loop through list from start to finish looking for a bracket
yi = 0;
for ipt = 2:npts
    goofbelow = xi-x(ipt-1);
    goofabove = xi-x(ipt);
    %  Test for a bracket
    if goofbelow*goofabove < 0
        if method2use == 'n'
            if abs(goofbelow) < abs(goofabove)
                yi = y(ipt-1);
            else
                yi = y(ipt);
            end
        else  % default to linear interpolation
            yi = y(ipt) + ((xi-x(ipt))*(y(ipt)-y(ipt-1))/(x(ipt)-x(ipt-1)));
        end
        %  Break out of the loop if there is a bracket.
        %  This makes sure you only find the slowest wind speed that brackets
        break
    end
end
%  If the backscatter is higher than "possible" assign the wind for the highest backscatter
[xmax,xmaxi] = max(x);
if (yi == 0) & (xi > xmax)
    yi = y(xmaxi);
end
%  Return from function




%  CMOD5 forward model - JHUAPL - Nathaniel Winstead - July 11, 2007
function sigma = cmod5hyb1c(u10, wdir, incid, r_look, alpha)
%
% CMOD 5 Modified for H-pol (using alpha param following Thompson
% Inputs are u10 (m/s), wind direction (degreed), incid. angle (deg)
% radar look angle (deg), alpha parameter
% result is a sigma_0 value
%
thetm = 40.0;
thethr = 25.0;
zpow = 1.6;
%
% c-coefficients
%
c(1) = -0.688;
c(2) = -0.793;
c(3) = 0.338;
c(4) = -0.173;
c(5) = 0.0;
c(6) = 0.004;
c(7) = 0.111;
c(8) = 0.0162;
c(9) = 6.34;
c(10) = 2.57;
c(11) = -2.18;
c(12) = 0.4;
c(13) = -0.6;
c(14) = 0.045;
c(15) = 0.007;
c(16) = 0.33;
c(17) = 0.012;
c(18) = 22.0;
c(19) = 1.95;
c(20) = 3.0;
c(21) = 8.39;
c(22) = -3.44;
c(23) = 1.36;
c(24) = 5.35;
c(25) = 1.99;
c(26) = 0.29;
c(27) = 3.80;
c(28) = 1.53;
y0 = c(19);
pn = c(20);
a = y0 - (y0 - 1.) / pn;
b = 1. / (pn * (y0 - 1.)^(pn - 1.0));
% Compute the difference between radar look and wind direction
% convert to radians and calculate the cosine
csfi = cos( deg2rad(wdir - r_look));
x = (incid - thetm) / thethr;
a0 = ((c(4)*x + c(3))*x + c(2))*x + c(1);
a1 = c(5) + c(6) * x;
a2 = c(7) + c(8) * x;
gam = (c(11) * x + c(10))*x + c(9);
s1 = c(12) + c(13)*x;
s2 = a2 * u10;
a3 = s2;
w = find(s2 < s1);
if size(w) > 0
  a3(w) = s1; 
end
a3 = 1./(1. + exp(-a3));
if size(w) > 0
  a3(w) = (s2(w)/s1).^(s1*(1.-a3(w))).*a3(w);
end
b0=(a3.^gam).*10.^(a0+a1*u10);
b1=c(15)*u10.*(0.5+x-tanh(4.*(x+c(16)+c(17)*u10)));
b1=(c(14)*(1.+x)-b1)./(exp(0.34.*(u10-c(18)))+1.0);
v0=(c(23)*x+c(22))*x+c(21);
d1=(c(26)*x+c(25))*x+c(24);
d2=c(27)+c(28)*x;
v2=(u10/v0+1.);
w=find(v2 < y0);
if size(w) > 0
  v2(w) = a+b*(v2(w)-1.0).^pn;
end
b2 = (-d1+d2*v2).*exp(-v2);
sigma = b0.*(1.0+b1.*csfi+b2.*(2.0*csfi.*csfi-1.0)).^zpow;
%
% Make into h-pol using hybrid model 
%
%ratio = 1.0;
al = alpha;
tan_theta= tan(deg2rad(incid));
%%tan_theta = tan(pi*theta/180);
tan_theta2 = tan_theta.^2;
ratio = 1 + al * tan_theta2;
ratio = ratio ./ (1.0 + 2.0*tan_theta2);
ratio = ratio.*ratio;
sigma = sigma.*ratio;






%  CMOD4 forward model - JHUAPL - Nathaniel Winstead - July 17, 2007
function sigma = cmod4hyb(wind_speed, wind_dir, theta, r_look, alpha)
x = (theta-40.)./25.;
p0 = 1.0;
p1 = x;
p2 = 0.5*(3.*x.^2-1.);
%
% c-coefficients
%
c(1) = -2.301523;
c(2) = -1.632686;
c(3) =  0.761210;
c(4) =  1.156619;
c(5) =  0.595955;
c(6) = -0.293819;
c(7) = -1.015244;
c(8) =  0.342175;
c(9) = -0.500786;
c(10)=  0.014430;
c(11)=  0.002484;
c(12)=  0.074450;
c(13)=  0.004023;
c(14)=  0.148810;
c(15)=  0.089286;
c(16)= -0.006667;
c(17)=  3.000000;
c(18)= -10.00000;
%
%br-coefficients
%
br(1:3)= 1.075;
br(4)=   1.072;
br(5)=   1.069;
br(6)=   1.066;
br(7)=   1.056;
br(8)=   1.030;
br(9)=   1.004;
br(10)=  0.979;
br(11)=  0.967;
br(12)=  0.958;
br(13)=  0.949;
br(14)=  0.941;
br(15)=  0.934;
br(16)=  0.927;
br(17)=  0.923;
br(18)=  0.930;
br(19)=  0.937;
br(20)=  0.944;
br(21)=  0.955;
br(22)=  0.967;
br(23)=  0.978;
br(24)=  0.988;
br(25)=  0.998;
br(26)=  1.009;
br(27)=  1.021;
br(28)=  1.033;
br(29)=  1.042;
br(30)=  1.050;
br(31)=  1.054;
br(32)=  1.053;
br(33)=  1.052;
br(34)=  1.047;
br(35)=  1.038;
br(36)=  1.028;
br(37)=  1.016;
br(38)=  1.002;
br(39)=  0.989;
br(40)=  0.941;
br(41:44)=0.929;
%th0=(16:60);
tint=(theta-16.0);
nlow=floor(tint)+1;
nhigh=ceil(tint)+1;
if nlow <= 0
    nlow = 1;
end
if nhigh <= 0
    nhigh = 1;
end
b_lookup=(br(nlow)+((br(nhigh)-br(nlow))*(tint-nlow)));
alpha0=  c(1)*p0+c(2)*p1+c(3)*p2;
gamma=   c(4)*p0+c(5)*p1+c(6)*p2;
beta=    c(7)*p0+c(8)*p1+c(9)*p2;
yy=wind_speed+beta;
y=yy;
ff1=find(y>0 & y<5.);
ff2=find(y>5.);
f1=y;
if size(ff1) > 0 
    f1(ff1)=log10(y(ff1));
end
if size(ff2) > 0
    f1(ff2)=sqrt(y(ff2))./3.2;
end
b0 = b_lookup*10.^(alpha0+gamma*f1);
f2 = tanh(2.5*(x+0.35))-0.61*(x+0.35);
b1 = c(10)*p0+c(11)*wind_speed+(c(12)*p0+c(13)*wind_speed)*f2;
b2 = c(14)*p0+c(15)*(1.+p1)*wind_speed;
b3 = 0.42*(1.+c(16)*(c(17)+x)*(c(18)+wind_speed));
%
phi = deg2rad(wind_dir - r_look);

sigma = b0 .* (1. + b1*cos(phi)+b3.*tanh(b2).*cos(2.*phi)).^1.6;
%
% Make into h-pol using hybrid model 
%
%ratio = 1.0;
al = alpha;
tan_theta= tan(deg2rad(theta));
tan_theta2 = tan_theta.^2;
ratio = 1 + al * tan_theta2;
ratio = ratio ./ (1.0 + 2.0*tan_theta2);
ratio = ratio*ratio;
sigma = sigma*ratio;


function xrad = deg2rad(xdeg)
xrad = pi*xdeg/180.0;

