writeFitsMap

function to write to fits file for the Descart Mapper
fitsfilenam - Name of Fits file to write data to
d -standard C-BASS data struct
batchStartMJD- MJD start date of the batch that will be processed (i.e
more than one fits file
batchEndMJD- as above
thisScanNumber - which scan is this of the total data set
totalScans- how many scans are there in the entire batch

0001 function scan = writeFitsMap(fitsFilename,d,batchStartMJD,batchEndMJD,thisScanNumber,totalScans,hdf_mode,varargin)
0002 
0003 %function to write to fits file for the Descart Mapper
0004 %fitsfilenam - Name of Fits file to write data to
0005 %d -standard C-BASS data struct
0006 %batchStartMJD- MJD start date of the batch that will be processed (i.e
0007 %more than one fits file
0008 %batchEndMJD- as above
0009 %thisScanNumber - which scan is this of the total data set
0010 %totalScans- how many scans are there in the entire batch
0011 
0012 % 4-Feb-2011 ACT added extra bit to allow both 6 column and 8 column calibrated.data to be
0013 % written out:
0014 %
0015 % 15-Apr-2014 Paddy: replaced call to sunUpDown (which defines sunrise as
0016 % solar elevation = 0 in a byzantine way) with d.flags.DayNight set in
0017 % flagMask using the measured observatory horizon. Value = 0 for night (as
0018 % before, 2 for day, 1 for solar centre within 1 degree of estimated
0019 % horizon (i.e. on horizon within systematic uncertainty).
0020 % Also added new output stream solarDist: distance in degrees from pointing
0021 % position to Sun.
0022 % Also changed sunBuffer to 0 mins, i.e. rely on the flag.
0023 %
0024 %
0025 
0026 
0027 
0028 
0029 
0030 optargin = size(varargin,2);
0031 if(optargin>0)
0032     flagsIn=varargin{1};
0033 else
0034     flagsIn=zeros(length(d.index.skydip.fast),1);
0035 end
0036 
0037 
0038 disp('writeFitsMap:: Starting writeFitsMap')
0039 sunBuffer=0;
0040 
0041 %lat=37.233; %ovro Latitude in degrees
0042 %lon= 118.282; %ovro longitude West positive in degrees
0043   lon=d.antenna0.tracker.siteActual(2,1);
0044   lat=d.antenna0.tracker.siteActual(2,2);
0045   
0046 
0047 startTimeMJD = d.antenna0.receiver.utc(1);
0048 endTimeMJD = d.antenna0.receiver.utc(length(d.antenna0.receiver.utc));
0049 
0050 disp('writeFitsMap:: EndTimeMJD writeFitsMap')
0051 %create the Fits Data bundle
0052 relativeTime = d.antenna0.receiver.utc-startTimeMJD;
0053 equatorialPosn = d.antenna0.servo.equa(:,1:2);
0054 horizontalPosn = d.antenna0.servo.apparent(:,1:2);
0055 correctedData = d.antenna0.receiver.data; %6 columns of data
0056 flags = logical(flagsIn)|logical(d.index.skydip.fast)|logical(d.index.noise.fast);
0057 sunPos = double(d.flags.dayNight);
0058 disp('writeFitsMap:: Flags writeFitsMap')
0059 try
0060 
0061 solarDist = d.antenna0.servo.solarDist;
0062 [galacticPosn,TotRot]=coco(equatorialPosn,'j2000.0','g','r','r');
0063 dataBundle = [relativeTime equatorialPosn horizontalPosn correctedData flags sunPos solarDist];
0064 %dataBundle = [relativeTime equatorialPosn horizontalPosn correctedData flags sunPos ];
0065 catch
0066 disp('writeFitsMap:: Failed galacticPosn,databundle in writeFitsMap')
0067 end
0068 
0069     % Determine number of columns 8 or 4 (CLASSIC or FILTERED)
0070     % and pass number down to cbass_write_data as the integer datatype
0071     % 8-columns gives [I Q1 U1 Q2 U2 Q3 U3 V]
0072     % 4-columns gives [I Q U V]
0073     datatype = size(d.antenna0.receiver.data,2);
0074     fprintf('datatype = %d\n',datatype)
0075 
0076     fits_header = prepareHeader(d);
0077     size_fits_header = size(fits_header);
0078     no_fits_keywords = size_fits_header(2);
0079     
0080   if (hdf_mode == 1 )
0081     disp('writeFitsMap:: Calling cbass_write_hdf');
0082     cbass_write_hdf(fitsFilename,dataBundle,0,0,startTimeMJD,endTimeMJD,thisScanNumber,totalScans,batchStartMJD,batchEndMJD,sunBuffer,datatype);
0083   else
0084     disp('writeFitsMap:: Calling cbass_write_data');  
0085     cbass_write_data(fitsFilename,dataBundle,0,0,startTimeMJD,endTimeMJD,thisScanNumber,totalScans,batchStartMJD,batchEndMJD,sunBuffer,datatype,fits_header,no_fits_keywords);
0086   end
0087 %catch
0088 %disp('writeFitsMap:: Failed to write the FITS file with cbass_write_data'
0089 %end
0090 scan = thisScanNumber;
0091 disp('writeFitsMap:: End writeFitsMap.')
0092 end
0093 
0094 
0095 function fits_header = prepareHeader(d)
0096 
0097 % Function to populate "dynamic" fits keywords from the d structure, and also
0098 % static hardwired headers.
0099 % output is a 3 x num_keywords array,
0100 % col1 = HEADER, col2 = VALUE col3 = COMMENT col4 = VALUE_TYPE (string double, etc.)
0101 [antenna_no antenna_name] = identifyTelescope(d);
0102 telescope_constants;  % load data about telescopes
0103 epoch_label= d.antenna0.servo.epoch
0104 
0105 num_keywords = 18;
0106 
0107 fits_header=cell(4,num_keywords);
0108 
0109 i=1;
0110 fits_header(:,i)={'TELESCOP'   antenna_name    'Name of Telescope'    'String'};
0111 i=i+1;
0112 
0113 fits_header(:,i)={'ORIGIN  '  char(antenna_observatory(antenna_no))  'Origin of file' 'String'};
0114 i=i+1;
0115 
0116 fits_header(:,i)={'---- x,y,z triplet for CBASS relative to centre of earth ----' 'blah' 'blah' 'Comment'};
0117 i=i+1;
0118 
0119 fits_header(:,i)={'ALT-OBS ' antenna_altitude(antenna_no) '[m] Height of observatory above sea level' 'Double'};
0120 i=i+1;
0121 
0122 
0123 fits_header(:,i)={'LAT-OBS ' antenna_latitude(antenna_no) '[deg] Latitude of observatory' 'Double'};
0124 i=i+1;
0125 
0126 fits_header(:,i)={'LONG-OBS' antenna_longitude(antenna_no) '[deg] Longitude of observatory' 'Double'};
0127 i=i+1;
0128 
0129 fits_header(:,i)={'RADECSYS' 'FK5     ' 'Reference frame for RA/DEC values' 'String'};
0130 i=i+1;
0131 
0132 fits_header(:,i)={'EQUINOX' epoch_label 'Epoch of reference equinox' 'String'};
0133 i=i+1;
0134 
0135 fits_header(:,i)={'VERSION ' 2 'CBASS FITS format version No.' 'Int'};
0136 i=i+1;
0137 
0138 fits_header(:,i)={'WORKING ' 1 'Is CBASS working?' 'Bool'};
0139 i=i+1;
0140 
0141 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0142 % JL This needs to be fixed, as regarless of filtering datatype
0143 % is comin through as 6 now.
0144 % this really should be set from the appropriate parameter in redScript.red
0145 
0146 %datatype = size(d.antenna0.receiver.data,2);
0147 %if (datatype==8)
0148 %  filter_mode=0;
0149 %elseif (datatype==4)
0150 %  filter_mode=1;
0151 %else
0152 %  error('writeFitsMap: prepareHeader  - size of d.antenna0.receiver.data neither 8 nor 4, as expected');
0153 %end
0154 %fits_header(:,i)={'FILTERED' filter_mode 'Was backend filtering used?' 'Bool'};
0155 %i=i+1;
0156 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0157 
0158 azmin=min(d.antenna0.servo.apparent(:,1));
0159 azmax=max(d.antenna0.servo.apparent(:,1));
0160 fits_header(:,i)={'AZMIN   ' azmin '[deg] Minimum azimuth value in the observation' 'Double'};
0161 i=i+1;
0162 fits_header(:,i)={'AZMAX   ' azmax '[deg] Maximum azimuth value in the observation' 'Double'};
0163 i=i+1;
0164 
0165 
0166 % unused as yet;
0167 
0168 elmin=min(d.antenna0.servo.apparent(:,2));
0169 elmax=max(d.antenna0.servo.apparent(:,2));
0170 fits_header(:,i)={'ELMIN   ' elmin '[deg] Minimum elevation value in the observation' 'Double'};
0171 i=i+1;
0172 fits_header(:,i)={'ELMAX   ' elmax '[deg] Maximum elevation value in the observation' 'Double'};
0173 i=i+1;
0174 
0175 ramin = min(d.antenna0.servo.equa(:,1));
0176 ramax = max(d.antenna0.servo.equa(:,1));
0177 fits_header(:,i)={'RAMIN   ' ramin '[rad] Min R.A. value in the obs (J2000)' 'Double'};
0178 i=i+1;
0179 fits_header(:,i)={'RAMAX   ' ramax '[rad] Max R.A. value in the obs (J2000)' 'Double'};
0180 i=i+1;
0181 
0182 decmin = min(d.antenna0.servo.equa(:,2));
0183 decmax = max(d.antenna0.servo.equa(:,2));
0184 fits_header(:,i)={'DECMIN   ' decmin '[rad] Min R.A. value in the obs (J2000)' 'Double'};
0185 i=i+1;
0186 fits_header(:,i)={'DECMAX   ' decmax '[rad] Max R.A. value in the obs (J2000)' 'Double'};
0187 i=i+1;
0188 
0189 
0190 
0191 
0192 return;
0193 end