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 ; select the first image of cont0 cont0 = reform(cont0(*,*,0)) ; eliminate NaNs, changing them to zero values whrnan = where(cont0 ne cont0,ncount) if (ncount ne 0) then cont0(whrnan) = 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. whr = where(azm gt 360.,kntazm) if kntazm ne 0 then azm(whr) = azm(whr) - 360. ; 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 = 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 = hdr.CRPIX1 - 1. yctr = hdr.CRPIX2 - 1. ; if these data are SHARP cutouts (active region cutouts rather ; than full-disk data) then CRPIX1,2 are offsets from the disk center ; position of the full disk data, the latter being given in the ; IMCRPIX1,2 header variables. Check for SHARP data by examining the ; size of the arrays (4096x4096 for full disk data) sz = size(fld) sharps = 0 if(sz(1) lt 4096 or sz(2) lt 4096) then begin xdim = sz(1) & ydim = sz(2) fdxctr = hdr.IMCRPIX1 - 1. fdyctr = hdr.IMCRPIX2 - 1. sharps = 1 ; for SHARPS, it seems that the azimuth is off by 180 consistently azm = azm - 180. whr = where(azm lt 0.,kntazm) if kntazm ne 0 then azm(whr) = azm(whr) + 360. endif ; find the rotation parameter, presume CLOCKWISE from solar N, degrees 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) cdelt1 = hdr.CDELT1 cdelt2 = hdr.CDELT2 xx = fld & xx(*,*) = 0. yy = fld & yy(*,*) = 0. ; for full disk data if sharps eq 0 then begin xctr = float(nxx)-1-xctr yctr = float(nyy)-1-yctr ; get arrays for x,y positions, arcsec for ii = 0,nxx-1 do xx(ii,*) = (float(ii) - xctr)*cdelt1 for jj = 0,nyy-1 do yy(*,jj) = (float(jj) - yctr)*cdelt2 endif else begin ; for SHARPS cutouts ;print,'sharps dimension: ',nxx,nyy,', dist from ctr: ',xctr,yctr ;stop for ii = 0,nxx-1 do xx(ii,*) = (float(ii) + xctr - nxx + 1)*cdelt1 for jj = 0,nyy-1 do yy(*,jj) = (float(jj) + yctr - nyy + 1)*cdelt2 endelse ; 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