% Calc for current force on a surface of changing area. 
% 
% This is based on the needs for FOCE "windmill" over a range of current
% 
% copyright (c) 2004 by W.J.Kirkwood of MBARI
% November 23, 2004
% Last Modified 12/2/04
%
% ///////////// ALL CALCS IN Kg - Meter - Sec \\\\\\\\\\\\\\\\\
% 
% formula based on F = 1/2 V**2 rho Cd A
%
% clean up from last activity 
%
clear all ; clf;
%
% set seawater density
%
rho = 1025;  % kg/m3 density of sea water 
%
% Setup of parameters
% ----- Current velocity ranges from 1 cm/sec to 50 cm/sec ------
%
Current = .01:((.5-.01)/49):.5;  % creating an equally spaced 50 element array
%
% ----------- Area estimated as 10 cm2 to .5 m2 -----------------
%
Area = .001:((.5-.001)/49):.5; % creating an equally spaced 50 element array
%
% Adjustment for angle as windmill changes position over 0 to 90 from current
% Both drag and area alter as windmill changes position, changing the force
% Based on estimated shape minimum drag is .04 and minimum area is .001 m2
% 
% Defining the angle mulitplier based on Cosine of 0 to 90 over 50 steps
%
Theta = 0:(pi/98):(pi/2); 
ThetaC = cos(Theta);
%
% Current;  % test statement 
% Area;     % test statement
% ThetaC;   % test statement 
%
% --------- Modified area equaiton becomes A_total = (.001 + (Area*theta)) --- 
%
for i = 1:50
    for j = 1:50
AreaT(i,j) = .001 + (Area(i) .* ThetaC(j));
    end
end
% AreaT;    % test statement 
%
% ------- Create the changing drag based on angle -----------------------
%
% drag of a flat plat is 2.00, faired body is approx .04
% assume chanign drag based on angle (this is inaccurate but gives a
% thumbnail estimation of the drag 
%
Cdrag = .04:((2.0-.04)/49):2.0; % drag range broken into 50 equal parts 
Cdrag; 
% 
% ---------- Add the Drag component of equation in to AreaT based on ------ 
%  ****  NOTE: cos part of math included in area calc as multiplier ****
%  *****   noted that this is not accurate but close enough   *******
%
AreaD = AreaT; % create a working drag area matrices
%
% Set up for loops to modify AreaD
%
for m = 1:50
    for n = 1:50
AreaD(m,n) = Cdrag(n) .* AreaT(m,n);
    end
end
% AreaD;  % test statement
%
% ---- Create flipped matrices for Force Calcs to make plotting easy---
%
AreaI = flipud (AreaD); 
%
% --------- Calc of force as a function of the area ---------------
%
% substituting the values into the formula except for Current
% creating the Force Multiplier Matrices 
for a = 1:50
    for b = 1:50
FMultiple(a,b) = .5 * rho .* AreaI(a,b);
    end
end
% FMultiple; % test statement 
%
% ****  Now the square of current across the elements ****
%       
% Input of current to generate force diagram for given areas 
VelCurrent = input ('Enter current in cm/sec : ');
VelTrans = (VelCurrent/100); 
%
for e = 1:50
    for f = 1:50
Force(e,f) = (VelTrans^2) .* FMultiple(e,f);
    end
end
Force; % test statement 
% testplot
meshgrid([0:50/49:50]);
Z = Force;
plot3(Current, AreaT, Force) 
xlabel('Current cm/sec')
ylabel('Area cm2')
zlabel('Newtons')
grid on