function [xar,yar,zz] = dune_plot(h,alpha,flagexpi)

% Plots the shape of a dune of height h.
%
% Michel Louge 7/26/09
% MYL Feb 12, 2014
% Error corrected by Jiannong Fang on line 96, March 6, 2014
%
% Reference:	G. Sauermann, P. Rognon, A. Poliakov, H.J. Herrmann : "The
% shape of the barchan dunes of Southern Morocco", 
% Geomorphology 36, 47 - 62 (2000)  
% available as "sauermann_barcan_shape.pdf"
%
% All distances in meters.
%
%	input	-	h       =   dune height (m)
%           -   alpha   =   inclination of the slip face (rad)
%                               ie, angle of horizontal and slip face
%           -   flagexpi=   flag == 0: same expression as Sauermann, et al.
%                                == 1: corrected with gaussian concavity
%
%	output	-	xar =   array of x-coord for 3d surface plot
%           -   yar =   array of y-coord for 3d surface plot
%           -   zz  =   upwind dune height (m)
%
%   global variables
%   
%           -   al  	=   dune length = al*h + bl (* dune length upwind from origin al = 5.4  bl = 10.7*);
%           -   bl      =   see above
%           -   aw      =   dune width  = aw*h + bw (* dune width aw = 5.8 bw = 1.9 *);
%           -   bw      =   see above
%           -   ab0     =   b0[h_]= ab0*h + bb0 (* brink curve coefficient ab0 = -2.4  bb0 = 7.9 *);
%           -   bb0     =   see above
%           -   bb2     =   b2[h] = bb2*l[h]/w[h]^2 (* brink curve coefficient bb2 = -1.46 *);
%           -   w	=   dune width (m)
%           -   l   =   dune length (m)
%           -   b0  =   brink curve coefficient (m)
%           -   b2  =   brink curve coefficient (m^-1)
%           -   tanalpha =  tangent of slip inclination (-)
%           -   flagexp  =  flagexpi
%           
%   to access the proper directory, type
%   cd ~/Documents/Michel_files/Research/Sahara/dune_life_cycle
%
% Example:
%
% clear all;[x,y,z]=dune_plot(5,pi/4,1);
%

global  al bl aw bw ab0 bb0 bb2 w l b0 b2 tanalpha flagexp

%   internal variables
%
%           -   ym4 =   y-coord of the intersection of brink, ground and windward face (m)
%           -   xm  =   x-coord of that point (m)
%           -   num =   number of points in each direction to define the plotting mesh
%           -   incx=   incremental distance in the x-direction
%           -   incy=   incremental distance in the y-direction

al  = 5.4   ;	% constants adopted by Sauermann, et al.
bl  = 10.7  ; 
aw  = 5.8   ;  
bw  = 1.9   ; 
ab0 = -2.4  ; 
bb0 = 7.9   ;
bb2 = -1.46 ;
tanalpha = tan(alpha) ;
flagexp  =  flagexpi  ;

num =   120 ;	% number of points in each direction to define the plotting mesh
%num =   40 ;	% number of points in each direction to define the plotting mesh

l       =   al*h + bl           ;	% dune length upwind from origin al = 5.4  bl = 10.7
w       =   aw*h + bw           ;	% dune width aw = 5.8  bw = 1.9
b0      =   ab0*h + bb0         ;	% brink curve coefficient ab0 = -2.4  bb0 = 7.9
b2      =   bb2*l/w^2           ;	% brink curve coefficient bb2 = -1.46

ym4	=	sqrt(-((1/(2*bb2^2*(bl + al*h)))*((bw + aw*h)^2*...
		(2*bb0*bb2 + bl + al*h + 2*ab0*bb2*h - ...
		sqrt(bl + al*h)*...
		sqrt(4*bb0*bb2 + bl + 4*bb2^2*bl + (al + 4*ab0*bb2 + ...
        4*al*bb2^2)*h)))))      ;   % y coord of the intersection of brink, ground and windward face
            
xm	=	b0 + b2*ym4^2           ;   % x-coord of brink's intersection with ground and windward face

incx=   abs(l - xm)/(num-1)     ;   % incremental distance in the x-direction
incy=   2*w/(num-1)             ;   % incremental distance in the y-direction

[xar,yar]	=	meshgrid(xm:incx:l,-w:incy:w)	;   % prepare the mesh for 3-D surface plot

for ix = 1:num  % loop along x
    for iy = 1:num   % loop along y

% build dune height array. Jiannong detected error on March 6, 2014:
% zz(ix,iy)	=   height(xar(1,ix),abs(yar(iy,1)),h)	;
        zz(iy,ix)	=   height(xar(1,ix),abs(yar(iy,1)),h)	;
    end
