clc
clear all
close all

set(groot, 'defaultAxesTickLabelInterpreter','latex');
set(groot,'defaulttextinterpreter','latex');
set(groot, 'defaultLegendInterpreter','latex');


an=365.25*24*3600;
hour=24;

fprintf('Analysis Bilma .... \n')

law=2;                  % Transport law (1) 42 classic; (2) Kok; (3) Pahtz; (else) classic - threshold et 25
dif=0;

fig=0;     		% Plotting test: 0 figure all, 1 figure flux, else no figure.
fsave=2;                % Printing test: 0 print, else no print.

% Name of the figures
iy=1;
nom1=strcat('wind_rose_',num2str(iy),'.eps');
nom2=strcat('flux_rose_',num2str(iy),'.eps');
nom3=strcat('flux_rose_',num2str(iy),'b.eps');
nom4=strcat('mode_rose_',num2str(iy),'.eps');
nom5=strcat('results_rose_',num2str(iy),'.txt');

%%%%%%%%%%%%%%%%% Loading the wind data
%a1=load('wd_+085.7_+37.8.txt');
%a1=load('wd_+085.8_+37.9.txt');
%a1=load('wd_+085.6_+37.9.txt');
%a1=load('wd_+085.6_+38.0.txt');
% The reference point
a1=load('wd_+085.7_+37.9.txt');
%%%%%%%%%%%%%%%%% Selecting a time period

deb_all=datenum(1981,1,1,1,0,0);
fin_all=datenum(2020,12,31,23,0,0);

debt=datenum(2010,1,1,0,0,0);
fint=datenum(2020,12,31,23,0,0);

indeb=round((debt-deb_all)*24+1);
infin=round((fint-deb_all)*24+1);

a=a1(indeb:infin,:);
%a=a1;

the=a(:,1);
u=a(:,2);
%in=find(the>135 & the<270);
%u(in)=0.0;
%in=find(u>2 | u<1);
%u(in)=0.0;
dt=(ones(length(a),1)/length(a))*an;

mu=tand(35);		% bed firction coefficient 
kappa = 0.4;            % Von Karman constant
z_0 = 1e-5;             % Roughness of the sediment layer [m]
g = 9.81;               % Gravitational acceleration [m/s^2]

d= 160e-6;               % Grain size [m] (from Lu, 2009)
%d= 332.5e-6;               % Grain size [m] (from Lu, 2009)
%d= 180e-6;               % Grain size [m] (from Lu, 2009)
%d= 400e-6;               % Grain size [m] (from Lu, 2009)

rhoair = 1.05; %1.293;          % Air density [kg/m^3]
rhosed = 2.55e3;        % sediment density [kg/m^3]
z = 10.0;                 % Height of wind speed measurment [m]

%phi = mod(-90-phi,360);
%Shear velocity
u_cis = u*kappa/log(z/z_0);

%Threshold shear velocity
switch law
  case {1,2,3}
	u_c = 0.082*sqrt((rhosed/rhoair)*g*d);
	%u_c = 0.140*sqrt((rhosed/rhoair)*g*d);
  otherwise
	u_c = 0.1*sqrt((rhosed/rhoair)*g*d);
end

iin=find(u_cis>=u_c);
m_u_cis=mean(u_cis(iin));
rap_u_cis=m_u_cis/u_c;
rap_t=length(iin)/length(u_cis);

%% FLUX
switch law
  case 1
        q=max(0,42*(rhoair/rhosed)*sqrt(d/g)*(u_cis.^2-u_c.^2));
  case 2
        q=max(0,(5/0.6)*(rhoair/rhosed)*(1/g)*u_c*(u_cis.^2-u_c.^2));
  case 3
        q=max(0,(2/(kappa*mu))*(rhoair/rhosed)*(1/g)*u_c*(u_cis.^2-u_c^2).*(1+(rhoair/rhosed)*(1/(g*d))*(1.7/mu)*(u_cis.^2-u_c^2)));
  otherwise
        q=max(0,25*(rhoair/rhosed)*sqrt(d/g)*(u_cis.^2-u_c.^2));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


