clc;
clear all;
clf;
CDir = 'C:\Users\Desktop\plotting\';

Point_Num = 1000; Curve_Num = 5;
f = logspace(6,8, Point_Num); % Hz
%r = linspace(15*10^-9,500*10^-9,Curve_Num); % meter, median of the radius of exosome
r = [14 25 55 100 200].*1e-9; % meter, median of the radius of exosome
%r= 50* 10^-9;

cmap  = [   1 0 1; % m
            0 0 1; % b
            0 0 0; % k
            1 0 0; % r
            0.4940 0.1840 0.5560; % purple
            ];      

for n2 = 1:Curve_Num
    r_cache = r(n2);

omega = 2*pi*f; % rad/s
epsilon_vacuum = 8.85*10^-12; % permitivity of vacuum, F/m
sigma_m = 2*10^-3; % estimated,conductivity of medium, S/m
sigma_p_bulk = 10^-14; % estimated, conductivity of polystyrene, S/m
sigma_p_surface = 1.8*10^-9; % a guess, surface conductance of particle, latex, S/m^2
sigma_p_composite = sigma_p_bulk + 2*sigma_p_surface / r_cache; % surface conductance dominates for sub-micron beads 
epsilon_m = 80 * epsilon_vacuum; % absolute permitivity, F/m
epsilon_p = 2.55 * epsilon_vacuum; % estimated, absolute permitivity, F/m


%initialization
ReCM = zeros(Curve_Num, Point_Num);


for n1 = 1: Point_Num;
    omega_cache = omega(n1);

 
    
        % frequency dependent
        sigma_m_cpx = sigma_m - 1i* omega_cache *epsilon_m;
        sigma_p_composite_cpx = sigma_p_composite - 1i * omega_cache * epsilon_p;
        epsilon_m_cpx = epsilon_m + 1i*(sigma_m / omega_cache);
        epsilon_p_cpx = epsilon_p + 1i*(sigma_p_composite / omega_cache);

        CM_eps = (epsilon_p_cpx - epsilon_m_cpx) / (epsilon_p_cpx +  epsilon_m_cpx); % latex sphere % here is different
        CM_sig = (sigma_p_composite_cpx - sigma_m_cpx)/(sigma_p_composite_cpx +  sigma_m_cpx);

        CM = CM_eps;

        ReCM (n2, n1) = real (CM);
        
     

    
end

% plotting
aet(n2) = semilogx(f./10^6, ReCM(n2,:), 'Color', cmap(n2,:,:), 'LineWidth', 2.5);xlabel('Frequency (MHz)');ylabel('Re(CM)');
hold on;
semilogx(f, zeros(1, Point_Num), 'k', 'LineWidth', 1.5)
hold on;
grid on;

end

legend ([aet(:)], '28 nm', '50 nm' , '100 nm', '200 nm', '400 nm');
legend boxoff;
legend('FontSize', 16);
xlim([1 100]);