end

scrsz   = get(0,'ScreenSize')                           ;   % finds the screen size of the current machine
figure('Position',...
    [1 scrsz(4) scrsz(3) scrsz(4)/2])                   ;	% reserves space on the screen for the figure

colormap jet
surf(xar,yar,zz+h,'FaceColor','interp',...
 'EdgeColor',[.7 .7 .7],...
 'FaceLighting','phong')	% plot dune surface elevated by the height h
%surfc(xar,yar,zz)	% plot dune surface
%colormap hot
hold on
contour(xar,yar,zz)	% plot dune surface
%axis equal
hold off


% --------------------------------------- %
function [ybrink] = brink_point(xx,yy,h)

% Calculates the y-coord of the brink just above a point at (x,y) on the
% avalanche face.
%
% Michel Louge 7/25/09
%
% Reference:	G. Sauermann, P. Rognon, A. Poliakov, H.J. Herrmann : "The
% shape of the barchan duns of Southern Morocco", 
% Geomorphology 36, 47 - 62 (2000)  sauermann_barcan_shape.pdf
%
%	input	-	xx      =	x-coord pointing upwind (m)
%           -   yy      =   y-coord pointing crosswind (m)
%           -   h       =   dune height (m)
%
%	output	-	ybrink  =   y-coord on brink above pt (x,y) on slip face
%
%   global variables
%   
%           -   see above
%           

global  al bl aw bw ab0 bb0 bb2 w l b0 b2 tanalpha flagexp

%   internal variables
%
%           -   none        =   none

