% Define instrument VISA address. The VISA address of the instrument
% may be obtained from the instrument's user interface or your VISA
% configuration utility
instrumentVISAAddress = 'TCPIP0::127.0.0.1::inst0::INSTR';
% Define frequency range of 2.3GHz to 2.6GHz
frequencyRange = [2.3e9 2.6e9];
% Number of points in measurement
numPoints = 401;

% Create a VISA connection to interface with instrument
instrObj = visa('agilent',instrumentVISAAddress);
% Set up connection parameters for transfer of measurement data from the
% instrument
instrObj.InputBufferSize = 10e6;
instrObj.ByteOrder = 'littleEndian';
% Open connection to the instrument and clear hardware buffer of instrument
fopen(instrObj);
clrdevice(instrObj);

% Display information about instrument
IDNString = query(instrObj,'*IDN?');
fprintf('Connected to: %s\n',IDNString);

% Prompt user to ask if they want to calibrate the instrument
doCalibration = questdlg('Do you want to calibrate your PNA?', 'ExampleMeasureSparameter','yes','no','no');
if strcmpi(doCalibration,'yes')
    % Launch wizard-guided calibration on the instrument and prompt user to confirm completion
    fprintf(instrObj,'SYSTem:CORR:WIZ MAIN');
    hMsgBox = msgbox('Press OK when calibration is complete', 'ExampleMeasureSparameter','help','modal');
    % Wait till user has confirmed that the calibration is complete
    uiwait(hMsgBox);
end

% Perform a System Preset
fprintf(instrObj,'SYSTem:FPReset');
fprintf(instrObj,'*CLS');

% Wait till system is ready as Preset could take time
opcStatus = 0;
while(~opcStatus)
    opcStatus = str2double(query(instrObj, '*OPC?'));
end

% Define a measurement name and parameter
fprintf(instrObj,'CALCulate:PARameter:DEFine:EXT ''SParamMeasurementS11'',S11');
fprintf(instrObj,'CALCulate:PARameter:DEFine:EXT ''SParamMeasurementS12'',S12');
fprintf(instrObj,'CALCulate:PARameter:DEFine:EXT ''SParamMeasurementS21'',S21');
fprintf(instrObj,'CALCulate:PARameter:DEFine:EXT ''SParamMeasurementS22'',S22');

% Create a new display window and turn it on
fprintf(instrObj,'DISPlay:WINDow1:STATE ON');

% Associate the measurements to WINDow1
fprintf(instrObj,'DISPlay:WINDow1:TRACe1:FEED ''SParamMeasurementS11''');
fprintf(instrObj,'DISPlay:WINDow1:TRACe2:FEED ''SParamMeasurementS12''');
fprintf(instrObj,'DISPlay:WINDow1:TRACe3:FEED ''SParamMeasurementS21''');
fprintf(instrObj,'DISPlay:WINDow1:TRACe4:FEED ''SParamMeasurementS22''');

% Turn ON the Title, Frequency, and Trace Annotation to allow for
% visualization of the measurements on the instrument display
fprintf(instrObj,'DISPlay:WINDow1:TITLe:STATe ON');
fprintf(instrObj,'DISPlay:ANNotation:FREQuency ON');
fprintf(instrObj,'DISPlay:WINDow1:TRACe1:STATe ON');
fprintf(instrObj,'DISPlay:WINDow1:TRACe2:STATe ON');
fprintf(instrObj,'DISPlay:WINDow1:TRACe3:STATe ON');
fprintf(instrObj,'DISPlay:WINDow1:TRACe4:STATe ON');

% Turn OFF averaging
fprintf(instrObj,'SENSe1:AVERage:STATe OFF');

% Set the number of points
fprintf(instrObj, sprintf('SENSe:SWEep:POINts %s',num2str(numPoints)));

% Set the frequency ranges
fprintf(instrObj, sprintf('SENSe:FREQuency:STARt %sHz',num2str(frequencyRange(1))));
fprintf(instrObj, sprintf('SENSe:FREQuency:STOP %sHz',num2str(frequencyRange(2))));

% Select measurements and set measurement trigger to immediate
fprintf(instrObj,'CALCulate:PARameter:SELect ''SParamMeasurementS11''');
fprintf(instrObj,'CALCulate:PARameter:SELect ''SParamMeasurementS12''');
fprintf(instrObj,'CALCulate:PARameter:SELect ''SParamMeasurementS21''');
fprintf(instrObj,'CALCulate:PARameter:SELect ''SParamMeasurementS22''');
fprintf(instrObj,'TRIG:SOURce IMMediate');

% Autoscale display
fprintf(instrObj, 'DISPlay:WIND:Y:AUTO');

% Select a single sweep across the frequency range to trigger a measurement
fprintf(instrObj,':SENSe:SWEep:MODE SINGLE');

% Since the instrument may take time to make the measurement, wait until it
% is done before requesting measurement data
opcStatus = 0;
while(~opcStatus)
    opcStatus = str2double(query(instrObj, '*OPC?'));
end

% Set instrument to return the data back using binblock format
fprintf(instrObj, 'FORMat REAL,64');

% Set byte order to swapped (little-endian) format. SWAPped is required
% when using IBM compatible computers
fprintf(instrObj, 'FORMat:BORDer SWAP');

% Request 2-port measurement data from instrument
fprintf(instrObj, 'CALC:DATA:SNP:PORTs? ''1,2''');

% Read the measured data
rawDataDB = binblockread(instrObj, 'double'); 
fread(instrObj,1);

% Read back the number of points in the measurement and reshape the
% measurement data
numPoints = str2double(query(instrObj,' SENSe:SWEep:POINts?'));
% Reshape measurement data to [frequency, real, imag] array
rawDataDB = reshape(rawDataDB, numPoints, 9);

% Loop through the error queue and display all errors
systemError = '';
while isempty(strfind(lower(systemError),'no error'))
    systemError = query(instrObj,'SYSTem:ERRor?');
    fprintf('System Error(s): %s',systemError);
end

% Close, delete, and clear instrument connections.
fclose(instrObj);
delete(instrObj);
clear instrObj;

% Store frequency range of the measurements
freqRange = rawDataDB(:,1);

% Convert retrieved magnitude info from dB
sparamMag = 10.^((1/20).*rawDataDB(:,2:2:8));

% Convert retrieved phase info from degrees to radians
sparamPhase = rawDataDB(:,3:2:9)*(pi/180);

% Extract S-Parameter vectors
rawDataRI = sparamMag.*(cos(sparamPhase)+1i*sin(sparamPhase));
S11 = reshape(rawDataRI(:,1),1,1,numPoints);
S12 = reshape(rawDataRI(:,3),1,1,numPoints);
S21 = reshape(rawDataRI(:,2),1,1,numPoints);
S22 = reshape(rawDataRI(:,4),1,1,numPoints);

% Assemble into a 3D matrix to be consumed by the RF Toolbox
SParameter3Ddata = [S11 S12; S21 S22];

rfwrite(SParameter3Ddata,freqRange,'LNASParams.s2p');

% Create the sparameters object
hSParamData = sparameters(SParameter3Ddata,freqRange);
% Plot the magnitude v/s frequency plot for the S-parameters of the LNA
rfplot(hSParamData);
title('2-Port S-Parameter Measurement Values');