%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% PURPOSE: Compute buoyancy and CB/CG displacement
% AUTHOR:  Rob McEwen
% DATE:    01/6/1
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Units and Fundamental Constants:
  units
  PbDensity    = 11340;        % Density of solid lead; kg/m^3
  dia          = 21*in2m;      % Maximum diameter of vehicle; m
  rhoK         = 1025;         % Density of sea water; kg/m^3; Knut's calc.
  vel          = 1.5433;       % 3 kts
  Zdelta       = 533.4159;     %.5*rho*(LinCtrlArea*dCl + RingArea*dClR);
  deltaE       = 5*d2r;        % Max allowable elevator trim
%
% Zdelta = Ydelta defined in my notes, because of symmetry.
%
% Load the measurements:
%
  wbData 
%
% Find the distance to the lift points:
%
  lf =  len/2 - 31*in2m;  %Measured by hand; dist to fwd ring.
  la = -len/2 + 52*in2m;  %Measured by hand; dist to aft ring.
%
% Start screen log:
  if exist([caseName,'.txt']) == 2,
    eval(['!rm run', caseName,'.txt'])
    disp([' Removed ', caseName,'.txt; starting another.'])
  end
  diary([caseName,'.txt'])
%
% Compute conversion from Pb mass, in kg, to weight-in-harbor-water, Newtons:
%
  PbKg2ww = (1 - rho/PbDensity)*g;    %See notes, Vol 4, p.31.
%
% Compute the forward and aft weight in water:
%
  Wf = Mf*PbKg2ww;
  Wa = Ma*PbKg2ww;
%  
% Residual buoyant force
%
  W = Wa + Wf;
%
% The total buoyancy of the vehicle is:
%
  B = mass*g + W;
%
% Compute W1 and W2 as defined on p29, Vol 4.  W1 is defined with Wf to
% STARBOARD, and Wa to PORT.  W2 assumes both are to STARBOARD.
%
  W1 = (Mf - Ma)*PbKg2ww;
  W2 = (Mf + Ma)*PbKg2ww;
%
% Compute c.m. in the z-y plane.  See notes, Vol 4, p.29
% 
  Acm = [ B*cos(phi1)    W1-B*sin(phi1);
          B*cos(phi2)    W2-B*sin(phi2)  ];
%
  Bcm = -[ W1  W2 ]'*dia/2;  
%
  cg_23 = Acm\Bcm;
%
  yGc = cg_23(1);
  zGc = cg_23(2);
%
% From notes, p. 33
% 
  xGc = (Wf*lf - Wa*la  +  ( (Wf + Wa)*(zGc + dia/2) - B*zGc )*tan(theta))/...
        (B + Wf + Wa);
%
  fprintf('\n')
  fprintf('\n')
  fprintf('\n')
  fprintf('  ****STATIC MEASUREMENT TEST REPORT****\n')
  fprintf('\n')
  fprintf('\n')
  fprintf('  Vehicle:        %s \n',vehicle)
  fprintf('  Date:           %s \n',testdate)
  fprintf('  Participants:   %s \n',crew)
  fprintf('  Comments:         \n')
  fprintf(comments)
  fprintf('\n')
  fprintf('\n')
  fprintf('  The density of the water in the test area was %8.1f kg/m^3.\n',...
	     rho);
  fprintf('\n')
  fprintf('  The weight-in-water added to the forward lifting ring \n')
  fprintf('  was %8.3f N (%8.3f lbf) \n',Wf, Wf/lbf2n)
  fprintf('\n')
  fprintf('  The weight-in-water added to the after lifting ring \n')
  fprintf('  was %8.3f N (%8.3f lbf) \n',Wa, Wa/lbf2n)
  fprintf('\n')
  fprintf('  The vehicle had a positive buoyant force of\n')
  fprintf('  %8.3f N (%8.3f lbf) \n',W, W/lbf2n)
  fprintf('\n')
%
% Compute the limits of GC travel given deltaE = 5 Deg.
%
%  xMaxG = (len/2 + xG)*Zdelta*deltaE/(B - Zdelta*deltaE);
%  fprintf(['  An x axis c.g. offset of %8.3f mm will result in an \n',...
%	   '  elevator trim angle of 5 degrees. \n'], xMaxG*1000);
  fprintf('\n')
% fprintf('  The x axis component of the c.g. was %8.3f mm. \n',xGc*1000);
  fprintf('  Including the entrained water;\n');
  fprintf('  The y axis component of the c.g. was %8.3f mm. \n',yGc*1000);
  fprintf('  The z axis component of the c.g. was %8.3f mm. \n',zGc*1000);
  
  diary off