ybrink = real(-(2*6^(1/3)*bb0*bb2^3*bl^3*bw^2 + 6^(1/3)*bb2^2*bl^2*bw^4 + 4*6^(1/3)*aw*bb0*bb2^3*bl^3*bw*h + 6*6^(1/3)*al*bb0*bb2^3*bl^2*bw^2*h + ...
       2*6^(1/3)*ab0*bb2^3*bl^3*bw^2*h + 4*6^(1/3)*aw*bb2^2*bl^2*bw^3*h + 2*6^(1/3)*al*bb2^2*bl*bw^4*h + 2*6^(1/3)*aw^2*bb0*bb2^3*bl^3*h^2 + 12*6^(1/3)*al*aw*bb0*bb2^3*bl^2*bw*h^2 + ...
       4*6^(1/3)*ab0*aw*bb2^3*bl^3*bw*h^2 + 6*6^(1/3)*al^2*bb0*bb2^3*bl*bw^2*h^2 + 6*6^(1/3)*aw^2*bb2^2*bl^2*bw^2*h^2 + 6*6^(1/3)*ab0*al*bb2^3*bl^2*bw^2*h^2 + ...
       8*6^(1/3)*al*aw*bb2^2*bl*bw^3*h^2 + 6^(1/3)*al^2*bb2^2*bw^4*h^2 + 6*6^(1/3)*al*aw^2*bb0*bb2^3*bl^2*h^3 + 2*6^(1/3)*ab0*aw^2*bb2^3*bl^3*h^3 + ...
       12*6^(1/3)*al^2*aw*bb0*bb2^3*bl*bw*h^3 + 4*6^(1/3)*aw^3*bb2^2*bl^2*bw*h^3 + 12*6^(1/3)*ab0*al*aw*bb2^3*bl^2*bw*h^3 + 2*6^(1/3)*al^3*bb0*bb2^3*bw^2*h^3 + ...
       12*6^(1/3)*al*aw^2*bb2^2*bl*bw^2*h^3 + 6*6^(1/3)*ab0*al^2*bb2^3*bl*bw^2*h^3 + 4*6^(1/3)*al^2*aw*bb2^2*bw^3*h^3 + 6*6^(1/3)*al^2*aw^2*bb0*bb2^3*bl*h^4 + ...
       6^(1/3)*aw^4*bb2^2*bl^2*h^4 + 6*6^(1/3)*ab0*al*aw^2*bb2^3*bl^2*h^4 + 4*6^(1/3)*al^3*aw*bb0*bb2^3*bw*h^4 + 8*6^(1/3)*al*aw^3*bb2^2*bl*bw*h^4 + ...
       12*6^(1/3)*ab0*al^2*aw*bb2^3*bl*bw*h^4 + 6*6^(1/3)*al^2*aw^2*bb2^2*bw^2*h^4 + 2*6^(1/3)*ab0*al^3*bb2^3*bw^2*h^4 + 2*6^(1/3)*al^3*aw^2*bb0*bb2^3*h^5 + ...
       2*6^(1/3)*al*aw^4*bb2^2*bl*h^5 + 6*6^(1/3)*ab0*al^2*aw^2*bb2^3*bl*h^5 + 4*6^(1/3)*al^2*aw^3*bb2^2*bw*h^5 + 4*6^(1/3)*ab0*al^3*aw*bb2^3*bw*h^5 + 6^(1/3)*al^2*aw^4*bb2^2*h^6 + ...
       2*6^(1/3)*ab0*al^3*aw^2*bb2^3*h^6 - 2*6^(1/3)*bb2^3*bl^3*bw^2*xx - 4*6^(1/3)*aw*bb2^3*bl^3*bw*h*xx - 6*6^(1/3)*al*bb2^3*bl^2*bw^2*h*xx - 2*6^(1/3)*aw^2*bb2^3*bl^3*h^2*xx - ...
       12*6^(1/3)*al*aw*bb2^3*bl^2*bw*h^2*xx - 6*6^(1/3)*al^2*bb2^3*bl*bw^2*h^2*xx - 6*6^(1/3)*al*aw^2*bb2^3*bl^2*h^3*xx - 12*6^(1/3)*al^2*aw*bb2^3*bl*bw*h^3*xx - ...
       2*6^(1/3)*al^3*bb2^3*bw^2*h^3*xx - 6*6^(1/3)*al^2*aw^2*bb2^3*bl*h^4*xx - 4*6^(1/3)*al^3*aw*bb2^3*bw*h^4*xx - 2*6^(1/3)*al^3*aw^2*bb2^3*h^5*xx - ...
       (9*bb2^4*bl^4*bw^4*yy + 36*aw*bb2^4*bl^4*bw^3*h*yy + 36*al*bb2^4*bl^3*bw^4*h*yy + 54*aw^2*bb2^4*bl^4*bw^2*h^2*yy + 144*al*aw*bb2^4*bl^3*bw^3*h^2*yy + ...
         54*al^2*bb2^4*bl^2*bw^4*h^2*yy + 36*aw^3*bb2^4*bl^4*bw*h^3*yy + 216*al*aw^2*bb2^4*bl^3*bw^2*h^3*yy + 216*al^2*aw*bb2^4*bl^2*bw^3*h^3*yy + 36*al^3*bb2^4*bl*bw^4*h^3*yy + ...
         9*aw^4*bb2^4*bl^4*h^4*yy + 144*al*aw^3*bb2^4*bl^3*bw*h^4*yy + 324*al^2*aw^2*bb2^4*bl^2*bw^2*h^4*yy + 144*al^3*aw*bb2^4*bl*bw^3*h^4*yy + 9*al^4*bb2^4*bw^4*h^4*yy + ...
         36*al*aw^4*bb2^4*bl^3*h^5*yy + 216*al^2*aw^3*bb2^4*bl^2*bw*h^5*yy + 216*al^3*aw^2*bb2^4*bl*bw^2*h^5*yy + 36*al^4*aw*bb2^4*bw^3*h^5*yy + 54*al^2*aw^4*bb2^4*bl^2*h^6*yy + ...
         144*al^3*aw^3*bb2^4*bl*bw*h^6*yy + 54*al^4*aw^2*bb2^4*bw^2*h^6*yy + 36*al^3*aw^4*bb2^4*bl*h^7*yy + 36*al^4*aw^3*bb2^4*bw*h^7*yy + 9*al^4*aw^4*bb2^4*h^8*yy + ...
         sqrt(3)*sqrt(bb2^6*(bl + al*h)^6*(bw + aw*h)^6*(2*(bw^2 + 2*ab0*bb2*bl*h + 2*aw*bw*h + aw^2*h^2 + 2*ab0*al*bb2*h^2 + 2*bb0*bb2*(bl + al*h) - 2*bb2*bl*xx - 2*al*bb2*h*xx)^3 + ...
             27*bb2^2*(bl + al*h)^2*(bw + aw*h)^2*yy^2)))^(2/3))/(6^(2/3)*bb2^2*(bl + al*h)^2*(9*bb2^4*bl^4*bw^4*yy + 36*aw*bb2^4*bl^4*bw^3*h*yy + 36*al*bb2^4*bl^3*bw^4*h*yy + ...
        54*aw^2*bb2^4*bl^4*bw^2*h^2*yy + 144*al*aw*bb2^4*bl^3*bw^3*h^2*yy + 54*al^2*bb2^4*bl^2*bw^4*h^2*yy + 36*aw^3*bb2^4*bl^4*bw*h^3*yy + 216*al*aw^2*bb2^4*bl^3*bw^2*h^3*yy + ...
        216*al^2*aw*bb2^4*bl^2*bw^3*h^3*yy + 36*al^3*bb2^4*bl*bw^4*h^3*yy + 9*aw^4*bb2^4*bl^4*h^4*yy + 144*al*aw^3*bb2^4*bl^3*bw*h^4*yy + 324*al^2*aw^2*bb2^4*bl^2*bw^2*h^4*yy + ...
        144*al^3*aw*bb2^4*bl*bw^3*h^4*yy + 9*al^4*bb2^4*bw^4*h^4*yy + 36*al*aw^4*bb2^4*bl^3*h^5*yy + 216*al^2*aw^3*bb2^4*bl^2*bw*h^5*yy + 216*al^3*aw^2*bb2^4*bl*bw^2*h^5*yy + ...
        36*al^4*aw*bb2^4*bw^3*h^5*yy + 54*al^2*aw^4*bb2^4*bl^2*h^6*yy + 144*al^3*aw^3*bb2^4*bl*bw*h^6*yy + 54*al^4*aw^2*bb2^4*bw^2*h^6*yy + 36*al^3*aw^4*bb2^4*bl*h^7*yy + ...
        36*al^4*aw^3*bb2^4*bw*h^7*yy + 9*al^4*aw^4*bb2^4*h^8*yy + sqrt(3)*sqrt(bb2^6*(bl + al*h)^6*(bw + aw*h)^6*...
           (2*(bw^2 + 2*ab0*bb2*bl*h + 2*aw*bw*h + aw^2*h^2 + 2*ab0*al*bb2*h^2 + 2*bb0*bb2*(bl + al*h) - 2*bb2*bl*xx - 2*al*bb2*h*xx)^3 + 27*bb2^2*(bl + al*h)^2*(bw + aw*h)^2*yy^2)))^(1/3)));

