MakeGroundTemplate

 function which reads all the HEALPix fits file names in the text file string name 'file' 
 and makes daily Ground Templates. If write = 1 templates get written out to
'name'I/Q/U.txt

0001 function MakeGroundTemplate(files, write, name, nside, elevation_list)
0002 % function which reads all the HEALPix fits file names in the text file string name 'file'
0003 % and makes daily Ground Templates. If write = 1 templates get written out to
0004 %'name'I/Q/U.txt
0005 
0006 [home,installeddir]=where_am_i();
0007 location = [home,'/',installeddir]; 
0008 
0009 locationDes = [location,'/descart/map_making/descart_cbass'];
0010 
0011 %TO DO - do not hardwire this - make it react to the appropraite value in the redscript.red
0012 locationGS = [location,'/webpage_logging/my_combinedDS.ini'];
0013 
0014 % JL HACK - TO PREVENT A DESCART link error on elephant/mammoth
0015 current_ld_lib_path = getenv('LD_LIBRARY_PATH');
0016 %setenv('LD_LIBRARY_PATH',['/usr/local/shared/cfitsio/lib/:',current_ld_lib_path]);
0017 setenv('LD_LIBRARY_PATH',['/usr/local/shared/cfitsio-3340-intel/lib/:',current_ld_lib_path])
0018 
0019 % Makes the fits files by running descart - your ini file needs to be reading in files.txt
0020 unix([locationDes,' ',locationGS,' output_filename=',name, 'Ground.fits']);
0021 
0022 setenv('LD_LIBRARY_PATH',current_ld_lib_path);
0023 
0024 % read in the fits files that was just generated by descart...
0025 
0026 fitsfile = strcat(name, 'Ground.fits');
0027 data = fitsread(fitsfile, 'BinTable'); 
0028 sI = data{1}';
0029 sQ = data{2}';
0030 sU = data{3}'; %'
0031 sI = reshape(sI, 1, 12*nside*nside); % Healpix formula for Npix from nside.
0032 sQ = reshape(sQ, 1, 12*nside*nside);
0033 sU = reshape(sU, 1, 12*nside*nside);
0034 runner = 1; 
0035 
0036 
0037 % cannot deal with healpix bad values, will use NaNs instead
0038 
0039 outI = find(sI < -1000.);
0040 sI(outI) = NaN;
0041 outQ = find(sQ < -1000.);
0042 sQ(outQ) = NaN;
0043 outU = find(sU < -1000.);
0044 sU(outU) = NaN;
0045 
0046 % Determine Az and El
0047 % This needs to MEALPix pix2ang function to be installed.
0048 tp_cells = pix2ang(nside, 'nest', false);
0049 
0050 % convert cells array into regular array for convenience
0051 
0052 %tp =zeros(2,length(tp_cells));
0053 
0054 theta=zeros(1,length(tp_cells));
0055 phi=zeros(1,length(tp_cells));
0056 
0057 for i=1:length(tp_cells)
0058   theta(i)=tp_cells{i}(1);
0059   phi(i)=tp_cells{i}(2);
0060   %  tp(1,i) = the_cells{i}(1);
0061   %  tp(2,i) = the_cells{i}(2);
0062 end;
0063 
0064 %theta = tp(1,:);
0065 %phi = tp(2,:);
0066 
0067 az = phi * 180./pi;
0068 el = -theta * 180./pi + 90;
0069 
0070 
0071 % Make a string version of the elevation list.
0072 for i = 1:length(elevation_list)
0073   elevation_string_list(i) =cellstr(sprintf('%i',elevation_list(i))); 
0074 end
0075 
0076 
0077 %>> class(elevation_string_list{1})
0078 %char
0079 %
0080 %>> class(elevation_string_list(1))
0081 %cell
0082 
0083 
0084 %elevation_slop = 1.0; % Get all elevations within 1 degree of each element in the elevation list.
0085 
0086 % JL 16 Jan 2014 elevation value of 1.0 does not pick up el =67 values (does for 77 values (increased).
0087 % Should think of a better solution for this.
0088 
0089 elevation_slop = 1.5; % Get all elevations within 1 degree of each element in the elevation list.
0090 
0091 for i = 1:length(elevation_list)
0092 
0093 
0094   el_pix = find(el >  (elevation_list(i) - elevation_slop/2) & el < (elevation_list(i) + elevation_slop/2) );
0095   disp( ['MakeGroundTemplate:: Found ',num2str(length(el_pix)),' pixels at elevation ',elevation_string_list{i}]);
0096 
0097   if (length(el_pix)==0)
0098     disp(['MakeGroundTemplate:: WARN:: NO SUITABLE PIXELS FOUND AT THIS ELEVATION, GROUND SUB WILL FAIL.']);
0099   end
0100  
0101   pixIs = sI(:,el_pix);
0102   pixQs = sQ(:,el_pix); 
0103   pixUs = sU(:,el_pix);
0104   pixIs = pixIs - nanmedian(pixIs);
0105   pixQs = pixQs - nanmedian(pixQs);
0106   pixUs = pixUs - nanmedian(pixUs);
0107   az_s = az(el_pix);
0108   if write == 1
0109 
0110     post_fix_stub = [elevation_string_list{i},'gs.txt'] 
0111     disp(['MakeGroundTemplate:: Preparing to write ground templates to ',location,'/constants/',name,'_I',post_fix_stub]);
0112 
0113     post_fix_stub = [elevation_string_list{i},'gs.txt'] 
0114     nameI = strcat(location,'/constants/',name, '_I',post_fix_stub);
0115     nameQ = strcat(location,'/constants/',name, '_Q',post_fix_stub);
0116     nameU = strcat(location,'/constants/',name, '_U',post_fix_stub);
0117 
0118     fidI = fopen(nameI,'w');
0119     for s = 1:length(pixIs)
0120       fprintf(fidI, '%8.4f %7.4f \n', az_s(s), pixIs(s));
0121     end
0122     fclose(fidI);   
0123 
0124     fidQ = fopen(nameQ,'w');
0125     for s = 1:length(pixQs)
0126       fprintf(fidQ, '%8.4f %7.4f \n', az_s(s), pixQs(s));
0127     end
0128     fclose(fidQ);   
0129 
0130     fidU = fopen(nameU,'w');
0131     for s = 1:length(pixUs)
0132       fprintf(fidU, '%8.4f %7.4f \n', az_s(s), pixUs(s));
0133     end
0134     fclose(fidU);   
0135   
0136   end
0137 
0138 end
0139