function [] = RmonitoringPcyclingD1_v20170802()
%Resistance monitoring in power cycling
% Power cycling and monitor electrical resistance with 34420A nanovoltmeter
% for multiple samples. The connection is controlled with a Keithley 7700 
% multiplexer cards installed on a Keithley 2700 multimeter. The example 
% code below is for cycling and monitoring 10 samples. The resistance is
% visualized in real time on PC screen. Temperature is monitored for one
% sample in each cycle. The code works with a custom designed PCB board
% described in my thesis. Modification is needed for different PCB designs.
% By Di Erick Xu, University of Waterloo, 2017.Aug.02

% Hardware: GPIB-USB cable, a Vantek programmable current source, an 
% Agilent 34420A, a Keithley 2700, a Keithley 7700 card.

% Tags: GPIB, Digital multimeter, Multiplexer, Resistance monitoring,
% Power cycling, Temperature monitoring,

%% initialzation
close all;
clear;
clc;
instrreset;
drawnow;

%% Variable definations
totalCycleNum = 720;
heatingSec = 2*60;
coolingSec = 4*60;

sampleName = 'Demo4MathworksCentral'; 
targetVolt = ['11V', num2str(heatingSec), 's_0V', num2str(coolingSec), 's']; 
fname = [datestr(now, 30) '_' sampleName '_' num2str(targetVolt)];

% Card 1 channels for temperature measurement
fresSwitchCH_Card1 = 141;
TempMeterCHs = '108, 128';
Pt1000CHs = {'101, 121', '102, 122', '103, 123', ...
    '104, 124', '105, 125', '106, 126', '107, 127'};

% Card 2 channels for resistance measurement and power cycling
fresSwitchCH_Card2 = 223;
nanoVoltCH = 209;
DCIsourceCH = 210;
sampleICHs = [201:208, 208, 208];
sampleVCHs = 211:220;
ResistorCHs = '202:208';

% temperature and resistance variables
temp = [];
tempTime = [];
xAxis = 0:totalCycleNum; % xAxis is the number of cycles
Datalength = length(xAxis);
resistance = nan(Datalength, length(sampleVCHs)); 

% parameters to calculate temperature from the recorded resistance of Pt1000
A = 3.9083e-3;
B = -5.775e-7;
RTD0 = 1000;

% gpib
nanovoltmeter = gpib('agilent', 7, 22);
multiplexer = gpib('agilent', 7, 25);
tempMeter = gpib('agilent', 7, 4);

% Open communication
fopen(multiplexer);
fopen(tempMeter);
fopen(nanovoltmeter);

%% instruments resets and configurations
% Temperature: Agilent 34401
fprintf(tempMeter, '*rst');
fprintf(tempMeter, '*cls');
fprintf(tempMeter, '*opc?');
fscanf(tempMeter);
fprintf(tempMeter, 'CONFigure:FRES');

% Multiplexer: Keithley 2700 + two 770x cards
fprintf(multiplexer, '*rst');
fprintf(multiplexer, '*cls');
fprintf(multiplexer, '*opc?');
fscanf(multiplexer);
fprintf(multiplexer, 'rout:open:all');

% Nanovoltmeter: Agilent 34420A
fprintf(nanovoltmeter, '*rst');
fprintf(nanovoltmeter, '*cls');
fprintf(nanovoltmeter, '*opc?');
fscanf(nanovoltmeter);
fprintf(nanovoltmeter, 'conf:fres');
fprintf(nanovoltmeter, 'sense:fres:ocom on');
fprintf(nanovoltmeter, 'sense:fres:nplc 10');

figure;
set(gcf,'Position',[50 50 1200 850]);


