function ScallopCharacterization

ScallopLengths = [0.055, 0.07, 0.10, 0.12];

for iScallop = 1: length(ScallopLengths)
    
    SL = ScallopLengths(iScallop);
    fprintf(1,'Scallop %d Length=%f\n',iScallop,SL);
    
    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.6;
    nTbins = 100;
    deltaT = T/nTbins;
    Tthrust = T/3;
    Tglide = T - Tthrust;
    Tglide = floor(Tglide/deltaT);
    Tthrust = floor(Tthrust/deltaT);

    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 = rho*power(LambdaA*SL,3)*Av*power(DBetaDt,2)/(2*LambdaAj);

    alpha = 8*pi/180.0;
    CL = 0.0565 + (1.21 + 3.8*SL)*alpha;
    CD = 0.06 + (1.47*CL^2)/(pi*AR);
   
    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,'CL = %f\n',CL);
    fprintf(1,'CD = %f\n',CD);
    %fprintf(1,'alpha = %f\n',alpha);
    
    
end
