pro rdhmifits,directory ; procedure input the inversion results from HMI, save ; in a common area, including UT ; INPUTS: directory = AZAM data directory ; common area for outputs common com_invertresults, $ fld, psi, azm, eta0, dmp, dop, bzero, b1mu, cen1, alpha, $ cct, chisq, occt, fld_er, psi_er, azm_er, $ cen1_er, alpha_er, xx, yy, slat, slong, $ yr, mnth, dy, hr, minit, sec, cdelt1, cdelt2,rconv ; HMI inversion data ; establish relevant part of file name for each parameter strhmiparams = [ $ 'alpha_mag.fits', 'azimuth.fits', 'chisq.fits', 'damping.fits', $ 'dop_width.fits', 'eta_0.fits', 'field.fits', 'inclination.fits', $ 'src_continuum.fits', 'src_grad.fits', 'vlos_mag.fits', $ '.alpha_err.fits','.azimuth_err.fits','.field_err.fits', $ '.inclination_err.fits','.vlos_err.fits' $ ] strvarname = [ $ 'alpha', 'azm', 'chisq', 'dmp', $ 'dop', 'eta0', 'fld', 'psi', $ 'bzero', 'b1mu', 'cen1', $ 'alpha_er','azm_er','fld_er', $ 'psi_er','cen1_er' $ ] ; number of inversion parameters ninvpar = n_elements(strhmiparams) ; Append directory path with '/'. if(strmid(directory, strlen(directory)-1, 1) ne '/') then begin directory = directory + '/' endif ; loop over parameters, reading in the FITS files for kk = 0,ninvpar-1 do begin derf = file_search(directory + '*' + strhmiparams(kk)) filnam = derf[0] ; tdat = readfits(filnam,hdr) ; use new IDL HMI read package read_sdo,derf,hdr,tdat ; eliminate NaNs, changing them to zero values whrnan = where(tdat ne tdat,ncount) if (ncount ne 0) then tdat(whrnan) = 0. ; assign names to the variable void = execute(strvarname[kk] + '= tdat') ; for first variable, assign RCONV, unity for fitted points if kk eq 0 then rconv = alpha endfor ; read in the continuum intensity from observed data derfi0 = file_search(directory + 'continuum.fits') ;cont0 = fitsio_read_image(string(derfi0(0)),hdrc0) ;cont0 = fitsio_read_image(derfi0(0),hdrc0) ; use new IDL HMI read package read_sdo,derfi0,hdrc0,cont0 whr = where(cont0 ne cont0) cont0(whr) = 0. occt = float(cont0) cct = occt ; constrain field > 0 fld = fld>0. ; rotate azimuth angle so that zero azimuth is solar west. ; The inversion code places zero azimuth along pixel column, ; increasing counter-clockwise azm = azm + 90. ; save HMI FITS header in the AZAM directory as an ; IDL save-file ; ???NOTE: read_sdo saves header as a structure, so sxpar does not work ; below we extract the desired header values from the structure instead) ; Therefore, header is not saved here ;header = hdr ;save,filename=directory+'HMI_FITSheader.save',header ; get the time information from the FITS file header ;timez = sxpar(hdr,'T_OBS') timez = hdr.T_OBS yr=fix(strmid(timez,0,4)) mnth=fix(strmid(timez,5,2)) dy=fix(strmid(timez,8,2)) hr=fix(strmid(timez,11,2)) minit=fix(strmid(timez,14,2)) sec=float(strmid(timez,17,6)) ; find the sun center point from the header values ; note: first pixel is 1, not 0! Remove this offset for IDL ;xctr = sxpar(hdr,'CRPIX1') - 1. ;yctr = sxpar(hdr,'CRPIX2') - 1. xctr = hdr.CRPIX1 - 1. yctr = hdr.CRPIX2 - 1. ; find the rotation parameter, presume CLOCKWISE from solar N, degrees ;angrot = sxpar(hdr,'CROTA2') angrot = hdr.CROTA2 ; if the rotation parameter is close to 180 degrees, then the HMI ; images will be rotated by exactly 180 degrees because we will eventually ; need to specify the azimuth for the original data presentation, not ; the rotated one. Must do the rotation otherwise the field angles ; will not be presented to the user correctly if(abs(angrot-180.) le 4.0) then begin for kk=0,ninvpar-1 do begin void = execute(strvarname[kk] + '= rotate(' + strvarname[kk] + ',2)' ) endfor endif else begin print,'CROTA2 = ',angrot,' > 4 deg. Cannot display in AZAM correctly' stop endelse ; rotate the continuum image occt = rotate(occt,2) cct = occt ; with 180 deg rotation, now must find new sun center positions ; get some dimensions npix = n_elements(fld) sz = size(fld) nxx = sz(1) nyy = sz(2) xctr = float(nxx)-1-xctr yctr = float(nyy)-1-yctr ; get arrays for x,y positions, arcsec ;cdelt1 = sxpar(hdr,'CDELT1') cdelt1 = hdr.CDELT1 xx = fld & xx(*,*) = 0. ;cdelt2 = sxpar(hdr,'CDELT2') cdelt2 = hdr.CDELT2 for ii = 0,nxx-1 do xx(ii,*) = (float(ii) - xctr)*cdelt1 yy = fld & yy(*,*) = 0. for jj = 0,nyy-1 do yy(*,jj) = (float(jj) - yctr)*cdelt2 ; find the solar latitude,longitude from XX, YY utimed = double(float(hr) + float(minit)/60. + sec/3600.) ; expand time arrays to full array size, not just per slit position iyr = intarr(nxx,nyy) & imnth = intarr(nxx,nyy) & idy = intarr(nxx,nyy) ihr = intarr(nxx,nyy) & iminit = intarr(nxx,nyy) secc = fltarr(nxx,nyy) & uutimed =fltarr(nxx,nyy) iyr(*,*) = yr imnth(*,*) = mnth idy(*,*) = dy ihr(*,*) = hr iminit(*,*) = minit secc(*,*) = sec uutimed = utimed azam_sunframe, iyr, imnth, idy, uutimed, xx, yy, $ slong, slat, valuemu, $ ras, dcs, gsdt, b0ang, peeang, ctrlong, r0 ; convert longitude,latitude from radians to degrees slong = float(slong/!dtor) slat = float(slat/!dtor) return end