% --------------------------------------- %
function [zz] = height(x,y,h)

% Calculates dune height. (Note y >0).
%
% Michel Louge 7/25/09
%
% Reference:	G. Sauermann, P. Rognon, A. Poliakov, H.J. Herrmann : "The
% shape of the barchan duns of Southern Morocco", 
% Geomorphology 36, 47 - 62 (2000)  sauermann_barcan_shape.pdf
%
%	input	-	x       =	x-coord pointing upwind (m)
%           -   y       =   y-coord pointing crosswind (m)
%           -   h       =   dune height (m)
%
%	output	-	zz      =   dune height upwind (m)
%
%   global variables
%   
%           -   see above
%           
%

global  al bl aw bw ab0 bb0 bb2 w l b0 b2 tanalpha flagexp

%   internal variables
%
%           -   bx      =   x-coord of the brink curve (m)
%           -   bz      =   height of the brink curve at coord y (m)
%           -   derbx   =   derivative of the brink curve at coord y
%           -   bx0     =   x-intersect. of slip face with ground starting from location y at brink
%           -   ybrink  =   y-coord of brink above point at y on slip face
%           -   bxbr    =   x-projection of brink point above y on the ground
%           -   bzbr    =   height of the brink curve above coord y (m)
%           -   derbr   =   derivative of the brink curve at coord ybrink
%           -   bxbr0   =   x-intersect. of slip face with ground starting from ybrink

bx      =   b0 + b2*y^2                         ;	% x-projection of the brink curve on the ground
bz      =   h*(1 - (bx/l)^2 - (y/w)^2)          ;	% height of the brink curve at coord y (m)
derbx	=	(2*bb2*(bl + al*h)*y)/(bw + aw*h)^2	;   % derivative of the brink curve at coord y
bx0     =	bx - bz/tanalpha/sqrt(1 + derbx^2)  ;   % x-intersect. of slip face with ground starting from location y at brink


    if x < bx   % in case x belongs to the slip face

        ybrink	=   brink_point(x,y,h)                  ;   % y-coord of brink above point at y on slip face
        bxbr    =   b0 + b2*ybrink^2                    ;	% x-projection of the brink curve on the ground
        bzbr    =   h*(1 - (bxbr/l)^2 - (ybrink/w)^2)	;	% height of the brink curve at ybrink
        derbr	=	(2*bb2*(bl + al*h)*ybrink)/(bw + aw*h)^2	;   % derivative of brink curve at ybrink
        bxbr0   =	bxbr - bzbr/tanalpha/sqrt(1 + derbr^2)      ;	% x-intersect. of slip face with ground starting from ybrink

zz      =   max(0,bzbr*((x - bxbr0)/(bxbr - bxbr0)))    ;   % dune height on slip face (m)

    else

zz      =	max(0,h*(1 - (x/l)^2 - (y/w)^2))            ;  % upwind dune height

    end
    
    if flagexp == 1 && x > 0
        zz =  zz*exp(-((x - max(0,b0))/l)^2)  ; % gaussian smoothing of toe
    end

% --------------------------------------- %
