%======================================================================

%                    D R A G M O D E L . M 

%                    doc: Wed Jan  7 16:23:52 2009

%                    dlm: Wed Jan  7 16:24:07 2009

%                    (c) 2009 A.M. Thurnherr

%                    uE-Info: 11 0 NIL 0 0 72 0 2 4 NIL ofnI

%======================================================================

% CHANGES BY ANT:

%   Jan  7, 2009: - tightened use of exist()

% 

%	dragmodel

%       predict the position of the CTD given the ships position and CTD depth

% 

% DEFINE DRAG MODEL constants

%

if ~exist('p','var')

 p.name='dragmodel';

end

% CTD/LADCP parameter

%  Weight of CTD/LADCP system in kg

p=setdefv(p,'dmodel_package_weight',400);

%  Horizontal area of CTD/LADCP system in m^2

p=setdefv(p,'dmodel_package_area_h',.2);

%  Vertical area of CTD/LADCP system in m^2

p=setdefv(p,'dmodel_package_area_v',0.1);

% drag coefficient for CTD/LADCP system (assume cylinder)

p=setdefv(p,'dmodel_package_cd',1.3);

% Cable parameter

% Cable weight in kg per meter length

p=setdefv(p,'dmodel_wire_weight',0.5);

% Cable wire diameter fraction of inch

p=setdefv(p,'dmodel_wire_diameter_inch',0.5);

% Cable thickness in m is equal to area 

p=setdefv(p,'dmodel_wire_area',0.0254*p.dmodel_wire_diameter_inch);

% Cable drag coefficient

p=setdefv(p,'dmodel_wire_cd',1.3);

% Model parameter:

% length of wire segments

p=setdefv(p,'dmodel_wire_length',10);

% Model design:

%

% each element has an angle with the vertical: alpha 

%    where 0 means pointing straight down

%

% each element had an horizontal angle: beta

%   where 0 means pointing to the north and positive angles point east

% each element has a u,v,w velocity defined at its center

%

% At each junction between elements the forces are calculated

%  in the north/east/up direction units are N=kg*m/s^2

%

% Acceleration of gravity  m/s^2

g=9.81;

% Density of water

rho0=1028;

if ~exist('al','var')

 % initialize model

 ncable=10;

 cable_l=ones(ncable,1)*p.dmodel_wire_length;

 % angle to the vertical

 al=zeros(ncable+1,1);

 % angle in horizontal

 be=zeros(ncable+1,1);

 % node velocity

 u=zeros(ncable+1,1);

 v=zeros(ncable+1,1);

 % node acceleration

 dudt=zeros(ncable+1,1);

 dvdt=zeros(ncable+1,1);

 % node position

 X=zeros(ncable+1,1);

 Y=zeros(ncable+1,1);

 Z=zeros(ncable+1,1);

 % node horizontal distance

 H=zeros(ncable+1,1);

 % node vertical distance

 D=zeros(ncable+1,1);

end

if ~exist('dr','var')

 dr.z=[0 3000];

 dr.u=[0 0];

 dr.v=[0 0];

 d.wm=[0 0];

end

% time step

nt=1;

dt=60;

% precalculate 

sal=sin(al*pi/180);

cal=cos(al*pi/180);

sbe=sin(be*pi/180);

cbe=cos(be*pi/180);

% get depth of each node

Z=al*0;

for i=1:length(cable_l)

 in=length(cable_l)-i+1;

 Z(in)=Z(in+1)-p.dmodel_wire_length*cos(al(in));

end

% velocity of nodes

vel=u+sqrt(-1)*v;

% interpolate ocean velocity to drag model points

uo=interp1(-dr.z,dr.u,Z-p.dmodel_wire_length/2,'linear','extrap');

vo=interp1(-dr.z,dr.v,Z-p.dmodel_wire_length/2,'linear','extrap');

velo=uo+sqrt(-1)*vo;

wctd=d.wm(nt);

%  Compute forces from bottom of CTD up to the ship

% 

% CTD/LADCP first

% gravity

Fv(1)=-p.dmodel_package_weight*g;

% drag

% normal drag Force in direction of flow

