function cmodtest12
%  cmodtest - a simple driver program calling to Matlab CMOD4 or CMOD5
%  functions for use as a portable SAR wind speed retrieval system


%  Housekeeping
clear all
close all
clc



%  Global variable definitions - used to communicate with some of the
%  functions *** double-check the neccesity of each of these ***
global sardir sarrcs sartheta speed sarspeed sarlookdir
global fileName speed ni nj
global cmodversion


%  <--------------------  USER INPUT REQUIRED HERE  ---------------------->
%  Set the program options
%  Specify whether to use CMOD4 or CMOD5: possible values = 4 or 5
cmodversion = 5;
%  Set test case: possible values = 1, 2, or 3
testcase = 3;
%  Set decimation factor - how many pixels to step between samples - 1
%  uses all pixels, 2 uses every other pixel, etc.
%  This parametercan be tuned as needed to look at more data or achieve faster runs
decfac = 1;
%  Change this to the relative or absolute path of the folder in which
%  you've stored your data files. 
filePath = 'Test_Cases/';
%  <-----------------------  END OF USER INPUT  -------------------------->

%   Set the file name - this can be changed as needed to read any file
casedates = ['0303UTC01Oct2004';'0326UTC01Jan2007';'0538UTC30APR2007']
fileName = [casedates(testcase,:),'_cmod',num2str(cmodversion)]


%   Load real data
inputFile = strcat(filePath, fileName, '.mat');
if (exist(inputFile) == 2)  % If the JHU/APL SAR file exists, load it
    display 'Loading an existing input data file'
    load(inputFile)
else % Note that there are no JHU/APL SAR data file so the program can't work
    display 'Input data file does not exist'
    return
end
%  See what variables have been loaded from the data file
who


%  Sub-sample the SAR variables to cut memory consumption
[isize,jsize] = size(sardir);
sarlat = sarlat(1:decfac:isize);
sarlon = sarlon(1:decfac:isize);
sardir = sardir(1:decfac:isize,1:decfac:jsize);
sarspeed = sarspeed(1:decfac:isize,1:decfac:jsize);
sarrcs = sarrcs(1:decfac:isize,1:decfac:jsize);
sartheta = sartheta(1:decfac:isize,1:decfac:jsize);
[iblock,jblock] = size(sarspeed);


% %   Display the initial fields
disp('About to start first figure')
fimage = figure

subplot(2,3,1)
size(sarrcs);
imagesc(sarrcs)
title('backscatter')
xlabel('X')
ylabel('Y')

subplot(2,3,4)
imagesc(-log(sarrcs))
title('-log of backscatter')
xlabel('X')
ylabel('Y')

subplot(2,3,2)
image(sarspeed)
title('Input wind speed')
xlabel('X')
ylabel('Y')

subplot(2,3,5)
imagesc(sardir)
title('input wind directon')
xlabel('X')
ylabel('Y')


%   Call CMOD - everything goes in through global variables
CMODcall;


%   Do thresholding checks on the data so that statistics and plots include
%   only data for which the algorithms are valid (theta 20 to 40 degrees and
%   wind speed 2 to 35 m/s).
isvalid = find((sartheta >= 20) & (sartheta <=40) & (sarspeed >= 2) & (sarspeed <= 35));


%   Display output fields
subplot(2,3,3)
image(speed)
xlabel('X')
ylabel('Y')
title('Output wind speed')


%  Display input vs output speeds
fscat = figure;
plot(sarspeed(isvalid),speed(isvalid),'.b',[0,35],[0,35],'-k')
title(['Output vs Input wind speed using CMOD',num2str(cmodversion),' on ',fileName])
xlabel('Input wind speed')
ylabel('Output wind speed')
saveas(fscat,['scatterplot_for_',fileName,'.png'])


%  Compute the comparison statistics statistics
clc
disp(fileName)
mean_diference = mean(sarspeed(isvalid)-speed(isvalid))
std_difference = std(sarspeed(isvalid)-speed(isvalid))
r_squared = corrcoef(sarspeed(isvalid),speed(isvalid));
r_squared = r_squared(1,2)
n = length(isvalid)


