function ScallopProject

AlphaAttack = [6,8,7,8,5,6,6,4,6,6,6,5,6,5,6];
ncontractions = length(AlphaAttack)-2;
ScallopLengths = [0.055, 0.06, 0.07, 0.10, 0.12]; % shell length m
xvals = zeros(ncontractions,length(ScallopLengths));
zvals = zeros(ncontractions,length(ScallopLengths));
tvals = zeros(ncontractions,length(ScallopLengths));

for iScallop = 1: length(ScallopLengths)
    SL = ScallopLengths(iScallop);
    fprintf(1,'Scallop %d Length=%f\n',iScallop,SL);
    
    rho = 1020; % density of sea water kg/m3
    g = 9.8; % acceleration due to gravity m/s2
    AR = 1.117 + 3.6*(SL - 0.04); % aspect ratio
    SW = AR*SL; % shell width m 
    RelativeThickness = 0.2; % height to length ratio
    ST = RelativeThickness*SL; % shell height
    Av = pi*SL*SW/4.0; % projected area of top valve m2
    Vs = Av*ST; % estimated closed volume
    Mb = Av*(100.0*SL - 3.5); % buoyant weight fit to data
                              % by dadsweil and weihs 1989
    M = Mb + rho*Vs; % mass of closed scallop
    JetAngle = 12*pi/180.0; % angle of jets cheng and demont 1996
    
    % time in seconds
    % subdividing the scallop open-close-glide cycle
    T = 0.5;
    nTbins = 1000;
    deltaT = T/nTbins;
    Tthrust = T/3;
    Tglide = T - Tthrust;
    Tthrust = floor(Tthrust/deltaT);
    Tglide = floor(Tglide/deltaT);
    T = 0.5;
    nTbins = 300;
    deltaT = T/nTbins;
    Tthrust = T/3;
    Tglide = T - Tthrust;
    Tthrust = floor(Tthrust/deltaT);
    Tglide = floor(Tglide/deltaT);
    
    %
    % cheng and demont's propulsive thrust model - gives advantage 
    % to larger scallops although their analysis concludes contrary
    %
    Am = 0.054*Av;
    LambdaA = 1.0 - Am/Av;
    BetaMax = (0.3*SL)*pi/180.0;
    DBetaDt = BetaMax/(Tthrust*deltaT);
    LambdaAj = (power((0.065244+0.0068523/SL),2))/3.0;
    LambdaAj = 1.0/(-4.1094*power(SL*100,2) + 81.592*(SL*100.0) - 243.46);
    LambdaAj = 0.0003*power(SL,-1.226);  
    F = 50000*rho*power(LambdaA*SL,2)*Av*power(DBetaDt,2)/(2*LambdaAj);

    %
    % Propulsive Force taken as nearly proportional to volume
    % and weighted to more or less match observations
    %   F =1500*SL^(3.21);
    F =2500*SL^(3.21);
    
    fprintf(1,'Scallop Dimensions: SL=%f SW=%f ST=%f\n',SL,SW,ST);
    fprintf(1,'                    Av=%f Volume=%f\n',Av,Vs);
    fprintf(1,'Scallop Mass: Mb=%f M=%f\n',Mb,M);
    fprintf(1,'Force of thrust: %f\n',F);
    fprintf(1,'Lambda Aj %f  %f\n',LambdaAj,LambdaAj*Av);
    
    X0 = 0; % start position 
    Z0 = 0;
    V0x = 0; % start velocity
    V0z = 0;
    V = sqrt(V0x^2 + V0z^2);
    for icontraction = 1:ncontractions+2
        
        alpha = pi*AlphaAttack(icontraction)/180.0;
        %alphaT = pi*AlphaTrajectory(icontraction)/180.0;
        
        if (icontraction == 3)
            X0 = 0;
            Z0 = 0;
        end
        
        if (icontraction > 2)
            xvals(icontraction-2,iScallop) = X0;
            zvals(icontraction-2,iScallop) = Z0;
            tvals(icontraction-2,iScallop) = (Tglide+Tthrust)*deltaT*(icontraction-2);
        end
        
        AlphaAttackVar = pi*(AlphaAttack(icontraction)+2*randn(Tthrust,1))/180.0;
        for i = 1:Tthrust
            alpha = AlphaAttackVar(i);
            CL = 0.0565 + (1.21 + 3.8*SL)*alpha; % Lift Coef (Cheng&Demont 1996)
            CD = 0.06 + (1.47*CL^2)/(pi*AR); % Drag Coef (Cheng&Demont 1996)
            CLp = 0.5*rho*CL*Av/M;
            CDp = 0.5*rho*CD*Av/M;
            X1 = X0 + V0x*deltaT;
            Z1 = Z0 + V0z*deltaT;
            if (V0x < 0)
                V1z = V0z + ( F*cos(JetAngle)*sin(alpha)/M - CDp*V*V0z - Mb*g/M)*deltaT;
                V1x = V0x + ( F*cos(JetAngle)*cos(alpha)/M - CDp*V*V0x)*deltaT;
            else
                V1z = V0z + ( F*cos(JetAngle)*sin(alpha)/M + CLp*V*V0x - CDp*V*V0z - Mb*g/M)*deltaT;
                V1x = V0x + ( F*cos(JetAngle)*cos(alpha)/M - CLp*V*V0z - CDp*V*V0x)*deltaT;
            end
            V0x = V1x;
            V0z = V1z;
            V = sqrt(V0x^2 + V0z^2);
            X0 = X1;
            Z0 = Z1;
      end
      AlphaAttackVar = pi*(AlphaAttack(icontraction)+5*randn(Tglide,1))/180.0;

        for i = 1:Tglide
            alpha = AlphaAttackVar(i);
            CL = 0.0565 + (1.21 + 3.8*SL)*alpha;
            CD = 0.06 + (1.47*CL^2)/(pi*AR);
            CLp = 0.5*rho*CL*Av/M;
            CDp = 0.5*rho*CD*Av/M;
            X1 = X0 + V0x*0.001;
            Z1 = Z0 + V0z*0.001;
            if (V0x < 0)
                V1z = V0z + ( -CDp*V*V0z - Mb*g/M)*deltaT;
                V1x = V0x + ( -CDp*V*V0x)*deltaT;
            else
                V1z = V0z + (  CLp*V*V0x - CDp*V*V0z - Mb*g/M)*deltaT;
                V1x = V0x + ( -CLp*V*V0z - CDp*V*V0x)*deltaT;
            end
            V0x = V1x;
            V0z = V1z;
            V = sqrt(V0x^2 + V0z^2);
            X0 = X1;
            Z0 = Z1;
        end
        %
        % The 'I'm getting tired!" factor
        %if (icontraction > 3 && icontraction < 7)
        if (icontraction > 3 && icontraction < 9)
            
            Tglide = 1.2*Tglide
        end
    end
end
figure(1);
hold on;
colorvals = ['r', 'g', 'b', 'k', 'm'];
for iScallop = 1: length(ScallopLengths)
    ds = 0;
    for icontraction = 2: ncontractions
        dx = xvals(icontraction,iScallop)-xvals(icontraction-1,iScallop);
        dz = zvals(icontraction,iScallop)-zvals(icontraction-1,iScallop);
        ds = ds + sqrt(dx*dx + dz*dz);
        fprintf(1,'ds = %f v=%f\n',ds,ds/((icontraction-1)*(Tglide+Tthrust)*deltaT));
    end
    plot(xvals(:,iScallop),zvals(:,iScallop),'color',colorvals(iScallop),'LineWidth',3);
    fprintf('average velocity = %f\n',ds/((ncontractions-1)*(Tglide+Tthrust)*deltaT));
end
plot([0 12],[0 0],'color',[.5 .5 .1],'LineWidth',5);
legend('0.055 m', '0.06 m', '0.07 m','0.10 m', '0.12 m');
xlabel('Length along ocean floor [m]','fontsize',14);
ylabel('Height above ocean floor [m]','fontsize',14);
title('Scallop Trajectories 13 Contractions','fontsize',14);
