Description of reduceData

0001 function d = reduceData(d, fn, interactive, field,sched_name)
0002 
0003 % function dcal = reduceData(d, fn, interactive, field,sched_name)
0004 %
0005 % VERSION 2.0
0006 %
0007 % Reduction Routine Using Matlab.
0008 %
0009 %  Inputs:
0010 %  d - raw/partially calibrated data structure
0011 %  output from read_arc.m
0012 %
0013 % Optional:
0014 %  fn - reduceData auto script, *.red
0015 %  interactive - 1/0 for user input during the pipeline
0016 %  field - fieldname you want the plots to be saved to.  if you leave it
0017 %          blank it will default to "other"
0018 %  sched_name - Name of the schedule if known (this will usually be passed in by
0019 %                           %               day_reduce)
0020 %
0021 % Stephen Muchovej
0022 %
0023 % Other contributors:
0024 %
0025 %
0026 %  HISTORY:
0027 %  version 0.0:  had many different sizes for the flags
0028 %  version 1.0:  all flagging will be done on basis of I1, all polarization,
0029 %            and I2    (implemented December 8, 2010)
0030 %  version 1.1:  changed step 'rfi1' to 'rfi'.  Removed step 'rfi2'.
0031 %                (MAS Jan 7, 2011).
0032 %  version 1.2:  routine list is now an input in .red file; alpha correction
0033 %  can be specified in that file, r-factor no longer combines the data. (sjcm, 2/24/2011)
0034 %  version 1.3:  gen, field now in script; removed 'intro' as an input
0035 %  (3/29/2010)
0036 %  version 1.4:  interactive now as input.  Major mods to subfunctions
0037 %  version 1.5:  plot in only 3 colors.  blah blah blah.  see sjcm''s email
0038 %  april 25,2011.
0039 %  version 1.6: sjcm:  changed rfactor to stokes.
0040 %  version 2.0: sjcm:  major modifications in how the noise diode
0041 %  temperature and opacities are calculated.  No longer backwards compatible
0042 %  with version 1.X
0043 % define version number
0044 REV='2.0';  
0045 % Add in a history field for history logging purposes.
0046 %d.history={'Began Data Reduction...'};
0047 
0048 % turn warnings off
0049 warning off
0050 
0051 % welcome!
0052 welcome(REV);
0053 
0054 olddir=pwd;
0055 
0056 is_reloaded=0; 
0057 
0058 % Check the inputs
0059 if (exist('fn'))
0060   if (~isempty(fn))
0061     disp('reduceData:: Reading script')
0062     pause(1)
0063     
0064     parm=readScript(fn);
0065     if (~matchRev(parm,REV))
0066       return
0067     end
0068     pause(1)
0069   else
0070     parm = [];
0071   end
0072 else
0073   parm = [];
0074 end
0075 
0076 if(length(d.array.frame.utc) < 90)
0077   display('reduceData:: Are you kidding me?');
0078   display('reduceData:: Are you seriously trying to reduce no data');
0079   display('reduceData:: GO AWAY!!! GO AWAY!! GO AWAY!! NOW!');
0080   return;
0081 end
0082 
0083 
0084 
0085 if (~exist('sched_name'))
0086  sched_name = 'not supplied';
0087 end;
0088 
0089 % check the inputs, and set defaults
0090 if(nargin==1)
0091   display('reduceData:: Interactive mode with only one input');
0092   display('reduceData:: Setting default values');
0093   parm = [];
0094   interactive = 1;
0095   field = 'other';
0096 elseif(nargin==2)
0097   % check if script was read properly
0098   
0099   if(isempty(parm))
0100     display('reduceData:: WARNING: script bad');
0101     display('reduceData:: Giving you interactive mode');
0102     interactive = 1;
0103     field = 'other';
0104   else
0105     if(isfield(parm, 'interactive'));
0106       interactive = parm.interactive.flag;
0107     else
0108       interactive = 1;
0109       parm.interactive.flag = 1;
0110     end
0111     if(isfield(parm, 'field'))
0112       field = parm.field.name;
0113     else
0114       field = 'other';
0115       parm.field.name = 'other';
0116     end
0117   end
0118 elseif(nargin==3)
0119   display('reduceData:: interactive variable in function call will supersede one in script');
0120   parm.interactive.flag = interactive;
0121   if(isfield(parm, 'field'))
0122     field = parm.field.name;
0123   else
0124     field = 'other';
0125     parm.field.name = field;
0126   end
0127 elseif(nargin==4)
0128   display('reduceData:: interactive  and field variable in function call will supersed one in script');
0129   parm.interactive.flag = interactive;
0130   if(isempty(field))
0131     field = 'other';
0132   end
0133   parm.field.name = field;
0134 end
0135 % done checking inputs
0136 
0137 % if it iss a calibrator schedule from before August 2012, we should check the
0138 % previous 10 minutes and following 10 minutes for some good skydips
0139 if(d.array.frame.utc(1) < date2mjd(2012, 08, 08))
0140   display('reduceData:: Reading in 10 minutes');
0141   display('reduceData:: So we could have a good skydip');
0142   timeRange = 10; % minutes
0143   timeRange = timeRange/60/24; % days
0144   % given that this is a scan cal observation, sometimes in the past our
0145   % skydips were not useful.  BUT -- we did have a skydip before and after the
0146   % schedule because survey scheds start and end with a skydip. so let us read
0147   % in some more data
0148   start1 = mjd2string(d.array.frame.utc(1) - timeRange);
0149   stop1  = mjd2string(d.array.frame.utc(1) - 0.5/60/60/24);
0150   dpre   = read_arc(start1, stop1);
0151   indpre = bitand(dpre.array.frame.features, 2^5)>0;
0152   if(any(indpre))
0153     display('reduceData:: Found a skydip before the track');
0154     [si ei] = findStartStop(indpre);
0155     start1  = mjd2string(dpre.array.frame.utc(last(si)));
0156     stop1   = mjd2string(dpre.array.frame.utc(last(ei)));
0157     dpre    = read_arc(start1, stop1);
0158     d       = catstruct(1, [dpre d]);
0159   end
0160 
0161   start2 = mjd2string(last(d.array.frame.utc)+0.5/60/60/24);
0162   stop2  = mjd2string(last(d.array.frame.utc) + timeRange);
0163   dpost  = read_arc(start2, stop2);
0164   indpost= bitand(dpost.array.frame.features, 2^5)>0;
0165   if(any(indpost))
0166     display('reduceData:: Found a skydip after the track');
0167     [si ei] = findStartStop(indpost);
0168     start2  = mjd2string(dpost.array.frame.utc(si(1)));
0169     stop2   = mjd2string(dpost.array.frame.utc(ei(1)));
0170     dpost   = read_arc(start2, stop2);
0171     d = catstruct(1, [d dpost]);
0172   end    
0173 end
0174 
0175 d.par = parm;
0176 
0177 % tag version to data
0178 d.rev=REV;
0179 
0180 % make sure the alpha database is uptodate
0181 isup2date = checkAlphaDatabase((d.array.frame.utc(1)));
0182 if(~isup2date)
0183   display('reduceData:: Your alpha database is not up to date');
0184   if(~isempty(parm))
0185     if(parm.interactive.flag)
0186       display('reduceData:: You will update the database');
0187       display('reduceData:: Don''t fudge it up');
0188       updateAlphaDatabase_DD(last(d.array.frame.utc));
0189     else
0190       display('reduceData:: Proceeding with last valid value');
0191     end
0192   end
0193 end
0194 
0195 
0196 % check run flags
0197 [plotparams, autorun]=setRunFlags(parm);  % from the file
0198 
0199 
0200 % check if we actually have any data.
0201 if(isempty(d.array.frame.utc) | all(d.array.frame.features==0))
0202   display('reduceData:: You have no data in this structure');
0203   display('reduceData:: Why would you try to calibrate it?');
0204   return;
0205 end
0206 
0207 % % Even if we're not saving plots, it's still necessary to run these two
0208 % % commands for the 'fits' step.  MAS 27-Feb-2012
0209   createDir(d,field, parm);
0210   maindir = getMainDir(d, field);
0211 
0212 % Grab a few things for later
0213 % thedate has format =  02-Nov-2013:14:35:52
0214 % thedate_day_portion is e.g. 02-Nov-2013
0215 % thedate_day_portion is passed further down to the individual wrappers
0216 % which is then used to forn the .txt filenames to avoid race conditions..
0217 
0218 [home,installeddir]=where_am_i();
0219 
0220 s2 = regexp(maindir, '/', 'split');
0221 thedate=char(s2(end));
0222 s3 = regexp(thedate, ':', 'split');
0223 thedate_day_portion=char(s3(1));
0224 
0225 
0226 
0227 % perform the first level of flagging
0228 if(plotparams.save)
0229   fn =  strcat(maindir, '/flags.txt');
0230   diary(fn);
0231   diary on;
0232   d = checkFlags(d, parm,thedate_day_portion);
0233   diary off;
0234 else
0235   d = checkFlags(d, parm,thedate_day_portion);
0236 end
0237 
0238 
0239 % if there's no noise source diode events, you can't calibrate.
0240 if(isempty(find(d.index.noise_event.fast)))
0241   display('reduceData:: There are no Noise Source diode events in your data');
0242   display('reduceData:: Can not calibrate this data set');
0243   return;
0244 end
0245 
0246 
0247 % check for the integration shortfall correction
0248 %[d,q] = logcal(d, 'shortfall');
0249 %if(~q)
0250 %  display('reduceData:: Correcting for the Integration Shortfall');
0251 %  d = shortfallCorr(d);
0252 %end
0253 
0254 
0255 
0256 % log start time
0257 time.startTime=d.array.frame.utc(1);
0258 time.stopTime=d.array.frame.utc(length(d.array.frame.utc));
0259 
0260 % plot the intro plots
0261 introPlots(d, plotparams, field);
0262 
0263 % set autoflag false
0264 autoflag=0;
0265 par=[];
0266 lastpar='auto';
0267 disp(' ')
0268 
0269 % Grab a few things for later
0270 % thedate has format =  02-Nov-2013:14:35:52
0271 % thedate_day_portion is e.g. 02-Nov-2013
0272 % thedate_day_portion is passed further down to the individual wrappers
0273 % which is then used to forn the .txt filenames to avoid race conditions..
0274  
0275 [home,installeddir]=where_am_i();
0276 
0277 s2 = regexp(maindir, '/', 'split');
0278 thedate=char(s2(end));
0279 s3 = regexp(thedate, ':', 'split');
0280 thedate_day_portion=char(s3(1));
0281 
0282 encoder_switch_filename = ['encoder_switch_',thedate_day_portion,'.txt'];
0283 
0284 % now the meat of the function
0285 while (1)
0286   if (autorun)
0287     [par,lastpar,autorun]=go(par,lastpar,autorun, parm);
0288   else
0289     [par,lastpar,autoflag]=prompt(par,lastpar,autoflag, parm);
0290   end
0291   
0292   switch par
0293     case 'help'
0294       menu;
0295       
0296     case 'intro'
0297       introPlots(d, plotparams, field);
0298       
0299     case 'pointing'
0300       display('reduceData:: Checking the pointing');
0301       d = pointingWrapper(d, plotparams, parm, field,thedate_day_portion);
0302 
0303       case 'deglitch'
0304       display('reduceData:: Deglitching data');
0305       d = deglitchWrapper(d, plotparams, parm, field);
0306       
0307     case 'rfi'
0308       display('reduceData:: RFI Flagging');
0309       d = rfiWrapper(d, plotparams, parm, field,thedate_day_portion); 
0310 
0311     case 'rfio'
0312       display('reduceData:: Old RFI flagging');
0313       d = rfioWrapper(d, plotparams, parm, field); 
0314     
0315       case 'mains'
0316       display('reduceData:: 60Hz Mains correction');
0317       d = mainsWrapper(d, plotparams, parm, field);  
0318 
0319     case 'alpha'
0320       display('reduceData:: Alpha Correction');
0321       d = alphaWrapper(d, plotparams, parm, field);        
0322       
0323     case 'stokes'
0324       display('reduceData:: Converting to Stokes parameters');
0325       d = stokesWrapper(d, plotparams, parm, field);      
0326       
0327     case 'load'
0328       display('reduceData:: Cold Load fluctuation Removal');
0329       d = loadWrapper(d, plotparams, parm, field);        
0330       
0331     case 'tau'
0332       display('Calculating the atmospheric opacity');
0333       d = tauCalWrapper(d, plotparams, parm, field,thedate_day_portion);   
0334       
0335     case 'save_d'
0336       display('reduceData:: About to save d structure...')  
0337       if(strcmp(parm.save_d.filename,'default') ~= 0)
0338         save_filename = [home,'/',installeddir,'/',thedate,'.mat']; 
0339       else
0340         save_filename = [home,'/',installeddir,'/',parm.save_d.filename]; 
0341       end
0342       disp(['reduceData:: Saving to ',save_filename]);
0343       save(save_filename,'d');
0344       disp('reduceData:: Save complete.');
0345 
0346      case 'load_d'
0347       display('reduceData:: About to load d structure...')  
0348       if(strcmp(parm.load_d.filename,'default') ~= 0)
0349         load_filename = [home,'/',installeddir,'/',thedate,'.mat']; 
0350       else
0351         load_filename = [home,'/',installeddir,'/',parm.load_d.filename]; 
0352       end
0353       disp(['reduceData:: Loading from ',load_filename]);
0354       load(load_filename,'d');
0355       disp('reduceData:: Load complete.');     
0356       is_reloaded=1; 
0357        
0358     case 'ground'
0359       display('reduceData:: Performing ground removal');     
0360       d = groundWrapper(d, plotparams, parm, field,thedate_day_portion);   
0361  
0362     case 'astro'
0363       display('reduceData:: Applying Astronomical Calibration');
0364       d = astroCalWrapper(d, plotparams, parm, field);     
0365       
0366     case 'power'
0367       display('reduceData:: Power Spectrum Plot');
0368       powerSpecWrapper(d, plotparams, parm, field);       
0369       
0370     case 'overf'
0371       display('reduceData:: 1/f plot');
0372       overfWrapper(d, plotparams, parm, field);  
0373       
0374     case 'rms'
0375       display('reduceData:: RMS calculation');
0376       rmsWrapper(d, plotparams, parm, field);                
0377       
0378     case 'tsys'
0379       display('Displaying system temperature of data');
0380       d = tsysWrapper(d, plotparams, parm, field,thedate_day_portion);
0381 
0382     case 'noise'
0383       display('reduceData:: Calculating the noise diode temperature');
0384       d = noiseWrapper(d, plotparams, parm, field,thedate_day_portion);                 
0385       
0386     case 'gain'
0387       display('reduceData:: Plotting quantities related to gain');
0388       d = gainWrapper(d, plotparams, parm, field);      
0389       
0390     case 'summary'
0391       display('reduceData:: Displaying the summary of all the flags');
0392       flagSummary(d, plotparams, parm, field,thedate_day_portion);
0393       
0394     case 'fits'
0395     
0396     if(is_reloaded == 1)
0397        reload_file_post_fix='_reload';
0398     else
0399        reload_file_post_fix='';
0400     end
0401 
0402       if(searchStruct(d, {'correction', 'tsys', 'allbad'}))
0403     if(d.correction.tsys.allbad == 1)
0404       display('reduceData:: Not writing out this crap data');
0405     else
0406       display('reduceData:: Writing out data to fits file');
0407       try
0408         writeOutFits(maindir, d, parm, reload_file_post_fix);
0409       catch exception
0410        disp('reduceData:: (1) Failed to write to fits data')
0411            rep = getReport(exception, 'basic');
0412            disp(rep);
0413            disp(['reduceData:: Returning to directory ',olddir]);
0414            cd(olddir);
0415         
0416       end
0417     end
0418       else
0419     try
0420       writeOutFits(maindir, d, parm, reload_file_post_fix);
0421     catch exception
0422       disp('reduceData:: (2) Failed to write to fits data')
0423           rep = getReport(exception, 'basic');
0424           disp(rep);         
0425        end
0426       end
0427 
0428       
0429 
0430       
0431     case 'auto'
0432       disp(' ')
0433       disp('reduceData:: ********Auto mode***********')
0434       disp('reduceData:: After each step any standard')
0435       disp('reduceData:: command will work')
0436       disp(' ')
0437       disp('reduceData:: To continue with calibration')
0438       disp('reduceData:: Hit RETURN without typing anything')
0439       disp(' ')
0440       autoflag=1;
0441       close all;
0442       
0443     case 'key'
0444       keyboard;
0445     
0446     case 'close'
0447       close all
0448       
0449     case 'plot off'
0450       plotflag=0;
0451       
0452     case 'plot on'
0453       plotflag=1;
0454 
0455     case 'save'
0456       save dataRed.mat
0457       
0458     case 'exit'
0459       disp(sprintf('reduceData:: Thank you for using Cbass Data Reduction. REDUCTION COMPLETE.\n'))
0460       cd(olddir);
0461       break;
0462 
0463     otherwise
0464       disp('reduceData:: Unknown command');
0465       disp(' ')
0466 
0467   end
0468 end
0469 
0470 % next, if we generated the plots, we should generate the pages.
0471 if(plotparams.web)
0472 
0473 
0474   display('reduceData:: Generating webpages of all previous plots');
0475   webpipe(d, field,0,sched_name,thedate_day_portion);
0476 
0477 end
0478 
0479 % make sure your plots appear from now on
0480 if(plotparams.plot==0)
0481   close all;
0482 end
0483 
0484   
0485 
0486 return;
0487   
0488   
0489     
0490 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0491 
0492 function welcome(rev)
0493 
0494 disp(' ')
0495 disp(' ')
0496 disp(' ')
0497 disp('***************************************************')
0498 disp('                 CBASS DATA REDUCTION ')
0499 disp(sprintf('                      rev %s',rev))
0500 disp('***************************************************')
0501 disp('Welcome to the Reduce Data Program');
0502 disp(' ')
0503 disp('Type ''help'' for menu')
0504 disp('  ')
0505 
0506 return;
0507 
0508 
0509 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0510 function q=matchRev(parm,rev)
0511 
0512 
0513 if (searchStruct(parm,{'rev','version'}))
0514   % compare only major number
0515   version=floor(parm.rev.version);
0516 else
0517   version=NaN;
0518 end
0519 
0520 if (isnan(version))
0521   disp(' ')
0522   disp('reduceData:: Script require unknown REV')
0523   disp('reduceData:: Continuing....')
0524   q=1;
0525   pause(1)
0526   return
0527 else
0528   q=version==floor(str2num(rev));
0529   if (~q)
0530     disp(sprintf('reduceData:: Script requires Redrum Major Version %s',...
0531         num2str(version)));
0532     pause(1)
0533   end
0534 end
0535 
0536 return;
0537 
0538 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0539 function d=checkFlags(d, parm,thedate_day_portion)
0540 
0541 % Identify the telecope we are using - used later in various sections
0542 antenna_no=identifyTelescope(d);
0543 
0544 parFlag = checkpar(parm, 'bitmask');
0545 
0546 if(parFlag)
0547   flagParams = parm.bitmask.bitVals;
0548 end
0549 
0550 encoder_switch_filename = ['encoder_switch_',thedate_day_portion,'.txt'];
0551 
0552 if(isfield(d.antenna0.servo, 'az_ccw_soft_limit'))
0553   perScan  = length(find(abs(deriv(d.antenna0.servo.slow_az_pos)) > ...
0554       3.5))./length(d.antenna0.servo.slow_az_pos);
0555   if(perScan > 0.5)
0556     % it is a survey track, we want to record it
0557     limVals = ...
0558     d.antenna0.servo.slow_az_pos(d.antenna0.servo.az_ccw_soft_limit>0);
0559     if(length(limVals)>0)
0560       % save [utcstart utcstop mean(az) rms(az) numSamps];
0561       thisData = [d.array.frame.utc(1) last(d.array.frame.utc) mean(limVals) ...
0562         sqrt(var(limVals)) length(limVals)];
0563       [home,installeddir]=where_am_i();
0564       %allEnc = load([home,'/',installeddir,'/constants/encoder_switch.txt']);
0565       %allEnc = [allEnc; thisData];
0566 
0567       save([home,'/',installeddir,'/constants/encoder_switch.txt'], 'thisData', '-ascii','-single','-append');
0568       save([home,'/',installeddir,'/constants/',encoder_switch_filename], 'thisData', '-ascii','-single','-append'); 
0569 
0570     end
0571     
0572     % we also want to correct it if there is a problem.
0573     if(length(limVals)>0)
0574 
0575       [switchmean switchstd] = check_encoder(d);
0576       meanVals = mean(limVals);
0577       offset   = meanVals - switchmean;
0578       stdVal   = switchstd;
0579       if(abs(offset)./stdVal > 3)
0580     % fix it.
0581     display('reduceData::checkFlags:: WARNING: ENCODER CALIBRATION INCORRECT');
0582     display('reduceData::checkFlags:: FIXING IT');
0583     d.antenna0.servo.slow_az_pos = d.antenna0.servo.slow_az_pos - offset;
0584     d.antenna0.servo.fast_az_pos = d.antenna0.servo.fast_az_pos - offset;
0585     d.antenna0.servo.az          = d.antenna0.servo.az - offset;
0586       end
0587     end
0588   end
0589 end
0590 
0591 
0592 
0593 if (~isfield(d, 'flags'))
0594   disp('reduceData::checkFlags::Flags not detected. Flagging data')
0595   pause(0.5)        
0596   disp(' ')    
0597   % Identify which antenna and pass to flagMask (needed for horizon
0598   % profile)
0599   
0600   d = flagMask(d,antenna_no);
0601   if(parFlag)
0602     d = flagData(d, flagParams);
0603   else
0604     d = flagData(d);
0605   end
0606 else
0607   if(~isfield(d.flags, 'mask'))
0608     disp('reduceData::checkFlags:: Flags not detected. Flagging data')
0609     pause(0.5)        
0610     disp(' ')    
0611     d = flagMask(d,antenna_no);
0612     if(parFlag)
0613       d = flagData(d, flagParams);
0614     else
0615       d = flagData(d);
0616     end    
0617   else
0618     disp('reduceData::checkFlags:: Flags found!')
0619     pause(0.5)
0620   end
0621 end
0622 
0623 if(~isfield(d, 'index'))
0624   display('reduceData::checkFlags:: Determining Indices');
0625   d = determineIndices(d);
0626 end
0627 
0628 if(~isfield(d.antenna0.servo, 'apparent'))
0629   % calculating apparent azimuth and elevation
0630   display('reduceData::checkFlags:: Calculating apparent Az/El');
0631   d = apparentAzEl(d, 0);
0632 end
0633 
0634 % apply nonliniarity correction
0635 %d = calibrate_linearity(d);
0636 
0637 return;
0638 
0639 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0640 function [plotparams, autorun]=setRunFlags(par)
0641 
0642 if (checkpar(par,'plotflag'))
0643   plotparams.plot = par.plotflag.flag;
0644 else
0645   plotparams.plot = 1;
0646 end
0647 
0648 if (checkpar(par,'autorun'))
0649   autorun=par.autorun.flag;
0650 else
0651   autorun=0;
0652 end
0653 
0654 if(checkpar(par, 'saveplot'))
0655   plotparams.save = par.saveplot.flag;
0656   plotparams.web = par.saveplot.web;
0657 else
0658   plotparams.save = 0;
0659   plotparams.web = 0;
0660 end
0661 
0662 if(checkpar(par, 'interactive'))
0663   plotparams.interactive = par.interactive.flag;
0664 else
0665   plotparams.interactive = 1;
0666 end
0667 
0668 if(plotparams.interactive == 1 & plotparams.plot == 0)
0669   display('reduceData::setRunFlags:: Warning: must be plot enabled in interactive mode');
0670   display('reduceData::setRunFlags:: Setting plotting to true');
0671   plotparams.plot = 1;
0672 end
0673 
0674 if(plotparams.interactive==1 & plotparams.save ==1)
0675   display('reduceData::setRunFlags:: Warning:  can not save plots in interactive mode');
0676   plotparams.save = 0;
0677 end
0678 
0679 
0680 return;
0681 
0682 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0683 function [par,lastpar,autorun]=go(par,lastpar,autorun, parm);
0684 
0685 if (autorun)
0686   par=proceed(lastpar, parm);
0687   lastpar=par;
0688   if(strcmp(par, 'exit'))
0689     autorun=0;
0690   end
0691 else
0692   par=input('CBASS> ', 's');
0693 end
0694 
0695 return;
0696 
0697 
0698 
0699 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0700 function par=proceed(par, parm)
0701 
0702 if(~isfield(parm, 'routines'))
0703   calstep={'deglitch', 'rfi', 'mains', 'alpha', 'load', 'pointing', 'tsys', ...
0704       'noise', 'gain', 'tau', 'stokes', ...
0705        'load_d','save_d','ground', 'astro', 'power', 'overf', 'rms', 'summary', 'fits'};
0706 else
0707   calstep = parm.routines.order;
0708 end
0709 
0710 
0711 
0712 if (strcmp(par,'auto'))
0713   close all
0714   par=calstep{1};
0715   return
0716 end
0717 
0718 f=find(strcmp(calstep,par));
0719 if(f==length(calstep))
0720   par = 'exit';
0721 else
0722   par=calstep{f(1)+1};
0723 end
0724 
0725 return;
0726 
0727 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0728 function [par,lastpar,autoflag]=prompt(par,lastpar,autoflag, parm);
0729 
0730 if (autoflag)
0731   autopar=input('Auto> ', 's');
0732   if isempty(autopar)
0733     par=proceed(lastpar, parm);
0734     lastpar=par;
0735     if (strcmp(par,'fits'))
0736       autoflag=0;
0737     end
0738   else
0739     par=autopar;
0740   end
0741 else
0742   par=input('CBASS> ', 's');
0743 end
0744 
0745 return;
0746 
0747 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0748 function menu()
0749 
0750 % function menu()
0751 %
0752 disp(' ');
0753 disp('  intro    - intro plots');
0754 disp('  pointing - checks any pointing crosses in the data');
0755 disp('  deglitch - runs the deglitching routine');
0756 disp('  rfi      - rfi flags');
0757 disp('  mains    - removes the effect of the mains');
0758 disp('  alpha    - splits load/sky, removes LHC/RHC phase offset, noise');
0759 disp('  load     - calibrates off the cold load variation');
0760 disp('  tau      - calculates atmospheric opacity');
0761 disp('  tsys     - calculates system temperatures');
0762 disp('  stokes   - converts to stokes (performs the r-factor on CLASSIC data');
0763 disp('  noise    - calculates the temperature of the noise source');
0764 disp('  gain     - calculates the time-varying gains');
0765 disp('  ground   - does the ground removal');
0766 disp('  astro    - applies opacity and gain fluctuation corrections');
0767 disp('  power    - displays power spectrum of data');
0768 disp('  overf    - calculates/displays 1/f plots');
0769 disp('  rms      - calculates noise in data');
0770 disp('  summary  - plots summary of why all the data was flagged');
0771 disp('  fits     - writes out fits file of on-source data');
0772 disp('  auto     - autocycle through calibration routine')
0773 disp('  key      - keyboard debugging mode')
0774 disp('  close    - close all figures')
0775 disp('  exit     - exits function')
0776 disp('  save     - saves the current sesssion to dataRed.mat');
0777 disp(' ');
0778 
0779 return;
0780 
0781 
0782 function flag=checkpar(par,type)
0783 
0784 % function flag=checkpar(par,type)
0785 %
0786 % check parameter structure for calibration "type"
0787 
0788 flag=0;
0789 
0790 if (~isempty(par))
0791   if (isfield(par,type))
0792     flag=1;
0793   end
0794 end
0795 
0796 return;
0797 
0798 
0799 function isUp2date = checkAlphaDatabase(time)
0800 
0801 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0802 %
0803 %  function isUp2date = checkAlphaDataBase(time)
0804 %
0805 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0806 
0807 t = readAlphaDatabase;
0808 
0809 
0810 if(last(t{1})>time)
0811   isUp2date = true;
0812 else
0813   isUp2date = false;
0814 end
0815 
0816 return;
0817 
0818 
0819 function olddir = writeOutFits(maindir, d, parm,  reload_file_post_fix)
0820 
0821 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0822 %
0823 %  function olddir = writeOutFits(maindir, d, parm)
0824 %
0825 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0826 %%%%%%%%%%%%%
0827 
0828 disp(['reduceData::writeOutFits:: reload postfix ',reload_file_post_fix])
0829 
0830 % Select the how you want to save data fits/hdf
0831 olddir =pwd;
0832 if(length(parm.fits.hdf) == 1)
0833   hdf=parm.fits.hdf; 
0834 else
0835   hdf = 0; % Defaults to .fits
0836 end
0837 display(['reduceData::writeOutFits:: hdf is ' , hdf])
0838 % Select the data you want to save as fits
0839 if(length(parm.fits.source) ~= 0)
0840  sourcename = parm.fits.source; 
0841 else
0842  sourcename ='source' % Defaults to just on-source data
0843 end
0844 
0845 display(['reduceData::writeOutFits:: Will cut data on ',sourcename])
0846 
0847 % Cut the data on required observation type, ready for saving to .fits
0848 % First check if the special cases of 'all' or 'allastro' have been
0849 % selected
0850 
0851 if (strcmpi(sourcename,'all'))
0852     df = d; % ie all the data will be writtn out as fits including skydips, noise events
0853 elseif (strcmpi(sourcename,'allastro'))
0854     df = cutObs(d, {'source', 'ncp','radio_point_scan','calibrator','beammap','abscal','radio_point_cross'},'or');
0855 else
0856     l = ['d.index.',char(sourcename),'.fast'];
0857     df = framecut(d,eval(l)); % Now only have the selected data in df
0858 end
0859 
0860 if (size(df.antenna0.receiver.utc,1)==0) 
0861     display('reduceData::writeOutFits:: Zero-sized data passed to writeOutFits.  Save will fail. Not going to make any maps');
0862   
0863 else
0864   
0865 
0866 root_name = strcat(maindir,'/fits/','fitstable');
0867 
0868 % JL TO DO check if file exists, if so delete.
0869 flags = df.flags.fast(:,1);
0870 startTimeMJD = df.antenna0.receiver.utc(1);
0871 endTimeMJD = df.antenna0.receiver.utc(length(df.antenna0.receiver.utc));
0872 
0873 
0874 [df q] = logcal(df, 'iv');
0875 if(~q)
0876   % JL 24/4/12 Everything shold be 8 column now irrespective of filter mode
0877   % parm.filtered.flag no longer exists.
0878   % df = getIV(df,parm.filtered.flag);
0879   disp('reduceData::writeOutFits:: Calling getIV');
0880   df = getIV(df,1); 
0881 end
0882 
0883 
0884 % need the equatorial coordinates
0885 
0886 % need the equatorial coordinates in J2000 co-ordinates
0887 % historically :  we do not need to calculate the ra/dec until we get ready to write things
0888 % out.
0889 
0890 disp('reduceData::writeOutFits:: Getting Equatorial co-ordinates');
0891 display('reduceData::writeOutFits:: Calculating RA/DEC'); 
0892 
0893 df = calcRADECJ2000(df);
0894 
0895  %
0896  %%%%%%%
0897  % uncomment this if in fact you want to keep RA and DEC in current epoch
0898  % form ( e.g. in case you want to compare the effect of precession)
0899  %
0900  % df = calcRADECcurrent(df);
0901  %
0902  %%%%%
0903  
0904 disp('Will now write out the .fits file');
0905 
0906 
0907 if (hdf==1)
0908     display('reduceData::writeOutFits:: Calling writefits map (hdf)');
0909     writeFitsMap([root_name,'.hdf'],df,startTimeMJD,endTimeMJD,1,1,hdf,flags); 
0910     clear df;
0911 else
0912     display('reduceData::writeOutFits:: Calling writefits map (fits)');
0913     writeFitsMap([root_name,'.fits'],df,startTimeMJD,endTimeMJD,1,1,hdf,flags);
0914     clear df;
0915 end
0916 
0917 disp('reduceData::writeOutFits:: Returned from writeFitsMap');
0918       
0919 % this -- exist('parm.fits.domapflag') returns zero for some reason
0920 % despite the fact that it exists!!!
0921 
0922 [home,installeddir]=where_am_i();
0923 location = [home,'/',installeddir];        
0924 %%%%%%%%%%%%%%%%%%%%%%%%%%
0925 % Get the date/time portion of the sirectory -will be useful for renaming
0926 % map and raw files later
0927 s2 = regexp(maindir, '/', 'split');
0928 thedate=char(s2(end));
0929 %%%%%%%%%%%%%%%%%%%%%%%%%%
0930 % Plot out the raw data for the fits file which has just been created...
0931 disp(['reduceData::writeOutFits:: Attempting to plot raw data to ',root_name,'_raw_plot.png...']); 
0932 olddir=pwd;
0933 disp(olddir);
0934 fitsdir=[maindir,'/fits/']; 
0935 disp(fitsdir);
0936 cd(fitsdir);
0937 %disp('reduceData::writeOutFits:: Changed dir');
0938 [aa, bb, cc, dd, hostnum] = whichHost();
0939 if(hostnum==11)
0940   raw_plot_command=['python2.7 ',location,'/webpage_logging/','plot_cbass_raw.py ',root_name,'.fits ',thedate,reload_file_post_fix,'_raw_plot.png'];
0941 else
0942   raw_plot_command=['python ',location,'/webpage_logging/','plot_cbass_raw.py ',root_name,'.fits ',thedate,reload_file_post_fix,'_raw_plot.png'];
0943 end
0944 unix(raw_plot_command);
0945 %disp('...done.');
0946 cd(olddir);
0947 %disp('Changed dir back');
0948 %%%%%%%%%%%%%%%%%%%%%%%%%%
0949 %Check if we want to make automated maps from the data
0950 
0951 if (strcmpi(parm.fits.domapflag,'yes')) % strcmpi() instead of '==' (MAS).
0952   
0953   disp('reduceData::writeOutFits:: Attempting to make a map from the fits file...');
0954 
0955   if(length(parm.fits.frame) ~= 0)
0956     frame = parm.fits.frame;
0957   else
0958     frame = 'heal';
0959   end
0960   disp('reduceData::writeOutFits:: Got frame');
0961 
0962   if(length(parm.fits.projection) ~= 0)
0963     projection = parm.fits.projection;
0964   else
0965     projection = 'moll';
0966   end
0967   disp('reduceData::writeOutFits:: Got projection');
0968 
0969 
0970 % Check if there is a specified .ini file for mapping
0971 
0972   if(length(parm.fits.ini) ~=0)
0973     descart_ini_file = parm.fits.ini;
0974   else
0975     if (strcmp(frame,'heal') & strcmp(projection,'moll'))
0976       descart_ini_file = 'mapping_params_heal_moll.ini';  
0977     elseif (strcmp(frame,'radec') & strcmp(projection,'tan'))
0978      name = df.antenna0.tracker.source(10)
0979      descart_ini_file = strcat(['mapping_params_radec_tan_',char(name),'.ini']); % AUTOMATICALLY PICKS UP SOURCE NAME
0980     elseif (strcmp(frame,'azel') & strcmp(projection,'plate'))
0981      descart_ini_file = 'mapping_params_azel_plate.ini';   
0982     else
0983      disp('reduceData::writeOutFits::  Could not parse frame and projection, defaulting to heal / moll');
0984      descart_ini_file = 'mapping_params_heal_moll.ini';
0985     end
0986   end
0987 
0988 % Check if there a designated python plotting file
0989     if(length(parm.fits.pyplot)~=0)
0990       png_plot_script = parm.fits.pyplot;
0991     else
0992       if (strcmp(frame,'heal') & strcmp(projection,'moll'))
0993     png_plot_script = 'plot_healpix_map.py';
0994       elseif (strcmp(frame,'radec') & strcmp(projection,'tan'))
0995     png_plot_script = 'plot_radec_map.py';
0996       elseif (strcmp(frame,'azel') & strcmp(projection,'plate'))
0997     png_plot_script = 'plot_cart_map.py';
0998       else
0999     disp('reduceData::writeOutFits::  Could not parse frame and projection, defaulting to heal / moll');
1000     png_plot_script = 'plot_healpix_map.py';
1001       end
1002     end
1003 
1004   disp(['reduceData::writeOutFits:: Will use ini file ',descart_ini_file,'  -- plot script ',png_plot_script]);
1005 
1006   %disp('reduceData::writeOutFits:: Setting Paths and Pythonpath...');
1007   % The following assumes that you have descart_cbass
1008   % and   plot_healpix_map.py installed in the appropriate locations.
1009   olddir=pwd;
1010   fitsdir=[maindir,'/fits/']; 
1011   cd(fitsdir);   
1012   current_pythonpath = getenv('PYTHONPATH');
1013   if strfind(current_pythonpath,'descart/data_selection')
1014      disp('reduceData::writeOutFits:: Pythonpath already set.');
1015   else
1016       if (isempty(current_pythonpath))
1017         new_pythonpath =[location,'/descart/data_selection/'];
1018       else  
1019         new_pythonpath =[location,'/descart/data_selection/:',current_pythonpath];
1020       end
1021       setenv('PYTHONPATH',new_pythonpath);
1022   end
1023 
1024   current_path = getenv('PATH'); 
1025   if isempty(strfind(current_path,'webpage_logging')) || ...
1026           isempty(strfind(current_path,'descart/map_making'))
1027       current_descart = getenv('DESCART');
1028       if (isempty(current_descart))
1029         loc2 = location;
1030       else
1031          loc2 = current_descart;
1032       end
1033       new_path=[loc2,'/webpage_logging/:',loc2,'/descart/map_making/:',current_path];
1034       setenv('PATH',new_path);
1035   else
1036       disp('reduceData::writeOutFits:: Path already set.');
1037   end
1038   disp('reduceData::writeOutFits:: ...done');
1039   %disp(['reduceData::writeOutFits:: New pythonpath',new_pythonpath])
1040 
1041   
1042   % delete the chunks
1043   disp('reduceData::writeOutFits:: Deleting chunks ...');
1044   if(hostnum~=11)
1045     delete_chunks_command=['python -m cbass.delete_chunks ',root_name,'.fits'];
1046   else
1047     delete_chunks_command=['python2.7 -m cbass.delete_chunks ',root_name,'.fits'];
1048   end
1049   unix(delete_chunks_command);
1050   disp('reduceData::writeOutFits:: ...done');
1051 
1052   % data selection stage
1053   do_select=1; 
1054   
1055   if(do_select==1)
1056     if(length(parm.dataselect)~=0)  
1057       split_flags = parm.dataselect.dflags;
1058       split_name  = parm.dataselect_name.fname;
1059       disp(['python -m cbass.split_data ',split_flags,' -N ', split_name])
1060       if(hostnum~=11)
1061     select_command=['python -m cbass.split_data ',split_flags,' -N ', split_name,' ',root_name,'.fits'];
1062       else
1063     select_command=['python2.7 -m cbass.split_data ',split_flags,' -N ', split_name,' ',root_name,'.fits'];
1064       end
1065     else
1066       if(hostnum~=11)
1067     select_command=['python -m cbass.split_data -zkg -N BASIC  ',root_name,'.fits'];
1068       else
1069     select_command=['python2.7 -m cbass.split_data -zkg -N BASIC  ',root_name,'.fits'];
1070       end      
1071     end
1072     % Override for ground subtraction.
1073     
1074     if(strcmp(split_name,'NoGalNo42'))
1075       
1076       disp(['reduceData::writeOutFits:: The split name is a ground subtraction one, splitting with split name',split_name]); 
1077       select_command=['python -m cbass.split_data ',split_flags,' -N ', split_name,' ',root_name,'.fits'];
1078       %select_command=['python -m cbass.split_data -N',split_name,' -zkg -G 10.0 -c 78@88:-10@0 ',root_name,'.fits'];
1079    
1080     end
1081  
1082     disp(['reduceData::writeOutFits:: Selecting the data with: ',select_command]);
1083     
1084     unix(select_command);
1085     disp('reduceData::writeOutFits:: ...done');
1086   end
1087 
1088   % prepare file list then called descart mapper.
1089   if (hdf==1)
1090     unix('echo fitstable.hdf > files.txt');
1091   else
1092     unix('echo fitstable.fits > files.txt');
1093   end
1094 
1095   % JL HACK - TO PREVENT A DESCART link error on mammoth / elephant
1096   % 17 = mammoth 8 = elephant
1097   current_ld_lib_path = getenv('LD_LIBRARY_PATH');
1098 
1099   if (hostnum == 17 | hostnum == 8)
1100       % elephant = '/usr/local/shared/cfitsio/lib/:'
1101       % mammoth = /usr/local/shared/cfitsio-3340-intel/lib/:
1102       % unix('echo $LD_LIBRARY_PATH');
1103 
1104       if (hostnum == 17)
1105         % mammoth
1106         setenv('LD_LIBRARY_PATH',['/usr/local/shared/cfitsio-3340-intel/lib/:',current_ld_lib_path])
1107       end
1108 
1109       if (hostnum == 8)
1110         % elephant
1111         setenv('LD_LIBRARY_PATH',['/usr/local/shared/cfitsio/lib/:',current_ld_lib_path]);
1112       end
1113 
1114   end
1115 
1116   disp(['reduceData::writeOutFits:: Calling descart_cbass to make map using ',descart_ini_file]);
1117   unix(['descart_cbass ',location,'/webpage_logging/',descart_ini_file,' cut_extension_name=',split_name]);
1118   
1119   % JL HACK - Now reset the LD_LIB_PATH to what it was before.
1120   setenv('LD_LIBRARY_PATH',current_ld_lib_path);
1121   %unix('echo $LD_LIBRARY_PATH');
1122 
1123   % Make the .png previews from the map fits file
1124   if (exist('map.fits')==2)
1125     disp(['reduceData::writeOutFits:: Making png previews using the dates ',png_plot_script]);
1126     if(hostnum==11)
1127       unix(['/home/sjcm/jaz/Python-2.6.8/python ',location,'/webpage_logging/',png_plot_script,' map.fits']);  
1128     else
1129       unix(['python ',location,'/webpage_logging/',png_plot_script,' map.fits']);        
1130     end
1131   
1132   else
1133     disp(['reduceData::writeOutFits:: map.fits does not exist will not make preview']);
1134   end;
1135 
1136   % Now rename relevant files
1137   disp('reduceData::writeOutFits:: Renaming files...');
1138   fitsfile =  strcat(thedate,reload_file_post_fix,'.fits');
1139   fitsmapfile= strcat(thedate,reload_file_post_fix,'_map.fits');
1140    hitsfile= strcat(thedate,reload_file_post_fix,'.hits');
1141   
1142   movefile('fitstable.fits',fitsfile);      
1143   movefile('map.fits',fitsmapfile);
1144   movefile('map.hits',hitsfile); 
1145 
1146   
1147 
1148 if (strcmp(frame,'heal') & strcmp(projection,'moll'))
1149 
1150   tmap= strcat(thedate,reload_file_post_fix,'_T.png');
1151   qmap= strcat(thedate,reload_file_post_fix,'_Q.png');
1152   umap= strcat(thedate,reload_file_post_fix,'_U.png');
1153   
1154   movefile('map_T.png',tmap);
1155   movefile('map_Q.png',qmap);
1156   movefile('map_U.png',umap);
1157   % Make Figs for the web page.
1158   copyfile(tmap,'fig1.png');
1159   copyfile(qmap,'fig2.png');
1160   copyfile(umap,'fig3.png');  
1161 
1162 elseif (strcmp(frame,'azel') & strcmp(projection,'plate'))
1163 
1164   tmap= strcat(thedate,reload_file_post_fix,'_azel_T.png');
1165   movefile('map_T.png',tmap);
1166   copyfile(tmap,'fig1.png');
1167 
1168   %%%
1169   % Angela added to make small raster maps from the radec,tan projections
1170   %%%
1171 elseif (strcmp(frame,'radec') & strcmp(projection,'tan'))
1172     disp('reduceData::writeOutFits:: Making radec tan .png maps')
1173     %date = thedate;
1174     name = df.antenna0.tracker.source(10);
1175     %d = script_pol_maps(d,name,date);
1176     if(hostnum==11)
1177       unix(['python2.7 ',location,'/webpage_logging/plot_radec_tan_map.py ', fitsmapfile]);
1178     else
1179       unix(['python ',location,'/webpage_logging/plot_radec_tan_map.py ', fitsmapfile]);
1180     end
1181     disp('reduceData::writeOutFits:: Copying fits file to named file');
1182 
1183     %disp(['cp ',fitsfile,' ',name,'_',fitsfile]); % For santiy in finding similar files later
1184     fitsfile
1185     movefile(fitsfile,[char(name),'_',fitsfile])
1186 end
1187 
1188   delete('files.txt');
1189   disp('reduceData::writeOutFits:: ...map making completed.');
1190 else
1191   disp('reduceData::writeOutFits:: Will not attempt to make a map...');
1192 end
1193   end
1194 disp(['reduceData::writeOutFits:: Returning to directory ',olddir]);
1195 cd(olddir);
1196 return;
1197
reduceData

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
