Description of decomposeSky

0001 function [Ivec Pvec Ivar Pvar] = decomposeSky(d, Nside, type)
0002 
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 %  function [Ivec Pvec Ivar Pvar] = decomposeSky(d, Nside, type)
0006 %
0007 %
0008 %    calculates the time-stream data for the sky-map.
0009 %
0010 %     d is our data stream
0011 %     Nside is the desired resolution of the map (right now 1024 or 256)
0012 %           default is 1024
0013 %     type:  1 -- closest pixel
0014 %            2 -- interpolation from 4 nearest pixels.  (default)
0015 %
0016 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0017 
0018 if(nargin<2)
0019   Nside = 1024;
0020   type = 2;
0021 end
0022 if(nargin<3)
0023   type = 2;
0024 end
0025 
0026 Ivec = [];
0027 Pvec = [];
0028 Ivar = [];
0029 Pvar = [];
0030 
0031 switch Nside
0032   case 1024
0033     load skymap_nside1024_20140101.mat
0034 
0035   case 256
0036     load skymap_nside256_20140101.mat
0037     
0038   otherwise
0039     error('decomposeSky:: Sky Map of desired Nside not available, use makeMapArrays.m');
0040 end
0041 
0042 
0043 % what we really want is to switch our data into galactic coordinates.  so
0044 % first we go from az/el to radec.
0045 display('decomposeSky:: Calculating RA/DEC');
0046 long=d.antenna0.tracker.siteFiducial(1,1);
0047 lat= d.antenna0.tracker.siteFiducial(1,2);
0048 az = d.antenna0.servo.apparent(:,1);
0049 el = d.antenna0.servo.apparent(:,2);
0050 jd=mjd2jd(d.antenna0.receiver.utc);
0051 [radec_d] = horiz_coo([pi/180*(az) (pi/180*(el))],jd,[pi/180*(long) ...
0052       pi/180*(lat)],'e');
0053 yearoff = (mean(d.array.frame.utc) - date2mjd(2000, 01, ...
0054     01,0,0,0))/365.25;
0055 yearoff = round( (2000+yearoff)*10)/10;
0056 txt = sprintf('[gal_d] = coco(radec_d, ''j%4.1f'', ''g'', ''r'', ''d'');', yearoff);
0057 eval(txt);
0058 [gal_d, TotRot]=coco(radec_d,'j2000.0','g','r','d');
0059 latd = gal_d(:,2);
0060 longd = gal_d(:,1);
0061 longd(longd<0) = longd(longd<0) + 360;
0062 
0063 
0064 Ivec = zeros(size(latd));
0065 Pvec = zeros(size(latd));
0066 
0067 newVal = 10;
0068 switch type
0069   case 1
0070     %% FASTEST METHOD.  JUST FIND THE CLOSEST  34.7s for 719900 samples
0071     tic
0072     for m=1:length(latd)
0073       amtdone = m./length(latd)*100;
0074       if(amtdone > newVal)
0075     display(sprintf('decomposeSky:: done with %d%% of data', newVal));
0076     newVal = newVal + 10;
0077       end
0078       f = binsearch(uniquelat, latd(m));
0079       ff = binsearch(longSplit{f}, longd(m));
0080       Ivec(m) = Isplit{f}(ff);
0081       Pvec(m) = Psplit{f}(ff);
0082     end
0083     display('decomposeSky::  done with 100% of data');
0084     toc
0085 
0086   case 2
0087     tic
0088     %% FIND THE NEAREST VALUES AS WELL:
0089     Ivar = zeros(size(latd));
0090     Pvar = Ivar;
0091     Ivals = Ivar;
0092     Pvals = Ivar;
0093     allI = zeros(length(latd),4);
0094     allP = zeros(length(latd),4);
0095     varI = zeros(length(latd),4);
0096     varP = zeros(length(latd),4);
0097     
0098     longVec = zeros(length(latd),4);
0099     latVec = zeros(length(latd),4);
0100     display('decomposeSky:: Decomposing sky model into time stream');
0101     for m=1:length(latd)
0102       amtdone = m./length(latd)*100;
0103       if(amtdone > newVal)
0104     display(sprintf('decomposeSky:: done with %d%% of data', newVal));
0105     newVal = newVal + 10;
0106       end
0107       f = binsearch(uniquelat, latd(m));
0108       if(latd(m) < uniquelat(f))
0109     f2 = f-1;
0110     if(f2==0)
0111       f2 = length(uniquelat);
0112     end
0113       else
0114     f2 = f+1;
0115     if(f2>length(uniquelat))
0116       f2 = 1;
0117     end
0118       end
0119       % toc
0120       %  tic
0121       ff = binsearch(longSplit{f}, longd(m));
0122       if(longd(m) < longSplit{f}(ff))
0123     ff_2 = ff - 1;
0124     if(ff_2==0)
0125       ff_2 = length(longSplit{f});
0126     end
0127       else
0128     ff_2 = ff + 1;
0129     if(ff_2 > length(longSplit{f}))
0130       ff_2 = 1;
0131     end
0132       end
0133       ff2 = binsearch(longSplit{f2}, longd(m));
0134       if(longd(m) < longSplit{f2}(ff2))
0135     ff2_2 = ff2 - 1;
0136     if(ff2_2==0)
0137       ff2_2 = length(longSplit{f2});
0138     end
0139       else
0140     ff2_2 = ff2 + 1;
0141     if(ff2_2 > length(longSplit{f2}))
0142       ff2_2 = 1;
0143     end
0144       end
0145       allI(m,:) = [Isplit{f}(ff); Isplit{f}(ff_2); Isplit{f2}(ff2); Isplit{f2}(ff2_2)];
0146       allP(m,:) = [Psplit{f}(ff); Psplit{f}(ff_2); Psplit{f2}(ff2); Psplit{f2}(ff2_2)];
0147       varI(m,:) = [Ivarsplit{f}(ff); Ivarsplit{f}(ff_2); Ivarsplit{f2}(ff2); Ivarsplit{f2}(ff2_2)];
0148       varP(m,:) = [Pvarsplit{f}(ff); Pvarsplit{f}(ff_2); Pvarsplit{f2}(ff2); Pvarsplit{f2}(ff2_2)];
0149       longVec(m,:) = [longSplit{f}(ff); longSplit{f}(ff_2); ...
0150     longSplit{f2}(ff2); longSplit{f2}(ff2_2)];
0151       latVec(m,:)  = [uniquelat(f); uniquelat(f); uniquelat(f2); uniquelat(f2)];
0152     end
0153     deltaLong = repmat(longd, [1,4]) - longVec;
0154     sinTerm   = repmat(sind(latd), [1,4]).*sind(latVec);
0155     cosTerm   = repmat(cosd(latd), [1,4]).*cosd(latVec);
0156     dist = acos(sinTerm + cosTerm.*cosd(deltaLong));
0157     gauss_dist = exp(-dist.^2./(2*(2/180*pi)^2));
0158     Ivec = sum(allI.*gauss_dist*1./varI,2)./(sum(gauss_dist./varI,2));
0159     Pvec = sum(allP.*gauss_dist./varP,2)./(sum(gauss_dist./varP,2));
0160     %% STILL NEED TO PROPAGATE THE ERRORS HERE.
0161     varFactor = sum( gauss_dist.^2,2) ./( (sum(gauss_dist,2)).^3 - ...
0162     sum(gauss_dist,2).*sum(gauss_dist.^2,2));
0163     Ivar  = varFactor.*sum(gauss_dist.*(repmat(Ivals, [1 4]) - ...
0164     allI).^2,2);% + sum(varI,2);
0165     Pvar  = varFactor.*sum(gauss_dist.*(repmat(Ivals, [1 4]) - ...
0166     allP).^2,2);% + sum(varP,2);
0167     display('decomposeSky:: done with 100% of data');
0168     toc
0169 end    
0170 
0171 return;
0172
decomposeSky

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SOURCE CODE
