How to plot contour ?

clc
clear
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define material properties and simulation parameters

n0 = 1;         % Refractive index of air
ns = 1.46;      % Refractive index of subtrate

% Parameters
lambda0=5e-6;        % Wavelength of light in micrometers

frequency = linspace(58,67,1000);

c=3e8;

transmission=zeros(1,1000);
n01 = 1.578;    % Refractive index of PS
n02 = 1.484;    % Refractive index of PMMA
n=25;





    
    nA = n01;
    dA = lambda0/(4*nA); %% Thickness of First Layer in meters

    

        
        nB = n02;
        dB = lambda0/(4*nB); %% Thickness of Second Layer in meters

        for i = 1:length(frequency)
            f = frequency(i);

            w=2*pi*f*1e12;        %%% Angular frequency by frequency

            DAA=dA * nA * (w/c);
            DBB=dB * nB * (w/c);

            %%% Transfer Matrix elements of first layer
            ma11=cos(DAA); ma12=-1i*sin(DAA)/nA; ma21=-1i*nA*sin(DAA); ma22=cos(DAA);
            MA=[ma11 ma12; ma21 ma22];

            %%% Transfer MAtrix elements of Second layer
            lb11=cos(DBB); lb12=-1i*sin(DBB)/nB; lb21=-1i*nB*sin(DBB); lb22=cos(DBB);
            MB=[lb11 lb12; lb21 lb22];            

            M_total_P = (MB*MA)^n*(MA*MB)^n;

            T1 = (2*ns/(ns*M_total_P(1,1)+ns*n0*M_total_P(1,2)+M_total_P(2,1)+n0*M_total_P(2,2)));
            T=(n0 / ns) * abs(T1)^2;
            transmission(i) = T;
        
        end
    

plot(frequency,transmission)

this code is transmission spectra

I need a contour plot of this by varying the n01 and n02 with respect to Q factor

I did a program but is not working

can anyone help me

clc
clear
close all
tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define material properties and simulation parameters
P = 0;          %%% Hydrostatic pressure in MPa
n0 = 1;         % Refractive index of air
ns = 1.46;      % Refractive index of subtrate

% Parameters
lambda0=5e-6;        % Wavelength of light in micrometers



%%% frequency Range
% f = 61.25;                   %% frequency in THz
frequency = linspace(58,67,1000);

c=3e8;

transmission=zeros(1,1000);
n01_ = linspace(1.3,4,25);    % Refractive index of PS
n02_ = linspace(1.3,4,25);    % Refractive index of PMMA


n=25;



% Q factor array
Q_factor_ = zeros(length(n01_), length(n02_));

for ii = 1:length(n01_)
    n01 = n01_(ii);

    
    nA = n01;
    dA = lambda0/(4*nA); %% Thickness of First Layer in meters (varying with respect to n01)

    for jj = 1:length(n02_)
        n02 = n02_(jj);

        
        nB = n02;
        dB = lambda0/(4*nB); %% Thickness of Second Layer in meters (varying with respect to n02)



        for i = 1:length(frequency)
            f = frequency(i);

            w=2*pi*f*1e12;        %%% Angular frequency by frequency

            DAA=dA * nA * (w/c); %%% (varying with respect to n01)
            DBB=dB * nB * (w/c); %%% (varying with respect to n02)




            %%% Transfer Matrix elements of first layer
            ma11=cos(DAA); ma12=-1i*sin(DAA)/nA; ma21=-1i*nA*sin(DAA); ma22=cos(DAA);
            MA=[ma11 ma12; ma21 ma22];
            



            %%% Transfer MAtrix elements of Second layer
            lb11=cos(DBB); lb12=-1i*sin(DBB)/nB; lb21=-1i*nB*sin(DBB); lb22=cos(DBB);
            MB=[lb11 lb12; lb21 lb22];
            

            

            M_total_P = (MB*MA)^n*(MA*MB)^n;

            T1 = (2*ns/(ns*M_total_P(1,1)+ns*n0*M_total_P(1,2)+M_total_P(2,1)+n0*M_total_P(2,2)));
            T=(n0 / ns) * abs(T1)^2;
            transmission(i) = T;
        

        
        
        desired_freq = 60;
        peak_idx = find((frequency) == (desired_freq));
        
        peak_val = transmission(peak_idx);
        half_max = peak_val / 2;
        peak_left_edge = find(transmission(1:peak_idx) < half_max, 1, 'last');
        peak_right_edge = find(transmission(peak_idx:end) < half_max, 1, 'first') + peak_idx - 1;
        fwhm = frequency(peak_right_edge) - frequency(peak_left_edge); %%% full width at half maximum

        % Calculate Q factor
        Q = desired_freq ./ fwhm;

        % Store Q factor in array
        Q_factor_(:,:,i) = Q;
        end
    end
