Lab 2b - gas sensor

clear
close all
clc
    % PROCESS_GRAPHENE_DATA
    % 1) Loads data (time, R, concentration)
    % 2) Splits data into on/off intervals
    % 3) Fits 'on' intervals with a log function
    % 4) Fits 'off' intervals with an exponential function
    % 5) Stores results in a cell array and plots them

    %% 1. Load data
    % data.txt must contain columns: time(sec), resistance(Ohm), concentration(ppm)
    data = load('data.txt');
    t = data(:,1)./1;    % time in seconds (or use t/60 if you prefer minutes)
    R = data(:,2)./1000;    % resistance
    C = data(:,3);    % concentration

    %% 2. Identify gas "on" vs. "off" intervals
    onThreshold = 0.1;      % user-defined threshold for "gas on"
    gasOn = (C > onThreshold);

    % Find indices where we transition from on->off or off->on
    switchPoints = find(diff(gasOn) ~= 0);

    % We'll store intervals in a cell array. Each element is a struct with:
    %   startIndex, endIndex, isOn
    intervals = {};
    startIdx = 1;
    currentState = gasOn(1);

    for i = 1:length(switchPoints)
        endIdx = switchPoints(i);
        intervals{end+1} = struct('startIndex', startIdx, ...
                                  'endIndex',   endIdx, ...
                                  'isOn',       currentState);
        % Next interval will begin after the switch
        startIdx = endIdx + 1;
        currentState = ~currentState;
    end
    % The last interval goes from startIdx to the end of the data
    if startIdx <= length(t)
        intervals{end+1} = struct('startIndex', startIdx, ...
                                  'endIndex',   length(t), ...
                                  'isOn',       currentState);
    end

    %% 3. Prepare a cell array to store fits
    % We'll store, for each interval k:
    %   fits{k}.type       = 'growth' or 'decay'
    %   fits{k}.startIndex
    %   fits{k}.endIndex
    %   fits{k}.fitParams
    %   fits{k}.fitModel   = 'A + B ln(t+1)' or 'C + D exp(-k t)'
    %   fits{k}.timeFit, fits{k}.RFit   (for plotting)
    fits = cell(length(intervals),1);

    %% 4. Fit each interval
    figure('Name','Graphene Growth/Decay Fits','NumberTitle','off');
    hold on; grid on;
    bluef = [0 0 1 0.1];
    plot(t, R,'color',bluef, 'DisplayName','Raw data');  % raw data in background

    for kInt = 1:length(intervals)
        sIdx = intervals{kInt}.startIndex;
        eIdx = intervals{kInt}.endIndex;

        % Extract the subset of data for this interval
        subset_t = t(sIdx:eIdx);
        subset_R = R(sIdx:eIdx);

        % We'll shift time so the start of the interval is t=0
        t0 = subset_t(1);
        shifted_t = subset_t - t0;  % so it goes from 0 up to (end - start)

        % We'll do a finer grid for plotting the fit
        fine_t = linspace(0, shifted_t(end), 200);

        if intervals{kInt}.isOn
            % "Growth" interval => fit a log function
            % R(t) = A + B*ln( (t - t0) + 1 ), but since shifted_t = t - t0,
            % effectively: R(t) = A + B * ln( shifted_t + 1 )

            % Use MATLAB's 'fit' with a custom fittype
            growthModel = fittype(@(A,B,x) A + B*log(x + 1), ...
                                  'independent','x','coefficients',{'A','B'});
            % Initial guesses
            A0 = subset_R(1);
            % Avoid dividing by zero if the interval is very short
            if shifted_t(end) == 0
                B0 = 0;
            else
                B0 = (subset_R(end) - subset_R(1)) / log(shifted_t(end) + 1 + eps);
            end

            fitGrowth = fit(shifted_t, subset_R, growthModel, 'StartPoint',[A0, B0]);

            % Evaluate the fit for plotting
            Rfit = fitGrowth.A + fitGrowth.B * log(fine_t + 1);

            % Store in fits cell
            fits{kInt}.type       = 'growth';
            fits{kInt}.startIndex = sIdx;
            fits{kInt}.endIndex   = eIdx;
            fits{kInt}.fitModel   = 'R(t) = A + B ln(t+1)';
            fits{kInt}.fitParams  = [fitGrowth.A, fitGrowth.B];
            fits{kInt}.timeFit    = t0 + fine_t;   % shift back to original time scale
            fits{kInt}.RFit       = Rfit;

            % Plot the fit
            plot(fits{kInt}.timeFit, fits{kInt}.RFit, '-', 'LineWidth',1.5, ...
                 'DisplayName', sprintf('Log fit (interval %d)', kInt));

        else
            % "Decay" interval => fit an exponential function
            % R(t) = C + D * exp(-k * (t - t0))
            % with shifted_t = t - t0, that becomes: R(t) = C + D * exp(-k * shifted_t)

            decayModel = fittype(@(C,D,k,x) C + D.*exp(-k.*x), ...
                                 'independent','x','coefficients',{'C','D','k'});
            % Initial guesses
            C0 = subset_R(end);             % final (plateau) value
            D0 = subset_R(1) - subset_R(end);
            k0 = 0.01;  % guess a small positive rate

            fitDecay = fit(shifted_t, subset_R, decayModel, 'StartPoint',[C0, D0, k0]);

            % Evaluate the fit for plotting
            Rfit = fitDecay.C + fitDecay.D * exp(-fitDecay.k * fine_t);

            % Store in fits cell
            fits{kInt}.type       = 'decay';
            fits{kInt}.startIndex = sIdx;
            fits{kInt}.endIndex   = eIdx;
            fits{kInt}.fitModel   = 'R(t) = C + D exp(-k t)';
            fits{kInt}.fitParams  = [fitDecay.C, fitDecay.D, fitDecay.k];
            fits{kInt}.timeFit    = t0 + fine_t;
            fits{kInt}.RFit       = Rfit;

            % Plot the fit
            plot(fits{kInt}.timeFit, fits{kInt}.RFit, '-', 'LineWidth',1.5, ...
                 'DisplayName', sprintf('Exp fit (interval %d)', kInt));
        end
    end

    xlabel('Time (s)');
    ylabel('Resistance (k\Omega)');
    legend('Location','best');
    title('Graphene Oxide Resistance Growth/Decay Fits');
    hold off;

    %% 5. Display final results in the Command Window
    fprintf('\n======= Fit Results by Interval =======\n');
    for kInt = 1:length(intervals)
        switch fits{kInt}.type
            case 'growth'
                A_(kInt) = fits{kInt}.fitParams(1);
                B_(kInt) = fits{kInt}.fitParams(2);
                fprintf('Interval %d (Growth): A = %.4g, B = %.4g\n', kInt, A_(kInt), B_(kInt));
            case 'decay'
                C_(kInt) = fits{kInt}.fitParams(1);
                D_(kInt) = fits{kInt}.fitParams(2);
                k_(kInt) = fits{kInt}.fitParams(3);
                fprintf('Interval %d (Decay): C = %.4g, D = %.4g, k = %.4g\n', kInt, C_(kInt), D_(kInt), k_(kInt));
        end
    end
ftype()
    %% calculating stuff for lab report

for i = 1:length(fits)
    R0(i) = fits{i}.RFit(1);
    Rend(i) = fits{i}.RFit(end);
    delta(i) = Rend(i) - R0(i);
    R90(i) = R0(i) + delta(i).*0.9;
    if mod(i,2) == 1
        index(i) = find(fits{i}.RFit > R90(i), 1);
    else
        index(i) = find(fits{i}.RFit < R90(i), 1);
    end
    T0(i) = fits{i}.timeFit(1);
    T90(i) = fits{i}.timeFit(index(i))-T0(i);

    S(i) = abs(delta(i))./min([R0(i) Rend(i)]).*100;
end

T90
S

%%
function ftype()
set(gca,'FontSize',13,'fontWeight','bold')
set(findall(gcf,'type','text'),'FontSize',13,'fontWeight','bold')
end