function [speed_t,speed,meanspeed]=model(h_rest,d_rest,squish,h_range,d_range,...
    f,t_c,t_r)
%% Hydrodynamic model information %%
%This hydrodynamic model of jetting jellyfish swimming can be used to
%calculate swimming speeds.

%For more information, see sources:

% N. W. Xu, J. P. Townsend, J. H. Costello, S. P. Colin, B. J. Gemmell, 
%     J. O. Dabiri. (2020). Field testing of biohybrid robotic jellyfish 
%     to demonstrate enhanced swimming speeds. Biomimetics, 5(4), 64. 
%     https://doi.org/10.3390/biomimetics5040064
% N. W. Xu, J. O. Dabiri, Low-power microelectronics embedded in live 
%     jellyfish enhance propulsion. Sci. Adv. 6(5), eaaz3194 (2020).
% T. L. Daniel, Mechanics and energetics of medusan jet propulsion. 
%     Can. J. Zool. 61, 1406–1420 (1983).

%% Input parameters
%morphological and time-dependent jellyfish parameters
%assume that the jellyfish is a semi-hollow hemisphere, i.e., concentric
%hemispheres to account for the subumbrellar volume

%h_rest: resting bell height (cm)
%d_rest: resting bell diameter (cm)
%squish: manubrium tissue height (cm)
%h_range: change in bell height (cm)
%d_range: change in bell diameter (cm)
%f: swimming frequency (Hz)
%t_c: contraction time (s)
%t_r: relaxation time (s)

%% Water parameters
rho=1024*100^(-3);  %saltwater density (kg/cm^3) @ 35 ppt, 21 deg C

%% Heavyside sinusoid function k(t)
TS=30;  %number of timesteps in each second (equiv to fps in video data)
period=1/f;    %time (s)
fallowperiod=period-t_c-t_r;  %jellyfish is motionless during this time

k_c=sin(linspace(0,pi/2,t_c*TS));   %contraction
k_r=sin(linspace(pi/2,0,t_r*TS));   %relaxation

rep=10; %calculate for rep number of periods
if fallowperiod>=0 %case when f<=1/(t_c+t_r)
    k_fallow=zeros(1,round(fallowperiod*TS));
    k_t=cat(2,k_c,k_r,k_fallow);
    k_t=repmat(k_t,1,rep);
else
    %case when f>1/(t_c+t_r)
    %jellyfish cannot fully relax before the next contraction
    partial=round((period-t_c)*TS);
    k_r=k_r(1:partial);
    k_t1=cat(2,k_c,k_r);
    k_c_new=linspace(k_r(partial),1,t_c*TS);
    k_t=cat(2,k_c_new,k_r);
    k_t=repmat(k_t,1,rep-1);
    k_t=cat(2,k_t1,k_t);
end
k_t=smooth(k_t);

%% Plot k_t over time
% figure(1)
% plot((1:size(k_t))/TS,k_t)
% xlabel('Time (s)','FontWeight','Bold')
% ylabel('k(t)','FontWeight','Bold')

%% Calculated jellyfish parameters (instantaneous)
h_t=h_range*k_t+h_rest; %height (cm)
d_t=d_range*(-k_t)+d_rest; %diameter (cm)
A_t=pi*(d_t/2).^2;   %bell area (cm^2)
hcav_t=h_t-squish;  %subumbrellar cavity height (cm)
V_t=2/3*pi*hcav_t.*(d_t/2).^2; %subumbrellar cavity volume (cm^3)
V_jelly=2/3*pi*squish.*(d_t/2).^2;

%check to ensure that all volume calcs are nonnegative
V_change=V_t(1)-(2/3*pi*hcav_t.*(d_t/2).^2);
for chk=1:length(V_change)
    if V_change(chk)<0
        disp('ERROR: difference in volumes is negative')
%         return
    end
end

mass=rho*1.0001*V_jelly;   %mass of jellyfish in kg
%assumption that jellyfish is SLIGHTLY denser than saltwater

%% Plot parameters
% figure(2)
% hold on
% plot((1:size(k_t))/TS,h_t)
% plot((1:size(k_t))/TS,d_t)
% plot((1:size(k_t))/TS,hcav_t)
% legend('jellyfish height','jellyfish diameter','subumbrellar cavity height')
% xlabel('Time (s)','FontWeight','Bold')
% ylabel('Length (cm)','FontWeight','Bold')
% hold off
% 
% figure(3)
% hold on
% plot((1:size(k_t))/TS,A_t)
% legend('subumbrellar area')
% xlabel('Time (s)','FontWeight','Bold')
% ylabel('Area (cm^2)','FontWeight','Bold')
% hold off
% 
% figure(4)
% hold on
% plot((1:size(k_t))/TS,V_t)
% plot((1:size(k_t))/TS,V_jelly)
% plot((1:size(k_t))/TS,V_change)
% legend('subumbrellar volume','jellyfish volume','change in subumbrellar volume')
% xlabel('Time (s)','FontWeight','Bold')
% ylabel('Volume (cm^3)','FontWeight','Bold')
% hold off

%% Set upn onlinear ODE with IC
syms u(t)   %speed

%% Calculate Thrust, T
T=rho./A_t.*([0; diff(V_t).^2]);

%% Calculate Drag, D
cd=0.42;

%% Calculate acceleration rxn force, AR
addmass=(2*h_t./d_t).^1.4;    %added mass coeff based on bell shape

%% ODE eqn
speed_t=zeros(1,length(addmass));
for swift=1:length(addmass)
    Du=diff(u);
    cond=Du(0)==0;
    ode=T(swift)==.5*cd*rho.*A_t(swift).*u.^2+addmass(swift).*...
        rho.*V_jelly(swift).*diff(u,t)+mass(swift).*diff(u,t);
    uSol(t)=dsolve(ode,cond);
%     disp(swift/length(addmass)) %display what percentage is complete
    %calculate the nonnegative speed
    value=abs(double(uSol(swift))); 
    speed_t(swift)=double(value(1));    %speed is at each time step
end

%% %%  Calculate speed
speed=speed_t*TS;   %speed in cm/s
% speed_nondim=speed_t*TS/d_rest;   %speed in body diameters per sec

%% Calculate mean speed
%first period is different from the other periods (which are identical) 

%calculate mean speed based on the last half of the periods only
% meanspeed=mean(speed(length(speed)/2+1:length(speed)));

%calculate mean speed based on all periods
meanspeed=mean(speed);