end
contourf(n01_,n02_,Q_factor_)

Answer 1

[Manikanta Aditya]

Hi @GULZAR,

The issue with your contour plot attempt lies in the incorrect handling of the Q_factor_ array. You're attempting to store the Q factor in a 3D array (Q_factor_(:,:,i) = Q;) inside a loop that iterates over frequency, but your intention seems to be to create a 2D matrix where each element corresponds to a combination of n01 and n02 indices. The Q factor calculation should be based on the transmission spectrum characteristics for each combination of n01 and n02, not for each frequency.

Check this code with plot of contour:

clc;
clear;
close all;
tic;

% Define material properties and simulation parameters
n0 = 1;         % Refractive index of air
ns = 1.46;      % Refractive index of substrate
lambda0 = 5e-6; % Wavelength of light in micrometers
frequency = linspace(58,67,1000); % Frequency range
c = 3e8; % Speed of light
n01_ = linspace(1.3,4,25); % Refractive index of PS range
n02_ = linspace(1.3,4,25); % Refractive index of PMMA range
n = 25; % Number of layers
Q_factor_ = zeros(length(n01_), length(n02_)); % Initialize Q factor matrix

for ii = 1:length(n01_)
    for jj = 1:length(n02_)
        transmission = zeros(1,length(frequency)); % Reset transmission for each n01, n02 combination
        n01 = n01_(ii);
        n02 = n02_(jj);
        dA = lambda0 / (4 * n01); % Thickness of First Layer in meters
        dB = lambda0 / (4 * n02); % Thickness of Second Layer in meters
        
        for i = 1:length(frequency)
            f = frequency(i);
            w = 2 * pi * f * 1e12; % Angular frequency
            DAA = dA * n01 * (w / c);
            DBB = dB * n02 * (w / c);
            % Transfer Matrix for first layer
            ma11 = cos(DAA); ma12 = -1i * sin(DAA) / n01;
            ma21 = -1i * n01 * sin(DAA); ma22 = cos(DAA);
            MA = [ma11 ma12; ma21 ma22];
            % Transfer Matrix for second layer
            lb11 = cos(DBB); lb12 = -1i * sin(DBB) / n02;
            lb21 = -1i * n02 * sin(DBB); lb22 = cos(DBB);
            MB = [lb11 lb12; lb21 lb22];
            % Total Transfer Matrix
            M_total_P = (MB * MA)^n * (MA * MB)^n;
            T1 = (2 * ns / (ns * M_total_P(1,1) + ns * n0 * M_total_P(1,2) + M_total_P(2,1) + n0 * M_total_P(2,2)));
            T = (n0 / ns) * abs(T1)^2;
            transmission(i) = T;
        end
        
        % Simplified Q factor calculation (for demonstration)
        % Ideally, you should find the peak and its FWHM more accurately
        [maxT, idx] = max(transmission);
        desired_freq = frequency(idx);
        fwhm = 1; % Placeholder for FWHM calculation
        Q = desired_freq / fwhm;
        Q_factor_(ii, jj) = Q; % Store Q factor
    end
end

