Fake_Raster

 Input matrix should be of the form
 (I,Q,U)(110000   = (I1,I2,Q1,Q2,U1,U2)
        001100
        000011)

0001 function [dfake dobs]= Fake_Raster(input_mangle_matrix,mangle_name,dfake)
0002 
0003 % Input matrix should be of the form
0004 % (I,Q,U)(110000   = (I1,I2,Q1,Q2,U1,U2)
0005 %        001100
0006 %        000011)
0007 
0008 % Make fake TauA data
0009 
0010 [home,installeddir] = where_am_i();
0011 addpath([home,'/',installeddir,'/Angelas_Raster_Code/']);
0012 maindir= [home,'/',installeddir,'/'];
0013 %datadir = 'n_raster_elliptical_newcal';
0014 %datadir_fits = 'n_raster_elliptical_newcal_fits';
0015 %datadir_fits = 'Raster_Tests_Jul14'
0016 datadir='Raster_Tests_Jul14'
0017 datadir_fits = ['n_raster_elliptical_newcal_',mangle_name,'_jul14']
0018 
0019 %Load in the pipelined data as a template
0020 %TauA
0021 input_date='17-Jul-2013';
0022 pipelined_file='17-Jul-2013:13:53:34_taua_I_pipelined';
0023 if(~exist('dfake'));
0024 load(['./',datadir,'/',pipelined_file,'.mat'])
0025 cal_source='TauA';
0026 source_name='TauA';
0027 %%
0028 % Need to calculate distance from TauA at each time point
0029 % Assume flat sky for now as small angles
0030 % need the equatorial coordinates in J2000 (TauA stored in J2000 coords)
0031 
0032 disp('Getting Equatorial co-ordinates');
0033 d = calcRADECJ2000(d)
0034 
0035 %   %if(~isfield(d.antenna0.servo, 'equa'))
0036 %
0037 %   % we do not need to calculate the ra/dec until we get ready to write things
0038 %   % out.
0039 %   display('Calculating RA/DEC');
0040 %   long=-118.2822;
0041 %   lat=37.2339;
0042 %
0043 %   az = d.antenna0.servo.apparent(:,1);
0044 %   el = d.antenna0.servo.apparent(:,2);
0045 %   jd=mjd2jd(d.antenna0.receiver.utc);
0046 %   [equa] = horiz_coo([pi/180*(az) pi/180*(el)],jd,[pi/180*(long) ...
0047 %     pi/180*(lat)],'e');
0048 %   d.antenna0.servo.equa=equa;
0049 %
0050 %   clear az;
0051 %   clear el;
0052 %   clear equa;
0053 %
0054 % end
0055 
0056 RA = rad2deg(d.antenna0.servo.equa(:,1)); % J2000
0057 DEC = rad2deg(d.antenna0.servo.equa(:,2)); % J2000
0058 
0059 [RA_tau DEC_tau] = sourcePos(source_name); % degrees
0060 
0061 RA_off = RA-RA_tau;
0062 DEC_off = DEC-DEC_tau;
0063 RA_off = RA_off.*(cosd(DEC)); % scale RA to declination
0064 
0065 distance = sqrt(RA_off.^2 + DEC_off.^2);
0066 
0067 %
0068 % Now make the fake data
0069 % get model of source
0070 
0071 mjd=tstr2mjd('17-Jul-2014:13:53');
0072 [sourceNum sourceFlux fluxErr polang polfrac] = ...
0073     sourceCorrespondance({cal_source}, mjd);
0074 
0075 % Want to keep things in K (rather than Jy) so convert this sourceFlux:
0076 % read in telescope parameters from /constants/telescop_constants.m
0077 [antenna_no antenna_name] = identifyTelescope(d);
0078 telescope_constants;
0079 Tsrc = sourceFlux.*(dishArea*apEff)./(2*k_b*1e26);
0080 model = [Tsrc polfrac polang]
0081 
0082 % Need a Gaussian beam
0083 
0084 fwhm = 0.74 % degrees
0085 sigma  = fwhm/(2*sqrt(2*log(2)))
0086 beam = exp((-distance.^2)./(2*sigma^2));
0087 
0088 
0089 % Calculate PA uses DEC at current epoch
0090 % Also need JD
0091     
0092 LAT_D = antenna_latitude(antenna_no)  % degrees
0093 LON_D = antenna_longitude(antenna_no) % degrees E
0094 
0095 JD = mjd2jd(d.antenna0.receiver.utc);
0096 LST=lst(JD,deg2rad(LON_D),'m'); % Uses East longitude in radians
0097                                 % Gives LST in fraction of days
0098 
0099 % Temporarily calculate current RA DEC
0100 d = calcRADECcurrent(d); 
0101 RA_current = d.antenna0.servo.equa(:,1);  %radians
0102 DEC_current = d.antenna0.servo.equa(:,2); %radians
0103 
0104 [PA]=parallactic_angle((RA_current),(DEC_current),LST,deg2rad(LAT_D)); % Everything must be in radians, output is in radians
0105 PA_D = rad2deg(PA); % gets definition in same way as real data. %% CHeck sign here - 13/8/2014 - ACT
0106 
0107 % Make sure that we store J2000 in the data structure
0108 d = calcRADECJ2000(d);
0109 
0110 %%
0111 % Now overwrite .data(:,:)
0112 %d.data has [I1,Q1,U1,Q2,U2,Q3,U3,I2]
0113 %   I = I
0114 %   Q = I * p * cos(2(theta + phi))
0115 %   U = I * p * sin(2(theta + phi))
0116 %   theta is position angle from model, phi is parallactic angle
0117 
0118 dfake = d;
0119 
0120 dfake.antenna0.receiver.data(:,1) = Tsrc*beam;%I1
0121 dfake.antenna0.receiver.data(:,8) = Tsrc*beam;%I2
0122 
0123 dfake.antenna0.receiver.data(:,2) = Tsrc*beam*polfrac.*cosd(2*(polang+PA_D)); %Q1
0124 dfake.antenna0.receiver.data(:,4) = Tsrc*beam*polfrac.*cosd(2*(polang+PA_D)); %Q2
0125 dfake.antenna0.receiver.data(:,6) = Tsrc*beam*polfrac.*cosd(2*(polang+PA_D)); %Q3
0126 
0127 dfake.antenna0.receiver.data(:,3) = Tsrc*beam*polfrac.*sind(2*(polang+PA_D)); %U1
0128 dfake.antenna0.receiver.data(:,5) = Tsrc*beam*polfrac.*sind(2*(polang+PA_D)); %U2
0129 dfake.antenna0.receiver.data(:,7) = Tsrc*beam*polfrac.*sind(2*(polang+PA_D)); %U3
0130 %%
0131 % Now write this out as a new fits file for descart
0132 %datadir_fits = 'n_raster_elliptical_newcal_fake_jul14';
0133 %datadir_fits = 'Raster_Tests_Jul14'
0134 %  For old data, need to add day night flag
0135 % Set up a daynight flag
0136 % sun within  sunlim degrees
0137 [dist az el]      = calcSourceDist(dfake, 'sun');
0138 dfake.antenna0.servo.solarDist = dist;
0139 az = wrap360(az);
0140 raz=round(60*az);
0141 raz(raz >= 21600) = 0;
0142 dfake.flags.dayNight = ones(size(az),'uint8');
0143 %datadir_fits = ['n_raster_elliptical_newcal_fake0_jul14']
0144 %dred = reduceData(dfake,'redScript_newcal_raster_2014_fits.red',0,datadir_fits)
0145 datadir_fits = ['n_raster_elliptical_newcal_',mangle_name,'_jul14']
0146 end
0147 %%
0148 % Now mangle the data and calibrate
0149 % remember how we calibrate in polcal2.m
0150 %(I,Q,U)(110000   = (I1,I2,Q1,Q2,U1,U2)
0151 %        001100
0152 %        000011)
0153 
0154 % Uncomment below or give mangle matrix and name as a variable
0155 
0156 %mangle3
0157 %mangle3 = [  0.9,  0.9,  0.0,  0.0, 0.0, 0.0;...
0158 %           -0.2, -0.3,  0.8,  0.7, 0.5, 0.6;...
0159 %           -0.3, -0.2, -0.4, -0.5, 0.9, 0.8];
0160 
0161 %fake0
0162 %mangle = [1,1,0,0,0,0;...
0163 %           0,0,1,1,0,0;...
0164 %           0,0,0,0,1,1];      %
0165 %
0166 %mangle1
0167 %mangle = [0.5,0.5,0,0,0,0;...
0168 %          0,0,1,1,0,0;...
0169 %          0,0,0,0,1,1];      %
0170 
0171 %mangle2a
0172 % mangle = [1,1,0,0,0,0;...
0173 %           0,0,0.5,0.5,0,0;...
0174 %           0,0,0,0,0.5,0.5];      %
0175 %
0176 %mangle4
0177 % mangle = [1,1,0,0,0,0;...
0178 %           0,0,0.3,0.3,0,0;...
0179 %           0,0,0,0,0.8,0.8];      %
0180  %mangle5
0181  %mangle = [1,1,0,0,0,0;...
0182  %          0,0,0.3,0.3,0.4,0.4;...
0183  %          0,0,0.5,0.5,0.8,0.8];      %
0184 %
0185 
0186 mangle=input_mangle_matrix;
0187 % Re-fill dfake with 'observed' mangled data (I,Q,U) above
0188 
0189 trueS = [dfake.antenna0.receiver.data(:,[1,2,3])];
0190     
0191 obsS = trueS*mangle; % (I1,I2,Q1,Q2,U1,U2)
0192 
0193 % Now copy this into dobs and write out as fits ready to go into descart
0194 dobs=dfake;
0195 dobs.antenna0.receiver.data(:,1) = obsS(:,1); %I1
0196 dobs.antenna0.receiver.data(:,2) = obsS(:,3); %Q1
0197 dobs.antenna0.receiver.data(:,3) = obsS(:,5); %U1
0198 dobs.antenna0.receiver.data(:,4) = obsS(:,4); %Q2
0199 dobs.antenna0.receiver.data(:,5) = obsS(:,6); %U2
0200 dobs.antenna0.receiver.data(:,6) = (obsS(:,3)+obsS(:,4))/2; %Q3
0201 dobs.antenna0.receiver.data(:,7) = (obsS(:,5)+obsS(:,6))/2; %U3
0202 dobs.antenna0.receiver.data(:,8) = obsS(:,2); %I2
0203 %%
0204 % Now write this out as a new fits file for descart
0205 %datadir_fits = ['n_raster_elliptical_newcal_',mangle_name,'_jul14'];
0206 
0207 % Save the mangeld observation in a form that can then be read in by the
0208 % raster_analysis and calibration code:
0209 
0210 %  For old data, need to add day night flag
0211 % Set up a daynight flag
0212 % sun within  sunlim degrees
0213 [dist az el]      = calcSourceDist(dobs, 'sun');
0214 dobs.antenna0.servo.solarDist = dist;
0215 az = wrap360(az);
0216 raz=round(60*az);
0217 raz(raz >= 21600) = 0;
0218 dobs.flags.dayNight = ones(size(az),'uint8');
0219 %d=dobs;
0220 dred = reduceData(dobs,'redScript_newcal_raster_2014_fits.red',0,datadir_fits)
0221 
0222 reduced_data_file = strcat([pipelined_file,'_',mangle_name,'.mat'])
0223 save([maindir,datadir,'/',reduced_data_file],'dobs','-v7.3')
0224 
0225 % Save the mangled matrix in the correct fom for descart
0226 %[I1 I2 Q1 Q2 U1 U2]
0227 
0228 descart_file=mangle(:,[1 2 3 5 4 6])';
0229 
0230 % Everything (fits and leakage matrix gets saved into
0231 % cbass_analysis/data_dir
0232 save([maindir,datadir,'/leakage_',mangle_name,'.txt'],'descart_file','-ascii','-v7.3')
0233 fitspath=[maindir,'reduc/',datadir_fits,'/',input_date,'*/fits/*.fits'] ;
0234 newfitsfile=[maindir,datadir,'/',input_date,'_',mangle_name,'_nodetrend.fits']
0235 system(['cp ',fitspath,' ',newfitsfile ])
0236 %whos
0237 
0238 % Make descart map with no calibration
0239 cd ([datadir])
0240 system(['ls ',input_date,'_',mangle_name,'_nodetrend.fits > files.txt'])
0241 % JL HACK - TO PREVENT A DESCART link error on elephant
0242   current_ld_lib_path = getenv('LD_LIBRARY_PATH');
0243   setenv('LD_LIBRARY_PATH',['/usr/local/shared/cfitsio/lib/:',current_ld_lib_path]);
0244 
0245 system(['descart_cbass ./default_params.ini do_polarization=T naive_mode=T nside=1024 subtractMeans=F leakage_matrix=leakage_fake0.txt  output_filename=',mangle_name,'_map_noC.fits'])
0246 % make map with leakage matrix applied
0247 system(['descart_cbass ./default_params.ini do_polarization=T naive_mode=T nside=1024 subtractMeans=F leakage_matrix=leakage_',mangle_name,'.txt  output_filename=',mangle_name,'_map.fits'])
0248 % Make the images using python
0249 system(['python plot_maps.py ',mangle_name,'_map.fits ',mangle_name,'_map.png'])
0250 system(['python plot_maps.py ',mangle_name,'_map_noC.fits ',mangle_name,'_map_noC.png'])
0251 cd ([maindir])