analyzePassbandMeasurement_Pol

 analyzePassbandMeasurement_Pol(dirname,prefix)
 OR
 analyzePassbandMeasurement_Pol(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]CableTransmission.Txt
 [prefix]0dBmSweep.Txt
 [prefix]InjectedPower.txt
 [prefix]SweepPol.txt

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

 Oliver, 18 Nov 2011

0001 function analyzePassbandMeasurement_Pol(varargin)
0002 % analyzePassbandMeasurement_Pol(dirname,prefix)
0003 % OR
0004 % analyzePassbandMeasurement_Pol(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]CableTransmission.Txt
0013 % [prefix]0dBmSweep.Txt
0014 % [prefix]InjectedPower.txt
0015 % [prefix]SweepPol.txt
0016 %
0017 % The function will produce some plots and write relevant ones as pdf files
0018 % to directory dirname
0019 %
0020 % Oliver, 18 Nov 2011
0021 %
0022 
0023 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0024 % The default file names are:
0025 fn_cable = 'CableTransmission.txt';
0026 fn_0dBmSweep = '0dBmSweep.txt';
0027 fn_InjectedPower = 'InjectedPower.txt';
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 % now add the prefix to the dirname string:
0048 dirname = [dirname prefix];
0049 
0050 % The cable loss:
0051 dc = importdata([dirname fn_cable],'\t',1);
0052 fcable = dc.data(:,2);
0053 Pcable = dc.data(:,4)-dc.data(:,3);
0054 
0055 dcn = importdata([dirname fn_0dBmSweep],'\t',1);
0056 fnocable = dcn.data(:,2);
0057 Pnocable = dcn.data(:,4)-dcn.data(:,3);
0058 
0059 Gcable = Pcable-Pnocable;
0060 
0061 figure
0062 plot(fcable,Gcable,'k-')
0063 xlabel('Frequency [GHz]')
0064 ylabel('Cable transmission [dB]')
0065 
0066 % The actual power being put out by the sig gen:
0067 dsg = importdata([dirname fn_InjectedPower],'\t',1);
0068 
0069 fsiggen = dsg.data(:,2);
0070 SetPower = dsg.data(:,3);
0071 MeasuredPower = dsg.data(:,4);
0072 
0073 figure
0074 subplot(2,1,1)
0075 plot(fsiggen,SetPower,'k.',fsiggen,MeasuredPower,'r.')
0076 xlabel('Frequency [GHz]')
0077 ylabel('Sig gen power [dB]')
0078 legend('Set Power','Measured Power')
0079 subplot(2,1,2)
0080 plot(fsiggen,SetPower-MeasuredPower,'k.')
0081 xlabel('Frequency [GHz]')
0082 ylabel('Sig gen power error [dB]')
0083 
0084 saveas(gcf,[dirname 'SignalGeneratorBehavior.pdf'])
0085 
0086 % So, the actual power arriving at the input to the cryostat is:
0087 Pin = MeasuredPower+Gcable-30; % assume a flat -30dB of coupling
0088 
0089 figure
0090 plot(fcable,Pin,'k-')
0091 xlabel('Frequency')
0092 ylabel('Power in [dBm]')
0093 
0094 saveas(gcf,[dirname 'PowerInProfile.pdf'])
0095 
0096 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0097 % Now load the actual data from the telescope and the signal generator:
0098 
0099 % The signal generator data:
0100 dpol = importdata([dirname fn_PolSweep],'\t',1);
0101 tpol = dpol.data(:,1)+OVRO_time_zone/24;
0102 fpol = dpol.data(:,2);
0103 Ppol = dpol.data(:,3);
0104 
0105 % The data to read in from the archive are:
0106 tstart = min(tpol)-2/60/24; % 2 minute buffer
0107 tend = max(tpol)+2/60/24;
0108 
0109 d = pipe_read(mjd2string(tstart),mjd2string(tend));
0110 
0111 % Calibrate the linearity of the measurements:
0112 d = calibrate_linearity(d);
0113 tdata = d.antenna0.receiver.utc;
0114 
0115 % The total power received at each diode:
0116 Pdata = d.antenna0.receiver.switchData(:,1:2:23) + d.antenna0.receiver.switchData(:,2:2:24);
0117 
0118 % For measuring the alpha value:
0119 P1 = d.antenna0.receiver.switchData(:,1)+d.antenna0.receiver.switchData(:,4); % "sky"
0120 P2 = d.antenna0.receiver.switchData(:,2)+d.antenna0.receiver.switchData(:,3); % "load"
0121 
0122 P11 = d.antenna0.receiver.switchData(:,21)+d.antenna0.receiver.switchData(:,24); % "sky"
0123 P12 = d.antenna0.receiver.switchData(:,22)+d.antenna0.receiver.switchData(:,23); % "load"
0124 
0125 % there are some throw-away events at the beginning and the end
0126 tlist = 5:(length(tleft)-4);
0127 
0128 % Declare some vectors we want the measured diode response to the injection
0129 % into LCP, and then into RCP.
0130 
0131 PmeasOL = zeros(length(tlist),12);
0132 PmeasOR = zeros(length(tlist),12);
0133 PmeasOP = zeros(length(tlist),12);
0134 
0135 P1L = zeros(length(tlist),1);
0136 P1R = zeros(length(tlist),1);
0137 P2L = zeros(length(tlist),1);
0138 P2R = zeros(length(tlist),1);
0139 
0140 % P3L = zeros(length(tlist),1);
0141 % P3R = zeros(length(tlist),1);
0142 % P4L = zeros(length(tlist),1);
0143 % P4R = zeros(length(tlist),1);
0144 %
0145 % P5L = zeros(length(tlist),1);
0146 % P5R = zeros(length(tlist),1);
0147 % P6L = zeros(length(tlist),1);
0148 % P6R = zeros(length(tlist),1);
0149 
0150 P11L = zeros(length(tlist),1);
0151 P11R = zeros(length(tlist),1);
0152 P12L = zeros(length(tlist),1);
0153 P12R = zeros(length(tlist),1);
0154 
0155 
0156 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0157 % The time keeping is not perfect, so match up the injected power times and
0158 % the data times. We are looking for a maximum correlation between a test
0159 % signal set we generate and the data.
0160 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0161 
0162 disp('analyzePassbandMeasurement_Pol:: Identifying power injection events')
0163 tDel = 0.35; % the duration of the power events
0164 tOffsetList = -1.5:0.05:1.5;
0165 for k=1:length(tOffsetList)
0166     % tdata: the data of the time samples
0167     % tlist: the list of time events
0168     IselectionL = zeros(size(tdata));
0169     IselectionR = zeros(size(tdata));
0170     IselectionP = zeros(size(tdata));
0171     for m=1:2:(length(tlist)-1)
0172         % Every two events: the on events
0173         ttL = tleft(tlist(m)); % start time of left event
0174         ttR = tright(tlist(m)); % start time of right event
0175         ttP = tpol(tlist(m)); % start time of pol event
0176         IselectionL = IselectionL | (tdata-ttL) > tOffsetList(k)/60/60/24 & (tdata-ttL) < (tOffsetList(k)+tDel)/60/60/24;
0177         IselectionR = IselectionR | (tdata-ttR) > tOffsetList(k)/60/60/24 & (tdata-ttR) < (tOffsetList(k)+tDel)/60/60/24;
0178         IselectionP = IselectionP | (tdata-ttP) > tOffsetList(k)/60/60/24 & (tdata-ttP) < (tOffsetList(k)+tDel)/60/60/24;
0179     end
0180     
0181     % sum up all the data
0182     SumR(k) = sum(Pdata(IselectionR,1),1);
0183     SumL(k) = sum(Pdata(IselectionL,11),1);
0184     SumP(k) = sum(Pdata(IselectionP,5),1);
0185 end
0186 
0187 % what is the best offset for each?
0188 [y,IR] = max(SumR);
0189 [y,IL] = max(SumL);
0190 [y,IP] = max(SumP);
0191 
0192 offsetR = tOffsetList(IR);
0193 offsetL = tOffsetList(IL);
0194 offsetP = tOffsetList(IP);
0195 
0196 disp('analyzePassbandMeasurement_Pol:: Reading in the data')
0197 % and now read in the data with these offsets
0198 IsL = zeros(size(tdata));
0199 IsR = zeros(size(tdata));
0200 IsP = zeros(size(tdata));
0201 for k=1:length(tlist)
0202     ttL = tleft(tlist(k)); % start time of event
0203     ttR = tright(tlist(k)); % start time of event
0204     
0205     IselL = (tdata-ttL) > offsetL/60/60/24 & (tdata-ttL) < (offsetL+tDel)/60/60/24;
0206     IselR = (tdata-ttR) > offsetR/60/60/24 & (tdata-ttR) < (offsetR+tDel)/60/60/24;
0207     IselP = (tdata-ttP) > offsetP/60/60/24 & (tdata-ttP) < (offsetP+tDel)/60/60/24;
0208     
0209     % for plotting purposes
0210     IsL = IsL | IselL;
0211     IsR = IsR | IselR;
0212     IsP = IsP | IselP;
0213     
0214     PmeasOL(k,:) = mean(Pdata(IselL,:),1);
0215     PmeasOR(k,:) = mean(Pdata(IselR,:),1);
0216     PmeasOP(k,:) = mean(Pdata(IselP,:),1);
0217     
0218     P1L(k,:) = mean(P1(IselL,:),1);
0219     P1R(k,:) = mean(P1(IselR,:),1);
0220     
0221     P2L(k,:) = mean(P2(IselL,:),1);
0222     P2R(k,:) = mean(P2(IselR,:),1);
0223     
0224     P11L(k,:) = mean(P11(IselL,:),1);
0225     P11R(k,:) = mean(P11(IselR,:),1);
0226     
0227     P12L(k,:) = mean(P12(IselL,:),1);
0228     P12R(k,:) = mean(P12(IselR,:),1);
0229     
0230 end
0231 
0232 figure
0233 subplot(2,1,1)
0234 plot(tdata,Pdata(:,1),'k-',tdata(IsR),Pdata(IsR,1),'r.',tdata(IsP),Pdata(IsP,1),'g.')
0235 xlabel('Time [mjd]')
0236 ylabel('Channel 1')
0237 legend('Data','Identified events','Pol')
0238 subplot(2,1,2)
0239 plot(tdata,Pdata(:,11),'k-',tdata(IsL),Pdata(IsL,11),'r.',tdata(IsP),Pdata(IsP,11),'g.')
0240 xlabel('Time [mjd]')
0241 ylabel('Channel 11')
0242 legend('Data','Identified events','Pol')
0243 saveas(gcf,[dirname 'IdentifiedEvents.pdf'])
0244 
0245 % Now process the data to extract passbands: (on - off)
0246 PmeasuredL = PmeasOL(1:2:end-1,:) - PmeasOL(2:2:end,:);
0247 PmeasuredR = PmeasOR(1:2:end-1,:) - PmeasOR(2:2:end,:);
0248 PmeasuredP = PmeasOP(1:2:end-1,:) - PmeasOP(2:2:end,:);
0249 
0250 dP1L = P1L(1:2:end-1,:) - P1L(2:2:end,:);
0251 dP1R = P1R(1:2:end-1,:) - P1R(2:2:end,:);
0252 
0253 dP2L = P2L(1:2:end-1,:) - P2L(2:2:end,:);
0254 dP2R = P2R(1:2:end-1,:) - P2R(2:2:end,:);
0255 
0256 dP11L = P11L(1:2:end-1,:) - P11L(2:2:end,:);
0257 dP11R = P11R(1:2:end-1,:) - P11R(2:2:end,:);
0258 
0259 dP12L = P12L(1:2:end-1,:) - P12L(2:2:end,:);
0260 dP12R = P12R(1:2:end-1,:) - P12R(2:2:end,:);
0261 
0262 Pnom = Pleft(tlist);
0263 Pnominal = Pnom(1:2:end-1);
0264 
0265 fE = fleft(1:2:end-1); % the frequency of the event.
0266 
0267 % Now exclude measurements 51, 52, 153, and 154
0268 Iselect = [1:50 53:152 155:length(Pnominal)];
0269 
0270 mV2mW = 10^(-23/10)/3.6; % approximate conversion factor
0271 DU2mV = 2.5/10^4*1000/6.5; % convert DU to mV at the diode output
0272 dPL = PmeasuredL(Iselect,:)*DU2mV*mV2mW; % diode input in mW
0273 dPR = PmeasuredR(Iselect,:)*DU2mV*mV2mW; % diode input in mW
0274 dPP = PmeasuredP(Iselect,:)*DU2mV*mV2mW;
0275 
0276 dP1L = dP1L(Iselect,:)*DU2mV*mV2mW;
0277 dP2L = dP2L(Iselect,:)*DU2mV*mV2mW;
0278 dP1R = dP1R(Iselect,:)*DU2mV*mV2mW;
0279 dP2R = dP2R(Iselect,:)*DU2mV*mV2mW;
0280 
0281 dP11L = dP11L(Iselect,:)*DU2mV*mV2mW;
0282 dP12L = dP12L(Iselect,:)*DU2mV*mV2mW;
0283 dP11R = dP11R(Iselect,:)*DU2mV*mV2mW;
0284 dP12R = dP12R(Iselect,:)*DU2mV*mV2mW;
0285 
0286 Psg = 10.^(Pin/10); % the actualy signal generator power, in mW
0287 
0288 % And hence the gains are:
0289 GL = (abs(dPL./repmat(Psg(:),1,12)));
0290 GR = (abs(dPR./repmat(Psg(:),1,12)));
0291 GPol = (abs(dPP./repmat(Psg(:),1,12)));
0292 
0293 % If you want to exclude a particular data point from the plots because of
0294 % RFI, you can use the Ip vector to do so.
0295 [y,Ir] = min(abs(fcable-4.78));
0296 
0297 Ip = true(size(fcable));
0298  Ip(50) = false;
0299  Ip(150) = false;
0300 
0301 % figure
0302 % plot(fcable(Ip),GdL(Ip,2),'k-',fcable(Ip),GdR(Ip,2),'r-')
0303 % xlabel('Frequency [GHz]')
0304 % ylabel('Gain at output 2')
0305 % legend('From L','From R')
0306 
0307 % the alpha values are
0308 % aL1 = (dP1L./dP2L-1)./(dP1L./dP2L+1);
0309 aR1 = (dP1R./dP2R-1)./(dP1R./dP2R+1);
0310 %aR1 = (dP2R./dP1R-1)./(dP2R./dP1R+1);
0311 
0312 aL2 = (dP11L./dP12L-1)./(dP11L./dP12L+1);
0313 %aL2 = (dP12L./dP11L-1)./(dP12L./dP11L+1);
0314 % aR2 = (dP11R./dP12R-1)./(dP11R./dP12R+1);
0315 
0316 D3 = dPR(:,3)./dPL(:,3);
0317 
0318 gainR = (GR(:,1)+GR(:,2));
0319 gainL = (GL(:,11)+GL(:,12));
0320 % What is the average alpha value?
0321 aRm = sum(gainR.*aR1)/sum(gainR);
0322 aLm = sum(gainL.*aL2)/sum(gainL);
0323 
0324 Isel = gainR>max(gainR/8) & Ip(:);
0325 
0326 f = fcable(Ip);
0327 GI1 = (GR((Ip),1)+GR((Ip),2)+GL((Ip),1)+GL((Ip),2))/4;
0328 GI2 = (GL((Ip),11)+GL((Ip),12)+GR((Ip),11)+GR((Ip),12))/4;
0329 GP = (GR((Ip),3)+GR((Ip),4)+GL((Ip),3)+GL((Ip),4))/4;
0330 % BWI1 = mean(GI1/max(GI1))*(max(f)-min(f));
0331 BWI1v2 = (max(f)-min(f))/length(f)*(sum(GI1))^2/(sum(GI1.^2));
0332 % BWI2 = mean(GI2/max(GI2))*(max(f)-min(f));
0333 BWI2v2 = (max(f)-min(f))/length(f)*(sum(GI2))^2/(sum(GI2.^2));
0334 % BWP = mean(GP/max(GP))*(max(f)-min(f));
0335 BWPv2 = (max(f)-min(f))/length(f)*(sum(GP))^2/(sum(GP.^2));
0336 fcI1 = sum(f.*GI1)/(sum(GI1));
0337 fcI2 = sum(f.*GI2)/(sum(GI2));
0338 fcP = sum(f.*GP)/(sum(GP));
0339 
0340 
0341 % fprintf('I1: f_centre = %4.3f, BW = %4.3f\n',fcI1,BWI1)
0342 % fprintf('I2: f_centre = %4.3f, BW = %4.3f\n',fcI2,BWI2)
0343 % fprintf('P: f_centre = %4.3f, BW = %4.3f\n',fcP,BWP)
0344 fprintf('I1: f_centre = %4.3f, BW = %4.3f\n',fcI1,BWI1v2)
0345 fprintf('I2: f_centre = %4.3f, BW = %4.3f\n',fcI2,BWI2v2)
0346 fprintf('P: f_centre = %4.3f, BW = %4.3f\n',fcP,BWPv2)
0347 
0348 % Now fit to aR1 and aL2.
0349 disp('analyzePassbandMeasurement_Pol:: Fitting to alpha spectrum')
0350 % Try a number of initial values for the lag error:
0351 Loff = 30*(1:10);
0352 Lstart = [-fliplr(Loff) 0 Loff]/1000;
0353 
0354 % Fit for phase too?
0355 % 0 = fit lag and amplitude ratios separately
0356 % 1 = fit lag, phase, and amplitude ratios separately
0357 % 2 = fit lag, amplitude ratios, and diode 3 simultaneously
0358 fitFlag = 0;
0359 
0360 cR = {fcable(Isel)*1e9 aR1(Isel)};
0361 cL = {fcable(Isel)*1e9 aL2(Isel)};
0362 cS = {fcable(Isel)*1e9 aR1(Isel) aL2(Isel) D3(Isel)};
0363 xR = [];
0364 xL = [];
0365 xS = [];
0366 fLbest = 1e10;
0367 fRbest = 1e10;
0368 fSbest = 1e10;
0369 for k=1:length(Lstart)
0370     switch fitFlag
0371         case 1
0372             [xkR,fvalR] = fminsearch(@(x) ssqerrP(x,cR),[1;Lstart(k);0]);
0373             [xkL,fvalL] = fminsearch(@(x) ssqerrP(x,cL),[1;Lstart(k);0]);
0374         case 0
0375             [xkR,fvalR] = fminsearch(@(x) ssqerr(x,cR),[1;Lstart(k)]);
0376             [xkL,fvalL] = fminsearch(@(x) ssqerr(x,cL),[1;Lstart(k)]);
0377         case 2
0378             % x(1) = xR
0379             % x(2) = LR
0380             % x(3) = xL
0381             % x(4) = LL
0382             % x(5) = xP
0383             [xk,fval] = fminsearch(@(x) ssqerrSim(x,cS),[1;Lstart(k);1;Lstart(k);1]);
0384     end
0385     if fitFlag == 2
0386         if fval < fSbest
0387             xS = xk;
0388             fSbest = fval;
0389         end
0390     else
0391         if fvalR < fRbest
0392             xR = xkR;
0393             fRbest = fvalR;
0394         end
0395         if fvalL < fLbest
0396             xL = xkL;
0397             fLbest = fvalL;
0398         end
0399     end
0400 end
0401 
0402 % Now assign the derived best fits to xR and xL if we are in the
0403 % simultaneous case:
0404 switch fitFlag
0405     case 2
0406         fitPhase = 0;
0407         xR = xS(1:2);
0408         xL = xS(3:4);
0409     case 1
0410         fitPhase = 1;
0411     case 0
0412         fitPhase = 0;
0413 end
0414 
0415 if fitPhase
0416     fprintf('I1: Path error = %3.1f mm, Gain ratio = %3.2f, Phase error = %3.2f deg\n',xR(2)*1000,xR(1),xR(3)*180/pi)
0417     fprintf('I2: Path error = %3.1f mm, Gain ratio = %3.2f, Phase error = %3.2f deg\n',xL(2)*1000,xL(1),xL(3)*180/pi)
0418     aRfit = alpha(xR(1),xR(2),fcable*1e9,xR(3));
0419     aLfit = alpha(xL(1),xL(2),fcable*1e9,xL(3));
0420 else
0421     fprintf('I1: Path error = %3.1f mm, Gain ratio = %3.2f\n',xR(2)*1000,xR(1))
0422     fprintf('I2: Path error = %3.1f mm, Gain ratio = %3.2f\n',xL(2)*1000,xL(1))
0423     aRfit = alpha(xR(1),xR(2),fcable*1e9,0);
0424     aLfit = alpha(xL(1),xL(2),fcable*1e9,0);
0425 end
0426 
0427 
0428 figure
0429 plot(fcable(Isel),aR1(Isel),'k.',fcable(Isel),ones(size(fcable(Isel)))*aRm,'k--',...
0430     fcable(Isel),aL2(Isel),'r.',fcable(Isel),ones(size(fcable(Isel)))*aLm,'r--',...
0431     fcable,aRfit,'k-',fcable,aLfit,'r-',...
0432     fcable,aR1,'k-.',fcable,aL2,'r-.')
0433 xlabel('Frequency [GHz]')
0434 ylabel('\alpha [goal = 1]')
0435 legend('\alpha_{I1}','mean','\alpha_{I2}','mean','Location','East')
0436 ylim([-1.05 1.05])
0437 if fitPhase
0438     sI1 = sprintf('I1: L=%3.1fmm, x=%3.2f, theta=%3.2fdeg',xR(2)*1000,xR(1),xR(3)*180/pi);
0439     sI2 = sprintf('I2: L=%3.1fmm, x=%3.2f, theta=%3.2fdeg',xL(2)*1000,xL(1),xL(3)*180/pi);
0440 else
0441     sI1 = sprintf('I1: L=%3.1fmm, x=%3.2f',xR(2)*1000,xR(1));
0442     sI2 = sprintf('I2: L=%3.1fmm, x=%3.2f',xL(2)*1000,xL(1));
0443 end
0444 
0445 title([sI1 '; ' sI2])
0446 
0447 saveas(gcf,[dirname 'alpha_Spectrum.pdf'])
0448 
0449 % and write the CW alpha values to disk:
0450 fn = [dirname 'alpha_cw.csv'];
0451 f1 = fopen(fn,'w');
0452 fprintf(f1,'# f [GHz], alpha_cw(I1), alpha_cw(I2), gain_I1 [dB], gain_I2 [dB]\n');
0453 fprintf(f1,'%4.3f, %4.3f, %4.3f, %4.3f, %4.3f\n',[fcable aR1 aL2 10*log10(gainR) 10*log10(gainL)]');
0454 fclose(f1);
0455 
0456 % Write the passbands to disk:
0457 fpoints = fcable(Ip);
0458 G1 = (GR(Ip,1)+GR(Ip,2)+GL(Ip,1)+GL(Ip,2))/4;
0459 G2 = (GR(Ip,3)+GR(Ip,4)+GL(Ip,3)+GL(Ip,4))/4;
0460 G3 = (GR(Ip,5)+GR(Ip,6)+GL(Ip,5)+GL(Ip,6))/4;
0461 G4 = (GR(Ip,7)+GR(Ip,8)+GL(Ip,7)+GL(Ip,8))/4;
0462 G5 = (GR(Ip,9)+GR(Ip,10)+GL(Ip,9)+GL(Ip,10))/4;
0463 G6 = (GL(Ip,11)+GL(Ip,12)+GR(Ip,11)+GR(Ip,12))/4;
0464 fn = [dirname 'Passband.csv'];
0465 f1 = fopen(fn,'w');
0466 fprintf(f1,'# f [GHz], Gain I1, Gain Q1, Gain U1, Gain Q2, Gain U2, Gain I2\n');
0467 fprintf(f1,'%4.3f, %4.3f, %4.3f, %4.3f, %4.3f, %4.3f, %4.3f\n',[fpoints G1 G2 G3 G4 G5 G6]');
0468 fclose(f1);
0469 
0470 % Plot the data!
0471 ylims = [20 90];
0472 
0473 figure
0474 plot(fcable(Ip),10*log10((GR(Ip,1)+GR(Ip,2)+GL(Ip,1)+GL(Ip,2))/4),'k-',...
0475     fcable(Ip),10*log10((GR(Ip,3)+GR(Ip,4)+GL(Ip,3)+GL(Ip,4))/4),'r-',...
0476     fcable(Ip),10*log10((GR(Ip,5)+GR(Ip,6)+GL(Ip,5)+GL(Ip,6))/4),'b-',...
0477     fcable(Ip),10*log10((GR(Ip,7)+GR(Ip,8)+GL(Ip,7)+GL(Ip,8))/4),'r--',...
0478     fcable(Ip),10*log10((GR(Ip,9)+GR(Ip,10)+GL(Ip,9)+GL(Ip,10))/4),'b--',...
0479     fcable(Ip),10*log10((GL(Ip,11)+GL(Ip,12)+GR(Ip,11)+GR(Ip,12))/4),'k--')
0480 ylim(ylims)
0481 xlabel('Frequency [GHz]')
0482 ylabel('Gain [dB]')
0483 legend('Channel 1','Channel 2','Channel 3','Channel 4','Channel 5',...
0484     'Channel 6','Location','SouthEast')
0485 title([sprintf('f_c [I1] = %4.3f GHz, BW = %4.3f GHz',fcI1,BWI1v2) ';' sprintf('f_c [I2] = %4.3f GHz, BW = %4.3f GHz',fcI2,BWI2v2) ';' sprintf('f_c [P] = %4.3f GHz, BW = %4.3f GHz',fcP,BWPv2)])
0486 grid on
0487 grid minor
0488 
0489 
0490 
0491 figure
0492 plot(fcable(Ip),10*log10((GPol(Ip,1)+GPol(Ip,2))/2),'k-',...
0493     fcable(Ip),10*log10((GPol(Ip,3)+GPol(Ip,4))/2),'r-',...
0494     fcable(Ip),10*log10((GPol(Ip,5)+GPol(Ip,6))/2),'b-',...
0495     fcable(Ip),10*log10((GPol(Ip,7)+GPol(Ip,8))/2),'r--',...
0496     fcable(Ip),10*log10((GPol(Ip,9)+GPol(Ip,10))/2),'b--',...
0497     fcable(Ip),10*log10((GPol(Ip,11)+GPol(Ip,12))/2),'k--')
0498 ylim(ylims)
0499 xlabel('Frequency [GHz]')
0500 ylabel('Gain [dB]')
0501 legend('Channel 1','Channel 2','Channel 3','Channel 4','Channel 5',...
0502     'Channel 6','Location','SouthEast')
0503 title('Pure U input')
0504 grid on
0505 grid minor
0506 
0507 saveas(gcf,[dirname 'C-BASS_Passband.pdf'])
0508 
0509 % Plot the response of all the diodes:
0510 figure
0511 subplot(3,4,1)
0512 plot(fcable(Ip),10*log10(GR(Ip,1)),'k-',fcable(Ip),10*log10(GL(Ip,1)),'r-')
0513 ylabel('Amplitude [dB]')
0514 ylim(ylims)
0515 title('Diode 1')
0516 subplot(3,4,2)
0517 plot(fcable(Ip),10*log10(GR(Ip,2)),'k-',fcable(Ip),10*log10(GL(Ip,2)),'r-')
0518 ylim(ylims)
0519 title('Diode 2')
0520 subplot(3,4,3)
0521 plot(fcable(Ip),10*log10(GR(Ip,11)),'k-',fcable(Ip),10*log10(GL(Ip,11)),'r-')
0522 ylim(ylims)
0523 title('Diode 11')
0524 subplot(3,4,4)
0525 plot(fcable(Ip),10*log10(GR(Ip,12)),'k-',fcable(Ip),10*log10(GL(Ip,12)),'r-')
0526 ylim(ylims)
0527 title('Diode 12')
0528 
0529 subplot(3,4,5)
0530 plot(fcable(Ip),10*log10(GR(Ip,3)),'k-',fcable(Ip),10*log10(GL(Ip,3)),'r-')
0531 ylim(ylims)
0532 ylabel('Amplitude [dB]')
0533 title('Diode 3')
0534 subplot(3,4,6)
0535 plot(fcable(Ip),10*log10(GR(Ip,4)),'k-',fcable(Ip),10*log10(GL(Ip,4)),'r-')
0536 ylim(ylims)
0537 title('Diode 4')
0538 subplot(3,4,7)
0539 plot(fcable(Ip),10*log10(GR(Ip,7)),'k-',fcable(Ip),10*log10(GL(Ip,7)),'r-')
0540 ylim(ylims)
0541 title('Diode 7')
0542 subplot(3,4,8)
0543 plot(fcable(Ip),10*log10(GR(Ip,8)),'k-',fcable(Ip),10*log10(GL(Ip,8)),'r-')
0544 ylim(ylims)
0545 title('Diode 8')
0546 
0547 subplot(3,4,9)
0548 plot(fcable(Ip),10*log10(GR(Ip,5)),'k-',fcable(Ip),10*log10(GL(Ip,5)),'r-')
0549 ylim(ylims)
0550 ylabel('Amplitude [dB]')
0551 xlabel('Frequency [GHz]')
0552 title('Diode 5')
0553 subplot(3,4,10)
0554 plot(fcable(Ip),10*log10(GR(Ip,6)),'k-',fcable(Ip),10*log10(GL(Ip,6)),'r-')
0555 ylim(ylims)
0556 xlabel('Frequency [GHz]')
0557 title('Diode 6')
0558 subplot(3,4,11)
0559 plot(fcable(Ip),10*log10(GR(Ip,9)),'k-',fcable(Ip),10*log10(GL(Ip,9)),'r-')
0560 ylim(ylims)
0561 xlabel('Frequency [GHz]')
0562 title('Diode 9')
0563 subplot(3,4,12)
0564 plot(fcable(Ip),10*log10(GR(Ip,10)),'k-',fcable(Ip),10*log10(GL(Ip,10)),'r-')
0565 ylim(ylims)
0566 xlabel('Frequency [GHz]')
0567 legend('RCP input','LCP input','Location','South')
0568 title('Diode 10')
0569 
0570 saveas(gcf,[dirname 'C-BASS_Passband_Diode_Responses.pdf'])
0571 
0572 end
0573 
0574 function err = ssqerr(x,c)
0575 % Function to fit the alpha value curve to:
0576 v = c{1}; % frequencies of the alpha points in Hz
0577 a = c{2}; % the alpha values
0578 
0579 % x(1) is ratio of voltage gains, x(2) is path error
0580 afit = alpha(x(1),x(2),v,0);
0581 err = sum((afit-a).^2);
0582 end
0583 
0584 function err = ssqerrP(x,c)
0585 % Function to fit the alpha value curve to:
0586 v = c{1}; % frequencies of the alpha points in Hz
0587 a = c{2}; % the alpha values
0588 
0589 % x(1) is ratio of voltage gains, x(2) is path error, x(3) is phase offset
0590 afit = alpha(x(1),x(2),v,x(3));
0591 err = sum((afit-a).^2);
0592 end
0593 
0594 function err = ssqerrSim(x,c)
0595 % Simultaneous fit to left, right, and pol channels
0596 v = c{1};
0597 aR = c{2};
0598 aL = c{3};
0599 D3 = c{4};
0600 
0601 % x(1) = xR
0602 % x(2) = LR
0603 % x(3) = xL
0604 % x(4) = LL
0605 % x(5) = xP
0606 
0607 aRfit = alpha(x(1),x(2),v,0);
0608 aLfit = alpha(x(3),x(4),v,0);
0609 D3fit = D3ratio(x(5),x(3),x(1),x(4),x(2),v);
0610 
0611 err = sum((aR-aRfit).^2+(aL-aLfit).^2+(D3-D3fit).^2);
0612 
0613 end
0614 
0615 function y = alpha(x,L,v,toff)
0616 % x = |g1|/|g2|
0617 % L = path length difference in m
0618 % v = frequency in Hz
0619 c = 3e8;
0620 Lc = L/c;
0621 
0622 y = 2*x*cos(2*pi*v*Lc+toff)/(1+x^2);
0623 end
0624 
0625 function y = D3ratio(xP,xL,xR,LL,LR,v)
0626 % the ratio of the diode 3 output powers
0627 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));
0628 
0629 end