% Plotting the contour
figure;
contourf(n01_, n02_, Q_factor_', 20); % Transpose Q_factor_ for correct orientation
colorbar;
xlabel('n01');
ylabel('n02');
title('Q Factor Contour');

toc;

Elapsed time is 10.374777 seconds.

Hope this helps you!

[GULZAR]

I need the Q factor according to full width at half maximum

[Manikanta Aditya]

Check this script for Q factor according to full width at half maximum, if it helps you can accept the answer:

The Q factor based on the full width at half maximum (FWHM) of the transmission peak for each combination of n01 and n02, and then plots a contour map of these Q factors.

clc;
clear;
close all;
tic;

% Define material properties and simulation parameters
n0 = 1;         % Refractive index of air
ns = 1.46;      % Refractive index of substrate
lambda0 = 5e-6; % Wavelength of light in micrometers
frequency = linspace(58,67,1000); % Frequency range
c = 3e8; % Speed of light
n01_ = linspace(1.3,4,25); % Refractive index of PS range
n02_ = linspace(1.3,4,25); % Refractive index of PMMA range
n = 25; % Number of layers
Q_factor_ = zeros(length(n01_), length(n02_)); % Initialize Q factor matrix

for ii = 1:length(n01_)
    for jj = 1:length(n02_)
        transmission = zeros(1,length(frequency)); % Reset transmission for each n01, n02 combination
        n01 = n01_(ii);
        n02 = n02_(jj);
        dA = lambda0 / (4 * n01); % Thickness of First Layer in meters
        dB = lambda0 / (4 * n02); % Thickness of Second Layer in meters
        
        for i = 1:length(frequency)
            f = frequency(i);
            w = 2 * pi * f * 1e12; % Angular frequency
            DAA = dA * n01 * (w / c);
            DBB = dB * n02 * (w / c);
            % Transfer Matrix for first layer
            ma11 = cos(DAA); ma12 = -1i * sin(DAA) / n01;
            ma21 = -1i * n01 * sin(DAA); ma22 = cos(DAA);
            MA = [ma11 ma12; ma21 ma22];
            % Transfer Matrix for second layer
            lb11 = cos(DBB); lb12 = -1i * sin(DBB) / n02;
            lb21 = -1i * n02 * sin(DBB); lb22 = cos(DBB);
            MB = [lb11 lb12; lb21 lb22];
            % Total Transfer Matrix
            M_total_P = (MB * MA)^n * (MA * MB)^n;
            T1 = (2 * ns / (ns * M_total_P(1,1) + ns * n0 * M_total_P(1,2) + M_total_P(2,1) + n0 * M_total_P(2,2)));
            T = (n0 / ns) * abs(T1)^2;
            transmission(i) = T;
        end
        
        % Find peak and FWHM
        [peakT, peakIdx] = max(transmission);
        peakFreq = frequency(peakIdx);
        
        % Find left and right half max points
        halfMax = peakT / 2;
        leftIdx = find(transmission(1:peakIdx) <= halfMax, 1, 'last');
        rightIdx = find(transmission(peakIdx:end) <= halfMax, 1, 'first') + peakIdx - 1;
        
        % Ensure fwhmFreq is defined before use
        if isempty(leftIdx) || isempty(rightIdx)
            fwhmFreq = NaN; % Assign NaN or another placeholder if FWHM can't be computed
        else
            fwhmFreq = frequency(rightIdx) - frequency(leftIdx);
        end
        
        % Calculate Q factor only if fwhmFreq is not NaN
        if ~isnan(fwhmFreq) && fwhmFreq > 0
            Q = peakFreq / fwhmFreq;
            Q_factor_(ii, jj) = Q; % Store Q factor
        else
            Q_factor_(ii, jj) = NaN; % Assign NaN or a placeholder value indicating calculation was not possible
        end
    end
end

% Plotting the contour
figure;
contourf(n01_, n02_, Q_factor_', 20); % Transpose Q_factor_ for correct orientation
colorbar;
xlabel('n01');
ylabel('n02');
title('Q Factor Contour');

toc;

Elapsed time is 9.891317 seconds.