automaticNorthColdLoad

% When do you want to read from, till and which schedules do you want to pull out
arcDir = '/data/arc'

0001 %% When do you want to read from, till and which schedules do you want to pull out
0002 %arcDir = '/data/arc'
0003 arcDir = '/elephant/CBASS_ARC/arc'
0004 date_start = '20-Jun-2013'%'01-Jan-2013';
0005 date_end = '07-Jan-2015'%'03-Jun-2014';
0006 path_to_redscript='reduc/redScript.red';
0007 sched_to_reduce = 'coldLoadStep';
0008 %sched_to_reduce = 'source_raster_';
0009 
0010 start_date_num = datenum(date_start);
0011 end_date_num = datenum(date_end);
0012 
0013 %% Select schedules which fall on the correct date
0014 [home,installeddir]=where_am_i();
0015 location = [home,'/',installeddir];
0016 %file = [location,'/','log/obs_log.html'];
0017 %file = '/home/cbassuser/cbass_analysis/log/obs_logSa.html'; %CBASS DELL
0018 file = '/elephant/cbass_analysis/log/obs_log.html'
0019 fid = fopen(file,'r');
0020 
0021 if fid < 0
0022     disp(['File failed to open: ',file]);
0023 end
0024 
0025 start_dates = {};
0026 end_dates = {};
0027 schedule = {};
0028 source_name = {};
0029 temp = 0;
0030 
0031 for j=0:(end_date_num - start_date_num)
0032     
0033     d_string = datestr(start_date_num + j);
0034     date = d_string;
0035     disp( date )
0036     fid = fopen(file,'r');
0037     
0038     i=0;
0039     
0040     while 1
0041         
0042         tline = fgetl(fid);
0043         if(~ischar(tline))
0044             %disp(['Closing file']);
0045             ST = fclose(fid);
0046             break;
0047         end;
0048         
0049         % Modified to ensure that the date appears in the START time.
0050         dloc = regexp(tline,date);
0051         if isempty(dloc)
0052             continue
0053         end
0054         
0055         % Well, is it?
0056         if (dloc(1) < 10)
0057             
0058             if ( length(regexp(tline,date)) ~=0)
0059                 i=i+1;
0060                 % disp(tline);
0061                 splitline = regexp(tline, '\s+', 'split');
0062                 % If the second element contains a / it is a pathname rather than a time
0063                 % i.e. no end time was written becuase the control system crashed.
0064                 if (~isempty(strfind(char(splitline(2)), '/')) )
0065                     start_dates(end+1) =  {'Error'};
0066                     end_dates(end+1) =  {'Error'};
0067                     schedule(end+1)  =  {'Error'};
0068                     source_name(end+1)={'Error'};
0069                 else
0070                     start_dates(end+1) =  splitline(1);
0071                     end_dates(end+1)  =  splitline(2);
0072                     schedule(end+1)  =  splitline(3);
0073                     % Grab source name if it is a raster or source scan
0074                     %                     if((strfind(char(schedule(end)),'source')))
0075                     %                         source_name(end+1) = splitline(6);
0076                     %                     else
0077                     %                         source_name(end+1) = schedule(end);% Just use sched name
0078                     %                     end
0079                     temp = i+1;
0080                 end
0081             end
0082             
0083             
0084         end
0085     end
0086 end
0087 %% Sort out empty cells (If there are any)
0088 % find empty cells
0089 emptyCells = cellfun(@isempty,schedule);
0090 % remove empty cells
0091 schedule(emptyCells) = [];
0092 start_dates(emptyCells) = [];
0093 end_dates(emptyCells) = [];
0094 
0095 %% Select which of these schedules are the correct ones
0096 %sched_to_reduce = 'coldLoadStep';
0097 
0098 start_string = {};
0099 end_string = {};
0100 sched_string = {};
0101 
0102 for i=1:length(start_dates)
0103     
0104     sched_bits = regexp(char(schedule(i)), '/', 'split');
0105     sched_name = char(sched_bits(length(sched_bits)));
0106     % strip out possible open bracket if there is one.
0107     sched_name = strrep(sched_name, '(','');
0108     
0109     % Check for error state.
0110     if strcmp(char(start_dates(i)),'Error')
0111         continue
0112     end
0113     
0114     if (exist('sched_to_reduce','var'))
0115         stripped_sched_to_reduce=sched_to_reduce;
0116         if (~isempty(regexp(sched_to_reduce,'^~')))
0117             stripped_sched_to_reduce = regexprep(sched_to_reduce,'^~','');
0118         end
0119         
0120         does_sched_match = ~isempty(regexp(sched_name,stripped_sched_to_reduce,'match'));
0121         if (does_sched_match) % This is a case sensitive regex style match
0122             start_string{end+1} =  strrep(char(start_dates(i)),'/','');
0123             end_string{end+1} =  strrep(char(end_dates(i)),'/','');
0124             sched_string{end+1} = char(schedule(i));
0125             continue;
0126         end;
0127     end;
0128     
0129 end
0130 
0131 start_string
0132 
0133 
0134 %% Begin finding the balance temp...
0135 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0136 %% Cold Load Step Balance
0137 % Get params we want
0138 imax = length( start_string );
0139 nchan = 2;
0140 Ttemp1 = zeros( imax, nchan );
0141 Ttemp2 = zeros( imax, nchan );
0142 for i=1:imax
0143     % Read in the data
0144     if ( tstr2mjd(char(end_string(i))) - tstr2mjd(char(start_string(i))) )<0.4
0145         a = char( [start_string(i),end_string(i)] )
0146         d = read_arc( a(1,:),a(2,:))%,'',arcDir )
0147         I1 = [ d.antenna0.receiver.switchData(:,1) fliplr(d.antenna0.receiver.switchData(:,4)) ];
0148         I2 = [ d.antenna0.receiver.switchData(:,21) fliplr(d.antenna0.receiver.switchData(:,24)) ];
0149         load1 = [ d.antenna0.receiver.switchData(:,2) fliplr(d.antenna0.receiver.switchData(:,3)) ];
0150         load2 = [ d.antenna0.receiver.switchData(:,22) fliplr(d.antenna0.receiver.switchData(:,23)) ];
0151         loadT = d.antenna0.thermal.ccTemperatureLoad;
0152         
0153         indMax = length( loadT );
0154         Ilen = length( I1 );
0155         interp_ratio = Ilen/indMax;
0156         x = [ 1:1:indMax ]';
0157         
0158         % I1 - Fit a straight line to the data
0159         for chan=1:nchan
0160             difference = I1(1:interp_ratio:Ilen,chan)-load1(1:interp_ratio:Ilen,chan);
0161             p = polyfit( x,difference,1 );
0162             ind = -p(2)/p(1);
0163             ind = round( ind );
0164             if ind<indMax & ind>0
0165                 temp1(chan) = loadT(ind);
0166             else
0167                 temp1(chan) = NaN;
0168             end
0169         end
0170         
0171         % I2 - Fit a straight line to the data
0172         for chan=1:nchan
0173             difference = I2(1:interp_ratio:Ilen,chan)-load2(1:interp_ratio:Ilen,chan);
0174             p = polyfit( x,difference,1 );
0175             ind = -p(2)/p(1);
0176             ind = round( ind );
0177             if ind<indMax & ind>0
0178                 temp2(chan) = loadT(ind);
0179             else
0180                 temp2(chan) = NaN;
0181             end
0182         end
0183         
0184         % Ttemp1/2 how they change with mjd
0185         Ttemp1(i,:) = temp1;
0186         Ttemp2(i,:) = temp2;
0187         mjd(i) = tstr2mjd( a(1,:) );
0188         
0189     else
0190         disp('Rejected schedule, too long')
0191         Ttemp1(i,:) = NaN;
0192         Ttemp2(i,:) = NaN;
0193         mjd(i) = tstr2mjd( a(1,:) );
0194     end
0195 end
0196 
0197 %% Prepare To Plot The Data
0198 Ttemp1(Ttemp1==0)=NaN;
0199 Ttemp2(Ttemp2==0)=NaN;
0200 [year1 month1 day1] = mjd2date(mjd);
0201 dateno = datenum( year1, month1, day1 );
0202 chpl = [1:1:nchan];
0203 
0204 X = repmat( chpl, 1, length(dateno) );
0205 tempY = repmat( dateno', 1, length(chpl) );
0206 Y = reshape( tempY',1, numel(tempY) );
0207 Z1 = reshape(Ttemp1',1,numel(Ttemp1) );
0208 Z2 = reshape(Ttemp2',1,numel(Ttemp2) );
0209 
0210 %% Save values as .txt file
0211 dlmwrite( 'northLoadBalTemp.txt', [year1',month1',day1',Ttemp1,Ttemp2], 'delimiter', '\t' )
0212 %% Newer Plotting
0213 figure()
0214 plot( mjd, Ttemp1(:,1),'.k' )
0215 hold()
0216 plot( mjd, Ttemp1(:,2),'.c' )
0217 plot( mjd, Ttemp2(:,1),'.r' )
0218 plot( mjd, Ttemp2(:,2),'.g' )
0219 %% New Plotting
0220 
0221 figure
0222 scatter3(Y, X, isnan(Z1'), 5, 'fill','sk')
0223 ylabel( 'Channel No' )
0224 xlabel( 'Date' )
0225 zlabel( 'Temp' )
0226 title( 'I1 and Load 1 Balance' )
0227 datetick('x',29,'keepticks','keeplimits')
0228 ylim( [0 nchan] )
0229 view(2)
0230 hold
0231 scatter3(Y, X, Z1', 20, Z1','fill','s')
0232 colorbar
0233 %caxis([8 30])
0234 xlim( [min(Y) max(Y)] )
0235 
0236 figure
0237 scatter3(Y, X, isnan(Z2'), 5, 'fill','sk')
0238 ylabel( 'Channel No' )
0239 xlabel( 'Date' )
0240 zlabel( 'Temp' )
0241 title( 'I2 and Load 2 Balance' )
0242 datetick('x',29,'keepticks','keeplimits')
0243 ylim( [0 nchan] )
0244 view(2)
0245 hold
0246 scatter3(Y, X, Z2', 20, Z2','fill','s')
0247 colorbar
0248 %caxis([8 30])
0249 xlim( [min(Y) max(Y)] )
0250 
0251 %% Saving
0252 save(['northLoadBalTemp_' date_start '_' date_end '.mat'],'X','Y','Z1','Z2')
0253 print 'SAVED'