; ----------------------------------------------------------------------------- ; ; Copyright (C) 1997-2005 Krzysztof M. Gorski, Eric Hivon, Anthony J. Banday ; ; ; ; ; ; This file is part of HEALPix. ; ; HEALPix is free software; you can redistribute it and/or modify ; it under the terms of the GNU General Public License as published by ; the Free Software Foundation; either version 2 of the License, or ; (at your option) any later version. ; ; HEALPix is distributed in the hope that it will be useful, ; but WITHOUT ANY WARRANTY; without even the implied warranty of ; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ; GNU General Public License for more details. ; ; You should have received a copy of the GNU General Public License ; along with HEALPix; if not, write to the Free Software ; Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA ; ; For more information about HEALPix see http://healpix.jpl.nasa.gov ; ; ----------------------------------------------------------------------------- pro list_all_ttypes, file, ttype, extn, nbefore case datatype(file) of 'STR': begin ; fits file fits_info, file, /silent, n_ext = n_ext ; junk first slot, to match structure generated by read_fits_s,/merge ttype = 'JUNK' extn = 0 nbefore = 0 nc_old = 0 ; image, if there is one hdr = headfits(file,exten=0) szimg = long(sxpar(hdr, 'NAXIS')) if (szimg gt 0) then begin ttype = ['JUNK', 'IMG'] extn = [-1,-1] nbefore = [0,0] nc_old = 1 endif for i=1,n_ext do begin hdr = headfits (file, exten=i) nc = long(sxpar(hdr, 'TFIELDS')) ttype1 = sxpar(hdr, 'TTYPE*', count = nt) if (nt lt nc) then ttype1 = strtrim(string(intarr(nc)+nc_old, form='(i2)'),2)+'s' ttype = [ttype, ttype1] extn = [extn, replicate(i,nc) - 1] nbefore = [nbefore, replicate(nc_old, nc)] nc_old = nc_old + nc endfor end 'STC': begin ; structure ttype = tag_names(file) nt = n_tabs(file) extn = replicate(0,nt) end endcase return end function extract_map_from_stc, stc, select_in, select=select, value=select_name, error=error, tonan=tonan ; ----------------------------------------------------------- ; stc : IN, structure ; select_in, IN, integer or string: column identifier ; select, OUT kwd integer actual column number ; value OUT kwd string actual column name ; error OUT kwd integer error code ; tonan IN kwd integer if set, map bad pixels to NaN ; ----------------------------------------------------------- routine = 'extract_map_from_stc' header = stc.(0) if (error ne 0) then return,0 init_healpix bad_value = !healpix.bad_value ; -1.63750e+30 junk = getsize_fits(header,/header,type = fits_type, nside = nside_head, tags=tag_names(stc)) ; coverage = strtrim(strupcase(sxpar(header,'OBJECT',Count=flag_cov)),2) ; nside_head = long(sxpar(header,'NSIDE',Count=flag_nside)) if (fits_type eq 3) then begin ; partial sky coverage: create a full sky map, setting missing pixels to BAD select_pix = index_word(tag_names(stc),'PIXEL',value=select_name,err=error2) if undefined(select_in) then begin select_map = index_word(tag_names(stc),'SIGNAL',value=select_name,err=error1) endif else begin select_map = index_word(tag_names(stc),select_in,value=select_name,err=error1) endelse if (error1 ne 0 or error2 ne 0 or select_map le 0 or select_pix le 0) then begin print,'wrong select in '+routine print,select_in,tag_names(stc) error = 1 return,0 endif npix_head = nside2npix(nside_head) if keyword_set(tonan) then bad_val = !values.f_nan else bad_val = bad_value map = make_array(npix_head,/float,value=bad_val) map[stc.(select_pix)] = stc.(select_map) select = select_map endif else if (fits_type eq 0 or fits_type eq 2) then begin ; full sky coverage (in image or binary table) if undefined(select_in) then begin select = index_word(tag_names(stc),1 ,value=select_name,err=error) endif else begin select = index_word(tag_names(stc),select_in,value=select_name,err=error) endelse if (error ne 0 or select eq -1) then begin print,'wrong select in '+routine print,select_in,tag_names(stc) error = 1 ; return,0 message,'Abort' endif if (keyword_set(tonan)) then begin bad = where ( stc.(select) le (bad_value*.9), nbad) map = stc.(select) if (nbad gt 0) then map[bad] = !values.f_nan endif else begin map = stc.(select) endelse endif else begin ; print,'Unrecognised format' ; stop message,'Unrecognised format' endelse return, map end ;================================================================================================== pro loaddata_healpix, file_in, select_in,$ data, pol_data, pix_type, pix_param, do_conv, do_rot, coord_in, coord_out, eul_mat, title_display, sunits, $ SAVE=save,ONLINE=online,NESTED=nested_online,UNITS=units,COORD=coord,ROT=rot_ang,QUADCUBE=quadcube,LOG=log, $ ERROR=error, FLIP=flip, POLARIZATION=polarization, FACTOR=factor_u, OFFSET = offset_u ;+ ; LOADDATA_HEALPIX ; ; IN: ; file_in, select_in ; OUT: ; data, pol_data, pix_type, pix_param, do_conv, do_rot, coord_in, coord_out, eul_mat, ; title_display, sunits ; KEYWORDS ; save, online, nested, units, coord, rot, quadcube, log, polarization, ; flip, factor ; OPTIONAL OUPUT ; error ; ; Aug 2003 : accepts ONLINE polarized data ; Feb 2006: edited to deal with multi-extension polarized cut sky files ; Apr 2006: edited to deal with arrays (online option) with more than 3 rows ;- factor = defined(factor_u) ? factor_u : 1. offset = defined(offset_u) ? offset_u : 0. ; ---------------------------------------------- ; check consistency of parameters ; ---------------------------------------------- error = 0 routine='loaddata_healpix' kw_save = KEYWORD_SET(save) & kw_online = KEYWORD_SET (online) if ((kw_save) and (kw_online)) then begin print,routine+' : /ONLINE and /SAVE can not be used together' error=1 return endif if (datatype(file_in) NE 'STR' AND kw_save) then begin print,routine+' : an external file is expected with /SAVE' error=1 return endif if (datatype(file_in) NE 'STR' AND NOT kw_online) then begin print,routine+' : a file name is expected or /ONLINE is missing' error=1 return endif if (datatype(file_in) EQ 'STR' AND kw_online) then begin print,routine+' : /ONLINE can not be used with an external file' error=1 return endif if (datatype(file_in) EQ 'STR') then begin ; looking for a file junk = FINDFILE(file_in, count=countfile) ; check its existence if (countfile eq 0) then begin print,file_in+' not found' error=1 return endif endif if keyword_set(flip) then flipconv=1 else flipconv = -1 ; longitude increase leftward by default (astro convention) do_polarization = keyword_set(polarization) if not do_polarization then polarization = 0 online_array = (kw_online and datatype(file_in) ne 'STC') online_structure = (kw_online and datatype(file_in) EQ 'STC') from_file = (not kw_online and datatype(file_in) EQ 'STR') ; ---------------------------------------------- ; reads in the FITS file or the SAVESET file ; ---------------------------------------------- print,'reading the file ...' ; in file_in, kw_save, kw_online, select_in, routine ; out data, title_display, header select_name = '' ;if undefined(select_in) then select_in = 1 t_keys = ['T','I'] q_keys = ['Q-P','Q_P','QS','Q'] u_keys = ['U-P','U_P','US','U'] do_rescale = (factor ne 1. or offset ne 0.) pol_data = 0. ;-------------------------------------------------------------- ; read the saveset file ;-------------------------------------------------------------- if (kw_save) then begin icolumn = (defined(select_in)) ? (select_in - 1) : 0 RESTORE, file_in, /VERBOSE dim2 = n_elements(data[0,*]) if (icolumn ge dim2) then message,'selected column not available in array read from saveset' data = data[*,icolumn] if (do_rescale) then begin bad_pixels = where(data le (!healpix.bad_value*0.9) or finite(data,/nan), nbad) if (factor ne 1.) then data = temporary(data) * factor if (offset ne 0.) then data = temporary(data) + factor*offset if (nbad gt 0) then data[bad_pixels] = !values.f_nan endif title_display = file_in if (do_polarization) then begin print,'No plotting of polarisation available with this kind of data set' error = 1 return endif if (keyword_set(quadcube) and defined(sky_cube)) then sixpack2vector,sky_cube, data goto, data_are_in endif ;-------------------------------------------------------------- ; non structure ONLINE data ;-------------------------------------------------------------- if (online_array) then begin title_display = ' on line processing ' dim1 = n_elements(file_in[*,0]) dim2 = n_elements(file_in[0,*]) if (do_polarization) then begin if (dim2 lt 3) then begin ; 3 rows or less -> crash print,'No plotting of polarisation available with this kind of data set' print,' expected array of shape [Npix,3], got : ',dim1,dim2 error = 1 return endif else begin if (dim2 gt 3) then begin ; more than 3 rows -> complain formdim2 = '(a,i12," x",i2)' print,'Input array has shape ',dim1,dim2,form=formdim2 print,'Expected ',dim1,3, form=formdim2 print,'WARNING: Will assume that first 3 rows are I,Q,U' endif if (do_rescale) then begin ; flag bad pixels with NaN bad_pixels = where(file_in le (!healpix.bad_value*0.9), nbad) if (nbad gt 0) then file_in[bad_pixels] = !values.f_nan endif case polarization of 1: begin select_name = 'Polarisation Amplitude' data = SQRT(file_in[*,1]^2 + file_in[*,2]^2)*factor ; amplitude = sqrt(U^2+Q^2) end 2: begin select_name = 'Polarisation Direction' ; change sign of U if astro convention data = .5*ATAN(file_in[*,2]*flipconv, file_in[*,1]) ; angle phi = 1/2 atan(U/Q) in radians end 3: begin select_name = 'Temperature + Polarisation' ; amplitude = sqrt(U^2+Q^2) & angle phi = 1/2 atan(U/Q) in radians pol_data = [[SQRT(file_in[*,1]^2 + file_in[*,2]^2)*factor], $ [.5*ATAN(file_in[*,2]*flipconv, file_in[*,1])]] ; temperature data = file_in[*,0]*factor if (offset ne 0.) then data = temporary(data) + factor*offset end endcase endelse endif else begin nbad = 0 icolumn = (defined(select_in)) ? (select_in - 1) : 0 if (icolumn ge dim2) then message,'selected column not available in array provided' if (do_rescale) then begin bad_pixels = where(file_in[*,icolumn] le (!healpix.bad_value*0.9) or finite(file_in[*,icolumn],/nan), nbad) endif data = file_in[*,icolumn] * factor if (offset ne 0.) then data = temporary(data) + factor*offset if (nbad gt 0) then data[bad_pixels] = !values.f_nan endelse goto, data_are_in endif ;-------------------------------------------------------------- ; FITS file ;-------------------------------------------------------------- if (from_file) then begin ; fits file title_display = file_in ; find out which extension to read list_all_ttypes, file_in, all_ttype_tags, ext_ttype_tag, ntags_before sel_tmp = (defined(select_in)) ? select_in : 1 iw = index_word(all_ttype_tags, sel_tmp, err= error) if (error ne 0 or iw lt 0) then begin print,'Invalid choice for field to plot: ',sel_tmp print,'Available choices: ',all_ttype_tags message,'Abort' endif ext2read = ext_ttype_tag[iw] ntags_before = ntags_before[iw] select = iw - ntags_before ; column number within extension select_name = all_ttype_tags[iw] ; name based on FITS column name ; junk = getsize_fits(file_in, type = fits_type) if (fits_type eq 3) then begin ; cut sky -> do structure file_keep = file_in if (polarization ge 1 and polarization le 3) then begin ; polarised data read_fits_s, file_keep, st_t, /merge, ext = 0 read_fits_s, file_keep, st_q, /merge, ext = 1 read_fits_s, file_keep, st_u, /merge, ext = 2 iw_signal = (index_word(all_ttype_tags, 'SIGNAL', err= error))[0] iw_pixel = (index_word(all_ttype_tags, 'PIXEL', err= error))[0] file_in = {JUNK:st_t.(0), PIXEL:st_t.(iw_pixel), TEMPERATURE:st_t.(iw_signal), $ Q_POLARISATION:st_q.(iw_signal), $ U_POLARISATION:st_u.(iw_signal)} if (n_elements(file_in.(2)) ne n_elements(file_in.(3)) or $ n_elements(file_in.(2)) ne n_elements(file_in.(4))) then begin message,'Unconsistent I, Q & U fields in cut-sky FITS file:'+file_keep endif endif else begin ; temperature only read_fits_s, file_keep, file_in, /merge, ext = ext2read select_in = select endelse online_structure = 1 goto, do_online_stc endif selectread, file_in, data, pol_data, header=header, extension = ext2read, pol = polarization, tonan=tonan, factor=factor, offset=offset,col=select ; pol_err_mess = 'can not find POLARISATION in file '+file_in case polarization of 1: select_name = 'Polarisation Amplitude' 2: select_name = 'Polarisation Direction' 3: select_name = 'Temperature + Polarisation' else: select_name = all_ttype_tags[iw] endcase goto, data_are_in endif ;-------------------------------------------------------------- ; STRUCTURE ONLINE data ;-------------------------------------------------------------- ; extract data from structure do_online_stc: if (online_structure) then begin ; structure header = file_in.(0) title_display = ' on line processing ' pol_err_mess = 'can not find POLARISATION in selected online data set' case polarization of 1: begin data_q = extract_map_from_stc(file_in, q_keys, select=select, value=select_name, error=error, /tonan) data = extract_map_from_stc(file_in, u_keys, select=select, value=select_name, error=error, /tonan) if (error ne 0) then message,pol_err_mess select_name = 'Polarisation Amplitude' data = temporary(SQRT(temporary(data)^2 + data_q^2))*factor ; amplitude = sqrt(U^2+Q^2) data_q = 0. end 2: begin data_q = extract_map_from_stc(file_in, q_keys, select=select, value=select_name, error=error, /tonan) data = extract_map_from_stc(file_in, u_keys, select=select, value=select_name, error=error, /tonan) if (error ne 0) then message,pol_err_mess select_name = 'Polarisation Direction' if (flipconv lt 0) then data = temporary(data) * flipconv ; change sign of U if astro convention data = .5*temporary(ATAN(temporary(data), data_q)) ; angle phi = 1/2 atan(U/Q) in radians data_q = 0. end 3: begin data_q = extract_map_from_stc(file_in, q_keys, select=select, value=select_name, error=error, /tonan) data_u = extract_map_from_stc(file_in, u_keys, select=select, value=select_name, error=error, /tonan) if (error ne 0) then message,pol_err_mess if (flipconv lt 0) then data_u = temporary(data_u) * flipconv ; change sign of U if astro convention ; amplitude = sqrt(U^2+Q^2) & angle phi = 1/2 atan(U/Q) in radians pol_data = [[SQRT(data_u^2 + data_q^2)*factor],[.5*(ATAN(data_u, data_q))]] if (offset ne 0.) then pol_data[*,0] = pol_data[*,0] + factor*offset data_q = 0. & data_u = 0. data = extract_map_from_stc(file_in, t_keys, select=select, value=select_name, error=error, /tonan) data = temporary(data)*factor if (offset ne 0.) then data = temporary(data) + factor*offset select_name = 'Temperature + Polarisation' end else: begin ; temperature alone data = extract_map_from_stc(file_in, select, select=select, value=select_name, error=error, tonan=do_rescale) if (datatype(data) eq 'STR') then message,'Abort: Wrong format' if (do_rescale) then begin bad_pixels = where(finite(data,/nan), nbad) data = temporary(data)*factor if (offset ne 0.) then data = temporary(data) + factor*offset if (nbad gt 0) then data[bad_pixels] = !values.f_nan endif end endcase if (error ne 0) then message,'error while reading structure' goto, data_are_in endif ;-------------------------------------------------------------- ;-------------------------------------------------------------- data_are_in: if defined(file_keep) then file_in = file_keep if (error ne 0) then return title_display = title_display + ': '+select_name if defined(header) then begin if n_elements(header) eq 1 then header =[header,' '] endif ; ---- Pixel scheme ---- ; Healpix : Nested or Ring scheme ; or Quadcube pix_type = 'R' ; Healpix ring by default IF KEYWORD_SET(nested_online) THEN pix_type = 'N' IF DEFINED(header) THEN BEGIN ordering = strtrim(SXPAR(header,'ORDERING',COUNT=flag_order),2) if (flag_order ne 0 and STRUPCASE(STRMID(ordering,0,4)) eq 'NEST') then pix_type = 'N' ENDIF if keyword_set(quadcube) then pix_type = 'Q' ; ---- units ---- sunits = '' if (datatype(units) eq 'STR') then begin ; if the units given is a string, take it sunits = units endif else begin ; otherwise find it in the file if defined(header) then begin key_unit = 'TUNIT'+STRTRIM(STRING(select,form='(i4)'),2) ; TUNITi sunits_read = SXPAR(header,key_unit,COUNT=flag_units) if (flag_units NE 0) then sunits = STRTRIM(sunits_read, 2) endif endelse if keyword_set(log) then sunits = 'Log ('+sunits+')' ;---------------------------------------------------- ; plot parameters defined by the user ;---------------------------------------------------- ; ---- input and output coord system ---- flag_coord = 0 & coord_in = 'G' IF DEFINED(header) THEN coord_in = STRUPCASE(STRMID(strtrim(SXPAR(header,'COORDSYS',COUNT=flag_coord),2),0,1)) if (flag_coord eq 0 or strtrim(coord_in) eq '') then coord_in = 'G' ; nothing in the header -> assume galactic coord_out = coord_in if N_elements(coord) EQ 1 then coord_out = STRUPCASE(STRMID(coord(0),0,1)) if N_elements(coord) EQ 2 then begin coord_in = STRUPCASE(STRMID(coord(0),0,1)) coord_out = STRUPCASE(STRMID(coord(1),0,1)) endif if coord_in EQ 'C' then coord_in = 'Q' ; cgis skyconv coding convention for celestial/equatorial if coord_out EQ 'C' then coord_out = 'Q' print,'input file : ',decode_coord(coord_in)+' coordinates' print,'plot coord : ',decode_coord(coord_out)+' coordinates' do_conv = (coord_in NE coord_out) ; ---- extra rotation ---- if defined(rot_ang) then rot_ang = ([rot_ang,0.,0.])(0:2) else rot_ang = [0., 0., 0.] ; eul_mat = euler_matrix(-rot_ang(0), -rot_ang(1), -rot_ang(2), /Deg, /ZYX) eul_mat = euler_matrix_new(rot_ang(0), -rot_ang(1), rot_ang(2), /Deg, /ZYX) do_rot = (TOTAL(ABS(rot_ang)) GT 1.e-5) ; ---- resolution parameter ---- npix = (size(data))[1] if (keyword_set(quadcube)) then begin pix_param = ROUND(ALOG(npix/6.)/ALOG(4.)) + 1 ; resolution parameter endif else begin nside = npix2nside(npix) pix_param = nside endelse return end