I = input("Body moment of inertia? in kgm^2");
l = input("Length of body? in m ");
b = input("Breadth of body? in m ");
w = input("Height of body? in m ");
d = input("cross-section diameter of towed array? in m ");
h = input("length of towed array? in m ");
s = input("span of hydrofoil? in m ");
c = input("chord of hydrofoil? in m ");
yf = input("y-coordinate of application of hydrofoil force? in m ");
xf = input("x-coordinate of application of hydrofoil force? in m ");
yc = input("y-coordinate of application of towed array force? in m ");
xc = input("x-coordinate of application of towed array force? in m ");
yt = input("y-coordinate of application of towing force? in m");
xt = input("x-coordinate of application of towing force? in m");
yb = input("y-coordinate of application of body force? in m");
xb = input("x-coordinate of application of body force? in m");
xcm = input("center of gravity x-coordinate?");
ycm = input("center of gravity y-coordinate?");
xcb = input("center of buoyancy x-coordinate?");
ycb = input("center of buoyancy y-coordinate?");
beta = input("Hydrofoil angle of attack? in rad ");
k = input("Initial tilt? in m/s ");
Cf = input("Skin friction coefficient");
Cdp = input("profile drag coefficient");
Clc = Cf*sind(k)+Cdp*(sind(k)^2);
Cd = input("Drag coefficient body? ");
Cl = input("Lift coefficient body? ");
Clf = 2*pi*(AR/(AR+3))*beta;
rhow = input("Fluid density? in kg/m^3 ");
rhos = input("Stabiliser density? in kg/m^3");
rhoc = input("Towed array density? in kg/m^3 ");
xk = xcm*cosd(k)+ycm*sind(k);
yk = -xcm*sind(k)+ycm*cosd(k);
Vf = 2*5*s*0.12*c*(0.1979*(c^1.5)-0.063*(c^2)-0.1172*(c^3)+0.0711*(c^4)-0.0203*(c^5));
B = (Mwb+2*Mwf+2*Mwc)*9.81;
Dc = 0.5*Cdc*rhow*(v^2)*((d)*h);
Lc = 0.5*Clc*rhow*(v^2)*((d)*h);
Db = 0.5*Cd*rhow*(v^2)*l*b;
Lb = 0.5*Cl*rhow*(v^2)*w*b;
Lf = 0.5*Clf*rhow*(v^2)*s*c*cos(beta);
Df = 0.5*Cdf*rhow*(v^2)*s*c*sin(beta);
theta = atand((-2*Lc-Lb-2*Lf-B+W)/(Db+2*Dc+2*Df));
T = (2*Dc+Db+2*Df)/cosd(theta);
A(v) = ( Lc*lc+Dc*bc-Lf*lf-Df*bf-Lb*lb-Db*bb+T*sind(theta)*lt+T*cosd(theta)*bt+B*lB);
xlabel('Towed array force');
ylabel('Variation of A');