function dsd = dsd_fun(inpars, dstep, nbins, kind)

% Calculates the Particle Size distribution at 'nbins' bins wide 'dstep [mm]'
% for a 'kind' of scatterer and a parametrization given by 'inpar'
% kind  scatterer
% 3     rain
% 37    rain
% 4     ice
% 5     snow
% 6     graupel
% 
%
% inpar(1)
% 
%
%
%

inpars = inpars * 1.0;
dstep = dstep * 1.0;

dmm = dstep*(1:nbins);

switch inpars(1)
case 0 % Exponential with free parameters
   mode = 0;
   nnot = inpars(2);
   lam = inpars(3);
case 1 % Gamma with free parameters
   mode = 1;
   nnot = inpars(2);
   lam = inpars(3);
   mu = inpars(4);
case 2 % 4 parameter Gamma with free parameters
   mode = 2;
   nnot = inpars(2);
   lam = inpars(3);
   mu = inpars(4);
   eps = inpars(5);
case 3 %
   mode = 0;
   disp('Exponential with fixed lambda not ready')
case 4 % Marshall Palmer : Ref. Doviak, p.213
   mode = 0;
   nnot = 8.0e3;
   lam = 4.1*(inpars(2)^(-0.21));
case 5 % R,D",S" - Ziad Haddad 1997 - TOGA/COARE dataset
   mode = 1;
   R = inpars(2);
   Ddp = inpars(3);
   S = inpars(4);
   if R < 0.001 
      mu = 1.0;
      lam = 1.0;
      nnot = 0.0;
   else
      mu = (1 / ((S^2)*(Ddp^0.33)*(R^0.074))) - 4;
      lam = 1 / ((S^2)*(Ddp^1.33)*(R^0.23));
      nnot = 55 * R * ( (lam^(mu+4)) / (gamma(mu+4) * (1- (1+0.53/lam)^(-mu-4) ) ) );
   end
case 6 % R,D',S' - Ziad Haddad 1996 - Darwin dataset
   mode = 1;
   M = inpars(2);
   D = inpars(3);
   S = inpars(4);
   % disp('Darwin parametrization Not ready')
   mu = 1/S^2 - 4;
   lam = 1 / ((S^2)* D);
   nnot = M / (pi/6 * gamma(1/S^2) * (D*S^2)^(1/S^2) * 1e-3);
case 7 % Meagher 2001 (?)
   mode = 3;
   disp('Joss-Meagher parametrization Not ready')
case 9  % Exponential with svetla's N0 and lam(M)
   mode = 0;
   switch kind 
   case 7
      nnot = 8.0e3;
      lam = 4.1*(inpars(2)^(-0.21)); % Rain Rate
   case 6 % Graupel
      nnot = 10.0e3;
      lam = (pi * 0.6e6 * nnot * 1.0e-9 / inpars(2)) ^ 0.25;
   case 5 % Snow
      nnot = 3.0e3;
      lam = (pi * 0.1e6 * nnot * 1.0e-9 / inpars(2)) ^ 0.25;
   case 4 % Cloud Ice, uniform 0-0.165mm
      nnot = inpars(2) * 24 / (pi * (0.92 * 1.0e-3) * max(d)^4); 
      lam = 0;
   case 37 % Rain
      nnot = 22.0e3;
      lam = (pi * 1.0e6 * nnot * 1.0e-9 / inpars(2)) ^ 0.25;
   case 2 % Cloud Water, uniform 0-0.33mm
      nnot = inpars(2) * 24 / (pi * (1.0e-3) * max(d)^4);   
      lam = 0;
   end
case 10 % Exponential with free N0 and lam(M) or lam(R)
   mode = 0
   nnot = inpars(3)
   switch kind 
   case 7, lam = 0.4*(inpars(2)^(-0.23))*(nnot*10)^0.21;	% Rain Self Made with Vt 9.65[1-exp(-5.3D)]
   case 6, lam = (pi * 0.6e6 * nnot * 1.0e-9 / inpars(2)) ^ 0.25; % Graupel
   case 5, lam = (pi * 0.1e6 * nnot * 1.0e-9 / inpars(2)) ^ 0.25; % Snow
   case 4, lam = 0; % Cloud Ice
   case 37, lam = (pi * 1.0e6 * nnot * 1.0e-9 / inpars(2)) ^ 0.25; % Rain
   case 2, lam = 0; % Cloud Water
   end
case 11 % Exponential with free lam and N0(M)
   mode = 0;
   lam = inpars(3);
   switch kind
   case 7, nnot = 8000;	% DUMMY - Rain Rate
   case 6, nnot = lam^4 * inpars(2) / (pi * 0.6e6 * 1.0e-9 ); % Graupel
   case 5, nnot = lam^4 * inpars(2) / (pi * 0.1e6 * 1.0e-9 ); % Snow
   case 4, nnot = 0; % Cloud Ice
   case 37, nnot = lam^4 * inpars(2) / (pi * 1.0e6 * 1.0e-9 ); % Rain
   case 2, nnot = 0; % Cloud Water
   end
otherwise, disp('Not available')
end

switch mode
case 0, dsd = nnot.*exp(-lam*dmm) * dstep;
case 1, dsd = nnot.*(dmm.^mu).*exp(-lam*dmm) * dstep;
case 2, dsd = nnot.*(dmm.^mu).*exp(-lam*dmm.^ooooooooooooes) * dstep;
case 3, disp('Not ready')
otherwise, disp('Not Available')
end ~isempty(h3)