%  Danalysis: perform a dimensional analysis.
%    Find non-dimensional variables via the basis set of the  
%  null space of the dimension matrix, D.  Several examples, 
%  starting with finding the drag on a sphere moving at a known  
%  speed through a viscous fluid, are shown below. The lines 
%  having comments denoted by the symbol %%% are the
%  lines that you will need to change to define a new problem.  
%  Enter the variables by calling Din; enter one dependent 
%  variable first, and then as many independent variables  
%  as appear in the problem.  The dimensions are entered as a 
%  row vector of the exponents on the dimensions Mass, Length, 
%  time and Temperature, i.e, velocity is [0 1 -1 0].  The
%  order of these dimensions is arbitrary, as is the choice 
%  of the fundamental dimensions themselves.  You must be  
%  consistent within any one problem but otherwise the choice
%  is yours. 
%
%    The raw output is a set of exponents, nsb, that are to be 
%  associated with each of the input variables.  The specific 
%  order of the exponents, and thus the specific form of the 
%  non-dimensional variables, is determined by the arbitrary order 
%  in which you entered the independent variables.  Only by 
%  accident would the form be optimal.  Your part of this job is to
%  manipulate the non-d variables into a useful form. For example,
%  in the first problem of drag on a moving sphere, you might 
%  guess that inertial drag is the zero order solution, and hence 
%  you may want to divide Pi(1) by Pi(2) to show this.  Note that
%  Pi is a symbolic variable, and to do this you could 
%  enter >>Pi(1) = Pi(1)/Pi(2), and then >>pretty(simplify(Pi))
%  (assuming that you have the symbolic toolbox). If you do not
%  have the symbolic toolbox (you should get it!), you will still
%  get the exponents on each variable returned through nsb.
% 
%    By Jim Price, July, 2001. jprice@whoi.edu, 508-289-2526. 
%  Please also let me know if you find a problem on which this 
%  method fails or if you have suggestions or comments that
%  would help improve the script.
%  
%   Public domain for all personal, educational purposes. 
%

function Danalysis
%  this has been set up as a function so that all the necessary 
%  subroutines can be included within one file, Din, Dclear,
%  Dnullspace

global D Ds jp nsb Pi
%  if you enter 'global Pi nsb' in the command window before 
%  running Danalysis, the exponent array nsb and the symbolic 
%  variable array Pi will be available in the command window
%  for further use (or at least the last computed values will
%  be available)

format compact
format rational

disp(' ')
disp('   ***  drag on a sphere moving at a known speed  ***')  %%%
Dclear                        %   always start with this
Din('Force',  [1  1 -2] )     %%% drag (force) = [m*l/t*t];  the dependent variable 
Din('rho',    [1 -3  0] )     %%% density of the fluid = [m/l*l*l]
Din('radius', [0  1  0] )     %%% radius of the sphere = [l]
Din('mu',     [1 -1 -1] )     %%% viscosity = [m/l*t]
Din('speed',  [0  1 -1] )     %%% the known speed of the sphere = [l/t]
D;                            %   print out D (or not) 
[nsb Pi] = Dnullspace(Ds, D); %   now perform the analysis on D
% 
%  You may prefer other forms of these Pis, as noted above. 
%   
disp(' hit any key to continue'), disp(' '), disp('  '), pause

%  now try with viscosity omitted
disp('   ***  drag on a sphere with viscosity omitted  ***')  %%%
Dclear                        %   always start with this
Din('Force',  [1  1 -2] )     %%% drag (force) = [m*l/t*t];  the dependent variable 
Din('rho',    [1 -3  0] )     %%% density of the fluid = [m/l*l*l]
Din('radius', [0  1  0] )     %%% radius of the sphere = [l]
%  Din('mu',     [1 -1 -1] )     %%% viscosity = [m/l*t]
Din('speed',      [0  1 -1] ) %%% the known speed of the sphere = [l/t]
D;                            %   print out D (or not) 
[nsb Pi] = Dnullspace(Ds, D); %   now perform the analysis on D
disp(' hit any key to continue'), disp(' '), disp('  '), pause

%  now try with density of the fluid omitted
disp('   ***  drag on a sphere with density omitted  ***')  %%%
Dclear                        %   always start with this
Din('Force',  [1  1 -2] )     %%% drag (force) = [m*l/t*t];  the dependent variable 
% Din('rho',    [1 -3  0] )     %%% density of the fluid = [m/l*l*l]
Din('radius', [0  1  0] )     %%% radius of the sphere = [l]
Din('mu',     [1 -1 -1] )     %%% viscosity = [m/l*t]
Din('speed',  [0  1 -1] )     %%% the known speed of the sphere = [l/t]
D;                            %   print out D (or not) 
[nsb Pi] = Dnullspace(Ds, D); %   now perform the analysis on D
disp(' hit any key to continue'), disp(' '), disp('  '), pause