Fdn(1)=(-sal(1)*wctd+cal(1)*abs(vel(1)-velo(1))) * ...

       abs(sal(1)*wctd+cal(1)*abs(vel(1)-velo(1))) * ...

      p.dmodel_package_area_h*p.dmodel_package_cd*rho0;

% along drag Force

Fda(1)=(-cal(1)*wctd+sal(1)*abs(vel(1)-velo(1))) * ...

       abs(cal(1)*wctd+sal(1)*abs(vel(1)-velo(1))) * ...

      p.dmodel_package_area_h*p.dmodel_package_cd*rho0;

% project in earth coordinates

% vertical drag

Fv(1)=Fv(1)+Fdn(1)*sal(1)+Fda(1)*cal(1);

% tension

Ft(1)=Fda(1)*sal(1)+Fv(1)*cal(1);

% horizontal drag

Fh(1)=Fdn(1)*cal(1)+Ft(1)*sal(1);

Fhx(1)=Fh(1)*cos(angle(vel(1)-velo(1)));

Fhy(1)=Fh(1)*sin(angle(vel(1)-velo(1)));

% loop over the wire elements to calculate Forces and Tension

for i=1:length(cable_l)

 i1=i+1;

 % gravity

 Fv(i1)=Fv(i)-p.dmodel_wire_weight*p.dmodel_wire_length*g;

 % drag

 % normal drag Force in direction of flow

 Fdn(i1)=(-sal(i1)*wctd+cal(i1)*abs(vel(i1)-velo(i1))) * ...

         abs(sal(i1)*wctd+cal(i1)*abs(vel(i1)-velo(i1))) * ...

      p.dmodel_wire_area*p.dmodel_wire_cd*p.dmodel_wire_length;

 % along drag Force

 Fda(i1)=(-cal(i1)*wctd+sal(i1)*abs(vel(i1)-velo(i1))) * ...

         abs(cal(i1)*wctd+sal(i1)*abs(vel(i1)-velo(i1))) * ...

      p.dmodel_wire_area*p.dmodel_wire_cd*p.dmodel_wire_length;

 % project in earth coordinates

 % vertical drag

 Fv(i1)=Fv(i1)+Fdn(i1)*sal(i1)+Fda(i1)*cal(i1);

 % tension

 Ft(i1)=Fda(i1)*sal(i1)+Fv(i1)*cal(i1);

 % horizontal drag

 Fh(i1)=Fdn(i1)*cal(i1)+Ft(i1)*sal(i1);

 Fhx(i1)=Fh(i1)*cos(angle(vel(i1)-velo(i1)));

 Fhy(i1)=Fh(i1)*sin(angle(vel(i1)-velo(i1)));

end

% loop over wire elements to calculate accelation from F = m*a

%  a = F/m

% use "dynamical mass" Fv/g

dudt(1) = Fhx(1)/abs(Fv(1)/g);

dvdt(1) = Fhy(1)/abs(Fv(1)/g);

for i=1:(length(cable_l))

  i1=i+1;

  dudt(i1) = 0.5*(Fhx(i1)+Fhx(i))/abs(Fv(i1)/g);

  dvdt(i1) = 0.5*(Fhy(i1)+Fhy(i))/abs(Fv(i1)/g);

end 

% loop to update velocity of nodes

for i=0:(length(cable_l))

  i1=i+1;

  u(i1) = u(i1) + dudt(i1)*dt/1000;

  v(i1) = v(i1) + dvdt(i1)*dt/1000;

end 

% loop to calculate horizontal position of nodes

for i=0:(length(cable_l))

  i1=i+1;

  X(i1)=X(i1)+u(i1)*dt;

  Y(i1)=Y(i1)+v(i1)*dt;

end

% loop to calculate vertical position of nodes

Z(length(cable_l)+1)=0;

for i=1:(length(cable_l))

  in=length(cable_l)+1-i;

  in1=in+1;

  H(in1)=sqrt((X(in1)-X(in)).^2 + (Y(in1)-Y(in)).^2);

  D(in1)=sqrt(p.dmodel_wire_length.^2 -  H(in).^2);

  Z(in)=Z(in1)-D(in1);

  al(in)=atan2(H(in),D(in))*180/pi;

  be(in)=atan2(Y(in1)-Y(in),X(in1)-X(in))*180/pi;

end