      %  Analyze sieve analysis of grain size to estimate
      %  permeability, porosity, and tortuosity.
      %  Data are analyzed for the median and standard deviation of a
      %  lognormal grain size distribution.
      %  Case 1: (1) millimeter, (2) cumulative wt %
      clear all
      fname=input('File name of input data: ','s'); % extension .dat expected
      clf reset
%     fname2='dum';
      eval(['load ',fname,'.dat'])
      dia=eval([fname ,'(:,2)']);  % diameter is expected in column 2
      wt =eval([fname ,'(:,7)']);  % cmulative % is expected in column 7
      origd=dia;
      origw=wt;
      % Eleminate data with cumulative wt < 1%
      drop=0;
      for i=1:length(wt)
        if wt(i) < 1.0
           drop=i;
        end
      end
      lwt=length(wt);
      wt(1:lwt-drop)=wt(drop+1:lwt);
      dia(1:lwt-drop)=dia(drop+1:lwt);
      if drop > 0
        wt(lwt-drop+1:lwt)=[];
        dia(lwt-drop+1:lwt)=[];
      end
      % Eleminate data with cumulative wt >  99%
      for i=length(wt):-1:1
        if wt(i) > 99.0
           wt(i)=[];
           dia(i)=[];
        end
      end
      %  Transform variables to y=ln(dia), P=wt/100
      y=log(dia);
      P=0.01*wt(:);
      u=sqrt(2)*erfinv(1-2*P);
      % Linear regression
      unit=ones(size(y));
      s=unit'*unit;
      sx=u'*unit;
      sy=y'*unit;
      sxx=u'*u;
      sxy=u'*y;
      del=s*sxx-sx^2;
      mu=(sxx*sy-sx*sxy)/del;
      sig=(s*sxy-sx*sy)/del;
      % Plot linear fit of data
      ul(1)=u(1);
      ul(2)=u(length(u));
      yl(1)=mu+sig*ul(1);
      yl(2)=mu+sig*ul(2);
      plot(u,y,'kx',ul,yl,'k-')
      title('Linear Regression: y = mu + sig u')
      xlabel('u')
      ylabel('y')
      input('Hit enter to continue')
      % Plot grain size distribution and fit
      clf reset
      semilogx(origd,origw,'xk')
      axis([.01 10 0 100])
      xlabel('Grain diameter, mm')
      ylabel('% Coarser by weight')
      hold
      fP=linspace(.001,0.999)';
      fy=zeros(size(fP));
      fy=mu+sig*sqrt(2)*erfinv(1-2*fP);
      fwt=100*fP;
      fdia=zeros(size(fP));
      fdia=exp(fy);
      semilogx(fdia,fwt,'-k')
      % Calculate porosity, tortuosity, and permeability
      mdia=exp(mu);
      por=0.428-0.0998*sig;
      tort=2.464+7.724*sig;
      perm=14070*por^3*mdia^2/((1-por)^2*tort);
      h=text(.015,  25,['mdia = ',num2str(mdia),' mm']);
      h=text(.015,  20,['por  = ',num2str(por)]);
      h=text(.015,  15,['tort = ',num2str(tort)]);
      h=text(.015,  10,['perm = ',num2str(perm),' d']);
      title(fname)