%  now add g, as if near the sea surface, but w/ viscosity omitted
disp('   ***  drag on a sphere, with g but not viscosity  ***')  %%%
Dclear                        %   always start with this
Din('Force',  [1  1 -2] )     %%% drag (force) = [m*l/t*t];  the dependent variable 
Din('rho',    [1 -3  0] )     %%% density of the fluid = [m/l*l*l]
Din('radius', [0  1  0] )     %%% radius of the sphere = [l]
Din('depth',  [0  1  0] )     %%% depth = [l]
Din('g', [0 1 -2])            %%% accel of gravity
%  Din('mu',     [1 -1 -1] )     %%% viscosity = [m/l*t]
Din('speed',  [0  1 -1] )     %%% the known speed of the sphere = [l/t]
D;                            %   print out D (or not) 
[nsb Pi] = Dnullspace(Ds, D); %   now perform the analysis on D
disp(' hit any key to continue'), disp(' '), disp('  '), pause

%  now add everything back in  
disp('   ***  drag on a moving sphere with all variables  ***')  %%%
Dclear                        %   always start with this
Din('Force',  [1  1 -2] )     %%% drag (force) = [m*l/t*t];  the dependent variable 
Din('rho',    [1 -3  0] )     %%% density of the fluid = [m/l*l*l]
Din('radius', [0  1  0] )     %%% radius of the sphere = [l]
Din('g', [0 1 -2])            %%% accel of gravity
Din('depth',  [0  1  0] )     %%% depth = [l]
Din('mu',     [1 -1 -1] )     %%% viscosity = [m/l*t]
Din('speed',  [0  1 -1] )     %%% the known speed of the sphere = [l/t]
D;                            %   print out D (or not) 
[nsb Pi] = Dnullspace(Ds, D); %   now perform the analysis on D
disp(' hit any key to continue'), disp(' '), disp('  '), pause


disp('   ***  GI Taylor`s famous analysis of the first A bomb  ***  ')
%  grim, but physically very interesting and historical
Dclear                    % reinitialize a counter and memory 
Din('Energy', [1 2 -2]) %%%  energy released (the unknown)
Din('radius', [0 1 0])  %%%  observed radius of the blast wave
Din('time', [0 0 1])    %%%  time elasped since the blast
Din('rho', [1 -3 0])    %%%  density of the ambient fluid
% D;                      % print out D (or not)
[nsb Pi] = Dnullspace(Ds, D);  % perform the analysis on D
%
%  some data read off photos in Sedov (1959), Fig 59-61
%  R is accurate to no better than about 10-20%
t = 1.e-3*[0.1 0.24 0.38 0.80 1.36 1.93 15.0 127.0]; % time, sec
R = 0.5*[24 36 50 68 80 90 220 340];  %  radius of blast wave
E = 8.0e13;   %  (kinetic) energy released, Joules, found by
%   an eyeball fit of the observed and non-d curves.  This
%   fitting would be plausible if the form (slope) of the
%   non-d function were consistent with the data (it seems 
%   to be in this case).  The best fit E corresponds roughly 
%   to 12 kilotons of TNT (perhaps more than you wanted to know?)

rho = 1.25;   %  air density
tmod = 1.e-4*[1:10:1000];
Rmod = (E/rho)^0.2*(tmod.^0.4);

%  some graphics settings
set(0,'DefaultLineLineWidth',1.0)
set(0,'DefaultTextFontSize',14)
set(0,'DefaultAxesLineWidth',1.1)
set(0,'DefaultAxesFontSize',12)
figure(9)
clf reset
loglog(t, R, '*', tmod, Rmod)
legend('observed radius', 'fit of r(t), E = 8x10^{13} J')
xlabel('elapsed time, sec')
ylabel('radius of the blast wave, m')
title('radius(t) of a very intense blast wave')
disp(' hit any key to continue'), disp(' '), disp('  '), pause


disp (' ')
disp('  ***  period of planetary orbits a la Kepler  *** ')  %%%  title
Dclear                  % reinitialize a counter and memory 
Din('T', [0 0 1])       %%%  Period, the unknown 
Din('R', [0 1 0])       %%%  radius of the orbit 
Din('m', [1 0 0])       %%%  mass of the planet
Din('M', [1 0 0])       %%%  mass of the sun  
Din('G', [-1 3 -2])     %%%  universal gravitational constant
%  D;                      % print out D (or not)
[nsb Pi] = Dnullspace(Ds, D);   % perform the nondimensional analysis of D
%
%  Kepler's data
%  planets are: Mercury, Venus, Earth, Mars, Jupiter, Saturn
Rp = [0.389 0.724 1.0 1.524 5.200 9.510];  %  AU 
Tp = [87.8 224.7 365.25 687.0 4332.6 10759.2];  % days 
mp = [0.054 0.815 1.0 0.108 317.8 95.2];  %  mass in units of Earth
Ms = 3.33e5*5.97e24;   %  mass of the sun
%
Rp = Rp*1.49e11;  % convert AU to meters
Tp = Tp*8.64e4;   % days to sec
mp = mp*5.97e24;  % earth masses to kg

