%This program uses two signal generators and a spectrum analyzer to get the
%OIP2.  
%Dependencies:
%   SigGen_init.m
%   SigGen_setlevel.m
%   sensitiveReadSpectrum.m
%Notes:
%You can either have user-input for the GPIB addresses, frequencies, and 
% resolution bandwidth or hard-code them, depending on what's commented.  
% Right now, they're hard-coded.
%The values for OIP calculated aren't right, they're just ballpark
%estimates

close all %Delivers errors on first call - comment this line out the first time program is run
clear all

%comment one of the following two sections:
%hard-coded initial info (i.e. addresses, frequencies):
analyzer_address = 20;      %#ok<NASGU>     (integers ...)
RFgen1_address = 28;        %#ok<NASGU>
RFgen2_address = 30;        %#ok<NASGU>
f1 = 0.5;                     %#ok<NASGU>     (GHz)
f2 = 0.505;                 %#ok<NASGU>     (GHz)
resBW = 30;                %#ok<NASGU>     (Hz) 
secondsweepcount = 5;
fundsweepcount = 3;
span = 500;

%user-input initial info:
%analyzer_address = input('Spectrum analyzer GPIB address (1-30):     ');
%RFgen1_address = input('RF generator 1 GPIB address (1-30):     ');
%RFgen2_address = input('RF generator 2 GPIB address (1-30):     ');
%f1 = input('frequency 1 (GHz):     ');
%f2 = input('frequency 2 (GHz):     ');
%resBW = input('resolution bandwidth (kHz):     ');

%error-check initial info (some of these ranges are made-up - they're just
%there to make sure the program doesn't try to set the frequency to
%anything other than a number
while (analyzer_address<0)||(analyzer_address>30)
    analyzer_address = input('Address invalid.  Please enter an integer 1-30:     ');
end
while ((analyzer_address<0)||(analyzer_address>30)||(RFgen1_address==analyzer_address))
    analyzer_address = input('Address invalid.  Please enter an integer 1-30:     ');
end
while ((analyzer_address<0)||(analyzer_address>30)||(RFgen2_address==analyzer_address)||(RFgen2_address==RFgen1_address))
    analyzer_address = input('Address invalid.  Please enter an integer 1-30:     ');
end
while (f1<0)||(f1>50)
    f1 = input('f1 invalid.  Retry:     ');
end
while (f2<0)||(f2>50)
    f2 = input('f2 invalid.  Retry:     ');
end

%get start power and number of power level steps from user:
start_power1 = input('start power level at f1 GHz (bottom of range, in dBm):     ');
start_power2 = input('start power level at f2 GHz (bottom of range, in dBm):     ');
num_steps = input('number of (power level) steps:     ');

%initialize RF generators (function sets these in MHz; here they're in GHz, 
%hence the *1e3)
SigGen_init(RFgen1_address, f1*1e3);
SigGen_init(RFgen2_address, f2*1e3);

%make RF magnitude vector:
RFmag_vec1 = start_power1:1:(start_power1+num_steps-1);
RFmag_vec2 = start_power2:1:(start_power2+num_steps-1);

cont = 1;                       %control variable: 1: continue; 2: quit
num_measurements = 0;           %start number of measurements counter
OIP = struct('f1', 1);

while cont == 1
    
    num_measurements = num_measurements+1;      %increment iteration counter
    
    %set RF magnitudes and notify user
    SigGen_setlevel(RFgen1_address, RFmag_vec1(num_measurements));
    SigGen_setlevel(RFgen2_address, RFmag_vec2(num_measurements));
    
    %Get first fundamental trace data - ignore spellcheck error "str2double
    %is faster - that's great, but it also doesn't translate properly
    fund1(num_measurements) = sensitiveReadSpectrum(analyzer_address, f1, span, resBW, fundsweepcount); %#ok<AGROW>
    sprintf('read at %f done', f1)
    
    %Get second fundamental trace data
    fund2(num_measurements) = sensitiveReadSpectrum(analyzer_address, f2, span, resBW, fundsweepcount); %#ok<AGROW>
    sprintf('read at %f done', f2)
    
    %Get first second intermod trace data
    CFreq = abs(f1-f2);
    second_imod1(num_measurements) = sensitiveReadSpectrum(analyzer_address, CFreq, span, resBW, secondsweepcount);
    sprintf('read at %f done', CFreq)
    
    %Get second second intermod trace data
    CFreq = f1+f2;
   	second_imod2(num_measurements) = sensitiveReadSpectrum(analyzer_address, CFreq, span, resBW, secondsweepcount);
    sprintf('read at %f done', CFreq)
    
%     %store data
%     OIP3(num_measurements) = OIP;
%     clear OIP;
    
    if(num_measurements==num_steps)
        cont = 2;
    end
end

%make data plottable
RFmag_vec = (RFmag_vec1+RFmag_vec2)/2;
for i = 1:num_measurements
    b11(i) = fund1(i)-RFmag_vec(i);
    b12(i) = fund2(i)-RFmag_vec(i);
    b21(i) = second_imod1(i)-2*RFmag_vec(i);
    b22(i) = second_imod2(i)-2*RFmag_vec(i);
end

%calculate OIP3s
B11 = sum(b11)/length(b11);
B12 = sum(b12)/length(b12);
B21 = sum(b21)/length(b21);
B22 = sum(b22)/length(b22);

OIP21 = 2*B11-B21;
OIP22 = 2*B12-B22;

f1slope = polyfit(RFmag_vec1, fund1, 1);
f2slope = polyfit(RFmag_vec2, fund2, 1);
imod1slope = polyfit(RFmag_vec2, second_imod1, 1);
imod2slope = polyfit(RFmag_vec2, second_imod2, 1);

%plot data
figure(2)
plot(RFmag_vec, fund1, 'bx', RFmag_vec, fund2, 'rx')
hold on
plot(RFmag_vec, second_imod1, 'go', RFmag_vec, second_imod2, 'co')
hold off
title(strcat('OIP2 = ', num2str(OIP21), ' or ', num2str(OIP22), ' dBm'));
key1 = strcat(num2str(f1), ' GHz');
key2 = strcat(num2str(f2), ' GHz');
key3 = strcat(num2str(2*f1-f2), ' GHz');
key4 = strcat(num2str(2*f2-f1), ' GHz');
legend(key1, key2, key3, key4, 'Location', 'Best');

%spit out slopes
f1slope(1)
f2slope(1)
imod1slope(1)
imod2slope(1)

%save everything
uisave();