function render_dune(filename,flag_tt)
% RENDER_DUNE(FILENAME,FLAG_TT)
%    Loads a CSP file with DUN model, extracts and smooths the topography, and display it with phong shading.
%    GRJ cells are rendered in red color.
%    If FLAG_TT is set, displays a turn table and removes sand outside of it.

close all
%clear all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% cellular space
%filename='Regen_DUN0900.csp';
%filename='dun55_DUN1000.csp';
[a,H,L,P] = read_csp(filename);

%b=(a==0 | a==6);% | a==4);
b=(a==0);
c=(a==1);

disp('topography')
for k=1:P
  for i=1:L
    alt(i,k)=H-1-min([H find(b(i,:,k))]);
  end
end

figure('PaperUnits','normalized', 'PaperPosition',[ 0 0 1 1 ],'color','white')
imagesc(alt)
daspect([1 1 1])
nom1=strcat('surf_2d','.png');
eval(['print -dpng ' nom1]);

%return

figure('PaperUnits','normalized', 'PaperPosition',[ 0 0 1 1 ],'color','white')
surf(alt)
shading interp
daspect([1 1 H/L*2])
nom1=strcat('surf_3d','.png');
eval(['print -dpng ' nom1]);

figure
%figure('PaperUnits','normalized', 'PaperPosition',[ 0 0 1 1 ],'color','white')
%subplot(2,2,3)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% turn table
if exist('flag_tt','var') && (flag_tt ~= 0)
  disp('turn table')
  tt=zeros(L,H,P);
  %% center
  cx=L/2;
  cy=P/2;
  %% radius
  r=L/2-2;
  %% height
  htab=1;

  for k=1:P
    for i=1:L
      if sqrt((i-cx)^2+(k-cy)^2)< r
        tt(i,H+1-htab:H,k)=1;
      end
    end
  end

  %% remove all cells outside the turn table
  for k=1:P
    for i=1:L
      if sqrt((i-cx)^2+(k-cy)^2)>= r
        b(i,:,k)=0;
        c(i,:,k)=0;
      end
    end
  end

  col_tt = [0.25 0.25 0.25];
  %col_tt = [0 0 0];
  data_tt = smooth3(tt,'box',1);
  p_tt = patch(isosurface(data_tt,0.15), 'FaceColor', col_tt, 'EdgeColor', 'none');
  isonormals(data_tt,p_tt)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% GR cells
disp('GR cells')
col = [.7 .4 0];
%col = [.8 .6 0];
data = smooth3(b,'box',3);
p = patch(isosurface(data,0.30), 'FaceColor', col, 'EdgeColor', 'none');
%q = patch(isocaps(data,0.30), 'FaceColor', 'interp', 'EdgeColor', 'none');
isonormals(data,p)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% GRJ cells
disp('GRJ cells')
col2 = [1 .2 .2];
%col2 = [0 .2 1];
data2 = smooth3(c,'box',1);
%p2 = patch(isosurface(data2,0.35), 'FaceColor', col2, 'EdgeColor', 'none');
p2 = patch(isosurface(data2,0.15), 'FaceColor', col2, 'EdgeColor', 'none');
%q2 = patch(isocaps(data2,0.35), 'FaceColor', 'interp', 'EdgeColor', 'none');
isonormals(data2,p2)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
axis([0 L 0 H 0 L])
%daspect([1 3*H/L 3])
daspect([1 1 1])
%view(3)
axis vis3d tight

%view(-112,-42)
%camlight(-112,-42)
%camlight(-12,-42)
%camlight(32,-42)
%camlight(32,42)
% %view(-67,-58)
% view(-135,-65)
% camlight(-75,-75)
% camlight(-75,-5)
% camlight(-5,-5)
% camlight(55,-5)

% view(90,0)
% camroll(-135)
% camorbit(0,-45)
% camproj('perspective')
% 
% camlight(75,75)
lighting phong
set(gca,'visible','off','box','off')

nom1=strcat('surf_smooth','.png');
eval(['print -dpng ' nom1]);

view(-90,0)
camroll(180)
camorbit(0,-45)
zoom(1.75)

camlight(-112,-42)
camlight(-12,-42)
camlight(32,-42)
camlight(32,42)

%% strong light towards east direction
%camlight(-90,45)
%camlight(-90,68)

%view(230,15)

%bonne vue couloir
%view(-121,-35)


% %%%%%% Vue plein est au ras du sol avec le parallelepipede couche sur le cote Nord.
% view(180,0)
% %print -djpeg99 fig_isole_180_0.jpg
% print -depsc fig_isole_180_0.eps
% 
% %%%%%% Vue plein ouest au ras du sol avec le parallelepipede couche sur le cote Sud.
% view(180,180)
% %print -djpeg99 fig_isole_180_180.jpg
% print -depsc fig_isole_180_180.eps
% 
% %%%%%% Vue de haut a l'est avec le parallelepipede couche sur le cote Nord.
 %view(230,15)
% %print -djpeg99 fig_isole_230_15.jpg
% print -depsc fig_isole_230_15.eps
% 
% %%%%%% Vue symmetrique
% view(320,15)
% %print -djpeg99 fig_isole_320_15.jpg
% print -depsc fig_isole_320_15.eps
% 

return

