read_skymodel

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  function  [...] = read_skymodel(xxx)

   Dec 11, 2012 (MAS) - Created.

   Reads de Oliveira Costa's (2008) sky model for a given
   ra and dec time series.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function [sky_timeseries] = ...
0002     read_skymodel(equa)
0003 
0004 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0005 %  function  [...] = read_skymodel(xxx)
0006 %
0007 %   Dec 11, 2012 (MAS) - Created.
0008 %
0009 %   Reads de Oliveira Costa's (2008) sky model for a given
0010 %   ra and dec time series.
0011 %
0012 %
0013 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0014 
0015 nside = 512;
0016 
0017 % Read in the JAVA Healpix library.
0018 % h = gov.fnal.eag.healpix.PixTools;
0019 
0020 
0021 % Galactic coords.
0022 [galactic,~]=coco(equa,'j2000.0','g','r','r');
0023 
0024 
0025 % Initialize the output array.
0026 sky_timeseries = zeros(size(equa,1),1);
0027 
0028 
0029 [home installeddir] = where_am_i();
0030 skymodel_healpix = fitsread([home '/' installeddir '/rfi_tuning/deoliveiracosta_5ghz.fits'],'BinTable');
0031 skymodel_healpix = transpose(skymodel_healpix{1});
0032 skymodel_healpix = skymodel_healpix(:);
0033 
0034 
0035 % Loop and grab the data...
0036 %disp(length(sky_timeseries));
0037 %input('adfasdf');
0038 % pix = zeros(1,length(sky_timeseries));
0039 % for i=1:length(sky_timeseries)
0040 % %    disp(size(galactic));
0041 % %    disp(i);
0042 % %    disp(h.ang2pix_ring(512,pi/2-galactic(i,2),galactic(i,1)));
0043 % %      sky_timeseries(i) = ...
0044 % %          skymodel_healpix(1+h.ang2pix_ring(512,pi/2-galactic(i,2),galactic(i,1)));
0045 %     pix(i) = 1+h.ang2pix_ring(512,pi/2-galactic(i,2),galactic(i,1));
0046 % end
0047 
0048 
0049 i_vec = zeros(1,length(sky_timeseries));
0050 j_vec = zeros(1,length(sky_timeseries));
0051 p_vec = zeros(1,length(sky_timeseries));
0052 
0053 z_vec = cos(pi/2 - galactic(:,2));
0054 
0055 
0056 north_pix = z_vec > 2/3;
0057 equa_pix = abs(z_vec) <= 2/3;
0058 south_pix = z_vec < -2/3;
0059 
0060 i_vec(north_pix) = round(sqrt(3*(1-z_vec(north_pix))) * nside);
0061 i_vec(equa_pix) = round(3/2*nside*(4/3 - z_vec(equa_pix)));
0062 i_vec(south_pix) = 4*nside - round(sqrt(3*(1-abs(z_vec(south_pix)))) * nside);
0063 
0064 % i_vec_lower = i_vec - 1;
0065 % i_vec_higher = i_vec + 1;
0066 
0067 
0068 i_vec(i_vec < 1) = 1;
0069 i_vec(i_vec > 4*nside-1) = 4*nside-1;
0070 % i_vec_lower(i_vec_lower < 1) = 1;
0071 % i_vec_lower(i_vec_lower > 4*nside-1) = 4*nside-1;
0072 % i_vec_higher(i_vec_higher < 1) = 1;
0073 % i_vec_higher(i_vec_higher > 4*nside-1) = 4*nside-1;
0074 
0075 
0076 
0077 
0078 north_i = i_vec < nside;
0079 equa_n_i = i_vec >= nside & i_vec <= 2*nside;
0080 equa_s_i = i_vec > 2*nside & i_vec <= 3*nside;
0081 south_i = i_vec > 3*nside;
0082 
0083 equa_i = equa_s_i | equa_n_i;
0084 j_vec(north_i) = ...
0085     round(2/pi * i_vec(north_i) .* mod(galactic(north_i,1)',2*pi) + 1/2);
0086 j_vec(equa_n_i) = ...
0087     round(2/pi * nside * mod(galactic(equa_n_i,1)',2*pi) + mod(i_vec(equa_n_i)-nside+1,2)/2);
0088 j_vec(equa_s_i) = ...
0089     round(2/pi * nside * mod(galactic(equa_s_i,1)',2*pi) + mod(3*nside-i_vec(equa_s_i)+1,2)/2);
0090 j_vec(south_i) = ...
0091     round(2/pi * (4*nside - i_vec(south_i)) .* mod(galactic(south_i,1)',2*pi) + 1/2);
0092 
0093 
0094 p_vec(north_i) = j_vec(north_i)-1 + 2*i_vec(north_i) .* (i_vec(north_i)-1);
0095 p_vec(equa_i) = 2*nside*(nside-1) + 4*nside*(i_vec(equa_i)-nside) + (j_vec(equa_i)-1);
0096 p_vec(south_i) = j_vec(south_i)-1 + ...
0097     12*nside^2 - 2*(4*nside - i_vec(south_i)) .* ((4*nside - i_vec(south_i))-1) - 4*(4*nside-i_vec(south_i));
0098 
0099 p_vec = p_vec + 1;
0100 
0101 
0102 
0103 % figure;
0104 % plot(pix,p_vec-pix,'.');
0105 %
0106 % input('adfads');
0107 %
0108 % figure;
0109 % plot(j_vec);
0110 %
0111 % figure;
0112 % plot(z_vec);
0113 %
0114 % figure;
0115 % plot(i_vec);
0116 % hold on;
0117 % plot(north_pix*100,'r');
0118 % plot(equa_pix*100+150,'r');
0119 % plot(south_pix*100+300,'r');
0120 
0121 sky_timeseries = skymodel_healpix(p_vec);
0122 
0123 end