function ScallopTrajectory


AlphaTrajectory = [7,3,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
AlphaTrajectory = [35,35,35,25,10,5,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
AlphaAttack = [6,6,6,7,8,9,10,11,12,15,10,8,6,6,6,6,6];

AlphaAttack = 8+1.5*randn(16,1);
AlphaAttack = [12,12,12,12,11,10,9,9,9,8,8,8,6,6,6,6,6];
AlphaAttack = [12,12,12,12,12,12,6,6,6,6,6,6,6,6,6,6,6];
AlphaAttack = 8+1.5*randn(16,1);
AlphaAttack = [7,8,9,10,11,12,15,10,8,18,16,12,16,6];
AlphaAttack = [12,10,8,6,5,2,2,1,1,1,2,3,2,1,0,0];
AlphaAttack = [12,12,12,10,8,7,6];
AlphaAttack = [8,8,8,8,9,10,12,15,15,15,15,15];
AlphaAttack = [12,12,12,12,12,12,12,12,12,12,12]; %,10,8,6,5,2,2,1,1,1,2,3,2,1,0,0];
AlphaAttack = [6,6,6,7,8,9,10,11,12,15,10,8,6,6,6,6,6];
AlphaAttack = 10+1.5*randn(16,1);
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];
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);
    alpha = 8;
    alpha = alpha*pi/180.0;
    
    rho = 1020;
    g = 9.8;
    AR = 1.117 + 3.6*(SL - 0.04);
    SW = AR*SL;
    RelativeThickness = 0.2;
    ST = RelativeThickness*SL;
    Av = pi*SL*SW/4.0;
    Vs = Av*ST;
    Mb = Av*(100.0*SL - 3.5);
    M = Mb + rho*Vs;
    JetAngle = 12*pi/180.0;
    
    % time in seconds
    T = 0.5;
    nTbins = 1000;
    deltaT = T/nTbins;
    Tthrust = T/3;
    Tglide = T - Tthrust;
    Tthrust = floor(Tthrust/deltaT);
    Tglide = floor(Tglide/deltaT);
    
    
    Am = 0.054*Av;
    LambdaA = 1.0 - Am/Av;
    LambdaAj = (power((0.065244+0.0068523/SL),2))/3.0;
    BetaMax = (0.3*SL)*pi/180.0;
    DBetaDt = BetaMax/(Tthrust*deltaT);
    
    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);
%     F = 10.2*SL;
%     F = 2.8*sqrt(SL);
%  % F = 15*SL^(3/2);
  F =1000*SL^(3.5);
     F =1500*SL^(3.2);
    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);
    %fprintf(1,'alpha = %f\n',alpha);
    %fprintf(1,'CL=%f CD=%f CLalpha=%f\n',CL,CD,(1.21 + 3.8*SL));
    
    X0 = 0;
    Z0 = 0;
    V0x = 0;
    V0z = 0;
    V = sqrt(V0x^2 + V0z^2);
    alphaT = 0.0;
    alphaT = 30.0*pi/180.0;
    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 > 4)
            %alphaT = 0.0;
        end
        
        if (icontraction > 2)
            xvals(icontraction-2,iScallop) = X0;
            zvals(icontraction-2,iScallop) = Z0;
            tvals(icontraction-2,iScallop) = (Tglide+Tthrust)*deltaT*(icontraction-2);
            fprintf(1,'%d %d trajectory velocity = %f %f t=%f position=(%f,%f)\n',iScallop,icontraction-2,X0/((icontraction-2)*(Tglide+Tthrust)*deltaT),V,(Tglide+Tthrust)*deltaT*(icontraction-2),X0,Z0);
        end
        
      %  V0x = 0;
       % V0z = 0;
        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;
            CD = 0.06 + (1.47*CL^2)/(pi*AR);
            CLp = 0.5*rho*CL*Av/M;
            CDp = 0.5*rho*CD*Av/M;
            %alphaT = AlphaTrajectory(icontraction)*pi/180.0;
            %    fprintf(1,'Thrust Velocity = %f (%f,%f) Position=(%f,%f)\n', V,V0x,V0z,X0,Z0);
            X1 = X0 + V0x*deltaT;
            Z1 = Z0 + V0z*deltaT;
            if (V0x < 0)
                V1z = V0z + ( F*cos(JetAngle)*sin(alpha+alphaT)/M - CDp*V*V0z - Mb*g/M)*deltaT;
                V1x = V0x + ( F*cos(JetAngle)*cos(alpha+alphaT)/M - CDp*V*V0x)*deltaT;
            else
                V1z = V0z + ( F*cos(JetAngle)*sin(alpha+alphaT)/M + CLp*V*V0x - CDp*V*V0z - Mb*g/M)*deltaT;
                V1x = V0x + ( F*cos(JetAngle)*cos(alpha+alphaT)/M - CLp*V*V0z - CDp*V*V0x)*deltaT;
            end
            % fprintf(1,'%f thrust direction = %f %f\n',(alpha_thrust_rad+alphaT),sin(alpha_thrust_rad+alphaT),cos(alpha_thrust_rad+alphaT));
            V0x = V1x;
            V0z = V1z;
            V = sqrt(V0x^2 + V0z^2);
            X0 = X1;
            Z0 = Z1;
            
%                 alphaT = abs(asin(V0z/V));
%                 if (V0x < 0 && V0z > 0)
%                     alphaT = 0;
%                 end
%                 if (V0x > 0 && V0z < 0)
%                     alphaT = 0;
%                 end
%                 if (alphaT > 2.0*pi/180)
%                     alphaT = 2.0*pi/180;
%                 end
%               
%              if (alphaT < 0.0*pi/180)
%                  alphaT = 0.0*pi/180;
%              end
%             
      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;
            
            %    fprintf(1,'Coast Velocity = %f (%f,%f) Position=(%f,%f)\n', V,V0x,V0z,X0,Z0);
            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)
            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));
    fprintf('average velocity = %f\n',ds/((ncontractions-1)*(Tglide+Tthrust)*deltaT));
end
