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 ; bypass reading data etc, already saved ;goto,jump44 ; 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 tdat = readfits(filnam,hdr) ; eliminate NaNs, changing them to zero values whr = where(tdat ne tdat,ncount) if (ncount ne 0) then tdat(whr) = 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 ;; !!!!temporarily use the computed continuum for continuum ;cct = bzero+b1mu ;occt = cct ;; read in the extreme raw intensity images, average them ;; search data directory for I0 and I5 ;goto, jump55 ;derfi0 = file_search(directory + '*' + 'I0' + '*' + '.fits') ;; read in I0 continuum intensity ;cont0 = fitsio_read_image(derfi0(0),hdrc0) ;whr = where(cont0 ne cont0) ;cont0(whr) = 0. ;; read in I5 continuum intensity ;derfi5 = file_search(directory + '*' + 'I5' + '*' + '.fits') ;cont5 = fitsio_read_image(derfi5(0),hdrc5) ;whr = where(cont5 ne cont5) ;cont5(whr) = 0. ;; average the images ;occt = float(0.5*(cont0 + cont5)) ;jump55: ; read in the continuum intensity from observed data derfi0 = file_search(directory + 'continuum.fits') cont0 = fitsio_read_image(derfi0(0),hdrc0) 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 header = hdr save,filename=directory+'HMI_FITSheader.save',header ; get the time information from the FITS file header timez = sxpar(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. ; find the rotation parameter, presume CLOCKWISE from solar N, degrees angrot = sxpar(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') xx = fld & xx(*,*) = 0. cdelt2 = sxpar(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 ; xy2lonlat calculates its own ephemeris values based on times ; make a loop to try to speed up this really slow calculation slong = fltarr(nxx,nyy) & slat = fltarr(nxx,nyy) ; process by column arcvec = fltarr(2,nxx) time1d = strarr(nxx) time1d(*) = timez print,' NOTICE: beginning computation of solar latitude,longitude.' print,' The SolarSoft routine xy2lonlat is VERY SLOW.....have patience!' for jj = 0,nyy-1 do begin ; calculate modulus, and print if 100 jjmod = jj mod 100 if jjmod eq 0 then print,'starting lat,long computation for row ',jj arcvec(0,*) = xx(*,jj) arcvec(1,*) = yy(*,jj) lonlat = xy2lonlat(arcvec,time1d) slong(*,jj) = lonlat(0,*) slat(*,jj) = lonlat(1,*) endfor ;save,filename='rdhmifits.save', $ ; 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 ;jump44: restore,'rdhmifits.save' return end