function stretch(filename)
% STRETCH(PARAM_FILENAME)
%    Reads ASCII parameter file PARAM_FILENAME, loads a CSP file (DUN
%    model), finds the contours of all dunes, calcultes the first and
%    second order moments for stretching analysis and saves the values.

%clear all
close all

%read_params_from_file=0;

% more on
% echo on


%%%%%%%% Parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%if read_params_from_file
if exist('filename')
    fprintf('reading %s ...\n',filename)
    %fid=fopen('PARAMS.txt','r');
    fid=fopen(filename,'r');
    C=textscan(fid,'%s');
    for i=1:size(C{1})
        eval([cell2mat(C{1}(i))]);
    end
    fclose(fid);
    disp('done')
else
    datafile='DUN0018.csp';
    h_lim=1;      %height threshold for countour detection
    size_min=10;  %minimal size of a dune (length and width)
    flag_find_dunes=1; %finds the exact contours of the dunes (0|1)
    show_figs=1;  %show figures (0|1)
    save_figs=0;  %save figures (0|1)
    save_data=0;  %save data (0|1)
end


%%%%%%%% Read CSP file %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ddcsp=datafile;
disp('loading ...')
[a,h,lx,ly] = read_csp(ddcsp);
disp('done')

%a=reshape(csp,lx,h,ly);


%%%%%%%% Topography %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%b=(a==0 | a==6);
b=(a==0);
alt=zeros(lx,ly); %memory allocation
dd=datafile(1:length(datafile)-4);

for k=1:ly
    for i=1:lx
        alt(i,k)=h-min([h find(b(k,:,i))]);
    end
end

altb=(alt>h_lim);
%min(min(alt))
%max(max(alt))

% altb=[altb altb
%     altb altb];
% 
% alt=[alt alt
%     alt alt];
% 
% lx=2*lx;

if show_figs
    disp('display topography ...')
    figure
    imagesc(alt)
    set(gca,'plotboxaspectratio',[1 1 1]);
end

if save_figs
    disp('saving figure 1 ...')
    %ddfig=strcat('TE',d,'_alt.fig');
    %saveas(gcf,ddfig);
    ddjpg=strcat(dd,'_alt.jpg');
    print('-djpeg','-r100',ddjpg);
    disp('done')
end

%figure
%imagesc(altb)


%%%%%%%% Contours of patterns %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if show_figs
    disp('display contours ...')
    figure
    [c,ho]=contour(altb,1);
    colormap(gray);
    set(gca,'plotboxaspectratio',[1 1 1]);
else
    c=contourc(altb,1);    
end

%c(:,find(c(1,:)==0.5))


%figure

if ~size(c,1)
    disp('no patterns !')
    close all;
    return
end

ind=(find(c(1,:)==0.5));

%%% number of patterns
ll=length(ind);
fprintf('number of patterns: %d\n',ll);


%%%%%%%% Detection of dunes, characteristic points and moments of inertia  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

jk=0;
st=zeros(ll,30);
Mst = zeros(3,3);
        