Gp = Rp.^3./(Ms*Tp.^2);
format short; 
Gavg = mean(Gp);
G = 6.67e-11;  %  the universal gravitational constant

figure(33); clf reset
loglog(Rp, Tp, Rp, Tp, 'r*')
xlabel('radius, m');  ylabel('period, sec')
Tmod = 2*3.14159*(Rp.^1.5)/sqrt(G*Ms);
hold on
loglog(Rp, Tmod, 'b')
title('orbital period of the inner planets and Jupiter`s moons')

figure(34);  clf reset
Tnd = Tp./((Rp.^1.5)/sqrt(Gavg*Ms));
loglog(mp/5.97e24, Tnd); 
axis([0.01 1000 0.1 10.])
xlabel('mass planet/mass Earth'); ylabel('period, non-d')

%  the moons of Jupiter; Io, Europa, Ganymede, Callisto
%
Rm = [5.57 8.87 14.15 24.90]*69.9e6;   % m
Tm = [1.53 3.07 6.19 14.5]*1.e5;   %  sec
Mj = 317.8*5.97e24;   %  mass of Jupiter in kg
G = 6.67e-11;
figure(33); hold on
loglog(Rm, Tm, 'g*')
Tmod = 2*3.14159*(Rm.^1.5)/sqrt(G*Mj);
hold on
loglog(Rm, Tmod, 'b')
hold on
% loglog(Rm, Tm/sqrt(3.33e5/317.8), 'm')
axis([1.e8 1.e14 1.e3 1.e9])
axis('square')
legend('fit of p(r); G = 6.67 x10^{-11}',4)


function [nsb, Pi] = Dnullspace(Ds, D)
%
%  Dnullspace: non-dimensional analysis of a dimension
%  matrix D.  Ds is a cell array of variable names. The
%  output is a matrix nsb of exponents on the dimensional
%  variables that will give a basis set of the null space 
%  of D. If symbolic toolbox is available, then the symbolic
%  variable Pi is a printable version of the non-d 
%  variables.  The resulting basis set is not unique, and is
%  unlikely to conform to your preferred set of non-d
%  variables.  Jim Price, August, 2001. 

disp(' ')
nrank = rank(D);
fprintf(' the rank of the dimension matrix is %g', nrank)
disp(' ')

[nv nd] = size(D);
numndv = nv - nrank;
fprintf(' the number of non-d variables is N = %g', numndv)
disp(' ')

%  The next two lines are the key steps. First compute the null space
%  basis of the dimensional matrix D by use of null, then put the
%  result into reduced row echelon form by rref.  The latter insures   
%  that the first dimensional variable (the dependent variable) will 
%  appear in one non-d variable only (the first one), and that 
%  it will have an exponent of 1. Beyond that, it is hard to say
%  how nsb (the exponents of the basis set) will be arranged. 

nsb = null(D.','r').';  %  find a basis set of non-d variables.
                        %  null seems to want the transpose of our D.

nsb = rref(nsb);        %  sort nsb so that pi(1) has an exponent 
                        %    =1 on the dependent variable.

%  all that follows is an attempt to make a useful display of 
%    the result computed just above

[nr nc] = size(nsb);

%  list the exponents on the input variables 
disp(' ')
disp(' the dimensional variables and the exponents required to make')
disp('    one possible basis set of non-d variables are:')
disp(' ')
for j = 1:nc
fprintf('    %4s  ', char(Ds(j)))
end
disp(' '), disp(' ')
for jr=1:nr
exponents = [sprintf('%9.4g ', nsb(jr,:))]
end
disp(' ')

%  if the symbolic toolbox is available, then we can 
%    make a more easily interpreted display of the results
%    (if not, you should still get the exponents)

Pi = 0;   % just in case you can't do the following:

%  check to see if the symbolic toolbox is (really) available
if exist('maplemex', 'file') == 3  
   
clear Q QQ Pi
eval('syms Q QQ Pi')

%  store the nc variable names 
for m=1:nc
   QQ(m) = char(Ds(1,m)); 
end

%  assemble the nr non-d variables
for j1=1:nr
  Pi(j1) = prod(QQ.^nsb(j1,:));
end

disp(' one possible set of non-d variables (Pi(1) ... Pi(N)) is:')
pretty(Pi)

end  % if on symbolic toolbox

disp(' ')

function Dclear
%  Dclear: call this once at the start of a non-d 
%  analysis problem to reset some variables.

global D Ds jp nsb
D = [];
Ds = [];
nsb = [];
jp = 0;

function Din(x, y)
%  Din: accept and store input to D and Ds. x is
%  a string with the variable name, y is a row matrix
%  of the exponents on each of the fundamental dimensions
%  that appear in x. The length of y must be the same for
%  each call to Din (within a given problem).

global D Ds jp nsb
jp = jp + 1;

if jp == 1
   Ds = {x};       %  store x in a cell array
else
   Ds = [Ds {x}];  %  add on to Ds with each call 
end

ny = length(y);    %  the number of fundamental dimensions
D(jp, 1:ny) = y(1:ny);  %  store the exponents