automaticNoiseStatsFromStares

% Empty arrays to fill with fit values

0001 %% Empty arrays to fill with fit values
0002 MJD = [];
0003 Ivals = [];
0004 Iupper = [];
0005 Ilower =[];
0006 Irms = [];
0007 
0008 Qvals = [];
0009 Qupper = [];
0010 Qlower =[];
0011 Qrms = [];
0012 
0013 Uvals = [];
0014 Uupper = [];
0015 Ulower =[];
0016 Urms = [];
0017 
0018 %% When do you want to read from, till and which schedules do you want to pull out
0019 %arcDir = '/data/arc'
0020 arcDir = '/elephant/CBASS_ARC/arc'
0021 date_start = '01-Oct-2012'%'01-Jan-2013';
0022 date_end = '19-Jan-2014'%'03-Jun-2014';
0023 path_to_redscript='reduc/redScript.red';
0024 sched_to_reduce = 'stare';
0025 
0026 start_date_num = datenum(date_start);
0027 end_date_num = datenum(date_end);
0028 %% Select schedules which fall on the correct date
0029 display( 'Selecting correct date range' )
0030 [home,installeddir]=where_am_i();
0031 location = [home,'/',installeddir];
0032 %file = [location,'/','log/obs_log.html'];
0033 %file = '/home/cbassuser/cbass_analysis/log/obs_logSa.html'; %CBASS DELL
0034 file = '/elephant/cbass_analysis/log/obs_log.html'
0035 fid = fopen(file,'r');
0036 
0037 if fid < 0
0038     disp(['File failed to open: ',file]);
0039 end
0040 
0041 start_dates = {};
0042 end_dates = {};
0043 schedule = {};
0044 source_name = {};
0045 temp = 0;
0046 
0047 for j=0:(end_date_num - start_date_num)
0048     d_string = datestr(start_date_num + j);
0049     date = d_string;
0050     disp( date )
0051     fid = fopen(file,'r');
0052     i=0;
0053     while 1
0054         tline = fgetl(fid);
0055         if(~ischar(tline))
0056             %disp(['Closing file']);
0057             ST = fclose(fid);
0058             break;
0059         end;
0060         % Modified to ensure that the date appears in the START time.
0061         dloc = regexp(tline,date);
0062         if isempty(dloc)
0063             continue
0064         end
0065         % Well, is it?
0066         if (dloc(1) < 10)
0067             if ( length(regexp(tline,date)) ~=0)
0068                 i=i+1;
0069                 % disp(tline);
0070                 splitline = regexp(tline, '\s+', 'split');
0071                 % If the second element contains a / it is a pathname rather than a time
0072                 % i.e. no end time was written becuase the control system crashed.
0073                 if (~isempty(strfind(char(splitline(2)), '/')) )
0074                     start_dates(end+1) =  {'Error'};
0075                     end_dates(end+1) =  {'Error'};
0076                     schedule(end+1)  =  {'Error'};
0077                     source_name(end+1)={'Error'};
0078                 else
0079                     start_dates(end+1) =  splitline(1);
0080                     end_dates(end+1)  =  splitline(2);
0081                     schedule(end+1)  =  splitline(3);
0082                     % Grab source name if it is a raster or source scan
0083                     %                     if((strfind(char(schedule(end)),'source')))
0084                     %                         source_name(end+1) = splitline(6);
0085                     %                     else
0086                     %                         source_name(end+1) = schedule(end);% Just use sched name
0087                     %                     end
0088                     temp = i+1;
0089                 end
0090             end
0091         end
0092     end
0093 end
0094 %% Sort out empty cells (If there are any)
0095 % find empty cells
0096 emptyCells = cellfun(@isempty,schedule);
0097 % remove empty cells
0098 schedule(emptyCells) = [];
0099 start_dates(emptyCells) = [];
0100 end_dates(emptyCells) = [];
0101 %% Select which of these schedules are the correct ones
0102 %sched_to_reduce = 'coldLoadStep';
0103 display( 'Selecting correct schedules' )
0104 
0105 start_string = {};
0106 end_string = {};
0107 sched_string = {};
0108 
0109 for i=1:length(start_dates)
0110     
0111     sched_bits = regexp(char(schedule(i)), '/', 'split');
0112     sched_name = char(sched_bits(length(sched_bits)));
0113     % strip out possible open bracket if there is one.
0114     sched_name = strrep(sched_name, '(','');
0115     
0116     % Check for error state.
0117     if strcmp(char(start_dates(i)),'Error')
0118         continue
0119     end
0120     
0121     if (exist('sched_to_reduce','var'))
0122         stripped_sched_to_reduce=sched_to_reduce;
0123         if (~isempty(regexp(sched_to_reduce,'^~')))
0124             stripped_sched_to_reduce = regexprep(sched_to_reduce,'^~','');
0125         end
0126         
0127         does_sched_match = ~isempty(regexp(sched_name,stripped_sched_to_reduce,'match'));
0128         if (does_sched_match) % This is a case sensitive regex style match
0129             start_string{end+1} =  strrep(char(start_dates(i)),'/','');
0130             end_string{end+1} =  strrep(char(end_dates(i)),'/','');
0131             sched_string{end+1} = char(schedule(i));
0132             continue;
0133         end;
0134     end;
0135     
0136 end
0137 
0138 start_string
0139 
0140 
0141 
0142 %%
0143 
0144 
0145 
0146 %% LOOP OVER SCHEDULES HERE
0147 jmax = length( start_string );
0148 for j=1:jmax
0149     % Read in the data
0150     if ( tstr2mjd(char(end_string(j))) - tstr2mjd(char(start_string(j))) )<0.4
0151         a = char( [start_string(j),end_string(j)] )
0152         d = read_arc( a(1,:),a(2,:));
0153         d = reduceData( d, './reduc/redScript_Luke.red' );
0154         d = determineIndices( d );
0155         % Find the start and stop indicies of each stare
0156         BLANK = d.index.blank.fast;
0157         len = length( BLANK );
0158         sIND = []; % start of a blank period
0159         eIND = []; % end of a blank period
0160         for i=1:len;
0161             if (i<len) & (BLANK(i)==0) & (BLANK(i+1)==1)
0162                 sIND(end+1)=i;
0163             end
0164             if (i>1) & (BLANK(i)==0) & (BLANK(i-1)==1)
0165                 eIND(end+1)=i;
0166             end
0167         end
0168         numOfStares = length(sIND)
0169         
0170         for i=1:numOfStares
0171             if length(eIND)~=length(sIND);
0172                 display( 'length(eIND) ~= length(sIND)' )
0173                 continue
0174             end
0175             idx = 1:(eIND(i)-sIND(i)+1);
0176             % Cut the relevant stare out
0177             I1 = d.antenna0.receiver.data(sIND(i):eIND(i),1);
0178             I1 = I1 - polyval(polyfit( idx', I1, 1 ),idx)';
0179             Q1 = d.antenna0.receiver.data(sIND(i):eIND(i),2);
0180             Q1 = Q1 - polyval(polyfit( idx', Q1, 1 ),idx)';
0181             U1 = d.antenna0.receiver.data(sIND(i):eIND(i),3);
0182             U1 = U1 - polyval(polyfit( idx', U1, 1 ),idx)';
0183             % Calculate the Power Spectra
0184             Fs = 100;
0185             L = eIND(i)-sIND(i)+1;
0186             NFFT = 2^nextpow2(L);
0187             x = Fs/2*linspace(0,1,NFFT/2+1); % Frequency
0188             % I
0189             y= fft(I1,NFFT);
0190             y = y.*conj(y)/(L*Fs);
0191             I1p = (y(1:NFFT/2+1)); % I1 Power
0192             % Q
0193             y= fft(Q1,NFFT);
0194             y = y.*conj(y)/(L*Fs);
0195             Q1p = (y(1:NFFT/2+1)); % Q1 Power
0196             % U
0197             y= fft(U1,NFFT);
0198             y = y.*conj(y)/(L*Fs);
0199             U1p = (y(1:NFFT/2+1)); % U1 Power
0200             % Fit model to the data
0201             fitLength = 4.;
0202             ft = fittype( 'noiseModel( x, sigma_w, f_knee, alpha )' );
0203             x_ = x(2:end/fitLength);
0204             I1p_ = I1p(2:end/fitLength);
0205             Q1p_ = Q1p(2:end/fitLength);
0206             U1p_ = U1p(2:end/fitLength);
0207             % I
0208             Iguess = [std(I1p(end/fitLength:end)),0.1,1];
0209             [Iparams, Igof, output] = fit( x_', real(log(I1p_)), ft, 'StartPoint', Iguess, 'Lower', [0.,0.,0.], 'Exclude', (x_<1.3)&(x_>1.1) ,'TolFun', 1e-20, 'TolX',1e-10 );
0210             Iconf = confint(Iparams,0.95);
0211             %Q
0212             Qguess = [std(Q1p(end/fitLength:end)),0.01,1];
0213             [Qparams, Qgof, output] = fit( x_', real(log(Q1p_)), ft, 'StartPoint', Qguess, 'Lower', [0.,0.,0.], 'Exclude', (x_<1.3)&(x_>1.1) ,'TolFun', 1e-20, 'TolX',1e-10 );
0214             Qconf = confint(Qparams,0.95);
0215             %U
0216             Uguess = [std(U1p(end/fitLength:end)),0.01,1];
0217             [Uparams, Ugof, output] = fit( x_', real(log(U1p_)), ft, 'StartPoint', Uguess, 'Lower', [0.,0.,0.], 'Exclude', (x_<1.3)&(x_>1.1) ,'TolFun', 1e-20, 'TolX',1e-10 );
0218             Uconf = confint(Uparams,0.95);
0219             % Now add these values to the arrays
0220             Ivals(end+1,1:3) = [Iparams.alpha, Iparams.f_knee, Iparams.sigma_w];
0221             Ilower(end+1,1:3) = [Iconf(1,1),Iconf(1,2),Iconf(1,3)];
0222             Iupper(end+1,1:3) = [Iconf(2,1),Iconf(2,2),Iconf(2,3)];
0223             Irms(end+1) = Igof.rmse;
0224             Qvals(end+1,1:3) = [Qparams.alpha, Qparams.f_knee, Qparams.sigma_w];
0225             Qlower(end+1,1:3) = [Qconf(1,1),Qconf(1,2),Qconf(1,3)];
0226             Qupper(end+1,1:3) = [Qconf(2,1),Qconf(2,2),Qconf(2,3)];
0227             Qrms(end+1) = Qgof.rmse;
0228             Uvals(end+1,1:3) = [Uparams.alpha, Uparams.f_knee, Uparams.sigma_w];
0229             Ulower(end+1,1:3) = [Uconf(1,1),Uconf(1,2),Uconf(1,3)];
0230             Uupper(end+1,1:3) = [Uconf(2,1),Uconf(2,2),Uconf(2,3)];
0231             Urms(end+1) = Ugof.rmse;
0232             MJD(end+1) = d.antenna0.receiver.utc(sIND(i));
0233         end
0234     end
0235 end
0236 
0237 %% Save Variables
0238 save( ['skyStareNoiseStats_' date_start '_' date_end '.mat'] , 'Ivals', 'Qvals', 'Uvals', 'Irms', 'Qrms', 'Urms', 'Ilower', 'Qlower', 'Ulower', 'Iupper', 'Qupper', 'Uupper', 'MJD' )
0239 display( 'Saved' )
0240 
0241 %% Convert Date
0242 [year month day] = mjd2date(MJD);
0243 dateno = datenum( year, month, day );
0244 
0245 % %% Side By Side Plotting
0246 % % KNEE FREQ
0247 % figure
0248 % subplot(131)
0249 % semilogy( dateno(Ivals(:,2)<10.), Ivals(Ivals(:,2)<10.,2), '.k', 'markerSize',10 )
0250 % datetick('x','mmm-YY','keepticks')%,'keeplimits')
0251 % xlabel( 'Date' )
0252 % ylabel( 'f_{knee}' )
0253 % title( 'I' )
0254 % grid()
0255 %
0256 % subplot(132)
0257 % semilogy( dateno(Qvals(:,2)<10.), Qvals(Qvals(:,2)<10.,2), '.k', 'markerSize',10 )
0258 % datetick('x','mmm-YY','keepticks')%,'keeplimits')
0259 % xlabel( 'Date' )
0260 % title( 'Q' )
0261 % grid()
0262 %
0263 % subplot(133)
0264 % semilogy( dateno(Uvals(:,2)<10.), Uvals(Uvals(:,2)<10.,2), '.k', 'markerSize',10 )
0265 % datetick('x','mmm-YY','keepticks')%,'keeplimits')
0266 % xlabel( 'Date' )
0267 % title( 'U' )
0268 % grid()
0269 %
0270 % % ALPHA
0271 % figure
0272 % subplot(131)
0273 % plot( dateno(Ivals(:,1)<10.), Ivals(Ivals(:,1)<10.,1), '.k', 'markerSize',10 )
0274 % datetick('x','mmm-YY','keepticks')%,'keeplimits')
0275 % xlabel( 'Date' )
0276 % ylabel( '\alpha' )
0277 % title( 'I' )
0278 % grid()
0279 %
0280 % subplot(132)
0281 % plot( dateno(Qvals(:,1)<10.), Qvals(Qvals(:,1)<10.,1), '.k', 'markerSize',10 )
0282 % datetick('x','mmm-YY','keepticks')%,'keeplimits')
0283 % xlabel( 'Date' )
0284 % title( 'Q' )
0285 % grid()
0286 %
0287 % subplot(133)
0288 % plot( dateno(Uvals(:,1)<10.), Uvals(Uvals(:,1)<10.,1), '.k', 'markerSize',10 )
0289 % datetick('x','mmm-YY','keepticks')%,'keeplimits')
0290 % xlabel( 'Date' )
0291 % title( 'U' )
0292 % grid()
0293 %
0294 % % SIGMA_W
0295 % figure
0296 % subplot(131)
0297 % semilogy( dateno(Ivals(:,3)<10.), Ivals(Ivals(:,3)<10.,3), '.k', 'markerSize',10 )
0298 % datetick('x','mmm-YY','keepticks')%,'keeplimits')
0299 % xlabel( 'Date' )
0300 % ylabel( '\sigma_w' )
0301 % title( 'I' )
0302 % grid()
0303 %
0304 % subplot(132)
0305 % semilogy( dateno(Qvals(:,3)<10.), Qvals(Qvals(:,3)<10.,3), '.k', 'markerSize',10 )
0306 % datetick('x','mmm-YY','keepticks')%,'keeplimits')
0307 % xlabel( 'Date' )
0308 % title( 'Q' )
0309 % grid()
0310 %
0311 % subplot(133)
0312 % semilogy( dateno(Uvals(:,3)<10.), Uvals(Uvals(:,3)<10.,3), '.k', 'markerSize',10 )
0313 % datetick('x','mmm-YY','keepticks')%,'keeplimits')
0314 % xlabel( 'Date' )
0315 % title( 'U' )
0316 % grid()
0317 %
0318 %
0319 % %% On Top Plotting
0320 % figure
0321 % semilogy( dateno(Ivals(:,2)<10.), Ivals(Ivals(:,2)<10.,2), '.r', 'markerSize',10 )
0322 % hold
0323 % semilogy( dateno(Qvals(:,2)<10.)+1, Qvals(Qvals(:,2)<10.,2), '.g', 'markerSize',10 )
0324 % semilogy( dateno(Uvals(:,2)<10.)+2, Uvals(Uvals(:,2)<10.,2), '.b', 'markerSize',10 )
0325 %
0326 % datetick('x','mmm-YY','keepticks')%,'keeplimits')
0327 % xlabel( 'Date' )
0328 % ylabel( 'f_{knee}' )
0329 % title( 'I-red, Q-green, U-blue' )
0330 % grid()
0331 %
0332 % % ALPHA
0333 % figure
0334 % plot( dateno(Ivals(:,1)<10.), Ivals(Ivals(:,1)<10.,1), '.r', 'markerSize',10 )
0335 % hold
0336 % plot( dateno(Qvals(:,1)<10.)+1, Qvals(Qvals(:,1)<10.,1), '.g', 'markerSize',10 )
0337 % plot( dateno(Uvals(:,1)<10.)+2, Uvals(Uvals(:,1)<10.,1), '.b', 'markerSize',10 )
0338 %
0339 % datetick('x','mmm-YY','keepticks')%,'keeplimits')
0340 % xlabel( 'Date' )
0341 % ylabel( '\alpha' )
0342 % title( 'I-red, Q-green, U-blue' )
0343 % grid()
0344 %
0345 % % SIGMA_W
0346 % figure
0347 % semilogy( dateno(Ivals(:,3)<10.), Ivals(Ivals(:,3)<10.,3), '.r', 'markerSize',10 )
0348 % hold
0349 % semilogy( dateno(Qvals(:,3)<10.)+1, Qvals(Qvals(:,3)<10.,3), '.g', 'markerSize',10 )
0350 % semilogy( dateno(Uvals(:,3)<10.)+2, Uvals(Uvals(:,3)<10.,3), '.b', 'markerSize',10 )
0351 %
0352 % datetick('x','mmm-YY','keepticks')%,'keeplimits')
0353 % xlabel( 'Date' )
0354 % ylabel( '\sigma_w' )
0355 % title( 'I-red, Q-green, U-blue' )
0356 % grid()
0357 %
0358 %
0359 % %% RMS Plotting
0360 % figure
0361 % plot( dateno, Irms,'.' )
0362 %
0363 % %% OldPlotting
0364 % hold
0365 % errorbar( dateno(Ivals(:,2)<10.), Ivals(Ivals(:,2)<10.,2), Ilower(Ivals(:,2)<10.,2), Iupper(Ivals(:,2)<10.,2),'ok' )
0366 %
0367 % %I
0368 % figure
0369 % errorbar( dateno, Ivals(:,2), Ilower(:,2), Iupper(:,2),'ok' )
0370 % ylabel( 'F_{knee} Hz' )
0371 % datetick('x','mmm-YY','keepticks','keeplimits')
0372 % xlabel( 'Date' )
0373 % title( 'I' )
0374 %
0375 % figure
0376 % errorbar( MJD, Ivals(:,1), Ilower(:,1), Iupper(:,1),'ok' )
0377 % ylabel( '\alpha' )
0378 % xlabel( 'MJD' )
0379 % title( 'I' )
0380 %
0381 % figure
0382 % errorbar( MJD, Ivals(:,3), Ilower(:,3), Iupper(:,3),'ok' )
0383 % ylabel( '\sigma_w K' )
0384 % xlabel( 'MJD' )
0385 % title( 'I' )
0386 %
0387 % % Q
0388 % figure
0389 % errorbar( MJD, Qvals(:,2), Qlower(:,2), Qupper(:,2),'ok' )
0390 % ylabel( 'F_{knee} Hz' )
0391 % xlabel( 'MJD' )
0392 % title( 'Q' )
0393 %
0394 % figure
0395 % errorbar( MJD, Qvals(:,1), Qlower(:,1), Qupper(:,1),'ok' )
0396 % ylabel( '\alpha' )
0397 % xlabel( 'MJD' )
0398 % title( 'Q' )
0399 %
0400 % figure
0401 % errorbar( MJD, Qvals(:,3), Qlower(:,3), Qupper(:,3),'ok' )
0402 % ylabel( '\sigma_w K' )
0403 % xlabel( 'MJD' )
0404 % title( 'Q' )
0405 %
0406 % % U
0407 % figure
0408 % errorbar( MJD, Uvals(:,2), Ulower(:,2), Uupper(:,2),'ok' )
0409 % ylabel( 'F_{knee} Hz' )
0410 % xlabel( 'MJD' )
0411 % title( 'U' )
0412 %
0413 % figure
0414 % errorbar( MJD, Uvals(:,1), Ulower(:,1), Uupper(:,1),'ok' )
0415 % ylabel( '\alpha' )
0416 % xlabel( 'MJD' )
0417 % title( 'U' )
0418 %
0419 % figure
0420 % errorbar( MJD, Uvals(:,3), Ulower(:,3), Uupper(:,3),'ok' )
0421 % ylabel( '\sigma_w K' )
0422 % xlabel( 'MJD' )
0423 % title( 'U' )
0424 %
0425 %
0426 % %figure;
0427 % %loglog( x(2:end), I1p(2:end) )
0428 % %hold()
0429 % %loglog( x_, exp(noiseModel(x_,Iparams.sigma_w,Iparams.f_knee,Iparams.alpha)),'r' )
0430 % %figure()
0431 % %errorbar( MJD, Ivals(:,2), Ilower(:,2), Iupper(:,2),'o' )
0432 %
0433 %
0434 %