for i=1:ll

    %%%%% bounding box %%%%%

    i1 = ind(i)+1;
    if i<ll
        i2 = ind(i+1)-1;
    else
        i2 = size(c,2);
    end
    
    minlar = min(c(1,i1:i2));
    minlon = min(c(2,i1:i2));
    maxlar = max(c(1,i1:i2));
    maxlon = max(c(2,i1:i2));
    
    %flag = minlar~=1 && minlon~=1 && maxlar~=lx && maxlon~=ly && (maxlar - minlar >= 5) && (maxlon - minlon >= 5);
    flag1 = (maxlar - minlar >= size_min) && (maxlon - minlon >= size_min);
    flag2 = minlar~=1 && minlon~=1 && maxlar~=lx && maxlon~=ly;
    if flag1 && flag2
        if show_figs
            hold on;
            if ~flag_find_dunes
                fill([minlar maxlar maxlar minlar], [minlon minlon maxlon maxlon],'r');
            else
                fill([minlar maxlar maxlar minlar], [minlon minlon maxlon maxlon],'w');
            end
        end
    
        %%%%% summit %%%%%
        minhau=h_lim-1;
        maxhau=0;
        jk=jk+1;
        for ii=floor(minlar)+1:floor(maxlar)
            for ij=floor(minlon)+1:floor(maxlon)
                aa=alt(ij,ii);
                if aa >= maxhau
                    loxhau=ii;
                    loyhau=ij;
                end
                maxhau=max(maxhau,aa);
            end
        end
        
        %%%%% localization of the dune %%%%%
        if flag_find_dunes
            dune=find_dune(altb(floor(minlon):floor(maxlon),floor(minlar):floor(maxlar)), loxhau-floor(minlar)+1, loyhau-floor(minlon)+1);
            if show_figs
                hold on;
                imagesc([floor(minlar) floor(maxlar)],[floor(minlon) floor(maxlon)],dune);
                colormap([[1;1] [1;0.2] [1;0]]);
                if exist('alpha') % check for octave compatibility
                    alpha(0.5);
                end
            end
        end
        
        %%%%% volume and moments of inertia %%%%%
        vol=0.0;
        mix=0.0;
        miy=0.0;
        miz=0.0;
        for ii=floor(minlar):floor(maxlar)
            for ij=floor(minlon):floor(maxlon)
                x=ii-floor(minlar);
                y=ij-floor(minlon);
                if ~flag_find_dunes || dune(y+1,x+1)
                    aa=alt(ij,ii);
                    mix=mix+aa*x;
                    miy=miy+aa*y;
                    miz=miz+aa*(aa+1)/2;
                    vol=vol+aa;
                end
            end
        end
        mix=mix/vol;
        miy=miy/vol;
        miz=miz/vol;

        mixx=0.0;
        miyy=0.0;
        mizz=0.0;
        mixy=0.0;
        mixz=0.0;
        miyz=0.0;
        for ii=floor(minlar):floor(maxlar)
            for ij=floor(minlon):floor(maxlon)
                x=ii-floor(minlar);
                y=ij-floor(minlon);
                if ~flag_find_dunes || dune(y+1,x+1)
                    aa=alt(ij,ii);
                    mixx=mixx+aa*(x-mix)^2;
                    miyy=miyy+aa*(y-miy)^2;
                    mizz=mizz+aa*(aa+1)*(2*aa+1)/6; % sum of the first n squares = n(n+1)(2n+1)/6
                    mixy=mixy+aa*(x-mix)*(y-miy);
                    mixz=mixz+(x-mix)*aa*(aa+1)/2;
                    miyz=miyz+(y-miy)*aa*(aa+1)/2;
                end
            end
        end
        mixx=mixx/vol;
        miyy=miyy/vol;
        mizz=(mizz/vol)-(miz^2);
        mixy=mixy/vol;
        mixz=(mixz/vol);%-mix*miz;
        miyz=(miyz/vol);%-miy*miz;

        %%%%% central hole %%%%%
        % approximation : we go straight to the east from the summit
        crx=loxhau;
        cry=loyhau;
        while crx<lx && alt(cry,crx+1)>h_lim
            crx=crx+1;
        end

        
        %%%%% horns %%%%%
        %cx1 = min(c(1,find(c(1,i1:i2)<loxhau)))
        %cy1 = min(c(2,find(c(1,i1:i2)<loxhau)))
        %c2 = minlon;
        cx1 = 0;
        cy1 = 0;
        cx2 = 0;
        cy2 = 0;
        for ii=i1:i2
            if c(2,ii) < loyhau && c(1,ii) > cx1
                cx1 = c(1,ii);
                cy1 = c(2,ii);
            end
            if c(2,ii) > loyhau && c(1,ii) > cx2
                cx2 = c(1,ii);
                cy2 = c(2,ii);
            end
        end
        
        if show_figs
            hold on;
            plot(cx1,cy1,'o');
            hold on;
            plot(cx2,cy2,'o');
        end

        %%%%% stretching matrix %%%%%
        Mst(1,1) = mixx;
        Mst(1,2) = mixy;
        Mst(1,3) = mixz;
        Mst(2,1) = mixy;
        Mst(2,2) = miyy;
        Mst(2,3) = miyz;
        Mst(3,1) = mixz;
        Mst(3,2) = miyz;
        Mst(3,3) = mizz;
        
        %%% dimensions
        st(jk,1)=maxlar-minlar;
        st(jk,2)=maxlon-minlon;
        st(jk,3)=maxhau-minhau;
        %%% various values
        st(jk,4)=i2-i1+1;
        st(jk,5)=minlar;
        st(jk,6)=maxlar;
        st(jk,7)=minlon;
        st(jk,8)=maxlon;
        st(jk,9)=minhau;
        st(jk,10)=maxhau;
        %%% summit
        st(jk,11)=loxhau;
        st(jk,12)=loyhau;
        %%% volume
        st(jk,13)=vol;
        %%% central hole
        st(jk,14)=crx;
        st(jk,15)=cry;
        %%% flag 'dune truncated' (by boundaries)
        st(jk,16)=flag2;
        %%% ends of the horns
        st(jk,17)=cx1;
        st(jk,18)=cy1;
        st(jk,19)=cx2;
        st(jk,20)=cy2;
        %%% moments of inertia: first order (center of mass)
        st(jk,21)=mix;
        st(jk,22)=miy;
        st(jk,23)=miz;
        %%% moments of inertia: second order (stretching)
        st(jk,24)=mixx;
        st(jk,25)=miyy;
        st(jk,26)=mizz;
        st(jk,27)=mixy;
        st(jk,28)=mixz;
        st(jk,29)=miyz;
    
    end
end

fprintf('number of dunes: %d\n',jk);



%%%%%%%% Save data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if save_data
    disp('saving data ...')
    ddmat=strcat(dd,'_stat.mat');
    save (ddmat,'st');
    ddtxt=strcat(dd,'_mi.txt');
    mi=st(:,21:29);
    save (ddtxt,'mi','-ASCII');    
    disp('done')
end


if show_figs
    hold on;
    plot(st(:,11), st(:,12), '+');
    hold on;
    plot(st(:,14), st(:,15), 'x');
    xlim([1 lx]);
    ylim([1 ly]);
    set(gca,'plotboxaspectratio',[1 1 1]);
end

if save_figs
    disp('saving figure 2 ...')
    %ddfig=strcat(dd,'_contour.fig');
    %saveas(gcf,ddfig);
    ddpng=strcat(dd,'_contour.png');
    print('-dpng','-r100',ddpng)
    disp('done')
end

return