%  Save the new wind results to to a different file
sarspeed = speed;
save(['unthresholded_output_for_',fileName],'sardir','sarlat','sarlon','sarlookdir','sarrcs','sarspeed','sartheta')
%  Zero out the wind speed for out-of-bounds pixels and do the same
notvalid = find((sartheta < 20) | (sartheta > 40) | (sarspeed < 2) | (sarspeed > 35));
sarspeed(notvalid) = 0;
save(['thresholded_output_for_',fileName],'sardir','sarlat','sarlon','sarlookdir','sarrcs','sarspeed','sartheta')


end




function CMODcall
%   The first level of CMOD calling code - does a bit of housekeeping and
%   then calls rhe rcs to wind conversion function.  That function does a
%   bit more housekeeping and then calls CMOD.

global sardir sarlookdir ni nj nimage w sarrcs sartheta speed sarspeed
% global sardir sarlookdir ni nj nimage w sarrcs sartheta sigmavec theta alpha rlook wind_dir wind_speed speed
global fileName speed ni nj


%  Prepare filename for output
outputFile = strcat('cmod_results_', fileName, '.mat')


%   <------------- Start of CMOD call
%  Locate the pixels with non-zero rcs (i.e. those that are actually in the swath)
[ni,nj] = find(sarrcs > 0.0);
disp('ni and nj sizes')
size(ni)
size(nj)
w = find(sarrcs > 0.0);
nimage=size(w);


%  Call the first of a linked sequence of two functions that make up CMOD
rcs2wind;   
%   <------------- End of real CMOD call

end




% ****************** Begin second level of cmod calling code
function rcs2wind
%   The second level of CMOD calling code - does a
%   bit more housekeeping and then calls CMOD.


global sardir sarlookdir ni nj nimage w sarrcs sartheta sigmavec theta alpha rlook wind_dir wind_speed speed
global cmodversion


%  Set the maximum wind to be retrieved
maxwind = 35.0;
%  Create a list of winds to be retrieved - linear interpolation of wind
%  speed given rcs will done between values spaced every m/s from 0 to
%  maxwind 
wind_speed = (0:1.0:maxwind); 
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);


%  Start timer
tic
%  Loop over all of the available pixels, calling CMOD once per pixel
for j = 1:nimage(1)
    %  Call either CMOD4 or CMOD5 - make sure to use the right test data
    %  file for the algotithm!
    if cmodversion == 4
        sigmavec = cmod4hyb(wind_speed, sardir(ni(j),nj(j)), sartheta(ni(j),nj(j)), sarlookdir, alpha);
    elseif cmodversion == 5
        sigmavec = cmod5hyb1c(wind_speed, sardir(ni(j),nj(j)), sartheta(ni(j),nj(j)), sarlookdir, alpha);
    else
        disp('Illegal CMOD version specified');
        crash
    end
    %  Store the sigma vector
    sigma(:,ni(j),nj(j)) = sigmavec;
end
% save('sigmafilenew.mat','sigma')
cmodtime = toc;


%  Use linear interpolation to look up the right wind in the sima table.
%  Start timer
tic;
%  Pre-dimension the matrix
speed = zeros(size(sartheta));
for j = 1:nimage(1)
    temp = sigma(:,ni(j),nj(j));
    utemp = wind_speed(:);
    xi = sarrcs(ni(j),nj(j));
    %  Because Matlab's interpolation function can fail, particularly on
    %  non-monotonic functions, we use our own one-dimensional linear
    %  interpolation function
    speed(ni(j),nj(j)) = interp1gsy(temp,utemp,xi,'l');
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.  There are other ways to approach this
% problem, but this is the approach used at JHUAPL and CSU-Meteorology
% *************************** Pete, you need to read this and respond to it  ************************************************
ww=find(sarrcs > max(max(max(sigma))));  %gsy *** - This implementationseems to lose the dependence of max sigma on theta.  
%  Is this a problem?  It seems to me that there may be too many maxes here.  I supect we need to do this pixel-by-pixel.
%  In my interpolation, I've handled this problem by by assigning the wind at max rcs rather than the max wind.  Which should it be?
speed(ww) = maxwind;
end




% --- 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
end



%  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;

end




%  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;
end

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