qdt=q.*dt;                % Here we calculate the flux vector for each time step.

temps_tot=sum(dt);

ref=[0:1:360];
ll=length(ref);

fluxx=qdt.*cosd(the);     % Components of the individual flux vectors on a flat sand bed
fluxy=qdt.*sind(the);

dp=sum(qdt)/temps_tot;    % DP = Drift Potential
rdp=sqrt(sum(fluxx)^2+sum(fluxy)^2)/temps_tot; % RDP = Resultant Drif Potential.
ano=mod(atan2(sum(fluxy),sum(fluxx))*180/pi,360);  % RDD = Resultant Drift Direction
rdp_dp=rdp/dp;  % RDP/DP ratio

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%  Figures - Roses
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

h1=figure('PaperType','A4','PaperUnits','normalized',...
'PaperPosition',[0 0 1 1],'PaperOrientation','portrait');


hw=all_wind_rose_clem(the,u,'cmap',[parula],'n',36,'legtype',2,'di',[0:1:10]);

set(gca,'color','none')

if fsave==0
   eval(['print -depsc ' nom1])
 end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%  Figures - Histogramme des flux
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

b=[];

for ii=1:ll-1
	if ii==1
 		fa=find(round(the)==0 | round(the)==360);
 	else
		fa=find(round(the)==ref(ii));
 	end
 if 1-isempty(fa)
   qa(ii,1)=ref(ii);
   qa(ii,2)=length(fa);
   qa(ii,3)=sum(qdt(fa));
        if qa(ii,2)~=0
            b=[b;ones(floor(1e5*qa(ii,3)./sum(qdt)),1)*qa(ii,1)];
        end
%   qa(ii,1)=ref(ii);
%   qa(ii,2)=sum(dt(fa));
%    if qa(ii,2)~=0
%       b=[b;ones(floor(10000*qa(ii,2)./temps_tot),1)*qa(ii,1)];
%    end
 end
end

if fig<=1
 h2=figure('PaperType','A4','PaperUnits','normalized',...
        'PaperPosition',[0 0 1 1],'PaperOrientation','portrait');

 hf=all_flux_rose_clem(b,ones(length(b),1),'n',36,'legtype',0);

 set(gca,'color','none')

 if fsave==0
   eval(['print -depsc ' nom2])
 end

 h2b=figure('PaperType','A4','PaperUnits','normalized',...
	'PaperPosition',[0 0 1 1],'PaperOrientation','portrait');

 pas=1;
 dec=0;
 [n,x]=hist(b,[0:pas:359]);

 subplot(211)
 plot(mod(x-dec,360),n/(pas*sum(n)),'r.')
 xlabel('Direction','fontsize',24);
 ylabel('PDF','fontsize',24);
 xlabelstr={'East';'North';'West';'South'};
 set(gca,'XLim',[0 360],'plotboxaspectratio',[2.5 1 1],'fontsize',20,'Xtick',[0 90 180 270 360],'Xticklabel',xlabelstr,'fontname','latex','color','none');

 subplot(212)
 plot(mod(x-dec,360),cumsum(n/(pas*sum(n))),'r.')
 xlabel('Direction','fontsize',24);
 ylabel('PDF','fontsize',24);
 xlabelstr={'East';'North';'West';'South'};
 set(gca,'XLim',[0 360],'plotboxaspectratio',[2.5 1 1],'fontsize',20,'Xtick',[0 90 180 270 360],'Xticklabel',xlabelstr,'fontname','latex','color','none');

end

if fsave==0
  eval(['print -depsc ' nom3])
end

%ngauss=3;
%dec=160;
%[mu_save,sig_save,w_save]=em_analyse_new(b,ngauss,dec)

%pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Choice of the modes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%figrot=figure('PaperType','A4','PaperUnits','normalized',...
%        'PaperPosition',[0 0 1 1],'PaperOrientation','portrait');

figpath=figure('PaperType','A4','PaperUnits','normalized',...
        'PaperPosition',[0 0 1 1],'PaperOrientation','portrait');
minqx=0;
maxqx=0;
minqy=0;
maxqy=0;

figtr=figure('PaperType','A4','PaperUnits','normalized',...
        'PaperPosition',[0 0 1 1],'PaperOrientation','portrait');
%plot([deb_all:1/24:fin_all],rdp*[1:length(a)]'/length(a)*an-(cumsum(fluxx).*cosd(ano)+cumsum(fluxy).*sind(ano)))
plot([debt:1/24:fint],rdp*[1:length(a)]'/length(a)*an-(cumsum(fluxx).*cosd(ano)+cumsum(fluxy).*sind(ano)))
datetick('x','yy')
xlabel('Time','FontName', 'latex', 'FontSize', 20)
ylabel('Transport Deficit $[{\rm m}^{2}]$','FontName', 'latex', 'FontSize', 20)
%title('$N_{\rm C_1},\; <t_{To}$','FontName', 'latex', 'FontSize', 20)
set(gca,'FontSize', 20, 'FontName','latex','linewidth',1.5,'plotboxaspectratio',[2 1 1])


prec=0.1;

k=0;

for gamma=[0 0.5];   % We explore different speed-up ratio values.

  k=k+1;            

%  ang(k)=transverse_instability_gamma_q_latex(the,qdt,prec,gamma,fig,dif,law);               % Estimation of crest orientation in the Bed Instability mode
  ang(k)=transverse_instability_gamma_u_latex(the,u_cis,u_c,dt,g,d,rhoair,rhosed,prec,gamma,fig,law);
%  if fig==0
%     pause
%     close all
%  end

 
  ucrest=u_cis.*(1+gamma*abs(sind(the-ang(k))));                                     %% shear velocity at the crest according to its orientation "alpha" [m/s]

  figure('PaperType','A4','PaperUnits','normalized',...
        'PaperPosition',[0 0 1 1],'PaperOrientation','portrait');

 all_wind_rose_clem(the,ucrest*log(z/z_0)/kappa,'cmap',[parula],'n',36,'legtype',2,'di',[0:1:10]);

switch law
  case 1
	qcrest=max(0,42*(rhoair/rhosed)*sqrt(d/g)*(ucrest.^2-u_c.^2).*dt);
  case 2
  	qcrest=max(0,(5/0.6)*(rhoair/rhosed)*(1/g)*u_c*(ucrest.^2-u_c.^2).*dt);
  case 3
	qcrest=max(0,(2/(kappa*mu))*(rhoair/rhosed)*(1/g)*u_c*(ucrest.^2-u_c^2).*(1+(rhoair/rhosed)*(1/(g*d))*(1.7/mu)*(ucrest.^2-u_c^2)).*dt);
  otherwise
	qcrest=max(0,25*(rhoair/rhosed)*sqrt(d/g)*(ucrest.^2-u_c.^2).*dt);
end

if gamma~=0

                b=[];

                                for ii=1:ll-1
                                        if ii==1
                                                fa=find(round(the)==0 | round(the)==360);
                                        else
                                                fa=find(round(the)==ref(ii));
                                        end
                                   if 1-isempty(fa)
                                     qa(ii,1)=ref(ii);
                                     qa(ii,2)=length(fa);
                                     qa(ii,3)=sum(qcrest(fa));
                                        if qa(ii,2)~=0
                                            b=[b;ones(floor(1e5*qa(ii,3)./sum(qcrest)),1)*qa(ii,1)];
                                        end
                                   end
                                end

                if fig<=1

                        hang=figure('PaperType','A4','PaperUnits','normalized',...
                              'PaperPosition',[0 0 1 1],'PaperOrientation','portrait');

                        hfang=all_flux_rose_clem(b,ones(length(b),1),'n',36,'legtype',0);

                        set(gca,'color','none')

			h2ang=figure('PaperType','A4','PaperUnits','normalized',...
        			     'PaperPosition',[0 0 1 1],'PaperOrientation','portrait');

 			pas=1;
 			dec=0;
 			[n,x]=hist(b,[0:pas:359]);
 			plot(mod(x-dec,360),n/(pas*sum(n)),'r.')
 			xlabel('Direction','fontsize',24);
 			ylabel('PDF','fontsize',24);
 			xlabelstr={'East';'North';'West';'South'};
 			set(gca,'XLim',[0 360],'plotboxaspectratio',[2.5 1 1],'fontsize',20,'Xtick',[0 90 180 270 360],...
			        'Xticklabel',xlabelstr,'fontname','latex','color','none');

                end

end

  figure(figpath)
  plot(cumsum(qcrest.*cosd(the)),cumsum(qcrest.*sind(the)))

  minqx=min(minqx,min(cumsum(qcrest.*cosd(the))));
  maxqx=max(maxqx,max(cumsum(qcrest.*cosd(the))));
  minqy=min(minqy,min(cumsum(qcrest.*sind(the))));
  maxqy=max(maxqy,max(cumsum(qcrest.*sind(the))));

  hold on

  angbi=atan2(sum(qcrest.*sind(the)),sum(qcrest.*cosd(the)))*180/pi;
  rotqcrest=[cumsum(qcrest.*cosd(the)) cumsum(qcrest.*sind(the))]*[cosd(angbi) -sind(angbi);sind(angbi) cosd(angbi)];

%  figure(figrot)
%  plot(rotqcrest(:,1)/rotqcrest(end,1),rotqcrest(:,2))
%  hold on

  Qy=sum(qcrest.*sind(the));
  Qx=sum(qcrest.*cosd(the));

  Qi(k)=sqrt(Qx^2+Qy^2)/temps_tot;      % Mean flux at the crest of dunes in the Bed Instability mode.
  aQi(k)=(atan2(Qy,Qx)*180/pi);         % Orientation of the mean flux at the crest of dunes in the Bed Instability mode.
                                        % This is the direction of the propagation of dunes in the Bed Instability mode.

  sig_i(k)=sum(qcrest.*abs(sind(the-ang(k))))/(temps_tot);  % Dune height growth rate of dunes in the Bed Instability mode.

  Qi_para(k)=sum(qcrest.*cosd(ang(k)-the))/(temps_tot);
  Qi_perp(k)=sum(qcrest.*sind(ang(k)-the))/(temps_tot);

%  fprintf('\n Transverse mode done .......\n');

%pause
  ang1(k)=fingering_instability_gamma_u_latex(the,u_cis,u_c,dt,g,d,rhoair,rhosed,prec,gamma,fig,law);      %  Estimation of crest orientation in the fingering mode
%  if fig==0
%    pause
%    close all
%  end

  ucrest=u_cis.*(1+gamma*abs(sind(the-ang1(k))));                                     %% shear velocity at the crest according to its orientation "alpha" [m/s]

 figure('PaperType','A4','PaperUnits','normalized',...
        'PaperPosition',[0 0 1 1],'PaperOrientation','portrait');

 all_wind_rose_clem(the,ucrest*log(z/z_0)/kappa,'cmap',[parula],'n',36,'legtype',2,'di',[0:1:10]);

switch law
  case 1
        qcrest=max(0,42*(rhoair/rhosed)*sqrt(d/g)*(ucrest.^2-u_c.^2).*dt);
  case 2
        qcrest=max(0,(5/0.6)*(rhoair/rhosed)*(1/g)*u_c*(ucrest.^2-u_c.^2).*dt);
  case 3
        qcrest=max(0,(2/(kappa*mu))*(rhoair/rhosed)*(1/g)*u_c*(ucrest.^2-u_c^2).*(1+(rhoair/rhosed)*(1/(g*d))*(1.7/mu)*(ucrest.^2-u_c^2)).*dt);
  otherwise
        qcrest=max(0,25*(rhoair/rhosed)*sqrt(d/g)*(ucrest.^2-u_c.^2).*dt);
end

  rotqcrest=[cumsum(qcrest.*cosd(the)) cumsum(qcrest.*sind(the))]*[cosd(ang1(k)) -sind(ang1(k));sind(ang1(k)) cosd(ang1(k))];

  if gamma~=0

  figure(figpath)
  plot(cumsum(qcrest.*cosd(the)),cumsum(qcrest.*sind(the)),'linewidth',1.5)
 
  minqx=min(minqx,min(cumsum(qcrest.*cosd(the))));
  maxqx=max(maxqx,max(cumsum(qcrest.*cosd(the))));
  minqy=min(minqy,min(cumsum(qcrest.*sind(the))));
  maxqy=max(maxqy,max(cumsum(qcrest.*sind(the))));

%  angbi=atan2(sum(qcrest.*sind(the)),sum(qcrest.*cosd(the)))*180/pi;
%  rotqcrest=[cumsum(qcrest.*cosd(the)) cumsum(qcrest.*sind(the))]*[cosd(angbi) -sind(angbi);sind(angbi) cosd(angbi)];
%  figure(figrot)
%  plot(rotqcrest(:,1)/rotqcrest(end,1),rotqcrest(:,2))
  
		b=[];

				for ii=1:ll-1
				        if ii==1
				                fa=find(round(the)==0 | round(the)==360);
				        else
				                fa=find(round(the)==ref(ii));
				        end
				   if 1-isempty(fa)
				     qa(ii,1)=ref(ii);
				     qa(ii,2)=length(fa);
				     qa(ii,3)=sum(qcrest(fa));
				        if qa(ii,2)~=0
				            b=[b;ones(floor(1e5*qa(ii,3)./sum(qcrest)),1)*qa(ii,1)];
				        end
				   end
				end

		if fig<=1
			hang1=figure('PaperType','A4','PaperUnits','normalized',...
        		      'PaperPosition',[0 0 1 1],'PaperOrientation','portrait');
			hfang1=all_flux_rose_clem(b,ones(length(b),1),'n',36,'legtype',0);
			set(gca,'color','none')

			h2ang1=figure('PaperType','A4','PaperUnits','normalized',...
                                     'PaperPosition',[0 0 1 1],'PaperOrientation','portrait');
                        pas=1;
                        dec=0;
                        [n,x]=hist(b,[0:pas:359]);
                        plot(mod(x-dec,360),n/(pas*sum(n)),'r.')
                        xlabel('Direction','fontsize',24);
                        ylabel('PDF','fontsize',24);
                        xlabelstr={'East';'North';'West';'South'};
                        set(gca,'XLim',[0 360],'plotboxaspectratio',[2.5 1 1],'fontsize',20,'Xtick',[0 90 180 270 360],...
                                'Xticklabel',xlabelstr,'fontname','latex','color','none');

		end

  end

  Qy=sum(qcrest.*sind(the));
  Qx=sum(qcrest.*cosd(the));

  Qf(k)=sqrt(Qx^2+Qy^2)/temps_tot;  % Mean flux at the crest of dunes in the fingering mode.
  aQf(k)=atan2(Qy,Qx)*180/pi;           % Orientation of the mean flux at the crest of dunes in the Fingering mode.
                                        % By definition, it is the direction of the propagation of dunes in the Fingering mode,
                                        % except if the sand source is mobile.

  sig_f(k)=sum(qcrest.*abs(sind(the-ang1(k))))/(temps_tot); %  Dune height growth rate of dunes in the Fingering mode.

%  fprintf('\n Fingering linear mode done .......\n');


end

figure(figpath)
xlabel('$Q_x$','fontsize',24,'fontname','latex');
ylabel('$Q_y$','fontsize',24,'fontname','latex');
axis([minqx maxqx minqy maxqy])
set(gca,'plotboxaspectratio',[maxqx-minqx maxqy-minqy 1],'fontsize',20,'fontname','latex','color','none');


    qtr1=[0 0 cosd(aQi(1)) sind(aQi(1))];
    qtr2=[0 0 cosd(aQi(2)) sind(aQi(2))];

    tl=[cosd(ang(1)) sind(ang(1)) -cosd(ang(1)) -sind(ang(1))];

    tnl=[cosd(ang(2)) sind(ang(2)) -cosd(ang(2)) -sind(ang(2))];

    fl=[-cosd(ang1(1)) -sind(ang1(1)) cosd(ang1(1)) sind(ang1(1))];

    fnl=[-cosd(ang1(2)) -sind(ang1(2)) cosd(ang1(2)) sind(ang1(2))];

    fl0=[ 0 0 cosd(ang1(1)) sind(ang1(1))];

    fnl0=[ 0 0 cosd(ang1(2)) sind(ang1(2))];

    fo=[-cosd(ano) -sind(ano) cosd(ano) sind(ano)];

    gg=[cos([0:0.01:2*pi])' sin([0:0.01:2*pi])'];

    eps=1-1.5e-1;
    headw=0.06;
    headh=0.08;

if fig<=1
    h3=figure('PaperType','A4','PaperUnits','normalized',...
        'PaperPosition',[0 0 1 1],'PaperOrientation','portrait');


    plot(gg(:,1),gg(:,2),'b--','linewidth',2)
    hold on
    plot_arrow(tl(1),tl(2),tl(3),tl(4),'Color','r','linewidth',2,'linestyle','-','facecolor','r','edgecolor','r','headwidth',headw,'headheight',headh);
    plot_arrow(tl(3),tl(4),tl(1),tl(2),'Color','r','linewidth',2,'linestyle','-','facecolor','r','edgecolor','r','headwidth',headw,'headheight',headh);
    plot_arrow(tnl(1),tnl(2),tnl(3),tnl(4),'Color','r','linewidth',2,'linestyle','--','facecolor','w','edgecolor','r','headwidth',headw,'headheight',headh);
    plot_arrow(eps*tnl(1),eps*tnl(2),tnl(1),tnl(2),'Color','w','linewidth',2,'linestyle','--','facecolor','w','edgecolor','r','headwidth',1,'headheight',1);
    plot_arrow(fl0(1),fl0(2),fl(3),fl(4),'Color','k','linewidth',2,'linestyle','-','facecolor','k','edgecolor','k','headwidth',headw,'headheight',headh);
    plot_arrow(fnl0(1),fnl0(2),fnl(3),fnl(4),'Color','k','linewidth',2,'linestyle','--','facecolor','w','edgecolor','k','headwidth',headw,'headheight',headh);
    plot_arrow(qtr2(1),qtr2(2),qtr2(3),qtr2(4),'Color','g','linewidth',2,'linestyle','-','facecolor','g','edgecolor','g','headwidth',headw,'headheight',headh);
    set(gca,'visible','off','plotboxaspectratio',[1 1 1],'color','none')
    if fsave==0
      eval(['print -depsc ' nom4])
    end
end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  if fsave==0
    eval(['fid=fopen(''' nom5 ''',''w'');'])

    fprintf(fid,'\n------------------------------------------------------------------- \n');
    fprintf(fid,'-------------------- Wind/flux propertie --------------------------- \n');
    fprintf(fid,'-------------------------------------------------------------------- \n');
    fprintf(fid,'\n flux total (Drift Potential)   =   %12.6g   m2/yr \n',dp);
    fprintf(fid,'\n Relative drift potential       =   %12.6g   m2/yr \n',rdp);
    fprintf(fid,'\n RDP/DP                         =   %12.6g \n',rdp_dp);

    fprintf(fid,'\n------------------------------------------------------------------- \n');
    fprintf(fid,'-------------------- Potential crest orientations ------------------ \n');
    fprintf(fid,'-------------------------------------------------------------------- \n');
    fprintf(fid,'\n Linear bed instability    =   %8.4g   degree clockwise from east \n',ang(1));
    fprintf(fid,'\n Nonlinear bed instability =   %8.4g   degree clockwise from east \n',ang(2));
    fprintf(fid,'\n Linear fingering          =   %8.4g   degree clockwise from east \n',ang1(1));
    fprintf(fid,'\n Nonlinear fingering       =   %8.4g   degree clockwise from east \n\n',ang1(2));
  end


fprintf('Done ... \n',i);

fprintf('\n     d = %10.6f       zs = %8.4f       z = %8.4f\n',d,z_0,z)
%fprintf('\n')

fprintf('\n     uc = %8.4f \n',u_c)
%fprintf('\n')

fprintf('\n Law = %i ---------------  delta = %8.4f\n',law,gamma)
%fprintf('\n')

fprintf('\n     <u*> = %8.4f \n',m_u_cis)
%fprintf('\n')

fprintf('\n     <u*>/uc = %8.4f \n',rap_u_cis)
%fprintf('\n')

fprintf('\n     Proportion temps = %8.4f \n',rap_t)
%fprintf('\n')

fprintf('\n     DP = %8.4f \n',dp*an)
%fprintf('\n')

fprintf('\n     RDP = %8.4f \n',rdp*an)
%fprintf('\n')

fprintf('\n     RDP/DP = %8.4f \n',rdp_dp)
%fprintf('\n')

fprintf('\n     RDD = %8.4f \n',ano)
%fprintf('\n')

fprintf('\n     alpha_I = %8.4g degree clockwise from east \n',ang(2))
%fprintf('\n')

fprintf('\n     alpha_F = %8.4g degree clockwise from east \n',ang1(2))
%fprintf('\n')

fprintf('\n     Q_I = %8.4f \n',Qi(2)*an)
%fprintf('\n')

fprintf('\n     Q_F = %8.4f \n',Qf(2)*an)
%fprintf('\n')

fprintf('\n     Q_M = %8.4f \n',Qi(2)*abs(sind(aQi(2)-ang(2))*an))
%fprintf('\n')

fprintf('\n     Direction Q_I = %8.4f \n',aQi(2))
%fprintf('\n')

fprintf('\n     Sigma_F/Sigma_I = %6.4f \n',sig_f(2)/sig_i(2))
%fprintf('\n')

fprintf('\n     Q_F/Q_I = %6.4f \n',Qf(2)/Qi(2))
%fprintf('\n')

fprintf('\n     Q_F/Q_M = %6.4f \n',Qf(2)/(Qi(2)*abs(sind(aQi(2)-ang(2)))))
%fprintf('\n')

fprintf('\n     Direction of Q_I with respect to Q_F = %6.4f \n',aQi(2)-aQf(2))
%fprintf('\n')

fprintf('\n     Qi_para/Qi_perp =%6.4f \n', Qi_para(2)/Qi_perp(2))

fprintf('\n     Q_M/Sigma_I = %6.4f \n', (Qi(2)*abs(sind(aQi(2)-ang(2))))/sig_i(2))

fprintf('\n     Q_F/Sigma_F = %6.4f \n', Qf(2)/sig_f(2))

%pause(3)

%for i=1:18
%	figure(i)
%	nom1=sprintf('figure_%02d',i);
%	eval(['print -depsc ' nom1]);
%end