%% Power cycling, temperature monitoring, resistance monitoring start
tic;
currentCycle = 0;
while ~(currentCycle > totalCycleNum)
    %% Power cycling with temperature monitoring of one sample in each cycle
    if currentCycle ~= 0 
        % connect temperature multimeter (potentially to a Pt1000)
        fprintf(multiplexer, ['rout:mult:clos (@', ...
            num2str(fresSwitchCH_Card1), ',', TempMeterCHs, ')']);
        fprintf(multiplexer, '*opc?');
        fscanf(multiplexer);
        % connect resistors to DC power supply
        fprintf(multiplexer, ['rout:mult:clos (@', ResistorCHs, ',' , ...
             num2str(DCIsourceCH), ')']);
        fprintf(multiplexer, '*opc?');
        fscanf(multiplexer);

        % Choose the sample to monitor the temperature in this cycle, and
        % connect the corresponding Pt1000 to the temperature multimeter
        Pt1000_ThisCycle = mod(currentCycle, length(Pt1000CHs));
        if Pt1000_ThisCycle == 0
            Pt1000_ThisCycle = length(Pt1000CHs);
        end
        fprintf(multiplexer, ['rout:mult:close (@', Pt1000CHs{Pt1000_ThisCycle}, ')']);
        fprintf(multiplexer, '*opc?');
        fscanf(multiplexer);

        % start heating all resistors with the DC power supply
        [~, ~] = dos('Vantek.exe -on -com12');
        ThisCycleStartTime = toc;
        
        % monitor the temperature until heating finishes
        while toc-ThisCycleStartTime < heatingSec
            % ask for a reading from temperature measuring multimeter
            tempR = [];
            while isempty(tempR)
                try
                    fprintf(tempMeter, 'read?');
                    tempR = str2double(fscanf(tempMeter));
                catch ME
                end
            end
            % check if the reading is in a valid range, otherwise, take NaN
            if tempR > 1000 && tempR < 2000
                RR = roots([B*RTD0, A*RTD0, RTD0-tempR]);
                temp = [temp, RR(RR<300)];
            else
                temp = [temp, nan];
            end
            tempTime = [tempTime, toc];
            % visualize the temperature measurement on PC screen
            subplot(1, 2, 1)
            plot(tempTime(max(1, length(tempTime)-1000):end)/60, temp(max(1, length(temp)-1000):end));
            xlabel('Time [min]');
            ylabel('Temperature [°C]');
            drawnow;
        end % end of while
        
        % stop heating all resistors with the DC power supply
        [~, ~] = dos('Vantek.exe -off -com12');
        
        % monitor the temperature until cooling finishes
        while toc-ThisCycleStartTime < heatingSec + coolingSec
            % ask for a reading from temperature measuring multimeter
            tempR = [];
            while isempty(tempR)
                try
                    fprintf(tempMeter, 'read?');
                    tempR = str2double(fscanf(tempMeter));
                catch ME
                end
            end
            % check if the reading is in a valid range, otherwise, take NaN
            if tempR > 1000 && tempR < 2000
                RR = roots([B*RTD0, A*RTD0, RTD0-tempR]);
                temp = [temp, RR(RR<300)];
            else
                temp = [temp, nan];
            end
            tempTime = [tempTime, toc];
            % visualize the temperature measurement on PC screen
            subplot(1, 2, 1)
            plot(tempTime(max(1, length(tempTime)-1000):end)/60, temp(max(1, length(temp)-1000):end));
            xlabel('Time [min]');
            ylabel('Temperature [°C]');
            drawnow;
        end % end of while
        
        % Open all channels
        fprintf(multiplexer, 'rout:open:all');
        fprintf(multiplexer, '*opc?');
        fscanf(multiplexer);
    end % end of if (i.e. Power cycling and Temperature monitoring)
        
    %% recording the resistances
    for i = 1:length(sampleVCHs)
        % Close channels for resistance measurement, one at a time.
        fprintf(multiplexer, ['rout:mult:clos (@', num2str(sampleICHs(i)), ...
            ',', num2str(sampleVCHs(i)), ',', num2str(fresSwitchCH_Card2), ...
            ',', num2str(nanoVoltCH), ')', ]);
        fprintf(multiplexer, '*opc?');
        fscanf(multiplexer);

        % ask for a reading from the nanovoltmeter
        rawStr = [];
        while isempty(rawStr)
            try
                fprintf(nanovoltmeter, 'read?');
                rawStr = fscanf(nanovoltmeter);
            catch ME
            end
        end
        Rindex = currentCycle+1;
        resistance(Rindex, i) = str2double(rawStr);
        
        % Open all channels
        fprintf(multiplexer, 'rout:open:all');
        fprintf(multiplexer, '*opc?');
        fscanf(multiplexer);
    end
    
    %% visualize the resistance vs. number of cycles on PC screen
    subplot(1, 2, 2)
    plot(xAxis(max(1, Rindex-100):Rindex),...
        resistance(max(1, Rindex-100):Rindex, 2:end));
    xlabel('Number of Cycles');
    ylabel('Resistance [Ohm]');
    drawnow;

    currentCycle = currentCycle + 1;
    
    %% Data save
    save([fname '.mat']);
end

%% disconnect communications with the instruments
fclose(multiplexer);
fclose(tempMeter);
fclose(nanovoltmeter);

%% Plot result and save figure
close all
plot(xAxis,resistance);
xlabel('Number of Cycles');
ylabel('Resistance [Ohm]');
figwind=[30          64        1493         729];
set(gcf,'Position',figwind);
saveas(gcf,[fname '_fig'],'fig');
