analyzePassbandMeasurement_Pol_method2

 analyzePassbandMeasurement_Pol_method2(dirname,prefix)
 OR
 analyzePassbandMeasurement_Pol_method2(dirname,prefix,timeZone)

 This function analyzes the data from passband measurements. It reads in a
 standard set of files from the directory given in dirname. It also reads
 in the appropriate data from the archive. The prefix to the standard file
 names can be blank: ''

 The files expected in directory dirname are:
 [prefix]SweepPol.txt

 This function assumes that you have used the single-sweep method, with
 the splitter and attenuator tree in the cone of the dish.

 The function will produce some plots and write relevant ones as pdf files
 to directory dirname

 Oliver, 2 Feb 2012

0001 function analyzePassbandMeasurement_Pol_method2(varargin)
0002 % analyzePassbandMeasurement_Pol_method2(dirname,prefix)
0003 % OR
0004 % analyzePassbandMeasurement_Pol_method2(dirname,prefix,timeZone)
0005 %
0006 % This function analyzes the data from passband measurements. It reads in a
0007 % standard set of files from the directory given in dirname. It also reads
0008 % in the appropriate data from the archive. The prefix to the standard file
0009 % names can be blank: ''
0010 %
0011 % The files expected in directory dirname are:
0012 % [prefix]SweepPol.txt
0013 %
0014 % This function assumes that you have used the single-sweep method, with
0015 % the splitter and attenuator tree in the cone of the dish.
0016 %
0017 % The function will produce some plots and write relevant ones as pdf files
0018 % to directory dirname
0019 %
0020 % Oliver, 2 Feb 2012
0021 %
0022 
0023 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0024 % The default file names are:
0025 fn_I2J2 = 'Input_to_J2.s2p';
0026 fn_I2J3 = 'Input_to_J3.s2p';
0027 fn_I2PM = 'Input_to_power_meter.s2p';
0028 fn_PolSweep = 'SweepPol.txt';
0029 
0030 if nargin == 2
0031     % current OVRO time zone?
0032     OVRO_time_zone = 7; % hours behind UTC
0033     dirname = varargin{1};
0034     prefix = varargin{2};
0035 else
0036     dirname = varargin{1};
0037     prefix = varargin{2};
0038     OVRO_time_zone = abs(varargin{3});
0039 end
0040 
0041 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0042 
0043 if ~strcmp(dirname(end),'/')
0044     dirname = [dirname '/'];
0045 end
0046 
0047 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0048 % The cable loss:
0049 
0050 dI2J2 = importdata(fn_I2J2,' ',5);
0051 dI2J3 = importdata(fn_I2J3,' ',5);
0052 dI2PM = importdata(fn_I2PM,' ',5);
0053 f_VNA = dI2J2.data(:,1)/10^9; % in GHz
0054 
0055 SJ2I = dI2J2.data(:,4); % in dB
0056 SJ3I = dI2J3.data(:,4); % in dB
0057 SPMI = dI2PM.data(:,4); % in dB
0058 
0059 % figure
0060 % plot(f_VNA,SJ2I,'k-',f_VNA,SJ3I,'r-',f_VNA,SPMI,'b-')
0061 % xlabel('Frequency [GHz]')
0062 % ylabel('Amplitude [dB]')
0063 % title('Signal coupling structure')
0064 % legend('Input to LCP','Input to RCP','Input to power meter')
0065 
0066 
0067 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0068 % Now load the actual data from the telescope and the signal generator:
0069 
0070 % now add the prefix to the dirname string:
0071 dirname = [dirname prefix];
0072 
0073 % The signal generator data:
0074 dpol = importdata([dirname fn_PolSweep],'\t',1);
0075 tpol = dpol.data(:,1)+OVRO_time_zone/24;
0076 fpol = dpol.data(:,2);
0077 Ppol = dpol.data(:,3); % the set power
0078 % is Pmeas in units of dBm or mW? If is it in units of dBm, then
0079 % some of the measurements will be negative
0080 if any(dpol.data(:,4) < 0)
0081     Pmeas = dpol.data(:,4);
0082 else
0083     Pmeas= 10*log10(dpol.data(:,4)*1000); % the actual measured power
0084 end
0085 
0086 % The data to read in from the archive are:
0087 tstart = min(tpol)-2/60/24; % 2 minute buffer
0088 tend = max(tpol)+2/60/24;
0089 
0090 d = pipe_read(mjd2string(tstart),mjd2string(tend));
0091 
0092 % Calibrate the linearity of the measurements:
0093 d = calibrate_linearity(d);
0094 tdata = d.antenna0.receiver.utc;
0095 
0096 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0097 % Identify the events that we want to find the powers of:
0098 % We want to throw away these events:
0099 % - events at 5 GHz with set powers of either -28 or -25 dB
0100 % - events with only one frequency point (this means that their partner was
0101 % lost for some reason)
0102 
0103 % Ireject = [1,2,3,4,103,104,105,106,307,308,309,310,409,410,411,412]; % old way of doing it
0104 Ireject = (fpol == 5 & Ppol == -28) | (fpol == 5 & Ppol == -25);
0105 % now sort by frequency
0106 f_unique = unique(fpol);
0107 for k=1:length(f_unique)
0108     if sum(f_unique(k)==fpol)==1
0109         Ireject = Ireject | (fpol==f_unique(k));
0110     end
0111 end
0112 
0113 disp('analyzePassbandMeasurement_Pol_method2:: Rejected samples are:')
0114 disp('Time  Freq  Power')
0115 fprintf('%7.5f   %3.2f   %3.2f\n',[tpol(Ireject) fpol(Ireject) Ppol(Ireject)]')%'
0116 
0117 Ikeep = ones(size(tpol)) == 1;
0118 Ikeep(Ireject) = 0;
0119 tpol = tpol(Ikeep);
0120 fpol = fpol(Ikeep);
0121 Ppol = Ppol(Ikeep);
0122 Pmeas = Pmeas(Ikeep);
0123 
0124 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0125 % Assemble the total power data vector:
0126 swd = d.antenna0.receiver.switchData;
0127 TotalPower = swd(:,1:2:23)+swd(:,2:2:24);
0128 
0129 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0130 % The time keeping is not perfect, so match up the injected power times and
0131 % the data times. We are looking for a maximum correlation between a test
0132 % signal set we generate and the data.
0133 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0134 
0135 disp('analyzePassbandMeasurement_Pol_method2:: Identifying power injection events')
0136 tDel = 0.35; % the duration of the power events
0137 tOffsetList = -2.5:0.1:2.5;
0138 SumP = zeros(size(tOffsetList));
0139 for k=1:length(tOffsetList)
0140     % tdata: the data of the time samples
0141     % tpol: the list of time events
0142     IselectionPon = zeros(size(tdata));
0143     IselectionPoff = zeros(size(tdata));
0144     for m=1:2:(length(tpol)-1)
0145         % Every two events: the on events
0146         ttP = tpol(m); % start time of pol event
0147         IselectionPon = IselectionPon | (tdata-ttP) > tOffsetList(k)/60/60/24 & (tdata-ttP) < (tOffsetList(k)+tDel)/60/60/24;
0148         ttP = tpol(m+1); % start time of pol event
0149         IselectionPoff = IselectionPoff | (tdata-ttP) > tOffsetList(k)/60/60/24 & (tdata-ttP) < (tOffsetList(k)+tDel)/60/60/24;
0150     end
0151     
0152     % sum up all the data
0153     SumP(k) = sum(TotalPower(IselectionPon,1),1)-sum(TotalPower(IselectionPoff,1),1);
0154 end
0155 
0156 % what is the best offset for each?
0157 [y,IP] = max(SumP);
0158 offsetP = tOffsetList(IP);
0159 
0160 % and now read in the data with these offsets
0161 disp('analyzePassbandMeasurement_Pol_method2:: Reading in the data')
0162 IsP = zeros(size(tdata)); % index vector for plotting identified events
0163 
0164 swdMeans = zeros(size(tpol,1),size(swd,2));
0165 swdSigma = zeros(size(tpol,1),size(swd,2));
0166 
0167 for k=1:length(tpol)
0168     ttP = tpol(k); % start time of event
0169     % make an index vector which selects samples for the event:
0170     IselP = (tdata-ttP) > offsetP/60/60/24 & (tdata-ttP) < (offsetP+tDel)/60/60/24;
0171     % for plotting purposes
0172     IsP = IsP | IselP;
0173     % and now calculate the mean and stddev of the swd vector for the event:
0174     swdMeans(k,:) = mean(swd(IselP,:),1);
0175     swdSigma(k,:) = std(swd(IselP,:),1);
0176 end
0177 
0178 figure
0179 subplot(2,1,1)
0180 plot(tdata,TotalPower(:,1),'k-',tdata(IsP),TotalPower(IsP,1),'r.')
0181 xlabel('Time [mjd]')
0182 ylabel('Channel 1')
0183 legend('Data','Identified events')
0184 subplot(2,1,2)
0185 plot(tdata,TotalPower(:,11),'k-',tdata(IsP),TotalPower(IsP,11),'r.')
0186 xlabel('Time [mjd]')
0187 ylabel('Channel 11')
0188 legend('Data','Identified events')
0189 saveas(gcf,[dirname 'IdentifiedEvents.pdf'])
0190 
0191 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0192 % calculate the power on - power off matrix, after converting to units of
0193 % mW:
0194 mV2mW = 10^(-23/10)/3.6; % approximate conversion factor
0195 DU2mV = 2.5/10^4*1000/6.5; % convert DU to mV at the diode output
0196 swdMeans = swdMeans*DU2mV*mV2mW;
0197 swdSigma = swdSigma*DU2mV*mV2mW;
0198 
0199 swdDelta = swdMeans(1:2:(end-1),:) - swdMeans(2:2:end,:); % units of mW
0200 swdDeltaS = sqrt(swdSigma(1:2:(end-1),:).^2 + swdSigma(2:2:end,:).^2);
0201 
0202 finput = fpol(Pmeas~=-150);
0203 I1sky = swdDelta(:,2)+swdDelta(:,3); sI1sky = sqrt(swdDeltaS(:,2).^2+swdDeltaS(:,3).^2);
0204 I1load =swdDelta(:,1)+swdDelta(:,4); sI1load = sqrt(swdDeltaS(:,1).^2+swdDeltaS(:,4).^2);
0205 I2sky = swdDelta(:,22)+swdDelta(:,23); sI2sky = sqrt(swdDeltaS(:,22).^2+swdDeltaS(:,23).^2);
0206 I2load =swdDelta(:,21)+swdDelta(:,24); sI2load= sqrt(swdDeltaS(:,21).^2+swdDeltaS(:,24).^2);
0207 
0208 figure
0209 subplot(2,1,1)
0210 errorbar(unique(finput),I1load,sI1load,'k.')
0211 hold on
0212 errorbar(unique(finput),I1sky,sI1sky,'r.')
0213 xlabel('Frequency [GHz]')
0214 ylabel('I1 (on-off) [mW]')
0215 legend('Load state','Sky state')
0216 subplot(2,1,2)
0217 errorbar(unique(finput),I2load,sI2load,'k.')
0218 hold on
0219 errorbar(unique(finput),I2sky,sI2sky,'r.')
0220 xlabel('Frequency [GHz]')
0221 ylabel('I2 (on-off) [mW]')
0222 legend('Load state','Sky state')
0223 saveas(gcf,[dirname 'TotalIntensityResponses.pdf'])
0224 
0225 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0226 % Calculate the gains of the channels. How do we define gain?
0227 % Total intensity channel: sky/input
0228 % Polarization channel: response/input
0229 
0230 % first, interpolate the VNA data to the power meter data:
0231 
0232 % the correction, in dB, to be applied to the power meter data
0233 S_corr = interp1(f_VNA,SJ2I - SPMI,finput,'nearest');
0234 
0235 finjected = finput; % the same as the finput vector
0236 Pinjected = 10.^((Pmeas(Pmeas~=-150)+S_corr-30)/10); % the actualy signal generator power, in mW
0237 
0238 % demodulate the polarization channels:
0239 % Polarization combination set 1:
0240 Q1 =  1/2*(swdDelta(:,5) -swdDelta(:,6) + swdDelta(:,8) -swdDelta(:,7));
0241 sQ1= 1/2*sqrt(swdDeltaS(:,5).^2+swdDeltaS(:,6).^2+swdDeltaS(:,8).^2+swdDeltaS(:,7).^2);
0242 U1 = -1/2*(swdDelta(:,10)-swdDelta(:,9) + swdDelta(:,11)-swdDelta(:,12));
0243 sU1= 1/2*sqrt(swdDeltaS(:,10).^2+swdDeltaS(:,9).^2+swdDeltaS(:,11).^2+swdDeltaS(:,12).^2);
0244 % Polarization combination set 2:
0245 Q2 =  1/2*(swdDelta(:,13)-swdDelta(:,14)+swdDelta(:,16)-swdDelta(:,15));
0246 sQ2 = 1/2*sqrt(swdDeltaS(:,13).^2+swdDeltaS(:,14).^2+swdDeltaS(:,16).^2+swdDeltaS(:,15).^2);
0247 U2 = -1/2*(swdDelta(:,17)-swdDelta(:,18)+swdDelta(:,20)-swdDelta(:,19));
0248 sU2 = 1/2*sqrt(swdDeltaS(:,17).^2+swdDeltaS(:,18).^2+swdDeltaS(:,20).^2+swdDeltaS(:,19).^2);
0249 % Polarization combination set 3:
0250 Q3 = (Q1+Q2)/2; sQ3 = 1/2*sqrt(sQ1.^2+sQ2.^2);
0251 U3 = (U1+U2)/2; sU3 = 1/2*sqrt(sU1.^2+sU2.^2);
0252 
0253 GI1 = (I1sky-I1load)./Pinjected; sGI1 = sqrt(sI1sky.^2+sI1load.^2)./abs(Pinjected);
0254 GI2 = (I2sky-I2load)./Pinjected; sGI2 = sqrt(sI2sky.^2+sI2load.^2)./abs(Pinjected);
0255 GP = sqrt(Q3.^2+U3.^2)./Pinjected;
0256 sGP = 1./(sqrt(Pinjected.^2.*(Q3.^2+U3.^2))).*sqrt(Q3.^2.*sQ3.^2+U3.^2.*sU3.^2);
0257 
0258 G1 = GI1; sG1 = sGI1;
0259 G2 = Q1./Pinjected; sG2 = sQ1./Pinjected;
0260 G3 = U1./Pinjected; sG3 = sU1./Pinjected;
0261 G4 = Q2./Pinjected; sG4 = sQ2./Pinjected;
0262 G5 = U2./Pinjected; sG5 = sU2./Pinjected;
0263 G6 = GI2; sG6 = sGI2;
0264 
0265 
0266 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0267 % Interesting plots:
0268 % 1) Alpha values for I1 and I2 channels
0269 % 2) Passband plots and effective bandwidths
0270 % 3) Mueller matrix parameters
0271 
0272 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0273 % 1) Alpha values for I1 and I2 channels
0274 % the alpha values are
0275 a1 = (I1sky./I1load-1)./(I1sky./I1load+1);
0276 sa1 = 2./(I1sky+I1load).^2.*sqrt(I1load.^2.*sI1sky.^2+I1sky.^2.*sI1load.^2);
0277 a2 = (I2sky./I2load-1)./(I2sky./I2load+1);
0278 sa2 = 2./(I2sky+I2load).^2.*sqrt(I2load.^2.*sI2sky.^2+I2sky.^2.*sI2load.^2);
0279 
0280 % Now fit to a1 and a2.
0281 disp('analyzePassbandMeasurement_Pol_method2:: Fitting to alpha spectrum')
0282 % Try a number of initial values for the lag error:
0283 Loff = 30*(1:10);
0284 Lstart = [-fliplr(Loff) 0 Loff]/1000;
0285 
0286 % Only fit to the in-band point:
0287 Isel = finjected > 4.6 & finjected < 4.8 | finjected > 5.35 & finjected < 5.5;
0288 
0289 c1 = {finjected(Isel)*1e9 a1(Isel) sa1(Isel)};
0290 c2 = {finjected(Isel)*1e9 a2(Isel) sa2(Isel)};
0291 
0292 f1best = 1e10;
0293 f2best = 1e10;
0294 for k=1:length(Lstart)
0295     [xk1,fval1] = fminsearch(@(x) ssqerr(x,c1),[1;Lstart(k)]);
0296     [xk2,fval2] = fminsearch(@(x) ssqerr(x,c2),[1;Lstart(k)]);
0297         
0298     if fval1 < f1best
0299        x1 = xk1;
0300        f1best = fval1;
0301     end
0302     if fval2 < f2best
0303        x2 = xk2;
0304        f2best = fval2;
0305     end
0306 end
0307 
0308 fprintf('I1: Path error = %3.1f mm, Gain ratio = %3.2f\n',x1(2)*1000,x1(1))
0309 fprintf('I2: Path error = %3.1f mm, Gain ratio = %3.2f\n',x2(2)*1000,x2(1))
0310 a1fit = alpha(x1(1),x1(2),finjected*1e9,0);
0311 a2fit = alpha(x2(1),x2(2),finjected*1e9,0);
0312 
0313 figure
0314 hold on
0315 errorbar(finjected(Isel),a1(Isel),sa1(Isel),'k.')
0316 errorbar(finjected(Isel),a2(Isel),sa2(Isel),'r.')
0317 plot(finjected,a1fit,'k-',finjected,a2fit,'r-',...
0318     finjected,a1,'k-.',finjected,a2,'r-.')
0319 xlabel('Frequency [GHz]')
0320 ylabel('\alpha [goal = 1]')
0321 legend('\alpha_{I1}','\alpha_{I2}','Location','East')
0322 ylim([-1.05 1.05])
0323 
0324 sI1 = sprintf('I1: L=%3.1fmm, x=%3.2f',x1(2)*1000,x1(1));
0325 sI2 = sprintf('I2: L=%3.1fmm, x=%3.2f',x2(2)*1000,x2(1));
0326 
0327 title([sI1 '; ' sI2])
0328 
0329 saveas(gcf,[dirname 'alpha_Spectrum.pdf'])
0330 
0331 % and write the CW alpha values to disk:
0332 fn = [dirname 'alpha_cw.csv'];
0333 f1 = fopen(fn,'w');
0334 fprintf(f1,'# f [GHz], alpha_cw(I1), alpha_cw(I2), gain_I1 [lin], gain_I2 [lin], uncertainties\n');
0335 fprintf(f1,'%4.3f, %4.3f, %4.3f, %4.3f, %4.3f, %4.3f, %4.3f, %4.3f, %4.3f\n',[finjected a1 a2 GI1 GI2 sa1 sa2 sGI1 sGI2]');
0336 fclose(f1);
0337 
0338 
0339 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0340 % 2) Passband plots and effective bandwidths
0341 f = finjected;
0342 BWI1 = (max(f)-min(f))/length(f)*(sum(GI1))^2/(sum(GI1.^2));
0343 BWI2 = (max(f)-min(f))/length(f)*(sum(GI2))^2/(sum(GI2.^2));
0344 BWP = (max(f)-min(f))/length(f)*(sum(GP))^2/(sum(GP.^2));
0345 fcI1 = sum(f.*GI1)/(sum(GI1));
0346 fcI2 = sum(f.*GI2)/(sum(GI2));
0347 fcP = sum(f.*GP)/(sum(GP));
0348 
0349 fprintf('I1: f_centre = %4.3f, BW = %4.3f\n',fcI1,BWI1)
0350 fprintf('I2: f_centre = %4.3f, BW = %4.3f\n',fcI2,BWI2)
0351 fprintf('P: f_centre = %4.3f, BW = %4.3f\n',fcP,BWP)
0352 
0353 % Write the passbands to disk:
0354 fn = [dirname 'Passband.csv'];
0355 f1 = fopen(fn,'w');
0356 fprintf(f1,'# f [GHz], Gain I1, Gain Q1, Gain U1, Gain Q2, Gain U2, Gain I2, uncertainties\n');
0357 fprintf(f1,'%4.3f, %4.3f, %4.3f, %4.3f, %4.3f, %4.3f, %4.3f, %4.3f, %4.3f, %4.3f, %4.3f, %4.3f, %4.3f\n',...
0358     [f G1 G2 G3 G4 G5 G6 sG1 sG2 sG3 sG4 sG5 sG6]');
0359 fclose(f1);
0360 
0361 % Plot the data!
0362 ylims = [20 90];
0363 
0364 % figure
0365 % plot(finjected,10*log10(G1),'k-',...
0366 %     finjected,10*log10(G2),'r-',...
0367 %     finjected,10*log10(G3),'b-',...
0368 %     finjected,10*log10(G4),'r--',...
0369 %     finjected,10*log10(G5),'b--',...
0370 %     finjected,10*log10(G6),'k--')
0371 % ylim(ylims)
0372 % xlabel('Frequency [GHz]')
0373 % ylabel('Gain [dB]')
0374 % legend('Channel 1','Channel 2','Channel 3','Channel 4','Channel 5',...
0375 %     'Channel 6','Location','SouthEast')
0376 % title([sprintf('f_c [I1] = %4.3f GHz, BW = %4.3f GHz',fcI1,BWI1) ';' ...
0377 %     sprintf('f_c [I2] = %4.3f GHz, BW = %4.3f GHz',fcI2,BWI2) ';' ...
0378 %     sprintf('f_c [P] = %4.3f GHz, BW = %4.3f GHz',fcP,BWP)])
0379 % grid on
0380 % grid minor
0381 
0382 figure
0383 hold on
0384 plot(finjected,10*log10(abs(GI1)),'k-',...
0385     finjected,10*log10(abs(GP)),'r-',...
0386     finjected,10*log10(abs(GI2)),'b-')
0387 ylim(ylims)
0388 %set(gca,'YScale','log')
0389 xlabel('Frequency [GHz]')
0390 ylabel('Gain [dB]')
0391 legend('I1','P','I2','Location','SouthEast')
0392 title([sprintf('f_c [I1] = %4.3f GHz, BW = %4.3f GHz',fcI1,BWI1) ';' ...
0393     sprintf('f_c [I2] = %4.3f GHz, BW = %4.3f GHz',fcI2,BWI2) ';' ...
0394     sprintf('f_c [P] = %4.3f GHz, BW = %4.3f GHz',fcP,BWP)])
0395 grid on
0396 grid minor
0397 
0398 saveas(gcf,[dirname 'C-BASS_Passband.pdf'])
0399 
0400 
0401 % Plot the response of all the diodes:
0402 figure
0403 subplot(3,4,1)
0404 plot(finjected,10*log10(abs(swdDelta(:,1))./Pinjected),'k-',finjected,10*log10(abs(swdDelta(:,2))./Pinjected),'r-')
0405 ylabel('Amplitude [dB]')
0406 ylim(ylims)
0407 title('Diode 1')
0408 subplot(3,4,2)
0409 plot(finjected,10*log10(abs(swdDelta(:,3))./Pinjected),'k-',finjected,10*log10(abs(swdDelta(:,4))./Pinjected),'r-')
0410 ylim(ylims)
0411 title('Diode 2')
0412 subplot(3,4,3)
0413 plot(finjected,10*log10(abs(swdDelta(:,21))./Pinjected),'k-',finjected,10*log10(abs(swdDelta(:,22))./Pinjected),'r-')
0414 ylim(ylims)
0415 title('Diode 11')
0416 subplot(3,4,4)
0417 plot(finjected,10*log10(abs(swdDelta(:,23))./Pinjected),'k-',finjected,10*log10(abs(swdDelta(:,24))./Pinjected),'r-')
0418 ylim(ylims)
0419 title('Diode 12')
0420 
0421 subplot(3,4,5)
0422 plot(finjected,10*log10(abs(swdDelta(:,5))./Pinjected),'k-',finjected,10*log10(abs(swdDelta(:,6))./Pinjected),'r-')
0423 ylim(ylims)
0424 ylabel('Amplitude [dB]')
0425 title('Diode 3')
0426 subplot(3,4,6)
0427 plot(finjected,10*log10(abs(swdDelta(:,7))./Pinjected),'k-',finjected,10*log10(abs(swdDelta(:,8))./Pinjected),'r-')
0428 ylim(ylims)
0429 title('Diode 4')
0430 subplot(3,4,7)
0431 plot(finjected,10*log10(abs(swdDelta(:,13))./Pinjected),'k-',finjected,10*log10(abs(swdDelta(:,14))./Pinjected),'r-')
0432 ylim(ylims)
0433 title('Diode 7')
0434 subplot(3,4,8)
0435 plot(finjected,10*log10(abs(swdDelta(:,15))./Pinjected),'k-',finjected,10*log10(abs(swdDelta(:,16))./Pinjected),'r-')
0436 ylim(ylims)
0437 title('Diode 8')
0438 
0439 subplot(3,4,9)
0440 plot(finjected,10*log10(abs(swdDelta(:,9))./Pinjected),'k-',finjected,10*log10(abs(swdDelta(:,10))./Pinjected),'r-')
0441 ylim(ylims)
0442 ylabel('Amplitude [dB]')
0443 xlabel('Frequency [GHz]')
0444 title('Diode 5')
0445 subplot(3,4,10)
0446 plot(finjected,10*log10(abs(swdDelta(:,11))./Pinjected),'k-',finjected,10*log10(abs(swdDelta(:,12))./Pinjected),'r-')
0447 ylim(ylims)
0448 xlabel('Frequency [GHz]')
0449 title('Diode 6')
0450 subplot(3,4,11)
0451 plot(finjected,10*log10(abs(swdDelta(:,17))./Pinjected),'k-',finjected,10*log10(abs(swdDelta(:,18))./Pinjected),'r-')
0452 ylim(ylims)
0453 xlabel('Frequency [GHz]')
0454 title('Diode 9')
0455 subplot(3,4,12)
0456 plot(finjected,10*log10(abs(swdDelta(:,19))./Pinjected),'k-',finjected,10*log10(abs(swdDelta(:,20))./Pinjected),'r-')
0457 ylim(ylims)
0458 xlabel('Frequency [GHz]')
0459 legend('State 1','State 2','Location','South')
0460 title('Diode 10')
0461 
0462 saveas(gcf,[dirname 'C-BASS_Passband_Diode_Responses.pdf'])
0463 
0464 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0465 % 3) Mueller matrix parameters
0466 
0467 % Normalize the gains using a broadband source:
0468 G1n = G1/sum(G1); % I1
0469 sG1n = sqrt(1./sum(G1).^2.*(sG1.^2 + G1.^2.*sum(sG1.^2)));
0470 G2n = G2/sum(sqrt(G2.^2+G3.^2)); % Q1
0471 sG2n = sG2/sum(sqrt(G2.^2+G3.^2));
0472 G3n = G3/sum(sqrt(G2.^2+G3.^2)); % U1
0473 sG3n = sG3/sum(sqrt(G2.^2+G3.^2));
0474 G4n = G4/sum(sqrt(G4.^2+G5.^2)); % Q2
0475 sG4n = sG4/sum(sqrt(G4.^2+G5.^2));
0476 G5n = G5/sum(sqrt(G4.^2+G5.^2)); % U2
0477 sG5n = sG5/sum(sqrt(G4.^2+G5.^2));
0478 G6n = G6/sum(G6); % I2
0479 sG6n = sqrt(1./sum(G6).^2.*(sG6.^2 + G6.^2.*sum(sG6.^2)));
0480 
0481 GI1n = GI1/sum(GI1);
0482 sGI1n = sqrt(1./sum(GI1).^2.*(sGI1.^2 + GI1.^2.*sum(sGI1.^2)));
0483 GI2n = GI2/sum(GI2);
0484 sGI2n = sqrt(1./sum(GI2).^2.*(sGI2.^2 + GI2.^2.*sum(sGI2.^2)));
0485 GPn = GP/sum(GP);
0486 sGPn = sqrt(1./sum(GP).^2.*(sGP.^2 + GP.^2.*sum(sGP.^2)));
0487 
0488 Gnorm = max(max([GI1n GI2n GPn]));
0489 
0490 MII = (G1n+G6n)/(2*Gnorm); sMII = sqrt((sG1n.^2+sG2n.^2)./Gnorm^2);
0491 MP = GPn/Gnorm; sMP = sGPn/Gnorm;
0492 MQU = (G2n+G4n)/2/Gnorm; sMQU = sqrt((sG2n.^2+sG4n.^2)./(2*Gnorm)^2);
0493 MUU = (G3n+G5n)/2/Gnorm; sMUU = sqrt((sG3n.^2+sG5n.^2)./(2*Gnorm)^2);
0494 
0495 theta_P = atan2(MQU,MUU);
0496 stheta_P = sqrt((MUU.^2.*sMQU.^2 + MQU.^2.*sMUU.^2)./(MQU.^2+MUU.^2).^2);
0497 
0498 % we have limited information about the other channels. so synthesize them.
0499 MUQ = -MP.*sin(theta_P); sMUQ = sqrt(sin(theta_P).^2.*sMP.^2+MP.^2.*cos(theta_P).^2.*stheta_P.^2);
0500 MQQ = MP.*cos(theta_P); sMQQ = sqrt(cos(theta_P).^2.*sMP.^2+MP.^2.*sin(theta_P).^2.*stheta_P.^2);
0501 
0502 
0503 % Write the passbands to disk:
0504 fn = [dirname 'MuellerMatrixParams.csv'];
0505 f1 = fopen(fn,'w');
0506 fprintf(f1,'# f [GHz], MII, MQU, MUU, MUQ, MQQ, MP, theta_P [rad], uncertainties\n');
0507 fprintf(f1,'%4.3f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %4.3f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %10.9f, %4.3f\n',...
0508     [f (MII) (MQU) (MUU) (MUQ) (MQQ) (MP) theta_P (sMII) (sMQU) (sMUU) (sMUQ) (sMQQ) (sMP) stheta_P]');
0509 fclose(f1);
0510 
0511 ylims = [-30 0];
0512 figure
0513 subplot(2,2,1)
0514 plot(finjected,10*log10(abs(MII)),'k-')
0515 ylabel('Amplitude [dB]')
0516 title('M_{II}')
0517 ylim(ylims)
0518 subplot(2,2,2)
0519 plot(finjected,10*log10(abs(MP)),'k-')
0520 title('M_{P}')
0521 ylim(ylims)
0522 subplot(2,2,3)
0523 plot(finjected,(MQU),'k-')
0524 xlabel('Frequency [GHz]')
0525 ylabel('Amplitude [linear]')
0526 title('M_{QU}')
0527 ylim([-1 1])
0528 subplot(2,2,4)
0529 plot(finjected,(MUU),'k-')
0530 title('M_{UU}')
0531 xlabel('Frequency [GHz]')
0532 ylim([-1 1])
0533 saveas(gcf,[dirname 'Mueller_Matrix.pdf'])
0534 
0535 
0536 b1l = 4.51;
0537 b1h = 4.82;
0538 b2l = 5.33;
0539 b2h = 5.52;
0540 Ihigh = (finjected >= b1l & finjected <= b1h) | (finjected >= b2l & finjected <= b2h);
0541 % We have a pure U input. Output Q is then proportional to sin(w*L/c) and
0542 % the U output is proportional to cos(w*L/c).
0543 
0544 
0545 % and try fitting directly to theta_P
0546 c = {theta_P(Ihigh), stheta_P(Ihigh), finjected(Ihigh)};
0547 f_low = inf;
0548 % x_low = [0;0];
0549 x_low = [0];
0550 for k=1:10
0551 %     x0 = [randn(1,1)*2; randn(1,1)*0.5]; % randomly generate a starting point
0552     x0 = [randn(1,1)*0.5]; % randomly generate a starting point
0553     [x,fval] = fminsearch(@(x) ErrFun(x,c), x0);
0554     if fval < f_low
0555         x_low = x;
0556         f_low = fval;
0557     end
0558 end
0559 x = x_low;
0560 fval = f_low;
0561 xold = [0;x(1)];
0562 fprintf('Best fit is L = %4.1f mm, phi = %3.1f deg, and fval = %f \n',xold(1)*10,xold(2)*180/pi,fval)
0563 
0564 xlims = [4 6];
0565 ylims = [-1 1];
0566 figure
0567 subplot(3,2,1)
0568 hold on
0569 errorbar(finjected,MQU,sMQU,'k.')
0570 plot(finjected,Qfunc(x,finjected,MP),'k-')
0571 plot([b1l b1l],ylims,'b--',[b1h b1h],ylims,'b--',...
0572     [b2l b2l],ylims,'m--',[b2h b2h],ylims,'m--')
0573 xlabel('Frequency [GHz]')
0574 ylabel('M_{QU}')
0575 xlim(xlims)
0576 ylim(ylims)
0577 legend('Data','Fit')
0578 subplot(3,2,2)
0579 hold on
0580 errorbar(finjected,MUU,sMUU,'k.')
0581 plot(finjected,Ufunc(x,finjected,MP),'k-')
0582 plot([b1l b1l],ylims,'b--',[b1h b1h],ylims,'b--',...
0583     [b2l b2l],ylims,'m--',[b2h b2h],ylims,'m--')
0584 xlabel('Frequency [GHz]')
0585 ylabel('M_{UU}')
0586 xlim(xlims)
0587 ylim(ylims)
0588 legend('Data','Fit')
0589 
0590 subplot(3,2,3:4)
0591 hold on
0592 ylims = [-180 180];
0593 errorbar(finjected(~Ihigh),theta_P(~Ihigh)*180/pi,stheta_P(~Ihigh)*180/pi,'k.')
0594 errorbar(finjected(Ihigh),theta_P(Ihigh)*180/pi,stheta_P(Ihigh)*180/pi,'r.')
0595 plot(finjected,phs(x,finjected)*180/pi,'r-')
0596 plot([b1l b1l],ylims,'b--',[b1h b1h],ylims,'b--',...
0597     [b2l b2l],ylims,'m--',[b2h b2h],ylims,'m--')
0598 ylim(ylims)
0599 xlim(xlims)
0600 xlabel('Frequency [GHz]')
0601 ylabel('Polarization angle [deg]')
0602 title(sprintf('Best-fit: Path error = %3.1f mm, phase error = %3.1f deg',xold(1)*10,xold(2)*180/pi))
0603 legend('All data','In-band data','Fit','Location','NorthWest')
0604 
0605 subplot(3,2,5:6)
0606 hold on
0607 ylims = [-30 30];
0608 errorbar(finjected(~Ihigh),(theta_P(~Ihigh)-phs(x,finjected(~Ihigh)))*180/pi,stheta_P(~Ihigh)*180/pi,'k.')
0609 errorbar(finjected(Ihigh),(theta_P(Ihigh)-phs(x,finjected(Ihigh)))*180/pi,stheta_P(Ihigh)*180/pi,'r.')
0610 plot(xlims,[0 0],'k--')
0611 plot([b1l b1l],ylims,'b--',[b1h b1h],ylims,'b--',...
0612     [b2l b2l],ylims,'m--',[b2h b2h],ylims,'m--')
0613 ylim(ylims)
0614 xlim(xlims)
0615 xlabel('Frequency [GHz]')
0616 ylabel('Residual [deg]')
0617 legend('Data-Fit','Location','NorthWest')
0618 
0619 saveas(gcf,[dirname 'Linear_Mueller_Fit.pdf'])
0620 
0621 disp('analyzePassbandMeasurement_Pol_method2:: Finished')
0622 end
0623 
0624 function y = phs(x,f)
0625 % L is in cm
0626 % y = atan(tan(  2*pi*f*10^9/3e8*(x(1)/100) ));
0627 
0628 % y = atan2(sin(  2*pi*f*10^9/3e8*(x(1)/100) +x(2) ),cos(  2*pi*f*10^9/3e8*(x(1)/100)+x(2) ));
0629 y = atan2(sin(  x ),cos(  x ));
0630 
0631 end
0632 function y = Qfunc(x,f,MP)
0633 % y = -MP.*sin( -2*pi*f*10^9/3e8*(x/100) );
0634 y = MP.*sin(phs(x,f));
0635 end
0636 function y = Ufunc(x,f,MP)
0637 % y = -MP.*cos( -2*pi*f*10^9/3e8*(x/100) );
0638 y = MP.*cos(phs(x,f));
0639 end
0640 
0641 function y = ErrFun(x,c)
0642 theta_P = c{1};
0643 stheta_P = c{2};
0644 f = c{3};
0645 y = sum( (( phs(x,f) - theta_P).^2) ./ (stheta_P.^2) );
0646 end
0647 
0648 function chi2 = ssqerr(x,c)
0649 % Function to fit the alpha value curve to:
0650 v = c{1}; % frequencies of the alpha points in Hz
0651 a = c{2}; % the alpha values
0652 sa = c{3}; % the errors in the alpha values
0653 
0654 % x(1) is ratio of voltage gains, x(2) is path error
0655 afit = alpha(x(1),x(2),v,0);
0656 % we want to minimize the chi^2:
0657 chi2 = sum((afit-a).^2./sa.^2);
0658 end
0659 
0660 function err = ssqerrP(x,c)
0661 % Function to fit the alpha value curve to:
0662 v = c{1}; % frequencies of the alpha points in Hz
0663 a = c{2}; % the alpha values
0664 
0665 % x(1) is ratio of voltage gains, x(2) is path error, x(3) is phase offset
0666 afit = alpha(x(1),x(2),v,x(3));
0667 err = sum((afit-a).^2);
0668 end
0669 
0670 function err = ssqerrSim(x,c)
0671 % Simultaneous fit to left, right, and pol channels
0672 v = c{1};
0673 aR = c{2};
0674 aL = c{3};
0675 D3 = c{4};
0676 
0677 % x(1) = xR
0678 % x(2) = LR
0679 % x(3) = xL
0680 % x(4) = LL
0681 % x(5) = xP
0682 
0683 aRfit = alpha(x(1),x(2),v,0);
0684 aLfit = alpha(x(3),x(4),v,0);
0685 D3fit = D3ratio(x(5),x(3),x(1),x(4),x(2),v);
0686 
0687 err = sum((aR-aRfit).^2+(aL-aLfit).^2+(D3-D3fit).^2);
0688 
0689 end
0690 
0691 function y = alpha(x,L,v,toff)
0692 % x = |g1|/|g2|
0693 % L = path length difference in m
0694 % v = frequency in Hz
0695 c = 3e8;
0696 Lc = L/c;
0697 
0698 y = -2*x*cos(2*pi*v*Lc+toff)/(1+x^2);
0699 end
0700 
0701 function y = D3ratio(xP,xL,xR,LL,LR,v)
0702 % the ratio of the diode 3 output powers
0703 y = xP^2*(1+xR^2+2*xR*cos(2*pi*v*LR/3e8))./(1+xL^2+2*xL*cos(2*pi*v*LL/3e8));
0704 
0705 end