% execute the file by typing the file name on the octave command prompt without the ".m" extension

% this script plots the drag coefficient of a sphere depending on the Kn number

clear; % clears all variables
%clc;
clf; % clears the plots
format long E;

%compute the way points for CD_sphere
CD_Ma_gt_1 = 1.0 - (1.5/18.5) % gives 0.91(891) while the nominal value is 0.96 (Ma = 1 - (1.5cm from 18.5cm on the electronic version));
CD_Ma_le_1 = CD_Ma_gt_1/2.0;
delta1 = 0.01;

Kn_way_point_0 = 1.0e-7;
CD_Kn_pair_0 = [Kn_way_point_0; CD_Ma_le_1];

Kn_way_point_1_coeff = 10^(8.6/10.25); %(6.75mm from 8mm on the paper version or 8.6cm from 10.25cm on the electronic version).
Kn_way_point_1 = (Kn_way_point_1_coeff - delta1)* 1.0e-6
CD_Kn_pair_1 = [Kn_way_point_1; CD_Ma_le_1];
   
Kn_way_point_2 = (Kn_way_point_1_coeff + delta1)* 1.0e-6
CD_Kn_pair_2 = [Kn_way_point_2; CD_Ma_gt_1];

Kn_way_point_3_coeff = 1.0; 
Kn_way_point_3 = (Kn_way_point_3_coeff)* 1.0e-3
CD_Kn_pair_3 = [Kn_way_point_3; CD_Ma_gt_1];

Kn_way_point_4_coeff = 10^(2.7/13.7/2.0); %(2.7cm from 13.7/2.0cm on the electronic version).
Kn_way_point_4 = (Kn_way_point_4_coeff)* 1.0e-3
CD_way_point_4 = 1.0 - (0.6/18.5) %= 0.9(675) = (Ma = 1 - (6mm from 185mm of 1 on the electronic version)
CD_Kn_pair_4 = [Kn_way_point_4; CD_way_point_4];

Kn_way_point_5_coeff = 10^((13.7/2.0-(13.7-7.05))/13.7/2.0); %(((13.7/2.0-(13.7-7.05))cm from 13.7/2 on the paper version).
Kn_way_point_5 = (Kn_way_point_4_coeff)* 1.0e-2
CD_way_point_5 = 1.0 + (1.4/18.5) %= 1.3(837) = (Ma = 1 + (1.4cm from 18.5cm of 1 on the electronic version)
CD_Kn_pair_5 = [Kn_way_point_5; CD_way_point_5];

Kn_way_point_6_coeff = 10^((13.7/2.0-(13.7-12.1))/13.7/2.0); %(((13.7/2.0-(13.7-12.1))cm from 13.7/2.0cm on the electronic version).
Kn_way_point_6 = (Kn_way_point_6_coeff)* 1.0e-2
CD_way_point_6 = 1.0 + (7.1/18.5) %= = (Ma = 1 + (7.1cm from 18.5cm of 1 on the electronic version)
CD_Kn_pair_6 = [Kn_way_point_6; CD_way_point_6];

Kn_way_point_7_coeff = 10^(3.8/13.7/2.0); %((3.8cm from 13.7/2.0cm on the electronic version).
Kn_way_point_7 = (Kn_way_point_7_coeff)* 1.0e-1
CD_way_point_7 = 1.0 + (14.2/18.5) %= 1.7(675) = (Ma = 1 + (14.2cm from 18.5cm of 1 on the electronic version)
CD_Kn_pair_7 = [Kn_way_point_7; CD_way_point_7];

Kn_way_point_7_coeff = 10^((13.7/2.0-(13.7-8.85))/13.7/2.0); %(((13.7/2.0-(13.7-8.85))cm from 13.7/2.0cm on the electronic version).
Kn_way_point_7 = (Kn_way_point_6_coeff)* 1.0e0
CD_way_point_7 = 1.0 + (17.8/18.5) %= 1.9(621) = (Ma = 1 + (17.8cm from 18.5cm of 1 on the electronic version)
CD_Kn_pair_7 = [Kn_way_point_7; CD_way_point_7];

Kn_way_point_8_coeff = 10^((13.7/2.0-(13.7-13.15))/13.7/2.0) %(((13.7/2.0-(13.7-13.15))cm from 13.7/2.0cm on the electronic version).
Kn_way_point_8 = (Kn_way_point_8_coeff)* 1.0e0
CD_way_point_8 = 1.0 + (19.8/18.5) %= 2.0(702)= (Ma = 1 + (19.8cm from 18.5cm of 1 on the electronic version)
CD_Kn_pair_8 = [Kn_way_point_8; CD_way_point_8];

Kn_way_point_9 = 1.0e1
CD_way_point_9 = 2.2
CD_Kn_pair_9 = [Kn_way_point_9; CD_way_point_9];

% put together the plotting data
CD_Kr = horzcat(CD_Kn_pair_0,
        CD_Kn_pair_1,
        CD_Kn_pair_2,
        CD_Kn_pair_3,
        CD_Kn_pair_4,
        CD_Kn_pair_5,
        CD_Kn_pair_6,
        CD_Kn_pair_7,
        CD_Kn_pair_8,
        CD_Kn_pair_9)
        
%plot the data
semilogx(CD_Kr(1,:), CD_Kr(2,:));
set(gca(), 
    'xlim', [1.0e-7,10.0],
    'ylim', [0.0,2.5],
    'xscale', 'log',
    'yscale', 'linear',
%    'xtick', [1.0e-7 1.0e-6 1.0e-5 1.0e-4 1.0e-3 1.0e-2 1.0e-1 1.0 10.0], % comment this out for correct result in v3.2.4
%    'xtick',[1.0e-7,1.0e-5,1.0e-3,1.0e-1,1.0e1], % comment this out for correct result in v3.2.4
    'xgrid', 'on',
    'xminorgrid', 'on',
    'xminortick', 'on',
%    'xtickmode', 'manual', % comment this out for correct result in v3.2.4
    'ygrid', 'on',
    'yminorgrid', 'on',
    'yminortick', 'on' 
    );

%label the plot
title ('C_D of a sphere over Kn number');
xlabel ('Knudsen number [-]');
ylabel ('Drag Coefficient [-]');
legend ({'C_d'}, 'location', 'north');

%debugging info           
get (gca(), 'xtick')

print -dpng drag_coefficient_sphere.png