clear
clc

ph_offset=38.1; %FILL IN THIS NUMBER
mag_offset=0.56; %FILL IN THIS NUMBER
User='Research'; %FILL IN THIS STRING

%--------------------------------------------------------------------------
% Setup Connection with the VNA
%--------------------------------------------------------------------------

Boardindex = 0;                                                             % Board Index of GPIB board
Gadd = 16;                                                                  % GPIB address of VNA
vna = gpib('ni', Boardindex, Gadd);
set(vna, 'InputBufferSize', 100000);                                        % Large input buffer
set(vna, 'Timeout', 200);
fopen(vna);                                                                 % Open communication

%--------------------------------------------------------------------------
% VNA Measurement Settings
%--------------------------------------------------------------------------
fc = 1956.25;                                                               % carrier frequency, MHz
TI=8;                                                                       % Sweep time, seconds (0.01s - 86400s)
npoints=1601;                                                               % Numbers of data points (3,11,26,51,101,201,401,801,1601)
sweepnum=4;                                                                 % # of sweeps. Total sweep time = sweepnum * TI (in secs)
ifbw=300;                                                                    % IF bandwidth (10, 30, 100, 300, 1000, 3000, 3700, and 6000 Hz)
sweepnumber = 4;


for rho=0:0.2:1
    for k1=0:4:20
        k2=0;
        while (k1-k2) >= 0
            mR1=[];
            mR2=[];
            mI1=[];
            mI2=[];
            pR1=[];
            pR2=[];
            pI1=[];
            pI2=[];

            dlmwrite(['C:\Documents and Settings\' User '\Desktop\Pol2ChParameters.txt'],[k1;k2;rho],'newline','pc');
            dos(['C:\Documents and Settings\' User '\My Documents\Visual Studio 2008\Projects\SR5500_Pol_2Chan\Debug\SR5500_Pol_2Chan.exe']);
            
            %--------------------------------------------------------------------------
            % Initialize VNA 
            %--------------------------------------------------------------------------

            fprintf(vna,'PRES;');                                                       % preset the VNA
            pause(2);                                                                   % pause for two seconds, wait for VNA to reset
            fprintf(vna,'INSMNETA;');                                                   % Set to tuned receiver mode. Uncommand to set to default Network Analyzer mode.
            fprintf(vna,'DUACON;');                                                     % Set dual channel
            fprintf(vna,'SOUP ON;');                                                   % Turn off the source power such that power at port A won't be larger than port B.
            fprintf(vna,'POWE -14.0 DB;');
            
            % Channel 1
            fprintf(vna,'CHAN1;');                                                      % Select Channel 1 as active channel. This operation will restart the current sweep.
            fprintf(vna,['IFBW' num2str(ifbw) 'Hz;']);                                  % set IF bandwidth
            fprintf(vna, ['CWFREQ' num2str(fc) 'MHz;']);                                % set CW frequency
            fprintf(vna,['MEASA; POIN' num2str(npoints) '; SWET' num2str(TI) ';']);     % 'MEASA' means measure the input port A.
            %         For port B, use 'MEASB'. R, 'MEASR'. A/R, 'AR'. B/R, 'BR'. A/B, 'AB'.
            %         S11, 'S11' and so on for other S parameters.

            fprintf(vna,'LOGM;');                                                       % set data in Log Magnitude
            fprintf(vna, 'FORM4');                                                      % 'FORM4' sets the data in ASCii format

            % % Channel 2
            fprintf(vna,'CHAN2;');
            fprintf(vna,['IFBW' num2str(ifbw) 'Hz;']);
            fprintf(vna,['MEASB; POIN' num2str(npoints) '; SWET' num2str(TI) ';']);
            fprintf(vna, ['CWFREQ' num2str(fc) 'MHz;']);
            fprintf(vna,'LOGM;');                                                       % set data as measured phase
            fprintf(vna,'FORM4');                                                       % set array-data format 4, ASCII floating-point format


            %--------------------------------------------------------------------------
            % Start the Measurement, and repeat
            %--------------------------------------------------------------------------

            k=1;
            location_number = ['Pol_TestK1_' num2str(k1) 'K2_' num2str(k2) 'rho_' num2str(rho) ];      % the filename used to save the result data files within the working folder.
            filename1 = ['R:\etc\www\poldiv\' location_number 'AandB.mat'];
            filename2 = ['R:\etc\www\poldiv\' location_number 'AoverB.mat'];

            for rounds = 1:1:sweepnumber

                % Start a single Sweep
                fprintf(vna,'OPC?SING;');                                   % Restart a sweep for both channels and wait till the sweep is complete.
                % The sweep will be stopped after it's done.
                swp_complete = 0;
                while(swp_complete == 0)
                    swp_complete = fscanf(vna)
                end


                %----------------------------------------------------------------------
                % After sweep is complete, transfer both the magnitude and phase data
                %----------------------------------------------------------------------
                fprintf(vna,'CHAN1;');
                fprintf(vna,'LOGM;AUTO;');                                  % select log magnitude and auto scale. Here the 'AUTO' won't restart the sweep, so it's ok in terms of saving data and time.
                fprintf(vna, 'OUTPFORM');                                   % output the formatted trace data from the current active channel(for LOG MAG, in 'dB')
                d=fscanf(vna);                                              % download data to controller
                m=strread(d,'%f','delimiter',',');                          % %f=fixed point notation
                mR1=vertcat(mR1,m(1:2:end));                                % Extract real component
                mI1=vertcat(mI1,m(2:2:end));                                % Extract imaginary component
                fprintf(vna, 'PHAS; AUTO;');                                % select phase and auto scale. Here the 'AUTO' won't restart the sweep, so it's ok in terms of saving data and time.
                fprintf(vna, 'OUTPFORM');                                   % output the formatted trace data from the current active channel(for Phase, in 'degree')
                f=fscanf(vna);
                m=strread(f,'%f','delimiter',',');                          % %f=fixed point notation
                pR1=vertcat(pR1,m(1:2:end));                                % Extract real component
                pI1=vertcat(pI1,m(2:2:end));                                % Extract imaginary component

                fprintf(vna,'CHAN2;');
                fprintf(vna, 'PHAS; AUTO;');                                % select phase and auto scale. Here the 'AUTO' won't restart the sweep, so it's ok in terms of saving data and time.
                fprintf(vna, 'OUTPFORM');                                   % output the formatted trace data from the current active channel(for Phase, in 'degree')
                f=fscanf(vna);
                m=strread(f,'%f','delimiter',',');                          % %f=fixed point notation
                pR2=vertcat(pR2,m(1:2:end));                                % Extract real component
                pI2=vertcat(pI2,m(2:2:end));                                % Extract imaginary component
                fprintf(vna,'LOGM;AUTO;');                                  % select log magnitude and auto scale. Here the 'AUTO' won't restart the sweep, so it's ok in terms of saving data and time.
                fprintf(vna, 'OUTPFORM');                                   % output the formatted trace data from the current active channel(for LOG MAG, in 'dB')
                d=fscanf(vna);                                              % download data to controller
                m=strread(d,'%f','delimiter',',');                          % %f=fixed point notation
                mR2=vertcat(mR2,m(1:2:end));                                % Extract real component
                mI2=vertcat(mI2,m(2:2:end));                                % Extract imaginary component

            end

            %--------------------------------------------------------------------------
            % Antenna Calibration -- Apply Amplitude and Phase Offset
            %--------------------------------------------------------------------------

            pR12 = pR1 - pR2;
            % Unwrap if necessary
            for i =1:1:length(pR12)
                if(pR12(i)>180)
                    pab(i) = pR12(i)-360;
                elseif(pR12(i)<-180)
                    pab(i) = 360+pR12(i);
                else
                    pab(i) =  pR12(i);
                end
            end
            %

            % Aapply phase difference   % These values are obtained from the setting
            % both channels to static modualtion.
            % NOTE: THE OFFSET VALUE FOR PHASE HERE IS NOT DIRECTLY TAKEN FROM THE
            % MEASUREMENT OF A OVER B. IT MUST BE ADJUSTED! HERE THE DIFFERENCE BETWEEN
            % A OVER B'S OFFSET IS AROUND 31 DEGREES.

            % pab_cal = pab + 97.44;   % value for june 19
            % data_degvab = pab_cal;

            % pab_cal = pab - 110.6;   % value for june 20
            % data_degvab = pab_cal;

            pab_cal = pab + ph_offset;   % value for june 27
            data_degvab = pab_cal;

            % Apply magnitude difference
            % data_dBmvab = ((mR1-mR2)+0.6)';   % value for june 19
            % data_dBmvab = ((mR1-mR2)+0.2)';   % value for june 20
            data_dBmvab = ((mR1-mR2)+mag_offset)';   % value for june 27

            linear_magAB = 10.^(data_dBmvab./20);   % linear A/B in voltage

            % %----------------------------------------
            % % 10 here is the averaging number. 7 places to change
            avg_factor = 1;
            for i =1:1:floor(length(linear_magAB)/avg_factor)
                linear_pwr = [];
                linear_pwr = linear_magAB((i-1)*avg_factor+1:i*avg_factor);
                mean_linear_magAB(i) = mean(linear_pwr);
                gamma(i) = 180*atan(mean_linear_magAB(i))/pi;
            end


            for i =1:1:floor(length(data_degvab)/avg_factor)
                phase_deg = [];
                phase_deg = data_degvab((i-1)*avg_factor+1:i*avg_factor);
                delta(i) = mean(phase_deg);

            end


            % averaging with every avg_factor datasamples, and compute the (T,E) angle pair.
            for i =1:1:floor(length(data_degvab)/avg_factor)
                tau(i) = atan2(cos(delta(i)/180*pi),cot(2*gamma(i)/180*pi))/pi*180/2;
                E(i) = asin(sin(2*gamma(i)/180*pi).*sin(delta(i)/180*pi))/pi*180/2;
            end

            save(filename1);
            beep;

            %--------------------------------------------------------------------------
            % Take A and B measurements
            %--------------------------------------------------------------------------

            mR1=[];
            mR2=[];
            mI1=[];
            mI2=[];
            pR1=[];
            pR2=[];
            pI1=[];
            pI2=[];

            pause(2);   % pause for two seconds, wait for VNA to reset

            % Channel 1
            fprintf(vna,'CHAN1;');                                          % Select Channel 1 as active channel. This operation will restart the current sweep.
            fprintf(vna,['IFBW' num2str(ifbw) 'Hz;']);                      % set IF bandwidth
            fprintf(vna, ['CWFREQ' num2str(fc) 'MHz;']);                    % set CW frequency
            fprintf(vna,['AB; POIN' num2str(npoints) '; SWET' num2str(TI)...
                ';']);
            fprintf(vna,'LOGM;');                                           % set data in Log Magnitude
            fprintf(vna, 'FORM4');                                          % 'FORM4' sets the data in ASCii format

            % % Channel 2
            fprintf(vna,'CHAN2;');
            fprintf(vna,['IFBW' num2str(ifbw) 'Hz;']);
            fprintf(vna,['MEASA; POIN' num2str(npoints) '; SWET' ...
                num2str(TI) ';']);
            fprintf(vna, ['CWFREQ' num2str(fc) 'MHz;']);
            fprintf(vna,'LOGM;');                                           % set data as measured phase
            fprintf(vna,'FORM4');                                           % set array-data format 4, ASCII floating-point format


            %--------------------------------------------------------------------------
            % Starting the Measurement, and repeat it with sweep number of times
            %--------------------------------------------------------------------------

            for rounds = 1:1:sweepnumber

                % Start a single Sweep
                fprintf(vna,'OPC?SING;');                                   % Restart a sweep for both channel and wait till the sweep is complete.
                % The sweep will be stopped after it's done.
                a=5
                swp_complete = 0;
                while(swp_complete == 0)
                    swp_complete = fscanf(vna)
                end


                %----------------------------------------------------------------------
                % After sweep is complete, transfer both the magnitude and phase data
                %----------------------------------------------------------------------
                fprintf(vna,'CHAN1;');
                fprintf(vna,'LOGM;AUTO;');                                  % select log magnitude and auto scale. Here the 'AUTO' won't restart the sweep, so it's ok in terms of saving data and time.
                fprintf(vna, 'OUTPFORM');                                   % output the formatted trace data from the current active channel(for LOG MAG, in 'dB')
                d=fscanf(vna);                                              % download data to controller
                m=strread(d,'%f','delimiter',',');                          % %f=fixed point notation
                mR1=vertcat(mR1,m(1:2:end));                                % Extract real component
                mI1=vertcat(mR1,m(2:2:end));                                % Extract imaginary component
                fprintf(vna, 'PHAS; AUTO;');                                % select phase and auto scale. Here the 'AUTO' won't restart the sweep, so it's ok in terms of saving data and time.
                fprintf(vna, 'OUTPFORM');                                   % output the formatted trace data from the current active channel(for Phase, in 'degree')
                f=fscanf(vna);
                m=strread(f,'%f','delimiter',',');                          % %f=fixed point notation
                pR1=vertcat(pR1,m(1:2:end));                                % Extract real component
                pI1=vertcat(pI1,m(2:2:end));                                % Extract imaginary component

                fprintf(vna,'CHAN2;');
                fprintf(vna, 'PHAS; AUTO;');                                % select phase and auto scale. Here the 'AUTO' won't restart the sweep, so it's ok in terms of saving data and time.
                fprintf(vna, 'OUTPFORM');                                   % output the formatted trace data from the current active channel(for Phase, in 'degree')
                f=fscanf(vna);
                m=strread(f,'%f','delimiter',',');                          % %f=fixed point notation
                pR2=vertcat(pR2,m(1:2:end));                                 % Extract real component
                pI2=vertcat(pR2,m(2:2:end));                                 % Extract imaginary component

                fprintf(vna,'LOGM;AUTO;');                                  % select log magnitude and auto scale. Here the 'AUTO' won't restart the sweep, so it's ok in terms of saving data and time.
                fprintf(vna, 'OUTPFORM');                                   % output the formatted trace data from the current active channel(for LOG MAG, in 'dB')
                d=fscanf(vna);                                              % download data to controller
                m=strread(d,'%f','delimiter',',');                          % %f=fixed point notation
                mR2=vertcat(mR2,m(1:2:end));                                % Extract real component
                mI2=vertcat(mI2,m(2:2:end));                                % Extract imaginary component

            end

            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            % Antenna Calibration -- Apply Amplitude and Phase Offset
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %
            pR12 = pR1;
            % Unwrap if necessary
            for i =1:1:length(pR12)
                if(pR12(i)>180)
                    pab(i) = pR12(i)-360;
                elseif(pR12(i)<-180)
                    pab(i) = 360+pR12(i);
                else
                    pab(i) =  pR12(i);
                end
            end
            %

            % Aapply phase difference   % These values are obtained from the setting
            % both channels to static modualtion.

            % pab_cal = pab + 128.44;   % value for june 19
            % data_degvab = pab_cal;

            % pab_cal = pab -79.5;   % value for june 20
            % data_degvab = pab_cal;

            pab_cal = pab - ph_offset;   % value for june 20
            % data_degvab = pab_cal;

            % Apply magnitude difference
            % data_dBmvab = (mR1+0.6)';   % value for june 19
            data_dBmvab = (mR1+mag_offset)';   % value for june 20

            linear_magAB = 10.^(data_dBmvab./20);   % linear A/B in voltage

            % %----------------------------------------
            % % 10 here is the averaging number. 7 places to change
            % avg_factor = 1;
            for i =1:1:floor(length(linear_magAB)/avg_factor)
                linear_pwr = [];
                linear_pwr = linear_magAB((i-1)*avg_factor+1:i*avg_factor);
                mean_linear_magAB(i) = mean(linear_pwr);
                gamma(i) = 180*atan(mean_linear_magAB(i))/pi;
            end


            for i =1:1:floor(length(data_degvab)/avg_factor)
                phase_deg = [];
                phase_deg = data_degvab((i-1)*avg_factor+1:i*avg_factor);
                delta(i) = mean(phase_deg);

            end


            % averaging with every avg_factor datasamples, and compute the (T,E) angle pair.
            for i =1:1:floor(length(data_degvab)/avg_factor)

                tau(i) = atan2(cos(delta(i)/180*pi),cot(2*gamma(i)/180*pi))/pi*180/2;
                E(i) = asin(sin(2*gamma(i)/180*pi).*sin(delta(i)/180*pi))/pi*180/2;

            end

            save(filename2);

            k2=k2+4;
        end
    end
end

fclose(vna);

  
